Tophat Pair-End Rnaseq Statistics

Question

[Xiaomin Yu]

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.

Answer 1

[Damian Kao]

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'
else
   strand = 'forward'

Answer 2

[Wen.Huang]

perhaps use a pipe to filter your bam file first?

samtools view -b -f 0x0010 input.bam | samtools flagstat

[pengchy]

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