Question: Tophat Pair-End Rnaseq Statistics
gravatar for Xiaomin Yu
7.8 years ago by
Xiaomin Yu20
Xiaomin Yu20 wrote:

Hi, I am analyzing pair-end RNAseq with tophat, I want to know how many fragments (paired) mapped to sense/antisense of the reference sequencing, and how many single reads (not paired) mapped to sense/antisense of the reference, what's tools I could use for this purpose? Thanks.

tophat • 2.8k views
ADD COMMENTlink written 7.8 years ago by Xiaomin Yu20
gravatar for Damian Kao
7.8 years ago by
Damian Kao15k
Damian Kao15k wrote:

You can get forward/reverse strand information by looking at the bitwise flag of the BAM/SAM file. [?]Here[?] is a good blog post on the bitwise flags. The post list some relavant flags:

0×0001 1 the read is paired in sequencing, no matter whether it is mapped in a pair
0×0002 2 the read is mapped in a proper pair (depends on the protocol, normally inferred during alignment)
0×0004 4 the query sequence itself is unmapped
0×0008 8 the mate is unmapped
0×0010 16 strand of the query (0 for forward; 1 for reverse strand)
0×0020 32 strand of the mate
0×0040 64 the read is the first read in a pair
0×0080 128 the read is the second read in a pair
0×0100 256 the alignment is not primary (a read having split hits may have multiple primary alignment records)

You can get the bitwise flag in python or perl with the bitwise operator '&'. For example in python:

if flag & 16:
   strand = 'reverse'
   strand = 'forward'
ADD COMMENTlink written 7.8 years ago by Damian Kao15k
gravatar for Wen.Huang
7.7 years ago by
Wen.Huang1.2k wrote:

perhaps use a pipe to filter your bam file first?

samtools view -b -f 0x0010 input.bam | samtools flagstat
ADD COMMENTlink written 7.7 years ago by Wen.Huang1.2k

Because the tophat will output multiple records for one read, the "samtools view -f " will not work.

ADD REPLYlink written 7.2 years ago by pengchy410
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 932 users visited in the last hour