Getting read counts per loci
1
1
Entering edit mode
8.4 years ago
JameB11 ▴ 10

How can I get the number of mapped reads for a particular region?

samtools view -c -F 4 my.bam gives me count in the entire bam file but I can't just add -r Chr1:0:1000 to get reads in that region only.

bam RNA-Seq samtools • 3.1k views
ADD COMMENT
3
Entering edit mode
8.4 years ago
samtools view in.bam chr1:0-1000 | wc -l

but better way is to do it with bedtools if you have many regions, for e.g like exon coordinates or peak coordinates, which is more efficient way for counting purpose.

bedtools multicov -bams aln1.bam [ aln2.bam aln3.bam . . ] -bed ivls-of-interest.bed
ADD COMMENT
0
Entering edit mode

Just to confirm, the number of lines in a BAM file give the info on the number of reads it has? Apologies, I am completely new to this area.

ADD REPLY
0
Entering edit mode

Pls check the link in updated answer.

ADD REPLY
0
Entering edit mode

Just confirming that for any given region, they both should give the same number?

ADD REPLY

Login before adding your answer.

Traffic: 2975 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6