78 lines
2.7 KiB
Python
78 lines
2.7 KiB
Python
#!/usr/bin/env python
|
||
# -*- coding: utf-8 -*-
|
||
"""
|
||
01_explore_raw.py — GSE292269 原始矩阵初步探索
|
||
- 加载 3 个样本的 raw feature-barcode 矩阵
|
||
- 计算每个 barcode 的 UMI/基因数分布，绘制 knee plot 辅助确定空液滴过滤阈值
|
||
- 输出：output/01_qc/*.png（300 ppi）+ 汇总打印
|
||
"""
|
||
import gzip
|
||
import numpy as np
|
||
import pandas as pd
|
||
import scanpy as sc
|
||
import scipy.io as sio
|
||
import matplotlib
|
||
matplotlib.use("Agg")
|
||
import matplotlib.pyplot as plt
|
||
from pathlib import Path
|
||
|
||
DATA = Path("data/GSE292269")
|
||
OUT = Path("output/01_qc")
|
||
OUT.mkdir(parents=True, exist_ok=True)
|
||
SAMPLES = ["WT", "Opa1V291D_S1", "Opa1V291D_S2"]
|
||
|
||
plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300, "font.size": 9})
|
||
|
||
|
||
def load_raw(sample: str) -> sc.AnnData:
|
||
d = DATA / sample
|
||
with gzip.open(d / "matrix.mtx.gz", "rb") as f:
|
||
X = sio.mmread(f).tocsr().T # barcodes x genes
|
||
bc = pd.read_csv(d / "barcodes.tsv.gz", header=None, sep="\t")[0].astype(str).values
|
||
feat = pd.read_csv(d / "features.tsv.gz", header=None, sep="\t")
|
||
adata = sc.AnnData(X=X)
|
||
adata.obs_names = bc
|
||
adata.var["gene_ids"] = feat[0].values
|
||
adata.var_names = feat[1].astype(str).values
|
||
adata.var_names_make_unique()
|
||
return adata
|
||
|
||
|
||
summary = []
|
||
fig, axes = plt.subplots(1, 3, figsize=(11, 3.2), sharey=False)
|
||
for ax, s in zip(axes, SAMPLES):
|
||
adata = load_raw(s)
|
||
umi = np.asarray(adata.X.sum(axis=1)).ravel()
|
||
n_genes = np.asarray((adata.X > 0).sum(axis=1)).ravel()
|
||
|
||
# knee plot：UMI 降序，log-log
|
||
umi_sorted = np.sort(umi)[::-1]
|
||
ax.plot(np.arange(1, len(umi_sorted) + 1), umi_sorted, lw=0.7)
|
||
ax.set_xscale("log"); ax.set_yscale("log")
|
||
ax.set_xlabel("Barcode rank"); ax.set_ylabel("UMI counts")
|
||
ax.set_title(s)
|
||
for thr, c in [(100, "grey"), (500, "orange"), (1000, "red")]:
|
||
n_pass = int((umi >= thr).sum())
|
||
ax.axhline(thr, ls="--", lw=0.6, color=c)
|
||
ax.text(5, thr * 1.15, f"≥{thr}: {n_pass:,}", fontsize=7, color=c)
|
||
|
||
for thr in (100, 500, 1000):
|
||
m = umi >= thr
|
||
summary.append({
|
||
"sample": s, "umi_threshold": thr,
|
||
"n_barcodes": int(m.sum()),
|
||
"median_umi": float(np.median(umi[m])) if m.any() else np.nan,
|
||
"median_genes": float(np.median(n_genes[m])) if m.any() else np.nan,
|
||
"total_genes_detected": int((np.asarray(adata.X[m].sum(axis=0)).ravel() > 0).sum()) if m.any() else 0,
|
||
})
|
||
print(f"[{s}] raw barcodes={adata.n_obs:,}, genes={adata.n_vars:,}, nnz={adata.X.nnz:,}")
|
||
|
||
fig.suptitle("Knee plots — GSE292269 raw matrices")
|
||
fig.tight_layout()
|
||
fig.savefig(OUT / "knee_plots.png")
|
||
plt.close(fig)
|
||
|
||
df = pd.DataFrame(summary)
|
||
df.to_csv(OUT / "barcode_threshold_summary.csv", index=False)
|
||
print("\n", df.to_string(index=False))
|