Scanpy 多样本分析完全指南：从数据读取到整合报告

多样本单细胞分析是科研项目里最常见的需求：多组对照、多时间点、多患者……单个样本分析还好说，多样本整合才是真正考验分析能力的地方。

---

一、项目结构规划

project/
├── data/
│   ├── sample1/filtered_feature_bc_matrix.h5
│   ├── sample2/filtered_feature_bc_matrix.h5
│   └── sample3/filtered_feature_bc_matrix.h5
├── metadata/
│   └── sample_info.csv   # sample_id, condition, patient_id, etc.
├── scripts/
│   └── analysis.py
└── results/
    ├── figures/
    └── tables/

---

二、批量读取和预处理

import scanpy as sc
import anndata as ad
import pandas as pd
import numpy as np
​
# 读取样本元数据
sample_meta = pd.read_csv("metadata/sample_info.csv")
print(sample_meta)
# sample_id  condition  patient_id  batch
# sample1    control    P01         batch1
# sample2    treatment  P01         batch1
# sample3    control    P02         batch2
​
# 批量读取并添加元数据
adatas = {}
for _, row in sample_meta.iterrows():
    sid = row["sample_id"]
    h5_path = f"data/{sid}/filtered_feature_bc_matrix.h5"
    tmp = sc.read_10x_h5(h5_path)
    tmp.var_names_make_unique()
    
    # 添加样本元数据
    for col in sample_meta.columns:
        tmp.obs[col] = row[col]
    
    # 确保 barcode 全局唯一
    tmp.obs_names = [f"{sid}_{bc}" for bc in tmp.obs_names]
    
    adatas[sid] = tmp
    print(f"读取 {sid}: {tmp.shape}")
​
# 合并所有样本
adata = ad.concat(adatas.values(), join="outer", fill_value=0)
print(f"\n合并后: {adata.shape}")

---

三、统一 QC 和过滤

# 标记线粒体基因（人类）
adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True)
​
# 按样本查看 QC 分布
import matplotlib.pyplot as plt
fig, axes = plt.subplots(3, 1, figsize=(12, 10))
for i, qc_var in enumerate(["n_genes_by_counts", "total_counts", "pct_counts_mt"]):
    import seaborn as sns
    sns.violinplot(data=adata.obs, x="sample_id", y=qc_var, ax=axes[i])
    axes[i].set_title(qc_var)
plt.tight_layout()
plt.savefig("results/figures/qc_per_sample.pdf")
​
# 过滤（基于各样本分布，可以设置不同阈值）
# 方法一：统一阈值
mask_cells = (
    (adata.obs["n_genes_by_counts"] >= 500) &
    (adata.obs["n_genes_by_counts"] <= 6000) &
    (adata.obs["total_counts"] >= 1000) &
    (adata.obs["pct_counts_mt"] <= 20)
)
​
# 方法二：MAD-based 自适应阈值（每个样本单独计算）
def mad_filter(values, n_mads=3):
    median = np.median(values)
    mad = np.median(np.abs(values - median))
    return (values >= median - n_mads * mad) & (values <= median + n_mads * mad)
​
sample_masks = []
for sid in sample_meta["sample_id"]:
    sample_idx = adata.obs["sample_id"] == sid
    sample_adata = adata[sample_idx]
    mask = (
        mad_filter(sample_adata.obs["n_genes_by_counts"].values) &
        mad_filter(np.log1p(sample_adata.obs["total_counts"].values)) &
        (sample_adata.obs["pct_counts_mt"] <= 20)
    )
    sample_masks.append(pd.Series(mask, index=sample_adata.obs_names))
​
adaptive_mask = pd.concat(sample_masks)[adata.obs_names].values
adata = adata[adaptive_mask].copy()
print(f"过滤后: {adata.shape}")

---

四、批次校正整合

# 保存原始 counts
adata.layers["counts"] = adata.X.copy()
​
# 标准化
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
​
# 高变基因（基于每个样本独立计算，取交集）
sc.pp.highly_variable_genes(
    adata,
    n_top_genes=3000,
    batch_key="sample_id",  # 关键参数！按样本独立计算
    flavor="seurat_v3",
    layer="counts"
)
print(f"高变基因数: {adata.var['highly_variable'].sum()}")
​
# PCA
sc.pp.scale(adata, max_value=10)
sc.pp.pca(adata, n_comps=50, use_highly_variable=True)
​
# Harmony 批次校正
import harmonypy as hm
ho = hm.run_harmony(
    adata.obsm["X_pca"],
    adata.obs,
    vars_use=["batch"],       # 按技术批次校正
    max_iter_harmony=50
)
adata.obsm["X_harmony"] = ho.Z_corr.T
print("Harmony 校正完成")

---

五、聚类和可视化

# 用校正后的表征构建邻居图
sc.pp.neighbors(adata, use_rep="X_harmony", n_neighbors=15)
sc.tl.umap(adata, random_state=42)
sc.tl.leiden(adata, resolution=0.5, random_state=42)
​
# 关键检查图
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
color_keys = ["leiden", "sample_id", "condition", "patient_id", 
              "total_counts", "pct_counts_mt"]
for ax, key in zip(axes.flatten(), color_keys):
    sc.pl.umap(adata, color=key, ax=ax, show=False, 
               frameon=False, title=key)
plt.tight_layout()
plt.savefig("results/figures/umap_overview.pdf", bbox_inches="tight")

---

六、差异分析：跨条件比较

# 伪批次差异表达（Pseudobulk，比直接做更可靠）
import pandas as pd
from scipy.sparse import issparse
​
def pseudobulk(adata, groupby_cols, layer="counts"):
    """将单细胞合并成伪批次，用于差异分析"""
    X = adata.layers[layer] if layer in adata.layers else adata.X
    if issparse(X):
        X = pd.DataFrame.sparse.from_spmatrix(
            X, index=adata.obs_names, columns=adata.var_names
        )
    else:
        X = pd.DataFrame(X, index=adata.obs_names, columns=adata.var_names)
    
    meta = adata.obs[groupby_cols]
    pb = X.join(meta).groupby(groupby_cols).sum()
    return pb
​
# 为每种细胞类型分别做伪批次差异分析
for cell_type in adata.obs["cell_type"].unique():
    ct_adata = adata[adata.obs["cell_type"] == cell_type]
    if ct_adata.n_obs < 50:  # 细胞数太少，跳过
        continue
    
    # 伪批次矩阵
    pb = pseudobulk(ct_adata, ["patient_id", "condition"])
    print(f"{cell_type}: {pb.shape}")
    # 后续可用 pydeseq2 做差异分析

---

七、输出标准化报告

# 保存完整分析结果
adata.write_h5ad("results/integrated_analysis.h5ad", compression="gzip")
​
# 导出细胞类型组成表
cell_composition = (
    adata.obs.groupby(["sample_id", "condition", "cell_type"])
    .size()
    .reset_index(name="count")
)
cell_composition.to_csv("results/tables/cell_composition.csv", index=False)
​
# 导出 UMAP 坐标
umap_df = pd.DataFrame(
    adata.obsm["X_umap"],
    index=adata.obs_names,
    columns=["UMAP1", "UMAP2"]
).join(adata.obs[["sample_id", "condition", "leiden", "cell_type"]])
umap_df.to_csv("results/tables/umap_coordinates.csv")
​
print("分析完成！结果已保存到 results/ 目录")

---

多样本分析的核心是：一致的预处理 + 合理的批次校正 + 严格的质量验证。每一步都要可视化检查，而不是"跑通了就行"。

Run2AI 运智（https://run2ai.open2ai.cn）专注多样本整合分析，支持 2-20 个样本的批量处理，提供完整的批次校正验证报告和标准化代码交付。