bio-expression-matrix-counts-ingest

What this skill does

## Version Compatibility

Reference examples tested with: pandas 2.2+, numpy 1.26+, scanpy 1.10+, anndata 0.10+, tximport 1.30+, tximeta 1.20+, GenomicFeatures 1.54+, Subread/featureCounts 2.0.6+ (post-v2.0.2 API), STAR 2.7.10+, Salmon 1.10+, kallisto 0.48+, RSEM 1.3.3+, HTSeq 2.0+

Before using code patterns, verify installed versions match. If versions differ:
- Python: `pip show <package>` then `help(module.function)` to check signatures
- R: `packageVersion('<pkg>')` then `?function_name` to verify parameters
- CLI: `<tool> --version` then `<tool> --help` to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.

# Count Matrix Ingestion

**"Load my featureCounts / Salmon / STAR output into a count matrix"** -> Parse the per-sample quantification, strip metadata, choose the correct strandedness or NumReads column, optionally apply length-bias correction via tximport, and return a gene-by-sample matrix appropriate for the downstream DE tool.

## The Single Most Important Modern Insight -- `lengthScaledTPM` is NOT TPM; the naming has fooled many

`tximport(..., countsFromAbundance='lengthScaledTPM')` returns a **count-scale matrix** (sum to library size) with the gene-length bias removed, intended as input to DE tools that cannot accept offsets (limma-voom). The word "TPM" in the option name has misled users into reporting these values as normalized abundance -- they are not.

The decision tree for `countsFromAbundance` (Soneson, Love, Robinson 2015 *F1000Res* 4:1521):

| Option | Returns | Use with |
|--------|---------|----------|
| `'no'` (default) | Raw NumReads; length matrix passed as DESeq2/edgeR offset | DESeq2 via `DESeqDataSetFromTximport`, edgeR via `DGEList` -- the correct path for DGE |
| `'scaledTPM'` | TPM scaled to library size | DGE when downstream tool cannot accept offsets (rare) |
| `'lengthScaledTPM'` | TPM scaled by avg transcript length AND library size | limma-voom for DGE; recommended modern default when offsets unavailable |
| `'dtuScaledTPM'` | Per-transcript scaling by median isoform length | Differential Transcript Usage (DRIMSeq, DEXSeq) ONLY; requires `txOut=TRUE` |

Two adjacent traps with the same flavor of silent miscount:

1. **featureCounts `-p` since v2.0.2 (March 2021) needs `--countReadPairs`**. Pre-v2.0.2, `-p` flagged paired-end AND counted pairs as fragments. Post-v2.0.2, `-p` only flags paired-end -- pairs are NOT counted as one fragment unless `--countReadPairs` is added. Old scripts run on new installs produce ~2x the expected counts (one count per mate). No warning. Always pass `-p --countReadPairs` together for paired-end.

2. **STAR `ReadsPerGene.out.tab` column choice depends on library protocol**. Column 2 is unstranded; column 3 is forward-stranded; column 4 is reverse-stranded. Illumina TruSeq Stranded (the dominant kit, dUTP-based) is reverse-stranded -- column 4. Reading column 3 for TruSeq throws away ~95% of reads. Reading column 2 conflates antisense expression. Verify with RSeQC `infer_experiment.py` on a few BAMs before assembling the matrix.

## Algorithmic Taxonomy

| Source | Output structure | Read into | Caveats |
|--------|------------------|-----------|---------|
| featureCounts (Liao 2014) | TSV with 6 metadata cols + N sample cols | `read.delim` (R), `pd.read_csv(comment='#')` (Python) | -p/--countReadPairs v2.0.2 break; -O double-counts; -s strandedness |
| HTSeq-count (Anders 2015) | Per-sample 2-col TSV with `__` summary lines | Per-sample read, drop `__` rows, concat | `--mode` matters (union vs intersection-strict vs intersection-nonempty); `-s` strandedness |
| STAR `--quantMode GeneCounts` | Per-sample 4-col `ReadsPerGene.out.tab` | Skip first 4 summary rows; pick column by strandedness | Column 4 = TruSeq Stranded |
| Salmon `quant.sf` | Per-sample 5-col TSV (Name, Length, EffectiveLength, TPM, NumReads) | tximport (recommended) or manual NumReads sum (length-biased) | Selective alignment is the default since Salmon 1.0.0 (Srivastava 2020); the `--validateMappings` flag is now a no-op |
| kallisto `abundance.tsv` / `.h5` | Per-sample 5-col TSV; bootstraps in `.h5` | tximport (gene-level) or sleuth (transcript-level with bootstrap variance) | `kallisto quant -b 100` for sleuth |
| RSEM `*.genes.results` / `*.isoforms.results` | Per-sample TSV; expected_count is non-integer | tximport (`type='rsem'`) with `round()` via DESeq2; or manual | Zero-length transcripts cause `lengths > 0 is not TRUE` |
| 10X Genomics CellRanger | filtered_feature_bc_matrix/ MTX dir OR `.h5` | scanpy `read_10x_mtx` / `read_10x_h5`; Seurat `Read10X` | Single-cell convention: cells in rows |
| AnnData `.h5ad` | scverse single-file binary | scanpy `read_h5ad`; in R via zellkonverter | `.X` vs `.layers['counts']` vs `.raw.X` semantics |
| RDS (from R) | R-serialized object | `pyreadr` (Python); base R `readRDS` | For Seurat objects, convert first |

## Decision Tree by Scenario

| Scenario | Recommended approach |
|----------|---------------------|
| Bulk RNA-seq with featureCounts paired-end | `featureCounts -p --countReadPairs -s 2 -t exon -g gene_id` (TruSeq Stranded) |
| Bulk RNA-seq with STAR | `--quantMode GeneCounts`; read column 4 for TruSeq Stranded |
| Salmon/kallisto -> gene-level DGE with DESeq2 | `tximport(type='salmon', tx2gene)` + `DESeqDataSetFromTximport()` -- offsets handled |
| Salmon -> gene-level DGE with limma-voom | `tximport(..., countsFromAbundance='lengthScaledTPM')` -- no offset |
| Salmon/kallisto -> DTE (differential transcript expression) | `tximport(..., txOut=TRUE, countsFromAbundance='dtuScaledTPM')` for DRIMSeq/DEXSeq; OR `edgeR::catchSalmon()` for the Baldoni 2024 framework |
| kallisto with bootstraps -> sleuth (uncertainty-aware DE) | `sleuth_prep()` directly; do not go through tximport |
| RSEM -> DESeq2 | `tximport(files, type='rsem', txIn=FALSE)` + `DESeqDataSetFromTximport()` |
| Want automatic annotation provenance | `tximeta` (Love 2020 *PLoS Comp Biol* 16:e1007664) |
| 3'-tagged library (10x bulk, QuantSeq) | `countsFromAbundance='no'` WITHOUT length offset -- length bias negligible |
| 10X single-cell | `sc.read_10x_h5()` or `Read10X_h5()` |
| Strandedness unknown | RSeQC `infer_experiment.py` on 1-2 BAMs BEFORE re-running quantification |

## featureCounts -- The `-p` and `-O` Traps

**Goal:** Run featureCounts correctly on paired-end stranded RNA-seq, then read the output without inheriting the metadata columns.

**Approach:** CLI invocation with `-p --countReadPairs -s <strandedness>`, plus optional `-O --fraction` for overlapping-gene handling; parse with pandas/read.delim stripping the 6 metadata columns.

```bash
featureCounts -T 8 -p --countReadPairs -s 2 \
    -t exon -g gene_id \
    -a annotation.gtf \
    -o featurecounts.txt \
    sample1.bam sample2.bam sample3.bam
```

```python
import pandas as pd

fc = pd.read_csv('featurecounts.txt', sep='\t', comment='#')
counts = fc.set_index('Geneid').iloc[:, 5:]
counts.columns = [c.replace('.bam', '').split('/')[-1] for c in counts.columns]
```

```r
fc <- read.delim('featurecounts.txt', comment.char = '#', row.names = 1)
counts <- fc[, 6:ncol(fc)]
colnames(counts) <- gsub('.*/|\\.bam$', '', colnames(counts))
```

| Flag | Meaning | Default | When to flip |
|------|---------|---------|--------------|
| `-p` | Input is paired-end | off | Always for paired-end -- AND add `--countReadPairs` |
| `--countReadPairs` | Count fragment (pair) as one | off in v2.0.2+ | Always for paired-end (post-v2.0.2 API change) |
| `-s 0|1|2` | Strandedness | 0 (unstranded) | `-s 2` for TruSeq Stranded; `-s 1` for forward kits |
| `-O` | Allow multi-overlap | off | Off for typical DGE; on creates double-counts in overlapping genes |
| `-M` | Count multi-mappers | off | Off for DGE; on with `--fraction` for fractional counting |
| `-t` | Feature type | `exon` | Almost always `exon` |
| `-g` | Group attribute |