单细胞 RNA+ATAC 联合分析：从双模态对象到 WNN

2022 年 9 月，我留下过一份使用 Seurat 与 Signac 处理单细胞 RNA+ATAC 数据的脚本。它完成了从 10x 多组学矩阵、双模态质控，到 RNA、ATAC 降维和 WNN 聚类的主链路。

这次复盘聚焦联合分析中每一步为什么存在，以及哪些参数不能机械照搬。

工作流概览

原始流程可以概括为：

读取每个样本的 RNA 与 ATAC 计数
→ 创建 RNA 对象
→ 添加 ATAC ChromatinAssay
→ 分样本质控
→ 合并样本
→ RNA 降维
→ ATAC 降维
→ WNN 联合邻居图
→ UMAP 与聚类

RNA 和 ATAC 先各自建立低维表示，再通过加权最近邻整合，而不是在最开始把两种矩阵简单拼接。

1. 从同一批细胞建立双模态对象

RNA 计数进入 Seurat assay，ATAC peak 计数进入 ChromatinAssay。ATAC 对象还需要：

peak 的染色体坐标格式；
fragments 文件；
与参考基因组一致的注释；
统一的染色体命名风格。

原脚本使用小鼠 mm10 和 EnsDb.Mmusculus.v79。参考版本必须和建库、比对及 peak calling 阶段一致，否则基因注释和坐标可能错位。

2. 质控应按样本观察，而不是复制阈值

原始脚本同时根据 RNA counts、ATAC counts 和线粒体比例过滤细胞。这是合理的起点，但阈值来自当时的数据分布，不能直接用于其他项目。

更稳妥的顺序是：

分样本绘制 nCount_RNA、nFeature_RNA、percent.mt；
检查 nCount_ATAC、nFeature_ATAC；
增加 TSS enrichment、nucleosome signal 等 ATAC 质量指标；
观察不同样本分布后再确定阈值；
记录每一步过滤前后的细胞数。

如果先合并再使用同一阈值，可能会把样本质量差异误当作生物差异。

3. RNA 模态：归一化和表达空间

原流程使用：

DefaultAssay(object) <- "RNA"
object <- SCTransform(object) |>
  RunPCA() |>
  RunUMAP(
    dims = 1:50,
    reduction.name = "umap.rna",
    reduction.key = "rnaUMAP_"
  )

PCA 维数不应固定照搬。需要结合样本规模、方差结构、批次处理方式和下游稳定性选择。

4. ATAC 模态：TF-IDF 与 LSI

ATAC peak 矩阵更稀疏，原流程使用 TF-IDF、特征筛选和 SVD：

DefaultAssay(object) <- "ATAC"
object <- RunTFIDF(object)
object <- FindTopFeatures(object, min.cutoff = "q0")
object <- RunSVD(object)
object <- RunUMAP(
  object,
  reduction = "lsi",
  dims = 2:50,
  reduction.name = "umap.atac",
  reduction.key = "atacUMAP_"
)

这里从第二个 LSI 维度开始，是因为第一个维度经常与测序深度高度相关。但这仍需要通过相关性和可视化确认，而不是当作永恒规则。

5. WNN：让每个细胞选择更有信息的模态

WNN 使用 RNA 的 PCA 空间和 ATAC 的 LSI 空间共同构建邻居关系：

object <- FindMultiModalNeighbors(
  object,
  reduction.list = list("pca", "lsi"),
  dims.list = list(1:50, 2:50)
)

object <- RunUMAP(
  object,
  nn.name = "weighted.nn",
  reduction.name = "wnn.umap",
  reduction.key = "wnnUMAP_"
)

object <- FindClusters(
  object,
  graph.name = "wsnn",
  algorithm = 3
)

它不是简单地给 RNA 和 ATAC 固定相同权重，而是利用局部邻域信息估计不同模态对每个细胞的贡献。

6. 至少并排检查三套结果

联合分析完成后，我会同时查看：

RNA UMAP；
ATAC UMAP；
WNN UMAP。

并分别按样本、批次、实验组和聚类标签着色。若 WNN 的主要分离仍由批次驱动，或者某个模态几乎完全决定邻居关系，需要回到质控、降维或批次校正阶段检查。

原脚本需要改进的地方

这份历史脚本还不是可直接复用的软件流程，主要缺少：

Seurat、Signac 和参考注释版本；
输入文件 manifest；
fragments 与 barcode 一致性检查；
双细胞检测；
分样本质控图和阈值依据；
批次效应评估；
随机种子与 session 信息；
每个检查点对应的细胞数量。

此外，工作流更适合使用命名列表、配置文件和循环函数管理样本，避免路径与对象散落在代码中。代码复用前还应固定软件版本、参考注释、随机种子和输入 manifest。