Answers.Body: <p>This could be organism-specific. We don't have a lot of info so far, so I would check a few more things:</p>
  <p>Run something like <a href="https://github.com/vezzi/FRC_align" rel="nofollow noreferrer">FRC_align</a>. Check if there's a clear signal between regions flagged as suspicious by it and your coverage graph.</p>
  <p>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.</p>
  <p>Worried about contamination? Run <a href="https://ccb.jhu.edu/software/kraken2/" rel="nofollow noreferrer">kraken2</a> on your raw data and see what you get there.</p>
  <p>Run <a href="http://www.repeatmasker.org/" rel="nofollow noreferrer">RepeatMasker</a> and check if repetive regions contribute to the effect you see.</p>
  <p>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 <a href="http://mummer.sourceforge.net/" rel="nofollow noreferrer">MUMmer</a> 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.</p>
  <p>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:
  <a href="https://i.sstatic.net/APLhU.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/APLhU.png" alt="GC distributions" /></a></p>
  <p>What does GC vs Length look like? Here's the one from the 45 genomes mentioned above:
  <a href="https://i.sstatic.net/Yme7I.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/Yme7I.png" alt="GC vs Length" /></a></p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-11T08:33:12.090Z
Answers.Text: Thank you for the tips I will try them and let you know. It is a genome from a male Hymenoptera (so haploid genome).
Answers.UserDisplayName: chippycentra
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-11T17:54:11.603Z
Answers.Text: 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?
Answers.UserDisplayName: chippycentra
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-11T18:29:45.240Z
Answers.Text: 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](https://hymenoptera.elsiklab.missouri.edu/) genomes.
Answers.UserDisplayName: Bastian Schiffthaler
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-12T09:57:40.313Z
Answers.Text: 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 ?
Answers.UserDisplayName: chippycentra
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-12T10:13:02.210Z
Answers.Text: 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.
Answers.UserDisplayName: Bastian Schiffthaler
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-12T10:59:45.047Z
Answers.Text: Yes sure: here is a room where we can talk  https://chat.stackoverflow.com/rooms/224463/chippycentra
Answers.UserDisplayName: chippycentra
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-11T08:14:21.997Z
Answers.LastActivityDate: 2020-11-11T18:47:37.073Z
Answers.LastEditDate: 2020-11-11T18:47:37.073Z
Answers.UserDisplayName: Bastian Schiffthaler
Answers.Body: <p>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 <a href="https://github.com/DRL/blobtools" rel="nofollow noreferrer">blobtools</a>, which creates a similar plot to the one you made manually, but also adds taxonomical annotations. It can look like this:</p>
  <p><a href="https://i.sstatic.net/wDDMW.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/wDDMW.jpg" alt="enter image description here" /></a></p>
  <p>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.</p>
  <p>-- edit --</p>
  <p>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 <a href="https://www.nature.com/articles/s41467-020-14998-3" rel="nofollow noreferrer">kind of problems in data</a>.</p>
  <p>-- edit 2 --</p>
  <p>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.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-13T09:28:46.267Z
Answers.Text: 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.
Answers.UserDisplayName: chippycentra
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-13T10:55:30.537Z
Answers.Text: That sounds like a lot. Are the scaffolds long? What is the sum of the 283 812 scaffolds?
Answers.UserDisplayName: Kamil S Jaron
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-13T11:01:27.327Z
Answers.Text: I also added one more thought to the answer. Is it possible you deal with a polyploid?
Answers.UserDisplayName: Kamil S Jaron
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-13T11:20:42.553Z
Answers.Text: 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)
Answers.UserDisplayName: juanjo75es
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-13T11:32:23.337Z
Answers.Text: 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.
Answers.UserDisplayName: Kamil S Jaron
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-13T13:18:06.820Z
Answers.Text: I added summary about assembly metrics on the post. And I'm sure it is an haploïd genome since it comes from an Hymenoptera male. I tought it would be contamination but as I said, even if the genome is quite fragmented it still found 82% BUSCO genes. And when I blast the BUSCO in both high and low (two pics) coverage scaffolds, I found an Hymenoptera hit...
Answers.UserDisplayName: chippycentra
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-13T14:56:50.017Z
Answers.Text: Just wanted to show my skepticism on the suitability of methods for assembly evaluation based on evaluating the genes obtained or in the GC percent. Especially with non-robust data. I am not sure how these methods work but seems to me that there are many requirements that could be missed leading to misleading data. But that's an interesting thread and I'm interested to learn.  Please go ahead...
Answers.UserDisplayName: juanjo75es
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-14T14:09:03.623Z
Answers.Text: I did the assembly with another assembler (MEGAHIT) and it seems that it does not come from here since I see the same coverage patterns (see figure in edit post)...
Answers.UserDisplayName: chippycentra
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-15T18:27:56.657Z
Answers.Text: Ha, having a hymenopteran male is really useful for ploidy doubts :-). Did you use a single male for the library prep? I thinking if you could have contamination of different lab strains? Is the BUSCO evenly distributed between high and low coverage scaffolds?
Answers.UserDisplayName: Kamil S Jaron
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-12T20:58:35.833Z
Answers.LastActivityDate: 2020-11-15T18:39:00.427Z
Answers.LastEditDate: 2020-11-15T18:39:00.427Z
Answers.UserDisplayName: Kamil S Jaron
Body: <p>I had a question for the community.</p>
  <p>I have a genome of a new species that has been sequenced via 150pb Illumina paired-end.</p>
  <p>To verify the quality of the assembly I used the <code>BUSCO pipeline</code> which allows me to search for 99% of the BUSCO genes (gene not duplicated) present in a taxon in the genome of my species.</p>
  <p>Globally I only miss <code>18%</code> of the BUSCO genes, which means that the assembly is not too bad.</p>
  <p>Now when I plot the <strong>G+C%</strong> and the <strong>coverage</strong> of all the scaffolds containing BUSCO genes (so a priori non-contaminating) I should find a fairly homogeneous cloud.</p>
  <p>Instead I end up with two different coverage profiles (noted red and blue below):</p>
  <p><a href="https://ibb.co/6X34kfq" rel="nofollow noreferrer"></a></p>
  <p>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...</p>
  <p>Does anyone have any idea what these are?</p>
  <p>here is an output of the <code>contigsTable.csv</code> where the two first row are scaffold with low coverage in the figure, and two last are scaffold with high coverage in the figure</p>
  <p><a href="https://i.sstatic.net/Dm85F.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/Dm85F.png" alt="enter image description here" /></a></p>
  <p>FRC output</p>
  <p><a href="https://i.sstatic.net/Z9pTV.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/Z9pTV.png" alt="enter image description here" /></a></p>
  <p>Coverage distribution</p>
  <p><a href="https://i.sstatic.net/jDU6S.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/jDU6S.png" alt="enter image description here" /></a></p>
  <p>Here is an Icarus summary of the assembly :</p>
  <pre><code>Assembly
  # contigs (&gt;= 0 bp) 1038631
  # contigs (&gt;= 1000 bp) 120132
  # contigs (&gt;= 5000 bp) 22245
  # contigs (&gt;= 10000 bp) 9986
  # contigs (&gt;= 25000 bp) 1322
  # contigs (&gt;= 50000 bp) 97
  Total length (&gt;= 0 bp) 655530043
  Total length (&gt;= 1000 bp) 451497861
  Total length (&gt;= 5000 bp) 260894112
  Total length (&gt;= 10000 bp) 174183054
  Total length (&gt;= 25000 bp) 45006012
  Total length (&gt;= 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
  </code></pre>
  <p>Here is the BUSCO summary as well :</p>
  <pre><code>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
  </code></pre>
  <p>EDIT</p>
  <p>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 :</p>
  <p><a href="https://i.sstatic.net/6cq6m.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/6cq6m.png" alt="enter image description here" /></a></p>
  
Tags: <assembly><read-mapping><coverage>
Title: Contamination on genome assembly
ThreadText: # Contamination on genome assembly
  
  (Asked by: username_0 on 2020-11-10)
  
  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
  
  [![enter image description here](https://i.sstatic.net/Dm85F.png)](https://i.sstatic.net/Dm85F.png)
  
  FRC output
  
  [![enter image description here](https://i.sstatic.net/Z9pTV.png)](https://i.sstatic.net/Z9pTV.png)
  
  Coverage distribution
  
  [![enter image description here](https://i.sstatic.net/jDU6S.png)](https://i.sstatic.net/jDU6S.png)
  
  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 :
  
  [![enter image description here](https://i.sstatic.net/6cq6m.png)](https://i.sstatic.net/6cq6m.png)
  
  **Tags:** assembly, read-mapping, coverage
  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_1 on 2020-11-11. Score: 5
  
  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](https://github.com/vezzi/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](https://ccb.jhu.edu/software/kraken2/) on your raw data and see what you get there.
  
  Run [RepeatMasker](http://www.repeatmasker.org/) 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](http://mummer.sourceforge.net/) 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:
  [![GC distributions](https://i.sstatic.net/APLhU.png)](https://i.sstatic.net/APLhU.png)
  
  What does GC vs Length look like? Here's the one from the 45 genomes mentioned above:
  [![GC vs Length](https://i.sstatic.net/Yme7I.png)](https://i.sstatic.net/Yme7I.png)
  
  > **Comments:**
  > **username_0**: Thank you for the tips I will try them and let you know. It is a genome from a male Hymenoptera (so haploid genome).  
  > **username_0**: 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?  
  > **username_1**: 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](https://hymenoptera.elsiklab.missouri.edu/) genomes.  
  > **username_0**: 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 ?  
  > **username_1**: 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.  
  > **username_0**: Yes sure: here is a room where we can talk  https://chat.stackoverflow.com/rooms/224463/username_0  
  
  ---
  
  ### Answer by username_2 on 2020-11-12. Score: 0
  
  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](https://github.com/DRL/blobtools), which creates a similar plot to the one you made manually, but also adds taxonomical annotations. It can look like this:
  
  [![enter image description here](https://i.sstatic.net/wDDMW.jpg)](https://i.sstatic.net/wDDMW.jpg)
  
  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](https://www.nature.com/articles/s41467-020-14998-3).
  
  -- 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.
  
  > **Comments:**
  > **username_0**: 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.  
  > **username_2**: That sounds like a lot. Are the scaffolds long? What is the sum of the 283 812 scaffolds?  
  > **username_2**: I also added one more thought to the answer. Is it possible you deal with a polyploid?  
  > **username_3**: 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)  
  > **username_2**: 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.  
  > **username_0**: I added summary about assembly metrics on the post. And I'm sure it is an haploïd genome since it comes from an Hymenoptera male. I tought it would be contamination but as I said, even if the genome is quite fragmented it still found 82% BUSCO genes. And when I blast the BUSCO in both high and low (two pics) coverage scaffolds, I found an Hymenoptera hit...  
  > **username_3**: Just wanted to show my skepticism on the suitability of methods for assembly evaluation based on evaluating the genes obtained or in the GC percent. Especially with non-robust data. I am not sure how these methods work but seems to me that there are many requirements that could be missed leading to misleading data. But that's an interesting thread and I'm interested to learn.  Please go ahead...  
  > **username_0**: I did the assembly with another assembler (MEGAHIT) and it seems that it does not come from here since I see the same coverage patterns (see figure in edit post)...  
  > **username_2**: Ha, having a hymenopteran male is really useful for ploidy doubts :-). Did you use a single male for the library prep? I thinking if you could have contamination of different lab strains? Is the BUSCO evenly distributed between high and low coverage scaffolds?  
  
  ---
  
  