3_ArchR_Exploration.Rmd

---
title: "R Notebook"
---

```{r}
library(ArchR)
library(knitr)
library(TFBSTools)
library(ComplexHeatmap)
library(ggplot2)
library(parallel)

addArchRThreads(threads = 62) 
set.seed(1)
projNC <- loadArchRProject(path = "/data/Austin/10XscATAC/ArchR_data_final_backup/")
projNC@projectMetadata$outputDirectory <- "./Figure_Exports"
```

Identification of Marker Genes
```{r}
markersGS <- getMarkerFeatures(
    ArchRProj = projNC, 
    useMatrix = "GeneScoreMatrix", 
    groupBy = "Clusters",
    bias = c("TSSEnrichment", "log10(nFrags)"),
    testMethod = "wilcoxon")

GS_df <- as.data.frame(getMarkers(markersGS, cutOff = "FDR <= 0.01 & Log2FC >= 1"))
write.csv(GS_df, file = "Figure_Exports/GeneScores_Results_FDR0.01_Log2FC-1.csv")

# Can we use an imputed gene matrix to clean things up a bit?

# It's important to look at local regions when considering genescores.
# Let's look at a few locations.
p1 <- plotBrowserTrack(
    ArchRProj = projNC, 
    groupBy = "Clusters", 
    geneSymbol = c("SOX10"), 
    upstream = 50000,
    downstream = 50000,
    loops = getCoAccessibility(projNC, corCutOff = 0.7))

## Initiate writing to PDF file
pdf("./Figure_Exports/Plots/SOX10_CoAccessibility.pdf", height = 6, width = 8)
grid::grid.draw(p1$SOX10)
dev.off()
```
Let's see if there are any highly correlated motifs to genescores
```{r}

```


```{r}

# markersPeaks <- getMarkerFeatures(
#     ArchRProj = projNC, 
#     useMatrix = "PeakMatrix", 
#     groupBy = "Clusters",
#   bias = c("TSSEnrichment", "log10(nFrags)"),
#   testMethod = "wilcoxon")
# 
# # This checks binary presence of a TF in a given peakset, but isn't chromvAR
# motifsUp <- peakAnnoEnrichment(
#     seMarker = markersPeaks,
#     ArchRProj = projNC,
#     peakAnnotation = "Motif",
#     cutOff = "FDR <= 0.1 & Log2FC >= 0.5"
#   )
# 
# df <- data.frame(TF = rownames(motifsUp), mlog10Padj = assay(motifsUp)[,1])
# df <- df[order(df$mlog10Padj, decreasing = TRUE),]
# df$rank <- seq_len(nrow(df))
# 
# 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
# 


# toplotMotifs <- getFeatures(projNC, select = paste(motifs, collapse="|"), useMatrix = "MotifMatrix")
# 
# toplotMotifs_z <- toplotMotifs[1:23]
# 
# assayNames <- names(SummarizedExperiment::assays(Cluster_Motifs))
# mat <- as.matrix(SummarizedExperiment::assays(Cluster_Motifs)[["Mean"]])
# 
# heatmapMotifs <- plotMarkerHeatmap(
#   seMarker = Cluster_Motifs, 
#   labelMarkers = toplotMotifs_z,
#   transpose = TRUE
# )
# 
# ComplexHeatmap::draw(heatmapMotifs, heatmap_legend_side = "bot", annotation_legend_side = "bot")
# 
# plotVarDev <- getVarDeviations(projNC, name = "MotifMatrix", plot = TRUE)
# var_tfs <- plotVarDev$data
# 
# grey_red <- ArchRPalettes$whiteRed
# grey_red[1] <- "grey"
# 
# motif_cutofss <- plotEmbedding(projNC, embedding = "UMAP", colorBy = "MotifMatrix",
#               name = toplotMotifs_z, pal = grey_red, plotAs = "points", sampleCells = NULL,
#               rastr = F,  quantCut = c(0.64999999, 0.65))
# 
# plotEmbedding(projNC, embedding = "UMAP", colorBy = "MotifMatrix",
#               name = "z:Pou5f1..Sox2_68", pal = ArchRPalettes$greenBlue,
#               rastr = F)
```

Identification of Marker Peaks
```{r}
# getMarkerFeatures on Peaks is broken until 1.0.2
markersPeaks_C1 <- getMarkerFeatures(
    ArchRProj = projNC, 
    useMatrix = "PeakMatrix", 
    groupBy = "Clusters",
    bias = c("TSSEnrichment", "log10(nFrags)"),
    testMethod = "wilcoxon", useGroups = "C1")

getMarkers(markersPeaks_C1)

p <- plotEmbedding(
    ArchRProj = projNC, 
    colorBy = "PeakMatrix", 
    name = getMarkers(markersPeaks_C1), 
    embedding = "UMAP",
    quantCut = c(0.01, 0.95),
    imputeWeights = NULL
)


GS_df <- as.data.frame(getMarkers(markersGS, cutOff = "FDR <= 0.01 & Log2FC >= 1"))
write.csv(GS_df, file = "temp_GeneScores.csv")

# It's important to look at local regions when considering genescores.
# Let's look at a few locations.
p1 <- plotBrowserTrack(
    ArchRProj = projNC, 
    groupBy = "Clusters", 
    geneSymbol = c("DLX4","SLC26A9","CLDN9"), 
    upstream = 50000,
    downstream = 50000,
    loops = getCoAccessibility(projNC, corCutOff = 0.7))

grid::grid.newpage()
grid::grid.draw(p1$DLX4)

```


# Plotting some motifs
```{r}
motifs <- c("SOX10","PAX3")

grey_red <- ArchRPalettes$whiteRed
grey_red[1] <- "grey"

plotting_motifs <- getFeatures(projNC, select = paste(motifs, collapse="|"), useMatrix = "MotifMatrix")

plotting_motifs_z <- paste0("z:",plotting_motifs)

# Problem with threads > 1...
motif_plot_nocutoff <- plotEmbedding(projNC,
                                     embedding = "UMAP",
                                     colorBy = "MotifMatrix",
                                     name = plotting_motifs_z,
                                     plotAs = "points",
                                     sampleCells = NULL,
                                     rastr = F,
                                     keepAxis = F,
                                     threads = 1)
motif_plot_nocutoff
```