上一章我們介紹了有關(guān)于scATAC-seq的相關(guān)內(nèi)容，包括scATAC-seq的優(yōu)缺點(diǎn)、分析軟件已經(jīng)分析方法的介紹。從本章開(kāi)始，開(kāi)始介紹使用ArchR來(lái)進(jìn)行scATAC-seq下游分析的代碼以及相關(guān)算法。有興趣的朋友可以官網(wǎng)查看官方操作手冊(cè)。

1 ArchR對(duì)象

首先，我們需要介紹的就是ArchRProject，也就是ArchR對(duì)象，在上一章中我們提到，ArchR對(duì)象是用來(lái)關(guān)聯(lián)Arrow文件的，除此之外，ArchR還儲(chǔ)存了一些其他信息，下面我們分別進(jìn)行介紹。

1.1 輸入文件

我們需要建立先建立一個(gè)ArchRProject。ArchR的輸入文件很多樣化，最常見(jiàn)的輸入文件格式有兩種，分別是fragment文件和BAM文件。Fragment文件是一個(gè)經(jīng)過(guò)tabix排序的文本文件，它包含了每一個(gè)scATAC-seq的fragment信息以及相對(duì)應(yīng)的細(xì)胞ID，其中每一行為一個(gè)fragment。BAM文件是一個(gè)經(jīng)過(guò)tabix排序的二進(jìn)制文件，它包含了每一個(gè)scATAC-seq的fragment信息、原始序列信息、細(xì)胞的barcode ID和其他的一些信息。這里，我們使用的數(shù)據(jù)是官網(wǎng)提供的數(shù)據(jù)。輸入文件是運(yùn)行完Cellranger-atac之后結(jié)果文件（fragment文件和BAM文件任選其一）。首先，我們需要基于這些數(shù)據(jù)文件創(chuàng)建Arrow文件，并且設(shè)置一些分析用到的參數(shù)。

library(ArchR)
addArchRThreads(threads = 3)
addArchRGenome("hg38")
#Create Arrow files
#/path 是存放fragment文件或者bam文件的路徑，如果是bam文件，請(qǐng)將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 細(xì)胞質(zhì)量過(guò)濾

去除一些低質(zhì)量的細(xì)胞有助于提高數(shù)據(jù)的準(zhǔn)確性。在ArchR中，作者考慮了3個(gè)方面的質(zhì)量特征：

• 唯一fragment片段的數(shù)量（不包含比對(duì)到線粒體DNA上的片段）；
• single-to-background比率。低的single-to-background比率通常是有凋亡或者正在凋亡的細(xì)胞引起的，這種狀態(tài)的下的細(xì)胞，它們的脫氧核苷酸可以在全基因組范圍內(nèi)進(jìn)行隨機(jī)轉(zhuǎn)位。該數(shù)值也是后面分析時(shí)提到的TSS富集值，我們?cè)谟^察TSS中心區(qū)域可及性的時(shí)候發(fā)現(xiàn)，相較于兩側(cè)（距離中心1900-2000bp的區(qū)域）的情況，中心區(qū)域會(huì)出現(xiàn)富集的情況，這些富集峰與兩側(cè)區(qū)域的比值代表了TSS的富集分?jǐn)?shù)；
• fragment片段大小分布。由于核小體的周期性，我們希望看到包裹核小體的DNA的在長(zhǎng)度上的逐漸消逝（大概在147bp的位置）。

1.3 計(jì)算雙細(xì)胞值

雙細(xì)胞數(shù)據(jù)是由于在微流控過(guò)程中產(chǎn)生的一個(gè)油包水中包含兩個(gè)細(xì)胞的情況，這樣的數(shù)據(jù)會(huì)影響我們數(shù)據(jù)的準(zhǔn)確性，因此需要將這部分?jǐn)?shù)據(jù)進(jìn)行過(guò)濾，提高數(shù)據(jù)的質(zhì)量。

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

在計(jì)算雙細(xì)胞值的過(guò)程中會(huì)產(chǎn)生一個(gè)R^2的值，如果該值小于0.9則表明數(shù)據(jù)中包含有少量的異質(zhì)性，這會(huì)導(dǎo)致我們對(duì)雙細(xì)胞預(yù)測(cè)的準(zhǔn)確性降低，因此需要停止雙細(xì)胞的預(yù)測(cè)與過(guò)濾。

1.4 創(chuàng)建ArchR對(duì)象

ArchR對(duì)象和Seurat對(duì)象一樣是S4類的格式，想要理解S4類請(qǐng)查看官方文檔。里面儲(chǔ)存了關(guān)于分析樣本的各種信息，比如Arrow文件的路徑、細(xì)胞分組信息、基因信息、peak信息、降維聚類信息等等。是進(jìn)行后續(xù)分析的一個(gè)基礎(chǔ)，所有的分析都是圍繞ArchR對(duì)象展開(kāi)的。

#創(chuàng)建ArchRProject
archr_obj <- ArchRProject(
  ArrowFiles = ArrowFiles,
  outputDirectory = output_dir,
  copyArrows = TRUE
)
#基于doubscore過(guò)濾雙細(xì)胞，如果R^2 < 0.9，這一步可以跳過(guò)
archr_obj <- filterDoublets(archr_obj, filterRatio = filterRatio)

創(chuàng)建完ArchR對(duì)象以及過(guò)濾掉雙細(xì)胞數(shù)據(jù)以后，可以觀察一下ATAC-seq常用的指標(biāo)來(lái)觀察一下我們得到的數(shù)據(jù)的質(zhì)量：比如Fragment大小分布和TSS富集情況：

fsd <- plotFragmentSizes(archr_obj)
fsd

Fragment大小分布圖

tep <- plotTSSEnrichment(archr_obj, groupBy="Sample")
tep

TSS富集譜

在富集譜中正常情況下可以看到在中心區(qū)域有一個(gè)比較明顯的峰，在這個(gè)峰的右邊大約147bp的位置還會(huì)有一個(gè)由于核小體影響而形成的小峰。