Issues using tophat2 to align paired sequences

[EReckase]

I work on software supporting targeted-amplicon sequencing where the insert size is typically shorter than the read size. One of my tasks was to evaluate tophat2 as an aligner for some of our more difficult cases, but I'm getting a really lousy alignment rate (granted, this is NextSeq data that I'm looking at, but even that can't explain the low percentages.)

Is there a way to debug this sort of thing, a log I can look at that will explain why a given read did not align? I get alignments with bwa, but tophat2 frequently will miss aligning R2. I suspect that my settings are not quite right, but I've varied the edit distance parameters as well as the inner distance metrics and nothing seems to help.

[Brian Bushnell]

If your insert size is shorter than read length, you should adapter-trim the reads prior to mapping. Are your reads paired, and do you know the adapter sequence?

[EReckase]

The reads are paired, and the adapters have already been trimmed off.

[Brian Bushnell]

Possibly, the error rate in R2 is too high. Since this is NextSeq, the quality scores won't give very useful information about the error rates; you'll have to determine that empirically by mapping. You could try BBMap, which is more error-tolerant than Tophat2 or bwa, and allows you to plot the error-rate histogram across both read1 and read2 (with the "mhist" flag), which may be a useful analytic tool in this case. Posting the base frequency histogram across the reads may also be useful.

Also, depending on the read lengths, you could consider error-correcting the data if you want to increase the mapping rates.

[dpryan]

I'm not sure why you would bother with tophat2 for amplicon sequencing. It has fixed parameters at various steps that you can't actually change from the command line (this is likely causing some of your difficulties). You'd probably be better off with a standard DNA aligner (BBmap, bwa, bowtie2, etc.).

Also, I'm curious why you bothered using paired-end reads if the read lengths are normally longer than the insert. That would seem to just be a waste of reads.

[EReckase]

Can BBMap handle fusion/split alignments?

[Brian Bushnell]

It can handle very long deletions, but not arbitrary rearrangements like inter-chromosomal fusions. Or, to be more precise, it only produces one sam line per read, so a read with part on one chromosome and part on another will map to the chromosome from which it got the majority of bases (if the "local" flag is enabled). But if it is a fusion created by, say, skipping 100kbp within a chromosome, that will be reported entirely in a single alignment (as Tophat2 would). To find these you should set the "maxindel" flag, which defaults to 16000. Depending on read length, the sensitivity drops eventually; I don't recommend setting it much above 200000 for 100bp reads.

[EReckase]

I'm not sure why you would bother with tophat2 for amplicon sequencing. It has fixed parameters at various steps that you can't actually change from the command line (this is likely causing some of your difficulties). You'd probably be better off with a standard DNA aligner (BBmap, bwa, bowtie2, etc.).

Also, I'm curious why you bothered using paired-end reads if the read lengths are normally longer than the insert. That would seem to just be a waste of reads.

We're performing targeted sequencing of an RNA sample that will cross exon-exon boundaries. Furthermore, we're looking for fusions in this data, which is why we were interested in tophat.

As far as using paired-end reads, the read lengths are not always longer than the amplicon sequence. Our chemistry is unique in that the amplicons are not always the same length, but they are anchored at one end (Archer AMP).