443 lines
21 KiB
R
443 lines
21 KiB
R
##################################################################################################
|
|
# 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 |