scATAC分析神器ArchR初探-簡(jiǎn)介（1）
scATAC分析神器ArchR初探-ArchR進(jìn)行doublet處理（2）
scATAC分析神器ArchR初探-創(chuàng)建ArchRProject（3）
scATAC分析神器ArchR初探-使用ArchR降維（4）
scATAC分析神器ArchR初探--使用ArchR進(jìn)行聚類（5）
scATAC分析神器ArchR初探-單細(xì)胞嵌入（6）
scATAC分析神器ArchR初探-使用ArchR計(jì)算基因活性值和標(biāo)記基因（7）
scATAC分析神器ArchR初探-scRNA-seq確定細(xì)胞類型（8）
scATAC分析神器ArchR初探-ArchR中的偽批次重復(fù)處理（9）
scATAC分析神器ArchR初探-使用ArchR-peak-calling（10）
scATAC分析神器ArchR初探-使用ArchR識(shí)別標(biāo)記峰（11）
scATAC分析神器ArchR初探-使用ArchR進(jìn)行主題和功能豐富(12)
scATAC分析神器ArchR初探-利用ArchR豐富ChromVAR偏差(13)
scATAC分析神器ArchR初探-使用ArchR進(jìn)行足跡（14）
scATAC分析神器ArchR初探-使用ArchR進(jìn)行整合分析（15）
scATAC分析神器ArchR初探-使用ArchR進(jìn)行軌跡分析（16）
15-使用ArchR進(jìn)行整合分析
ArchR的主要優(yōu)勢(shì)之一是它能夠集成多個(gè)級(jí)別的信息以提供新穎的見(jiàn)解驮吱。這可以采取僅ATAC-seq的分析形式，例如識(shí)別峰的可共同接近性以預(yù)測(cè)調(diào)節(jié)相互作用萧吠，或采用整合scRNA-seq數(shù)據(jù)的分析左冬，例如通過(guò)峰-基因連鎖分析預(yù)測(cè)增強(qiáng)子活性。無(wú)論哪種情況纸型，ArchR都可以輕松地從scATAC-seq數(shù)據(jù)中獲得更深刻的見(jiàn)解拇砰。
15.1創(chuàng)建單元格的低重疊聚合
ArchR促進(jìn)了許多涉及特征相關(guān)性的綜合分析。在稀疏的單細(xì)胞數(shù)據(jù)中執(zhí)行這些計(jì)算會(huì)導(dǎo)致這些相關(guān)分析中的大量噪聲绊袋。為了克服這一挑戰(zhàn)毕匀，我們采用了Cicero引入的方法在這些分析之前創(chuàng)建單個(gè)細(xì)胞的低重疊聚集體铸鹰。為了減少偏差癌别，我們會(huì)過(guò)濾與其他任何聚合具有大于80％重疊的聚合。為了提高這種方法的速度蹋笼，我們使用“ Rcpp”包開(kāi)發(fā)了一種優(yōu)化的迭代重疊檢查例程的實(shí)現(xiàn)展姐，以及C ++中快速特征關(guān)聯(lián)的實(shí)現(xiàn)。這些優(yōu)化的方法用于ArchR中剖毯，用于計(jì)算峰的可及性圾笨，峰與基因的鏈接以及其他鏈接分析。這些低重疊的聚合的使用是在引擎蓋下進(jìn)行的逊谋，但是為了清楚起見(jiàn)擂达，我們?cè)诖颂峒啊?/p>
15.2與ArchR的可共訪問(wèn)性
共可及性是許多單個(gè)細(xì)胞中兩個(gè)峰之間可及性的關(guān)聯(lián)。換句話說(shuō)胶滋，當(dāng)在單個(gè)單元中可以訪問(wèn)峰A時(shí)板鬓，通常也可以訪問(wèn)峰B。我們?cè)谙旅嬉砸曈X(jué)方式說(shuō)明了這一概念究恤，表明Enhancer E3通常與Promoter P可以同時(shí)訪問(wèn)俭令。
關(guān)于協(xié)同訪問(wèn)分析，要注意的一件事是部宿，它經(jīng)常將特定于細(xì)胞類型的峰鑒定為協(xié)同訪問(wèn)。這是因?yàn)檫@些峰通常在單個(gè)單元格類型中都可以一起訪問(wèn)棉胀，而在所有其他單元格類型中通常都不能訪問(wèn)斤葱。這會(huì)產(chǎn)生很強(qiáng)的相關(guān)性，但不一定意味著這些峰之間存在調(diào)節(jié)關(guān)系绵患。
要計(jì)算ArchR中的可訪問(wèn)性，我們使用addCoAccessibility()將峰值可訪問(wèn)性信息存儲(chǔ)在中的功能ArchRProject悟耘。
projHeme5 <- addCoAccessibility(
ArchRProj = projHeme5,
reducedDims = "IterativeLSI"
)
我們可以ArchRProject通過(guò)getCoAccessibility()函數(shù)從中檢索此可訪問(wèn)性信息藏雏，該函數(shù)返回DataFrame對(duì)象if returnLoops = FALSE。
cA <- getCoAccessibility(
ArchRProj = projHeme5,
corCutOff = 0.5,
resolution = 1,
returnLoops = FALSE
)
其中DataFrame包含一些重要的信息作煌。的queryHits和subjectHits列表示被發(fā)現(xiàn)相關(guān)聯(lián)掘殴，這兩個(gè)峰的索引。該correlation列給出了這兩個(gè)峰之間可及性的數(shù)字相關(guān)性粟誓。
cA
這種可共訪問(wèn)性DataFrame還具有包含GRanges相關(guān)峰對(duì)象的元數(shù)據(jù)組件奏寨。上面提到的queryHits和的索引subjectHits適用于此GRanges對(duì)象。
metadata(cA)[[1]]
如果我們?cè)O(shè)置returnLoops = TRUE鹰服，則將getCoAccessibility()以循環(huán)跟蹤的形式返回可訪問(wèn)性信息病瞳。在此GRanges對(duì)象中，IRanges每次互動(dòng)的圖譜的起點(diǎn)和終點(diǎn)分別指向兩個(gè)不同的可同時(shí)訪問(wèn)的峰悲酷。該resolution參數(shù)設(shè)置這些循環(huán)的堿基對(duì)分辨率套菜。當(dāng)為時(shí)resolution = 1，這將創(chuàng)建連接每個(gè)峰中心的循環(huán)设易。
cA <- getCoAccessibility(
ArchRProj = projHeme5,
corCutOff = 0.5,
resolution = 1,
returnLoops = TRUE
)
我們可以將此GRanges對(duì)象與DataFrame上面生成的對(duì)象進(jìn)行比較逗柴。
cA[[1]]
相反，如果我們降低到的循環(huán)的分辨率resolution = 1000顿肺，則可能有助于過(guò)度繪制協(xié)同訪問(wèn)交互戏溺。在下面，我們看到GRanges對(duì)象中的總條目比上面的少屠尊。
cA <- getCoAccessibility(
ArchRProj = projHeme5,
corCutOff = 0.5,
resolution = 1000,
returnLoops = TRUE
)
cA[[1]]
同樣旷祸，如果我們進(jìn)一步降低分辨率resolution = 10000，我們將確定更少的協(xié)同訪問(wèn)交互讼昆。
cA <- getCoAccessibility(
ArchRProj = projHeme5,
corCutOff = 0.5,
resolution = 10000,
returnLoops = TRUE
)
cA[[1]]
15.2.1繪制瀏覽器的可共通性軌跡
在向我們添加了可訪問(wèn)性信息后托享，我們ArchRProject就可以在繪制瀏覽器軌跡時(shí)將其用作循環(huán)軌跡。我們通過(guò)函數(shù)的loops參數(shù)來(lái)實(shí)現(xiàn)plotBrowserTrack()浸赫。在這里闰围，我們使用的是默認(rèn)參數(shù)getCoAccessibility()，其中包括corCutOff = 0.5掺炭，resolution = 1000辫诅，和returnLoops = TRUE。
markerGenes <- c(
"CD34", #Early Progenitor
"GATA1", #Erythroid
"PAX5", "MS4A1", #B-Cell Trajectory
"CD14", #Monocytes
"CD3D", "CD8A", "TBX21", "IL7R" #TCells
)
p <- plotBrowserTrack(
ArchRProj = projHeme5,
groupBy = "Clusters2",
geneSymbol = markerGenes,
upstream = 50000,
downstream = 50000,
loops = getCoAccessibility(projHeme5)
)
為了繪制瀏覽器軌跡涧狮，我們使用該grid.draw函數(shù)炕矮，并使用$訪問(wèn)器按名稱選擇特定的標(biāo)記基因么夫。
grid::grid.newpage()
grid::grid.draw(p$CD14)
要保存此圖的可編輯矢量化版本，請(qǐng)使用plotPDF()肤视。
plotPDF(plotList = p,
name = "Plot-Tracks-Marker-Genes-with-CoAccessibility.pdf",
ArchRProj = projHeme5,
addDOC = FALSE, width = 5, height = 5)
15.3 Peak2Gene與ArchR的鏈接
與協(xié)同訪問(wèn)相似档痪，ArchR也可以識(shí)別所謂的“峰到基因鏈接”。峰與基因之間的鏈接與共同訪問(wèn)之間的主要區(qū)別在于邢滑，共同訪問(wèn)是僅ATAC-seq的分析腐螟，用于尋找兩個(gè)峰之間的可訪問(wèn)性的相關(guān)性，而峰到基因的鏈接則利用整合的scRNA-seq數(shù)據(jù)尋找峰可及性與基因表達(dá)之間的相關(guān)性困后。這些代表解決類似問(wèn)題的正交方法乐纸。但是，由于峰與基因之間的聯(lián)系將scATAC-seq和scRNA-seq數(shù)據(jù)相關(guān)聯(lián)摇予，因此我們經(jīng)常認(rèn)為這些聯(lián)系與基因調(diào)控相互作用更為相關(guān)汽绢。
為了識(shí)別ArchR中的峰到基因鏈接，我們使用該addPeak2GeneLinks()函數(shù)侧戴。
projHeme5 <- addPeak2GeneLinks(
ArchRProj = projHeme5,
reducedDims = "IterativeLSI"
)
然后宁昭，我們可以使用類似于使用getPeak2GeneLinks()函數(shù)檢索共通性交互的方式來(lái)檢索這些峰-基因鏈接。如我們先前所見(jiàn)酗宋，此功能允許用戶指定鏈接的關(guān)聯(lián)和分辨率的截止值积仗。
p2g <- getPeak2GeneLinks(
ArchRProj = projHeme5,
corCutOff = 0.45,
resolution = 1,
returnLoops = FALSE
)
當(dāng)returnLoops設(shè)置為false時(shí)，此函數(shù)將返回與所返回的DataFrame對(duì)象類似的DataFrame對(duì)象getCoAccessibility()蜕猫。主要區(qū)別在于寂曹，scATAC-seq峰idxATAC的索引存儲(chǔ)在列中，而scRNA-seq基因的索引存儲(chǔ)在idxRNA列中丹锹。
p2g
這種峰到基因的鏈接DataFrame還具有包含GRanges相關(guān)峰對(duì)象的元數(shù)據(jù)組件稀颁。idxATAC上面提到的索引適用于此GRanges對(duì)象芬失。
metadata(p2g)[[1]]
如果設(shè)置returnLoops = TRUE楣黍，getPeak2GeneLinks()則將返回GRanges連接峰和基因的循環(huán)跟蹤對(duì)象。至于共通性棱烂，IRanges對(duì)象的開(kāi)始和結(jié)束代表峰的位置和被鏈接的基因租漂。當(dāng)為時(shí)resolution = 1，這會(huì)將峰的中心鏈接到基因的單堿基TSS颊糜。
p2g <- getPeak2GeneLinks(
ArchRProj = projHeme5,
corCutOff = 0.45,
resolution = 1,
returnLoops = TRUE
)
p2g[[1]]
我們可以通過(guò)設(shè)置來(lái)降低這些鏈接的分辨率resolution = 1000哩治。這在將鏈接繪制為瀏覽器軌跡時(shí)非常有用，因?yàn)樵谀承┣闆r下衬鱼，許多附近的峰都鏈接到同一基因业筏，這可能很難可視化。
p2g <- getPeak2GeneLinks(
ArchRProj = projHeme5,
corCutOff = 0.45,
resolution = 1000,
returnLoops = TRUE
)
p2g[[1]]
降低分辨率甚至?xí)M(jìn)一步減少所識(shí)別的峰到基因鏈接的總數(shù)鸟赫。
p2g <- getPeak2GeneLinks(
ArchRProj = projHeme5,
corCutOff = 0.45,
resolution = 10000,
returnLoops = TRUE
)
p2g[[1]]
15.3.1使用峰到基因鏈接繪制瀏覽器軌跡
要將這些峰到基因鏈接繪制為瀏覽器軌跡蒜胖，我們使用上一部分中所示的可共訪問(wèn)性相同的工作流程消别。在這里我們使用plotBrowserTrack()功能
markerGenes <- c(
"CD34", #Early Progenitor
"GATA1", #Erythroid
"PAX5", "MS4A1", #B-Cell Trajectory
"CD14", #Monocytes
"CD3D", "CD8A", "TBX21", "IL7R" #TCells
)
p <- plotBrowserTrack(
ArchRProj = projHeme5,
groupBy = "Clusters2",
geneSymbol = markerGenes,
upstream = 50000,
downstream = 50000,
loops = getPeak2GeneLinks(projHeme5)
)
為了繪制瀏覽器軌跡，我們使用該grid.draw函數(shù)台谢，并使用$訪問(wèn)器按名稱選擇特定的標(biāo)記基因寻狂。
grid::grid.newpage()
grid::grid.draw(p$CD14)
To save an editable vectorized version of this plot, we use plotPDF().
plotPDF(plotList = p,
name = "Plot-Tracks-Marker-Genes-with-Peak2GeneLinks.pdf",
ArchRProj = projHeme5,
addDOC = FALSE, width = 5, height = 5)
15.3.2 Plotting a heatmap of peak-to-gene links
To visualize the correspondence of all of our peak-to-gene links, we can plot a peak-to-gene heatmap which contains two side-by-side heatmaps, one for our scATAC-seq data and one for our scRNA-seq data. To do this, we use the plotPeak2GeneHeatmap()
p <- plotPeak2GeneHeatmap(ArchRProj = projHeme5, groupBy = "Clusters2")
The heatmap rows are clustered using k-means clustering based on the value passed to the parameter k, which defaults to 25 as shown below.
p
15.4確定正TF調(diào)節(jié)器
ATAC-seq可以無(wú)偏鑒定TF，這些TF在包含其DNA結(jié)合基序的位點(diǎn)上朋沮，染色質(zhì)可及性表現(xiàn)出很大的變化蛇券。但是，當(dāng)通過(guò)位置權(quán)重矩陣（PWM）進(jìn)行匯總時(shí)樊拓，TF家族（例如GATA因子）在其綁定圖案中具有相似的功能纠亚。
這種基序相似性使其難以鑒定可能導(dǎo)致觀察到的染色質(zhì)可及性在其預(yù)測(cè)的結(jié)合位點(diǎn)發(fā)生變化的特定TF。為了規(guī)避這一挑戰(zhàn)筋夏，我們先前使用了ATAC-seq和RNA-seq來(lái)鑒定其基因表達(dá)與其相應(yīng)基序可及性變化呈正相關(guān)的TF菜枷。我們將這些TF稱為“積極監(jiān)管者”。但是叁丧，該分析依賴于并非在所有實(shí)驗(yàn)中都容易獲得的匹配基因表達(dá)數(shù)據(jù)啤誊。為了克服這種依賴性，ArchR可以識(shí)別其推斷的基因得分與其chromVAR TF偏差z得分相關(guān)的TF拥娄。為此蚊锹，ArchR將TF基序的chromVAR偏差z分?jǐn)?shù)與來(lái)自低重疊細(xì)胞聚集體的TF基因的基因活性得分相關(guān)聯(lián)。將scRNA-seq與ArchR結(jié)合使用時(shí)稚瘾，
15.4.1第1步牡昆。確定偏差TF圖案
鑒定陽(yáng)性TF調(diào)節(jié)劑的第一部分是鑒定異常TF基序。我們?cè)谏弦徽轮羞M(jìn)行了此分析摊欠，MotifMatrix為所有圖案創(chuàng)建了一個(gè)chromVAR偏差和z分?jǐn)?shù)偏差丢烘。我們可以使用getGroupSE()返回a 的函數(shù)來(lái)獲得按聚類平均的數(shù)據(jù)SummarizedExperiment。
seGroupMotif <- getGroupSE(ArchRProj = projHeme5, useMatrix = "MotifMatrix", groupBy = "Clusters2")
由于此SummarizedExperiment對(duì)象來(lái)自些椒，因此MotifMatrix具有兩個(gè)序列名稱-“偏差”和“ z”-對(duì)應(yīng)于chromVAR的原始偏差和z分?jǐn)?shù)偏差播瞳。
seGroupMotif
我們可以將其子集化為SummarizedExperiment偏差z得分。
seZ <- seGroupMotif[rowData(seGroupMotif)$seqnames=="z",]
然后免糕，我們可以確定所有群集之間z得分的最大增量赢乓。這將有助于根據(jù)跨群集觀察到的變化程度對(duì)主題進(jìn)行分層。
rowData(seZ)$maxDelta <- lapply(seq_len(ncol(seZ)), function(x){
rowMaxs(assay(seZ) - assay(seZ)[,x])
}) %>% Reduce("cbind", .) %>% rowMaxs
15.4.2步驟2.識(shí)別相關(guān)的TF基序和TF基因得分/表達(dá)
為了鑒定其基序可及性與其自身基因活性相關(guān)的TF（通過(guò)基因評(píng)分或基因表達(dá)）石窑，我們使用該correlateMatrices()功能并提供我們感興趣的兩個(gè)矩陣牌芋，在這種情況下為GeneScoreMatrix和MotifMatrix。如前所述松逊，這些相關(guān)性是在reducedDims參數(shù)中指定的較低維空間中確定的許多低重疊單元聚合中確定的躺屁。
corGSM_MM <- correlateMatrices(
ArchRProj = projHeme5,
useMatrix1 = "GeneScoreMatrix",
useMatrix2 = "MotifMatrix",
reducedDims = "IterativeLSI"
)
此函數(shù)返回一個(gè)DataFrame對(duì)象，其中包含來(lái)自每個(gè)矩陣的元素以及低重疊單元格聚合之間的相關(guān)性经宏。
corGSM_MM
我們可以使用GeneIntegrationMatrix代替進(jìn)行相同的分析GeneScoreMatrix犀暑。
corGIM_MM <- correlateMatrices(
ArchRProj = projHeme5,
useMatrix1 = "GeneIntegrationMatrix",
useMatrix2 = "MotifMatrix",
reducedDims = "IterativeLSI"
)
corGIM_MM
步驟3.將最大Delta偏差添加到相關(guān)數(shù)據(jù)幀
對(duì)于這些相關(guān)分析中的每一個(gè)熄捍，我們都可以使用在步驟1中計(jì)算出的聚類之間觀察到的最大增量來(lái)注釋每個(gè)基序。
corGSM_MM$maxDelta <- rowData(seZ)[match(corGSM_MM$MotifMatrix_name, rowData(seZ)$name), "maxDelta"]
corGIM_MM$maxDelta <- rowData(seZ)[match(corGIM_MM$MotifMatrix_name, rowData(seZ)$name), "maxDelta"]
步驟4.確定正TF調(diào)節(jié)器
我們可以使用所有這些信息來(lái)確定正TF調(diào)節(jié)器母怜。在以下示例中余耽，我們將正調(diào)節(jié)劑視為其基序與基因得分（或基因表達(dá)）之間的相關(guān)性大于0.5，且調(diào)整后的p值小于0.01苹熏，且z值偏差的最大簇間差異為正的那些TF在前四分之一碟贾。
我們應(yīng)用這些選擇標(biāo)準(zhǔn)，并進(jìn)行一些文字處理以隔離TF名稱轨域。
corGSM_MM <- corGSM_MM[order(abs(corGSM_MM$cor), decreasing = TRUE), ]
corGSM_MM <- corGSM_MM[which(!duplicated(gsub("\\-.*","",corGSM_MM[,"MotifMatrix_name"]))), ]
corGSM_MM$TFRegulator <- "NO"
corGSM_MM$TFRegulator[which(corGSM_MM$cor > 0.5 & corGSM_MM$padj < 0.01 & corGSM_MM$maxDelta > quantile(corGSM_MM$maxDelta, 0.75))] <- "YES"
sort(corGSM_MM[corGSM_MM$TFRegulator=="YES",1])
從基因得分和基序偏差z分?jǐn)?shù)中識(shí)別出這些陽(yáng)性TF調(diào)節(jié)劑后袱耽，我們可以在點(diǎn)圖中突出顯示它們。
p <- ggplot(data.frame(corGSM_MM), aes(cor, maxDelta, color = TFRegulator)) +
geom_point() +
theme_ArchR() +
geom_vline(xintercept = 0, lty = "dashed") +
scale_color_manual(values = c("NO"="darkgrey", "YES"="firebrick3")) +
xlab("Correlation To Gene Score") +
ylab("Max TF Motif Delta") +
scale_y_continuous(
expand = c(0,0),
limits = c(0, max(corGSM_MM$maxDelta)*1.05)
)
p
我們可以對(duì)從得出的相關(guān)性執(zhí)行相同的分析GeneIntegrationMatrix干发。
corGIM_MM <- corGIM_MM[order(abs(corGIM_MM$cor), decreasing = TRUE), ]
corGIM_MM <- corGIM_MM[which(!duplicated(gsub("\\-.*","",corGIM_MM[,"MotifMatrix_name"]))), ]
corGIM_MM$TFRegulator <- "NO"
corGIM_MM$TFRegulator[which(corGIM_MM$cor > 0.5 & corGIM_MM$padj < 0.01 & corGIM_MM$maxDelta > quantile(corGIM_MM$maxDelta, 0.75))] <- "YES"
sort(corGIM_MM[corGIM_MM$TFRegulator=="YES",1])
p <- ggplot(data.frame(corGIM_MM), aes(cor, maxDelta, color = TFRegulator)) +
geom_point() +
theme_ArchR() +
geom_vline(xintercept = 0, lty = "dashed") +
scale_color_manual(values = c("NO"="darkgrey", "YES"="firebrick3")) +
xlab("Correlation To Gene Expression") +
ylab("Max TF Motif Delta") +
scale_y_continuous(
expand = c(0,0),
limits = c(0, max(corGIM_MM$maxDelta)*1.05)
)
p
