第1章: ArchR基礎(chǔ)入門
這一章將會(huì)介紹如何導(dǎo)入數(shù)據(jù)，如何構(gòu)建Arrow文件算谈，這是后續(xù)ArchR分析的基礎(chǔ)。
1.1 ATAC-seq術(shù)語(yǔ)介紹
"fragment"是ATAC-seq實(shí)驗(yàn)中的一個(gè)重要概念瞳秽，它指的是通過(guò)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í)的"開(kāi)放"位置的中心在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)過(guò)偏移校正后的單堿基起始位置，以及經(jīng)過(guò)偏移校正后單堿基結(jié)束位置千元，每個(gè)fragment都會(huì)對(duì)應(yīng)唯一的細(xì)胞條形碼苫昌。類似的，"Insertions"這得是偏移校正后的單堿基位置幸海，它位于開(kāi)放位置的正中心祟身。
"fragment"和"insertions"這兩個(gè)詞我并沒(méi)有將其翻譯成中文奥务，我覺(jué)得這兩個(gè)單詞可能就和PCR一樣，當(dāng)說(shuō)到它們的時(shí)候袜硫，腦中會(huì)有一個(gè)畫(huà)面描述它們氯葬，而不是局限一個(gè)詞。
1.2 為什么是用ArchR
現(xiàn)在其實(shí)已經(jīng)有一些工具能夠處理單細(xì)胞ATAC-seq數(shù)據(jù)父款，我們?yōu)槭裁匆~外造一個(gè)輪子溢谤，開(kāi)發(fā)一個(gè)新的項(xiàng)目呢放案？主要是ArchR提供了其他工具目前尚不能實(shí)現(xiàn)的功能
不僅如此禀苦，ArchR通過(guò)優(yōu)化數(shù)據(jù)結(jié)構(gòu)降低了內(nèi)存消耗，使用并行提高了運(yùn)行速度乳幸，因此保證其性能優(yōu)于其他同類型工具肝集。在超過(guò)70,000個(gè)細(xì)胞的分析項(xiàng)目中瞻坝，一些軟件需要超過(guò)128G的可用內(nèi)存。
ArchR最初就是根據(jù)Unix的筆記本進(jìn)行設(shè)計(jì)（我覺(jué)得他說(shuō)的是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并沒(méi)有針對(duì)Windows進(jìn)行優(yōu)化厕倍。這句話的意思是指，ArchR的并行策略是基于Unix系統(tǒng)而非Windows系統(tǒng)贩疙，因此上述說(shuō)的性能提升不包括Windows讹弯。
1.3 什么是Arrow文件/ArchRProject
正如開(kāi)頭所說(shuō)，ArchR分析項(xiàng)目的基礎(chǔ)是Arrow文件这溅。每個(gè)Arrow文件記錄著每個(gè)獨(dú)立樣本的所有相關(guān)信息（例如元信息闸婴、開(kāi)放的fragment和數(shù)據(jù)矩陣）。這里說(shuō)的獨(dú)立樣本最好是最詳盡的分析單元（例如芍躏，一種特定條件下的單個(gè)重復(fù)）。在創(chuàng)建Arrow文件以及一些附加分析中降狠，ArchR會(huì)編輯和更新相應(yīng)的Arrow文件对竣，在其中添加額外的信息層庇楞。值得注意的是，Arrow文件實(shí)際指的是磁盤上的文件路徑否纬。更確切的說(shuō)吕晌，Arrow文件并不是一個(gè)存放在內(nèi)存中的R語(yǔ)言對(duì)象，而是存放在磁盤上HDF5文件临燃。正因如此睛驳，我們使用ArchRProject對(duì)象用來(lái)關(guān)聯(lián)這些Arrow文件，將其關(guān)聯(lián)到單個(gè)分析框架中膜廊，從而保證在R中能高效訪問(wèn)它們乏沸。而這個(gè)ArchRProject對(duì)象占用內(nèi)存不多，因此才是存放在內(nèi)存中的R語(yǔ)言對(duì)象爪瓜。
有一些操作會(huì)直接修改Arrow文件蹬跃，而一些操作會(huì)先作用于ArchRproject，接著反過(guò)來(lái)更新每個(gè)相關(guān)Arrow文件铆铆。因?yàn)锳rrow文件是以非常大的HDF5格式存放蝶缀，所以ArchR的get-er函數(shù)通過(guò)和ArchRProject進(jìn)行交互獲取數(shù)據(jù)，而add-er函數(shù)既能直接在Arrow文件中添加數(shù)據(jù)薄货，也能直接在ArchRpoject里添加數(shù)據(jù)翁都，或者通過(guò)ArchRpoject向Arrow文件添加數(shù)據(jù)。
1.4 ArchR的輸入文件類型
ArchR主要以scATAC-seq原始數(shù)據(jù)經(jīng)上游處理后的兩種常見(jiàn)輸出文件(BAM, fragment)作為輸入谅猾。Fragment文件記錄著scATAC-seq的fragment以及對(duì)應(yīng)的細(xì)胞ID柄慰，每一行都是一條記錄，該文件需要是tabix(見(jiàn)注1)排序并建立索引保證能被高效讀取赊瞬。BAM文件則是二進(jìn)制格式下的tabix排序文件先煎，記錄著scATAC-seq的fragment、原始數(shù)據(jù)巧涧、細(xì)胞條形碼和其他信息薯蝎。具體使用何種文件作為輸入取決于你的上游處理流程。以10XGenomics的CellRanger軟件為例谤绳，它的輸出文件是fragments占锯，而sci-ATAC-seq流程則輸出BAM文件。
ArchR使用scanTabix函數(shù)讀取fragment文件缩筛，使用scanBAM讀取BAM文件消略。輸入過(guò)程并不會(huì)直接讀取所有數(shù)據(jù)，而是每次讀取一大塊(chunks)瞎抛，然后將這一塊數(shù)據(jù)轉(zhuǎn)換成fragment格式(見(jiàn)注2)艺演，經(jīng)過(guò)壓縮先暫時(shí)以HDF5格式保存到本地磁盤上，避免消耗過(guò)多的內(nèi)存。最后胎撤，等所有數(shù)據(jù)都加載完成之后晓殊，該數(shù)據(jù)相關(guān)的所有臨時(shí)HDF5文件會(huì)被重新讀取，經(jīng)過(guò)組織之后以單個(gè)HDF5形式保存為Arrow文件伤提。ArchR之所以能以較小的內(nèi)存量高效地讀取大文件巫俺，就是因?yàn)樗捎玫氖沁@種分塊處理的方法，由于每一塊數(shù)據(jù)的處理都互不干擾肿男，因此也就能夠并行計(jì)算介汹。
- 注1: 當(dāng)以制表符記錄基因組位置信息時(shí)，我們可以通過(guò)壓縮讓其體積變小(bgzip/gzip)舶沛，通過(guò)建立tabix索引高效訪問(wèn)給定位置的信息嘹承。
- 注2: fragments一共有5列，分別是chromosome, offset-adjusted chromosome start position, offset-adjusted chromosome end position, and cellular barcode ID, duplicate count
1.5 開(kāi)始設(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上，因此無(wú)法直接用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ù)舞痰，線程多多益善土榴，但是不要超過(guò)總線程的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ù)無(wú)縫地為ArchR提供這些信息准夷。該函數(shù)告訴ArchR，在之后的所有分析中莺掠，它應(yīng)該使用定義在ArchRgenome的相關(guān)genomeAnnotation和geneAnnotation衫嵌。每一個(gè)原生支持的基因組都包括四個(gè)對(duì)象(object)
- BSgenome: 記錄 染色體坐標(biāo)信息和染色體序列信息
- GRanges: 記錄blacklist, 即對(duì)分析沒(méi)有用設(shè)置可能產(chǎn)生干擾的區(qū)域
- TxDb: 定義所有基因的位置信息
- OrgDb: 提供基因編號(hào)，以及不同基因編號(hào)之間的轉(zhuǎn)換
如下是ArchR自帶的物種注釋的數(shù)據(jù)來(lái)源
ArchR的預(yù)編譯的hg19基因組用的是BSgenome.Hsapiens.UCSC.hg19, TxDb.Hsapiens.UCSC.hg19.knownGene, org.Hs.eg.db彻秆。而黑名單(記錄著一直打開(kāi)的區(qū)域楔绞，對(duì)后續(xù)分析沒(méi)有幫助的位置）則來(lái)源于 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唇兑。而黑名單(記錄著一直打開(kāi)的區(qū)域酒朵，對(duì)后續(xù)分析沒(méi)有幫助的位置）則來(lái)源于 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扎附。而黑名單(記錄著一直打開(kāi)的區(qū)域蔫耽，對(duì)后續(xù)分析沒(méi)有幫助的位置）則來(lái)源于 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帕棉。而黑名單(記錄著一直打開(kāi)的區(qū)域针肥，對(duì)后續(xù)分析沒(méi)有幫助的位置）則來(lái)源于 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: 包括基因組的序列信息，你可以通過(guò)谷歌查找指定物種的Bioconductor包
- GRanges: 記錄基因組中的黑名單區(qū)域即纲，用來(lái)過(guò)濾下游分析中不需要的區(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
如果你沒(méi)有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ì)做一些必要的檢查碌识，來(lái)避免你做一些ArchR覺(jué)得"不符合常理"的操作碾篡。其中就是檢查你的基因組注釋的seqnames都需要以"chr"開(kāi)頭。大部分情況下這都沒(méi)有問(wèn)題丸冕，除了一些基因組的染色體并不都是以"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)建過(guò)程中瞻佛，一些基本的元信息和矩陣會(huì)增加到各個(gè)Arrow文件中脱篙，包括TileMatrix和GeneScoreMatrix. 其中TileMarix記錄的以500-bp作為分窗統(tǒng)計(jì)染色體上各個(gè)位置上是否有insertion（具體見(jiàn)addTileMatrix）, GeneScoreMatrix則是基于鄰近基因啟動(dòng)子的insetion數(shù)推斷的基因表達(dá)量（見(jiàn)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)在開(kāi)始疯兼，我們將會(huì)用10-15分鐘的時(shí)間創(chuàng)建Arrow文件然遏。對(duì)每一個(gè)樣本，都會(huì)有如下操作
- 從給定文件中讀取開(kāi)放信息(fragments)
- 為每個(gè)細(xì)胞計(jì)算質(zhì)控信息(TSS富集得分和核小體信息)
- 根據(jù)質(zhì)控參數(shù)過(guò)濾細(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ù)不需要過(guò)高吧彪，后續(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)出來(lái)就是信噪比低赡艰。
- fragment大小分布( fragment size distribution)售淡。 由于核小體周期性，我們期望看到在147-bp附近出現(xiàn)一個(gè)低谷慷垮，因?yàn)槔p繞一個(gè)核小體的DNA序列大約為147bp揖闸。
第一個(gè)參數(shù), 唯一比對(duì)數(shù)非常的直觀，如果一個(gè)細(xì)胞中的唯一比對(duì)太少料身，顯然在后續(xù)分析中也沒(méi)有太多可用價(jià)值汤纸，可以直接剔除掉。
第二個(gè)參數(shù)惯驼，信噪比是根據(jù)TSS富集分?jǐn)?shù)進(jìn)行計(jì)算蹲嚣。傳統(tǒng)的混池ATAC-seq分析中递瑰，會(huì)使用TSS富集得分作為標(biāo)準(zhǔn)流程的一部分，用于確定信噪比（如 ENCODE計(jì)劃）隙畜。我們和其他人通過(guò)混池ATAC-seq和scATAC-seq分析發(fā)現(xiàn)抖部，TSS富集得分在大部分的細(xì)胞類型中都具有代表性。TSS富集得分的背后思想是议惰，由于大蛋白復(fù)合體會(huì)結(jié)合在啟動(dòng)子區(qū)域慎颗，ATAC-seq數(shù)據(jù)更多的富集在基因的TSS區(qū)域而不基因組其他區(qū)域。通過(guò)檢查這些TSS區(qū)域中心的開(kāi)放水平言询，我們發(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è)堿基的開(kāi)放性夫啊，之后這個(gè)譜會(huì)被用于確定TSS富集得分。在scATAC-seq中通過(guò)該方法計(jì)算每個(gè)細(xì)胞的TSS富集得分速度較慢辆憔，并且需要很強(qiáng)的算力撇眯。為了提高計(jì)算效率，同時(shí)還能得到和傳統(tǒng)計(jì)算接近的結(jié)果虱咧，我們以TSS位置為中心熊榛，以50-bp作為分窗計(jì)算兩邊的平均開(kāi)放程度，之后該值除以TSS兩側(cè)位置(+-1900-200bp)的平均開(kāi)放程度腕巡。 這種計(jì)算方法和原來(lái)相對(duì)比玄坦，兩者的相關(guān)性大于0.99，并且結(jié)果幾乎一樣绘沉。
第三個(gè)參數(shù)fragment大小分布并不是特別重要煎楣，最好是人工檢查下。由于DNA纏繞核小體的模式车伞，我們預(yù)期在fragment大小分布中看到核小體的周期性分布转质。這些山峰和低谷的出現(xiàn)正是由于DNA纏繞0，1帖世，2個(gè)核小體的結(jié)果(Tn5無(wú)法切割纏繞在核小體的DNA序列休蟹，只能切割兩邊）
ArchR默認(rèn)會(huì)過(guò)濾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, 虛線表示過(guò)濾閾值。第二圖注釋fragment大小分布圖珊膜。
對(duì)我們的教程數(shù)據(jù)容握，我們的三個(gè)樣本表現(xiàn)如下
BMMC:
CD34+ BMMC:
PBMC:
我們現(xiàn)在整理完畢我們的Arrow文件，接下來(lái)就是創(chuàng)建ArchRProject了车柠。
轉(zhuǎn)載來(lái)自：
作者：xuzhougeng
鏈接：http://www.reibang.com/p/f7975da8e1e8
來(lái)源：簡(jiǎn)書(shū)
著作權(quán)歸作者所有剔氏。商業(yè)轉(zhuǎn)載請(qǐng)聯(lián)系作者獲得授權(quán)，非商業(yè)轉(zhuǎn)載請(qǐng)注明出處竹祷。
