前面的學(xué)習(xí)筆記里，對(duì)于scATAC-seq的分析，使用的是signac R包進(jìn)行的。前幾天在一次偶然的機(jī)會(huì)里，有幸和一個(gè)專門做生信的老外聊了聊（紐約某校的PI），從他那里聽說ArchR是一個(gè)非常強(qiáng)大且值得學(xué)習(xí)的軟件。那么既然專家都這么說了，那必須趕快學(xué)起來。我知道網(wǎng)上有關(guān)于ArchR的系列教程已經(jīng)有人寫過了，但還是決定自己跟著英文教程過一遍，加深印象~
ArchR是一個(gè)用于分析單細(xì)胞染色質(zhì)可接近性數(shù)據(jù)（scATAC-seq）的功能齊全的R包套件。它能夠用于處理數(shù)十萬個(gè)單細(xì)胞，無需很大的內(nèi)存或強(qiáng)大的電腦配置，與商業(yè)平臺(tái)如10x Genomics Chromium系統(tǒng)可實(shí)現(xiàn)的實(shí)驗(yàn)規(guī)模保持同步。
這個(gè)官網(wǎng)（here）是一個(gè)完整的用戶使用手冊(cè)，是今年5月才更新過的。作者通過所能想到的所有例子來解釋如何使用ArchR。還包括了一些關(guān)于scATAC -seq分析過程中重點(diǎn)的文檔，例如降維和聚類。所有的章節(jié)都將使用相同造血細(xì)胞數(shù)據(jù)集Granja* et al. Nature Biotechnology 2019。每一章都將建立在前面的基礎(chǔ)上，所以理想情況下你應(yīng)該從頭開始并按順序執(zhí)行所有的分析。需要說明的是，這個(gè)教程是在本地運(yùn)行ArchR。（非服務(wù)器）
（一）ArchR的安裝
安裝ArchR:
# install ArchR
# make sure you already installed devtools and BiocManager
devtools::install_github("GreenleafLab/ArchR", ref="master", repos = BiocManager::repositories())
然后再安裝ArchR的依賴包：
# install all of the ArchR dependencies
library(ArchR)
ArchR::installExtraPackages()
如果安裝有問題，請(qǐng)參考：https://github.com/GreenleafLab/ArchR
（二）ArchR使用的初次了解
這一章將介紹如何將數(shù)據(jù)導(dǎo)入到ArchR里，并講解如何創(chuàng)建ArrowFiles，這個(gè)是ArchR分析的基本單位。
（1）ATAC-seq的簡(jiǎn)要介紹
任何ATAC-seq實(shí)驗(yàn)最基本的組成部分都是一個(gè)“片段”。在ATAC-seq中，一個(gè)片段指的是由兩次換位事件產(chǎn)生的可測(cè)序的DNA分子。該片段的每一端都使用paired-end進(jìn)行測(cè)序。根據(jù)Tn5的插入偏移調(diào)整片段開始和結(jié)束的單堿基位置。就像已經(jīng)報(bào)道的文獻(xiàn)說的一樣，Tn5轉(zhuǎn)座酶以一種同源二聚體與DNA結(jié)合，在兩個(gè)Tn5分子之間有9bp的DNA。因此，每個(gè)Tn5同源二聚體結(jié)合事件都會(huì)產(chǎn)生兩個(gè)插入，間隔9bp。因此，“可接近”位點(diǎn)的實(shí)際中心點(diǎn)在Tn5二聚體的正中心，而不是每個(gè)Tn5插入的位置。考慮到這點(diǎn)，我們對(duì)每個(gè)Tn5插入應(yīng)用一個(gè)偏移量，將正鏈插入事件調(diào)整+4 bp，將負(fù)鏈插入事件調(diào)整-5 bp。這與ATAC-seq最初描述時(shí)提出的定義是一致的。因此，在ArchR中，“fragments”指的是一組對(duì)象：包含染色體、偏移調(diào)整的單堿基染色體起始位置、偏移調(diào)整的單堿基染色體結(jié)束位置以及每個(gè)唯一細(xì)胞barcode ID對(duì)應(yīng)的測(cè)序片段。類似地，“insertions”指的是可接近位點(diǎn)最中心的偏移調(diào)整后的單堿基位置。
（2）為什么要選擇ArchR?
目前有很多的scATAC-seq分析軟件可供使用，那么為什么選擇ArchR呢？最主要的是ArchR提供了很多其他軟件沒有的功能：
此外，ArchR比其他可用工具更快，占用的內(nèi)存更少，這是由于對(duì)數(shù)據(jù)結(jié)構(gòu)和并行化方法進(jìn)行了大量?jī)?yōu)化，而這些方法構(gòu)成了ArchR軟件的基礎(chǔ)。當(dāng)分析超過70,000個(gè)細(xì)胞時(shí)，其他工具需要高性能計(jì)算環(huán)境，需要內(nèi)存超過128 GB 。
ArchR最初基于unix的筆記本電腦被設(shè)計(jì)出來的。對(duì)于中等大小的實(shí)驗(yàn)(少于100,000個(gè)細(xì)胞)，ArchR的速度足夠快，執(zhí)行分析并可視化結(jié)果，使其能夠以更深入和更具生物學(xué)意義的方式與數(shù)據(jù)交互。當(dāng)然，對(duì)于較高的細(xì)胞數(shù)或喜歡使用服務(wù)器的用戶，ArchR提供了導(dǎo)出可在服務(wù)器上生成后下載和使用的plot和projects。
目前，ArchR還沒有在Windows上優(yōu)化。它可以工作，但是ArchR中的并行化還沒有在Windows中啟用，所以上面提到的性能還無法轉(zhuǎn)換。
（3）什么是Arrow file / ArchRProject?
在ArchR中，一個(gè)可分析項(xiàng)目的基本單元稱為Arrow file。每個(gè)Arrow file存儲(chǔ)與單個(gè)樣本相關(guān)的所有數(shù)據(jù)(即元數(shù)據(jù)、可接近的片段和數(shù)據(jù)矩陣)。在創(chuàng)建期間和執(zhí)行其他分析時(shí)，ArchR更新并編輯每個(gè)Arrow file，以包含額外的信息層。值得注意的是，對(duì)于ArchR來說，Arrow file實(shí)際上只是存儲(chǔ)在磁盤上的外部文件的路徑。更明確地說，Arrow file不是存儲(chǔ)在內(nèi)存中的R語言對(duì)象，而是存儲(chǔ)在磁盤上的HDF5格式文件。因此，我們使用ArchRProject對(duì)象將這些Arrow file關(guān)聯(lián)到一個(gè)單獨(dú)的分析框架中，該框架可以在R中被快速訪問。這個(gè)ArchRProject對(duì)象很小，并且存儲(chǔ)在內(nèi)存中。
某些操作可以直接在Arrow files上執(zhí)行，而有些操作則需要在ArchRProject上執(zhí)行，反過來再更新每個(gè)相關(guān)的Arrow files。因?yàn)锳rrow files存儲(chǔ)為一個(gè)大的hdf5格式文件，ArchR中的get-er函數(shù)通過與ArchRProject交互來獲取數(shù)據(jù)，而add-er函數(shù)可以(i)直接向Arrow files添加數(shù)據(jù)，(ii)直接向ArchRProject添加數(shù)據(jù)，或者(iii)通過與ArchRProject交互向Arrow files添加數(shù)據(jù)。
（4）ArchR的輸入文件格式
ArchR可以利用scATAC-seq數(shù)據(jù)的多種輸入格式，最常見的格式是fragment文件和BAM文件。
fragment文件是tabix排序的文本文件，包含每個(gè)scATAC-seq片段和相應(yīng)的細(xì)胞ID，每行一個(gè)片段。
BAM文件是二進(jìn)制的tabix排序文件，包含每個(gè)scATAC-seq片段、原始序列、細(xì)胞barcode id和其他信息。
所使用的輸入格式將取決于所使用的預(yù)處理流程使用的格式。例如，10x Genomics Cell Ranger軟件返回fragment文件，而sci-ATAC-seq應(yīng)用程序?qū)⑹褂肂AM文件。ArchR使用“scanTabix”讀取fragment文件，使用“scanBam”讀取BAM文件。在此輸入過程中，inputs被切分成chunk，每個(gè)input chunk 轉(zhuǎn)換為壓縮表格，包含每個(gè)片段的染色體，offset-adjusted（偏移調(diào)整的）染色體開始位置，offset-adjusted染色體末端的位置，和細(xì)胞barcode ID。這些被切分的片段存儲(chǔ)在臨時(shí)的HDF5-formatted文件里，保護(hù)內(nèi)存使用量，同時(shí)保持快速訪問每一個(gè)chunk。最后，所有與每個(gè)染色體相關(guān)聯(lián)的chunks都被讀取、組織和重寫到一個(gè)名為“fragments”的單個(gè)HDF5組中的Arrow file中。這種預(yù)分塊過程使ArchR能夠高效地處理超大的輸入文件，并降低內(nèi)存使用率，從而實(shí)現(xiàn)并行處理。
（5）Getting Set Up
我們要做的第一件事是更改到我們想要的工作目錄，設(shè)置我們想要使用的線程數(shù)量，并加載我們的基因和基因組注釋。根據(jù)本地環(huán)境的配置，您可能需要修改addArchRThreads()中使用的線程數(shù)。默認(rèn)情況下，ArchR使用了可用線程總數(shù)的一半，但你可以根據(jù)需要手動(dòng)調(diào)整。如果你使用的是Windows，可用的線程將自動(dòng)設(shè)置為1，因?yàn)锳rchR中的并行處理是基于unix的操作系統(tǒng)構(gòu)建的。
首先，調(diào)用ArchR包：
> library(ArchR)
接下來，我們?yōu)锳rchR設(shè)置默認(rèn)的線程數(shù)。在每次新的R程序中必須要做的事情。我們建議將線程設(shè)置為可用內(nèi)核總數(shù)的1/2到3/4。ArchR中的內(nèi)存使用通常會(huì)隨著線程數(shù)量的增加而增加，因此允許ArchR使用更多線程也就意味著要使用更多的內(nèi)存。
> addArchRThreads(threads = 16) #比如我現(xiàn)在用的是windows，我想設(shè)置16線程，但是ArchR自動(dòng)識(shí)別設(shè)備，并設(shè)置線程數(shù)為1
Detected windows OS, setting threads to 1.
Setting default number of Parallel threads to 1.
然后，我們?cè)O(shè)置基因組（用于注釋）。如上所述，這也是你在每次新的R會(huì)話中必須要做的。當(dāng)然，這個(gè)基因組版本必須與用于比對(duì)的基因組版本相匹配。對(duì)于本教程中使用的數(shù)據(jù)，我們將使用hg19基因組。
> addArchRGenome("hg19")
Setting default genome to Hg19.
ArchR需要基因和基因組注釋來做一些事情，比如計(jì)算TSS富集分?jǐn)?shù)、核苷酸含量和基因活性分?jǐn)?shù)。因?yàn)槲覀兊慕坛讨械膁ataset使用已經(jīng)與hg19參考基因組比對(duì)的scATAC-seq數(shù)據(jù)。但是，ArchR支持“hg19”、“hg38”、“mm9”和“mm10”，你可以使用createGeneAnnotation()和createGenomeAnnotation()函數(shù)創(chuàng)建自己的基因組和基因注釋。
通過addArchRGenome()函數(shù)向ArchR提供信息。這個(gè)函數(shù)告訴ArchR，對(duì)于當(dāng)前會(huì)話中的所有分析，它應(yīng)該使用與已定義的ArchRGenome相關(guān)的genomeAnnotation和geneAnnotation。每個(gè)支持的基因組由以下內(nèi)容組成：BSgenome對(duì)象（定義了基因組坐標(biāo)和每個(gè)染色體的序列），一個(gè)GRanges對(duì)象（包含一組“黑名單”區(qū)域），一個(gè)TxDb對(duì)象（定義了所有基因的位置和結(jié)構(gòu)），和一個(gè)OrgDb對(duì)象（提供了一個(gè)中央基因標(biāo)識(shí)符，包含這個(gè)標(biāo)識(shí)符之間和其他類型的標(biāo)識(shí)符的mapping）。
在ArchR中，hg19基因組的預(yù)編版本使用BSgenome.Hsapiens.UCSC，hg19 TxDb.Hsapiens.UCSC.hg19，knownGene org.Hs.eg，以及使用ArchR::mergeGR()合并的黑名單（來自hg19 v2 blacklist regions和來自mitochondrial regions that show high mappability to the hg19 nuclear genome）。
（6）創(chuàng)建ArchRGenome
如上所述，ArchRGenome由基因組注釋和基因注釋組成。
要?jiǎng)?chuàng)建自定義基因組注釋，我們使用createGenomeAnnotation()。為此，你將需要以下信息:
1.一個(gè)
BSgenome對(duì)象，它包含一個(gè)基因組的序列信息
(例如BSgenome.Hsapiens.UCSC.hg38)。
2.一個(gè)GRanges基因組范圍對(duì)象，包含一組黑名單區(qū)域，將用于過濾不需要的區(qū)域。這不是必需的，但建議這樣做。有關(guān)如何創(chuàng)建黑名單的信息，請(qǐng)參閱publication on the ENCODE blacklists。
這里，教程展示了如何從Drosophila melanogaster創(chuàng)建一個(gè)自定義基因組注釋，有需要的童鞋可以自行按照教程構(gòu)建，這里我就不翻譯了~因?yàn)槲抑皇褂眯∈蠛腿祟惢蚪M。官網(wǎng)代碼見：here
（7）創(chuàng)建Arrow Files
再次提醒~這個(gè)教程使用的數(shù)據(jù)來自文章Granja* et al. Nature Biotechnology 2019里的造血細(xì)胞數(shù)據(jù)集。包括骨髓單核細(xì)胞(BMMC)，外周血單核細(xì)胞(PBMC)和CD34+造血干細(xì)胞和祖細(xì)胞(CD34 BMMC)。
一旦我們有了fragment文件，給createArrowFiles()提供文件的路徑。在創(chuàng)建過程中，
創(chuàng)建期間,一些基本的元數(shù)據(jù)和矩陣被添加到每個(gè)Arrow file中，包括“TileMatrix”（包含在全基因組500bp bins的insertion counts），和“GeneScoreMatrix”（基于鄰近基因啟動(dòng)子的insetion counts推斷的基因表達(dá)量）。
示例數(shù)據(jù)使用getTutorialData()下載：
> 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"
現(xiàn)在我們開始創(chuàng)建Arrow Files，大概要花10-15分鐘。對(duì)于每一個(gè)樣品來說：
1.從input文件里讀取可接近性 fragments
2.計(jì)算每個(gè)細(xì)胞的質(zhì)量控制信息（如TSS富集分?jǐn)?shù)和核小體信息）
3.根據(jù)質(zhì)量控制參數(shù)過濾細(xì)胞
4.創(chuàng)建基因組范圍內(nèi)的TileMatrix（500bp的bin）
5.創(chuàng)建GeneScoreMatrix
> ArrowFiles <- createArrowFiles(
inputFiles = inputFiles,
sampleNames = names(inputFiles),
filterTSS = 4, #Don't set this too high because you can always increase later
filterFrags = 1000,
addTileMat = TRUE,
addGeneScoreMat = TRUE
)
#運(yùn)行時(shí)會(huì)彈出很多東西：
#從彈出信息的前兩行可以看出，上面的代碼里有兩個(gè)參數(shù)的名字已經(jīng)更改了，但是也可以運(yùn)行
filterFrags is no longer a valid input. Please use minFrags! Setting filterFrags value to minFrags!
filterTSS is no longer a valid input. Please use minTSS! Setting filterTSS value to minTSS!
Using GeneAnnotation set by addArchRGenome(Hg19)!
Using GeneAnnotation set by addArchRGenome(Hg19)!
ArchR logging to : ArchRLogs\ArchR-createArrows-3364659e77d1-Date-2020-11-17_Time-22-13-26.log
If there is an issue, please report to github with logFile!
2020-11-17 22:13:27 : Batch Execution w/ safelapply!, 0 mins elapsed.
Attempting to index HemeFragments/scATAC_BMMC_R1.fragments.tsv.gz as tabix..
2020-11-17 22:13:37 : (scATAC_BMMC_R1 : 1 of 3) Reading In Fragments from inputFiles (readMethod = tabix), 0.171 mins elapsed.
2020-11-17 22:13:37 : (scATAC_BMMC_R1 : 1 of 3) Tabix Bed To Temporary File, 0.174 mins elapsed.
......
2020-11-17 22:36:24 : (scATAC_PBMC_R1 : 3 of 3) Adding GeneScoreMatrix!, 22.962 mins elapsed.
2020-11-17 22:37:43 : (scATAC_PBMC_R1 : 3 of 3) Finished Creating Arrow File, 24.271 mins elapsed.
ArchR logging successful to : ArchRLogs\ArchR-createArrows-3364659e77d1-Date-2020-11-17_Time-22-13-26.log
現(xiàn)在我們看一下ArrowFiles對(duì)象，確認(rèn)一下它確實(shí)是一個(gè)Arrow file路徑的向量：
> ArrowFiles
[1] "scATAC_BMMC_R1.arrow" "scATAC_CD34_BMMC_R1.arrow" "scATAC_PBMC_R1.arrow"
（8）每個(gè)細(xì)胞的質(zhì)量控制
scATAC-seq數(shù)據(jù)的嚴(yán)格質(zhì)量控制是很重要的，首先就是要去掉低質(zhì)量細(xì)胞。在ArchR里，考慮下面三個(gè)因素：
1.唯一核fragments的數(shù)量（沒有比對(duì)到線粒體上的DNA）
2.信噪比。低信噪比通常是死細(xì)胞或者快死了的細(xì)胞，細(xì)胞會(huì)含有去染色質(zhì)化的DNA，使得Tn5在基因組上隨便切。
3.fragment大小的分布，由于核小體的周期性，我們希望fragment在147bp左右消失。
第一個(gè)指標(biāo)是唯一的核fragments，這很簡(jiǎn)單：只有很少的可用的fragments的細(xì)胞不能提供足夠的數(shù)據(jù)來進(jìn)行有用的分析，因此應(yīng)該被排除在外。
第二個(gè)度量標(biāo)準(zhǔn)，信噪比，被計(jì)算為TSS富集分?jǐn)?shù)。傳統(tǒng)的ATAC-seq分析使用這個(gè)TSS富集分?jǐn)?shù)作為確定signal-to-background的標(biāo)準(zhǔn)工作流的一部分。我們和其他人發(fā)現(xiàn)TSS富集在大部分ATAC-seq和scac -seq測(cè)試的細(xì)胞類型中具有代表性。TSS富集評(píng)分指標(biāo)的idea是，與其他基因組區(qū)域相比，由于與啟動(dòng)子結(jié)合的大型蛋白復(fù)合物，ATAC-seq數(shù)據(jù)在基因TSS區(qū)域普遍富集。通過觀察以這些TSS區(qū)域?yàn)橹行牡拿總€(gè)bp的可接近性，我們可以看到相對(duì)于側(cè)翼區(qū)域(在兩個(gè)方向上1900-2000 bp遠(yuǎn)端)的局部富集。峰的富集(以TSS為中心)與側(cè)翼區(qū)域的比值代表TSS富集分?jǐn)?shù)。
傳統(tǒng)上，對(duì)每個(gè)ATAC-seq樣本計(jì)算每個(gè)堿基對(duì)的可接近性，然后使用此profile確定TSS富集分?jǐn)?shù)。在scATAC-seq中，按細(xì)胞執(zhí)行此操作相對(duì)較慢，且計(jì)算成本較高。為了精確地評(píng)估每個(gè)細(xì)胞的TSS富集評(píng)分，我們計(jì)算以每個(gè)單堿基TSS位置為中心的50 bp區(qū)域內(nèi)的平均可接近性，并將其除以TSS側(cè)翼位置的平均可接近性(+/- 1900 - 2000 bp)。這個(gè)近似值(R > 0.99)與原始方法高度相關(guān)，并且數(shù)值在量級(jí)上非常接近。
第三個(gè)度量標(biāo)準(zhǔn)，fragment大小分布，通常不那么重要，但是手工檢查總是好的。由于DNA包裹核小體的方式，我們期望在我們的數(shù)據(jù)里看到fragment大小分布的周期性。出現(xiàn)這些“峰”和“谷”是因?yàn)槠伪仨毧缭?,1,2個(gè)核小體(Tn5不能切斷緊緊包裹在核小體周圍的DNA)。
默認(rèn)情況下，在ArchR中，通過過濾細(xì)胞被識(shí)別為TSS富集分?jǐn)?shù)大于4和超過1000個(gè)唯一核片段的細(xì)胞。需要注意的是，TSS富集分?jǐn)?shù)的實(shí)際數(shù)值取決于所使用的TSSs集。ArchR中的默認(rèn)值是為人類數(shù)據(jù)設(shè)計(jì)的，在運(yùn)行createArrowFiles()時(shí)更改默認(rèn)閾值是非常重要的。
在運(yùn)行完上面的代碼后，你的當(dāng)前工作文件夾里會(huì)多出一個(gè)文件夾，叫QualityControl。包含每個(gè)樣品的兩個(gè)相關(guān)的圖，第一個(gè)點(diǎn)圖顯示的是log10(unique nuclear fragments) vs TSS富集分?jǐn)?shù)，并顯示閾值；第二個(gè)圖是fragment大小分布圖：
比如我們打開第一個(gè)BMMC的結(jié)果來看：
現(xiàn)在我們已經(jīng)準(zhǔn)備整理這些Arrow files，后面將展示如何創(chuàng)建一個(gè)ArchRProject。
