Commit 274fe372 authored by Perrine Woringer's avatar Perrine Woringer
Browse files

Add network plots, and many improvements

parent 1a37f5fa
options(download.file.method = 'wget', bitmapType='cairo')
#options(download.file.method = 'wget', bitmapType='cairo')
options(bitmapType='cairo')
if (!require("Rcpp")){
install.packages("Rcpp", repos="http://cran.univ-paris1.fr/")
......@@ -83,7 +84,8 @@ library(biomformat)
}
if (!require(scatterD3)) {
devtools::install_github('aghozlane/scatterD3')
#devtools::install_github('aghozlane/scatterD3')
install.packages("scatterD3")
library(scatterD3)
}
......@@ -109,7 +111,7 @@ if (!require(shinyjs)) {
}
if(!require(d3vennR)){
install_github("timelyportfolio/d3vennR")
devtools::install_github("timelyportfolio/d3vennR")
library(d3vennR)
}
......@@ -129,8 +131,7 @@ if (!require(gplots)) {
}
if (!require(DESeq2)) {
source("https://bioconductor.org/biocLite.R")
biocLite("DESeq2")
devtools::install_github("aghozlane/DESeq2shaman")
library(DESeq2)
}
......@@ -230,8 +231,8 @@ if (!require("rAmCharts")){
if(!require("colourpicker")){
install.packages("colourpicker")
library(coulourpicker)
}
library(colourpicker)
}
if(!require("data.table")){
install.packages("data.table")
......@@ -246,4 +247,19 @@ if(!require("UpSetR")){
if(!require("ggrepel")){
install.packages("ggrepel")
library(ggrepel)
}
\ No newline at end of file
}
# if(!require("networkD3")){
# install.packages("networkD3")
# library(networkD3)
# }
if(!require("igraph")){
install.packages("igraph")
library(igraph)
}
if(!require("visNetwork")){
install.packages("visNetwork")
library(visNetwork)
}
......@@ -46,7 +46,7 @@ Plot_Visu_Heatmap_FC <- function(input, BaseContrast, resDiff, ContrastListDebou
##############################
## P VALUE DENSITY PLOT
##############################
Plot_pValue_Density <- function(input, BaseContrast, resDiff, ContrastListDebounce){
Plot_pValue_Density <- function(input, BaseContrast, resDiff, ContrastListDebounce, alphaVal){
res = NULL
#SelContrast = input$ContrastList_table_FC
SelContrast = ContrastListDebounce()
......@@ -64,7 +64,8 @@ Plot_pValue_Density <- function(input, BaseContrast, resDiff, ContrastListDeboun
data <- rbind(data, data_cont)
}
data$contrast <- factor(data$contrast, levels = SelContrast)
p <- ggplot(data, aes(x = padj, color = contrast, fill = contrast)) + geom_density(alpha = input$fillOpacity, size = input$lineWidth)
# uncomment to zoom (X axis between 0 and input$alphaVal) (use "+ xlim(0,as.numeric(alphaVal))" instead, to redraw geom_density with only data under threshold)
p <- ggplot(data, aes(x = padj, color = contrast, fill = contrast)) + geom_density(alpha = input$fillOpacity, size = input$lineWidth) + theme_minimal() #+ coord_cartesian(xlim = c(0,as.numeric(alphaVal)))
p <- p + theme(axis.title = element_text(size = input$FontSizepValueDensity),
axis.text = element_text(size = input$FontSizepValueDensity),
legend.title = element_text(size = input$FontSizepValueDensity),
......@@ -187,7 +188,7 @@ Plot_Comp_Logit <- function(input, BaseContrast, resDiff, SelectTaxoPlotCompDebo
if(export){
data <- data.frame(logit_x, logit_y, color_var)
plot <- ggplot(data, aes(x = logit_x, y = logit_y))
plot <- ggplot(data, aes(x = logit_x, y = logit_y)) + theme_minimal()
plot <- plot + geom_point(aes(color = color_var), size = (input$pointSizeLogit %/% 30), alpha = input$pointOpacityLogit)
plot <- plot + scale_x_continuous(paste("logit p value", input$Contrast1))
plot <- plot + scale_y_continuous(paste("logit p value", input$Contrast2))
......@@ -216,20 +217,23 @@ Plot_Comp_Logit <- function(input, BaseContrast, resDiff, SelectTaxoPlotCompDebo
##############################
## VENN DIAGRAM
##############################
Plot_Visu_Venn <- function(input,BaseContrast,resDiff, ContrastListDebounce, export=FALSE){
Plot_Visu_Venn <- function(input,BaseContrast,resDiff, ContrastListVennDebounce, export=FALSE){
res = NULL
#SelContrast = input$ContrastList_table_FC
SelContrast = ContrastListDebounce()
#SelContrast = input$ContrastList_table_FCVenn
SelContrast = ContrastListVennDebounce()
if(length(SelContrast)>=2 & length(SelContrast)<=4){
data = GetData_venn(input,SelContrast,BaseContrast,resDiff)$res
res = venn_tooltip(d3vennR(data=data))
}
return(res)
}
##############################
## CONTRASTS COMPAIR
## UPSET
##############################
Plot_UpSet <- function(input,BaseContrast, resDiff, ContrastListDebounce, export=FALSE){
plot = NULL
df = NULL
#SelContrast = input$ContrastList_table_FC
SelContrast = ContrastListDebounce()
if(length(SelContrast)>=2){
......@@ -240,16 +244,38 @@ Plot_UpSet <- function(input,BaseContrast, resDiff, ContrastListDebounce, export
new <- list(na.omit(data[,i]))
listInput <- c(listInput, new)
}
print(input$pointSizeContrastsCompair)
print(class(input$pointSizeContrastsCompair))
names(listInput) <- colnames(data)
# for plot
plot <- upset(fromList(listInput), nsets = n, nintersects = NA, set_size.show = input$showNumbers,
show.numbers = if(input$showNumbers){"yes"}else{"no"}, sets = SelContrast, keep.order = !(input$orderBySize),
point.size = input$pointSizeContrastsCompair, line.size = 1.5, order.by = if(input$orderBySize){"freq"}else{"degree"},
point.size = input$pointSizeUpSet, line.size = 1.5, order.by = input$orderByUpset,
mainbar.y.label = "Number of differential features in common", sets.x.label = "Number of differential features",
text.scale = c(2, 1.6, 1.6, 1.6, 2, 1.5))
text.scale = c(2, 1.6, 1.6, 1.6, 2, 1.5), decreasing = input$decreasingUpset)
# for table
combos <- Reduce(c,lapply(rev(1:length(listInput)), function(x) combn(rev(colnames(data)),x,simplify=FALSE) ))
listSignif <- Reduce(function(a,b) {unique(c(as.character(a),as.character(b)))}, listInput)
n_col = length(combos)
n_row = length(listSignif)
df <- matrix(rep(listSignif, times = n_col), ncol = n_col)
lst <- list()
for(j in 1:n_col){for(i in 1:n_row){group_contrasts <- unlist(combos[j])
if (is.element(df[i,j], lst)){df[i,j] <- NA}
else{if (!is.element(df[i,j], Reduce(intersect, listInput[group_contrasts]))){df[i,j] <- NA}
else{lst <- c(lst,df[i,j])}
}
}
}
colnames(df) <- lapply(combos, function(lstCont) paste(lstCont, collapse = " & "))
df = as.data.frame(apply(df,2,Go_data_top))
maxRow = max(apply(df,2,FUN=function(x) length(which(!is.na(x)))))
df = df[1:max(maxRow,1),]
df = df[,which(apply(!is.na(df),2,any))]
}
return(plot)
return(list(plot=plot,table=df))
}
......@@ -463,10 +489,24 @@ GetData_venn <-function(input,SelContrast,BaseContrast,resDiff)
for(i in 1:(ncont))
{
for(j in i:ncont)
{
if(i!=j) res[[cmp]] = list(sets=list(names.df[i],names.df[j]),size= length(which(!is.na(intersect(df[,i],df[,j])))))
if(i==j) res[[cmp]] = list(sets=list(names.df[i]),size= length(which(!is.na(intersect(df[,i],df[,i])))))
cmp=cmp+1
# For Venn diagram (until intersection of 4 contrasts)
{ if(i==j) {res[[cmp]] = list(sets=list(names.df[i]),size= length(which(!is.na(intersect(df[,i],df[,i])))))
cmp=cmp+1}
if(i!=j) {res[[cmp]] = list(sets=list(names.df[i],names.df[j]),size= length(which(!is.na(intersect(df[,i],df[,j])))))
cmp=cmp+1
for(k in j:ncont){
if(k!=j){
res[[cmp]] = list(sets=list(names.df[i],names.df[j],names.df[k]), size = length(which(!is.na(intersect(intersect(df[,i],df[,j]),df[,k])))))
cmp=cmp+1
for(l in k:ncont){
if(l!=k){
res[[cmp]] = list(sets=list(names.df[i],names.df[j],names.df[k],names.df[l]), size = length(which(!is.na(intersect(intersect(intersect(df[,i],df[,j]),df[,k]),df[,l])))))
cmp=cmp+1}
}
}
}
}
# For 'Multiple Venn Compair'
if(i!=j) {
if (is.null(res_multiple_venn))
{res_multiple_venn <- data.frame(name = c(rownames_multiple_venn, paste(names.df[i],names.df[j], sep = " vs ")),
......
......@@ -47,22 +47,22 @@ Get_dds_object <- function(input,counts,target,design,normFactorsOTU,CT_noNorm,C
# }
# countsNorm = counts(dds, normalized = TRUE)
if(as.numeric(R.Version()$major)>=3 && as.numeric(R.Version()$minor) >=1.3){
if(nrow(counts)*nrow(target)>=50000)
{
try(DESeq(dds,fitType=input$fitType,parallel = TRUE,minReplicatesForReplace = Inf) ->ddstmp,silent = TRUE)
if(!is.null(ddstmp)) dds = ddstmp
}
if(nrow(counts)*nrow(target)<50000 || is.null(ddstmp))
{
dds <- estimateDispersions(dds, fitType=input$fitType)
dds <- nbinomWaldTest(dds)
}
}else{
# if(as.numeric(R.Version()$major)>=3 && as.numeric(R.Version()$minor) >=1.3){
# if(nrow(counts)*nrow(target)>=50000)
# {
# try(DESeq(dds,fitType=input$fitType,parallel = TRUE,minReplicatesForReplace = Inf) ->ddstmp,silent = TRUE)
# if(!is.null(ddstmp)) dds = ddstmp
# }
# if(nrow(counts)*nrow(target)<50000 || is.null(ddstmp))
# {
# dds <- estimateDispersions(dds, fitType=input$fitType)
# dds <- nbinomWaldTest(dds)
# }
# }else{
dds <- estimateDispersions(dds, fitType=input$fitType)
dds <- nbinomWaldTest(dds,modelMatrixType = "expanded")
# dds <- DESeq(dds,fitType=input$fitType,modelMatrixType = "expanded",parallel = TRUE)
}
# }
countsNorm = counts(dds, normalized = TRUE)
......
......@@ -73,15 +73,13 @@ Bar_Chart_Tables <- function(input, data, export = FALSE) {
creditsPosition = "bottom-right"
)
plot <- setProperties(plot, fontSize = input$fontSize)
return(plot)
}
if(export){
print(input$fontSize)
else{
data <- data_table
data$Id <- factor(data$Id, levels = rev(data$Id))
plot <- ggplot(data, aes(x = Id, y = log2FoldChange))
plot <- plot + geom_bar(stat = "identity", aes(color = Group, fill = Group))
plot <- ggplot(data, aes(x = Id, y = log2FoldChange)) + theme_minimal()
plot <- plot + geom_bar(stat = "identity", aes(color = Group, fill = Group), width = 0.7)
plot <- plot + coord_flip()
plot <- plot + scale_color_manual(values = c("Down" = input$colour1,"Not significant" = input$colour2,"Up" = input$colour3))
plot <- plot + scale_fill_manual(values = c("Down" = input$colour1,"Not significant" = input$colour2,"Up" = input$colour3))
......@@ -90,11 +88,10 @@ Bar_Chart_Tables <- function(input, data, export = FALSE) {
axis.text = element_text(size = input$fontSize),
legend.title = element_text(size = input$fontSize),
legend.text = element_text(size = input$fontSize))
if(input$removeLabelsBarChart){print("ici")
plot <- plot + theme(axis.ticks.y = element_blank(),
axis.text.y = element_blank())}
return(plot)
if(input$removeLabelsBarChart){plot <- plot + theme(axis.ticks.y = element_blank(),
axis.text.y = element_blank())}
}
return(plot)
}
###########################
......@@ -186,8 +183,7 @@ Volcano_Plot <- function(input, data, export = FALSE) {
}
else{NULL},
col_var = color_var,
col_lab = NA #if(input$showLegend){"Legend title"}else{NA}
,
col_lab = NA, #if(input$showLegend){"Legend title"}else{NA}
colors = c(
"Down" = input$colour1,
"Not significant" = input$colour2,
......@@ -205,14 +201,14 @@ Volcano_Plot <- function(input, data, export = FALSE) {
menu = FALSE,
height = input$heightVolcanoPlot,
width=if(input$modifwidthVolcano){input$widthVolcanoPlot},
disable_wheel = TRUE
disable_wheel = FALSE
)
return(plot)
}
if(export){
data <- data.frame(x,y_log, color_var)
plot <- ggplot(data, aes(x = x, y = y_log))
plot <- ggplot(data, aes(x = x, y = y_log)) + theme_minimal()
plot <- plot + geom_point(aes(color = color_var), size = (input$pointSize %/% 30), alpha = input$pointOpacity)
plot <- plot + scale_x_continuous("Log 2 Fold Change", limits = xlimits)
plot <- plot + scale_y_continuous("- Log 10 p value adjusted", limits = ylimits)
......
......@@ -449,17 +449,17 @@ Plot_Visu_Diversity <- function(input,resDiff,ForScatter=FALSE){
ci.alpha.down = pmax(alpha - 1.96*tapply(TaxoNumber(counts_tmp_combined), targetInt$AllVar, sd)/sqrt.nb,0)
ci.alpha.up = alpha + 1.96*tapply(TaxoNumber(counts_tmp_combined), targetInt$AllVar, sd)/sqrt.nb
shan <- tapply(diversity(counts_tmp_combined, index = "shannon"), targetInt$AllVar, mean)
ci.shan.down = pmax(shan - 1.96*tapply(diversity(counts_tmp_combined, index = "shannon"), targetInt$AllVar, sd)/sqrt.nb,0)
ci.shan.up = shan + 1.96*tapply(diversity(counts_tmp_combined, index = "shannon"), targetInt$AllVar, sd)/sqrt.nb
shan <- tapply(vegan::diversity(counts_tmp_combined, index = "shannon"), targetInt$AllVar, mean)
ci.shan.down = pmax(shan - 1.96*tapply(vegan::diversity(counts_tmp_combined, index = "shannon"), targetInt$AllVar, sd)/sqrt.nb,0)
ci.shan.up = shan + 1.96*tapply(vegan::diversity(counts_tmp_combined, index = "shannon"), targetInt$AllVar, sd)/sqrt.nb
simpson <- tapply(diversity(counts_tmp_combined, index = "simpson"), targetInt$AllVar, mean)
ci.simpson.down = pmax(simpson - 1.96*tapply(diversity(counts_tmp_combined, index = "simpson"), targetInt$AllVar, sd)/sqrt.nb,0)
ci.simpson.up = simpson + 1.96*tapply(diversity(counts_tmp_combined, index = "simpson"), targetInt$AllVar, sd)/sqrt.nb
simpson <- tapply(vegan::diversity(counts_tmp_combined, index = "simpson"), targetInt$AllVar, mean)
ci.simpson.down = pmax(simpson - 1.96*tapply(vegan::diversity(counts_tmp_combined, index = "simpson"), targetInt$AllVar, sd)/sqrt.nb,0)
ci.simpson.up = simpson + 1.96*tapply(vegan::diversity(counts_tmp_combined, index = "simpson"), targetInt$AllVar, sd)/sqrt.nb
invsimpson <- tapply(diversity(counts_tmp_combined, index = "invsimpson"), targetInt$AllVar, mean)
ci.invsimpson.down = pmax(invsimpson - 1.96*tapply(diversity(counts_tmp_combined, index = "invsimpson"), targetInt$AllVar, sd)/sqrt.nb,0)
ci.invsimpson.up = invsimpson + 1.96*tapply(diversity(counts_tmp_combined, index = "invsimpson"), targetInt$AllVar, sd)/sqrt.nb
invsimpson <- tapply(vegan::diversity(counts_tmp_combined, index = "invsimpson"), targetInt$AllVar, mean)
ci.invsimpson.down = pmax(invsimpson - 1.96*tapply(vegan::diversity(counts_tmp_combined, index = "invsimpson"), targetInt$AllVar, sd)/sqrt.nb,0)
ci.invsimpson.up = invsimpson + 1.96*tapply(vegan::diversity(counts_tmp_combined, index = "invsimpson"), targetInt$AllVar, sd)/sqrt.nb
gamma <- TaxoNumber(counts_tmp_combined, targetInt$AllVar)
beta = gamma/alpha - 1
......@@ -605,8 +605,12 @@ expand.grid2.list <- function(listInput)
## Put the data in the right format to be plot
GetDataToPlot <- function(input,resDiff,VarInt,ind_taxo,aggregate=TRUE,rarefy=FALSE)
GetDataToPlot <- function(input,resDiff,VarInt,ind_taxo,sec_variable = NULL, aggregate=TRUE,rarefy=FALSE)
{
sec_variable_added_to_VarInt <- FALSE
if(!is.null(sec_variable)){if(!is.element(sec_variable,VarInt)){VarInt <- c(VarInt,sec_variable)
sec_variable_added_to_VarInt <- TRUE}}
dds = resDiff$dds
val = c()
list.val = list()
......@@ -640,7 +644,10 @@ GetDataToPlot <- function(input,resDiff,VarInt,ind_taxo,aggregate=TRUE,rarefy=FA
expr=parse(text=Tinput)
## All the modalities for all the var of interest
val = c(val,eval(expr))
list.val[[i]] = eval(expr)
if(sec_variable_added_to_VarInt){mod <- eval(expr)
if(is.null(mod)){mod <- as.character(unique(as.factor(target[,VarInt[i]])))}
list.val[[i]] = mod}
else{list.val[[i]] = eval(expr)}
}
if (!is.null(val) && !is.null(list.val))
{
......@@ -699,9 +706,7 @@ GetDataToPlot <- function(input,resDiff,VarInt,ind_taxo,aggregate=TRUE,rarefy=FA
}
}
}
return(list(counts = counts_tmp_combined,targetInt=targetInt,prop=prop_tmp_combined,namesCounts=namesCounts,levelsMod=levelsMod,prop_all=prop_all))
return(list(counts=counts_tmp_combined,targetInt=targetInt,prop=prop_tmp_combined,namesCounts=namesCounts,levelsMod=levelsMod,prop_all=prop_all))
}
......@@ -830,3 +835,78 @@ Plot_Visu_Tree <- function(input,resDiff,CT_Norm_OTU,taxo_table)
return(res)
}
#############################
## NETWORK
#############################
Plot_network <- function(input,resDiff,availableTaxo, ind_taxo, qualiVariable, export = FALSE){
plot = NULL
dataVN = NULL
VarInt = input$VisuVarInt
if(isolate(input$colorCorr)){sec_variable = isolate(input$sec_variable)}
else{sec_variable = NULL}
data <- GetDataToPlot(input,resDiff,VarInt,availableTaxo, sec_variable = sec_variable, aggregate = FALSE)
if(!is.null(data)){
counts_tmp_combined <- data$counts
dataVariables <- as.matrix(data$targetInt)
if(isolate(input$colorCorr && qualiVariable())){dataVariables[,sec_variable] <- sapply(dataVariables[,sec_variable], function(x) if(is.element(x,isolate(input$values1))){1}else{0})}
if(!is.null(counts_tmp_combined)){
countsMatrix <- as.matrix(counts_tmp_combined)
n <- ncol(countsMatrix)
resCorrTest <- corr.test(countsMatrix, ci = FALSE)
cor <- resCorrTest$r
pval <- resCorrTest$p
pval_bool <- t(apply(pval, 1, function(v) {sapply(v, function(x){x < 0.05})}))
cor_sgn <- t(apply(cor, 1, function(v) {sapply(v, sign)}))
adjacency <- matrix(mapply(function(a,b) {mapply(function(x,y){x*y}, x=a, y=b)}, a=pval_bool, b=cor_sgn), nrow = n)
rownames(adjacency) <- colnames(countsMatrix)
colnames(adjacency) <- colnames(countsMatrix)
# ### Remove rows and columns with only NA
# adjacency <- adjacency[apply(adjacency, 1, function(y) !all(is.na(y))),]
# adjacency <- t(adjacency)
# adjacency <- adjacency[apply(adjacency, 1, function(y) !all(is.na(y))),]
# adjacency <- t(adjacency)
### Replace NA by zeros (ie "no correlation")
adjacency[is.na(adjacency)] <- 0
adjacency <- adjacency[,ind_taxo]
adjacency <- adjacency[ind_taxo,]
igraphGraph <- graph_from_adjacency_matrix(adjacency, diag = FALSE, mode = "upper" , weighted = TRUE) # "upper" for adjusted p-value, lower for p-value not adjusted
list_to_label <- isolate(input$ToLabelNetwork)
dataVN <- toVisNetworkData(igraphGraph)
dataVN$nodes$title <- paste0("<b>", dataVN$nodes$id,"</b>")
dataVN$nodes$label <- sapply(dataVN$nodes$id, function(x)if(is.element(x, list_to_label)){x}else{""})
dataVN$edges$color <- sapply(dataVN$edges$weight, function(x)if(x==1){isolate(input$edgeColorPositive)}else{isolate(input$edgeColorNegative)})
if(!is.null(sec_variable)){
cor <- sapply(dataVN$nodes$id, function(x)cor(as.numeric(dataVariables[,sec_variable]), countsMatrix[,x]))
dataVN$nodes$cor <- round(cor, digits = 5)
scale <- if(isolate(input$scaleFree)){max(c(max(dataVN$nodes$cor),-min(dataVN$nodes$cor)))}else{1}
dataVN$nodes$color.background <- sapply(dataVN$nodes$cor, function(x) colorRampPalette(rev(brewer.pal(9, isolate(input$colorPalette))))(200)[round(x, digits = 2)*100/scale+100])
dataVN$nodes$color.highlight.background <- dataVN$nodes$color.background
}
else{dataVN$nodes$color.background <- isolate(input$colorBackground)
dataVN$nodes$color.highlight.background <- isolate(input$colorHighlightBackground)}
dataVN$nodes$color.border <- isolate(input$colorBorder)
dataVN$nodes$color.highlight.border <- isolate(input$colorHighlightBorder)
plot <- visNetwork(nodes = dataVN$nodes, edges = dataVN$edges)
plot <- visIgraphLayout(plot, layout = "layout_nicely", physics = FALSE, smooth = FALSE)
plot <- visNodes(plot, size = 20) #, scaling = list(label = list(min = 30, max = 30, maxVisible = 30)))
plot <- visEdges(plot, width = 1)
plot <- visOptions(plot, width = if(isolate(input$modifwidthVisu)){isolate(input$widthVisu)}, height = isolate(input$heightVisu), autoResize = FALSE)
#plot <- visLegend(plot, addEdges = data.frame(color = c("red", "blue"), label = c("Positive correlation","Negative correlation")))
#plot <- visExport(plot, type = "pdf", name = "network_SHAMAN.pdf", float="bottom")
}
}
return(list(plot = plot, data = dataVN))
}
This diff is collapsed.
This diff is collapsed.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment