Setting the file. One moment.
Qc Analysis · Single Cell Rna Qc · anthropics/knowledge-work-plugins · Skills Docs
Back to the top of the page 22
Validate Data
Tech Debt
62
Recruiting Pipeline
71
Vendor Check
125
Zoom Meeting SDK Web
88
Vendor Review
181
Create An Asset
Video Sdk/web
(opens in a new tab)
scripts/ qc_analysis.py
Python · 232 lines · 9 KB
import
os
15 import argparse
16
17 # Import our modular utilities
18 from 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 )
26 from qc_plotting import (
27 plot_qc_distributions,
28 plot_filtering_thresholds,
29 plot_qc_after_filtering
30 )
31
32 print ( "=" * 80 )
33 print ( "Single-Cell RNA-seq Quality Control Analysis" )
34 print ( "=" * 80 )
35
36 # Default parameters (single source of truth)
37 DEFAULT_MAD_COUNTS = 5
38 DEFAULT_MAD_GENES = 5
39 DEFAULT_MAD_MT = 3
40 DEFAULT_MT_THRESHOLD = 8
41 DEFAULT_MIN_CELLS = 20
42 DEFAULT_MT_PATTERN = 'mt-,MT-'
43 DEFAULT_RIBO_PATTERN = 'Rpl,Rps,RPL,RPS'
44 DEFAULT_HB_PATTERN = '^Hb[^(p)]|^HB[^(P)]'
45
46 # Parse command-line arguments
47 parser = argparse.ArgumentParser(
48 description = 'Quality Control Analysis for Single-Cell RNA-seq Data' ,
49 formatter_class = argparse.RawDescriptionHelpFormatter,
50 epilog = """
51 Examples:
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 )
59
60 parser.add_argument( 'input_file' , help = 'Input .h5ad or .h5 file (10X Genomics format)' )
61 parser.add_argument( '--output-dir' , type = str , help = 'Output directory (default: <input_basename>_qc_results)' )
62 parser.add_argument( '--mad-counts' , type = float , default = DEFAULT_MAD_COUNTS , help = f 'MAD threshold for total counts (default: { DEFAULT_MAD_COUNTS } )' )
63 parser.add_argument( '--mad-genes' , type = float , default = DEFAULT_MAD_GENES , help = f 'MAD threshold for gene counts (default: { DEFAULT_MAD_GENES } )' )
64 parser.add_argument( '--mad-mt' , type = float , default = DEFAULT_MAD_MT , help = f 'MAD threshold for mitochondrial percentage (default: { DEFAULT_MAD_MT } )' )
65 parser.add_argument( '--mt-threshold' , type = float , default = DEFAULT_MT_THRESHOLD , help = f 'Hard threshold for mitochondrial percentage (default: { DEFAULT_MT_THRESHOLD } )' )
66 parser.add_argument( '--min-cells' , type = int , default = DEFAULT_MIN_CELLS , help = f 'Minimum cells for gene filtering (default: { DEFAULT_MIN_CELLS } )' )
67 parser.add_argument( '--mt-pattern' , type = str , default = DEFAULT_MT_PATTERN , help = f 'Comma-separated mitochondrial gene prefixes (default: " { DEFAULT_MT_PATTERN } ")' )
68 parser.add_argument( '--ribo-pattern' , type = str , default = DEFAULT_RIBO_PATTERN , help = f 'Comma-separated ribosomal gene prefixes (default: " { DEFAULT_RIBO_PATTERN } ")' )
69 parser.add_argument( '--hb-pattern' , type = str , default = DEFAULT_HB_PATTERN , help = f 'Hemoglobin gene regex pattern (default: " { DEFAULT_HB_PATTERN } ")' )
70
71 args = parser.parse_args()
72
73 # Verify input file exists
74 if not os.path.exists(args.input_file):
75 print ( f " \n Error: File ' { args.input_file } ' not found!" )
76 sys.exit( 1 )
77
78 input_file = args.input_file
79 base_name = os.path.splitext(os.path.basename(input_file))[ 0 ]
80
81 # Set up output directory
82 if args.output_dir:
83 output_dir = args.output_dir
84 else :
85 output_dir = f " { base_name } _qc_results"
86
87 os.makedirs(output_dir, exist_ok = True )
88 print ( f " \n Output directory: { output_dir } " )
89
90 # Display parameters
91 print ( f " \n Parameters:" )
92 print ( f " MAD thresholds: counts= { args.mad_counts } , genes= { args.mad_genes } , MT%= { args.mad_mt } " )
93 print ( f " MT hard threshold: { args.mt_threshold } %" )
94 print ( f " Min cells for gene filtering: { args.min_cells } " )
95 print ( f " Gene patterns: MT= { args.mt_pattern } , Ribo= { args.ribo_pattern } " )
96
97 # Load the data
98 print ( " \n [1/5] Loading data..." )
99 file_ext = os.path.splitext(input_file)[ 1 ].lower()
100
101 if file_ext == '.h5ad' :
102 adata = ad.read_h5ad(input_file)
103 print ( f "Loaded .h5ad file: { adata.n_obs } cells × { adata.n_vars } genes" )
104 elif 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()
109 else :
110 print ( f " \n Error: Unsupported file format ' { file_ext } '. Expected .h5ad or .h5" )
111 sys.exit( 1 )
112
113 # Store original counts for comparison
114 n_cells_original = adata.n_obs
115 n_genes_original = adata.n_vars
116
117 # Calculate QC metrics
118 print ( " \n [2/5] Calculating QC metrics..." )
119 calculate_qc_metrics(adata, mt_pattern = args.mt_pattern,
120 ribo_pattern = args.ribo_pattern,
121 hb_pattern = args.hb_pattern,
122 inplace = True )
123
124 print ( f " Found { adata.var[ 'mt' ].sum() } mitochondrial genes (pattern: { args.mt_pattern } )" )
125 print ( f " Found { adata.var[ 'ribo' ].sum() } ribosomal genes (pattern: { args.ribo_pattern } )" )
126 print ( f " Found { adata.var[ 'hb' ].sum() } hemoglobin genes (pattern: { args.hb_pattern } )" )
127
128 print_qc_summary(adata, label = 'QC Metrics Summary (before filtering)' )
129
130 # Create before-filtering visualizations
131 print ( " \n [3/5] Creating QC visualizations..." )
132 before_plot = os.path.join(output_dir, 'qc_metrics_before_filtering.png' )
133 plot_qc_distributions(adata, before_plot, title = 'Quality Control Metrics - Before Filtering' )
134 print ( f " Saved: { before_plot } " )
135
136 # Apply MAD-based filtering
137 print ( " \n [4/5] Applying MAD-based filtering thresholds..." )
138
139 # Detect outliers for each metric
140 adata.obs[ 'outlier_counts' ] = detect_outliers_mad(adata, 'total_counts' , args.mad_counts)
141 adata.obs[ 'outlier_genes' ] = detect_outliers_mad(adata, 'n_genes_by_counts' , args.mad_genes)
142 adata.obs[ 'outlier_mt' ] = detect_outliers_mad(adata, 'pct_counts_mt' , args.mad_mt)
143
144 # Apply hard threshold for mitochondrial content
145 print ( f " \n Applying hard threshold for mitochondrial content (> { args.mt_threshold } %):" )
146 high_mt_mask = apply_hard_threshold(adata, 'pct_counts_mt' , args.mt_threshold, operator = '>' )
147
148 # Combine MT filters (MAD + hard threshold)
149 adata.obs[ 'outlier_mt' ] = adata.obs[ 'outlier_mt' ] | high_mt_mask
150
151 # Overall filtering decision
152 adata.obs[ 'pass_qc' ] = ~ (
153 adata.obs[ 'outlier_counts' ] |
154 adata.obs[ 'outlier_genes' ] |
155 adata.obs[ 'outlier_mt' ]
156 )
157
158 print ( f " \n Total cells failing QC: { ( ~ adata.obs[ 'pass_qc' ]).sum() } ( { ( ~ adata.obs[ 'pass_qc' ]).sum() / adata.n_obs * 100 :.2f} %)" )
159 print ( f " Cells passing QC: { adata.obs[ 'pass_qc' ].sum() } ( { adata.obs[ 'pass_qc' ].sum() / adata.n_obs * 100 :.2f} %)" )
160
161 # Visualize filtering thresholds
162 outlier_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 }
167
168 thresholds = {
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 }
173
174 threshold_plot = os.path.join(output_dir, 'qc_filtering_thresholds.png' )
175 plot_filtering_thresholds(adata, outlier_masks, thresholds, threshold_plot)
176 print ( f " \n Saved: { threshold_plot } " )
177
178 # Apply filtering
179 print ( " \n [5/5] Applying filters..." )
180 adata_filtered = filter_cells(adata, adata.obs[ 'pass_qc' ].values, inplace = False )
181 print ( f " Cells after filtering: { adata_filtered.n_obs } (removed { n_cells_original - adata_filtered.n_obs } )" )
182
183 # Filter genes
184 print ( f " \n Filtering genes detected in < { args.min_cells } cells..." )
185 filter_genes(adata_filtered, min_cells = args.min_cells, inplace = True )
186 print ( f " Genes after filtering: { adata_filtered.n_vars } (removed { n_genes_original - adata_filtered.n_vars } )" )
187
188 # Generate summary statistics
189 print ( " \n " + "=" * 80 )
190 print ( "QC Summary" )
191 print ( "=" * 80 )
192
193 print ( " \n Before filtering:" )
194 print ( f " Cells: { n_cells_original } " )
195 print ( f " Genes: { n_genes_original } " )
196
197 print ( " \n After filtering:" )
198 print ( f " Cells: { adata_filtered.n_obs } ( { adata_filtered.n_obs / n_cells_original * 100 :.1f} % retained)" )
199 print ( f " Genes: { adata_filtered.n_vars } ( { adata_filtered.n_vars / n_genes_original * 100 :.1f} % retained)" )
200
201 print_qc_summary(adata_filtered, label = ' \n Filtered data QC metrics' )
202
203 # Create after-filtering visualizations
204 after_plot = os.path.join(output_dir, 'qc_metrics_after_filtering.png' )
205 plot_qc_after_filtering(adata_filtered, after_plot)
206 print ( f " \n Saved: { after_plot } " )
207
208 # Save filtered data
209 print ( " \n Saving filtered data..." )
210 output_filtered = os.path.join(output_dir, f ' { base_name } _filtered.h5ad' )
211 output_with_qc = os.path.join(output_dir, f ' { base_name } _with_qc.h5ad' )
212 adata_filtered.write(output_filtered)
213 print ( f " Saved: { output_filtered } " )
214
215 # Also save the unfiltered data with QC annotations
216 adata.write(output_with_qc)
217 print ( f " Saved: { output_with_qc } (original data with QC annotations)" )
218
219 print ( " \n " + "=" * 80 )
220 print ( "Quality Control Analysis Complete!" )
221 print ( "=" * 80 )
222 print ( f " \n All results saved to: { output_dir } /" )
223 print ( " \n Generated files:" )
224 print ( " 1. qc_metrics_before_filtering.png - Initial QC visualizations" )
225 print ( " 2. qc_filtering_thresholds.png - MAD-based threshold visualization" )
226 print ( " 3. qc_metrics_after_filtering.png - Post-filtering QC visualizations" )
227 print ( f " 4. { base_name } _filtered.h5ad - Filtered dataset" )
228 print ( f " 5. { base_name } _with_qc.h5ad - Original dataset with QC annotations" )
229 print ( " \n Next steps:" )
230 print ( " - Consider ambient RNA correction (SoupX)" )
231 print ( " - Consider doublet detection (scDblFinder)" )
232 print ( " - Proceed with normalization and downstream analysis" )