Pydeseq2 Differential Expression

pip install pydeseq2 matplotlib seaborn

decisions:
  - id: D1
    param: design
    kind: required
    source: data
    ask: "Which column of the sample sheet separates the groups you want to compare?"
    default: null

  - id: D2
    param: contrast
    kind: required
    source: user
    depends_on: [D1]
    ask: "Within that column, which group is the baseline the others are measured against?"
    default: null

  - id: D3
    param: design
    kind: required
    source: data
    depends_on: [D1]
    ask: "Is there a nuisance variable to adjust for (batch, donor, sex, sequencing run)?"
    default: "none - compare groups without covariates"
    skip_if: "sample sheet carries no column besides the grouping variable"

  - id: D4
    param: alpha
    kind: optional
    source: user
    ask: "How strong must the evidence be before a gene counts as changed (false-discovery rate)?"
    default: 0.05

  - id: D5
    param: lfc_cutoff
    kind: optional
    source: user
    ask: "How large a fold change is worth reporting, on top of the statistical call?"
    default: "|log2FC| > 0.58 (1.5x change)"

  - id: D6
    param: min_total_counts
    kind: optional_conditional
    source: data
    ask: "Genes with very few reads cannot be tested - where should the floor sit?"
    default: 10

  - id: D7
    param: lfc_shrink
    kind: optional
    source: user
    ask: "Shrink fold-change estimates so low-count genes stop dominating the rankings and plots?"
    default: "applied for ranking and visualization, not for significance calls"

  - id: D8
    param: refit_cooks, cooks_filter, independent_filter
    kind: optional
    source: user
    ask: "Should DESeq2 retain its default outlier replacement and adaptive filtering behavior?"
    default: true

  - id: D9
    param: n_cpus
    kind: never_ask
    source: data
    reason: "Affects runtime only, not the result"
    default: 4

import pandas as pd

# Load data — typical CSV has genes as rows, samples as columns
counts_raw = pd.read_csv("counts.csv", index_col=0)
metadata = pd.read_csv("metadata.csv", index_col=0)

# Transpose if needed: PyDESeq2 requires samples x genes
if counts_raw.shape[0] > counts_raw.shape[1]:
    counts_df = counts_raw.T  # genes x samples → samples x genes
else:
    counts_df = counts_raw

# Validate alignment
common_samples = counts_df.index.intersection(metadata.index)
counts_df = counts_df.loc[common_samples]
metadata = metadata.loc[common_samples]

print(f"Samples: {counts_df.shape[0]}, Genes: {counts_df.shape[1]}")
print(f"Metadata columns: {list(metadata.columns)}")
print(f"Condition counts:\n{metadata['condition'].value_counts()}")

# Filter genes with total counts below threshold
min_total_counts = 10
gene_counts = counts_df.sum(axis=0)
genes_to_keep = gene_counts[gene_counts >= min_total_counts].index
counts_df = counts_df[genes_to_keep]

# Optional: require minimum counts in a minimum number of samples
min_count_per_sample = 5
min_samples = 3
genes_expressed = (counts_df >= min_count_per_sample).sum(axis=0) >= min_samples
counts_df = counts_df.loc[:, genes_expressed]

print(f"Genes after filtering: {counts_df.shape[1]}")

from pydeseq2.dds import DeseqDataSet

dds = DeseqDataSet(
    counts=counts_df,
    metadata=metadata,
    design="~condition",   # Wilkinson-style formula
    refit_cooks=True,      # Refit after Cook's outlier removal
    n_cpus=4               # Parallel threads
)

# Run: size factors → dispersions → trend → MAP shrinkage → LFC fitting
dds.deseq2()

# Inspect normalization
print(f"Size factors (first 5): {dds.obsm['size_factors'][:5]}")
print(f"Size factor range: {dds.obsm['size_factors'].min():.2f} - {dds.obsm['size_factors'].max():.2f}")

from pydeseq2.ds import DeseqStats

ds = DeseqStats(
    dds,
    contrast=["condition", "treated", "control"],
    alpha=0.05,              # FDR threshold
    cooks_filter=True,       # Filter Cook's outliers
    independent_filter=True  # Independent filtering for power
)

ds.summary()

# Access full results
results = ds.results_df
print(f"Total genes tested: {len(results)}")
print(f"Significant (padj < 0.05): {(results.padj < 0.05).sum()}")

# Apply shrinkage — modifies results_df.log2FoldChange in place
ds.lfc_shrink()

# Compare pre/post shrinkage effect
print(f"Max |LFC| after shrinkage: {results.log2FoldChange.abs().max():.2f}")
print(f"Genes with |LFC| > 2: {(results.log2FoldChange.abs() > 2).sum()}")

import numpy as np

# Significance + effect size filter
significant = results[
    (results.padj < 0.05) &
    (results.log2FoldChange.abs() > 1.0)
].copy()

# Separate up/down-regulated
up = significant[significant.log2FoldChange > 0].sort_values("padj")
down = significant[significant.log2FoldChange < 0].sort_values("padj")

print(f"Upregulated: {len(up)}, Downregulated: {len(down)}")

# Export
results.to_csv("deseq2_all_results.csv")
significant.to_csv("deseq2_significant.csv")
up.to_csv("deseq2_upregulated.csv")
down.to_csv("deseq2_downregulated.csv")
print("Results exported to CSV files")

# Volcano input = log2FoldChange vs padj; saved in Step 6 as deseq2_all_results.csv
volcano_input = results[["log2FoldChange", "padj"]].dropna()
print(f"Volcano input: {len(volcano_input)} genes -> deseq2_all_results.csv")
# Render with the omics-plotting SKILL (`skills/data-visualization/omics-plotting/SKILL.md`) "Volcano" recipe.

# MA input = log2FoldChange vs mean normalized count (baseMean); from deseq2_all_results.csv
ma_input = results[["baseMean", "log2FoldChange", "padj"]].dropna(subset=["padj"])
print(f"MA input: {len(ma_input)} genes -> deseq2_all_results.csv")
# Render with the omics-plotting SKILL (`skills/data-visualization/omics-plotting/SKILL.md`) "MA plot" recipe.