Setting the file. One moment. Qc Plotting · Single Cell Rna Qc · anthropics/knowledge-work-plugins · Skills Docs22
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
scripts/qc_plotting.py
scripts/qc_plotting.py
Python·235 lines·8 KB
'Quality Control Metrics'
):
14 """
15 Create comprehensive QC distribution plots.
16
17 Parameters
18 ----------
19 adata : AnnData
20 Annotated data matrix with QC metrics
21 output_path : str
22 Path to save the figure
23 title : str
24 Figure title (default: 'Quality Control Metrics')
25 """
26 fig, axes = plt.subplots(3, 3, figsize=(15, 12))
27 fig.suptitle(title, fontsize=16, y=0.995)
28
29 # Row 1: Histograms
30 axes[0, 0].hist(adata.obs['total_counts'], bins=100, color='steelblue', edgecolor='black')
31 axes[0, 0].set_xlabel('Total counts per cell')
32 axes[0, 0].set_ylabel('Number of cells')
33 axes[0, 0].set_title('Distribution of Total Counts')
34 axes[0, 0].axvline(adata.obs['total_counts'].median(), color='red', linestyle='--', label='Median')
35 axes[0, 0].legend()
36
37 axes[0, 1].hist(adata.obs['n_genes_by_counts'], bins=100, color='forestgreen', edgecolor='black')
38 axes[0, 1].set_xlabel('Genes per cell')
39 axes[0, 1].set_ylabel('Number of cells')
40 axes[0, 1].set_title('Distribution of Detected Genes')
41 axes[0, 1].axvline(adata.obs['n_genes_by_counts'].median(), color='red', linestyle='--', label='Median')
42 axes[0, 1].legend()
43
44 axes[0, 2].hist(adata.obs['pct_counts_mt'], bins=100, color='coral', edgecolor='black')
45 axes[0, 2].set_xlabel('Mitochondrial %')
46 axes[0, 2].set_ylabel('Number of cells')
47 axes[0, 2].set_title('Distribution of Mitochondrial Content')
48 axes[0, 2].axvline(adata.obs['pct_counts_mt'].median(), color='red', linestyle='--', label='Median')
49 axes[0, 2].legend()
50
51 # Row 2: Violin plots
52 axes[1, 0].violinplot([adata.obs['total_counts']], positions=[0], showmeans=True, showmedians=True)
53 axes[1, 0].set_ylabel('Total counts')
54 axes[1, 0].set_title('Total Counts per Cell')
55 axes[1, 0].set_xticks([])
56
57 axes[1, 1].violinplot([adata.obs['n_genes_by_counts']], positions=[0], showmeans=True, showmedians=True)
58 axes[1, 1].set_ylabel('Genes detected')
59 axes[1, 1].set_title('Genes per Cell')
60 axes[1, 1].set_xticks([])
61
62 axes[1, 2].violinplot([adata.obs['pct_counts_mt']], positions=[0], showmeans=True, showmedians=True)
63 axes[1, 2].set_ylabel('Mitochondrial %')
64 axes[1, 2].set_title('Mitochondrial Content')
65 axes[1, 2].set_xticks([])
66
67 # Row 3: Scatter plots
68 scatter1 = axes[2, 0].scatter(
69 adata.obs['total_counts'],
70 adata.obs['n_genes_by_counts'],
71 c=adata.obs['pct_counts_mt'],
72 cmap='viridis',
73 alpha=0.5,
74 s=10
75 )
76 axes[2, 0].set_xlabel('Total counts')
77 axes[2, 0].set_ylabel('Genes detected')
78 axes[2, 0].set_title('Counts vs Genes (colored by MT%)')
79 plt.colorbar(scatter1, ax=axes[2, 0], label='MT %')
80
81 axes[2, 1].scatter(
82 adata.obs['total_counts'],
83 adata.obs['pct_counts_mt'],
84 alpha=0.5,
85 s=10,
86 color='coral'
87 )
88 axes[2, 1].set_xlabel('Total counts')
89 axes[2, 1].set_ylabel('Mitochondrial %')
90 axes[2, 1].set_title('Total Counts vs Mitochondrial %')
91
92 axes[2, 2].scatter(
93 adata.obs['n_genes_by_counts'],
94 adata.obs['pct_counts_mt'],
95 alpha=0.5,
96 s=10,
97 color='forestgreen'
98 )
99 axes[2, 2].set_xlabel('Genes detected')
100 axes[2, 2].set_ylabel('Mitochondrial %')
101 axes[2, 2].set_title('Genes vs Mitochondrial %')
102
103 plt.tight_layout()
104 plt.savefig(output_path, dpi=300, bbox_inches='tight')
105 plt.close()
106
107
108def plot_filtering_thresholds(adata, outlier_masks, thresholds, output_path):
109 """
110 Visualize filtering thresholds overlaid on distributions.
111
112 Parameters
113 ----------
114 adata : AnnData
115 Annotated data matrix with QC metrics
116 outlier_masks : dict
117 Dictionary mapping metric names to boolean outlier masks
118 Example: {'total_counts': mask1, 'n_genes_by_counts': mask2, 'pct_counts_mt': mask3}
119 thresholds : dict
120 Dictionary with threshold information for each metric
121 Example: {'total_counts': {'n_mads': 5}, 'pct_counts_mt': {'n_mads': 3, 'hard': 8}}
122 output_path : str
123 Path to save the figure
124 """
125 fig, axes = plt.subplots(1, 3, figsize=(15, 4))
126 fig.suptitle('MAD-Based Filtering Thresholds', fontsize=16)
127
128 # Helper function to plot with thresholds
129 def plot_with_threshold(ax, metric, outlier_mask, n_mads, hard_threshold=None):
130 data = adata.obs[metric]
131 median = np.median(data)
132 mad = median_abs_deviation(data)
133 lower = median - n_mads * mad
134 upper = median + n_mads * mad
135
136 ax.hist(data[~outlier_mask], bins=100, alpha=0.7, label='Pass QC', color='steelblue')
137 ax.hist(data[outlier_mask], bins=100, alpha=0.7, label='Fail QC', color='coral')
138 ax.axvline(lower, color='red', linestyle='--', linewidth=2, label=f'Thresholds ({n_mads} MADs)')
139 ax.axvline(upper, color='red', linestyle='--', linewidth=2)
140
141 if hard_threshold is not None:
142 ax.axvline(hard_threshold, color='darkred', linestyle=':', linewidth=2,
143 label=f'Hard threshold ({hard_threshold})')
144
145 ax.set_xlabel(metric.replace('_', ' ').title())
146 ax.set_ylabel('Number of cells')
147 ax.legend()
148
149 # Plot each metric
150 metrics = [
151 ('total_counts', 'Total Counts'),
152 ('n_genes_by_counts', 'Genes Detected'),
153 ('pct_counts_mt', 'Mitochondrial %')
154 ]
155
156 for idx, (metric, label) in enumerate(metrics):
157 if metric in outlier_masks and metric in thresholds:
158 hard = thresholds[metric].get('hard', None)
159 plot_with_threshold(axes[idx], metric, outlier_masks[metric],
160 thresholds[metric]['n_mads'], hard)
161
162 plt.tight_layout()
163 plt.savefig(output_path, dpi=300, bbox_inches='tight')
164 plt.close()
165
166
167def plot_qc_after_filtering(adata, output_path):
168 """
169 Create QC plots for filtered data (simplified version without outlier overlay).
170
171 Parameters
172 ----------
173 adata : AnnData
174 Filtered annotated data matrix with QC metrics
175 output_path : str
176 Path to save the figure
177 """
178 fig, axes = plt.subplots(2, 3, figsize=(15, 8))
179 fig.suptitle('Quality Control Metrics - After Filtering', fontsize=16, y=0.995)
180
181 # Row 1: Histograms
182 axes[0, 0].hist(adata.obs['total_counts'], bins=100, color='steelblue', edgecolor='black')
183 axes[0, 0].set_xlabel('Total counts per cell')
184 axes[0, 0].set_ylabel('Number of cells')
185 axes[0, 0].set_title('Distribution of Total Counts')
186
187 axes[0, 1].hist(adata.obs['n_genes_by_counts'], bins=100, color='forestgreen', edgecolor='black')
188 axes[0, 1].set_xlabel('Genes per cell')
189 axes[0, 1].set_ylabel('Number of cells')
190 axes[0, 1].set_title('Distribution of Detected Genes')
191
192 axes[0, 2].hist(adata.obs['pct_counts_mt'], bins=100, color='coral', edgecolor='black')
193 axes[0, 2].set_xlabel('Mitochondrial %')
194 axes[0, 2].set_ylabel('Number of cells')
195 axes[0, 2].set_title('Distribution of Mitochondrial Content')
196
197 # Row 2: Scatter plots
198 scatter1 = axes[1, 0].scatter(
199 adata.obs['total_counts'],
200 adata.obs['n_genes_by_counts'],
201 c=adata.obs['pct_counts_mt'],
202 cmap='viridis',
203 alpha=0.5,
204 s=10
205 )
206 axes[1, 0].set_xlabel('Total counts')
207 axes[1, 0].set_ylabel('Genes detected')
208 axes[1, 0].set_title('Counts vs Genes (colored by MT%)')
209 plt.colorbar(scatter1, ax=axes[1, 0], label='MT %')
210
211 axes[1, 1].scatter(
212 adata.obs['total_counts'],
213 adata.obs['pct_counts_mt'],
214 alpha=0.5,
215 s=10,
216 color='coral'
217 )
218 axes[1, 1].set_xlabel('Total counts')
219 axes[1, 1].set_ylabel('Mitochondrial %')
220 axes[1, 1].set_title('Total Counts vs Mitochondrial %')
221
222 axes[1, 2].scatter(
223 adata.obs['n_genes_by_counts'],
224 adata.obs['pct_counts_mt'],
225 alpha=0.5,
226 s=10,
227 color='forestgreen'
228 )
229 axes[1, 2].set_xlabel('Genes detected')
230 axes[1, 2].set_ylabel('Mitochondrial %')
231 axes[1, 2].set_title('Genes vs Mitochondrial %')
232
233 plt.tight_layout()
234 plt.savefig(output_path, dpi=300, bbox_inches='tight')
235 plt.close()