TF enrichment of linked peaks

[waltno]

I followed the methods in

https://www.sciencedirect.com/science/article/pii/S0092867421009429?via%3Dihub

Which is a Greenleaf paper :) If I did something wrong feel free to correct me!

`
#get linked peaks

#I have three kmeans clusters
p <- plotPeak2GeneHeatmap(ArchRProj = m2, groupBy = "cell_line", k = 3, returnMatrices = T)
link_df<- data.frame(p@listData$Peak2GeneLinks)
names(p@listData$RNA@listData$kmeansId)<- rownames(p@listData[["RNA"]]@ListData[["matrix"]])
link_df$kmean <- p@listData$RNA@listData$kmeansId[match(names(p@listData$RNA@listData$kmeansId), rownames(link_df))]

link_df$peak <- gsub(":", "-", link_df$peak)
k1<- link_df[link_df$kmean == "1",]
k2<- link_df[link_df$kmean == "2",]
k3<- link_df[link_df$kmean == "3",]

peak<-getPeakSet(m2)
names(peak)<-NULL

list_peak<-GRangesList(StringToGRanges(k1$peak), StringToGRanges(k2$peak), StringToGRanges(k3$peak))

##I went here and manually downloaded a LOLA Core database http://big.databio.org/regiondb/
##not sure which one is the best core database, I picked the most recent one

library(LOLA)

regionDB = loadRegionDB("~/nm/t1/resources/regions/LOLACore/hg38")
locResults = runLOLA(userSets = list_peak, userUniverse = peak, regionDB, cores=1)
writeCombinedEnrichment(locResults, outFolder= "lolaResults", includeSplits=TRUE)
`

[rcorces]

@waltno - I'm not sure how what you are looking to accomplish is different from what is enabled by peakAnnoEnrichment() in ArchR. Can you look through https://www.archrproject.com/bookdown/archr-enrichment.html and clarify what you need?

Perhaps you're requesting how to do peak enrichment when you dont have a SummarizedExperiment object from ArchR?

[waltno]

Right. Im just trying to get motif enrichment for linked peaks. If I can make a summarized experiment of just my peaks of interest, with the necessary stats (Log2FC, AUC, etc) for peakAnnoEnrichment() that'd be great.

[rcorces]

I dont believe there is a plug and play way to do this but its a reasonable feature request. In the meantime, you should be able to hack this by adding the lola annotations to your whole peakset via addArchRAnnotations() and then creating a summarized experiment from scratch that mimics what is expected as input to peakAnnoEnrichment and making sure the values you set for FDR and Log2FC pass whatever you give to cutoff

[Zepeng-Mu]

I think this is very similar to this discussion: #1331

[waltno]

Yea! I dont see how the "de_res"part was made?

[Zepeng-Mu]

Oh, that's just some testing results, similar to what you get from marker peaks analysis in ArchR. You should replace with a matrix that contains statistics from your analysis.

[waltno]

Thank you all for your input. This is still not making total sense to me...

If return the peak-to-gene matrix, I get Correlation and FDR values for each link per kmeans cluster. Would that be sufficient to construct a Summarized Experiment for peakAnnoEnrichment()? It still would not have as many peaks as my ArchR project.

[Zepeng-Mu]

If I'm understanding your needs correctly, you have peaks that are linked to genes (enhancers) and they belong to different clusters from ArchR kmeans clustering strategy. You want to test the enrichment of TF binding in enhancers from each of the K clusters.

I think it would make sense to construct a matrix M for example, which has enhancers as rows and k-clusters as columns. The tricky thing is each cluster have different number of enhancers, so to construct a matrix you need to used the union of all enhancers from the K clusters. In this M matrix, I would set M_ij = 1 if enhancer i belongs to cluster j, otherwise set to 0.

Then construct a SummarizedExperiment using M as the assay. When doing enrichment analysis, use cutOff = 'M == 1'. By doing so, AcrhR evaluate M==1 on each element in the matrix using eval(parse(text="M==1")) under the hood, which is used to subset the matrix, effectively getting enhancers in each of K clusters from M and then do enrichment analysis only on those peaks. In this way, the output is also a matrix that can be nicely plotted using ArchR functions.

If you want to use FDR and correlation as cutoff, you can build a SummarizedExperiment by having to matrices similar to M. You can set FDR=1 and correlation=0 for peaks that do not belong to a certain K, for example. Then you can use cutOff = "FDR < 0.1 & corr > 0.8" in your enrichment call.

This is the best explanation that I can provide based on my understanding of your question. Let me know if I'm misunderstanding what you are trying to do.

[waltno]

Okay I think I got it! Here's a giant chunk of code

p <- plotPeak2GeneHeatmap(ArchRProj = m2, groupBy = "cell_line", k = 3, returnMatrices = T)

names(p@listData$RNA@listData$kmeansId)<- rownames(p@listData[["RNA"]]@listData[["matrix"]])

link_df$kmean <- p@listData$RNA@listData$kmeansId[match(names(p@listData$RNA@listData$kmeansId), rownames(link_df))]

link_df<- data.frame(p@listData$Peak2GeneLinks)

mat<-matrix(nrow = length(link_df$peak), ncol = 3)
rownames(mat)<-link_df$peak
colnames(mat)<-c("1", "2", "3")

for(x in colnames(mat)){
  mat[, x]<-rep(0, length(rownames(mat)))
  hit<- which(rownames(mat) %in% link_df[link_df$kmean == x,]$peak)
  mat[,x][hit]<-1
}

link_df$peak <- gsub(":", "-", link_df$peak)

library(Signac)
peak<-StringToGRanges(link_df$peak)

mySe <- SummarizedExperiment(
  assays = list(m = mat),
  rowData = as.data.frame(peak),
  metadata = list(Params = list(useMatrix = "PeakMatrix"))
)
enhancer_project<-addPeakSet(m2,  peakSet = peak,
  force = T)

rowData(mySe)<-sort(rowData(mySe))

enrichMotifs <- peakAnnoEnrichment(
  seMarker = mySe,
  ArchRProj = enhancer_project,
  peakAnnotation = "Motif",
  cutOff = "m == 1"
)
heatmapEM <- plotEnrichHeatmap(enrichMotifs, n = 10, transpose = TRUE,   returnMatrix = F)
plotPDF(heatmapEM, name = "Motifs-Linked-Heatmap", width = 5, height = 5, ArchRProj = m2, addDOC = FALSE)

I got this heatmap!!!

Motifs-Linked-Heatmap.pdf