Commit 8ec80231 authored by avandenheuvel1's avatar avandenheuvel1

Updated Ref in app.R

parent 49a252b0
################################### START ######################################
####################### Get Packages ##########################
lib.location <- "Packages"
# get packages
library("ggplot2", lib.loc = lib.location)
library("reshape2", lib.loc = lib.location)
library("irlba", lib.loc = lib.location)
library("monocle", lib.loc = lib.location) # http://bioconductor.org/packages/release/bioc/html/monocle.html
library("cellrangerRkit", lib.loc = lib.location)
# source("https://s3-us-west-2.amazonaws.com/10x.files/supp/cell-exp/rkit-install-1.1.0.R") # https://github.com/hb-gitified/cellrangerRkit.git
library("shiny", lib.loc = lib.location)
source("https://raw.githubusercontent.com/obigriffith/biostar-tutorials/master/Heatmaps/heatmap.3.R")
#
# library("Rmisc")
# https://github.com/hb-gitified/cellrangerRkit/blob/master/scripts/rkit-install-1.1.0.R
# library("cellrangerRkit", lib.loc = "Package/")
# install.packages("Package/cellrangerRkit-1.1.0.tar.gz", type="source", repos=NULL)
# install_github('hb-gitified/cellrangerRkit')
# library("stringr")
# library("colorspace")
# library("RBGL")
# library("dplyr")
# library("rhdf5")
# library("stringi")
# library("VGAM")
# library("irlba")
# library("DDRTree")
# library("slam")
# library("combinat")
# library("limma")
# library("jsonlite")
# library("gridExtra")
# library("grid")
# library("RGraphics")
# library("digest")
# library("bit64")
# library("shiny")
# library("BiocGenerics")
# library("Biobase")
# library("FNN")
####################### Get Packages ##########################
################################### END ######################################
################################### START ######################################
####################### Get Data and Settings ##########################
################### get data ready
### Get data PrimaryData
Primary.dataname <- load("PrimaryData/TrajectoryPlot.Rdata")
Primary.data <- get(Primary.dataname)
pData(Primary.data)$SampleID <- colsplit(pData(Primary.data)$barcode, pattern = "-", names = c("a", "b"))$b
gbm.Pr <- load_cellranger_matrix("PrimaryData")
fData(gbm.Pr)$gene_short_name <- fData(gbm.Pr)$symbol
analysis_results.original.Pr <- load_cellranger_analysis_results("PrimaryData")
Table.data.Pr <- fData(Primary.data)[,c("id", "gene_short_name", "num_cells_expressed", "Pooled_mRNAs", "myotube.expression", "FSHDpop.expression", "Fold.change")]
################### get data ready
####################### Get Data and Settings ##########################
################################### END ######################################
################################### START ######################################
#################### Defining new Functions #######################
##### function --- PrimaryData analysis
prepping.data.for.plotting.Primary <- function(x,
num.cells.state1 = c(1,25), num.cells.state2 = c(1,25), num.cells.state3 = c(1,25),
geneID.or.gene_short_name = "gene_short_name", genes = geneset,
order.on1 = "Pseudotime", order.on2 = NA, order.on3 = NA) {
##### data.processing ##########################################################
State1branch <- x[,pData(x)$State == 1]
State1branch <- State1branch[,order(pData(State1branch)$Pseudotime)]
State1branch.cells.to.inlcude <- State1branch[, rownames(pData(State1branch)[c(num.cells.state1[1]:num.cells.state1[2]),])]
State2branch <- x[,pData(x)$State == 2]
State2branch <- State2branch[,order(pData(State2branch)$Pseudotime)]
State2branch.cells.to.inlcude <- State2branch[, rownames(pData(State2branch)[c(num.cells.state2[1]:num.cells.state2[2]),])]
State3branch <- x[,pData(x)$State == 3]
State3branch <- State3branch[,order(pData(State3branch)$Pseudotime)]
State3branch.cells.to.inlcude <- State3branch[, rownames(pData(State3branch)[c(num.cells.state3[1]:num.cells.state3[2]),])]
data.to.plot <<- x[,pData(x)$barcode %in% pData(State1branch.cells.to.inlcude)$barcode | pData(x)$barcode %in% pData(State2branch.cells.to.inlcude)$barcode | pData(x)$barcode %in% pData(State3branch.cells.to.inlcude)$barcode]
# select only wanted genes
my_genes <- row.names(subset(fData(data.to.plot),
fData(data.to.plot)[,paste(geneID.or.gene_short_name)] %in% genes[,1]))
data.to.plot.genefiltered <- data.to.plot[fData(data.to.plot)$id %in% my_genes,]
# Reorder this data.subset for Pseudotime order
ordering.list <- c(pData(data.to.plot.genefiltered)[paste(order.on1)])
if(!is.na(order.on3)){
index <- with(exprs(data.to.plot.genefiltered), order(pData(data.to.plot.genefiltered)[paste(order.on1)], pData(data.to.plot.genefiltered)[paste(order.on2)], pData(data.to.plot.genefiltered)[paste(order.on3)]))}
else if(!is.na(order.on2)){
index <- with(exprs(data.to.plot.genefiltered), order(pData(data.to.plot.genefiltered)[paste(order.on1)], pData(data.to.plot.genefiltered)[paste(order.on2)]))}
else
index <- with(exprs(data.to.plot.genefiltered), order(pData(data.to.plot.genefiltered)[paste(order.on1)]))
data.to.plot.genefiltered.reordered <<- data.to.plot.genefiltered[,index]
rownames(data.to.plot.genefiltered.reordered) <<- fData(data.to.plot.genefiltered.reordered)$gene_short_name
#reorder genes
if(geneID.or.gene_short_name == "gene_short_name"){
genes.included <- subset(genes, genes[,1] %in% fData(data.to.plot.genefiltered.reordered)$gene_short_name)
rownames(data.to.plot.genefiltered.reordered) <<- fData(data.to.plot.genefiltered.reordered)$gene_short_name
row.index <- with(exprs(data.to.plot.genefiltered.reordered), genes.included[,1])
data.to.plot.genefiltered.reordered <<- data.to.plot.genefiltered.reordered[row.index,]}
if(geneID.or.gene_short_name == "id"){
genes.included <- subset(genes, genes[,1] %in% fData(data.to.plot.genefiltered.reordered)$id)
rownames(data.to.plot.genefiltered.reordered) <<- fData(data.to.plot.genefiltered.reordered)$id
row.index <- with(exprs(data.to.plot.genefiltered.reordered), genes.included[,1])
data.to.plot.genefiltered.reordered <<- data.to.plot.genefiltered.reordered[row.index,]}
# Make a matrix & Log2.transformed matrix for the final heatmap
data.to.plot.matrix <<- as.matrix(exprs(data.to.plot.genefiltered.reordered))
data.to.plot.matrix.log2 <<- log((data.to.plot.matrix + 0.1), base = 2)
##### data.processing ##########################################################
##### legend.info ##############################################################
# Make legends for the heatmap
# CellType legend
CellType.colors <- c()
CellType.table <- data.frame("CellType" = c("Unknown","FSHDpop","Healthy.compare.cells"), "CellType.colors" = c("blue", "darkred", "green"))
for (i in 1:length(pData(data.to.plot.genefiltered.reordered)$CellType)){
for(j in 1:length(CellType.table$CellType)){
if(pData(data.to.plot.genefiltered.reordered)$CellType[i] == CellType.table$CellType[j]) {CellType.colors[i] <- as.character(CellType.table$CellType.colors[j])}
}
}
CellType.table <<- CellType.table
# State legend
State.colors <- c()
State.table <- data.frame("State.nr" = c(1,2,3), "State.colors" = c("dodgerblue", "indianred3", "darkgreen"))
for (i in 1:length(pData(data.to.plot.genefiltered.reordered)$State)){
for(j in 1:length(State.table$State.nr)){
if(pData(data.to.plot.genefiltered.reordered)$State[i] == State.table$State.nr[j]) {State.colors[i] <- as.character(State.table$State.colors[j])}
}
}
State.table <<- State.table
# SampleID legend
# first add SampleID column to data
SampleID.color <- c()
SampleID.color.renaming.table <- data.frame("SampleID" = c(1,2,3,4,5,6), "Sample.color" = c("red","green","pink","darkgreen","lightblue","blue"))
for (i in 1:length(pData(data.to.plot.genefiltered.reordered)$SampleID)){
for(j in 1:length(SampleID.color.renaming.table$SampleID)){
if(pData(data.to.plot.genefiltered.reordered)$SampleID[i] == SampleID.color.renaming.table$SampleID[j]) {SampleID.color[i] <- as.character(SampleID.color.renaming.table$Sample.color[j])}
}
}
SampleID.color.renaming.table <<- SampleID.color.renaming.table
color.gradient <- function(x, colors=c("darkblue", "white"), colsteps=75) {
return( colorRampPalette(colors) (colsteps) [ findInterval(x, seq(min(x),max(x), length.out=colsteps)) ] )
}
Pseudotimecol <- color.gradient(x = pData(data.to.plot.genefiltered.reordered)$Pseudotime)
# legends table
clab <- cbind(CellType.colors,State.colors, SampleID.color, Pseudotimecol)
colnames(clab)=c("CellType","State", "SampleID", "Pseudotime")
clab <<- clab
##### legend.info ##############################################################
}
#################### Defining new Functions #######################
################################### END ######################################
################################### START ######################################
#################### Define UI #######################
#### Define UI settings
ui <- fluidPage(
tags$head(tags$style(HTML("label {font-size:75%;}"))),
titlePanel("SingleCellData Analysis"),
sidebarLayout(sidebarPanel(selectInput(inputId = "GeneName.or.GeneID", label = "Used 'gene_short_names' (eg DUX4) or 'geneID' (eg 'ENSG00000260596') ", choices = c("gene_short_name", "id"), selected = c("gene_short_name"),multiple = FALSE, selectize = TRUE),
textAreaInput(inputId = "Genes", label = "Genes.to.plot", value = c("MYOG")),
p(span(strong("Select which cells to visualize (per State)"))),
sliderInput(inputId = "Pr.State1", label = "1.Pr.State1", value = c(1,length(subset(pData(Primary.data)$State, pData(Primary.data)$State == 1))), min =0, max = length(subset(pData(Primary.data)$State, pData(Primary.data)$State == 1)),step = 1),
sliderInput(inputId = "Pr.State2", label = "2.Pr.State2", value = c(1,length(subset(pData(Primary.data)$State, pData(Primary.data)$State == 2))), min =0, max = length(subset(pData(Primary.data)$State, pData(Primary.data)$State == 2)),step = 1),
sliderInput(inputId = "Pr.State3", label = "3.Pr.State3", value = c(1,length(subset(pData(Primary.data)$State, pData(Primary.data)$State == 3))), min =0, max = length(subset(pData(Primary.data)$State, pData(Primary.data)$State == 3)),step = 1),
p(span(strong("Select additional settings"))),
sliderInput(inputId = "Pr.limit.settings", label = "Heatmap color limits", value = c(5), min =1, max = 50,step = 1),
selectInput(inputId = "Pr.GeneOrder",width = "100%", label = "Cluster genes or not (no cluster-tree shown)", choices = c("TRUE", "FALSE"), selected = c("TRUE"),multiple = FALSE),
selectInput(inputId = "Pr.CellOrder",width = "100%", label = "re-order cells in heatmap based on: (sequential)", choices = c("State", "Pseudotime", "CellType","SampleID", "GeneExprs.Gene1"), selected = c("State", "Pseudotime"),multiple = TRUE, selectize = TRUE),
width = 3),
mainPanel(tabsetPanel(
tabPanel("Tool Info",
h1("General info on the use of this tool"),
br(),
p("This ShinyApp is specifically designed to visualize and browse the scRNAseq data of the 'Primary mononuclear myocyte' samples generated by the",span(strong("Van der Maarel-lab, LUMC, Dep. Human Genetics, Leiden, The Netherlands:")), "This data has been published in 'A van den Heuvel. et al. Human Molecular Genetics Nov 16, 2018'"),
p("Our scRNAseq data is analyzed to generate a Pseudotime Trajectory. For this, using each cell's own transcripome, we re-ordered all cells in the sample in such a way that they form a Branched Trajectory-like structure", span(strong("as depicted below"))),
p("This tool allows you to check the expression of any of your gene(s)-of-interest in the scRNAseq in general, but also to check what happens with this expression when a cell would progress in pseudotime."),
p("On the left you find some interactive options with which you can select which cells you want to plot, which genes you want to analyze, etc."),
br(),
h4("In the 'Heatmap'-tab you find:"),
p("A heatmap for the gene expression per cell", span(em("Cells ordered on the Horizontal axis, genes ordered on the vertical axis", style = "color:grey")),
br(),
"With the 're-order cells in heatmap based on:'-option you can re-order cells based on different criteria. You can add mutliple sequential criteria here.",
br(),
"With the 'Cluster genes or not'-option you can select whether the genes (rows) in the heatmap will be clustered based on expression pattern or not.",
br(),
"The heatmap is row-normalized"),
br(),
h4("In the 'Trajectory'-tab you find:"),
p("A plot as depicted below, but then highlighting in", span(strong("blue"),style = "color:blue"), "the cells included in the heatmap."),
br(),
h4("In the 'Expr.in.Trajectory'-tab you find:"),
p("A plot as depicted below, but then depicting the expression level of your gene(s)-of-interest by a color-gradient.",
br(),
em("!! Keep in mind that the gradient is set for all genes you visualize. So one gene with very high expression levels will skewing the gradient for others", style = "color:grey"),
br(),
em("!! Do not insert too many genes when using this Tab...", style = "color:grey")),
br(),
h4("In the 'Gene.Expression.Table'-tab you find:"),
p("An interactive table per selected geneset-of-interest, including info on Read Counts, number of cells expressing the gene, foldChanges for the DUX4-affected cells versus non-affected stage-matched myocytes, etc."),
br(),
h4("In the 'Expression.in.tSNE'-tab you find:"),
p("A tSNE plot including all cells of the siz samples tested ('Aggr Dataset in paper'), clustered by the CellRanger analysis pipeline of 10XGenomics.",
"In this tSNE plot, you can plot the expression of any gene of interest by color-gradient. Using the 'Expression.limits'-option you can change the color-gradient max limit.",
em("If the max limit is set higher than the actual maximum value in the data, it will not rescale any higher", style = "color:grey")),
fluidRow(
h3("Example of Pseudotime trajectory",
br(),
em("Primary Data Trajectory", style = "color:grey")),
column(width = 5,
p(strong("Colored by CellType")),
plotOutput("Trajectory.example.CellType")),
column(width = 5,
p(strong("Colored by Branch/State")),
plotOutput("Trajectory.example.State"))),
br(),
em("Generated by dr. Anita van den Heuvel, Human Genetics LUMC. Last update 24-01-2019"),
br(),
em("For Questions or Comments, please email to A.van_den_heuvel.HG@lumc.nl")),
tabPanel("Expression.in.tSNE",
sliderInput(inputId = "Expression.limits", label = "Expression.limits", value = c(1.5), min =0.5, max = 50,step = 0.5),
h3("Expression in CellRanger-based tSNE"),
plotOutput("Expression.in.tSNE")),
tabPanel("Heatmap.Primary",
plotOutput("Heatmap.Primary",height = "800px")),
tabPanel("Trajectory",
h3("Pseudotime trajectory plots"),
fluidRow(
column(width = 5,
p(strong("Colored for Cells Included")),
plotOutput("Trajectory.check")),
column(width = 5,
p(strong("Colored by CellType.or.State")),
plotOutput("Trajectory.CellType.or.State")))),
tabPanel("Expr.In.Trajectory",
h3("Dot.color is expression level"),
plotOutput("Trajectory.expression")),
tabPanel("Gene.Expression.Table", dataTableOutput("fData.FSHD.vs.LateMyos"))),
width = 9)))
#################### Define UI #######################
################################### END ######################################
################################### START ######################################
#################### Define Server #######################
server <- function(input,output) {
########
output$Trajectory.example.CellType <- renderPlot({
plot_cell_trajectory(Primary.data, color_by = "CellType", show_branch_points = FALSE, cell_size = 3)})
output$Trajectory.example.State <- renderPlot({
plot_cell_trajectory(Primary.data, color_by = "State",show_branch_points = FALSE, cell_size = 3)})
########
output$Expression.in.tSNE <- renderPlot({
data.for.expression.in.tSNE.plotting <- gbm.Pr
analysis_results <- analysis_results.original.Pr
tsne_proj <- analysis_results$tsne
geneset.full <- data.frame( "x" = strsplit(input$Genes, split = "\n"))
geneset <- subset(geneset.full, geneset.full[,1] %in% fData(data.for.expression.in.tSNE.plotting)[,paste(input$GeneName.or.GeneID)])
visualize_gene_markers(data.for.expression.in.tSNE.plotting,geneset[,1],tsne_proj[c("TSNE.1","TSNE.2")],limits=c(0,input$Expression.limits), marker_size = 2, title = input$Sample.to.be.Analyzed)
})
########
output$Heatmap.Primary <- renderPlot({
data.for.my.own.heatmap <- Primary.data
geneset.exclDUX4 <- data.frame( "x" = strsplit(input$Genes, split = "\n"))
geneset.dux4 <- data.frame( "V1" = c("DUX4", "ENSG00000260596")) # "ENSG00000260596"
colnames(geneset.dux4) <- colnames(geneset.exclDUX4)
geneset <- rbind(geneset.exclDUX4, geneset.dux4)
gene.to.sort.on <- geneset.exclDUX4[1,1]
geneID.to.sort.on <- rownames(subset(fData(data.for.my.own.heatmap), fData(data.for.my.own.heatmap)[,paste(input$GeneName.or.GeneID)] == gene.to.sort.on))
pData(data.for.my.own.heatmap)$GeneExprs.Gene1 <- c(NA)
pData(data.for.my.own.heatmap)$GeneExprs.Gene1 <- as.numeric(t(exprs(data.for.my.own.heatmap[geneID.to.sort.on,])))
prepping.data.for.plotting.Primary(x = data.for.my.own.heatmap,
num.cells.state1 = input$Pr.State1,
num.cells.state2 = input$Pr.State2,
num.cells.state3 = input$Pr.State3,
genes = geneset, geneID.or.gene_short_name = input$GeneName.or.GeneID,
order.on1 = input$Pr.CellOrder[1], order.on2 = input$Pr.CellOrder[2], order.on3 = input$Pr.CellOrder[3])
Rowordering.yes.or.no <- input$Pr.GeneOrder[1]
limits <-input$Pr.limit.settings
if(input$Pr.GeneOrder == "TRUE"){
test <<- heatmap.3(data.to.plot.matrix,
dendrogram = "none",
Rowv = TRUE, Colv = FALSE,
labRow = rownames(data.to.plot.matrix), labCol = "",
trace = "none", density.info = "none",
scale = "row",breaks = seq(-limits,limits,length.out = 76),
col = colorRampPalette(c("blue", "white","red"))(n=75), na.color = "grey",
ColSideColors = clab,
ColSideColorsSize = 4,
cexRow = 1,
main = paste("Primary.data"),
keysize = 1)
legend("topleft", inset = c(0.0,0.8),
legend = CellType.table$CellType,
col = c("blue", "darkred", "green"),
lty = 1,
lwd = 5,
cex= 0.8,
bty = "n",
title = "CellType")
legend("topleft", inset = c(0.01,0.55) ,
legend = State.table$State.nr,
col = c("dodgerblue", "indianred3", "darkgreen"),
lty = 1,
lwd = 5,
cex= 0.8,
bty = "n",
title = "State/Branch")
legend("topleft", inset = c(0.01,0.15) ,
legend = SampleID.color.renaming.table$SampleID,
col = c("red","green","pink","darkgreen","lightblue","blue"),
lty = 1,
lwd = 5,
cex= 0.8,
bty = "n",
title = "SampleID")
}
else
if(input$Pr.GeneOrder == "FALSE"){
test <<- heatmap.3(data.to.plot.matrix,
dendrogram = "none",
Rowv = FALSE, Colv = FALSE,
labRow = rownames(data.to.plot.matrix), labCol = "",
trace = "none", density.info = "none",
scale = "row",breaks = seq(-limits,limits,length.out = 76),
col = colorRampPalette(c("blue", "white","red"))(n=75), na.color = "grey",
ColSideColors = clab,
ColSideColorsSize = 4,
cexRow = 1,
main = paste("Primary.data"),
keysize = 1)
legend("topleft", inset = c(0.0,0.8),
legend = CellType.table$CellType,
col = c("blue", "darkred", "green"),
lty = 1,
lwd = 5,
cex= 0.8,
bty = "n",
title = "CellType")
legend("topleft", inset = c(0.01,0.55) ,
legend = State.table$State.nr,
col = c("dodgerblue", "indianred3", "darkgreen"),
lty = 1,
lwd = 5,
cex= 0.8,
bty = "n",
title = "State/Branch")
legend("topleft", inset = c(0.01,0.15) ,
legend = SampleID.color.renaming.table$SampleID,
col = c("red","green","pink","darkgreen","lightblue","blue"),
lty = 1,
lwd = 5,
cex= 0.8,
bty = "n",
title = "SampleID")
}
})
########
output$fData.FSHD.vs.LateMyos <- renderDataTable({
data.for.my.own.table <- Table.data.Pr
geneset.exclDUX4 <- data.frame( "x" = strsplit(input$Genes, split = "\n"))
geneset.dux4 <- data.frame( "V1" = c("DUX4", "ENSG00000260596")) # "ENSG00000260596"
colnames(geneset.dux4) <- colnames(geneset.exclDUX4)
geneset.full <- rbind(geneset.exclDUX4, geneset.dux4)
geneset <- subset(geneset.full, geneset.full[,1] %in% data.for.my.own.table[,paste(input$GeneName.or.GeneID)])
if(input$Genes == ""){
data.for.table <- data.for.my.own.table}
else
if(input$GeneName.or.GeneID == "gene_short_name"){
data.for.table <- subset(data.for.my.own.table, data.for.my.own.table$gene_short_name %in% geneset[,1])}
else
if(input$GeneName.or.GeneID == "id"){
data.for.table <- subset(data.for.my.own.table, data.for.my.own.table$id %in% geneset[,1])}
as.matrix(data.for.table)
},options = list(scrollX = TRUE))
########
output$Trajectory.check <- renderPlot({
data.for.my.own.heatmap <- Primary.data
geneset.exclDUX4 <- data.frame( "x" = strsplit(input$Genes, split = "\n"))
geneset.dux4 <- data.frame( "V1" = c("DUX4", "ENSG00000260596")) # "ENSG00000260596"
colnames(geneset.dux4) <- colnames(geneset.exclDUX4)
geneset <- rbind(geneset.exclDUX4, geneset.dux4)
prepping.data.for.plotting.Primary(x = data.for.my.own.heatmap,
num.cells.state1 = input$Pr.State1,
num.cells.state2 = input$Pr.State2,
num.cells.state3 = input$Pr.State3,
genes = geneset, geneID.or.gene_short_name = input$GeneName.or.GeneID,
order.on1 = input$Pr.CellOrder[1], order.on2 = input$Pr.CellOrder[2], order.on3 = input$Pr.CellOrder[3])
## plot cells in trajectory -- color_by In.data.or.not
pData(data.for.my.own.heatmap)$InHeatmap <- c(NA)
for(i in 1:length(pData(data.for.my.own.heatmap)$InHeatmap)){
if(pData(data.for.my.own.heatmap)$barcode[i] %in% pData(data.to.plot.genefiltered.reordered)$barcode){
pData(data.for.my.own.heatmap)$InHeatmap[i] <- "1.Included"} else pData(data.for.my.own.heatmap)$InHeatmap[i] <- "2.Excluded"}
plot_cell_trajectory(data.for.my.own.heatmap, color_by = "InHeatmap",show_branch_points = TRUE, cell_size = 3)+
scale_color_manual(values = c("lightblue", "darkred"))+
ggtitle("Primary.data")
})
########
output$Trajectory.CellType.or.State <- renderPlot({
data.for.my.own.heatmap <- Primary.data
plot_cell_trajectory(data.for.my.own.heatmap, color_by = "CellType",show_branch_points = TRUE, cell_size = 3)+
ggtitle("Primary.data")
})
########
output$Trajectory.expression <- renderPlot({
data.for.my.own.heatmap <- Primary.data
geneset.full <- data.frame( "x" = strsplit(input$Genes, split = "\n"))
geneset <- subset(geneset.full, geneset.full[,1] %in% fData(data.for.my.own.heatmap)[,paste(input$GeneName.or.GeneID)])
plot_cell_trajectory(data.for.my.own.heatmap, markers = geneset[,1], use_color_gradient = TRUE, show_branch_points = FALSE, show_backbone = FALSE, show_tree = FALSE, cell_size = 2, markers_linear = FALSE)+
ggtitle("Primary.data")
})
}
#################### Define Server #######################
################################### END ######################################
################################### START ######################################
#################### Launch tool #######################
#### Launch tool
shinyApp(ui = ui,server = server)
#################### Launch tool #######################
################################### END ######################################
################################### START ######################################
####################### Get Packages ##########################
lib.location <- "Packages"
# get packages
library("ggplot2", lib.loc = lib.location)
library("reshape2", lib.loc = lib.location)
library("irlba", lib.loc = lib.location)
library("monocle", lib.loc = lib.location) # http://bioconductor.org/packages/release/bioc/html/monocle.html
library("cellrangerRkit", lib.loc = lib.location)
# source("https://s3-us-west-2.amazonaws.com/10x.files/supp/cell-exp/rkit-install-1.1.0.R") # https://github.com/hb-gitified/cellrangerRkit.git
library("shiny", lib.loc = lib.location)
source("https://raw.githubusercontent.com/obigriffith/biostar-tutorials/master/Heatmaps/heatmap.3.R")
#
# library("Rmisc")
# https://github.com/hb-gitified/cellrangerRkit/blob/master/scripts/rkit-install-1.1.0.R
# library("cellrangerRkit", lib.loc = "Package/")
# install.packages("Package/cellrangerRkit-1.1.0.tar.gz", type="source", repos=NULL)
# install_github('hb-gitified/cellrangerRkit')
# library("stringr")
# library("colorspace")
# library("RBGL")
# library("dplyr")
# library("rhdf5")
# library("stringi")
# library("VGAM")
# library("irlba")
# library("DDRTree")
# library("slam")
# library("combinat")
# library("limma")
# library("jsonlite")
# library("gridExtra")
# library("grid")
# library("RGraphics")
# library("digest")
# library("bit64")
# library("shiny")
# library("BiocGenerics")
# library("Biobase")
# library("FNN")
####################### Get Packages ##########################
################################### END ######################################
################################### START ######################################
####################### Get Data and Settings ##########################
################### get data ready
### Get data PrimaryData
Primary.dataname <- load("PrimaryData/TrajectoryPlot.Rdata")
Primary.data <- get(Primary.dataname)
pData(Primary.data)$SampleID <- colsplit(pData(Primary.data)$barcode, pattern = "-", names = c("a", "b"))$b
gbm.Pr <- load_cellranger_matrix("PrimaryData")
fData(gbm.Pr)$gene_short_name <- fData(gbm.Pr)$symbol
analysis_results.original.Pr <- load_cellranger_analysis_results("PrimaryData")
Table.data.Pr <- fData(Primary.data)[,c("id", "gene_short_name", "num_cells_expressed", "Pooled_mRNAs", "myotube.expression", "FSHDpop.expression", "Fold.change")]
################### get data ready
####################### Get Data and Settings ##########################
################################### END ######################################
################################### START ######################################
#################### Defining new Functions #######################
##### function --- PrimaryData analysis
prepping.data.for.plotting.Primary <- function(x,
num.cells.state1 = c(1,25), num.cells.state2 = c(1,25), num.cells.state3 = c(1,25),
geneID.or.gene_short_name = "gene_short_name", genes = geneset,
order.on1 = "Pseudotime", order.on2 = NA, order.on3 = NA) {
##### data.processing ##########################################################
State1branch <- x[,pData(x)$State == 1]
State1branch <- State1branch[,order(pData(State1branch)$Pseudotime)]
State1branch.cells.to.inlcude <- State1branch[, rownames(pData(State1branch)[c(num.cells.state1[1]:num.cells.state1[2]),])]
State2branch <- x[,pData(x)$State == 2]
State2branch <- State2branch[,order(pData(State2branch)$Pseudotime)]
State2branch.cells.to.inlcude <- State2branch[, rownames(pData(State2branch)[c(num.cells.state2[1]:num.cells.state2[2]),])]
State3branch <- x[,pData(x)$State == 3]
State3branch <- State3branch[,order(pData(State3branch)$Pseudotime)]
State3branch.cells.to.inlcude <- State3branch[, rownames(pData(State3branch)[c(num.cells.state3[1]:num.cells.state3[2]),])]
data.to.plot <<- x[,pData(x)$barcode %in% pData(State1branch.cells.to.inlcude)$barcode | pData(x)$barcode %in% pData(State2branch.cells.to.inlcude)$barcode | pData(x)$barcode %in% pData(State3branch.cells.to.inlcude)$barcode]
# select only wanted genes
my_genes <- row.names(subset(fData(data.to.plot),
fData(data.to.plot)[,paste(geneID.or.gene_short_name)] %in% genes[,1]))
data.to.plot.genefiltered <- data.to.plot[fData(data.to.plot)$id %in% my_genes,]
# Reorder this data.subset for Pseudotime order
ordering.list <- c(pData(data.to.plot.genefiltered)[paste(order.on1)])
if(!is.na(order.on3)){
index <- with(exprs(data.to.plot.genefiltered), order(pData(data.to.plot.genefiltered)[paste(order.on1)], pData(data.to.plot.genefiltered)[paste(order.on2)], pData(data.to.plot.genefiltered)[paste(order.on3)]))}
else if(!is.na(order.on2)){
index <- with(exprs(data.to.plot.genefiltered), order(pData(data.to.plot.genefiltered)[paste(order.on1)], pData(data.to.plot.genefiltered)[paste(order.on2)]))}
else
index <- with(exprs(data.to.plot.genefiltered), order(pData(data.to.plot.genefiltered)[paste(order.on1)]))
data.to.plot.genefiltered.reordered <<- data.to.plot.genefiltered[,index]
rownames(data.to.plot.genefiltered.reordered) <<- fData(data.to.plot.genefiltered.reordered)$gene_short_name
#reorder genes
if(geneID.or.gene_short_name == "gene_short_name"){
genes.included <- subset(genes, genes[,1] %in% fData(data.to.plot.genefiltered.reordered)$gene_short_name)
rownames(data.to.plot.genefiltered.reordered) <<- fData(data.to.plot.genefiltered.reordered)$gene_short_name
row.index <- with(exprs(data.to.plot.genefiltered.reordered), genes.included[,1])
data.to.plot.genefiltered.reordered <<- data.to.plot.genefiltered.reordered[row.index,]}
if(geneID.or.gene_short_name == "id"){
genes.included <- subset(genes, genes[,1] %in% fData(data.to.plot.genefiltered.reordered)$id)
rownames(data.to.plot.genefiltered.reordered) <<- fData(data.to.plot.genefiltered.reordered)$id
row.index <- with(exprs(data.to.plot.genefiltered.reordered), genes.included[,1])
data.to.plot.genefiltered.reordered <<- data.to.plot.genefiltered.reordered[row.index,]}
# Make a matrix & Log2.transformed matrix for the final heatmap
data.to.plot.matrix <<- as.matrix(exprs(data.to.plot.genefiltered.reordered))
data.to.plot.matrix.log2 <<- log((data.to.plot.matrix + 0.1), base = 2)
##### data.processing ##########################################################
##### legend.info ##############################################################
# Make legends for the heatmap
# CellType legend
CellType.colors <- c()
CellType.table <- data.frame("CellType" = c("Unknown","FSHDpop","Healthy.compare.cells"), "CellType.colors" = c("blue", "darkred", "green"))
for (i in 1:length(pData(data.to.plot.genefiltered.reordered)$CellType)){
for(j in 1:length(CellType.table$CellType)){
if(pData(data.to.plot.genefiltered.reordered)$CellType[i] == CellType.table$CellType[j]) {CellType.colors[i] <- as.character(CellType.table$CellType.colors[j])}
}
}
CellType.table <<- CellType.table
# State legend
State.colors <- c()
State.table <- data.frame("State.nr" = c(1,2,3), "State.colors" = c("dodgerblue", "indianred3", "darkgreen"))
for (i in 1:length(pData(data.to.plot.genefiltered.reordered)$State)){
for(j in 1:length(State.table$State.nr)){
if(pData(data.to.plot.genefiltered.reordered)$State[i] == State.table$State.nr[j]) {State.colors[i] <- as.character(State.table$State.colors[j])}
}
}
State.table <<- State.table
# SampleID legend
# first add SampleID column to data
SampleID.color <- c()
SampleID.color.renaming.table <- data.frame("SampleID" = c(1,2,3,4,5,6), "Sample.color" = c("red","green","pink","darkgreen","lightblue","blue"))
for (i in 1:length(pData(data.to.plot.genefiltered.reordered)$SampleID)){
for(j in 1:length(SampleID.color.renaming.table$SampleID)){
if(pData(data.to.plot.genefiltered.reordered)$SampleID[i] == SampleID.color.renaming.table$SampleID[j]) {SampleID.color[i] <- as.character(SampleID.color.renaming.table$Sample.color[j])}
}
}
SampleID.color.renaming.table <<- SampleID.color.renaming.table
color.gradient <- function(x, colors=c("darkblue", "white"), colsteps=75) {
return( colorRampPalette(colors) (colsteps) [ findInterval(x, seq(min(x),max(x), length.out=colsteps)) ] )
}
Pseudotimecol <- color.gradient(x = pData(data.to.plot.genefiltered.reordered)$Pseudotime)
# legends table
clab <- cbind(CellType.colors,State.colors, SampleID.color, Pseudotimecol)
colnames(clab)=c("CellType","State", "SampleID", "Pseudotime")
clab <<- clab
##### legend.info ##############################################################
}
#################### Defining new Functions #######################
################################### END ######################################
################################### START ######################################
#################### Define UI #######################
#### Define UI settings
ui <- fluidPage(
tags$head(tags$style(HTML("label {font-size:75%;}"))),
titlePanel("SingleCellData Analysis"),
sidebarLayout(sidebarPanel(selectInput(inputId = "GeneName.or.GeneID", label = "Used 'gene_short_names' (eg DUX4) or 'geneID' (eg 'ENSG00000260596') ", choices = c("gene_short_name", "id"), selected = c("gene_short_name"),multiple = FALSE, selectize = TRUE),
textAreaInput(inputId = "Genes", label = "Genes.to.plot", value = c("MYOG")),
p(span(strong("Select which cells to visualize (per State)"))),
sliderInput(inputId = "Pr.State1", label = "1.Pr.State1", value = c(1,length(subset(pData(Primary.data)$State, pData(Primary.data)$State == 1))), min =0, max = length(subset(pData(Primary.data)$State, pData(Primary.data)$State == 1)),step = 1),
sliderInput(inputId = "Pr.State2", label = "2.Pr.State2", value = c(1,length(subset(pData(Primary.data)$State, pData(Primary.data)$State == 2))), min =0, max = length(subset(pData(Primary.data)$State, pData(Primary.data)$State == 2)),step = 1),
sliderInput(inputId = "Pr.State3", label = "3.Pr.State3", value = c(1,length(subset(pData(Primary.data)$State, pData(Primary.data)$State == 3))), min =0, max = length(subset(pData(Primary.data)$State, pData(Primary.data)$State == 3)),step = 1),
p(span(strong("Select additional settings"))),
sliderInput(inputId = "Pr.limit.settings", label = "Heatmap color limits", value = c(5), min =1, max = 50,step = 1),
selectInput(inputId = "Pr.GeneOrder",width = "100%", label = "Cluster genes or not (no cluster-tree shown)", choices = c("TRUE", "FALSE"), selected = c("TRUE"),multiple = FALSE),
selectInput(inputId = "Pr.CellOrder",width =