This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Calculating number of reads in a specific region from sam file

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

rna-seq sequencing
• 5,327 views • link
updated by GouthamAtla GouthamAtla • written by ag1194 • ag1194

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

• 1 views • link
written by ag1194 • ag1194

bedtools coverage with a -d option should report that

• 1 views • link
written by microfuge microfuge

great, let me check. Thanks!

• 1 views • link
written by ag1194 • ag1194

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

• 1 views • link
written by microfuge microfuge

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

• 1 views • link
written by vakul.mohanty vakul.mohanty

1 answer

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

• 1 views • link
written by GouthamAtla GouthamAtla

Log in to answer this question.

More posts like this