First version with the server part.

TCGA_request/app.R

@@ -1,5 +1,14 @@
 library(shiny)
 
+## Library loadings
+library(tidyverse)
+library(GSVA)
+library(rstatix)
+library(ggpubr)
+library(coxphf)
+library(survival)
+library(survminer)
+
 ui <- fluidPage(
 
     # Application title
@@ -15,18 +24,339 @@ ui <- fluidPage(
 
         # Show a plot of the generated distribution
         mainPanel(
-           plotOutput("distPlot"),
-           plotOutput("plotEndo")
+           uiOutput("distPlot")
+           # plotOutput("plotEndo")
         )
     )
 )
 
 # Define server logic
 server <- function(input, output) {
+    vars<-reactiveValues()
+    vars$height<-900
+    
+    output$distPlot<- renderUI({
+      observeEvent(input$goButton, {})
+      height=vars$height
+      width=900
+      plotOutput("dsPlot", width=paste0(width,"px"), height = paste0(height, "px"))
+    })
+    
+    grafic<-eventReactive(input$goButton, {
+      progress <- shiny::Progress$new(min=0, max=6)
+      progress$set(message = "Start", value = 1)
+      
+      if (input$ds_input == "endo"){
+        vars$height=1500
+      }else{
+        vars$height=800
+      }
+      
+      ## Parameters definition
+      ds.list<-input$ds_input
+      gene<-input$gene_input
+      
+      for (ds in ds.list){
+        ## Loading pathways
+        load("../pathways.Rdata")
+        Pathways
+        
+        ## Loading dataset
+        datasets<-c("endo"="../data/endo_tcga_2020.Rdata",
+                    "uveal"="../data/uveal_tcga_MCG.Rdata",
+                    "prostate"="../data/prostate_tcga_MCG.Rdata")
+        load(datasets[ds])
+        
+        ## Gene expression
+        gex<-ex %>%
+          as.data.frame() %>%
+          add_column("Gene"=rownames(.)) %>% 
+          gather(SampleID, GEX, -Gene) %>% 
+          filter(Gene %in% gene) %>% 
+          mutate(MedExp=factor(GEX > median(GEX, na.rm=T), levels=c(FALSE,TRUE), labels=c("Low","High")))
+        
+        ### ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
+        ### Pathways Analysis
+        ### ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        
+        ## Runned once for each dataset and stored in an Object
+        # Perform GSVA function
+        # GSVA <- gsva(expr=as.matrix(ex), gset.idx.list=Pathways, 
+        #              mx.diff=TRUE, verbose=TRUE)
+        # summary(reshape2::melt(GSVA)$value)
+        # 
+        # # Vector to change pathways names
+        # pathways <- c("PD1 Signaling", "VEGF (Angiogenesis)",  
+        #               "Extracellular matrix interaction", "Epitelial Mesenchimal Transition", 
+        #               "Adherent junctions", "TGFB Signaling", "Inflamatory Signaling", 
+        #               "B cell activation","CTLA4 Tcell Signaling","Antigen Processing",
+        #               "IL6 Signaling", "CD8 T cell receptor",  "Antigen Presentation",
+        #               "IFN Gamma Signaling", "T cell activation")
+        # 
+        # rownames(GSVA) <- pathways
+        # colnames(GSVA) <- gsub("\\.", "-", colnames(GSVA))
+        
+        load("../GSVAs/GSVAs.Rdata")
+        if (ds == "endo"){GSVA<-GSVAendo}
+        if (ds == "uveal"){GSVA<-GSVAuveal}
+        if (ds == "prostate"){GSVA<-GSVAprostate}
+        
+        progress$set(message = "datasets loaded", value = 2)
+        
+        GSVAsummary<-GSVA %>%
+          as.data.frame %>% 
+          add_column("Path"=rownames(.)) %>% 
+          gather(SampleID, GSVA, -Path) %>%
+          merge(gex %>% select(-Gene)) %>% 
+          group_by(Path, MedExp) %>% 
+          summarise(Value=median(GSVA))
+        
+        GSVAstats<- GSVA %>%
+          as.data.frame %>% 
+          add_column("Path"=rownames(.)) %>% 
+          gather(SampleID, GSVA, -Path) %>%
+          merge(gex %>% select(-Gene)) %>% 
+          group_by(Path) %>% 
+          wilcox_test(GSVA~MedExp)
+        
+        GSVAstats["signif"]<-ifelse(GSVAstats$p < 0.05, "*", "")
+        
+        GSVAplot<-ggplot(GSVAsummary, aes(x=MedExp, y=Path))+
+          geom_tile(aes(fill=Value), color="white", size=1.75) +
+          scale_fill_gradientn(colours=c("darkblue", "blue3","snow","red3", "firebrick"),
+                               limits=c(summary(GSVAsummary$Value)[1],summary(GSVAsummary$Value)[6])) +
+          scale_x_discrete(position = "top", limits=c("High", "Low"))+
+          scale_y_discrete(limits=rev(rownames(GSVA)))+
+          guides(fill = guide_colorbar(barwidth = 6, barheight = 0.5))+
+          geom_text(data=GSVAstats,aes(x=2.9, label=signif), size=6, hjust=2, vjust=0.65, colour="black")+
+          # theme_minimal()+
+          theme(legend.justification = c(1, 2),
+                legend.position = "bottom",
+                legend.direction = "horizontal",
+                legend.text=element_text(size=8),
+                plot.title = element_text(hjust = 0.5, lineheight=.8, face="bold"),
+                text = element_text(size=12),
+                axis.ticks = element_blank(),
+                axis.title.x = element_blank(),
+                axis.title.y = element_blank(),
+                axis.text.x=element_text(face="bold",size=12),
+                axis.text.y=element_text(size=12),
+                panel.grid.major = element_blank(),
+                panel.border = element_blank(),
+                panel.background = element_blank(),
+                plot.margin=unit(c(1,2,1,1.5),"cm"))
+        progress$set(message = "Heatmap", value = 3)
+        ### ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        ### Survival Analysis
+        ### ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        
+        if (ds == "endo"){
+          labels<-merge(clin, gex, by.x = "PATIENT_ID", by.y="SampleID")
+          labels$OS <- as.numeric(as.character(factor(labels$OS_STATUS, levels=c("LIVING", "DECEASED"), labels=c(0,1))))
+          labels$OS.Months <- as.numeric(labels$OS_MONTHS)
+        }
+        if (ds == "uveal"){
+          labels<-merge(clin, gex, by.x = "Patient.ID", by.y="SampleID")
+        }
+        if (ds == "prostate"){
+          labels<-merge(clin, gex, by.x = "id", by.y="SampleID")
+          labels$OS <- labels$os
+          labels$OS.Months <- labels$OS_MONTHS
+        }
+        
+        fit_os <- survfit(Surv(OS.Months, OS) ~ MedExp, data=labels)
+        fit_os
+        
+        ### DFS
+        if (ds == "endo"){
+          labels$DFS <- as.numeric(as.character(factor(labels$DFS_STATUS, levels=c("DiseaseFree", "Recurred/Progressed"), labels=c(0,1))))
+          labels$DFS.Months <- as.numeric(labels$DFS_MONTHS)
+        }
+        if (ds == "prostate"){
+          labels$DFS <- labels$dfs
+          labels$DFS.Months <- as.numeric(labels$DFS_MONTHS)
+        }
+        
+        fit_dfs <- survfit(Surv(as.numeric(DFS.Months), DFS) ~ MedExp, data=labels)
+        fit_dfs
+        
+        # plots
+        km1 <- ggsurvplot(fit=fit_os, data=labels, 
+                          font.main = c(14, "plain", "black"),
+                          font.x = c(12, "plain", "black"),
+                          font.y = c(12, "plain", "black"),
+                          pval = TRUE, 
+                          legend.title = gene,
+                          legend.labs = c("Low", "High"),
+                          #legend = c(0.2, 0.2),
+                          palette=c("dodgerblue4","firebrick3"),
+                          font.legend = c(12, "plain", "black"), 
+                          ylab="Overall survival probability")
+        
+        
+        km2 <- ggsurvplot(fit=fit_dfs, data=labels, 
+                          font.main = c(14, "plain", "black"),
+                          font.x = c(12, "plain", "black"),
+                          font.y = c(12, "plain", "black"),
+                          pval = TRUE, 
+                          legend.title = gene,
+                          legend.labs = c("Low", "High"),
+                          #legend = c(0.2, 0.2),
+                          palette=c("dodgerblue4","firebrick3"),
+                          font.legend = c(12, "plain", "black"), 
+                          ylab="Progression free probability")
+        
+        km<-ggarrange(km1$plot, km2$plot, ncol=2, nrow=1, common.legend = T)
+        
+        progress$set(message = "Survival", value = 4)
+        # ggsave("",height=3.5, width=7)
+        
+        ### ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        ### Boxplot Recurrences
+        ### ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        
+        if (ds == "endo"){
+          bp2 <- ggplot(labels %>% filter(!is.na(DFS_STATUS)), aes(x=DFS_STATUS, y=GEX)) +
+            geom_boxplot(color="black", outlier.shape = NA) +
+            geom_jitter(position=position_jitter(0.2), size=1.5, 
+                        aes(x=DFS_STATUS, y=GEX, color=DFS_STATUS)) +
+            labs(title=gene, 
+                 x="", y = "log2(gene expression)") +
+            scale_color_manual(values=c("dodgerblue3", "firebrick3")) +
+            theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), panel.spacing = unit(1, "lines"),
+                  axis.text.x = element_text(size=12),
+                  axis.text.y = element_text(size=8),
+                  strip.text.x = element_text(size = 12))+
+            stat_compare_means(method = "wilcox.test",
+                               size=4, label="p.format", label.x=2.1)
+        }
+        if (ds == "uveal"){
+          bp2 <- ggplot(labels %>% 
+                          filter(!is.na(DFS)) %>% 
+                          mutate(DFS=factor(DFS, labels = c("DiseaseFree","Recurred/Progressed"))), 
+                        aes(x=DFS, y=GEX)) +
+            geom_boxplot(color="black", outlier.shape = NA) +
+            geom_jitter(position=position_jitter(0.2), size=1.5, 
+                        aes(x=DFS, y=GEX, color=DFS)) +
+            labs(title=gene, 
+                 x="", y = "log2(gene expression)") +
+            scale_color_manual(values=c("dodgerblue3", "firebrick3")) +
+            theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), panel.spacing = unit(1, "lines"),
+                  axis.text.x = element_text(size=12),
+                  axis.text.y = element_text(size=8),
+                  strip.text.x = element_text(size = 12))+
+            stat_compare_means(method = "wilcox.test",
+                               size=4, label="p.format", label.x=2.1)
+        }
+        if (ds == "prostate"){
+          bp2 <- ggplot(labels %>% 
+                          filter(!is.na(DFS_STATUS)) %>% 
+                          filter(DFS_STATUS != "[Not Available]"), 
+                        aes(x=DFS_STATUS, y=GEX)) +
+            geom_boxplot(color="black", outlier.shape = NA) +
+            geom_jitter(position=position_jitter(0.2), size=1.5, 
+                        aes(x=DFS_STATUS, y=GEX, color=DFS_STATUS)) +
+            labs(title=gene, 
+                 x="", y = "log2(gene expression)") +
+            scale_color_manual(values=c("dodgerblue3", "firebrick3")) +
+            theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), panel.spacing = unit(1, "lines"),
+                  axis.text.x = element_text(size=12),
+                  axis.text.y = element_text(size=8),
+                  strip.text.x = element_text(size = 12))+
+            stat_compare_means(method = "wilcox.test",
+                               size=4, label="p.format", label.x=2.1)
+        }
+        progress$set(message = "Boxplot", value = 5)
+        # jpeg("/shared/users/Sandra/01_collaborations/genes_piulats/HDAC6/uveal_HDAC6_boxplot1.jpg", 
+        # #     units="in", width=4,5, height=3.5, res=300)
+        # bp2
+        # dev.off()
+        
+        ### ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        ### Plot part 1
+        ### ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        
+        plotAll<-ggarrange(km, ggarrange(GSVAplot, bp2), ncol = 1)
+        plotAll<-annotate_figure(plotAll, top = text_grob(paste0(gene,"-",ds), 
+                                                          color = "black", face = "bold", size = 14))
+        # ggsave(paste0("result/",gene,"_",ds,"_results.png"), plotAll,width = 11, height = 9)
+      }
+      
+      if ("endo" %in% ds.list){
+        # -------------------------------------------------------------------------------
+        # by histological grade
+        
+        table(labels$hist_grade)
+        labels$hist_grade <- factor(labels$hist_grade, 
+                                    levels = c("G1", "G2", "G3", "High Grade"))
+        labels <- labels[order(labels$hist_grade),]
+        
+        
+        hg <- ggplot(labels, aes(x=hist_grade, y=GEX))+
+          geom_boxplot(color="black", outlier.shape = NA, aes(fill=hist_grade))+
+          scale_fill_brewer(palette = "Pastel2")+
+          geom_jitter(position=position_jitter(0.2), size=1.5,
+                      aes(x=hist_grade, y=GEX), alpha=0.5) +
+          guides(fill="none")+
+          labs(title=paste0(gene,", Histological grade"), 
+               x="", y = "log2(gene expression)") +
+          theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), panel.spacing = unit(1, "lines"),
+                axis.text.x = element_text(size=12),
+                axis.text.y = element_text(size=8),
+                strip.text.x = element_text(size = 12))+ 
+          stat_compare_means(method = "kruskal.test", size=4, label="p.format", label.x=4)
+        
+        table(labels$cluster)
+        labels$cluster <- factor(labels$cluster, levels = c("CN_low", "CN_high", "MSI", "POLE"))
+        labels <- labels[order(labels$cluster),]
+        
+        bp.cluster <- ggplot(labels, aes(x=cluster, y=GEX)) +
+          geom_boxplot(color="black", outlier.shape = NA, aes(fill=cluster)) +
+          scale_fill_brewer(palette = "Pastel2")+
+          scale_x_discrete(limits=levels(labels$cluster))+
+          guides(fill="none")+
+          geom_jitter(position=position_jitter(0.2), size=1.5, 
+                      aes(x=cluster, y=GEX), alpha=0.5) +
+          labs(title=paste0(gene," Molecular classification"), 
+               x="", y = "log2(gene expression)") +
+          theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), panel.spacing = unit(1, "lines"),
+                axis.text.x = element_text(size=12),
+                axis.text.y = element_text(size=8),
+                strip.text.x = element_text(size = 12))+
+          stat_compare_means(method = "kruskal.test", size=4, label="p.format", label.x=4)
+        
+        table(labels$histology, exclude=NULL)
+        labels$histology <- factor(labels$histology)
+        
+        bp.hg <- ggplot(labels, aes(x=histology, y=GEX)) +
+          geom_boxplot(color="black", outlier.shape = NA, aes(fill=histology)) +
+          scale_fill_brewer(palette = "Pastel2")+
+          guides(fill="none")+
+          geom_jitter(position=position_jitter(0.2), size=1.5, 
+                      aes(x=histology, y=GEX), alpha=0.5) +
+          labs(title=paste0(gene, ", Histology"), 
+               x="", y = "log2(gene expression)") +
+          theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), panel.spacing = unit(1, "lines"),
+                axis.text.x = element_text(size=12),
+                axis.text.y = element_text(size=8),
+                strip.text.x = element_text(size = 12))+
+          stat_compare_means(method = "kruskal.test", size=4, label="p.format", label.x=2)
+        progress$set(message = "Ploting", value = 6)
+        progress$close()
+        return(ggarrange(
+          plotAll,
+          ggarrange(hg, bp.cluster, bp.hg),
+          ncol=1
+        ))
+      }else{
+        return(plotAll)
+      }
+    })
+    
+    output$dsPlot <- renderPlot({
+      grafic()
+    })
 
-    # output$distPlot <- renderPlot({
-    #     
-    # })
 }
 
 # Run the application