Calculating number of reads in a specific region from sam file

Question

[ag1194]

Hi,

Hi I have a mapped sequencing data in sam format (RNA). I want to calculate number of reads are on one region (let's say exon 10), then number of gene are on another region (lets say exon 11), then compare these two read number for each gene.

So basically I have a table who contains two region for each gene, and I want to calculate reads on each region to compare.

Is there any tool I can do such calculation?

Thanks a lot!

-A

[bioinforesearchquestions]

Hi,

You can use HTseq count tool for getting the number of reads assigned to genes or exons.

Refer this link (http://www-huber.embl.de/HTSeq/doc/overview.html)

[ag1194]

I want to check the reads falls into specific region though. Not a defined exon.

[microfuge]

bedtools coverage with a -d option should report that

[ag1194]

great, let me check. Thanks!

[microfuge]

Goutham's answer seems most appropriate. Bedtools will give the depth and not the number of reads. I commented too soon.

[vakul.mohanty]

If your comfortable with python pysam will let you do this pretty eaily

Answer 1

[GouthamAtla]

Use featureCounts.

Make an SAF file with regions of interest.

GeneID  Chr Start   End Strand
497097  chr1    3204563 3207049 -
497097  chr1    3411783 3411982 -
497097  chr1    3660633 3661579 -

.

featureCounts -F SAF -T <int> -a regions.saf -o counts.txt <in.bam>

You could also count paired-end reads as one fragment using -p