Hello everyone,
I've been assembling yeast genomes for the past 3 years. My data consist of 301 bp paired-end reads with an average coverage of ~350×. The genomes belong to novel yeast species and range from 12–19 Mb, depending on the genus. For this particular genome, I have approximately 16 million reads available (8M forward + 8M reverse).
I usually perform de novo assemblies, and although I have tested several assemblers (ABySS, MaSuRCA, SOAPdenovo, among others), SPAdes consistently produces the best results. For most genomes, I obtain assemblies with N50 values >200 kb, although assemblies with >1000 contigs/scaffolds are still common.
A few months ago, however, I encountered a genome that seems impossible to assemble properly.
The reads have good quality, and I additionally trimmed them with fastp. I ran GenomeScope, Jellyfish, and KmerGenie to investigate genome characteristics and optimize k-mer selection. I also tested more than 14 different SPAdes parameter combinations, in addition to trying other assemblers, but none produced a satisfactory assembly (I get extremelly big gensizes and super small N50).
One approach that partially "worked" was aligning the reads against a closely related species that I had previously assembled (>1000 contigs, N50 ~130 kb) and then performing a de novo assembly using the --trusted-contigs option in SPAdes. Even then, the result was poor (>3000 contigs and N50 of only 12 kb), although BUSCO completeness was >90%.
My current idea is to use the assembly of a closely related species as a reference, map the reads, call variants, and generate a consensus assembly. However, I am concerned that this approach could introduce a strong reference bias. Since I am also performing ANI analyses, I worry that the resulting assembly would artificially inflate similarity to the reference genome. In fact, when I tested this approach, the resulting assembly showed >99% ANI to the reference genome. I am fairly confident this does not reflect the true relationship between these species, as phylogenomic analyses based on BUSCO genes place them much farther apart than other species pairs in my dataset that share only ~97% ANI.
At this point, I would appreciate any suggestions on alternative strategies or analyses that might help explain why this genome is so difficult to assemble. It is hard to accept having this much data and still being unable to obtain an assembly with an N50 above ~90 kb. Even more so beacause of the fasct that its 'sisters' were easly assembled.
I'm happy to provide additional details if needed. Thanks very much!
0 answers
No answers yet.
Log in to answer this question.
You seem to have a total of 1.4 terabases (4.8 billion reads with 301 bp each, assuming this is adding up both R1/R2 reads). With a 19 Mb genome this would be 76000x coverage (if you actually have 4.8 billion reads each for R1/R2 then this coverage becomes 152,000x).
That is simply too much data for a genome of that size. You need to normalize the data down, if you are not actually using 350x to do the assemblies.
Ideally you should get your hands on some long read data, which would be the way to go to solve the genome.
Oops, that was a typo.
I have approximately 8 million reads in each fastq file, totaling around 16 million reads overall. For this particular genome, each fastq file is approximately 5GB in size.
Regarding the long-read data, we are trying to obtain it, but unfortunately those are very expensive in my country. Since I already have high-quality Illumina data, I was hoping to find a way to improve the assembly using the resources currently available to me.
I'll check for the normalization tho. I have seen that even 300x covered can be too much for some assemblers.
Thank you for you answer :D
Even 300x? I would subsample this to circa 60X and retry.
Long reads are likely to be very helpful as others have said. Maybe you can collaborate with another local user in your country who has a Minion established, or send internationally.
Sticking with why the assembly is poor, I guess some of my first instincts would say that what there was something odd about the sample you extracted DNA from...I few things below to maybe try get at an answer.
When you align the reads against the closely related species, what percentage of the reads align?
For the BUSCOs do you get 90% complete and single copy? or are many duplicated?
Are all your other species homozygote and is it possible that you are dealing with a heavily heterozygote species? You would see that with genomescope.
Also you could check out the spades assembly graph to see what that structure looks like.
Did you check for contamination? Kraken can help to see if there are other DNAs messing with the assembly.