scATAC分析神器ArchR初探-使用ArchR-peak-calling(10)

scATAC分析神器ArchR初探-簡介(1)
scATAC分析神器ArchR初探-ArchR進行doublet處理(2)
scATAC分析神器ArchR初探-創(chuàng)建ArchRProject(3)
scATAC分析神器ArchR初探-使用ArchR降維(4)
scATAC分析神器ArchR初探--使用ArchR進行聚類(5)
scATAC分析神器ArchR初探-單細胞嵌入(6)
scATAC分析神器ArchR初探-使用ArchR計算基因活性值和標記基因(7)
scATAC分析神器ArchR初探-scRNA-seq確定細胞類型(8)
scATAC分析神器ArchR初探-ArchR中的偽批次重復處理(9)
scATAC分析神器ArchR初探-使用ArchR-peak-calling(10)
scATAC分析神器ArchR初探-使用ArchR識別標記峰(11)
scATAC分析神器ArchR初探-使用ArchR進行主題和功能豐富(12)
scATAC分析神器ArchR初探-利用ArchR豐富ChromVAR偏差(13)
scATAC分析神器ArchR初探-使用ArchR進行足跡(14)
scATAC分析神器ArchR初探-使用ArchR進行整合分析(15)
scATAC分析神器ArchR初探-使用ArchR進行軌跡分析(16)

10-使用ArchR-peak-calling

peak-calling是ATAC-seq數(shù)據(jù)分析中最基本的過程之一。因為每單元scATAC-seq數(shù)據(jù)本質上是二進制的(可訪問或不可訪問),所以我們不能在單個單元上peak-calling。因此,我們在上一章中定義了單元格組,通常是群集。此外,我們創(chuàng)建了偽大量復制品,以使我們能夠評估峰調用的可重復性。

10.1迭代重疊峰合并過程

我們首先在2018年引入了迭代重疊峰合并的策略。其他峰值合并策略也遇到一些我們在下面概述的關鍵問題。

10.1.1固定寬度與可變寬度峰

我們使用501-bp的固定寬度峰,因為它們不需要歸一化峰長,因此使下游計算更容易。此外,ATAC-seq中的絕大多數(shù)峰寬都小于501-bp。使用可變寬度的峰也很難合并來自多個樣本的峰調用。通常,我們認為使用可變寬度峰帶來的潛在收益并不超過成本。更廣泛地說,大多數(shù)分析相對于使用的峰集或峰樣式是穩(wěn)定的。

下面,我們使用具有幾個不同峰的幾種細胞類型的相同玩具示例,來說明這些常用峰合并方法之間的差異。

10.1.2使用bedtools合并原始峰重疊

原始峰重疊包括將任何彼此重疊的峰合并為一個更大的峰。在此方案中,菊花鏈成為一個大問題,因為彼此不直接重疊的峰包含在同一較大的峰中,因為它們被共享的內部峰橋接。這種方法的另一個問題是,如果您要跟蹤峰頂,則必須為每個新的合并峰選擇一個新峰,或者跟蹤適用于每個新峰的所有峰。 。通常,這種峰合并方法是使用bedtools merge命令實現(xiàn)的。

10.1.3使用bedtools群集的群集重疊

聚集的重疊峰聚集在一起,并選出一個獲勝者。這通常是通過使用bedtools cluster命令并在每個群集中保留最重要的峰值來完成的。根據(jù)我們的經驗,這最終會到達呼入峰,而錯過附近的較小峰。


10.1.4 ArchR中的迭代重疊

迭代重疊刪除可避免上述問題。峰首先按其重要性排序。保留最高有效峰,并將與最高有效峰直接重疊的任何峰從進一步分析中除去。然后,在剩余的峰中,重復此過程,直到不再存在任何峰為止。這樣可以避免菊花鏈,并且仍然允許使用固定寬度的峰。


10.1.5高峰調用方式比較

比較所有這些方法產生的峰調用,可以直接看出最終峰集中的明顯差異。我們認為,迭代重疊峰合并過程可產生具有最少警告的最佳峰集。


10.1.6那么這一切在ArchR中如何工作?

迭代重疊峰合并過程以分層方式執(zhí)行,以最佳地保留特定于細胞類型的峰。

想象一下,您有3種細胞類型(A,B和C),每種細胞類型都有3個偽大量復制品。ArchR使用稱為addReproduciblePeakSet()執(zhí)行此迭代重疊峰合并過程。首先,ArchR會分別為每個偽批量復制調用峰。然后,ArchR將分析單個單元格類型的所有偽批量復制,一起執(zhí)行迭代重疊移除的第一次迭代。重要的是要注意,ArchR對峰使用歸一化的顯著性度量來比較在不同樣本中調用的峰的顯著性。這是因為報告的MACS2顯著性與測序深度成正比,因此峰的顯著性在樣品之間無法立即比較。在迭代重疊移除的第一次迭代之后,ArchR檢查以查看偽批量重復中每個峰的重現(xiàn)性,并且僅保留通過該峰指示的閾值的峰。reproducibility參數(shù)。在此過程結束時,我們將為A,B和C這3種像元中的每一種都有一個合并的峰集。

然后,我們將重復此過程以合并A,B和C峰集。為此,我們將不同單元格類型的峰顯著性重新歸一化,并執(zhí)行迭代重疊去除。最終結果是固定寬度峰的單個合并峰集。

10.1.7如果我不喜歡這種迭代的重疊峰合并過程,該怎么辦?

迭代重疊峰合并過程由ArchR通過來實現(xiàn),addReproduciblePeakSet()但您始終可以通過使用自己的峰集ArchRProj <- addPeakSet()。

10.2使用 Macs2調用峰

如上所述,我們使用該addReproduciblePeakSet()函數(shù)在ArchR中生成了可重現(xiàn)的峰集。默認情況下,ArchR嘗試使用MACS2調用峰;但是,ArchR還實現(xiàn)了自己的本機峰值調用程序,當無法安裝MACS2(例如,我們尚未在Windows上成功安裝MACS2)時,可以使用該調用程序-下一節(jié)將介紹這種替代的峰值調用方法。

要使用MACS2調用峰,ArchR必須能夠找到MACS2可執(zhí)行文件。首先,ArchR查找您的PATH環(huán)境變量。如果失敗,則ArchR嘗試確定您是否已使用pip或來安裝MACS2 pip3。如果這些都不成功,則ArchR放棄并提供錯誤消息。如果已安裝MACS2,但ArchR找不到它,則應addReproduciblePeakSet()通過pathToMacs2參數(shù)提供該功能的路徑。

pathToMacs2 <- findMacs2()

確定了通往MACS2的路徑后,我們便可以創(chuàng)建帶有MACS2的可重現(xiàn)合并峰集(約5-10分鐘)。為了避免單元數(shù)很少的偽批量復制產生偏差,我們可以通過peaksPerCell參數(shù)為每個單元調用的峰數(shù)上限提供一個截止值。這樣可以防止單元格很少的群集將大量低質量峰貢獻給合并的峰集。還有許多其他參數(shù)可以調整addReproduciblePeakSet()-嘗試?addReproduciblePeakSet獲取更多信息。

每個ArchRProject對象只能包含一個峰集。因此,我們將輸出分配給addReproduciblePeakSet()期望的ArchRProject。如果要嘗試使用不同的峰集,則必須保存您的副本,ArchRProject從而也復制Arrow文件。盡管這確實占用了更多的磁盤存儲空間,但是鑒于Arrow文件的結構以及在Arrow文件中存儲峰矩陣信息,這是不可避免的。

projHeme4 <- addReproduciblePeakSet(
    ArchRProj = projHeme4, 
    groupBy = "Clusters2", 
    pathToMacs2 = pathToMacs2
)

要檢索此峰集作為GRanges對象,我們使用getPeakSet()函數(shù)。此峰集包含每個峰起源的組的注釋。但是,這些注釋并不固有地意味著僅在該組中調用給定的峰,而是帶注釋的組對該峰調用具有最高的歸一化意義。

getPeakSet(projHeme4)
10.3使用 TileMatrix調用峰

如前所述,ArchR還實現(xiàn)了其自己的本機峰值調用程序。盡管我們已針對MACS2對該峰值調用者進行了基準測試,并注意到性能非常相似,但除非絕對必要,否則我們不建議您使用此本地峰值調用者。

ArchR本機峰值調用者在500 bp處調用峰TileMatrix,我們表明addReproduciblePeakSet()我們想通過peakMethod參數(shù)使用該峰值調用者。請注意,我們沒有將輸出存儲到projHeme4對象中,因為我們不打算保持此峰值不變,并且此分析僅用于說明目的。存儲到ArchRProject對象中將覆蓋已經存儲在中的先前可重現(xiàn)峰集projHeme4。

projHemeTmp <- addReproduciblePeakSet(
    ArchRProj = projHeme4, 
    groupBy = "Clusters2",
    peakMethod = "Tiles",
    method = "p"
)

我們可以類似地檢查此合并的峰集。

getPeakSet(projHemeTmp)
10.4添加峰矩陣

現(xiàn)在,我們可以projHeme4使用saveArchRProject()函數(shù)保存原始圖像。這ArchRProject包含MACS2衍生的合并峰集。

saveArchRProject(ArchRProj = projHeme4, outputDirectory = "Save-ProjHeme4", load = FALSE)

要為下游分析做準備,我們可以創(chuàng)建一個新的ArchRProject名為projHeme5并添加新的基質含有成為我們新合并的波峰集合內插入計數(shù)。

projHeme5 <- addPeakMatrix(projHeme4)

現(xiàn)在我們可以看到已將一個新矩陣添加到projHeme5“ PeakMatrix”中。這是另一個類似于GeneScoreMatrix和的保留名稱矩陣TileMatrix。如前所述,每個ArchRProject對象只能有一個峰集和一個峰集PeakMatrix。當然,您可以創(chuàng)建數(shù)量不限的不同名稱的自定義特征矩陣,但是PeakMatrix保留給從存儲在中的峰集派生的插入計數(shù)矩陣ArchRProject。

getAvailableMatrices(projHeme5)
參考材料:

https://www.archrproject.com/

最后編輯于
?著作權歸作者所有,轉載或內容合作請聯(lián)系作者
【社區(qū)內容提示】社區(qū)部分內容疑似由AI輔助生成,瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明:文章內容(如有圖片或視頻亦包括在內)由作者上傳并發(fā)布,文章內容僅代表作者本人觀點,簡書系信息發(fā)布平臺,僅提供信息存儲服務。

友情鏈接更多精彩內容