ArchR官网教程学习笔记12:Motif和Feature富集

系列回顾：
ArchR官网教程学习笔记1：Getting Started with ArchR
ArchR官网教程学习笔记2:基于ArchR推测Doublet
ArchR官网教程学习笔记3:创建ArchRProject
ArchR官网教程学习笔记4:ArchR的降维
ArchR官网教程学习笔记5:ArchR的聚类
ArchR官网教程学习笔记6:单细胞嵌入（Single-cell Embeddings）
ArchR官网教程学习笔记7:ArchR的基因评分和Marker基因
ArchR官网教程学习笔记8:定义与scRNA-seq一致的聚类
ArchR官网教程学习笔记9:ArchR的伪批量重复
ArchR官网教程学习笔记10:ArchR的call peak
ArchR官网教程学习笔记11:鉴定Marker峰

在确定了一个峰集之后，我们通常想要预测哪些转录因子可能会介导产生这些可接近的染色质位点的结合事件。这有助于评估marker峰或差异峰，以了解这些峰是否富集于特定转录因子的结合位点。例如，我们经常发现，在细胞类型特异的可接近染色质区域中，定义关键谱系的TFs得到了富集。以类似的方式，我们希望测试不同的峰组，以富集其他已知特征。例如，我们想知道A细胞的细胞类型特异的ATAC-seq峰是否富集于另一组基因组区域，如ChIP-seq峰。本章详细说明了如何在ArchR中执行这些富集。

（一）差异峰的motif富集

继续我们对前一章的差异峰的分析，我们可以在各种细胞类型中寻找在上调的峰上或下调的峰富集的motif。要做到这一点，我们必须首先将这些motif注解添加到我们的ArchRProject中。创建一个二进制矩阵，其中存在于在每个峰的motif是用数字表示的。我们使用addMotifAnnotations()函数来实现这一点，该函数决定了motif是否在ArchRProject中存储的峰集中存在。

> devtools::install_github("GreenleafLab/chromVARmotifs")
> projHeme5 <- addMotifAnnotations(ArchRProj = projHeme5, motifSet = "cisbp", name = "Motif")

然后，我们可以使用在前一章中生成的SummarizedExperiment对象markerTest的差异测试来定义我们感兴趣的显著差异峰的集合测试motif富集。在本例中，我们寻找FDR <= 0.1和Log2FC >= 0.5的峰。在markerTest的差异比较中，这些峰在“Erythroid”细胞中比在“Progenitor”中更容易接近。我们可以使用peakAnnoEnrichment()函数来测试这些不同的可接近的峰，以富集各种motif。这个函数可以用于许多不同的富集测试，我们将在本章中演示。

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

peakAnnoEnrichment()这个函数的输出是一个SummarizedExperiment对象，包含了多个assays，储存了超几何检验富集测试的结果：

> 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):

为了准备使用ggplot绘制这个数据，我们可以创建一个简化的data.frame对象，其中包含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))

正如预期的那样，在“Erythroid”细胞中更容易接近的峰中最丰富的motif对应于GATA转录因子，这与已被充分研究的GATA1在erythroid分化中的作用一致：

> head(df)
           TF mlog10Padj rank
383 GATA1_383   593.4009    1
388 GATA2_388   575.7300    2
385 GATA5_385   451.0765    3
384 GATA3_384   449.8240    4
386 GATA4_386   341.9064    5
387 GATA6_387   232.8139    6

使用ggplot，我们可以绘制排序的TF motifs，并根据它们的富集的显著性为它们着色。这里我们使用ggrepel来标记每个TF motif:

> ggUp <- ggplot(df, aes(rank, mlog10Padj, color = mlog10Padj)) + 
  geom_point(size = 1.5) +
  ggrepel::geom_label_repel(
    data = df[rev(seq_len(30)), ], aes(x = rank, y = mlog10Padj, label = TF), 
    size = 2.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

我们还可以画出Progenitor细胞中更易接近的峰的motif（Log2FC <= -0.5）：

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

> motifsDo
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):

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

> head(df)
           TF mlog10Padj rank
326  ELF2_326   88.94726    1
56   TCF12_56   83.25016    2
733 RUNX1_733   70.43406    3
843 ASCL1_843   61.00876    4
801  CBFB_801   59.59376    5
336  SPIB_336   54.85216    6

在这个例子里，在Progenitor细胞中更易接近的峰富集的motif有RUNX, ELF,和CBFB。

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

（二）Marker峰的motif富集

与上一节中对差异峰进行的motif富集分析类似，我们也可以对marker峰使用getMarkerFeatures()进行motif富集。

为此，我们将marker peak SummarizedExperiment (markersPeaks)传递给peakAnnotationEnrichment()函数。

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

peakAnnoEnrichment()的输出是SummarizedExperiment对象，包含了多个assays，储存了超几何检验的富集测试结果。

> 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()函数直接绘制所有细胞群的motif富集。在这个函数中，我们使用参数n限制每个细胞群显示的motif的总数：

> heatmapEM <- plotEnrichHeatmap(enrichMotifs, n = 7, transpose = TRUE)
> ComplexHeatmap::draw(heatmapEM, heatmap_legend_side = "bot", annotation_legend_side = "bot")

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

（三）ArchR的富集

除了测试峰的motif富集，ArchR还可以定制更个性化的富集。为了促进这一层次的数据探索，我们已经挑选了几个不同的特征集，可以很容易地测试你感兴趣峰区域的富集。我们将在下面描述每一个经过挑选的特性集。

（１）Encode TF Binding Sites

ENCODE已经比对许多细胞类型和因子的TF结合位点。我们可以使用这些TFBS collections来更好的理解我们的clusters。例如，在一个完全不知道的细胞类型里，这些富集可以帮助我们鉴定细胞类型。为了能够用这些ENCODE TFBS分析特征集，我们使用addArchRAnnotations()函数，并使用collection = "EncodeTFBS"。在使用addPeakAnnotations()时，它会创建以二进制代表所有marker peaks和所有ENCODE TFBS的overlap。

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

我们可以用我们的峰集来测试ENCODE TFBSs的富集，使用peakAnnoEnrichment()函数：

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

> 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):

画热图：

> heatmapEncode <- plotEnrichHeatmap(enrichEncode, n = 7, transpose = TRUE)
> ComplexHeatmap::draw(heatmapEncode, heatmap_legend_side = "bot", annotation_legend_side = "bot")
#保存
> plotPDF(heatmapEncode, name = "EncodeTFBS-Enriched-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme5, addDOC = FALSE)

（２）Bulk ATAC-seq

与ENCODE TF结合位点类似，我们也挑选了可用于重叠富集测试的bulk ATAC-seq实验的peak calls。我们通过设置collection = "ATAC"来访问这些 bulk ATAC-seq峰集。

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

> enrichATAC <- peakAnnoEnrichment(
    seMarker = markersPeaks,
    ArchRProj = projHeme5,
    peakAnnotation = "ATAC",
    cutOff = "FDR <= 0.1 & Log2FC >= 0.5"
  )
> 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):

画热图：

> heatmapATAC <- plotEnrichHeatmap(enrichATAC, n = 7, transpose = TRUE)
> ComplexHeatmap::draw(heatmapATAC, heatmap_legend_side = "bot", annotation_legend_side = "bot")
> plotPDF(heatmapATAC, name = "ATAC-Enriched-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme5, addDOC = FALSE)

（３）Codex TFBS

我们也可以使用CODEX TFBSs来富集我们峰集的motif。

> projHeme5 <- addArchRAnnotations(ArchRProj = projHeme5, collection = "Codex")
> enrichCodex <- peakAnnoEnrichment(
    seMarker = markersPeaks,
    ArchRProj = projHeme5,
    peakAnnotation = "Codex",
    cutOff = "FDR <= 0.1 & Log2FC >= 0.5"
  )
> enrichCodex
class: SummarizedExperiment 
dim: 189 11 
metadata(0):
assays(10): mlog10Padj mlog10p ... CompareFrequency
  feature
rownames(189): 1.STAT5-No_drug_(DMSO)...
  2.RUNX3-GM12878_cell_fr... ...
  188.TP53-codex_Embryonic... 189.TP53-codex_Embryonic...
rowData names(0):
colnames(11): B CD4.M ... PreB Progenitor
colData names(0):

画热图：

> heatmapCodex <- plotEnrichHeatmap(enrichCodex, n = 7, transpose = TRUE)
> ComplexHeatmap::draw(heatmapCodex, heatmap_legend_side = "bot", annotation_legend_side = "bot")
> plotPDF(heatmapCodex, name = "Codex-Enriched-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme5, addDOC = FALSE)

（四）自定义富集

除了所有这些精心挑选的注释集，ArchR还能够接受用户定义的注释来执行自定义富集。在下面的示例中，我们将说明如何基于select ENCODE ChIP-seq实验创建自定义注释。

首先，我们将定义将要使用的数据集，并提供下载链接。本地文件也可以以同样的方式使用。

> 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"
)

接下来我们添加自定义注释到ArchRProject对象，使用addPeakAnnotation()函数。这里，我们把这个自定义注释叫做“ChIP”：

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

与前面一样，我们使用这个自定义注释通过peakAnnoEnrichment()执行peak注释富集，并按照相同的步骤创建我们的注释热图：

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

> enrichRegions
class: SummarizedExperiment 
dim: 4 11 
metadata(0):
assays(10): mlog10Padj mlog10p ... CompareFrequency
  feature
rownames(4): Encode_K562_GATA1 Encode_GM12878_CEBPB
  Encode_K562_Ebf1 Encode_K562_Pax5
rowData names(0):
colnames(11): B CD4.M ... PreB Progenitor
colData names(0):

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

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