diff --git a/TCGA_request/app.R b/TCGA_request/app.R index 88627ef..9b52d28 100644 --- a/TCGA_request/app.R +++ b/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