ArchR 单细胞 ATAC 分析 — 高性能单细胞染色质可及性全流程分析¶

# 安装 ArchR（需要 R 4.1+ 和 Bioconductor 3.13+）
devtools::install_github("GreenleafLab/ArchR", ref="master",
                          repos = BiocManager::repositories())

# 安装所有额外依赖
ArchR::installExtraPackages()

# 安装 Cairo（图形后端，用于 PDF/PNG 输出）
install.packages("Cairo")

# 设置并行核心数（推荐设置为 CPU 核数 - 2）
library(ArchR)
addArchRThreads(threads = 8)           # 设置 8 个并行线程

# 设置基因组（hg38 人类，mm10 小鼠）
addArchRGenome("hg38")                 # 下载并设置 hg38 基因组注释

# 验证安装
packageVersion("ArchR")               # 应显示 "1.0.3" 或以上

library(ArchR)
set.seed(1)                            # 固定随机种子，保证结果可重复

# ── 创建 ArchRProject ─────────────────────────────────
# 准备片段文件列表（每个样本一个 fragments.tsv.gz 文件）
inputFiles <- c(
  'sample1' = 'sample1/fragments.tsv.gz',   # Cellranger ATAC 输出
  'sample2' = 'sample2/fragments.tsv.gz'
)

# 创建箭头文件（Arrow Files）：ArchR 的核心数据存储格式
# 每个样本生成一个 .arrow 文件，存储压缩的 tile 矩阵和 peak 矩阵
ArrowFiles <- createArrowFiles(
  inputFiles = inputFiles,             # 输入片段文件
  sampleNames = names(inputFiles),     # 样本名
  filterTSS = 4,                       # TSS 富集分数过滤阈值（TSS score > 4）
  filterFrags = 1000,                  # 最少片段数阈值（过滤低质量细胞）
  addTileMat = TRUE,                   # 添加 tile 矩阵（500bp 窗口）
  addGeneScoreMat = TRUE               # 添加基因得分矩阵（用于与 RNA 比较）
)

# 创建 ArchRProject 对象（包含所有样本的统一分析空间）
proj <- ArchRProject(
  ArrowFiles = ArrowFiles,             # 箭头文件列表
  outputDirectory = 'ArchR_output/',   # 输出目录
  copyArrows = TRUE                    # 复制箭头文件到输出目录
)

# 查看项目基本信息
proj                                   # 打印项目摘要
print(paste0("细胞数：", nCells(proj)))

library(ArchR)
addArchRThreads(8)                     # 8 线程并行
addArchRGenome("hg38")

# ── 完整 ArchR 分析流程 ──────────────────────────────

# 假设已创建 proj（从上面的 ArchRProject 继续）

# 1. 质控可视化
# 双细胞检测（预测 doublet）
proj <- addDoubletScores(
  input = proj,                        # 输入 ArchRProject
  k = 10,                              # KNN 近邻数
  knnMethod = "UMAP",                  # 用 UMAP 空间判断 doublet
  LSIMethod = 1                        # LSI 方法（1 = TF-IDF 后 PCA）
)
proj <- filterDoublets(proj)           # 过滤预测的 doublet 细胞

# 2. 降维（LSI = Latent Semantic Indexing）
# LSI 是 ATAC-seq 的标准降维方法（相当于 scRNA-seq 的 PCA）
proj <- addIterativeLSI(
  ArchRProj = proj,
  useMatrix = "TileMatrix",            # 基于 500bp tile 矩阵
  name = "IterativeLSI",              # 降维结果的键名
  iterations = 2,                      # 迭代次数（2 轮足够）
  clusterParams = list(
    resolution = c(0.2),               # 第一轮聚类分辨率
    sampleCells = 10000,               # 第一轮用于聚类的细胞数
    n.start = 10                       # 重启次数
  ),
  varFeatures = 25000,                 # 使用方差最大的 25000 个 feature
  dimsToUse = 1:30                     # 使用前 30 个 LSI 维度
)

# 3. 批次校正（可选，多样本整合）
proj <- addHarmony(
  ArchRProj = proj,
  reducedDims = "IterativeLSI",        # 输入 LSI 嵌入
  name = "Harmony",                    # 输出 Harmony 校正嵌入的键名
  groupBy = "Sample"                   # 按样本批次校正
)

# 4. 聚类
proj <- addClusters(
  input = proj,
  reducedDims = "Harmony",             # 用 Harmony 校正后的嵌入
  method = "Seurat",                   # 聚类方法（Seurat = Leiden 算法）
  name = "Clusters",                   # 聚类结果列名
  resolution = 0.8                     # 聚类分辨率
)

# 5. UMAP 可视化
proj <- addUMAP(
  ArchRProj = proj,
  reducedDims = "Harmony",
  name = "UMAP",
  nNeighbors = 30,                     # KNN 近邻数
  minDist = 0.5                        # UMAP 最小距离
)

# 绘制 UMAP
p1 <- plotEmbedding(
  ArchRProj = proj,
  colorBy = "cellColData",             # 按元数据着色
  name = "Clusters",                   # 显示聚类信息
  embedding = "UMAP"
)
p2 <- plotEmbedding(
  ArchRProj = proj,
  colorBy = "cellColData",
  name = "TSSEnrichment",              # 显示 TSS 富集分数（质控指标）
  embedding = "UMAP"
)
plotPDF(p1, p2, name = "UMAP_plots.pdf",
        ArchRProj = proj, addDOC = FALSE, width = 5, height = 5)

# 6. Peak Calling（调用 MACS3 检测 peaks）
proj <- addGroupCoverages(
  ArchRProj = proj,
  groupBy = "Clusters"                 # 按 cluster 合并生成覆盖度轨迹
)
proj <- addReproduciblePeakSet(
  ArchRProj = proj,
  groupBy = "Clusters",               # 按 cluster 分别 call peaks
  pathToMacs2 = "macs3",             # MACS3 路径（如果安装了 MACS3）
  genomeSize = 'hs',                  # 基因组大小（hs=人类，mm=小鼠）
  maxPeaks = 150000                   # 每组最多 150000 个 peak
)
proj <- addPeakMatrix(proj)           # 计算 peak 矩阵（细胞 × peak）

# 7. 差异可及性分析
markersPeaks <- getMarkerFeatures(
  ArchRProj = proj,
  useMatrix = "PeakMatrix",           # 基于 peak 矩阵分析
  groupBy = "Clusters",               # 按 cluster 比较
  testMethod = "wilcoxon",            # Wilcoxon 检验（比 T 检验更稳健）
  bias = c("TSSEnrichment", "log10(nFrags)")  # 校正 TSS 富集和片段数的偏差
)

# 提取显著差异 peaks
markerList <- getMarkers(
  markersPeaks,
  cutOff = "FDR <= 0.01 & Log2FC >= 1"  # FDR<1% 且 log2倍变 >1
)
print(markerList)

# 8. 基序富集分析（鉴定驱动开放染色质的 TF）
proj <- addMotifAnnotations(
  ArchRProj = proj,
  motifSet = "cisbp",                 # CisBP 基序数据库（人类/小鼠）
  name = "Motif"
)
enrichMotifs <- peakAnnoEnrichment(
  seMarker = markersPeaks,
  ArchRProj = proj,
  peakAnnotation = "Motif",
  cutOff = "FDR <= 0.1 & Log2FC >= 0.5"
)
heatmapEM <- plotEnrichHeatmap(
  enrichMotifs, n = 7, transpose = TRUE   # 每个 cluster 展示 top7 基序
)
plotPDF(heatmapEM, name = "MotifEnrichment.pdf", ArchRProj = proj)

# 9. 保存项目
saveArchRProject(ArchRProj = proj, outputDirectory = 'ArchR_output/')

# 设置
addArchRGenome("hg38")
addArchRThreads(8)

# 核心流程（6步）
ArrowFiles <- createArrowFiles(inputFiles, sampleNames, filterTSS=4, filterFrags=1000)
proj <- ArchRProject(ArrowFiles, outputDirectory="out/")
proj <- addIterativeLSI(proj, useMatrix="TileMatrix", name="IterativeLSI")
proj <- addClusters(proj, reducedDims="IterativeLSI")
proj <- addUMAP(proj, reducedDims="IterativeLSI")
proj <- addGroupCoverages(proj, groupBy="Clusters")
proj <- addReproduciblePeakSet(proj, groupBy="Clusters")