Question: Rna-Seq : How To Check For Very High Coverage Region
2
gravatar for Nicolas Rosewick
6.4 years ago by
Belgium, Brussels
Nicolas Rosewick8.3k wrote:

Hi,

How can I check for very high coverage region in a bam file. I've several samples with good alignment rate (~90%) but when I count the reads for known genes (ensembl), I don't have good results compared to previous data. I think it's rRNA or other contamination. So how can I check for high coverage region in my data (more than 10000 reads per position) ?

Thanks a lot,

N.

coverage rna-seq • 2.0k views
ADD COMMENTlink modified 6.4 years ago by Ian5.5k • written 6.4 years ago by Nicolas Rosewick8.3k

Various relevant methods are discussed here: What is the fastest method to determine the number of positions in a BAM file with >N coverage?

ADD REPLYlink written 6.4 years ago by Malachi Griffith18k
1
gravatar for bruce.moran
6.4 years ago by
bruce.moran650
Ireland
bruce.moran650 wrote:

Hi,

I was going to suggest 'samtools depth' but that doesn't appear in the manual anymore, might still be functioning though (works on my cluster)... Essentially, depth used to give number of reads aligning to a region. So you would give it a BAM and the region and reference.fa and get your (per base) 'depth', then awk out those higher than 10000.

Bruce.

ADD COMMENTlink modified 6.4 years ago • written 6.4 years ago by bruce.moran650
1
gravatar for Dan Gaston
6.4 years ago by
Dan Gaston7.1k
Canada
Dan Gaston7.1k wrote:

You could use the GATK DepthOfCoverage walker with a provided BED file. Usually these are per exon but any sort of intervals will do. I believe you can get per position coverage information from that tool as well. Otherwise for intervals you can set what sort of percentage/depth cutoffs you want to report.

ADD COMMENTlink written 6.4 years ago by Dan Gaston7.1k
1
gravatar for Damian Kao
6.4 years ago by
Damian Kao15k
USA
Damian Kao15k wrote:

You can generate a mpileup output with samtools. Then use a script to scan through the file to get the top X bases with the highest read mapping. You can use this top X list to narrow down the regions with highest coverage.

To get counts for each base, you can simply count the number of '.' or ',' in the 4th column of the mpileup output. A '.' (period) stands for read mapping to forward strand and ',' (comma) stands for read mapping to reverse strand.

ADD COMMENTlink modified 6.4 years ago • written 6.4 years ago by Damian Kao15k
1
gravatar for Ian
6.4 years ago by
Ian5.5k
University of Manchester, UK
Ian5.5k wrote:

I would use bedtools coverage:

Summary: Returns the depth and breadth of coverage of features from A on the intervals in B.

bedtools coverage -abam reads.bam -b interesting_intervals.bed > output

other flags can help you get the format you want.

ADD COMMENTlink written 6.4 years ago by Ian5.5k
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: 1354 users visited in the last hour