748. 空间转录组细胞解卷积¶

10x Visium等平台：
  每个spot直径55μm → 通常包含1-10个细胞
  测到的基因表达 = 所有细胞表达的混合

问题：
  Spot A的表达谱 = 30% T细胞 + 50% 巨噬细胞 + 20% 成纤维细胞
  但我们只看到了混合后的总表达

解卷积就是：
  输入：每个spot的混合表达谱 + 参考scRNA-seq（知道各细胞类型的"特征表达"）
  输出：每个spot中各细胞类型的比例

# 安装cell2location
pip install cell2location  # 安装主包（需要PyTorch）

# 如果有GPU
pip install cell2location[gpu]

# 验证
python -c "import cell2location; print(cell2location.__version__)"

# ===== 准备输入数据 =====
import scanpy as sc  # 导入scanpy
import cell2location  # 导入cell2location
import numpy as np

# 1. 加载空间转录组数据
adata_vis = sc.read_visium("visium_output/")  # 读取10x Visium数据
adata_vis.var_names_make_unique()  # 确保基因名唯一

# 基本质控
sc.pp.filter_genes(adata_vis, min_cells=3)  # 至少3个spot表达的基因
adata_vis.var["mt"] = adata_vis.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata_vis, qc_vars=["mt"], inplace=True)

# 2. 加载参考scRNA-seq数据
adata_ref = sc.read_h5ad("reference_scrna.h5ad")  # 读取参考数据
# 确保adata_ref.obs["cell_type"]列有细胞类型注释

# 过滤参考数据
sc.pp.filter_genes(adata_ref, min_cells=10)  # 至少10个细胞表达

# 取交集基因
shared_genes = adata_vis.var_names.intersection(adata_ref.var_names)
adata_vis = adata_vis[:, shared_genes]
adata_ref = adata_ref[:, shared_genes]
print(f"共有基因数: {len(shared_genes)}")

# ===== 第一步：从scRNA-seq学习细胞类型特征 =====
from cell2location.models import RegressionModel

# 设置模型
RegressionModel.setup_anndata(
    adata_ref,
    batch_key="sample",  # 参考数据的批次（如果有）
    labels_key="cell_type"  # 细胞类型标签
)

# 创建回归模型
mod_ref = RegressionModel(adata_ref)

# 训练
mod_ref.train(
    max_epochs=250,  # 最大训练轮数
    batch_size=2500,  # 批次大小
    train_size=1,     # 使用全部数据训练
    lr=0.002          # 学习率
)

# 查看训练历史
mod_ref.plot_history(20)  # 查看最后20个epoch的loss

# 导出每种细胞类型的特征表达谱
adata_ref = mod_ref.export_posterior(
    adata_ref,
    sample_kwargs={"num_samples": 1000, "batch_size": 2500}
)

# 保存参考特征
ref_sig = adata_ref.varm["means_per_cluster_mu_fg"]  # 细胞类型特征矩阵

# ===== 第二步：在空间数据上做解卷积 =====
from cell2location.models import Cell2location

# 准备空间数据
# 将参考特征添加到空间数据
inf_aver = adata_ref.varm["means_per_cluster_mu_fg"].copy()
adata_vis.uns["mod"] = {"inf_aver": inf_aver}

# 设置Cell2location模型
Cell2location.setup_anndata(
    adata_vis,
    batch_key=None  # Visium数据通常不分批次
)

# 创建模型
mod_vis = Cell2location(
    adata_vis,
    cell_state_df=inf_aver,  # 参考细胞类型特征
    N_cells_per_location=10,  # 每个spot预期的细胞数
    detection_alpha=200       # 检测灵敏度参数
)

# 训练
mod_vis.train(
    max_epochs=30000,  # 需要更多轮次
    batch_size=None,   # 使用全部数据
    train_size=1,
    lr=0.002
)

# 导出结果
adata_vis = mod_vis.export_posterior(
    adata_vis,
    sample_kwargs={"num_samples": 1000, "batch_size": mod_vis.adata.n_obs}
)

# 每种细胞类型在每个spot中的估计细胞数
# 存储在 adata_vis.obsm["q05_cell_abundance_w_sf"]
cell_counts = adata_vis.obsm["q05_cell_abundance_w_sf"]
print(f"Cell type columns: {cell_counts.columns.tolist()}")

# ===== 可视化解卷积结果 =====
import matplotlib.pyplot as plt

# 在组织切片上显示各细胞类型的空间分布
cell_types = cell_counts.columns.tolist()

fig, axes = plt.subplots(2, 4, figsize=(20, 10))
axes = axes.flatten()

for i, ct in enumerate(cell_types[:8]):  # 显示前8种细胞类型
    adata_vis.obs[ct] = cell_counts[ct].values
    sc.pl.spatial(
        adata_vis,
        color=ct,
        cmap="Reds",
        size=1.3,
        ax=axes[i],
        show=False,
        title=ct
    )

plt.tight_layout()
plt.savefig("cell_type_spatial_distribution.png", dpi=300)
plt.show()

# 饼图可视化（每个spot一个饼图显示细胞类型比例）
# 转为比例
cell_proportions = cell_counts.div(cell_counts.sum(axis=1), axis=0)
adata_vis.obsm["cell_proportions"] = cell_proportions

# ===== RCTD解卷积 =====
library(spacexr)  # 安装：devtools::install_github("dmcable/spacexr")
library(Seurat)

# 准备参考数据
reference <- Reference(
  counts = ref_counts,  # 参考scRNA-seq的count矩阵
  cell_types = ref_cell_types,  # 细胞类型向量
  nUMI = ref_nUMI  # 每个细胞的UMI总数
)

# 准备空间数据
coords <- data.frame(
  x = spatial_coords$x,
  y = spatial_coords$y
)
puck <- SpatialRNA(
  coords = coords,  # 空间坐标
  counts = spatial_counts,  # 空间count矩阵
  nUMI = colSums(spatial_counts)  # 每个spot的UMI数
)

# 创建RCTD对象
myRCTD <- create.RCTD(
  puck,
  reference,
  max_cores = 8  # CPU核心数
)

# 运行解卷积
myRCTD <- run_RCTD(
  myRCTD,
  doublet_mode = "full"  # full=允许每个spot多种细胞类型
)

# 提取结果
results <- myRCTD@results
weights <- results$weights  # 每个spot中各细胞类型的比例

# ===== 空间解卷积速查 =====

# Cell2location (Python)
pip install cell2location
# Step 1: 训练参考模型
mod_ref = RegressionModel(adata_ref)
mod_ref.train(max_epochs=250)
# Step 2: 空间解卷积
mod_vis = Cell2location(adata_vis, cell_state_df=inf_aver)
mod_vis.train(max_epochs=30000)
# 结果: adata_vis.obsm["q05_cell_abundance_w_sf"]

# RCTD (R)
library(spacexr)
reference <- Reference(counts, cell_types, nUMI)
puck <- SpatialRNA(coords, counts, nUMI)
myRCTD <- create.RCTD(puck, reference)
myRCTD <- run_RCTD(myRCTD, doublet_mode="full")
# 结果: myRCTD@results$weights

# 参考数据要求：
# 同组织同物种 + 每种细胞类型≥25个 + 注释准确