本文首發(fā)于我的個(gè)人博客， http://xuzhougeng.top/
第1章: ArchR基礎(chǔ)入門
這一章將會(huì)介紹如何導(dǎo)入數(shù)據(jù)，如何構(gòu)建Arrow文件，這是后續(xù)ArchR分析的基礎(chǔ)。
1.1 ATAC-seq術(shù)語介紹
"fragment"是ATAC-seq實(shí)驗(yàn)中的一個(gè)重要概念，它指的是通過Tn5轉(zhuǎn)座酶對(duì)DNA分子進(jìn)行酶切，然后經(jīng)由雙端測(cè)序得到的序列。
我們會(huì)根據(jù)Tn5插入導(dǎo)致的偏移從read比對(duì)得到的位置推斷出fragment的起始和結(jié)束位置。根據(jù)之前的報(bào)道，Tn5轉(zhuǎn)座酶以同源二聚體的形式結(jié)合到DNA上，在兩個(gè)Tn5分子間隔著9-bp的DNA序列。根據(jù)這個(gè)情況，每個(gè)Tn5同源二聚體的結(jié)合事件會(huì)產(chǎn)生2個(gè)Insertions，中間隔著9bp。因此，真實(shí)的"開放"位置的中心在Tn5二聚體的正中間，而不是Tn5的插入位置。為了盡可能的還原真實(shí)情況，我們對(duì)Tn5的Insertions進(jìn)行了校正，即正鏈的插入結(jié)果往右移動(dòng)4bp(+4 bp), 負(fù)鏈的插入結(jié)果往左偏移5bp(-5 bp)。這個(gè)和最早提出的ATAC-seq里的描述是一致的。因此，在ArchR中，"fragment"指的是一個(gè)table或genomic ranges對(duì)象, 記錄在染色體上，經(jīng)過偏移校正后的單堿基起始位置，以及經(jīng)過偏移校正后單堿基結(jié)束位置，每個(gè)fragment都會(huì)對(duì)應(yīng)唯一的細(xì)胞條形碼。類似的，"Insertions"這得是偏移校正后的單堿基位置，它位于開放位置的正中心。
"fragment"和"insertions"這兩個(gè)詞我并沒有將其翻譯成中文，我覺得這兩個(gè)單詞可能就和PCR一樣，當(dāng)說到它們的時(shí)候，腦中會(huì)有一個(gè)畫面描述它們，而不是局限一個(gè)詞。
1.2 為什么是用ArchR
現(xiàn)在其實(shí)已經(jīng)有一些工具能夠處理單細(xì)胞ATAC-seq數(shù)據(jù)，我們?yōu)槭裁匆~外造一個(gè)輪子，開發(fā)一個(gè)新的項(xiàng)目呢？主要是ArchR提供了其他工具目前尚不能實(shí)現(xiàn)的功能
不僅如此，ArchR通過優(yōu)化數(shù)據(jù)結(jié)構(gòu)降低了內(nèi)存消耗，使用并行提高了運(yùn)行速度，因此保證其性能優(yōu)于其他同類型工具。在超過70,000個(gè)細(xì)胞的分析項(xiàng)目中，一些軟件需要超過128G的可用內(nèi)存。
ArchR最初就是根據(jù)Unix的筆記本進(jìn)行設(shè)計(jì)（我覺得他說的是MacBook Pro），因此對(duì)于中等大小的實(shí)驗(yàn)（小于100,000個(gè)細(xì)胞），ArchR能夠保證一些特殊分析的運(yùn)行速度，并能實(shí)時(shí)的展示結(jié)果，讓我們能夠更深入的和數(shù)據(jù)進(jìn)行互動(dòng)，給出有意義的生物學(xué)解釋。當(dāng)然，如果細(xì)胞數(shù)更多，你最好使用服務(wù)器進(jìn)行分析。ArchR提供了方便的圖形導(dǎo)出功能，在服務(wù)器處理完項(xiàng)目之后，可以直接下載到本地進(jìn)行使用。
目前，ArchR并沒有針對(duì)Windows進(jìn)行優(yōu)化。這句話的意思是指，ArchR的并行策略是基于Unix系統(tǒng)而非Windows系統(tǒng)，因此上述說的性能提升不包括Windows。
1.3 什么是Arrow文件/ArchRProject
正如開頭所說，ArchR分析項(xiàng)目的基礎(chǔ)是Arrow文件。每個(gè)Arrow文件記錄著每個(gè)獨(dú)立樣本的所有相關(guān)信息（例如元信息、開放的fragment和數(shù)據(jù)矩陣）。這里說的獨(dú)立樣本最好是最詳盡的分析單元（例如，一種特定條件下的單個(gè)重復(fù)）。在創(chuàng)建Arrow文件以及一些附加分析中，ArchR會(huì)編輯和更新相應(yīng)的Arrow文件，在其中添加額外的信息層。值得注意的是，Arrow文件實(shí)際指的是磁盤上的文件路徑。更確切的說，Arrow文件并不是一個(gè)存放在內(nèi)存中的R語言對(duì)象，而是存放在磁盤上HDF5文件。正因如此，我們使用ArchRProject對(duì)象用來關(guān)聯(lián)這些Arrow文件，將其關(guān)聯(lián)到單個(gè)分析框架中，從而保證在R中能高效訪問它們。而這個(gè)ArchRProject對(duì)象占用內(nèi)存不多，因此才是存放在內(nèi)存中的R語言對(duì)象。
有一些操作會(huì)直接修改Arrow文件，而一些操作會(huì)先作用于ArchRproject，接著反過來更新每個(gè)相關(guān)Arrow文件。因?yàn)锳rrow文件是以非常大的HDF5格式存放，所以ArchR的get-er函數(shù)通過和ArchRProject進(jìn)行交互獲取數(shù)據(jù)，而add-er函數(shù)既能直接在Arrow文件中添加數(shù)據(jù)，也能直接在ArchRpoject里添加數(shù)據(jù)，或者通過ArchRpoject向Arrow文件添加數(shù)據(jù)。
1.4 ArchR的輸入文件類型
ArchR主要以scATAC-seq原始數(shù)據(jù)經(jīng)上游處理后的兩種常見輸出文件(BAM, fragment)作為輸入。Fragment文件記錄著scATAC-seq的fragment以及對(duì)應(yīng)的細(xì)胞ID，每一行都是一條記錄，該文件需要是tabix(見注1)排序并建立索引保證能被高效讀取。BAM文件則是二進(jìn)制格式下的tabix排序文件，記錄著scATAC-seq的fragment、原始數(shù)據(jù)、細(xì)胞條形碼和其他信息。具體使用何種文件作為輸入取決于你的上游處理流程。以10XGenomics的CellRanger軟件為例，它的輸出文件是fragments，而sci-ATAC-seq流程則輸出BAM文件。
ArchR使用scanTabix函數(shù)讀取fragment文件，使用scanBAM讀取BAM文件。輸入過程并不會(huì)直接讀取所有數(shù)據(jù)，而是每次讀取一大塊(chunks)，然后將這一塊數(shù)據(jù)轉(zhuǎn)換成fragment格式(見注2)，經(jīng)過壓縮先暫時(shí)以HDF5格式保存到本地磁盤上，避免消耗過多的內(nèi)存。最后，等所有數(shù)據(jù)都加載完成之后，該數(shù)據(jù)相關(guān)的所有臨時(shí)HDF5文件會(huì)被重新讀取，經(jīng)過組織之后以單個(gè)HDF5形式保存為Arrow文件。ArchR之所以能以較小的內(nèi)存量高效地讀取大文件，就是因?yàn)樗捎玫氖沁@種分塊處理的方法，由于每一塊數(shù)據(jù)的處理都互不干擾，因此也就能夠并行計(jì)算。
- 注1: 當(dāng)以制表符記錄基因組位置信息時(shí)，我們可以通過壓縮讓其體積變小(bgzip/gzip)，通過建立tabix索引高效訪問給定位置的信息。
- 注2: fragments一共有5列，分別是chromosome, offset-adjusted chromosome start position, offset-adjusted chromosome end position, and cellular barcode ID, duplicate count
1.5 開始設(shè)置
在后續(xù)分析之前，請(qǐng)先設(shè)置好我們的工作目錄，設(shè)置將要使用的線程數(shù)，加載我們的基因和基因組注釋。由于每個(gè)人的環(huán)境都不太一樣，所以你后續(xù)可能需要用addArchRThreads()修改線程數(shù)。默認(rèn)情況下，ArchR使用系統(tǒng)一半的可用線程，你可以手動(dòng)進(jìn)行調(diào)整。如果你用的是Windows，那么默認(rèn)都是1，這是因?yàn)锳rchR的多線程是基于Unix操作系統(tǒng)。
第零步，安裝ArchR。目前ArchR托管在GitHub上，因此無法直接用CRAN或者Bioconductor中直接下載安裝。
方法1：適用于網(wǎng)絡(luò)狀態(tài)好的情況
if (!requireNamespace("devtools", quietly = TRUE)) install.packages("devtools")
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
devtools::install_github("GreenleafLab/ArchR", ref="master", repos = BiocManager::repositories())
方法2: 適用于上述方法多次失敗的情況
先從Github上下載項(xiàng)目到本地, 大約有262Mb
git clone https://github.com/GreenleafLab/ArchR.git
然后在 R里面進(jìn)行安裝
BiocManager::install(c("nabor","motifmatchr","chromVAR","ComplexHeatmap"))
install.packages("./ArchR", repo=NULL)
接著安裝一些額外包
ArchR::installExtraPackages()
如果安裝失敗，就手動(dòng)安裝 Seurat, immunogenomics/harmony, immunogenomics/presto, Cairo，shiny, shinythemes, rhandsontable
第一步，我們加載ArchR包。
library(ArchR)
第二步, 我們需要設(shè)置ArchR函數(shù)的默認(rèn)線程數(shù)。你需要在每個(gè)新的R session中都設(shè)置該參數(shù)，線程多多益善，但是不要超過總線程的3/4。因?yàn)榫€程會(huì)和內(nèi)存掛鉤，所以線程數(shù)越多，內(nèi)存相應(yīng)使用的也越多。
addArchRThreads(threads = 16)
第三步: 我們?cè)O(shè)置需要使用基因和基因組注釋。同樣，這也是每個(gè)新的R session需要設(shè)置的參數(shù)。當(dāng)然，我們需要使用和序列比對(duì)階段相同的基因組版本。對(duì)于本教程使用的數(shù)據(jù)，我們會(huì)使用hg19參考基因組。ArchR支持多種基因組注釋，并且允許自定義基因組注釋。
addArchRGenome("hg19")
基因和基因組注釋信息是后續(xù)計(jì)算TSS富集得分，核小體含量和基因活躍度得分所必需的。同樣，你得保證這里選擇的基因組版本得和你上游數(shù)據(jù)處理時(shí)用到的基因組版本一致。我們這里設(shè)置"hg19"就是因?yàn)樯蠑?shù)據(jù)處理用的是hg19作為參考基因組。除了hg19外，ArchR還提供了"hg38", "mm9", "mm10", 并且允許用戶使用createGeneAnnotation和createGenomeAnnotation函數(shù)創(chuàng)建其他物種注釋。
我們使用addArchRGenome函數(shù)無縫地為ArchR提供這些信息。該函數(shù)告訴ArchR，在之后的所有分析中，它應(yīng)該使用定義在ArchRgenome的相關(guān)genomeAnnotation和geneAnnotation。每一個(gè)原生支持的基因組都包括四個(gè)對(duì)象(object)
- BSgenome: 記錄 染色體坐標(biāo)信息和染色體序列信息
- GRanges: 記錄blacklist, 即對(duì)分析沒有用設(shè)置可能產(chǎn)生干擾的區(qū)域
- TxDb: 定義所有基因的位置信息
- OrgDb: 提供基因編號(hào)，以及不同基因編號(hào)之間的轉(zhuǎn)換
如下是ArchR自帶的物種注釋的數(shù)據(jù)來源
ArchR的預(yù)編譯的hg19基因組用的是BSgenome.Hsapiens.UCSC.hg19, TxDb.Hsapiens.UCSC.hg19.knownGene, org.Hs.eg.db。而黑名單(記錄著一直打開的區(qū)域，對(duì)后續(xù)分析沒有幫助的位置）則來源于 hg19 v2 blacklist regions 和mitochondrial regions that show high mappability to the hg19 nuclear genome from Caleb Lareau and Jason Buenrostro，用ArchR::mergeGR()函數(shù)進(jìn)行合并
ArchR的預(yù)編譯的hg38基因組用的是BSgenome.Hsapiens.UCSC.hg38, TxDb.Hsapiens.UCSC.hg38.knownGene, org.Hs.eg.db。而黑名單(記錄著一直打開的區(qū)域，對(duì)后續(xù)分析沒有幫助的位置）則來源于 hg19 v2 blacklist regions 和mitochondrial regions that show high mappability to the hg38 nuclear genome from Caleb Lareau and Jason Buenrostro，用ArchR::mergeGR()函數(shù)進(jìn)行合并
ArchR的預(yù)編譯的mm9基因組用的是BSgenome.Mmusculus.UCSC.mm9, TxDb.Mmusculus.UCSC.mm9.knownGene, org.Mm.eg.db。而黑名單(記錄著一直打開的區(qū)域，對(duì)后續(xù)分析沒有幫助的位置）則來源于 mm9 v1 blacklist regions 和mitochondrial regions that show high mappability to the mm9 nuclear genome from Caleb Lareau and Jason Buenrostro，用ArchR::mergeGR()函數(shù)進(jìn)行合并
ArchR的預(yù)編譯的mm10基因組用的是BSgenome.Mmusculus.UCSC.mm10, TxDb.Mmusculus.UCSC.mm10.knownGene, org.Mm.eg.db。而黑名單(記錄著一直打開的區(qū)域，對(duì)后續(xù)分析沒有幫助的位置）則來源于 mm10 v2 blacklist regions 和mitochondrial regions that show high mappability to the mm10 nuclear genome from Caleb Lareau and Jason Buenrostro，用ArchR::mergeGR()函數(shù)進(jìn)行合并
1.5.1: 自定義ArchRGenome
如上所述，一個(gè)ArchRGenome需要包括基因組注釋和基因注釋
為了構(gòu)建自定義的基因組注釋，我們會(huì)使用createGenomeAnnotation. 他需要如下2個(gè)信息
- BSgenome: 包括基因組的序列信息，你可以通過谷歌查找指定物種的Bioconductor包
- GRanges: 記錄基因組中的黑名單區(qū)域，用來過濾下游分析中不需要的區(qū)間。這不是必須的，畢竟不是所有物種都能有這個(gè)文件。 publication on the ENCODE blacklists提供黑名單列表的制作方法。
例如，我們?nèi)绻枰獮镈rosophila melanogaster創(chuàng)建一個(gè)基因組注釋，那么我們需要先下載對(duì)應(yīng)的BSgenome
if (!requireNamespace("BSgenome.Dmelanogaster.UCSC.dm6", quietly = TRUE)){
BiocManager::install("BSgenome.Dmelanogaster.UCSC.dm6")
}
library(BSgenome.Dmelanogaster.UCSC.dm6)
接著，我們從BSgenome對(duì)象中創(chuàng)建基因組注釋
genomeAnnotation <- createGenomeAnnotation(genome = BSgenome.Dmelanogaster.UCSC.dm6)
查看這個(gè)對(duì)象，你可以看到一個(gè)ArchR基因組注釋的組成內(nèi)容
genomeAnnotation
## List of length 3
## names(3): genome chromSizes blacklist
為了構(gòu)造一個(gè)自定義基因注釋，我們需要用到createGeneAnnotation()，它需要你提供如下兩個(gè)信息
- TxDb: 記錄gene/transcript的坐標(biāo)信息
- OrgDb: 用于基因名和其他基因編號(hào)的轉(zhuǎn)換
繼續(xù)之前的Drosophila melanogaster案例，我們需要先安裝并加載相關(guān)的TxDb和OrgDb對(duì)象
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)
library(org.Dm.eg.db)
記著，我們構(gòu)建基因注釋對(duì)象
geneAnnotation <- createGeneAnnotation(TxDb = TxDb.Dmelanogaster.UCSC.dm6.ensGene, OrgDb = org.Dm.eg.db)
檢查該對(duì)象，可以了解ArchR關(guān)于基因注釋的存放形式
geneAnnotation
## List of length 3
## names(3): genes exons TSS
如果你沒有TxDb和OrgDb對(duì)象，你也可以直接根據(jù)如下信息創(chuàng)建geneAnnotation對(duì)象
- "genes": GRanges對(duì)象，記錄基因的坐標(biāo)信息，起始位置和結(jié)束。必須要有一列是和"exons"對(duì)象中其中一列匹配
- "exons": 記錄每個(gè)基因外顯子的坐標(biāo)。必須要有一列和"genes"對(duì)象中的一列匹配
- "GRanges": 記錄TSS(轉(zhuǎn)錄起始位點(diǎn))的坐標(biāo)
geneAnnotation <- createGeneAnnotation(
TSS = geneAnnotation$TSS,
exons = geneAnnotation$exons,
genes = geneAnnotation$genes
)
1.5.2 使用非標(biāo)準(zhǔn)基因組
ArchR會(huì)做一些必要的檢查，來避免你做一些ArchR覺得"不符合常理"的操作。其中就是檢查你的基因組注釋的seqnames都需要以"chr"開頭。大部分情況下這都沒有問題，除了一些基因組的染色體并不都是以"chr"作為前綴（例如玉米）。如果用ArchR分析這些基因組，你需要告知ArchR忽略染色體名前綴，否則會(huì)報(bào)錯(cuò)停止。你需要在創(chuàng)建Arrow文件前，先運(yùn)行addArchRChrPrefix(chrPrefix = FALSE). 它會(huì)當(dāng)前的R session里停止所有對(duì)染色體名前綴的檢查操作。
此外，ArchR默認(rèn)會(huì)將染色體名/seqnames轉(zhuǎn)成字符串，因此如果你的seqnames都是數(shù)字，你需要以字符串形式提供這seqnames, 例如，你需要執(zhí)行useSeqnames = c("1", "2", "3")，而不是useSeqnames = c(1, 2, 3)
你可以用getArchRChrPrefix隨時(shí)檢查當(dāng)前R session是否需要染色體前綴。
1.6 創(chuàng)建Arrow文件
在本教程的后續(xù)部分，我們會(huì)使用Granja* et al. Nature Biotechnology 2019文章的數(shù)據(jù)進(jìn)行展示，當(dāng)然不是所有的數(shù)據(jù)集，我們會(huì)使用降抽樣的造血細(xì)胞數(shù)據(jù)集，保證大部分人的電腦都能帶動(dòng)。該數(shù)據(jù)集包括了骨髓單核細(xì)胞(BMMC)和外周血單核細(xì)胞(PBMC)，以及CD34+造血干細(xì)胞和骨髓前體細(xì)胞(CD34 BMMC)。
我們下載的數(shù)據(jù)以fragment格式存放，記錄每個(gè)比對(duì)序列在基因組上的位置。Fragments文件是10X Genomics分析平臺(tái)的其中一個(gè)輸出文件，或者你也可以從BAM文件進(jìn)
一旦我們得到了fragment文件，我們將它們的路徑記錄在一個(gè)向量中，然后作為參數(shù)傳給createArrowFiles(). 在構(gòu)建過程中，一些基本的元信息和矩陣會(huì)增加到各個(gè)Arrow文件中，包括TileMatrix和GeneScoreMatrix. 其中TileMarix記錄的以500-bp作為分窗統(tǒng)計(jì)染色體上各個(gè)位置上是否有insertion（具體見addTileMatrix）, GeneScoreMatrix則是基于鄰近基因啟動(dòng)子的insetion數(shù)推斷的基因表達(dá)量（見addGeneScoreMatrix()）。
教程用到的數(shù)據(jù)，可以用getTutorialData進(jìn)行下載，大約為0.5G。
inputFiles <- getTutorialData("Hematopoiesis")
當(dāng)然這一步實(shí)際是在本地建立一個(gè)HemeFragments目錄，并下載指定數(shù)據(jù)，因此如果網(wǎng)絡(luò)太差，可以自己嘗試下載。
mkdir HemeFragments && cd HemeFragments
wget 'https://jeffgranja.s3.amazonaws.com/ArchR/TestData/HemeFragments/scATAC_BMMC_R1.fragments.tsv.gz'
wget 'https://jeffgranja.s3.amazonaws.com/ArchR/TestData/HemeFragments/scATAC_CD34_BMMC_R1.fragments.tsv.gz'
wget 'https://jeffgranja.s3.amazonaws.com/ArchR/TestData/HemeFragments/scATAC_PBMC_R1.fragments.tsv.gz'
在啟動(dòng)項(xiàng)目前，我們必須得設(shè)置ArchRGenome并設(shè)置線程數(shù)threads
addArchRGenome("hg19") # hg38, mm9, mm10
addArchRThreads(threads = 16)
從現(xiàn)在開始，我們將會(huì)用10-15分鐘的時(shí)間創(chuàng)建Arrow文件。對(duì)每一個(gè)樣本，都會(huì)有如下操作
- 從給定文件中讀取開放信息(fragments)
- 為每個(gè)細(xì)胞計(jì)算質(zhì)控信息(TSS富集得分和核小體信息)
- 根據(jù)質(zhì)控參數(shù)過濾細(xì)胞
- 以500bp為窗口構(gòu)建全基因組范圍的TileMatrix
- 使用自定義的
geneAnnotaiton創(chuàng)建GeneScoreMatrix, 其中geneAnnotation在我們調(diào)用addArchRGenome時(shí)定義
ArrowFiles <- createArrowFiles(
inputFiles = inputFiles,
sampleNames = names(inputFiles),
filterTSS = 4, # 這個(gè)參數(shù)不需要過高，后續(xù)可以調(diào)整
filterFrags = 1000,
addTileMat = TRUE,
addGeneScoreMat = TRUE
)
我們可以看下ArrowFiles對(duì)象，檢查它是否真的是一個(gè)存放Arrow文件路徑的向量
ArrowFiles
# [1] "scATAC_BMMC_R1.arrow" "scATAC_CD34_BMMC_R1.arrow"
# [3] "scATAC_PBMC_R1.arrow"
1.7 細(xì)胞質(zhì)控
scATAC-seq的嚴(yán)格質(zhì)控對(duì)于剔除低質(zhì)量細(xì)胞至關(guān)重要。在ArchR中，我們考慮數(shù)據(jù)的三種信息
- 每個(gè)細(xì)胞中唯一比對(duì)數(shù)(number of unique nuclear fragments)，不包括比對(duì)到線粒體的部分(unique nuclear fragments, 類似于單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)的表達(dá)的基因量)
- 信噪比(signal-to-background ratio)。如果是死細(xì)胞或者快死的細(xì)胞，那么DNA傾向于去染色質(zhì)化，就會(huì)導(dǎo)致全局轉(zhuǎn)座酶隨機(jī)切割，體現(xiàn)出來就是信噪比低。
- fragment大小分布( fragment size distribution)。 由于核小體周期性，我們期望看到在147-bp附近出現(xiàn)一個(gè)低谷，因?yàn)槔p繞一個(gè)核小體的DNA序列大約為147bp。
第一個(gè)參數(shù), 唯一比對(duì)數(shù)非常的直觀，如果一個(gè)細(xì)胞中的唯一比對(duì)太少，顯然在后續(xù)分析中也沒有太多可用價(jià)值，可以直接剔除掉。
第二個(gè)參數(shù)，信噪比是根據(jù)TSS富集分?jǐn)?shù)進(jìn)行計(jì)算。傳統(tǒng)的混池ATAC-seq分析中，會(huì)使用TSS富集得分作為標(biāo)準(zhǔn)流程的一部分，用于確定信噪比（如 ENCODE計(jì)劃）。我們和其他人通過混池ATAC-seq和scATAC-seq分析發(fā)現(xiàn)，TSS富集得分在大部分的細(xì)胞類型中都具有代表性。TSS富集得分的背后思想是，由于大蛋白復(fù)合體會(huì)結(jié)合在啟動(dòng)子區(qū)域，ATAC-seq數(shù)據(jù)更多的富集在基因的TSS區(qū)域而不基因組其他區(qū)域。通過檢查這些TSS區(qū)域中心的開放水平，我們發(fā)現(xiàn)中心相對(duì)于兩側(cè)(兩邊的1900-2000 bp處)存在富集現(xiàn)象。因此，我們定義富集中心(TSS中心)相對(duì)于兩側(cè)區(qū)域的比值為TSS富集得分(TSS enrichment score)
傳統(tǒng)上，我們會(huì)計(jì)算每個(gè)混池ATAC-seq樣本中每個(gè)堿基的開放性，之后這個(gè)譜會(huì)被用于確定TSS富集得分。在scATAC-seq中通過該方法計(jì)算每個(gè)細(xì)胞的TSS富集得分速度較慢，并且需要很強(qiáng)的算力。為了提高計(jì)算效率，同時(shí)還能得到和傳統(tǒng)計(jì)算接近的結(jié)果，我們以TSS位置為中心，以50-bp作為分窗計(jì)算兩邊的平均開放程度，之后該值除以TSS兩側(cè)位置(+-1900-200bp)的平均開放程度。 這種計(jì)算方法和原來相對(duì)比，兩者的相關(guān)性大于0.99，并且結(jié)果幾乎一樣。
第三個(gè)參數(shù)fragment大小分布并不是特別重要，最好是人工檢查下。由于DNA纏繞核小體的模式，我們預(yù)期在fragment大小分布中看到核小體的周期性分布。這些山峰和低谷的出現(xiàn)正是由于DNA纏繞0，1，2個(gè)核小體的結(jié)果(Tn5無法切割纏繞在核小體的DNA序列，只能切割兩邊）
ArchR默認(rèn)會(huì)過濾TSS富集得分低于4或唯一比對(duì)數(shù)小于1000（也就是保留TSS富集得分大于4且唯一比對(duì)數(shù)大于1000的細(xì)胞）。切記，ArchR默認(rèn)的參數(shù)最初是用于人類數(shù)據(jù)，不能直接用于其他物種。每個(gè)數(shù)據(jù)都有他的獨(dú)特性，切勿生搬硬套，我們需要按照實(shí)際情況設(shè)置參數(shù)。一定要在createArrowFiles()運(yùn)行前設(shè)置好該參數(shù)。
創(chuàng)建Arrow文件會(huì)在當(dāng)前目錄下生成一個(gè)"QualityControl"目錄，這里面包括2個(gè)和樣本相關(guān)的圖。第一個(gè)圖展示log10(unique nuclear fragments) vs TSS enrichment score, 虛線表示過濾閾值。第二圖注釋fragment大小分布圖。
對(duì)我們的教程數(shù)據(jù)，我們的三個(gè)樣本表現(xiàn)如下
BMMC:
CD34+ BMMC:
PBMC:
我們現(xiàn)在整理完畢我們的Arrow文件，接下來就是創(chuàng)建ArchRProject了。
