What does the Pydeseq2 AI skill do?

Differential gene expression analysis for bulk RNA-seq with PyDESeq2, including formulaic designs, Wald tests, FDR correction, LFC shrinkage, and result visualization.

Why use Pydeseq2 on TypingMind?

Because you install it once and use it with any model. Pydeseq2 is plain Markdown rather than provider-specific code, so the same skill runs on GPT-5, Claude, Gemini, Grok, or a local model — and you can switch model mid-chat without it breaking. TypingMind runs on your own API keys, so you pay providers directly instead of a per-seat subscription, and your skills and chats stay in your own storage.

How do I install Pydeseq2 in TypingMind?

Open Plugins → Skills → Install from GitHub in TypingMind and paste https://github.com/dralkh/seerai/tree/main/skills/pydeseq2. TypingMind reads its SKILL.md and bundles its files and installs it as a skill you can enable per chat.

Which AI models can use Pydeseq2?

Any model you connect in TypingMind. AI skills are plain Markdown instructions rather than provider-specific code, so GPT, Claude, Gemini, Grok, and local models can all load this skill when a request matches it.

How many AI models can I use with Pydeseq2?

As many as you like. As long as a model supports skills, you can use Pydeseq2 with it — GPT, Claude, Gemini, Grok, DeepSeek, Mistral, Llama and more — all on TypingMind with your own API keys.

Is the Pydeseq2 AI skill free?

Yes. It is published on GitHub by dralkh under the MIT license. You only pay your own AI provider for the tokens you use.

An AI skill is a reusable instruction bundle that teaches an AI model how to do one specific task. It follows the open Agent Skills format: a SKILL.md file with a name and description, plus any scripts, templates or reference files the model may need. The model reads the instructions only when your request matches the skill, so an installed skill costs nothing until it is used.

How are AI skills different from plugins or MCP servers?

A plugin or MCP server gives a model new tools to call — code that runs somewhere and returns a result. An AI skill gives the model knowledge and process instead: how to approach a task, which steps to follow, what good output looks like. Skills are plain Markdown, so they need no server, no API key and no runtime, and they work with any model.

Pydeseq2 by dralkh — AI Skill for Any LLM

PyDESeq2

Overview

PyDESeq2 is a Python implementation of DESeq2 for differential expression analysis with bulk RNA-seq data. Design and execute complete workflows from data loading through result interpretation, including formulaic single-factor and multi-factor designs, Wald tests with multiple testing correction, optional apeGLM shrinkage, and integration with pandas and AnnData.

When to Use This Skill

This skill should be used when:

Analyzing bulk RNA-seq count data for differential expression
Comparing gene expression between experimental conditions (e.g., treated vs control)
Performing multi-factor designs accounting for batch effects or covariates
Converting R-based DESeq2 workflows to Python
Integrating differential expression analysis into Python-based pipelines
Users mention "DESeq2", "differential expression", "RNA-seq analysis", or "PyDESeq2"

Quick Start Workflow

For users who want to perform a standard differential expression analysis:

python
import pandas as pd
from pydeseq2.dds import DeseqDataSet
from pydeseq2.default_inference import DefaultInference
from pydeseq2.ds import DeseqStats

# 1. Load data
counts_df = pd.read_csv("counts.csv", index_col=0).T  # Transpose to samples × genes
metadata = pd.read_csv("metadata.csv", index_col=0)

# 2. Filter low-count genes
genes_to_keep = counts_df.columns[counts_df.sum(axis=0) >= 10]
counts_df = counts_df[genes_to_keep]

# 3. Make the reference level explicit and fit DESeq2
metadata["condition"] = pd.Categorical(
    metadata["condition"], categories=["control", "treated"]
)
inference = DefaultInference(n_cpus=4)
dds = DeseqDataSet(
    counts=counts_df,
    metadata=metadata,
    design="~condition",
    refit_cooks=True,
    inference=inference,
)
dds.deseq2()

# 4. Perform statistical testing
ds = DeseqStats(
    dds,
    contrast=["condition", "treated", "control"],
    inference=inference,
)
ds.summary()

# 5. Access results
results = ds.results_df
significant = results[results.padj < 0.05]
print(f"Found {len(significant)} significant genes")

Core Workflow Steps

Step 1: Data Preparation

Input requirements:

Count matrix: Samples × genes DataFrame with non-negative integer read counts
Metadata: Samples × variables DataFrame with experimental factors

Common data loading patterns:

python
# From CSV (typical format: genes × samples, needs transpose)
counts_df = pd.read_csv("counts.csv", index_col=0).T
metadata = pd.read_csv("metadata.csv", index_col=0)

# From TSV
counts_df = pd.read_csv("counts.tsv", sep="\t", index_col=0).T

# From AnnData
import anndata as ad
adata = ad.read_h5ad("data.h5ad")
counts_df = pd.DataFrame(adata.X, index=adata.obs_names, columns=adata.var_names)
metadata = adata.obs

Data filtering:

python
# Remove low-count genes
genes_to_keep = counts_df.columns[counts_df.sum(axis=0) >= 10]
counts_df = counts_df[genes_to_keep]

# Remove samples with missing metadata
samples_to_keep = ~metadata.condition.isna()
counts_df = counts_df.loc[samples_to_keep]
metadata = metadata.loc[samples_to_keep]

Step 2: Design Specification

The design formula specifies how gene expression is modeled.

Single-factor designs:

python
design = "~condition"  # Simple two-group comparison

Multi-factor designs:

python
design = "~batch + condition"  # Control for batch effects
design = "~age + condition"     # Include continuous covariate
design = "~group + condition + group:condition"  # Interaction effects

Design formula guidelines:

Use formulaic/Wilkinson formula notation (R-style)
Put adjustment variables (e.g., batch) before the main variable of interest
Ensure variables exist as columns in the metadata DataFrame
Use appropriate data types; continuous variables are detected from the formula, and categorical variables can be forced with C(variable) or a pandas categorical dtype
Do not use deprecated design_factors, continuous_factors, or ref_level arguments in new workflows

Step 3: DESeq2 Fitting

Initialize the DeseqDataSet and run the complete pipeline:

python
from pydeseq2.dds import DeseqDataSet
from pydeseq2.default_inference import DefaultInference

inference = DefaultInference(n_cpus=4)
dds = DeseqDataSet(
    counts=counts_df,
    metadata=metadata,
    design="~condition",
    refit_cooks=True,  # Refit after removing outliers
    inference=inference,
    low_memory=False,
)

# Run the complete DESeq2 pipeline
dds.deseq2()

What deseq2() does:

Computes size factors (normalization)
Fits genewise dispersions
Fits dispersion trend curve
Computes dispersion priors
Fits MAP dispersions (shrinkage)
Fits log fold changes
Calculates Cook's distances (outlier detection)
Refits if outliers detected (optional)

Step 4: Statistical Testing

Perform Wald tests to identify differentially expressed genes:

python
from pydeseq2.ds import DeseqStats

ds = DeseqStats(
    dds,
    contrast=["condition", "treated", "control"],  # Test treated vs control
    alpha=0.05,                # Significance threshold
    cooks_filter=True,         # Filter outliers
    independent_filter=True    # Filter low-power tests
)

ds.summary()

Contrast specification:

Format: [variable, test_level, reference_level]
Example: ["condition", "treated", "control"] tests treated vs control
Use a numeric contrast vector for continuous variables or complex coefficients
Default contrasts are no longer supported in PyDESeq2 0.5.x; always provide contrast

Result DataFrame columns:

baseMean: Mean normalized count across samples
log2FoldChange: Log2 fold change between conditions
lfcSE: Standard error of LFC
stat: Wald test statistic
pvalue: Raw p-value
padj: Adjusted p-value (FDR-corrected via Benjamini-Hochberg)

Step 5: Optional LFC Shrinkage

Apply shrinkage to reduce noise in fold change estimates:

python
ds.lfc_shrink(coeff="condition[T.treated]")  # Applies apeGLM shrinkage

When to use LFC shrinkage:

For visualization (volcano plots, heatmaps)
For ranking genes by effect size
When prioritizing genes for follow-up experiments

Important: Shrinkage affects only the log2FoldChange values, not the statistical test results (p-values remain unchanged). Use shrunk values for visualization but report unshrunken p-values for significance.

Step 6: Result Export

Save results and intermediate objects:

python
# Export results as CSV
ds.results_df.to_csv("deseq2_results.csv")

# Save significant genes only
significant = ds.results_df[ds.results_df.padj < 0.05]
significant.to_csv("significant_genes.csv")

# Save a portable AnnData object for later inspection
dds.to_picklable_anndata().write_h5ad("dds_result.h5ad")

Avoid loading pickle files from untrusted sources. For exchange between agents, pipelines, or collaborators, prefer CSV results and .h5ad AnnData files.

Common Analysis Patterns

Two-Group Comparison

Standard case-control comparison:

python
dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~condition")
dds.deseq2()

ds = DeseqStats(dds, contrast=["condition", "treated", "control"])
ds.summary()

results = ds.results_df
significant = results[results.padj < 0.05]

Multiple Comparisons

Testing multiple treatment groups against control:

python
dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~condition")
dds.deseq2()

treatments = ["treatment_A", "treatment_B", "treatment_C"]
all_results = {}

for treatment in treatments:
    ds = DeseqStats(dds, contrast=["condition", treatment, "control"])
    ds.summary()
    all_results[treatment] = ds.results_df

    sig_count = len(ds.results_df[ds.results_df.padj < 0.05])
    print(f"{treatment}: {sig_count} significant genes")

Accounting for Batch Effects

Control for technical variation:

python
# Include batch in design
dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~batch + condition")
dds.deseq2()

# Test condition while controlling for batch
ds = DeseqStats(dds, contrast=["condition", "treated", "control"])
ds.summary()

Continuous Covariates

Include continuous variables like age or dosage:

python
# Ensure continuous variable is numeric
metadata["age"] = pd.to_numeric(metadata["age"])

dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~age + condition")
dds.deseq2()

ds = DeseqStats(dds, contrast=["condition", "treated", "control"])
ds.summary()

Using the Analysis Script

This skill includes a complete command-line script for standard analyses:

bash
# Basic usage
python scripts/run_deseq2_analysis.py \
  --counts counts.csv \
  --metadata metadata.csv \
  --design "~condition" \
  --contrast condition treated control \
  --output results/

# With additional options
python scripts/run_deseq2_analysis.py \
  --counts counts.csv \
  --metadata metadata.csv \
  --design "~batch + condition" \
  --contrast condition treated control \
  --output results/ \
  --min-counts 10 \
  --alpha 0.05 \
  --n-cpus 4 \
  --shrink-coeff "condition[T.treated]" \
  --plots

Script features:

Automatic data loading and validation
Gene and sample filtering
Complete DESeq2 pipeline execution
Statistical testing with customizable parameters
Result export (CSV and portable AnnData/H5AD)
Explicit LFC shrinkage coefficient support for PyDESeq2 0.5.x
Optional visualization (volcano and MA plots)

Refer users to scripts/run_deseq2_analysis.py when they need a standalone analysis tool or want to batch process multiple datasets.

Result Interpretation

Identifying Significant Genes

python
# Filter by adjusted p-value
significant = ds.results_df[ds.results_df.padj < 0.05]

# Filter by both significance and effect size
sig_and_large = ds.results_df[
    (ds.results_df.padj < 0.05) &
    (abs(ds.results_df.log2FoldChange) > 1)
]

# Separate up- and down-regulated
upregulated = significant[significant.log2FoldChange > 0]
downregulated = significant[significant.log2FoldChange < 0]

print(f"Upregulated: {len(upregulated)}")
print(f"Downregulated: {len(downregulated)}")

Ranking and Sorting

python
# Sort by adjusted p-value
top_by_padj = ds.results_df.sort_values("padj").head(20)

# Sort by absolute fold change (use shrunk values)
ds.lfc_shrink(coeff="condition[T.treated]")
ds.results_df["abs_lfc"] = abs(ds.results_df.log2FoldChange)
top_by_lfc = ds.results_df.sort_values("abs_lfc", ascending=False).head(20)

# Sort by a combined metric
ds.results_df["score"] = -np.log10(ds.results_df.padj) * abs(ds.results_df.log2FoldChange)
top_combined = ds.results_df.sort_values("score", ascending=False).head(20)

Quality Metrics

python
# Check normalization (size factors should be close to 1)
print("Size factors:", dds.obs["size_factors"])

# Examine dispersion estimates
import matplotlib.pyplot as plt
plt.hist(dds.var["dispersions"], bins=50)
plt.xlabel("Dispersion")
plt.ylabel("Frequency")
plt.title("Dispersion Distribution")
plt.show()

# Check p-value distribution (should be mostly flat with peak near 0)
plt.hist(ds.results_df.pvalue.dropna(), bins=50)
plt.xlabel("P-value")
plt.ylabel("Frequency")
plt.title("P-value Distribution")
plt.show()

Visualization Guidelines

Volcano Plot

Visualize significance vs effect size:

python
import matplotlib.pyplot as plt
import numpy as np

results = ds.results_df.copy()
results["-log10(padj)"] = -np.log10(results.padj)

plt.figure(figsize=(10, 6))
significant = results.padj < 0.05

plt.scatter(
    results.loc[~significant, "log2FoldChange"],
    results.loc[~significant, "-log10(padj)"],
    alpha=0.3, s=10, c='gray', label='Not significant'
)
plt.scatter(
    results.loc[significant, "log2FoldChange"],
    results.loc[significant, "-log10(padj)"],
    alpha=0.6, s=10, c='red', label='padj < 0.05'
)

plt.axhline(-np.log10(0.05), color='blue', linestyle='--', alpha=0.5)
plt.xlabel("Log2 Fold Change")
plt.ylabel("-Log10(Adjusted P-value)")
plt.title("Volcano Plot")
plt.legend()
plt.savefig("volcano_plot.png", dpi=300)

MA Plot

Show fold change vs mean expression:

python
plt.figure(figsize=(10, 6))

plt.scatter(
    np.log10(results.loc[~significant, "baseMean"] + 1),
    results.loc[~significant, "log2FoldChange"],
    alpha=0.3, s=10, c='gray'
)
plt.scatter(
    np.log10(results.loc[significant, "baseMean"] + 1),
    results.loc[significant, "log2FoldChange"],
    alpha=0.6, s=10, c='red'
)

plt.axhline(0, color='blue', linestyle='--', alpha=0.5)
plt.xlabel("Log10(Base Mean + 1)")
plt.ylabel("Log2 Fold Change")
plt.title("MA Plot")
plt.savefig("ma_plot.png", dpi=300)

Troubleshooting Common Issues

Data Format Problems

Issue: "Index mismatch between counts and metadata"

Solution: Ensure sample names match exactly

python
print("Counts samples:", counts_df.index.tolist())
print("Metadata samples:", metadata.index.tolist())

# Take intersection if needed
common = counts_df.index.intersection(metadata.index)
counts_df = counts_df.loc[common]
metadata = metadata.loc[common]

Issue: "All genes have zero counts"

Solution: Check if data needs transposition

python
print(f"Counts shape: {counts_df.shape}")
# If genes > samples, transpose is needed
if counts_df.shape[1] < counts_df.shape[0]:
    counts_df = counts_df.T

Design Matrix Issues

Issue: "Design matrix is not full rank"

Cause: Confounded variables (e.g., all treated samples in one batch)

Solution: Remove confounded variable or add interaction term

python
# Check confounding
print(pd.crosstab(metadata.condition, metadata.batch))

# Either simplify design or add interaction
design = "~condition"  # Remove batch
# OR
design = "~condition + batch + condition:batch"  # Model interaction

No Significant Genes

Diagnostics:

python
# Check dispersion distribution
plt.hist(dds.var["dispersions"], bins=50)
plt.show()

# Check size factors
print(dds.obs["size_factors"])

# Look at top genes by raw p-value
print(ds.results_df.nsmallest(20, "pvalue"))

Possible causes:

Small effect sizes
High biological variability
Insufficient sample size
Technical issues (batch effects, outliers)

Reference Documentation

For comprehensive details beyond this workflow-oriented guide:

API Reference (references/api_reference.md): Complete documentation of PyDESeq2 classes, methods, and data structures. Use when needing detailed parameter information or understanding object attributes.
Workflow Guide (references/workflow_guide.md): In-depth guide covering complete analysis workflows, data loading patterns, multi-factor designs, troubleshooting, and best practices. Use when handling complex experimental designs or encountering issues.

Load these references into context when users need:

Detailed API documentation: Read references/api_reference.md
Comprehensive workflow examples: Read references/workflow_guide.md
Troubleshooting guidance: Read references/workflow_guide.md (see Troubleshooting section)

Key Reminders

Data orientation matters: Count matrices typically load as genes × samples but need to be samples × genes. Always transpose with .T if needed.
Sample filtering: Remove samples with missing metadata before analysis to avoid errors.
Gene filtering: Filter low-count genes (e.g., < 10 total reads) to improve power and reduce computational time.
Design formula order: Put adjustment variables before the variable of interest (e.g., "~batch + condition" not "~condition + batch").
LFC shrinkage timing: Apply shrinkage after statistical testing and only for visualization/ranking purposes. P-values remain based on unshrunken estimates.
Result interpretation: Use padj < 0.05 for significance, not raw p-values. The Benjamini-Hochberg procedure controls false discovery rate.
Contrast specification: The format is [variable, test_level, reference_level] where test_level is compared against reference_level.
Save intermediate objects: Prefer dds.to_picklable_anndata().write_h5ad("dds_result.h5ad") for portable outputs. Only load pickle files that you created yourself and trust.

Installation and Requirements

bash
uv pip install pydeseq2==0.5.4

System requirements:

Python 3.11+
PyDESeq2 0.5.4
pandas 2.2.0+
numpy 2.0.0+
scipy 1.12.0+
scikit-learn 1.4.0+
anndata 0.11.0+
formulaic 1.0.2+ and formulaic-contrasts 0.2.0+

Optional for visualization:

matplotlib
seaborn

Additional Resources

Official Documentation: https://pydeseq2.readthedocs.io
GitHub Repository: https://github.com/scverse/PyDESeq2
Publication: Muzellec et al. (2023) Bioinformatics, DOI: 10.1093/bioinformatics/btad547
Original DESeq2 (R): Love et al. (2014) Genome Biology, DOI: 10.1186/s13059-014-0550-8

Publisher	dralkh
Repository	`seerai`
Skill name	`pydeseq2`
Stars	96
Forks	5
Bundled files	3
License	MIT
Links	Skill source Repository

Overview

Installation

Install in TypingMind

Install in another agent

Use it in TypingMind

Works with any AI model

Skill instructions

PyDESeq2

Overview

When to Use This Skill

Quick Start Workflow

Core Workflow Steps

Step 1: Data Preparation

Step 2: Design Specification

Step 3: DESeq2 Fitting

Step 4: Statistical Testing

Step 5: Optional LFC Shrinkage

Step 6: Result Export

Common Analysis Patterns

Two-Group Comparison

Multiple Comparisons

Accounting for Batch Effects

Continuous Covariates

Using the Analysis Script

Result Interpretation

Identifying Significant Genes

Ranking and Sorting

Quality Metrics

Visualization Guidelines

Volcano Plot

MA Plot

Troubleshooting Common Issues

Data Format Problems

Design Matrix Issues

No Significant Genes

Reference Documentation

Key Reminders

Installation and Requirements

Additional Resources

Bundled files

Frequently asked questions

Related AI Skills

Set up your own AI workspace now