ATAC + RNA 多组学：基于 ArchR 的基础分析

作者: SeekGene

时长: 28 分钟

字数: 7.0k 字

更新: 2026-07-07

阅读: 0 次

1. 教程简介

本教程基于 ArchR 分析框架，对 SeekArc 单细胞多组学数据（RNA + ATAC） 进行 数据导入与预处理、多样本联合分析、降维聚类、细胞类型注释。

该流程面向同一细胞同时获得 转录组表达 与 染色质开放性 的 SeekArc 数据：以 ATAC fragments 构建 Arrow 项目，将 seekARC 基因表达矩阵挂载到 ArchRProject，在多样本场景下完成从注释准备、质控过滤到联合分析与调控机制解析的完整流程，实现以下核心分析目标：

多样本多组学联合分析：在同一样本或跨多个 SeekArc 样本中，同时利用 RNA 表达与 ATAC 开放染色质信息，更全面地刻画细胞状态与生物学异质性。
数据导入与预处理：基于 GTF 构建基因坐标注释；从 filtered_feature_bc_matrix 读取 RNA 并整合为 SummarizedExperiment；由 fragments 文件创建 Arrow（含 Tile Matrix 与 Gene Score Matrix）；初始化 ArchRProject 并完成 RNA/ATAC 条形码匹配与基因表达矩阵挂载，为后续分析提供统一的项目对象。
数据质控与过滤：依据 TSS 富集、fragment 数量等指标对 ATAC 细胞进行质控，并结合 RNA 挂载时的严格匹配策略，保留双模态一致的可靠细胞。
单模态与多模态聚类分析：分别基于 peak 开放（ATAC）与 Gene Expression（RNA）开展降维聚类，并通过 ArchR 的联合降维策略整合两种模态，用于比较不同模态下的细胞分群结果。
联合降维与可视化：融合 RNA 与 ATAC 信息构建更稳定的低维表示，并通过 UMAP 展示细胞分布与样本间整合效果。
细胞类型注释：结合 marker 基因的表达与基因组开放模式，对聚类结果进行生物学解释和细胞类型标注。

library(ArchR)
library(rtracklayer)
library(GenomicFeatures)
library(dplyr)
library(DropletUtils)
library(SummarizedExperiment)

2. 输入文件准备

输入文件要求

请保证每个样本均按照统一的目录结构进行整理。每个样本使用一个独立文件夹，文件夹名称建议直接使用样本 ID，例如 pbmc_1、pbmc_2 等。每个样本目录下应包含以下内容：

filtered_feature_bc_matrix：scRNA-seq 表达矩阵目录，包含 barcodes.tsv.gz、features.tsv.gz 和 matrix.mtx.gz 文件。
filtered_peaks_bc_matrix：scATAC-seq peak 开放矩阵目录，包含 barcodes.tsv.gz、features.tsv.gz 和 matrix.mtx.gz 文件。
{样本ID}_A_fragments.tsv.gz：ATAC 片段文件，用于后续构建 ChromatinAssay、计算 TSS 富集和核小体信号等分析。
{样本ID}_A_fragments.tsv.gz.tbi：ATAC 片段文件对应的索引文件。
genes.gtf：参考基因组注释文件，与比对所用基因组版本一致。

目录结构示例

text

pbmc_1/
├── filtered_feature_bc_matrix/
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz
├── filtered_peaks_bc_matrix/
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz
├── pbmc_1_A_fragments.tsv.gz
└── pbmc_1_A_fragments.tsv.gz.tbi

pbmc_2/
├── filtered_feature_bc_matrix/
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz
├── filtered_peaks_bc_matrix/
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz
├── pbmc_2_A_fragments.tsv.gz
└── pbmc_2_A_fragments.tsv.gz.tbi

addArchRGenome("hg38") #按照自己物种信息选择
samples = c('pbmc_1', 'pbmc_2')
rna_dirs <- c(
  "pbmc_1" = "/path/to/pbmc_1/filtered_feature_bc_matrix",
  "pbmc_2" = "/path/to/pbmc_2/filtered_feature_bc_matrix"
)

frag_files <- c(
  "pbmc_1" = "/path/to/pbmc_1/pbmc_1_A_fragments.tsv.gz",
  "pbmc_2" = "/path/to/pbmc_2/pbmc_2_A_fragments.tsv.gz"
)
gtf_file <- "/path/to/genes.gtf"   #人参考基因组里面的genes.gtf注释文件

3. 数据加载与预处理

本部分从 ATAC 的 fragment 文件和 RNA 的矩阵文件构建初始分析对象，为多组学整合提供数据结构基础。

3.1 初始化 ArchRProject

首先利用 ATAC fragment 文件构建 ArchR 专用的 Arrow 文件格式，这是 ArchR 处理 ATAC 数据的标准输入格式。通过 createArrowFiles 函数从每个样本的片段文件生成 Arrow 文件，包含基因组 tile 矩阵和基因打分矩阵。随后通过 ArchRProject 函数初始化项目对象，为后续分析建立容器。

# --------------------------
# 4. 构建 ArchR 项目
# --------------------------
message("[3/5] 创建 Arrow files ...")

ArrowFiles <- createArrowFiles(
  inputFiles = frag_files,
  sampleNames = names(frag_files),
  minTSS = 1,
  minFrags = 200,
  addTileMat = TRUE,
  addGeneScoreMat = TRUE
)

message("[4/5] 初始化 ArchRProject ...")
projMulti <- ArchRProject(ArrowFiles = ArrowFiles, outputDirectory = "Save-ProjMulti", copyArrows = TRUE)

💡 Note说明：ArchR 默认会过滤 TSS 富集得分低于 4 和唯一比对数小于 1000（也就是保留 TSS 富集得分大于等于 4 且唯一比对数大于等于 1000）的细胞。在实际分析过程中，建议根据数据的实际分布情况，合理设置 minTSS 和 minFrags 这两个参数来调整过滤标准。创建 Arrow 文件会在当前目录下生成一个 QualityControl 目录，这里面包括 2 个和样本相关的质控图。第一个图展示 log10(unique nuclear fragments) 与 TSS enrichment score 的散点图，虚线表示过滤阈值。第二图为 fragment 大小分布图。

3.2 构建 SummarizedExperiment 对象

由于 seekARC 数据中 RNA 与 ATAC 矩阵是分离的，而 ArchR 的 import10xFeatureMatrix 函数仅支持 Cell Ranger ARC 输出的多模态整合文件格式，因此我们需要自己手动构建 RNA 数据的 SummarizedExperiment 对象。核心步骤包括：从 GTF 文件提取基因坐标信息，读取 RNA 表达矩阵，匹配基因坐标并过滤无位置信息的基因，将行名转换为唯一基因符号，最终生成与 ArchR 的 addGeneExpressionMatrix 函数兼容的数据结构。

# --------------------------
# 2. 注释准备
# --------------------------
message("[1/5] 载入基因注释 ...")
gtf <- rtracklayer::import(gtf_file)
genes_gtf <- gtf[gtf$type == "gene"]
gene_id_name <- unique(data.frame(
  gene_id = genes_gtf$gene_id,
  gene_name = genes_gtf$gene_name,
  stringsAsFactors = FALSE
))

txdb <- makeTxDbFromGFF(gtf_file, format = "gtf", organism = "Mus musculus", taxonomyId = 10090)
genes_gr <- genes(txdb)
gene_ids <- mcols(genes_gr)$gene_id
names(genes_gr) <- gene_ids  # 确保后续可通过 gene_id 直接索引

# --------------------------
# 3. 构建 SummarizedExperiment
# --------------------------
load_rna_se <- function(sample_id) { 
  sce <- read10xCounts(rna_dirs[[sample_id]], col.names = TRUE) 
  rowData(sce)$gene_id <- rowData(sce)$ID 
  rowData(sce)$gene_name <- rowData(sce)$Symbol 
  
  se <- SummarizedExperiment( 
    assays = list(counts = counts(sce)), 
    rowData = rowData(sce), 
    colData = colData(sce) 
  ) 
  # 内部计算仍然使用 gene_id，因为它最可靠
  rownames(se) <- rowData(se)$gene_id 
  
  keep <- rowData(se)$gene_id %in% names(genes_gr) 
  if (any(!keep)) { 
    message(sprintf("样本 %s 有 %d 个基因缺少基因坐标，已自动剔除。", sample_id, sum(!keep))) 
  } 
  se <- se[keep, ] 
  
  matched_gene_ids <- rowData(se)$gene_id 
  gr <- genes_gr[matched_gene_ids] 
  matched_gene_names <- gene_id_name$gene_name[match(matched_gene_ids, gene_id_name$gene_id)] 
  mcols(gr)$gene_name <- matched_gene_names 
  
  rowRanges(se) <- gr 
  rowData(se)$symbol <- matched_gene_names 
  colnames(se) <- paste(sample_id, colnames(se), sep = "#") 
  
  # --- 修改部分：在这里将 rownames 替换为 gene name ---
  # 使用 make.unique 确保 rownames 唯一，否则在后续 ArchR 整合时会报错
  rownames(se) <- make.unique(rowData(se)$symbol)
  # ----------------------------------------------
  
  rm(sce) 
  gc() 
  se 
}

RNA_list <- lapply(samples, load_rna_se)
se_combined <- do.call(cbind, RNA_list)

3.3 将 RNA 的 SummarizedExperiment 添加到 ArchRProject 中

通过匹配 ATAC 和 RNA 数据的共同 Barcode ，筛选保留两模态共有的细胞。使用 addGeneExpressionMatrix 函数将 RNA 表达矩阵整合到 ArchRProject 中，创建基因表达矩阵，实现 RNA 与 ATAC 数据的联合存储，为后续多模态联合分析建立基础。

# --------------------------
# 5. 条形码匹配 & RNA 挂载
# --------------------------
message("[5/5] 匹配 RNA / ATAC 细胞并加入表达矩阵 ...")
# 匹配细胞ID
common_cells <- intersect(getCellNames(projMulti), colnames(se_combined))


if (length(common_cells) == 0) {
  stop("没有找到可以同时匹配的细胞，请检查条形码格式或白名单。")
}

projMulti <- subsetArchRProject(
  ArchRProj = projMulti,
  cells = common_cells,
  outputDirectory = "Save-ProjMulti2",
  force = TRUE
)

se_combined <- se_combined[, common_cells, drop = FALSE]

projMulti <- addGeneExpressionMatrix(
  input = projMulti,
  seRNA = se_combined,
  strictMatch = TRUE,
  force = TRUE
)

4. 数据质量控制

4.1 质量指标

在完成数据读取和项目初始化后，需要对每个细胞进行全面的质量控制，以去除低质量细胞，保证后续分析的可靠性。ArchR 在构建 Arrow 文件过程中会自动计算多种质量指标，涵盖 ATAC 和 RNA 两个模态。

RNA 质控指标：

Gex_nUMI：每个细胞中检测到的 RNA 总分子数（UMI 总数），反映测序深度。过低可能提示细胞捕获效率低，过高可能提示双胞。
Gex_nGenes：每个细胞中检测到的唯一基因数，反映转录组复杂度。过低通常伴随高线粒体比例，提示低质量细胞；过高可能提示双胞。
Gex_MitoRatio：线粒体基因表达比例，用于评估细胞状态。正常细胞通常 < 0.2，过高（> 0.2-0.3）可能提示受损或低质量细胞。

ATAC 质控指标：

TSSEnrichment：转录起始位点富集分数，评估开放染色质信号在基因启动子区域的富集程度。高值表明数据质量良好，通常 > 2 为合格，> 4 为优质。
nFrags：每个细胞的 ATAC 总片段数，衡量染色质开放信号的强度。过低表明测序覆盖度不足，过高可能提示多重捕获或双胞。
NucleosomeRatio：核小体信号比率，反映片段长度分布特征。正常细胞通常 < 2，过高（> 4）可能提示凋亡细胞或低质量细胞。

# 1. ATAC质量分布
plotQC <- plotGroups(
  ArchRProj = projMulti,
  groupBy = "Sample",
  colorBy = "cellColData",
  name = c("TSSEnrichment", "nFrags", "NucleosomeRatio"),
  alpha = 0.4,
  plotAs = "violin"
)

# 2. RNA质量分布
plotQC_RNA <- plotGroups(
  ArchRProj = projMulti,
  groupBy = "Sample",
  colorBy = "cellColData",
  name = c("Gex_nUMI", "Gex_MitoRatio","Gex_nGenes"),
  alpha = 0.4,
  plotAs = "violin"
)

plotQC$`nFrags`

plotQC$NucleosomeRatio

plotQC$TSSEnrichment

plotQC_RNA$Gex_MitoRatio

plotQC_RNA$Gex_nUMI

plotQC_RNA$Gex_nGenes

Note 建议重点关注以下内容：
是否存在明显的异常值、长尾分布或双峰分布：在质控指标的可视化中，注意检查是否有极端高值或低值的离群细胞，这些可能代表低质量细胞、双胞体或技术伪影。长尾或双峰分布可能提示样本中存在不同状态的细胞亚群。
不同样本之间质控指标分布是否一致：比较各样本在关键质控指标（如 TSS 富集分数、片段数、线粒体比例等）上的分布差异。显著的批次差异可能需要在降维前进行批次矫正。
根据可视化结果合理调整过滤阈值：避免使用固定的绝对阈值，而应根据数据的实际分布设置过滤标准，在保留足够细胞数的同时去除明显的低质量细胞，以提升下游聚类和 UMAP 可视化的清晰度与可靠性。

4.2 低质量细胞过滤

根据前面小提琴图展示的质控指标，对低质量细胞进行过滤，以去除异常细胞并提高后续分析结果的可靠性。具体过滤阈值需要结合样本实际情况进行设置，通常可参考质控指标的小提琴图和散点图分布结果进行调整。

projMulti_filtered <- subsetArchRProject(
  ArchRProj = projMulti,
  cells = getCellNames(projMulti)[
    projMulti$TSSEnrichment > 1 &
    projMulti$nFrags > 500 &
    projMulti$nFrags < 50000 &
    projMulti$NucleosomeRatio < 1 &
    projMulti$Gex_nUMI > 500 &
    projMulti$Gex_nUMI < 30000 &
    projMulti$Gex_nGenes > 200 &
    projMulti$Gex_nGenes < 8000 &
    projMulti$Gex_MitoRatio < 0.20
  ],
  outputDirectory = "Filtered_Project",
  force = TRUE
)

4.3 双胞识别

低质量细胞过滤后，还可以利用 ArchR 提供的自动化的双胞检测和过滤功能，进一步过滤掉潜在的双胞（Doublets）。

通过 addDoubletScores() 函数计算每个细胞的双胞评分，该函数基于细胞的染色质开放特征（模拟双胞）预测其为双胞的可能性。
使用 filterDoublets() 函数根据双胞评分过滤细胞。函数中关键的 filterRatio 参数控制了过滤的严格程度，其定义为：基于通过质控的细胞数计算的预测双胞最大过滤比例。

projMulti_filtered <- addDoubletScores(projMulti_filtered, force = TRUE)

projMulti_filtered <- filterDoublets(projMulti_filtered)

5. 线性降维分析

在构建多模态数据对象后，分别对 ATAC 和 RNA 数据进行降维处理，以提取核心数据结构并减少计算维度。ArchR 采用迭代隐语义索引（LSI）算法，可有效处理高维稀疏的染色质开放数据和基因表达数据。

5.1 迭代 LSI 降维

针对两种数据类型，采用不同的特征选择策略：

ATAC 数据：基于全基因组 500-bp 窗口构建 tile 矩阵，通过迭代 LSI 逐步聚焦于变异最大的染色质开放区域。
RNA 数据：基于高变异基因构建表达矩阵，通过 LSI 识别主要的转录组模式。

5.2 多模态数据整合

将 ATAC 和 RNA 的 LSI 降维结果通过 addCombinedDims 函数进行整合，生成联合降维结果 LSI_Combined，为后续的多模态联合分析提供统一的基础。

#对ATAC数据进行lsi降维分析
projMulti_filtered <- addIterativeLSI(
  ArchRProj = projMulti_filtered, 
  clusterParams = list(
    resolution = 0.2, 
    sampleCells = 10000,
    n.start = 10
  ),
  saveIterations = FALSE,
  useMatrix = "TileMatrix", 
  depthCol = "nFrags",
  name = "LSI_ATAC"
)

#对RNA数据进行lsi降维分析
projMulti_filtered <- addIterativeLSI(
  ArchRProj = projMulti_filtered, 
  clusterParams = list(
    resolution = 0.2, 
    sampleCells = 10000,
    n.start = 10
  ),
  saveIterations = FALSE,
  useMatrix = "GeneExpressionMatrix", 
  depthCol = "Gex_nUMI",
  varFeatures = 2500,
  firstSelection = "variable",
  binarize = FALSE,
  name = "LSI_RNA"
)

#将 ATAC 和 RNA 的 LSI 降维结果进行整合
projMulti_filtered <- ArchR::addCombinedDims(projMulti_filtered, 
                                             reducedDims = c("LSI_ATAC", "LSI_RNA"),   
                                             name = "LSI_Combined")

6. 非线性降维与聚类分析

在完成 LSI 降维与多模态整合后，我们通过非线性降维方法（UMAP）将高维数据映射到二维空间，以可视化细胞间的复杂关系。UMAP 能够有效保留高维空间中的局部与全局结构，便于直观观察细胞亚群分布。

6.1 三种降维策略

整合模态降维：基于整合后的多模态 LSI 结果（LSI_Combined）进行 UMAP 可视化，体现 ATAC 和 RNA 数据的联合信号。
ATAC 单模态降维：基于纯 ATAC 的 LSI 结果（LSI_ATAC）进行 UMAP 可视化，突出染色质可及性模式。
RNA 单模态降维：基于纯 RNA 的 LSI 结果（LSI_RNA）进行 UMAP 可视化，反映基因表达谱相似性。

projMulti_filtered <- addUMAP(projMulti_filtered, reducedDims = "LSI_Combined", name = "UMAP_Combined", minDist = 0.8, force = TRUE)
projMulti_filtered <- addUMAP(projMulti_filtered, reducedDims = "LSI_ATAC", name = "UMAP_ATAC", minDist = 0.8, force = TRUE)
projMulti_filtered <- addUMAP(projMulti_filtered, reducedDims = "LSI_RNA", name = "UMAP_RNA", minDist = 0.8, force = TRUE)
projMulti_filtered <- addClusters(projMulti_filtered, reducedDims = "LSI_Combined", name = "Clusters_Combined", resolution = 0.4, force = TRUE)
projMulti_filtered <- addClusters(projMulti_filtered, reducedDims = "LSI_ATAC", name = "Clusters_ATAC", resolution = 0.4, force = TRUE)
projMulti_filtered <- addClusters(projMulti_filtered, reducedDims = "LSI_RNA", name = "Clusters_RNA", resolution = 0.4, force = TRUE)

Note
参数说明
minDist = 0.8：控制 UMAP 中点的最小间距，较低值可分离紧密的细胞亚群，较高值可呈现更连续的细胞状态过渡。通常设置为 0.1~0.8。
resolution = 0.4：控制聚类的粒度，较低值产生较粗的聚类（较少簇），较高值产生更精细的亚群划分（通常基于 Seurat 的 FindClusters 算法）。

6.2 降维聚类可视化

p1 <- plotEmbedding(projMulti_filtered, name = "Clusters_ATAC", embedding = "UMAP_ATAC", size = 1, labelAsFactors=F, labelMeans=F)
p2 <- plotEmbedding(projMulti_filtered, name = "Clusters_RNA", embedding = "UMAP_RNA", size = 1, labelAsFactors=F, labelMeans=F)
p3 <- plotEmbedding(projMulti_filtered, name = "Clusters_Combined", embedding = "UMAP_Combined", size = 1, labelAsFactors=F, labelMeans=F)

options(repr.plot.width = 14, repr.plot.height = 7)
p1+p2+p3

p4 <- plotEmbedding(projMulti_filtered, name = "Sample", embedding = "UMAP_ATAC", size = 1, labelAsFactors=F, labelMeans=F)
p5 <- plotEmbedding(projMulti_filtered, name = "Sample", embedding = "UMAP_RNA", size = 1, labelAsFactors=F, labelMeans=F)
p6 <- plotEmbedding(projMulti_filtered, name = "Sample", embedding = "UMAP_Combined", size = 1, labelAsFactors=F, labelMeans=F)

options(repr.plot.width = 14, repr.plot.height = 7)
p4+p5+p6

Note 在 UMAP 可视化中观察到样本间存在明显的批次分离，表明 LSI 迭代降维虽能缓解部分由实验技术引入的批次效应，但对于显著的样本间技术变异，其矫正能力可能有限。若批次效应仍干扰下游生物学信号的解析，建议考虑在 LSI 降维后额外应用专门设计的批次矫正算法（如 Harmony）进行进一步处理。

6.3 用 Harmony 分别对 ATAC 和 RNA 进行批次矫正

在多样本数据分析中，若 LSI 降维后仍观察到明显的样本间批次分离，说明需要额外的批次矫正步骤。ArchR 内置了 Harmony 算法，该算法专为单细胞数据设计，通过主成分的迭代聚类与线性校正，能够有效整合来自不同样本或实验批次的数据。

Harmony 在 ArchR 中通过 addHarmony() 函数调用，用户需指定 groupBy 参数来定义需要矫正的批次变量（通常为样本来源）。

本次分析分别对 ATAC 和 RNA 数据进行独立的 Harmony 批次矫正：

ATAC 数据矫正：基于 LSI_ATAC 降维结果，生成矫正后的嵌入 Harmony_ATAC。
RNA 数据矫正：基于 LSI_RNA 降维结果，生成矫正后的嵌入 Harmony_RNA。

最后，将两个矫正后的降维结果通过 addCombinedDims() 进行整合，生成联合的批次矫正降维表示 LSI_Combined_Harmony，用于后续的多模态联合分析。

# 1. 对ATAC数据进行Harmony批次矫正
projMulti_filtered <- addHarmony(
    ArchRProj = projMulti_filtered,
    reducedDims = "LSI_ATAC",       # ✔️ 输入：使用您刚刚创建的ATAC LSI结果
    name = "Harmony_ATAC",           # ✔️ 输出：新名称，避免覆盖
    groupBy = "Sample",              # 按样本批次进行矫正
    theta = 5, 
    sigma = 0.01,
    force = TRUE
)

# 2. 对RNA数据进行Harmony批次矫正
projMulti_filtered <- addHarmony(
    ArchRProj = projMulti_filtered,
    reducedDims = "LSI_RNA",         # ✔️ 输入：使用您刚刚创建的RNA LSI结果
    name = "Harmony_RNA",            # ✔️ 输出：新名称
    groupBy = "Sample",
    theta = 5,
    sigma = 0.01,
    force = TRUE
)

projMulti_filtered <- ArchR::addCombinedDims(projMulti_filtered, 
                                             reducedDims = c("Harmony_ATAC", "Harmony_RNA"),   
                                             name = "LSI_Combined_Harmony")

6.4 基于批次矫正结果的降维与聚类

完成 Harmony 批次矫正后，基于矫正后的降维结果重新执行 UMAP 非线性降维和细胞聚类，以评估批次矫正效果并获取更准确的细胞分群。

降维与聚类策略

ATAC 批次矫正结果：基于 Harmony_ATAC 进行 UMAP 降维和聚类，得到 UMAP_Harmony_ATAC 和 Clusters_Harmony_ATAC，反映去除批次效应后的染色质开放模式。
RNA 批次矫正结果：基于 Harmony_RNA 进行 UMAP 降维和聚类，得到 UMAP_Harmony_RNA 和 Clusters_Harmony_RNA，反映去除批次效应后的转录组特征。
联合批次矫正结果：基于整合的 LSI_Combined_Harmony 进行 UMAP 降维和聚类，得到 UMAP_Combined_Harmony 和 Clusters_Combined_Harmony，综合考量去批次后的多模态信息。

# 基于Harmony矫正结果进行UMAP降维
projMulti_filtered <- addUMAP(projMulti_filtered, reducedDims = "Harmony_ATAC", 
                              name = "UMAP_Harmony_ATAC", minDist = 0.8, force = TRUE)
projMulti_filtered <- addUMAP(projMulti_filtered, reducedDims = "Harmony_RNA", 
                              name = "UMAP_Harmony_RNA", minDist = 0.8, force = TRUE)
projMulti_filtered <- addUMAP(projMulti_filtered, reducedDims = "LSI_Combined_Harmony", 
                              name = "UMAP_Combined_Harmony", minDist = 0.8, force = TRUE)

# 基于Harmony矫正结果进行聚类
projMulti_filtered <- addClusters(projMulti_filtered, reducedDims = "Harmony_ATAC", 
                                  name = "Clusters_Harmony_ATAC", resolution = 0.4, force = TRUE)
projMulti_filtered <- addClusters(projMulti_filtered, reducedDims = "Harmony_RNA", 
                                  name = "Clusters_Harmony_RNA", resolution = 0.4, force = TRUE)
projMulti_filtered <- addClusters(projMulti_filtered, reducedDims = "LSI_Combined_Harmony", 
                                  name = "Clusters_Combined_Harmony", resolution = 0.4, force = TRUE)

# 样本分布可视化（评估批次效应）
p7 <- plotEmbedding(projMulti_filtered, name = "Sample", embedding = "UMAP_Harmony_ATAC", size = 1, labelAsFactors=F, labelMeans=F)
p8 <- plotEmbedding(projMulti_filtered, name = "Sample", embedding = "UMAP_Harmony_RNA", size = 1, labelAsFactors=F, labelMeans=F)
p9 <- plotEmbedding(projMulti_filtered, name = "Sample", embedding = "UMAP_Combined_Harmony", size = 1, labelAsFactors=F, labelMeans=F)

options(repr.plot.width = 14, repr.plot.height = 7)
p7+p8+p9

# 加载拼图包
library(cowplot)

# 1. 获取所有样本名称
sample_list <- unique(projMulti_filtered$Sample)

# 2. 循环绘制：每个样本单独高亮，背景为灰色
plot_list <- lapply(sample_list, function(sample_i) {
  plotEmbedding(
    ArchRProj = projMulti_filtered,
    embedding = "UMAP_Combined_Harmony",      # 使用矫正后的UMAP
    colorBy = "cellColData",
    name = "Clusters_Combined_Harmony",       # 使用矫正后的聚类
    highlightCells = getCellNames(projMulti_filtered)[projMulti_filtered$Sample == sample_i],
    pal = c("Non.Highlighted" = "lightgrey"), # 非当前样本设为灰色
    size = 0.8,
    baseSize = 8,
    title = sample_i                          # 标题设为样本名
  )
})

# 3. 组合分面图 (调整 ncol 控制每行显示数量)
p_faceted <- plot_grid(plotlist = plot_list, ncol = 3)
print(p_faceted)

Note: 批次矫正效果评估

批次矫正效果的评估需基于可视化结果进行判断：
混合度检查：对比批次矫正前后（UMAP_Combined vs. UMAP_Combined_Harmony）的样本分布，理想情况下矫正后不同样本的细胞在UMAP空间中应更均匀地混合，而非按样本分离。
生物学信号保留：检查批次矫正是否保留了清晰的细胞类型分离。过度矫正会导致细胞亚群界限模糊，而矫正不足则难以有效整合样本。
聚类一致性：观察同一细胞类型在不同样本中是否被分配到相同的聚类，确保批次矫正未引入新的技术偏差。
下游分析选择：若矫正效果理想，推荐使用批次矫正后的聚类结果（如 Clusters_Combined_Harmony）进行后续的差异分析和功能注释，以获得更稳定、可重复的生物学发现。

# 批次矫正后的可视化（使用Harmony矫正结果）
p10 <- plotEmbedding(projMulti_filtered, name = "Clusters_Harmony_ATAC", embedding = "UMAP_Harmony_ATAC", size = 1, labelAsFactors=F, labelMeans=F)
p11 <- plotEmbedding(projMulti_filtered, name = "Clusters_Harmony_RNA", embedding = "UMAP_Harmony_RNA", size = 1, labelAsFactors=F, labelMeans=F)
p12 <- plotEmbedding(projMulti_filtered, name = "Clusters_Combined_Harmony", embedding = "UMAP_Combined_Harmony", size = 1, labelAsFactors=F, labelMeans=F)

options(repr.plot.width = 14, repr.plot.height = 7)
p10+p11+p12

7. 分组 Peak Calling 与矩阵构建

在完成细胞聚类（Clusters_Combined_Harmony）分析后，基于细胞类型分群结果，我们采用分组策略重新进行染色质开放区域（Peak）识别，并构建细胞 × Peak 的计数矩阵，为后续的调控元件分析提供高分辨率输入。

7.1 分组覆盖度计算

addGroupCoverages() 函数根据指定的细胞分组（groupBy = "Clusters_Combined_Harmony"），将同一聚类内所有细胞的 ATAC-seq 片段进行合并，生成各细胞亚群的聚合覆盖度文件。此步骤旨在为每个细胞类型积累充分的信号深度，确保后续 Peak Calling 的灵敏度，尤其有助于捕获稀有细胞亚群中开放度较弱但生物学意义重要的区域。

addArchRThreads(4)

projMulti_filtered <- addGroupCoverages(ArchRProj = projMulti_filtered, groupBy = "Clusters_Combined_Harmony", verbose = FALSE)

7.2 分组 Peak Calling

addReproduciblePeakSet() 是核心执行步骤，采用分组独立识别策略：

独立识别：在每个细胞亚群内部独立运行 MACS2（需指定 pathToMacs2 路径）进行 Peak Calling，从而识别出该细胞类型特有的开放区域。
合并去重：将所有亚群识别出的 Peak 坐标合并，生成一个非冗余的统一 Peak 集合（Union Peak Set）。相比在全数据集上一次 Calling，此策略能更稳健地捕获细胞类型特异的调控元件，避免信号被稀释或掩盖。

pathToMacs2 <- "/PROJ/development/changyanxia/software/bin/macs2"
projMulti_filtered <- addReproduciblePeakSet(ArchRProj = projMulti_filtered, groupBy = "Clusters_Combined_Harmony", pathToMacs2 = pathToMacs2)

7.3 Peak 矩阵构建

addPeakMatrix() 函数基于整合后的 Peak 集合，为每个细胞计算在这些区域内的片段计数，构建细胞 × Peak 的开放度矩阵。此矩阵是后续差异 Peak 分析、Motif 富集、Peak‑to‑Gene 关联等表观遗传分析的核心数据基础。

projMulti_filtered <- addPeakMatrix(ArchRProj = projMulti_filtered)

8. 细胞类型注释

8.1 基因表达差异分析

识别各细胞簇中特异性高表达的基因，基于表达比例和差异倍数筛选代表性标记基因，为细胞类型鉴定提供初步依据。可通过热图（Heatmap）查看不同 cluster 中 Top 差异基因的表达情况，从而辅助判断各个 cluster 对应的细胞类型。

#getCellColData(projMulti_filtered)[1:6,]
#查看可以用于差异分析的矩阵有哪些
addArchRThreads(threads = 1) 
getAvailableMatrices(projMulti_filtered)
DEG <- getMarkerFeatures(
    ArchRProj = projMulti_filtered,
    groupBy = "Clusters_Combined_Harmony",
    #groupBy = "Clusters_Combined",
    useMatrix = "GeneExpressionMatrix",
    bias = c("TSSEnrichment", "log10(nFrags)", "log10(Gex_nUMI)")  # 包含RNA偏差
)

markerList <- getMarkers(DEG, cutOff = "FDR < 0.05 & Log2FC > 0.5")
#markerList <- markerList[c(1:11,14)]
# 2. 提取每个cluster的Top 4差异基因
top_n <- 4
top_markers_list <- lapply(markerList, function(df) {
  df_sorted <- df[order(df$Log2FC, decreasing = TRUE), ]
  return(head(df_sorted, top_n))
})
for(cluster in names(top_markers_list)){
  top_markers_list[[cluster]]$Cluster <- cluster
}
all_top_markers <- do.call(rbind, top_markers_list)

heatmap_plot1 <- plotMarkerHeatmap(
  seMarker = DEG,
  cutOff = "FDR <= 0.1 & Log2FC >= 1",
  binaryClusterRows = TRUE,
  labelMarkers = unique(all_top_markers$name),  # 标记这些基因
  clusterCols = TRUE,
  transpose = FALSE
)

heatmap_plot1

8.2 在 scRNA-seq 数据层面查看 marker 基因的表达情况

本示例数据为 PBMC，因此整理了常见免疫细胞类型及其对应的 marker 基因集。实际分析中，应根据样本来源和组织类型，自定义相应的 marker 基因。后续可通过热图查看不同 cluster 中 marker 基因的表达情况，从而辅助判断各个 cluster 对应的细胞类型。

pbmc_marker_integrated <- list(
  "Plasma Cells" = c("IGHG1", "IGKC"),
  "Monocytes" = c("VCAN", "FCN1", "TREM1"),
  "B cells" = c("MS4A1", "BACH2", "PAX5"),
  "CMP" = c("MPO"),
  "Pro B cells" = c("DNTT"),
  "Erythroblast" = c("ALAS2", "NFIA", "SOX6"),
  "T cells" = c("PLXNA4", "THEMIS", "ITK", "BCL11B"),
  "NK cells" = c("NCAM1", "KLRC3"),
  "Dividing B cells" = c("MKI67", "TOP2A"),
  "pDC" = c("RUNX2")
)

heatmap_plot <- plotMarkerHeatmap(
  seMarker = DEG,  # 这里se是你的差异分析结果，如getMarkerFeatures的输出
  cutOff = "FDR <= 0.1 & Log2FC >= 1",  # 设置显著性阈值
  labelMarkers = unlist(pbmc_marker_integrated),  # 标记这些基因
  binaryClusterRows = TRUE,
  clusterCols = TRUE,
  transpose = FALSE
)

heatmap_plot

8.3 在 scATAC-seq 数据层面查看 marker 基因组区域开放情况

除 RNA 表达外，还可以从 scATAC-seq 数据层面查看 marker 基因附近调控区域的开放情况，以进一步辅助细胞类型注释。通常可结合基因附近的 peak 信号、染色质开放程度以及不同 cluster 间的可及性差异，判断特定 marker 基因是否在对应细胞群中具有活跃的调控状态。

my_clusters <- sort(unique(projMulti_filtered$Clusters_Combined_Harmony))

p1 <- plotBrowserTrack(
  ArchRProj = projMulti_filtered,
  geneSymbol  = c("VCAN", "RUNX2", "MPO", "ALAS2", "THEMIS", "KLRC3"),  # 基因区域
  groupBy = "Clusters_Combined_Harmony",  # 按细胞类型分组
  #groupBy = "Clusters_Combined",  
  useGroups = my_clusters,  # 确保只显示这些，且顺序固定
  tileSize = 100,
  upstream = 5000,
  downstream = 5000,
  title = "VCAN Region"
)

options(repr.plot.width = 16, repr.plot.height = 8)
patchwork::wrap_plots(p1)

8.4 细胞类型标注及其可视化

将已推断出的细胞类型名称，通过 addCellColData 映射并写入到 ArchRProject 中，方便在后续的可视化与差异分析中直接按照实际生物学意义的分类（如 T cells、B cells）进行展示。

cat('开始细胞类型重命名...\n')

# 1. 先查看实际的聚类名称
print("当前聚类标签：")
print(table(projMulti_filtered$UMAP_Combined_Harmony))

# 2. 建立实际的映射关系
# 假设你的聚类名是 "C1", "C2" 等格式
celltype_mapping <- c(
  "C1" = "T cells", "C2" = "Monocytes", "C3" = "NK cells",
  "C4" = "Erythroblast", "C5" = "B cells", "C6" = "T cells",
  "C7" = "B cells", "C8" = "Dividing B cells", "C9" = "Pro B cells",
  "C10" = "CMP", "C11" = "T cells", "C12" = "pDC",
  "C13" = "T cells", "C14" = "B cells", "C15" = "Monocytes",
  "C16" = "Plasma Cells", "C17" = "NK cells", "C18" = "NA","C19" = "NA","C20" = "NA"
)

# 3. 执行映射
celltype_annotations <- celltype_mapping[as.character(projMulti_filtered$UMAP_Combined_Harmony)]

# 4. 添加到项目
projMulti_filtered <- addCellColData(
  ArchRProj = projMulti_filtered,
  data = celltype_annotations,
  name = "CellType",
  cells = getCellNames(projMulti_filtered)
)

print(table(projMulti_filtered$CellType))

p1=plotEmbedding(projMulti_filtered, name = "CellType", embedding = "UMAP_Combined_Harmony", size = 1, labelAsFactors=F, labelMeans=F)

options(repr.plot.width = 10, repr.plot.height = 10)
p1

9. 保存结果

projMulti_filtered <- saveArchRProject(ArchRProj = projMulti_filtered, outputDirectory = "Save-ProjMulti2", overwrite = TRUE, load = TRUE)

ATAC + RNA 多组学：基于 ArchR 的基础分析 ​

1. 教程简介 ​

2. 输入文件准备 ​

输入文件要求 ​

目录结构示例 ​

3. 数据加载与预处理 ​

3.1 初始化 ArchRProject ​

3.2 构建 SummarizedExperiment 对象 ​

3.3 将 RNA 的 SummarizedExperiment 添加到 ArchRProject 中 ​

4. 数据质量控制 ​

4.1 质量指标 ​

4.2 低质量细胞过滤 ​

4.3 双胞识别 ​

5. 线性降维分析 ​

5.1 迭代 LSI 降维 ​

5.2 多模态数据整合 ​

6. 非线性降维与聚类分析 ​

6.1 三种降维策略 ​

参数说明 ​

6.2 降维聚类可视化 ​

6.3 用 Harmony 分别对 ATAC 和 RNA 进行批次矫正 ​

6.4 基于批次矫正结果的降维与聚类 ​

7. 分组 Peak Calling 与矩阵构建 ​

7.1 分组覆盖度计算 ​

7.2 分组 Peak Calling ​

7.3 Peak 矩阵构建 ​

8. 细胞类型注释 ​

8.1 基因表达差异分析 ​

8.2 在 scRNA-seq 数据层面查看 marker 基因的表达情况 ​

8.3 在 scATAC-seq 数据层面查看 marker 基因组区域开放情况 ​

8.4 细胞类型标注及其可视化 ​

9. 保存结果 ​