scATAC分析神器ArchR初探-简介（1）

简介

近来单细胞技术发展迅速，不仅scRNA技术应用逐渐成熟，scATAC技术也在开始规模应用，分析工具也多种多样，这里介绍由ATAC发明者greenleaf实验室开发的单细胞ATAC分析工具ArchR使用说明，ArchR：对单细胞染色质可及性数据的可靠且可扩展的分析，ArchR是功能齐全的软件套件，用于分析单细胞染色质可访问性数据。它设计用于处理成千上万个单细胞，而没有大的内存或计算需求，与10x Genomics Chromium系统等商业平台可达到的实验规模保持同步。这部分介绍如何将数据导入ArchR，以及如何创建ArrowFiles（ArchR分析的基本单元）

1.1关于ATAC-seq术语的简要入门

任何ATAC-seq实验的最基本组成部分是“ 片段”。在ATAC-seq中，片段是指通过两次转座事件产生的可测序DNA分子。使用配对末端测序对片段的每个末端进行测序。基于Tn5的插入偏移量调整推断的片段起点和终点的单碱基位置。如先前报道，Tn5转座酶以同源二聚体的形式与DNA结合，两个Tn5分子之间具有9 bp的DNA。因此，每个Tn5同型二聚体结合事件都会产生两个插入，间隔为9 bp。因此，“可及”位点的实际中心点位于Tn5二聚体的正中心，而不是每个Tn5插入的位置。为了解决这个问题，我们将偏移量应用于单个Tn5插入，将正链插入事件调整为+4 bp，将负链插入事件调整为-5 bp。ATAC-seq的原始描述。因此，在ArchR中，“ 片段 ”是指一个表或基因组范围对象，其中包含染色体，偏移量调整后的单碱基染色体起始位置，偏移量调整后的单碱基染色体终点位置以及与每个序列化片段相对应的唯一细胞条形码ID 。类似地，“ 插入 ”是指在可到达部位的正中央的偏移调整后的单碱基位置。

1.2为什么要使用ArchR？

那里有多种用于单细胞ATAC序列分析的工具，那么为什么要使用ArchR？最重要的是，ArchR提供功能并支持其他工具无法做到的分析：

此外，由于对数据结构和并行化方法进行了重大优化，构成了ArchR软件的基础，因此与其他可用工具相比，ArchR更快并且使用的内存更少。当分析超过70,000个单元时，某些工具需要高性能的计算环境，超过128 GB的可用内存（OoM =内存不足）。

ArchR被设计用于基于Unix的笔记本电脑上。对于中等规模的实验（少于100,000个细胞），ArchR足够快，可以执行临时分析并实时显示结果，从而可以更深入，更有意义的方式与数据进行交互。当然，对于更高的单元数或更喜欢基于服务器的分析的用户，ArchR提供了可轻松导出的图和项目，可在服务器上生成后下载并使用。
当前，尚未将ArchR优化为在Windows上运行。它应该可以工作，但是尚未为Windows启用ArchR中的并行化功能，因此上述性能提升将无法转化。

1.3什么是箭头文件/ ArchRProject？

ArchR中分析项目的基本单位称为箭头文件。每个Arrow文件都存储与单个样本关联的所有数据（即，元数据，可访问的片段和数据矩阵）。在此，“单个样本”将是所需的最详细的分析单位（例如，特定条件的单个重复）。在创建过程中以及执行附加分析时，ArchR更新和编辑每个Arrow文件以包含附加信息层。值得注意的是，对于ArchR，Arrow文件实际上只是磁盘上存储的外部文件的路径。更明确地说，箭头文件不是存储在内存中的R语言对象，而是存储在磁盘上的HDF5格式的文件。因此，我们使用ArchRProject对象，将这些Arrow文件关联在一起，成为一个可以在R中快速访问的单个分析框架。此ArchRProject对象尺寸很小，并存储在内存中。

某些操作可以直接在Arrow文件上执行，而其他操作则可以在上执行ArchRProject，从而依次更新每个关联的Arrow文件。由于Arrow文件存储为大型HDF5格式文件，因此ArchR中的“ get-er”功能通过与ArchRProjectwhile“ add-er”功能进行交互来检索数据，或者（i）将数据直接添加到Arrow文件，（ii）直接添加数据到ArchRProject，或（iii）通过与交互数据添加到箭头文件ArchRProject。

1.4 ArchR中的输入文件类型

ArchR可以利用scATAC-seq数据的多种输入格式，这种格式最常见的是片段文件和BAM文件的格式。片段文件是由tabix排序的文本文件，包含每个scATAC-seq片段和相应的单元ID，每行一个片段。BAM文件是按Tabix排序的二进制文件，包含每个scATAC-seq片段，原始序列，蜂窝条形码ID和其他信息。使用的输入格式将取决于使用的预处理管道。例如，10x Genomics Cell Ranger软件返回片段文件，而sci-ATAC-seq应用程序将使用BAM文件。ArchR使用“ scanTabix”读取片段文件，使用“ scanBam”读取BAM文件。在此输入过程中，对输入进行分块，然后将每个输入块转换为基于压缩表的片段表示形式，其中包含每个片段染色体，偏移校正后的染色体起始位置，偏移校正后的染色体终止位置和细胞条形码ID。然后，将这些逐块片段存储在HDF5格式的临时文件中，以保留内存使用率，同时保持对每个块的快速访问。最后，与每个染色体相关的所有块均被读取，组织并重新写入一个称为“片段”的HDF5组内的Arrow文件。此预分块过程使ArchR可以有效地处理非常大的输入文件，并且内存使用率较低，从而可以充分利用并行处理。与每个染色体相关的所有块均被读取，组织并重写为一个称为“片段”的HDF5组中的Arrow文件。此预分块过程使ArchR可以有效地处理非常大的输入文件，并且内存使用率较低，从而可以充分利用并行处理。与每个染色体相关的所有块均被读取，组织并重写为一个称为“片段”的HDF5组中的Arrow文件。此预分块过程使ArchR可以有效地处理非常大的输入文件，并且内存使用率较低，从而可以充分利用并行处理。

1.5设置

我们要做的第一件事是更改到所需的工作目录，设置我们要使用的线程数，并加载我们的基因和基因组注释。根据您本地环境的配置，您可能需要在中修改threads以下使用的数量addArchRThreads()。默认情况下，ArchR使用threads可用总数的一半，但是您可以根据需要手动进行调整。如果使用Windows，则threads因为ArchR中的并行处理是针对基于Unix的操作系统而构建的，所以可用将自动设置为1。

首先，我们加载ArchR库。如果失败，则说明您尚未正确安装ArchR，应重新访问安装说明。

library(ArchR)

接下来，我们设置ArchR函数的默认线程数。这是您在每个新的R会话期间都必须执行的操作。我们建议将threads可用内核总数设置为1/2到3/4。ArchR中的内存使用量通常会随所用线程的数量而定，因此允许ArchR使用更多线程也将导致更高的内存使用量。

addArchRThreads(threads = 16)

将并行线程的默认数量设置为16。

然后，我们设置用于基因和基因组注释的基因组。如上所述，这是您在每个新的R会话期间都必须执行的操作。当然，该基因组版本必须与用于比对的基因组版本匹配。对于本教程中使用的数据，我们将使用hg19参考基因组，但是ArchR本机支持其他基因组注释和自定义基因组注释，如下一节所述。

addArchRGenome("hg19")

将默认基因组设置为Hg19。##设置基因组和基因注释

ArchR需要基因和基因组注释来执行诸如计算TSS富集得分，核苷酸含量和基因活性得分之类的事情。由于我们的教程数据集使用已经与hg19参考基因组对齐的scATAC-seq数据，因此在上一节中我们将“ hg19”设置为默认基因组。但是，ArchR本地支持“ hg19”，“ hg38”，“ mm9”和“ mm10”，您可以使用createGeneAnnotation()和createGenomeAnnotation()函数创建自己的基因组和基因注释。

通过该addArchRGenome()功能可以简化向ArchR提供此信息的过程。此函数告诉ArchR，对于当前会话中的所有分析，它都应使用genomeAnnotation和geneAnnotation与define 关联ArchRGenome。每个由本地支持的基因组均由BSgenome定义每个染色体的基因组坐标和序列的GRanges对象，包含一组列入黑名单的区域的TxDb对象，定义所有基因的位置和结构的OrgDb对象以及提供中心基因的对象组成基因标识符，并包含此标识符与其他种类的标识符之间的映射。

以下是如何为本地支持的基因组加载基因和基因组注释的示例，以及有关其BSgenome和黑名单组件的信息。

在预编译的版本hg19基因组ArchR使用BSgenome.Hsapiens.UCSC.hg19，TxDb.Hsapiens.UCSC.hg19.knownGene，org.Hs.eg.db，和用合并黑名单ArchR::mergeGR()从hg19 V2黑名单区和显示高mappability到hg19核基因组线粒体地区从迦勒Lareau和Jason Buenrostro。要将全局基因组默认设置为预编译的hg19基因组，请执行以下操作：

addArchRGenome("hg19")

将默认基因组设置为Hg19。

在预编译的版本hg38基因组ArchR使用BSgenome.Hsapiens.UCSC.hg38，TxDb.Hsapiens.UCSC.hg38.knownGene，org.Hs.eg.db，和用合并黑名单ArchR::mergeGR()从hg38 V2黑名单区和显示高mappability到hg38核基因组线粒体地区从迦勒Lareau和Jason Buenrostro。要将全局基因组默认设置为预编译的hg38基因组，请执行以下操作：

addArchRGenome("hg38")

将默认基因组设置为Hg38。

在预编译的版本MM9在ArchR用途的基因组BSgenome.Mmusculus.UCSC.mm9，TxDb.Mmusculus.UCSC.mm9.knownGene，org.Mm.eg.db，并使用合并，这是一个黑名单ArchR::mergeGR()从MM9 V1黑名单地区从Anshul Kundaje和显示高mappability至MM9核基因组线粒体地区从迦勒Lareau和Jason Buenrostro。要将全局基因组默认设置为预编译的mm9基因组，请执行以下操作：

addArchRGenome("mm9")

将默认基因组设置为Mm9。

在预编译的版本MM10基因组ArchR使用BSgenome.Mmusculus.UCSC.mm10，TxDb.Mmusculus.UCSC.mm10.knownGene，org.Mm.eg.db，和用合并黑名单ArchR::mergeGR()从MM10 V2黑名单区和显示高mappability到MM10核基因组线粒体地区从迦勒Lareau和Jason Buenrostro。要将全局基因组默认设置为预编译的mm10基因组，请执行以下操作：

addArchRGenome("mm10")

将默认基因组设置为Mm10。

1.5.1创建自定义`ArchRGenome`

如上所述，ArchRGenome由基因组注释和基因注释组成。

要创建自定义基因组注释，我们使用createGenomeAnnotation()。为此，您将需要以下信息：

甲BSgenome其中包含一个基因组序列信息的对象。这些通常是Bioconductor软件包（例如，BSgenome.Hsapiens.UCSC.hg38），可通过google轻松找到。
一个GRanges基因组范围对象，其中包含一组列入黑名单的区域，这些区域将用于从下游分析中过滤掉不需要的区域。这不是必需的，但建议这样做。有关如何创建黑名单的信息，请参阅ENCODE黑名单上的出版物。

例如，如果我们想从果蝇中创建自定义基因组注释，我们将首先识别并安装并加载相关BSgenome对象。

if (!requireNamespace("BSgenome.Dmelanogaster.UCSC.dm6", quietly = TRUE)){
  BiocManager::install("BSgenome.Dmelanogaster.UCSC.dm6")
}
library(BSgenome.Dmelanogaster.UCSC.dm6)

加载所需的软件包：BSgenome
加载所需的软件包：Biostrings
加载所需的软件包：XVector
附加软件包：'Biostrings'
以下对象从'package：base'中屏蔽：strsplit

然后，我们从该BSgenome对象创建基因组注释。

genomeAnnotation <- createGenomeAnnotation(genome = BSgenome.Dmelanogaster.UCSC.dm6)

正在获取基因组..
正在获取chromSizes ..
正在获取黑名单..
未下载黑名单！如果没有继续，请注意下游的偏差。

检查该对象显示了ArchR基因组注释的组成部分。

genomeAnnotation

长度为3的列表
名称（3）：基因组chromSizes黑名单

要创建自定义基因注释，我们使用createGeneAnnotation()。为此，您将需要以下信息：

甲TxDb从Bioconductor的对象（转录数据库），其包含用于基因/转录物的坐标的信息。
OrgDb来自Bioconductor 的对象（生物数据库），提供了一个统一的框架来映射基因名称和各种基因标识符。

继续以果蝇（Drosophila melanogaster）为例，我们首先安装并加载相关对象TxDb和OrgDb对象。

if (!requireNamespace("TxDb.Dmelanogaster.UCSC.dm6.ensGene", quietly = TRUE)){
  BiocManager::install("TxDb.Dmelanogaster.UCSC.dm6.ensGene")
}
if (!requireNamespace("org.Dm.eg.db", quietly = TRUE)){
  BiocManager::install("org.Dm.eg.db")
}
library(TxDb.Dmelanogaster.UCSC.dm6.ensGene)

加载所需的软件包：GenomicFeatures
加载所需的软件包：AnnotationDbi

library(org.Dm.eg.db)

然后我们创建基因注释对象。

geneAnnotation <- createGeneAnnotation(TxDb = TxDb.Dmelanogaster.UCSC.dm6.ensGene, OrgDb = org.Dm.eg.db)

正在获取基因..
确定的注释样式= ENSEMBL
正在获取外显子..
正在获取TSS ..

检查该对象显示了ArchR基因注释的组成部分。

geneAnnotation

长度为3的列表
名称（3）：基因外显子TSS

另外，如果您没有TxDb和OrgDb对象，则可以geneAnnotation根据以下信息创建一个对象：

“基因”对象- GRanges包含基因坐标的对象（开始到结束）。该对象的符号列必须与下面描述的“外显子”对象的符号列匹配。
“外显子”对象- GRanges包含基因外显子坐标的对象。必须具有与上述“基因”对象的符号列匹配的符号列。
一个GRanges包含standed转录起始位点（对象TSS）的坐标。

geneAnnotation <- createGeneAnnotation(
  TSS = geneAnnotation$TSS, 
  exons = geneAnnotation$exons, 
  genes = geneAnnotation$genes
)

1.5.2在ArchR中使用非标准基因组

ArchR实施了一些检查，以防止您执行我们认为“超出正常范围”的事情。这些检查之一会强制seqnames您的基因组注释以“ chr”开头。在大多数应用中都是如此，但是有些基因组（例如玉米）不使用“ chr”作为不同染色体的前缀。要对此类基因组进行任何ArchR分析，您必须告诉ArchR忽略染色体前缀。为此，您必须addArchRChrPrefix(chrPrefix = FALSE)在创建Arrow文件之前运行。这将在当前R会话中全局关闭对seqnames的“ chr”前缀的要求。请记住，ArchR默认将染色体/序列名称转换为字符。因此，如果你seqnames仅仅是整数，你需要提供他们的角色，即 useSeqnames = c("1", "2", "3")替代useSeqnames = c(1, 2, 3)。您始终可以使用来检查当前R会话中是否需要染色体前缀getArchRChrPrefix()。

1.6创建箭头文件

在本教程的其余部分中，我们将使用来自Granja *等人的造血细胞降采样数据集的数据。自然生物技术2019。这包括来自骨髓单核细胞（BMMC），外周血单核细胞（PBMC）和来自骨髓的CD34 +造血干细胞和祖细胞（CD34 BMMC）的数据。

该数据以片段文件的形式下载，其中包含所有比对的测序片段的起始和终止基因组坐标。片段文件是10x Genomics分析平台（和其他平台）的基本文件类型之一，可以从任何BAM文件轻松创建。有关创建自己的片段文件以输入到ArchR的信息，请参见10x Genomics网站。

有了片段文件后，我们会将它们的路径作为字符向量提供给createArrowFiles()。在创建过程中，一些基本的元数据和矩阵被添加到每个Arrow文件中，其中包括一个“ TileMatrix”（包含跨全基因组范围500 bp的条带的插入计数addTileMatrix()）（请参见参考资料）和“ GeneScoreMatrix”，该存储基于在图块中加权插入计数的预测基因表达在基因启动子附近（请参阅参考资料addGeneScoreMatrix()）。

可以使用该getTutorialData()功能下载教程数据。教程数据的大小约为0.5 GB。如果您已经在当前工作目录中下载了该教程，则ArchR将绕过下载。

library(ArchR)

inputFiles <- getTutorialData("Hematopoiesis")
inputFiles

scATAC_BMMC_R1
“HemeFragments / scATAC_BMMC_R1.fragments.tsv.gz”
scATAC_CD34_BMMC_R1
“HemeFragments / scATAC_CD34_BMMC_R1.fragments.tsv.gz”
scATAC_PBMC_R1
“HemeFragments / scATAC_PBMC_R1.fragments.tsv.gz”

与往常一样，在开始项目之前，我们必须为并行化设置ArchRGenome和默认值threads。

addArchRGenome("hg19")

将默认基因组设置为Hg19。

addArchRThreads(threads = 16)

将并行线程的默认数量设置为16。

现在，我们将创建耗时10-15分钟的箭头文件。对于每个样本，此步骤将：

从提供的输入文件中读取可访问的片段。
计算每个细胞的质量控制信息（即TSS富集得分和核小体信息）。
根据质量控制参数过滤单元。
使用500 bp的条带创建全基因组的TileMatrix。
使用geneAnnotation调用时定义的自定义创建GeneScoreMatrix addArchRGenome()。

ArrowFiles <- createArrowFiles(
  inputFiles = inputFiles,
  sampleNames = names(inputFiles),
  filterTSS = 4, #Dont set this too high because you can always increase later
  filterFrags = 1000, 
  addTileMat = TRUE,
  addGeneScoreMat = TRUE
)

使用addArchRGenome（Hg19）设置的GeneAnnotation！
使用addArchRGenome（Hg19）设置的GeneAnnotation！
ArchR记录到：ArchRLogs / ArchR-createArrows-dfa159ddbf6e-Date-2020-04-15_Time-09-21-27.log
如果有问题，请使用logFile向github报告！
清理临时文件
2020-04-15 09:21:28：批量执行w / safelapply !，已过去0分钟。
ArchR记录成功到：ArchRLogs / ArchR-createArrows-dfa159ddbf6e-Date-2020-04-15_Time-09-21-27.log

我们可以检查ArrowFiles对象，以查看它实际上只是Arrow文件路径的字符向量。

ArrowFiles

“ scATAC_BMMC_R1.arrow”“ scATAC_CD34_BMMC_R1.arrow”
“ scATAC_PBMC_R1.arrow”

1.7每单元质量控制

scATAC-seq数据的严格质量控制（QC）对于消除低质量细胞的作用至关重要。在ArchR中，我们考虑数据的三个特征：

独特核片段的数量（即不映射到线粒体DNA）。
信噪比。低的信噪比通常归因于死亡或垂死的细胞，这些细胞具有脱色的DNA，可以在全基因组范围内进行随机转座。
片段大小分布。由于核小体的周期性，我们预计会看到碎片消失，这些碎片是包裹在核小体周围的DNA长度（大约147 bp）。

第一个度量标准是唯一的核片段，很简单-可用片段很少的细胞将无法提供足够的数据来做出有用的解释，因此应排除在外。

第二个指标，即信噪比，被计算为TSS富集得分。传统的批量ATAC-seq分析已使用此TSS富集得分作为确定信号到背景（例如，ENCODE项目）的标准工作流程的一部分）。我们和其他人已经发现，在批量ATAC-seq和scATAC-seq中测试的大多数细胞类型中，TSS富集具有代表性。TSS富集得分度量标准背后的想法是，由于与启动子结合的大型蛋白质复合物，与其他基因组区域相比，ATAC-seq数据普遍在基因TSS区域富集。通过查看以这些TSS区域为中心的每个碱基对的可及性，我们发现相对于侧翼区域（两个方向的远端1900-2000 bp）而言，局部富集。相对于这些侧翼区域，该富集峰（以TSS为中心）之间的比率表示TSS富集得分。

传统上，针对每个批量ATAC-seq样本计算每个碱基对的可访问性，然后使用此配置文件确定TSS富集得分。在scATAC-seq中以每个单元为单位执行此操作相对较慢且计算量很大。为了准确估算每个单细胞的TSS富集得分，我们计算以每个单碱基TSS位置为中心的50 bp区域内的平均可访问性，并将其除以TSS侧翼位置的平均可访问性（+/- 1900 – 2000 bp ）。该近似值与原始方法高度相关（R> 0.99），并且值在大小上非常接近。

第三个度量标准，片段大小分布通常不那么重要，但始终可以手动检查。由于DNA包裹核小体的模式化方式，我们希望在我们的数据中看到片段大小分布的核小体周期性。之所以出现这些丘陵和山谷，是因为片段必须跨越0、1、2等。核小体（Tn5无法切割紧紧包裹在核小体周围的DNA。

默认情况下，在ArchR中，将通过过滤器细胞识别为那些TSS富集得分大于4且具有1000个以上独特核片段的细胞。重要的是要注意，TSS富集得分的实际数值取决于所使用的TSS的集合。ArchR中的默认值是为人类数据设计的，在运行时更改默认阈值可能很重要createArrowFiles()。

创建Arrow文件将在当前工作目录中创建一个名为“ QualityControl”的文件夹，该文件夹将包含与每个样本关联的2个图。第一张图显示了log10(unique nuclear fragments)vs TSS富集得分，并指出了虚线所用的阈值。第二个显示片段大小分布。

对于我们的教程数据，我们有三个示例，如下所示：

对于BMMC：

对于CD34 BMMC：

对于PBMC：

参考材料：

https://www.archrproject.com/