ffl: Identifies Feed Forward Gene Networks

#file: DEm-6-summary-plot-fn-v2-percent-denom-upordown-genes-not-all-DE-genes

# 7. summary plot ---------------------------------------------------------

#added to imports
#library(ggplot2)
#library(ggpubr)

# summary plot function ---------------------------------------------------


#' @title DEm_summary_plot
#' @description summarizes results in FFL plot
#'
#' @param arm1 DEm_arm1_miRNA_gene output
#' @param arm2 DEm_arm2_tf_gene output
#' @param ge DEm_ge_analysis output
#' @param ffl_fdr_cutoff FDR cutoff for FFLs (FDR from gene enrichment analysis)
#'
#' @return list of data.frame (modified ge) and plots
#'

DEm_summary_plot <- function(arm1, arm2, ge, ffl_fdr_cutoff){

  #####ffl plot (bubble plot)
  #convert miRNA and TF columns to factors
  ge$miRNA <- as.factor(ge$miRNA)
  ge$TF <- as.factor(ge$TF)

  #add fdr_threshold column
  ge$fdr_threshold <- as.character(ge$fdr < ffl_fdr_cutoff)
  ge$fdr_threshold[ge$fdr_threshold == "TRUE"] <- paste("FDR less than", ffl_fdr_cutoff)
  ge$fdr_threshold[ge$fdr_threshold == "FALSE"] <- paste("FDR greater/equal to", ffl_fdr_cutoff)

  #bubble plot: ffls
  bubble <- ggplot(ge, aes(x = reorder(TF, -percent_upordowngenes_tftargets), y = reorder(miRNA, percent_upordowngenes_mirtargets))) +
    geom_point(aes(size = gene_count, color = fdr_threshold), alpha = 0.3) +
    geom_text(aes(label = gene_count, size = 5)) +
    #geom_text(aes(label = paste("FDR =", formatC(fdr, format = "e", digits = 1.5)), size = 1.5, vjust = 2.5)) +
    scale_size(range = c(0.5, 11)) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
    scale_color_manual(values = c("#E7B800", "#9370DB")) +
    theme(legend.position = "top", plot.margin = unit(c(0.5, 0.5, 0.5, 0.5), "lines")) +
    labs(x = "TFs in FFL", y = "miRNAs in FFL") + guides(size = FALSE)

  #####mirna plot (bar plot)
  #for each miRNA in ge, indicate whether it is up or downregulated in the case group
  mirna_log2fc <- arm1$db_hits[arm1$db_hits$miRNA %in% ge$miRNA, c("miRNA", "log-ratio(miRNA)")]
  mirna_log2fc <- mirna_log2fc[!duplicated(mirna_log2fc), ]
  ge <- merge(x = ge, y = mirna_log2fc, by = "miRNA", all.x = TRUE)
  ge$miRNA_fc_direction <- as.character(ge$`log-ratio(miRNA)` > 0)
  ge$miRNA_fc_direction[ge$miRNA_fc_direction == "TRUE"] <- "increase"
  ge$miRNA_fc_direction[ge$miRNA_fc_direction == "FALSE"] <- "decrease"

  #subset ge and remove duplicate rows (otherwise, get wrong bar plots)
  mirna_percent <- ge[ , c("miRNA", "percent_DEgenes_mirtargets", "percent_upordowngenes_mirtargets", "miRNA_fc_direction")]
  mirna_percent <- mirna_percent[!duplicated(mirna_percent), ]

  #barplot: miRNA
  bar_mirna <- ggplot(mirna_percent, aes(x = reorder(miRNA, percent_upordowngenes_mirtargets), y = percent_upordowngenes_mirtargets, fill = miRNA_fc_direction)) +
    geom_bar(stat = "identity") +
    #geom_text(aes(label = paste(round(percent_DEgenes_mirtargets, digits = 1), "%")), size = 3, hjust = 1.2) +
    theme_minimal() +
    coord_flip() +
    theme(legend.position = "none") +
    theme(axis.title.y = element_blank()) +
    labs(y = "percent up/down genes explained by targets", fill = "FC direction")

  #####tf plot (bar plot)
  #for each TF in ge, indicate whether it is up or downregulated
  tf_log2fc <- arm2$db_hits[arm2$db_hits$TF %in% ge$TF, c("TF", "log-ratio(TF)")]
  tf_log2fc <- tf_log2fc[!duplicated(tf_log2fc), ]
  ge <- merge(x = ge, y = tf_log2fc, by = "TF", all.x = TRUE)
  ge$TF_fc_direction <- as.character(ge$`log-ratio(TF)` > 0)
  ge$TF_fc_direction[ge$TF_fc_direction == "TRUE"] <- "increase"
  ge$TF_fc_direction[ge$TF_fc_direction == "FALSE"] <- "decrease"

  #subset ge and remove duplicate rows (otherwise, get wrong bar plots)
  tf_percent <- ge[ , c("TF", "percent_DEgenes_tftargets", "percent_upordowngenes_tftargets", "TF_fc_direction")]
  tf_percent <- tf_percent[!duplicated(tf_percent), ]

  #barplot: TF
  bar_tf <- ggplot(tf_percent, aes(x = reorder(TF, -percent_upordowngenes_tftargets), y = percent_upordowngenes_tftargets, fill = TF_fc_direction)) +
    geom_bar(stat = "identity") +
    #geom_text(aes(label = paste(round(percent_DEgenes_tftargets, digits = 1), "%")), size = 2, hjust = 1.2) +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 90)) +
    theme(legend.position = "none") +
    theme(axis.title.x = element_blank()) +
    labs(y = "percent up/down genes explained by targets", fill = "FC direction")

  #####combine the three plots ()
  combined <- ggarrange(bar_tf, NULL, bubble, bar_mirna,
                        ncol = 2, nrow = 2,  align = "hv",
                        widths = c(2.5, 1), heights = c(1, 2.5),
                        common.legend = TRUE, legend = "top")

  return(list(new_ge = ge, bubble_ffl = bubble, bar_mirna = bar_mirna, bar_tf = bar_tf, combined_plot = combined))
}