bio-expression-matrix-metadata-joins

What this skill does

## Version Compatibility

Reference examples tested with: pandas 2.2+, DESeq2 1.42+, edgeR 4.0+, limma 3.58+, variancePartition / dream 1.32+, somalier 0.2.18+ (CLI), pyensembl 2.3+, anndata 0.10+

Before using code patterns, verify installed versions match. If versions differ:
- R: `packageVersion('<pkg>')` then `?function_name` to verify parameters
- Python: `pip show <package>` then `help(module.function)` to check signatures
- 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.

# Metadata Joins

**"Align my sample metadata with my count matrix and build a design"** -> Reconcile sample identifiers, validate the join, set factor levels explicitly to control fold-change direction, encode the experimental structure (paired, repeated measures, interaction) in the design formula, and detect swaps before downstream DE.

## The Single Most Important Modern Insight -- Alphabetical reference levels invert fold changes silently

DESeq2 picks the reference level alphabetically if not told otherwise. With `condition = c('Treated', 'Untreated')`, `T < U`, so `Treated` becomes the reference. The reported `log2FoldChange` is then `Untreated vs Treated` -- **the opposite of what the methods section says**. No error; the volcano plot looks plausible; the gene list is correct but with reversed sign.

```r
dds$condition <- relevel(dds$condition, ref = 'Untreated')
coldata$condition <- factor(coldata$condition, levels = c('Untreated', 'Treated'))
```

Set BEFORE `DESeq()`; relevel-then-DESeq again to take effect.

A second insight: in interaction designs `~ genotype * treatment`, `results(dds, name='treatment_drug_vs_vehicle')` returns the drug effect IN THE WT REFERENCE only, not the marginal/average effect. The interaction coefficient `genotypeKO.treatmentdrug` is the DIFFERENCE in drug effect between KO and WT (the difference of differences), NOT the drug effect in KO. The cleanest fix is to use `~ 0 + group` with `group = paste(genotype, treatment, sep='_')` so every comparison of interest is a single contrast.

## Algorithmic Taxonomy

| Pattern | Encoding | Tests | Caveat |
|---------|----------|-------|--------|
| Simple two-group | `~ condition` | `results(name='condition_treated_vs_control')` | Set reference level explicitly |
| With known batch | `~ batch + condition` | `results(name='condition_...')` | Variable of interest LAST is the convention; not required |
| Paired (pre/post per subject) | `~ subject + condition` | `results(name='condition_...')` | Pairing FIRST; subject absorbs baseline variability |
| Repeated measures (>2 time per subject) | DREAM `~ condition + (1|subject)` | `topTable(coef='condition')` | Mixed model; uses voom + lmer |
| Interaction (2x2) | `~ A * B` (expands to A + B + A:B) | `results(name='A.B')` for diff-in-diff; combined factor for cleaner contrasts | resultsNames trap |
| Multi-group all pairwise | `~ 0 + group` + `makeContrasts` | Contrasts named directly | DESeq2 needs intercept; works for edgeR/limma |
| Multi-level "any change" | LRT `reduced = ~ 1` | `results(dds)` from LRT fit | padj is omnibus; LFC is one specific level |

## Decision Tree by Scenario

| Scenario | Recommended approach |
|----------|---------------------|
| Tumor / normal paired by patient | `~ patient + tissue`; pairing FIRST -- absorbs subject variability |
| Pre / post drug, same patient (2 time points) | `~ subject + condition` |
| Longitudinal, 4+ time points per subject | DREAM mixed model with `(1 \| subject)` random effect |
| Known batch (sequencing run, library prep date) | `~ batch + condition`; do NOT subtract |
| Many batches (>10 levels, n=20-40) | `~ condition + (1 \| batch)` via DREAM; or aggregate batches |
| Continuous covariate (age, RIN) | Center first; include linearly OR via `ns(x, df=3)` if non-linear |
| Mixed-sex cohort | Include sex unless sex-specific; report sex-stratified sensitivity |
| Multi-group, want pairwise contrasts | `~ 0 + group` + `makeContrasts` |
| Interaction question (does effect of A differ across B?) | Either `~ A * B` (handle resultsNames trap) or combined factor `~ 0 + group` |
| Many technical reps within bio reps | `duplicateCorrelation` (limma) OR aggregate via `collapseReplicates` (DESeq2) |
| Cohort >=20 samples | Run somalier or NGSCheckMate genotype check at the matrix-build step |

## Basic Join

**Goal:** Align count matrix columns with metadata rows; remove samples present in only one source; verify alignment before downstream use.

**Approach:** Intersect sample identifiers; reorder both data sources; assert match.

```python
import pandas as pd

counts = pd.read_csv('counts.tsv', sep='\t', index_col=0)
metadata = pd.read_csv('metadata.csv', index_col=0)

common = counts.columns.intersection(metadata.index)
only_counts = set(counts.columns) - set(metadata.index)
only_meta = set(metadata.index) - set(counts.columns)
if only_counts: print(f'In counts not metadata: {only_counts}')
if only_meta: print(f'In metadata not counts: {only_meta}')

counts = counts[common]
metadata = metadata.loc[common]
assert all(counts.columns == metadata.index)
```

```r
common <- intersect(colnames(counts), rownames(coldata))
counts <- counts[, common]
coldata <- coldata[common, , drop = FALSE]
stopifnot(all(colnames(counts) == rownames(coldata)))
```

When sample names differ by formatting (underscore vs dash, BAM suffix, case), try systematic transformations -- replace `_` <-> `-`, strip `.bam`, lower-case, take prefix before `_` -- before giving up.

## Reference Level (Critical -- Set Before DESeq)

**Goal:** Pick the baseline against which fold changes are reported, controlling sign direction.

**Approach:** `relevel()` or construct the factor with explicit `levels=`. Set BEFORE `DESeq()` runs.

```r
coldata$condition <- factor(coldata$condition, levels = c('control', 'treated'))

coldata$condition <- relevel(coldata$condition, ref = 'control')
```

```python
metadata['condition'] = pd.Categorical(metadata['condition'],
                                        categories=['control', 'treated'],
                                        ordered=True)
```

With factor levels explicit, the LFC reads `treated / control` -- treated up means LFC > 0.

A reminder for edgeR / limma users: `makeContrasts(Treated - Control)` is explicit and immune to the reference-level trap. Same caution applies less because the user names the contrast.

## Paired Designs

```r
design = ~ patient + tissue
```

Pairing variable FIRST (convention). Patient absorbs inter-subject baseline variability, dramatically increasing power for the tissue effect.

For DESeq2:
```r
dds <- DESeqDataSetFromMatrix(counts, coldata, design = ~ patient + tissue)
dds <- DESeq(dds)
res <- results(dds, name = 'tissue_tumor_vs_normal')
```

For edgeR:
```r
design <- model.matrix(~ patient + tissue, coldata)
y <- estimateDisp(y, design, robust = TRUE)
fit <- glmQLFit(y, design, robust = TRUE)
qlf <- glmQLFTest(fit, coef = 'tissuetumor')
```

Common mistake: writing `~ tissue + patient`. Numerically the model is the same; the convention of pairing-first improves readability and matches the natural mental model.

## Interaction Terms -- the resultsNames Trap

```r
design = ~ genotype + treatment + genotype:treatment
dds <- DESeqDataSetFromMatrix(counts, coldata, design = design)
dds <- DESeq(dds)
resultsNames(dds)
```

Output names (after relevel):
```
"Intercept"
"genotype_KO_vs_WT"
"treatment_drug_vs_vehicle"
"genotypeKO.treatmentdrug"
```

| Question | Wrong answer | Right answer |
|----------|--------------|--------------|
| Drug effect averaged over genotypes | `results(name='treatment_drug_vs_vehicle')` -- this is drug effect in WT only | Combined factor `~ 0 + group`; or contrast that explicitly averages |
| Is drug effect different between genotypes? | n/a | `results(name='genotypeKO.trea