################################################################################################## # Author : E. Gerrits # Date : 11-dec-2018 # Dataset : Microglia from NAWM, NAGM and LESIONS from MS donors # Purpose : Make an application to visualize the data and make it more accessible for others # Updated on 06-01-2019. Run by " shiny::runApp() " ################################################################################################## setwd("/media/astrid/My Passport/Analysis6_SP2-16_NEW/Shiny/") # CHANGE Packages <- c("shiny", "shinythemes", "dplyr", "readr", "ggplot2", "plotly", "data.table", "feather", "DT") no <- 14 # install packages (if not already installed) inst <- Packages %in% installed.packages() if(length(Packages[!inst]) > 0) install.packages(Packages[!inst], repos='http://cran.us.r-project.org') # load packages into session lapply(Packages, require, character.only=TRUE) # load donor info table # ik heb er een tabelletje met info over de doners in gezet # load datasets datafile <- feather("seu.all.UMAP_datafile.feather") umapdata <- read_feather("seu.all.UMAP_umapdata.feather") pcadata <- read_feather("seu.all.UMAP_pcadata.feather") barplot <- read.csv("seu.all.UMAP_barplot.csv", row.names = 1) markers <- read.csv("Markers.csv", row.names = 1) #-------------------------------------------------------------------------------------------------# #-------------------------------------------------------------------------------------------------# #-------------------------------------------------------------------------------------------------# # define UI ui <- fluidPage(theme = shinytheme("simplex"), # lumen is ook mooi titlePanel("Single-cell RNA-sequencing of microglia from nondemented elderly and AD donors"), sidebarLayout( sidebarPanel( h4("Select variables for the plots"), # select dataset # in dit geval is elke donor een andere dataset, maar jij hebt maar 1 dataset dus misschien moet je dit veranderen # selectInput("donor", "Choose donor", c("2018-115" ,"2018-116", "2018-127"), selected = "2018-115"), # select how to group and color the data # zeg maar de "group.by" die je bij seurat gebruikt selectInput("groups", "Group/color by", choices = c("DonorGroup", "Cluster", "DonorID", "Superpool", "InitialDiagnosis", "expression"), selected = "DonorGroup"), # select gene to plot textInput("ingene", "Enter the name of a gene (in capitals)", value = "ITGAX"), # more info hr(), print("Developed by Emma Gerrits, Last updated: 18-04-2019 by Astrid Alsema") ), mainPanel( tabsetPanel(tabPanel("Information", # eerste tablad gebruik ik voor algemene info van het experiment hr(), # een lijn print("Microglia were mechanically isolated from post-mortem brain tissue and FACS-sorted in plates. An adapted, barcoded version of SmartSeq2 was used for RNA-seq library construction.") ), tabPanel("PCA", plotlyOutput("pca", height = "600px"), # plot output van de output$pca variabele die je in de server definieert hr(), print("This plot shows a linear dimensionality reduction of the data. The principal component analysis was computed on the highly variable genes. Each dot is a cell and the distances between the cells indicate how different these cells are based on their gene expression. You can choose to color the dots by DonorGroup, cluster or gene expression. Also, you can hover over the plot to get specific information about each cell, and you can zoom in by scrolling with your mouse.") ), tabPanel("UMAP", plotOutput("umap", height = "600px"), hr(), print("This plot shows a non-linear dimensionality reduction of the data. Each dot is a cell, but since this method is non-linear the distances between the dots do NOT indicate how different the cells are. This plot is very useful to visualize the clustering of the cells. In order to generate this clustering, the effects of number of UMIs, percentage of mitochondrial and ribosomal genes were removed. The dimensionality reduction was calculated on the first 10 PCs.") ), # tabPanel("Barplot", # plotlyOutput("barplot", # height = "600px"), # hr(), # print("This plot shows the relative presence of each cluster per DonorGrup You can hover over the plot to see the percentages.") # ), tabPanel("Cluster markers", DT::dataTableOutput("markers", height = "600px"), hr(), print("Differentially expressed genes per cluster compared to all other clusters. These can be seen as cluster markers. 'pct.1' is the percentage of cells expressing the gene in the specific cluster, 'pct.2' is the percentage of cells expressing that gene in all the other clusters together.") ), tabPanel("LogFC plot", plotOutput("logFCplot", height = "600px"), hr(), print("DE analysis was performed for each cluster compared to all other cells together. In this plot you see the logFC values and the colorscale represents the percentage of cells in a cluster that expresses the gene. ") ), tabPanel("Violin plot", print(tags$h2("Violin Plot")), plotOutput("violins", height = "600px"), hr(), print("In this plot you visualize the normalized expression values of a selected gene. Each dot is a cell and the violin represents a density distribution. The horizontal position of the dots is random and changes everytime you make the plot."), print(tags$h2("Violin Plot2")), plotOutput("violins2", height = "600px"), hr(), hr(), print(tags$h2("Density plot")), plotOutput("density", height = "600px"), hr(), print("This plot shows the distribution of the gene expression within a group. On the x-axis you see the expression values, and on the y-axis you see the density. So the heigher the line the more cells have that specific expression value in that group. Think of the plot as a variation of a histogram where the distribution is smoothened.") ), tabPanel("Test", DT::dataTableOutput("test") ) )# end tabsetPanel )# end MainPanel ) # end SidebarLayout ) # end fluidPage #-------------------------------------------------------------------------------------------------# #-------------------------------------------------------------------------------------------------# #-------------------------------------------------------------------------------------------------# server <- function(input, output, session) { # define colorscheme cols <- reactive( if (input$groups == "DonorGroup" || input$groups == "expression"){ # default to DonorGroup colors c("limegreen", "navyblue", "red3") } else if (input$groups == "Cluster"){ # dependent on number of clusters which is dependent on donor #c("#FFE167", "#86CCB3", "#FBA080", "#A0AED2", "#E9AECB", "#B4E178", "#E8CDA5") c("#A0AED2", "#FFE167", "#86CCB3", "#FBA080", "#E9AECB", "#B4E178", "#E8CDA5", "#FFE167", "#8567FF", "darkgoldenrod", "azure4", "cadetblue","bisque1","darkorchid3","aquamarine") } else if (input$groups == "expression"){ NULL } ) # define dataset based on donor # aangezien jij denk ik maar 1 dataset hebt (denk ik), kun je dit stuk heel erg verkorten, want je laad altijd dezelfde datafiles # dus dan kun je all ifelsjes vermijden datafile <- reactive({ datafile} ) pcadata <- reactive({ pcadata} ) umapdata <- reactive({ umapdata} ) barplotdata <- reactive({ barplot} ) Markers <- reactive({ markers} ) #-------------------------------------------DONOR INFO--------------------------------------------# #-----------------------------------------------PCA----------------------------------------------# # make dataframe PCAdata <- reactive({ if (input$groups == "Cluster" || input$groups == "DonorGroup"){ as.data.frame(cbind(pcadata()[,1:2], pcadata()[, input$groups])) } else if (input$groups == "expression"){ validate(need(try(input$ingene %in% colnames(datafile())), "Gene not detectable in single-cell data...")) as.data.frame(cbind(pcadata()[,1:2], as.numeric(as.character(unlist(datafile()[,input$ingene]))), pcadata()[,"DonorGroup"], pcadata()[,"Cluster"])) } }) # output$test <- DT::renderDataTable({ # logFCdata() # }) # # make the plot output$pca <- renderPlotly({ if (input$groups == "DonorGroup"){ plot_ly(x = PCAdata()[,1], y = PCAdata()[,2], mode = "markers", color = ~factor(PCAdata()[,3]) ) } else if (input$groups == "Cluster"){ plot_ly(x = PCAdata()[,1], y = PCAdata()[,2], mode = "markers", color = ~factor(PCAdata()[,3]) ) } else if (input$groups == "expression"){ plot_ly(x = PCAdata()[,1], y = PCAdata()[,2], mode = "markers", colors = colorRamp(colors = c("lightgrey", "lemonchiffon1", "tan1", "firebrick3")) ) # size = 0.7, # alpha = 1, # hoverinfo = 'text' # text = ~paste("DonorGroup: ", PCAdata()[,5], # "\n Cluster: ", PCAdata()[,6], # "\n LogTPM: ", format(round(PCAdata()[,4], 2), nsmall = 2)) # ) %>% # layout(scene = list( # xaxis = list(title = "PC1", showspikes = FALSE), # yaxis = list(title = "PC2", showspikes = FALSE) # ) # ) %>% # config(displayModeBar = FALSE) } }) #----------------------------------------------UMAP-----------------------------------------------# # make dataframe UMAPdata <- reactive({ if (input$groups == "Cluster" || input$groups == "DonorGroup"){ as.data.frame(cbind(umapdata()[,1:2], umapdata()[,input$groups])) } else if (input$groups == "expression"){ validate(need(try(input$ingene %in% colnames(datafile())), "Gene not found...")) as.data.frame(cbind(umapdata()[,1:2], as.numeric(as.character(unlist(datafile()[,input$ingene]))))) } }) # make the plot output$umap <- renderPlot({ if (input$groups == "Cluster" || input$groups == "DonorGroup"){ ggplot(UMAPdata(), aes(x = UMAPdata()[,1], y = UMAPdata()[,2], colour = as.factor(UMAPdata()[,3]))) + geom_point(size = 0.7) + scale_colour_manual(values = cols()) + xlab("UMAP1") + ylab("UMAP2") + labs(colour = input$groups, size = 11) + theme_classic() + theme(legend.text=element_text(size=rel(1.1))) + guides(colour = guide_legend(override.aes = list(size=6))) } else if (input$groups == "expression"){ ggplot(UMAPdata(), aes(x = UMAPdata()[,1], y = UMAPdata()[,2])) + geom_point(aes(colour = UMAPdata()[,3]), size = 0.7) + scale_colour_gradientn(colours = c("lightgrey", "lemonchiffon1", "tan1", "firebrick3")) + xlab("UMAP1") + ylab("UMAP2") + labs(colour = "LogTPM", size = 11) + theme_classic() } }) # #---------------------------------------------BARPLOT---------------------------------------------# # output$barplot <- renderPlotly({ # if (input$donor == "2018-115" || input$donor == "2018-116"){ # plot_ly(barplotdata(), # x = ~factor(barplotdata()[,2], levels = c("NAGM", "NAWM", "LESIONS")), # y = ~barplotdata()[,3], # type = 'bar', # name = ~barplotdata()[,1], # color = ~as.factor(barplotdata()[,1]), # colors = c("#A0AED2", "#FFE167", "#86CCB3", "#FBA080", "#E9AECB", "#B4E178", "#E8CDA5"), # #colors = c("#F8766D", "#CD9600", "#7CAE00", "#00BE67", "#00BFC4", "#00A9FF", "#C77CFF"), # hoverinfo = 'text', # text = ~paste('Region: ', barplotdata()[,2], # '\n Cluster: ', barplotdata()[,1], # '\n Percentage: ', format(round(barplotdata()[,3], 2), nsmall = 2), "%") # ) %>% # layout(xaxis = list(title = ""), # yaxis = list(title = 'Percentage (%)'), # barmode = 'stack' # ) %>% # config(displayModeBar = FALSE) # } else if (input$donor == "2018-127"){ # plot_ly(barplotdata(), # x = ~factor(barplotdata()[,2], levels = c("NAGM", "NAWM", "LESIONS")), # y = ~barplotdata()[,3], # type = 'bar', # name = ~barplotdata()[,1], # color = ~as.factor(barplotdata()[,1]), # colors = c("#86CCB3", "#FBA080", "#A0AED2"), # hoverinfo = 'text', # text = ~paste('Region: ', barplotdata()[,2], # '\n Cluster: ', barplotdata()[,1], # '\n Percentage: ', format(round(barplotdata()[,3], 2), nsmall = 2), "%") # ) %>% # layout(xaxis = list(title = ""), # yaxis = list(title = 'Percentage (%)'), # barmode = 'stack' # ) %>% # config(displayModeBar = FALSE) # } # }) #-------------------------------------------CLUSTERMARKERS----------------------------------------# output$markers = DT::renderDataTable({ t <- as.data.frame(Markers()) t <- t[,-c(7)] d <- t[,c(6,1,2,3,4,5)] d }) #----------------------------------------------LOGFCPLOT------------------------------------------# logFCdata <- reactive({ validate(need(try(input$ingene %in% Markers()[,7]), "Gene not found or not significantly differentially expressed between clusters..")) dat <- Markers()[Markers()[,7] == as.character(input$ingene), ] all_clusters <- unique(Markers()[,6]) d <- dat[,6] toAdd <- setdiff(all_clusters, d) Add <- data.frame(matrix(NA, nrow = length(toAdd), ncol = ncol(Markers()))) for (i in 1:length(toAdd)){ new <- c(NA, 0, 0, NA, NA, toAdd[i], input$ingene) Add[i,] <- new } colnames(Add) <- colnames(Markers()) dat <- rbind(dat, Add) dat[,6] <- factor(as.character(dat[,6]), levels = seq(from = 0, to = no)) dat[,2] <- as.numeric(dat[,2]) dat[,3] <- as.numeric(dat[,3]) dat <- dat[order(dat[,6]),] dat }) output$logFCplot <- renderPlot({ ggplot(logFCdata(), aes(x = logFCdata()[,6], y = logFCdata()[,2], fill = logFCdata()[,3])) + geom_bar(colour = "white", stat = "identity") + ylim(-3,3) + xlab("Cluster") + ylab("Average logFC") + ggtitle(input$ingene) + geom_hline(yintercept = 0, color = "black") + scale_fill_gradientn(colours = c("lightgrey", "lemonchiffon1", "tan1", "firebrick3"), name = paste(" ratio of cells \n where", input$ingene, "\n is expressed \n in the cluster"), space = "lab", guide = "colourbar", aesthetics = "fill", limits = c(0,1)) + theme_classic() }) #----------------------------------------------VIOLINS--------------------------------------------# # subset data subsetdata <- reactive({ validate(need(try(input$ingene %in% colnames(datafile())), "Gene not detectable in single-cell data...")) as.data.frame(unlist(datafile()[,input$ingene])) }) # add metadata subsetwithmeta <- reactive({ validate(need(try(input$ingene %in% colnames(datafile())), "Gene not detectable in single-cell data...")) if (input$groups == "expression"){ Sub <- as.data.frame(cbind(as.numeric(as.character(subsetdata()[,1])), factor(pcadata[,5]))) # place storing DonorGroup } else if (input$groups == "Cluster"){ Sub <- as.data.frame(cbind(as.numeric(as.character(subsetdata()[,1])), factor(pcadata[,4]))) # place storing Cluster } } ) output$test <- DT::renderDataTable({ head(subsetwithmeta()) # subsetwithmeta() is null }) # make the violin plot output$violins <- renderPlot({ ggplot(subsetwithmeta(), aes(x = subsetwithmeta()[,2], y = subsetwithmeta()[,1], fill = as.factor(subsetwithmeta()[,2]))) + scale_fill_manual(values = cols()) + geom_violin(aes(x = as.factor(subsetwithmeta()[,2]), y = subsetwithmeta()[,1]), scale = "width", show.legend = FALSE) + geom_jitter(size = 0.1, show.legend = FALSE) + xlab("") + ylab("cpm") + theme_classic() }) #version 2 output$violins2 <- renderPlot({ggplot(subsetwithmeta(), aes(x = subsetwithmeta()[,2], y = subsetwithmeta()[,1], fill = as.factor(subsetwithmeta()[,2]))) + geom_violin(scale="width", trim = FALSE) + scale_fill_manual(values=c("limegreen", "navyblue", "red3")) + guides(fill=FALSE, color = FALSE) + geom_boxplot(width=0.2, aes(group=subsetwithmeta()[,2]), fill = NA, outlier.shape = NA) + geom_jitter(color = "black", size = 1, alpha = 0.1) + #aes(group=interaction(as.factor(dataset[[xCol]]), as.factor(dataset[,colCondition]))), #position = position_dodge(width = .9)) #+ labs(title = paste("Expression of", as.character(input$ingene))) + ylab("log counts - Seurat2 ") }) # #-------------------------------------------DENSITY PLOT-----------------------------------------# output$density <- renderPlot({ if (input$groups == "DonorGroup" || input$groups == "expression"){ ggplot(subsetwithmeta(), aes(x = subsetwithmeta()[,1], colour = as.factor(subsetwithmeta()[,2]), fill = as.factor(subsetwithmeta()[,2]))) + geom_density(size = 0.8, alpha = 0.05, trim = TRUE) + scale_colour_manual(values = cols()) + scale_fill_manual(values = cols()) + xlab("cpm") + ylab("Density") + theme_classic() + guides(fill=FALSE) + labs(colour = "DonorGroup") + theme(legend.text=element_text(size=rel(1.1))) } else if (input$groups == "Cluster"){ ggplot(subsetwithmeta(), aes(x = subsetwithmeta()[,1], colour = as.factor(subsetwithmeta()[,2]), fill = as.factor(subsetwithmeta()[,2]))) + geom_density(size = 0.8, alpha = 0.05, trim = TRUE) + scale_colour_manual(values = cols()) + scale_fill_manual(values = cols()) + xlab("cpm") + ylab("Density") + theme_classic() + guides(fill=FALSE) + labs(colour = input$groups) + theme(legend.text=element_text(size=rel(1.1))) } }) #--------------------------------------------CLOSE APP------------------------------------------# session$onSessionEnded(function() { # dit stuk code zegt eigenlijk: sluit de app ook echt af wanneer je het scherm weg klikt, is vooral nuttig als je de app gaat delen stopApp() }) } #-------------------------------------------------------------------------------------------------# #-------------------------------------------------------------------------------------------------# #-------------------------------------------------------------------------------------------------# # create shiny app shinyApp(ui = ui, server = server) # dit maakt daadwerkelijk een app van alle code die hier boven staat