inst/sake/server.r
In naikai/sake: Single-cell RNA-Seq Analysis and Klustering Evaluation
library(shiny)
library(shinyBS)
library(ggplot2)
library(data.table)
library(magrittr)
library(corrplot)
library(biomaRt)
library(BiocParallel)
library(AnnotationHub)
library(RColorBrewer)
library(matrixStats)
library(tools)
library(snowfall)
library(NMF)
library(DT)
library(gage)
library(gageData)
library(pathview)
shinyServer(function(input, output, session) {
fileinfo <- reactive({
if (input$selectfile == 'upload'){
inFile <- input$file1
req(inFile)
validate(
need(!is.null(inFile), "Please upload a gene count data set")
)
filename <- basename(inFile$name)
filepath <- inFile$datapath
}else if (input$selectfile == "preload"){
req(input$inputdata)
validate(
need(input$inputdata!="", "Please select a gene count data set")
)
filename <- basename(input$inputdata)
filepath <- file.path(raw_fd, input$inputdata)
}else if (input$selectfile == "saved"){
inFile <- input$rda
req(inFile)
validate(
need(!is.null(inFile), "Please upload a .rda file")
)
filename <- basename(inFile$name)
filepath <- inFile$datapath
}else{
return()
}
fileinfo <- list()
fileinfo[['name']] <- filename
fileinfo[['path']] <- filepath
return(fileinfo)
})
### Output Prefix ###
File <- reactive({
if (input$selectfile == 'upload'){
File <- input$file1$name
}else if (input$selectfile == "preload"){
File <- input$inputdata
}else if (input$selectfile == "saved"){
File <- input$rda$name
}else{
return()
}
return(File)
})
File_ext <- reactive({
file_ext(File())
})
file_prefix <- reactive({
if(input$selectfile == "saved"){
gsub(paste0(".", File_ext()), "", File())
}else{
paste( gsub(paste0(".", File_ext()), "", File()),
paste0("k", input$num_cluster),
paste0(input$top.num,input$algorithm),
paste0("nrun", input$nrun),
paste0("predict_", input$predict),
sep=".")
}
})
rda <- reactive({
req(fileinfo())
req(input$selectfile == "saved")
filepath <- fileinfo()[['path']]
validate(
need(File_ext() == "rda" || File_ext() == "RData",
message = paste("file extension", File_ext(), "unknown, please try again"))
)
rda <- local({load(filepath); environment()})
# start checking for all the required variables
validate(
need(!is.null(rda$rawdata), "Variable 'rawdata' does not exist in your .rda file, please check it and upload again") %then%
need(nrow(rda$rawdata) > 0, "Variable 'rawdata' in your .rda file is empty, please check it and upload again") %then%
need(!is.null(rda$nmfres), "Variable 'nmfres' does not exist in your .rda file, please check it and upload again")
# need(!is.null(rda$tsne_2d), "Variable 'tsne_2d' does not exist in your .rda file, please check it and upload again") %then%
# need(!is.null(rda$tsne_3d), "Variable 'tsne_3d' does not exist in your .rda file, please check it and upload again") #%then%
# need(!is.null(rda$dds), "Variable 'dds' does not exist in your .rda file, please check it and upload again")
)
return(rda)
})
observeEvent(rda(), {
tsne_2d$data <- rda()$tsne_2d
tsne_3d$data <- rda()$tsne_3d
})
rawdata <- reactive({
req(fileinfo())
filepath <- fileinfo()[['path']]
n <- 2
withProgress(message = 'Loading data', value = 0, {
incProgress(1/n, detail = "Usually takes ~10 seconds")
if(input$selectfile == "saved"){
rawdata <- rda()$rawdata
}else{
rawdata <- myfread.table(filepath, check.platform=T, sep=input$sep, detect.file.ext=FALSE)
}
})
colnames(rawdata) <- gsub("\\.", "-", colnames(rawdata))
return (rawdata)
})
#' Raw data
output$rawdatatbl = DT::renderDataTable({
rawdata <- rawdata()
validate(
need(ncol(rawdata)>1,
"Please modify the parameters on the left and try again!") %then%
need(!any(duplicated(t(rawdata))),
"There are columns with the exact same value in your data, Please check and try again!")
)
withProgress(message = 'Displying datatable', value = NULL, {
DT::datatable(rawdata[1:20, 1:min(500, ncol(rawdata))], rownames= TRUE,
options = list(scroll = TRUE,
scrollX = TRUE,
scrollY = TRUE,
pageLength = 8
)
) %>% formatRound(1:ncol(rawdata), 3)
})
}, server=TRUE)
#' Info box
output$nameBox <- renderInfoBox({
infoBox(
"FILENAME", fileinfo()[['name']], icon = icon("list"),
color = "green"
)
})
output$sampleBox <- renderInfoBox({
rawdata <- rawdata()
infoBox(
"SAMPLE", ncol(rawdata), icon = icon("list"),
color = "blue"
)
})
output$featureBox <- renderInfoBox({
rawdata <- rawdata()
infoBox(
"FEATURE", nrow(rawdata), icon = icon("list"),
color = "yellow"
)
})
# transform_data <- reactiveValues(data=NULL, rawdata=NULL)
transform_data <- reactiveValues(data=data.frame(a=NULL, b=NULL),
rawdata=data.frame(a=NULL, b=NULL))
run_selsamp <- reactiveValues(go=FALSE)
run_trans <- function(rawdata){
transdata <- rawdata
n <- 2
withProgress(message = 'Transforming data', value = 0, {
incProgress(1/n, detail = "Takes around 5~10 seconds")
if(input$normdata=="takeRPM"){
transdata <- rpm(transdata, colSums(transdata))
}else if(input$normdata=="takesizeNorm"){
sf <- norm_factors(transdata)
transdata <- t(t(transdata)/sf)
}else if(input$normdata=="takeuq"){
transdata <- uq(transdata)
}
transdata <- rmv_constant_0(transdata, pct=0.95, minimum=0)
if(input$transformdata=="takevst"){
transdata <- vst(round(transdata))
}else if(input$transformdata=="takelog"){
validate(
need(!any(transdata<0), "Please make sure all the numbers in the data are positive values")
)
transdata <- log2(transdata + 1)
}
})
return(transdata)
}
observeEvent(rawdata(), {
transdata <- run_trans(rawdata())
transform_data$data <- transdata
transform_data$rawdata <- rawdata()
})
observeEvent(input$run_transf, {
transdata <- run_trans(rawdata())
transform_data$data <- transdata
transform_data$rawdata <- rawdata()
})
output$dl_transformdata <- downloadHandler(
filename <- function() {
paste(file_prefix(), "filtered.txt", sep=".")
},
content = function(file) {
data <- transform_data$data
colnames(data) <- gsub("\\.", "-", colnames(data))
write.table(data, file = file, quote=F, row.names=T, sep="\t")
}
)
output$readDistrib <- renderPlotly({
validate(
need(dim(transform_data$data)[1] > 1, "Did you upload your file or select preloaded data yet?")
)
total <- colSums(transform_data$data)
genecov <- colSums(transform_data$rawdata > 0)
run_selsamp$go <- TRUE
textlab <- paste(names(total),
paste("Total:", signif(total/1000000,3), "Million"),
paste("Gene Cov:", genecov),
sep="<br>")
plot_ly(x=genecov, y=total, text=textlab, source="readdist", #type="scattergl",
mode="markers", hoverinfo="text",
themes="Catherine", marker = list(size=10, opacity=0.65)) %>%
layout(yaxis = list(title = "Total transcripts counts", zeroline=TRUE),
xaxis = list(title = "Gene Coverage", zeroline=TRUE),
dragmode = "select",
showlegend=FALSE)
})
sel_samp <- reactive({
# Get subset based on selection
event.data <- event_data("plotly_selected", source="readdist")
# If NULL dont do anything
if(is.null(event.data) == T) return(NULL)
if(length(event.data$curveNumber) == 0) return(NULL)
# Get selected samples
sel_idx <- subset(event.data, curveNumber == 0)$pointNumber + 1
non_sel_idx <- subset(event.data, curveNumber == 0)$pointNumber + 1
if(max(sel_idx) > ncol(transform_data$data)) return(NULL)
res <- data.frame("Total_Reads"=colSums(transform_data$data)[sel_idx],
"Gene_Coverage"=colSums(transform_data$rawdata>0)[sel_idx])
return(res)
})
output$selsamp_tb <- DT::renderDataTable({
DT::datatable(sel_samp(),
rownames= TRUE,
options = list(scrollX = TRUE,
scrollY = TRUE,
dom = 'tipr',
pageLength = 5
)
)
}, server=TRUE)
output$selsamp_txt <- renderPrint({
cat("Selected samples\n\n")
idx <- input$selsamp_tb_rows_selected
if(length(idx) > 0){
cat(paste(rownames(sel_samp())[idx], sep = ', '))
}else{
cat(NULL)
}
})
### Filter selected samples from user
observeEvent(input$run_filtsamp, {
idx <- input$selsamp_tb_rows_selected
filtsamp <- rownames(sel_samp())[idx]
if(!is.null(filtsamp)){
idx <- match(filtsamp, colnames(transform_data$data))
if(length(idx) > 0){
temp <- transform_data$data[, -idx]
transform_data$data <- temp
temp <- transform_data$rawdata[, -idx]
transform_data$rawdata <- temp
}else{
warning("Can not find selected sample names in transform_data")
}
}
})
merged <- reactive({
rawdata <- transform_data$data
genelist <- NULL
if(!is.null(input$selected_samples) && sum(input$selected_samples %in% colnames(rawdata))==length(input$selected_samples)){
rawdata <- rawdata[, input$selected_samples]
}
if(input$list_type=='Top Ranks'){
select.genes <- select_top_n(apply(rawdata,1,input$math), n=input$top.num, bottom=ifelse(input$orders=="bottom", T, F))
genelist <- names(select.genes)
}else if(input$list_type=='Upload gene list'){
validate(
need(!is.null(input$genefile), "Please upload a file with list of genes, with column name 'Gene' ")
)
genefile <- read.table(input$genefile$datapath, header=T, sep="\t")
ext <- file_ext(input$genefile[1])
n <- 2
withProgress(message = 'Loading gene list', value = 0, {
incProgress(1/n, detail = "Usually takes ~10 seconds")
if (ext=="csv"){
genefile <- read.table(input$genefile$datapath, header=T, sep =',', stringsAsFactors = FALSE)
}else if (ext=="out" || ext=="tsv" || ext=="txt"){
genefile <- read.table(input$genefile$datapath, header=T, sep ='\t', stringsAsFactors = FALSE)
}else{
print("File format doesn't support, please try again")
return(NULL)
}
})
validate(
need(sum(colnames(genefile)=='Gene'), "Please make sure genefile contains column name 'Gene'")
)
validate(
need(length(intersect(genefile$Gene, rownames(rawdata)))>=2, "We did not find enough genes that overlap with the expression data (<2). Please check again and provide a new list! ")
)
genelist <- genefile$Gene
}else if(input$list_type == 'Whole transcriptome'){
genelist <- rownames(rawdata)
}
idx <- unique(match(genelist, rownames(rawdata)))
merged <- rawdata[idx[!is.na(idx)], ]
return(merged)
})
output$filterdatatbl = DT::renderDataTable({
filter_data <- merged()
DT::datatable(filter_data, rownames= TRUE,
selection = 'single',
options = list(scrollX = TRUE,
scroller = TRUE,
pageLength = 5
)
) %>% formatRound(1:ncol(filter_data), 2)
}, server=TRUE)
callModule(feature, "sample", reactive({ merged() }))
callModule(network, "gene", reactive({ merged() }))
#' Scatter plot
observeEvent(transform_data$data, {
transform_data <- transform_data$data
updateSelectizeInput(session, 'cor_sample1',
server = TRUE,
choices = sort(as.character(colnames(transform_data))),
selected = sort(as.character(colnames(transform_data)))[1]
)
updateSelectizeInput(session, 'cor_sample2',
server = TRUE,
choices = sort(as.character(colnames(transform_data))),
selected = sort(as.character(colnames(transform_data)))[2]
)
})
scatterData <- reactive({
transform_data <- transform_data$data
validate(
need(input$cor_sample1 != "", "Please select one sample for X axis"),
need(input$cor_sample2 != "", "Please select one sample for Y axis")
)
x <- transform_data[, input$cor_sample1]
y <- transform_data[, input$cor_sample2]
if(input$scatter_log){
if(input$transformdata == "takelog"){
message("Already did log transformation")
}else{
x <- log2(x+1)
y <- log2(y+1)
}
}
res <- data.frame(x=x, y=y)
return(res)
})
output$scatterPlot <- renderPlotly({
transform_data <- transform_data$data
x <- scatterData()$x %>% as.numeric
y <- scatterData()$y %>% as.numeric
plot_range <- c(0, min(x,y):max(x,y), max(x,y) + 0.2)
reg = lm(y ~ x)
modsum = summary(reg)
r <- signif(cor(x,y), 3)
R2 = signif(summary(reg)$r.squared, 3)
textlab <- paste(paste0("x = ", signif(x,3)), paste0("y = ", signif(y,3)), rownames(transform_data), sep="<br>")
n <- 2
withProgress(message = 'Generating Scatter Plot', value = 0, {
# p <- (x=~x, y=~y, text=~textlab, type="scattergl", #source = "scatter",
# mode="markers", hoverinfo="text",
# #themes="Catherine",
# marker = list(size=8, opacity=0.65)) %>%
p <- plot_ly() %>%
add_trace(x=~x, y=~y, text=~textlab, type="scattergl",
mode="markers", hoverinfo="text",
marker = list(size=8, opacity=0.65)) %>%
add_trace(x=~plot_range, y=~plot_range, type="scattergl", mode = "lines", name = "Ref", line = list(width=2)) %>%
layout(title = paste("cor:", r),
xaxis = list(title = input$cor_sample1, zeroline=TRUE),
yaxis = list(title = input$cor_sample2, zeroline=TRUE),
showlegend=FALSE)
})
# if(input$show_r2){
# p <- p %>% layout( annotations = list(x = x, y = y, text=paste0("R2=", R2), showarrow=FALSE) )
# }
if(input$show_fc){
df <- data.frame(x = plot_range[plot_range>=1], y=plot_range[plot_range>=1]-1)
p <- p %>%
add_trace(data=df, x= ~x, y= ~y, type = "scattergl", mode = "lines", name = "2-FC", line = list(width=1.5, color="#fec44f")) %>%
add_trace(data=df, x= ~y, y= ~x, type = "scattergl", mode = "lines", name = "2-FC", line = list(width=1.5, color="#fec44f"))
}
p #%>% toWebGL()
})
output$scatterdatatbl = DT::renderDataTable({
x <- scatterData()$x
y <- scatterData()$y
if(input$transformdata == "takelog"){
logFC = x - y
}else{
esp <- 0.001
logFC = log2((x+esp)/(y+esp))
}
GeneCard <- paste0("<a href='http://www.genecards.org/cgi-bin/carddisp.pl?gene=",
rownames(transform_data$data), "'>", rownames(transform_data$data), "</a>")
scatter_data <- data.frame(Gene=GeneCard,
X = x,
Y = y,
logFC = logFC,
meanFC = mean(x+y)*logFC) %>%
dplyr::arrange(desc(abs(logFC)))
DT::datatable(scatter_data, rownames= FALSE,
options = list(scrollX = TRUE,
pageLength = 8),
escape = FALSE
) %>% formatRound(2:5, 3)
}, server=TRUE)
### Main part of the program ###
nmf_seed <- eventReactive(input$runNMF, {
nmf_seed <- as.numeric(input$nmf_seed)
if(nmf_seed == 0){
set.seed(as.integer(Sys.time()))
sample(1:999999, 1)
}
})
nmf_res <- eventReactive(input$runNMF, {
merged <- merged()
validate(
need(is.numeric(input$num_cluster) & is.whole(input$num_cluster) & input$num_cluster>1, "Please provide a valid positive numeric integer value for 'Num of clusters(K)'") %then%
need(input$num_cluster <= ncol(merged), "Number of clusters(K) must be smaller than number of samples\nPlease Try selecting another K") %then%
need(is.numeric(input$nrun) & is.whole(input$nrun) & input$nrun>1, "Please provide a valid numeric integer value for 'nrun'")
)
if(input$selectfile == "saved"){
nmfres <- rda()$nmfres
}else{
if(ncol(merged) >= 300){
# toggleModal(session, "nmfModal1", toggle = "open")
# Run NMF through YABI
# test_yabi <- try(system("yabish login yabi buffalo! backends"))
# Will remove yabi later, right now just don't let yabi command start and run NMF on the shiny01 server
test_yabi <- try(system("naikai"))
if(test_yabi == 0){
createAlert(session, "YabiAlert", "YabiAlert1", title = "Running NMF on Amazon Cloud", style = "success",
content = paste0("Sample size is ", ncol(merged),
". We are running tiis job on the cloud and will notify you when the results are ready.<br>"),
append = FALSE)
# pop up window asking for email
data_path <- "/mnt/sake-uploads"
# data_path <- "~/Desktop/sake-uploads"
output <- file.path(data_path, paste0(file_prefix(), ".txt"))
write.table(merged, output, sep="\t", quote=F)
command <- paste("su - centos -c '/home/centos/yabish_NMF_email.sh",
"-d", output,
"-t", nrow(merged),
"-k", input$num_cluster,
"-r", input$nrun,
"-m", input$mode,
"-n", "FALSE",
"-a", input$algorithm,
"-s", input$nmf_seed,
"-c", 8,
"-x", "FALSE",
"-f", 0,
"-q", "FALSE'")
# send email to the user and stop sake
# t1 <- isolate(try(system(command, intern = TRUE)))
print("Sending command to yabi")
t1 <- try(system(command, wait = FALSE))
print(paste0("Command sent:", command))
if(t1 == 0){
Sys.sleep(5)
closeAlert(session, "YabiAlert1")
Sys.sleep(5)
print("closing YabiAlert and stopping sake")
# stopApp(11)
stop(paste("sent jobs to yabi and stop sake:", command))
}else{
createAlert(session, "YabiAlert", "YabiAlert2", title = "Error message from running Yabi", style = "danger",
content = t1,
append = FALSE)
Sys.sleep(5)
stop(paste("command creates erros:", command))
}
}else{
createAlert(session, "NMFAlert", "NMFAlert1", title = "WARNING", style = "warning",
content = paste0("Sample size is ", ncol(merged), ". This might take long to run. <br>"),
append = FALSE)
}
}
# Run NMF on local server
ptm <- proc.time()
withProgress(message = 'Running NMF', value = 0, {
incProgress(1/2, detail = "Takes around 30~60 seconds for each K")
nmfres <- myNMF(merged,
mode=input$mode,
cluster=input$num_cluster,
nrun=input$nrun,
ncores = cores(),
algorithm = input$algorithm,
seed = nmf_seed()
)
b <- proc.time() - ptm
print("### Time to run NMF ### ")
print(b)
})
closeAlert(session, "NMFAlert1")
closeAlert(session, "YabiAlert2")
}
return(nmfres)
})
### Download NMF ###
observeEvent(input$runNMF, {
output$dlNMF_UI <- renderUI({
tagList(
column(width=2, br(), downloadButton("dlNMF", "Download NMF result", class = 'dwnld'))
)
})
})
output$dlNMF <- downloadHandler(
filename <- function() {
paste( file_prefix(), "nmf_results.zip", sep=".")
},
content = function(file) {
tmpdir <- tempdir()
current_dir <- getwd()
setwd(tmpdir)
filenames <- sapply(c("nmf_Groups", "nmf_Features", "Original_plus_NMF"),
function(x) paste(file_prefix(), x, "csv", sep="."))
write.csv(nmf_groups(), filenames[1], quote=F, row.names=F)
write.csv(nmf_features(), filenames[2], quote=F, row.names=F)
write.csv(ori_plus_nmfResult()[["Total"]], filenames[3], quote=F)
for(i in 1:(length(ori_plus_nmfResult())-1)){
filenames <- c(filenames, paste(file_prefix(), paste0("onlyNMF", i), "csv", sep="."))
write.csv(ori_plus_nmfResult()[[paste("NMF", i)]], filenames[length(filenames)], quote=F)
}
nmfres <- nmf_res()
rawdata <- rawdata()
nmfres_filename <- paste(file_prefix(), "nmfres", "rda", sep=".")
save(rawdata, nmfres, file=nmfres_filename)
zip(zipfile=file, files=c(filenames, nmfres_filename))
setwd(as.character(current_dir))
},
contentType = "application/zip"
)
### Estimate K ###
runSummary <- reactive({
req(class(nmf_res()) == "NMF.rank")
nmf_res <- nmf_res()
summary <- nmf_summary(nmf_res)
if(input$mode=="estim"){
summary <- t(summary)
}
return(summary)
})
output$estimSummary = DT::renderDataTable({
runSummary <- runSummary()
DT::datatable(runSummary, rownames= TRUE,
extensions = 'Buttons',
options = list(dom = 'Bfrtip',
buttons =
list('copy', 'print', list(
extend = 'collection',
buttons = list(list(extend='csv',
filename = filename),
list(extend='excel',
filename = filename),
list(extend='pdf',
filename= filename)),
text = 'Download'
)),
scrollX = TRUE,
pageLength = 12
)
) %>% formatRound(1:ncol(runSummary), 2)
}, server=FALSE)
output$estimPlot <- renderPlot({
validate(
need(class(nmf_res()) == "NMF.rank", "Seems like you haven't run NMF 'estim run' yet\nOr the loaded .rda file is from 'real run' result")
)
nmf_res <- nmf_res()
if(input$estimtype=="consensus"){
consensusmap(nmf_res)
}else if(input$estimtype=="quality"){
plot(nmf_res)
}
}, height=777)
### NMF plots ###
output$nmfplot <- renderPlot({
validate(
need(class(nmf_res()) == "NMFfitX1", "Seems like you haven't run NMF 'real' run yet\nOr the loaded .rda file is from 'estim run' result")
)
nmf_res <- nmf_res()
nmf_plot(nmf_res, type=input$plottype, silorder=input$nmfplot_silhouette,
subsetRow = ifelse(input$select_feature_num>0, as.numeric(input$select_feature_num), TRUE ))
},height = 777)
output$silhouetteplot <- renderPlot({
validate(
need(class(nmf_res()) == "NMFfitX1", "Seems like you haven't run NMF 'real' run yet")
)
nmf_res <- nmf_res()
nmf_silhouette_plot(nmf_res, type=input$plottype)
},height = 777)
output$dl_nmf_estimplot <- downloadHandler(
filename <- function() {
paste(file_prefix(), input$mode, "nmf.pdf", sep=".")
},
content = function(file) {
nmf_res <- nmf_res()
pdf(file, width=input$estim_pdf_w, height=input$estim_pdf_h)
if(input$mode == "estim"){
consensusmap(nmf_res)
print(plot(nmf_res))
}
dev.off()
}
)
output$dl_nmf_realplot <- downloadHandler(
filename <- function() {
paste(file_prefix(), input$mode, "nmf.pdf", sep=".")
},
content = function(file) {
nmf_res <- nmf_res()
pdf(file, width=input$real_pdf_w, height=input$real_pdf_h)
if(input$mode == "real"){
nmf_plot(nmf_res, type="samples", silorder=input$nmfplot_silhouette)
nmf_plot(nmf_res, type="features", silorder=input$nmfplot_silhouette,
subsetRow = ifelse(input$select_feature_num>0, as.numeric(input$select_feature_num), TRUE ))
nmf_plot(nmf_res, type="consensus")
if(!is.null(tsne_2d$data)){
print(nmfplottsne())
}
}
dev.off()
}
)
nmf_groups <- reactive({
validate(
need(class(nmf_res()) == "NMFfitX1", "Seems like you haven't run NMF 'real' run yet")
)
withProgress(message = 'Extracting Groups info', value = NULL, {
nmf_extract_group(nmf_res(), type=input$predict)
})
})
output$nmfGroups <- DT::renderDataTable({
nmf_groups <- nmf_groups()
filename <- "nmf_groups"
colnames(nmf_groups)[which(colnames(nmf_groups)=="nmf_subtypes")] <- "Group"
# Get subset based on selection from scatter plot
event.data <- event_data("plotly_selected", source = "nmfGroups")
if(is.null(event.data) == TRUE){
sel_idx <- NULL
}else{
sel_idx <- subset(event.data, curveNumber == 0)$pointNumber + 1
}
DT::datatable(nmf_groups, rownames=FALSE, escape=-1,
selection = list(selected = sel_idx),
options = list(dom = 'rtip',
scrollX = TRUE,
pageLength = 8,
order=list(list(1,'asc'))
)
) %>% formatRound(3:ncol(nmf_groups), 3)
}, server=FALSE)
# Add t-SNE plot next to NMF group
observeEvent(input$run_nmftSNE, {
plot_data <- plot_data()[['plot_data']]
nsamples <- ncol(plot_data()[['plot_data']])
if(input$nmftsne_perplexity*3 > (nsamples-1)) {
createAlert(session, "perplexityAlert", "perplexityAlert1", title = "WARNING", style = "warning",
content = paste("Perpleixty", input$nmftsne_perplexity,
"is too large compared to num of samples", nsamples),
append = FALSE)
}else{
closeAlert(session, "perplexityAlert1")
withProgress(message = 'Running t-SNE', value = NULL, {
incProgress(1/3, detail = "For t-SNE 2D")
try(
tsne_2d$data <- run_tsne(plot_data, iter=input$tsne_iter, dims=2,
initial_dims = input$tsne_pca_num, theta = input$tsne_theta,
perplexity=input$nmftsne_perplexity, cores=cores())
)
incProgress(2/3, detail = "For t-SNE 3D")
try(
tsne_3d$data <- run_tsne(plot_data, iter=input$tsne_iter, dims=3,
initial_dims = input$tsne_pca_num, theta = input$tsne_theta,
perplexity=input$nmftsne_perplexity, cores=cores())
)
})
}
})
nmfplottsne <- reactive({
plot_data <- plot_data()[['plot_data']]
nsamples <- ncol(plot_data)
validate(
need(input$mode=="real", "This part won't work for 'estim' module\nPlease Try NMF 'Real' run") %then%
need(!is.null(tsne_2d$data), "Please hit 'Run t-SNE' button")
)
color <- NULL
# Will add an option to select by default column names from samples
if(length(nmf_groups()$nmf_subtypes)>0){
idx <- match(colnames(plot_data), nmf_groups()$Sample_ID)
nmf_subtypes <- nmf_groups()$nmf_subtypes[idx]
# color <- create.brewer.color(nmf_subtypes, length(unique(nmf_subtypes)), "naikai")
color <- nmf_subtypes
}else if(length(ColSideColors()[['name']][,1])>0){
color <- ColSideColors()[['name']][,1]
}else{
stop("Error: color scheme undefined")
}
# # default by user selection
point_size <- input$plot_point_size - 5
if(input$nmf_prob){
if(input$predict=="consensus"){
point_size <- point_size + point_size * (1-nmf_groups()$sil_width)
}else if(input$predict == "samples"){
point_size <- point_size + point_size * (1-nmf_groups()$prob)
}
}
nmf_tsne_res <- tsne_2d$data
naikai <- plot_tsne(nmf_tsne_res, color=color, alpha=input$plot_point_alpha,
add.label = input$plot_label, save.plot=F, real.plot=F,
add.legend = input$plot_legend,
point.size=point_size, label.size=input$plot_label_size)
### Check if user select rows (samples)
s = input$nmfGroups_rows_selected
if (length(s)){
if(!input$nmf_prob){
sub_color <- create.brewer.color(nmf_groups()$nmf_subtypes, length(unique(color)), "naikai")
sub_tsne_res <- parse_tsne_res(nmf_tsne_res) %>% as.data.frame %>% dplyr::slice(s)
naikai <- naikai + geom_point(data=sub_tsne_res, aes(x=x, y=y), size=point_size+2, colour=sub_color[s], alpha=input$plot_point_alpha)
}
}
p <- ggplotly(naikai, source="nmfGroups") %>% layout(dragmode = "select")
return(p)
})
output$nmftsneplot <- renderPlotly({
withProgress(message = 'Genrating t-SNE plot', value = NULL, {
nmfplottsne()
})
})
output$nmfgrp_var <- renderPlotly({
nmf_groups <- nmf_groups()
plot_data <- plot_data()[['plot_data']]
idx <- match(colnames(plot_data), nmf_groups$Sample_ID)
colnames(plot_data) <- paste0("NMF", nmf_groups$nmf_subtypes[idx])
dd <- plot_data %>% t %>% reshape2::melt(.) %>% group_by(Var2, Var1) %>% summarise(var=var(value))
num_clus <- dd$Var1 %>% levels %>% length
mycolor <- create.brewer.color(1:num_clus, num=num_clus, name="naikai")
dd$Var1 <- factor(dd$Var1, levels = paste0("NMF", 1:num_clus))
a <- ggplot(dd, aes(var, fill=Var1, colour=Var1)) +
geom_density(alpha=input$cl_alpha) + scale_x_log10() +
scale_colour_manual(values = mycolor) +
scale_fill_manual(values = mycolor) +
labs(x="Log (variance)", y="Density", colour="", fill="") +
theme_bw()
p <- plotly_build(a)
# remove duplicated name
p$data <- lapply(p$data, FUN = function(x){
x$name <- sub("\\((\\S+\\d),.*", "\\1", x$name, perl=TRUE)
x$line$width <- 1.2
return(x)
})
p$layout$legend$font$size <- as.numeric(input$cl_lgsize)
p$layout$xaxis$tickfont$size <- input$cl_tickfont
p$layout$xaxis$titlefont$size <- input$cl_axisfont
p$layout$xaxis$titlefont$family <- "Arial"
p$layout$yaxis$tickfont$size <- input$cl_tickfont
p$layout$yaxis$titlefont$size <- input$cl_axisfont
p$layout$yaxis$titlefont$family <- "Arial"
p
})
output$nmfgrp_expgene <- renderPlotly({
rawdata <- rawdata()
nmf_groups <- nmf_groups()
idx <- match(colnames(rawdata), nmf_groups$Sample_ID)
rawdata <- rawdata[, !is.na(idx)]
colnames(rawdata) <- paste0("NMF", nmf_groups$nmf_subtypes[idx[!is.na(idx)]])
dd <- data.frame(Group=colnames(rawdata), Genes=colSums(rawdata>0))
num_clus <- dd$Group %>% levels %>% length
mycolor <- create.brewer.color(1:num_clus, num=num_clus, name="naikai")
a <- ggplot(dd, aes(Genes, fill=Group, colour=Group)) +
geom_density(alpha=input$cl_alpha) +
scale_fill_manual(values = mycolor, guide="none") +
scale_color_manual(values = mycolor) +
labs(x="No. of expressed genes", y="Density", colour="", fill="") +
theme_bw()
p <- plotly_build(a)
# remove duplicated name
p$data <- lapply(p$data, FUN = function(x){
x$name <- sub("\\((\\S+\\d),.*", "\\1", x$name, perl=TRUE)
x$line$width <- 1.2
return(x)
})
p$layout$legend$font$size <- as.numeric(input$cl_lgsize)
p$layout$xaxis$tickfont$size <- input$cl_tickfont
p$layout$xaxis$titlefont$size <- input$cl_axisfont
p$layout$xaxis$titlefont$family <- "Arial"
p$layout$yaxis$tickfont$size <- input$cl_tickfont
p$layout$yaxis$titlefont$size <- input$cl_axisfont
p$layout$yaxis$titlefont$family <- "Arial"
p
})
output$nmfgrp_coef <- renderPlotly({
nmf_groups <- nmf_groups()
plot_data <- plot_data()[['plot_data']]
idx <- match(colnames(plot_data), nmf_groups$Sample_ID)
nmf_groups <- nmf_groups[idx, ]
colnames(plot_data) <- paste0("NMF", nmf_groups$nmf_subtypes)
pair_cor <- function(x) {
res <- WGCNA::cor(x, nThreads = 4)
return(res[lower.tri(res)])
}
num_clus <- nmf_groups$nmf_subtypes %>% unique %>% length
mycolor <- create.brewer.color(1:num_clus, num=num_clus, name="naikai")
num_poss_cor <- nmf_groups$nmf_subtypes %>% table %>% apply(., 1, function(x) x * (x-1) / 2) %>% sum
group_cor <- data.frame(matrix(NA_real_, nrow=num_poss_cor, ncol=2))
j <- 1
for(i in 1:num_clus){
idx <- nmf_groups$nmf_subtypes == i
res <- pair_cor(plot_data[, idx])
group_cor[j:(j+length(res)-1), ] <- cbind(paste0("NMF", i), res)
j <- j + length(res)
}
colnames(group_cor) <- c("NMF", "Cor")
group_cor$Cor <- as.numeric(group_cor$Cor)
a <- ggplot(data=group_cor, aes(x=NMF, y=Cor, color=NMF)) +
geom_boxplot(aes(color = NMF)) +
stat_summary(fun.y=mean, shape=19, col='black', geom='point') +
theme_bw() +
scale_colour_manual(values = mycolor) +
labs(x="", y="Cor", colour="")
# ggplotly(a)
p <- plotly_build(a)
# remove outliers for plotly boxplot
p$data <- lapply(p$data, FUN = function(x){
if(x$type == "box"){
x$marker = list(opacity = 0)
}
return(x)
})
p$layout$legend$font$size <- as.numeric(input$cl_lgsize)
p$layout$xaxis$tickfont$size <- input$cl_tickfont
p$layout$xaxis$titlefont$size <- input$cl_axisfont
p$layout$xaxis$titlefont$family <- "Arial"
p$layout$yaxis$tickfont$size <- input$cl_tickfont
p$layout$yaxis$titlefont$size <- input$cl_axisfont
p$layout$yaxis$titlefont$family <- "Arial"
p
})
nmf_features <- reactive({
validate(
need(input$mode=="real", "This part won't work for eastimating k\nPlease Try NMF 'Realrun'") %then%
need(class(nmf_res())=="NMFfitX1", "Please rerun NMF with 'real run'")
)
nmf_res <- nmf_res()
# add more options to select features here. then compare with ClaNC
# Add rank by MAD expression
if(input$select_method=="total"){
res <- nmf_extract_feature(nmf_res, method="total") %>% unique
}
else if(input$select_method=="default"){
res <- nmf_extract_feature(nmf_res, method=input$select_method, manual.num = as.numeric(input$select_feature_num)) %>% unique
validate(
need(nrow(res)!=1 & !is.na(res$Gene), "None of the features passed through default NMF filtering criteria, please manually specify number of genes from each group or use 'select features by Rank'") %then%
need(nrow(res)>1, "Not enought features are selected by default! \nMaybe try with more input features or try select features by Rank")
)
}else if(input$select_method=="rank"){
res <- nmf_extract_feature(nmf_res, rawdata = merged(), method=input$select_method, manual.num = as.numeric(input$select_feature_num),
FScutoff=input$select_FScutoff, math=input$select_math) %>% unique
validate(
need(nrow(res)>1, "Not enought features! \nMaybe try lower the FeatureScore cutoff")
)
}
return(res)
})
nmf_features_annot <- reactive({
data <- nmf_features() %>% as.data.table
merged <- NULL
n <- 2
withProgress(message = 'Extracting features info', value = 0, {
incProgress(1/n, detail = "Takes around 5~10 seconds")
ext.func.data <- ext_gene_function(data$Gene, ensembl = ensembl())
setkey(data, "Gene")
setkey(ext.func.data, "hgnc_symbol")
merged <- ext.func.data[data, nomatch=0, allow.cartesian = TRUE]
merged <- as.data.frame(merged) %>%
dplyr::arrange(Group, desc(featureScore))
merged <- merged[, c(4,2,7,6,8)]
merged <- unique(merged)
merged$featureScore <- round(merged$featureScore, 3)
merged$prob <- round(merged$prob, 3)
})
return(merged)
})
output$nmfFeatures <- DT::renderDataTable({
filename <- "nmf_features"
DT::datatable(nmf_features_annot(), rownames=FALSE, escape=-1,
extensions = 'Buttons', selection = 'multiple',
options = list(dom = 'Bfrtip',
buttons =
list('copy', 'print', list(
extend = 'collection',
buttons = list(list(extend='csv',
filename = filename),
list(extend='excel',
filename = filename),
list(extend='pdf',
filename= filename)),
text = 'Download'
)),
pageLength = 50,
autoWidth = TRUE
)
)
}, server = FALSE)
output$nmf_vioplot <- renderPlotly({
req(nmf_features_annot())
validate(
need(!is.null(input$nmfFeatures_rows_selected), "Please select a gene")
)
withProgress(message = 'Generating box/violin plot', value = NULL, {
gene <- nmf_features_annot()[input$nmfFeatures_rows_selected, "GeneCard"] %>%
gsub(".*'>(.*)</a>", "\\1", .)
gene.data <- transform_data$rawdata[gene, ]
gene.data <- gene.data %>%
tibble::rownames_to_column() %>%
dplyr::rename(genename = rowname) %>%
tidyr::gather(Sample, Expr, -genename) %>%
dplyr::mutate(NMF = factor(paste0("NMF", rep(nmf_groups()$nmf_subtypes, each=length(gene)))))
if(input$sel_vioscale == "raw"){
print('Use original scale')
}else if(input$sel_vioscale == "log2"){
gene.data$Expr <- log2(gene.data$Expr+1)
}else if(input$sel_vioscale == "log10"){
gene.data$Expr <- log10(gene.data$Expr+1)
}
# match the coloring from PCA, t-SNE, and heatmap
num_clus <- gene.data$NMF %>% levels %>% length
mycolor <- create.brewer.color(1:num_clus, num=num_clus, name="naikai")
if(input$sel_grp == "Manual"){
col.num <- as.numeric(factor(mycolor, levels = mycolor))
col.lvl <- as.character(mycolor)
col.lvl <- col.lvl[as.numeric(input$sel_grp_order)]
mycolor <- col.lvl[col.num]
}
# add option to select either boxplot or violin plot
if(input$sel_vioplot == "violin"){
a <- ggplot(data=gene.data, aes(x=NMF, y=Expr, color=NMF)) +
geom_violin(aes(color = NMF), scale="width", width=0.6, show.legend = FALSE) +
geom_jitter(aes(color=NMF), alpha=0.6, width=0.1, show.legend = FALSE) +
theme_bw() +
facet_wrap(~ genename, ncol = input$sel_ncol) +
scale_colour_manual(values = mycolor) +
labs(x="", y="", colour="")
}else if(input$sel_vioplot == "box"){
a <- ggplot(data=gene.data, aes(x=NMF, y=Expr, color=NMF)) +
geom_boxplot(aes(color = NMF), show.legend = FALSE) +
geom_jitter(aes(color=NMF), alpha=0.6, width=0.1, show.legend = FALSE) +
theme_bw() +
facet_wrap(~ genename, ncol = input$sel_ncol) +
scale_colour_manual(values = mycolor) +
labs(x="", y="", colour="")
}else{
warning(paste("Unknown plot type:", input$sel_vioplot, "Please check again"))
}
# custom range for Y-axis is only being triggered when ymax > 0
if(!is.na(input$sel_ymax) & input$sel_ymax>0){
validate(
need(input$sel_ymax > input$sel_ymin, "Please make sure 'ymax' is greater than 'ymin'")
)
a <- a + ylim(c(input$sel_ymin, input$sel_ymax))
}
p <- plotly_build(a)
# remove outliers for plotly boxplot
p[[1]]$data <- lapply(p[[1]]$data, FUN = function(x){
if(x$type == "box"){
x$marker = list(outliercolor = 0)
}
return(x)
})
p[[1]]$layout$margin$l <- p[[1]]$layout$margin$l + 5
p[[1]]$layout$legend$font$size <- input$sel_legend_size
p[[1]]$layout$showlegend = input$sel_show_legend
# loop through ncol(Xaxis) and gene (Yaxis)
num_xaxis <- grep("xaxis", names(p[[1]]$layout)) %>% length
num_yaxis <- grep("yaxis", names(p[[1]]$layout)) %>% length
for(i in 1:num_xaxis){
p[[1]]$layout[[paste0("xaxis", ifelse(i>1, i, ""))]]$tickfont$size <- input$sel_xfont_size
}
for(i in 1:num_yaxis){
p[[1]]$layout[[paste0("yaxis", ifelse(i>1, i, ""))]]$tickfont$size <- input$sel_yfont_size
}
for(i in 1:length(gene)){
p[[1]]$layout$annotations[[i]]$font$size <- input$sel_title_size
}
p
})
})
ori_plus_nmfResult <-reactive({
rawdata <- rawdata()
validate(
need(!is.null(nmf_groups()), "Please run NMF and get the result first")
)
res <- list()
idx <- match(colnames(rawdata), nmf_groups()$Sample_ID)
newdata <- rawdata
colnames(newdata) <- paste(paste0("NMF",nmf_groups()$nmf_subtypes[idx]), colnames(newdata), sep="_")
res[["Total"]] <- newdata
# add separate files for each NMF group
for(i in 1:length(unique(nmf_groups()$nmf_subtypes))){
subidx <- nmf_groups()$nmf_subtypes == i
res[[paste("NMF", i)]] <- newdata[, idx[subidx]]
}
return(res)
})
### Statistics ###
output$runStatSummary <- renderTable({
runSummary()
}, sanitize.text.function = function(x) x, options = list(orderClasses = TRUE, lengthMenu = c(10,25,50,100, NROW(merged)), pageLength=25)
)
output$runStatPlot <- renderPlot({
nmf_res <- nmf_res()
plot(nmf_res)
}, height=666)
ensembl <- reactive({
# Need to add a testing here when the BioMart is down
# Or find offline annotation db
# ensembl = useMart(biomart="ENSEMBL_MART_ENSEMBL", host="grch37.ensembl.org", path="/biomart/martservice" ,dataset="hsapiens_gene_ensembl")
ensembl = useMart(biomart="ENSEMBL_MART_ENSEMBL", host="jul2016.archive.ensembl.org", path="/biomart/martservice" ,dataset="hsapiens_gene_ensembl")
})
#' Heatmap setting
geneListInfo <- reactive({
rawdata <- transform_data$data
gene.list <- NULL
gene.groups <- NULL
if(input$heatmapGene=='nmf'){
nmf_features <- nmf_features()
validate(
need(nrow(nmf_features)>0, "NMF not run yet? \nPlease run NMF classification and then try again!")
)
gene.list <- nmf_features$Gene
}
else if(input$heatmapGene=='rank'){
select.genes <- select_top_n(apply(rawdata,1,input$heat.math), n=as.numeric(input$heat.top.num), bottom=F)
gene.list <- names(select.genes)
}
else if(input$heatmapGene=='preload'){
validate(
need(!is.null(input$predefined_list), "Please select at least one predefined gene list")
)
for(i in 1:length(input$predefined_list)){
gene.list <- c(gene.list, pre_genelist[[ input$predefined_list[i] ]]$Gene)
# if file contains 'Group' info in one of the columns, use that for RowColors
if("Group" %in% colnames(pre_genelist[[ input$predefined_list[i] ]])){
gene.groups <- c(gene.groups, pre_genelist[[ input$predefined_list[i] ]]$Group)
}
}
}
else if(input$heatmapGene=='manual'){
validate(
need(!is.null(input$manual_genes), "You didn't select any genes. Please select genes from the data") %then%
need(length(input$manual_genes)>=2, "Please select at least two genes! ")
)
gene.list <- input$manual_genes
}
else if(input$heatmapGene=='upload'){
validate(
need(!is.null(input$heatmapGeneFile), "Please upload a file with list of genes, with column name 'Gene' ")
)
ext <- file_ext(input$heatmapGeneFile)[1]
if (ext=="csv"){
genefile <- read.csv(input$heatmapGeneFile$datapath, header=T, stringsAsFactors = FALSE)
}else if (ext=="out" || ext=="tsv" || ext=="txt"){
genefile <- read.table(input$heatmapGeneFile$datapath, header=T, sep="\t", stringsAsFactors = FALSE)
}else{
print("File format doesn't support, please use '.csv', '.txt', '.out', '.tsv' and try again")
return(NULL)
}
validate(
need(sum(colnames(genefile)=='Gene'), "Please make sure genefile contains column name 'Gene'")
)
gene.list <- genefile$Gene
if("Group" %in% colnames(genefile)){
gene.groups <- genefile$Group
}
}
validate(
need(length(intersect(gene.list, rownames(rawdata)))>=2, "We did not find enough genes that overlap with the expression data (<2). \nPlease check again and provide a new list! \nMaybe this is from species other than Human (default gene list are from Human)")
)
# only select unique gene name on rawdata
idx <- unique(match(gene.list, rownames(rawdata)))
na.genes <- gene.list[is.na(idx)]
if(length(na.genes)>0){
print(c("These genes cannot be found in the data:", paste(na.genes, collapse=", ")))
}
data.rownames <- rownames(rawdata)[idx[!is.na(idx)]]
# match back those gene names to genelist to get idx for gene.groups
ii <- match(data.rownames, gene.list)
if(!is.null(gene.groups)){
gene.groups <- gene.groups[ii]
}
gene.list <- data.rownames
# finish loading or selecting gene list file
res <- list()
res[['list']] <- gene.list
res[['group']] <- gene.groups
return(res)
})
OrdColopts <- reactive({
c("Default", "Hierarchical", "Filename", "GeneExpr")
})
ColClrByopts <- reactive({
c("Filename")
})
observeEvent(rawdata(), {
updateSelectizeInput(session, 'OrdCol',
server = TRUE,
choices = as.character(OrdColopts()),
selected = "Hierarchical"
)
updateSelectizeInput(session, 'ColClrBy',
server = TRUE,
choices = as.character(ColClrByopts()),
selected = "Filename"
)
})
observeEvent(input$runNMF, {
updateSelectizeInput(session, 'OrdCol',
server = TRUE,
choices = as.character(c(OrdColopts(), "NMF Group")),
selected = "NMF Group"
)
updateSelectizeInput(session, 'ColClrBy',
server = TRUE,
choices = as.character(c(ColClrByopts(), "NMF Group")),
selected = "NMF Group"
)
})
# separate out this for PCA and t-SNE plots
plot_data <- reactive({
rawdata <- transform_data$data
genelist <- geneListInfo()[['list']]
idx <- match(genelist, rownames(rawdata))
plot_data <- rawdata[idx,]
res <- list()
res[["plot_data"]] <- plot_data
res[["genelist"]] <- genelist
return(res)
})
heatmap_data <- reactive({
rawdata <- transform_data$data
genelist <- geneListInfo()[['list']]
gene.groups <- geneListInfo()[['group']]
# only select unique gene name on rawdata
idx <- match(genelist, rownames(rawdata))
heatmap_data <- rawdata[idx,]
if(input$OrdCol == 'Filename'){ # Add options to sort by column groups
column.names <- strsplit(colnames(heatmap_data), split="\\_")
validate(
need(input$sortcolumn_num <= length(column.names[[1]]), "Group num is too big, select a smaller number!")
)
idx <- order(sapply(column.names, function(x) x[[input$sortcolumn_num]]))
}else if(input$OrdCol == 'GeneExpr'){ # Add options to sort by expression value of specific gene
validate(
need(input$sortgene_by!="", "Please select a gene!")
)
idx <- order(subset(heatmap_data, rownames(heatmap_data)==input$sortgene_by))
}
else if(input$OrdCol == 'NMF Group') { # Add options to sort by nmf groups
idx <- nmf_groups()$nmf_subtypes %>% order
}else{
idx <- colnames(heatmap_data)
}
heatmap_data <- heatmap_data[, idx]
res <- list()
res[["heatmap_data"]] <- heatmap_data
res[["genelist"]] <- genelist
res[["gene.groups"]] <- gene.groups
return(res)
})
color <- reactive({
color <- switch(input$heatColorScheme,
"RdBu_Brewer" = colorRampPalette(rev(brewer.pal(input$color_num, "RdBu")))(51),
"RdYlBu_Brewer" = colorRampPalette(rev(brewer.pal(input$color_num, "RdYlBu")))(51),
"Red_Green" = colorRampPalette(c("red2", "white", "green"))(input$color_num),
"Red_Blue" = colorRampPalette(c("blue", "white", "red"))(input$color_num),
"Rainbow" = colorRampPalette(c("lightblue", "blue", "white", "gold2", "red2"))(input$color_num),
"Black_White" = colorRampPalette(c("black", "white"))(input$color_num)
)
return(color)
})
ColScheme <- reactive({
ColScheme <- list()
ColScheme <- c(input$ColScheme1, input$ColScheme2, input$ColScheme3, input$ColScheme4, input$ColScheme5, input$ColScheme6)
return(ColScheme)
})
ColSideColors <- reactive({
heatmap_data <- heatmap_data()[['heatmap_data']]
res <- list()
if(input$ColClrBy=='NMF Group'){
req(nmf_groups)
# rematch the sample order to the sorted order based on selection from input$OrdCol
idx <- match(colnames(heatmap_data), nmf_groups()$Sample_ID)
nmf_subtypes <- nmf_groups()$nmf_subtypes[idx]
res[['color']] <- create.brewer.color(nmf_subtypes, length(unique(nmf_subtypes)), "naikai") %>% as.matrix
res[['name']] <- paste0("NMF", nmf_subtypes) %>% as.matrix
}else if(input$ColClrBy == 'Filename'){
res <- name_to_color(colnames(heatmap_data), split_pattern ="\\_",
num_color = 6,
ColScheme = ColScheme()[1:6] )
# num_color = input$ColSideColorsNum,
# ColScheme = ColScheme()[1:input$ColSideColorsNum] )
# add options to select which column to show
if(is.null(input$SelColSideColors)){
res[['color']] <- NULL
res[['name']] <- NULL
}else{
res[['color']] <- res[['color']][, as.numeric(input$SelColSideColors)] %>% as.matrix()
res[['name']] <- res[['name']][, as.numeric(input$SelColSideColors)] %>% as.matrix()
}
}
return(res)
})
RowSideColors <- reactive({
groups <- heatmap_data()[['gene.groups']]
row.color <- NULL
res <- list()
res[['color']] <- NULL
res[['name']] <- NULL
if(!is.null(groups)){
row.color <- rbind(row.color, create.brewer.color(groups, input$RowColScheme.num, input$RowColScheme))
res[["color"]] <- as.matrix(row.color)
res[["name"]] <- as.matrix(groups)
}
return(res)
})
Colv <- reactive({
Colv <- NULL
validate(
need(nrow(heatmap_data()[['heatmap_data']] > 2), "Did not get enough genes that overlapped with your data\nPlease try another gene list or use rank from data set")
)
if(input$OrdCol=="Hierarchical"){
Colv <- heatmap_data()[['heatmap_data']] %>% t %>% dist(method=input$distance) %>% hclust(method=input$linkage) %>% as.dendrogram
ColMean <- colMeans(heatmap_data()[['heatmap_data']], na.rm=T)
Colv <- reorderfun(Colv, ColMean)
}
return(Colv)
})
Rowv <- reactive({
Rowv <- NULL
validate(
need(nrow(heatmap_data()[['heatmap_data']] > 2), "Did not get enough genes that overlapped with your data\nPlease try another gene list or use rank from data set")
)
if(input$OrdRow=="hier"){
Rowv <- heatmap_data()[['heatmap_data']] %>% dist(method=input$distance) %>% hclust(method=input$linkage) %>% as.dendrogram
RowMean <- rowMeans(heatmap_data()[['heatmap_data']], na.rm=T)
Rowv <- reorderfun(Rowv, RowMean)
}
return(Rowv)
})
dendrogram <- reactive({
if(is.null(Rowv()) & is.null(Colv())){
return("none")
}else if(is.null(Rowv()) & !is.null(Colv())){
return("column")
}else if(!is.null(Rowv()) & is.null(Colv())){
return("row")
}else if(!is.null(Rowv()) & !is.null(Colv())){
return("both")
}
})
### Observe and update input selection ###
observe({
transform_data <- transform_data$data
if(!is.null(transform_data)){
if (length(rownames(transform_data)) > 1){
# update the render function for selectize - update manual if user choose to select its own markers
updateSelectizeInput(session, 'manual_genes',
server = TRUE,
choices = as.character(rownames(transform_data))
)
}
}
})
observe({
if (!is.null(geneListInfo()[['list']])){
updateSelectizeInput(session, 'sortgene_by',
server = TRUE,
choices = as.character(geneListInfo()[['list']])
)
}
})
output$d3heatmap <- renderD3heatmap({
heatmap_data <- heatmap_data()[['heatmap_data']]
ColSideColors <- ColSideColors()[["color"]]
if(!is.null(ColSideColors)){
ColSideColors <- t(ColSideColors[, ncol(ColSideColors):1])
}
myd3Heatmap(
as.matrix(heatmap_data),
type="heatmap",
dist.m=input$distance,
hclust.m=input$linkage,
scale="row",
scale.method=input$scale_method,
scale.first=input$scale_first,
RowSideColors = RowSideColors()[["color"]],
ColSideColors = ColSideColors,
color = color(),
Colv = Colv(),
Rowv = Rowv(),
cexCol=input$cexRow,
cexRow=input$cexCol,
na.rm = T,
show_grid = FALSE,
dendrogram = dendrogram(),
dendro.ord = "manual",
anim_duration = 0
)
})
plotInput <- function(){
heatmap_data <- heatmap_data()[['heatmap_data']]
filename <- paste(input$type, input$math, paste0(input$orders, paste0(nrow(merged()), "genes")), sep=".")
myHeatmap.3(as.matrix(heatmap_data),
type="heatmap",
color=color(),
RowSideColors = RowSideColors()[["color"]],
RowSideColors.name=RowSideColors()[["name"]],
ColSideColors = ColSideColors()[["color"]],
ColSideColorsSize=input$ColSideColorsSize,
ColSideColors.name=ColSideColors()[["name"]],
cor.method=input$cor.method,
dist.m=input$distance,
hclust.m=input$linkage,
title=input$title,
Colv = Colv(),
Rowv = Rowv(),
scale="row",
scale.method=input$scale_method,
scale.first=input$scale_first,
dendro.ord="manual",
dendrogram=dendrogram(),
cexCol=input$cexCol,
cexRow=input$cexRow,
col.legend = input$col_legend,
row.legend = input$row_legend,
file.prefix=filename,
save.image=F
)
}
output$dl_heatmap<- downloadHandler(
filename <- function() {
paste(file_prefix(), "heatmap.pdf", sep="_")
},
content = function(file) {
pdf(file, width = input$heat_pdf_wd, height = input$heat_pdf_ht)
plotInput()
dev.off()
}
)
plot_colopts <- reactive({
c("Default", "Filename", "GeneExpr", "ReadDepth", "NumExpresGenes")
})
observeEvent(rawdata(), {
updateSelectizeInput(session, 'pt_col',
server = TRUE,
choices = as.character(plot_colopts()),
selected = "Filename"
)
choices <- rawdata() %>% colnames() %>% strsplit("_") %>% .[[1]] %>% length %>% seq(1, .) %>% as.character()
updateSelectizeInput(session, 'pt_file_grp',
server = TRUE,
choices = choices,
selected = choices[1]
)
})
observeEvent(input$pt_file_grp, {
# this should change with respect to the selection of filename group
choices <- rawdata() %>% colnames() %>% strsplit("_") %>%
sapply(., function(x) x[as.numeric(input$pt_file_grp)]) %>% unique %>% length %>%
seq(1, .) %>% as.character()
updateSelectizeInput(session, 'pt_grp_order',
server = TRUE,
choices = choices,
selected = choices[1]
)
})
observeEvent(transform_data$data, {
updateSelectizeInput(session, 'pt_allgene',
server = TRUE,
choices = as.character(rownames(transform_data$data)),
selected = as.character(rownames(transform_data$data)[1])
)
})
observeEvent(input$runNMF,{
updateSelectizeInput(session, 'pt_col',
server = TRUE,
choices = as.character(c(plot_colopts(), "NMF Group", "NMF Feature")),
selected = 'NMF Group'
)
updateSelectizeInput(session, 'pt_nmfgene',
server = TRUE,
choices = as.character(nmf_features()$Gene),
selected = as.character(nmf_features()$Gene[1])
)
})
observeEvent(nmf_groups()$nmf_subtypes, {
choices <- nmf_groups()$nmf_subtypes %>% unique %>% length %>% seq(1, .) %>% as.character()
updateSelectizeInput(session, 'sel_grp_order',
server = TRUE,
choices = choices,
selected = choices[1]
)
})
point_col <- reactive({
# col <- toRGB("steelblue")
col <- "steelblue"
plot_data <- plot_data()[['plot_data']]
group <- NULL
if(input$pt_col == "Default"){
group <- "Default"
}else if(input$pt_col == "NMF Group"){
req(nmf_groups)
idx <- match(colnames(plot_data), nmf_groups()$Sample_ID)
nmf_subtypes <- nmf_groups()$nmf_subtypes[idx]
col <- create.brewer.color(nmf_subtypes, length(unique(nmf_subtypes)), "naikai")
group <- paste0("NMF", nmf_subtypes)
}else if(input$pt_col == "NMF Feature"){
req(input$pt_nmfgene)
idx <- match(colnames(plot_data), colnames(transform_data$data))
gene_data <- transform_data$data[input$pt_nmfgene, idx] %>% as.numeric
col <- create.brewer.color(gene_data, num = 9, name="YlOrRd")
group <- input$pt_nmfgene
}else if(input$pt_col == "Filename"){
filename <- colnames(plot_data)
filename <- sapply(strsplit(filename, "_"), function(x) paste0(x[as.numeric(input$pt_file_grp)], collapse = "_"))
col <- create.brewer.color(filename, length(unique(filename)), "naikai") %>% as.matrix
group <- filename %>% as.matrix
# add option to manual select the order of coloring
if(input$pt_sel_grp_order == "Manual"){
col <- col[, 1]
col.num <- as.numeric(factor(col, levels = unique(col) %>% as.character))
col.lvl <- unique(col)
col.lvl <- col.lvl[as.numeric(input$pt_grp_order)]
col <- col.lvl[col.num] %>% as.matrix
}
}else if(input$pt_col == "GeneExpr"){
req(input$pt_allgene)
idx <- match(colnames(plot_data), colnames(transform_data$data))
gene_data <- transform_data$data[input$pt_allgene, idx] %>% as.numeric
col <- create.brewer.color(gene_data, num = 9, name="YlOrRd")
group <- input$pt_allgene
}else if(input$pt_col == "ReadDepth"){
req(rawdata())
idx <- match(colnames(plot_data), colnames(rawdata()))
gene_data <- rawdata()[, idx] %>% colSums() %>% as.numeric
col <- create.brewer.color(gene_data, num = 9, name="YlOrRd")
group <- "Read Depth"
}else if(input$pt_col == "NumExpresGenes"){
req(rawdata())
idx <- match(colnames(plot_data), colnames(rawdata()))
gene_data <- colSums(rawdata()[, idx]>0) %>% as.numeric
col <- create.brewer.color(gene_data, num = 9, name="YlOrRd")
group <- "Num Expressed Genes"
}else{
warning("Wrong point color assigning method!")
}
res <- list()
res[['color']] <- col
res[['group']] <- group
return(res)
})
pca_data <- reactive({
data <- plot_data()[['plot_data']]
withProgress(message = 'Calculating PCs', value = NULL, {
pc <- prcomp(t(data), center = TRUE)
vars<- pc$sdev^2
vars<- signif(vars/sum(vars) * 100, 2)
title = paste(paste0("PC", 1:length(vars)), paste0("(", vars, "% Var)"))
projection <- as.data.frame(pc$x)
})
res <- list()
res[['proj']] <- projection
res[['title']] <- title
return(res)
})
#' PCA plot
output$pcaPlot_2D <- renderPlotly({
projection <- pca_data()[['proj']]
title <- pca_data()[['title']]
validate(
need(as.numeric(input$pca_x) != as.numeric(input$pca_y), "PC on the X-axis is the same as PC on the Y-axis.\nPlease select different number for these two axes")
)
idx <- c(as.numeric(input$pca_x), as.numeric(input$pca_y))
withProgress(message = 'Generating 2d PCA plot', value = NULL, {
projection <- projection %>%
dplyr::select(idx) %>%
set_colnames(c("pca_x", "pca_y")) %>%
mutate(color = point_col()[['color']] %>% as.character(),
group = point_col()[['group']] %>% as.character(),
Sample = colnames(plot_data()[['plot_data']])) %>%
arrange(group)
if(input$plot_label){
t <- list( size=input$plot_label_size, color=toRGB("grey50") )
p <- plot_ly(projection, x = ~pca_x, y = ~pca_y, type="scatter", mode="markers+text",
color = ~group, colors = ~unique(color),
text = ~Sample, hoverinfo="text", textfont=t, textposition="top middle",
marker = list(color = ~color, size=input$plot_point_size+3, opacity=input$plot_point_alpha))
}else{
p <- plot_ly(projection, x = ~pca_x, y = ~pca_y, type="scatter", mode="markers",
color = ~group, colors = ~unique(color),
text = ~Sample, hoverinfo="text",
marker = list(color = ~color, size=input$plot_point_size+3, opacity=input$plot_point_alpha)) #%>%
# toWebGL()
}
p %>% layout(showlegend = input$plot_legend,
xaxis = list(title = title[idx[1]], zeroline=FALSE),
yaxis = list(title = title[idx[2]], zeroline=FALSE))
})
})
output$pcaPlot_3D <- renderPlotly({
projection <- pca_data()[['proj']]
title <- pca_data()[['title']]
idx <- c(1,2,3)
withProgress(message = 'Generating 3d PCA plot', value = NULL, {
projection <- projection %>%
dplyr::select(idx) %>%
set_colnames(c("pca_x", "pca_y", "pca_z")) %>%
mutate(color = point_col()[['color']] %>% as.character(),
group = point_col()[['group']] %>% as.character(),
Sample = colnames(plot_data()[['plot_data']])) %>%
arrange(group)
if(input$plot_label){
p <- plot_ly(data=projection, x = ~pca_x, y = ~pca_y, z = ~pca_z, type="scatter3d", mode="markers+text",
color = ~group, colors = ~unique(color), text = ~Sample, hoverinfo="text",
marker = list(color = ~color, size=input$plot_point_size, opacity=input$plot_point_alpha))
}else{
p <- plot_ly(data=projection, x = ~pca_x, y = ~pca_y, z = ~pca_z, type="scatter3d", mode="markers",
color = ~group, colors = ~unique(color), text = ~Sample, hoverinfo="text",
marker = list(color = ~color, size=input$plot_point_size, opacity=input$plot_point_alpha)) #%>%
# toWebGL()
}
p %>% layout(showlegend = input$plot_legend,
scene = list(
xaxis = list(title = title[1]),
yaxis = list(title = title[2]),
zaxis = list(title = title[3]))
)
})
})
#' t-SNE plot
cores <- reactive({
avail_cores <- parallel::detectCores()
if(avail_cores==1){
return(1)
}else if(as.numeric(input$ncores) <= avail_cores){
return(as.numeric(input$ncores))
}else{
return(avail_cores)
}
})
tsne_2d <- reactiveValues(data=NULL)
tsne_3d <- reactiveValues(data=NULL)
observeEvent(input$runtSNE, {
nsamples <- ncol(plot_data()[['plot_data']])
if(input$tsne_perplexity*3 > (nsamples-1)) {
createAlert(session, "visualperplexityAlert", "visualperplexityAlert1", title = "WARNING", style = "warning",
content = paste("Perpleixty", input$tsne_perplexity,
"is too large compared to num of samples", nsamples),
append = FALSE)
}else {
closeAlert(session, "visualperplexityAlert1")
plot_data <- plot_data()[['plot_data']]
withProgress(message = 'Running t-SNE', value = NULL, {
incProgress(1/3, detail = "For t-SNE 2D")
tsne_2d$data <- run_tsne(plot_data, iter=input$tsne_iter, dims=2,
perplexity=input$tsne_perplexity, cores=cores())
incProgress(2/3, detail = "For t-SNE 3D")
tsne_3d$data <- run_tsne(plot_data, iter=input$tsne_iter, dims=3,
perplexity=input$tsne_perplexity, cores=cores())
})
}
})
plot_tsne_2d <- reactive({
nsamples <- ncol(plot_data()[['plot_data']])
validate(
need(!is.null(tsne_2d$data), "Please click 'Run t-SNE' button")
)
tsne_out <- isolate(tsne_2d$data)
withProgress(message = 'Genrating t-SNE plot', value = 0, {
incProgress(2/3, detail = "Usually takes 15~20 seconds")
color <- "steelblue"
projection <- parse_tsne_res(tsne_out) %>%
data.frame %>%
mutate(group = point_col()[['group']] %>% as.character(),
color = point_col()[['color']] %>% as.character(),
Sample = rownames(.)) %>%
arrange(group)
min.cost <- signif(tsne_out$itercosts[length(tsne_out$itercosts)], 2)
if(input$plot_label){
t <- list( size=input$plot_label_size, color=toRGB("grey50") )
p <- plot_ly(projection, x = ~x, y = ~y, type="scatter", mode="markers+text",
color = ~group, colors = ~unique(color),
text= ~Sample, hoverinfo="text", textposition="top middle", textfont=t,
marker = list(color = ~color, size=input$plot_point_size+3, opacity=input$plot_point_alpha))
}else{
p <- plot_ly(projection, x = ~x, y = ~y, type="scatter", mode="markers", text= ~Sample, hoverinfo="text",
color = ~group, colors = ~unique(color),
marker = list(color = ~color, size=input$plot_point_size+3, opacity=input$plot_point_alpha)) #%>%
# toWebGL()
}
p %>% layout(showlegend = input$plot_legend,
title = title,
xaxis = list(title = "Component 1", zeroline=FALSE),
yaxis = list(title = "Component 2", zeroline=FALSE))
})
})
output$tsneplot_2d <- renderPlotly({
plot_tsne_2d()
})
output$tsneplot_3d <- renderPlotly({
nsamples <- ncol(plot_data()[['plot_data']])
if(is.null(tsne_3d$data)) return ()
projection <- isolate(as.data.frame(tsne_3d$data$Y)) %>%
set_colnames(c("tsne_1", "tsne_2", "tsne_3")) %>%
mutate(group = point_col()[['group']] %>% as.character(),
color = point_col()[['color']] %>% as.character(),
Sample = rownames(.)) %>%
arrange(group)
if(input$plot_label){
p <- plot_ly(data=projection, x = ~tsne_1, y = ~tsne_2, z = ~tsne_3, type="scatter3d", mode="markers+text",
color = ~group, colors = ~unique(color), text= ~Sample, hoverinfo="text",
marker = list(color = ~color, size=input$plot_point_size, opacity=input$plot_point_alpha))
}else{
p <- plot_ly(data=projection, x = ~tsne_1, y = ~tsne_2, z = ~tsne_3, type="scatter3d", mode="markers",
color = ~group, colors = ~unique(color), text= ~Sample, hoverinfo="text",
marker = list(color = ~color, size=input$plot_point_size, opacity=input$plot_point_alpha)) %>%
toWebGL()
}
p %>% layout(showlegend = input$plot_legend,
scene = list(
xaxis = list(title = "Component 1"),
yaxis = list(title = "Component 2"),
zaxis = list(title = "Component 3"))
)
})
### DESeq2 ###
observeEvent(input$runNMF, {
validate(
need(is.numeric(input$num_cluster), "num of cluster is not a valid numeric number")
)
updateSelectizeInput(session, 'de_group1',
server = TRUE,
choices = as.character(paste0("NMF", sort(unique(nmf_groups()$nmf_subtypes)))),
selected = "NMF1"
)
updateSelectizeInput(session, 'de_group2',
server = TRUE,
choices = as.character(paste0("NMF", sort(unique(nmf_groups()$nmf_subtypes)))),
selected = "NMF2"
)
})
deseq_res <- eventReactive(input$runDESeq, {
if(input$selectfile == "saved"){
if(!is.null(rda()$dds)){
return(rda()$dds)
}
}
rawdata <- rawdata()
if(input$de_conv2int)
rawdata <- round(rawdata)
if(!all(sapply(rawdata, is.whole))){
createAlert(session, "alert", "exampleAlert", title = "WARNING", style = "warning",
content = "Data contains non-integer value, Please use raw count data for DESeq2. <br>
Or Use the option above to 'Force convert to integer' <b>(Results are not statistically rigorous)</b>",
append = FALSE)
return(NULL)
}else{
closeAlert(session, "exampleAlert")
}
# since the user might have filtered samples during the QC steps,
# rematch sample_ID and only kept the ones that are presented in NMF res
idx <- match(nmf_groups()$Sample_ID, colnames(rawdata))
rawdata <- rawdata[, idx]
register(BiocParallel::MulticoreParam(workers = cores()))
colData <- data.frame(Group = paste0("NMF", nmf_groups()$nmf_subtypes))
ddsfeatureCounts <- DESeq2::DESeqDataSetFromMatrix(countData = rawdata,
colData = colData,
design = ~ Group)
### Set betaPrior=FALSE to go with MLE LFC to get simple LFC = (avg in group2/ avg in group1)
# dds <- DESeq(ddsfeatureCounts, parallel=T, betaPrior=FALSE)
withProgress(message = 'Running DESeq2', value = NULL, {
dds <- DESeq2::DESeq(ddsfeatureCounts, parallel=T)
})
return(dds)
})
filt_deseq_res <- reactive({
req(deseq_res())
withProgress(message = 'Summarizing DESeq2 result', value = NULL, {
register(BiocParallel::MulticoreParam(workers = 8))
DESeq2::results(deseq_res(), parallel = TRUE,
contrast = c("Group", input$de_group1, input$de_group2)) %>%
as.data.frame %>%
subset(padj <= input$de_alpha) %>%
.[order(.$padj), ]
})
})
output$deseq_table <- DT::renderDataTable({
validate(
need(nrow(filt_deseq_res()) > 0, "No significant DE results found, Maybe try lower your alpha threshold?")
)
filename <- fileinfo()[['name']]
GeneCard <- paste0("<a href='http://www.genecards.org/cgi-bin/carddisp.pl?gene=",
rownames(filt_deseq_res()), "'>", rownames(filt_deseq_res()), "</a>")
filt_res <- data.frame(Gene=GeneCard,
filt_deseq_res())
DT::datatable(filt_res, selection='single', rownames=FALSE,
extensions = 'Buttons',
escape = FALSE,
options = list(dom = 'Bfrtip',
buttons =
list('colvis', 'copy', 'print', list(
extend = 'collection',
buttons = list(list(extend='csv',
filename = filename),
list(extend='excel',
filename = filename)),
text = 'Download'
)),
scrollY = TRUE,
pageLength = 12,
autoWidth = TRUE
)
) %>% formatRound(2:ncol(filt_res), 3)
}, server = FALSE)
output$deseq_boxplot <- renderPlotly({
req(deseq_res())
validate(
need(!is.null(input$deseq_table_rows_selected), "Please select a gene")
)
withProgress(message = 'Generating boxplot', value = NULL, {
norm.data <- DESeq2::counts(deseq_res(), normalized=TRUE)
gene <- rownames(filt_deseq_res())[input$deseq_table_rows_selected]
gene.data <- norm.data[gene, ]
gene.data <- data.frame(Sample = names(gene.data),
Expr = gene.data,
NMF = factor(paste0("NMF", nmf_groups()$nmf_subtypes)))
# match the coloring from PCA, t-SNE, and heatmap
num_clus <- gene.data$NMF %>% levels %>% length
mycolor <- create.brewer.color(1:num_clus, num=num_clus, name="naikai")
a <- ggplot(data=gene.data, aes(x=NMF, y=Expr, color=NMF)) +
geom_boxplot(aes(color = NMF)) +
geom_jitter(aes(color=NMF, text=Sample), alpha=0.6, width=0.1) +
theme_bw() +
scale_colour_manual(values = mycolor) +
labs(title=gene, x="", y="Expression", colour="")
p <- plotly_build(a)
# remove outliers for plotly boxplot
p$data <- lapply(p$data, FUN = function(x){
if(x$type == "box"){
x$marker = list(opacity = 0)
}
return(x)
})
p$layout$margin$l <- p$layout$margin$l + 15
p$layout$annotations[[1]]$x <- -0.1
p$layout$xaxis$tickfont$size <- 8
p$layout$annotations[[1]]$font$size <- 11
p$layout$yaxis$tickfont$size <- 10
p
})
})
### Pathway analysis ###
nmf_group_feature_rank <- reactive({
validate(
need(is.numeric(input$num_cluster), "num of cluster is not a numeric value") %then%
need(!is.null(nmf_groups()), "Please run NMF first")
)
rawdata <- transform_data$data
# rawdata <- rawdata[rowMeans(rawdata) >= input$min_rowMean, ]
nmf_group <- nmf_groups()$nmf_subtypes
num_nmf_subtypes <- nmf_group %>% unique %>% length
if(input$rank_method=="featureScore"){
print('Do we want to use featureScore, there are only few genes')
}else if(input$rank_method=="logFC"){
nmf_feature_rank <- data.frame(matrix(NA_real_, nrow=nrow(rawdata), ncol=num_nmf_subtypes))
for(i in 1:num_nmf_subtypes){
idx <- nmf_group==i
nmf_feature_rank[, i] <- rawdata %>% as.data.frame %>%
dplyr::mutate(logFC = log(rowMeans(.[idx])/rowMeans(.[!idx]))) %>%
dplyr::select(logFC)
}
colnames(nmf_feature_rank) <- paste0("NMF", 1:num_nmf_subtypes)
rownames(nmf_feature_rank) <- rownames(rawdata)
}
return(nmf_feature_rank)
})
geneset.file <- reactive({
validate(
need(!is.null(input$predefined_list), "Please select at least one predefined gene sets")
)
file <- pre_geneset[[ input$predefined_list ]]
return (file)
})
ncpus <- reactive({
cpus <- c(16,8,6,4,2,1)
cpus[min(which(sapply(cpus, function(x) input$nPerm%%x)==0))]
})
observeEvent(input$runNMF, {
updateSelectizeInput(session, 'pathway_group',
server = TRUE,
choices = as.character(paste0("NMF", sort(unique(nmf_groups()$nmf_subtypes))))
)
})
pathview.species <- reactive({
species <- input$species
if(species == "Human" | species == "human" | species == "hg19" | species == "hsa"){
pathview.species <- "hsa"
}else if(species == "Mouse" | species == "mouse" | species == "mm9" | species == "mmu"){
pathview.species <- "mmu"
}else if(species == "Drosophila" | species == "dm3" | species == "dm6"){
pathview.species <- "dm3"
}
})
id.org <- reactive({
species <- input$species
if(species == "Human" | species == "human" | species == "hg19" | species == "hsa"){
id.org <- "Hs"
}else if(species == "Mouse" | species == "mouse" | species == "mm9" | species == "mmu"){
id.org <- "Mm"
}
})
kegg.gs <- reactive({
species <- input$species
withProgress(message = 'Extracting KEGG pathway information', value = NULL, {
if(species == "Human" | species == "human" | species == "hg19" | species == "hsa"){
kg.human<- kegg.gsets("human")
kegg.gs<- kg.human$kg.sets[kg.human$sigmet.idx]
}else if(species == "Mouse" | species == "mouse" | species == "mm9" | species == "mmu"){
kg.mouse<- kegg.gsets("mouse")
kegg.gs<- kg.mouse$kg.sets[kg.mouse$sigmet.idx]
}
})
return(kegg.gs)
})
go.gs <- reactive({
species <- input$species
withProgress(message = 'Extracting GO term information', value = NULL, {
if(species == "Human" | species == "human" | species == "hg19" | species == "hsa"){
go.gs <- go.gsets(species="human", pkg.name=NULL, id.type = "eg")
}else if(species == "Mouse" | species == "mouse" | species == "mm9" | species == "mmu"){
go.gs <- go.gsets(species="mouse", pkg.name=NULL, id.type = "eg")
}
})
return(go.gs)
})
gage.exp.fc <- eventReactive(input$run_enrich, {
rankdata <- nmf_group_feature_rank()
num_groups <- ncol(rankdata)
emp <- matrix(NA, nrow(rankdata), 1)
res <- rep(list(emp), num_groups)
names(res) <- colnames(rankdata)
withProgress(message = 'Calculating LogFC between groups', value = NULL, {
for (i in 1:num_groups){
# In logFC, error might come from the fact that some of them are Inf or -Inf
idx <- is.finite(rankdata[, i])
sub_rankdata <- rankdata[idx, i, drop=F]
# ID conversion
id.map.sym2eg <- id2eg(ids=rownames(sub_rankdata), category = "SYMBOL", org=id.org())
gene.entrez <- mol.sum(mol.data = sub_rankdata, id.map = id.map.sym2eg, sum.method = "mean")
res[[i]] <- gene.entrez[, 1]
res[[i]] <- res[[i]][!is.na(res[[i]])]
}
})
return(res)
})
go.Res <- reactive({
gage.exp.fc <- gage.exp.fc()
gsets <- NULL
res <- rep(list(NA), length(gage.exp.fc))
names(res) <- names(gage.exp.fc)
if(input$EnrichType == 'GO'){
gsets <- go.gs()
gsets <- gsets$go.sets[ gsets$go.subs[[input$goTerm]] ]
}else if(input$EnrichType == 'KEGG'){
gsets <- kegg.gs()
}
validate(
need(length(gsets)>0, "Did not find any Gene Sets, maybe your gene count data is not from 'Human'?")
)
withProgress(message = 'Running enrichment analysis', value = NULL, {
for (i in 1:length(gage.exp.fc)){
exp.fc <- gage.exp.fc[[i]]
goRes <- gage(exp.fc, gsets = gsets, ref = NULL, samp = NULL, same.dir=ifelse(input$samedir, TRUE, FALSE))
res[[i]] <- goRes
}
})
return(res)
})
go_summary <- reactive({
goRes <- go.Res()
genes <- goRes[[1]]$greater %>% rownames
res <- list()
num_genes <- 5
hi_res <- lo_res <- matrix(NA, num_genes * length(goRes), length(goRes))
colnames(hi_res) <- colnames(lo_res) <- names(goRes)
total_hi_names <- total_lo_names <- rep(NA, num_genes * length(goRes))
for(i in 1:length(goRes)){
range <- ((i-1) * num_genes + 1) : (i * num_genes)
hi_genenames <- goRes[[i]]$greater %>% rownames %>% head(n = num_genes)
hi_res[range, ] <- sapply(goRes, function(x) x$greater[hi_genenames, "p.val"])
total_hi_names[range] <- hi_genenames
if(!is.na(goRes[[i]]$less[1])){
lo_genenames <- goRes[[i]]$less %>% rownames %>% head(n = num_genes)
lo_res[range, ] <- sapply(goRes, function(x) x$less[lo_genenames, "p.val"])
total_lo_names[range] <- lo_genenames
}
}
rownames(hi_res) <- total_hi_names
rownames(lo_res) <- total_lo_names
# save only the unique GO:Term
res[["hi_sub_pval"]] <- hi_res[!duplicated(rownames(hi_res)), ]
res[["lo_sub_pval"]] <- lo_res[!duplicated(rownames(lo_res)), ]
# save total results
res[["hi_total_pval"]] <- sapply(goRes, function(x) x$greater[genes, "p.val"])
res[["hi_total_qval"]] <- sapply(goRes, function(x) x$greater[genes, "q.val"])
res[["lo_total_pval"]] <- sapply(goRes, function(x) x$less[genes, "p.val"])
res[["lo_total_qval"]] <- sapply(goRes, function(x) x$less[genes, "q.val"])
return(res)
})
output$goplot_hi <- renderPlotly({
req(go_summary())
validate(
need(!is.na(go_summary()[["hi_sub_pval"]]), "Did not find any GO terms for greater, maybe you select the wrong 'species'?")
)
withProgress(message = 'Generating summary plot for GO term (Hi)', value = NULL, {
aa <- isolate(
melt(go_summary()[["hi_sub_pval"]]) %>%
set_colnames(c("GO", "NMF", "P")) %>%
mutate(P = replace(P, is.na(P), 1))
)
if(input$go_logscale){
p <- plot_ly(aa, x = ~NMF, y = ~GO, z = ~log10(P), type = "heatmap", colors = input$go_colorscale)
}else{
p <- plot_ly(aa, x = ~NMF, y = ~GO, z = ~P, type = "heatmap", colors = input$go_colorscale)
}
m <- list( l = input$go_leftmargin, r = 10, b = 40, t = 20)
l <- list(margin = m,
xaxis = list(title = "", tickfont=list(family = "Arial", size = input$go_xfontsize), zeroline=FALSE),
yaxis = list(title = "", tickfont=list(family = "Arial", size = input$go_yfontsize), zeroline=FALSE))
colorbar_format <- list(family = "Arial", size = 8)
naikai <- plotly_build(p)
naikai$x$data[[1]]$colorbar <- list(title = "P", thickness = input$go_barwidth, len=0.5, titlefont=colorbar_format, tickfont=colorbar_format)
naikai$x$layout <- l
naikai
})
})
output$goplot_lo <- renderPlotly({
req(go_summary())
validate(
need(!is.na(go_summary()[["lo_sub_pval"]]), "Did not find any GO terms for less, maybe you select the wrong 'species'?")
)
withProgress(message = 'Generating summary plot for GO term (Low)', value = NULL, {
aa <- isolate(
melt(go_summary()[["lo_sub_pval"]]) %>%
set_colnames(c("GO", "NMF", "P")) %>%
mutate(P = replace(P, is.na(P), 1))
)
if(input$go_logscale){
p <- plot_ly(aa, x = ~NMF, y = ~GO, z = ~log10(P), type = "heatmap", colors = input$go_colorscale)
}else{
p <- plot_ly(aa, x = ~NMF, y = ~GO, z = ~P, type = "heatmap", colors = input$go_colorscale)
}
m <- list( l = input$go_leftmargin, r = 10, b = 40, t = 20)
l <- list(margin = m,
xaxis = list(title = "", tickfont=list(family = "Arial", size = input$go_xfontsize), zeroline=FALSE),
yaxis = list(title = "", tickfont=list(family = "Arial", size = input$go_yfontsize), zeroline=FALSE))
colorbar_format <- list(family = "Arial", size = 8)
naikai <- plotly_build(p)
naikai$x$data[[1]]$colorbar <- list(title = "P", thickness = input$go_barwidth, len=0.5, titlefont=colorbar_format, tickfont=colorbar_format)
naikai$x$layout <- l
naikai
})
})
go_table <- reactive({
goRes <- go.Res()
validate(
need(input$pathway_group!="", "Please select the enrichment result for the group you are interest in")
)
goRes <- goRes[[input$pathway_group]]
if (is.null(goRes$greater)){
stop("No significant gene sets found using this q.val cutoff")
}else{
return(goRes$greater)
}
})
output$go_summary <- DT::renderDataTable({
filename <- paste(input$EnrichType, input$pathway_group, input$goTerm, sep = "_")
idx <- go_table()[, "q.val"] <= input$qval_cutoff
go_table <- go_table()[which(!is.na(idx) & idx), ]
if(nrow(go_table) == 0){
stop("No significant gene sets found using this q.val cutoff")
}
go_table <- go_table[, c(3,4,2,5)]
datatable(go_table, selection = 'single',
extensions = c('Buttons', 'Responsive'),
options = list(dom = 'Bfrtip',
buttons =
list('colvis', 'copy', 'print', list(
extend = 'collection',
buttons = list(list(extend='csv',
filename = filename),
list(extend='excel',
filename = filename)),
text = 'Download'
)),
pageLength = 6,
autoWidth = TRUE
# columnDefs = list(list(width = '10px', targets = "_all"))
)
) %>% formatRound(1:3, 3)
}, server = FALSE)
output$pathviewImg <- renderImage({
go_table <- go_table()
s.idx <- input$go_summary_rows_selected[length(input$go_summary_rows_selected)]
validate(
need(length(s.idx)>0, "Please select one KEGG pathway from above.")
)
withProgress(message = 'Rendering pathview map', value = NULL, {
path.pid <- substr(rownames(go_table)[s.idx], 1, 8)
out.suffix <- "nmf"
pathview(gene.data=gage.exp.fc()[[input$pathway_group]], pathway.id=path.pid, species=pathview.species(), out.suffix=out.suffix)
})
filename <- path.pid
outfile <- paste(filename, out.suffix, 'png', sep=".")
# Return a list containing the filename
list(src = outfile,
contentType = 'image/png',
width = 1000,
height = 777,
alt = "This is alternate text")
}, deleteFile = FALSE)
nmf_feature_rank_for_gsea <- reactive({
species <- input$species
nmf_feature_rank <- nmf_group_feature_rank()
# check if genenames contains all CAPITAL LETTERS (after removing all numeric characters)
cap_genes <- grepl("^[[:upper:]]+$", gsub('[[:digit:]]+', '', rownames(nmf_feature_rank))) %>% sum
cap_pct <- cap_genes/ nrow(nmf_feature_rank)
# some of the genes might contains '.', '-', 'orf', etc, so we use percentage of genes qualified as filtering criteria
# if(species == "Mouse" | species == "mouse" | species == "mm9" | species == "mmu"){
if(cap_pct < 0.1){
n <- 2
withProgress(message = 'Assuming the input gene symbols are from Mouse, Converting to Human..', value = 0, {
incProgress(1/n, detail = "Takes around 10~15 seconds")
mouse_human_convert_genes <- gene_convert_mouse_to_human(rownames(nmf_feature_rank))
## After conversion, some of the mouse genes might point to the same human gene
## So select thoide mouse genes that has highest mean expression across groups?
idx <- duplicated(mouse_human_convert_genes$HGNC.symbol)
if (sum(idx)>0){
warning("Some of the mouse genes are mapped back to the same human genes..")
}
max_dup_mean_lfc_genes <-
cbind(mouse_human_convert_genes, nmf_feature_rank) %>%
mutate(avg_lfc=apply(nmf_feature_rank, 1, mean),
abs_lfc=apply(nmf_feature_rank, 1, function(x) max(abs(x)))
) %>%
dplyr::group_by(HGNC.symbol) %>%
dplyr::filter(min_rank(desc(avg_lfc))==1) %>%
dplyr::slice(which.max(abs_lfc)) %>%
ungroup # this will give Error: corrupt 'grouped_df', workaround is to ungroup it below
print(dim(max_dup_mean_lfc_genes))
## Return final de-duplicated genes only
nmf_feature_rank <- max_dup_mean_lfc_genes %>%
dplyr::select(starts_with("NMF")) %>%
as.data.frame
rownames(nmf_feature_rank) <- max_dup_mean_lfc_genes$HGNC.symbol
})
}
return(nmf_feature_rank)
})
output$downloadPathData <- downloadHandler(
filename <- function() {
paste( file_prefix(), "pathway.zip", sep=".")
},
content = function(file) {
tmpdir <- tempdir()
current_dir <- getwd()
setwd(tempdir())
### Save GO Term results
go_summary <- go_summary()
total_idx <- grep("total", names(go_summary))
go_filenames <- paste(file_prefix(), paste0("Enrichment_", names(go_summary)[total_idx]), "txt", sep=".")
for(i in 1:length(total_idx)){
write.table(go_summary[[ total_idx[i] ]], go_filenames[i], quote=F, row.names = TRUE, col.names=NA, sep="\t")
}
### Rank file for GSEA
# modify it to save one file for each rank ###
nmf_group_feature_rank <- nmf_feature_rank_for_gsea()
num_samples <- ncol(nmf_group_feature_rank)
filenames <- paste(file_prefix(), paste0("group_feature_rank", 1:num_samples), "rnk", sep=".")
for (i in 1:num_samples){
sub.data <- data.frame(Gene=rownames(nmf_group_feature_rank), LogFC=nmf_group_feature_rank[,i])
idx <- is.finite(sub.data$LogFC)
write.table(sub.data[idx, ], filenames[i], quote=F, row.names = F, sep="\t")
}
zip(zipfile=file, files=c(go_filenames, filenames))
setwd(as.character(current_dir))
},
contentType = "application/zip"
)
})
naikai/sake documentation built on Feb. 15, 2023, 11 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(shiny)
library(shinyBS)
library(ggplot2)
library(data.table)
library(magrittr)
library(corrplot)
library(biomaRt)
library(BiocParallel)
library(AnnotationHub)
library(RColorBrewer)
library(matrixStats)
library(tools)
library(snowfall)
library(NMF)
library(DT)
library(gage)
library(gageData)
library(pathview)
shinyServer(function(input, output, session) {
fileinfo <- reactive({
if (input$selectfile == 'upload'){
inFile <- input$file1
req(inFile)
validate(
need(!is.null(inFile), "Please upload a gene count data set")
)
filename <- basename(inFile$name)
filepath <- inFile$datapath
}else if (input$selectfile == "preload"){
req(input$inputdata)
validate(
need(input$inputdata!="", "Please select a gene count data set")
)
filename <- basename(input$inputdata)
filepath <- file.path(raw_fd, input$inputdata)
}else if (input$selectfile == "saved"){
inFile <- input$rda
req(inFile)
validate(
need(!is.null(inFile), "Please upload a .rda file")
)
filename <- basename(inFile$name)
filepath <- inFile$datapath
}else{
return()
}
fileinfo <- list()
fileinfo[['name']] <- filename
fileinfo[['path']] <- filepath
return(fileinfo)
})
### Output Prefix ###
File <- reactive({
if (input$selectfile == 'upload'){
File <- input$file1$name
}else if (input$selectfile == "preload"){
File <- input$inputdata
}else if (input$selectfile == "saved"){
File <- input$rda$name
}else{
return()
}
return(File)
})
File_ext <- reactive({
file_ext(File())
})
file_prefix <- reactive({
if(input$selectfile == "saved"){
gsub(paste0(".", File_ext()), "", File())
}else{
paste( gsub(paste0(".", File_ext()), "", File()),
paste0("k", input$num_cluster),
paste0(input$top.num,input$algorithm),
paste0("nrun", input$nrun),
paste0("predict_", input$predict),
sep=".")
}
})
rda <- reactive({
req(fileinfo())
req(input$selectfile == "saved")
filepath <- fileinfo()[['path']]
validate(
need(File_ext() == "rda" || File_ext() == "RData",
message = paste("file extension", File_ext(), "unknown, please try again"))
)
rda <- local({load(filepath); environment()})
# start checking for all the required variables
validate(
need(!is.null(rda$rawdata), "Variable 'rawdata' does not exist in your .rda file, please check it and upload again") %then%
need(nrow(rda$rawdata) > 0, "Variable 'rawdata' in your .rda file is empty, please check it and upload again") %then%
need(!is.null(rda$nmfres), "Variable 'nmfres' does not exist in your .rda file, please check it and upload again")
# need(!is.null(rda$tsne_2d), "Variable 'tsne_2d' does not exist in your .rda file, please check it and upload again") %then%
# need(!is.null(rda$tsne_3d), "Variable 'tsne_3d' does not exist in your .rda file, please check it and upload again") #%then%
# need(!is.null(rda$dds), "Variable 'dds' does not exist in your .rda file, please check it and upload again")
)
return(rda)
})
observeEvent(rda(), {
tsne_2d$data <- rda()$tsne_2d
tsne_3d$data <- rda()$tsne_3d
})
rawdata <- reactive({
req(fileinfo())
filepath <- fileinfo()[['path']]
n <- 2
withProgress(message = 'Loading data', value = 0, {
incProgress(1/n, detail = "Usually takes ~10 seconds")
if(input$selectfile == "saved"){
rawdata <- rda()$rawdata
}else{
rawdata <- myfread.table(filepath, check.platform=T, sep=input$sep, detect.file.ext=FALSE)
}
})
colnames(rawdata) <- gsub("\\.", "-", colnames(rawdata))
return (rawdata)
})
#' Raw data
output$rawdatatbl = DT::renderDataTable({
rawdata <- rawdata()
validate(
need(ncol(rawdata)>1,
"Please modify the parameters on the left and try again!") %then%
need(!any(duplicated(t(rawdata))),
"There are columns with the exact same value in your data, Please check and try again!")
)
withProgress(message = 'Displying datatable', value = NULL, {
DT::datatable(rawdata[1:20, 1:min(500, ncol(rawdata))], rownames= TRUE,
options = list(scroll = TRUE,
scrollX = TRUE,
scrollY = TRUE,
pageLength = 8
)
) %>% formatRound(1:ncol(rawdata), 3)
})
}, server=TRUE)
#' Info box
output$nameBox <- renderInfoBox({
infoBox(
"FILENAME", fileinfo()[['name']], icon = icon("list"),
color = "green"
)
})
output$sampleBox <- renderInfoBox({
rawdata <- rawdata()
infoBox(
"SAMPLE", ncol(rawdata), icon = icon("list"),
color = "blue"
)
})
output$featureBox <- renderInfoBox({
rawdata <- rawdata()
infoBox(
"FEATURE", nrow(rawdata), icon = icon("list"),
color = "yellow"
)
})
# transform_data <- reactiveValues(data=NULL, rawdata=NULL)
transform_data <- reactiveValues(data=data.frame(a=NULL, b=NULL),
rawdata=data.frame(a=NULL, b=NULL))
run_selsamp <- reactiveValues(go=FALSE)
run_trans <- function(rawdata){
transdata <- rawdata
n <- 2
withProgress(message = 'Transforming data', value = 0, {
incProgress(1/n, detail = "Takes around 5~10 seconds")
if(input$normdata=="takeRPM"){
transdata <- rpm(transdata, colSums(transdata))
}else if(input$normdata=="takesizeNorm"){
sf <- norm_factors(transdata)
transdata <- t(t(transdata)/sf)
}else if(input$normdata=="takeuq"){
transdata <- uq(transdata)
}
transdata <- rmv_constant_0(transdata, pct=0.95, minimum=0)
if(input$transformdata=="takevst"){
transdata <- vst(round(transdata))
}else if(input$transformdata=="takelog"){
validate(
need(!any(transdata<0), "Please make sure all the numbers in the data are positive values")
)
transdata <- log2(transdata + 1)
}
})
return(transdata)
}
observeEvent(rawdata(), {
transdata <- run_trans(rawdata())
transform_data$data <- transdata
transform_data$rawdata <- rawdata()
})
observeEvent(input$run_transf, {
transdata <- run_trans(rawdata())
transform_data$data <- transdata
transform_data$rawdata <- rawdata()
})
output$dl_transformdata <- downloadHandler(
filename <- function() {
paste(file_prefix(), "filtered.txt", sep=".")
},
content = function(file) {
data <- transform_data$data
colnames(data) <- gsub("\\.", "-", colnames(data))
write.table(data, file = file, quote=F, row.names=T, sep="\t")
}
)
output$readDistrib <- renderPlotly({
validate(
need(dim(transform_data$data)[1] > 1, "Did you upload your file or select preloaded data yet?")
)
total <- colSums(transform_data$data)
genecov <- colSums(transform_data$rawdata > 0)
run_selsamp$go <- TRUE
textlab <- paste(names(total),
paste("Total:", signif(total/1000000,3), "Million"),
paste("Gene Cov:", genecov),
sep="<br>")
plot_ly(x=genecov, y=total, text=textlab, source="readdist", #type="scattergl",
mode="markers", hoverinfo="text",
themes="Catherine", marker = list(size=10, opacity=0.65)) %>%
layout(yaxis = list(title = "Total transcripts counts", zeroline=TRUE),
xaxis = list(title = "Gene Coverage", zeroline=TRUE),
dragmode = "select",
showlegend=FALSE)
})
sel_samp <- reactive({
# Get subset based on selection
event.data <- event_data("plotly_selected", source="readdist")
# If NULL dont do anything
if(is.null(event.data) == T) return(NULL)
if(length(event.data$curveNumber) == 0) return(NULL)
# Get selected samples
sel_idx <- subset(event.data, curveNumber == 0)$pointNumber + 1
non_sel_idx <- subset(event.data, curveNumber == 0)$pointNumber + 1
if(max(sel_idx) > ncol(transform_data$data)) return(NULL)
res <- data.frame("Total_Reads"=colSums(transform_data$data)[sel_idx],
"Gene_Coverage"=colSums(transform_data$rawdata>0)[sel_idx])
return(res)
})
output$selsamp_tb <- DT::renderDataTable({
DT::datatable(sel_samp(),
rownames= TRUE,
options = list(scrollX = TRUE,
scrollY = TRUE,
dom = 'tipr',
pageLength = 5
)
)
}, server=TRUE)
output$selsamp_txt <- renderPrint({
cat("Selected samples\n\n")
idx <- input$selsamp_tb_rows_selected
if(length(idx) > 0){
cat(paste(rownames(sel_samp())[idx], sep = ', '))
}else{
cat(NULL)
}
})
### Filter selected samples from user
observeEvent(input$run_filtsamp, {
idx <- input$selsamp_tb_rows_selected
filtsamp <- rownames(sel_samp())[idx]
if(!is.null(filtsamp)){
idx <- match(filtsamp, colnames(transform_data$data))
if(length(idx) > 0){
temp <- transform_data$data[, -idx]
transform_data$data <- temp
temp <- transform_data$rawdata[, -idx]
transform_data$rawdata <- temp
}else{
warning("Can not find selected sample names in transform_data")
}
}
})
merged <- reactive({
rawdata <- transform_data$data
genelist <- NULL
if(!is.null(input$selected_samples) && sum(input$selected_samples %in% colnames(rawdata))==length(input$selected_samples)){
rawdata <- rawdata[, input$selected_samples]
}
if(input$list_type=='Top Ranks'){
select.genes <- select_top_n(apply(rawdata,1,input$math), n=input$top.num, bottom=ifelse(input$orders=="bottom", T, F))
genelist <- names(select.genes)
}else if(input$list_type=='Upload gene list'){
validate(
need(!is.null(input$genefile), "Please upload a file with list of genes, with column name 'Gene' ")
)
genefile <- read.table(input$genefile$datapath, header=T, sep="\t")
ext <- file_ext(input$genefile[1])
n <- 2
withProgress(message = 'Loading gene list', value = 0, {
incProgress(1/n, detail = "Usually takes ~10 seconds")
if (ext=="csv"){
genefile <- read.table(input$genefile$datapath, header=T, sep =',', stringsAsFactors = FALSE)
}else if (ext=="out" || ext=="tsv" || ext=="txt"){
genefile <- read.table(input$genefile$datapath, header=T, sep ='\t', stringsAsFactors = FALSE)
}else{
print("File format doesn't support, please try again")
return(NULL)
}
})
validate(
need(sum(colnames(genefile)=='Gene'), "Please make sure genefile contains column name 'Gene'")
)
validate(
need(length(intersect(genefile$Gene, rownames(rawdata)))>=2, "We did not find enough genes that overlap with the expression data (<2). Please check again and provide a new list! ")
)
genelist <- genefile$Gene
}else if(input$list_type == 'Whole transcriptome'){
genelist <- rownames(rawdata)
}
idx <- unique(match(genelist, rownames(rawdata)))
merged <- rawdata[idx[!is.na(idx)], ]
return(merged)
})
output$filterdatatbl = DT::renderDataTable({
filter_data <- merged()
DT::datatable(filter_data, rownames= TRUE,
selection = 'single',
options = list(scrollX = TRUE,
scroller = TRUE,
pageLength = 5
)
) %>% formatRound(1:ncol(filter_data), 2)
}, server=TRUE)
callModule(feature, "sample", reactive({ merged() }))
callModule(network, "gene", reactive({ merged() }))
#' Scatter plot
observeEvent(transform_data$data, {
transform_data <- transform_data$data
updateSelectizeInput(session, 'cor_sample1',
server = TRUE,
choices = sort(as.character(colnames(transform_data))),
selected = sort(as.character(colnames(transform_data)))[1]
)
updateSelectizeInput(session, 'cor_sample2',
server = TRUE,
choices = sort(as.character(colnames(transform_data))),
selected = sort(as.character(colnames(transform_data)))[2]
)
})
scatterData <- reactive({
transform_data <- transform_data$data
validate(
need(input$cor_sample1 != "", "Please select one sample for X axis"),
need(input$cor_sample2 != "", "Please select one sample for Y axis")
)
x <- transform_data[, input$cor_sample1]
y <- transform_data[, input$cor_sample2]
if(input$scatter_log){
if(input$transformdata == "takelog"){
message("Already did log transformation")
}else{
x <- log2(x+1)
y <- log2(y+1)
}
}
res <- data.frame(x=x, y=y)
return(res)
})
output$scatterPlot <- renderPlotly({
transform_data <- transform_data$data
x <- scatterData()$x %>% as.numeric
y <- scatterData()$y %>% as.numeric
plot_range <- c(0, min(x,y):max(x,y), max(x,y) + 0.2)
reg = lm(y ~ x)
modsum = summary(reg)
r <- signif(cor(x,y), 3)
R2 = signif(summary(reg)$r.squared, 3)
textlab <- paste(paste0("x = ", signif(x,3)), paste0("y = ", signif(y,3)), rownames(transform_data), sep="<br>")
n <- 2
withProgress(message = 'Generating Scatter Plot', value = 0, {
# p <- (x=~x, y=~y, text=~textlab, type="scattergl", #source = "scatter",
# mode="markers", hoverinfo="text",
# #themes="Catherine",
# marker = list(size=8, opacity=0.65)) %>%
p <- plot_ly() %>%
add_trace(x=~x, y=~y, text=~textlab, type="scattergl",
mode="markers", hoverinfo="text",
marker = list(size=8, opacity=0.65)) %>%
add_trace(x=~plot_range, y=~plot_range, type="scattergl", mode = "lines", name = "Ref", line = list(width=2)) %>%
layout(title = paste("cor:", r),
xaxis = list(title = input$cor_sample1, zeroline=TRUE),
yaxis = list(title = input$cor_sample2, zeroline=TRUE),
showlegend=FALSE)
})
# if(input$show_r2){
# p <- p %>% layout( annotations = list(x = x, y = y, text=paste0("R2=", R2), showarrow=FALSE) )
# }
if(input$show_fc){
df <- data.frame(x = plot_range[plot_range>=1], y=plot_range[plot_range>=1]-1)
p <- p %>%
add_trace(data=df, x= ~x, y= ~y, type = "scattergl", mode = "lines", name = "2-FC", line = list(width=1.5, color="#fec44f")) %>%
add_trace(data=df, x= ~y, y= ~x, type = "scattergl", mode = "lines", name = "2-FC", line = list(width=1.5, color="#fec44f"))
}
p #%>% toWebGL()
})
output$scatterdatatbl = DT::renderDataTable({
x <- scatterData()$x
y <- scatterData()$y
if(input$transformdata == "takelog"){
logFC = x - y
}else{
esp <- 0.001
logFC = log2((x+esp)/(y+esp))
}
GeneCard <- paste0("<a href='http://www.genecards.org/cgi-bin/carddisp.pl?gene=",
rownames(transform_data$data), "'>", rownames(transform_data$data), "</a>")
scatter_data <- data.frame(Gene=GeneCard,
X = x,
Y = y,
logFC = logFC,
meanFC = mean(x+y)*logFC) %>%
dplyr::arrange(desc(abs(logFC)))
DT::datatable(scatter_data, rownames= FALSE,
options = list(scrollX = TRUE,
pageLength = 8),
escape = FALSE
) %>% formatRound(2:5, 3)
}, server=TRUE)
### Main part of the program ###
nmf_seed <- eventReactive(input$runNMF, {
nmf_seed <- as.numeric(input$nmf_seed)
if(nmf_seed == 0){
set.seed(as.integer(Sys.time()))
sample(1:999999, 1)
}
})
nmf_res <- eventReactive(input$runNMF, {
merged <- merged()
validate(
need(is.numeric(input$num_cluster) & is.whole(input$num_cluster) & input$num_cluster>1, "Please provide a valid positive numeric integer value for 'Num of clusters(K)'") %then%
need(input$num_cluster <= ncol(merged), "Number of clusters(K) must be smaller than number of samples\nPlease Try selecting another K") %then%
need(is.numeric(input$nrun) & is.whole(input$nrun) & input$nrun>1, "Please provide a valid numeric integer value for 'nrun'")
)
if(input$selectfile == "saved"){
nmfres <- rda()$nmfres
}else{
if(ncol(merged) >= 300){
# toggleModal(session, "nmfModal1", toggle = "open")
# Run NMF through YABI
# test_yabi <- try(system("yabish login yabi buffalo! backends"))
# Will remove yabi later, right now just don't let yabi command start and run NMF on the shiny01 server
test_yabi <- try(system("naikai"))
if(test_yabi == 0){
createAlert(session, "YabiAlert", "YabiAlert1", title = "Running NMF on Amazon Cloud", style = "success",
content = paste0("Sample size is ", ncol(merged),
". We are running tiis job on the cloud and will notify you when the results are ready.<br>"),
append = FALSE)
# pop up window asking for email
data_path <- "/mnt/sake-uploads"
# data_path <- "~/Desktop/sake-uploads"
output <- file.path(data_path, paste0(file_prefix(), ".txt"))
write.table(merged, output, sep="\t", quote=F)
command <- paste("su - centos -c '/home/centos/yabish_NMF_email.sh",
"-d", output,
"-t", nrow(merged),
"-k", input$num_cluster,
"-r", input$nrun,
"-m", input$mode,
"-n", "FALSE",
"-a", input$algorithm,
"-s", input$nmf_seed,
"-c", 8,
"-x", "FALSE",
"-f", 0,
"-q", "FALSE'")
# send email to the user and stop sake
# t1 <- isolate(try(system(command, intern = TRUE)))
print("Sending command to yabi")
t1 <- try(system(command, wait = FALSE))
print(paste0("Command sent:", command))
if(t1 == 0){
Sys.sleep(5)
closeAlert(session, "YabiAlert1")
Sys.sleep(5)
print("closing YabiAlert and stopping sake")
# stopApp(11)
stop(paste("sent jobs to yabi and stop sake:", command))
}else{
createAlert(session, "YabiAlert", "YabiAlert2", title = "Error message from running Yabi", style = "danger",
content = t1,
append = FALSE)
Sys.sleep(5)
stop(paste("command creates erros:", command))
}
}else{
createAlert(session, "NMFAlert", "NMFAlert1", title = "WARNING", style = "warning",
content = paste0("Sample size is ", ncol(merged), ". This might take long to run. <br>"),
append = FALSE)
}
}
# Run NMF on local server
ptm <- proc.time()
withProgress(message = 'Running NMF', value = 0, {
incProgress(1/2, detail = "Takes around 30~60 seconds for each K")
nmfres <- myNMF(merged,
mode=input$mode,
cluster=input$num_cluster,
nrun=input$nrun,
ncores = cores(),
algorithm = input$algorithm,
seed = nmf_seed()
)
b <- proc.time() - ptm
print("### Time to run NMF ### ")
print(b)
})
closeAlert(session, "NMFAlert1")
closeAlert(session, "YabiAlert2")
}
return(nmfres)
})
### Download NMF ###
observeEvent(input$runNMF, {
output$dlNMF_UI <- renderUI({
tagList(
column(width=2, br(), downloadButton("dlNMF", "Download NMF result", class = 'dwnld'))
)
})
})
output$dlNMF <- downloadHandler(
filename <- function() {
paste( file_prefix(), "nmf_results.zip", sep=".")
},
content = function(file) {
tmpdir <- tempdir()
current_dir <- getwd()
setwd(tmpdir)
filenames <- sapply(c("nmf_Groups", "nmf_Features", "Original_plus_NMF"),
function(x) paste(file_prefix(), x, "csv", sep="."))
write.csv(nmf_groups(), filenames[1], quote=F, row.names=F)
write.csv(nmf_features(), filenames[2], quote=F, row.names=F)
write.csv(ori_plus_nmfResult()[["Total"]], filenames[3], quote=F)
for(i in 1:(length(ori_plus_nmfResult())-1)){
filenames <- c(filenames, paste(file_prefix(), paste0("onlyNMF", i), "csv", sep="."))
write.csv(ori_plus_nmfResult()[[paste("NMF", i)]], filenames[length(filenames)], quote=F)
}
nmfres <- nmf_res()
rawdata <- rawdata()
nmfres_filename <- paste(file_prefix(), "nmfres", "rda", sep=".")
save(rawdata, nmfres, file=nmfres_filename)
zip(zipfile=file, files=c(filenames, nmfres_filename))
setwd(as.character(current_dir))
},
contentType = "application/zip"
)
### Estimate K ###
runSummary <- reactive({
req(class(nmf_res()) == "NMF.rank")
nmf_res <- nmf_res()
summary <- nmf_summary(nmf_res)
if(input$mode=="estim"){
summary <- t(summary)
}
return(summary)
})
output$estimSummary = DT::renderDataTable({
runSummary <- runSummary()
DT::datatable(runSummary, rownames= TRUE,
extensions = 'Buttons',
options = list(dom = 'Bfrtip',
buttons =
list('copy', 'print', list(
extend = 'collection',
buttons = list(list(extend='csv',
filename = filename),
list(extend='excel',
filename = filename),
list(extend='pdf',
filename= filename)),
text = 'Download'
)),
scrollX = TRUE,
pageLength = 12
)
) %>% formatRound(1:ncol(runSummary), 2)
}, server=FALSE)
output$estimPlot <- renderPlot({
validate(
need(class(nmf_res()) == "NMF.rank", "Seems like you haven't run NMF 'estim run' yet\nOr the loaded .rda file is from 'real run' result")
)
nmf_res <- nmf_res()
if(input$estimtype=="consensus"){
consensusmap(nmf_res)
}else if(input$estimtype=="quality"){
plot(nmf_res)
}
}, height=777)
### NMF plots ###
output$nmfplot <- renderPlot({
validate(
need(class(nmf_res()) == "NMFfitX1", "Seems like you haven't run NMF 'real' run yet\nOr the loaded .rda file is from 'estim run' result")
)
nmf_res <- nmf_res()
nmf_plot(nmf_res, type=input$plottype, silorder=input$nmfplot_silhouette,
subsetRow = ifelse(input$select_feature_num>0, as.numeric(input$select_feature_num), TRUE ))
},height = 777)
output$silhouetteplot <- renderPlot({
validate(
need(class(nmf_res()) == "NMFfitX1", "Seems like you haven't run NMF 'real' run yet")
)
nmf_res <- nmf_res()
nmf_silhouette_plot(nmf_res, type=input$plottype)
},height = 777)
output$dl_nmf_estimplot <- downloadHandler(
filename <- function() {
paste(file_prefix(), input$mode, "nmf.pdf", sep=".")
},
content = function(file) {
nmf_res <- nmf_res()
pdf(file, width=input$estim_pdf_w, height=input$estim_pdf_h)
if(input$mode == "estim"){
consensusmap(nmf_res)
print(plot(nmf_res))
}
dev.off()
}
)
output$dl_nmf_realplot <- downloadHandler(
filename <- function() {
paste(file_prefix(), input$mode, "nmf.pdf", sep=".")
},
content = function(file) {
nmf_res <- nmf_res()
pdf(file, width=input$real_pdf_w, height=input$real_pdf_h)
if(input$mode == "real"){
nmf_plot(nmf_res, type="samples", silorder=input$nmfplot_silhouette)
nmf_plot(nmf_res, type="features", silorder=input$nmfplot_silhouette,
subsetRow = ifelse(input$select_feature_num>0, as.numeric(input$select_feature_num), TRUE ))
nmf_plot(nmf_res, type="consensus")
if(!is.null(tsne_2d$data)){
print(nmfplottsne())
}
}
dev.off()
}
)
nmf_groups <- reactive({
validate(
need(class(nmf_res()) == "NMFfitX1", "Seems like you haven't run NMF 'real' run yet")
)
withProgress(message = 'Extracting Groups info', value = NULL, {
nmf_extract_group(nmf_res(), type=input$predict)
})
})
output$nmfGroups <- DT::renderDataTable({
nmf_groups <- nmf_groups()
filename <- "nmf_groups"
colnames(nmf_groups)[which(colnames(nmf_groups)=="nmf_subtypes")] <- "Group"
# Get subset based on selection from scatter plot
event.data <- event_data("plotly_selected", source = "nmfGroups")
if(is.null(event.data) == TRUE){
sel_idx <- NULL
}else{
sel_idx <- subset(event.data, curveNumber == 0)$pointNumber + 1
}
DT::datatable(nmf_groups, rownames=FALSE, escape=-1,
selection = list(selected = sel_idx),
options = list(dom = 'rtip',
scrollX = TRUE,
pageLength = 8,
order=list(list(1,'asc'))
)
) %>% formatRound(3:ncol(nmf_groups), 3)
}, server=FALSE)
# Add t-SNE plot next to NMF group
observeEvent(input$run_nmftSNE, {
plot_data <- plot_data()[['plot_data']]
nsamples <- ncol(plot_data()[['plot_data']])
if(input$nmftsne_perplexity*3 > (nsamples-1)) {
createAlert(session, "perplexityAlert", "perplexityAlert1", title = "WARNING", style = "warning",
content = paste("Perpleixty", input$nmftsne_perplexity,
"is too large compared to num of samples", nsamples),
append = FALSE)
}else{
closeAlert(session, "perplexityAlert1")
withProgress(message = 'Running t-SNE', value = NULL, {
incProgress(1/3, detail = "For t-SNE 2D")
try(
tsne_2d$data <- run_tsne(plot_data, iter=input$tsne_iter, dims=2,
initial_dims = input$tsne_pca_num, theta = input$tsne_theta,
perplexity=input$nmftsne_perplexity, cores=cores())
)
incProgress(2/3, detail = "For t-SNE 3D")
try(
tsne_3d$data <- run_tsne(plot_data, iter=input$tsne_iter, dims=3,
initial_dims = input$tsne_pca_num, theta = input$tsne_theta,
perplexity=input$nmftsne_perplexity, cores=cores())
)
})
}
})
nmfplottsne <- reactive({
plot_data <- plot_data()[['plot_data']]
nsamples <- ncol(plot_data)
validate(
need(input$mode=="real", "This part won't work for 'estim' module\nPlease Try NMF 'Real' run") %then%
need(!is.null(tsne_2d$data), "Please hit 'Run t-SNE' button")
)
color <- NULL
# Will add an option to select by default column names from samples
if(length(nmf_groups()$nmf_subtypes)>0){
idx <- match(colnames(plot_data), nmf_groups()$Sample_ID)
nmf_subtypes <- nmf_groups()$nmf_subtypes[idx]
# color <- create.brewer.color(nmf_subtypes, length(unique(nmf_subtypes)), "naikai")
color <- nmf_subtypes
}else if(length(ColSideColors()[['name']][,1])>0){
color <- ColSideColors()[['name']][,1]
}else{
stop("Error: color scheme undefined")
}
# # default by user selection
point_size <- input$plot_point_size - 5
if(input$nmf_prob){
if(input$predict=="consensus"){
point_size <- point_size + point_size * (1-nmf_groups()$sil_width)
}else if(input$predict == "samples"){
point_size <- point_size + point_size * (1-nmf_groups()$prob)
}
}
nmf_tsne_res <- tsne_2d$data
naikai <- plot_tsne(nmf_tsne_res, color=color, alpha=input$plot_point_alpha,
add.label = input$plot_label, save.plot=F, real.plot=F,
add.legend = input$plot_legend,
point.size=point_size, label.size=input$plot_label_size)
### Check if user select rows (samples)
s = input$nmfGroups_rows_selected
if (length(s)){
if(!input$nmf_prob){
sub_color <- create.brewer.color(nmf_groups()$nmf_subtypes, length(unique(color)), "naikai")
sub_tsne_res <- parse_tsne_res(nmf_tsne_res) %>% as.data.frame %>% dplyr::slice(s)
naikai <- naikai + geom_point(data=sub_tsne_res, aes(x=x, y=y), size=point_size+2, colour=sub_color[s], alpha=input$plot_point_alpha)
}
}
p <- ggplotly(naikai, source="nmfGroups") %>% layout(dragmode = "select")
return(p)
})
output$nmftsneplot <- renderPlotly({
withProgress(message = 'Genrating t-SNE plot', value = NULL, {
nmfplottsne()
})
})
output$nmfgrp_var <- renderPlotly({
nmf_groups <- nmf_groups()
plot_data <- plot_data()[['plot_data']]
idx <- match(colnames(plot_data), nmf_groups$Sample_ID)
colnames(plot_data) <- paste0("NMF", nmf_groups$nmf_subtypes[idx])
dd <- plot_data %>% t %>% reshape2::melt(.) %>% group_by(Var2, Var1) %>% summarise(var=var(value))
num_clus <- dd$Var1 %>% levels %>% length
mycolor <- create.brewer.color(1:num_clus, num=num_clus, name="naikai")
dd$Var1 <- factor(dd$Var1, levels = paste0("NMF", 1:num_clus))
a <- ggplot(dd, aes(var, fill=Var1, colour=Var1)) +
geom_density(alpha=input$cl_alpha) + scale_x_log10() +
scale_colour_manual(values = mycolor) +
scale_fill_manual(values = mycolor) +
labs(x="Log (variance)", y="Density", colour="", fill="") +
theme_bw()
p <- plotly_build(a)
# remove duplicated name
p$data <- lapply(p$data, FUN = function(x){
x$name <- sub("\\((\\S+\\d),.*", "\\1", x$name, perl=TRUE)
x$line$width <- 1.2
return(x)
})
p$layout$legend$font$size <- as.numeric(input$cl_lgsize)
p$layout$xaxis$tickfont$size <- input$cl_tickfont
p$layout$xaxis$titlefont$size <- input$cl_axisfont
p$layout$xaxis$titlefont$family <- "Arial"
p$layout$yaxis$tickfont$size <- input$cl_tickfont
p$layout$yaxis$titlefont$size <- input$cl_axisfont
p$layout$yaxis$titlefont$family <- "Arial"
p
})
output$nmfgrp_expgene <- renderPlotly({
rawdata <- rawdata()
nmf_groups <- nmf_groups()
idx <- match(colnames(rawdata), nmf_groups$Sample_ID)
rawdata <- rawdata[, !is.na(idx)]
colnames(rawdata) <- paste0("NMF", nmf_groups$nmf_subtypes[idx[!is.na(idx)]])
dd <- data.frame(Group=colnames(rawdata), Genes=colSums(rawdata>0))
num_clus <- dd$Group %>% levels %>% length
mycolor <- create.brewer.color(1:num_clus, num=num_clus, name="naikai")
a <- ggplot(dd, aes(Genes, fill=Group, colour=Group)) +
geom_density(alpha=input$cl_alpha) +
scale_fill_manual(values = mycolor, guide="none") +
scale_color_manual(values = mycolor) +
labs(x="No. of expressed genes", y="Density", colour="", fill="") +
theme_bw()
p <- plotly_build(a)
# remove duplicated name
p$data <- lapply(p$data, FUN = function(x){
x$name <- sub("\\((\\S+\\d),.*", "\\1", x$name, perl=TRUE)
x$line$width <- 1.2
return(x)
})
p$layout$legend$font$size <- as.numeric(input$cl_lgsize)
p$layout$xaxis$tickfont$size <- input$cl_tickfont
p$layout$xaxis$titlefont$size <- input$cl_axisfont
p$layout$xaxis$titlefont$family <- "Arial"
p$layout$yaxis$tickfont$size <- input$cl_tickfont
p$layout$yaxis$titlefont$size <- input$cl_axisfont
p$layout$yaxis$titlefont$family <- "Arial"
p
})
output$nmfgrp_coef <- renderPlotly({
nmf_groups <- nmf_groups()
plot_data <- plot_data()[['plot_data']]
idx <- match(colnames(plot_data), nmf_groups$Sample_ID)
nmf_groups <- nmf_groups[idx, ]
colnames(plot_data) <- paste0("NMF", nmf_groups$nmf_subtypes)
pair_cor <- function(x) {
res <- WGCNA::cor(x, nThreads = 4)
return(res[lower.tri(res)])
}
num_clus <- nmf_groups$nmf_subtypes %>% unique %>% length
mycolor <- create.brewer.color(1:num_clus, num=num_clus, name="naikai")
num_poss_cor <- nmf_groups$nmf_subtypes %>% table %>% apply(., 1, function(x) x * (x-1) / 2) %>% sum
group_cor <- data.frame(matrix(NA_real_, nrow=num_poss_cor, ncol=2))
j <- 1
for(i in 1:num_clus){
idx <- nmf_groups$nmf_subtypes == i
res <- pair_cor(plot_data[, idx])
group_cor[j:(j+length(res)-1), ] <- cbind(paste0("NMF", i), res)
j <- j + length(res)
}
colnames(group_cor) <- c("NMF", "Cor")
group_cor$Cor <- as.numeric(group_cor$Cor)
a <- ggplot(data=group_cor, aes(x=NMF, y=Cor, color=NMF)) +
geom_boxplot(aes(color = NMF)) +
stat_summary(fun.y=mean, shape=19, col='black', geom='point') +
theme_bw() +
scale_colour_manual(values = mycolor) +
labs(x="", y="Cor", colour="")
# ggplotly(a)
p <- plotly_build(a)
# remove outliers for plotly boxplot
p$data <- lapply(p$data, FUN = function(x){
if(x$type == "box"){
x$marker = list(opacity = 0)
}
return(x)
})
p$layout$legend$font$size <- as.numeric(input$cl_lgsize)
p$layout$xaxis$tickfont$size <- input$cl_tickfont
p$layout$xaxis$titlefont$size <- input$cl_axisfont
p$layout$xaxis$titlefont$family <- "Arial"
p$layout$yaxis$tickfont$size <- input$cl_tickfont
p$layout$yaxis$titlefont$size <- input$cl_axisfont
p$layout$yaxis$titlefont$family <- "Arial"
p
})
nmf_features <- reactive({
validate(
need(input$mode=="real", "This part won't work for eastimating k\nPlease Try NMF 'Realrun'") %then%
need(class(nmf_res())=="NMFfitX1", "Please rerun NMF with 'real run'")
)
nmf_res <- nmf_res()
# add more options to select features here. then compare with ClaNC
# Add rank by MAD expression
if(input$select_method=="total"){
res <- nmf_extract_feature(nmf_res, method="total") %>% unique
}
else if(input$select_method=="default"){
res <- nmf_extract_feature(nmf_res, method=input$select_method, manual.num = as.numeric(input$select_feature_num)) %>% unique
validate(
need(nrow(res)!=1 & !is.na(res$Gene), "None of the features passed through default NMF filtering criteria, please manually specify number of genes from each group or use 'select features by Rank'") %then%
need(nrow(res)>1, "Not enought features are selected by default! \nMaybe try with more input features or try select features by Rank")
)
}else if(input$select_method=="rank"){
res <- nmf_extract_feature(nmf_res, rawdata = merged(), method=input$select_method, manual.num = as.numeric(input$select_feature_num),
FScutoff=input$select_FScutoff, math=input$select_math) %>% unique
validate(
need(nrow(res)>1, "Not enought features! \nMaybe try lower the FeatureScore cutoff")
)
}
return(res)
})
nmf_features_annot <- reactive({
data <- nmf_features() %>% as.data.table
merged <- NULL
n <- 2
withProgress(message = 'Extracting features info', value = 0, {
incProgress(1/n, detail = "Takes around 5~10 seconds")
ext.func.data <- ext_gene_function(data$Gene, ensembl = ensembl())
setkey(data, "Gene")
setkey(ext.func.data, "hgnc_symbol")
merged <- ext.func.data[data, nomatch=0, allow.cartesian = TRUE]
merged <- as.data.frame(merged) %>%
dplyr::arrange(Group, desc(featureScore))
merged <- merged[, c(4,2,7,6,8)]
merged <- unique(merged)
merged$featureScore <- round(merged$featureScore, 3)
merged$prob <- round(merged$prob, 3)
})
return(merged)
})
output$nmfFeatures <- DT::renderDataTable({
filename <- "nmf_features"
DT::datatable(nmf_features_annot(), rownames=FALSE, escape=-1,
extensions = 'Buttons', selection = 'multiple',
options = list(dom = 'Bfrtip',
buttons =
list('copy', 'print', list(
extend = 'collection',
buttons = list(list(extend='csv',
filename = filename),
list(extend='excel',
filename = filename),
list(extend='pdf',
filename= filename)),
text = 'Download'
)),
pageLength = 50,
autoWidth = TRUE
)
)
}, server = FALSE)
output$nmf_vioplot <- renderPlotly({
req(nmf_features_annot())
validate(
need(!is.null(input$nmfFeatures_rows_selected), "Please select a gene")
)
withProgress(message = 'Generating box/violin plot', value = NULL, {
gene <- nmf_features_annot()[input$nmfFeatures_rows_selected, "GeneCard"] %>%
gsub(".*'>(.*)</a>", "\\1", .)
gene.data <- transform_data$rawdata[gene, ]
gene.data <- gene.data %>%
tibble::rownames_to_column() %>%
dplyr::rename(genename = rowname) %>%
tidyr::gather(Sample, Expr, -genename) %>%
dplyr::mutate(NMF = factor(paste0("NMF", rep(nmf_groups()$nmf_subtypes, each=length(gene)))))
if(input$sel_vioscale == "raw"){
print('Use original scale')
}else if(input$sel_vioscale == "log2"){
gene.data$Expr <- log2(gene.data$Expr+1)
}else if(input$sel_vioscale == "log10"){
gene.data$Expr <- log10(gene.data$Expr+1)
}
# match the coloring from PCA, t-SNE, and heatmap
num_clus <- gene.data$NMF %>% levels %>% length
mycolor <- create.brewer.color(1:num_clus, num=num_clus, name="naikai")
if(input$sel_grp == "Manual"){
col.num <- as.numeric(factor(mycolor, levels = mycolor))
col.lvl <- as.character(mycolor)
col.lvl <- col.lvl[as.numeric(input$sel_grp_order)]
mycolor <- col.lvl[col.num]
}
# add option to select either boxplot or violin plot
if(input$sel_vioplot == "violin"){
a <- ggplot(data=gene.data, aes(x=NMF, y=Expr, color=NMF)) +
geom_violin(aes(color = NMF), scale="width", width=0.6, show.legend = FALSE) +
geom_jitter(aes(color=NMF), alpha=0.6, width=0.1, show.legend = FALSE) +
theme_bw() +
facet_wrap(~ genename, ncol = input$sel_ncol) +
scale_colour_manual(values = mycolor) +
labs(x="", y="", colour="")
}else if(input$sel_vioplot == "box"){
a <- ggplot(data=gene.data, aes(x=NMF, y=Expr, color=NMF)) +
geom_boxplot(aes(color = NMF), show.legend = FALSE) +
geom_jitter(aes(color=NMF), alpha=0.6, width=0.1, show.legend = FALSE) +
theme_bw() +
facet_wrap(~ genename, ncol = input$sel_ncol) +
scale_colour_manual(values = mycolor) +
labs(x="", y="", colour="")
}else{
warning(paste("Unknown plot type:", input$sel_vioplot, "Please check again"))
}
# custom range for Y-axis is only being triggered when ymax > 0
if(!is.na(input$sel_ymax) & input$sel_ymax>0){
validate(
need(input$sel_ymax > input$sel_ymin, "Please make sure 'ymax' is greater than 'ymin'")
)
a <- a + ylim(c(input$sel_ymin, input$sel_ymax))
}
p <- plotly_build(a)
# remove outliers for plotly boxplot
p[[1]]$data <- lapply(p[[1]]$data, FUN = function(x){
if(x$type == "box"){
x$marker = list(outliercolor = 0)
}
return(x)
})
p[[1]]$layout$margin$l <- p[[1]]$layout$margin$l + 5
p[[1]]$layout$legend$font$size <- input$sel_legend_size
p[[1]]$layout$showlegend = input$sel_show_legend
# loop through ncol(Xaxis) and gene (Yaxis)
num_xaxis <- grep("xaxis", names(p[[1]]$layout)) %>% length
num_yaxis <- grep("yaxis", names(p[[1]]$layout)) %>% length
for(i in 1:num_xaxis){
p[[1]]$layout[[paste0("xaxis", ifelse(i>1, i, ""))]]$tickfont$size <- input$sel_xfont_size
}
for(i in 1:num_yaxis){
p[[1]]$layout[[paste0("yaxis", ifelse(i>1, i, ""))]]$tickfont$size <- input$sel_yfont_size
}
for(i in 1:length(gene)){
p[[1]]$layout$annotations[[i]]$font$size <- input$sel_title_size
}
p
})
})
ori_plus_nmfResult <-reactive({
rawdata <- rawdata()
validate(
need(!is.null(nmf_groups()), "Please run NMF and get the result first")
)
res <- list()
idx <- match(colnames(rawdata), nmf_groups()$Sample_ID)
newdata <- rawdata
colnames(newdata) <- paste(paste0("NMF",nmf_groups()$nmf_subtypes[idx]), colnames(newdata), sep="_")
res[["Total"]] <- newdata
# add separate files for each NMF group
for(i in 1:length(unique(nmf_groups()$nmf_subtypes))){
subidx <- nmf_groups()$nmf_subtypes == i
res[[paste("NMF", i)]] <- newdata[, idx[subidx]]
}
return(res)
})
### Statistics ###
output$runStatSummary <- renderTable({
runSummary()
}, sanitize.text.function = function(x) x, options = list(orderClasses = TRUE, lengthMenu = c(10,25,50,100, NROW(merged)), pageLength=25)
)
output$runStatPlot <- renderPlot({
nmf_res <- nmf_res()
plot(nmf_res)
}, height=666)
ensembl <- reactive({
# Need to add a testing here when the BioMart is down
# Or find offline annotation db
# ensembl = useMart(biomart="ENSEMBL_MART_ENSEMBL", host="grch37.ensembl.org", path="/biomart/martservice" ,dataset="hsapiens_gene_ensembl")
ensembl = useMart(biomart="ENSEMBL_MART_ENSEMBL", host="jul2016.archive.ensembl.org", path="/biomart/martservice" ,dataset="hsapiens_gene_ensembl")
})
#' Heatmap setting
geneListInfo <- reactive({
rawdata <- transform_data$data
gene.list <- NULL
gene.groups <- NULL
if(input$heatmapGene=='nmf'){
nmf_features <- nmf_features()
validate(
need(nrow(nmf_features)>0, "NMF not run yet? \nPlease run NMF classification and then try again!")
)
gene.list <- nmf_features$Gene
}
else if(input$heatmapGene=='rank'){
select.genes <- select_top_n(apply(rawdata,1,input$heat.math), n=as.numeric(input$heat.top.num), bottom=F)
gene.list <- names(select.genes)
}
else if(input$heatmapGene=='preload'){
validate(
need(!is.null(input$predefined_list), "Please select at least one predefined gene list")
)
for(i in 1:length(input$predefined_list)){
gene.list <- c(gene.list, pre_genelist[[ input$predefined_list[i] ]]$Gene)
# if file contains 'Group' info in one of the columns, use that for RowColors
if("Group" %in% colnames(pre_genelist[[ input$predefined_list[i] ]])){
gene.groups <- c(gene.groups, pre_genelist[[ input$predefined_list[i] ]]$Group)
}
}
}
else if(input$heatmapGene=='manual'){
validate(
need(!is.null(input$manual_genes), "You didn't select any genes. Please select genes from the data") %then%
need(length(input$manual_genes)>=2, "Please select at least two genes! ")
)
gene.list <- input$manual_genes
}
else if(input$heatmapGene=='upload'){
validate(
need(!is.null(input$heatmapGeneFile), "Please upload a file with list of genes, with column name 'Gene' ")
)
ext <- file_ext(input$heatmapGeneFile)[1]
if (ext=="csv"){
genefile <- read.csv(input$heatmapGeneFile$datapath, header=T, stringsAsFactors = FALSE)
}else if (ext=="out" || ext=="tsv" || ext=="txt"){
genefile <- read.table(input$heatmapGeneFile$datapath, header=T, sep="\t", stringsAsFactors = FALSE)
}else{
print("File format doesn't support, please use '.csv', '.txt', '.out', '.tsv' and try again")
return(NULL)
}
validate(
need(sum(colnames(genefile)=='Gene'), "Please make sure genefile contains column name 'Gene'")
)
gene.list <- genefile$Gene
if("Group" %in% colnames(genefile)){
gene.groups <- genefile$Group
}
}
validate(
need(length(intersect(gene.list, rownames(rawdata)))>=2, "We did not find enough genes that overlap with the expression data (<2). \nPlease check again and provide a new list! \nMaybe this is from species other than Human (default gene list are from Human)")
)
# only select unique gene name on rawdata
idx <- unique(match(gene.list, rownames(rawdata)))
na.genes <- gene.list[is.na(idx)]
if(length(na.genes)>0){
print(c("These genes cannot be found in the data:", paste(na.genes, collapse=", ")))
}
data.rownames <- rownames(rawdata)[idx[!is.na(idx)]]
# match back those gene names to genelist to get idx for gene.groups
ii <- match(data.rownames, gene.list)
if(!is.null(gene.groups)){
gene.groups <- gene.groups[ii]
}
gene.list <- data.rownames
# finish loading or selecting gene list file
res <- list()
res[['list']] <- gene.list
res[['group']] <- gene.groups
return(res)
})
OrdColopts <- reactive({
c("Default", "Hierarchical", "Filename", "GeneExpr")
})
ColClrByopts <- reactive({
c("Filename")
})
observeEvent(rawdata(), {
updateSelectizeInput(session, 'OrdCol',
server = TRUE,
choices = as.character(OrdColopts()),
selected = "Hierarchical"
)
updateSelectizeInput(session, 'ColClrBy',
server = TRUE,
choices = as.character(ColClrByopts()),
selected = "Filename"
)
})
observeEvent(input$runNMF, {
updateSelectizeInput(session, 'OrdCol',
server = TRUE,
choices = as.character(c(OrdColopts(), "NMF Group")),
selected = "NMF Group"
)
updateSelectizeInput(session, 'ColClrBy',
server = TRUE,
choices = as.character(c(ColClrByopts(), "NMF Group")),
selected = "NMF Group"
)
})
# separate out this for PCA and t-SNE plots
plot_data <- reactive({
rawdata <- transform_data$data
genelist <- geneListInfo()[['list']]
idx <- match(genelist, rownames(rawdata))
plot_data <- rawdata[idx,]
res <- list()
res[["plot_data"]] <- plot_data
res[["genelist"]] <- genelist
return(res)
})
heatmap_data <- reactive({
rawdata <- transform_data$data
genelist <- geneListInfo()[['list']]
gene.groups <- geneListInfo()[['group']]
# only select unique gene name on rawdata
idx <- match(genelist, rownames(rawdata))
heatmap_data <- rawdata[idx,]
if(input$OrdCol == 'Filename'){ # Add options to sort by column groups
column.names <- strsplit(colnames(heatmap_data), split="\\_")
validate(
need(input$sortcolumn_num <= length(column.names[[1]]), "Group num is too big, select a smaller number!")
)
idx <- order(sapply(column.names, function(x) x[[input$sortcolumn_num]]))
}else if(input$OrdCol == 'GeneExpr'){ # Add options to sort by expression value of specific gene
validate(
need(input$sortgene_by!="", "Please select a gene!")
)
idx <- order(subset(heatmap_data, rownames(heatmap_data)==input$sortgene_by))
}
else if(input$OrdCol == 'NMF Group') { # Add options to sort by nmf groups
idx <- nmf_groups()$nmf_subtypes %>% order
}else{
idx <- colnames(heatmap_data)
}
heatmap_data <- heatmap_data[, idx]
res <- list()
res[["heatmap_data"]] <- heatmap_data
res[["genelist"]] <- genelist
res[["gene.groups"]] <- gene.groups
return(res)
})
color <- reactive({
color <- switch(input$heatColorScheme,
"RdBu_Brewer" = colorRampPalette(rev(brewer.pal(input$color_num, "RdBu")))(51),
"RdYlBu_Brewer" = colorRampPalette(rev(brewer.pal(input$color_num, "RdYlBu")))(51),
"Red_Green" = colorRampPalette(c("red2", "white", "green"))(input$color_num),
"Red_Blue" = colorRampPalette(c("blue", "white", "red"))(input$color_num),
"Rainbow" = colorRampPalette(c("lightblue", "blue", "white", "gold2", "red2"))(input$color_num),
"Black_White" = colorRampPalette(c("black", "white"))(input$color_num)
)
return(color)
})
ColScheme <- reactive({
ColScheme <- list()
ColScheme <- c(input$ColScheme1, input$ColScheme2, input$ColScheme3, input$ColScheme4, input$ColScheme5, input$ColScheme6)
return(ColScheme)
})
ColSideColors <- reactive({
heatmap_data <- heatmap_data()[['heatmap_data']]
res <- list()
if(input$ColClrBy=='NMF Group'){
req(nmf_groups)
# rematch the sample order to the sorted order based on selection from input$OrdCol
idx <- match(colnames(heatmap_data), nmf_groups()$Sample_ID)
nmf_subtypes <- nmf_groups()$nmf_subtypes[idx]
res[['color']] <- create.brewer.color(nmf_subtypes, length(unique(nmf_subtypes)), "naikai") %>% as.matrix
res[['name']] <- paste0("NMF", nmf_subtypes) %>% as.matrix
}else if(input$ColClrBy == 'Filename'){
res <- name_to_color(colnames(heatmap_data), split_pattern ="\\_",
num_color = 6,
ColScheme = ColScheme()[1:6] )
# num_color = input$ColSideColorsNum,
# ColScheme = ColScheme()[1:input$ColSideColorsNum] )
# add options to select which column to show
if(is.null(input$SelColSideColors)){
res[['color']] <- NULL
res[['name']] <- NULL
}else{
res[['color']] <- res[['color']][, as.numeric(input$SelColSideColors)] %>% as.matrix()
res[['name']] <- res[['name']][, as.numeric(input$SelColSideColors)] %>% as.matrix()
}
}
return(res)
})
RowSideColors <- reactive({
groups <- heatmap_data()[['gene.groups']]
row.color <- NULL
res <- list()
res[['color']] <- NULL
res[['name']] <- NULL
if(!is.null(groups)){
row.color <- rbind(row.color, create.brewer.color(groups, input$RowColScheme.num, input$RowColScheme))
res[["color"]] <- as.matrix(row.color)
res[["name"]] <- as.matrix(groups)
}
return(res)
})
Colv <- reactive({
Colv <- NULL
validate(
need(nrow(heatmap_data()[['heatmap_data']] > 2), "Did not get enough genes that overlapped with your data\nPlease try another gene list or use rank from data set")
)
if(input$OrdCol=="Hierarchical"){
Colv <- heatmap_data()[['heatmap_data']] %>% t %>% dist(method=input$distance) %>% hclust(method=input$linkage) %>% as.dendrogram
ColMean <- colMeans(heatmap_data()[['heatmap_data']], na.rm=T)
Colv <- reorderfun(Colv, ColMean)
}
return(Colv)
})
Rowv <- reactive({
Rowv <- NULL
validate(
need(nrow(heatmap_data()[['heatmap_data']] > 2), "Did not get enough genes that overlapped with your data\nPlease try another gene list or use rank from data set")
)
if(input$OrdRow=="hier"){
Rowv <- heatmap_data()[['heatmap_data']] %>% dist(method=input$distance) %>% hclust(method=input$linkage) %>% as.dendrogram
RowMean <- rowMeans(heatmap_data()[['heatmap_data']], na.rm=T)
Rowv <- reorderfun(Rowv, RowMean)
}
return(Rowv)
})
dendrogram <- reactive({
if(is.null(Rowv()) & is.null(Colv())){
return("none")
}else if(is.null(Rowv()) & !is.null(Colv())){
return("column")
}else if(!is.null(Rowv()) & is.null(Colv())){
return("row")
}else if(!is.null(Rowv()) & !is.null(Colv())){
return("both")
}
})
### Observe and update input selection ###
observe({
transform_data <- transform_data$data
if(!is.null(transform_data)){
if (length(rownames(transform_data)) > 1){
# update the render function for selectize - update manual if user choose to select its own markers
updateSelectizeInput(session, 'manual_genes',
server = TRUE,
choices = as.character(rownames(transform_data))
)
}
}
})
observe({
if (!is.null(geneListInfo()[['list']])){
updateSelectizeInput(session, 'sortgene_by',
server = TRUE,
choices = as.character(geneListInfo()[['list']])
)
}
})
output$d3heatmap <- renderD3heatmap({
heatmap_data <- heatmap_data()[['heatmap_data']]
ColSideColors <- ColSideColors()[["color"]]
if(!is.null(ColSideColors)){
ColSideColors <- t(ColSideColors[, ncol(ColSideColors):1])
}
myd3Heatmap(
as.matrix(heatmap_data),
type="heatmap",
dist.m=input$distance,
hclust.m=input$linkage,
scale="row",
scale.method=input$scale_method,
scale.first=input$scale_first,
RowSideColors = RowSideColors()[["color"]],
ColSideColors = ColSideColors,
color = color(),
Colv = Colv(),
Rowv = Rowv(),
cexCol=input$cexRow,
cexRow=input$cexCol,
na.rm = T,
show_grid = FALSE,
dendrogram = dendrogram(),
dendro.ord = "manual",
anim_duration = 0
)
})
plotInput <- function(){
heatmap_data <- heatmap_data()[['heatmap_data']]
filename <- paste(input$type, input$math, paste0(input$orders, paste0(nrow(merged()), "genes")), sep=".")
myHeatmap.3(as.matrix(heatmap_data),
type="heatmap",
color=color(),
RowSideColors = RowSideColors()[["color"]],
RowSideColors.name=RowSideColors()[["name"]],
ColSideColors = ColSideColors()[["color"]],
ColSideColorsSize=input$ColSideColorsSize,
ColSideColors.name=ColSideColors()[["name"]],
cor.method=input$cor.method,
dist.m=input$distance,
hclust.m=input$linkage,
title=input$title,
Colv = Colv(),
Rowv = Rowv(),
scale="row",
scale.method=input$scale_method,
scale.first=input$scale_first,
dendro.ord="manual",
dendrogram=dendrogram(),
cexCol=input$cexCol,
cexRow=input$cexRow,
col.legend = input$col_legend,
row.legend = input$row_legend,
file.prefix=filename,
save.image=F
)
}
output$dl_heatmap<- downloadHandler(
filename <- function() {
paste(file_prefix(), "heatmap.pdf", sep="_")
},
content = function(file) {
pdf(file, width = input$heat_pdf_wd, height = input$heat_pdf_ht)
plotInput()
dev.off()
}
)
plot_colopts <- reactive({
c("Default", "Filename", "GeneExpr", "ReadDepth", "NumExpresGenes")
})
observeEvent(rawdata(), {
updateSelectizeInput(session, 'pt_col',
server = TRUE,
choices = as.character(plot_colopts()),
selected = "Filename"
)
choices <- rawdata() %>% colnames() %>% strsplit("_") %>% .[[1]] %>% length %>% seq(1, .) %>% as.character()
updateSelectizeInput(session, 'pt_file_grp',
server = TRUE,
choices = choices,
selected = choices[1]
)
})
observeEvent(input$pt_file_grp, {
# this should change with respect to the selection of filename group
choices <- rawdata() %>% colnames() %>% strsplit("_") %>%
sapply(., function(x) x[as.numeric(input$pt_file_grp)]) %>% unique %>% length %>%
seq(1, .) %>% as.character()
updateSelectizeInput(session, 'pt_grp_order',
server = TRUE,
choices = choices,
selected = choices[1]
)
})
observeEvent(transform_data$data, {
updateSelectizeInput(session, 'pt_allgene',
server = TRUE,
choices = as.character(rownames(transform_data$data)),
selected = as.character(rownames(transform_data$data)[1])
)
})
observeEvent(input$runNMF,{
updateSelectizeInput(session, 'pt_col',
server = TRUE,
choices = as.character(c(plot_colopts(), "NMF Group", "NMF Feature")),
selected = 'NMF Group'
)
updateSelectizeInput(session, 'pt_nmfgene',
server = TRUE,
choices = as.character(nmf_features()$Gene),
selected = as.character(nmf_features()$Gene[1])
)
})
observeEvent(nmf_groups()$nmf_subtypes, {
choices <- nmf_groups()$nmf_subtypes %>% unique %>% length %>% seq(1, .) %>% as.character()
updateSelectizeInput(session, 'sel_grp_order',
server = TRUE,
choices = choices,
selected = choices[1]
)
})
point_col <- reactive({
# col <- toRGB("steelblue")
col <- "steelblue"
plot_data <- plot_data()[['plot_data']]
group <- NULL
if(input$pt_col == "Default"){
group <- "Default"
}else if(input$pt_col == "NMF Group"){
req(nmf_groups)
idx <- match(colnames(plot_data), nmf_groups()$Sample_ID)
nmf_subtypes <- nmf_groups()$nmf_subtypes[idx]
col <- create.brewer.color(nmf_subtypes, length(unique(nmf_subtypes)), "naikai")
group <- paste0("NMF", nmf_subtypes)
}else if(input$pt_col == "NMF Feature"){
req(input$pt_nmfgene)
idx <- match(colnames(plot_data), colnames(transform_data$data))
gene_data <- transform_data$data[input$pt_nmfgene, idx] %>% as.numeric
col <- create.brewer.color(gene_data, num = 9, name="YlOrRd")
group <- input$pt_nmfgene
}else if(input$pt_col == "Filename"){
filename <- colnames(plot_data)
filename <- sapply(strsplit(filename, "_"), function(x) paste0(x[as.numeric(input$pt_file_grp)], collapse = "_"))
col <- create.brewer.color(filename, length(unique(filename)), "naikai") %>% as.matrix
group <- filename %>% as.matrix
# add option to manual select the order of coloring
if(input$pt_sel_grp_order == "Manual"){
col <- col[, 1]
col.num <- as.numeric(factor(col, levels = unique(col) %>% as.character))
col.lvl <- unique(col)
col.lvl <- col.lvl[as.numeric(input$pt_grp_order)]
col <- col.lvl[col.num] %>% as.matrix
}
}else if(input$pt_col == "GeneExpr"){
req(input$pt_allgene)
idx <- match(colnames(plot_data), colnames(transform_data$data))
gene_data <- transform_data$data[input$pt_allgene, idx] %>% as.numeric
col <- create.brewer.color(gene_data, num = 9, name="YlOrRd")
group <- input$pt_allgene
}else if(input$pt_col == "ReadDepth"){
req(rawdata())
idx <- match(colnames(plot_data), colnames(rawdata()))
gene_data <- rawdata()[, idx] %>% colSums() %>% as.numeric
col <- create.brewer.color(gene_data, num = 9, name="YlOrRd")
group <- "Read Depth"
}else if(input$pt_col == "NumExpresGenes"){
req(rawdata())
idx <- match(colnames(plot_data), colnames(rawdata()))
gene_data <- colSums(rawdata()[, idx]>0) %>% as.numeric
col <- create.brewer.color(gene_data, num = 9, name="YlOrRd")
group <- "Num Expressed Genes"
}else{
warning("Wrong point color assigning method!")
}
res <- list()
res[['color']] <- col
res[['group']] <- group
return(res)
})
pca_data <- reactive({
data <- plot_data()[['plot_data']]
withProgress(message = 'Calculating PCs', value = NULL, {
pc <- prcomp(t(data), center = TRUE)
vars<- pc$sdev^2
vars<- signif(vars/sum(vars) * 100, 2)
title = paste(paste0("PC", 1:length(vars)), paste0("(", vars, "% Var)"))
projection <- as.data.frame(pc$x)
})
res <- list()
res[['proj']] <- projection
res[['title']] <- title
return(res)
})
#' PCA plot
output$pcaPlot_2D <- renderPlotly({
projection <- pca_data()[['proj']]
title <- pca_data()[['title']]
validate(
need(as.numeric(input$pca_x) != as.numeric(input$pca_y), "PC on the X-axis is the same as PC on the Y-axis.\nPlease select different number for these two axes")
)
idx <- c(as.numeric(input$pca_x), as.numeric(input$pca_y))
withProgress(message = 'Generating 2d PCA plot', value = NULL, {
projection <- projection %>%
dplyr::select(idx) %>%
set_colnames(c("pca_x", "pca_y")) %>%
mutate(color = point_col()[['color']] %>% as.character(),
group = point_col()[['group']] %>% as.character(),
Sample = colnames(plot_data()[['plot_data']])) %>%
arrange(group)
if(input$plot_label){
t <- list( size=input$plot_label_size, color=toRGB("grey50") )
p <- plot_ly(projection, x = ~pca_x, y = ~pca_y, type="scatter", mode="markers+text",
color = ~group, colors = ~unique(color),
text = ~Sample, hoverinfo="text", textfont=t, textposition="top middle",
marker = list(color = ~color, size=input$plot_point_size+3, opacity=input$plot_point_alpha))
}else{
p <- plot_ly(projection, x = ~pca_x, y = ~pca_y, type="scatter", mode="markers",
color = ~group, colors = ~unique(color),
text = ~Sample, hoverinfo="text",
marker = list(color = ~color, size=input$plot_point_size+3, opacity=input$plot_point_alpha)) #%>%
# toWebGL()
}
p %>% layout(showlegend = input$plot_legend,
xaxis = list(title = title[idx[1]], zeroline=FALSE),
yaxis = list(title = title[idx[2]], zeroline=FALSE))
})
})
output$pcaPlot_3D <- renderPlotly({
projection <- pca_data()[['proj']]
title <- pca_data()[['title']]
idx <- c(1,2,3)
withProgress(message = 'Generating 3d PCA plot', value = NULL, {
projection <- projection %>%
dplyr::select(idx) %>%
set_colnames(c("pca_x", "pca_y", "pca_z")) %>%
mutate(color = point_col()[['color']] %>% as.character(),
group = point_col()[['group']] %>% as.character(),
Sample = colnames(plot_data()[['plot_data']])) %>%
arrange(group)
if(input$plot_label){
p <- plot_ly(data=projection, x = ~pca_x, y = ~pca_y, z = ~pca_z, type="scatter3d", mode="markers+text",
color = ~group, colors = ~unique(color), text = ~Sample, hoverinfo="text",
marker = list(color = ~color, size=input$plot_point_size, opacity=input$plot_point_alpha))
}else{
p <- plot_ly(data=projection, x = ~pca_x, y = ~pca_y, z = ~pca_z, type="scatter3d", mode="markers",
color = ~group, colors = ~unique(color), text = ~Sample, hoverinfo="text",
marker = list(color = ~color, size=input$plot_point_size, opacity=input$plot_point_alpha)) #%>%
# toWebGL()
}
p %>% layout(showlegend = input$plot_legend,
scene = list(
xaxis = list(title = title[1]),
yaxis = list(title = title[2]),
zaxis = list(title = title[3]))
)
})
})
#' t-SNE plot
cores <- reactive({
avail_cores <- parallel::detectCores()
if(avail_cores==1){
return(1)
}else if(as.numeric(input$ncores) <= avail_cores){
return(as.numeric(input$ncores))
}else{
return(avail_cores)
}
})
tsne_2d <- reactiveValues(data=NULL)
tsne_3d <- reactiveValues(data=NULL)
observeEvent(input$runtSNE, {
nsamples <- ncol(plot_data()[['plot_data']])
if(input$tsne_perplexity*3 > (nsamples-1)) {
createAlert(session, "visualperplexityAlert", "visualperplexityAlert1", title = "WARNING", style = "warning",
content = paste("Perpleixty", input$tsne_perplexity,
"is too large compared to num of samples", nsamples),
append = FALSE)
}else {
closeAlert(session, "visualperplexityAlert1")
plot_data <- plot_data()[['plot_data']]
withProgress(message = 'Running t-SNE', value = NULL, {
incProgress(1/3, detail = "For t-SNE 2D")
tsne_2d$data <- run_tsne(plot_data, iter=input$tsne_iter, dims=2,
perplexity=input$tsne_perplexity, cores=cores())
incProgress(2/3, detail = "For t-SNE 3D")
tsne_3d$data <- run_tsne(plot_data, iter=input$tsne_iter, dims=3,
perplexity=input$tsne_perplexity, cores=cores())
})
}
})
plot_tsne_2d <- reactive({
nsamples <- ncol(plot_data()[['plot_data']])
validate(
need(!is.null(tsne_2d$data), "Please click 'Run t-SNE' button")
)
tsne_out <- isolate(tsne_2d$data)
withProgress(message = 'Genrating t-SNE plot', value = 0, {
incProgress(2/3, detail = "Usually takes 15~20 seconds")
color <- "steelblue"
projection <- parse_tsne_res(tsne_out) %>%
data.frame %>%
mutate(group = point_col()[['group']] %>% as.character(),
color = point_col()[['color']] %>% as.character(),
Sample = rownames(.)) %>%
arrange(group)
min.cost <- signif(tsne_out$itercosts[length(tsne_out$itercosts)], 2)
if(input$plot_label){
t <- list( size=input$plot_label_size, color=toRGB("grey50") )
p <- plot_ly(projection, x = ~x, y = ~y, type="scatter", mode="markers+text",
color = ~group, colors = ~unique(color),
text= ~Sample, hoverinfo="text", textposition="top middle", textfont=t,
marker = list(color = ~color, size=input$plot_point_size+3, opacity=input$plot_point_alpha))
}else{
p <- plot_ly(projection, x = ~x, y = ~y, type="scatter", mode="markers", text= ~Sample, hoverinfo="text",
color = ~group, colors = ~unique(color),
marker = list(color = ~color, size=input$plot_point_size+3, opacity=input$plot_point_alpha)) #%>%
# toWebGL()
}
p %>% layout(showlegend = input$plot_legend,
title = title,
xaxis = list(title = "Component 1", zeroline=FALSE),
yaxis = list(title = "Component 2", zeroline=FALSE))
})
})
output$tsneplot_2d <- renderPlotly({
plot_tsne_2d()
})
output$tsneplot_3d <- renderPlotly({
nsamples <- ncol(plot_data()[['plot_data']])
if(is.null(tsne_3d$data)) return ()
projection <- isolate(as.data.frame(tsne_3d$data$Y)) %>%
set_colnames(c("tsne_1", "tsne_2", "tsne_3")) %>%
mutate(group = point_col()[['group']] %>% as.character(),
color = point_col()[['color']] %>% as.character(),
Sample = rownames(.)) %>%
arrange(group)
if(input$plot_label){
p <- plot_ly(data=projection, x = ~tsne_1, y = ~tsne_2, z = ~tsne_3, type="scatter3d", mode="markers+text",
color = ~group, colors = ~unique(color), text= ~Sample, hoverinfo="text",
marker = list(color = ~color, size=input$plot_point_size, opacity=input$plot_point_alpha))
}else{
p <- plot_ly(data=projection, x = ~tsne_1, y = ~tsne_2, z = ~tsne_3, type="scatter3d", mode="markers",
color = ~group, colors = ~unique(color), text= ~Sample, hoverinfo="text",
marker = list(color = ~color, size=input$plot_point_size, opacity=input$plot_point_alpha)) %>%
toWebGL()
}
p %>% layout(showlegend = input$plot_legend,
scene = list(
xaxis = list(title = "Component 1"),
yaxis = list(title = "Component 2"),
zaxis = list(title = "Component 3"))
)
})
### DESeq2 ###
observeEvent(input$runNMF, {
validate(
need(is.numeric(input$num_cluster), "num of cluster is not a valid numeric number")
)
updateSelectizeInput(session, 'de_group1',
server = TRUE,
choices = as.character(paste0("NMF", sort(unique(nmf_groups()$nmf_subtypes)))),
selected = "NMF1"
)
updateSelectizeInput(session, 'de_group2',
server = TRUE,
choices = as.character(paste0("NMF", sort(unique(nmf_groups()$nmf_subtypes)))),
selected = "NMF2"
)
})
deseq_res <- eventReactive(input$runDESeq, {
if(input$selectfile == "saved"){
if(!is.null(rda()$dds)){
return(rda()$dds)
}
}
rawdata <- rawdata()
if(input$de_conv2int)
rawdata <- round(rawdata)
if(!all(sapply(rawdata, is.whole))){
createAlert(session, "alert", "exampleAlert", title = "WARNING", style = "warning",
content = "Data contains non-integer value, Please use raw count data for DESeq2. <br>
Or Use the option above to 'Force convert to integer' <b>(Results are not statistically rigorous)</b>",
append = FALSE)
return(NULL)
}else{
closeAlert(session, "exampleAlert")
}
# since the user might have filtered samples during the QC steps,
# rematch sample_ID and only kept the ones that are presented in NMF res
idx <- match(nmf_groups()$Sample_ID, colnames(rawdata))
rawdata <- rawdata[, idx]
register(BiocParallel::MulticoreParam(workers = cores()))
colData <- data.frame(Group = paste0("NMF", nmf_groups()$nmf_subtypes))
ddsfeatureCounts <- DESeq2::DESeqDataSetFromMatrix(countData = rawdata,
colData = colData,
design = ~ Group)
### Set betaPrior=FALSE to go with MLE LFC to get simple LFC = (avg in group2/ avg in group1)
# dds <- DESeq(ddsfeatureCounts, parallel=T, betaPrior=FALSE)
withProgress(message = 'Running DESeq2', value = NULL, {
dds <- DESeq2::DESeq(ddsfeatureCounts, parallel=T)
})
return(dds)
})
filt_deseq_res <- reactive({
req(deseq_res())
withProgress(message = 'Summarizing DESeq2 result', value = NULL, {
register(BiocParallel::MulticoreParam(workers = 8))
DESeq2::results(deseq_res(), parallel = TRUE,
contrast = c("Group", input$de_group1, input$de_group2)) %>%
as.data.frame %>%
subset(padj <= input$de_alpha) %>%
.[order(.$padj), ]
})
})
output$deseq_table <- DT::renderDataTable({
validate(
need(nrow(filt_deseq_res()) > 0, "No significant DE results found, Maybe try lower your alpha threshold?")
)
filename <- fileinfo()[['name']]
GeneCard <- paste0("<a href='http://www.genecards.org/cgi-bin/carddisp.pl?gene=",
rownames(filt_deseq_res()), "'>", rownames(filt_deseq_res()), "</a>")
filt_res <- data.frame(Gene=GeneCard,
filt_deseq_res())
DT::datatable(filt_res, selection='single', rownames=FALSE,
extensions = 'Buttons',
escape = FALSE,
options = list(dom = 'Bfrtip',
buttons =
list('colvis', 'copy', 'print', list(
extend = 'collection',
buttons = list(list(extend='csv',
filename = filename),
list(extend='excel',
filename = filename)),
text = 'Download'
)),
scrollY = TRUE,
pageLength = 12,
autoWidth = TRUE
)
) %>% formatRound(2:ncol(filt_res), 3)
}, server = FALSE)
output$deseq_boxplot <- renderPlotly({
req(deseq_res())
validate(
need(!is.null(input$deseq_table_rows_selected), "Please select a gene")
)
withProgress(message = 'Generating boxplot', value = NULL, {
norm.data <- DESeq2::counts(deseq_res(), normalized=TRUE)
gene <- rownames(filt_deseq_res())[input$deseq_table_rows_selected]
gene.data <- norm.data[gene, ]
gene.data <- data.frame(Sample = names(gene.data),
Expr = gene.data,
NMF = factor(paste0("NMF", nmf_groups()$nmf_subtypes)))
# match the coloring from PCA, t-SNE, and heatmap
num_clus <- gene.data$NMF %>% levels %>% length
mycolor <- create.brewer.color(1:num_clus, num=num_clus, name="naikai")
a <- ggplot(data=gene.data, aes(x=NMF, y=Expr, color=NMF)) +
geom_boxplot(aes(color = NMF)) +
geom_jitter(aes(color=NMF, text=Sample), alpha=0.6, width=0.1) +
theme_bw() +
scale_colour_manual(values = mycolor) +
labs(title=gene, x="", y="Expression", colour="")
p <- plotly_build(a)
# remove outliers for plotly boxplot
p$data <- lapply(p$data, FUN = function(x){
if(x$type == "box"){
x$marker = list(opacity = 0)
}
return(x)
})
p$layout$margin$l <- p$layout$margin$l + 15
p$layout$annotations[[1]]$x <- -0.1
p$layout$xaxis$tickfont$size <- 8
p$layout$annotations[[1]]$font$size <- 11
p$layout$yaxis$tickfont$size <- 10
p
})
})
### Pathway analysis ###
nmf_group_feature_rank <- reactive({
validate(
need(is.numeric(input$num_cluster), "num of cluster is not a numeric value") %then%
need(!is.null(nmf_groups()), "Please run NMF first")
)
rawdata <- transform_data$data
# rawdata <- rawdata[rowMeans(rawdata) >= input$min_rowMean, ]
nmf_group <- nmf_groups()$nmf_subtypes
num_nmf_subtypes <- nmf_group %>% unique %>% length
if(input$rank_method=="featureScore"){
print('Do we want to use featureScore, there are only few genes')
}else if(input$rank_method=="logFC"){
nmf_feature_rank <- data.frame(matrix(NA_real_, nrow=nrow(rawdata), ncol=num_nmf_subtypes))
for(i in 1:num_nmf_subtypes){
idx <- nmf_group==i
nmf_feature_rank[, i] <- rawdata %>% as.data.frame %>%
dplyr::mutate(logFC = log(rowMeans(.[idx])/rowMeans(.[!idx]))) %>%
dplyr::select(logFC)
}
colnames(nmf_feature_rank) <- paste0("NMF", 1:num_nmf_subtypes)
rownames(nmf_feature_rank) <- rownames(rawdata)
}
return(nmf_feature_rank)
})
geneset.file <- reactive({
validate(
need(!is.null(input$predefined_list), "Please select at least one predefined gene sets")
)
file <- pre_geneset[[ input$predefined_list ]]
return (file)
})
ncpus <- reactive({
cpus <- c(16,8,6,4,2,1)
cpus[min(which(sapply(cpus, function(x) input$nPerm%%x)==0))]
})
observeEvent(input$runNMF, {
updateSelectizeInput(session, 'pathway_group',
server = TRUE,
choices = as.character(paste0("NMF", sort(unique(nmf_groups()$nmf_subtypes))))
)
})
pathview.species <- reactive({
species <- input$species
if(species == "Human" | species == "human" | species == "hg19" | species == "hsa"){
pathview.species <- "hsa"
}else if(species == "Mouse" | species == "mouse" | species == "mm9" | species == "mmu"){
pathview.species <- "mmu"
}else if(species == "Drosophila" | species == "dm3" | species == "dm6"){
pathview.species <- "dm3"
}
})
id.org <- reactive({
species <- input$species
if(species == "Human" | species == "human" | species == "hg19" | species == "hsa"){
id.org <- "Hs"
}else if(species == "Mouse" | species == "mouse" | species == "mm9" | species == "mmu"){
id.org <- "Mm"
}
})
kegg.gs <- reactive({
species <- input$species
withProgress(message = 'Extracting KEGG pathway information', value = NULL, {
if(species == "Human" | species == "human" | species == "hg19" | species == "hsa"){
kg.human<- kegg.gsets("human")
kegg.gs<- kg.human$kg.sets[kg.human$sigmet.idx]
}else if(species == "Mouse" | species == "mouse" | species == "mm9" | species == "mmu"){
kg.mouse<- kegg.gsets("mouse")
kegg.gs<- kg.mouse$kg.sets[kg.mouse$sigmet.idx]
}
})
return(kegg.gs)
})
go.gs <- reactive({
species <- input$species
withProgress(message = 'Extracting GO term information', value = NULL, {
if(species == "Human" | species == "human" | species == "hg19" | species == "hsa"){
go.gs <- go.gsets(species="human", pkg.name=NULL, id.type = "eg")
}else if(species == "Mouse" | species == "mouse" | species == "mm9" | species == "mmu"){
go.gs <- go.gsets(species="mouse", pkg.name=NULL, id.type = "eg")
}
})
return(go.gs)
})
gage.exp.fc <- eventReactive(input$run_enrich, {
rankdata <- nmf_group_feature_rank()
num_groups <- ncol(rankdata)
emp <- matrix(NA, nrow(rankdata), 1)
res <- rep(list(emp), num_groups)
names(res) <- colnames(rankdata)
withProgress(message = 'Calculating LogFC between groups', value = NULL, {
for (i in 1:num_groups){
# In logFC, error might come from the fact that some of them are Inf or -Inf
idx <- is.finite(rankdata[, i])
sub_rankdata <- rankdata[idx, i, drop=F]
# ID conversion
id.map.sym2eg <- id2eg(ids=rownames(sub_rankdata), category = "SYMBOL", org=id.org())
gene.entrez <- mol.sum(mol.data = sub_rankdata, id.map = id.map.sym2eg, sum.method = "mean")
res[[i]] <- gene.entrez[, 1]
res[[i]] <- res[[i]][!is.na(res[[i]])]
}
})
return(res)
})
go.Res <- reactive({
gage.exp.fc <- gage.exp.fc()
gsets <- NULL
res <- rep(list(NA), length(gage.exp.fc))
names(res) <- names(gage.exp.fc)
if(input$EnrichType == 'GO'){
gsets <- go.gs()
gsets <- gsets$go.sets[ gsets$go.subs[[input$goTerm]] ]
}else if(input$EnrichType == 'KEGG'){
gsets <- kegg.gs()
}
validate(
need(length(gsets)>0, "Did not find any Gene Sets, maybe your gene count data is not from 'Human'?")
)
withProgress(message = 'Running enrichment analysis', value = NULL, {
for (i in 1:length(gage.exp.fc)){
exp.fc <- gage.exp.fc[[i]]
goRes <- gage(exp.fc, gsets = gsets, ref = NULL, samp = NULL, same.dir=ifelse(input$samedir, TRUE, FALSE))
res[[i]] <- goRes
}
})
return(res)
})
go_summary <- reactive({
goRes <- go.Res()
genes <- goRes[[1]]$greater %>% rownames
res <- list()
num_genes <- 5
hi_res <- lo_res <- matrix(NA, num_genes * length(goRes), length(goRes))
colnames(hi_res) <- colnames(lo_res) <- names(goRes)
total_hi_names <- total_lo_names <- rep(NA, num_genes * length(goRes))
for(i in 1:length(goRes)){
range <- ((i-1) * num_genes + 1) : (i * num_genes)
hi_genenames <- goRes[[i]]$greater %>% rownames %>% head(n = num_genes)
hi_res[range, ] <- sapply(goRes, function(x) x$greater[hi_genenames, "p.val"])
total_hi_names[range] <- hi_genenames
if(!is.na(goRes[[i]]$less[1])){
lo_genenames <- goRes[[i]]$less %>% rownames %>% head(n = num_genes)
lo_res[range, ] <- sapply(goRes, function(x) x$less[lo_genenames, "p.val"])
total_lo_names[range] <- lo_genenames
}
}
rownames(hi_res) <- total_hi_names
rownames(lo_res) <- total_lo_names
# save only the unique GO:Term
res[["hi_sub_pval"]] <- hi_res[!duplicated(rownames(hi_res)), ]
res[["lo_sub_pval"]] <- lo_res[!duplicated(rownames(lo_res)), ]
# save total results
res[["hi_total_pval"]] <- sapply(goRes, function(x) x$greater[genes, "p.val"])
res[["hi_total_qval"]] <- sapply(goRes, function(x) x$greater[genes, "q.val"])
res[["lo_total_pval"]] <- sapply(goRes, function(x) x$less[genes, "p.val"])
res[["lo_total_qval"]] <- sapply(goRes, function(x) x$less[genes, "q.val"])
return(res)
})
output$goplot_hi <- renderPlotly({
req(go_summary())
validate(
need(!is.na(go_summary()[["hi_sub_pval"]]), "Did not find any GO terms for greater, maybe you select the wrong 'species'?")
)
withProgress(message = 'Generating summary plot for GO term (Hi)', value = NULL, {
aa <- isolate(
melt(go_summary()[["hi_sub_pval"]]) %>%
set_colnames(c("GO", "NMF", "P")) %>%
mutate(P = replace(P, is.na(P), 1))
)
if(input$go_logscale){
p <- plot_ly(aa, x = ~NMF, y = ~GO, z = ~log10(P), type = "heatmap", colors = input$go_colorscale)
}else{
p <- plot_ly(aa, x = ~NMF, y = ~GO, z = ~P, type = "heatmap", colors = input$go_colorscale)
}
m <- list( l = input$go_leftmargin, r = 10, b = 40, t = 20)
l <- list(margin = m,
xaxis = list(title = "", tickfont=list(family = "Arial", size = input$go_xfontsize), zeroline=FALSE),
yaxis = list(title = "", tickfont=list(family = "Arial", size = input$go_yfontsize), zeroline=FALSE))
colorbar_format <- list(family = "Arial", size = 8)
naikai <- plotly_build(p)
naikai$x$data[[1]]$colorbar <- list(title = "P", thickness = input$go_barwidth, len=0.5, titlefont=colorbar_format, tickfont=colorbar_format)
naikai$x$layout <- l
naikai
})
})
output$goplot_lo <- renderPlotly({
req(go_summary())
validate(
need(!is.na(go_summary()[["lo_sub_pval"]]), "Did not find any GO terms for less, maybe you select the wrong 'species'?")
)
withProgress(message = 'Generating summary plot for GO term (Low)', value = NULL, {
aa <- isolate(
melt(go_summary()[["lo_sub_pval"]]) %>%
set_colnames(c("GO", "NMF", "P")) %>%
mutate(P = replace(P, is.na(P), 1))
)
if(input$go_logscale){
p <- plot_ly(aa, x = ~NMF, y = ~GO, z = ~log10(P), type = "heatmap", colors = input$go_colorscale)
}else{
p <- plot_ly(aa, x = ~NMF, y = ~GO, z = ~P, type = "heatmap", colors = input$go_colorscale)
}
m <- list( l = input$go_leftmargin, r = 10, b = 40, t = 20)
l <- list(margin = m,
xaxis = list(title = "", tickfont=list(family = "Arial", size = input$go_xfontsize), zeroline=FALSE),
yaxis = list(title = "", tickfont=list(family = "Arial", size = input$go_yfontsize), zeroline=FALSE))
colorbar_format <- list(family = "Arial", size = 8)
naikai <- plotly_build(p)
naikai$x$data[[1]]$colorbar <- list(title = "P", thickness = input$go_barwidth, len=0.5, titlefont=colorbar_format, tickfont=colorbar_format)
naikai$x$layout <- l
naikai
})
})
go_table <- reactive({
goRes <- go.Res()
validate(
need(input$pathway_group!="", "Please select the enrichment result for the group you are interest in")
)
goRes <- goRes[[input$pathway_group]]
if (is.null(goRes$greater)){
stop("No significant gene sets found using this q.val cutoff")
}else{
return(goRes$greater)
}
})
output$go_summary <- DT::renderDataTable({
filename <- paste(input$EnrichType, input$pathway_group, input$goTerm, sep = "_")
idx <- go_table()[, "q.val"] <= input$qval_cutoff
go_table <- go_table()[which(!is.na(idx) & idx), ]
if(nrow(go_table) == 0){
stop("No significant gene sets found using this q.val cutoff")
}
go_table <- go_table[, c(3,4,2,5)]
datatable(go_table, selection = 'single',
extensions = c('Buttons', 'Responsive'),
options = list(dom = 'Bfrtip',
buttons =
list('colvis', 'copy', 'print', list(
extend = 'collection',
buttons = list(list(extend='csv',
filename = filename),
list(extend='excel',
filename = filename)),
text = 'Download'
)),
pageLength = 6,
autoWidth = TRUE
# columnDefs = list(list(width = '10px', targets = "_all"))
)
) %>% formatRound(1:3, 3)
}, server = FALSE)
output$pathviewImg <- renderImage({
go_table <- go_table()
s.idx <- input$go_summary_rows_selected[length(input$go_summary_rows_selected)]
validate(
need(length(s.idx)>0, "Please select one KEGG pathway from above.")
)
withProgress(message = 'Rendering pathview map', value = NULL, {
path.pid <- substr(rownames(go_table)[s.idx], 1, 8)
out.suffix <- "nmf"
pathview(gene.data=gage.exp.fc()[[input$pathway_group]], pathway.id=path.pid, species=pathview.species(), out.suffix=out.suffix)
})
filename <- path.pid
outfile <- paste(filename, out.suffix, 'png', sep=".")
# Return a list containing the filename
list(src = outfile,
contentType = 'image/png',
width = 1000,
height = 777,
alt = "This is alternate text")
}, deleteFile = FALSE)
nmf_feature_rank_for_gsea <- reactive({
species <- input$species
nmf_feature_rank <- nmf_group_feature_rank()
# check if genenames contains all CAPITAL LETTERS (after removing all numeric characters)
cap_genes <- grepl("^[[:upper:]]+$", gsub('[[:digit:]]+', '', rownames(nmf_feature_rank))) %>% sum
cap_pct <- cap_genes/ nrow(nmf_feature_rank)
# some of the genes might contains '.', '-', 'orf', etc, so we use percentage of genes qualified as filtering criteria
# if(species == "Mouse" | species == "mouse" | species == "mm9" | species == "mmu"){
if(cap_pct < 0.1){
n <- 2
withProgress(message = 'Assuming the input gene symbols are from Mouse, Converting to Human..', value = 0, {
incProgress(1/n, detail = "Takes around 10~15 seconds")
mouse_human_convert_genes <- gene_convert_mouse_to_human(rownames(nmf_feature_rank))
## After conversion, some of the mouse genes might point to the same human gene
## So select thoide mouse genes that has highest mean expression across groups?
idx <- duplicated(mouse_human_convert_genes$HGNC.symbol)
if (sum(idx)>0){
warning("Some of the mouse genes are mapped back to the same human genes..")
}
max_dup_mean_lfc_genes <-
cbind(mouse_human_convert_genes, nmf_feature_rank) %>%
mutate(avg_lfc=apply(nmf_feature_rank, 1, mean),
abs_lfc=apply(nmf_feature_rank, 1, function(x) max(abs(x)))
) %>%
dplyr::group_by(HGNC.symbol) %>%
dplyr::filter(min_rank(desc(avg_lfc))==1) %>%
dplyr::slice(which.max(abs_lfc)) %>%
ungroup # this will give Error: corrupt 'grouped_df', workaround is to ungroup it below
print(dim(max_dup_mean_lfc_genes))
## Return final de-duplicated genes only
nmf_feature_rank <- max_dup_mean_lfc_genes %>%
dplyr::select(starts_with("NMF")) %>%
as.data.frame
rownames(nmf_feature_rank) <- max_dup_mean_lfc_genes$HGNC.symbol
})
}
return(nmf_feature_rank)
})
output$downloadPathData <- downloadHandler(
filename <- function() {
paste( file_prefix(), "pathway.zip", sep=".")
},
content = function(file) {
tmpdir <- tempdir()
current_dir <- getwd()
setwd(tempdir())
### Save GO Term results
go_summary <- go_summary()
total_idx <- grep("total", names(go_summary))
go_filenames <- paste(file_prefix(), paste0("Enrichment_", names(go_summary)[total_idx]), "txt", sep=".")
for(i in 1:length(total_idx)){
write.table(go_summary[[ total_idx[i] ]], go_filenames[i], quote=F, row.names = TRUE, col.names=NA, sep="\t")
}
### Rank file for GSEA
# modify it to save one file for each rank ###
nmf_group_feature_rank <- nmf_feature_rank_for_gsea()
num_samples <- ncol(nmf_group_feature_rank)
filenames <- paste(file_prefix(), paste0("group_feature_rank", 1:num_samples), "rnk", sep=".")
for (i in 1:num_samples){
sub.data <- data.frame(Gene=rownames(nmf_group_feature_rank), LogFC=nmf_group_feature_rank[,i])
idx <- is.finite(sub.data$LogFC)
write.table(sub.data[idx, ], filenames[i], quote=F, row.names = F, sep="\t")
}
zip(zipfile=file, files=c(go_filenames, filenames))
setwd(as.character(current_dir))
},
contentType = "application/zip"
)
})
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.