# Debug a Failed RNA-seq Alignment and featureCounts Pipeline
End-to-end failure checklist for an RNA-seq pipeline that produces wrong, empty, or suspiciously low gene counts: FASTQ QC (FastQC flags) → aligner output health (STAR / HISAT2) → index-vs-annotation mismatch → strandedness → multimapping and MAPQ filtering. Work the checklist in order; each stage gates the next.
## Stage 1 — FASTQ QC: FastQC failure flags
Run FastQC on every input FASTQ and read `summary.txt`, which lists one PASS / WARN / FAIL per module:
```bash
fastqc SAMPLE_R1.fastq.gz SAMPLE_R2.fastq.gz -o FASTQC_OUT
grep FAIL FASTQC_OUT/SAMPLE_R1_fastqc/summary.txt
```
Act on these modules first (they are the ones that break alignment downstream):
- **FAIL "Per base sequence quality"** — low-quality tails. Trim or lower the aligner's expectations; untrimmed tails drive STAR "Number of reads unmapped: too short" up.
- **FAIL "Per sequence quality scores"** — a large fraction of reads are globally low quality. Suspect sequencing failure, not biology.
- **FAIL "Overrepresented sequences" + "Adapter Content"** — adapter read-through, especially in short-insert libraries. Trim adapters (e.g. with fastp or cutadapt) and re-run FastQC until "Adapter Content" passes.
- **FAIL "Per sequence GC content"** — a GC distribution with extra humps or sharp spikes (not the normal single bell) means contamination or rRNA over-representation.
- **WARN/FAIL "Sequence Duplication Levels"** — high duplication can be PCR duplicates from over-amplification. Consider duplicate-aware counting later (featureCounts `--ignoreDup` only if BAM duplicates were marked).
Red line: if the raw FASTQ is failing quality and adapter modules, do not debug the aligner — fix the reads first.
## Stage 2 — Aligner health: STAR and HISAT2 summary stats
### STAR: read Log.final.out
```bash
cat SAMPLE_Log.final.out
```
The fields that matter:
| Field in Log.final.out | Healthy range | Failure reading |
|---|---|---|
| `Number of input reads` | matches your FASTQ read count | if far lower, reads were lost upstream |
| `Uniquely mapped reads %` | ~70–90% for good bulk RNA-seq | <50% → suspect contamination, adapter read-through, or wrong reference |
| `Number of reads mapped to multiple loci` / `% of reads mapped to multiple loci` | a few percent typical | high → reads too short after trimming, rRNA contamination, or paralog-rich annotation |
| `Number of reads unmapped: too short` | low | high → reads over-trimmed or adapter-contaminated (ties back to Stage 1) |
| `Mismatch rate per base, %` | ~0.2–0.5% | >1–2% → wrong reference genome, high error rate, or bisulfite/untrimmed contamination |
STAR caps reported multimappers with `--outFilterMultimapNmax` (default 10); reads exceeding the cap are dropped, so a 0 in `Number of reads mapped to too many loci` does not mean zero multimapping pressure.
### HISAT2: read the alignment summary
```bash
hisat2 -p 8 -x INDEX_BASE -1 SAMPLE_R1.fastq.gz -2 SAMPLE_R2.fastq.gz --summary-file hisat2_summary.txt 2> SAMPLE_hisat2.log
cat hisat2_summary.txt
```
Look for the `overall alignment rate`. Red lines: overall alignment rate below ~70%, or a large share "aligned concordantly 0 times" → same suspects as STAR: wrong/contaminated reads or wrong reference.
### Stage 2 diagnostic one-liners
```bash
# BAM sanity: read counts, properly-paired rate, duplicates
samtools flagstat SAMPLE.bam
# unmapped share
samtools view -c -f 4 SAMPLE.bam
# low-MAPQ share (STAR writes 255 for uniquely mapped reads)
samtools view -c -q 10 SAMPLE.bam
```
## Stage 3 — Index vs annotation mismatch
This is the classic "alignment looks great, but every count is zero" failure. The aligner's reference and the annotation given to featureCounts must be the same genome build with identical chromosome names. Check chromosome-name agreement three ways:
```bash
# 1. Chromosome names in the GTF (column 1, minus comment lines)
grep -v '^#' ANNOT.gtf | cut -f1 | sort -u | head
# 2. Chromosome names in the reference FASTA
grep '^>' REF.fa | head
# 3. Chromosome names baked into the STAR index
head STAR_GENOME_DIR/chrName.txt
# 4. Chromosome names in the BAM you will count
samtools view -H SAMPLE.bam | grep '^@SQ' | cut -f2 | head
```
Failure modes:
- **STAR dies at genomeGenerate with `Fatal INPUT FILE error, no valid exon lines in the GTF file`** — almost always a chromosome-naming difference between the FASTA and the GTF (e.g. `chr1` vs `1`). Normalise one side to match the other; STAR cannot match them otherwise.
- **GTF has no usable `gene_id` attribute in column 9** — featureCounts fails with `ERROR: failed to find the gene identifier attribute in the 9th column of the provided GTF file`. This happens with some GFF-style annotations (e.g. NCBI GFF uses `ID=gene-LOC...;Dbxref=...` instead of `gene_id "..."`). Fix by using a proper GTF (Ensembl/GENCODE) or by setting `-g` to the attribute that actually exists, e.g. `-g Dbxref` — but downstream parsing gets messy, so prefer the real GTF.
- **Coordinates disagree** — alignment unique-mapping % is high yet `Unassigned_NoFeatures` dominates the featureCounts summary. The GTF's exon coordinates do not match the reference the reads were aligned to (different build or chromosome assembly). Re-align against the correct reference or fetch the matching GTF.
Confirm with featureCounts' own summary (see Stage 5).
## Stage 4 — Strandedness: the silent killer
featureCounts has three strand settings; the Subread documentation defines them exactly:
- `-s 0` — unstranded (the default)
- `-s 1` — stranded (forward)
- `-s 2` — reversely stranded (dUTP / Illumina TruSeq stranded, the most common stranded protocol)
The default `-s 0` silently absorbs small strand errors; the silent kill comes when you *think* the data is stranded and set `-s 1` or `-s 2` wrong, or the library is stranded and you leave `-s 0` and then wonder why antisense signal leaks. Wrong `-s` shows up as near-zero `Assigned` or a flood of `Unassigned_NoFeatures` in the summary — without any error message. **Never guess strandedness from the sample sheet. Measure it.**
### Detect it with RSeQC
```bash
# convert GTF to BED12 once (RSeQC ships gtf2bed)
gtf2bed ANNOT.gtf > ANNOT.bed12
infer_experiment.py -r ANNOT.bed12 -i SAMPLE.bam
```
Output interpretation (paired-end example):
```
Fraction of reads explained by "1++,1--,2+-,2-+": 0.9441
Fraction of reads explained by "1+-,1-+,2++,2--": 0.0487
```
| Fractions | Library | featureCounts setting |
|---|---|---|
| both ~0.5 / ~0.5 | unstranded | `-s 0` |
| "1++,1--,2+-,2-+" clearly dominant (~0.9+) | forward stranded | `-s 1` |
| "1+-,1-+,2++,2--" clearly dominant (~0.9+) | reverse stranded (dUTP / TruSeq stranded) | `-s 2` |
Single-end data reports `"++,--"` and `"+-,-+"` instead. For reads, strand is read-strand convention; for paired-end fragments featureCounts uses the first read's strand.
### Verify the setting against featureCounts itself
Run counts with the candidate `-s` value and compare the summary's `Assigned` row:
```bash
featureCounts -T 8 -p --countReadPairs -s 2 -a ANNOT.gtf -o counts_s2.txt SAMPLE.bam
featureCounts -T 8 -p --countReadPairs -s 1 -a ANNOT.gtf -o counts_s1.txt SAMPLE.bam
grep 'Assigned' counts_s2.txt.summary counts_s1.txt.summary
```
The correct `-s` assigns dramatically more reads. If both settings assign very little, the problem is upstream (Stages 1–3), not strandedness.
## Stage 5 — Multimapping and low-MAPQ filtering
featureCounts default behaviour: multimapping reads are NOT counted (`-M` off), reads overlapping multiple features are not counted (`-O` off), and no MAPQ filter is applied (`-Q` unset). Decide deliberately per project:
```bash
# Require unique mappers only (the DE-analysis default)
featureCounts -T 8 -p --countReadPairs -s 2 -a ANNOT.gtf -o counts.txt SAMPLE.bam
# Keep multimappers with fractional counts (uses the NH tag; needs -M or -O)
featureCounts -T 8 -p --countReadPairs -M --fraction -s 2 -a ANNOT.gtf -o counts_frac.txt SAMPLE.bam
# Count only high-confidence alignments
featureCounts -T 8 -p --countReadPairs -Q 10 -s 2 -a ANNOT.gtf -o counts_q10.txt SAMPLE.bam
# Ignore marked duplicates and chimeric fragments
featureCounts -T 8 -p --countReadPairs -B -C --ignoreDup -s 2 -a ANNOT.gtf -o counts_strict.txt SAMPLE.bam
# Primary alignments only
featureCounts -T 8 -p --countReadPairs --primary -s 2 -a ANNOT.gtf -o counts_primary.txt SAMPLE.bam
```
Flag reference:
- `-p` — input is paired-end; `--countReadPairs` counts fragments, not individual reads.
- `-B` — count only fragments with both ends aligned; `-C` — exclude chimeric fragments.
- `-M` — count multimapping reads (STAR reports all their alignments; featureCounts identifies them via the NH tag). `--fraction` assigns each a count of 1/x instead of 1.
- `-O` — allow reads overlapping more than one feature to be counted (turns `Unassigned_Ambiguity` into assigned counts; inflates genes in dense regions — usually leave off for DE).
- `-Q INT` — minimum mapping quality.
- `--primary` — count primary alignments only (STAR marks secondary alignments in the BAM).
- `--ignoreDup` — ignore reads marked as duplicates (requires a marked BAM, e.g. from Picard or samtools markdup).
To compare against STAR's own gene counts (`ReadsPerGene.out.tab`, produced with `--quantMode GeneCounts`), note STAR counts unstranded by default and its `N_multimapping` / `N_noFeature` / `N_ambiguous` rows are the same concepts as featureCounts' summary rows.
### Read the featureCounts summary like a diagnosis
Every run writes `counts.txt.summary`. Its rows are the diagnosis:
| Row | Meaning | Typical cause when dominant |
|---|---|---|
| `Assigned` | reads counted | the healthy number |
| `Unassigned_Unmapped` | reads unmapped in BAM | aligner failure (Stage 2) |
| `Unassigned_NoFeatures` | mapped but no exon under the read | wrong `-s`, wrong GTF build, chr-name mismatch (Stages 3–4) |
| `Unassigned_Ambiguity` | read overlaps multiple features | dense annotation; `-O` counts these at the cost of specificity |
| `Unassigned_MultiMapping` | read maps to multiple loci | high for short reads or rRNA; see `-M`/`--fraction` |
| `Unassigned_MappingQuality` | below `-Q` | too aggressive a MAPQ cutoff, or poor alignments |
A STAR Log.final.out can report ~85% uniquely mapped while featureCounts assigns only ~50%: the gap is almost always multimapping detail (STAR counts a read once; featureCounts sees all NH-tagged alignments and calls it multimapping) or strand/build mismatch. Work the summary rows in the order above before touching any parameters.
## Stage 6 — Per-read evidence (last resort)
To confirm exactly why a read was not assigned, ask featureCounts to report per-read decisions:
```bash
featureCounts -T 8 -p --countReadPairs -s 2 -R CORE -a ANNOT.gtf -o counts.txt SAMPLE.bam
# produces SAMPLE.bam.featureCounts with: read name, status (assigned or reason), target count, target list
grep -c 'Unassigned_NoFeatures' SAMPLE.bam.featureCounts
```
Sample the reads, inspect their CIGAR/strand/position with `samtools view`, and compare against the GTF exons. If reads land exactly where the GTF says exons are, your annotation and reference disagree — go back to Stage 3.
## Quick decision tree
1. `Unassigned_Unmapped` high → Stage 2 (aligner) and Stage 1 (FASTQ).
2. `Unassigned_NoFeatures` high but STAR unique-mapping % fine → Stage 4 (run `infer_experiment.py`), then Stage 3 (GTF build / chr names / `gene_id`).
3. `Unassigned_MultiMapping` high → Stage 5 (`-M --fraction`, or accept unique-only).
4. Counts exist but look wrong (e.g. antisense signal, near-ties between `-s 1` and `-s 2`) → strandedness is still suspect; re-measure with `infer_experiment.py` on the actual BAM.