omics.co.in logo omics.co.in

Module 8: Gene Expression Quantification with featureCounts (Subread Suite)

1 min read • Verified October 2026 • Peer-Reviewed Learning Resource
Tested Stacks: GATK 4.6.1 Samtools 1.21 Nextflow 24.10 Seurat v5.1 Python 3.12
Prerequisites: Basic Linux CLI & Sequencing Formats
Module 08 • Quantification Track 18 • Practical Transcriptomics Hands-On Code Included

Module 8: Gene Expression Quantification with featureCounts (Subread Suite)

Building the gene-level raw count matrix: GTF/GFF3 exon feature parsing, resolving multi-overlapping reads, stranded count flags (-s 1 / -s 2), paired-end fragment counting (-p -B), and formatting the raw counts matrix.

8. Gene Expression Quantification with featureCounts

Differential expression analysis tools (DESeq2, edgeR) require raw, unnormalized integer read counts. Do NOT feed normalized TPM, RPKM, or FPKM values into DESeq2 because its statistical model relies on the sampling variability inherent in discrete Poisson/Negative Binomial counts.

Running featureCounts Across All BAM Files

# Execute featureCounts on all BAM alignments
featureCounts \
  -T 16 \
  -p --countReadPairs \
  -B -C \
  -s 2 \
  -t exon \
  -g gene_id \
  -a ref/Homo_sapiens.GRCh38.110.gtf \
  -o 05_counts/gene_raw_counts.txt \
  04_alignment/*.sorted.bam

# Format matrix for direct import into R (GeneID + Sample count columns)
cut -f 1,7- 05_counts/gene_raw_counts.txt | sed '1d' > 05_counts/counts_matrix_clean.tsv
Parameter Cheat-Sheet:
  • -p --countReadPairs: Counts fragments instead of individual reads (prevents double-counting paired reads).
  • -B -C: Only counts read pairs that have both ends successfully mapped to the same chromosome.
  • -s 2: Reverse stranded counting (standard for Illumina TruSeq Stranded / dUTP protocols).
⚠️

Common Pitfalls & Troubleshooting Checklist

  • Java Heap Exhaustion (GATK / Picard): If processes terminate with java.lang.OutOfMemoryError, set your environment variable explicitly: export _JAVA_OPTIONS="-Xmx16g -Djava.io.tmpdir=/tmp".
  • Index Timestamp Mismatch: If tools complain that a BAM index is missing or out of date, regenerate with samtools index -@ 4 sample.bam.
  • Conda Channel Conflicts: Avoid ABI incompatibilities between bioconda and conda-forge by setting strict priority: conda config --set channel_priority strict.
  • Chromosome Naming Contamination: Always verify whether reference FASTA headers use chr1 or 1; divergent naming causes silent alignment drops.
Was this guide helpful to your research?
Help us maintain 100% accurate, peer-reviewed bioinformatics standards.
📝 Propose Edit on GitHub ↗