系列回顧：
ArchR官網(wǎng)教程學(xué)習(xí)筆記1：Getting Started with ArchR
ArchR官網(wǎng)教程學(xué)習(xí)筆記2:基于ArchR推測(cè)Doublet
ArchR官網(wǎng)教程學(xué)習(xí)筆記3:創(chuàng)建ArchRProject
ArchR官網(wǎng)教程學(xué)習(xí)筆記4:ArchR的降維
ArchR官網(wǎng)教程學(xué)習(xí)筆記5:ArchR的聚類
ArchR官網(wǎng)教程學(xué)習(xí)筆記6:單細(xì)胞嵌入（Single-cell Embeddings）
ArchR官網(wǎng)教程學(xué)習(xí)筆記7:ArchR的基因評(píng)分和Marker基因
ArchR官網(wǎng)教程學(xué)習(xí)筆記8:定義與scRNA-seq一致的聚類
ArchR官網(wǎng)教程學(xué)習(xí)筆記9:ArchR的偽批量重復(fù)
call peak是ATAC-seq數(shù)據(jù)分析中最基本的過(guò)程之一。因?yàn)槊總€(gè)細(xì)胞的scATAC-seq數(shù)據(jù)本質(zhì)上是二進(jìn)制的(可接近或不可接近)，所以我們不能在單個(gè)細(xì)胞的基礎(chǔ)上call peak。出于這個(gè)原因，我們?cè)谇耙徽轮卸x了細(xì)胞群，通常是clusters。此外，我們創(chuàng)建了偽批量重復(fù)，以允許我們?cè)u(píng)估call peak的再現(xiàn)性。
（一）迭代重疊峰合并
我們首先要介紹的是 a strategy for iterative overlap peak merging in 2018. 其他合并peak的策略里有些關(guān)鍵問(wèn)題我們也會(huì)在下面談到。
（1）固定寬度的峰 vs 可變寬度的峰
我們使用501-bp固定寬度的峰，因?yàn)樗鼈兪瓜掠斡?jì)算更加容易，因?yàn)榉彘L(zhǎng)度不需要標(biāo)準(zhǔn)化。此外，ATAC-seq的絕大多數(shù)峰寬度小于501-bp。使用可變寬度的峰也使得從多個(gè)樣本合并call peak變得困難。一般來(lái)說(shuō)，我們并不認(rèn)為使用可變寬度的峰帶來(lái)的潛在好處大于其成本。更廣泛地說(shuō)，大多數(shù)分析對(duì)于使用的峰集或峰樣式是穩(wěn)定的。
下面，我們使用具有不同峰的幾個(gè)細(xì)胞類型的示例來(lái)說(shuō)明這些峰合并方法之間的差異。
（2）使用bedtools merge合并原始peak overlap
原始峰重疊（Raw peak overlap）是指取任何相互重疊的峰，并將它們合并成一個(gè)更大的峰。在這個(gè)方案中，daisy-chaining成為了一個(gè)大問(wèn)題，因?yàn)闆](méi)有直接重疊的峰被包含在同一個(gè)更大的峰中，因?yàn)樗鼈儽灰粋€(gè)共享的內(nèi)部峰連接。這種方法的另一個(gè)問(wèn)題是，如果你想要跟蹤峰的頂點(diǎn)，你要么為每個(gè)新的合并的峰選擇一個(gè)新的頂點(diǎn)，要么用每個(gè)峰的所有頂點(diǎn)。通常，這種峰合并方法是使用bedtools合并命令實(shí)現(xiàn)的。
（3）使用bedtools cluster重疊clusters
clusters重疊是指從聚集在一起peaks中選出一個(gè)“贏家”。通常使用bedtools cluster命令來(lái)完成，然后在每個(gè)cluster中保留最顯著的峰。根據(jù)我們的經(jīng)驗(yàn)，這會(huì)導(dǎo)致忽略附近較小的峰。
（4）ArchR中的迭代重疊
迭代重疊消除避免了上述問(wèn)題。峰的重要性排在第一位。保留最顯著峰，并將與最顯著峰直接重疊的任何峰從進(jìn)一步分析中去除。然后，對(duì)于剩余的峰，重復(fù)這個(gè)過(guò)程，直到不再存在多于的峰為止。這就避免了daisy-chaining，并且仍然允許使用固定寬度的峰。
（5）比較peak calling的方法
比較所有這些peak calling的方法，可以直接顯示最終峰集的明顯差異。我們的觀點(diǎn)是，迭代重疊峰合并過(guò)程產(chǎn)生了最好的峰集。
（6）在ArchR中這些是如何工作的呢？
迭代重疊峰合并過(guò)程以分層的方式執(zhí)行，以最佳地保留細(xì)胞類型特異的峰。
想象這樣一種情況，你有3種細(xì)胞類型，A, B和C，每種細(xì)胞類型有3次偽批量重復(fù)。ArchR使用一個(gè)名為addReproduciblePeakSet()的函數(shù)來(lái)執(zhí)行這個(gè)迭代重疊峰值合并過(guò)程。首先，ArchR將分別為每個(gè)偽批量重復(fù)call peak。然后，ArchR將分析來(lái)自單個(gè)細(xì)胞類型的所有偽批量重復(fù)，執(zhí)行迭代重疊去除的第一次迭代。值得注意的是，ArchR使用峰的標(biāo)準(zhǔn)化顯著性度量來(lái)比較不同樣本間peak的顯著性。這是因?yàn)榈腗ACS2顯著性與測(cè)序深度成比例，所以峰顯著性不能立即在樣本間進(jìn)行比較。在第一次迭代之后，ArchR檢查偽批量重復(fù)中每個(gè)峰的重現(xiàn)性，并且只保留那些通過(guò)重現(xiàn)性參數(shù)閾值的峰。在這個(gè)過(guò)程的最后，我們將為3種細(xì)胞類型(A、B和C)分別設(shè)置一個(gè)合并的峰集。
然后，我們重復(fù)這個(gè)步驟來(lái)合并A、B和C的峰集。為了做到這一點(diǎn)，我們重新normalize不同細(xì)胞類型的峰的顯著性，并執(zhí)行迭代重疊去除。這樣做的最終結(jié)果是一個(gè)固定寬度峰的合并峰集。
（7）如果我不喜歡這個(gè)迭代重疊峰合并過(guò)程呢?
迭代重疊峰合并過(guò)程是由ArchR通過(guò)addReproduciblePeakSet()實(shí)現(xiàn)的，但是你也可以通過(guò)ArchRProj <- addPeakSet()使用你自己的峰集。
（二）使用Macs2進(jìn)行call peak
如上所述，我們使用addReproduciblePeakSet()函數(shù)在ArchR中生成一個(gè)可重復(fù)的峰集。默認(rèn)情況下，ArchR嘗試使用MACS2 call peaks。然而，ArchR也實(shí)現(xiàn)了它自己的本地peak caller，它可以在MACS2無(wú)法安裝時(shí)使用(例如，我們沒(méi)有成功地在Windows上安裝MACS2)——這種備選call peaks方法將在下一節(jié)中描述。
要使用MACS2 call peaks，ArchR必須能夠找到MACS2可執(zhí)行文件。首先，ArchR查看PATH環(huán)境變量。如果不成功，ArchR將嘗試確定你是否安裝了帶有pip或pip3的MACS2。如果這兩個(gè)都不成功，ArchR將放棄并提供一個(gè)錯(cuò)誤消息。如果你安裝了MACS2, ArchR找不到它，你應(yīng)該通過(guò)pathToMacs2參數(shù)給addReproduciblePeakSet()函數(shù)提供路徑。
這一步我是在Linux系統(tǒng)里運(yùn)行的，因?yàn)槲业膚indows系統(tǒng)實(shí)在是不給力，無(wú)法找到macs2的路徑，所以我只好在linux系統(tǒng)里重新運(yùn)行了一遍，包括前面的所有章節(jié)。
> pathToMacs2 <- findMacs2() #find the path of macs2
Searching For MACS2..
Not Found in $PATH
Not Found with pip
Not Found with pip3
Error in findMacs2() :
Could Not Find Macs2! Please install w/ pip, add to your $PATH, or just supply the macs2 path directly and avoid this function!
#這里顯示找不到MACS2路徑，實(shí)際上我是安裝了MACS2的
如果ArchR找不到路徑，你可以指定路徑給它：
> projHeme4 <- addReproduciblePeakSet(
ArchRProj = projHeme4,
groupBy = "Clusters2",
pathToMacs2 = "/home/yanfang/Downloads/Anaconda/Anaconda_install/bin/macs2"
)
#運(yùn)行成功會(huì)有下面的提示：
ArchR logging to : ArchRLogs/ArchR-addReproduciblePeakSet-58eb1b1ff90-Date-2020-11-24_Time-23-00-04.log
If there is an issue, please report to github with logFile!
2020-11-24 23:00:05 : Peak Calling Parameters!, 0.01 mins elapsed.
Group nCells nCellsUsed nReplicates nMin nMax maxPeaks
B B 432 428 2 174 254 150000
CD4.M CD4.M 639 610 2 110 500 150000
CD4.N CD4.N 1279 552 2 52 500 150000
CLP CLP 384 384 2 84 300 150000
Erythroid Erythroid 732 666 2 166 500 150000
GMP GMP 1201 849 2 349 500 150000
Mono Mono 2579 1000 2 500 500 150000
NK NK 874 786 2 286 500 150000
pDC pDC 300 290 2 137 153 145000
PreB PreB 358 358 2 40 318 150000
Progenitor Progenitor 1473 658 2 158 500 150000
2020-11-24 23:00:06 : Batching Peak Calls!, 0.011 mins elapsed.
2020-11-24 23:00:06 : Batch Execution w/ safelapply!, 0 mins elapsed.
2020-11-24 23:10:39 : Identifying Reproducible Peaks!, 10.577 mins elapsed.
2020-11-24 23:11:09 : Creating Union Peak Set!, 11.073 mins elapsed.
Converged after 6 iterations!
Plotting Ggplot!
2020-11-24 23:11:14 : Finished Creating Union Peak Set (140865)!, 11.147 mins elapsed.
找到MACS2的路徑后，我們就可以用MACS2創(chuàng)建一個(gè)可重復(fù)的合并峰集。為了避免細(xì)胞數(shù)量很少的偽批量重復(fù)產(chǎn)生的偏差，我們可以通過(guò)peaksPerCell參數(shù)為每個(gè)細(xì)胞call的峰數(shù)量的上限提供一個(gè)cutoff。這可以防止細(xì)胞數(shù)量很少的clusters里大量低質(zhì)量的峰影響合并峰集。還有許多其他參數(shù)可以在addReproduciblePeakSet()中調(diào)整，請(qǐng)嘗試?addReproduciblePeakSet以獲得更多信息。
每個(gè)ArchRProject對(duì)象只能包含一個(gè)峰集。因此，我們將addReproduciblePeakSet()的輸出分配給我們想要的ArchRProject。如果您想嘗試使用不同的峰集，則必須保存ArchRProject的一個(gè)副本，從而也復(fù)制Arrow files。雖然這會(huì)使用更多的磁盤存儲(chǔ)空間，但考慮到Arrow files的結(jié)構(gòu)和Arrow files中峰值矩陣信息的存儲(chǔ)，這是不可避免的。
要將這個(gè)峰集作為GRanges對(duì)象，我們使用getPeakSet()函數(shù)。這個(gè)峰集包含每個(gè)峰產(chǎn)生的組的注釋。但是，這些注釋并不意味著只在該組中被call，而是說(shuō)明注釋的組對(duì)該call peak具有最高的標(biāo)準(zhǔn)化意義。
> getPeakSet(projHeme4)
GRanges object with 140865 ranges and 12 metadata columns:
seqnames ranges strand | score replicateScoreQuantile
<Rle> <IRanges> <Rle> | <numeric> <numeric>
Mono chr1 752495-752995 * | 27.3149 0.865
CD4.N chr1 757871-758371 * | 5.55037 0.835
CD4.N chr1 762690-763190 * | 17.9242 0.977
GMP chr1 773670-774170 * | 4.90925 0.513
B chr1 801006-801506 * | 15.2793 0.685
... ... ... ... . ... ...
NK chrX 154807254-154807754 * | 4.84969 0.403
Mono chrX 154840907-154841407 * | 3.66286 0.376
Mono chrX 154841994-154842494 * | 12.6809 0.747
NK chrX 154862043-154862543 * | 13.8961 0.525
GMP chrX 154996991-154997491 * | 4.90925 0.513
groupScoreQuantile Reproducibility GroupReplicate distToGeneStart
<numeric> <numeric> <character> <integer>
Mono 0.737 2 Mono._.scATAC_BMMC_R1 10156
CD4.N 0.541 2 CD4.N._.scATAC_PBMC_R1 4780
CD4.N 0.917 2 CD4.N._.scATAC_PBMC_R1 30
GMP 0.176 2 GMP._.scATAC_CD34_BMMC_R1 10948
B 0.418 2 B._.scATAC_PBMC_R1 10925
... ... ... ... ...
NK 0.13 2 NK._.scATAC_BMMC_R1 35117
Mono 0.081 2 Mono._.scATAC_PBMC_R1 1464
Mono 0.512 2 Mono._.scATAC_BMMC_R1 377
NK 0.253 2 NK._.scATAC_PBMC_R1 19670
GMP 0.176 2 GMP._.scATAC_CD34_BMMC_R1 154618
nearestGene peakType distToTSS nearestTSS GC idx
<character> <character> <integer> <character> <numeric> <integer>
Mono LINC00115 Distal 10156 uc001aau.3 0.485 1
CD4.N LINC00115 Distal 4780 uc001aau.3 0.5609 2
CD4.N LINC01128 Promoter 30 uc021oeh.1 0.6966 3
GMP LINC01128 Intronic 10741 uc021oeh.1 0.4571 4
B FAM41C Distal 10925 uc021oei.1 0.4371 5
... ... ... ... ... ... ...
NK TMLHE Intronic 35117 uc004cin.3 0.499 3405
Mono TMLHE Intronic 1464 uc004cin.3 0.4391 3406
Mono TMLHE Intronic 377 uc004cin.3 0.6148 3407
NK TMLHE Distal 19670 uc004cin.3 0.4251 3408
GMP TMLHE Distal 209 uc004cin.3 0.5449 3409
-------
seqinfo: 23 sequences from an unspecified genome; no seqlengths
（三）使用TileMatrix進(jìn)行call peak
如前所述，ArchR還實(shí)現(xiàn)了它自己的peak caller。雖然我們已經(jīng)用MACS2對(duì)這個(gè)peak caller進(jìn)行了基準(zhǔn)測(cè)試，并且性能非常相似，但除非絕對(duì)必要，我們不建議使用這個(gè)peak caller。
ArchR自身的peak caller call的是500-bp TileMatrix上的峰，我們通過(guò)peakMethod參數(shù)向addReproduciblePeakSet()說(shuō)明我們希望使用這個(gè)peak caller。請(qǐng)注意，我們沒(méi)有將輸出存儲(chǔ)到projHeme4對(duì)象中，因?yàn)槲覀儾淮蛩惚Ａ暨@個(gè)峰集，并且此分析僅用于說(shuō)明目的。存儲(chǔ)到ArchRProject對(duì)象中會(huì)覆蓋之前已經(jīng)存儲(chǔ)在projHeme4中的峰集。
> projHemeTmp <- addReproduciblePeakSet(
ArchRProj = projHeme4,
groupBy = "Clusters2",
peakMethod = "Tiles",
method = "p"
)
> getPeakSet(projHemeTmp)
GRanges object with 271847 ranges and 9 metadata columns:
seqnames ranges strand | mlog10p Group distToGeneStart
<Rle> <IRanges> <Rle> | <numeric> <Rle> <integer>
[1] chr1 752500-752999 * | 6.939 Mono 10152
[2] chr1 758000-758499 * | 3.185 CD4.N 4652
[3] chr1 758500-758999 * | 1.03 NK 4152
[4] chr1 762000-762499 * | 2.409 NK 652
[5] chr1 762500-762999 * | 23.286 NK 152
... ... ... ... . ... ... ...
[271843] chrX 154842000-154842499 * | 6.121 CLP 372
[271844] chrX 154842500-154842999 * | 5.545 CD4.M 126
[271845] chrX 154862000-154862499 * | 2.325 NK 19626
[271846] chrX 154862500-154862999 * | 1.362 Progenitor 20126
[271847] chrX 154997000-154997499 * | 2.075 Progenitor 154626
nearestGene peakType distToTSS nearestTSS GC idx
<character> <character> <integer> <character> <numeric> <integer>
[1] LINC00115 Distal 10152 uc001aau.3 0.484 1
[2] LINC00115 Distal 4652 uc001aau.3 0.552 2
[3] LINC00115 Distal 4152 uc001aau.3 0.56 3
[4] LINC00115 Promoter 652 uc001aau.3 0.574 4
[5] LINC00115 Promoter 152 uc001aau.3 0.684 5
... ... ... ... ... ... ...
[271843] TMLHE Intronic 372 uc004cin.3 0.612 6496
[271844] TMLHE Promoter 126 uc004cin.3 0.554 6497
[271845] TMLHE Distal 19626 uc004cin.3 0.43 6498
[271846] TMLHE Distal 20126 uc004cin.3 0.42 6499
[271847] TMLHE Distal 201 uc004cin.3 0.542 6500
-------
seqinfo: 23 sequences from an unspecified genome; no seqlengths
（四）比較兩種peak calling的方法
將上面的兩種peak calling進(jìn)行比較。
首先，我們檢查macs2的call peak與tilematrix的call peak重疊的百分比:
> length(subsetByOverlaps(getPeakSet(projHeme4), getPeakSet(projHemeTmp))) / length(getPeakSet(projHeme4))
[1] 0.9710006
相反，我們看tilematrix的call peak與macs2的call peak重疊的百分比：
> length(subsetByOverlaps(getPeakSet(projHemeTmp), getPeakSet(projHeme4))) / length(getPeakSet(projHemeTmp))
[1] 0.7472328
如果我們將峰的寬度擴(kuò)大(1000-bp的峰而不是500-bp的峰)，macs2的峰重疊的百分比不會(huì)發(fā)生太大變化：
> length(subsetByOverlaps(resize(getPeakSet(projHeme4), 1000, "center"), getPeakSet(projHemeTmp))) / length(getPeakSet(projHeme4))
[1] 0.9752671
但是tilematrix的call peak與MACS2重疊的比例增加了：
> length(subsetByOverlaps(getPeakSet(projHemeTmp), resize(getPeakSet(projHeme4), 1000, "center"))) / length(getPeakSet(projHemeTmp))
[1] 0.8220985
（五）添加Peak矩陣
現(xiàn)在我們可以保存projHeme4了。這個(gè)ArchRProject包含了MACS2合并峰集。
> saveArchRProject(ArchRProj = projHeme4, outputDirectory = "Save-ProjHeme4", load = FALSE)
為了為下游分析做準(zhǔn)備，我們創(chuàng)建一個(gè)名為projHeme5的新ArchRProject，并向其添加一個(gè)新的矩陣，其中包含在新的合并峰集中的insertion counts。
> projHeme5 <- addPeakMatrix(projHeme4)
現(xiàn)在我們可以看到，projHeme5中添加了一個(gè)名為“PeakMatrix”的新矩陣。這是另一個(gè)保留名稱的矩陣，類似于GeneScoreMatrix和TileMatrix。如前所述，每個(gè)ArchRProject對(duì)象只能有一個(gè)峰集和一個(gè)峰矩陣。當(dāng)然，你可以創(chuàng)建數(shù)量不限、名稱不同的自定義特征矩陣，但PeakMatrix是保留給ArchRProject中存儲(chǔ)的峰集派生的insertion counts矩陣的。
> getAvailableMatrices(projHeme5)
[1] "GeneIntegrationMatrix" "GeneScoreMatrix" "PeakMatrix"
[4] "TileMatrix"
