Question: Removing contigs with improbable coverage using R?
0
gravatar for Lesley Sitter
4.5 years ago by
Lesley Sitter490
Netherlands
Lesley Sitter490 wrote:

Hi, 

I have a question involving assemblies and checking the quality of the contigs in those assemblies.

The data i have consists of Pair ended reads, i merged them using PEAR. The reads that were left unassembled by PEAR because they did not have enough overlap were fed to SOAP denovo for assembling. The two assemblies were merged. This gave me a large set of contigs and to check their quality i mapped my reads back to them. The data then looks something like this;

contig_name    contig_coverage
1                     14
2                     9
3                     8
4                     15
5                     14
6                     13
7                     3.000 

The distribution looks something like this

 distribution

On the left is the distribution with the contig that has 3000x coverage, and on the right i used boxplot to find outliers and removed them. It is pretty apparent that contig 7 is an outlier, but the way i removed it is not a proper way because the boxplot method is only useful for a normally distributed dataset and this looks more like log normally distributed.

I have tried using fitdistr() to find the exact type of distribution but it does not produce a definitive answer.

Our statistician also mentioned using generalized linear models, then take the deviance and divide by degrees of freedom to get the distribution. and the closest to 1 is the distribution which i can use to determine the type of outlier detection i need... but i can't get this to work properly (see example down below)

glm1 <- glm(contig_coverage~1,data=data_set,  family = gaussian(log))
glm2 <- glm(contig_coverage~1,data=data_set,  family = poisson(log))
glm3 <- glm(contig_coverage~1,data=data_set,  family = Gamma(log))

 

Can anybody help me with a rule/method/function with which i can find out which contigs have improbable coverage in R so that i can remove them from the final assembly?  

ADD COMMENTlink modified 4.5 years ago by Michael Dondrup47k • written 4.5 years ago by Lesley Sitter490

Some doubts:

  • is it at all a valid approach to make a decision based only on the coverage? you'd remove the contig with the largest amount of evidence in that case.
  • this contig might contain a large stretch of repeats, should you check for that?
  • I guess it's genomic dna and reads are also genomic dna? If not then other concerns would come up.
  • If you don't know the distribution, you could go for non-parametric trimming e.g. trim based on percentile (highest 5% and lowest?)
  • Did you try to fit negative binomial, that seems to be the standard approach today?
ADD REPLYlink modified 4.5 years ago • written 4.5 years ago by Michael Dondrup47k

Thanks for reading my question :D to answer your points

  • The goal is a pipeline that can do a complete assembly, quality filtering, and eventually SNP calling so that the output gets standardized and the laboratory analists don't have to worry about the data handling. We would rather sacrifce a small number of possibly problematic contigs than having to do a manual check each time we want to run this pipeline.
  • I have checked the top ten outliers and they contain contaminations like TE's so it's not valuable in this test and we expect these types of contigs coming up regularly that is why we want to remove them.
  • You are correct, the reads are genomic DNA andour assemblies therefore are also genomic DNA
  • We have tried some non-parametric approaches, but the problem is that in this case the 3000x coverage is highly improbable, but if for example we had a perfect run and the max coverage we have is for example 250, this would fall within our predicted distribution but when trimming for example 5% we would remove our best contigs even though they would be considered to be accurate according to the distribution.
  • Have not yet tried that, thank you for suggesting it, i will look into it right away :D 
ADD REPLYlink written 4.5 years ago by Lesley Sitter490

How did you merge the reads? If I understood you correctly, you assembled using only the unmerged reads, is this correct? Is this genome or transcriptome? Why do you want to remove them?

Anyway, if your outliers are like contig 7, finding them should be straightforward.

ADD REPLYlink written 4.5 years ago by h.mon28k

Sorry i wasn't clear on that part. We merged reads using PEAR (Pair Ended reAd mergeR), that way we end up with lots of contigs. The read pairs that couldn't be merged by PEAR were fed to SOAP and assembled. Both assemblies combined results in final assembly.

Contig 7 is straightforward, but this is just and example, as you can see in the graph there isn't a definitive threshold where correct contigs end and improbable ones start so we want some function that can set that threshold for us, preferabbly one that takes into account the distribution so that it not just removes contigs based on a predetermined static value.

ADD REPLYlink modified 4.5 years ago • written 4.5 years ago by Lesley Sitter490

PEAR is not an assembler as in "genome assembler", I would rather call it "read merger".

ADD REPLYlink written 4.5 years ago by h.mon28k
2
gravatar for Michael Dondrup
4.5 years ago by
Bergen, Norway
Michael Dondrup47k wrote:

I think you shouldn't base the decision on coverage alone, because high coverage is an indication of presence of repetitive elements as I suspected, but they might be surrounded by 'good' sequence. Instead, I would try to mask repeats using RepeatMasker and attempt to align reads to non-masked regions only. Imagine you have a region containing a long tandem repeat, this region would yield abnormal coverage as a peak, while you couldn't see the TR in the assembly, however the surrounding regions might be perfectly fine.

Also, variant calling pipelines might have a parameter for maximum local coverage to avoid similar problems.

 

Does this make sense?

Wrt your response:

  • We have tried some non-parametric approaches, but the problem is that in this case the 3000x coverage is highly improbable, but if for example we had a perfect run and the max coverage we have is for example 250, this would fall within our predicted distribution but when trimming for example 5% we would remove our best contigs even though they would be considered to be accurate according to the distribution.

In this case fitting a neg. binomial using fitdistr and using a p value of 0.05 or 0.01 as the cutoff, - or using  negative binomial regression - should yield better results, meaning that trimming would be done conditionally only.

ADD COMMENTlink modified 4.5 years ago • written 4.5 years ago by Michael Dondrup47k

We have a step where we use dustmasker to check low-complexity regions. But haven't thought about actually masking those regions we just used it to determine the percentage of the contig that is considered low-complexity. This is definitely something i will look into, might solve the coverage problem before it becomes a problem.

I am also going to check the negative binomial approach you mentioned, thank you for the tips :D

ADD REPLYlink written 4.5 years ago by Lesley Sitter490
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1696 users visited in the last hour