Why Tophat calls as unique reads which clearly map to many locations?

[Daniel]

Hello,
I am starting to use Tophat (latest build) to map RNA seq reads (human) and I am trying to understand some of the results which I am seeing. I see many reads which Tophat calls as "unique" (based on both max MAPQ score of 50, and also NH flag =1) - yet when I simply BLAT the sequences - I see equal or better alignment to many (>20) locations in human genome (hg19) to which I am aligning. Example of two ends of paired reads are:

HWI-ST1220:175:C9Q45TEXX:1:1303:19652:99794 163 chr20 25165769 50 51M = 25165847 129 TTTTCTTTAAGAATGTTAAATATTGGCCCCCACTCTCTTCTGGCTTGTAGG CCCFFFFFHHHHHJJHHJJJJIJJIJJJJJJJJJJJJJJJJJJJJJGHGGE AS:i:-11 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:23C26A0 YT:Z:UU XS:A:- NH:i:1

HWI-ST1220:175:C9Q45TEXX:1:1303:19652:99794 83 chr20 25165847 50 51M = 25165769 -129 CTGATGGGCTTCCCGTTGTGGGTAACCCGACCTTTCTCTCTGGCTGCCCTT HHHIGJJJIJJJJHHJJJJJJJIIHJJIIHJJJJIJJJHHHHHFFFFFCCB AS:i:-10 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:14T31G4 YT:Z:UU XS:A:- NH:i:1

It does not appear that mapping of both ends together would produce any clearer unique result, and I see many examples of regions where I see many apparently unique aligned reads - yet mappability of these regions is very low, and when I go back and look at individual sequences - I see that they in fact align (using BLAT) to many locations. I understand that Tophat is going to place multiply-mapped reads somewhere - but is there some better method to determine confidence level of correct placement for each read? (I am using default pms for alignment.)

Thanks for everyone's help...

[dpryan]

Are the other locations in the blat hits within the transcriptome (and if they are, would the resulting insert sizes be reasonable)? If not, there's your answer.

[Daniel]

Thanks for the quick reply. Yes - I see what you mean. The (uniquely) mapped location does lie on an annotated UCSC element, while the others do not appear to map to any known transcript. However, this behavior does not provide me with any great confidence that this particular transcript is necessarily the "true" origin of my RNA. We know that there are a lot of as yet unannotated transcripts out there, and if reads are by default assigned to only known locations (and designated as unique)- it both (possibly incorrectly) substantiates the existence of these previously annotated transcripts, and misses out on the very likely alternative origin of these RNAs.

I understand that Tophat is using known transcriptome to try and more accurately (and quickly) place reads - but is there nevertheless any way to see whether these transcriptome mapped reads are in fact unique (and thus correct)? Or does Tophat not even bother to look at entire genome once read is found in transcriptome?. I guess that I could map reads to genome without including annotated transcriptome - is this the best/only way of confirming the validity of mapping?

I am looking by the way at a nonstandard cell type - so I expect many of reads to not necessarily originate from known, annotated transcripts.

Thanks

[dpryan]

Tophat won't bother with a read/pair once it maps to the transcriptome. If you're looking at a weird cell type where you expect expression of novel species then I would strongly encourage you to not supply a GTF file during alignment. Also, use STAR, it's much faster (tophat is incredibly slow).

[Daniel]

Thanks for pointing me to STAR. It looks like not only is it superfast, but it would solve my issue of (as with Tophat) first aligning to transcriptome, and then only aligning the non-transcriptome aligned reads to genome. I look forward to giving this a try!