Setting the file. One moment.
Qc Core · 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
Next
Script Qc Plotting
scripts/ qc_core.py
Python · 233 lines · 7 KB
scipy.stats
import
median_abs_deviation
14
15
16 def calculate_qc_metrics (adata, mt_pattern = 'mt-,MT-' , ribo_pattern = 'Rpl,Rps,RPL,RPS' ,
17 hb_pattern = '^Hb[^(p)]|^HB[^(P)]' , inplace = True ):
18 """
19 Calculate QC metrics for single-cell RNA-seq data.
20
21 Parameters
22 ----------
23 adata : AnnData
24 Annotated data matrix
25 mt_pattern : str
26 Comma-separated mitochondrial gene prefixes (default: 'mt-,MT-')
27 ribo_pattern : str
28 Comma-separated ribosomal gene prefixes (default: 'Rpl,Rps,RPL,RPS')
29 hb_pattern : str
30 Regex pattern for hemoglobin genes (default: '^Hb[^(p)]|^HB[^(P)]')
31 inplace : bool
32 Modify adata in place (default: True)
33
34 Returns
35 -------
36 AnnData or None
37 If inplace=False, returns modified AnnData. Otherwise modifies in place.
38 """
39 if not inplace:
40 adata = adata.copy()
41
42 # Identify gene categories
43 mt_prefixes = tuple (mt_pattern.split( ',' ))
44 adata.var[ 'mt' ] = adata.var_names.str.startswith(mt_prefixes)
45
46 ribo_prefixes = tuple (ribo_pattern.split( ',' ))
47 adata.var[ 'ribo' ] = adata.var_names.str.startswith(ribo_prefixes)
48
49 adata.var[ 'hb' ] = adata.var_names.str.match(hb_pattern)
50
51 # Calculate QC metrics
52 sc.pp.calculate_qc_metrics(
53 adata,
54 qc_vars = [ 'mt' , 'ribo' , 'hb' ],
55 percent_top = None ,
56 log1p = False ,
57 inplace = True
58 )
59
60 if not inplace:
61 return adata
62
63
64 def detect_outliers_mad (adata, metric, n_mads, verbose = True ):
65 """
66 Detect outliers using Median Absolute Deviation (MAD).
67
68 Parameters
69 ----------
70 adata : AnnData
71 Annotated data matrix with QC metrics
72 metric : str
73 Column name in adata.obs to use for outlier detection
74 n_mads : float
75 Number of MADs to use as threshold
76 verbose : bool
77 Print outlier statistics (default: True)
78
79 Returns
80 -------
81 np.ndarray
82 Boolean mask where True indicates outliers
83 """
84 metric_values = adata.obs[metric]
85 median = np.median(metric_values)
86 mad = median_abs_deviation(metric_values)
87
88 # Calculate bounds
89 lower = median - n_mads * mad
90 upper = median + n_mads * mad
91
92 # Identify outliers
93 outlier_mask = (metric_values < lower) | (metric_values > upper)
94
95 if verbose:
96 print ( f " { metric } :" )
97 print ( f " Median: { median :.2f} , MAD: { mad :.2f} " )
98 print ( f " Bounds: [ { lower :.2f} , { upper :.2f} ] ( { n_mads } MADs)" )
99 print ( f " Outliers: { outlier_mask.sum() } cells ( { outlier_mask.sum() / len (metric_values) * 100 :.2f} %)" )
100
101 return outlier_mask
102
103
104 def apply_hard_threshold (adata, metric, threshold, operator = '>' , verbose = True ):
105 """
106 Apply a hard threshold filter.
107
108 Parameters
109 ----------
110 adata : AnnData
111 Annotated data matrix
112 metric : str
113 Column name in adata.obs to filter on
114 threshold : float
115 Threshold value
116 operator : str
117 Comparison operator: '>', '<', '>=', '<=' (default: '>')
118 verbose : bool
119 Print filtering statistics (default: True)
120
121 Returns
122 -------
123 np.ndarray
124 Boolean mask where True indicates cells to filter out
125 """
126 metric_values = adata.obs[metric]
127
128 if operator == '>' :
129 mask = metric_values > threshold
130 elif operator == '<' :
131 mask = metric_values < threshold
132 elif operator == '>=' :
133 mask = metric_values >= threshold
134 elif operator == '<=' :
135 mask = metric_values <= threshold
136 else :
137 raise ValueError ( f "Invalid operator: { operator } . Use '>', '<', '>=', or '<='" )
138
139 if verbose:
140 print ( f " { metric } { operator } { threshold } :" )
141 print ( f " Cells filtered: { mask.sum() } ( { mask.sum() / len (metric_values) * 100 :.2f} %)" )
142
143 return mask
144
145
146 def filter_cells (adata, mask, inplace = False ):
147 """
148 Filter cells based on a boolean mask.
149
150 Parameters
151 ----------
152 adata : AnnData
153 Annotated data matrix
154 mask : np.ndarray or pd.Series
155 Boolean mask where True indicates cells to KEEP
156 inplace : bool
157 Modify adata in place (default: False)
158
159 Returns
160 -------
161 AnnData
162 Filtered AnnData object
163 """
164 if inplace:
165 # This is actually a bit tricky - AnnData doesn't support true inplace filtering
166 # Return filtered copy which caller should reassign
167 return adata[mask].copy()
168 else :
169 return adata[mask].copy()
170
171
172 def filter_genes (adata, min_cells = 20 , min_counts = None , inplace = True ):
173 """
174 Filter genes based on detection thresholds.
175
176 Parameters
177 ----------
178 adata : AnnData
179 Annotated data matrix
180 min_cells : int
181 Minimum number of cells a gene must be detected in (default: 20)
182 min_counts : int, optional
183 Minimum total counts across all cells
184 inplace : bool
185 Modify adata in place (default: True)
186
187 Returns
188 -------
189 AnnData or None
190 If inplace=False, returns filtered AnnData
191 """
192 if not inplace:
193 adata = adata.copy()
194
195 if min_cells is not None :
196 sc.pp.filter_genes(adata, min_cells = min_cells)
197
198 if min_counts is not None :
199 sc.pp.filter_genes(adata, min_counts = min_counts)
200
201 if not inplace:
202 return adata
203
204
205 def print_qc_summary (adata, label = '' ):
206 """
207 Print summary statistics for QC metrics.
208
209 Parameters
210 ----------
211 adata : AnnData
212 Annotated data matrix with QC metrics
213 label : str
214 Label to prepend to output (e.g., 'Before filtering', 'After filtering')
215 """
216 if label:
217 print ( f " \n{ label } :" )
218 print ( f " Cells: { adata.n_obs } " )
219 print ( f " Genes: { adata.n_vars } " )
220
221 if 'total_counts' in adata.obs:
222 print ( f " Mean counts per cell: { adata.obs[ 'total_counts' ].mean() :.0f} " )
223 print ( f " Median counts per cell: { adata.obs[ 'total_counts' ].median() :.0f} " )
224
225 if 'n_genes_by_counts' in adata.obs:
226 print ( f " Mean genes per cell: { adata.obs[ 'n_genes_by_counts' ].mean() :.0f} " )
227 print ( f " Median genes per cell: { adata.obs[ 'n_genes_by_counts' ].median() :.0f} " )
228
229 if 'pct_counts_mt' in adata.obs:
230 print ( f " Mean mitochondrial %: { adata.obs[ 'pct_counts_mt' ].mean() :.2f} %" )
231
232 if 'pct_counts_ribo' in adata.obs:
233 print ( f " Mean ribosomal %: { adata.obs[ 'pct_counts_ribo' ].mean() :.2f} %" )