This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How do you remove coverage outliers in sequencing count data?

I've recently looked at methylation coverage (bisulfite seq data) in IGV, and there are obvious coverage outliers which I'm interpreting as mismappings of repetitive sequences. How do you identify these outliers from methylation calls/count data (I don't have BAM/SAM files) and remove them?

bisulfite sequencing qc r quality control
• 3,716 views • link
updated by dariober dariober • written by biohack92 biohack92

2 answers

The quantile function in R is handy in these cases. To discard datapoints in the top 1% (and keep the bottom 99%) do:

xv<- rnbinom(10000, 1, 0.1)     ## Test data
qqcut<- quantile(xv, p= 0.99)   ## Cut off to use in R or awk
xvcut<- xv[xv < qqcut]          ## Keep only points in lower 99%
• 0 views • link
written by dariober dariober

Load the coverage into R, plot the distribution, and then use awk to filter the files given some reasonable threshold given the coverage distribution.

• 0 views • link
written by Devon Ryan Devon Ryan

Thanks @Devon Ryan. I followed your advice and this is the plot I created. X-axis shows the # of reads/counts and Y-axis is the frequency of each count. Is there a way to determine which threshold is 'reasonable'?

• 0 views • link
written by biohack92 biohack92

You could make a qq-plot, but it's generally OK to just eye-ball things. From the distribution, it looks like a threshold around 200 would be reasonable.

• 0 views • link
written by Devon Ryan Devon Ryan

Log in to answer this question.

More posts like this