Debug failed RNA-seq alignment and featureCounts: FASTQ QC, index mismatch, strandedness, multimapping

TL;DR

For "Fatal INPUT FILE error, no valid exon lines in the GTF file": 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.

Fatal INPUT FILE error, no valid exon lines in the GTF file

When to use

You are seeing this: 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. Use this skill when you run into "Debug failed RNA-seq alignment and featureCounts: FASTQ QC, index mismatch, strandedness, multimapping".

When not to use

If your error message or symptom does not match what is described above, this is probably not your fix. Search for your exact error text instead of forcing this one to fit.

Versions

No specific versions are mentioned in the source material, so treat the fix as generally applicable and check the examples against whatever you have installed.

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:

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

cat SAMPLE_Log.final.out

The fields that matter:

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

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

# 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:

# 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

# 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

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:

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:

# 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:

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:

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.

Subread/Rsubread Users Guide (Bioconductor): -s values 0 unstranded (default), 1 stranded, 2 reversely stranded; -p/--countReadPairs, -B, -C, -M, --fraction, --primary, -O, -Q, --ignoreDup flag semantics; -t exon/-g geneid defaults; -R CORE per-read report. RSeQC manual: inferexperiment.py output 'Fraction of reads explained by 1++,1--,2+-,2-+ / 1+-,1-+,2++,2--' and unstranded/forward/reverse interpretation. STAR: Log.final.out field names (Uniquely mapped reads %, Number of reads mapped to multiple loci, Number of reads unmapped: too short, Mismatch rate per base), Fatal INPUT FILE error no valid exon lines caused by FASTA/GTF chromosome-name mismatch, --outFilterMultimapNmax default behavior. BioStars and Galaxy help threads: coordinate/build mismatch produces high unique-mapping with high UnassignedNoFeatures; NCBI GFF lacks geneid attribute causing featureCounts 9th-column error.