MIC functional enrichment

Last update: 24-10-18

Last markdown compiled: 24-10-18

knitr::opts_chunk$set(warning = FALSE, message = FALSE, results = FALSE) 

# dirs
## Sebastian @Linux
dir_moa <- "/home/isar/Dropbox/scn_moa/"
## Sebastian @MAC
# dir_moa <- "/Users/home_sh/Library/CloudStorage/Dropbox/scn_moa"

# data, revised 2023.11.29
# exp
load(file = paste0(dir_moa, "/data/final/", "gnp.pro.RData"))
load(file = paste0(dir_moa, "/data/final/", "gnp.mrn.RData"))
load(file = paste0(dir_moa, "/data/final/", "gnp.mic.RData"))

# results
load(file = paste0(dir_moa, "/data/final/", "dex.all.RData"))
load(file = paste0(dir_moa, "/data/final/", "timeomics.singles.RData"))
load(file = paste0(dir_moa, "/data/final/", "res.wgcna.RData"))


## packages
source(paste0(dir_moa, "/scripts/", "scn_moa_fct.R"))

# load packages
sapply(packages_general, require, character.only = TRUE)
sapply(packages_enrichment, require, character.only = TRUE)
#sapply(packages_pca, require, character.only = TRUE)
#sapply(packages_cluster, require, character.only = TRUE)
#sapply(packages_heatmap, require, character.only = TRUE)

library(enrichMiR)
library(STRINGdb)
library(rbioapi)

# graphs
global_size = 21

# WGCAN results dont have the correct mic data
load(file = paste0(dir_moa, "/data/final/", "res_wgcna_mic.RData")) # inset into wgcna!


fan <- gnp.mrn$fan

fct

format import data

add_mic.Fc <- function(comp, eMIR){
  fc.mic_comp <- fc.mic %>% filter(comparison == comp)
  eMIR$logFC <- fc.mic_comp$logFC[match(eMIR$miRNA, fc.mic_comp$fid)]
  return(eMIR)
}

# modify colnames
edit_colnames_eMIR <- function(eMIR){
  names(eMIR)[1] <- "motif"
  #names(eMIR)[3] <- "miRNA"
  eMIR$members <- eMIR$miRNA
  eMIR$miRNA <- gsub("\\;.*","",eMIR$miRNA)
  eMIR$members <- gsub(";([A-Za-z])", "; \\1", eMIR$members)
  return(eMIR)
}

bubble plots

#function for reverse log-scale the y-axis
revlog_trans <- function(base = exp(1)){
  ## Define the desired transformation.
  trans <- function(x){
  -log(x, base)
                }
  ## Define the reverse of the desired transformation
  inv <- function(x){
  base^(-x)
                }
  ## Creates the transformation
    trans_new(paste("revlog-", base, sep = ""),
  trans, ## The transformation function (can be defined using anonymous functions)
  inv,  ## The reverse of the transformation
              log_breaks(base = base), ## default way to define the scale breaks
  domain = c(1e-100, Inf) ## The domain over which the transformation is valued
             )
}



plot_bouble.mic.negPos <- function(eMIR, neg = FALSE) {
  plot <- eMIR %>%
    ggplot(aes(x = enrichment, y = pvalue)) + # pvalue
    geom_point(aes(size = set_size, color = logFC)) +
    scale_y_reverse() +
    scale_color_gradient2(low = 'purple', mid = 'orange', high =  '#E01631', midpoint = 0) +
    scale_y_continuous(
      name = "p-value",
      labels = scales::scientific,
      #breaks = scales::pretty_breaks(breaks),
      trans = revlog_trans(base = 10)) +
    geom_text_repel(aes(label = ifelse(FDR < 0.05, stringr::str_wrap(members, 10), NA)), size = 6, color = 'black', fontface = "bold") +
    geom_text_repel(aes(label = ifelse(FDR > 0.05 & pvalue < 0.05, stringr::str_wrap(members, 10), NA)), size = 6, color = 'black') +
    #ggtitle(title) +
    xlab("Enrichment") +
    guides(size = "none") +
    theme(legend.title = element_text(size = 2.7), 
          legend.text = element_text(size = 3))
  return(plot)
}

plot_bouble.mic.neg <- function(eMIR, neg = FALSE) {
    plot <- eMIR %>%
      ggplot(aes(x = enrichment, y = pvalue)) + # pvalue
      geom_point(aes(size = set_size, color = logFC)) +
      scale_y_reverse() +
      scale_color_gradient2(low = 'blue', mid = 'lightblue', high = 'purple') +
      scale_y_continuous(
        name = "p-value",
        labels = scales::scientific,
        #breaks = scales::pretty_breaks(breaks),
        trans = revlog_trans(base = 10)) +
      geom_text_repel(aes(label = ifelse(FDR < 0.05, stringr::str_wrap(members, 10), NA)), size = 6, color = 'black', fontface = "bold") +
      geom_text_repel(aes(label = ifelse(FDR > 0.05 & pvalue < 0.05, stringr::str_wrap(members, 10), NA)), size = 6, color = 'black') +
      #ggtitle(title) +
      xlab("Enrichment") +
     guides(size = "none") +
    theme(legend.title = element_text(size = 2.7), 
          legend.text = element_text(size = 3))
  return(plot)
}




# 
# ######### below unnneded? 
# 
# # labels
# plot_bouble.mic.label <- function(eMIR, title, neg = FALSE) {
#   minFC <- min(eMIR$logFC, na.rm = T) # 0.7
#   
#   if(minFC < 0) {
#     plot <- eMIR %>%
#       ggplot(aes(x = enrichment, y = pvalue)) + # pvalue
#       geom_point(aes(size = set_size, color = logFC)) +
#       scale_y_reverse() +
#       scale_color_gradient2(low = 'blue', mid = 'blue', high =  '#E01631') +
#       scale_y_continuous(
#           name = "p-value",
#           labels = scales::scientific,
#           #breaks = scales::pretty_breaks(breaks),
#           trans = revlog_trans(base = 10)) +
#       geom_label_repel(aes(label = ifelse(FDR < 0.05, stringr::str_wrap(members, 10), NA)), size = 3.5, color = 'black', fontface = "bold", max.overlaps = Inf) +
#       geom_label_repel(aes(label = ifelse(FDR > 0.05 & pvalue < 0.05, stringr::str_wrap(members, 10), NA)), size = 3.5, color = 'black', max.overlaps = Inf) +
#       ggtitle(title)}
#   
#     if(minFC > 0) {
#       plot <- eMIR %>%
#         ggplot(aes(x = enrichment, y = pvalue)) + # pvalue
#         geom_point(aes(size = set_size, color = logFC)) +
#         scale_y_reverse() +
#         scale_color_gradient2(low = 'blue', mid = 'blue', high =  '#E01631') +
#         scale_y_continuous(
#             name = "p-value",
#             labels = scales::scientific,
#             #breaks = scales::pretty_breaks(breaks),
#             trans = revlog_trans(base = 10)) +
#         geom_label_repel(aes(label = ifelse(FDR < 0.05, stringr::str_wrap(members, 10), NA)), size = 3.5, color = 'black', fontface = "bold", max.overlaps = Inf) +
#         geom_label_repel(aes(label = ifelse(FDR > 0.05 & pvalue < 0.05, stringr::str_wrap(members, 10), NA)), size = 3.5, color = 'black', max.overlaps = Inf) +
#         ggtitle(title)}
#   return(plot)
# }
# 
# 
# get_SYM_of_mic_targets <- function(eMIR){
#   SYMdn <- eMIR %>% filter(FDR < 0.05) %>% pull(genes.down)
#   SYMdn <- unlist(strsplit(SYM, ","))
#   return(SYM)
# }
# 
# 
# get_SYM_of_mic_targets.pv <- function(eMIR){
#   SYM <- eMIR %>% filter(pvalue < 0.05) %>% pull(genes.down)
#   SYM <- unlist(strsplit(SYM, ","))
#   return(SYM)
# }

dose response

# redefine colors
defcols <- enrichMiR:::.siteTypeColors()
defcols['8mer'] <- 'blue'
defcols['7mer-1a'] <- 'tan1'
defcols['7mer-m8'] <- 'firebrick1'


CDplotWrapper.own <- function(eMIR.dea, eMIR.TS){
  # fixed
  by = "auto"
  pvals = FALSE
  line.size = 1.2
  point.size = 0.8
  # dependent
  setName = eMIR.TS$set[1]
  # change data (not sure how)
  dea <- enrichMiR:::.homogenizeDEA(eMIR.dea)
  sets <- enrichMiR:::.list2DF(eMIR.TS)
  sets <- sets[sets$set == eMIR.TS$set[1], ] # im not 100% sure what this does
  sets <- enrichMiR:::.addBestType(sets) # im not 100% sure what this does
  by <- head(intersect(c("best_stype", "score", "sites"),
           colnames(sets)), 1)
  dea <- enrichMiR:::.dea2sig(enrichMiR:::.homogenizeDEA(dea), "logFC")
  # if statements
  if (by == "best_stype" | by == "type") {
    if (by == "best_stype") {
        by <- as.character(sets[["best_stype"]])
        names(by) <- sets[["feature"]]
        by2 <- by[names(dea)]
        by2[is.na(by2)] <- "no site"
    }
    
    bylvls <- levels(as.factor(by2))
    defcols <- enrichMiR:::.siteTypeColors()
    defcols['8mer'] <- 'blue'
    if (all(bylvls %in% names(defcols))) {
        by2 <- factor(by2, intersect(names(defcols), bylvls))
    }
    # if (addN)
    #     bylvls <- levels(by2) <- paste0(levels(by2), " (n=",
    #         as.integer(table(by2)), ")")
    p <- CDplot(dea, by = by2, k = length(bylvls), sameFreq = sameFreq,
        size = line.size, pvals = pvals)
    if (all(bylvls %in% names(defcols))) {
        p <- p + scale_colour_manual(values = defcols[levels(by2)])
    }}
  p.fin <- p + xlab("logFC") + ggtitle(setName)
  # results
  res <- list(
    by2 = by2,
    p.fin = p.fin
  )
  return(res)
}

add.numbers2base <- function(plot, eMIRobject, setSize){
  p.fin <- plot + theme(
    legend.position = c(0.99, 0.01),
    legend.justification = c("right", "bottom")) +
    scale_color_manual(labels = c(paste0("no site (n=", setSize$noSite, ")") ,
                                  paste0("7mer-1a (n=", setSize$seven1a, ")") ,
                                  paste0("7mer-m8 (n=", setSize$sevenM8, ")") ,
                                  paste0("8mer (n=", setSize$eight, ")")),
                       values = defcols[levels(eMIRobject$by2)])
  return(p.fin)
}



### PROTEINS

# all pro features as symbols
fid.pro <- unique(dex.all$dex.pro$dex_stage_pro_vsn$SYMBOL)

# change DEA from mrn to pro
switch.DEA.mrnTOpro <- function(eMIR.DEA_mrn, dex.pro){
  eMIR.DEA_pro <- eMIR.DEA_mrn[rownames(eMIR.DEA_mrn) %in% fid.pro,]
  eMIR.DEA_pro$logFC <- dex.pro$logFC[match(rownames(eMIR.DEA_pro), dex.pro$Gene)]
  eMIR.DEA_pro$PValue <- dex.pro$PValue[match(rownames(eMIR.DEA_pro), dex.pro$Gene)]
  eMIR.DEA_pro$FDR <- dex.pro$FDR[match(rownames(eMIR.DEA_pro), dex.pro$Gene)]
  return(eMIR.DEA_pro)
}




add.numbers2base.pro <- function(plot, eMIRobject, setSize){
  non.reg = 3156 - (setSize$seven1a + setSize$sevenM8 + setSize$eight)
  p.fin <- plot + theme(
    legend.position = c(0.99, 0.01),
    legend.justification = c("right", "bottom")) +
    scale_color_manual(labels = c(paste0("no site (n=", non.reg, ")") ,
                                  paste0("7mer-1a (n=", setSize$seven1a, ")") ,
                                  paste0("7mer-m8 (n=", setSize$sevenM8, ")") ,
                                  paste0("8mer (n=", setSize$eight, ")")),
                       values = defcols[levels(eMIRobject$by2)])
  return(p.fin)
}


# add numbers
# add.numbers2base <- function(plot, eMIRobject, setSize){
#   p.fin <- plot + theme(
#     legend.position = c(0.95, 0.15),
#     legend.justification = c("right", "bottom")) +
#     scale_color_manual(labels = c(paste0("no site (n=", setSize$noSite, ")") ,
#                                   paste0("7mer-1a (n=", setSize$seven1a, ")") ,
#                                   paste0("7mer-m8 (n=", setSize$sevenM8, ")") ,
#                                   paste0("8mer (n=", setSize$eight, ")")),
#                                   values = defcols[levels(eMIRobject$by2)])
#   return(p.fin)
# }


# # CDplotWrapper: modified above 

#define: because package is broken
#sets = TS
#dea = dea
# setName = TS$set[1]

# CDplotWrapper <- function (dea, sets, setName, k = 3, addN = FALSE, by = c("auto", 
#     "sites", "score", "best_stype", "type"), sameFreq = NULL, 
#     line.size = 1.2, point.size = 0.8, checkSynonyms = TRUE, 
#     pvals = FALSE) {
#     message("Preparing inputs...")
#     dea <- enrichMiR:::.homogenizeDEA(dea)
#     # if (checkSynonyms) 
#     #     dea <- enrichMiR:::.applySynonyms(dea, sets)
#     sets <- enrichMiR:::.list2DF(sets)
#     sets <- sets[sets$set == setName, ]
#     #if (is.null(sets$sites)) 
#     #    sets$sites <- 1
#    # by <- match.arg(by)
#     #if (by == "auto" || by == "best_stype") 
#         sets <- enrichMiR:::.addBestType(sets)
#    # if (by == "auto") 
#        by <- head(intersect(c("best_stype", "score", "sites"),
#            colnames(sets)), 1)
#   #  if (by != "type" && !(by %in% colnames(sets))) 
#    #     stop("`by` not available in `sets`.")
#    # if (nrow(sets) == 0) 
#   #      stop("setName not found in `sets`.")
#  #   if (is.null(dim(dea))) {
#   #      if (!is.numeric(dea) || is.null(names(dea))) 
#     #        stop("`dea` should be a named numeric vector or a data.frame containing ", 
#    #             "the results of a differential expression analysis.")
#   #  }
#   #  else {
#        dea <- enrichMiR:::.dea2sig(enrichMiR:::.homogenizeDEA(dea), "logFC")
#   #  }
#     # if (!any(sets$feature %in% names(dea))) 
#     #     stop("There seems to be no overlap between the rows of `dea` and the ", 
#     #         "features of `sets`.")
#     # if (is.null(sameFreq)) 
#     #     sameFreq <- by == "score"
#     if (by == "best_stype" | by == "type") {
#         if (by == "best_stype") {
#             by <- as.character(sets[["best_stype"]])
#             names(by) <- sets[["feature"]]
#             by2 <- by[names(dea)]
#             by2[is.na(by2)] <- "no site"
#         }
#         
#         bylvls <- levels(as.factor(by2))
#         defcols <- enrichMiR:::.siteTypeColors()
#         defcols['8mer'] <- 'blue'
#         if (all(bylvls %in% names(defcols))) {
#             by2 <- factor(by2, intersect(names(defcols), bylvls))
#         }
#         # if (addN)
#         #     bylvls <- levels(by2) <- paste0(levels(by2), " (n=",
#         #         as.integer(table(by2)), ")")
#         p <- CDplot(dea, by = by2, k = length(bylvls), sameFreq = sameFreq,
#             size = line.size, pvals = pvals)
#         if (all(bylvls %in% names(defcols))) {
#             p <- p + scale_colour_manual(values = defcols[levels(by2)])
#         }}
#     #     else {
#     #         bylvls2 <- gsub(" \\(n=[0-9]+\\)", "", bylvls)
#     #         if (all(bylvls2 %in% names(defcols))) {
#     #             p <- p + scale_colour_manual(values = setNames(defcols[bylvls2],
#     #               bylvls))
#     #         }
#     #     }
#     # }
#     # else{
#     #     by <- rowsum(sets[[by]], sets$feature)[, 1]
#     #     by2 <- by[names(dea)]
#     #     by2[is.na(by2)] <- 0
#     #     if (k == 2) {
#     #         ll <- split(dea, by2 != 0)
#     #         names(ll) <- c("non-targets", "targets")
#     #         p <- CDplot(ll, size = line.size, pvals = pvals, 
#     #             ...)
#     #         p <- p + scale_colour_manual(values = c("#000004FF", 
#     #             "#E65D2FFF"))
#     #     }
#     #     else {
#     #         p <- CDplot(dea, by = by2, k = k, sameFreq = sameFreq, 
#     #             size = line.size, pvals = pvals, ...)
#     #         if (k == 3) {
#     #             p <- p + scale_color_manual(values = c("#F5DB4BFF", 
#     #               "#A92E5EFF", "#000004FF"))
#     #         }
#     #     }
#     # }
#     # if (point.size > 0) 
#     #     p <- p + geom_point(size = point.size)
#     p.fin <- p + xlab("logFC") + ggtitle(setName)
# }


get.setSizes.pro <- function(sets){
  seven1a = sum(sets$Sites_7mer_1a)
  sevenM8 = sum(sets$Sites_7mer_m8)
  eight = sum(sets$Sites_8mer)
  noSite = 3156 - (seven1a + sevenM8 + eight)
  res <- list(seven1a = seven1a,
              sevenM8 = sevenM8,
              eight = eight,
              noSite = noSite)
  return(res)
}

get.setSizes.mrn <- function(sets){
  seven1a = sum(sets$Sites_7mer_1a)
  sevenM8 = sum(sets$Sites_7mer_m8)
  eight = sum(sets$Sites_8mer)
  noSite = 8432 - (seven1a + sevenM8 + eight)
  res <- list(seven1a = seven1a,
              sevenM8 = sevenM8,
              eight = eight,
              noSite = noSite)
  return(res)
}

enrichment

cluster profiler

get_SYM_of_mic_targets.all <- function(eMIR, FIDmic){
  SYMdn <- eMIR %>% filter(miRNA == FIDmic) %>% pull(genes.down)
  SYMdn <- unlist(strsplit(SYMdn, ","))
  SYMup <- eMIR %>% filter(miRNA == FIDmic) %>% pull(genes.up)
  SYMup <- unlist(strsplit(SYMup, ","))
  SYM <- c(SYMdn, SYMup)
  return(SYM)
}

get_SYM_of_mic_targets.dn <- function(eMIR, FIDmic){
  SYM <- eMIR %>% filter(miRNA == FIDmic) %>% pull(genes.down)
  SYM <- unlist(strsplit(SYM, ","))
  return(SYM)
}



prep.eMIR.4enrich.mrn <- function(SYMBOLS, COMP){
  # create df with ENTREZ & Fc
  mirTar <- data.frame(ENSG = fan$ENSG[match(SYMBOLS, fan$SYMBOL)])
  mirTar$ENTREZ <- fan$ENTREZ[match(mirTar$ENSG, fan$ENSG)]
  FCs <- dex.all$dex.mrn$dex_stage_mrn_vsn %>% filter(comparison == COMP)
  mirTar$Fc <- FCs$logFC[match(mirTar$ENSG, FCs$fid)]
  # turn into ordered vector
  geneList <- setNames(mirTar$Fc, mirTar$ENTREZ)
  geneList <- sort(geneList, decreasing = TRUE)
  
  return(geneList)
}


prep.eMIR.4enrich.pro <- function(SYMBOLS, COMP){
  # create df with ENTREZ & Fc
  mirTar <- data.frame(ENSG = fan$ENSG[match(SYMBOLS, fan$SYMBOL)])
  mirTar$ENTREZ <- fan$ENTREZ[match(mirTar$ENSG, fan$ENSG)]
  mirTar$SYMBOL <- fan$SYMBOL[match(mirTar$ENSG, fan$ENSG)]
  FCs <- dex.all$dex.pro$dex_stage_pro_vsn %>% filter(comparison == COMP)
  mirTar$Fc <- FCs$logFC[match(mirTar$SYMBOL, FCs$SYMBOL)]
  # turn into ordered vector
  geneList <- setNames(mirTar$Fc, mirTar$ENTREZ)
  geneList <- sort(geneList, decreasing = TRUE)
  
  return(geneList)
}

enrichR

# try enrichR
library(enrichR)

# websiteLive <- getOption("enrichR.live")
# if (websiteLive) {
#     listEnrichrSites()
#     setEnrichrSite("Enrichr") # Human genes   
# }
# 
# if (websiteLive) dbs <- listEnrichrDbs()
# if (websiteLive) (dbs)

dbs <- c(
'Reactome_2016',
'WikiPathways_2016',
'KEGG_2016',
'PPI_Hub_Proteins',
'NCI-Nature_2016',
'Jensen_COMPARTMENTS',
'GO_Molecular_Function_2018',
'GO_Cellular_Component_2018',
'GO_Biological_Process_2018')

# check colnames
# purrr::map(rs.enR_mic.DN_pro, colnames) # all the same!

example

data prep

exp data

# FCs of mic (for bubble plots)
fc.mic <- dex.all$dex.mic$dex_stage_mic_vsn


##### MICs

# load
load(file = paste0(dir_moa, "/data/final/preps/timeomics/", "centralities.RData"))

# reduce to 90%
centrality.mic.90 <- centrality.mic %>% select(ENSG, traj, centrality) %>% 
  group_by(traj) %>% 
  arrange(traj, desc(centrality)) %>% 
  filter(centrality > quantile(centrality, .10))

# compare 90% to unfiltered (I feel that we need more IDs as the webiste filters out a lot again)
centrality.mic.all <- centrality.mic %>% select(ENSG, traj, centrality) %>% 
  group_by(traj) %>% 
  arrange(traj, desc(centrality))

setdiff(centrality.mic.all[centrality.mic.all$traj == 'up',]$ENSG, centrality.mic.90[centrality.mic.90$traj == 'up',]$ENSG) # 7

setdiff(centrality.mic.all[centrality.mic.all$traj == 'dn',]$ENSG, centrality.mic.90[centrality.mic.90$traj == 'dn',]$ENSG) # 8


# copy to clip for website . use all becuase trajectories re very clean in mics
# centrality.mic.all %>% filter(traj == 'dn') %>% pull(ENSG) %>% write_clip()
# centrality.mic.all %>% filter(traj == 'up') %>% pull(ENSG) %>% write_clip()
# centrality.mic.all %>% filter(traj == 'dnUp') %>% pull(ENSG) %>% write_clip()
# centrality.mic.all %>% filter(traj == 'upDn') %>% pull(ENSG) %>% write_clip()
# centrality.mic.all %>% filter(traj == 'upDn late') %>% pull(ENSG) %>% write_clip()

export

fct

##### mrn & pro

### DEX
dex.mrn <- dex.all$dex.mrn$dex_stage_mrn_vsn
dex.pro <- dex.all$dex.pro$dex_stage_pro_vsn

# add FC & pval
create.dex.4enrichMIR <- function(dex, COMP){
  dex.comp <- dex %>% filter(comparison == COMP)
  dex.comp %<>% select(Gene = 'SYMBOL', logFC, PValue= 'P.Value', FDR= 'adj.P.Val')
  return(dex.comp)
}


### TRAJ
centrality.mrn.90 <- centrality.mrn %>% select(ENSG, traj, centrality) %>% 
  group_by(traj) %>% 
  arrange(traj, desc(centrality)) %>% 
  filter(centrality > quantile(centrality, .10))

centrality.pro.90 <- centrality.pro %>% select(UNIPROT = ENSG, traj, centrality) %>% 
  group_by(traj) %>% 
  arrange(traj, desc(centrality)) %>% 
  filter(centrality > quantile(centrality, .10))
centrality.pro.90 <- add.fan2UNIPROT(centrality.pro.90)


# add FC & pval from dex to traj
addFC.2traj.4enrichMIR_mrn <- function(TRAJ, COMP){
  # tra
  foi <- centrality.mrn.90 %>% ungroup %>% filter(traj == TRAJ) %>% select(ENSG)
  foi <- add.fan2ENSG(foi)
  # dex
  dex.comp <- dex.mrn %>% filter(comparison == COMP)
  # add to foi
  foi$logFC <- dex.comp$logFC[match(foi$ENSG, dex.comp$fid)]
  foi$PValue <- dex.comp$P.Value[match(foi$ENSG, dex.comp$fid)]
  foi$FDR <- dex.comp$adj.P.Val[match(foi$ENSG, dex.comp$fid)]
  foi$Gene <- foi$SYMBOL
  # finalize
  res <- foi %>% select(Gene,logFC,PValue,FDR)
  return(res)
}

# add FC & pval from dex to traj
addFC.2traj.4enrichMIR_pro <- function(TRAJ, COMP){
  # tra
  foi <- centrality.pro.90 %>% ungroup %>% filter(traj == TRAJ)
  # dex
  dex.comp <- dex.pro %>% filter(comparison == COMP)
  # add to foi
  foi$logFC <- dex.comp$logFC[match(foi$UNIPROT, dex.comp$fid)]
  foi$PValue <- dex.comp$P.Value[match(foi$UNIPROT, dex.comp$fid)]
  foi$FDR <- dex.comp$adj.P.Val[match(foi$UNIPROT, dex.comp$fid)]
  foi$Gene <- foi$SYMBOL
  # finalize
  res <- foi %>% select(Gene,logFC,PValue,FDR)
  return(res)
}






# 
# 
# ##### dump
# 
# create.traj.4enrichMIR_mrn <- function(timeOmics, dex, COMP){
#   # dex
#   dex.comp <- dex %>% filter(comparison == COMP)
#   # tra
#   tra <- data.frame(ENSG = timeOmics)
#   tra <- add.fan2ENSG(tra)
#   tra$Gene <- tra$SYMBOL
#   tra$logFC <- dex.comp$logFC[match(tra$SYMBOL, dex.comp$SYMBOL)]
#   #tra$logFC <- tra$logFC * -1
#   tra$PValue <- dex.comp$P.Value[match(tra$SYMBOL, dex.comp$SYMBOL)]
#   tra$FDR <- dex.comp$adj.P.Val[match(tra$SYMBOL, dex.comp$SYMBOL)]
#   tra$ENSG <- NULL
#   tra$SYMBOL <- NULL
#   tra$UNIPROT <- NULL
#   return(tra)
# }
# 
# create.traj.4enrichMIR_pro <- function(timeOmics, dex, COMP){
#   # dex
#   dex.comp <- dex %>% filter(comparison == COMP)
#   # tra
#   tra <- data.frame(UNIPROT = timeOmics)
#   tra <- add.fan2UNIPROT(tra)
#   tra$Gene <- tra$SYMBOL
#   tra$logFC <- dex.comp$logFC[match(tra$SYMBOL, dex.comp$SYMBOL)]
#   #tra$logFC <- tra$logFC * -1
#   tra$PValue <- dex.comp$P.Value[match(tra$SYMBOL, dex.comp$SYMBOL)]
#   tra$FDR <- dex.comp$adj.P.Val[match(tra$SYMBOL, dex.comp$SYMBOL)]
#   res <- tra %>% select(Gene, logFC, PValue, FDR)
#   return(res)
# }

up

# mic
# centrality.mic.90 %>% filter(traj == 'up') %>% pull(ENSG) %>% write_clip() # recognized: 19
# up -> diminishing stem cell potential during granulopoiesis 


# mrn
dex.mrn_MBtoPMN <- create.dex.4enrichMIR(dex.mrn, 'MBtoPMN')
dex.mrn_MBtoPMN %>% filter(Gene == "S100A8") # up: ok

write_csv(dex.mrn_MBtoPMN, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_MBtoPMN", ".csv"))
# check: 2x FDR. perfect curves



## pro
# MB to PMN
dex.pro_MBtoPMN <- create.dex.4enrichMIR(dex.pro, 'MBtoPMN')
dex.pro_MBtoPMN %>% filter(Gene == "S100A8") # up: wrong! should be down bc mics are inhibiting!

write_csv(dex.pro_MBtoPMN, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.pro_MBtoPMN", ".csv"))
# res: 1 pval, shitty curve


# PM to MM
dex.pro_PMtoMM <- create.dex.4enrichMIR(dex.pro, 'PMtoMM')
dex.pro_PMtoMM %>% filter(Gene == "S100A8") # up: wrong! should be down bc mics are inhibiting!

write_csv(dex.pro_PMtoMM, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.pro_PMtoMM", ".csv"))
# res: 2 pval, good curve

upDn

# mic
# centrality.mic.90 %>% filter(traj == 'upDn') %>% pull(ENSG) %>% write_clip() # recognized: 11
# upDn -> short inhibition of features in MC

# mrn MB -> MC
dex.mrn_MBtoMC <- create.dex.4enrichMIR(dex.mrn, 'MBtoMC')
dex.mrn_MBtoMC %>% filter(Gene == "S100A8") # up: ok?

write_csv(dex.mrn_MBtoMC, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_MBtoMC", ".csv"))
# check: nothing

# try vs
dex.mrn_MBvsMC <- dex.mrn_MBtoMC
dex.mrn_MBvsMC$logFC <- dex.mrn_MBvsMC$logFC * -1
dex.mrn_MBvsMC %>% filter(Gene == "S100A8") # ok! 

write_csv(dex.mrn_MBvsMC, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_MBvsMC", ".csv")) # 3871 sig
# check: 1 pval, good curve


# pro
dex.pro_MBtoMC <- create.dex.4enrichMIR(dex.pro, 'MBtoMC')
dex.pro_MBvsMC <- dex.pro_MBtoMC
dex.pro_MBvsMC$logFC <- dex.pro_MBvsMC$logFC * -1

write_csv(dex.pro_MBvsMC, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.pro_MBvsMC", ".csv")) # 1005 sig
# no sig

# MCtoMM
dex.pro_MCtoMM <- create.dex.4enrichMIR(dex.pro, 'MCtoMM')
write_csv(dex.pro_MCtoMM, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.pro_MCtoMM", ".csv")) # 1005 sig

dex.pro_MCvsMM <- dex.pro_MCtoMM
dex.pro_MCvsMM$logFC <- dex.pro_MCvsMM$logFC * -1
write_csv(dex.pro_MCvsMM, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.pro_MCvsMM", ".csv")) # 171 sig
# 1 pval, ok curve. saved


# MBtoMM
dex.pro_MBtoMM <- create.dex.4enrichMIR(dex.pro, 'MBtoMM')
write_csv(dex.pro_MBtoMM, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.pro_MBtoMM", ".csv")) # 1097 sig

dnUp

Shortly disinhibiting -> ribosomes, translation

# mic
# centrality.mic.90 %>% filter(traj == 'dnUp') %>% pull(ENSG) %>% write_clip() # recognized: 10
# down -> high in MB and falling in granulopoieis -> inhibiting maturation in MB!

# mrn
dex.mrn_MBtoMM <- create.dex.4enrichMIR(dex.mrn, 'MBtoMM')
write_csv(dex.mrn_MBtoMM, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_MBtoMM", ".csv")) # 5514

dex.mrn_MBvsMM <- dex.mrn_MBtoMM
dex.mrn_MBvsMM$logFC <- dex.mrn_MBvsMM$logFC * -1
write_csv(dex.mrn_MBvsMM, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_MBvsMM", ".csv"))


# pro
dex.pro_MBtoMM <- create.dex.4enrichMIR(dex.pro, 'MBtoMM')
dex.pro_MBvsMM <- dex.pro_MBtoMM
dex.pro_MBvsMM$logFC <- dex.pro_MBvsMM$logFC * -1
write_csv(dex.pro_MBvsMM, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.pro_MBvsMM", ".csv")) # 1097

# pro
dex.pro_MCtoB <- create.dex.4enrichMIR(dex.pro, 'MCtoB')
write_csv(dex.pro_MCtoB, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.pro_MCtoB", ".csv")) # 1097
dex.pro_MCvsB <- dex.pro_MCtoB
dex.pro_MCvsB$logFC <- dex.pro_MCvsB$logFC * -1
write_csv(dex.pro_MCvsB, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.pro_MCvsB", ".csv")) # 1097


dex.pro_MCtoMM <- create.dex.4enrichMIR(dex.pro, 'MCtoMM')
write_csv(dex.pro_MCtoMM, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.pro_MCtoMM", ".csv")) # 171 nothing
dex.pro_MCvsMM <- dex.pro_MCtoMM
dex.pro_MCvsMM$logFC <- dex.pro_MCvsMM$logFC * -1
write_csv(dex.pro_MCvsMM, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.pro_MCvsMM", ".csv")) # 171 nothing

upDn late

# mic
# centrality.mic.90 %>% filter(traj == 'upDn late') %>% pull(ENSG) %>% write_clip() # recognized: 2
# down -> high in MB and falling in granulopoieis -> inhibiting maturation in MB!

# mrn
dex.mrn_MMtoPMN <- create.dex.4enrichMIR(dex.mrn, 'MMtoPMN')
write_csv(dex.mrn_MMtoPMN, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_MMtoPMN", ".csv")) # ok suprisingly

dex.mrn_BtoPMN <- create.dex.4enrichMIR(dex.mrn, 'BtoPMN')
write_csv(dex.mrn_BtoPMN, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_BtoPMN", ".csv")) 

dex.mrn_StoPMN <- create.dex.4enrichMIR(dex.mrn, 'StoPMN')
write_csv(dex.mrn_StoPMN, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_StoPMN", ".csv")) 

dn

mics with dn traj = high in MB and falling in granulopoieis -> inhibiting maturation in MB Look for mics and features of neutrophil maturation

mrn

using conserved TS8

prep

# mics
# centrality.mic.all %>% filter(traj == 'dn') %>% pull(ENSG) %>% write_clip() 
# recognized: conserved: 19 / all: 57
# centrality.mic.90 %>% filter(traj == 'dn') %>% pull(ENSG) %>% write_clip() 
# recognized: 19



# MBtoPMN: all of granulopoiesis
dex.mrn_MBtoPMN <- create.dex.4enrichMIR(dex.mrn, 'MBtoPMN')
dex.mrn_MBtoPMN %>% filter(Gene == "S100A8") # up: wrong! should be down bc mics are inhibiting!

dex.mrn_MBvsPMN <- dex.mrn_MBtoPMN

dex.mrn_MBvsPMN$logFC <- dex.mrn_MBvsPMN$logFC * -1
dex.mrn_MBvsPMN %>% filter(Gene == "S100A8") # ok! 

write_csv(dex.mrn_MBvsPMN, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_MBvsPMN", ".csv"))

# check: 2x FDR with conserved

MBtoPMN

bubble

eMIR_mic.DN_mrn.MBtoPMN <- read_csv(paste0(dir_moa, "/data/final/ext/eMIR/in/mic_dn/mrn/MBvsPMN/", "EnrichMir_hits_siteoverlap.down_2024-05-09.csv"))

# prep
eMIR_mic.DN_mrn.MBtoPMN <- edit_colnames_eMIR(eMIR_mic.DN_mrn.MBtoPMN)
eMIR_mic.DN_mrn.MBtoPMN <- add_mic.Fc('MBtoPMN', eMIR_mic.DN_mrn.MBtoPMN)

# plot
range(eMIR_mic.DN_mrn.MBtoPMN$logFC) # -6 to 3
p.bub_mic.DN_mrn.MBtoPMN_pre <- plot_bouble.mic.neg(eMIR_mic.DN_mrn.MBtoPMN)
p.bub_mic.DN_mrn.MBtoPMN_pre

# format for export
p.bub_mic.DN_mrn.MBtoPMN_fin <- p.bub_mic.DN_mrn.MBtoPMN_pre + 
  theme_light(base_size=20) + 
  theme(legend.title = element_text(size = 9), legend.text = element_text(size = 9), legend.margin=margin(-5,-5,-5,-15)) + 
  labs(title = NULL)

p.bub_mic.DN_mrn.MBtoPMN_fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.bub_mic.DN_mrn.MBtoPMN_fin", ".svg"),
  plot = p.bub_mic.DN_mrn.MBtoPMN_fin,
  device = svg,
  width = 6,
  height = 5)

dose

mrn

# import data & prep
load(paste0(dir_moa, "/data/final/ext/eMIR/in/mic_dn/mrn/MBvsPMN/", "CDplot_AAAGUGC", ".RData"))
eMIR.TS_mic.DN_mrn <- TS
eMIR.DEA_mic.DN_mrn <- dea

# create
obj.dose_mic.DN_mrn <- CDplotWrapper.own(eMIR.DEA_mic.DN_mrn, eMIR.TS_mic.DN_mrn)
p.dose_mic.DN_mrn_pre <- obj.dose_mic.DN_mrn$p.fin + 
                          xlim(-4, 4) + 
                          theme_light(base_size=20)
p.dose_mic.DN_mrn_pre

# numbers
setSize_mic.DN_mrn <- get.setSizes.mrn(eMIR.TS_mic.DN_mrn)
p.dose_mic.DN_mrn.fin <- add.numbers2base(p.dose_mic.DN_mrn_pre, obj.dose_mic.DN_mrn,  setSize_mic.DN_mrn)
p.dose_mic.DN_mrn.fin <- p.dose_mic.DN_mrn.fin + labs(title = NULL)
p.dose_mic.DN_mrn.fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.dose_mic.DN_mrn.fin", ".svg"),
  plot = p.dose_mic.DN_mrn.fin,
  device = svg,
  width = 5.5,
  height = 5)

pro

dex.pro_MBtoPMN <- create.dex.4enrichMIR(dex.pro, 'MBtoPMN')
dex.pro_MBtoPMN %>% filter(Gene == "S100A8") 

dex.pro_MBvsPMN <- dex.pro_MBtoPMN
dex.pro_MBvsPMN$logFC <- dex.pro_MBvsPMN$logFC * -1

eMIR.DEA_mic.DN_pro <- switch.DEA.mrnTOpro(eMIR.DEA_mic.DN_mrn, dex.pro_MBvsPMN)

obj.dose_mic.DN_pro <- CDplotWrapper.own(eMIR.DEA_mic.DN_pro, eMIR.TS_mic.DN_mrn)
p.dose_mic.DN_pro_pre <- obj.dose_mic.DN_pro$p.fin + 
                          xlim(-4, 4) + 
                          theme_light(base_size=20)
p.dose_mic.DN_pro_pre

# numbers
# reduce to pro
eMIR.TS_mic.DN_pro <- eMIR.TS_mic.DN_mrn[eMIR.TS_mic.DN_mrn$feature %in% fid.pro,]

setSize_mic.DN_pro <- get.setSizes.mrn(eMIR.TS_mic.DN_pro)

p.dose_mic.DN_pro.fin <- add.numbers2base.pro(p.dose_mic.DN_pro_pre, obj.dose_mic.DN_mrn,  setSize_mic.DN_pro)

p.dose_mic.DN_pro.fin <- p.dose_mic.DN_pro.fin + labs(title = NULL)
p.dose_mic.DN_pro.fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/supp/", "p.dose_mic.DN_pro.fin", ".svg"),
  plot = p.dose_mic.DN_pro.fin,
  device = svg,
  width = 5.5,
  height = 5)

enrichment

prep

# extract gene SYMBOLS: not sure if up and down or only down???
SYM.mic.DN_mrn.19a <- get_SYM_of_mic_targets.dn(eMIR_mic.DN_mrn.MBtoPMN, 'hsa-miR-19a-3p')
sort(SYM.mic.DN_mrn.19a) # 309

SYM.mic.DN_mrn.106a <- get_SYM_of_mic_targets.dn(eMIR_mic.DN_mrn.MBtoPMN, 'hsa-miR-106a-5p')
sort(SYM.mic.DN_mrn.106a) # 329
                                 

# prep
geneList_mic.DN_mrn.19a <- prep.eMIR.4enrich.mrn(SYM.mic.DN_mrn.19a, "MBtoPMN")
length(geneList_mic.DN_mrn.19a) # 309

geneList_mic.DN_mrn.106a <- prep.eMIR.4enrich.mrn(SYM.mic.DN_mrn.106a, "MBtoPMN")
length(geneList_mic.DN_mrn.106a) # 329

hsa-miR-106a-5p

# GO
# rs.GO_mic.DN_mrn.106a <- enrichGO(gene         = names(geneList_mic.DN_mrn.106a),
#                 OrgDb         = org.Hs.eg.db,
#                 ont           = "ALL",
#                 pAdjustMethod = "BH",
#                 pvalueCutoff  = 0.05,
#                 readable = TRUE,
#                 universe = universe_mrn.ENTREZ)
# 
# save(rs.GO_mic.DN_mrn.106a,
#      file = paste0(dir_moa, "/data/final/preps/enrichment/", "rs.GO_mic.DN_mrn.106a.RData"))

load(file = paste0(dir_moa,  "/data/final/preps/enrichment/", "rs.GO_mic.DN_mrn.106a.RData"))

length(rs.GO_mic.DN_mrn.106a@result$Description) # 107
rs.GO_mic.DN_mrn.106a@result %>% filter(p.adjust < 0.05) %>% nrow() # 107
clusterProfiler::dotplot(rs.GO_mic.DN_mrn.106a, showCategory= 9) # 

rs.GO_mic.DN_mrn.106a@result %>% filter(p.adjust < 0.05) %>% pull(Description)
rs.GO_mic.DN_mrn.106a@result %>% filter(p.adjust < 0.05) %>% view()

# dot plot
terms.mrn.106.dot <- c(
  "cell growth",
  "positive regulation of autophagy",
  "myeloid cell differentiation",
  "regulation of JNK cascade",
  "Wnt signaling pathway",
  "histone modification"
)

p.GO.mic.DN_mrn.106 <- clusterProfiler::dotplot(rs.GO_mic.DN_mrn.106a, showCategory=terms.mrn.106.dot, font.size = 14) #
p.GO.mic.DN_mrn.106

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.GO.mic.DN_mrn.106", ".svg"),
  plot = p.GO.mic.DN_mrn.106,
  device = svg,
  width = 5,
  height = 5)

# cNet

terms.mrn.106.cnt <- c(
  "histone modification",
  "myeloid cell differentiation"
)

p.cnet.mic.DN_mrn.106a <- rs.GO_mic.DN_mrn.106a %>% 
  filter(Description %in% terms.mrn.106.cnt) %>% 
  cnetplot(foldChange=geneList_mic.DN_mrn.106a)

p.cnet.mic.DN_mrn.106a

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.cnet.mic.DN_mrn.106a", ".svg"),
  plot = p.cnet.mic.DN_mrn.106a,
  device = svg,
  width = 7,
  height = 7)



# KEGG
# rs.KG_mic.DN_mrn.106a <- enrichKEGG(gene         = names(geneList_mic.DN_mrn.106a),
#                  organism     = 'hsa',
#                  pvalueCutoff = 0.01,
#                  qvalueCutoff  = 0.05,
#                  pAdjustMethod = "BH",
#                  universe = universe_mrn.ENTREZ)


## check
# rs.KG_mic.DN_mrn.106a@result$Description
# 
# clusterProfiler::dotplot(rs.KG_mic.DN_mrn.106a, showCategory=20) # much better!
# 
# # automatic
# clusterProfiler::dotplot(rs.KG_mic.DN_mrn, showCategory=7) # perfect!






#browseKEGG(rs.KG_mic.DN_mrn, 'hsa05221') # AML
#browseKEGG(rs.KG_mic.DN_mrn, 'hsa05220') # CML


# p.cnet.ME1 <- cnetplot(rs.GO_ME1.mrn, foldChange=geneList.ME1.mrn, 
#          cex.params = list(category_node = 2, gene_node = 2, category_label = 2, gene_label = 2))
# 
# p.cnet.ME1 <- p.cnet.ME1 + guides(size = 'none')


# manual selection?
# rs.KG_mic.DN_mrn@result %>% filter(p.adjust < 0.05) %>% view()
# terms.mic.DN_mrn.KG <- c(GnRH signaling pathway,
#                          
#                           'PD-L1 expression and PD-1 checkpoint pathway in cancer',
#                          'mTOR signaling pathway',
#                          'FoxO signaling pathway')

hsa-miR-19

# GO
# rs.GO_mic.DN_mrn.19a <- enrichGO(gene         = names(geneList_mic.DN_mrn.19a),
#                 OrgDb         = org.Hs.eg.db,
#                 ont           = "ALL",
#                 pAdjustMethod = "BH",
#                 pvalueCutoff  = 0.05,
#                 readable = TRUE,
#                 universe = universe_mrn.ENTREZ)
# 
# save(rs.GO_mic.DN_mrn.19a,
#      file = paste0(dir_moa, "/data/final/preps/enrichment/", "rs.GO_mic.DN_mrn.19a.RData"))

load(file = paste0(dir_moa,  "/data/final/preps/enrichment/", "rs.GO_mic.DN_mrn.19a.RData"))



# total number
rs.GO_mic.DN_mrn.19a@result %>% filter(p.adjust < 0.05) %>% nrow() # 38
# top 9
clusterProfiler::dotplot(rs.GO_mic.DN_mrn.19a, showCategory= 9) # 

# overview
rs.GO_mic.DN_mrn.19a@result %>% filter(p.adjust < 0.05) %>% view()

rs.GO_mic.DN_mrn.19a@result %>% filter(p.adjust < 0.05) %>% pull(Description)


# dot plot
terms.mrn.106.dot <- c(
  "cell growth",
  "positive regulation of autophagy",
  "myeloid cell differentiation",
  "regulation of JNK cascade",
  "Wnt signaling pathway",
  "histone modification"
)

# clusterProfiler::dotplot(rs.GO_mic.DN_mrn.106a, showCategory= terms.mrn.106.dot) # 

p.GO.mic.DN_mrn.106 <- clusterProfiler::dotplot(rs.GO_mic.DN_mrn.106a, showCategory=terms.mrn.106.dot, font.size = 14) #

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.GO.mic.DN_mrn.106", ".svg"),
  plot = p.GO.mic.DN_mrn.106,
  device = svg,
  width = 5,
  height = 5)

# cNet

terms.mrn.106.cnt <- c(
  "histone modification",
  "myeloid cell differentiation"
)

p.cnet.mic.DN_mrn.106a <- rs.GO_mic.DN_mrn.106a %>% 
  filter(Description %in% terms.mrn.106.cnt) %>% 
  cnetplot(foldChange=geneList_mic.DN_mrn.106a)

p.cnet.mic.DN_mrn.106a

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.cnet.mic.DN_mrn.106a", ".svg"),
  plot = p.cnet.mic.DN_mrn.106a,
  device = svg,
  width = 7,
  height = 7)

pro

using all TS8

prep

# mics
# centrality.mic.all %>% filter(traj == 'dn') %>% pull(ENSG) %>% write_clip() 
# recognized: conserved: 19 
# centrality.mic.90 %>% filter(traj == 'dn') %>% pull(ENSG) %>% write_clip() 
# recognized: 19


# MBvsMC: might be better bc more


# MBtoPMN: 
dex.pro_MBtoPMN <- create.dex.4enrichMIR(dex.pro, 'MBtoPMN')
dex.pro_MBtoPMN %>% filter(Gene == "S100A8") # up: wrong! should be down bc mics are inhibiting!

dex.pro_MBvsPMN <- dex.pro_MBtoPMN
dex.pro_MBvsPMN$logFC <- dex.pro_MBvsPMN$logFC * -1
dex.pro_MBvsPMN %>% filter(Gene == "S100A8") # ok! 

write_csv(dex.pro_MBvsPMN, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.pro_MBvsPMN", ".csv"))

# check: 2x pval with all TS8!

MBvsPMN

bubble

eMIR_mic.DN_pro.MBtoPMN <- read_csv(paste0(dir_moa, "/data/final/ext/eMIR/in/mic_dn/pro/MBvsPMN/", "EnrichMir_hits_siteoverlap.down_2024-05-20.csv"))

# prep
eMIR_mic.DN_pro.MBtoPMN <- edit_colnames_eMIR(eMIR_mic.DN_pro.MBtoPMN)
eMIR_mic.DN_pro.MBtoPMN <- add_mic.Fc('MBtoPMN', eMIR_mic.DN_pro.MBtoPMN)

# plot
range(eMIR_mic.DN_pro.MBtoPMN$logFC) # -6 to 3
p.bub_mic.DN_pro.MBtoPMN_pre <- plot_bouble.mic.neg(eMIR_mic.DN_pro.MBtoPMN)

# format for export
p.bub_mic.DN_pro.MBtoPMN_fin <- p.bub_mic.DN_pro.MBtoPMN_pre + 
  theme_light(base_size=20) + 
  theme(legend.title = element_text(size = 9), legend.text = element_text(size = 9), legend.margin=margin(-5,-5,-5,-15)) + 
  labs(title = NULL)

p.bub_mic.DN_pro.MBtoPMN_fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/supp/mic/", "p.bub_mic.DN_pro.MBtoPMN_fin", ".svg"),
  plot = p.bub_mic.DN_pro.MBtoPMN_fin,
  device = svg,
  width = 6,
  height = 5)

dose

# import data & prep
load(paste0(dir_moa, "/data/final/ext/eMIR/in/mic_dn/pro/MBvsPMN/", "CDplot_GGAGGAA", ".RData"))
eMIR.TS_mic.DN_pro <- TS
eMIR.DEA_mic.DN_pro <- dea

# create
obj.dose_mic.DN_pro <- CDplotWrapper.own(eMIR.DEA_mic.DN_pro, eMIR.TS_mic.DN_pro)
p.dose_mic.DN_pro_pre <- obj.dose_mic.DN_pro$p.fin + 
                          xlim(-3, 3) + 
                          theme_light(base_size=20)
p.dose_mic.DN_pro_pre

# numbers
setSize_mic.DN_pro <- get.setSizes.pro(eMIR.TS_mic.DN_pro)
p.dose_mic.DN_pro.fin <- add.numbers2base(p.dose_mic.DN_pro_pre, obj.dose_mic.DN_pro,  setSize_mic.DN_pro)
p.dose_mic.DN_pro.fin <- p.dose_mic.DN_pro.fin + labs(title = NULL)
p.dose_mic.DN_pro.fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/supp/mic/", "p.dose_mic.DN_pro.fin", ".svg"),
  plot = p.dose_mic.DN_pro.fin,
  device = svg,
  width = 5.5,
  height = 5)

enrichment

prep

# extract gene SYMBOLS: not sure if up and down or only down???
SYM.mic.DN_pro.766 <- get_SYM_of_mic_targets.dn(eMIR_mic.DN_pro.MBtoPMN, 'hsa-miR-766-5p')
geneList_mic.DN_pro.766 <- prep.eMIR.4enrich.pro(SYM.mic.DN_pro.766, "MBtoPMN")
length(geneList_mic.DN_pro.766) # 157

hsa-miR-766-5p

# GO
# rs.GO_mic.DN_pro.766 <- enrichGO(gene = names(geneList_mic.DN_pro.766),
#                 OrgDb         = org.Hs.eg.db,
#                 ont           = "ALL",
#                 pAdjustMethod = "BH",
#                 pvalueCutoff  = 0.01,
#                 readable = TRUE,
#                 universe = universe_pro.ENTREZ)
# 
# save(rs.GO_mic.DN_pro.766,
#      file = paste0(dir_moa, "/data/final/preps/enrichment/", "rs.GO_mic.DN_pro.766.RData"))

load(file = paste0(dir_moa,  "/data/final/preps/enrichment/", "rs.GO_mic.DN_pro.766.RData"))


rs.GO_mic.DN_pro.766@result %>% filter(p.adjust < 0.05) %>% nrow() # 102
clusterProfiler::dotplot(rs.GO_mic.DN_pro.766, showCategory= 10) # 

rs.GO_mic.DN_pro.766@result %>% filter(p.adjust < 0.05) %>% view() # 239

rs.GO_mic.DN_pro.766@result %>% filter(p.adjust < 0.05) %>% pull(Description)


# dot plot
terms.pro.766.dot <- c(
'azurophil granule',
'specific granule',
'tertiary granule',
'myeloid leukocyte activation',
'leukocyte migration',
'leukocyte chemotaxis'
)

clusterProfiler::dotplot(rs.GO_mic.DN_pro.766, showCategory= terms.pro.766.dot) # 

p.GO.mic.DN_pro.766 <- clusterProfiler::dotplot(rs.GO_mic.DN_pro.766, showCategory=terms.pro.766.dot, font.size = 14) #

ggsave(
  filename = paste0(dir_moa, "/figures/final/supp/mic/", "p.GO.mic.DN_pro.766", ".svg"),
  plot = p.GO.mic.DN_pro.766,
  device = svg,
  width = 5.5,
  height = 5)

# cNet

terms.pro.766.cnt <- c(
  'azurophil granule',
  'leukocyte migration'
)

p.cnet.mic.DN_pro.766 <- rs.GO_mic.DN_pro.766 %>% 
  filter(Description %in% terms.pro.766.cnt) %>% 
  cnetplot(foldChange=geneList_mic.DN_pro.766)

p.cnet.mic.DN_pro.766

ggsave(
  filename = paste0(dir_moa, "/figures/final/supp/mic/", "p.cnet.mic.DN_pro.766", ".svg"),
  plot = p.cnet.mic.DN_pro.766,
  device = svg,
  width = 7,
  height = 7)


# 
# # KEGG
# rs.KG_mic.DN_pro.106a <- enrichKEGG(gene         = names(geneList_mic.DN_pro.106a),
#                  organism     = 'hsa',
#                  pvalueCutoff = 0.01,
#                  qvalueCutoff  = 0.05,
#                  pAdjustMethod = "BH",
#                  universe = universe_pro.ENTREZ)
# 
# 
# ## check
# # GO
# 
# 
# 
# # KEG
# rs.KG_mic.DN_pro.106a@result$Description
# 
# clusterProfiler::dotplot(rs.KG_mic.DN_pro.106a, showCategory=20) # much better!
# 
# 
# 
# # automatic
# clusterProfiler::dotplot(rs.KG_mic.DN_pro, showCategory=7) # perfect!
# 
# 
# 
# 
# 
# 
# #browseKEGG(rs.KG_mic.DN_pro, 'hsa05221') # AML
# #browseKEGG(rs.KG_mic.DN_pro, 'hsa05220') # CML
# 
# 
# # p.cnet.ME1 <- cnetplot(rs.GO_ME1.pro, foldChange=geneList.ME1.pro, 
# #          cex.params = list(category_node = 2, gene_node = 2, category_label = 2, gene_label = 2))
# # 
# # p.cnet.ME1 <- p.cnet.ME1 + guides(size = 'none')
# 
# 
# # manual selection?
# # rs.KG_mic.DN_pro@result %>% filter(p.adjust < 0.05) %>% view()
# # terms.mic.DN_pro.KG <- c(GnRH signaling pathway,
# #                          
# #                           'PD-L1 expression and PD-1 checkpoint pathway in cancer',
# #                          'mTOR signaling pathway',
# #                          'FoxO signaling pathway')

up

mics with up traj = rising in granulopoieis -> down regulating “unneeded” stuff

mrn

using conserved TS8

prep

# mics
# centrality.mic.all %>% filter(traj == 'up') %>% pull(ENSG) %>% write_clip() 
# recognized: conserved: 23 = used!
# all: 57


# MBtoPMN: all of granulopoiesis but to, not vs!
dex.mrn_MBtoPMN <- create.dex.4enrichMIR(dex.mrn, 'MBtoPMN')
dex.mrn_MBtoPMN %>% filter(Gene == "S100A8") # up: wrong! should be down bc mics are inhibiting!

# write_csv(dex.mrn_MBtoPMN, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_MBtoPMN", ".csv"))

# check: 1x FDR with conserved

MBtoPMN

bubble

eMIR_mic.UP_mrn.MBtoPMN <- read_csv(paste0(dir_moa, "/data/final/ext/eMIR/in/mic_up/mrn/MBtoPMN/", "EnrichMir_hits_siteoverlap.down_2024-05-20.csv"))

# prep
eMIR_mic.UP_mrn.MBtoPMN <- edit_colnames_eMIR(eMIR_mic.UP_mrn.MBtoPMN)
eMIR_mic.UP_mrn.MBtoPMN <- add_mic.Fc('MBtoPMN', eMIR_mic.UP_mrn.MBtoPMN)

# plot
range(eMIR_mic.UP_mrn.MBtoPMN$logFC) # -1 to 5
p.bub_mic.UP_mrn.MBtoPMN_pre <- plot_bouble.mic.negPos(eMIR_mic.UP_mrn.MBtoPMN)

# format for export
p.bub_mic.UP_mrn.MBtoPMN_fin <- p.bub_mic.UP_mrn.MBtoPMN_pre + 
  theme_light(base_size=20) + 
  theme(legend.title = element_text(size = 9), legend.text = element_text(size = 9), legend.margin=margin(-5,-5,-5,-15)) + 
  labs(title = NULL)

p.bub_mic.UP_mrn.MBtoPMN_fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.bub_mic.UP_mrn.MBtoPMN_fin", ".svg"),
  plot = p.bub_mic.UP_mrn.MBtoPMN_fin,
  device = svg,
  width = 6,
  height = 5)

dose

mrn

# import data & prep
load(paste0(dir_moa, "/data/final/ext/eMIR/in/mic_up/mrn/MBtoPMN/", "CDplot_GGGUCUU", ".RData"))
eMIR.TS_mic.UP_mrn <- TS
eMIR.DEA_mic.UP_mrn <- dea

# create
obj.dose_mic.UP_mrn <- CDplotWrapper.own(eMIR.DEA_mic.UP_mrn, eMIR.TS_mic.UP_mrn)
p.dose_mic.UP_mrn_pre <- obj.dose_mic.UP_mrn$p.fin + 
                          xlim(-3, 3) + 
                          theme_light(base_size=20)
p.dose_mic.UP_mrn_pre

# numbers
setSize_mic.UP_mrn <- get.setSizes.mrn(eMIR.TS_mic.UP_mrn)
p.dose_mic.UP_mrn.fin <- add.numbers2base(p.dose_mic.UP_mrn_pre, obj.dose_mic.UP_mrn,  setSize_mic.UP_mrn)
p.dose_mic.UP_mrn.fin <- p.dose_mic.UP_mrn.fin + labs(title = NULL)
p.dose_mic.UP_mrn.fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.dose_mic.UP_mrn.fin", ".svg"),
  plot = p.dose_mic.UP_mrn.fin,
  device = svg,
  width = 5.5,
  height = 5)

pro

eMIR.DEA_mic.UP_pro <- switch.DEA.mrnTOpro(eMIR.DEA_mic.UP_mrn, dex.pro_MBtoPMN)

obj.dose_mic.UP_pro <- CDplotWrapper.own(eMIR.DEA_mic.UP_pro, eMIR.TS_mic.UP_mrn)
p.dose_mic.UP_pro_pre <- obj.dose_mic.UP_pro$p.fin + 
                          xlim(-4, 4) + 
                          theme_light(base_size=20)
p.dose_mic.UP_pro_pre

# numbers
# reduce to pro
eMIR.TS_mic.UP_pro <- eMIR.TS_mic.UP_mrn[eMIR.TS_mic.UP_mrn$feature %in% fid.pro,]

setSize_mic.UP_pro <- get.setSizes.mrn(eMIR.TS_mic.UP_pro)

p.dose_mic.UP_pro.fin <- add.numbers2base.pro(p.dose_mic.UP_pro_pre, obj.dose_mic.UP_mrn,  setSize_mic.UP_pro)

p.dose_mic.UP_pro.fin <- p.dose_mic.UP_pro.fin + labs(title = NULL)
p.dose_mic.UP_pro.fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/supp/mic/", "p.dose_mic.UP_pro.fin", ".svg"),
  plot = p.dose_mic.UP_pro.fin,
  device = svg,
  width = 5.5,
  height = 5)

enrichment

prep

# genes down only
SYM.mic.UP_mrn.193 <- get_SYM_of_mic_targets.dn(eMIR_mic.UP_mrn.MBtoPMN, 'hsa-miR-193a-5p')
geneList_mic.UP_mrn.193 <- prep.eMIR.4enrich.mrn(SYM.mic.UP_mrn.193, "MBtoPMN")
length(geneList_mic.UP_mrn.193) # 42

# check
eMIR_mic.UP_mrn.MBtoPMN %>% filter(miRNA == 'hsa-miR-193a-5p') %>% pull(genes.down) %>% str_count(",")
eMIR_mic.UP_mrn.MBtoPMN %>% filter(miRNA == 'hsa-miR-193a-5p') %>% pull(genes.up)%>% str_count(",")

# 42+22 # 64: why are there less genes in the table than in the TS?

# 15+35+36 # 86 ???. 85 because one has two binding sites. ok


# overlaps TS and table
feature_TS <- droplevels(eMIR.TS_mic.UP_mrn$feature)

feature.tab.down <- eMIR_mic.UP_mrn.MBtoPMN %>% filter(miRNA == 'hsa-miR-193a-5p') %>% pull(genes.down) %>% strsplit(",") %>% unlist()

feature.tab.up <- eMIR_mic.UP_mrn.MBtoPMN %>% filter(miRNA == 'hsa-miR-193a-5p') %>% pull(genes.up) %>% strsplit(",") %>% unlist()

setdiff(feature_TS, c(feature.tab.down, feature.tab.up)) # 21
setdiff(c(feature.tab.down, feature.tab.up), feature_TS) # 0
# maybe the table contains all the genes and the TS only my specific genes?


eMIR_mic.UP_mrn.MBtoPMN %>% filter(miRNA == 'hsa-miR-15b-5p') %>% pull(genes.down) %>% strsplit(",") %>% unlist()

hsa-miR-193

# GO
# rs.GO_mic.UP_mrn.193 <- enrichGO(gene = names(geneList_mic.UP_mrn.193),
#                 OrgDb         = org.Hs.eg.db,
#                 ont           = "ALL",
#                 pAdjustMethod = "BH",
#                 pvalueCutoff  = 0.05,
#                 readable = TRUE)
# 
# save(rs.GO_mic.UP_mrn.193,
#      file = paste0(dir_moa, "/data/final/preps/enrichment/", "rs.GO_mic.UP_mrn.193.RData"))

load(file = paste0(dir_moa,  "/data/final/preps/enrichment/", "rs.GO_mic.UP_mrn.193.RData"))



rs.GO_mic.UP_mrn.193@result %>% filter(p.adjust < 0.05) %>% nrow() # 9

clusterProfiler::dotplot(rs.GO_mic.UP_mrn.193, showCategory= 9) # redundant

clusterProfiler::dotplot(rs.GO_mic.UP_mrn.193, showCategory= 5) # better

# dot plot
p.GO.mic.DN_mrn.193 <- clusterProfiler::dotplot(rs.GO_mic.UP_mrn.193, showCategory=5, font.size = 14) #

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.GO.mic.DN_mrn.193", ".svg"),
  plot = p.GO.mic.DN_mrn.193,
  device = svg,
  width = 5.5,
  height = 5)


# cNet
rs.GO_mic.UP_mrn.193@result %>% filter(p.adjust < 0.05) %>% pull(Description)

terms.mrn193.cnt <- c(
  "ATPase complex",
  "SWI/SNF superfamily-type complex"
)

p.cnet.mic.UP_mrn.193 <- rs.GO_mic.UP_mrn.193 %>% 
  filter(Description %in% terms.mrn193.cnt) %>% 
  cnetplot(foldChange=geneList_mic.UP_mrn.193)

p.cnet.mic.UP_mrn.193

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.cnet.mic.UP_mrn.193", ".svg"),
  plot = p.cnet.mic.UP_mrn.193,
  device = svg,
  width = 7,
  height = 7)

pro

using all TS8

prep

# mics
# centrality.mic.all %>% filter(traj == 'up') %>% pull(ENSG) %>% write_clip() 
# recognized: conserved: 19 

#centrality.mic.90 %>% filter(traj == 'dn') %>% pull(ENSG) %>% write_clip() 
# recognized: 19

# MBvsMC: might be better bc more


# MBtoPMN: 
dex.pro_MBtoPMN <- create.dex.4enrichMIR(dex.pro, 'MBtoPMN')
dex.pro_MBtoPMN %>% filter(Gene == "S100A8") # up: wrong! should be down bc mics are inhibiting!

# write_csv(dex.pro_MBtoPMN, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.pro_MBtoPMN", ".csv"))

# check: 2x pval with conserved TS8!

MBtoPMN

bubble

eMIR_mic.UP_pro.MBtoPMN <- read_csv(paste0(dir_moa, "/data/final/ext/eMIR/in/mic_up/pro/MBtoPMN/", "EnrichMir_hits_siteoverlap.down_2024-05-20.csv"))

# prep
eMIR_mic.UP_pro.MBtoPMN <- edit_colnames_eMIR(eMIR_mic.UP_pro.MBtoPMN)
eMIR_mic.UP_pro.MBtoPMN <- add_mic.Fc('MBtoPMN', eMIR_mic.UP_pro.MBtoPMN)

# plot
range(eMIR_mic.UP_pro.MBtoPMN$logFC) # -6 to 3
p.bub_mic.UP_pro.MBtoPMN_pre <- plot_bouble.mic.negPos(eMIR_mic.UP_pro.MBtoPMN)

# format for export
p.bub_mic.UP_pro.MBtoPMN_fin <- p.bub_mic.UP_pro.MBtoPMN_pre + 
  theme_light(base_size=20) + 
  theme(legend.title = element_text(size = 9), legend.text = element_text(size = 9), legend.margin=margin(-5,-5,-5,-15)) + 
  labs(title = NULL)

p.bub_mic.UP_pro.MBtoPMN_fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.bub_mic.UP_pro.MBtoPMN_fin", ".svg"),
  plot = p.bub_mic.UP_pro.MBtoPMN_fin,
  device = svg,
  width = 6,
  height = 5)

dose

# import data & prep
load(paste0(dir_moa, "/data/final/ext/eMIR/in/mic_up/pro/MBtoPMN/", "CDplot_CUCCUGC", ".RData"))
eMIR.TS_mic.UP_pro <- TS
eMIR.DEA_mic.UP_pro <- dea

# create
obj.dose_mic.UP_pro <- CDplotWrapper.own(eMIR.DEA_mic.UP_pro, eMIR.TS_mic.UP_pro)
p.dose_mic.UP_pro_pre <- obj.dose_mic.UP_pro$p.fin + 
                          xlim(-2, 3) + 
                          theme_light(base_size=20)
p.dose_mic.UP_pro_pre

# numbers
setSize_mic.UP_pro <- get.setSizes.pro(eMIR.TS_mic.UP_pro)
p.dose_mic.UP_pro.fin <- add.numbers2base(p.dose_mic.UP_pro_pre, obj.dose_mic.UP_pro,  setSize_mic.UP_pro)
p.dose_mic.UP_pro.fin <- p.dose_mic.UP_pro.fin + labs(title = NULL)
p.dose_mic.UP_pro.fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.dose_mic.UP_pro.fin", ".svg"),
  plot = p.dose_mic.UP_pro.fin,
  device = svg,
  width = 5.5,
  height = 5)

enrichment

prep

# extract gene SYMBOLS: not sure if up and down or only down???
SYM.mic.UP_pro.1976 <- get_SYM_of_mic_targets.dn(eMIR_mic.UP_pro.MBtoPMN, 'hsa-miR-1976')
geneList_mic.UP_pro.1976 <- prep.eMIR.4enrich.pro(SYM.mic.UP_pro.1976, "MBtoPMN")
length(geneList_mic.UP_pro.1976) # 107

hsa-miR-1976

# GO
# rs.GO_mic.UP_pro.1976 <- enrichGO(gene = names(geneList_mic.UP_pro.1976),
#                 OrgDb         = org.Hs.eg.db,
#                 ont           = "ALL",
#                 pAdjustMethod = "BH",
#                 pvalueCutoff  = 0.01,
#                 readable = TRUE)
# 
# save(rs.GO_mic.UP_pro.1976,
#      file = paste0(dir_moa, "/data/final/preps/enrichment/", "rs.GO_mic.UP_pro.1976.RData"))

load(file = paste0(dir_moa,  "/data/final/preps/enrichment/", "rs.GO_mic.UP_pro.1976.RData"))


rs.GO_mic.UP_pro.1976@result %>% filter(p.adjust < 0.05) %>% nrow() # 7
clusterProfiler::dotplot(rs.GO_mic.UP_pro.1976, showCategory= 10) # 

# rs.GO_mic.UP_pro.1976@result %>% filter(p.adjust < 0.05) %>% view() # 239

rs.GO_mic.UP_pro.1976@result %>% filter(p.adjust < 0.05) %>% pull(Description)


# # dot plot
terms.pro.766.dot <- c(
"SWI/SNF superfamily-type complex",
'specific granule',
'tertiary granule',
'myeloid leukocyte activation',
'leukocyte migration',
'leukocyte chemotaxis'
)

clusterProfiler::dotplot(rs.GO_mic.UP_pro.1976, showCategory= terms.pro.766.dot) # 

p.GO.mic.UP_pro.1976 <- clusterProfiler::dotplot(rs.GO_mic.UP_pro.1976, showCategory=10, font.size = 14) #

ggsave(
  filename = paste0(dir_moa, "/figures/final/supp/mic/", "p.GO.mic.UP_pro.1976", ".svg"),
  plot = p.GO.mic.UP_pro.1976,
  device = svg,
  width = 5.5,
  height = 5)

# cNet

terms.pro.1976.cnt <- c(
  "SWI/SNF superfamily-type complex",
  "histone deacetylase complex"
)

p.cnet.mic.UP_pro.1976 <- rs.GO_mic.UP_pro.1976 %>% 
  filter(Description %in% terms.pro.1976.cnt) %>% 
  cnetplot(foldChange=geneList_mic.UP_pro.1976)

p.cnet.mic.UP_pro.1976

ggsave(
  filename = paste0(dir_moa, "/figures/final/supp/mic/", "p.cnet.mic.UP_pro.1976", ".svg"),
  plot = p.cnet.mic.UP_pro.1976,
  device = svg,
  width = 7,
  height = 7)

dnUp

mics low in mid but high in MB and S/PMN

mrn

using conserved TS8

prep

# mics
# centrality.mic.all %>% filter(traj == 'dnUp') %>% pull(ENSG) %>% write_clip() 
# recognized: conserved: 10 / all: 57
#centrality.mic.90 %>% filter(traj == 'dn') %>% pull(ENSG) %>% write_clip() 
# recognized: 19



# MBtoMM
dex.mrn_MBtoMM <- create.dex.4enrichMIR(dex.mrn, 'MBtoMM')
dex.mrn_MBtoPMN %>% filter(Gene == "S100A8") # up: wrong! should be down bc mics are inhibiting!

dex.mrn_MBvsMM <- dex.mrn_MBtoPMN

dex.mrn_MBvsMM$logFC <- dex.mrn_MBvsMM$logFC * -1
dex.mrn_MBvsMM %>% filter(Gene == "S100A8") # ok! 

# write_csv(dex.mrn_MBvsMM, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_MBvsMM", ".csv"))

# check: 2x FDR with conserved

MBvsMM

bubble

eMIR_mic.DNuP_mrn.MBvsMM <- read_csv(paste0(dir_moa, "/data/final/ext/eMIR/in/mic_dnUp/mrn/MBvsMM/", "EnrichMir_hits_siteoverlap.down_2024-05-20.csv"))

# prep
eMIR_mic.DNuP_mrn.MBvsMM <- edit_colnames_eMIR(eMIR_mic.DNuP_mrn.MBvsMM)
eMIR_mic.DNuP_mrn.MBvsMM <- add_mic.Fc('MBtoMM', eMIR_mic.DNuP_mrn.MBvsMM)

# plot
range(eMIR_mic.DNuP_mrn.MBvsMM$logFC) # -6 to 3
p.bub_mic.DNuP_mrn.MBvsMM_pre <- plot_bouble.mic.neg(eMIR_mic.DNuP_mrn.MBvsMM)

# format for export
p.bub_mic.DNuP_mrn.MBvsMM_fin <- p.bub_mic.DNuP_mrn.MBvsMM_pre + 
  theme_light(base_size=20) + 
  theme(legend.title = element_text(size = 9), legend.text = element_text(size = 9), legend.margin=margin(-5,-5,-5,-15)) + 
  labs(title = NULL)

p.bub_mic.DNuP_mrn.MBvsMM_fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.bub_mic.DNuP_mrn.MBvsMM_fin", ".svg"),
  plot = p.bub_mic.DNuP_mrn.MBvsMM_fin,
  device = svg,
  width = 6,
  height = 5)

dose

# import data & prep
load(paste0(dir_moa, "/data/final/ext/eMIR/in/mic_dnUp/mrn/MBvsMM/", "CDplot_CCCUGAG", ".RData"))
eMIR.TS_mic.DNuP_mrn <- TS
eMIR.DEA_mic.DNuP_mrn <- dea

# create
obj.dose_mic.DNuP_mrn <- CDplotWrapper.own(eMIR.DEA_mic.DNuP_mrn, eMIR.TS_mic.DNuP_mrn)
p.dose_mic.DNuP_mrn_pre <- obj.dose_mic.DNuP_mrn$p.fin + 
                          xlim(-4, 5) + 
                          theme_light(base_size=20)
p.dose_mic.DNuP_mrn_pre

# numbers
setSize_mic.DNuP_mrn <- get.setSizes.mrn(eMIR.TS_mic.DNuP_mrn)
p.dose_mic.DNuP_mrn.fin <- add.numbers2base(p.dose_mic.DNuP_mrn_pre, obj.dose_mic.DNuP_mrn,  setSize_mic.DNuP_mrn)
p.dose_mic.DNuP_mrn.fin <- p.dose_mic.DNuP_mrn.fin + labs(title = NULL)
p.dose_mic.DNuP_mrn.fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.dose_mic.DNuP_mrn.fin", ".svg"),
  plot = p.dose_mic.DNuP_mrn.fin,
  device = svg,
  width = 5.5,
  height = 5)

enrichment

prep

# extract gene SYMBOLS: not sure if up and down or only down???
SYM.mic.DNuP_mrn.125a <- get_SYM_of_mic_targets.dn(eMIR_mic.DNuP_mrn.MBvsMM, 'hsa-miR-125a-5p')
geneList_mic.DNuP_mrn.125a <- prep.eMIR.4enrich.mrn(SYM.mic.DNuP_mrn.125a, "MBtoMM")
length(geneList_mic.DNuP_mrn.125a) # 225

hsa-miR-125a

# GO
# rs.GO_mic.DNuP_mrn.125a <- enrichGO(gene = names(geneList_mic.DNuP_mrn.125a),
#                 OrgDb         = org.Hs.eg.db,
#                 ont           = "ALL",
#                 pAdjustMethod = "BH",
#                 pvalueCutoff  = 0.05,
#                 readable = TRUE)
# 
# save(rs.GO_mic.DNuP_mrn.125a,
#      file = paste0(dir_moa, "/data/final/preps/enrichment/", "rs.GO_mic.DNuP_mrn.125a.RData"))

load(file = paste0(dir_moa,  "/data/final/preps/enrichment/", "rs.GO_mic.DNuP_mrn.125a.RData"))


rs.GO_mic.DNuP_mrn.125a@result %>% filter(p.adjust < 0.05) %>% nrow() # 9
clusterProfiler::dotplot(rs.GO_mic.DNuP_mrn.125a, showCategory= 7) # 

rs.GO_mic.DNuP_mrn.125a@result %>% filter(p.adjust < 0.05) %>% pull(Description)
# rs.GO_mic.DNuP_mrn.125a@result %>% filter(p.adjust < 0.05) %>% view()

# dot plot
# terms.mrn.106.dot <- c(
#   "cell growth",
#   "positive regulation of autophagy",
#   "myeloid cell differentiation",
#   "regulation of JNK cascade",
#   "Wnt signaling pathway",
#   "histone modification"
# )

clusterProfiler::dotplot(rs.GO_mic.DNuP_mrn.125a, showCategory= 7) # 

p.GO.mic.DNuP_mrn.125a <- clusterProfiler::dotplot(rs.GO_mic.DNuP_mrn.125a, showCategory=7, font.size = 14) #

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.GO.mic.DNuP_mrn.125a", ".svg"),
  plot = p.GO.mic.DNuP_mrn.125a,
  device = svg,
  width = 5,
  height = 5)

# cNet

terms.mrn.125a.cnt <- c(
  "cell growth",
  "myeloid cell differentiation"
)

p.cnet.mic.DNuP_mrn.125a <- rs.GO_mic.DNuP_mrn.125a %>% 
  filter(Description %in% terms.mrn.125a.cnt) %>% 
  cnetplot(foldChange=geneList_mic.DNuP_mrn.125a)

p.cnet.mic.DNuP_mrn.125a

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.cnet.mic.DNuP_mrn.125a", ".svg"),
  plot = p.cnet.mic.DNuP_mrn.125a,
  device = svg,
  width = 7,
  height = 7)

pro

using all TS8

prep

table(centrality.mic.all$traj)
# mics
# centrality.mic.all %>% filter(traj == 'dnUp') %>% pull(ENSG) %>% write_clip() 
# recognized: conserved: 19 
#centrality.mic.90 %>% filter(traj == 'dn') %>% pull(ENSG) %>% write_clip() 
# recognized: 19


# MBvsMC: might be better bc more


# MBtoPMN: 
dex.pro_MBtoMM <- create.dex.4enrichMIR(dex.pro, 'MBtoMM')
dex.pro_MBtoMM %>% filter(Gene == "S100A8") # up: wrong! should be down bc mics are inhibiting!

dex.pro_MBvsMM <- dex.pro_MBtoMM
dex.pro_MBvsMM$logFC <- dex.pro_MBvsMM$logFC * -1
dex.pro_MBvsMM %>% filter(Gene == "S100A8") # ok! 

# write_csv(dex.pro_MBvsMM, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.pro_MBvsPMN", ".csv"))

# check: 0


# MBtoPMN: 
dex.pro_MBtoB <- create.dex.4enrichMIR(dex.pro, 'MBtoB')
dex.pro_MBtoB %>% filter(Gene == "S100A8") # up: wrong! should be down bc mics are inhibiting!

dex.pro_MBvsB <- dex.pro_MBtoB
dex.pro_MBvsB$logFC <- dex.pro_MBvsB$logFC * -1
dex.pro_MBvsB %>% filter(Gene == "S100A8") # ok! 

# write_csv(dex.pro_MBvsB, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.pro_MBvsB", ".csv"))

# check: good

MBvsB

bubble

eMIR_mic.DNuP_pro.MBvsB <- read_csv(paste0(dir_moa, "/data/final/ext/eMIR/in/mic_dnUp/pro/MBvsB/", "EnrichMir_hits_siteoverlap.down_2024-05-21.csv"))

# prep
eMIR_mic.DNuP_pro.MBvsB <- edit_colnames_eMIR(eMIR_mic.DNuP_pro.MBvsB)
eMIR_mic.DNuP_pro.MBvsB <- add_mic.Fc('MBvsB', eMIR_mic.DNuP_pro.MBvsB)

# plot
range(eMIR_mic.DNuP_pro.MBvsB$logFC) # -6 to 3
p.bub_mic.DNuP_pro.MBvsB_pre <- plot_bouble.mic.neg(eMIR_mic.DNuP_pro.MBvsB)

# format for export
p.bub_mic.DNuP_pro.MBvsB_fin <- p.bub_mic.DNuP_pro.MBvsB_pre + 
  theme_light(base_size=20) + 
  theme(legend.title = element_text(size = 9), legend.text = element_text(size = 9), legend.margin=margin(-5,-5,-5,-15)) + 
  labs(title = NULL)

p.bub_mic.DNuP_pro.MBvsB_fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/supp/mic/", "p.bub_mic.DNuP_pro.MBvsB_fin", ".svg"),
  plot = p.bub_mic.DNuP_pro.MBvsB_fin,
  device = svg,
  width = 6,
  height = 5)

dose

# import data & prep
load(paste0(dir_moa, "/data/final/ext/eMIR/in/mic_dnUp/pro/MBvsB/", "CDplot_CCCUGAG", ".RData"))
eMIR.TS_mic.DNuP_pro <- TS
eMIR.DEA_mic.DNuP_pro <- dea

# create
obj.dose_mic.DNuP_pro <- CDplotWrapper.own(eMIR.DEA_mic.DNuP_pro, eMIR.TS_mic.DNuP_pro)
p.dose_mic.DNuP_pro_pre <- obj.dose_mic.DNuP_pro$p.fin + 
                          xlim(-3, 3) + 
                          theme_light(base_size=20)
p.dose_mic.DNuP_pro_pre

# numbers
setSize_mic.DNuP_pro <- get.setSizes.pro(eMIR.TS_mic.DNuP_pro)
p.dose_mic.DNuP_pro.fin <- add.numbers2base(p.dose_mic.DNuP_pro_pre, obj.dose_mic.DNuP_pro,  setSize_mic.DNuP_pro)
p.dose_mic.DNuP_pro.fin <- p.dose_mic.DNuP_pro.fin + labs(title = NULL)
p.dose_mic.DNuP_pro.fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.dose_mic.DNuP_pro.fin", ".svg"),
  plot = p.dose_mic.DNuP_pro.fin,
  device = svg,
  width = 5.5,
  height = 5)

enrichment

prep

# extract gene SYMBOLS: not sure if up and down or only down???
SYM.mic.DNuP_pro.125a <- get_SYM_of_mic_targets.dn(eMIR_mic.DNuP_pro.MBvsB, 'hsa-miR-125a-5p')
geneList_mic.DNuP_pro.125a <- prep.eMIR.4enrich.pro(SYM.mic.DNuP_pro.125a, "MBtoB")
length(geneList_mic.DNuP_pro.125a) # 40

hsa-miR-125a-5p

# GO
# rs.GO_mic.DNuP_pro.125a <- enrichGO(gene = names(geneList_mic.DNuP_pro.125a),
#                 OrgDb         = org.Hs.eg.db,
#                 ont           = "ALL",
#                 pAdjustMethod = "BH",
#                 pvalueCutoff  = 0.05,
#                 readable = TRUE)
# 
# save(rs.GO_mic.DNuP_pro.125a,
#      file = paste0(dir_moa, "/data/final/preps/enrichment/", "rs.GO_mic.DNuP_pro.125a.RData"))

load(file = paste0(dir_moa,  "/data/final/preps/enrichment/", "rs.GO_mic.DNuP_pro.125a.RData"))


rs.GO_mic.DNuP_pro.125a@result %>% filter(p.adjust < 0.05) %>% nrow() # 
clusterProfiler::dotplot(rs.GO_mic.DNuP_pro.125a, showCategory= 10) # 

# rs.GO_mic.DNuP_pro.125a@result %>% filter(p.adjust < 0.05) %>% view() # 239

rs.GO_mic.DNuP_pro.125a@result %>% filter(p.adjust < 0.05) %>% pull(Description)


# dot plot
terms.pro.125a.dot <- c(
'cell cortex',
'vesicle organization',
'nuclear membrane organization',
'myeloid cell differentiation',
'tertiary granule',
'growth hormone receptor signaling pathway via JAK-STAT'
)

clusterProfiler::dotplot(rs.GO_mic.DNuP_pro.125a, showCategory= terms.pro.125a.dot) # 

p.GO.mic.DNuP_pro.125a <- clusterProfiler::dotplot(rs.GO_mic.DNuP_pro.125a, showCategory=terms.pro.125a.dot, font.size = 14) #

ggsave(
  filename = paste0(dir_moa, "/figures/final/supp/mic/", "p.GO.mic.DNuP_pro.125a", ".svg"),
  plot = p.GO.mic.DNuP_pro.125a,
  device = svg,
  width = 5.5,
  height = 5)

# cNet

terms.pro.125a.cnt <- c(
'myeloid cell differentiation',
'tertiary granule'
)

p.cnet.mic.DNuP_pro.125a <- rs.GO_mic.DNuP_pro.125a %>% 
  filter(Description %in% terms.pro.125a.cnt) %>% 
  cnetplot(foldChange=geneList_mic.DNuP_pro.125a)

p.cnet.mic.DNuP_pro.125a

ggsave(
  filename = paste0(dir_moa, "/figures/final/supp/mic/", "p.cnet.mic.DNuP_pro.125a", ".svg"),
  plot = p.cnet.mic.DNuP_pro.125a,
  device = svg,
  width = 7,
  height = 7)

upDn

UpDn = temporarily repressing transcripts

mrn

using conserved TS8

prep

# mics
# centrality.mic.all %>% filter(traj == 'upDn') %>% pull(ENSG) %>% write_clip() 
# recognized: conserved: 12 = used!

# MBvsB! = final
# 1 FDR


# MBtoMM
# check: 1x pval

# MBvsMM
# 

# MCtoMM or B
dex.mrn_MCtoMM <- create.dex.4enrichMIR(dex.mrn, 'MCtoMM')
dex.mrn_MCtoMM %>% filter(Gene == "S100A8") # up: ok
# write_csv(dex.mrn_MCtoMM, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_MCtoMM", ".csv"))
# no FDR

# MCtoB
dex.mrn_MCtoB <- create.dex.4enrichMIR(dex.mrn, 'MCtoB')
dex.mrn_MCtoB %>% filter(Gene == "S100A8") # up: ok
# write_csv(dex.mrn_MCtoB, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_MCtoB", ".csv"))
# no FDR

# MBtoMC
# nothing

MBvsB

bubble

eMIR_mic.UpDn_mrn.MBvsB <- read_csv(paste0(dir_moa, "/data/final/ext/eMIR/in/mic_upDn/mrn/MBvsB/", "EnrichMir_hits_siteoverlap.down_2024-05-21.csv"))

# prep
eMIR_mic.UpDn_mrn.MBvsB <- edit_colnames_eMIR(eMIR_mic.UpDn_mrn.MBvsB)
eMIR_mic.UpDn_mrn.MBvsB <- add_mic.Fc('MBtoMC', eMIR_mic.UpDn_mrn.MBvsB)

# plot
range(eMIR_mic.UpDn_mrn.MBvsB$logFC) # -1 to 5
p.bub_mic.UpDn_mrn.MBvsB_pre <- plot_bouble.mic.negPos(eMIR_mic.UpDn_mrn.MBvsB)

# format for export
p.bub_mic.UpDn_mrn.MBvsB_fin <- p.bub_mic.UpDn_mrn.MBvsB_pre + 
  theme_light(base_size=20) + 
  theme(legend.title = element_text(size = 9), legend.text = element_text(size = 9), legend.margin=margin(-5,-5,-5,-15)) + 
  labs(title = NULL)

p.bub_mic.UpDn_mrn.MBvsB_fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.bub_mic.UpDn_mrn.MBvsB_fin", ".svg"),
  plot = p.bub_mic.UpDn_mrn.MBvsB_fin,
  device = svg,
  width = 6,
  height = 5)

dose

# import data & prep
load(paste0(dir_moa, "/data/final/ext/eMIR/in/mic_upDn/mrn/MBvsB/", "CDplot_AGUGCAA", ".RData"))
eMIR.TS_mic.UpDn_mrn <- TS
eMIR.DEA_mic.UpDn_mrn <- dea

# create
obj.dose_mic.UpDn_mrn <- CDplotWrapper.own(eMIR.DEA_mic.UpDn_mrn, eMIR.TS_mic.UpDn_mrn)
p.dose_mic.UpDn_mrn_pre <- obj.dose_mic.UpDn_mrn$p.fin + 
                          xlim(-4, 4) + 
                          theme_light(base_size=20)
p.dose_mic.UpDn_mrn_pre

# numbers
setSize_mic.UpDn_mrn <- get.setSizes.mrn(eMIR.TS_mic.UpDn_mrn)
p.dose_mic.UpDn_mrn.fin <- add.numbers2base(p.dose_mic.UpDn_mrn_pre, obj.dose_mic.UpDn_mrn,  setSize_mic.UpDn_mrn)
p.dose_mic.UpDn_mrn.fin <- p.dose_mic.UpDn_mrn.fin + labs(title = NULL)
p.dose_mic.UpDn_mrn.fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.dose_mic.UpDn_mrn.fin", ".svg"),
  plot = p.dose_mic.UpDn_mrn.fin,
  device = svg,
  width = 5.5,
  height = 5)

enrichment

prep

# extract gene SYMBOLS: not sure if up and down or only down???
SYM.mic.UpDn_mrn.454 <- get_SYM_of_mic_targets.dn(eMIR_mic.UpDn_mrn.MBvsB, 'hsa-miR-454-3p')
geneList_mic.UpDn_mrn.454 <- prep.eMIR.4enrich.mrn(SYM.mic.UpDn_mrn.454, "MBtoB")
length(geneList_mic.UpDn_mrn.454) # 257

hsa-miR-454

# GO
# rs.GO_mic.UpDn_mrn.454 <- enrichGO(gene = names(geneList_mic.UpDn_mrn.454),
#                 OrgDb         = org.Hs.eg.db,
#                 ont           = "ALL",
#                 pAdjustMethod = "BH",
#                 pvalueCutoff  = 0.05,
#                 readable = TRUE,
#                 universe = universe_mrn.ENTREZ)

# save(rs.GO_mic.UpDn_mrn.454,
#      file = paste0(dir_moa, "/data/final/preps/enrichment/", "rs.GO_mic.UpDn_mrn.454.RData"))

load(file = paste0(dir_moa,  "/data/final/preps/enrichment/", "rs.GO_mic.UpDn_mrn.454.RData"))


rs.GO_mic.UpDn_mrn.454@result %>% filter(p.adjust < 0.05) %>% nrow() # 4

clusterProfiler::dotplot(rs.GO_mic.UpDn_mrn.454, showCategory= 9) # redundant

clusterProfiler::dotplot(rs.GO_mic.UpDn_mrn.454, showCategory= 5) # better

# dot plot
p.GO.mic.UpDn_mrn.454 <- clusterProfiler::dotplot(rs.GO_mic.UpDn_mrn.454, showCategory=5, font.size = 14) #

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.GO.mic.UpDn_mrn.454", ".svg"),
  plot = p.GO.mic.UpDn_mrn.454,
  device = svg,
  width = 5.5,
  height = 5)


# cNet
rs.GO_mic.UpDn_mrn.454@result %>% filter(p.adjust < 0.05) %>% pull(Description)

terms.mrn454.cnt <- c(
  "negative regulation of transcription by RNA polymerase II"
)

p.cnet.mic.UpDn_mrn.454 <- rs.GO_mic.UpDn_mrn.454 %>% 
  filter(Description %in% terms.mrn454.cnt) %>% 
  cnetplot(foldChange=geneList_mic.UpDn_mrn.454)

p.cnet.mic.UpDn_mrn.454

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.cnet.mic.UpDn_mrn.454", ".svg"),
  plot = p.cnet.mic.UpDn_mrn.454,
  device = svg,
  width = 7,
  height = 7)

pro

prep

# mics
# centrality.mic.all %>% filter(traj == 'upDn') %>% pull(ENSG) %>% write_clip() 
# recognized: conserved: 19 

# MBtoB: 
dex.pro_MBtoB <- create.dex.4enrichMIR(dex.pro, 'MBtoB')
dex.pro_MBvsB <- dex.pro_MBtoB
dex.pro_MBvsB$logFC <- dex.pro_MBvsB$logFC * -1
dex.pro_MBvsB %>% filter(Gene == "S100A8") # ok! 

# write_csv(dex.pro_MBvsB, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.pro_MBvsB", ".csv"))
# check: 2x pval with conserved TS8!

MBvsB

bubble

eMIR_mic.upDn_pro.MBvsB <- read_csv(paste0(dir_moa, "/data/final/ext/eMIR/in/mic_upDn/pro/MBvsB/", "EnrichMir_hits_siteoverlap.down_2024-05-21.csv"))

# prep
eMIR_mic.upDn_pro.MBvsB <- edit_colnames_eMIR(eMIR_mic.upDn_pro.MBvsB)
eMIR_mic.upDn_pro.MBvsB <- add_mic.Fc('MBtoB', eMIR_mic.upDn_pro.MBvsB)

# plot
range(eMIR_mic.upDn_pro.MBvsB$logFC) # -6 to 3
p.bub_mic.upDn_pro.MBvsB_pre <- plot_bouble.mic.negPos(eMIR_mic.upDn_pro.MBvsB)

# format for export
p.bub_mic.upDn_pro.MBvsB_fin <- p.bub_mic.upDn_pro.MBvsB_pre + 
  theme_light(base_size=20) + 
  theme(legend.title = element_text(size = 9), legend.text = element_text(size = 9), legend.margin=margin(-5,-5,-5,-15)) + 
  labs(title = NULL)

p.bub_mic.upDn_pro.MBvsB_fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/supp/mic/", "p.bub_mic.upDn_pro.MBvsB_fin", ".svg"),
  plot = p.bub_mic.upDn_pro.MBvsB_fin,
  device = svg,
  width = 6,
  height = 5)

dose

# import data & prep
load(paste0(dir_moa, "/data/final/ext/eMIR/in/mic_upDn/pro/MBvsB/", "CDplot_GAGAUGA", ".RData"))
eMIR.TS_mic.upDn_pro <- TS
eMIR.DEA_mic.upDn_pro <- dea

# create
obj.dose_mic.upDn_pro <- CDplotWrapper.own(eMIR.DEA_mic.upDn_pro, eMIR.TS_mic.upDn_pro)
p.dose_mic.upDn_pro_pre <- obj.dose_mic.upDn_pro$p.fin + 
                          xlim(-2, 3) + 
                          theme_light(base_size=20)
p.dose_mic.upDn_pro_pre

# numbers
setSize_mic.upDn_pro <- get.setSizes.pro(eMIR.TS_mic.upDn_pro)
p.dose_mic.upDn_pro.fin <- add.numbers2base(p.dose_mic.upDn_pro_pre, obj.dose_mic.upDn_pro,  setSize_mic.upDn_pro)
p.dose_mic.upDn_pro.fin <- p.dose_mic.upDn_pro.fin + labs(title = NULL)
p.dose_mic.upDn_pro.fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/supp/mic/", "p.dose_mic.upDn_pro.fin", ".svg"),
  plot = p.dose_mic.upDn_pro.fin,
  device = svg,
  width = 5.5,
  height = 5)

enrichment

prep

# extract gene SYMBOLS: not sure if up and down or only down???
SYM.mic.upDn_pro.143 <- get_SYM_of_mic_targets.dn(eMIR_mic.upDn_pro.MBvsB, 'hsa-miR-143-3p')
geneList_mic.upDn_pro.143 <- prep.eMIR.4enrich.pro(SYM.mic.upDn_pro.143, "MBtoB")
length(geneList_mic.upDn_pro.143) # 33

hsa-miR-143

# # GO
# rs.GO_mic.upDn_pro.143 <- enrichGO(gene = names(geneList_mic.upDn_pro.143),
#                 OrgDb         = org.Hs.eg.db,
#                 ont           = "ALL",
#                 pAdjustMethod = "BH",
#                 pvalueCutoff  = 0.01,
#                 readable = TRUE)
# 
# save(rs.GO_mic.upDn_pro.143,
#      file = paste0(dir_moa, "/data/final/preps/enrichment/", "rs.GO_mic.upDn_pro.143.RData"))

load(file = paste0(dir_moa,  "/data/final/preps/enrichment/", "rs.GO_mic.upDn_pro.143.RData"))


rs.GO_mic.upDn_pro.143@result %>% filter(p.adjust < 0.05) %>% nrow() # 4
clusterProfiler::dotplot(rs.GO_mic.upDn_pro.143, showCategory= 10) # 

# rs.GO_mic.upDn_pro.143@result %>% filter(p.adjust < 0.05) %>% view() # 239

rs.GO_mic.upDn_pro.143@result %>% filter(p.adjust < 0.05) %>% pull(Description)


# # dot plot
# terms.pro.766.dot <- c(
# "SWI/SNF superfamily-type complex"
# 'specific granule',
# 'tertiary granule',
# 'myeloid leukocyte activation',
# 'leukocyte migration',
# 'leukocyte chemotaxis'
# )

clusterProfiler::dotplot(rs.GO_mic.upDn_pro.143, showCategory= 4) # 

p.GO.mic.upDn_pro.143 <- clusterProfiler::dotplot(rs.GO_mic.upDn_pro.143, showCategory=4, font.size = 14) #

ggsave(
  filename = paste0(dir_moa, "/figures/final/supp/mic/", "p.GO.mic.upDn_pro.143", ".svg"),
  plot = p.GO.mic.upDn_pro.143,
  device = svg,
  width = 5.5,
  height = 5)

# cNet

terms.pro.143.cnt <- c(
  "focal adhesion",
  "cell-substrate junction"
)

p.cnet.mic.upDn_pro.143 <- rs.GO_mic.upDn_pro.143 %>% 
  filter(Description %in% terms.pro.143.cnt) %>% 
  cnetplot(foldChange=geneList_mic.upDn_pro.143)

p.cnet.mic.upDn_pro.143

ggsave(
  filename = paste0(dir_moa, "/figures/final/supp/mic/", "p.cnet.mic.upDn_pro.143", ".svg"),
  plot = p.cnet.mic.upDn_pro.143,
  device = svg,
  width = 7,
  height = 7)

upDn late

UpDn = temporarily repressing transcripts

mrn

using conserved TS8

prep

table(centrality.mic.all$traj)
# mics
# centrality.mic.all %>% filter(traj == 'upDn late') %>% pull(ENSG) %>% write_clip() 
# recognized: 
# conserved: 2
# all = 10. USED

# MMtoB = final!
# 1 FDR, curve is meh

# MCtoPMN
dex.mrn_MCtoPMN <- create.dex.4enrichMIR(dex.mrn, 'MCtoPMN')
dex.mrn_MCtoPMN %>% filter(Gene == "S100A8") # up: ok
# write_csv(dex.mrn_MCtoPMN, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_MCtoPMN", ".csv"))
# nothing

dex.mrn_MCvsPMN <- dex.mrn_MCtoPMN
dex.mrn_MCvsPMN$logFC <- dex.mrn_MCvsPMN$logFC * -1
# write_csv(dex.mrn_MCvsPMN, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_MCvsPMN", ".csv"))
# 1x pval


# BtoPMN
dex.mrn_BtoPMN <- create.dex.4enrichMIR(dex.mrn, 'BtoPMN')
dex.mrn_BtoPMN %>% filter(Gene == "S100A8") # up: ok
# write_csv(dex.mrn_BtoPMN, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_BtoPMN", ".csv"))
# nothing

dex.mrn_BvsPMN <- dex.mrn_BtoPMN
dex.mrn_BvsPMN$logFC <- dex.mrn_BvsPMN$logFC * -1
# write_csv(dex.mrn_BvsPMN, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_BvsPMN", ".csv"))
# 1x pval


# MMtoPMN: 1 FDR, curve not great
# MMvsPMN: nothing

# MMtoS
dex.mrn_MMtoS <- create.dex.4enrichMIR(dex.mrn, 'MMtoS')
dex.mrn_MMtoS %>% filter(Gene == "S100A8") # up: ok
# write_csv(dex.mrn_MMtoS, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_MMtoS", ".csv"))
# nothing

dex.mrn_MMvsS <- dex.mrn_MMtoS
dex.mrn_MMvsS$logFC <- dex.mrn_MMvsS$logFC * -1
# write_csv(dex.mrn_MMvsS, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_MMvsS", ".csv"))
# 1x pval


# MCtoB
dex.mrn_MCtoS <- create.dex.4enrichMIR(dex.mrn, 'MCtoS')
dex.mrn_MCtoS %>% filter(Gene == "S100A8") # up: ok
# write_csv(dex.mrn_MCtoS, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_MCtoS", ".csv"))
# nothing


# MMtoS
dex.mrn_MMtoB <- create.dex.4enrichMIR(dex.mrn, 'MMtoB')
# write_csv(dex.mrn_MMtoB, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_MMtoB", ".csv"))
# 1x pval, bad curve

dex.mrn_MMvsB <- dex.mrn_MMtoB
dex.mrn_MMvsB$logFC <- dex.mrn_MMvsB$logFC * -1
# write_csv(dex.mrn_MMvsB, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_MMvsB", ".csv"))
# 1x pval


# MCtoB
dex.mrn_MCtoB <- create.dex.4enrichMIR(dex.mrn, 'MCtoB')
dex.mrn_MCtoB %>% filter(Gene == "S100A8") # up: ok
# write_csv(dex.mrn_MCtoB, paste0(dir_moa, "/data/final/ext/eMIR/out/", "dex.mrn_MCtoB", ".csv"))
# no FDR

# MBtoMC
# nothing

MMtoB

bubble

eMIR_mic.UpDnL_mrn.MMtoB <- read_csv(paste0(dir_moa, "/data/final/ext/eMIR/in/mic_upDnL/mrn/MMtoB/", "EnrichMir_hits_siteoverlap.down_2024-05-21.csv"))

# prep
eMIR_mic.UpDnL_mrn.MMtoB <- edit_colnames_eMIR(eMIR_mic.UpDnL_mrn.MMtoB)
eMIR_mic.UpDnL_mrn.MMtoB <- add_mic.Fc('MMtoB', eMIR_mic.UpDnL_mrn.MMtoB)

# plot
range(eMIR_mic.UpDnL_mrn.MMtoB$logFC) # -1 to 5
p.bub_mic.UpDnL_mrn.MMtoB_pre <- plot_bouble.mic.negPos(eMIR_mic.UpDnL_mrn.MMtoB)

# format for export
p.bub_mic.UpDnL_mrn.MMtoB_fin <- p.bub_mic.UpDnL_mrn.MMtoB_pre + 
  theme_light(base_size=20) + 
  theme(legend.title = element_text(size = 9), legend.text = element_text(size = 9), legend.margin=margin(-5,-5,-5,-15)) + 
  labs(title = NULL)

p.bub_mic.UpDnL_mrn.MMtoB_fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.bub_mic.UpDnL_mrn.MMtoB_fin", ".svg"),
  plot = p.bub_mic.UpDnL_mrn.MMtoB_fin,
  device = svg,
  width = 6,
  height = 5)

dose

# import data & prep
load(paste0(dir_moa, "/data/final/ext/eMIR/in/mic_upDnL/mrn/MMtoB/", "CDplot_UCCUGGG", ".RData"))
eMIR.TS_mic.UpDnL_mrn <- TS
eMIR.DEA_mic.UpDnL_mrn <- dea

# create
obj.dose_mic.UpDnL_mrn <- CDplotWrapper.own(eMIR.DEA_mic.UpDnL_mrn, eMIR.TS_mic.UpDnL_mrn)
p.dose_mic.UpDnL_mrn_pre <- obj.dose_mic.UpDnL_mrn$p.fin + 
                          xlim(-2, 1.5) + 
                          theme_light(base_size=20)
p.dose_mic.UpDnL_mrn_pre

# numbers
setSize_mic.UpDnL_mrn <- get.setSizes.mrn(eMIR.TS_mic.UpDnL_mrn)
p.dose_mic.UpDnL_mrn.fin <- add.numbers2base(p.dose_mic.UpDnL_mrn_pre, obj.dose_mic.UpDnL_mrn,  setSize_mic.UpDnL_mrn)
p.dose_mic.UpDnL_mrn.fin <- p.dose_mic.UpDnL_mrn.fin + labs(title = NULL)
p.dose_mic.UpDnL_mrn.fin

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.dose_mic.UpDnL_mrn.fin", ".svg"),
  plot = p.dose_mic.UpDnL_mrn.fin,
  device = svg,
  width = 5.5,
  height = 5)

enrichment

prep

# extract gene SYMBOLS: not sure if up and down or only down???
SYM.mic.UpDnL_mrn.1343 <- get_SYM_of_mic_targets.dn(eMIR_mic.UpDnL_mrn.MMtoB, 'hsa-miR-1343-3p')
geneList_mic.UpDnL_mrn.1343 <- prep.eMIR.4enrich.mrn(SYM.mic.UpDnL_mrn.1343, "MBtoB")
length(geneList_mic.UpDnL_mrn.1343) # 257

hsa-miR-1343

# GO
# rs.GO_mic.UpDnL_mrn.1343 <- enrichGO(gene = names(geneList_mic.UpDnL_mrn.1343),
#                 OrgDb         = org.Hs.eg.db,
#                 ont           = "ALL",
#                 pAdjustMethod = "BH",
#                 pvalueCutoff  = 0.05,
#                 readable = TRUE,
#                 universe = universe_mrn.ENTREZ)
# 
# save(rs.GO_mic.UpDnL_mrn.1343,
#      file = paste0(dir_moa, "/data/final/preps/enrichment/", "rs.GO_mic.UpDnL_mrn.1343.RData"))

load(file = paste0(dir_moa,  "/data/final/preps/enrichment/", "rs.GO_mic.UpDnL_mrn.1343.RData"))

rs.GO_mic.UpDnL_mrn.1343@result %>% filter(p.adjust < 0.05) %>% nrow() # 22

clusterProfiler::dotplot(rs.GO_mic.UpDnL_mrn.1343, showCategory= 9) # redundant

clusterProfiler::dotplot(rs.GO_mic.UpDnL_mrn.1343, showCategory= 8) # better

# rs.GO_mic.UpDnL_mrn.1343@result %>% filter(p.adjust < 0.05) %>% view()

# dot plot
terms.mrn1343.dot <- c(
  'translation',
  'ribosome',
  'mitochondrial inner membrane'
)


p.GO.mic.UpDnL_mrn.1343 <- clusterProfiler::dotplot(rs.GO_mic.UpDnL_mrn.1343, showCategory=terms.mrn1343.dot, font.size = 14) #

p.GO.mic.UpDnL_mrn.1343

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.GO.mic.UpDnL_mrn.1343", ".svg"),
  plot = p.GO.mic.UpDnL_mrn.1343,
  device = svg,
  width = 5.5,
  height = 5)


# cNet
rs.GO_mic.UpDnL_mrn.1343@result %>% filter(p.adjust < 0.05) %>% pull(Description)

terms.mrn1343.cnt <- c(
  'translation'
)

p.cnet.mic.UpDnL_mrn.1343 <- rs.GO_mic.UpDnL_mrn.1343 %>% 
  filter(Description %in% terms.mrn1343.cnt) %>% 
  cnetplot(foldChange=geneList_mic.UpDnL_mrn.1343)

p.cnet.mic.UpDnL_mrn.1343

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_4/", "p.cnet.mic.UpDnL_mrn.1343", ".svg"),
  plot = p.cnet.mic.UpDnL_mrn.1343,
  device = svg,
  width = 7,
  height = 7)