library(DESeq2)
library(magrittr)
library(tibble)
library(dplyr)
library(tidyr)
library(knitr)
library(ggplot2)
library(SummarizedExperiment)
library(DEGreport)
library(ggrepel)
library(EnhancedVolcano)
library(openxlsx)
library(RColorBrewer)
library(org.Hs.eg.db)
library(pheatmap)
library(clusterProfiler)
library(ComplexHeatmap)
library(DOSE)

## using FDR 0.05 ##
logfc.thr <- 0
adj.pval.thr <- 0.05

# Volcano plot
obj <- readRDS(file = "../20191109_120832_DGE.rds")
DGE.results <- as.data.frame(obj$deseq$dge_res$`HSCMPP_THAL-HSCMPP_ND`)
colnames(DGE.results) <- c("log2FoldChange","lfcSE","baseMean","PValue","padj")

DGE.results <- readRDS("Fig1F_input.rds") # DGE output with padj < 0.05
data <- data.frame(gene = row.names(DGE.results), pval = -log10(DGE.results$padj), lfc = DGE.results$log2FoldChange)
data <- na.omit(data)
data <- mutate(data, color = case_when(data$lfc > logfc.thr & data$pval > -log10(adj.pval.thr) ~ "Overexpressed",
                                       data$lfc < -logfc.thr & data$pval > -log10(adj.pval.thr) ~ "Underexpressed",
                                       abs(data$lfc) < logfc.thr | data$pval < -log10(adj.pval.thr) ~ "NonSignificant"))
data <- data[order(data$pval, decreasing = TRUE),]

selected_genes <-read.table(file = "Selected_genes", header = TRUE)
rownames(selected_genes) <- selected_genes$GENE

highl2 <- subset(data, subset = gene %in% selected_genes$GENE, color != "NonSignificant")

vol <- ggplot(data, aes(x = lfc, y = pval, colour = color)) +
  theme_bw(base_size = 12) +
  theme(legend.position = "right") +
  xlab("log2 Fold Change") + 
  ylab(expression(-log[10]("adjusted p-value"))) +
  geom_hline(yintercept = -log10(adj.pval.thr), colour = "darkgray")
if (logfc.thr > 0) {
  vol <- vol +
    geom_vline(xintercept = logfc.thr, colour = "darkgray") +
    geom_vline(xintercept = -logfc.thr, colour = "darkgray")
}
vol <- vol +
  geom_point(size = 2, alpha = 0.8, na.rm = T) +
  scale_color_manual(name = "Expression",
                     values = c(Overexpressed = "#CD4F39",
                                Underexpressed = "#0B5394",
                                NonSignificant = "darkgray")) +
  geom_label_repel(data = highl2, aes(label = gene), show.legend = FALSE, max.overlaps = 20)
ggsave(filename = "Fig1F.png",
       plot = vol,
       width = 9, height = 7, dpi = 600)