test

Test dataset

This directory contains data from simulated hybrid genome with 5% divergence between parentals and 40% of loss of heterozygosity (LOH).

To run the test example, execute:

./redundans.py -v -i test/*.fq.gz -f test/contigs.fa -o test/run1 

Table of contents

• Generation of test set
  • Simulations
    • Heterozygous genome
    • Short reads
  • De novo genome assembly
  • Redundans pipeline
    • Run statistics
      • De novo assembly
      • Parameters estimation
      • Reduction
      • Scaffolding
      • Gap closing
      • Summary statistics
• Accuracy estimation

Generation of test set

Below, you can find detailed description for regeneration of test dataset with any input genome.

Simulations

Test dataset is based on 100 Kb region from C. parapsilosis CDC317 [HE605202:100000-200000].

samtools faidx CANPA.fa HE605202:100000-200000 > ref.fa

Heterozygous genome

Heterozygous regions, 5% divergence and 40% LOH, have been simulated with fasta2diverged.py.

fasta2diverged.py -v -i ref.fa --loh 0.40 --dna -d 0.05 > ref.d05.l40.fa

Short reads

Two shotgun sequencing libraries including typical Illumina-related errors were simulated using GemSIM:

• 100 bp paired-end reads with 600 bp insert [200X coverage]:
• 600_1.fq.gz
• 600_2.fq.gz
• 50 bp mate-pair reads with 5 kb insert [20X coverage]:
• 5000_1.fq.gz
• 5000_2.fq.gz

# generate SNPs
~/src/GemSIM_v1.6/GemHaps.py -r ref.fa -g '.50,0 .50,6'

# 2x 50,000 reads
## 100 bp paired-end with 600 bp insert +/- 50 bp
~/src/GemSIM_v1.6/GemReads.py -r ref.fa -n 50000 -g ref.txt -l 100 -m /home/lpryszcz/src/GemSIM_v1.6/models/ill100v5_p.gzip -q 33 -o 600 -u 600 -s 50 -p
~/src/GemSIM_v1.6/GemReads.py -r ref.d05.l40.fa -n 50000 -g ref.txt -l 100 -m /home/lpryszcz/src/GemSIM_v1.6/models/ill100v5_p.gzip -q 33 -o 600.d05.l40 -u 600 -s 50 -p

# 2x 10,000 reads
## 50 bp mate-pairs with 5 kb insert +/- 500 bp
~/src/GemSIM_v1.6/GemReads.py -r ref.fa -n 10000 -g ref.txt -l 50 -m /home/lpryszcz/src/GemSIM_v1.6/models/ill100v5_p.gzip -q 33 -o 5000 -u 5000 -s 500 -p
~/src/GemSIM_v1.6/GemReads.py -r ref.d05.l40.fa -n 10000 -g ref.txt -l 50 -m /home/lpryszcz/src/GemSIM_v1.6/models/ill100v5_p.gzip -q 33 -o 5000.d05.l40 -u 5000 -s 500 -p

# combine both libraries in random order
## 600_
paste <(zcat 600_1.fastq.gz 600.d05.l40_1.fastq.gz) <(zcat 600_2.fastq.gz 600.d05.l40_2.fastq.gz) | paste - - - - | shuf | awk -F'\t' '{OFS="\n"; print $1,$3,$5,$7 | "gzip > random.600.d05.l40_1.fq.gz"; print $2,$4,$6,$8 | "gzip > random.600.d05.l40_2.fq.gz"}'
## 5000_
paste <(zcat 5000_1.fastq.gz 600.d05.l40_1.fastq.gz) <(zcat 5000_2.fastq.gz 5000.d05.l40_2.fastq.gz) | paste - - - - | shuf | awk -F'\t' '{OFS="\n"; print $1,$3,$5,$7 | "gzip > random.5000.d05.l40_1.fq.gz"; print $2,$4,$6,$8 | "gzip > random.5000.d05.l40_2.fq.gz"}'

# symlinks with handier names
for f in *.fq.gz; do ln -s $f `echo $f | sed 's/random.//;s/.d05.l40//'`; done

# mate-pairs contaminated with paired-end reads
#for i in 1 2; do zcat <(zcat test/5000_$i.fq.gz test/5000_$i.fq.gz | fastx2reverse_complement.py fastq | gzip) test/600_$i.fq.gz > test/mixed.$i.fq.gz; done

De novo genome assembly

Simulated reads were assembled with SPAdes v.3.5.0.

dipspades.py --only-assembler -t 4 -o 600 -1 600_1.fq.gz -2 600_2.fq.gz

Redundans pipeline

Finally, heterozygous genome assembly pipeline can be applied to simulated heterozygous contigs and paired/mate-pairs reads.

./redundans.py -v -i test/*.fq.gz -f test/contigs.fa -o test/run1 

Run statistics

De novo assembly

This is optional step performed only if no contigs (-f) are given. Redundans uses Platanus for contigs assembly, scaffolding and gap closing. In first order, Redundans uses libraries with reads between 70-150 bp. If no such libraries are available, libraries with other read lenghts are used, but only libraries with cummulative sequence length of at least 2/3 of the largest library will be used. Note, assembly is be performed in single-end mode and only resulting contigs (not scaffolds) are going to be used.
Subsequently, library statistics are estimated and if FR paired-end libraries are detected, Redundans will proceed with scaffolding and gap closing of de novo contigs using Platanus. This shouldn't be mixed with Redundans scaffolding and gap closing modules - those are using different software and set of parameters.

Parameters estimation

At the beginning, Redundans estimates number of parameters:

• number of reads that it's going to align (based on -l/--limit parameter)

 Aligning 32779 mates per library...

• library statistics:
• insert size: median, mean, stdev
• mates orientation: FF, FR, RF, RR

Table 1: Library statistics

FastQ files	read length	median	mean	stdev	FF	FR	RF	RR
test/5000_1.fq.gz test/5000_2.fq.gz	50	4998	4990.20	721.47	0	4,674	0	0
test/600_1.fq.gz test/600_2.fq.gz	100	599	598.63	47.68	0	10,000	0	0

Above, you can see, there are two libraries:

• test/5000_1.fq.gz & test/5000_2.fq.gz:
  • FR orientation
  • mate-pairs with 5000 bp insert size +/- 600 bp
• test/600_1.fq.gz & test/600_2.fq.gz:
  • FR orientation
  • paired-end with 600 bp insert size +/- 39 bp

Based on these statistics, the program will perform scaffolding and gap closing starting from libraries with the smallest insert size.
Also at this stage, the user is notified about libraries with poor statistics i.e. large stdev or not consisted orientation. For these cases, the statistics are updated before every scaffolding iteration. This is done, because for libraries with large insert size, the estimation of insert size may be highly affected by the fragmentation of the assembly.

Reduction

In the first step, Redundans will recognise and selectively remove heterozygous contigs from initial assembly. The details about all removed heterozygous contigs are stored in tab-delimited file contigs.reduced.fa.hetero.tsv. Additionally, the histogram with identities of reduced contigs are saved as contigs.reduced.fa.hist.png.

At the end of this step, the user will be notified about:

• initial assembly size and fragmentation
• what portion of the the initial assembly have been recognised as heterozygous
• identity between heterozygous regions
• size and fragmentation of reduced assembly

Table 2: Reduction statistics

file name	genome size	contigs	heterozygous size	[%]	heterozygous contigs	[%]	identity [%]	possible joins	homozygous size	[%]	homozygous contigs	[%]
test/run6/contigs.fa	163897	245	76363	46.59	222	90.61	94.837	0	97319	59.38	23	9.39

Above, you can see that:

• initial SPAdes assembly was fragmented (245 contigs) and larger (163.8 Kb) than reference sequence (100 Kb)
• Redundans:
• removed 220 heterozygous contigs, resulting in reduced assembly of 97.8 Kb (very similar to reference) in 28 contigs
• estimated 94.837% identity between heterozygous regions, which is very similar to the divergence that was used for simulations (5%)

Scaffolding

Redundans repeat scaffolding on each library several times (--iters). For each scaffolding iteration, the program reports number of:

• processed pairs
• pairs passing alignment quality
• pairs that aligned in different contigs

 iteration 1.1 ...
   19572 pairs. 18505 passed filtering [94.55%]. 1993 in different contigs [10.18%].
 iteration 1.2 ...
   19572 pairs. 18569 passed filtering [94.88%]. 340 in different contigs [1.74%].
 iteration 2.1 ...
   19572 pairs. 17009 passed filtering [86.90%]. 755 in different contigs [3.86%].
  closing gaps ...
 iteration 2.2 ...
   19572 pairs. 17261 passed filtering [88.19%]. 0 in different contigs [0.00%].
  closing gaps ...

Above, you can see that in iteration 1.1 above 10% of reads aligned in different contigs, while in subsequent iteration, less than 2% reads. This informs that in general it's enough to use 1-2 iterations of scaffolding per library.

Gap closing

For the time-being, only the iteration number is reported.

 iteration 1.1 ...
 iteration 1.2 ...
 iteration 2.1 ...
 iteration 2.2 ...

Summary statistics

At the end, the program reports details of each step:

• file name
• fragmentation and size of the assembly at this stage
• %GC content
• fragmentation and size of assembly in contigs >1 Kb
• N50 & N90
• cumulative size of gaps
• the size of the longest contig

fname	contigs	bases	GC [%]	contigs >1kb	bases in contigs >1kb	N50	N90	Ns	longest
test/run6/contigs.fa	245	163897	40.298	24	117391	3975	233	0	29603
test/run6/contigs.reduced.fa	23	97319	39.333	17	94157	7321	2195	0	29603
test/run6/_sspace.1.1.fa	4	97242	39.315	3	96924	87727	87727	473	87727
test/run6/_sspace.1.2.fa	3	97321	39.315	2	97003	87727	87727	552	87727
test/run6/_sspace.2.1.fa	1	100280	39.315	1	100280	100280	100280	3511	100280
test/run6/_sspace.2.1.filled.fa	1	100251	39.363	1	100251	100251	100251	2782	100251
test/run6/_sspace.2.2.fa	1	100251	39.363	1	100251	100251	100251	2782	100251
test/run6/_sspace.2.2.filled.fa	1	100123	39.360	1	100123	100123	100123	2677	100123
test/run6/scaffolds.fa	1	100123	39.360	1	100123	100123	100123	2677	100123
test/run6/_gapcloser.1.1.fa	1	100495	39.542	1	100495	100495	100495	1265	100495
test/run6/_gapcloser.1.2.fa	1	100498	39.579	1	100498	100498	100498	975	100498
test/run6/_gapcloser.2.1.fa	1	100502	39.679	1	100502	100502	100502	4	100502
test/run6/_gapcloser.2.2.fa	1	100506	39.679	1	100506	100506	100506	4	100506
test/run6/scaffolds.filled.fa	1	100506	39.679	1	100506	100506	100506	4	100506

Accuracy estimation

Accuracy of recovered contigs can be assessed by alignment of final scaffolds (run1/scaffolds.filled.fa) back onto reference (ref.fa).

NUCMER

This can be quickly achieved with nucmer (works for small genomes):

cd run1
# align
nucmer -p nucmer.ref ../ref.fa scaffolds.filled.fa
# generate pairwise plot
mummerplot --png --large nucmer.ref.delta -p nucmer.ref.plot
# open plot
eog nucmer.ref.plot.png

Figure 1: Pairwise alignment of reference sequence and redundans scaffolds

You may want to compare final redundans scaffolds (run1/scaffolds.filled.fa) with initial SPAdes contigs (contigs.fa).

cd run1
# align
nucmer -p nucmer.contigs scaffolds.filled.fa ../contigs.fa
# generate pairwise plot
mummerplot --png --large nucmer.contigs.delta -p nucmer.contigs.plot
# open plot
eog nucmer.contigs.plot.png

Figure 2: Pairwise alignment of redundans scaffolds and SPAdes contigs

Finally, comparison of consensus contigs from dipSPAdes (consensus_contigs.fa) and reference sequnece (ref.fa) informs about missing regions in dipSPAdes reconstruction.

# align
nucmer -p nucmer.ref_dipspades ref.fa consensus_contigs.fa
# generate pairwise plot
mummerplot --png --large nucmer.ref_dipspades.delta -p nucmer.ref_dipspades.plot
# open plot
eog nucmer.ref_dipspades.plot.png

Figure 3: Pairwise alignment of reference sequence and dipSPAdes consensus contigs

LAST

The same can be done using LAST (this will work also for large genomes):

# index reference genome (need to be done only once per each reference)
lastdb ../ref.fa ../ref.fa

# align & generate dotplot on the fly
lastal -f TAB ../ref.fa scaffolds.filled.fa | last-dotplot - scaffolds.filled.fa.png

# open plot
eog scaffolds.filled.fa.png

Figure 4: Pairwise alignment of reference sequence and Redundans contigs

All programs/scripts not disclosed in this repository (i.e. fasta2diverged.py), can be found in https://github.com/lpryszcz/bin.