|
| 1 | +import os |
| 2 | +import h5py |
| 3 | +import numpy as np |
| 4 | +import matplotlib.pyplot as plt |
| 5 | +from scipy.sparse import csr_matrix |
| 6 | + |
| 7 | + |
| 8 | +def diagnose_qc(h5ad_path, output_path, min_umis=500, min_genes=200, max_mt=0.15): |
| 9 | + print(f"Loading {h5ad_path}...") |
| 10 | + |
| 11 | + with h5py.File(h5ad_path, "r") as f: |
| 12 | + # Load barcodes and gene names |
| 13 | + barcodes = [ |
| 14 | + b.decode("utf-8") if isinstance(b, bytes) else b |
| 15 | + for b in f["obs"]["_index"][:] |
| 16 | + ] |
| 17 | + gene_names = [ |
| 18 | + g.decode("utf-8") if isinstance(g, bytes) else g |
| 19 | + for g in f["var"]["_index"][:] |
| 20 | + ] |
| 21 | + |
| 22 | + # Load expression matrix X |
| 23 | + X_group = f["X"] |
| 24 | + if isinstance(X_group, h5py.Group): |
| 25 | + data = X_group["data"][:] |
| 26 | + indices = X_group["indices"][:] |
| 27 | + indptr = X_group["indptr"][:] |
| 28 | + X = csr_matrix( |
| 29 | + (data, indices, indptr), shape=(len(barcodes), len(gene_names)) |
| 30 | + ) |
| 31 | + else: |
| 32 | + X = f["X"][:] |
| 33 | + |
| 34 | + # Spatial coordinates (usually in obsm/spatial) |
| 35 | + if "obsm" in f and "spatial" in f["obsm"]: |
| 36 | + coords = f["obsm"]["spatial"][:] |
| 37 | + else: |
| 38 | + # Fallback if spatial not found |
| 39 | + coords = np.zeros((len(barcodes), 2)) |
| 40 | + print("Warning: Spatial coordinates not found in obsm/spatial") |
| 41 | + |
| 42 | + # Calculate QC metrics |
| 43 | + n_counts = np.array(X.sum(axis=1)).flatten() |
| 44 | + n_genes = np.array((X > 0).sum(axis=1)).flatten() |
| 45 | + |
| 46 | + # MT fraction |
| 47 | + # Robust MT detection: check for mt- or MT- anywhere, but prioritize common patterns |
| 48 | + mt_genes = [ |
| 49 | + i |
| 50 | + for i, name in enumerate(gene_names) |
| 51 | + if "mt-" in name.lower() or "mt:" in name.lower() |
| 52 | + ] |
| 53 | + if mt_genes: |
| 54 | + print(f"Found {len(mt_genes)} mitochondrial genes.") |
| 55 | + mt_counts = np.array(X[:, mt_genes].sum(axis=1)).flatten() |
| 56 | + pct_counts_mt = mt_counts / (n_counts + 1e-9) |
| 57 | + else: |
| 58 | + pct_counts_mt = np.zeros_like(n_counts) |
| 59 | + print("No mitochondrial genes found.") |
| 60 | + |
| 61 | + # Filter mask |
| 62 | + pass_umi = n_counts >= min_umis |
| 63 | + pass_gene = n_genes >= min_genes |
| 64 | + pass_mt = pct_counts_mt <= max_mt |
| 65 | + |
| 66 | + keep_mask = pass_umi & pass_gene & pass_mt |
| 67 | + |
| 68 | + total_spots = len(barcodes) |
| 69 | + kept_spots = np.sum(keep_mask) |
| 70 | + |
| 71 | + print(f"Total spots: {total_spots}") |
| 72 | + print(f"Kept spots: {kept_spots} ({kept_spots/total_spots:.1%})") |
| 73 | + print(f"Filtered: {total_spots - kept_spots}") |
| 74 | + print(f" - Low UMI (<{min_umis}): {np.sum(~pass_umi)}") |
| 75 | + print(f" - Low Genes (<{min_genes}): {np.sum(~pass_gene)}") |
| 76 | + print(f" - High MT (>{max_mt:.1%}): {np.sum(~pass_mt)}") |
| 77 | + |
| 78 | + # Plotting |
| 79 | + fig, axes = plt.subplots(2, 2, figsize=(15, 12)) |
| 80 | + |
| 81 | + # 1. UMI vs Genes |
| 82 | + axes[0, 0].scatter(n_counts, n_genes, c=keep_mask, cmap="RdYlGn", alpha=0.5, s=10) |
| 83 | + axes[0, 0].axvline( |
| 84 | + min_umis, color="red", linestyle="--", label=f"Min UMI={min_umis}" |
| 85 | + ) |
| 86 | + axes[0, 0].axhline( |
| 87 | + min_genes, color="blue", linestyle="--", label=f"Min Genes={min_genes}" |
| 88 | + ) |
| 89 | + axes[0, 0].set_xlabel("Total UMI counts") |
| 90 | + axes[0, 0].set_ylabel("Number of detected genes") |
| 91 | + axes[0, 0].set_title("QC: UMI vs Genes") |
| 92 | + axes[0, 0].legend() |
| 93 | + |
| 94 | + # 2. MT Fraction distribution |
| 95 | + axes[0, 1].hist(pct_counts_mt, bins=50, color="gray", alpha=0.7) |
| 96 | + axes[0, 1].axvline( |
| 97 | + max_mt, color="red", linestyle="--", label=f"Max MT={max_mt:.0%}" |
| 98 | + ) |
| 99 | + axes[0, 1].set_xlabel("Mitochondrial Fraction") |
| 100 | + axes[0, 1].set_ylabel("Count") |
| 101 | + axes[0, 1].set_title("QC: MT Fraction Distribution") |
| 102 | + axes[0, 1].legend() |
| 103 | + |
| 104 | + # 3. Spatial: Before (Total) |
| 105 | + axes[1, 0].scatter( |
| 106 | + coords[:, 0], coords[:, 1], c="lightgray", s=15, label="All Spots" |
| 107 | + ) |
| 108 | + axes[1, 0].scatter( |
| 109 | + coords[keep_mask, 0], coords[keep_mask, 1], c="green", s=15, label="Pass QC" |
| 110 | + ) |
| 111 | + axes[1, 0].set_title("Spatial Distribution: Kept vs Filtered") |
| 112 | + axes[1, 0].set_aspect("equal") |
| 113 | + axes[1, 0].legend() |
| 114 | + |
| 115 | + # 4. Summary Table (as text) |
| 116 | + stats_text = ( |
| 117 | + f"Sample: {os.path.basename(h5ad_path)}\n\n" |
| 118 | + f"Total Spots: {total_spots}\n" |
| 119 | + f"Kept Spots: {kept_spots} ({kept_spots/total_spots:.1%})\n" |
| 120 | + f"Filtered: {total_spots - kept_spots}\n\n" |
| 121 | + f"Thresholds:\n" |
| 122 | + f"Min UMI: {min_umis}\n" |
| 123 | + f"Min Genes: {min_genes}\n" |
| 124 | + f"Max MT: {max_mt:.0%}" |
| 125 | + ) |
| 126 | + axes[1, 1].text(0.1, 0.5, stats_text, fontsize=14, family="monospace") |
| 127 | + axes[1, 1].axis("off") |
| 128 | + |
| 129 | + plt.tight_layout() |
| 130 | + plt.savefig(output_path) |
| 131 | + print(f"Plot saved to {output_path}") |
| 132 | + |
| 133 | + |
| 134 | +if __name__ == "__main__": |
| 135 | + import argparse |
| 136 | + |
| 137 | + parser = argparse.ArgumentParser() |
| 138 | + parser.add_argument("--sample", type=str, default="MEND29") |
| 139 | + args = parser.parse_args() |
| 140 | + |
| 141 | + sample_id = args.sample |
| 142 | + h5ad_file = f"A:\\hest_data\\st\\{sample_id}.h5ad" |
| 143 | + output_file = f"qc_diagnosis_{sample_id}.png" |
| 144 | + |
| 145 | + if not os.path.exists(h5ad_file): |
| 146 | + print(f"Error: {h5ad_file} not found.") |
| 147 | + else: |
| 148 | + diagnose_qc(h5ad_file, output_file) |
0 commit comments