From 642df836a9f602cd02514e5538b136286dc61fbd Mon Sep 17 00:00:00 2001
From: marcelcosta <marcelcosta@protonmail.com>
Date: Tue, 3 Nov 2020 07:34:13 +0100
Subject: [PATCH] add core files

---
 app.R        | 169 +++++++++++++++++++++++++++++++++++++++++++++++++++
 elispots.Rmd |  95 +++++++++++++++++++++++++++++
 funcions.R   | 137 +++++++++++++++++++++++++++++++++++++++++
 3 files changed, 401 insertions(+)
 create mode 100644 app.R
 create mode 100644 elispots.Rmd
 create mode 100644 funcions.R

diff --git a/app.R b/app.R
new file mode 100644
index 0000000..660ef66
--- /dev/null
+++ b/app.R
@@ -0,0 +1,169 @@
+library(shiny)
+library(openxlsx)
+library(readxl)
+library(ggplot2)
+library(reshape2)
+library(dplyr)
+library(ggbeeswarm)
+library(magrittr)
+library(flextable)
+source("funcions.R")
+
+ui <- fluidPage(
+    
+    # Application title
+    titlePanel("ELISPOTs"),
+    
+    sidebarLayout(
+        sidebarPanel(
+            fileInput(inputId = "file1", label = "Dades", multiple = F),
+            selectInput(inputId = "test", "Test Estadístic", selected = "Ttest", choices = c("T-test (adj Holm)","Wilcoxon (adj Holm)")),
+            sliderInput(inputId = "umbral_pos", "Mínimo para positivo:", min = 0, max=100, step = 5, value = 10),
+            checkboxInput(inputId = "positive", label = "Mostrar positivitat", value = F),
+            checkboxInput(inputId = "showstats", label = "Mostrar estadística", value = F),
+            downloadButton("downloadData", "Descarregar Informe")
+        ),
+        
+        mainPanel(
+            plotOutput("distPlot"),
+            uiOutput("flexstats")
+        )
+    )
+)
+
+# Define server logic required to draw a histogram
+server <- function(input, output) {
+    dades<-reactiveValues()
+    dades$taula<-NULL
+    dades$stats<-NULL
+    observe({
+        if (!is.null(input$file1)){ 
+            dades$taula<-read_xlsx(input$file1$datapath)
+        }
+    })
+    
+    output$distPlot <- renderPlot({
+        observeEvent(dades$taula, {})
+        if (!is.null(dades$taula)){
+            
+            ctrl<-"Ctrl+"
+            mock<-"Mock"
+            
+            table<-dades$taula[,colnames(dades$taula) != "Groups"]
+            t_mean<-dcast(melt(table, id="Mice"),Mice~variable, mean, na.rm=T)
+        
+            t_substr<-data.frame("Mice"=t_mean[,1],
+                                 as.data.frame(t(apply(t_mean[,2:ncol(table)], 1, function(x) x-x[mock])))
+            )
+            t_mean_group<-merge(t_mean, unique(dades$taula[c("Mice","Groups")]), id="Mice")
+            mock_mean<-dcast(t_mean_group, Groups~., value.var=mock, mean)
+            mock_mean[,2]<-mock_mean[,2]*2
+            mock_mean[mock_mean$. < input$umbral_pos,2]<-input$umbral_pos
+            
+            t_substr<-merge(t_substr, unique(dades$taula[c("Mice","Groups")]), id="Mice")
+            t_substr<-t_substr[,c(1, ncol(t_substr), 2:(ncol(t_substr)-1))]
+            t_substr<-t_substr[,c("Mice", "Groups", colnames(t_substr)[!colnames(t_substr) %in% c("Mice", "Groups")])]
+            t_substr[,3:ncol(t_substr)]<-apply(t_substr[,3:ncol(t_substr)],2, function(x) replace(x, which(x < 0),0))
+            colnames(t_substr)<-c("Mice", "Groups", colnames(t_mean)[2:ncol(t_mean)])
+            
+            t_substr_gp<-t_substr
+            t_substr_gp[3:ncol(t_substr)]<-apply(t_substr[3:ncol(t_substr)], 2, function(x) gsub(".",",",x, fixed=T))
+            t_substr_gp<-t_substr_gp[order(factor(t_substr_gp$Mice, levels = unique(table$Mice))),]
+            doc<-t(t_substr_gp)
+            
+            write.xlsx(doc, "data4graphpad.xlsx",rowNames=T)
+            
+            t<-melt(t_substr[,!colnames(t_substr) %in% c(ctrl, mock)])
+            if (input$test == "T-test (adj Holm)"){
+                t_stats<-multi_stats(t, "value", "variable", "Groups", stat.test = "ttest")
+            }
+            if (input$test == "Wilcoxon (adj Holm)"){
+                t_stats<-multi_stats(t, "value", "variable", "Groups", stat.test = "wilcox")
+            }
+            
+            dades$stats<<-t_stats
+            t_stats<-t_stats %>% filter(p.signif != "ns")
+            
+            t_maps<-generate_labstats(t_stats, t, "value", "variable", "Groups")
+            
+            if (input$showstats == F){
+                t_stats<-as.data.frame(matrix(nrow=0, ncol=6))
+                colnames(t_stats)<-c("variable", "group1", "group2", "p.adj", "p.signif", "Method")
+                t_maps<-list()
+                t_maps[["label"]]<-as.data.frame(matrix(nrow = 0, ncol=2))
+                colnames(t_maps$label)<-c("x", "y")
+                t_maps[["brackets"]]<-as.data.frame(matrix(nrow = 0, ncol=4))
+                colnames(t_maps$brackets)<-c("y1", "y2", "x1", "x2")
+            }
+            
+            set.seed(123)
+            if (input$positive == T){
+                ggplot(melt(t_substr, id=c("Mice", ctrl, "Groups")), aes(variable, value))+
+                    labs(x="", y="Spots/2.5*10^5 cells")+
+                    # geom_errorbar(stat="summary", position=position_dodge(width=0.9), width=0.5, aes(fill=Groups))+
+                    geom_hline(data=mock_mean, aes(color=Groups, yintercept = `.`))+
+                    # geom_bar(stat="summary", position="dodge", color="black", aes(fill=Groups))+
+                    geom_boxplot(color="black", aes(fill=Groups), alpha=0.4, outlier.alpha = 0)+
+                    geom_jitter(position=position_jitterdodge(jitter.width = 0.2), shape=21, aes(fill=Groups), size=3)+
+                    # geom_quasirandom(position = position_quasirandom(), shape=21)+
+                    scale_x_discrete(limits=colnames(t_substr)[!colnames(t_substr) %in% c("Mice", "Groups", ctrl, mock)])+
+                    geom_segment(data=t_maps$brackets, aes(x=x1, xend=x2, y=y1, yend=y2), color="black")+
+                    geom_text(data=t_stats, aes(t_maps$label$x, t_maps$label$y, label=p.signif), color="black")+
+                    theme_bw()+
+                    theme(axis.text.x=element_text(angle=45, hjust=1))
+            }else{
+                ggplot(melt(t_substr, id=c("Mice", ctrl, "Groups")), aes(variable, value))+
+                    labs(x="", y="Spots/2.5*10^5 cells")+
+                    # geom_errorbar(stat="summary", position=position_dodge(width=0.9), width=0.5, aes(fill=Groups))+
+                    # geom_bar(stat="summary", position="dodge", color="black", aes(fill=Groups))+
+                    geom_boxplot(color="black", aes(fill=Groups), alpha=0.4, outlier.alpha = 0)+
+                    geom_jitter(position=position_jitterdodge(jitter.width = 0.2), shape=21, aes(fill=Groups), size=3)+
+                    # geom_quasirandom(width=0.2, position=position_dodge(), shape=21)+
+                    scale_x_discrete(limits=colnames(t_substr)[!colnames(t_substr) %in% c("Mice", "Groups", ctrl, mock)])+
+                    geom_segment(data=t_maps$brackets, aes(x=x1, xend=x2, y=y1, yend=y2), color="black")+
+                    geom_text(data=t_stats, aes(t_maps$label$x, t_maps$label$y, label=p.signif), color="black")+
+                    theme_bw()+
+                    theme(axis.text.x=element_text(angle=45, hjust=1))
+            }
+        }
+        
+    })
+    
+    output$flexstats <- renderUI({
+      observeEvent(dades$stats, {})
+      t_stats<-dades$stats
+      if (!is.null(dades$stats)){
+        return(
+          t_stats %>%
+            flextable() %>%
+            theme_vanilla() %>%
+            fontsize(size=14, part="all") %>%
+            padding(padding=10, part="all") %>%
+            color(~ p.adj < 0.05, color = "red")%>%
+            autofit()
+        ) %>%
+          htmltools_value()
+      }
+    })
+    
+    output$downloadData <- downloadHandler(
+        filename = function() {
+            paste("elispot", ".zip", sep="")
+        },
+        content = function(file){
+            print(file)
+            # tempReport <- file.path(tempdir(), "elispots.Rmd")
+            # file.copy("elispots.Rmd", tempReport, overwrite = TRUE)
+            params=list(file=input$file1$datapath, positive=input$positive, showstats=input$showstats, test=input$test, umbral_pos=input$umbral_pos)
+            rmarkdown::render("elispots.Rmd", output_file="Results_elispot.html", params=params, envir = new.env(parent = globalenv()))
+            zip(file,
+                c("Results_elispot.html","data4graphpad.xlsx")
+            )
+            
+        },
+        contentType="application/zip"
+    )
+}
+
+# Run the application 
+shinyApp(ui = ui, server = server)
diff --git a/elispots.Rmd b/elispots.Rmd
new file mode 100644
index 0000000..07b2f70
--- /dev/null
+++ b/elispots.Rmd
@@ -0,0 +1,95 @@
+---
+title: "Elispots"
+output: html_document
+params:
+  file: "../Elispot-validation-spots.xlsx"
+  positive: F
+  showstats: F
+  test: "Ttest"
+  umbral_pos: 10
+---
+
+```{r, echo=F, warning=F, message=F}
+dades<-read_xlsx(params$file)
+
+ctrl<-"Ctrl+"
+mock<-"Mock"
+
+table<-dades[,colnames(dades) != "Groups"]
+t_mean<-dcast(melt(table, id="Mice"),Mice~variable, mean, na.rm=T)
+
+t_substr<-data.frame("Mice"=t_mean[,1],
+                     as.data.frame(t(apply(t_mean[,2:ncol(table)], 1, function(x) x-x[mock])))
+)
+t_mean_group<-merge(t_mean, unique(dades[c("Mice","Groups")]), id="Mice")
+mock_mean<-dcast(t_mean_group, Groups~., value.var=mock, mean)
+mock_mean[,2]<-mock_mean[,2]*2
+mock_mean[mock_mean$. < params$umbral_pos,2]<-params$umbral_pos
+
+t_substr<-merge(t_substr, unique(dades[c("Mice","Groups")]), id="Mice")
+t_substr<-t_substr[,c(1, ncol(t_substr), 2:(ncol(t_substr)-1))]
+t_substr<-t_substr[,c("Mice", "Groups", colnames(t_substr)[!colnames(t_substr) %in% c("Mice", "Groups")])]
+t_substr[,3:ncol(t_substr)]<-apply(t_substr[,3:ncol(t_substr)],2, function(x) replace(x, which(x < 0),0))
+colnames(t_substr)<-c("Mice", "Groups", colnames(t_mean)[2:ncol(t_mean)])
+
+t<-melt(t_substr[,!colnames(t_substr) %in% c(ctrl, mock)])
+if (params$test == "T-test (adj Holm)"){
+    t_stats<-multi_stats(t, "value", "variable", "Groups", stat.test = "ttest")
+}
+if (params$test == "Wilcoxon (adj Holm)"){
+    t_stats<-multi_stats(t, "value", "variable", "Groups", stat.test = "wilcox")
+}
+t_stats<-t_stats %>% filter(p.signif != "ns")
+
+t_maps<-generate_labstats(t_stats, t, "value", "variable", "Groups")
+
+if (params$showstats == F){
+    t_stats<-as.data.frame(matrix(nrow=0, ncol=6))
+    colnames(t_stats)<-c("variable", "group1", "group2", "p.adj", "p.signif", "Method")
+    t_maps<-list()
+    t_maps[["label"]]<-as.data.frame(matrix(nrow = 0, ncol=2))
+    colnames(t_maps$label)<-c("x", "y")
+    t_maps[["brackets"]]<-as.data.frame(matrix(nrow = 0, ncol=4))
+    colnames(t_maps$brackets)<-c("y1", "y2", "x1", "x2")
+}
+
+set.seed(123)
+if (params$positive == T){
+    ggplot(melt(t_substr, id=c("Mice", ctrl, "Groups")), aes(variable, value))+
+        labs(x="", y="Spots/2.5*10^5 cells")+
+        geom_errorbar(stat="summary", position=position_dodge(width=0.9), width=0.5, aes(fill=Groups))+
+        geom_hline(data=mock_mean, aes(color=Groups, yintercept = `.`))+
+        geom_bar(stat="summary", position="dodge", color="black", aes(fill=Groups))+
+        geom_jitter(position=position_jitterdodge(jitter.width = 0.2), shape=21, aes(fill=Groups))+
+        # geom_quasirandom(position = position_quasirandom(), shape=21)+
+        scale_x_discrete(limits=colnames(t_substr)[!colnames(t_substr) %in% c("Mice", "Groups", ctrl, mock)])+
+        geom_segment(data=t_maps$brackets, aes(x=x1, xend=x2, y=y1, yend=y2), color="black")+
+        geom_text(data=t_stats, aes(t_maps$label$x, t_maps$label$y, label=p.signif), color="black")+
+        theme_bw()+
+        theme(axis.text.x=element_text(angle=45, hjust=1))
+}else{
+    ggplot(melt(t_substr, id=c("Mice", ctrl, "Groups")), aes(variable, value))+
+        labs(x="", y="Spots/2.5*10^5 cells")+
+        geom_errorbar(stat="summary", position=position_dodge(width=0.9), width=0.5, aes(fill=Groups))+
+        geom_bar(stat="summary", position="dodge", color="black", aes(fill=Groups))+
+        geom_jitter(position=position_jitterdodge(jitter.width = 0.2), shape=21, aes(fill=Groups))+
+        # geom_quasirandom(width=0.2, position=position_dodge(), shape=21)+
+        scale_x_discrete(limits=colnames(t_substr)[!colnames(t_substr) %in% c("Mice", "Groups", ctrl, mock)])+
+        geom_segment(data=t_maps$brackets, aes(x=x1, xend=x2, y=y1, yend=y2), color="black")+
+        geom_text(data=t_stats, aes(t_maps$label$x, t_maps$label$y, label=p.signif), color="black")+
+        theme_bw()+
+        theme(axis.text.x=element_text(angle=45, hjust=1))
+}
+
+```
+
+```{r, echo=F, warning=F, message=F}
+t_stats %>%
+            flextable() %>%
+            theme_vanilla() %>%
+            fontsize(size=14, part="all") %>%
+            padding(padding=10, part="all") %>%
+            color(~ p.adj < 0.05, color = "red")%>%
+            autofit()
+```
+
diff --git a/funcions.R b/funcions.R
new file mode 100644
index 0000000..7e4b402
--- /dev/null
+++ b/funcions.R
@@ -0,0 +1,137 @@
+library(dplyr)
+library(dunn.test)
+library(tibble)
+library(reshape2)
+library(ggplot2)
+
+corr_multi<-function(table, genes1, genes2, font_size=5, ruta="./", exportar=F)
+{
+  cor.list<-list()
+  cont<-1
+  for (gene1 in genes1){
+    for (gene2 in genes2){
+      print(gene1)
+      print(gene2)
+      cor_temp<-cor.test(table[,gene1], table[,gene2])
+      if (!is.na(cor_temp$p.value)){
+        pval_plot<-paste("cor =",format(cor_temp$estimate, digits=2),if(cor_temp$p.value < 2.2e-16){" , p < 2.2e-16"}else{paste(", p=",format(cor_temp$p.value,digits=3))})
+      }
+      cor.list[[cont]]<-ggplotGrob(ggplot(data=table, aes(x=table[,gene1], y=table[,gene2]))+ 
+                                     geom_point()+
+                                     geom_smooth(method="lm")+
+                                     geom_text(aes(x=min(table[,gene1], na.rm=T), y=(max(table[,gene2],na.rm = T)+max(table[,gene2],na.rm = T)*0.2), label=pval_plot), hjust="inward", family="serif", size=font_size)+
+                                     labs(x=gene1, y=gene2)+
+                                     theme_bw())
+      if (exportar == T){
+        png(paste0(ruta,gene1,"vs",gene2,".png"), res=900, width=3000, height=3000)
+        plot(cor.list[[cont]])
+        dev.off()
+      }
+      cont<-cont+1
+    }
+  }
+  return(cor.list)
+}
+
+multi_stats<-function(table, value.var, x, group, stat.test, adjust="default", paired=F){
+  ## Requires dplyr and tibble packets
+  defaults=c("dunn"="none", "ttest"="holm", "wilcox"="holm")
+  funs<-c("ttest"="pairwise.t.test", "wilcox"="pairwise.wilcox.test")
+  if (adjust == "default"){adjust=defaults[stat.test]}
+  stat.def<-as.data.frame(matrix(nrow=0, ncol=5))
+  colnames(stat.def)<-c(x, "group1", "group2", "p.adj", "p.signif")
+  for (point in unique(table[,x])){
+    condition<-all(table %>% filter(table[,x] == point) %>% pull(value.var) == 0) == F
+    len_group<-length(unique(table %>% filter(table[,x] == point) %>% pull(group)))
+    if (condition == T & !is.na(condition) & len_group > 1){
+      if(stat.test == "dunn"){
+        test<-dunn.test(table %>% filter(table[,x] == point) %>% pull(value.var), table %>% filter(table[,x] == point) %>% pull(group), method=adjust)
+        comp<-strsplit(test$comparisons, " - ")
+        stat.temp<-data.frame(matrix(unlist(comp), nrow=length(comp), byrow=T), "p.adj"=test$P.adjusted, stringsAsFactors = F, check.names = F)
+        colnames(stat.temp)[1:2]<-c("group1", "group2")
+      }else if (stat.test %in% names(funs)){
+        test<-get(funs[stat.test])(table %>% filter(table[,x] == point) %>% pull(value.var), table %>% filter(table[,x] == point) %>% pull(group), method=adjust, paired=paired)$p.value
+        stat.temp<-melt(test)
+        colnames(stat.temp)<-c("group1", "group2","p.adj")
+      }
+      
+      stat.temp["p.signif"]<-case_when(
+        stat.temp$p.adj >= 0.05 ~ "ns",
+        stat.temp$p.adj < 0.0001 ~ "****",
+        stat.temp$p.adj < 0.001 ~ "***",
+        stat.temp$p.adj < 0.01 ~ "**",
+        stat.temp$p.adj < 0.05 ~ "*"
+      )
+      stat.temp<-stat.temp %>% add_column(x=point, .before=T)
+      colnames(stat.temp)[1]<-x
+      stat.def<-rbind(stat.def, stat.temp)
+    }
+  }
+  stat.def["Method"]<-stat.test
+  return(stat.def)
+}
+
+generate_labstats<-function(table_stat, table, value.var, x, group, y="max", bracket.offset=0.05, bracket.length=0.02){
+  table[,group]<-as.factor(table[,group])
+  table[,x]<-as.factor(table[,x])
+  se<-function(x, na.rm=F) sd(x, na.rm = na.rm)/sqrt(length(x))
+  if (y == "max"){ 
+    formula<-as.formula(paste0(colnames(table_stat)[1], "~."))
+    agg<-dcast(table, formula, value.var = value.var, fun.aggregate = max, na.rm=T)
+  }else if (y == "mean"){ 
+    formula<-as.formula(paste0(colnames(table_stat)[1], "~", group))
+    agg<-dcast(table, formula, value.var = value.var, fun.aggregate = mean, na.rm=T)
+    agg<- data.frame(x=agg[,1], "."=apply(agg[,2:ncol(agg)], 1, max, na.rm=T))
+    colnames(agg)[1]<-x
+  }else if (y == "mean+sd"){ 
+    formula<-as.formula(paste0(colnames(table_stat)[1], "~", group))
+    agg<- dcast(table, formula, value.var = value.var, fun.aggregate = function(x) mean(x,na.rm=T)+sd(x,na.rm=T))
+    agg<- data.frame(x=agg[,1], "."=apply(agg[,2:ncol(agg)], 1, max, na.rm=T))
+    colnames(agg)[1]<-x
+  }else if (y == "mean+se"){ 
+    formula<-as.formula(paste0(colnames(table_stat)[1], "~", group))
+    agg<- dcast(table, formula, value.var = value.var, fun.aggregate = function(x) mean(x,na.rm=T)+se(x,na.rm=T))
+    agg<- data.frame(agg[,1], "."=apply(agg[,2:ncol(agg)], 1, max, na.rm=T))
+    colnames(agg)[1]<-x
+  }
+  
+  t<-data.frame("y1"=merge(table_stat, agg ,sort=F)[,"."]+diff(range(table[value.var], na.rm = T))*bracket.offset,
+                "y2"=merge(table_stat, agg ,sort=F)[,"."]+diff(range(table[value.var], na.rm = T))*bracket.offset,
+                "x1"= match(table_stat[,x], unique(table[,x]))+
+                  0.75*((match(table_stat$group1, levels(table[,group]))-0.5)/length(levels(table[,group]))-0.5),
+                "x2"= match(table_stat[,x], unique(table[,x]))+
+                  0.75*((match(table_stat$group2, levels(table[,group]))-0.5)/length(levels(table[,group]))-0.5)
+  )
+  for (dia in unique(table_stat[,1])){
+    t[table_stat[,x] == dia,"y1"]<-seq(t[table_stat[,x] == dia,"y1"][1],
+                                           t[table_stat[,x] == dia,"y1"][1]+diff(range(table[,value.var], na.rm = T))*0.05*(nrow(table_stat[table_stat[,x] == dia,])-1), 
+                                           by=diff(range(table[,value.var], na.rm = T))*0.05)
+    t[table_stat[,x] == dia,"y2"]<-t[table_stat[,x] == dia,"y1"]
+  }
+  
+  t_def<-as.data.frame(matrix(ncol=4, nrow=0))
+  for (row in 1:nrow(t)){
+    t_def<-rbind(t_def, t[row,],
+                 c(t[row,"y1"]-diff(range(table[,value.var], na.rm = T))*bracket.length, t[row,"y1"], t[row,"x1"], t[row,"x1"]),
+                 c(t[row,"y1"]-diff(range(table[,value.var], na.rm = T))*bracket.length, t[row,"y1"], t[row,"x2"], t[row,"x2"]))
+  }
+  t_lab<-data.frame("x"=t$x1+(t$x2-t$x1)/2, "y"=t$y1+diff(range(table[,value.var], na.rm = T))*0.005, check.names = F)
+  return(list("label"=t_lab, "brackets"=t_def))
+}
+
+secfile<-function(file){
+  ext<-strsplit(file, ".", fixed = T)[[1]]
+  ext<-ext[length(ext)]
+  num<-1
+  while(file.exists(file) == T){
+    if(num == 1){
+      file_tmp<-strsplit(file, ".", fixed=T)[[1]]
+      file<-paste0(paste(file_tmp[1:(length(file_tmp)-1)], collapse = "."),"_",num,".",file_tmp[length(file_tmp)])
+    }else{
+      file_tmp<-paste(strsplit(file, "_", fixed=T)[[1]][-length(strsplit(file, "_", fixed=T)[[1]])], collapse = "_")
+      file<-paste0(file_tmp, "_", num,".",ext)
+    }
+    num<-num+1
+  }
+  return(file)
+}
\ No newline at end of file