单细胞多组学（scATAC-seq / CITE-seq / Multiome）¶

scRNA-seq = 只看你的"购物清单"（mRNA），猜你要做什么菜
scATAC-seq = 看你厨房里"哪些柜子是打开的"（染色质开放区域），推测你可能要用哪些食材
CITE-seq  = 看你身上"穿了什么工作服"（表面蛋白），判断你是厨师还是服务员
Multiome  = 同时看你的购物清单 + 哪些柜子开着（RNA + ATAC 联合测）

DNA 在细胞核里不是散乱的，而是缠绕在"组蛋白"上形成"染色质"。
- 染色质"松开"的区域 → 基因可以被转录 → 叫做"开放区域"(open chromatin)
- 染色质"紧缩"的区域 → 基因被锁住了 → 叫做"关闭区域"

scATAC-seq 的原理：
用一种叫 Tn5 的酶（转座酶），它只能切割"松开"的 DNA 区域。
切完之后测序，看哪些位置被切了 → 就知道哪些区域是开放的。

白话：想象一栋大楼有 1000 扇门（基因），scATAC-seq 就是派一个人去
      试着推每扇门，能推开的就记录下来 → 你就知道哪些门是"开着的"。

# ====== 安装 ArchR ======
# 需要 R >= 4.0
if (!requireNamespace("devtools", quietly = TRUE)) install.packages("devtools")
devtools::install_github("GreenleafLab/ArchR", ref="master", repos = BiocManager::repositories())

library(ArchR)          # 加载 ArchR 包
set.seed(1)             # 设置随机种子，保证可重复
addArchRThreads(threads = 8)  # 设置并行线程数（根据你的电脑 CPU 核心数调整）
addArchRGenome("hg38")  # 设置参考基因组（人类用 hg38，小鼠用 mm10）

# ====== Step 1: 创建 Arrow 文件（ArchR 的核心数据格式）======
# fragments.tsv.gz 是 10x CellRanger-ATAC 的输出文件
ArrowFiles <- createArrowFiles(
  inputFiles = "fragments.tsv.gz",   # 输入：片段文件
  sampleNames = "PBMC",              # 样本名
  filterTSS = 4,                     # 质控：TSS 富集分数 ≥ 4（太低说明质量差）
  filterFrags = 1000,                # 质控：至少 1000 个 unique fragments
  addTileMat = TRUE,                 # 创建 500bp bin 的 Tile 矩阵
  addGeneScoreMat = TRUE             # 创建基因活性评分矩阵（用来推断基因表达）
)

# ====== Step 2: 推断 doublet 并创建 ArchRProject ======
doubScores <- addDoubletScores(input = ArrowFiles, k = 10, knnMethod = "UMAP")
proj <- ArchRProject(ArrowFiles = ArrowFiles, outputDirectory = "output")
proj <- filterDoublets(ArchRProj = proj)  # 过滤掉 doublets（一个液滴里包了两个细胞）

# ====== Step 3: 降维（Iterative LSI，而不是 PCA）======
# LSI = Latent Semantic Indexing，借鉴自文本检索领域
# scATAC 数据太稀疏，PCA 效果不好，LSI 更适合
proj <- addIterativeLSI(ArchRProj = proj, useMatrix = "TileMatrix", name = "IterativeLSI")

# ====== Step 4: 聚类和 UMAP 可视化 ======
proj <- addClusters(input = proj, reducedDims = "IterativeLSI")
proj <- addUMAP(ArchRProj = proj, reducedDims = "IterativeLSI")

# ====== Step 5: 用 Gene Score 推断细胞类型 ======
# Gene Score = 根据基因附近的 ATAC peaks 信号推断该基因的"活性"
proj <- addImputeWeights(proj)  # 用 MAGIC 平滑 dropout 噪声
plotEmbedding(proj, colorBy = "GeneScoreMatrix", name = "CD14",
              embedding = "UMAP", imputeWeights = getImputeWeights(proj))

# ====== Step 6: Peak Calling（识别开放区域）======
proj <- addGroupCoverages(ArchRProj = proj, groupBy = "Clusters")
proj <- addReproduciblePeakSet(ArchRProj = proj, groupBy = "Clusters")
proj <- addPeakMatrix(proj)  # 创建 Peak × Cell 矩阵

# ====== Step 7: Motif 分析（找转录因子结合位点）======
proj <- addMotifAnnotations(ArchRProj = proj, motifSet = "cisbp", name = "Motif")
markerMotifs <- getMarkerFeatures(
  ArchRProj = proj, useMatrix = "MotifMatrix",
  groupBy = "Clusters", testMethod = "wilcoxon"
)

library(Signac)         # 加载 Signac
library(Seurat)         # 加载 Seurat（Signac 依赖 Seurat）
library(EnsDb.Hsapiens.v86)  # 基因注释（hg38 用这个）

# ====== Step 1: 创建 ChromatinAssay ======
counts <- Read10X_h5("filtered_peak_bc_matrix.h5")   # 读取 10x 输出
fragpath <- "fragments.tsv.gz"                        # 片段文件路径

# 获取基因注释信息
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotations) <- "UCSC"  # 统一染色体命名格式

# 创建 ChromatinAssay（Signac 特有的 assay 类型）
chrom_assay <- CreateChromatinAssay(
  counts = counts, sep = c(":", "-"),
  fragments = fragpath, annotation = annotations
)
pbmc <- CreateSeuratObject(counts = chrom_assay, assay = "peaks")

# ====== Step 2: 质控 ======
pbmc <- NucleosomeSignal(pbmc)   # 计算核小体信号（应 < 4）
pbmc <- TSSEnrichment(pbmc)      # 计算 TSS 富集分数（应 > 2）
pbmc <- subset(pbmc, nCount_peaks > 3000 & nCount_peaks < 30000 &
               nucleosome_signal < 4 & TSS.enrichment > 2)

# ====== Step 3: 降维（TF-IDF + SVD = LSI）======
pbmc <- RunTFIDF(pbmc)           # TF-IDF 归一化（替代 scRNA 的 NormalizeData）
pbmc <- FindTopFeatures(pbmc, min.cutoff = "q0")  # 选择特征
pbmc <- RunSVD(pbmc)             # SVD 降维（替代 PCA）

# ====== Step 4: 聚类和可视化 ======
pbmc <- RunUMAP(pbmc, reduction = "lsi", dims = 2:30)  # 注意跳过第 1 维（通常与测序深度相关）
pbmc <- FindNeighbors(pbmc, reduction = "lsi", dims = 2:30)
pbmc <- FindClusters(pbmc, algorithm = 3, resolution = 0.5)

# 用 Signac/Seurat 做跨模态整合
# 核心思想：用 Gene Activity Score（scATAC 的基因活性推断值）
#           与 scRNA-seq 的基因表达进行锚点匹配（anchor-based integration）

gene.activities <- GeneActivity(atac_obj)  # 从 ATAC peaks 推断基因活性
atac_obj[["RNA"]] <- CreateAssayObject(counts = gene.activities)
atac_obj <- NormalizeData(atac_obj, assay = "RNA")

# 找到两个数据集之间的锚点
transfer.anchors <- FindTransferAnchors(
  reference = rna_obj,         # scRNA-seq 作为参考
  query = atac_obj,            # scATAC-seq 作为查询
  reduction = "cca"            # 用 CCA 做跨模态匹配
)

# 把 scRNA-seq 的细胞类型标签转移到 scATAC-seq 上
predicted.labels <- TransferData(
  anchorset = transfer.anchors,
  refdata = rna_obj$celltype,  # 参考数据的细胞类型
  weight.reduction = atac_obj[["lsi"]]
)
atac_obj$predicted.celltype <- predicted.labels$predicted.id

CITE-seq (Cellular Indexing of Transcriptomes and Epitopes by Sequencing)

核心思想：
在做 scRNA-seq 的同时，给细胞表面的蛋白贴上"DNA 条形码抗体"。
测序时同时读到 mRNA 和这些"蛋白标签"的条形码。

白话类比：
传统 scRNA-seq = 给每个学生考语文
CITE-seq      = 语文考试的同时，还用扫码枪扫了他们身上的"学生证"
                （学生证上写着班级、年级等信息 = 表面蛋白类型）

具体流程：
1. 准备带 DNA 标签的抗体（antibody-derived tag, ADT）
   每种抗体识别一种表面蛋白（如 CD4, CD8, CD19 等）
2. 把抗体和细胞孵育 → 抗体结合到细胞表面对应蛋白上
3. 做常规 10x scRNA-seq → 液滴里同时包含 mRNA 和抗体的 DNA 标签
4. 测序 → 同时得到基因表达矩阵 (RNA) 和蛋白丰度矩阵 (ADT)

library(Seurat)

# ====== 加载 CITE-seq 数据（以 Seurat bmcite 数据集为例）======
# bmcite 数据包含 30672 个骨髓细胞，25 种表面蛋白
library(SeuratData)
InstallData("bmcite")     # 安装示例数据
bm <- LoadData(ds = "bmcite")  # 加载数据

# 数据结构：一个 Seurat 对象里有两个 assay
# bm[["RNA"]]  → 基因表达矩阵（~2 万基因）
# bm[["ADT"]]  → 蛋白丰度矩阵（25 种蛋白）

# ====== RNA 的标准处理 ======
DefaultAssay(bm) <- "RNA"
bm <- NormalizeData(bm) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA()

# ====== ADT 的特殊处理 ======
DefaultAssay(bm) <- "ADT"
VariableFeatures(bm) <- rownames(bm[["ADT"]])  # ADT 通道蛋白数少，全部用作特征
bm <- NormalizeData(bm, normalization.method = "CLR", margin = 2) %>%
      # CLR = Centered Log-Ratio 归一化（适合成分数据，和 RNA 不同！）
      # margin = 2 表示按特征（蛋白）方向归一化
      ScaleData() %>%
      RunPCA(reduction.name = "apca")  # 单独做 PCA，结果命名为 "apca"

WNN 的核心思想：
对于每个细胞，根据 RNA 和蛋白（或 ATAC）的信息含量，
自动学习每种模态的"权重"，然后用加权组合构建细胞间的邻居图。

白话类比：
假设你要给学生分班：
- 有些学生的"语文成绩"更能区分他们（比如文科生），那语文权重高
- 有些学生的"数学成绩"更能区分他们（比如理科生），那数学权重高
- WNN 就是根据每个学生的具体情况，动态调整两科的权重

实际效果：
- 祖细胞（progenitor cells）：RNA 权重高（因为抗体面板没有好的祖细胞标记）
- T 细胞亚群：ADT 权重高（因为 CD4/CD8 蛋白标记比 mRNA 更准确）

# ====== WNN 分析（Seurat v5 实现）======
# FindMultiModalNeighbors 是核心函数
bm <- FindMultiModalNeighbors(
  bm,
  reduction.list = list("pca", "apca"),  # RNA 用 "pca"，ADT 用 "apca"
  dims.list = list(1:30, 1:18),          # RNA 用前 30 个 PC，ADT 用前 18 个
  modality.weight.name = "RNA.weight"    # 权重存储在这个 metadata 列
)
# 该函数会生成：
# bm[["weighted.nn"]] → 加权最近邻图
# bm[["wknn"]]        → WNN 的 KNN 图
# bm[["wsnn"]]        → WNN 的 SNN 图（用于聚类）
# bm$RNA.weight       → 每个细胞的 RNA 权重（0~1，越高说明 RNA 信息越重要）

# ====== 基于 WNN 图做 UMAP 和聚类 ======
bm <- RunUMAP(bm, nn.name = "weighted.nn", reduction.name = "wnn.umap",
              reduction.key = "wnnUMAP_")  # WNN UMAP
bm <- FindClusters(bm, graph.name = "wsnn", algorithm = 3, resolution = 2)

# ====== 对比三种 UMAP ======
# 1) 只用 RNA 的 UMAP
bm <- RunUMAP(bm, reduction = "pca", dims = 1:30,
              reduction.name = "rna.umap", reduction.key = "rnaUMAP_")
# 2) 只用 ADT 的 UMAP
bm <- RunUMAP(bm, reduction = "apca", dims = 1:18,
              reduction.name = "adt.umap", reduction.key = "adtUMAP_")
# 3) WNN UMAP（已在上面生成）

# 对比结果：WNN UMAP 通常能更清晰地分开不同细胞亚群
# RNA UMAP 在区分祖细胞方面更好
# ADT UMAP 在区分 T 细胞亚群方面更好
# WNN UMAP 结合了两者的优势

10x Genomics Multiome ATAC + Gene Expression：
从同一个细胞核中同时获取 RNA 表达和 ATAC 信号。

关键区别：
- 分别做 scRNA-seq 和 scATAC-seq → 来自不同细胞，需要"计算整合"
- Multiome → 来自同一个细胞核，天然配对，不需要计算锚点匹配

白话：
分别做 = 两个班各考一科（语文班考语文，数学班考数学），
         事后按"相似"配对 → 可能配错
Multiome = 同一个班同时考语文和数学 → 每个学生两科成绩天然对应

library(Seurat)
library(Signac)
library(EnsDb.Hsapiens.v86)

# ====== Step 1: 读取 10x Multiome 数据 ======
# 10x Multiome 的 filtered_feature_bc_matrix.h5 里同时包含 RNA 和 ATAC 数据
inputdata.10x <- Read10X_h5("filtered_feature_bc_matrix.h5")

rna_counts <- inputdata.10x$`Gene Expression`  # 提取 RNA 计数矩阵
atac_counts <- inputdata.10x$Peaks              # 提取 ATAC peaks 矩阵

# ====== Step 2: 创建多模态 Seurat 对象 ======
pbmc <- CreateSeuratObject(counts = rna_counts)
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

# 添加 ATAC assay
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotations) <- "UCSC"
chrom_assay <- CreateChromatinAssay(
  counts = atac_counts, sep = c(":", "-"),
  genome = "hg38", fragments = "atac_fragments.tsv.gz",
  min.cells = 10, annotation = annotations
)
pbmc[["ATAC"]] <- chrom_assay  # 一个对象里同时有 RNA 和 ATAC

# ====== Step 3: 质控 ======
pbmc <- subset(pbmc,
  nCount_ATAC < 7e4 & nCount_ATAC > 5e3 &
  nCount_RNA < 25000 & nCount_RNA > 1000 &
  percent.mt < 20
)

# ====== Step 4: 分别预处理两个模态 ======
# RNA
DefaultAssay(pbmc) <- "RNA"
pbmc <- SCTransform(pbmc) %>% RunPCA() %>%
        RunUMAP(dims = 1:50, reduction.name = "umap.rna")

# ATAC
DefaultAssay(pbmc) <- "ATAC"
pbmc <- RunTFIDF(pbmc) %>% FindTopFeatures(min.cutoff = "q0") %>%
        RunSVD() %>% RunUMAP(reduction = "lsi", dims = 2:50,
                             reduction.name = "umap.atac")

# ====== Step 5: WNN 联合分析 ======
pbmc <- FindMultiModalNeighbors(pbmc,
  reduction.list = list("pca", "lsi"),   # RNA 用 PCA，ATAC 用 LSI
  dims.list = list(1:50, 2:50)           # ATAC 跳过第 1 维
)
pbmc <- RunUMAP(pbmc, nn.name = "weighted.nn",
                reduction.name = "wnn.umap")
pbmc <- FindClusters(pbmc, graph.name = "wsnn", algorithm = 3)

# ====== Step 6: 联合找转录因子调控子 ======
# 关键应用：同时看 TF 的 mRNA 表达 + motif 可及性
# 如果两者都在某个细胞类型中富集 → 该 TF 很可能是该细胞类型的关键调控因子

MOFA+ 做了什么：
把多种组学数据（可以是 RNA + ATAC + 蛋白 + 甲基化...任意组合）
分解成少数几个"潜在因子"(latent factors)。

白话：
假设你有一群学生，考了语文、数学、英语、物理四科。
MOFA+ 会找到几个"隐藏因子"：
- 因子 1 = "理科能力"（和数学、物理成绩高度相关）
- 因子 2 = "文科能力"（和语文、英语成绩高度相关）
- 因子 3 = "整体学习态度"（和所有科目都弱相关）

每个因子同时解释多种组学的变异 → 找到跨模态的共同信号。

# ====== MOFA+ 的 Python 接口（使用 mofapy2）======
# 安装：pip install mofapy2
from mofapy2.run.entry_point import entry_point
import numpy as np

# 准备数据：list of numpy arrays，每个元素是一个组学的矩阵
# 行 = 样本/细胞，列 = 特征
data = [
    [rna_matrix],    # 视图 1: RNA 数据 (cells × genes)
    [atac_matrix]    # 视图 2: ATAC 数据 (cells × peaks)
]

# 初始化 MOFA 模型
ent = entry_point()
ent.set_data_options(scale_groups=False, scale_views=True)  # 按视图标准化
ent.set_data_matrix(data, groups_names=["group1"], views_names=["RNA", "ATAC"])
ent.set_model_options(factors=10)        # 设置因子数为 10
ent.set_training_options(iter=1000, seed=42)
ent.build()
ent.run()
ent.save("mofa_model.hdf5")  # 保存模型

# ====== MultiVI 分析流程 ======
# 安装：pip install scvi-tools
import scvi
import anndata as ad

# 输入数据：
# - adata_multi: 配对的 Multiome 数据（同一细胞有 RNA + ATAC）
# - adata_rna:   只有 RNA 的数据
# - adata_atac:  只有 ATAC 的数据

# 合并所有数据
adata = ad.concat([adata_multi, adata_rna, adata_atac])

# 设置 MultiVI 模型
scvi.model.MULTIVI.setup_anndata(adata, batch_key="modality")

model = scvi.model.MULTIVI(adata, n_genes=2000, n_regions=2000)
model.train(max_epochs=500)  # 训练 VAE 模型

# 获取联合潜在空间
latent = model.get_latent_representation()  # 所有细胞的联合 embedding
adata.obsm["X_multivi"] = latent

# 在联合空间中做下游分析
import scanpy as sc
sc.pp.neighbors(adata, use_rep="X_multivi")  # 基于联合 embedding 找邻居
sc.tl.umap(adata)                             # UMAP 可视化
sc.tl.leiden(adata)                            # 聚类

import muon as mu

# MuData = 多模态 AnnData 的容器
# 每个模态是一个 AnnData，共享同一批细胞
mdata = mu.read_10x_h5("filtered_feature_bc_matrix.h5")  # 读取 10x Multiome 数据
# mdata.mod["rna"]  → RNA 的 AnnData
# mdata.mod["atac"] → ATAC 的 AnnData

# RNA 预处理
mu.pp.intersect_obs(mdata)  # 保留两个模态都有的细胞
sc.pp.normalize_total(mdata.mod["rna"])
sc.pp.log1p(mdata.mod["rna"])

# ATAC 预处理
mu.atac.pp.tfidf(mdata.mod["atac"])
mu.atac.tl.lsi(mdata.mod["atac"])

# WNN 分析（muon 也实现了 WNN）
mu.pp.neighbors(mdata, key_added="wnn")
mu.tl.umap(mdata, neighbors_key="wnn")

你想测什么？
├── 只测染色质开放性 → scATAC-seq（ArchR 或 Signac）
├── 基因表达 + 表面蛋白 → CITE-seq（Seurat WNN）
├── 基因表达 + 染色质开放性（同一细胞）→ 10x Multiome（Seurat+Signac WNN）
├── 基因表达 + 染色质开放性（不同细胞）→ 计算整合（Signac transfer / MultiVI）
└── 三模态或更多 → TEA-seq / MOFA+ / GLUE

你的数据类型？
├── CITE-seq (RNA + ADT)
│   ├── Seurat WNN（R，首选）
│   ├── totalVI（Python, scvi-tools）
│   └── muon WNN（Python）
├── 10x Multiome (RNA + ATAC，配对)
│   ├── Seurat + Signac WNN（R，首选）
│   ├── ArchR + scRNA 整合（R）
│   ├── MultiVI（Python）
│   └── muon（Python）
├── 非配对 scRNA + scATAC
│   ├── Signac 的 TransferData（R）
│   ├── MultiVI（Python, 推荐）
│   └── GLUE（Python）
└── 多个组学（≥3 种）
    ├── MOFA+（R/Python）
    └── GLUE（Python）