Version Compatibility

Reference examples tested with: Subread 2.0+, STAR 2.7.11+, HISAT2 2.2.1+, DESeq2 1.42+, edgeR 4.0+, pandas 2.2+

Before using code patterns, verify installed versions match. If versions differ:

• CLI: <tool> --version then <tool> --help to confirm flags
• Python: pip show <package> then help(module.function) to check signatures
• R: packageVersion('<pkg>') then ?function_name to verify parameters

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

featureCounts Counting

"Count reads per gene from my BAM files" -> Assign each aligned read to at most one gene by overlap with a GTF, discarding ambiguous reads, to produce an integer gene-by-sample matrix for differential expression.

• CLI: featureCounts -a genes.gtf -o counts.txt sample1.bam sample2.bam

featureCounts is bookkeeping, not inference: it tallies reads to genes and discards anything ambiguous. That is correct for gene-level DE, where almost every read's gene of origin is unambiguous even when its isoform is not. The two settings that silently corrupt the matrix if wrong are strandedness (-s) and, for paired-end data, fragment counting (--countReadPairs).

Basic Usage

# Multiple samples in one run -> a single aligned matrix (recommended)
featureCounts -a annotation.gtf -o counts.txt sample1.bam sample2.bam sample3.bam

# Defaults: -t exon -g gene_id (count reads over exons, aggregate by gene)

Decision 1: strandedness (-s) is the load-bearing setting

The dominant chemistry, dUTP / Illumina TruSeq Stranded / NEBNext Directional, is reverse, -s 2. Setting the wrong strand does not error; it silently destroys the matrix. For a truly stranded library, the correct -s assigns ~80-90% of reads while the opposite setting collapses to ~5-20% (counting only antisense background). Do not trust the kit name; determine it empirically:

# Method A: RSeQC reports the strand pattern fractions
infer_experiment.py -r genes.bed -i sample.bam

# Method B: run all three and pick the one that maximizes Assigned in the .summary
for s in 0 1 2; do featureCounts -s $s -a annotation.gtf -o counts_s$s.txt sample.bam; done

If -s 1 and -s 2 give wildly different Assigned fractions, the data are stranded (use the higher); if both are roughly equal and about half of -s 0, the data are unstranded. STAR --quantMode GeneCounts provides a free cross-check (see below).

Decision 2: paired-end fragment counting

# Subread >= 2.0.2: -p only declares paired input; --countReadPairs is REQUIRED to count fragments
featureCounts -p --countReadPairs -a annotation.gtf -o counts.txt *.bam

# Stricter: require both ends mapped, exclude chimeric/discordant pairs
featureCounts -p --countReadPairs -B -C -a annotation.gtf -o counts.txt *.bam

Omitting --countReadPairs on paired-end data counts each mate separately, roughly doubling counts and breaking the count model. -B requires both ends aligned; -C excludes pairs mapping across chromosomes or in the wrong orientation.

Decision 3: multi-mapping and multi-overlap reads

# Default (recommended for gene-level DE): discard both -> uniquely, unambiguously assigned reads only
featureCounts -a annotation.gtf -o counts.txt *.bam

# Count multimappers fractionally (1/N) or fully (1 each) -- NOT recommended for DE
featureCounts -M --fraction -a annotation.gtf -o counts.txt *.bam
featureCounts -M -a annotation.gtf -o counts.txt *.bam

# Count reads overlapping >1 gene in all of them
featureCounts -O -a annotation.gtf -o counts.txt *.bam

Discarding multimappers is the right default for gene-level DE. -M --fraction looks principled but biases exactly the genes where resolution matters: a read truly from gene A that also maps to paralog A' is split 0.5/0.5, diluting both. This is the regime where alignment-free EM quantifiers (rna-quantification/alignment-free-quant) outperform featureCounts, because they reassign by full likelihood rather than a flat split.

Quality and feature options

featureCounts -Q 10 -a annotation.gtf -o counts.txt *.bam      # min MAPQ (aligner-specific scale)
featureCounts --primary -a annotation.gtf -o counts.txt *.bam  # primary alignments only
featureCounts -t CDS -g gene_id -a annotation.gtf -o counts.txt *.bam  # count CDS instead of exon

-Q thresholds mapping quality, but MAPQ conventions are aligner-specific (STAR assigns 255 to unique reads, low values to multimappers), so confirm the scheme before choosing a cutoff. Do NOT add --ignoreDup for standard RNA-seq: high duplication is expected from highly expressed genes, and position-based deduplication discards real signal. Deduplicate only with UMIs. For exon-level usage testing (DEXSeq), use a flattened annotation rather than gene-level counting (alternative-splicing/isoform-switching).

Output

counts.txt:           Geneid Chr Start End Strand Length sample1.bam sample2.bam ...
counts.txt.summary:   Status              sample1.bam  sample2.bam
                      Assigned            1523456      1678234
                      Unassigned_NoFeatures 234567     245678

Reading the .summary is the primary QC step. A good poly-A library assigns ~70-90% of mapped reads (rRNA-depletion libraries run lower).

STAR cross-check

If aligned with STAR --quantMode GeneCounts, ReadsPerGene.out.tab gives a free independent count: column 2 = unstranded (≈ -s 0), column 3 = forward (≈ -s 1), column 4 = reverse (≈ -s 2). The larger of columns 3 vs 4 reveals the strand directly, and the per-gene counts should track featureCounts at the matching -s.

Extract the matrix

cut -f1,7- counts.txt | tail -n +2 > count_matrix.txt   # drop the 6 annotation columns

import pandas as pd
counts = pd.read_csv('counts.txt', sep='\t', comment='#')
mat = counts.set_index('Geneid').iloc[:, 5:]
mat.columns = [c.replace('.bam', '') for c in mat.columns]
mat.to_csv('count_matrix.csv')

Common Errors

Related Skills

• alignment-files/sam-bam-basics - Input BAM handling and filtering
• read-alignment/star-alignment - Producing BAMs and ReadsPerGene.out.tab strand cross-check
• genome-intervals/gtf-gff-handling - GTF/GFF annotation files
• rna-quantification/alignment-free-quant - EM-based alternative; better for multimappers/paralogs
• rna-quantification/count-matrix-qc - QC the resulting matrix before DE
• differential-expression/deseq2-basics - Gene-level DE from these counts

References

• Liao Y, Smyth GK, Shi W. 2014. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30(7):923-930. doi:10.1093/bioinformatics/btt656
• Liao Y, Smyth GK, Shi W. 2013. The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote. Nucleic Acids Res 41(10):e108. doi:10.1093/nar/gkt214