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

往期回顧:

第12章使用ArchR進行motif和特征富集分析

在鑒定到可靠的peak集之后，我們也會想預(yù)測有哪些轉(zhuǎn)錄因子參與了結(jié)合事件(binding events)祝钢，從而產(chǎn)生了這些染色質(zhì)開放位點。在分析標記peak或差異peak時，這能幫助我們更好的理解為什么某組的peak會富集某一類轉(zhuǎn)錄因子的結(jié)合位點投蝉。舉個例子，我們想在細胞特異的染色質(zhì)開放區(qū)域中找到定義譜系的關(guān)鍵轉(zhuǎn)錄因子征堪。同樣瘩缆，我們也想根據(jù)其他已知特征對不同的peak進行富集分析。比如說请契，我們像知道是否細胞類型A的細胞特異性ATAC-seq peak對于另一組基因組區(qū)域（如ChIP-seq peak）也富集咳榜。這一章會詳細介紹ArchR中的富集分析原理。

12.1 差異peak中的motif富集

繼續(xù)上一章的差異peak分析爽锥，我們可以尋找在不同類型細胞富集的peak中的motif涌韩。我們需要先將motif的注釋信息加入到我們的ArchRProject中。我們調(diào)用addMotifAnnotations()函數(shù)分析ArchRProject的peak中是否存在motif氯夷。運行結(jié)束后會在ArchRProject對象中加入一個新的二值矩陣臣樱，用于判斷peak是否包括motif。

projHeme5 <- addMotifAnnotations(ArchRProj = projHeme5, motifSet = "cisbp", name = "Motif")

接著我們使用上一章差異檢驗得到的markerTest分析motif的富集情況腮考，這是一個SummarizedExperiment對象雇毫。我們用peakAnnoEnrichments()函數(shù)分析這些差異開放peak是否富集某一類moitf〔任担可以設(shè)置cutOff來過濾peak棚放，例如PDR <= 0.1 & Log2FC >=0.5記錄的是"Erythroid"比"Progenitor"更開放的peak。

注: peakAnnoEnrichment()能用于多種差異富集檢驗馅闽，在后續(xù)章節(jié)還會介紹飘蚯。

motifsUp <- peakAnnoEnrichment(
    seMarker = markerTest,
    ArchRProj = projHeme5,
    peakAnnotation = "Motif",
    cutOff = "FDR <= 0.1 & Log2FC >= 0.5"
  )

輸出的peakAnnoEnrichment()是一個SummarizedExperiment對象，里面存放著多個assays, 記錄著超幾何檢驗的富集結(jié)果福也。

motifsUp
# class: SummarizedExperiment
# dim: 870 1
# metadata(0):
# assays(10): mlog10Padj mlog10p … CompareFrequency feature
# rownames(870): TFAP2B_1 TFAP2D_2 … TBX18_869 TBX22_870
# rowData names(0):
# colnames(1): Erythroid
# colData names(0):

然后局骤，我們創(chuàng)建一個data.frame對象用于ggplot作圖，包括motif名暴凑，矯正的p值和顯著性排序峦甩。

df <- data.frame(TF = rownames(motifsUp), mlog10Padj = assay(motifsUp)[,1])
df <- df[order(df$mlog10Padj, decreasing = TRUE),]
df$rank <- seq_len(nrow(df))

正如我們所預(yù)期的那樣，"Erythroid"里開放的peak富集的motif主要是GATA轉(zhuǎn)錄因子现喳，符合以往研究中"GATA1"在erythroid分化中發(fā)揮的作用凯傲。

head(df)

使用ggplot展示結(jié)果犬辰，以ggrepel來標識每個TF motif名。

ggUp <- ggplot(df, aes(rank, mlog10Padj, color = mlog10Padj)) + 
  geom_point(size = 1) +
  ggrepel::geom_label_repel(
        data = df[rev(seq_len(30)), ], aes(x = rank, y = mlog10Padj, label = TF), 
        size = 1.5,
        nudge_x = 2,
        color = "black"
  ) + theme_ArchR() + 
  ylab("-log10(P-adj) Motif Enrichment") + 
  xlab("Rank Sorted TFs Enriched") +
  scale_color_gradientn(colors = paletteContinuous(set = "comet"))

ggUp

Erythroid-vs-Progenitor-Markers-Motifs-Enriched_1

通過設(shè)置Log2FC <= 0.5我們可以挑選出在"Progenitor"里更加開放的peak泣洞，然后分析其中富集的motif忧风。

motifsDo <- peakAnnoEnrichment(
    seMarker = markerTest,
    ArchRProj = projHeme5,
    peakAnnotation = "Motif",
    cutOff = "FDR <= 0.1 & Log2FC <= -0.5"
  )
motifsDo

準備繪圖所需數(shù)據(jù)框

df <- data.frame(TF = rownames(motifsDo), mlog10Padj = assay(motifsDo)[,1])
df <- df[order(df$mlog10Padj, decreasing = TRUE),]
df$rank <- seq_len(nrow(df))

此時，我們會發(fā)現(xiàn)在"Progenitor"細胞更加開放的peak中球凰，更多富集RUNX, ELF和CBFB狮腿。

head(df)
# TF mlog10Padj rank
# 326 ELF2_326 88.68056 1
# 733 RUNX1_733 64.00586 2
# 801 CBFB_801 53.55426 3
# 732 RUNX2_732 53.14766 4
# 734 ENSG00000250096_734 53.14766 5
# 336 SPIB_336 52.79666 6

使用ggplot展示結(jié)果。

ggDo <- ggplot(df, aes(rank, mlog10Padj, color = mlog10Padj)) + 
  geom_point(size = 1) +
  ggrepel::geom_label_repel(
        data = df[rev(seq_len(30)), ], aes(x = rank, y = mlog10Padj, label = TF), 
        size = 1.5,
        nudge_x = 2,
        color = "black"
  ) + theme_ArchR() + 
  ylab("-log10(FDR) Motif Enrichment") +
  xlab("Rank Sorted TFs Enriched") +
  scale_color_gradientn(colors = paletteContinuous(set = "comet"))

ggDo

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

plotPDF(ggUp, ggDo, name = "Erythroid-vs-Progenitor-Markers-Motifs-Enriched", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)

12.2 標記Peak的motif富集分析

和之前利用差異peak的motif富集分析類似缘厢，我們同樣能用getMarkerFeatures()分析標記peak里富集的motif。

我們向函數(shù)peakAnnotationEnrichment()傳入存放標記peak的SummarizedExperiment對象甩挫，即markersPeaks

enrichMotifs <- peakAnnoEnrichment(
    seMarker = markersPeaks,
    ArchRProj = projHeme5,
    peakAnnotation = "Motif",
    cutOff = "FDR <= 0.1 & Log2FC >= 0.5"
  )

輸出的peakAnnoEnrichment()是一個SummarizedExperiment對象贴硫，里面存放著多個assays, 記錄著超幾何檢驗的富集結(jié)果。

enrichMotifs
# class: SummarizedExperiment
# dim: 870 11
# metadata(0):
# assays(10): mlog10Padj mlog10p … CompareFrequency feature
# rownames(870): TFAP2B_1 TFAP2D_2 … TBX18_869 TBX22_870
# rowData names(0):
# colnames(11): B CD4.M … PreB Progenitor
# colData names(0):

直接用plotEnrichHeatmap()函數(shù)繪制不同細胞組的富集的motif伊者。通過設(shè)置參數(shù)n限制每個細胞分組中展示的motif英遭。

heatmapEM <- plotEnrichHeatmap(enrichMotifs, n = 7, transpose = TRUE)

使用ComplexHeatmap::draw()函數(shù)展示結(jié)果

ComplexHeatmap::draw(heatmapEM, heatmap_legend_side = "bot", annotation_legend_side = "bot")

Motifs-Enriched-Marker-Heatmap

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

plotPDF(heatmapEM, name = "Motifs-Enriched-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme5, addDOC = FALSE)

12.3 ArchR富集分析

除了分析peak中富集motif, ArchR還能進行個性化的富集分析亦渗。為了方便這類數(shù)據(jù)探索挖诸，我們已人工確定了一些特征數(shù)據(jù)集，它們能比較容易地在我們感興趣的peak區(qū)間進行檢驗法精。我們接下來將會逐個介紹這些特征數(shù)據(jù)集多律。該分析最初受LOLA啟發(fā)。

12.3.1 Encode TF 結(jié)合位點

ENCODE協(xié)會已經(jīng)將TF結(jié)合位點(TFBS)匹配到多種細胞類型和因子中搂蜓。我們可以利用這些TFBS去更好地理解聚類結(jié)果狼荞。例如，我們可以根據(jù)富集結(jié)果去判斷未知細胞類型的可能類型帮碰。為了能夠使用ENCODE TFBS特征集進行分析相味，我們需要調(diào)用addArchRAnnotations()函數(shù)，設(shè)置collection = "EncodeTFBS". 和使用addPeakAnnotations()類似殉挽，這會創(chuàng)建一個二值矩陣丰涉，記錄我們的標記peak是否和ENCODE TFBS有重疊。

projHeme5 <- addArchRAnnotations(ArchRProj = projHeme5, collection = "EncodeTFBS")

我們接著使用peakAnnoEnrichment()函數(shù)分析這些 ENCODE TFBS是否在我們的peak中富集此再。

enrichEncode <- peakAnnoEnrichment(
    seMarker = markersPeaks,
    ArchRProj = projHeme5,
    peakAnnotation = "EncodeTFBS",
    cutOff = "FDR <= 0.1 & Log2FC >= 0.5"
  )

和之前一樣昔搂，該函數(shù)返回一個SummarizedExperiment對象玲销。

enrichEncode
# class: SummarizedExperiment
# dim: 689 11
# metadata(0):
# assays(10): mlog10Padj mlog10p … CompareFrequency feature
# rownames(689): 1.CTCF-Dnd41… 2.EZH2_39-Dnd41… …
# 688.CTCF-WERI_Rb_1… 689.CTCF-WI_38…
# rowData names(0):
# colnames(11): B CD4.M … PreB Progenitor
# colData names(0):

我們可以使用plotEnrichHeatmap函數(shù)從富集結(jié)果中創(chuàng)建熱圖输拇。

heatmapEncode <- plotEnrichHeatmap(enrichEncode, n = 7, transpose = TRUE)

然后用ComplexHeatmap::draw()繪制熱圖

ComplexHeatmap::draw(heatmapEncode, heatmap_legend_side = "bot", annotation_legend_side = "bot")

EncodeTFBS-Enriched-Marker-Heatmap

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

plotPDF(heatmapEncode, name = "EncodeTFBS-Enriched-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme5, addDOC = FALSE)

12.3.2 混池ATAC-seq

和ENCODE TFBS類似贤斜，我們還可以使用混池ATAC-seq實驗鑒定的peak策吠，分析兩者的重疊情況逛裤。通過設(shè)置collection="ATAC"來調(diào)用混池ATAC-seqpeak數(shù)據(jù)集。

projHeme5 <- addArchRAnnotations(ArchRProj = projHeme5, collection = "ATAC")

接著通過設(shè)置peakAnnotation = "ATAC"檢驗我們的標記peak是否富集了混池ATAC-seq的peak猴抹。

enrichATAC <- peakAnnoEnrichment(
    seMarker = markersPeaks,
    ArchRProj = projHeme5,
    peakAnnotation = "ATAC",
    cutOff = "FDR <= 0.1 & Log2FC >= 0.5"
  )

和之前一樣带族，該函數(shù)會輸出SummarizedExperiment對象，記錄著富集結(jié)果

enrichATAC
# class: SummarizedExperiment
# dim: 96 11
# metadata(0):
# assays(10): mlog10Padj mlog10p … CompareFrequency feature
# rownames(96): Brain_Astrocytes Brain_Excitatory_neurons … Heme_MPP
# Heme_NK
# rowData names(0):
# colnames(11): B CD4.M … PreB Progenitor
# colData names(0):

我們用plotEnrichHeatmap()函數(shù)基于SummarizedExperiment繪制富集熱圖

heatmapATAC <- plotEnrichHeatmap(enrichATAC, n = 7, transpose = TRUE)

使用ComplexHeatmap::draw()繪制結(jié)果

ComplexHeatmap::draw(heatmapATAC, heatmap_legend_side = "bot", annotation_legend_side = "bot")

ATAC-Enriched-Marker-Heatmap

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

plotPDF(heatmapATAC, name = "ATAC-Enriched-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme5, addDOC = FALSE)

12.3.3 Codex TFBS

相同類型的分析還能用于 CODEX TFBS蝙砌，只要設(shè)置collection = "Codex"即可

projHeme5 <- addArchRAnnotations(ArchRProj = projHeme5, collection = "Codex")
enrichCodex <- peakAnnoEnrichment(
    seMarker = markersPeaks,
    ArchRProj = projHeme5,
    peakAnnotation = "Codex",
    cutOff = "FDR <= 0.1 & Log2FC >= 0.5"
  )
heatmapCodex <- plotEnrichHeatmap(enrichCodex, n = 7, transpose = TRUE)
ComplexHeatmap::draw(heatmapCodex, heatmap_legend_side = "bot", annotation_legend_side = "bot")

Codex-Enriched-Marker-Heatmap

于是我們就保存圖片了

plotPDF(heatmapCodex, name = "Codex-Enriched-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme5, addDOC = FALSE)

12.4 自定義富集

除了之前這些經(jīng)過人工審核的注釋數(shù)據(jù)集，ArchR還能處理用戶自定義注釋信息來執(zhí)行富集分析跋理。接下來择克，我們會介紹如何根據(jù)ENCODE ChIP-seq實驗來創(chuàng)建自定義的注釋信息。

首先前普，我們先提供后續(xù)將被使用并下載的數(shù)據(jù)集肚邢，也可以提供本地文件。

EncodePeaks <- c(
  Encode_K562_GATA1 = "https://www.encodeproject.org/files/ENCFF632NQI/@@download/ENCFF632NQI.bed.gz",
  Encode_GM12878_CEBPB = "https://www.encodeproject.org/files/ENCFF761MGJ/@@download/ENCFF761MGJ.bed.gz",
  Encode_K562_Ebf1 = "https://www.encodeproject.org/files/ENCFF868VSY/@@download/ENCFF868VSY.bed.gz",
  Encode_K562_Pax5 = "https://www.encodeproject.org/files/ENCFF339KUO/@@download/ENCFF339KUO.bed.gz"
)

然后拭卿，我們用addPeakAnnotation()函數(shù)在ArchRProject函數(shù)中增加自定義注釋骡湖。我們這里將其命名為"ChIP"

projHeme5 <- addPeakAnnotations(ArchRProj = projHeme5, regions = EncodePeaks, name = "ChIP")

和之前一樣，我們使用peakAnnoEnrichment()函數(shù)根據(jù)自定義的注釋信息執(zhí)行peak注釋富集分析

enrichRegions <- peakAnnoEnrichment(
    seMarker = markersPeaks,
    ArchRProj = projHeme5,
    peakAnnotation = "ChIP",
    cutOff = "FDR <= 0.1 & Log2FC >= 0.5"
  )

并以相同的步驟生成注釋熱圖峻厚。

heatmapRegions <- plotEnrichHeatmap(enrichRegions, n = 7, transpose = TRUE)
ComplexHeatmap::draw(heatmapRegions, heatmap_legend_side = "bot", annotation_legend_side = "bot")

Regions-Enriched-Marker-Heatmap

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

plotPDF(heatmapRegions, name = "Regions-Enriched-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme5, addDOC = FALSE)

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

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

第12章使用ArchR進行motif和特征富集分析

12.1 差異peak中的motif富集

12.2 標記Peak的motif富集分析

12.3 ArchR富集分析

12.3.1 Encode TF 結(jié)合位點

12.3.2 混池ATAC-seq

12.3.3 Codex TFBS

12.4 自定義富集

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

第12章 使用ArchR進行motif和特征富集分析

12.1 差異peak中的motif富集

12.2 標記Peak的motif富集分析

12.3 ArchR富集分析

12.3.1 Encode TF 結(jié)合位點

12.3.2 混池ATAC-seq

12.3.3 Codex TFBS

12.4 自定義富集

第12章使用ArchR進行motif和特征富集分析