本文首發(fā)于我的個人博客， http://xuzhougeng.top/

往期回顧:

第11章: 使用ArchR鑒定標(biāo)記Peak

在討論基因得分(gene score)這一章中翠勉，我們已經(jīng)介紹了標(biāo)記特征的鑒定政冻。相同的函數(shù)getMakerFeatures()也能夠用于從ArchRProject任意矩陣中鑒定標(biāo)記特征阐斜。所謂的標(biāo)記特征指的是相對于其他細(xì)胞分組唯一的特征。這些特征能幫助我們理解類群或者細(xì)胞類型特異的生物學(xué)現(xiàn)象。在這一章中阐滩，我們會討論如何使用該功能鑒定標(biāo)記peak虑省。

11.1: 使用ArchR鑒定標(biāo)記Peak

通常而言匿刮，我們想知道哪些peak是某個聚類或者某一些聚類所特有的。在ArchR中探颈，這可以通過設(shè)置addMarkFeatures()函數(shù)的useMatrix="PeakMatrix"來實(shí)現(xiàn)（無需監(jiān)督）熟丸。

首先，我們需要再看一眼projHeme5中有哪些細(xì)胞類型伪节，以及它們的相對比例

table(projHeme5$Cluster2)

現(xiàn)在光羞，讓我們調(diào)用getMarkerFeatures參數(shù)，并設(shè)置useMatrix="PeakMatrix". 此外怀大，為了降低不同細(xì)胞組之間的數(shù)據(jù)質(zhì)量對結(jié)果的影響纱兑，我們可以設(shè)置bias參數(shù)，其中bias = c("TSSEnrichment", "log10(nFrags)")就是用來避免TSS富集和每個細(xì)胞的fragment數(shù)對結(jié)果的影響化借。

markersPeaks <- getMarkerFeatures(
    ArchRProj = projHeme5, 
    useMatrix = "PeakMatrix", 
    groupBy = "Clusters2",
  bias = c("TSSEnrichment", "log10(nFrags)"),
  testMethod = "wilcoxon"
)

getMarkerFeatures()函數(shù)返回一個SummarizedExperiment對象潜慎，該對象包含一些不同的assays

markerPeaks

接著，我們可以用getMarkers函數(shù)從輸出的SummarizedExperiment對象中提取我們感興趣的部分。默認(rèn)情況下铐炫，它會返回一個包含多個DataFrame的列表垒手，不同的DataFrame表示來自不同的細(xì)胞分組。

markerList <- getMarkers(markersPeaks, cutOff = "FDR <= 0.01 & Log2FC >= 1")
markerList

如果我們對特定的一個細(xì)胞分組感興趣倒信，我們可以用$進(jìn)行提取科贬。

markerList$Erythroid

除了返回一個包含多個DataFrame的列表外，我們還可以用getMarkers()返回一個GRangesList鳖悠，只要設(shè)置returnGR=TRUE即可榜掌。

markerList <- getMarkers(markersPeaks, cutOff = "FDR <= 0.01 & Log2FC >= 1", returnGR = TRUE)
markerList

這個GRangesList同樣可以用$提取特定細(xì)胞組的結(jié)果，返回的是一個GRanges對象

markerList$Erythroid

11.2 在ArchR中繪制Marker Peaks

ArchR提供了許多繪圖函數(shù)用于getMarkerfeatuers()返回的SummarizedExperiment對象的可視化竞穷。

11.2.1 Marker Peak Heatmap

markerHeatmap能以熱圖的形式展示標(biāo)記Marker Peak（或者其他getMarkerFeatures()輸出的特征）

heatmapPeaks <- markerHeatmap(
  seMarker = markersPeaks, 
  cutOff = "FDR <= 0.1 & Log2FC >= 0.5",
  transpose = TRUE
)

使用draw函數(shù)繪制結(jié)果

plotPDF(heatmapPeaks, name = "Peak-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme5, addDOC = FALSE)

Marker Peak Heatmap

plotFDF()函數(shù)能夠以可編輯的矢量版本保存圖片

plotPDF(heatmapPeaks, name = "Peak-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme5, addDOC = FALSE)

11.2.2 Marker Peak MA和火山圖

除了繪制熱圖唐责，我們也可以為每個細(xì)胞分組繪制MA或者火山圖(volcano)。這些圖可以用markerPlot()函數(shù)繪制瘾带。對于MA圖鼠哥，需要設(shè)置參數(shù)plotAs="MA". 我們以"Erythroid"細(xì)胞分組為例，設(shè)置參數(shù)name = "Erythroid"

pma <- markerPlot(seMarker = markersPeaks, name = "Erythroid", cutOff = "FDR <= 0.1 & Log2FC >= 1", plotAs = "MA")
pma

iMA圖

同樣的看政，只要設(shè)置plotAs="Volcano"就可以繪制火山圖

pv <- markerPlot(seMarker = markersPeaks, name = "Erythroid", cutOff = "FDR <= 0.1 & Log2FC >= 1", plotAs = "Volcano")
pv

volcano plot

plotFDF()函數(shù)能夠以可編輯的矢量版本保存圖片朴恳。

plotPDF(pma, pv, name = "Erythroid-Markers-MA-Volcano", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)

11.2.3 Browser Tracks的Marker Peak

此外，我們在基因組瀏覽器上檢查這些peak區(qū)域允蚣，只需要為plotBrowserTrack()函數(shù)的features參數(shù)傳入相應(yīng)的peak區(qū)間于颖。這會額外在我們的ArchR track圖的下方以BED形式展示marker peak區(qū)域。

p <- plotBrowserTrack(
    ArchRProj = projHeme5, 
    groupBy = "Clusters2", 
    geneSymbol = c("GATA1"),
    features =  getMarkers(markersPeaks, cutOff = "FDR <= 0.1 & Log2FC >= 1", returnGR = TRUE)["Erythroid"],
    upstream = 50000,
    downstream = 50000
)

我們使用grid::grid.draw()繪制結(jié)果

grid::grid.newpage()
grid::grid.draw(p$GATA1)

browser track

plotFDF()函數(shù)能夠以可編輯的矢量版本保存圖片嚷兔。

plotPDF(p, name = "Plot-Tracks-With-Features", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)

11.3 組間配對檢驗(yàn)

標(biāo)記特征鑒定是一種特別的差異表達(dá)檢驗(yàn)森渐。此外，使用相同的getMarkerFeatures()函數(shù)也能實(shí)現(xiàn)標(biāo)準(zhǔn)化的差異分析冒晰。我們只需要設(shè)置useGroup為一組細(xì)胞同衣，然后設(shè)置bgdGroup為另一組細(xì)胞即可。這就可以對給定兩組進(jìn)行差異分析壶运。在這些差異分析中耐齐，在useGroups比較高的peak的倍數(shù)變化值是正數(shù)，在bgdGroups比較高的peak則是由負(fù)的倍數(shù)變化值蒋情。

這里埠况，我們對"Erythroid"與"Progenitor"細(xì)胞組進(jìn)行配對檢驗(yàn)。

markerTest <- getMarkerFeatures(
  ArchRProj = projHeme5, 
  useMatrix = "PeakMatrix",
  groupBy = "Clusters2",
  testMethod = "wilcoxon",
  bias = c("TSSEnrichment", "log10(nFrags)"),
  useGroups = "Erythroid",
  bgdGroups = "Progenitor"
)

使用markerPlot()函數(shù)可以繪制MA或者火山圖棵癣。MA圖需要設(shè)置plotAs="MA"

pma <- markerPlot(seMarker = markerTest, name = "Erythroid", cutOff = "FDR <= 0.1 & abs(Log2FC) >= 1", plotAs = "MA")
pma

group marker peak MA plot

火山圖需要設(shè)置plotAs="Volcano"

pv <- markerPlot(seMarker = markerTest, name = "Erythroid", cutOff = "FDR <= 0.1 & abs(Log2FC) >= 1", plotAs = "Volcano")
pv

group marker peak volcano

plotFDF()函數(shù)能夠以可編輯的矢量版本保存圖片辕翰。

plotPDF(pma, pv, name = "Erythroid-vs-Progenitor-Markers-MA-Volcano", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)

在后續(xù)章節(jié)中，我們還會進(jìn)行差異分析狈谊，因?yàn)闀谖覀兊牟町愰_放的peak中搜索富集的motif金蜀。

使用ArchR分析單細(xì)胞ATAC-seq數(shù)據(jù)(第十一章)