Tophat: Proper-Paired Reads

Question

[ai]

Hi Everyone,

I tried Tophat v 2.0.6 and GSNAP v 2013-02-05 to align Illumina paired-end reads. However, the numbers of proper-paired reads were very different. I also used cufflinks for assembly, but got similar number of genes/isoforms. Is there anything wrong with my commands? Why did Tophat result in such low rates of proper-paired reads? Thanks!

To get proper-paired reads: samtools view -f 0x2 accepted_hits.bam | cut -f1 | sort | uniq | wc -l

Tophat comands:

tophat -G genes.gtf -o tophat --no-novel-juncs genome read_1 read_2

Results: 1,431,500 (6.34%) proper-paired reads; 3746 genes (cufflinks).

GSNAP command

gsnap -A sam -N 0 -D ~/Software/gmap-2013-02-05/ -d mm10 -s mm10.splicesites.iit read_1 read_2

Results: 20,543,759 (90.9%) proper-paired reads;3959 genes (cufflinks).

Answer 1

[swbarnes2]

Isn't cut f1 | uniq going to compress two properly paired reads into a single read name?

What does samtools flagstat accepted_hits.bam tell you is the # of properly paired reads?

Answer 2

[Ashutosh Pandey]

Well this is not an answer but you should try '-c' to count. For e.g. samtools view -c -f 0x2 accepted_hits.bam. It should give you the counts of the properly paired mapped reads.