Running TopHat for RNA-Seq, and What to Run Instead

By the end of this guide you will have a coordinate-sorted BAM file of spliced alignments from your own RNA-Seq FASTQs, a splice junction file, and a gene-level count matrix. You will also have a transcript-level abundance table in TPM, plus enough quality metrics to know whether any of it is trustworthy. You will have produced these two ways: once with TopHat2, because a large body of published pipelines and older lab protocols still specify it, and once with the aligner we would use today. You need paired-end FASTQ files, typically 20–60 million read pairs per sample for a standard poly-A library. You also need a reference genome FASTA and a matching gene annotation GTF. Finally, you need a Linux machine with at least 32 GB of RAM and a few hundred gigabytes of free disk, plus comfort at a shell prompt. All examples use human GRCh38 and GENCODE annotation.

1. Decide whether you need TopHat at all

TopHat solved a specific problem in 2009: aligning short reads to a genome when a read may be split across an intron tens of thousands of bases long, without knowing in advance where the introns are. It mapped reads with Bowtie, collected the unmapped ones, and split them into segments. It then used the segment alignments to infer splice sites ab initio and realigned reads across those inferred junctions.¹ TopHat2 rebuilt this on Bowtie2 and added handling of insertions, deletions, and gene fusions. It also gained the ability to align first against a transcriptome index built from a supplied GTF before falling back to the genome.² The TopHat/Cufflinks pair became the default teaching pipeline for genome-guided expression analysis and is the pipeline described in most of the protocol chapters from that era.³⁴

TopHat is no longer maintained. Its own authors have recommended HISAT2 as the successor for several years, and the practical consequences of using TopHat in 2026 are concrete. It is slow, taking hours per sample where HISAT2 takes tens of minutes. It depends on a Python 2 and Boost toolchain that most modern distributions no longer ship, and its indel and junction handling has been superseded. We would run TopHat2 for exactly two reasons. The first is reproducing a specific published result where the alignment step is part of what you are checking. The second is reanalyzing a dataset whose downstream files (the junctions.bed, the accepted_hits.bam read groups) another tool expects in TopHat’s exact format. For new analysis of your own data, skip to step 5.

2. Quality-control and trim the reads first

Alignment rate is the metric everyone looks at, and adapter contamination is the most common reason it is bad. Run fastp, which does adapter detection, quality trimming, and reporting in a single pass:

fastp \
  -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
  -o sample_R1.trim.fastq.gz -O sample_R2.trim.fastq.gz \
  --detect_adapter_for_pe \
  --cut_front --cut_tail --cut_mean_quality 20 \
  --length_required 36 \
  --thread 8 \
  --json sample.fastp.json --html sample.fastp.html

Read the JSON as well as the HTML. The HTML alone is not enough. Three fields matter. A duplication.rate above roughly 40–50% on a 30-million-pair library suggests you are sequencing a small number of distinct molecules, usually from a low-input or over-amplified prep. When adapter_cutting.adapter_trimmed_reads is a large fraction of total reads, the inserts are short, which will also show up later as a compressed fragment-length distribution. The read1_before_filtering.total_cycles field tells you your read length, which you need for the TopHat segment parameters below.

Do not aggressively quality-trim beyond this. Hard trimming to high Phred thresholds removes real signal and biases coverage at read ends, and modern spliced aligners soft-clip poor tails on their own.

3. Build the Bowtie2 and transcriptome indexes

TopHat2 sits on Bowtie2, so you need a Bowtie2 index whose basename matches the FASTA filename, in the same directory.

bowtie2-build --threads 16 GRCh38.primary_assembly.genome.fa GRCh38.primary_assembly.genome

This takes roughly an hour on 16 cores and produces about 4 GB of .bt2 files. Then build the transcriptome index once, rather than letting TopHat rebuild it on every run:

tophat2 \
  -G gencode.v44.annotation.gtf \
  --transcriptome-index=tx_index/gencode_v44 \
  -o /tmp/tophat_index_build \
  -p 16 \
  GRCh38.primary_assembly.genome

TopHat will build the index and then exit if you give it no reads. Two things break here reliably. First, your GTF chromosome names must match your FASTA chromosome names exactly. GENCODE uses chr1 and Ensembl uses 1, and mixing them yields an index with zero usable transcripts and no error message you would notice. Check with cut -f1 gencode.v44.annotation.gtf | grep -v '^#' | sort -u | head against grep '^>' GRCh38.primary_assembly.genome.fa | head. Second, TopHat requires Python 2.7 on PATH. The cleanest fix is a container rather than fighting your system Python:

docker run --rm -v "$PWD":/data -w /data \
  quay.io/biocontainers/tophat:2.1.1--py27h5e3d4d5_5 \
  tophat2 --version

4. Run TopHat2

For 2×100 bp paired-end human data with annotation available:

tophat2 \
  -o tophat_out/sample \
  -p 16 \
  --transcriptome-index=tx_index/gencode_v44 \
  --library-type fr-firststrand \
  --mate-inner-dist 150 --mate-std-dev 60 \
  --min-intron-length 50 --max-intron-length 500000 \
  --segment-length 25 --segment-mismatches 2 \
  --min-anchor-length 8 \
  --no-coverage-search \
  GRCh38.primary_assembly.genome \
  sample_R1.trim.fastq.gz sample_R2.trim.fastq.gz

The flags that change results, in order of how often people get them wrong:

--library-type must match your kit. Illumina TruSeq Stranded mRNA and most dUTP-based protocols are fr-firststrand. Unstranded libraries are fr-unstranded. Getting this wrong will silently invert strand assignment in every downstream quantification step. Alignment itself does not usually suffer much. If you do not know, align a few million reads and check with infer_experiment.py from RSeQC against a BED12 of your annotation.

--mate-inner-dist is the gap between mates. The fragment length is not what it asks for. If your library has a 350 bp mean fragment and 2×100 bp reads, the inner distance is about 150. TopHat2 is more forgiving of a bad estimate than TopHat1 was, but a value off by more than a couple of hundred bases will cost you properly-paired alignments.

--segment-length should be at most half your read length. The default is 25, which is right for reads of 50 bp and up. For 2×36 bp reads drop it to 18. This parameter controls the ab initio junction discovery: reads that fail to align whole are cut into segments of this length, and junctions are inferred where neighbouring segments land far apart.¹

--no-coverage-search disables the island-based junction search. Leave it off for annotated genomes with reasonable depth. Turn it on only for shallow libraries or poorly annotated organisms, and expect the run time to increase several-fold.

The output directory contains accepted_hits.bam (already coordinate-sorted), unmapped.bam, junctions.bed with a score column giving the supporting read count for each junction, deletions.bed, insertions.bed, and align_summary.txt. Read align_summary.txt before anything else. For good human poly-A data you should see a concordant pair alignment rate above 85%, with overall mapping above 90%.

5. The pipeline we would use instead: HISAT2 plus featureCounts

HISAT2 does the same job, produces a compatible BAM, and runs an order of magnitude faster on a human index that fits in about 4 GB of RAM. Extract known splice sites from the GTF and build the index once:

hisat2_extract_splice_sites.py gencode.v44.annotation.gtf > splicesites.txt
hisat2_extract_exons.py gencode.v44.annotation.gtf > exons.txt

hisat2-build -p 16 \
  --ss splicesites.txt --exon exons.txt \
  GRCh38.primary_assembly.genome.fa GRCh38_tran

The graph index build with --ss and --exon needs roughly 160–200 GB of RAM. If you do not have that, download the prebuilt grch38_tran index rather than building without annotation.

hisat2 \
  -p 16 -x GRCh38_tran \
  --rna-strandness RF \
  --dta \
  --known-splicesite-infile splicesites.txt \
  --summary-file sample.hisat2.summary \
  -1 sample_R1.trim.fastq.gz -2 sample_R2.trim.fastq.gz \
| samtools sort -@ 8 -m 2G -o sample.bam -
samtools index -@ 8 sample.bam

--rna-strandness RF is the HISAT2 spelling of fr-firststrand. --dta reports alignments tailored for transcript assembly, which you want if StringTie comes next. Then count reads per gene:

featureCounts \
  -T 16 -p --countReadPairs \
  -s 2 \
  -t exon -g gene_id \
  -a gencode.v44.annotation.gtf \
  -o counts.txt \
  sample.bam

-s 2 is reverse-stranded and must agree with --rna-strandness RF. --countReadPairs counts fragments rather than reads, which is what you want for paired data.

If your only goal is expression of known transcripts, we would skip genome alignment entirely and use salmon with selective alignment. It is faster, gives transcript-level TPM directly, and corrects sequence-specific and GC fragment bias:

salmon index -t gencode.v44.transcripts.fa -i salmon_idx --gencode -k 31

salmon quant -i salmon_idx -l ISR \
  -1 sample_R1.trim.fastq.gz -2 sample_R2.trim.fastq.gz \
  --validateMappings --seqBias --gcBias --numBootstraps 30 \
  -p 16 -o salmon_out/sample

The tradeoff is real. Salmon cannot tell you about novel junctions, structural events, or intronic signal. There is also no BAM to look at in IGV. Align to the genome when you want to see what the reads are doing. Quantify against the transcriptome when you want numbers.

6. Check the data before interpreting anything

RNA-Seq measures relative abundance with a dynamic range spanning several orders of magnitude, and its usefulness depends entirely on whether the library represents the RNA population you think it does.⁵ For a single personal sample, run RSeQC and look at four things.

read_distribution.py -i sample.bam -r gencode.v44.bed12
geneBody_coverage.py -i sample.bam -r housekeeping.hg38.bed -o sample
junction_saturation.py -i sample.bam -r gencode.v44.bed12 -o sample
inner_distance.py -i sample.bam -r gencode.v44.bed12 -o sample

For a poly-A whole-blood or PBMC library, expect 70% or more of assigned tags in coding exons plus UTRs, and a gene-body coverage curve that is roughly flat rather than sloping toward the 3′ end. A strong 3′ slope means degraded input RNA, and it will distort every long transcript in the sample. Junction saturation that is still climbing steeply at your full depth means you do not have enough reads to make claims about splicing. Whole blood adds one more problem: globin transcripts can consume 50–70% of your reads unless the prep used globin depletion, which collapses your effective depth for everything else.

7. Understand what one transcriptome can and cannot tell you

Expression is tissue-specific and time-specific. A blood transcriptome tells you about circulating immune cells and reflects the hour you drew the sample, your last meal, your last workout, and whether you were fighting a cold. That sensitivity is what makes RNA-Seq informative about state rather than just genotype, and it is why the technique became standard across fields that need to read condition-dependent biology.⁶

The consequence for a single person is that you should treat one sample as one observation in a series, not as a measurement of “your expression.” Differential expression tools assume replication, and with n=1 versus n=1 the fold change has no error estimate. What works is repeated sampling of yourself under recorded conditions, so that you accumulate a within-person distribution against which a later sample can be compared. Interpretation gets substantially stronger when transcript levels sit alongside genotype, protein, and metabolite data from the same person. Each layer constrains the others.⁷⁸

Two boundaries are firm. RNA-Seq of a blood sample is not a clinical test, and gene-level counts do not establish disease. If something in your data appears to concern your health, take the underlying question to a physician, who can order a validated assay in a certified laboratory. Nothing in an accepted_hits.bam file is a basis for changing what you take or do.

Common problems

Low alignment rate with high unmapped reads, typically under 60%. Check unmapped.bam by pulling a few thousand reads and BLASTing them. The usual culprits are adapter dimer that trimming missed, rRNA carryover from a failed depletion, or a species contaminant. A quick seqkit head -n 5000 followed by a look at the most frequent sequences resolves this in minutes.

TopHat reports almost no junctions. Nearly always a GTF and FASTA naming mismatch, or a --transcriptome-index path that TopHat silently rebuilt as empty. Confirm with wc -l tophat_out/sample/junctions.bed, which should be in the hundreds of thousands for human at 30 million pairs.

TopHat fails to build or crashes on startup. Do not patch the Python 2 dependency by hand. Use the biocontainer image shown above, or accept that the modern path is HISAT2. The features that justified keeping TopHat2 around, particularly its indel-aware and fusion-aware alignment, are now covered by better-maintained tools.²

Counts that disagree between featureCounts and salmon by more than a factor of two for specific genes. Expected for genes with high multi-mapping (paralog families, immunoglobulin loci, histones). featureCounts discards ambiguous fragments by default while salmon distributes them via its expectation-maximization model. Neither is wrong. Decide which behaviour you want before comparing samples.

Strandedness disagreement. If infer_experiment.py reports roughly 50/50 strand assignment, your library is unstranded and both -s 2 and --rna-strandness RF are wrong. Rerun with fr-unstranded and -s 0.

Memory failures during hisat2-build --ss --exon. The graph index is the memory-hungry step. Alignment is not where the memory goes. Build on a high-memory instance once, or use the prebuilt index and pass known splice sites at alignment time with --known-splicesite-infile, which recovers most of the benefit.

Footnotes

1. Cole Trapnell, Lior Pachter, Steven L. Salzberg. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics, 2009. https://doi.org/10.1093/bioinformatics/btp120 ↩ ↩²

2. Daehwan Kim, Geo Pertea, Cole Trapnell, et al. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biology, 2013. https://doi.org/10.1186/gb-2013-14-4-r36 ↩ ↩²

3. Jacob Pollier, Stephane Rombauts, Alain Goossens. Analysis of RNA-Seq Data with TopHat and Cufflinks for Genome-Wide Expression Analysis of Jasmonate-Treated Plants and Plant Cultures. Methods in Molecular Biology, 2013. https://doi.org/10.1007/978-1-62703-414-2_24 ↩

4. Sreya Ghosh, Chon-Kit Kenneth Chan. Analysis of RNA-Seq Data Using TopHat and Cufflinks. Methods in Molecular Biology, 2016. https://doi.org/10.1007/978-1-4939-3167-5_18 ↩

5. Zhong Wang, Mark Gerstein, Michael Snyder. RNA-Seq: a revolutionary tool for transcriptomics. Nature Reviews Genetics, 2009. https://doi.org/10.1038/nrg2484 ↩

6. Erica V. Todd, Michael A. Black, Neil J. Gemmell. The power and promise of RNA ‐seq in ecology and evolution. Molecular Ecology, 2016. https://doi.org/10.1111/mec.13526 ↩

7. Biswapriya B Misra, Carl Langefeld, Michael Olivier, et al. Integrated omics: tools, advances and future approaches. Journal of Molecular Endocrinology, 2019. https://doi.org/10.1530/jme-18-0055 ↩

8. Xiaofeng Dai, Li Shen. Advances and Trends in Omics Technology Development. Frontiers in Medicine, 2022. https://doi.org/10.3389/fmed.2022.911861 ↩