Skills › Knowledge Work Plugins › single-cell-rna-qc › qc_core.py
qc_core.py
233 lines · bio-research/skills/single-cell-rna-qc/scripts/qc_core.py at 2574e29
1#!/usr/bin/env python3
2"""
3Core utility functions for single-cell RNA-seq quality control.
5This module provides building blocks for metrics calculation and filtering
6while following scverse best practices from:
7https://www.sc-best-practices.org/preprocessing_visualization/quality_control.html
8"""
10import anndata as ad
11import scanpy as sc
12import numpy as np
13from scipy.stats import median_abs_deviation
16def 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.
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)
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()
42 # Identify gene categories
43 mt_prefixes = tuple(mt_pattern.split(','))
44 adata.var['mt'] = adata.var_names.str.startswith(mt_prefixes)
46 ribo_prefixes = tuple(ribo_pattern.split(','))
47 adata.var['ribo'] = adata.var_names.str.startswith(ribo_prefixes)
49 adata.var['hb'] = adata.var_names.str.match(hb_pattern)
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 )
60 if not inplace:
61 return adata
64def detect_outliers_mad(adata, metric, n_mads, verbose=True):
65 """
66 Detect outliers using Median Absolute Deviation (MAD).
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)
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)
88 # Calculate bounds
89 lower = median - n_mads * mad
90 upper = median + n_mads * mad
92 # Identify outliers
93 outlier_mask = (metric_values < lower) | (metric_values > upper)
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}%)")
101 return outlier_mask
104def apply_hard_threshold(adata, metric, threshold, operator='>', verbose=True):
105 """
106 Apply a hard threshold filter.
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)
121 Returns
122 -------
123 np.ndarray
124 Boolean mask where True indicates cells to filter out
125 """
126 metric_values = adata.obs[metric]
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 '<='")
139 if verbose:
140 print(f" {metric} {operator} {threshold}:")
141 print(f" Cells filtered: {mask.sum()} ({mask.sum()/len(metric_values)*100:.2f}%)")
143 return mask
146def filter_cells(adata, mask, inplace=False):
147 """
148 Filter cells based on a boolean mask.
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)
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()
172def filter_genes(adata, min_cells=20, min_counts=None, inplace=True):
173 """
174 Filter genes based on detection thresholds.
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)
187 Returns
188 -------
189 AnnData or None
190 If inplace=False, returns filtered AnnData
191 """
192 if not inplace:
193 adata = adata.copy()
195 if min_cells is not None:
196 sc.pp.filter_genes(adata, min_cells=min_cells)
198 if min_counts is not None:
199 sc.pp.filter_genes(adata, min_counts=min_counts)
201 if not inplace:
202 return adata
205def print_qc_summary(adata, label=''):
206 """
207 Print summary statistics for QC metrics.
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}")
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}")
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}")
229 if 'pct_counts_mt' in adata.obs:
230 print(f" Mean mitochondrial %: {adata.obs['pct_counts_mt'].mean():.2f}%")
232 if 'pct_counts_ribo' in adata.obs:
233 print(f" Mean ribosomal %: {adata.obs['pct_counts_ribo'].mean():.2f}%")