使用TopHat/Bowtie2做RNA-seq比对异常：55M保留reads仅1M用于比对

Let’s work through what’s happening with your RNA-seq alignment—your 13GB FASTQ (55M total reads) resulting in a tiny 50MB accepted_hits.bam file, paired with that confusing note about only 1M reads being used for alignment, is definitely worth digging into. Here’s how to diagnose the issue step by step:

1. Parse the TopHat Log Files First #

Your first stop should be the main tophat.log file and the align_summary.txt you mentioned. Look for these critical details:

• Initial read filtering stats: TopHat automatically filters low-quality reads, reads with too many Ns, or unpaired reads (if you’re working with paired-end data) before sending anything to Bowtie2. Check if a large chunk of your 55M reads was discarded here—this would explain the small BAM file.
• Clarify the "1M reads used" metric: That number might not mean what you think. It could refer to unique mapped reads used for novel junction discovery, not the total number of aligned reads. Compare it to the overall alignment rate (97%)—if 97% of 55M reads is ~53.3M, that’s way more than 1M, so this stat is likely a subset of aligned reads, not the total.

2. Audit the BAM File Directly #

Don’t rely solely on TopHat’s summaries—verify the actual content of your accepted_hits.bam with command-line tools:

• Count total reads in the BAM:

  samtools view -c accepted_hits.bam
  

• Get a detailed breakdown of alignment statuses:

  samtools flagstat accepted_hits.bam
  

This will tell you exactly how many reads are aligned, unaligned, multi-mapped, etc. If the total read count here is way lower than 55M, it confirms reads were filtered out before or during alignment.

3. Verify Your Reference Index and Input Data #

• Check Bowtie2 index completeness: Ensure you’re using the full set of iGenomes hg19 Bowtie2 index files (you should have genome.1.bt2, genome.2.bt2, genome.3.bt2, genome.4.bt2, genome.rev.1.bt2, genome.rev.2.bt2). Missing files can cause partial or failed alignments, leading to a tiny BAM.
• Confirm read length: A 13GB FASTQ with 55M reads works out to ~236bp per read (uncompressed), which is reasonable. But double-check with a quick peek:

  head -20 your_sample.fastq
  

Short reads would compress more, but this size doesn’t add up to ultra-short reads, so this is probably not the issue—but it’s worth ruling out.

4. Review Your TopHat Parameter Settings #

Did you use any restrictive parameters that might be filtering reads aggressively? For example:

• --max-multihits set to 1 (which discards all multi-mapped reads, a large portion of RNA-seq data)
• --read-mismatches or --read-gap-length set to extremely low values, rejecting reads with even minor mismatches
• --no-novel-juncs which limits alignment to known transcripts only—this can filter out reads from unannotated regions, but rarely to this extreme.

Quick Reality Check on File Sizes #

BAM is a compressed format, so it’s normal for it to be smaller than uncompressed FASTQ—but 50MB for 55M reads is way too small. A typical aligned BAM for 50M 200bp reads is usually several gigabytes, not megabytes. This confirms your BAM has far fewer reads than expected, so the earlier steps will help you find where they’re being lost.

Start with the log file and samtools flagstat—those two steps will almost always pinpoint whether reads are being filtered early, misaligned due to index issues, or if that "1M reads" stat was misinterpreted.

内容的提问来源于stack exchange，提问作者Alex Stuckel