Calculating number of reads in a specific region from sam file
1
0
Entering edit mode
6.2 years ago
ag1194 • 0

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 • 3.2k views
ADD COMMENT
0
Entering edit mode

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)

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

bedtools coverage with a -d option should report that

ADD REPLY
0
Entering edit mode

great, let me check. Thanks!

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
3
Entering edit mode
6.2 years ago

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

ADD COMMENT

Login before adding your answer.

Traffic: 1538 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