This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Varying assembly statistics: Abyss & SPades

Hi guys!

I am trying to assemble a genome of ~5MB and I have 250bp PE reads. I tried assembling with Abyss:

abyss-pe k=64 name=novo in='r1.fastq r2.fastq'

This resulted in a scaffold N50 of 50k and total assembly size of 5MB (and if you run with k=128, N50=60k and total assembly size=4.5MB). I wanted to see if I can improve this using SPades so I ran it like this:

spades.py -1 r1.fastq -2 r2.fastq --careful -k 21,33,55,77,99,127 -o spades_assembly

And then used QUAST to get assembly stats like so:

quast-5.0.2/quast.py scaffolds.fasta -o report

The surprising result is that the resulting stats are a lot worse (low N50 of 1455, very high total length of 18 million) and I have to think that something went wrong or maybe I am missing something. --careful flag runs error correcting which is not done in Abyss but I don't think this is the reason? Full output of QUAST is below:

Assembly                    scaffolds
# contigs (>= 0 bp)         32842    
# contigs (>= 1000 bp)      2918     
# contigs (>= 5000 bp)      170      
# contigs (>= 10000 bp)     85       
# contigs (>= 25000 bp)     50       
# contigs (>= 50000 bp)     34       
Total length (>= 0 bp)      25940092 
Total length (>= 1000 bp)   11070707 
Total length (>= 5000 bp)   6648867  
Total length (>= 10000 bp)  6094453  
Total length (>= 25000 bp)  5557460  
Total length (>= 50000 bp)  4920491  
# contigs                   14063    
Largest contig              446405   
Total length                18374693 
GC (%)                      42.66    
N50                         1455     
N75                         708      
L50                         1327     
L75                         6199     
# N's per 100 kbp           84.79

Thanks for any input!

assembly abyss spades
• 476 views • link
updated by h.mon h.mon • written by max_19 max_19

2 answers

Spades often has trouble with 250bp PE reads. I wrote Shovill to resolve this problem: https://github.com/tseemann/shovill In your case I think however you have 1000s of TINY contigs. You need to remove all the small contigs first to fix the N50 statistic. seqtk seq -L 500 contigs.fa > good_contigs.fa

• 0 views • link
written by Torst Torst

Thanks for your suggestion! I tried doing that and then running good_contigs.fa with QUAST, however the N50 did not change. Do you think I have to reassemble with SPades (i.e. remove all the reads which map to those contigs, then reassemble)

• 0 views • link
written by max_19 max_19

If you look at the total assembly size of contigs greater than 50k, you will see you have approximately the same assembly size as the abyss assemblies:

Total length (>= 50000 bp) 4920491

As Torst said, you have a lot of short contigs (SPAdes doesn't have a minimum contig length filter), probably caused by not so good quality towards the end of the reads. If you remove these contigs, SPAdes N50 will be a lot better than the abyss N50. What is the expected assembled genome size? Did you trim the adapters from the reads? Did you apply some read quality filtering? What is the sequencing coverage?

P.S.: you know you can feed quast more than one assembly simultaneously? Something like:

quast-5.0.2/quast.py -o report -l "SPAdes,ABYSS" abyss/contigs.fasta spades/scaffolds.fasta

Then you will have a nice side-by-side comparison

• 0 views • link
written by h.mon h.mon

That is good to know, thank you! Expected assembly size is ~4.8MB. Adapters have been trimmed and I applied a quality filtering of q=30 with trimmomatic. Sequencing coverage is very high (~1000X).

Once I remove those scaffolds that are very small from scaffolds.fasta, should I then remove the reads mapping to those scaffolds, AND then reassemble?

• 0 views • link
written by max_19 max_19

I suspect 1000x is too high coverage, you should reduce it with digital normalization or just plain down-sampling. After that, assemble again with either SPAdes or shovil, as suggested above.

Your approach may work, but seems too complicated and I am not sure it wouldn't introduce artifacts. For example, pieces of the original genome may be unrepresented in the larger contigs set. Also, better than filtering just based on just contig length, using contig length and coverage would be more efficient.

• 0 views • link
written by h.mon h.mon

I actually tried subsampling (Randomly select 15% of the data) and I assembled with SPades, see stats below:

Assembly                    scaffolds
# contigs (>= 0 bp)         12763    
# contigs (>= 1000 bp)      1672     
# contigs (>= 5000 bp)      87       
# contigs (>= 10000 bp)     39       
# contigs (>= 25000 bp)     20       
# contigs (>= 50000 bp)     14       
Total length (>= 0 bp)      10596879 
Total length (>= 1000 bp)   5513812  
Total length (>= 5000 bp)   2823624  
Total length (>= 10000 bp)  2520546  
Total length (>= 25000 bp)  2246504  
Total length (>= 50000 bp)  2025521  
# contigs                   5166     
Largest contig              278079   
Total length                7872888  
GC (%)                      38.46    
N50                         1939     
N75                         873      
L50                         492      
L75                         2092     
# N's per 100 kbp           112.00

How else would you suggest going about filtering the contings? (Sorry if this is a basic question)

• 0 views • link
written by max_19 max_19

Rather than plain down- or sub-sampling you should normalize the data. You could do that using BBNorm from BBMap suite. A guide to use bbnorm is available here.

• 0 views • link
written by GenoMax GenoMax

Will try that - thank you!

• 0 views • link
written by max_19 max_19

Log in to answer this question.

More posts like this