Skills › Scientific Agent Skills › pydeseq2 › run_deseq2_analysis.py
run_deseq2_analysis.py
466 lines · skills/pydeseq2/scripts/run_deseq2_analysis.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3PyDESeq2 Analysis Script
5This script performs a complete differential expression analysis using PyDESeq2.
6It can be used as a template for standard RNA-seq DEA workflows.
8Usage:
9 python run_deseq2_analysis.py --counts counts.csv --metadata metadata.csv \
10 --design "~condition" --contrast condition treated control \
11 --output results/
13Requirements:
14 - pydeseq2
15 - pandas
16 - matplotlib (installed by PyDESeq2; used for plots)
17"""
19import argparse
20import csv
21import json
22from importlib.metadata import version
23import sys
24from pathlib import Path
26import numpy as np
27import pandas as pd
29try:
30 from pydeseq2.dds import DeseqDataSet
31 from pydeseq2.default_inference import DefaultInference
32 from pydeseq2.ds import DeseqStats
33except ImportError:
34 print("Error: pydeseq2 not installed. Install with: uv pip install pydeseq2==0.5.4")
35 sys.exit(1)
38def read_indexed_csv(path):
39 """Reject ambiguous source headers before pandas can rename duplicates."""
40 with open(path, newline="", encoding="utf-8-sig") as handle:
41 reader = csv.reader(handle)
42 header = next(reader, [])
43 for line_number, row in enumerate(reader, start=2):
44 if len(row) != len(header):
45 raise ValueError(f"CSV row {line_number} has the wrong number of fields")
46 if len(header) < 2 or any(not item.strip() for item in header[1:]):
47 raise ValueError("CSV requires an identifier column and named data columns")
48 if len(set(header)) != len(header):
49 raise ValueError("CSV contains duplicate headers")
50 frame = pd.read_csv(path, dtype={0: str})
51 frame = frame.set_index(frame.columns[0])
52 for axis in (frame.index, frame.columns):
53 if axis.has_duplicates or axis.isna().any() or any(not str(x).strip() for x in axis):
54 raise ValueError("Sample, gene, and metadata identifiers must be nonempty and unique")
55 return frame
58def validate_counts(counts_df):
59 """Validate the count domain; integer-valued data do not prove raw-count provenance."""
60 if counts_df.empty:
61 raise ValueError("Count matrix must contain samples and genes")
62 try:
63 values = counts_df.to_numpy(dtype=float)
64 except (TypeError, ValueError) as exc:
65 raise ValueError("Counts must be numeric") from exc
66 if not np.isfinite(values).all():
67 raise ValueError("Counts must be finite without missing values")
68 if (values < 0).any():
69 raise ValueError("Count matrix contains negative values")
70 if (values != np.floor(values)).any() or (values > 2**53 - 1).any():
71 raise ValueError("Counts must be exact nonnegative integers no larger than 2**53 - 1")
72 if (values.sum(axis=1) == 0).any():
73 raise ValueError("Every sample must have positive total counts")
76def load_and_validate_data(counts_path, metadata_path, transpose_counts=True):
77 """Load explicitly oriented CSVs, require the same sample set, then align."""
78 counts_df = read_indexed_csv(counts_path)
79 if transpose_counts:
80 counts_df = counts_df.T
81 metadata = read_indexed_csv(metadata_path)
82 if set(counts_df.index) != set(metadata.index):
83 raise ValueError("Sample sets differ; resolve omissions explicitly before analysis")
84 metadata = metadata.loc[counts_df.index].copy()
85 validate_counts(counts_df)
86 print(f"[OK] Loaded {counts_df.shape[0]} samples x {counts_df.shape[1]} genes")
87 return counts_df.astype(np.int64), metadata
90def filter_data(counts_df, metadata, min_counts=10, condition_col=None):
91 """Apply a declared total-count gene prefilter without silently excluding samples."""
92 if isinstance(min_counts, bool) or not isinstance(min_counts, (int, np.integer)) or min_counts < 0:
93 raise ValueError("min_counts must be a nonnegative integer")
94 if not counts_df.index.equals(metadata.index):
95 raise ValueError("Counts and metadata sample order must match")
96 if condition_col:
97 if condition_col not in metadata:
98 raise ValueError(f"Missing contrast column: {condition_col}")
99 if metadata[condition_col].isna().any():
100 raise ValueError("Missing contrast annotations; resolve sample exclusions explicitly")
101 validate_counts(counts_df)
102 counts_df = counts_df.loc[:, counts_df.sum(axis=0) >= min_counts].copy()
103 validate_counts(counts_df)
104 print(f"[OK] Retained {counts_df.shape[1]} genes at total counts >= {min_counts}")
105 return counts_df, metadata
108def prepare_contrast_metadata(metadata, contrast):
109 """Set the CLI categorical reference before fitting, preserving all observed levels."""
110 variable, tested, reference = contrast
111 if variable not in metadata or metadata[variable].isna().any():
112 raise ValueError("Contrast variable must exist and have no missing annotations")
113 values = metadata[variable].astype(str)
114 levels = list(pd.unique(values))
115 if tested == reference or tested not in levels or reference not in levels:
116 raise ValueError("Contrast requires two distinct observed levels")
117 metadata = metadata.copy()
118 metadata[variable] = pd.Categorical(values, categories=[reference] + [x for x in levels if x != reference])
119 return metadata
122def run_deseq2(counts_df, metadata, design, n_cpus=1):
123 """Run DESeq2 normalization and fitting."""
124 print(f"\nInitializing DeseqDataSet with design: {design}")
126 validate_counts(counts_df)
127 if n_cpus < 1:
128 raise ValueError("n_cpus must be positive")
129 inference = DefaultInference(n_cpus=n_cpus)
130 dds = DeseqDataSet(
131 counts=counts_df,
132 metadata=metadata,
133 design=design,
134 refit_cooks=True,
135 inference=inference,
136 quiet=False
137 )
139 matrix = dds.obsm["design_matrix"]
140 if not matrix.index.equals(counts_df.index) or not np.isfinite(matrix.to_numpy()).all():
141 raise ValueError("Design matrix must retain every sample and contain finite values")
142 if np.linalg.matrix_rank(matrix.to_numpy()) != matrix.shape[1]:
143 raise ValueError("Design is not full rank; inspect confounding and redundant terms")
144 if matrix.shape[0] <= matrix.shape[1]:
145 raise ValueError("Design requires residual degrees of freedom and biological replication")
147 print("\nRunning DESeq2 pipeline...")
148 print(" Step 1/7: Computing size factors...")
149 print(" Step 2/7: Fitting genewise dispersions...")
150 print(" Step 3/7: Fitting dispersion trend curve...")
151 print(" Step 4/7: Computing dispersion priors...")
152 print(" Step 5/7: Fitting MAP dispersions...")
153 print(" Step 6/7: Fitting log fold changes...")
154 print(" Step 7/7: Calculating Cook's distances...")
156 dds.deseq2()
158 print("\n[OK] DESeq2 fitting complete")
160 return dds, inference
163def infer_shrink_coeff(dds, contrast, coeff=None):
164 """Require the tested contrast to equal one positive design coefficient."""
165 columns = list(dds.obsm["design_matrix"].columns)
166 vector = np.asarray(dds.contrast(column=contrast[0], group_to_compare=contrast[1], baseline=contrast[2]))
167 matches = [name for i, name in enumerate(columns)
168 if np.allclose(vector, np.eye(len(columns))[i], rtol=0, atol=1e-12)]
169 if len(matches) != 1 or (coeff is not None and coeff != matches[0]):
170 raise ValueError(
171 "LFC shrinkage needs a single positive coefficient matching the exact contrast. "
172 f"Contrast vector: {vector.tolist()}; columns: {columns}. "
173 "Relevel and refit, or use --no-shrink; --shrink-coeff cannot override a mismatch."
174 )
175 return matches[0]
178def run_statistical_tests(dds, contrast, alpha=0.05, shrink_lfc=True, inference=None, shrink_coeff=None):
179 """Perform Wald tests and compute p-values."""
180 if not np.isfinite(alpha) or not 0 < alpha < 1:
181 raise ValueError("alpha must be finite and between 0 and 1")
182 coeff = infer_shrink_coeff(dds, contrast, shrink_coeff) if shrink_lfc else None
183 print(f"\nPerforming statistical tests...")
184 print(f" Contrast: {contrast}")
185 print(f" Significance threshold: {alpha}")
187 ds = DeseqStats(
188 dds,
189 contrast=contrast,
190 alpha=alpha,
191 cooks_filter=True,
192 independent_filter=True,
193 inference=inference,
194 quiet=False
195 )
197 print("\n Running Wald tests...")
198 print(" Filtering outliers based on Cook's distance...")
199 print(" Applying independent filtering...")
200 print(" Adjusting p-values (Benjamini-Hochberg)...")
202 ds.summary()
203 ds.unshrunk_results_df = ds.results_df.copy(deep=True)
204 ds.shrink_coeff = coeff
206 print("\n[OK] Statistical testing complete")
208 # Optional LFC shrinkage
209 if shrink_lfc:
210 print("\nApplying LFC shrinkage for visualization...")
211 ds.lfc_shrink(coeff=coeff)
212 print("[OK] LFC shrinkage complete")
214 return ds
217def save_results(ds, dds, output_dir):
218 """Save results and intermediate objects."""
219 output_dir = Path(output_dir)
220 output_dir.mkdir(parents=True, exist_ok=True)
222 print(f"\nSaving results to {output_dir}/")
224 # Save statistical results
225 results_path = output_dir / "deseq2_results.csv"
226 ds.results_df.to_csv(results_path)
227 unshrunk = getattr(ds, "unshrunk_results_df", ds.results_df)
228 unshrunk.to_csv(output_dir / "deseq2_results_unshrunken.csv")
229 print(f" Saved: {results_path}")
231 # Save significant genes
232 alpha = ds.alpha
233 significant = ds.results_df[ds.results_df.padj < alpha]
234 sig_path = output_dir / "significant_genes.csv"
235 significant.to_csv(sig_path)
236 print(f" Saved: {sig_path} ({len(significant)} significant genes)")
238 # Save sorted results
239 sorted_results = ds.results_df.sort_values("padj")
240 sorted_path = output_dir / "results_sorted_by_padj.csv"
241 sorted_results.to_csv(sorted_path)
242 print(f" Saved: {sorted_path}")
244 # Save as AnnData/H5AD to avoid unsafe pickle interchange.
245 dds_path = output_dir / "deseq_dataset.h5ad"
246 dds.to_picklable_anndata().write_h5ad(dds_path)
247 print(f" Saved: {dds_path}")
249 manifest = {
250 "alpha": alpha,
251 "design": getattr(dds, "design", None),
252 "size_factors_fit_type_requested": getattr(dds, "size_factors_fit_type", None),
253 "retained_genes": len(ds.results_df),
254 "contrast": list(ds.contrast),
255 "shrink_coefficient": getattr(ds, "shrink_coeff", None),
256 "pvalues_and_stat": "unshrunken Wald test",
257 "versions": {name: version(name) for name in ("pydeseq2", "anndata", "numpy", "pandas")},
258 }
259 (output_dir / "analysis_manifest.json").write_text(json.dumps(manifest, indent=2) + "\n")
261 # Print summary
262 print(f"\n{'='*60}")
263 print("ANALYSIS SUMMARY")
264 print(f"{'='*60}")
265 print(f"Genes retained: {len(ds.results_df)}; finite p-values: {ds.results_df.pvalue.notna().sum()}")
266 print(f"Significant genes (padj < {alpha}): {len(significant)}")
267 print(f"Upregulated: {len(significant[significant.log2FoldChange > 0])}")
268 print(f"Downregulated: {len(significant[significant.log2FoldChange < 0])}")
269 print(f"{'='*60}")
271 # Show top genes
272 print("\nTop 10 most significant genes:")
273 print(sorted_results.head(10)[["baseMean", "log2FoldChange", "pvalue", "padj"]])
275 return results_path
278def create_plots(ds, output_dir):
279 """Create basic visualization plots."""
280 try:
281 import matplotlib.pyplot as plt
282 import numpy as np
283 except ImportError:
284 print("\nNote: matplotlib not installed. Skipping plot generation.")
285 return
287 output_dir = Path(output_dir)
288 results = ds.results_df.copy()
289 alpha = ds.alpha
290 output_dir.mkdir(parents=True, exist_ok=True)
292 print("\nGenerating plots...")
294 # Volcano plot
295 # Missing adjusted p-values remain missing; zero p-values are capped only for display.
296 results["-log10(padj)"] = -np.log10(results.padj.clip(lower=np.finfo(float).tiny))
298 plt.figure(figsize=(10, 6))
299 significant = results.padj < alpha
300 plt.scatter(
301 results.loc[~significant, "log2FoldChange"],
302 results.loc[~significant, "-log10(padj)"],
303 alpha=0.3, s=10, c='gray', label='Not significant'
304 )
305 plt.scatter(
306 results.loc[significant, "log2FoldChange"],
307 results.loc[significant, "-log10(padj)"],
308 alpha=0.6, s=10, c='red', label=f'Significant (padj < {alpha})'
309 )
310 plt.axhline(-np.log10(alpha), color='blue', linestyle='--', linewidth=1, alpha=0.5)
311 plt.axvline(1, color='gray', linestyle='--', linewidth=1, alpha=0.5)
312 plt.axvline(-1, color='gray', linestyle='--', linewidth=1, alpha=0.5)
313 plt.xlabel("Log2 Fold Change", fontsize=12)
314 plt.ylabel("-Log10(Adjusted P-value)", fontsize=12)
315 plt.title("Volcano Plot", fontsize=14, fontweight='bold')
316 plt.legend()
317 plt.tight_layout()
318 volcano_path = output_dir / "volcano_plot.png"
319 plt.savefig(volcano_path, dpi=300)
320 plt.close()
321 print(f" Saved: {volcano_path}")
323 # MA plot: distinguish absent adjusted tests from tested nonsignificant genes.
324 plt.figure(figsize=(10, 6))
325 tested_nonsignificant = (~significant) & results.padj.notna()
326 plt.scatter(
327 np.log10(results.loc[tested_nonsignificant, "baseMean"] + 1),
328 results.loc[tested_nonsignificant, "log2FoldChange"],
329 alpha=0.3, s=10, c='gray', label='Not significant'
330 )
331 plt.scatter(
332 np.log10(results.loc[significant, "baseMean"] + 1),
333 results.loc[significant, "log2FoldChange"],
334 alpha=0.6, s=10, c='red', label=f'Significant (padj < {alpha})'
335 )
336 missing = results.padj.isna()
337 plt.scatter(
338 np.log10(results.loc[missing, "baseMean"] + 1),
339 results.loc[missing, "log2FoldChange"],
340 alpha=0.4, s=12, c='silver', marker='x', label='No adjusted p-value'
341 )
342 plt.axhline(0, color='blue', linestyle='--', linewidth=1, alpha=0.5)
343 plt.xlabel("Log10(Base Mean + 1)", fontsize=12)
344 plt.ylabel("Log2 Fold Change", fontsize=12)
345 plt.title("MA Plot", fontsize=14, fontweight='bold')
346 plt.legend()
347 plt.tight_layout()
348 ma_path = output_dir / "ma_plot.png"
349 plt.savefig(ma_path, dpi=300)
350 plt.close()
351 print(f" Saved: {ma_path}")
354def main():
355 parser = argparse.ArgumentParser(
356 description="Run PyDESeq2 differential expression analysis",
357 formatter_class=argparse.RawDescriptionHelpFormatter,
358 epilog="""
359Examples:
360 # Basic analysis
361 python run_deseq2_analysis.py \\
362 --counts counts.csv \\
363 --metadata metadata.csv \\
364 --design "~condition" \\
365 --contrast condition treated control \\
366 --output results/
368 # Multi-factor analysis
369 python run_deseq2_analysis.py \\
370 --counts counts.csv \\
371 --metadata metadata.csv \\
372 --design "~batch + condition" \\
373 --contrast condition treated control \\
374 --output results/ \\
375 --n-cpus 4 \\
376 --shrink-coeff "condition[T.treated]"
377 """
378 )
380 parser.add_argument("--counts", required=True, help="Path to count matrix CSV file")
381 parser.add_argument("--metadata", required=True, help="Path to metadata CSV file")
382 parser.add_argument("--design", required=True, help="Design formula (e.g., '~condition')")
383 parser.add_argument(
384 "--contrast",
385 nargs=3,
386 required=True,
387 metavar=("VARIABLE", "TEST", "REFERENCE"),
388 help="Contrast specification: variable test_level reference_level",
389 )
390 parser.add_argument("--output", default="results", help="Output directory (default: results)")
391 parser.add_argument(
392 "--min-counts",
393 type=int,
394 default=10,
395 help="Minimum total counts for gene filtering (default: 10)",
396 )
397 parser.add_argument("--alpha", type=float, default=0.05, help="Significance threshold (default: 0.05)")
398 parser.add_argument(
399 "--no-transpose",
400 action="store_true",
401 help="Don't transpose count matrix (use if already samples × genes)",
402 )
403 parser.add_argument("--no-shrink", action="store_true", help="Skip LFC shrinkage")
404 parser.add_argument(
405 "--shrink-coeff",
406 help="Design-matrix coefficient to shrink (e.g., 'condition[T.treated]')",
407 )
408 parser.add_argument(
409 "--n-cpus",
410 type=int,
411 default=1,
412 help="Number of CPUs for parallel processing (default: 1)",
413 )
414 parser.add_argument("--plots", action="store_true", help="Generate volcano and MA plots")
416 args = parser.parse_args()
417 if not np.isfinite(args.alpha) or not 0 < args.alpha < 1:
418 parser.error("--alpha must be finite and between 0 and 1")
419 if args.n_cpus < 1 or args.min_counts < 0:
420 parser.error("--n-cpus must be positive and --min-counts nonnegative")
421 if args.no_shrink and args.shrink_coeff:
422 parser.error("--shrink-coeff cannot be combined with --no-shrink")
424 # Load data
425 counts_df, metadata = load_and_validate_data(
426 args.counts,
427 args.metadata,
428 transpose_counts=not args.no_transpose,
429 )
431 # Filter data
432 condition_col = args.contrast[0]
433 counts_df, metadata = filter_data(
434 counts_df,
435 metadata,
436 min_counts=args.min_counts,
437 condition_col=condition_col,
438 )
440 metadata = prepare_contrast_metadata(metadata, args.contrast)
442 # Run DESeq2
443 dds, inference = run_deseq2(counts_df, metadata, args.design, n_cpus=args.n_cpus)
445 # Statistical testing
446 ds = run_statistical_tests(
447 dds,
448 contrast=args.contrast,
449 alpha=args.alpha,
450 shrink_lfc=not args.no_shrink,
451 inference=inference,
452 shrink_coeff=args.shrink_coeff,
453 )
455 # Save results
456 save_results(ds, dds, args.output)
458 # Create plots if requested
459 if args.plots:
460 create_plots(ds, args.output)
462 print(f"\n[OK] Analysis complete! Results saved to {args.output}/")
465if __name__ == "__main__":
466 main()