/featurecounts-rna-counting

name: "featurecounts-rna-counting"
description: "Counts RNA-seq reads overlapping GTF gene features. Takes sorted STAR BAMs plus GTF; outputs a per-gene tab-delimited matrix across samples. Handles strandedness (0/1/2), paired-end, multi-sample batch counting in one command, and outputs assignment statistics. Use Salmon for alignment-free quantification; use featureCounts when STAR BAMs already exist."
license: "GPL-3.0"

featureCounts — RNA-seq Read Counting

Overview

featureCounts (part of the Subread package) assigns sequencing reads in BAM files to genomic features defined in a GTF/GFF annotation. It counts how many reads overlap each gene (or exon, intron, or custom feature), producing a gene × sample count matrix suitable for differential expression analysis with DESeq2 or edgeR. featureCounts processes multiple BAM files in a single command, reporting read assignment statistics (assigned, unassigned by category) alongside the count matrix. It is the standard counting step after STAR alignment in RNA-seq pipelines.

When to Use

• Generating gene-level count matrices from STAR-aligned BAM files for DESeq2 or edgeR
• Counting reads from multiple samples simultaneously in a single featureCounts command
• Handling stranded RNA-seq libraries where sense/antisense assignment matters
• Producing exon-level or custom-feature counts (e.g., for splicing analysis with DEXSeq)
• Verifying strandedness of an RNA-seq library when protocol documentation is unavailable
• Use **Salmon** instead when no BAM file exists and fast pseudoalignment is preferred
• Use **HTSeq-count** as an alternative with slower but more flexible counting modes

Prerequisites

• **Software**: Subread package (contains `featureCounts`)
• **Input**: Sorted BAM files from STAR or HISAT2, plus a matching GTF annotation file

> **Check before installing**: The tool may already be available in the current environment (e.g., inside a `pixi` / `conda` env). Run `command -v featureCounts` first and skip the install commands below if it returns a path. When running inside a pixi project, invoke the tool via `pixi run featureCounts` rather than bare `featureCounts`.

# Install with conda (recommended)
conda install -c bioconda subread

# Verify
featureCounts -v
# featureCounts v2.0.6

# Alternative: install via apt (Ubuntu/Debian)
sudo apt-get install subread

Pre-flight Interview

Settle these with the user before writing any analysis code.

decisions:
  - id: D1
    param: annotation
    kind: required
    source: upstream
    ask: "Which GTF/GFF annotation file should be used, and does its assembly exactly match the BAM alignment reference?"
    default: null

  - id: D2
    param: strandSpecific
    kind: required
    source: data
    ask: "Which strandedness did the representative-alignment check infer, and does it agree with the library preparation?"
    default: null

  - id: D3
    param: -t, -g
    kind: optional
    source: user
    ask: "Should expression be summarized from exons to genes, or from another feature type for a specific downstream question?"
    default: "-t exon -g gene_id"

  - id: D4
    param: -M, -O
    kind: optional
    source: user
    ask: "Should reads with multiple genomic mappings or overlapping gene assignments be included in the count matrix?"
    default: "exclude both multi-mappers and multi-overlap reads"

  - id: D5
    param: --minOverlap, --fracOverlap
    kind: optional_conditional
    source: user
    ask: "Does the assay or feature definition require a minimum overlap stricter than one aligned base?"
    default: "--minOverlap 1 --fracOverlap 0"

Incorrect library strandedness can substantially reduce or misassign counts; the magnitude varies by gene and annotation, so infer it from the data and treat an ambiguous result as a QC failure. Enabling multi-mapper counting can inflate counts for paralogs and repeat-associated genes. Neither setting necessarily produces a software error.

Quick Start

# Count reads for multiple samples (unstranded paired-end RNA-seq)
featureCounts \
    -a gencode.v47.annotation.gtf \
    -o counts/gene_counts.txt \
    -T 8 \
    -p --countReadPairs \
    results/sample1/Aligned.sortedByCoord.out.bam \
    results/sample2/Aligned.sortedByCoord.out.bam

echo "Count matrix: counts/gene_counts.txt"
head -3 counts/gene_counts.txt

Workflow

Step 1: Prepare BAM Files and GTF

Ensure BAM files are sorted and indexed, and the GTF matches the genome assembly.

# Verify BAM files are sorted
samtools view -H results/sample1/Aligned.sortedByCoord.out.bam | grep "SO:"
# Expected: SO:coordinate

# List BAMs to count
ls results/*/Aligned.sortedByCoord.out.bam | head -5

# Download GENCODE GTF (same version used for STAR indexing)
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_47/gencode.v47.primary_assembly.annotation.gtf.gz
gunzip gencode.v47.primary_assembly.annotation.gtf.gz

echo "GTF lines: $(wc -l < gencode.v47.primary_assembly.annotation.gtf)"

Step 2: Determine Library Strandedness

Test strandedness using a small read count to set the `-s` parameter correctly.

# Quick strandedness check: count 1 sample with all 3 modes
# Compare assigned rates: highest = correct mode
for strand in 0 1 2; do
    echo "=== Strandedness -s $strand ==="
    featureCounts \
        -a gencode.v47.primary_assembly.annotation.gtf \
        -o /tmp/test_s${strand}.txt \
        -T 4 \
        -p --countReadPairs \
        -s $strand \
        results/sample1/Aligned.sortedByCoord.out.bam 2>&1 \
        | grep "Successfully assigned"
done
# Rule: 0=unstranded if similar rates; 1 or 2 if one is much higher

Step 3: Count Unstranded Paired-End RNA-seq

Standard configuration for unstranded libraries (most polyA-selected RNA-seq).

mkdir -p counts

# Multi-sample batch counting: pass all BAMs as positional arguments
featureCounts \
    -a gencode.v47.primary