ATAC-seq领域超级大佬William J. Greenleaf团队开发的scATAC-seq分析软件：ArchR

生信技能树

发布于 2025-05-08 07:29:16

17900

代码可运行

文章被收录于专栏：生信技能树生信技能树

运行总次数：0

代码可运行

scATAC-seq数据分析除了经典的 Signac软件，还有一款也使用超多的软件 ArchR，官网给到你：

https://www.archrproject.com/index.html

这个软件于2021年2月25号发表在 Nature Genetics 杂志上，文献标题为《ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis》。作者认为 ArchR 是目前最全面且用户友好的单细胞ATAC-seq和多组学数据分析工具。开发团队将继续改进和扩展ArchR的功能，以保持其在该领域的领先地位。

分析流程：

大家看文章的通讯作者就知道是 ATAC-seq 这个领域的大佬了：William J. Greenleaf

ATAC-seq方法是 2013 年由斯坦福大学的 William J. Greenleaf 和 Howard Y. Chang 实验室共同开发的用于研究染色质可及性/开放性的方法。第一个scATAC-seq数据是 2015 年由Greenleaf(Buenrostro, Wu et al. 2015)和Shendure(Cusanovich, Daza et al. 2015)实验室的分别发布Nature和Science期刊上，他们通过修改ATAC-seq protocal获取了几百~上万个细胞。其中Greenleaf实验室Nature文章中是依赖物理隔离单细胞，而Shendure实验室避免了单细胞反应体积使用两步组合索引策略。

fragment概念

在ATAC-seq数据中，fragment 指由两次转座事件产生的可测序DNA分子，每个片段的两端通过配对末端测序进行测序。

在 ArchR 中, fragments 指的是一个表格或基因组范围对象，其中包含每个测序片段对应的染色体、经过偏移调整的单碱基染色体起始位置、经过偏移调整的单碱基染色体终止位置以及唯一的细胞条形码ID。

ATAC-seq更加详细的概念可以看这篇文献：https://pmc.ncbi.nlm.nih.gov/articles/PMC3959825/

软件安装

我这里安装的github上面的最新版：

## 使用西湖大学的 Bioconductor镜像
options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor")
options("repos"=c(CRAN="https://mirrors.westlake.edu.cn/CRAN/"))
# 安装
devtools::install_github("GreenleafLab/ArchR", ref="master", repos = BiocManager::repositories())

# 安装完后，设置一下环境
rm(list=ls())
library(ArchR)
set.seed(1)

addArchRThreads(threads = 60) 
## Setting default number of Parallel threads to 1.

addArchRGenome("hg19")
## Setting default genome to Hg19.
# addArchRGenome("hg38")

0.输入数据

可以接受两种格式的文件：fragment files and BAM files

Fragment files：Fragment files are tabix-sorted text files containing each scATAC-seq fragment and the corresponding cell ID, one fragment per line.

BAM files：BAM files are binarized tabix-sorted files that contain each scATAC-seq fragment, raw sequence, cellular barcode id and other information.

示例数据

来自文献：https://pubmed.ncbi.nlm.nih.gov/31792411/

可以在GEO上下载得到：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE139369

这个数据的分析代码也可以下载：https://github.com/GreenleafLab/MPAL-Single-Cell-2019.

数据包括三个样本：

BMMC（Bone Marrow Mononuclear Cells）：骨髓单个核细胞，是从骨髓中分离的单个核细胞群体，包括造血干细胞、祖细胞以及其他免疫细胞。
PBMC（Peripheral Blood Mononuclear Cells）：外周血单个核细胞，是从外周血中分离的单个核细胞群体，主要包括淋巴细胞、单核细胞等。
CD34+ BMMC（CD34+ Hematopoietic Stem and Progenitor Cells from Bone Marrow）：来自骨髓的CD34+造血干细胞和祖细胞，这类细胞具有自我更新和分化为多种血细胞类型的能力，是造血系统的重要组成部分。

####### 示例数据下载
library(parallel)
inputFiles <- getTutorialData("Hematopoiesis")
inputFiles
names(inputFiles)

1.创建Arrow Files

Arrow 文件是 ArchR 中用于存储单个样本数据的基本单位。存储内容：

元数据（Metadata）：样本的基本信息，如样本来源、细胞类型等；
可访问的片段（Accessible Fragments）：经过ATAC-seq等技术检测到的开放染色质区域的DNA片段；
数据矩阵（Data Matrices）：用于分析的数值化数据，如基因表达矩阵、片段计数矩阵等。

Arrow文件将样本的所有相关信息整合在一起，便于管理和分析，是ArchR进行单细胞数据分析的基础。Arrow Files的结构如下：

创建Arrow Files，这个步骤会对每一个样本做以下处理：

读取accessible fragments；
计算质控指标 (i.e. TSS enrichment scores and nucleosome info)；
过滤细胞；
使用500-bp bins窗口创建genome-wide TileMatrix；
创建GeneScoreMatrix。

这个过程会耗时~15 mins 左右

ArrowFiles <- createArrowFiles(
  inputFiles = inputFiles,
  sampleNames = names(inputFiles),
  minTSS = 4, # Dont set this too high because you can always increase later
  minFrags = 1000, 
  addTileMat = TRUE,
  addGeneScoreMat = TRUE
)
ArrowFiles
# [1] "scATAC_BMMC_R1.arrow"      "scATAC_CD34_BMMC_R1.arrow" "scATAC_PBMC_R1.arrow"

数据质控

scATAC-seq数据质控考虑三个指标：

unique nuclear fragments：细胞核内fragments去重数，即没有比对到线粒体DNA上的；
signal-to-background ratio：信噪比（signal-to-background ratio），低信噪比通常归因于死亡或濒死的细胞，这些细胞的DNA已经去染色质化，从而允许转座子在整个基因组范围内随机插入；计算为 the TSS enrichment score；
fragment size distribution：片段大小分布。由于核小体的周期性，我们期望看到与核小体包裹的DNA长度（大约147个碱基对）相当的片段出现缺失。

上面过滤的标准 minTSS = 4, minFrags = 1000 需要根据数据的情况进行调整。

运行完上面的步骤后目录下会生成一下文件：

TSS_by_Unique_Frags.pdf：log10(unique nuclear fragments) vs TSS enrichment score，过滤阈值使用虚线表示

the fragment size distribution：

2.Doublet预测

在几乎所有平台上生成的单细胞数据都容易出现双细胞（doublet）的情况。有单细胞转录组分析经验，doublets在这里也非常好理解，就是一个液滴中只有一个barcode但是进了两个以上的细胞。

ArchR这个软件可以预测双包体，前面的Signac软件没有这个环节。

双细胞检测方法：

通过计算机模拟（in silico）合成双细胞，即将不同单个细胞的读段（reads）混合，生成模拟的双细胞数据。
将合成的双细胞投影到UMAP（一种降维可视化工具）嵌入空间中。
识别合成双细胞的最近邻细胞，判断其与真实细胞的相似性。

######## 双细胞鉴定
doubScores <- addDoubletScores(
  input = ArrowFiles,
  k = 10, #Refers to how many cells near a "pseudo-doublet" to count.
  knnMethod = "UMAP", #Refers to the embedding to use for nearest neighbor search with doublet projection.
  LSIMethod = 1
)
doubScores

运行完后，在QualityControl/该文件夹中，每个样本都有3个相关的图表：

双细胞富集度（Doublet Enrichments）：衡量模拟双细胞在单细胞周围的富集程度。通过与均匀分布的预期值比较，判断是否存在双细胞聚集；
双细胞评分（Doublet Scores）：通过计算二项式校正后的P值（取负对数）来衡量双细胞的显著性。由于其稳定性较差，通常不用于双细胞识别；
双细胞密度（Doublet Density）：展示模拟双细胞在二维嵌入空间中的分布密度。用于可视化合成双细胞的投影位置，帮助评估双细胞检测的合理性。

如其中一个样本BMMC：

3.创建ArchRProject

ArchRProject的作用：将多个Arrow文件整合到一个项目中，便于统一管理和分析。

重要注意事项：不可合并多个ArchRProject：所有需要分析的样本必须在创建项目时以Arrow文件的形式包含进去，后续无法合并多个项目。

推荐设置：copyArrows = TRUE，用于在后续操作中保留Arrow文件的原始副本，以便将来使用。

####### 创建 ArchRProject
projHeme1 <- ArchRProject(
  ArrowFiles = ArrowFiles, 
  outputDirectory = "HemeTutorial",
  copyArrows = TRUE
)
projHeme1

ArchRProject的初始化属性：

输出目录（outputDirectory）：指定用于保存分析结果和图表的路径；
样本名称（samples）：从Arrow文件中提取的样本标识；
样本数据矩阵（sampleColData）：存储与样本相关的信息；
细胞数据矩阵（cellColData）：存储与每个细胞相关的信息，包括双细胞富集分数和双细胞评分（如果已计算）；
细胞总数：经过双细胞识别和去除后的细胞数量；
中位TSS富集分数（衡量转录起始位点的富集程度）；
中位片段数（每个细胞的DNA片段数量）。

是一个高级的S4对象，数据结构如下：

简单探索：

# which data matrices are available
getAvailableMatrices(projHeme1)

# the metadata
head(projHeme1@cellColData)

# access the cell names associated with each cell:
head(projHeme1$cellNames)

# access the sample names associated with each cell
head(projHeme1$Sample)

# access the TSS Enrichment Scores for each cell:
quantile(projHeme1$TSSEnrichment)

# 取子集
idxSample <- BiocGenerics::which(projHeme1$Sample %in% "scATAC_BMMC_R1")
idxSample
projSubset <- subsetArchRProject(
  ArchRProj = projHeme1,
  cells = projHeme1$cellNames[idxSample],
  outputDirectory = "ArchRSubset",
  dropCells = TRUE,
  force = TRUE
)

# 提取细胞指标的某一列
df <- getCellColData(projHeme1, select = "nFrags")
df

# 多列
df <- getCellColData(projHeme1, select = c("log10(nFrags)", "nFrags - 1"))
df

单独绘制QC指标图：前面在创建 ArrowFiles 的时候，函数一键自动化全给我们绘制并生成了文件，我们也可以自己提取数据进行单独绘制

# TSS_by_Unique_Frags.pdf：log10(unique nuclear fragments) vs TSS enrichment score
df <- getCellColData(projHeme1, select = c("log10(nFrags)", "TSSEnrichment"))
df

p <- ggPoint(
    x = df[,1], 
    y = df[,2], 
    colorDensity = TRUE,
    continuousSet = "sambaNight",
    xlabel = "Log10 Unique Fragments",
    ylabel = "TSS Enrichment",
    xlim = c(log10(500), quantile(df[,1], probs = 0.99)),
    ylim = c(0, quantile(df[,2], probs = 0.99))
) + geom_hline(yintercept = 4, lty = "dashed") + geom_vline(xintercept = 3, lty = "dashed")

p

结果如下：

4.绘制一些统计指标

在整合多个不同样本的数据集时，比较不同样本之间的各种指标是常见的需求。

ArchR提供的绘图工具：

山脊图（Ridge Plots）：用于展示不同组（如样本或聚类）中某个变量的分布情况，通过堆叠的密度曲线来直观比较；
小提琴图（Violin Plots）：结合了箱线图和核密度图，展示数据的分布和概率密度，适合比较不同组之间的数据分布。

绘制 the TSS enrichment scores

山脊图：

p1 <- plotGroups(
  ArchRProj = projHeme1, 
  groupBy = "Sample", 
  colorBy = "cellColData", 
  name = "TSSEnrichment",
  plotAs = "ridges",
  baseSize = 10
)
p1

绘制 the TSS enrichment scores 的小提琴图

p2 <- plotGroups(
  ArchRProj = projHeme1, 
  groupBy = "Sample", 
  colorBy = "cellColData", 
  name = "TSSEnrichment",
  plotAs = "violin",
  alpha = 0.4,
  baseSize = 10,
  addBoxPlot = TRUE,
)
p2

其他类似并保存到一个pdf中：

p3 <- plotGroups(
  ArchRProj = projHeme1, 
  groupBy = "Sample", 
  colorBy = "cellColData", 
  name = "log10(nFrags)",
  plotAs = "ridges",
  baseSize = 10
)

p4 <- plotGroups(
  ArchRProj = projHeme1, 
  groupBy = "Sample", 
  colorBy = "cellColData", 
  name = "log10(nFrags)",
  plotAs = "violin",
  alpha = 0.4,
  baseSize = 10,
  addBoxPlot = TRUE
)

plotPDF(p1,p2,p3,p4, name = "QC-Sample-Statistics.pdf", 
        ArchRProj = projHeme1, addDOC = FALSE, width = 4, height = 4)

绘制 Fragment size distributions

在ATAC-seq实验中，片段大小分布是分析数据质量的一个重要指标。不同样本、细胞类型和批次之间的片段大小分布可能会有很大差异。这种差异是常见的，不一定反映数据质量的差异。

使用plotFragmentSizes()函数可以绘制所有样本的片段大小分布图：

p1 <- plotFragmentSizes(ArchRProj = projHeme1)
p1

绘制 TSS enrichment profiles

TSS富集曲线应该在中心位置显示出一个明显的峰值，并且在中心右侧有一个较小的肩峰，这是由于定位良好的+1核小体所引起的。

一个清晰的中心峰值和较小的右侧肩峰表明数据质量良好，TSS区域的染色质开放性较高；
这种模式反映了转录起始位点附近染色质的典型结构，有助于验证实验的可靠性和数据的质量。

p2 <- plotTSSEnrichment(ArchRProj = projHeme1)
p2

# 保存
plotPDF(p1,p2, name = "QC-Sample-FragSizes-TSSProfile.pdf", ArchRProj = projHeme1, addDOC = FALSE, width = 5, height = 5)

附一个应用：https://pages.10xgenomics.com/rs/446-PBO-704/images/CN_10x_LIT055_RevB_AppNote_EpigeneticRegulationATAC_A4_digital.pdf

实验过程：https://pages.10xgenomics.com/rs/446-PBO-704/images/10x_LIT000039_Product_Sheet_Chromium_Single_Cell_ATAC_Letter_digital.pdf

5.保存ArchRProject对象

使用saveArchRProject()函数保存projHeme1项目，以便在后续章节中使用。

saveArchRProject()的功能：

复制Arrow文件，并在指定目录中保存一个.RDS格式的ArchRProject对象；
提供dropCells参数，用于在复制之前从Arrow文件中移除已从ArchRProject中删除的细胞。

# save
projHeme1 <- saveArchRProject(ArchRProj = projHeme1, outputDirectory = "Save-ProjHeme1", load = TRUE)
projHeme1

list.files(path = "./Save-ProjHeme1")
# load 
projHeme1 <- loadArchRProject(path = "./Save-ProjHeme1")
projHeme1

6.过滤Doublets

在不移除真正的单细胞的情况下，尽可能多地过滤掉双细胞。通过绘制DoubletEnrichment或DoubletScore的分布，识别出高分或高富集度的潜在双细胞群体。

双细胞过滤流程：

使用addDoubletScores()添加双细胞预测分数；
使用filterDoublets()移除预测的双细胞。

关键参数：

filterRatio：

用于确定基于通过过滤的细胞数量的最大双细胞过滤比例。
公式为filterRatio * (细胞数量)^2 / 100000。
适用于不同样本，即使它们的双细胞比例因细胞加载浓度不同而不同。
值越高，过滤掉的细胞越多。

cutEnrich：

表示模拟双细胞作为细胞最近邻的比例。
是双细胞富集度的阈值。

cutScore：

表示双细胞富集的-log10(pval)。
比cutEnrich更不准确，通常不推荐作为主要过滤依据。

前面已经运行过：addDoubletScores()

注意选择合适的阈值：

# 过滤
projHeme2 <- filterDoublets(projHeme1,filterRatio = 1.5)

# Filtering 614 cells from ArchRProject!
#   scATAC_BMMC_R1 : 364 of 4932 (7.4%)
#   scATAC_CD34_BMMC_R1 : 160 of 3275 (4.9%)
#   scATAC_PBMC_R1 : 90 of 2453 (3.7%)
未完待续~