added volcano plot

Paper_figures.Rmd

@@ -445,3 +445,102 @@ dev.off()
 
 ```
 
+
+
+volcano plot
+```{r echo=FALSE,echo=FALSE, warning=FALSE,message=FALSE} 
+
+ds <- c('ON Old vs ON Young' = '/beegfs/scratch/ric.dimicco/ric.dimicco/DiMiccoR_1653_RNAseq/results/05-DGE-ON/tables/DESeq-ON_Old-ON_Young-allgene.txt')
+```
+
+```{r echo=FALSE,echo=FALSE, warning=FALSE,message=FALSE} 
+library(EnhancedVolcano)
+library(dplyr)
+library(ggrepel)
+plotVolcano <- function(dge_res, adj.pval.thr, logfc.thr,titlex) {
+  data <- data.frame(gene = row.names(dge_res), pval = -log10(dge_res$FDR), lfc = dge_res$logFC)
+  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),]
+  highl <- head(subset(data, color == "Overexpressed"), 10)
+  lowl <-  head(subset(data, color == "Underexpressed"), 10)
+  tt <- rbind( highl,lowl)
+  vol <- ggplot(data, aes(x = lfc, y = pval, colour = color)) +
+    theme(legend.position = "right") +
+    xlab("log2 Fold Change") +
+    ylab(expression(-log[10]("adjusted p-value"))) +
+    geom_hline(yintercept = -log10(adj.pval.thr), colour = "darkgray") +
+    geom_point(size = 1, alpha = 0.8, na.rm = T) +
+    scale_color_manual(name = "Expression",
+                       values = c(Overexpressed = "#CD4F39",
+                                  Underexpressed = 'blue',#"#00488b",
+                                  NonSignificant = "darkgray")) +
+    geom_text_repel(data=tt,aes(x = lfc, y = pval,label=gene),show.legend= FALSE, box.padding = 0.3)+
+    theme_classic() +
+    theme(axis.text.x = element_text(colour="black",size=12),
+          axis.text.y = element_text(colour="black",size=12),
+          
+          legend.position = "bottom",
+          #legend.justification = 0.05,
+          panel.grid.major.x = element_line(colour = "grey80", linetype = "dashed"),
+          panel.grid.minor.y = element_blank(),
+          panel.grid.major.y = element_line(colour = "grey80", linetype = "dashed"),
+          panel.grid.minor.x = element_blank()) +
+    ggtitle(titlex)+theme(plot.title = element_text(size = 12, face = "bold"))
+  return(vol)
+}
+
+adj.pval.thr = 0.05
+logfc.thr = 0
+titlex = 'OLD vs Young'
+
+
+
+x <- list()
+xx <- list()
+folder = '/beegfs/scratch/ric.dimicco/ric.dimicco/DiMiccoR_1653_RNAseq/bulk_rnaseq'
+for (name in names(ds)) {
+  dt <- read.table(ds[[name]])
+  x[[name]] <- plotVolcano(dt, 0.05, 0,name) #myvolcano(dt,0.05,0,name)
+  xx[[name]] <-   EnhancedVolcano(dt,
+   lab = rownames(dt),
+   x = 'logFC',
+   y = 'PValue')+ ggtitle(name)
+  # 
+  svg(filename= paste(folder, paste0(name, '.svg'),sep = "/"),
+      width = 6, height = 6, pointsize=12)
+  plot(x[[name]])
+  dev.off()
+  
+  #svg(filename= paste(folder, paste0(name, '.svg'),sep = "/"),
+  #    width = 6, height = 5, pointsize=12)
+  #plot(xx[[name]])
+  #dev.off()
+}
+
+
+
+
+
+x <- list()
+xx <- list()
+folder = '/beegfs/scratch/ric.dimicco/ric.dimicco/DiMiccoR_1653_RNAseq/results'
+for (name in names(ds)) {
+  dt <- read.table(ds[[name]])
+  x[[name]] <- plotVolcano(dt, 0.05, 0,name) #myvolcano(dt,0.05,0,name)
+  xx[[name]] <-   EnhancedVolcano(dt,
+                                  lab = rownames(dt),
+                                  x = 'logFC',
+                                  y = 'PValue')+ ggtitle(name)
+  # 
+  svg(filename= paste(folder, paste0(name, '_noY2.svg'),sep = "/"),
+      width = 6, height = 6, pointsize=12)
+  plot(x[[name]])
+  dev.off()
+  
+
+}
+
+```