Skip to content

GitLab

Commit d3d4f423 authored by svolant's avatar svolant
Browse files

NMDS + permanova

parent 956df065
Hide whitespace changes
Inline Side-by-side
Rfunctions/DiagPlot.R
View file @ d3d4f423
......@@ -40,15 +40,53 @@ Plot_diag <- function(input,resDiff,getTable=FALSE)
if(input$DiagPlot=="SfactorsVStot") res = diagSFactors(input,normFactors,resDiff$raw_counts)
if(input$DiagPlot=="pcaPlot") res = PCAPlot_meta(input,dds, group, type.trans = input$TransType, col = colors)
if(input$DiagPlot=="pcoaPlot") res = PCoAPlot_meta(input,dds, group, col = colors)
if(input$DiagPlot=="nmdsPlot") res = NMDSPlot(input, dds, group, col = colors)
if(input$DiagPlot=="clustPlot") res = HCPlot(input,dds,group,type.trans=input$TransType,counts,col=colors)
}
if(getTable && input$DiagPlot=="pcaPlot") res = Get_pca_table(input,dds, group, type.trans = input$TransType)
if(getTable && input$DiagPlot=="pcoaPlot") res = Get_pcoa_table(input,dds, group)
if(getTable && input$DiagPlot=="pcoaPlot") res = Get_pcoa_table(input,dds, group)$table
if(getTable && input$DiagPlot=="nmdsPlot") res = Get_nmds_table(input,dds, group)$table
return(res)
}
##############################################################
##
## Permanova test
##
##############################################################
Perma_test_Diag <- function(input,resDiff)
{
VarInt = input$VarInt
dds = resDiff$dds
counts = resDiff$raw_counts
if(input$CountsType=="Normalized") counts = resDiff$countsNorm
target = resDiff$target
normFactors = resDiff$normFactors
## Counts at the OTU level
CT = resDiff$CT_noNorm
if(input$CountsType=="Normalized") CT = resDiff$CT_Norm
group = as.data.frame(target[,VarInt])
rownames(group) = rownames(target)
res = NULL
if(ncol(group)>0 && !is.null(dds))
{
if(input$DiagPlot=="pcoaPlot") res = Get_pcoa_table(input,dds, group)$test$aov.tab$`Pr(>F)`[1]
if(input$DiagPlot=="nmdsPlot") res = Get_nmds_table(input,dds, group)$test$aov.tab$`Pr(>F)`[1]
}
return(res)
}
##############################################################
##
## Plot functions
......@@ -463,9 +501,9 @@ PCAPlot_meta <-function(input,dds, group_init, n = min(500, nrow(counts(dds))),
prp <- round(prp, 2)
ncol1 <- ncol(group) == 1
abs = range(pca$x[, 1])
abs = range(proj[, as.numeric(gsub("PC","",input$PCaxe1))])
abs = abs(abs[2] - abs[1])/25
ord = range(pca$x[, 2])
ord = range(proj[, as.numeric(gsub("PC","",input$PCaxe2))])
ord = abs(ord[2] - ord[1])/25
......@@ -583,11 +621,12 @@ Get_pcoa_table <-function (input, dds, group_init)
if(input$DistClust!="sere") dist.counts.norm = vegdist(t(counts.norm), method = input$DistClust)
if(input$DistClust=="sere") dist.counts.norm = as.dist(SEREcoef(counts.norm))
permanova_test = adonis(dist.counts.norm~group)
print(permanova_test)
## Do PCoA
pco.counts.norm = dudi.pco(d = dist.counts.norm, scannf = FALSE,nf=ncol(counts.norm))
return(pco.counts.norm$li)
return(list(table = pco.counts.norm$li,test=permanova_test))
}
}
......@@ -637,7 +676,57 @@ Get_pca_table <-function(input,dds, group_init, n = min(500, nrow(counts(dds))),
}
Get_nmds_table <-function(input,dds, group_init)
{
## Var of interest
VarInt = input$VarInt
group = as.character(apply(group_init,1,paste, collapse = "-"))
group_init = group
## Keep only some sample
val = c()
for(i in 1:length(VarInt))
{
Tinput = paste("input$","Mod",VarInt[i],sep="")
expr=parse(text=Tinput)
## All the modalities for all the var of interest
val = c(val,eval(expr))
}
if(length(VarInt)>1) Kval = apply(expand.grid(val,val),1,paste, collapse = "-")
else Kval = val
ind_kept = which(as.character(group)%in%Kval)
## Get the group corresponding to the modalities
group = group[ind_kept]
nb = length(unique((group)))
group = as.factor(group)
## To select the colors
indgrp =as.integer(as.factor(group_init))[ind_kept]
if(nlevels(group)!=0)
{
## Get the norm data
counts.norm = as.data.frame(round(counts(dds)))
if(input$CountsType=="Normalized") counts.norm = as.data.frame(round(counts(dds, normalized = TRUE)))
# was removed
counts.norm = counts.norm[,ind_kept]
## Get the distance
if(input$DistClust!="sere") dist.counts.norm = vegdist(t(counts.norm), method = input$DistClust)
if(input$DistClust=="sere") dist.counts.norm = as.dist(SEREcoef(counts.norm))
permanova_test = adonis(dist.counts.norm~group)
## Do NMDS
nmds.counts.norm = metaMDS(dist.counts.norm,k=min(round((nrow(counts.norm)-1)/2-1),round(ncol(counts.norm)/2)), trymax = 1)
proj = nmds.counts.norm$points
return(list(table = proj,test=permanova_test))
}
}
......@@ -646,16 +735,11 @@ Get_pca_table <-function(input,dds, group_init, n = min(500, nrow(counts(dds))),
### NMDS
NMDSPlot <-function (input, dds, group_init, col = c("SpringGreen","dodgerblue","black","firebrick1"))
{
cval=c()
time_set = 0
# Set of shape
shape=c(19,17,15,18)
## Var of interest
VarInt = input$VarInt
## Group
group = as.character(apply(group_init,1,paste, collapse = "-"))
group_init = group
## Keep only some sample
val = c()
......@@ -675,8 +759,13 @@ NMDSPlot <-function (input, dds, group_init, col = c("SpringGreen","dodgerblue",
nb = length(unique((group)))
group = as.factor(group)
## To select the colors
indgrp =as.integer(as.factor(group_init))[ind_kept]
if(nlevels(group)!=0 && !is.null(input$PCaxe1) && !is.null(input$PCaxe2))
{
## Get the norm data
counts.norm = as.data.frame(round(counts(dds)))
if(input$CountsType=="Normalized") counts.norm = as.data.frame(round(counts(dds, normalized = TRUE)))
......@@ -687,71 +776,27 @@ NMDSPlot <-function (input, dds, group_init, col = c("SpringGreen","dodgerblue",
if(input$DistClust!="sere") dist.counts.norm = vegdist(t(counts.norm), method = input$DistClust)
if(input$DistClust=="sere") dist.counts.norm = as.dist(SEREcoef(counts.norm))
save(indgrp,counts.norm,dist.counts.norm,group,group_init,file="testNMDS.RData")
## Do NMDS
nmds.counts.norm = metaMDS(dist.counts.norm,k=round(ncol(counts.norm)/2), trymax = 25)
nmds.counts.norm = metaMDS(dist.counts.norm,k=min(round((nrow(counts.norm)-1)/2-1),round(ncol(counts.norm)/2)), trymax = 25)
dudi.pco(d = dist.counts.norm, scannf = FALSE,nf=ncol(counts.norm))
## Get eigen values
eigen=(pco.counts.norm$eig/sum(pco.counts.norm$eig))*100
proj = nmds.counts.norm$points
## xlim and ylim of the plot
min = min(pco.counts.norm$li)
max = max(pco.counts.norm$li)
## get condition set
condition_set=val[which(val %in% unique(group_init$condition))]
time_set=val[which(val %in% unique(group_init$time))]
## Colors
if(length(col)<length(condition_set) * length(time_set))# && !input$colorgroup)
{
col = rainbow(length(condition_set) * length(time_set))
}
#else if(length(col)<length(condition_set) * length(time_set) && input$colorgroup){
# col = rep(col[1:length(condition_set)], length(time_set))
#}
if (length(time_set) == 1 && length(condition_set) <= 4){
cval = apply(expand.grid(condition_set,time_set),1,paste, collapse = "-")
cval = sort(cval)
}
min = min(proj); max = max(proj)
abs = range(proj[, as.numeric(gsub("PC","",input$PCaxe1))])
abs = abs(abs[2] - abs[1])/25
ord = range(proj[, as.numeric(gsub("PC","",input$PCaxe2))])
ord = abs(ord[2] - ord[1])/25
# to reactivate
#pco.counts.norm$li = pco.counts.norm$li[ind_kept,]
if (plot == "pcoa"){
par(cex=input$cexTitleDiag,mar=c(6,6,4,5))
## Plot axis, label and circles
v_axes = c(as.numeric(gsub("PC","",input$PCaxe1)),as.numeric(gsub("PC","",input$PCaxe2)))
plot(pco.counts.norm$li[v_axes],
xlab=paste(input$PCaxe1, ": ",round(eigen[v_axes[1]],1),"%") ,
ylab=paste(input$PCaxe2, ": ",round(eigen[v_axes[2]],1),"%"),
xlim=c(min+0.25*min,max+0.25*max), ylim=c(min-0.1,max+0.1),
cex.axis=1, cex.lab=1,lwd=2, type="n",main='Principal Coordinates Analysis ')
# Set different shapes
if(input$labelPCOA == "Group"){
if(!is.null(cval)){
for (i in 1:length(cval)){
points(pco.counts.norm$li[which(group==cval[i]),v_axes],pch=shape[i],col=col[i], cex=input$cexpoint)
}
s.class(dfxy = pco.counts.norm$li[v_axes], fac = group, col = col, label = levels(group),
add.plot = TRUE, cpoint = 0, cell=input$cexcircle, clabel=input$cexLabelDiag, cstar = input$cexstar)
}else s.class(dfxy = pco.counts.norm$li[v_axes], fac = group, col = col, label = levels(group),
add.plot = TRUE, cpoint = input$cexpoint, cell=input$cexcircle, clabel=input$cexLabelDiag, cstar = input$cexstar)
}
else{
s.label(pco.counts.norm$li, clabel = input$cexLabelDiag,boxes=FALSE, add.plot = TRUE)
s.class(dfxy = pco.counts.norm$li, fac = group, col = col, label = levels(group), add.plot = TRUE, cpoint = 0, clabel = 0, cstar = input$cexstar, cell=input$cexcircle)
}
}else{
v_axes = c(as.numeric(gsub("PC","",input$PCaxe1)),as.numeric(gsub("PC","",input$PCaxe2)))
nbBar = max(7,max(v_axes))
col = rep("grey",nbBar)
col[v_axes] = "black"
barplot(eigen[1:nbBar], xlab="Dimensions", ylab="Eigenvalues (%)", names.arg = 1:nbBar, col = col, ylim=c(0,max(eigen)+5), cex.axis=1.2, cex.lab=1.4,cex.names=1.2)
}
plot(proj[,as.numeric(gsub("PC","",input$PCaxe1))],proj[,as.numeric(gsub("PC","",input$PCaxe2))], las = 1, cex = input$cexTitleDiag, col = col[indgrp],
pch = 16,
xlab = paste("MDS",gsub("PC","",input$PCaxe1)),
ylab = paste("MDS",gsub("PC","",input$PCaxe2)),
main = "Non-metric multidimensional scaling"
)
abline(h = 0, v = 0, lty = 2, col = "lightgray")
text(proj[,as.numeric(gsub("PC","",input$PCaxe1))] - ifelse(proj[,as.numeric(gsub("PC","",input$PCaxe1))] > 0, abs, -abs), proj[,as.numeric(gsub("PC","",input$PCaxe2))] - ifelse(proj[,as.numeric(gsub("PC","",input$PCaxe2))] > 0, ord, -ord), colnames(counts.norm), col = col[indgrp],cex=input$cexLabelDiag)
}
}
......
server.R
View file @ d3d4f423
......@@ -1451,6 +1451,31 @@ shinyServer(function(input, output,session) {
})
output$ResPermaTestBox <- renderUI({
resDiff = ResDiffAnal()
resTest = Perma_test_Diag(input,resDiff)
resBox = NULL
print(resTest)
if(!is.null(resDiff) && !is.null(resTest))
{
res = list()
## Title
res$title = paste("<center><b><font size='+3'>","Permanova test","</font></b></center><br/>")
## Subtitle
res$subtitle = paste("<center><em>","Analysis of variance using distance matrices","</em></center><br/>")
## Pvalue
res$ccl = paste("<center><b><font size='+1'>p-value :",round(resTest,5),"</font></b></center><br/>")
resBox = box(title="Permanova ",width = 6, status = "primary", solidHeader = TRUE,collapsible = TRUE,collapsed = FALSE,
HTML(unlist(res))
)
}
return(resBox)
})
output$SizeFactTable <- DT::renderDataTable(
SizeFactor_table(),
options = list(scrollX=TRUE,searching = FALSE, processing=FALSE
......
ui.R
View file @ d3d4f423
......@@ -445,24 +445,29 @@ body <- dashboardBody(
tags$style(type='text/css', "#ExportSizeFactor { margin-top: 37px;}")
)
),
fluidRow(
conditionalPanel(condition="input.DiagPlot=='pcaPlot'",
box(title = "Eigen values", width = 6, status = "primary", solidHeader = TRUE,collapsible = TRUE,collapsed= FALSE,
plotOutput("PlotEigen",height="100%")
)
),
conditionalPanel(condition="input.DiagPlot=='pcoaPlot'",
box(title = "Eigen values", width = 6, status = "primary", solidHeader = TRUE,collapsible = TRUE,collapsed= FALSE,
plotOutput("PlotpcoaEigen",height="100%")
)
),
conditionalPanel(condition="input.DiagPlot=='pcoaPlot' || input.DiagPlot=='nmdsPlot'",
uiOutput("ResPermaTestBox")
)
)
),
column(width=3,
box(title = "Select your plot", width = NULL, status = "primary", solidHeader = TRUE,collapsible = FALSE,collapsed= FALSE,
selectInput("DiagPlot","",c("Total barplot"="barplotTot","Nul barplot"="barplotNul",
"Maj. taxonomy"="MajTax","Boxplots" = "boxplotNorm", "Density"="densityPlot", "Dispersion" = "DispPlot",
"Size factors VS total"="SfactorsVStot", "PCA"="pcaPlot", "PCoA"="pcoaPlot","Clustering" = "clustPlot"))
"Size factors VS total"="SfactorsVStot", "PCA"="pcaPlot", "PCoA"="pcoaPlot","NMDS"="nmdsPlot","Clustering" = "clustPlot"))
),
box(title = "Options", width = NULL, status = "primary", solidHeader = TRUE,collapsible = TRUE,collapsed= FALSE,
conditionalPanel(condition="input.DiagPlot!='clustPlot' && input.DiagPlot!='pcaPlot' && input.DiagPlot!='SfactorsVStot' && input.DiagPlot!='DispPlot'",
......@@ -472,7 +477,7 @@ body <- dashboardBody(
checkboxInput("RemoveNullValue","Remove 0",value = TRUE)
),
conditionalPanel(condition="input.DiagPlot!='Sfactors' && input.DiagPlot!='SfactorsVStot' ",uiOutput("VarIntDiag")),
conditionalPanel(condition="input.DiagPlot=='pcoaPlot' || input.DiagPlot=='pcaPlot'",
conditionalPanel(condition="input.DiagPlot=='pcoaPlot' || input.DiagPlot=='pcaPlot' || input.DiagPlot=='nmdsPlot'",
h5(strong("Select the modalities")),
uiOutput("ModMat"),
fluidRow(
......@@ -481,7 +486,7 @@ body <- dashboardBody(
column(width=5,uiOutput("PC2_sel"))
)
),
conditionalPanel(condition="input.DiagPlot=='pcoaPlot' || input.DiagPlot=='SERE' || input.DiagPlot=='clustPlot' ",
conditionalPanel(condition="input.DiagPlot=='pcoaPlot' || input.DiagPlot=='nmdsPlot' || input.DiagPlot=='SERE' || input.DiagPlot=='clustPlot' ",
selectInput("DistClust","Distance",c("euclidean", "SERE"="sere", "canberra", "bray", "kulczynski", "jaccard",
"gower", "altGower", "morisita", "horn","mountford","raup","binomial",
"chao","cao","mahalanobis"),selected="canberra")
......@@ -523,7 +528,7 @@ body <- dashboardBody(
fluidRow(
column(width=12, p(strong("Size"))),
column(width=6,sliderInput("cexTitleDiag", h6("Axis"),min=0,max=5,value = 1,step =0.1)),
conditionalPanel(condition="input.DiagPlot=='SfactorsVStot' || input.DiagPlot=='pcaPlot' || input.DiagPlot=='pcoaPlot'",column(width=6,sliderInput("cexLabelDiag", h6("Points"),min=0,max=5,value = 1,step =0.1)))
conditionalPanel(condition="input.DiagPlot=='SfactorsVStot' || input.DiagPlot=='pcaPlot' || input.DiagPlot=='pcoaPlot' || input.DiagPlot=='nmdsPlot' ",column(width=6,sliderInput("cexLabelDiag", h6("Points"),min=0,max=5,value = 1,step =0.1)))
)
......
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