Skills › Knowledge Work Plugins › single-cell-rna-qc › qc_analysis.py
qc_analysis.py
232 lines · bio-research/skills/single-cell-rna-qc/scripts/qc_analysis.py at 2574e29
1#!/usr/bin/env python3
2"""
3Quality Control Analysis for Single-Cell RNA-seq Data
4Following scverse best practices from:
5https://www.sc-best-practices.org/preprocessing_visualization/quality_control.html
7This is a convenience script that runs a complete QC workflow using the
8modular functions from qc_core.py and qc_plotting.py.
9"""
11import anndata as ad
12import scanpy as sc
13import sys
14import os
15import argparse
17# Import our modular utilities
18from qc_core import (
19 calculate_qc_metrics,
20 detect_outliers_mad,
21 apply_hard_threshold,
22 filter_cells,
23 filter_genes,
24 print_qc_summary
25)
26from qc_plotting import (
27 plot_qc_distributions,
28 plot_filtering_thresholds,
29 plot_qc_after_filtering
30)
32print("=" * 80)
33print("Single-Cell RNA-seq Quality Control Analysis")
34print("=" * 80)
36# Default parameters (single source of truth)
37DEFAULT_MAD_COUNTS = 5
38DEFAULT_MAD_GENES = 5
39DEFAULT_MAD_MT = 3
40DEFAULT_MT_THRESHOLD = 8
41DEFAULT_MIN_CELLS = 20
42DEFAULT_MT_PATTERN = 'mt-,MT-'
43DEFAULT_RIBO_PATTERN = 'Rpl,Rps,RPL,RPS'
44DEFAULT_HB_PATTERN = '^Hb[^(p)]|^HB[^(P)]'
46# Parse command-line arguments
47parser = argparse.ArgumentParser(
48 description='Quality Control Analysis for Single-Cell RNA-seq Data',
49 formatter_class=argparse.RawDescriptionHelpFormatter,
50 epilog="""
51Examples:
52 python3 qc_analysis.py data.h5ad
53 python3 qc_analysis.py raw_feature_bc_matrix.h5
54 python3 qc_analysis.py data.h5ad --mad-counts 4 --mad-genes 4 --mad-mt 2.5
55 python3 qc_analysis.py data.h5ad --mt-threshold 10 --min-cells 10
56 python3 qc_analysis.py data.h5ad --mt-pattern "^mt-" --ribo-pattern "^Rpl,^Rps"
57 """
58)
60parser.add_argument('input_file', help='Input .h5ad or .h5 file (10X Genomics format)')
61parser.add_argument('--output-dir', type=str, help='Output directory (default: <input_basename>_qc_results)')
62parser.add_argument('--mad-counts', type=float, default=DEFAULT_MAD_COUNTS, help=f'MAD threshold for total counts (default: {DEFAULT_MAD_COUNTS})')
63parser.add_argument('--mad-genes', type=float, default=DEFAULT_MAD_GENES, help=f'MAD threshold for gene counts (default: {DEFAULT_MAD_GENES})')
64parser.add_argument('--mad-mt', type=float, default=DEFAULT_MAD_MT, help=f'MAD threshold for mitochondrial percentage (default: {DEFAULT_MAD_MT})')
65parser.add_argument('--mt-threshold', type=float, default=DEFAULT_MT_THRESHOLD, help=f'Hard threshold for mitochondrial percentage (default: {DEFAULT_MT_THRESHOLD})')
66parser.add_argument('--min-cells', type=int, default=DEFAULT_MIN_CELLS, help=f'Minimum cells for gene filtering (default: {DEFAULT_MIN_CELLS})')
67parser.add_argument('--mt-pattern', type=str, default=DEFAULT_MT_PATTERN, help=f'Comma-separated mitochondrial gene prefixes (default: "{DEFAULT_MT_PATTERN}")')
68parser.add_argument('--ribo-pattern', type=str, default=DEFAULT_RIBO_PATTERN, help=f'Comma-separated ribosomal gene prefixes (default: "{DEFAULT_RIBO_PATTERN}")')
69parser.add_argument('--hb-pattern', type=str, default=DEFAULT_HB_PATTERN, help=f'Hemoglobin gene regex pattern (default: "{DEFAULT_HB_PATTERN}")')
71args = parser.parse_args()
73# Verify input file exists
74if not os.path.exists(args.input_file):
75 print(f"\nError: File '{args.input_file}' not found!")
76 sys.exit(1)
78input_file = args.input_file
79base_name = os.path.splitext(os.path.basename(input_file))[0]
81# Set up output directory
82if args.output_dir:
83 output_dir = args.output_dir
84else:
85 output_dir = f"{base_name}_qc_results"
87os.makedirs(output_dir, exist_ok=True)
88print(f"\nOutput directory: {output_dir}")
90# Display parameters
91print(f"\nParameters:")
92print(f" MAD thresholds: counts={args.mad_counts}, genes={args.mad_genes}, MT%={args.mad_mt}")
93print(f" MT hard threshold: {args.mt_threshold}%")
94print(f" Min cells for gene filtering: {args.min_cells}")
95print(f" Gene patterns: MT={args.mt_pattern}, Ribo={args.ribo_pattern}")
97# Load the data
98print("\n[1/5] Loading data...")
99file_ext = os.path.splitext(input_file)[1].lower()
101if file_ext == '.h5ad':
102 adata = ad.read_h5ad(input_file)
103 print(f"Loaded .h5ad file: {adata.n_obs} cells × {adata.n_vars} genes")
104elif file_ext == '.h5':
105 adata = sc.read_10x_h5(input_file)
106 print(f"Loaded 10X .h5 file: {adata.n_obs} cells × {adata.n_vars} genes")
107 # Make variable names unique (10X data sometimes has duplicate gene names)
108 adata.var_names_make_unique()
109else:
110 print(f"\nError: Unsupported file format '{file_ext}'. Expected .h5ad or .h5")
111 sys.exit(1)
113# Store original counts for comparison
114n_cells_original = adata.n_obs
115n_genes_original = adata.n_vars
117# Calculate QC metrics
118print("\n[2/5] Calculating QC metrics...")
119calculate_qc_metrics(adata, mt_pattern=args.mt_pattern,
120 ribo_pattern=args.ribo_pattern,
121 hb_pattern=args.hb_pattern,
122 inplace=True)
124print(f" Found {adata.var['mt'].sum()} mitochondrial genes (pattern: {args.mt_pattern})")
125print(f" Found {adata.var['ribo'].sum()} ribosomal genes (pattern: {args.ribo_pattern})")
126print(f" Found {adata.var['hb'].sum()} hemoglobin genes (pattern: {args.hb_pattern})")
128print_qc_summary(adata, label='QC Metrics Summary (before filtering)')
130# Create before-filtering visualizations
131print("\n[3/5] Creating QC visualizations...")
132before_plot = os.path.join(output_dir, 'qc_metrics_before_filtering.png')
133plot_qc_distributions(adata, before_plot, title='Quality Control Metrics - Before Filtering')
134print(f" Saved: {before_plot}")
136# Apply MAD-based filtering
137print("\n[4/5] Applying MAD-based filtering thresholds...")
139# Detect outliers for each metric
140adata.obs['outlier_counts'] = detect_outliers_mad(adata, 'total_counts', args.mad_counts)
141adata.obs['outlier_genes'] = detect_outliers_mad(adata, 'n_genes_by_counts', args.mad_genes)
142adata.obs['outlier_mt'] = detect_outliers_mad(adata, 'pct_counts_mt', args.mad_mt)
144# Apply hard threshold for mitochondrial content
145print(f"\n Applying hard threshold for mitochondrial content (>{args.mt_threshold}%):")
146high_mt_mask = apply_hard_threshold(adata, 'pct_counts_mt', args.mt_threshold, operator='>')
148# Combine MT filters (MAD + hard threshold)
149adata.obs['outlier_mt'] = adata.obs['outlier_mt'] | high_mt_mask
151# Overall filtering decision
152adata.obs['pass_qc'] = ~(
153 adata.obs['outlier_counts'] |
154 adata.obs['outlier_genes'] |
155 adata.obs['outlier_mt']
156)
158print(f"\n Total cells failing QC: {(~adata.obs['pass_qc']).sum()} ({(~adata.obs['pass_qc']).sum()/adata.n_obs*100:.2f}%)")
159print(f" Cells passing QC: {adata.obs['pass_qc'].sum()} ({adata.obs['pass_qc'].sum()/adata.n_obs*100:.2f}%)")
161# Visualize filtering thresholds
162outlier_masks = {
163 'total_counts': adata.obs['outlier_counts'].values,
164 'n_genes_by_counts': adata.obs['outlier_genes'].values,
165 'pct_counts_mt': adata.obs['outlier_mt'].values
166}
168thresholds = {
169 'total_counts': {'n_mads': args.mad_counts},
170 'n_genes_by_counts': {'n_mads': args.mad_genes},
171 'pct_counts_mt': {'n_mads': args.mad_mt, 'hard': args.mt_threshold}
172}
174threshold_plot = os.path.join(output_dir, 'qc_filtering_thresholds.png')
175plot_filtering_thresholds(adata, outlier_masks, thresholds, threshold_plot)
176print(f"\n Saved: {threshold_plot}")
178# Apply filtering
179print("\n[5/5] Applying filters...")
180adata_filtered = filter_cells(adata, adata.obs['pass_qc'].values, inplace=False)
181print(f" Cells after filtering: {adata_filtered.n_obs} (removed {n_cells_original - adata_filtered.n_obs})")
183# Filter genes
184print(f"\n Filtering genes detected in <{args.min_cells} cells...")
185filter_genes(adata_filtered, min_cells=args.min_cells, inplace=True)
186print(f" Genes after filtering: {adata_filtered.n_vars} (removed {n_genes_original - adata_filtered.n_vars})")
188# Generate summary statistics
189print("\n" + "=" * 80)
190print("QC Summary")
191print("=" * 80)
193print("\nBefore filtering:")
194print(f" Cells: {n_cells_original}")
195print(f" Genes: {n_genes_original}")
197print("\nAfter filtering:")
198print(f" Cells: {adata_filtered.n_obs} ({adata_filtered.n_obs/n_cells_original*100:.1f}% retained)")
199print(f" Genes: {adata_filtered.n_vars} ({adata_filtered.n_vars/n_genes_original*100:.1f}% retained)")
201print_qc_summary(adata_filtered, label='\nFiltered data QC metrics')
203# Create after-filtering visualizations
204after_plot = os.path.join(output_dir, 'qc_metrics_after_filtering.png')
205plot_qc_after_filtering(adata_filtered, after_plot)
206print(f"\n Saved: {after_plot}")
208# Save filtered data
209print("\nSaving filtered data...")
210output_filtered = os.path.join(output_dir, f'{base_name}_filtered.h5ad')
211output_with_qc = os.path.join(output_dir, f'{base_name}_with_qc.h5ad')
212adata_filtered.write(output_filtered)
213print(f" Saved: {output_filtered}")
215# Also save the unfiltered data with QC annotations
216adata.write(output_with_qc)
217print(f" Saved: {output_with_qc} (original data with QC annotations)")
219print("\n" + "=" * 80)
220print("Quality Control Analysis Complete!")
221print("=" * 80)
222print(f"\nAll results saved to: {output_dir}/")
223print("\nGenerated files:")
224print(" 1. qc_metrics_before_filtering.png - Initial QC visualizations")
225print(" 2. qc_filtering_thresholds.png - MAD-based threshold visualization")
226print(" 3. qc_metrics_after_filtering.png - Post-filtering QC visualizations")
227print(f" 4. {base_name}_filtered.h5ad - Filtered dataset")
228print(f" 5. {base_name}_with_qc.h5ad - Original dataset with QC annotations")
229print("\nNext steps:")
230print(" - Consider ambient RNA correction (SoupX)")
231print(" - Consider doublet detection (scDblFinder)")
232print(" - Proceed with normalization and downstream analysis")