TopHat 1.3.0 Splicing Bug ?

[Dario1984]

Hello,

I was getting some strange results with cufflinks, where exons were being called as separate transcripts, so I delved deeper, and found that TopHat wasn't mapping across some splice junctions.

I decided to take the gap between exon 3 and exon 4 of GAPDH as an example (this intron is longer than 70 bp.). I had lots of reads mapping in exon 3 and in exon 4 but zero reads across the junction. So what I did was I made a 'seed' by pasting the last 15 bases of exon 3 to the first 15 bases of exon 4, and looking for this sequence in the fastq file, and there are lots of reads which have the seed sequence inside them (and extend out further, by visual inspection).

I made a small size fastq file containing reads spanning this junction and reads that fall within the two exons. The datafile is here [530 KB download].

I'm using 1.3.0 (the latest minor release from June 10) and hg18.
My search sequences were :
junction : GACCTCAACTACATGGTTTACATGTTCCAA
exon 3 : GCTTTTAACTCTGGTAAAGTGGATATTGC
exon 4 : GATTCCACCCATGGCAAATTCCATGGCAC

all searches allowed 1 mismatch.

Can the TopHat developers (or forum enthusiasts) run this example, to verify that it creates an empty junctions.bed file and misses all of the spliced mappings ? I tried -F 0 as a final idea, and it doesn't make it better.

[agent99]

I think that if your segment-length is longer than 1/2 your read length you won't get any mapped splice junctions.

[Dario1984]

But most of the reads are 75 bases long, some are shorter because of adapter trimming, but these are the minority, and the default segment length of TopHat is 25 bases, so this shouldn't be an issue ?

[agent99]

Ah, I just assumed you had short reads because you were pasting together such short segments of exons.

[spikesd17]

did you find an answer to this problem? I am having the same problem.(probably) I was using v1.3.1 though. It seems that segment_jcns is failing before it runs... I found I encountered the problem after getting sequencing data from a new machine/facility (singapore) that uses illumina 1.5 rather than 1.3 to encode quality scores. I verified there should be spliced reads via blat. It seems strange that running tophat in --phred64-quals mode would use intermediate files (found in the tophat tmp/ directory with the --keep-tmp option) which aren't in phred64 format:


== fastq lines

@1
AGCAATATAAAGTTCTGGTGCTA
+HWUSI-EAS1712_0042_FC:5:1:3971:1222#0/1
EE@BBHHHHHHHGHGHHHHH:;?
@2
GCTTTGGTGACCATCAGTAG
+HWUSI-EAS1712_0042_FC:5:1:1971:1230#0/1
HHHHHHHFHHGGDGG<GACG
@3
CTCCTTTTCCATTTCATGGAC
+HWUSI-EAS1712_0042_FC:5:1:8107:1223#0/1
HHEHHHHDHGGGEDDDEGGB@
@4
TCAGGATCAGGCTCTGGTTCAGCA
+HWUSI-EAS1712_0042_FC:5:1:3863:1223#0/1
IIIIIIIIIIIIIIIIIIII=;=@
==

Anyways, I'm still trying to fix it...

[Dario1984]

No, I didn't resolve it. I sent a report to the e-mail address on TopHat's website, but they seem to completely ignore bug reports about their software. It seems like there is no motivation to fix it once it is published in Nature Biotechnology

[spikesd17]

that's no good. How are you living without your spliced reads?! Should we pester Cole Trapnell?

[Dario1984]

It doesn't happen for me on all junctions. I manually made a file of junctions sequences based on RefSeq, and then mapped to this with bowtie. It rescued about 500000 reads that were being discarded by TopHat mapping. I have 18 million mapping, so I'm not fussed about getting those last few mappings.

I'm avoiding doing differential exon usage analysis for now, as our experimental protocol seems to only get a decent coverage of reads at the 3' end. There are very few reads around TSSs.

I think it is best if you try to get a response from [email protected]

[lamon]

I've had similar problems and worse because sometimes it misaligns the reads altogether if there is a second isoform with a slightly different exon length. I think is somewhat due to the how the program uses the binucleotide sequence immediately after and before exons. I have tried using GSNAP and I find it does better a job getting the junctions correct but sometimes it goes a little too far and splices rather distant exons from other genes and intergenic regions even when it could match to a neighboring exon. Decreasing the allowable local splice distance just results in an unspliced read with several mismatches (even though I set the allowable number to 2). For the same junction and set of reads, tophat won't align them at all. I would love to hear other users experiences and tips.
Thanks,
Lynn

[spikesd17]

Yeah, I've switched to using GSNAP as well and I do notice several long-distance spliced reads. I think that these might be worth it, given the increased coverage GSNAP gives, relative to tophat, of spliced reads. Moving to GSNAP now, tophat had too many problems and no user support.

[lamon]

I just figured out one way the decrease the long-distance alignments from gsnap. I built a splice site file as described in the README file and reran the alignment and so far, it seems that misaligned reads have aligned correctly. And, a read that I knew was really alternative spliced aligned also so I am feeling pretty good with this method.

[spikesd17]

Hmm, I glad I started using that splice sites file outright then. I didn't try any alignments without the .iit file.