11a. Add motif annotations

proj5 <- addMotifAnnotations(ArchRProj = proj5,
                              motifSet = "cisbp",
                              name = "Motif_LSI_Harmony")

11b. Motif Enrichment in Marker Peaks

Here, we're looking at our specific motifs of interest.

enrichMotifs_LSI_Harmony <- peakAnnoEnrichment(
    seMarker = markersPeaks,
    ArchRProj = proj5,
    peakAnnotation = "Motif_LSI_Harmony",
    cutOff = "FDR <= 0.1 & Log2FC >= 0.25"
  )
enrichMotifs_LSI_Harmony
heatmapEM <- plotEnrichHeatmap(enrichMotifs_LSI_Harmony,
                               n = 5, transpose = TRUE, rastr=FALSE)
svglite(paste0(sample_name,
    "_scATAC-seq_Manual-Annotation_enrichMotif_heatmap_",
    Sys.Date(), ".svg"))
ComplexHeatmap::draw(heatmapEM, heatmap_legend_side = "bot",
                     annotation_legend_side = "bot")
dev.off()

# add a set of background peaks
proj5 <- addBgdPeaks(proj5)
proj5 <- addDeviationsMatrix(
  ArchRProj = proj5, 
  peakAnnotation = "Motif_LSI_Harmony",
  force = TRUE
)
plotVarDev <- getVarDeviations(proj5,
                               name = "Motif_LSI_HarmonyMatrix",
                               plot = TRUE)
svglite(paste0(sample_name,
    "_scATAC-seq_Manual-Annotation_chromVar_deviations_",
    Sys.Date(), ".svg"))
plotVarDev
dev.off()

motifs  <- c("FOXD1", "NOTCH", "RBPJ", "SP1", "CITED1", "SIX2", "SMTN", "SP1",
             "SMAD1", "E2F1", "HOXB7", "NFI", "MEF2A", "MEF2B", "MEF2C", 
             "MEF2D", "SMARCC1", "BACH1", "BACH2", "TCF21", "SNAI2",
             "NFIC", "NFIX", "JUN", "JUNB", "JUND")
markerMotifs <- getFeatures(proj5, select = paste(motifs, collapse="|"),
                            useMatrix = "Motif_LSI_HarmonyMatrix")
markerMotifs

# Get just z scores
markerMotifs <- grep("z:", markerMotifs, value = TRUE)
markerMotifs <- markerMotifs[markerMotifs %ni% "z:Mesp1_58"]
markerMotifs <- markerMotifs[markerMotifs %ni% "z:Sp100_714"]
markerMotifs <- markerMotifs[markerMotifs %ni% "z:Sp140_717"]
markerMotifs <- markerMotifs[markerMotifs %ni% "z:Sp110_718"]
markerMotifs <- markerMotifs[markerMotifs %ni% "z:Nfil3_138"]
markerMotifs <- markerMotifs[markerMotifs %ni% "z:Nfil3_131"]
markerMotifs <- markerMotifs[markerMotifs %ni% "z:Nfil3_137"]
markerMotifs <- markerMotifs[markerMotifs %ni% "z:Nfia_862"]
markerMotifs <- markerMotifs[markerMotifs %ni% "z:Nfib_859"]
markerMotifs

p <- plotGroups(ArchRProj = proj5, 
  groupBy = "ManualClusters", 
  colorBy = "Motif_LSI_HarmonyMatrix", 
  name = markerMotifs,
  imputeWeights = getImputeWeights(proj5)
)

p2 <- lapply(seq_along(p), function(x){
  if(x != 1){
    p[[x]] + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6) +
    theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm")) +
    theme(
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank(),
        axis.title.y=element_blank()
    ) + ylab("")
  }else{
    p[[x]] + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6) +
    theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm")) +
    theme(
        axis.ticks.y=element_blank(),
        axis.title.y=element_blank()
    ) + ylab("")
  }
})

# This plot can show us how FOXD1, for example, changes in accessibility across time points by cluster
svglite(paste0(sample_name,
               "_scATAC-seq_Manual-Annotation_markerMotifs_",
               Sys.Date(), ".svg"),
        height=8, width=24)
do.call(cowplot::plot_grid,
    c(list(nrow = 1, rel_widths = c(2, rep(1, length(p2) - 1))),p2))
dev.off()

11c. Motif enrichment of reported Ren1 enhancer binding TFs

enhancerMotifs <- c("NFIX", "SP1", "SP3", "CREB", "CREM", "USF1",
                    "USF2", "RAR", "RXR", "EAR2", "NFYA")
markerEnhMotifs <- getFeatures(proj5,
    select = paste(enhancerMotifs, collapse="|"),
    useMatrix = "Motif_LSI_HarmonyMatrix")
markerEnhMotifs

# Get just z scores
markerEnhMotifs <- grep("z:", markerEnhMotifs, value = TRUE)
markerEnhMotifs <- markerEnhMotifs[markerEnhMotifs %ni% "z:Mesp1_58"]
markerEnhMotifs <- markerEnhMotifs[markerEnhMotifs %ni% "z:Sp100_714"]
markerEnhMotifs <- markerEnhMotifs[markerEnhMotifs %ni% "z:Sp140_717"]
markerEnhMotifs <- markerEnhMotifs[markerEnhMotifs %ni% "z:Sp110_718"]
markerEnhMotifs <- markerEnhMotifs[markerEnhMotifs %ni% "z:Creb3l3_116"]
markerEnhMotifs <- markerEnhMotifs[markerEnhMotifs %ni% "z:Rarb_663"]
markerEnhMotifs <- markerEnhMotifs[markerEnhMotifs %ni% "z:Rarg_654"]
markerEnhMotifs <- markerEnhMotifs[markerEnhMotifs %ni% "z:Creb3_870"]
markerEnhMotifs <- markerEnhMotifs[markerEnhMotifs %ni% "z:Crebzf_122"]
markerEnhMotifs <- markerEnhMotifs[markerEnhMotifs %ni% "z:Creb3l2_118"]
markerEnhMotifs

p <- plotGroups(ArchRProj = proj5, 
  groupBy = "ManualClusters", 
  colorBy = "Motif_LSI_HarmonyMatrix", 
  name = markerEnhMotifs,
  imputeWeights = getImputeWeights(proj5)
)

p2 <- lapply(seq_along(p), function(x){
  if(x != 1){
    p[[x]] + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6) +
    theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm")) +
    theme(
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank(),
        axis.title.y=element_blank()
    ) + ylab("")
  }else{
    p[[x]] + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6) +
    theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm")) +
    theme(
        axis.ticks.y=element_blank(),
        axis.title.y=element_blank()
    ) + ylab("")
  }
})

# This plot can show us how FOXD1, for example, changes in accessibility across time points by cluster
svglite(paste0(sample_name,
    "_scATAC-seq_Manual-Annotation_markerEnhMotifs_", Sys.Date(), ".svg"),
    height=8, width=24)
do.call(cowplot::plot_grid,
    c(list(nrow = 1, rel_widths = c(2, rep(1, length(p2) - 1))),p2))
dev.off()

p <- plotEmbedding(
    ArchRProj = proj5, 
    colorBy = "Motif_LSI_HarmonyMatrix", 
    name = sort(markerEnhMotifs), 
    embedding = "UMAP_LSI_Harmony",
    imputeWeights = getImputeWeights(proj5)
)
p2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
svglite(paste0(sample_name,
    "_scATAC-seq_Manual-Annotation_enhancerMotifs_UMAP_",
    Sys.Date(), ".svg"), height=8, width=24)
do.call(cowplot::plot_grid, c(list(ncol = 3),p2))
dev.off()

# Compare how deviations compare to inferred gene expression on UMAP
markerEnhRNA <- c("Sp3",  "Nfix", "Rxrb", "Rara", "Rxra",
                  "Rxrg", "Nfya", "Sp1",  "Crem", "Creb1",
                  "Usf2",  "Usf1")

p <- plotEmbedding(
    ArchRProj = proj5, 
    colorBy = "GeneScoreMatrix", 
    name = sort(markerEnhRNA), 
    embedding = "UMAP_LSI_Harmony",
    imputeWeights = getImputeWeights(proj5)
)

p2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
svglite(paste0(sample_name,
    "_scATAC-seq_Manual-Annotation_enhancer_motifs_inferredExpression_UMAP_",
    Sys.Date(), ".svg"), height=8, width=24)
do.call(cowplot::plot_grid, c(list(ncol = 3),p2))
dev.off()

11d. Motif enrichment of reported Ren1 promoter binding TFs

PREP1 = PKNOX1

promoterMotifs <- c("CBF1", "NIC", "LXR", "NFIX", "SP1", "SP3", "HOXA9",
                    "HOXA10", "PBX1", "PKNOX1")

markerPromMotifs <- getFeatures(proj5,
                                select = paste(promoterMotifs, collapse="|"),
                                useMatrix = "Motif_LSI_HarmonyMatrix")
markerPromMotifs

# Get just z scores
markerPromMotifs <- grep("z:", markerPromMotifs, value = TRUE)
markerPromMotifs <- markerPromMotifs[markerPromMotifs %ni% "z:Mesp1_58"]
markerPromMotifs <- markerPromMotifs[markerPromMotifs %ni% "z:Sp100_714"]
markerPromMotifs <- markerPromMotifs[markerPromMotifs %ni% "z:Sp140_717"]
markerPromMotifs <- markerPromMotifs[markerPromMotifs %ni% "z:Sp110_718"]
markerPromMotifs

p <- plotGroups(ArchRProj = proj5, 
  groupBy = "ManualClusters", 
  colorBy = "Motif_LSI_HarmonyMatrix", 
  name = markerPromMotifs,
  imputeWeights = getImputeWeights(proj5)
)

p2 <- lapply(seq_along(p), function(x){
  if(x != 1){
    p[[x]] + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6) +
    theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm")) +
    theme(
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank(),
        axis.title.y=element_blank()
    ) + ylab("")
  }else{
    p[[x]] + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6) +
    theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm")) +
    theme(
        axis.ticks.y=element_blank(),
        axis.title.y=element_blank()
    ) + ylab("")
  }
})

# plot the distribution of chromVAR deviation scores for each cluster.
# This plot can show us how FOXD1, for example, changes in accessibility across time points by cluster
svglite(paste0(sample_name,
               "_scATAC-seq_Manual-Annotation_markerPromMotifs_",
               Sys.Date(), ".svg"),
        height=8, width=24)
do.call(cowplot::plot_grid,
    c(list(nrow = 1, rel_widths = c(2, rep(1, length(p2) - 1))),p2))
dev.off()

# Overlay promoter motifs to our UMAP
p <- plotEmbedding(
    ArchRProj = proj5, 
    colorBy = "Motif_LSI_HarmonyMatrix", 
    name = sort(markerPromMotifs), 
    embedding = "UMAP_LSI_Harmony",
    imputeWeights = getImputeWeights(proj5)
)
p2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
svglite(paste0(sample_name,
               "_scATAC-seq_Manual-Annotation_promoterMotifs_UMAP_",
               Sys.Date(), ".svg"),
        height=8, width=24)
do.call(cowplot::plot_grid, c(list(ncol = 3),p2))
dev.off()

markerPromRNA <- c("Sp3", "Nfix", "Pbx1", "Hoxa9", "Pknox1", "Hoxa10", "Sp1")

p <- plotEmbedding(
    ArchRProj = proj5, 
    colorBy = "GeneScoreMatrix", 
    name = sort(markerPromRNA), 
    embedding = "UMAP_LSI_Harmony",
    imputeWeights = getImputeWeights(proj5)
)

p2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
svglite(paste0(sample_name,
    "_scATAC-seq_Manual-Annotation_promoter_motifs_inferredExpression_UMAP_",
    Sys.Date(), ".svg"), height=8, width=24)
do.call(cowplot::plot_grid, c(list(ncol = 3),p2))
dev.off()

save.image(paste0(sample_name, "_scATAC-seq_Motif-Annotation_",
           Sys.Date(), ".RData"))

11e. Motif footprinting

motifPositions <- getPositions(proj5)
motifPositions

# subset to TFs of interest
motifs <- c("Foxd1", "Notch", "Rbpj", "Sp1", "Sp3", "Cited1", "Six2", "Tagln",
            "Smtn", "Smad1", "E2f1", "Hoxb7", "Nfi", "Nfix", "Nfic", "Myh11",
            "Mef2a", "Mef2b", "Mef2c", "Mef2d", "Smarcc1")

markerMotifs2 <- unlist(lapply(motifs, function(x) grep(x, names(motifPositions), value = TRUE)))
markerMotifs2 <- markerMotifs2[markerMotifs2 %ni% "Sp100_714"]
markerMotifs2 <- markerMotifs2[markerMotifs2 %ni% "Sp140_717"]
markerMotifs2 <- markerMotifs2[markerMotifs2 %ni% "Sp110_718"]
markerMotifs2

seFoot <- getFootprints(
  ArchRProj = proj5, 
  positions = motifPositions[markerMotifs2], 
  groupBy = "ManualClusters"
)

plotFootprints(
  seFoot = seFoot,
  ArchRProj = proj5, 
  normMethod = "Subtract",
  plotName = paste0(sample_name,
                    "_scATAC-seq_Manual-Annotation_footprints-subtract-bias_",
                    Sys.Date()),
  addDOC = FALSE,
  smoothWindow = 5
)

plotFootprints(
  seFoot = seFoot,
  ArchRProj = proj5, 
  normMethod = "Divide",
  plotName = paste0(sample_name,
                    "_scATAC-seq_Manual-Annotation_footprints-divide-bias_",
                    Sys.Date()),
  addDOC = FALSE,
  smoothWindow = 5
)

11f. Feature footprinting

TSS footprints

seTSS <- getFootprints(
  ArchRProj = proj5, 
  positions = GRangesList(TSS = getTSS(proj5)), 
  groupBy = "ManualClusters",
  flank = 2000
)

plotFootprints(
  seFoot = seTSS,
  ArchRProj = proj5, 
  normMethod = "None",
  plotName = paste0(sample_name,
                    "_scATAC-seq_Manual-Annotation_TSS-No-Normalization_",
                    Sys.Date()),
  addDOC = FALSE,
  flank = 2000,
  flankNorm = 100
)