Question: How Many Reads In A Bam File Are Aligned To An Arbitrary Region
gravatar for David Quigley
10.0 years ago by
David Quigley11k
San Francisco
David Quigley11k wrote:

I have BAM files from an Illumina GA experiment. I'd like to know how many of my aligned reads fall in a particular region of a given chromosome. It's a bonus if I can filter this by alignment quality. I've checked Picard and SAMtools and don't see an obvious canned command. It looks like I can use PySam to code this myself (e.g. by adapting the example in the

next-gen sequence sequencing • 7.3k views
ADD COMMENTlink modified 9.5 years ago by Aaronquinlan11k • written 10.0 years ago by David Quigley11k
gravatar for Aaronquinlan
10.0 years ago by
United States
Aaronquinlan11k wrote:

Assuming you've sorted and indexed your bam file with "samtools sort" and samtools index", you can pull out reads from a given region as follows (for example from chr1:1000000-2000000)

samtools view <BAM> chr1:1000000-2000000 | wc -l

Now, if you also want to limit on mapping quality, use the -q option (here I am asking for MAPQ >= 10).

samtools view -q 10 <BAM> chr1:1000000-2000000 | wc -l

If you need to look for reads overlapping many regions, use BEDTools.

intersectBed -abam <BAM> -b regions.bed >
ADD COMMENTlink modified 22 months ago by RamRS27k • written 10.0 years ago by Aaronquinlan11k
Please log in to add an answer.


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