Count Matrix

A count matrix is a rectangular table whose rows are features (usually genes, sometimes transcripts or exons) and whose columns are samples, with each cell holding the number of sequencing reads assigned to that feature in that sample. It is the pivot point of RNA-seq: everything upstream is reads and alignment, everything downstream is statistics on this table.

How it works

Reads come off the instrument as FASTQ. Two routes produce counts.

The alignment route maps reads to the genome with a splice-aware aligner (STAR, HISAT2) and then counts overlaps against a GTF annotation with featureCounts or htseq-count. The relevant flags matter more than people expect: featureCounts -p --countReadPairs -t exon -g gene_id -s 2 -a gencode.v44.annotation.gtf. Getting -s wrong (strandedness) on a dUTP library silently halves your counts and skews antisense-overlapping genes. Check the .summary file for Unassigned_NoFeatures and Unassigned_Ambiguity; if either exceeds roughly 20 percent of assigned reads, something is misconfigured.

The lightweight route uses salmon quant or kallisto against a transcriptome index, producing transcript-level estimates that get summed to genes with tximport. We prefer this one. It handles multi-mapping reads across transcript isoforms probabilistically instead of discarding them, it runs in minutes on a laptop, and transcript-level estimates measurably improve gene-level inference because they correct for shifts in isoform usage that change effective gene length between samples ¹. Use tximport(files, type="salmon", tx2gene=t2g, countsFromAbundance="lengthScaledTPM") when feeding limma-voom, or the default when feeding DESeq2, which wants the length offset supplied separately.

The values are integers (or, from an EM quantifier, fractional estimates rounded at import). They are not expression levels. A count is proportional to abundance times transcript length times sequencing depth, which is why you never compare raw counts across genes or across samples.

In your own data

The file is usually a TSV or CSV with one header row, or an .rds/.h5ad if it came out of a pipeline. Load it and check four things before anything else.

Library size per column: colSums(counts). For bulk poly-A RNA-seq, 20 to 30 million assigned reads per sample is a working floor for gene-level work. A column at 4 million is a failed library, not a low-expression sample.

Detected genes per column: colSums(counts > 0). Expect 13,000 to 16,000 for whole blood, higher for tissue. A sharp drop in one sample usually tracks a low RIN.

RIN itself is not in the matrix but should be in your metadata. RIN runs 1 to 10. We treat 8 and above as clean for standard poly-A libraries, 7 to 8 as usable with degradation as a covariate, and below 7 as something where 3′ bias is real and you should either switch to a ribo-depletion or 3′-tag protocol or accept that long transcripts will read low. Put RIN in your design matrix if it varies across samples.

Row filtering: drop features with near-zero counts across all samples before testing. keep <- rowSums(counts >= 10) >= smallestGroupSize. This is not cosmetic. It reduces the multiple-testing burden and stabilizes dispersion estimation ².

Normalization: DESeq2’s median-of-ratios and edgeR’s TMM operate on the raw integer matrix and produce size factors, not a transformed matrix. Never hand TPM or FPKM to a count-based model. For visualization, use vst() or rlog(), which return log-scale values with variance roughly independent of the mean.

On TPM versus RPKM: RPKM divides by length then by total reads, so the resulting values do not sum to the same number in every sample, which means a proportion in one sample is not comparable to the same proportion in another. TPM reverses the order (length first, then normalize so each column sums to 1e6), which makes the units a genuine per-million-transcripts proportion and comparable across samples ³. Use TPM for cross-sample description, counts for testing.

Limitations

The matrix inherits every choice made upstream. Different annotation releases give different gene counts for the same reads. Different count-based methods disagree most on low-count genes and on small sample sizes, where dispersion is poorly estimated ⁴, and the ranking of methods depends on sample size and the size of the true effect ⁵.

A personal longitudinal matrix has n=1 per timepoint. You have no biological replicates, so the dispersion a package estimates is technical plus within-person temporal variance combined. Treat fold changes as descriptive trajectories rather than tested effects, and require a gene to move consistently across three or more timepoints before you take it seriously.

Single-cell matrices are a different object: sparse, UMI-deduplicated, cells rather than samples in columns, with 90 percent or more zeros. Quality control thresholds, normalization, and doublet handling all differ, and the sequence of steps changes the result ⁶.

Nothing in a count matrix is clinical. Expression changes are not diagnoses. Take anything that looks medically meaningful to a physician.

Woolf Software builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.