transcription factors

Last update: 24-10-05

ToDo: - format

# version: 22.03.24
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 2024.01.09
# 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 packagea
sapply(packages_general, require, character.only = TRUE)
sapply(packages_enrichment, require, character.only = TRUE)

library("viridis")


## graphics

# size
global_size = 21

thema_nix = theme(axis.title.x=element_blank(),
                  axis.title.y=element_blank())
thema_noY = theme(axis.title.y=element_blank())
thema_noX = theme(axis.title.x=element_blank())

# shapes (mrn & pro)
shps <- c(1, 16)

base data

medians

# extract data
fan <- gnp.mrn$fan

# pro: single values
exp_proL <- as.data.frame(t(gnp.pro$exp_pro_vsn))
exp_proL$UNIPROT <- rownames(exp_proL)
exp_proL <- exp_proL %>% pivot_longer(!UNIPROT, names_to = "id", values_to = "value")
exp_proL$stage <- gnp.pro$san_pro$stage[match(exp_proL$id, gnp.pro$san_pro$id)]
exp_proL$SYMBOL <- fan$SYMBOL[match(exp_proL$UNIPROT, fan$UNIPROT)]
exp_proL$ENSG <- fan$ENSG[match(exp_proL$UNIPROT, fan$UNIPROT)]

# pro: means
exp_proLc <- exp_proL %>%
  select(ENSG, stage, value) %>%
  group_by(ENSG, stage) %>%
  dplyr::summarise(mean= mean(value), median= median(value))
exp_proLc$SYMBOL <- fan$SYMBOL[match(exp_proLc$ENSG, fan$ENSG)]
exp_proLc$stage <- factor(exp_proLc$stage, levels = c('MB', 'PM', 'MC', 'MM', 'B', 'PMN'))

exp_proWc <- exp_proLc %>% select(SYMBOL, stage, median) %>%
  pivot_wider(names_from = stage, values_from = median)



# mrn: single values
exp_mrnL <- as.data.frame(t(gnp.mrn$exp_mrn_tpm_vsn))
exp_mrnL$ENSG <- rownames(exp_mrnL)
exp_mrnL <- exp_mrnL %>% pivot_longer(!ENSG, names_to = "id", values_to = "value")
exp_mrnL$stage <- gnp.mrn$san_mrn$stage[match(exp_mrnL$id, gnp.mrn$san_mrn$id)]
exp_mrnL$SYMBOL <- fan$SYMBOL[match(exp_mrnL$ENSG, fan$ENSG)]

# mrn: means
exp_mrnLc <- exp_mrnL %>%
  select(ENSG, stage, value) %>%
  group_by(ENSG, stage) %>%
  dplyr::summarise(median= median(value, na.rm = T))
exp_mrnLc$SYMBOL <- fan$SYMBOL[match(exp_mrnLc$ENSG, fan$ENSG)]
exp_mrnLc$stage <- factor(exp_mrnLc$stage, levels = c('MB', 'PM', 'MC', 'MM', 'B', 'S', 'PMN'))

# turn into wide to get direction 
exp_mrnWc <- exp_mrnLc %>% select(SYMBOL, stage, median) %>%
  pivot_wider(names_from = stage, values_from = median)



# mic: means

# single values
exp_micL <- as.data.frame(t(gnp.mic$exp_mic_vsn))
exp_micL$SYMBOL <- rownames(exp_micL)
exp_micL <- exp_micL %>% pivot_longer(!SYMBOL, names_to = "id", values_to = "value")
exp_micL$stage <- gnp.mic$san_mic$stage[match(exp_micL$id, gnp.mic$san_mic$id)]

exp_micLc <- exp_micL %>%
  select(SYMBOL, stage, value) %>%
  group_by(SYMBOL, stage) %>%
  dplyr::summarise(median= median(value, na.rm = T))
exp_micLc$stage <- factor(exp_micLc$stage, levels = c('MB', 'MC', 'MM', 'B', 'S', 'PMN'))

# turn into wide to get direction 
exp_micWc <- exp_micLc %>% select(SYMBOL, stage, median) %>%
  pivot_wider(names_from = stage, values_from = median)


#exp_micWc$direction.MBtoMC <- ifelse(exp_micWc$MB > exp_micWc$MC, 'down', 'up')


#exp_micWc$direction <- ifelse(exp_micWc$MB > exp_micWc$PMN, 'down', 'up')

# transfer
#exp_micLc$direction.MBtoMC <- exp_micWc$direction.MBtoMC[match(exp_micLc$SYMBOL, exp_micWc$SYMBOL)]

universes

# extract data
fan <- gnp.mrn$fan

# ENTREZ
dex.all$dex.mrn$dex_stage_mrn_vsn$ENTREZ <- fan$ENTREZ[match(dex.all$dex.mrn$dex_stage_mrn_vsn$fid, fan$ENSG)]
dex.all$dex.pro$dex_stage_pro_vsn$ENTREZ <- fan$ENTREZ[match(dex.all$dex.pro$dex_stage_pro_vsn$fid, fan$UNIPROT)]

# ENSG
dex.all$dex.mrn$dex_stage_mrn_vsn$ENSG <- dex.all$dex.mrn$dex_stage_mrn_vsn$fid
dex.all$dex.pro$dex_stage_pro_vsn$ENSG <- fan$ENSG[match(dex.all$dex.pro$dex_stage_pro_vsn$fid, fan$UNIPROT)]

# SYMBOL
dex.all$dex.mrn$dex_stage_mrn_vsn$SYMBOL <- fan$SYMBOL[match(dex.all$dex.mrn$dex_stage_mrn_vsn$fid, fan$ENSG)]
dex.all$dex.pro$dex_stage_pro_vsn$SYMBOL <- fan$SYMBOL[match(dex.all$dex.pro$dex_stage_pro_vsn$fid, fan$UNIPROT)]


#### backgrounds

universe_mrn.ENSG <- unique(dex.all$dex.mrn$dex_stage_mrn_vsn$ENSG)
universe_pro.ENSG <- unique(dex.all$dex.pro$dex_stage_pro_vsn$ENSG)
universe_combo.ENSG <- unique(c(universe_mrn.ENSG, universe_pro.ENSG))

universe_mrn.ENTREZ <- unique(dex.all$dex.mrn$dex_stage_mrn_vsn$ENTREZ)
universe_pro.ENTREZ <- unique(dex.all$dex.pro$dex_stage_pro_vsn$ENTREZ)
universe_combo.ENTREZ <- unique(c(universe_mrn.ENTREZ, universe_pro.ENTREZ))

universe_mrn.SYMBOL <- fan$SYMBOL[match(universe_mrn.ENSG, fan$ENSG)]
universe_pro.SYMBOL <- fan$SYMBOL[match(universe_pro.ENSG, fan$ENSG)]

# remove garbage
universe_mrn.SYMBOL <- universe_mrn.SYMBOL[-c(2612)]
# write_clip(universe_mrn.SYMBOL)

ChEA3

fct

library(httr)
library(jsonlite)

recieve_ChEA3 <- function(vec_SYM, name){
  # fixed variables
  url <-  "https://maayanlab.cloud/chea3/api/enrich/"
  encode <- "json"
  # data depended var
  payload <- list(query_name = name, gene_set = vec_SYM)
  response <- POST(url = url, body = payload, encode = encode)
  json <- content(response, "text")
  #results as list of R dataframes
  results <- fromJSON(json)
  return(results)
}


# dex results: Fc filtered as SYMBOLS vectors
get.vecSYM_dex.dir.FC <- function(dex, COMPoi, DIRECTION, FClevel){
  dex_comp <- dex %>% filter(comparison %in% COMPoi)
  if(DIRECTION == 'up'){
    foi <- dex_comp %>% filter(adj.P.Val < 0.05 & logFC >= FClevel) %>% arrange(desc(logFC)) %>% pull(SYMBOL)
  }
  if(DIRECTION == 'down'){
    foi <- dex_comp %>% filter(adj.P.Val < 0.05 & logFC <= -FClevel) %>% arrange(logFC) %>% pull(SYMBOL)
  }
  #fid <- fid[!is.na(fid)]
  return(foi)
}

prep feature lists

# upregulated transcripts as SYMBOLS

l.dexNfc1_SYM_gnp.mrn_UP <- list(
  'MB->PM' = get.vecSYM_dex.dir.FC(dex.all$dex.mrn$dex_stage_mrn_vsn, 'MBtoPM', 'up', 1),
  'PM->MC' = get.vecSYM_dex.dir.FC(dex.all$dex.mrn$dex_stage_mrn_vsn, 'PMtoMC', 'up', 1),
  'MC->MM' = get.vecSYM_dex.dir.FC(dex.all$dex.mrn$dex_stage_mrn_vsn, 'MCtoMM', 'up', 1),
  'MM->B' = get.vecSYM_dex.dir.FC(dex.all$dex.mrn$dex_stage_mrn_vsn, 'MMtoB', 'up', 1),
  'B->S' = get.vecSYM_dex.dir.FC(dex.all$dex.mrn$dex_stage_mrn_vsn, 'BtoS', 'up', 1),
  'S->PMN' = get.vecSYM_dex.dir.FC(dex.all$dex.mrn$dex_stage_mrn_vsn, 'StoPMN', 'up', 1)
)

l.dexNfc1_SYM_gnp.mrn_DN <- list(
  'MB->PM' = get.vecSYM_dex.dir.FC(dex.all$dex.mrn$dex_stage_mrn_vsn, 'MBtoPM', 'down', 1),
  'PM->MC' = get.vecSYM_dex.dir.FC(dex.all$dex.mrn$dex_stage_mrn_vsn, 'PMtoMC', 'down', 1),
  'MC->MM' = get.vecSYM_dex.dir.FC(dex.all$dex.mrn$dex_stage_mrn_vsn, 'MCtoMM', 'down', 1),
  'MM->B' = get.vecSYM_dex.dir.FC(dex.all$dex.mrn$dex_stage_mrn_vsn, 'MMtoB', 'down', 1),
  'B->S' = get.vecSYM_dex.dir.FC(dex.all$dex.mrn$dex_stage_mrn_vsn, 'BtoS', 'down', 1),
  'S->PMN' = get.vecSYM_dex.dir.FC(dex.all$dex.mrn$dex_stage_mrn_vsn, 'StoPMN', 'down', 1)
)

l.dexNfc1_SYM_gnp.pro_UP <- list(
  'MB->PM' = get.vecSYM_dex.dir.FC(dex.all$dex.pro$dex_stage_pro_vsn, 'MBtoPM', 'up', 1),
  'PM->MC' = get.vecSYM_dex.dir.FC(dex.all$dex.pro$dex_stage_pro_vsn, 'PMtoMC', 'up', 1),
  'MC->MM' = get.vecSYM_dex.dir.FC(dex.all$dex.pro$dex_stage_pro_vsn, 'MCtoMM', 'up', 1),
  'MM->B' = get.vecSYM_dex.dir.FC(dex.all$dex.pro$dex_stage_pro_vsn, 'MMtoB', 'up', 1),
  'B->PMN' = get.vecSYM_dex.dir.FC(dex.all$dex.pro$dex_stage_pro_vsn, 'StoPMN', 'up', 1)
)

l.dexNfc1_SYM_gnp.pro_DN <- list(
  'MB->PM' = get.vecSYM_dex.dir.FC(dex.all$dex.pro$dex_stage_pro_vsn, 'MBtoPM', 'down', 1),
  'PM->MC' = get.vecSYM_dex.dir.FC(dex.all$dex.pro$dex_stage_pro_vsn, 'PMtoMC', 'down', 1),
  'MC->MM' = get.vecSYM_dex.dir.FC(dex.all$dex.pro$dex_stage_pro_vsn, 'MCtoMM', 'down', 1),
  'MM->B' = get.vecSYM_dex.dir.FC(dex.all$dex.pro$dex_stage_pro_vsn, 'MMtoB', 'down', 1),
  'B->PMN' = get.vecSYM_dex.dir.FC(dex.all$dex.pro$dex_stage_pro_vsn, 'StoPMN', 'down', 1)
)


# write_clip(l.dexNfc1_SYM_gnp.mrn_UP$`MB->PM`)

import from ChEA3

# below there is a save & load file! Data could use better names 

# ## transcriptome
# # up
# res_MBtoPM_mrn_up <- recieve_ChEA3(l.dexNfc1_SYM_gnp.mrn_UP$`MB->PM`, "MBtoPM_mrn_up")
# res_PMtoMC_mrn_up <- recieve_ChEA3(l.dexNfc1_SYM_gnp.mrn_UP$`PM->MC`, "PMtoMC_mrn_up")
# res_MCtoMM_mrn_up <- recieve_ChEA3(l.dexNfc1_SYM_gnp.mrn_UP$`MC->MM`, "MCtoMM_mrn_up")
# res_MMtoB_mrn_up <- recieve_ChEA3(l.dexNfc1_SYM_gnp.mrn_UP$`MM->B`, "MMtoB_mrn_up")
# res_BtoS_mrn_up <- recieve_ChEA3(l.dexNfc1_SYM_gnp.mrn_UP$`B->S`, "BtoS_mrn_up")
# res_StoPMN_mrn_up <- recieve_ChEA3(l.dexNfc1_SYM_gnp.mrn_UP$`S->PMN`, "StoPMN_mrn_up")
# 
# # down
# res_MBtoPM_mrn_DN <- recieve_ChEA3(l.dexNfc1_SYM_gnp.mrn_DN$`MB->PM`, "MBtoPM_mrn_DN")
# res_PMtoMC_mrn_DN <- recieve_ChEA3(l.dexNfc1_SYM_gnp.mrn_DN$`PM->MC`, "PMtoMC_mrn_DN")
# res_MCtoMM_mrn_DN <- recieve_ChEA3(l.dexNfc1_SYM_gnp.mrn_DN$`MC->MM`, "MCtoMM_mrn_DN")
# res_MMtoB_mrn_DN <- recieve_ChEA3(l.dexNfc1_SYM_gnp.mrn_DN$`MM->B`, "MMtoB_mrn_DN")
# res_BtoS_mrn_DN <- recieve_ChEA3(l.dexNfc1_SYM_gnp.mrn_DN$`B->S`, "BtoS_mrn_DN")
# res_StoPMN_mrn_DN <- recieve_ChEA3(l.dexNfc1_SYM_gnp.mrn_DN$`S->PMN`, "StoPMN_mrn_DN")
# 
# ## proteome
# # up
# res_MBtoPM_pro_up <- recieve_ChEA3(l.dexNfc1_SYM_gnp.pro_UP$`MB->PM`, "MBtoPM_pro_up")
# res_PMtoMC_pro_up <- recieve_ChEA3(l.dexNfc1_SYM_gnp.pro_UP$`PM->MC`, "PMtoMC_pro_up")
# res_MCtoMM_pro_up <- recieve_ChEA3(l.dexNfc1_SYM_gnp.pro_UP$`MC->MM`, "MCtoMM_pro_up")
# res_MMtoB_pro_up <- recieve_ChEA3(l.dexNfc1_SYM_gnp.pro_UP$`MM->B`, "MMtoB_pro_up")
# # down
# res_MBtoPM_pro_DN <- recieve_ChEA3(l.dexNfc1_SYM_gnp.pro_DN$`MB->PM`, "MBtoPM_pro_DN")
# res_PMtoMC_pro_DN <- recieve_ChEA3(l.dexNfc1_SYM_gnp.pro_DN$`PM->MC`, "PMtoMC_pro_DN")
# res_MCtoMM_pro_DN <- recieve_ChEA3(l.dexNfc1_SYM_gnp.pro_DN$`MC->MM`, "MCtoMM_pro_DN")
# res_MMtoB_pro_DN <- recieve_ChEA3(l.dexNfc1_SYM_gnp.pro_DN$`MM->B`, "MMtoB_pro_DN")
# 
# 
# save(res_MBtoPM_mrn_up,
#       res_PMtoMC_mrn_up,
#       res_MCtoMM_mrn_up,
#       res_MMtoB_mrn_up,
#       res_BtoS_mrn_up,
#       res_StoPMN_mrn_up,
#       
#       res_MBtoPM_mrn_DN,
#       res_PMtoMC_mrn_DN,
#       res_MCtoMM_mrn_DN,
#       res_MMtoB_mrn_DN,
#       res_BtoS_mrn_DN,
#       res_StoPMN_mrn_DN,
#       
#       res_MBtoPM_pro_up,
#       res_PMtoMC_pro_up,
#       res_MCtoMM_pro_up,
#       res_MMtoB_pro_up,
#       
#       res_MBtoPM_pro_DN,
#       res_PMtoMC_pro_DN,
#       res_MCtoMM_pro_DN,
#       res_MMtoB_pro_DN,
#      file = paste0(dir_moa, "/data/final/preps/enrichment/", "res.ChEA3.RData"))

load(file = paste0(dir_moa,  "/data/final/preps/enrichment/", "res.ChEA3.RData"))

prep imported data

prep_4_TFe_plot <- function(res_ChEA3, DEX, COMP){
  # format to num
  res_ChEA3$`Integrated--topRank`$Score <- as.numeric(res_ChEA3$`Integrated--topRank`$Score)
  res_ChEA3$`Integrated--topRank`$Rank <- as.numeric(res_ChEA3$`Integrated--topRank`$Rank)
  # add number of overlapping genes
  res_ChEA3$`Integrated--topRank`$Overlapping_num <- str_count(res_ChEA3$`Integrated--topRank`$Overlapping_Genes, '\\w+')
  # add FC
  add_Fc_mrn <- function(res_ChEA3, DEX, COMP){
    dex <- DEX %>% filter(comparison == COMP)
    res_ChEA3$`Integrated--topRank`$logFC <- dex$logFC[match(res_ChEA3$`Integrated--topRank`$TF, dex$SYMBOL)]
    res_ChEA3$`Integrated--topRank`$AveExpr <- dex$AveExpr[match(res_ChEA3$`Integrated--topRank`$TF, dex$SYMBOL)]
    return(res_ChEA3)}
  res_ChEA3 <- add_Fc_mrn(res_ChEA3, DEX, COMP)
  # 1 - comp score
  res_ChEA3$`Integrated--topRank`$score_1 <- 1-  res_ChEA3$`Integrated--topRank`$Score
  return(res_ChEA3)
}

## mrn
# up
prp_MBtoPM_mrn_up <- prep_4_TFe_plot(res_MBtoPM_mrn_up, dex.all$dex.mrn$dex_stage_mrn_vsn, 'MBtoPM')
prp_PMtoMC_mrn_up <- prep_4_TFe_plot(res_PMtoMC_mrn_up, dex.all$dex.mrn$dex_stage_mrn_vsn, 'PMtoMC')
prp_MCtoMM_mrn_up <- prep_4_TFe_plot(res_MCtoMM_mrn_up, dex.all$dex.mrn$dex_stage_mrn_vsn, 'MCtoMM')
prp_MMtoB_mrn_up <- prep_4_TFe_plot(res_MMtoB_mrn_up, dex.all$dex.mrn$dex_stage_mrn_vsn, 'MMtoB')
prp_BtoS_mrn_up <- prep_4_TFe_plot(res_BtoS_mrn_up, dex.all$dex.mrn$dex_stage_mrn_vsn, 'BtoS')
prp_StoPMN_mrn_up <- prep_4_TFe_plot(res_StoPMN_mrn_up, dex.all$dex.mrn$dex_stage_mrn_vsn, 'StoPMN')

# down
prp_MBtoPM_mrn_DN <- prep_4_TFe_plot(res_MBtoPM_mrn_DN, dex.all$dex.mrn$dex_stage_mrn_vsn, 'MBtoPM')
prp_PMtoMC_mrn_DN <- prep_4_TFe_plot(res_PMtoMC_mrn_DN, dex.all$dex.mrn$dex_stage_mrn_vsn, 'PMtoMC')
prp_MCtoMM_mrn_DN <- prep_4_TFe_plot(res_MCtoMM_mrn_DN, dex.all$dex.mrn$dex_stage_mrn_vsn, 'MCtoMM')
prp_MMtoB_mrn_DN <- prep_4_TFe_plot(res_MMtoB_mrn_DN, dex.all$dex.mrn$dex_stage_mrn_vsn, 'MMtoB')
prp_BtoS_mrn_DN <- prep_4_TFe_plot(res_BtoS_mrn_DN, dex.all$dex.mrn$dex_stage_mrn_vsn, 'BtoS')
prp_StoPMN_mrn_DN <- prep_4_TFe_plot(res_StoPMN_mrn_DN, dex.all$dex.mrn$dex_stage_mrn_vsn, 'StoPMN')




## pro
# up
prp_MBtoPM_pro_up <- prep_4_TFe_plot(res_MBtoPM_pro_up, dex.all$dex.pro$dex_stage_pro_vsn, 'MBtoPM')
prp_PMtoMC_pro_up <- prep_4_TFe_plot(res_PMtoMC_pro_up, dex.all$dex.pro$dex_stage_pro_vsn, 'PMtoMC')
prp_MCtoMM_pro_up <- prep_4_TFe_plot(res_MCtoMM_pro_up, dex.all$dex.pro$dex_stage_pro_vsn, 'MCtoMM')
prp_MMtoB_pro_up <- prep_4_TFe_plot(res_MMtoB_pro_up, dex.all$dex.pro$dex_stage_pro_vsn, 'MMtoB')
# dn
prp_MBtoPM_pro_DN <- prep_4_TFe_plot(res_MBtoPM_pro_DN, dex.all$dex.pro$dex_stage_pro_vsn, 'MBtoPM')
prp_PMtoMC_pro_DN <- prep_4_TFe_plot(res_PMtoMC_pro_DN, dex.all$dex.pro$dex_stage_pro_vsn, 'PMtoMC')
prp_MCtoMM_pro_DN <- prep_4_TFe_plot(res_MCtoMM_pro_DN, dex.all$dex.pro$dex_stage_pro_vsn, 'MCtoMM')
prp_MMtoB_pro_DN <- prep_4_TFe_plot(res_MMtoB_pro_DN, dex.all$dex.pro$dex_stage_pro_vsn, 'MMtoB')

plots

fct

breaks <- c(-1, 1 , 3, 6)

# prep for colored labels
plot_TFe_up <- function(prp, TIT){
  p <- prp$`Integrated--topRank` %>%
    filter(Rank < 16) %>%
    filter(!is.na(AveExpr)) %>%
    #arrange(desc(Rank)) %>%
    ggplot(aes(x= score_1, y = reorder(TF, -Rank), size = AveExpr)) +
    geom_point(aes(colour = logFC)) +
    #geom_text(aes(label = Overlapping_num), size  = 4, hjust = -0.1) +
    coord_cartesian(xlim = c(0.995, 1)) +
    scale_size_area() +
    scale_colour_gradient2(low = "yellow", mid = "orange", high = "red", midpoint = 1, breaks = breaks) +
    scale_x_continuous(breaks=c(0.995,1)) +
    ggtitle(TIT) +
    xlab('Enrichment factor') +
    ylab('Transcription factor') +
    guides(size = "none") +    
    theme_light(base_size=18) +
    theme(
        text = element_text(size=18),
        axis.text.x = element_text(angle = 30, vjust = 0.9, hjust= 1,  size=16),
        #legend.position = "none",
        #axis.text.x=element_text(angle = 30, vjust = 0.9, hjust= 1, color="black", size=18),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(), 
        panel.border = element_blank(), 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.line.x = element_line(color="black", size = 1),
        axis.line.y = element_line(color="black", size = 1),
        axis.text.y = element_text(color="black", size=18),
        axis.ticks.y= element_line(color="black", size = 1),
        axis.ticks.x= element_line(color="black", size = 1))
  return(p)
}


plot_TFe_dn <- function(prp, TIT){
  p <- prp$`Integrated--topRank` %>%
    filter(Rank < 16) %>%
    filter(!is.na(AveExpr)) %>%
    #arrange(desc(Rank)) %>%
    ggplot(aes(x= score_1, y = reorder(TF, -Rank), size = AveExpr)) +
    geom_point(aes(colour = logFC)) +
    #geom_text(aes(label = Overlapping_num), size  = 4, hjust = -0.1) +
    coord_cartesian(xlim = c(0.995, 1)) +
    scale_size_area() +
    scale_colour_gradient2(low = "darkblue", mid = "yellow", high = "yellow", midpoint = 1, breaks = breaks) +
    ggtitle(TIT) +
    xlab('Enrichment factor') +
    ylab('Transcription factor') +
    guides(size = "none") +    
    theme_light(base_size=20) +
    theme(text = element_text(size=20),
          axis.text.x = element_text(angle = 30, vjust = 0.9, hjust= 1))
  return(p)
}

transcriptome

# up
p.TFe_MBtoPM_mrn_up <- plot_TFe_up(prp_MBtoPM_mrn_up, 'MB->PM')
p.TFe_PMtoMC_mrn_up <- plot_TFe_up(prp_PMtoMC_mrn_up, 'PM->MC')
p.TFe_MCtoMM_mrn_up <- plot_TFe_up(prp_MCtoMM_mrn_up, 'MC->MM')
p.TFe_MMtoB_mrn_up <- plot_TFe_up(prp_MMtoB_mrn_up, 'MM->B')
p.TFe_BtoS_mrn_up <- plot_TFe_up(prp_BtoS_mrn_up, 'B->S')
p.TFe_StoPMN_mrn_up <- plot_TFe_up(prp_StoPMN_mrn_up, 'S->PMN')

# fig
fig_TFe_mrn_up <- ggarrange(
  p.TFe_MBtoPM_mrn_up + theme_noX, 
  p.TFe_PMtoMC_mrn_up + theme_nada,
  p.TFe_MCtoMM_mrn_up+ theme_nada,
  p.TFe_MMtoB_mrn_up + theme_nada,
  p.TFe_BtoS_mrn_up + theme_nada,
  p.TFe_StoPMN_mrn_up  + theme_nada,
  common.legend = T,
  legend = 'right',
  nrow = 1,
  ncol = 6) 
#heights=c(0.95,1))
fig_TFe_mrn_up

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_5/", "fig_TFe_mrn_up", ".svg"),
  plot = fig_TFe_mrn_up,
  device = svg,
  width = 15.5,
  height = 5)


# DN
p.TFe_MBtoPM_mrn_DN <- plot_TFe_dn(prp_MBtoPM_mrn_DN, 'MB->PM')
p.TFe_PMtoMC_mrn_DN <- plot_TFe_dn(prp_PMtoMC_mrn_DN, 'PM->MC')
p.TFe_MCtoMM_mrn_DN <- plot_TFe_dn(prp_MCtoMM_mrn_DN, 'MC->MM')
p.TFe_MMtoB_mrn_DN <- plot_TFe_dn(prp_MMtoB_mrn_DN, 'MM->B')
p.TFe_BtoS_mrn_DN <- plot_TFe_dn(prp_BtoS_mrn_DN, 'B->S')
p.TFe_StoPMN_mrn_DN <- plot_TFe_dn(prp_StoPMN_mrn_DN, 'S->PMN')

# fig
fig_TFe_mrn_DN <- ggarrange(
  p.TFe_MBtoPM_mrn_DN + theme_noX, 
  p.TFe_PMtoMC_mrn_DN + theme_nada,
  p.TFe_MCtoMM_mrn_DN+ theme_nada,
  p.TFe_MMtoB_mrn_DN + theme_nada,
  p.TFe_BtoS_mrn_DN  + theme_nada,
  p.TFe_StoPMN_mrn_DN  + theme_nada,
  common.legend = T,
  legend = 'right',
  nrow = 1,
  ncol = 6)
fig_TFe_mrn_DN

ggsave(
  filename = paste0(dir_moa, "/figures/fig4/", "fig_TFe_mrn_DN", ".png"),
  plot = fig_TFe_mrn_DN,
  device = png,
  width = 23,
  height = 5)


## Proteome
# up
p.TFe_MBtoPM_pro_up <- plot_TFe_up(prp_MBtoPM_pro_up, 'MBtoPM')
p.TFe_PMtoMC_pro_up <- plot_TFe_up(prp_PMtoMC_pro_up, 'PMtoMC')
p.TFe_MCtoMM_pro_up <- plot_TFe_up(prp_MCtoMM_pro_up, 'MCtoMM')
p.TFe_MMtoB_pro_up <- plot_TFe_up(prp_MMtoB_pro_up, 'MMtoB')

fig_TFe_pro_up <- ggarrange(
  NULL,
  p.TFe_MBtoPM_pro_up + theme_noX, 
  p.TFe_PMtoMC_pro_up + theme_nada,
  NULL,
  p.TFe_MCtoMM_pro_up,
  p.TFe_MMtoB_pro_up+ theme_noY,
  common.legend = T,
  legend = 'right',
  nrow = 2,
  ncol = 3,
  heights=c(0.95,1))
fig_TFe_pro_up

ggsave(
  filename = paste0(dir_moa, "/figures/fig4/", "fig_TFe_pro_up", ".png"),
  plot = fig_TFe_pro_up,
  device = png,
  width = 11,
  height = 8.5)

# dn
p.TFe_MBtoPM_pro_DN <- plot_TFe_dn(prp_MBtoPM_pro_DN, 'MBtoPM')
p.TFe_PMtoMC_pro_DN <- plot_TFe_dn(prp_PMtoMC_pro_DN, 'PMtoMC')
p.TFe_MCtoMM_pro_DN <- plot_TFe_dn(prp_MCtoMM_pro_DN, 'MCtoMM')
p.TFe_MMtoB_pro_DN <- plot_TFe_dn(prp_MMtoB_pro_DN, 'MMtoB')

fig_TFe_pro_DN <- ggarrange(
  NULL,
  p.TFe_MBtoPM_pro_DN + theme_noX, 
  p.TFe_PMtoMC_pro_DN + theme_nada,
  NULL,
  p.TFe_MCtoMM_pro_DN,
  p.TFe_MMtoB_pro_DN+ theme_noY,
  common.legend = T,
  legend = 'right',
  nrow = 2,
  ncol = 3,
  heights=c(0.95,1))
fig_TFe_pro_DN

ggsave(
  filename = paste0(dir_moa, "/figures/fig4/", "fig_TFe_pro_DN", ".png"),
  plot = fig_TFe_pro_DN,
  device = png,
  width = 11,
  height = 8.5)

Stage specific TFs: upset

library(ComplexUpset)

l.TFs.mrn.up <- list(
  MBtoPM.up = prp_MBtoPM_mrn_up$`Integrated--topRank`,
  PMtoMC.up = prp_PMtoMC_mrn_up$`Integrated--topRank`,
  MCtoMM.up = prp_MCtoMM_mrn_up$`Integrated--topRank`,
  MMtoB.up = prp_MMtoB_mrn_up$`Integrated--topRank`,
  BtoS.up = prp_BtoS_mrn_up$`Integrated--topRank`,
  StoPMN.up = prp_StoPMN_mrn_up$`Integrated--topRank`
  #MBtoPM.dn = prp_MBtoPM_mrn_DN$`Integrated--topRank`,
  #PMtoMC.dn = prp_PMtoMC_mrn_DN$`Integrated--topRank`,
  #MCtoMM.dn = prp_MCtoMM_mrn_DN$`Integrated--topRank`,
  #MMtoB.dn = prp_MMtoB_mrn_DN$`Integrated--topRank`,
  #BtoS.dn = prp_BtoS_mrn_DN$`Integrated--topRank`,
  #StoPMN.dn = prp_StoPMN_mrn_DN$`Integrated--topRank`
)

clean_steps <- function(TF.df){
  res <- TF.df %>% filter(Score < 0.05)
  res$step <- sub("_.*","",res$`Query Name`)
  return(res)
}

TFs.mrn.filt <- purrr::map(l.TFs.mrn.up, clean_steps)

top_steps <- function(TF.df){
  res <- TF.df %>% filter(Rank < 15)
  res$step <- sub("_.*","",res$`Query Name`)
  return(res)
}

TFs.mrn.top <- purrr::map(l.TFs.mrn.up, top_steps)
TFs.mrn.top <- purrr::reduce(TFs.mrn.top, bind_rows)

# create TRUE/FALSE table
allTF <- bind_rows(TFs.mrn.filt, .id = "group")

TFs.4upset <- data.frame(TF = unique(allTF$TF))
TFs.4upset$MBtoPM <- ifelse(TFs.4upset$TF %in% TFs.mrn.filt$MBtoPM.up$TF, T, F)
TFs.4upset$PMtoMC <- ifelse(TFs.4upset$TF %in% TFs.mrn.filt$PMtoMC.up$TF, T, F)
TFs.4upset$MCtoMM <- ifelse(TFs.4upset$TF %in% TFs.mrn.filt$MCtoMM.up$TF, T, F)
TFs.4upset$MMtoB <- ifelse(TFs.4upset$TF %in% TFs.mrn.filt$MMtoB.up$TF, T, F)
TFs.4upset$BtoS <- ifelse(TFs.4upset$TF %in% TFs.mrn.filt$BtoS.up$TF, T, F)
TFs.4upset$StoPMN <- ifelse(TFs.4upset$TF %in% TFs.mrn.filt$StoPMN.up$TF, T, F)


# get the features to get top 3 per group
TFs.4upset.features <- TFs.4upset
# for(i in 1:ncol(TFs.4upset.features)){
# 
#     if(is.logical(TFs.4upset.features[, i]) == TRUE) TFs.4upset.features[, i] <- as.numeric(TFs.4upset.features[, i]) 
# 
# }
# TFs.4upset.features

TFs.4upset.features$MBtoPM <- ifelse(TFs.4upset.features$MBtoPM == T, 'MBtoPM', NA)
TFs.4upset.features$PMtoMC <- ifelse(TFs.4upset.features$PMtoMC == T, 'PMtoMC', NA)
TFs.4upset.features$MCtoMM <- ifelse(TFs.4upset.features$MCtoMM == T, 'MCtoMM', NA)
TFs.4upset.features$MMtoB <- ifelse(TFs.4upset.features$MMtoB == T, 'MMtoB', NA)
TFs.4upset.features$BtoS <- ifelse(TFs.4upset.features$BtoS == T, 'BtoS', NA)
TFs.4upset.features$StoPMN <- ifelse(TFs.4upset.features$StoPMN == T, 'StoPMN', NA)


intersecting.TFs <- TFs.4upset.features %>%
  unite(col = "intersection", -c("TF"), sep = ";", na.rm = TRUE)

table(intersecting.TFs$intersection)

intersecting.TFs %>%
  group_by(intersection) %>%
  summarise(n = n()) %>%
  arrange(desc(n))

TFs.all <- intersecting.TFs %>% filter(intersection == 'MBtoPM;PMtoMC;MCtoMM;MMtoB;BtoS;StoPMN') %>% 
  rownames_to_column('row') %>% 
  mutate(rank = as.numeric(row)) %>%
  top_n(rank, n=-5) %>% pull(TF)
TFs.notLast <- intersecting.TFs %>% filter(intersection == 'MBtoPM;PMtoMC;MCtoMM;MMtoB;BtoS') %>% 
  rownames_to_column('row') %>% 
  mutate(rank = as.numeric(row)) %>%
  top_n(rank, n=-5) %>% pull(TF)

TFs.MBtoPM <- intersecting.TFs %>% filter(intersection == 'MBtoPM') %>% 
  rownames_to_column('row') %>% 
  mutate(rank = as.numeric(row)) %>%
  top_n(rank, n=-5) %>% pull(TF)
TFs.PMtoMC <- intersecting.TFs %>% filter(intersection == 'PMtoMC') %>% 
  rownames_to_column('row') %>% 
  mutate(rank = as.numeric(row)) %>%
  top_n(rank, n=-5) %>% pull(TF)
TFs.MCtoMM <- intersecting.TFs %>% filter(intersection == 'MCtoMM') %>% 
  rownames_to_column('row') %>% 
  mutate(rank = as.numeric(row)) %>%
  top_n(rank, n=-5) %>% pull(TF)
TFs.MMtoB <- intersecting.TFs %>% filter(intersection == 'MMtoB') %>% 
  rownames_to_column('row') %>% 
  mutate(rank = as.numeric(row)) %>%
  top_n(rank, n=-5) %>% pull(TF)
TFs.BtoS <- intersecting.TFs %>% filter(intersection == 'BtoS') %>% 
  rownames_to_column('row') %>% 
  mutate(rank = as.numeric(row)) %>%
  top_n(rank, n=-3) %>% pull(TF)
TFs.StoPMN <- intersecting.TFs %>% filter(intersection == 'StoPMN') %>% 
  rownames_to_column('row') %>% 
  mutate(rank = as.numeric(row)) %>%
  top_n(rank, n=-5) %>% pull(TF)
TFs.MBtoMC <- intersecting.TFs %>% filter(intersection == 'MBtoPM;PMtoMC') %>% 
  rownames_to_column('row') %>% 
  mutate(rank = as.numeric(row)) %>%
  top_n(rank, n=-5) %>% pull(TF)
TFs.MBtoMM <- intersecting.TFs %>% filter(intersection == 'MBtoPM;PMtoMC;MCtoMM') %>% 
  rownames_to_column('row') %>% 
  mutate(rank = as.numeric(row)) %>%
  top_n(rank, n=-5) %>% pull(TF)
TFs.MBtoB <- intersecting.TFs %>% filter(intersection == 'MBtoPM;PMtoMC;MCtoMM;MMtoB') %>% 
  rownames_to_column('row') %>% 
  mutate(rank = as.numeric(row)) %>%
  top_n(rank, n=-5) %>% pull(TF)

TFs.foi <- c(TFs.all, TFs.notLast, TFs.MBtoPM, TFs.PMtoMC, TFs.MCtoMM, TFs.MMtoB, TFs.BtoS, TFs.StoPMN, TFs.MBtoMC, TFs.MBtoMM, TFs.MBtoB)


# plot prep
TF.groups <- colnames(TFs.4upset)[-1]

# plot
p.TF.upset <- ComplexUpset::upset(
  TFs.4upset, 
  TF.groups,
  name='Steps', 
  width_ratio=0.1,
  intersections=list(
    c("MBtoPM"),
    c("PMtoMC"),
    c("MCtoMM"),
    c("MMtoB"),
    c("BtoS"),
    c("StoPMN"),
    c("MBtoPM", "PMtoMC", "MCtoMM", "MMtoB" , "BtoS" ,"StoPMN"),
    c("MBtoPM", "PMtoMC", "MCtoMM", "MMtoB" , "BtoS"),
    c("MBtoPM", "PMtoMC", "MCtoMM", "MMtoB"),
    c("MBtoPM", "PMtoMC", "MCtoMM"),
    c("MBtoPM", "PMtoMC")),
  base_annotations=list(
    'Intersection size'=(
      intersection_size(
        bar_number_threshold=1,
        color='grey9',
        fill='grey80',
        text_mapping = (aes(fontface = "bold")))
  )))

p.TF.upset

ggsave(
  filename = paste0(dir_moa, "/figures/final/fig_5/", "p.TF.upset", ".svg"),
  plot = p.TF.upset,
  device = svg,
  width = 12,
  height = 7)



 
#       themes = upset_default_themes(text=element_text(size=20),
#         set_sizes=(upset_set_size() + 
#                  theme(axis.text.x=element_text(angle=35), text=element_text(size=15)))))
# 
# 
#     geom_text_repel((
#         mapping = aes(label = ifelse(TF %in% TFs.foi, TF, NA)),
#         position = position_stack(vjust = 90),
#         na.rm = TRUE,
#         direction = "y",
#         segment.coor = 'transparent'))
#      ))
# 


# p.TF.upset <- ComplexUpset::upset(
#   TFs.4upset, 
#   TF.groups,
#   name='Steps', 
#   width_ratio=0.1,
#   intersections=list(
#     c("MBtoPM"),
#     c("PMtoMC"),
#     c("MCtoMM"),
#     c("MMtoB"),
#     c("BtoS"),
#     c("StoPMN"),
#     c("MBtoPM", "PMtoMC", "MCtoMM", "MMtoB" , "BtoS" ,"StoPMN"),
#     c("MBtoPM", "PMtoMC", "MCtoMM", "MMtoB" , "BtoS"),
#     c("MBtoPM", "PMtoMC", "MCtoMM", "MMtoB"),
#     c("MBtoPM", "PMtoMC", "MCtoMM"),
#     c("MBtoPM", "PMtoMC")
#   ),
#   base_annotations=list(
#     'Intersection size'=(
#       intersection_size(
#         bar_number_threshold=1,
#         color='grey9',
#         fill='grey80',
#         text_mapping = (aes(fontface = "bold"))
#       ) +
#         geom_text_repel((
#           mapping = aes(label = ifelse(TF %in% TFs.foi, TF, NA)),
#           position = position_stack(vjust = 90),
#           na.rm = TRUE,
#           direction = "y",
#           segment.coor = 'transparent'))
#      ),
#     themes = upset_default_themes(text=element_text(size=20)),
#     set_sizes=(upset_set_size() + 
#                  theme(axis.text.x=element_text(angle=35), text=element_text(size=15)))
#   ))
#                     
# p.TF.upset
#                     
#