.Rhistory

}
else{
pca = plotPCA(normObject, intgroup = c(input$variable))
}
# plot
output$pca = renderPlot({pca})
# pca interactive brush info
output$pca_info = renderPrint({
brushedPoints(pca[["data"]], input$pca_brush)
})
#### VOLCANO PLOT ####
observeEvent(input$volcPlot, {
# in case user specified a different normalization method, display message that results table and volcano plots are still based on DESeq standard:
if(input$normMethod != "Size Factor Division"){
showNotification("NOTE: Volcano plot is based on Size Factor Divison.")
}
fc_bound = input$volcFcThr
# transform results data
res = na.omit(as.data.frame(ddsRes))
# significance
res$Expression = "NS"                                                                      # create new column "significance", initially all genes = NS
res[res$log2FoldChange >= fc_bound & res$padj < input$alpha,]$Expression = "UP"            # UP if padj > alpha & logFC >= fc threshold
res[res$log2FoldChange <= -fc_bound & res$padj < input$alpha,]$Expression = "DOWN"         # DOWN if padj > alpha & logFC <= - fc threshold
# new column with -log10(padj) - interactive shiny CAN'T handle if you change parameters in aes() of ggplot and behaves weird
res$neglog10_p_value = -log10(res$padj)
# plot:
vp = ggplot(res, aes(x = log2FoldChange, y = neglog10_p_value, color = Expression, tooltip = padj))
output$volc = renderPlot({
vp +
# scatter:
geom_point() +
# lines (alpha and logFC):
geom_hline(yintercept = (log10(input$alpha)/log10(10))*(-1), color = "darkgrey") +       # horizontal for alpha
geom_vline(xintercept = c(fc_bound, -fc_bound), color = "darkgrey") +                    # vertical for fc and -fc
# color:
scale_color_manual(values = c("red", "black", "blue")) +
# x-axis ticks:
scale_x_continuous(breaks = c(round(min(res$log2FoldChange))):round(max(res$log2FoldChange))) # integers from rounded minimum to maximum of the log2FC
}) # render volcano plot close
# volcano plot interactive brush info
output$volc_info = renderPrint({
brushedPoints(res, input$volc_brush)
})
}) # volcano plot button close
#### HEATMAP OF EXPERIMENTS ####
log2normCounts = normCounts
# if data was normalized by size factor division => data is not yet log transformed!
if(input$normMethod == "Size Factor Division"){
# log2normCounts = log2(log2normCounts)            # log of size factor counts, this produces -inf counts were counts were filtered (set to 0) by DESeq
# log2normCounts[log2normCounts == -Inf] = NA      # set these values to NA
# log2normCounts = na.omit(log2normCounts)         # omit NA rows
log2normCounts = logTransform(normCounts)
}
# calculate distance
samp_dist = dist(t(log2normCounts))
# render plot
plotExpr = pheatmap(as.matrix(samp_dist))
output$heatExp = renderPlot({plotExpr})
#### HEATMAP OF TOP GENES (BASED ON VARIANCE) ####
observeEvent(input$plotGeneHeat, {
# Selection of genes:
numberOfGenes = input$geneHeatNo                                                   # number of genes the user wants to display
highVarIndex = head(order(rowVars(log2normCounts), decreasing = TRUE), numberOfGenes)  # indexes of the [numberOfGenes] with highest variance
topVarGenes = log2normCounts[highVarIndex, ]                                       # subset log2 transformed dataset accordingly
topVarGenes = topVarGenes - rowMeans(topVarGenes)                                  # mean centering to acquire (log2-)deviation from the mean
# column annotation:
colAnno = data.frame(colData(dds)[, c(input$variable)])                            # annotation is based on the experimental variables the user chose before pressing the analyze button!
row.names(colAnno) = row.names(colData(dds))                                       # if only one variable is selected, R omits the rownames meaning the need to be re-specified!
colnames(colAnno) = input$variable                                                 # format column name
plotGenes = pheatmap(topVarGenes, annotation_col = colAnno)
# plot heatmap:
output$heatGene = renderPlot({plotGenes}, height = input$geneHeatHeight)
}) # gene Heatmap button close
#### BOXPLOTS ####
output$boxplot = renderPlot(boxplot(log2normCounts))
}) # analyze button close
}) # upload button close
}) # server close
#### ======== START APP ======== ####
shinyApp(ui, server)
# load
library(DESeq2)
library(ggplot2)
library(pheatmap)
library(shiny)
library(shinythemes)
library(rtracklayer)
library(preprocessCore)
sortThatData = function(rawCounts, infoData, gffData){
# Purpose of this function is to sort the info data
# and set the column names of the raw data so they
# match with the info data (e.g. remove .bam ending)
vec = numeric(0) #Character vector containing the correct order of column names. Will be used to sort info data
for(i in 1:ncol(rawCounts)){                   # outer loop sets column of raw data...
checkCol = colnames(rawCounts)[i]
for(j in 1:nrow(infoData)){                  # ...which will be compared with nested loop using grepl
if(grepl(infoData[j, 1], checkCol)){       # QBiC Code must be the first column of the info data for this to work
vec[i] = j                               # vec will contain the correct order the info data must be sorted with
colnames(rawCounts)[i] = infoData[j, 1]  # replace column names with their correpsonding QBiC Code
}
}
}
row.names(infoData) = infoData[,1]                  # set row names of info data to QBiC Code so it can be sorted by column names of count data
# row.names(rawCounts) = rawCounts[, 1]               # set row names of raw data to gene ID
rawCounts = rawCounts[,-(2:7)]                      # remove columns 2 to 6
infoData = infoData[colnames(rawCounts)[-1],]       # sort accordingly
# Change row names of raw counts to corrsponding gene name:
names = gffData[gffData$locus_tag %in% rawCounts$Geneid & gffData$gbkey == "Gene",]$Name  # Match locus_tag of gff with Geneid and get gene names
row.names(rawCounts) = make.names(names, unique = TRUE)  # Problem: Same gene name for > 1 locus tag => will be annotated like this: name, name.1, name.2 ...
return(list(rawCounts, infoData))
}
# Method to logtransform AND remove rows containing -Inf values
logTransform = function(dataset){
log2normCounts = log2(dataset)                   # log of size factor counts, this produces -inf counts were counts were filtered (set to 0) by DESeq
log2normCounts[log2normCounts == -Inf] = NA      # set these values to NA
log2normCounts = na.omit(log2normCounts)         # omit NA rows
return(log2normCounts)
}
#### =========== UI ============ ####
ui = fluidPage(
theme = shinytheme("slate"),
navbarPage("DESeq2 Analysis",
### First tab for uploading and displaying count data, design data, specifying DESeq Parameters ###
tabPanel("Data Upload & Analysis Parameters",
## sidebar (upload request, DESeq specifications) ##
sidebarPanel(
## Data Upload ##
h4("Data Upload"),
div(style = "margin-top: +10px"), # reduce/increase space
# file browsers for raw data and info data
fileInput("countFile", "Upload count data (.txt)", accept = ".txt"),
div(style = "margin-top: -20px"),
fileInput("infoFile", "Upload design data (.tsv)", accept = ".tsv"),
div(style = "margin-top: -15px"),
fileInput("gffFile", "Upload General feature format (.gff)", accept = ".gff"),
div(style = "margin-top: -15px"),
# Upload button
actionButton("upload", "Upload!", width = '100%',class = "btn-warning"),
## DESeq Design ##
div(style = "margin-top: +45px"),
h4("Experimental Design"),
# Dropdown menu to specify design variable for DESeq Analysis
selectInput("variable", "Experimental Variable:", "-"),
# Radio buttons to specify normalization method
radioButtons("normMethod",
"Normalization Method:",
c("Size Factor Division", "VST", "Quantile Normalization")),
# Slider to selecet significance level
sliderInput("alpha",
"Significance Level:",
min = 0.01,
max = 1,
step = 0.01,
value = 0.05),
# Button to start analysis
actionButton("analyze", "Analyze!", width = '100%',class = "btn-warning")
), # side bar close
## Main panel displaying data, results, plots ##
mainPanel(tabsetPanel(type = "tabs",
# Raw data and info table
tabPanel("Raw counts", tableOutput("countTable")),
tabPanel("Design", tableOutput("designTable")),
# Normalized data and results
tabPanel("Normalized Counts", tableOutput("normalizedTable")),
tabPanel("Results", textOutput("resText"), tableOutput("resTable"))
)
) # main panel close
), # tabPanel close
tabPanel("Plots",
tabsetPanel(type = "tabs",
# BOXPLOTS
tabPanel("Boxplots", plotOutput("boxplot")),
# PCA
tabPanel("PCA", plotOutput("pca", brush = "pca_brush"), verbatimTextOutput("pca_info")), # interactive PCA plot
# HEATMAPS
tabPanel("Heatmaps",
tabsetPanel(type = "tabs",
tabPanel("Experiments", plotOutput("heatExp")), # distance heatmap of experiment data
tabPanel("Genes",                               # tabPanel for heatmap of high variance genes
sidebarPanel(
h4("Select amount of Genes:"),
sliderInput("geneHeatNo",            # slider to select amount of genes
"Number",
min = 5,
max = 200,
step = 1,
value = 20),
h4("Plot height:"),
sliderInput("geneHeatHeight",
"Pixels",
min = 400,
max = 2000,
step = 10,
value = 500),
actionButton("plotGeneHeat", "Refresh Plot!", width = '100%', class = "btn-warning")
),
mainPanel(plotOutput("heatGene"))
) # tabPanel "Genes" close
) # tabsetPanel close
), # tabPanel "Heatmaps" close
# VOLCANO PLOT
tabPanel("Volcano",
sidebarPanel(
h4("Select LogFC Threshold:"),
sliderInput("volcFcThr",       # slider to select threshold for the logFC of the volcanoplot
"Threshold (absolute)",
min = 0,
max = 3,
step = 0.1,
value = 1),
actionButton("volcPlot", "Refresh Plot!", width = '100%', class = "btn-warning")
),
mainPanel(plotOutput("volc", brush = "volc_brush"), verbatimTextOutput("volc_info")))
)
) # tabPanel close
) # navBarPage close
) # ui close
#### ========= SERVER ========== ####
server = shinyServer(function(input, output, session){
options(shiny.sanitize.errors = TRUE)
observeEvent(input$analyze, {
## Error message if data has not been uploaded before ##
if(is.null(input$countFile) | is.null(input$infoFile) | is.null(input$gffFile)){
showNotification("Please upload count data, sample preparation info and .gff-file", type = "error")
}
})
##############
### UPLOAD ###
##############
## Display Data if user presses 'Upload!' button ##
observeEvent(input$upload,{
## Exception handling and display of error message if data is missing ##
if(is.null(input$countFile) | is.null(input$infoFile) | is.null(input$gffFile)){
showNotification("Please upload count data, sample preparation info and .gff-file", type = "error")
}
req(input$countFile)
req(input$infoFile)
req(input$gffFile)
## INPUT RAW DATA ##
countFile = input$countFile
rawdat = read.table(countFile$datapath, header = TRUE)
## INPUT INFO DATA ##
infoFile = input$infoFile
infodat = read.csv(infoFile$datapath, sep = '\t')
## INPUT GFF FILE ##
gffFile = input$gffFile
gffdat = as.data.frame(readGFF(gffFile$datapath))
## SORT DATA ##
dats = sortThatData(rawdat, infodat, gffdat)
## DISPLAY DATA ##
output$countTable = renderTable(dats[[1]], rownames = TRUE)
output$designTable = renderTable(dats[[2]], rownames = TRUE)
## UPDATE SELECINPUT BASED ON INFO DATA ##
updateSelectInput(session, "variable", choices = factor(colnames(dats[[2]])))
##################
### RUN DESEQ2 ###
##################
observeEvent(input$analyze, {
## DESIGN MATRIX ##
dsm = model.matrix(~dats[[2]][, c(input$variable)]) # selects column from info data and creates design matrix to be used in the following command
## CREATE DESEQ DATASET ##
dds = DESeqDataSetFromMatrix(countData = dats[[1]][-1], colData = dats[[2]], design = dsm)  # ignore GeneName Column when building DESeq Object
dds = DESeq(dds)                                    # vst, quantil, tpm
## DISPLAY NORMALIZED COUNTS ##
# normalization based on selected method (radio button)
if(input$normMethod == "Size Factor Division"){
normCounts = counts(dds, normalized = TRUE)
}
else if(input$normMethod == "VST"){
normObject = varianceStabilizingTransformation(dds)  # S4 object, will e.g. be used in PCA
normCounts = assay(normObject)
}
# else if(input$normMethod == "Quantile Normalization"){
#   # Quantile normalization is not included in DESeq2 => log-transform counts => remove -inf-values => normalize
#   logCount = logTransform(dats[[1]][,-1])
#   normCounts = normalize.quantiles(as.matrix(logCount), copy = TRUE)
# }
output$normalizedTable = renderTable(normCounts, rownames = TRUE)
## DISPLAY RESULTS ##
ddsRes = results(dds, alpha = input$alpha) # significance level is chosen by user via slider
output$resTable = renderTable(as.data.frame(ddsRes), rownames = TRUE)
# in case user specified a different normalization method, display message that results table and volcano plots are still based on DESeq standard:
if(input$normMethod != "Size Factor Division"){
output$resText = renderText("NOTE: Differential Expression Results will always be based on DESeq2's standard normalization method (Size Factor Division).")
}
##################
##### PLOTS ######
##################
#### PCA ####
# plot PCA based on chosen normalization method
if(input$normMethod == "Size Factor Division"){
pca = plotPCA(rlog(dds), intgroup = c(input$variable))
}
else{
pca = plotPCA(normObject, intgroup = c(input$variable))
}
# plot
output$pca = renderPlot({pca})
# pca interactive brush info
output$pca_info = renderPrint({
brushedPoints(pca[["data"]], input$pca_brush)
})
#### VOLCANO PLOT ####
observeEvent(input$volcPlot, {
# in case user specified a different normalization method, display message that results table and volcano plots are still based on DESeq standard:
if(input$normMethod != "Size Factor Division"){
showNotification("NOTE: Volcano plot is based on Size Factor Divison.")
}
fc_bound = input$volcFcThr
# transform results data
res = na.omit(as.data.frame(ddsRes))
# significance
res$Expression = "NS"                                                                      # create new column "significance", initially all genes = NS
res[res$log2FoldChange >= fc_bound & res$padj < input$alpha,]$Expression = "UP"            # UP if padj > alpha & logFC >= fc threshold
res[res$log2FoldChange <= -fc_bound & res$padj < input$alpha,]$Expression = "DOWN"         # DOWN if padj > alpha & logFC <= - fc threshold
# new column with -log10(padj) - interactive shiny CAN'T handle if you change parameters in aes() of ggplot and behaves weird
res$neglog10_p_value = -log10(res$padj)
# plot:
vp = ggplot(res, aes(x = log2FoldChange, y = neglog10_p_value, color = Expression, tooltip = padj))
output$volc = renderPlot({
vp +
# scatter:
geom_point() +
# lines (alpha and logFC):
geom_hline(yintercept = (log10(input$alpha)/log10(10))*(-1), color = "darkgrey") +       # horizontal for alpha
geom_vline(xintercept = c(fc_bound, -fc_bound), color = "darkgrey") +                    # vertical for fc and -fc
# color:
scale_color_manual(values = c("red", "black", "blue")) +
# x-axis ticks:
scale_x_continuous(breaks = c(round(min(res$log2FoldChange))):round(max(res$log2FoldChange))) # integers from rounded minimum to maximum of the log2FC
}) # render volcano plot close
# volcano plot interactive brush info
output$volc_info = renderPrint({
brushedPoints(res, input$volc_brush)
})
}) # volcano plot button close
#### HEATMAP OF EXPERIMENTS ####
log2normCounts = normCounts
# if data was normalized by size factor division => data is not yet log transformed!
if(input$normMethod == "Size Factor Division"){
# log2normCounts = log2(log2normCounts)            # log of size factor counts, this produces -inf counts were counts were filtered (set to 0) by DESeq
# log2normCounts[log2normCounts == -Inf] = NA      # set these values to NA
# log2normCounts = na.omit(log2normCounts)         # omit NA rows
log2normCounts = logTransform(normCounts)
}
# calculate distance
samp_dist = dist(t(log2normCounts))
# render plot
plotExpr = pheatmap(as.matrix(samp_dist))
output$heatExp = renderPlot({plotExpr})
#### HEATMAP OF TOP GENES (BASED ON VARIANCE) ####
observeEvent(input$plotGeneHeat, {
# Selection of genes:
numberOfGenes = input$geneHeatNo                                                   # number of genes the user wants to display
highVarIndex = head(order(rowVars(log2normCounts), decreasing = TRUE), numberOfGenes)  # indexes of the [numberOfGenes] with highest variance
topVarGenes = log2normCounts[highVarIndex, ]                                       # subset log2 transformed dataset accordingly
topVarGenes = topVarGenes - rowMeans(topVarGenes)                                  # mean centering to acquire (log2-)deviation from the mean
# column annotation:
colAnno = data.frame(colData(dds)[, c(input$variable)])                            # annotation is based on the experimental variables the user chose before pressing the analyze button!
row.names(colAnno) = row.names(colData(dds))                                       # if only one variable is selected, R omits the rownames meaning the need to be re-specified!
colnames(colAnno) = input$variable                                                 # format column name
plotGenes = pheatmap(topVarGenes, annotation_col = colAnno)
# plot heatmap:
output$heatGene = renderPlot({plotGenes}, height = input$geneHeatHeight)
}) # gene Heatmap button close
#### BOXPLOTS ####
output$boxplot = renderPlot(boxplot(log2normCounts))
}) # analyze button close
}) # upload button close
}) # server close
#### ======== START APP ======== ####
shinyApp(ui, server)
#### ======== START APP ======== ####
shinyApp(ui, server)
library(DESeq2)
library(pheatmap)
#----read/transform data ----
dat = read.table("./counts.txt", header = TRUE)
row.names(dat) = dat$Geneid
dat = dat[,-(1:7)] #not required
#----Info data
info = read.csv('./QVTLF_sample_preparations.tsv', sep = '\t')
row.names(info) = info[,1] #For boolean masking
#----provided info data is not structured according to colnames order :-(
vec = numeric(0)
#compare count columns to info data
for(i in 1:ncol(dat)){
checkCol = colnames(dat)[i]
for(j in 1:nrow(info)){
if(grepl(info[j, 1], checkCol)){
vec[i] = info[j, 1]
}
}
}
info = info[vec,] #sort accordingly
colnames(dat) = vec
#----Make new info data----
hours = as.factor(c(5, 2, 2, 5, 2, 2, 5, 2, 2, 5, 5, 5))    #just read it from the samples infos...ugly, aber passt schon für's Einlernen? :-)
type = info$Condition..strain_or_phenotype
newInf = data.frame(info$QBiC.Code, hours, type)
#----Start DESeq-Action----
dds = DESeqDataSetFromMatrix(countData = dat, colData = newInf, design = ~type) #this design made the most sense to me. Is this the same as a linear model?
#PCA
plotPCA(rlog(dds), intgroup = c(colnames(colData(dds))[2], "hours"))
# Trying to make my own PCA
pca_data = plotPCA(rlog(dds), intgroup = c("type", "hours"))
View(pca_data)
View(as.data.frame(pca_data))
ggplot(pca_data[["data"]], aes(x = PC1, y = PC2, color = type, shape = hours)) +
geom_point(size = 2) +
theme(legend.title = element_blank()) +
labs(xlab = pca_data[["labels"]][["x"]])
library(ggplot2)
ggplot(pca_data[["data"]], aes(x = PC1, y = PC2, color = type, shape = hours)) +
geom_point(size = 2) +
theme(legend.title = element_blank()) +
labs(xlab = pca_data[["labels"]][["x"]])
pca_data[["labels"]][["x"]]
ggplot(pca_data[["data"]], aes(x = PC1, y = PC2, color = type, shape = hours)) +
geom_point(size = 2) +
theme(legend.title = element_blank()) +
labs(xlab = "ji")
ggplot(pca_data[["data"]], aes(x = PC1, y = PC2, color = type, shape = hours)) +
geom_point(size = 2) +
theme(legend.title = element_blank()) +
labs(xlab = "ji")
?labs
ggplot(pca_data[["data"]], aes(x = PC1, y = PC2, color = type, shape = hours)) +
geom_point(size = 2) +
theme(legend.title = element_blank()) +
labs(x = pca_data[["labels"]][["x"]])
setwd("./Shiny_DESeq2_Tool")
source("./Components/ShinySeq_Packages.R")
source("./Components/ShinySeq_Functions.R")
source("./Components/ShinySeq_UI.R")
source("./Components/ShinySeq_Server.R")
shinyApp(ui, server)
source("./Components/ShinySeq_Packages.R")
source("./Components/ShinySeq_Functions.R")
source("./Components/ShinySeq_UI.R")
source("./Components/ShinySeq_Server.R")
shinyApp(ui, server)
#----Heatmap ----
deseq_norm[deseq_norm == -Inf] = 0
deseq_norm = log2(counts(dds, normalized = TRUE))
dds = DESeq(dds)
deseq_norm = log2(counts(dds, normalized = TRUE))
#----Heatmap ----
deseq_norm[deseq_norm == -Inf] = 0
samp_dist = dist(t(deseq_norm))
pheatmap(as.matrix(samp_dist))
dist(t(deseq_norm))
View(dist(t(deseq_norm)))
# heatmap of genes:
geneIndex = head(order(rowVars(deseq_norm), decreasing = TRUE), 20)
selGenes = deseq_norm[geneIndex,]
selGenes = selGenes - rowMeans(selGenes)
anno = data.frame(colData(dds)[, c("type")])
row.names(anno) = row.names(colData(dds))
pheatmap(selGenes, annotation_col = anno)
pheatmap(as.matrix(samp_dist))
?pheatmap
pheatmap(as.matrix(samp_dist), color = c("green", "red"))
makeColorRampPalette()
colorRampPalette()
colorRampPalette("green", "red")
colorRampPalette("green", "red", length = 100)
colorRampPalette(c("green", "red"), length = 100)
colorRampPalette(c("green", "red"))(length = 100)
colorRampPalette(c("green", "red"))(100)
pheatmap(as.matrix(samp_dist), color = colorRampPalette(c("green", "red"))(100))
# heatmap of genes:
geneIndex = head(order(rowVars(deseq_norm), decreasing = TRUE), 100)
pheatmap(as.matrix(samp_dist), color = colorRampPalette(c("white", "red"))(100))
pheatmap(as.matrix(samp_dist), color = colorRampPalette(c("white", "orange", "red"))(100))
pheatmap(as.matrix(samp_dist), color = colorRampPalette(c("white", "orange", "yellow" , "red"))(100))
pheatmap(as.matrix(samp_dist), color = colorRampPalette(c("white", "yellow" , "red"))(100))
pheatmap(as.matrix(samp_dist), color = colorRampPalette(c("black", "yellow" , "red"))(100))
source("./Components/ShinySeq_Server.R")
shinyApp(ui, server)
source("./Components/ShinySeq_Server.R")
shinyApp(ui, server)
source("./Components/ShinySeq_Server.R")
shinyApp(ui, server)
source("./Components/ShinySeq_Server.R")
shinyApp(ui, server)
source("./Components/ShinySeq_Server.R")
shinyApp(ui, server)
source("./Components/ShinySeq_Server.R")
shinyApp(ui, server)
source("./Components/ShinySeq_Server.R")
shinyApp(ui, server)
source("./Components/ShinySeq_Server.R")
shinyApp(ui, server)
source("./Components/ShinySeq_Server.R")
shinyApp(ui, server)
source("./Components/ShinySeq_Server.R")
shinyApp(ui, server)
source("./Components/ShinySeq_Server.R")
shinyApp(ui, server)
shinyApp(ui, server)
install.packages('rsconnect')
install.packages("rsconnect")
install.packages("Bioconductor")
install.packages("BiocManager")
install.packages("BiocManager")
source("./Components/ShinySeq_Packages.R")
source("./Components/ShinySeq_Functions.R")
source("./Components/ShinySeq_UI.R")
source("./Components/ShinySeq_Server.R")
shinyApp(ui, server)
source("./Components/ShinySeq_Packages.R")