This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Basic Question : #Reads Used By Tophat

Hey Guys

I have a basic question in trying to make sense of the number of reads mapped by tophat. My understanding is the following

#input_reads (to tophat) = #unmapped_reads(by tophat) + #primary_mapped reads

In my case the above equation is not working out as #primarymapped reads + #unmappedreads is coming out to be more than #input_reads. My hunch is that either the flag is not correctly set in the bam file by tophat or my equation/understanding of how the tophat uses & flags the reads is not right.

Also I am calculating the #primary reads by using the sam flag -F 256

Thanks! -Abhi

rna-seq cufflinks transcriptome
• 3,349 views • link
updated by Arun Arun • written by Abhi Abhi

Did you check if your equation works if you set --max-multihits to 1?

• 0 views • link
written by Federico Giorgi Federico Giorgi

No I dint try that until now

• 0 views • link
written by Abhi Abhi

1 answer

What I'd do is to get the total number of reads in the sam file corresponding to every fastq header (or read id). If you have a paired end read (PE), one would expect this count to be not more than 2( and <= 1 for SE). You can isolate those reads that are represented >2 and then check if it is due to multiple hits as Federico points out or something else.

• 0 views • link
written by Arun Arun

I see you point..for now I was trying to do it with bitwise flag from the bam file. Counting based on read head is a bit memory intensive but worth as a check. --Thanks

• 0 views • link
written by Abhi Abhi

Log in to answer this question.

More posts like this