Quantifying Gene Counts

In the previous lessons you turned FASTQ into coordinate-sorted BAMs with a splice-aware aligner. Differential expression tools need the next artifact: a matrix of counts per gene per sample. Two routes dominate, and this lesson runs both: count aligned reads over gene features with featureCounts, or quantify transcripts directly from FASTQ with Salmon.

Counting with featureCounts

featureCounts, part of the Subread package, assigns each aligned read or fragment to the gene it overlaps using the exon records in a GTF grouped by gene_id. It is fast and deterministic, and writes a plain count table plus a summary of how every read was categorized.

featureCounts \
  -T 4 \
  -p --countReadPairs \
  -s 2 \
  -t exon \
  -g gene_id \
  -a gencode.v44.annotation.gtf \
  -o counts.txt \
  sample.sorted.bam

The -p --countReadPairs pair counts one fragment per mate pair instead of one per mate, and -s sets strandedness: 0 unstranded, 1 forward, 2 reverse, where reverse matches the common dUTP stranded protocol. By default reads mapping to multiple loci or overlapping more than one gene are left unassigned, which is the safe choice for gene-level DE; add -M to count multimappers, --fraction to split their weight, and -O to allow overlap with several features. If strandedness is unknown, run -s 1 and -s 2 and compare their assignment rates: a large asymmetry points to the protocol strand (the higher one), while near-equal rates mean the library is unstranded, so use -s 0. RSeQC's infer_experiment.py or Salmon's -l A report the same information directly.

Status	sample.sorted.bam
Assigned	24853120
Unassigned_Unmapped	0
Unassigned_Read_Type	0
Unassigned_Singleton	412006
Unassigned_MappingQuality	0
Unassigned_Chimera	0
Unassigned_FragmentLength	0
Unassigned_Duplicate	0
Unassigned_MultiMapping	1830450
Unassigned_Secondary	0
Unassigned_NonSplit	0
Unassigned_NoFeatures	3120045
Unassigned_Overlapping_Length	0
Unassigned_Ambiguity	905233

Quantifying with Salmon

Salmon skips base-level alignment: it uses selective alignment against a transcriptome index and an expectation-maximization step to share reads across isoforms they map to equally well. It reads FASTQ directly, models multi-mapping probabilistically, and can correct fragment GC and sequence-composition bias once you enable --gcBias and --seqBias, which are both off by default. For production runs, build a decoy-aware index from the transcriptome plus the genome so reads from unannotated loci are not force-assigned to transcripts.

salmon index \
  -t gencode.v44.transcripts.fa.gz \
  -i salmon_index \
  -k 31 --gencode

salmon quant \
  -i salmon_index \
  -l A \
  -1 sample_R1.fastq.gz \
  -2 sample_R2.fastq.gz \
  -p 4 \
  -g gencode.v44.annotation.gtf \
  -o sample_salmon

Name	Length	EffectiveLength	TPM	NumReads
ENST00000456328.2	1657	1490.319	0.293921	9.000
ENST00000450305.2	632	465.117	0.000000	0.000
ENST00000488147.1	1351	1184.204	14.882733	362.114
ENST00000619216.1	68	24.000	0.000000	0.000
ENST00000473358.1	712	545.000	0.525728	5.887

Because -g was passed, Salmon also writes quant.genes.sf, summing transcript estimates to the gene level for you. When you need finer control you can import quant.sf into R with tximport and a two-column transcript-to-gene map, which also yields average transcript lengths that DESeq2 uses as offsets. Note that Salmon's NumReads are estimated fractional counts, so round or use tximport before feeding integer-count tools.