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
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 > reads.in.regions.bam