Skills › Scientific Agent Skills › scanpy › analysis_template.py
analysis_template.py
298 lines · skills/scanpy/assets/analysis_template.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Complete Single-Cell Analysis Template
5This template provides a complete workflow for single-cell RNA-seq analysis
6using scanpy, from data loading through clustering and cell type annotation.
8Illustrative template targeting Scanpy 1.12.4. Customize thresholds, metadata,
9and annotation; this file is not an automatically validated biological analysis.
10"""
12import scanpy as sc
13import pandas as pd
14import numpy as np
15import matplotlib.pyplot as plt
17# ============================================================================
18# CONFIGURATION
19# ============================================================================
21# File paths
22INPUT_FILE = 'data/raw_counts.h5ad' # Change to your input file
23OUTPUT_DIR = 'results/'
24FIGURES_DIR = 'figures/'
26# QC parameters
27MIN_GENES = 200 # Minimum genes per cell
28MIN_CELLS = 3 # Minimum cells per gene
29MT_THRESHOLD = 5 # Maximum mitochondrial percentage
31# Analysis parameters
32N_TOP_GENES = 2000 # Number of highly variable genes
33N_PCS = 40 # Number of principal components
34N_NEIGHBORS = 10 # Number of neighbors for graph
35LEIDEN_RESOLUTION = 0.5 # Clustering resolution
37# Scanpy settings
38sc.settings.verbosity = 3
39sc.set_figure_params(dpi=80, facecolor='white')
40sc.settings.figdir = FIGURES_DIR
41sc.settings.autosave = True
43# ============================================================================
44# 1. LOAD DATA
45# ============================================================================
47print("=" * 80)
48print("LOADING DATA")
49print("=" * 80)
51# Load data (adjust based on your file format)
52adata = sc.read_h5ad(INPUT_FILE)
53# adata = sc.read_10x_mtx('data/filtered_gene_bc_matrices/') # For 10X data
54# adata = sc.read_csv('data/counts.csv') # For CSV data
56print(f"Loaded: {adata.n_obs} cells x {adata.n_vars} genes")
58# ============================================================================
59# 2. QUALITY CONTROL
60# ============================================================================
62print("\n" + "=" * 80)
63print("QUALITY CONTROL")
64print("=" * 80)
66# Identify mitochondrial genes
67adata.var['mt'] = adata.var_names.str.startswith('MT-')
69# Calculate QC metrics
70sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None,
71 log1p=False, inplace=True)
73# Visualize QC metrics before filtering
74sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],
75 jitter=0.4, multi_panel=True, save='_qc_before_filtering')
77sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt', save='_qc_mt')
78sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts', save='_qc_genes')
80# Filter cells and genes
81print(f"\nBefore filtering: {adata.n_obs} cells, {adata.n_vars} genes")
83sc.pp.filter_cells(adata, min_genes=MIN_GENES)
84sc.pp.filter_genes(adata, min_cells=MIN_CELLS)
85adata = adata[adata.obs.pct_counts_mt < MT_THRESHOLD, :].copy()
87print(f"After filtering: {adata.n_obs} cells, {adata.n_vars} genes")
89# Optional: doublet detection (uncomment; run before normalization)
90# sc.pp.scrublet(adata)
91# adata = adata[~adata.obs['predicted_doublet'], :].copy()
92# print(f"After doublet removal: {adata.n_obs} cells")
94# ============================================================================
95# 3. NORMALIZATION
96# ============================================================================
98print("\n" + "=" * 80)
99print("NORMALIZATION")
100print("=" * 80)
102# Identified raw counts: preserve before any normalization
103adata.layers["counts"] = adata.X.copy()
104# Normalize to 10,000 counts per cell
105sc.pp.normalize_total(adata, target_sum=1e4)
107# Log-transform
108sc.pp.log1p(adata)
110# Store normalized data
111adata.raw = adata.copy()
113# ============================================================================
114# 4. FEATURE SELECTION
115# ============================================================================
117print("\n" + "=" * 80)
118print("FEATURE SELECTION")
119print("=" * 80)
121# Identify highly variable genes
122sc.pp.highly_variable_genes(adata, n_top_genes=N_TOP_GENES)
124# Visualize
125sc.pl.highly_variable_genes(adata, save='_hvg')
127print(f"Selected {sum(adata.var.highly_variable)} highly variable genes")
129# Keep all genes/counts in the exported object; PCA automatically uses the HVG mask.
131# ============================================================================
132# 5. SCALING AND REGRESSION
133# ============================================================================
135print("\n" + "=" * 80)
136print("SCALING AND REGRESSION")
137print("=" * 80)
139# Regress out unwanted sources of variation
140# sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])
142# Scale data
143# sc.pp.scale(adata, max_value=10) # Optional; may densify sparse data
145# ============================================================================
146# 6. DIMENSIONALITY REDUCTION
147# ============================================================================
149print("\n" + "=" * 80)
150print("DIMENSIONALITY REDUCTION")
151print("=" * 80)
153# PCA
154n_pcs = min(N_PCS, adata.n_obs - 1, int(adata.var.highly_variable.sum()) - 1)
155sc.pp.pca(adata, n_comps=n_pcs, svd_solver='arpack')
156sc.pl.pca_variance_ratio(adata, log=True, save='_pca_variance')
158# Compute neighborhood graph
159sc.pp.neighbors(adata, n_neighbors=N_NEIGHBORS, n_pcs=n_pcs, use_rep="X_pca")
161# UMAP
162sc.tl.umap(adata)
164# ============================================================================
165# 7. CLUSTERING
166# ============================================================================
168print("\n" + "=" * 80)
169print("CLUSTERING")
170print("=" * 80)
172# Leiden clustering
173sc.tl.leiden(adata, resolution=LEIDEN_RESOLUTION, flavor="igraph", directed=False, n_iterations=2)
175# Visualize
176sc.pl.umap(adata, color='leiden', legend_loc='on data', save='_leiden')
178print(f"Identified {len(adata.obs['leiden'].unique())} clusters")
180# ============================================================================
181# 8. MARKER GENE IDENTIFICATION
182# ============================================================================
184print("\n" + "=" * 80)
185print("MARKER GENE IDENTIFICATION")
186print("=" * 80)
188# Find marker genes (exploratory — pseudobulk + pydeseq2 for rigorous DE)
189sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon', use_raw=True,
190 reference='rest', corr_method='benjamini-hochberg', tie_correct=True)
192# Visualize top markers
193sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False, save='_markers')
194sc.pl.rank_genes_groups_heatmap(adata, n_genes=10, save='_markers_heatmap')
195sc.pl.rank_genes_groups_dotplot(adata, n_genes=5, save='_markers_dotplot')
197# Get top markers for each cluster
198for cluster in adata.obs['leiden'].unique():
199 print(f"\nCluster {cluster} top markers:")
200 markers = sc.get.rank_genes_groups_df(adata, group=cluster).head(10)
201 print(markers[['names', 'scores', 'pvals_adj']].to_string(index=False))
203# ============================================================================
204# 9. CELL TYPE ANNOTATION (CUSTOMIZE THIS SECTION)
205# ============================================================================
207print("\n" + "=" * 80)
208print("CELL TYPE ANNOTATION")
209print("=" * 80)
211# Example marker genes for common cell types (customize for your data)
212marker_genes = {
213 'T cells': ['CD3D', 'CD3E', 'CD3G'],
214 'B cells': ['MS4A1', 'CD79A', 'CD79B'],
215 'Monocytes': ['CD14', 'LYZ', 'S100A8'],
216 'NK cells': ['NKG7', 'GNLY', 'KLRD1'],
217 'Dendritic cells': ['FCER1A', 'CST3'],
218}
220# Visualize marker genes
221for cell_type, genes in marker_genes.items():
222 available_genes = [g for g in genes if g in adata.raw.var_names]
223 if available_genes:
224 sc.pl.umap(adata, color=available_genes, use_raw=True,
225 save=f'_{cell_type.replace(" ", "_")}')
227# Manual annotation based on marker expression (customize this mapping)
228cluster_to_celltype = {} # Fill only after reviewing markers, provenance and uncertainty.
230# Apply annotations
231adata.obs['cell_type'] = adata.obs['leiden'].astype(str).map(cluster_to_celltype)
232adata.obs['cell_type'] = adata.obs['cell_type'].fillna('Unknown').astype('category')
234# Visualize annotated cell types
235sc.pl.umap(adata, color='cell_type', legend_loc='on data', save='_celltypes')
237# ============================================================================
238# 10. ADDITIONAL ANALYSES (OPTIONAL)
239# ============================================================================
241print("\n" + "=" * 80)
242print("ADDITIONAL ANALYSES")
243print("=" * 80)
245# PAGA trajectory analysis (optional)
246# sc.tl.paga(adata, groups='leiden')
247# sc.pl.paga(adata, color='leiden', save='_paga')
249# Gene set scoring (optional)
250# example_gene_set = ['CD3D', 'CD3E', 'CD3G']
251# sc.tl.score_genes(adata, example_gene_set, score_name='T_cell_score')
252# sc.pl.umap(adata, color='T_cell_score', save='_gene_set_score')
254# ============================================================================
255# 11. SAVE RESULTS
256# ============================================================================
258print("\n" + "=" * 80)
259print("SAVING RESULTS")
260print("=" * 80)
262import os
263os.makedirs(OUTPUT_DIR, exist_ok=True)
265# Save processed AnnData object
266adata.write(f'{OUTPUT_DIR}/processed_data.h5ad')
267print(f"Saved processed data to {OUTPUT_DIR}/processed_data.h5ad")
269# Export metadata
270adata.obs.to_csv(f'{OUTPUT_DIR}/cell_metadata.csv')
271adata.var.to_csv(f'{OUTPUT_DIR}/gene_metadata.csv')
272print(f"Saved metadata to {OUTPUT_DIR}/")
274# Export marker genes
275for cluster in adata.obs['leiden'].unique():
276 markers = sc.get.rank_genes_groups_df(adata, group=cluster)
277 markers.to_csv(f'{OUTPUT_DIR}/markers_cluster_{cluster}.csv', index=False)
278print(f"Saved marker genes to {OUTPUT_DIR}/")
280# ============================================================================
281# 12. SUMMARY
282# ============================================================================
284print("\n" + "=" * 80)
285print("ANALYSIS SUMMARY")
286print("=" * 80)
288print(f"\nFinal dataset:")
289print(f" Cells: {adata.n_obs}")
290print(f" Genes: {adata.n_vars}")
291print(f" Clusters: {len(adata.obs['leiden'].unique())}")
293print(f"\nCell type distribution:")
294print(adata.obs['cell_type'].value_counts())
296print("\n" + "=" * 80)
297print("ANALYSIS COMPLETE")
298print("=" * 80)