omics.co.in logo omics.co.in

Gene & Isoform Quantification: featureCounts, HTSeq & RSEM

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
PLANT BIOTECHNOLOGY • TRACK 20 CHAPTER 04 OF 10
50 min read • 9 Code Examples • Subread & RSEM

Gene & Isoform Quantification: featureCounts, HTSeq & RSEM

Quantifying exon and gene-level reads using featureCounts and HTSeq; executing Bowtie2 + RSEM for transcript isoform estimation; comparing raw integer counts, FPKM, and TPM.

Phase: Day 2 • Expression Quantification

Key Learning Deliverables & Competencies

  • Running Subread featureCounts for paired-end plant RNA-Seq datasets
  • HTSeq-count configuration: union, intersection-strict, and strand-specific handling
  • Pipeline A: Bowtie2 + RSEM isoform-level quantification (rsem-calculate-expression)
  • Mathematical evaluation: raw counts vs FPKM vs TPM and why raw counts are mandatory for GLMs

1. Gene-Level Quantification: featureCounts vs HTSeq vs RSEM

Once reads are mapped, the next phase is summarizing aligned read fragments into a numeric count matrix $\mathbf{C} \in \mathbb{N}^{G \times S}$ (where $G$ is genes and $S$ is samples).

Quantification Metric Comparison

  • Raw Read Counts ($C_{ij}$): Exact integer counts of cDNA fragments assigned to gene $i$ in sample $j$. Mandatory input for DESeq2 and edgeR, because statistical power directly depends on the observed binomial count sampling depth.
  • FPKM / RPKM: Normalizes for gene length ($l_i$) and library size ($N_j$). Useful historically, but flawed for cross-sample comparison because the sum of FPKM varies between samples.
  • TPM (Transcripts Per Million): \[ \text{TPM}_i = \frac{\frac{C_{ij}}{l_i / 10^3}}{\sum_k \left(\frac{C_{kj}}{l_k / 10^3}\right)} \times 10^6 \] Guarantees $\sum_i \text{TPM}_i = 10^6$ in every sample. Best for comparing relative abundance of different genes within the same sample.

2. High-Throughput Counting with Subread featureCounts

#!/usr/bin/env bash
# Run featureCounts across all 18 coordinate-sorted BAM files
featureCounts -T 12 \
    -p --countReadPairs \
    -B -C \
    -t exon \
    -g gene_id \
    -a Solanum_tuberosum.DM_v6.1.gtf \
    -o potato_gene_counts.txt \
    bam_alignments/*.sorted.bam

3. Isoform Quantification: Bowtie2 + RSEM (Pipeline A)

When alternative splicing is the primary investigative target, RSEM (RNA-Seq by Expectation Maximization) estimates isoform-level abundance directly against transcript sequences:

#!/usr/bin/env bash
# Prepare RSEM reference
rsem-prepare-reference \
    --gtf Solanum_tuberosum.DM_v6.1.gtf \
    --bowtie2 \
    Solanum_tuberosum.DM_v6.1.dna.toplevel.fa \
    rsem_ref/potato_rsem

# Calculate expression with paired-end sequencing
rsem-calculate-expression --paired-end \
    --bowtie2 \
    -p 12 \
    trimmed_fastq/Tol_D3_R1_clean_R1.fastq.gz \
    trimmed_fastq/Tol_D3_R1_clean_R2.fastq.gz \
    rsem_ref/potato_rsem \
    rsem_out/Tol_D3_R1
⚠️

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 ↗