# Contamination on genome assembly

I had a question for the community.

I have a genome of a new species that has been sequenced via 150pb Illumina paired-end.

To verify the quality of the assembly I used the BUSCO pipeline which allows me to search for 99% of the BUSCO genes (gene not duplicated) present in a taxon in the genome of my species.

Globally I only miss 18% of the BUSCO genes, which means that the assembly is not too bad.

Now when I plot the G+C% and the coverage of all the scaffolds containing BUSCO genes (so a priori non-contaminating) I should find a fairly homogeneous cloud.

Instead I end up with two different coverage profiles (noted red and blue below):

I then thought that the low coverage sequencing profile around 4x might be a eukaryotic contaminant like a acarians, however when I blast the BUSCOs present on these low coverage scaffolds I get hits that are very close to the genome of my species...

Does anyone have any idea what these are?

here is an output of the contigsTable.csv where the two first row are scaffold with low coverage in the figure, and two last are scaffold with high coverage in the figure

FRC output

Coverage distribution

Here is an Icarus summary of the assembly :

Assembly
# contigs (>= 0 bp) 1038631
# contigs (>= 1000 bp) 120132
# contigs (>= 5000 bp) 22245
# contigs (>= 10000 bp) 9986
# contigs (>= 25000 bp) 1322
# contigs (>= 50000 bp) 97
Total length (>= 0 bp) 655530043
Total length (>= 1000 bp) 451497861
Total length (>= 5000 bp) 260894112
Total length (>= 10000 bp) 174183054
Total length (>= 25000 bp) 45006012
Total length (>= 50000 bp) 5956145
# contigs 214542
Largest contig 120620
Total length 519138286
GC (%) 32.94
N50 5062
N75 1624
L50 21982
L75 71030
# N's per 100 kbp 22.28


Here is the BUSCO summary as well :

2544 Complete BUSCOs (C)
2312 Complete and single-copy BUSCOs (S)
232 Complete and duplicated BUSCOs (D)
1060 Fragmented BUSCOs (F)
811 Missing BUSCOs (M)
4415 Total BUSCO groups searched
C:57.7%[S:52.4%,D:5.3%],F:24.0%,M:18.3%,n:4415


EDIT

I ran the fastq sequences with another aligner (MEGAHIT) to see if it was not because of an assembly issue and I got the same profils coverage see here :

This could be organism-specific. We don't have a lot of info so far, so I would check a few more things:

Run something like FRC_align. Check if there's a clear signal between regions flagged as suspicious by it and your coverage graph.

Is it a eukaryote? Plant? Check where mitchondria and chloroplasts are on the plot. They will have different GC/coverage signals than the rest of your assembly, which is fine, expected even.

Worried about contamination? Run kraken2 on your raw data and see what you get there.

Run RepeatMasker and check if repetive regions contribute to the effect you see.

Finally, unless your specimen is low heterozygosity, haploid, or inbred, you probably have some effect from ploidy. There are numerous ways to check that. A self-self dot plot with MUMmer and looking coverage histograms are a good start. The idea is that you should see two peaks (let's assume diploid). One for the split haplotigs, and one for collapsed.

I checked the GC distribution of 45 genomes in your order. It's kind of all over the place, but there are some bimodal ones:

What does GC vs Length look like? Here's the one from the 45 genomes mentioned above:

• Thank you for the tips I will try them and let you know. It is a genome from a male Hymenoptera (so haploid genome). – chippycentra Nov 11 '20 at 8:33
• I ran FRC_align but to be honest I do not really understand what I have to look into the FRC_align output (they are many). Is there one that gives you suspicious candidats? – chippycentra Nov 11 '20 at 17:54
• Go by OUTPUT_HEADER_Features.txt. It will have features for abnormal coverage as well as abnormal PE alignments - given your input size. Compression is where they align too close together, stretch is the opposite. I also edited my post with some more info on genome GC densities from these genomes. – Bastian Schiffthaler Nov 11 '20 at 18:29
• Thank you very much for that. I added the outputs I got, I do not have OUTPUT_HEADER_Features.txt... But I'm wondering is the issue does not come from the fact that there are 2 coverage profils rather than GC content ? I also added content of the file contigTable.csv to show you 4 scaffolds, maybe you will se what is going on there ? – chippycentra Nov 12 '20 at 9:57
• The OUTPUT_HEADER part is just a placeholder. In your image, the file is Features.txt. I suggest if there are any more concerns of that nature, we should move this to chat, as comments are discouraged from becoming extended discussions. – Bastian Schiffthaler Nov 12 '20 at 10:13

That is indeed puzzling. My best guess would be contamination, even if they carry some of the BUSCOs. I run all de-novo genomes through blobtools, which creates a similar plot to the one you made manually, but also adds taxonomical annotations. It can look like this:

Also, what is the proportion of the genome that is low-coverage? Does it have a significant span? I almost always see a bunch of tiny contigs with low coverage in de novo assemblies, but they usually represent a very small fraction of the genome.

-- edit --

One more thought, are you sure you don't have a tetraploid species? The two coverage peaks kind of correspond to 8 and 32x right? One is 1/4 of the other. We developed some tools to investigate this kind of problems in data.

-- edit 2 --

Just a few more comments. Coverage differences can have various reasons, but all the sequencing biases would cause a skew of the distribution, hence the apparent bimodality of your coverage (more apparent from the two 2d plots) suggests two distinct sources. In general, could be due to various ploidy levels or contamination either by bacteria, endosymbionts or even other species getting into DNA isolation. Now that polyploidy is ruled out, it boils down to contamination, and to figure out what contamination it is, I recommend the tool I liked above.

• Hello, I have 553 265 scaffolds with coverage >0 in total in which 283 812 have covdepth between [0.1-15] (so more than half of scaffolds), in fact when I plot the coverage I get the plot I added on the post. – chippycentra Nov 13 '20 at 9:28
• That sounds like a lot. Are the scaffolds long? What is the sum of the 283 812 scaffolds? – Kamil S Jaron Nov 13 '20 at 10:55
• I also added one more thought to the answer. Is it possible you deal with a polyploid? – Kamil S Jaron Nov 13 '20 at 11:01
• I don't have much experience but my common sense tells me I would focus on getting larger scaffolds with higher coverage before trying to make fireworks with the content of these scaffolds. If you can't do that for all the genome maybe you should focus on a region and analyze in-depth that region with better data. (addressed to the OP) – juanjo75es Nov 13 '20 at 11:20
• In my experience, sequencing data must make sense. If you have extra coverage peak that affects half of your genome, it's a rather big deal and I do think OP is on right track to investigate where it comes from. – Kamil S Jaron Nov 13 '20 at 11:32