Bio Expression Matrix Metadata Joins

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.

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

Decision Tree by Scenario

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.

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)

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.

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

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

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

design = ~ patient + tissue

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

For DESeq2:

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

For edgeR:

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

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"

The cleaner alternative for designs with many contrasts of interest:

coldata$group <- factor(paste(coldata$genotype, coldata$treatment, sep = '_'))
dds <- DESeqDataSetFromMatrix(counts, coldata, design = ~ 0 + group)
dds <- DESeq(dds)

res_drug_in_ko <- results(dds, contrast = c('group', 'KO_drug', 'KO_vehicle'))
res_drug_in_wt <- results(dds, contrast = c('group', 'WT_drug', 'WT_vehicle'))
res_diff       <- results(dds, contrast = list(c('groupKO_drug', 'groupWT_vehicle'),
                                                c('groupKO_vehicle', 'groupWT_drug')))

~ 0 + group parameterization is the long-standing edgeR / limma recommendation (Smyth and Robinson User's Guides) for any design with multiple pairwise contrasts of interest.

LRT and the Reduced Model

dds <- DESeq(dds, test = 'LRT', reduced = ~ batch)

Reduced model drops the term being tested. With design = ~ batch + condition and reduced = ~ batch, the LRT tests condition. With interaction designs:

dds <- DESeq(dds, test = 'LRT', reduced = ~ genotype + treatment)

(Tests the interaction.)

The reduced model must be NESTED in the full model (every term in reduced must appear in full).

Continuous Covariates

coldata$age_c <- coldata$age - mean(coldata$age)
design = ~ age_c + RIN + condition

library(splines)
design = ~ ns(age, df = 3) + condition

DO NOT include library size as a covariate -- it is handled by size factors / normalization factors internally.

Repeated Measures -- duplicateCorrelation vs dream

Goal: Correctly model within-subject correlation when the same subject contributes multiple samples.

Approach: For technical reps within bio reps OR paired pre/post: duplicateCorrelation (limma) is adequate. For >2 time points per subject or random slopes: DREAM (variancePartition).

library(limma)
library(edgeR)

v <- voom(y, design)
corfit <- duplicateCorrelation(v, design, block = coldata$donor)
v <- voom(y, design, block = coldata$donor, correlation = corfit$consensus)
corfit <- duplicateCorrelation(v, design, block = coldata$donor)
fit <- lmFit(v, design, block = coldata$donor, correlation = corfit$consensus)
fit <- eBayes(fit, robust = TRUE)

The double pass is intentional: estimate correlation, re-voom with correlation, re-estimate, fit. Limma's duplicateCorrelation assumes ONE within-subject correlation across all genes -- approximation.

For proper per-gene mixed models:

library(variancePartition)

form <- ~ condition + (1 | donor)
vobj <- voomWithDreamWeights(y, form, coldata)
fitmm <- dream(vobj, form, coldata)
fitmm <- eBayes(fitmm)
tt <- topTable(fitmm, coef = 'condition')

See differential-expression/timeseries-de for full longitudinal designs.

Before committing to a design, variancePartition::fitExtractVarPartModel(vobj, form, coldata) quantifies the fraction of expression variance explained by each covariate, gene-by-gene. If a candidate "nuisance" covariate explains <1% of variance across most genes, it can usually be dropped; if a known biological factor explains <5% and isn't of direct interest, model it as a random effect rather than fixed.

Sample Swap Detection (Mandatory for Cohort >= 20)

Sex check via XIST and chrY expression

Goal: Detect mislabeled samples by checking that gene-expression sex matches reported sex.

Approach: Compare XIST expression (high in XX, low/absent in XY) against chrY-gene expression (DDX3Y, RPS4Y1, UTY, KDM5D, EIF1AY -- high in XY, absent in XX).

import pandas as pd

def sex_check(counts, metadata, sex_column='sex'):
    y_genes = ['DDX3Y', 'RPS4Y1', 'UTY', 'KDM5D', 'EIF1AY']
    y_avail = [g for g in y_genes if g in counts.index]
    if 'XIST' not in counts.index or not y_avail:
        return None
    predicted = pd.Series('unknown', index=counts.columns)
    predicted[counts.loc[y_avail].sum() > counts.loc['XIST']] = 'M'
    predicted[counts.loc['XIST'] > counts.loc[y_avail].sum()] = 'F'
    if sex_column in metadata.columns:
        mis = predicted != metadata[sex_column]
        if mis.any():
            print(f'SEX MISMATCHES: {list(metadata.index[mis])}')
    return predicted

sex_check <- function(counts, coldata, sex_col = 'sex') {
    y_genes <- c('DDX3Y', 'RPS4Y1', 'UTY', 'KDM5D', 'EIF1AY')
    y_expr <- colSums(counts[intersect(y_genes, rownames(counts)), , drop = FALSE])
    predicted <- ifelse(y_expr > counts['XIST', ], 'M', 'F')
    mis <- predicted != coldata[[sex_col]]
    if (any(mis)) cat('Sex mismatches:', colnames(counts)[mis], '\n')
    predicted
}

CAVEAT: tumors with X loss, sex chromosome aneuploidies, HeLa (XXX with mixed inactivation) muddy this. Genotype-based methods are more robust.

Genotype-based: somalier and NGSCheckMate

somalier extract -d extracted/ --sites sites.GRCh38.vcf.gz \
    -f reference.fa sample.bam
somalier relate --infer extracted/*.somalier

ncm_fastq.py -l fastq_list.txt -O outdir -bed common_sites.bed

Somalier (Pedersen et al. 2020 Genome Med 12:62) extracts a few thousand SNP sketches per sample (sub-second per sample) and computes pairwise relatedness from BAM/CRAM/VCF. NGSCheckMate (Lee et al. 2017 NAR 45:e103) computes VAF correlation across a common-SNP panel; works on FASTQ/BAM/VCF including RNA-seq.

For any cohort >=20 samples, run one of these at the matrix-build step. Catching a swap in raw data is cheap; finding it after DE is expensive.

SABV -- Sex as Biological Variable

NIH 2016+ requires sex consideration in vertebrate animal and human studies. Mauvais-Jarvis F et al. 2020 Lancet 396:565 reviews effect-size differences across diseases.

Practical implication:

• Always include sex in the metadata, even when not in the model.
• For sex-balanced cohorts, ~ sex + condition rarely hurts and captures real biology.
• For sex-confounded cohorts (all-male disease cohort vs mixed-sex control), can't rescue but documents the limitation.
• For chrX/chrY analyses, sex MUST be in the model OR the analysis is uninterpretable.
• Report DE counts overall AND sex-stratified.

Simpson's Paradox -- Collapsing Technical Replicates

Goal: Aggregate technical replicates from the same subject correctly; never treat them as independent biological replicates.

Approach: Sum (not average) technical replicates of the same subject BEFORE downstream DE.

library(DESeq2)
dds_collapsed <- collapseReplicates(dds, groupby = dds$subject)

counts_per_subject = counts.T.groupby(metadata['subject']).sum().T
metadata_per_subject = metadata.drop_duplicates(subset='subject').set_index('subject')

Why sum and not average? Reads add. Two technical replicates yielding 1M reads each are equivalent to one library yielding 2M reads. Averaging would understate the effective library size.

Treating technical replicates as independent biological samples is the cardinal sin: it inflates the apparent sample size and deflates standard errors. With 4 patients x 3 tech reps = 12 samples, naive DE assumes 12 independent observations; the truth is closer to 4. p-values are compressed ~3x.