上一章我们介绍了有关于scATAC-seq的相关内容，包括scATAC-seq的优缺点、分析软件已经分析方法的介绍。从本章开始，开始介绍使用ArchR来进行scATAC-seq下游分析的代码以及相关算法。有兴趣的朋友可以官网查看官方操作手册。
1 ArchR对象
首先，我们需要介绍的就是ArchRProject，也就是ArchR对象，在上一章中我们提到，ArchR对象是用来关联Arrow文件的，除此之外，ArchR还储存了一些其他信息，下面我们分别进行介绍。
1.1 输入文件
我们需要建立先建立一个ArchRProject。ArchR的输入文件很多样化，最常见的输入文件格式有两种，分别是fragment文件和BAM文件。Fragment文件是一个经过tabix排序的文本文件，它包含了每一个scATAC-seq的fragment信息以及相对应的细胞ID，其中每一行为一个fragment。BAM文件是一个经过tabix排序的二进制文件，它包含了每一个scATAC-seq的fragment信息、原始序列信息、细胞的barcode ID和其他的一些信息。这里，我们使用的数据是官网提供的数据。输入文件是运行完Cellranger-atac之后结果文件（fragment文件和BAM文件任选其一）。首先，我们需要基于这些数据文件创建Arrow文件，并且设置一些分析用到的参数。
library(ArchR)
addArchRThreads(threads = 3)
addArchRGenome("hg38")
#Create Arrow files
#/path 是存放fragment文件或者bam文件的路径，如果是bam文件，请将gz$改为bam$
input <- list.files(path = "/path", pattern = "gz$")
#为输入的文件添加名字
names <- unlist(lapply(input , function(x){
name <- unlist(strsplit(x, "[.]"))
name <- name[1]
}))
input <- unlist(lapply(input , function(x){
x <- file.path("/path", x)
}))
names(input ) <- names
input
scATAC_BMMC_R1
"/path/scATAC_BMMC_R1.fragments.tsv.gz"
scATAC_CD34_BMMC_R1
"/path/scATAC_CD34_BMMC_R1.fragments.tsv.gz"
scATAC_PBMC_R1
"/path/scATAC_PBMC_R1.fragments.tsv.gz"
#Create Arrow files
ArrowFiles <- createArrowFiles(inputFiles=input,
sampleNames=names(input),
minTSS=4,
minFrags=1000,
addTileMat=TRUE,
addGeneScoreMat=TRUE)
1.2 细胞质量过滤
去除一些低质量的细胞有助于提高数据的准确性。在ArchR中，作者考虑了3个方面的质量特征：
- 唯一fragment片段的数量（不包含比对到线粒体DNA上的片段）；
- single-to-background比率。低的single-to-background比率通常是有凋亡或者正在凋亡的细胞引起的，这种状态的下的细胞，它们的脱氧核苷酸可以在全基因组范围内进行随机转位。该数值也是后面分析时提到的TSS富集值，我们在观察TSS中心区域可及性的时候发现，相较于两侧（距离中心1900-2000bp的区域）的情况，中心区域会出现富集的情况，这些富集峰与两侧区域的比值代表了TSS的富集分数；
- fragment片段大小分布。由于核小体的周期性，我们希望看到包裹核小体的DNA的在长度上的逐渐消逝（大概在147bp的位置）。
1.3 计算双细胞值
双细胞数据是由于在微流控过程中产生的一个油包水中包含两个细胞的情况，这样的数据会影响我们数据的准确性，因此需要将这部分数据进行过滤，提高数据的质量。
doubscore <- addDoubletScores(input=ArrowFiles,
k=10,
knnMethod='UMAP',
LSIMethod=1,
threads=1)
## If there is an issue, please report to github with logFile!
## 2020-04-15 09:28:44 : Batch Execution w/ safelapply!, 0 mins elapsed.
## 2020-04-15 09:28:44 : scATAC_BMMC_R1 (1 of 3) : Computing Doublet Statistics, 0.001 mins elapsed.
## scATAC_BMMC_R1 (1 of 3) : UMAP Projection R^2 = 0.9736
## scATAC_BMMC_R1 (1 of 3) : UMAP Projection R^2 = 0.9736
## 2020-04-15 09:31:15 : scATAC_CD34_BMMC_R1 (2 of 3) : Computing Doublet Statistics, 2.511 mins elapsed.
## scATAC_CD34_BMMC_R1 (2 of 3) : UMAP Projection R^2 = 0.99046
## scATAC_CD34_BMMC_R1 (2 of 3) : UMAP Projection R^2 = 0.99046
## 2020-04-15 09:32:40 : scATAC_PBMC_R1 (3 of 3) : Computing Doublet Statistics, 3.936 mins elapsed.
## scATAC_PBMC_R1 (3 of 3) : UMAP Projection R^2 = 0.97507
## scATAC_PBMC_R1 (3 of 3) : UMAP Projection R^2 = 0.97507
在计算双细胞值的过程中会产生一个R^2的值，如果该值小于0.9则表明数据中包含有少量的异质性，这会导致我们对双细胞预测的准确性降低，因此需要停止双细胞的预测与过滤。
1.4 创建ArchR对象
ArchR对象和Seurat对象一样是S4类的格式，想要理解S4类请查看官方文档。里面储存了关于分析样本的各种信息，比如Arrow文件的路径、细胞分组信息、基因信息、peak信息、降维聚类信息等等。是进行后续分析的一个基础，所有的分析都是围绕ArchR对象展开的。
#创建ArchRProject
archr_obj <- ArchRProject(
ArrowFiles = ArrowFiles,
outputDirectory = output_dir,
copyArrows = TRUE
)
#基于doubscore过滤双细胞，如果R^2 < 0.9，这一步可以跳过
archr_obj <- filterDoublets(archr_obj, filterRatio = filterRatio)
创建完ArchR对象以及过滤掉双细胞数据以后，可以观察一下ATAC-seq常用的指标来观察一下我们得到的数据的质量：比如Fragment大小分布和TSS富集情况：
fsd <- plotFragmentSizes(archr_obj)
fsd
tep <- plotTSSEnrichment(archr_obj, groupBy="Sample")
tep
在富集谱中正常情况下可以看到在中心区域有一个比较明显的峰，在这个峰的右边大约147bp的位置还会有一个由于核小体影响而形成的小峰。
