Commit de934b49 authored by svolant's avatar svolant
Browse files

Check target + modif filtering + comment tree

parent 0ad799ec
......@@ -89,6 +89,73 @@ CheckTaxoTable <- function(taxo,counts)
return(list(Error=Error,Warning=Warning))
}
CheckTargetModel <- function(input,target,labeled,CT)
{
Error = NULL
HowTo = NULL
InterVar = input$InterestVar
labels = rownames(target)
ind = which(colnames(CT)%in%labels)
## At least one variable selected
if(is.null(Error) && length(ind)<=1){
Error = "Less than two samples names fit with the counts table"
HowTo = "Check the samples names in the target file. They must be in the first column and must correspond EXACTLY to the names in the count table."
}
## At least one variable selected
if(is.null(Error) && length(InterVar)==0){
Error = "At least one variable must be selected for the model"
HowTo = "Add at least one variable in the 'Select the variables' widget"
}
## Names of samples correct ?
if(is.null(Error) && labeled==0){
Error = "The names of the samples in the target file do not fit the counts table"
HowTo = "Check the samples names in the target file. They must be in the first column and must correspond EXACTLY to the names in the count table."
}
## Number of columns
if(is.null(Error) && ncol(target)<2){
Error = "The number of columns of the target file must be at least 2"
HowTo = "Add at least one additional variable to describe your samples"
}
## Number of rows
if(is.null(Error) && nrow(target)<=1){
Error = "The number of rows if the target file must be at least 2"
HowTo = "Add information about more than 2 samples"
}
## NA values
if(is.null(Error) && (any(is.na(target)) || any(target ==""))){
Error = "NA's or missing values are not allowed in the target file"
HowTo = "Remove all the samples for which one or more variables are NA or missing"
}
## Full rank matrix
if(is.null(Error) && length(InterVar)>0)
{
design = GetDesign(input)
testRank = CheckMatrixRank(design,target)
if(!testRank){
Error = "The model matrix is not full rank. One or more variables or interaction terms
are linear combinations of the others and must be removed."
HowTo = "Remove variable(s) that provide the same information, i.e, if the value of a variable is totaly determine by an other variable remove one of them."
}
}
return(list(Error=Error,HowTo=HowTo))
}
## Get the percentage of annotated OTU
PercentAnnot <- function(counts,taxo)
{
......@@ -247,12 +314,10 @@ GetCountsMerge <- function(input,dataInput,taxoSelect,target,design)
# Only interesting OTU
# merged_table = merge(CT, taxo[order(rownames(CT)),], by="row.names")
save(CT,taxo,file="testMerge.RData")
merged_table = merge(CT, taxo, by="row.names")
## Why ??
CT = merged_table[,2: (dim(CT)[2]+1)]
taxo = merged_table[,(dim(CT)[2]+2):dim(merged_table)[2]]
## ???
rownames(CT) = merged_table[,1]
rownames(taxo) = merged_table[,1]
......@@ -407,7 +472,7 @@ plot_filter <- function(counts,th.samp,th.abund,type="Scatter")
## Initial plot for plotly
if(type == 'Abundance' || type == 'Samples'){
dataNull = data.frame(x=c(0,0),y=c(1,2),col=c("white","white"))
res = ggplotly(ggplot(dataNull,aes(x=x,y=y))+geom_point(aes(colour = col))+theme_bw()+ scale_color_manual(values = "white"))
res = ggplot(dataNull,aes(x=x,y=y))+geom_point(aes(colour = col))+theme_bw()+ scale_color_manual(values = "white")
}
res_filter = Filtered_feature(counts,th.samp,th.abund)
......@@ -430,12 +495,12 @@ plot_filter <- function(counts,th.samp,th.abund,type="Scatter")
## plot the results
gg = ggplot(df,aes(lab,y,fill=State)) + geom_bar(stat='identity') + theme_bw() +theme(axis.text.x = element_text(angle = 90, hjust = 1,vjust = 0.5))
gg = gg + geom_hline( yintercept = th.abund,linetype = "longdash")
gg = gg + geom_hline( yintercept = th.abund,linetype = "longdash") + xlab("")
gg = gg + scale_fill_manual(values = c('springgreen3','firebrick'))
if(!"Kept"%in%df$State ) gg = gg + scale_fill_manual(values = 'firebrick')
if(!"Removed"%in%df$State ) gg = gg + scale_fill_manual(values = 'springgreen3')
res = ggplotly(gg)
res = gg
}
if(type == 'Samples' && !is.null(th.samp) && !is.null(th.abund))
......@@ -455,12 +520,12 @@ plot_filter <- function(counts,th.samp,th.abund,type="Scatter")
## plot the results
gg = ggplot(df,aes(lab,y,fill=State)) + geom_bar(stat='identity') + theme_bw() +theme(axis.text.x = element_text(angle = 90, hjust = 1,vjust = 0.5))
gg = gg + geom_hline( yintercept = th.samp,linetype = "longdash")
gg = gg + geom_hline( yintercept = th.samp,linetype = "longdash") + xlab("")
gg = gg + scale_fill_manual(values = c('springgreen3','firebrick'))
if(!"Kept"%in%df$State ) gg = gg + scale_fill_manual(values = 'firebrick')
if(!"Removed"%in%df$State ) gg = gg + scale_fill_manual(values = 'springgreen3')
res = ggplotly(gg)
res = gg
}
......@@ -488,9 +553,7 @@ plot_filter <- function(counts,th.samp,th.abund,type="Scatter")
State = df$State
PointSize = 2
colors_scat = list(Kept="#00CD66",Removed="#b22222")
print(State)
print(colors_scat)
print(head(df))
res = scatterD3(x = x_var,
y = y_var,
lab = rownames(df),
......@@ -532,15 +595,17 @@ plot_filter <- function(counts,th.samp,th.abund,type="Scatter")
SelectThreshAb <- function(infile,lambda=500,graph=TRUE){
rs <- rowSums(infile)
lambda = max(rs)
test_Filtre <- sapply(c(min(rs):lambda),FUN=function(x) table(rs>x))
x <- c(min(rs):lambda)
reslm <- lm(test_Filtre[1,]~x)$coefficients
val <- which(test_Filtre[1,]>reslm[1])[1]
if (graph==TRUE){
plot(test_Filtre[1,],ylab="Nb especes avec moins de x comptages sur tous les echantillons",xlab="x")
abline(v=val,col="red")
}
return(val)
if(!is.list(test_Filtre))
{
x <- c(min(rs):lambda)
reslm <- lm(test_Filtre[1,]~x)$coefficients
val <- which(test_Filtre[1,]>reslm[1])[1]
if (graph==TRUE){
plot(test_Filtre[1,],ylab="Nb especes avec moins de x comptages sur tous les echantillons",xlab="x")
abline(v=val,col="red")
}
} else val = min(rs)
return(max(val, min(rs)+1))
}
......@@ -12,20 +12,22 @@
## Get the possible interaction for the statistical model
GetInteraction2 <- function(target)
GetInteraction2 <- function(target,VarInt)
{
res=c()
namesTarget = colnames(target)[2:ncol(target)]
k=1
for(i in 1:(length(namesTarget)-1))
{
for(j in (i+1):length(namesTarget))
if(length(VarInt)>=2)
{
k=1
for(i in 1:(length(VarInt)-1))
{
res[k] = paste(namesTarget[i],":",namesTarget[j],sep="")
k = k+1
for(j in (i+1):length(VarInt))
{
res[k] = paste(VarInt[i],":",VarInt[j],sep="")
k = k+1
}
}
}
return(res)
}
......@@ -33,7 +35,6 @@ GetInteraction2 <- function(target)
## Get the dds object of DESeq2
Get_dds_object <- function(input,counts,target,design,normFactorsOTU,CT_noNorm,CT_Norm)
{
dds <- DESeqDataSetFromMatrix(countData=counts, colData=target, design=design,ignoreRank = TRUE)
sizeFactors(dds) = normFactorsOTU
......@@ -53,10 +54,13 @@ Get_dds_object <- function(input,counts,target,design,normFactorsOTU,CT_noNorm,C
## Get the design according to the input
GetDesign <- function(input)
{
design = NULL
Interaction = NULL
InterVar = input$InterestVar
Interaction = input$Interaction2
if(length(InterVar)>1) Interaction = input$Interaction2
alltmp = c(InterVar,Interaction)
design = as.formula(paste("~", paste0(alltmp, collapse= "+")))
if(length(alltmp)>0) design = as.formula(paste("~", paste0(alltmp, collapse= "+")))
return(design)
}
......@@ -167,6 +171,23 @@ PrintContrasts <- function (coefs, contrasts,contnames)
}
## Check the rank of the matrix for the statistical model.
CheckMatrixRank <- function(design,target)
{
testRank = FALSE
if(!is.null(design) && length(as.character(design)))
{
tmp = strsplit(as.character(design)[2],"[ ]*[+][ ]*")
testDesign = any(unlist(tmp)%in%colnames(target))
}
if(nrow(target)>=1 && ncol(target)>=2 && nrow(target)>1 && !is.null(design) && testDesign){
modelMatrix <- model.matrix(design, data = as.data.frame(target))
rk = qr(modelMatrix)$rank
testRank = (ncol(modelMatrix)==rk)
}
return(testRank)
}
......@@ -692,7 +692,7 @@ Plot_Visu_Tree <- function(input,resDiff,CT_Norm_OTU,taxo_table)
{
tmp = CreateTableTree(input,resDiff,CT_Norm_OTU,taxo_table,VarInt)
if(nrow(tmp$counts)>0 && !is.null(tmp$counts))
if(nrow(tmp$counts)>0 && !is.null(tmp$counts) && !is.null(input$TaxoTree))
{
merge_dat = merge(taxo_table,round(t(tmp$counts)),by="row.names")
......
source('LoadPackages.R')
library(plotly)
library(treeWeightD3)
if(!require(plotly)){
install.packages("plotly")
library(plotly)
}
shinyServer(function(input, output,session) {
hide(id = "loading-content", anim = TRUE, animType = "fade",time=1.5)
......@@ -175,33 +177,39 @@ shinyServer(function(input, output,session) {
normFactors = NULL
CT_noNorm = NULL
CT_Norm = NULL
#labeled= NULL
ChTM = NULL
data = isolate(dataInput()$data)
target = isolate(dataInputTarget()$target)
labeled= isolate(dataInputTarget()$labeled)
taxo = isolate(input$TaxoSelect)
withProgress(
if(!is.null(data$counts) && !is.null(data$taxo) && nrow(data$counts)>0 && nrow(data$taxo)>0 && !is.null(taxo) && taxo!="..." && !is.null(target))
{
design = GetDesign(isolate(input))
tmp = isolate(GetCountsMerge(input,data,taxo,target,design))
counts = tmp$counts
ChTM = CheckTargetModel(input,target,labeled,data$counts)$Error
## Filtering the counts
if(isolate(input$AddFilter) && !is.null(isolate(input$SliderThSamp)) && !is.null(isolate(input$SliderThAb)))
if(!is.null(design) && is.null(ChTM))
{
ind.filter =Filtered_feature(counts,isolate(input$SliderThSamp),isolate(input$SliderThAb))$ind
counts = counts[-ind.filter,]
tmp = isolate(GetCountsMerge(input,data,taxo,target,design))
counts = tmp$counts
## Filtering the counts
if(isolate(input$AddFilter) && !is.null(isolate(input$SliderThSamp)) && !is.null(isolate(input$SliderThAb)))
{
ind.filter =Filtered_feature(counts,isolate(input$SliderThSamp),isolate(input$SliderThAb))$ind
counts = counts[-ind.filter,]
}
CheckTarget = tmp$CheckTarget
#target = tmp$target
#labeled = tmp$labeled
normFactors = tmp$normFactors
## OTU table, norm and no norm
CT_noNorm = tmp$CT_noNorm
CT_Norm = tmp$CT_Norm
}
CheckTarget = tmp$CheckTarget
#target = tmp$target
#labeled = tmp$labeled
normFactors = tmp$normFactors
## OTU table, norm and no norm
CT_noNorm = tmp$CT_noNorm
CT_Norm = tmp$CT_Norm
}
,message="Merging the counts ...")
return(list(counts=counts,CheckTarget=CheckTarget,normFactors=normFactors,CT_noNorm=CT_noNorm, CT_Norm=CT_Norm))
......@@ -277,7 +285,8 @@ shinyServer(function(input, output,session) {
res = NULL
counts = isolate(dataMergeCounts()$counts)
tot = rowSums(counts)
tmp = SelectThreshAb(counts,lambda=max(round(sum(counts)/nrow(counts)*0.05),min(tot)+1),graph=FALSE)
save(counts,tot,file="testFilter.RData")
withProgress({tmp = SelectThreshAb(counts,lambda=max(round(sum(counts)/nrow(counts)*0.05),min(tot)+1),graph=FALSE)},message="Loading...")
res = sliderInput("SliderThAb","Threshold on the total abundance (in log)",min=0,max=round(max(log(tot+1)),1),value = log(tmp+1))
return(res)
......@@ -305,14 +314,14 @@ shinyServer(function(input, output,session) {
# plot_filter(counts,th.samp,th.abund,type="Scatter")
output$Plot_ThAb <- renderPlotly({
output$Plot_ThAb <- renderPlot({
counts = dataMergeCounts()$counts
## output of plot_filter is ggplot class
plot_filter(counts,input$SliderThSamp,input$SliderThAb,type="Abundance")
})
output$Plot_ThSamp <- renderPlotly({
output$Plot_ThSamp <- renderPlot({
counts = dataMergeCounts()$counts
## output of plot_filter is ggplot class
plot_filter(counts,input$SliderThSamp,input$SliderThAb,type="Samples")
......@@ -480,17 +489,17 @@ shinyServer(function(input, output,session) {
ind = which(rownames(data)%in%colnames(counts))
data = as.data.frame(data[ind,])
colnames(data) = names
## Replace "-" by "."
ind_num = which(sapply(as.data.frame(data[,-1]),is.numeric)) + 1
if(length(ind_num)>0){
data_tmp =cbind( as.data.frame(apply(as.data.frame(data[,-ind_num]),2,gsub,pattern = "-",replacement = ".")),data[,ind_num])
colnames(data_tmp) = c(colnames(data)[-ind_num],colnames(data)[ind_num])
data = data_tmp
if(ncol(data)>1 && nrow(data)>1){
ind_num = which(sapply(as.data.frame(data[,-1]),is.numeric)) + 1
if(length(ind_num)>0){
data_tmp =cbind( as.data.frame(apply(as.data.frame(data[,-ind_num]),2,gsub,pattern = "-",replacement = ".")),data[,ind_num])
colnames(data_tmp) = c(colnames(data)[-ind_num],colnames(data)[ind_num])
data = data_tmp
}
if(length(ind_num)==0){data = as.data.frame(apply(data,2,gsub,pattern = "-",replacement = "."))}
}
if(length(ind_num)==0){data = as.data.frame(apply(data,2,gsub,pattern = "-",replacement = "."))}
target = data
target = as.data.frame(data)
# target = as.data.frame(apply(target,2,gsub,pattern = "-",replacement = "."))
......@@ -522,18 +531,20 @@ shinyServer(function(input, output,session) {
})
## Interactions
output$SelectInteraction2 <- renderUI({
target = dataInputTarget()$target
if(!is.null(target))
VarInt = input$InterestVar
res = NULL
if(!is.null(target) && length(input$InterestVar)>1)
{
Interac = GetInteraction2(target)
selectInput("Interaction2",h6(strong("Add interactions")),Interac,selected=NULL,multiple=TRUE)
Interac = GetInteraction2(target,VarInt)
res = selectInput("Interaction2",h6(strong("Add interactions")),Interac,selected=NULL,multiple=TRUE)
}
if(length(input$InterestVar)==1) res = NULL
return(res)
})
......@@ -1134,6 +1145,65 @@ shinyServer(function(input, output,session) {
})
# Infobox model/target
# output$InfoModel<- renderInfoBox({
# res = infoBox(h6(strong("Target file and variables")),
# subtitle = h6(strong("Your target file must contain at least 2 columns and 2 rows. NA's values are not allowed and the variables must not be collinear.")),
# icon = icon("book"),color = "green",width=NULL,fill=TRUE)
#
# target = dataInputTarget()$target
# taxo = input$TaxoSelect
# ChTM = NULL
#
# ## Return NULL if there is no error
# if(!is.null(target)) ChTM = CheckTargetModel(input,target)
#
# if(!is.null(ChTM)) res = infoBox(h6(strong("Error")), subtitle = h6(ChTM), icon = icon("thumbs-o-down"),color = "red",width=NULL,fill=TRUE)
#
# return(res)
#
# })
output$InfoModel<- renderUI({
CT = dataInput()$data$counts
target = dataInputTarget()$target
labeled = dataInputTarget()$labeled
taxo = input$TaxoSelect
ChTM = NULL
## Return NULL if there is no error
if(!is.null(target)) ChTM = CheckTargetModel(input,target,labeled,CT)
if(!is.null(ChTM$Error)) {
box(title = "Error", status = "danger",width = 6,
h6(strong(ChTM$Error)),
footer = em("Reminder: Your target file must contain at least 2 columns and 2 rows. NA's values are not allowed and the variables must not be collinear.")
)
}
})
output$InfoModelHowTo<- renderUI({
CT = dataInput()$data$counts
target = dataInputTarget()$target
labeled = dataInputTarget()$labeled
taxo = input$TaxoSelect
ChTM = NULL
## Return NULL if there is no error
if(!is.null(target)) ChTM = CheckTargetModel(input,target,labeled,CT)
if(!is.null(ChTM$HowTo)) {
box(title = "How To", status = "success",width = 6,
h6(strong(ChTM$HowTo))
)
}
})
#####################################################
##
......@@ -1503,9 +1573,18 @@ shinyServer(function(input, output,session) {
output$RunButton <- renderUI({
res = NULL
ChTM = "Error"
target = dataInputTarget()$target
labeled = dataInputTarget()$labeled
CT = dataInput()$data$counts
taxo = input$TaxoSelect
if(!is.null(target) && taxo!="...") res = actionButton("RunDESeq",strong("Run analysis"),icon = icon("caret-right"))
VarInt = input$InterestVar
## Return NULL if there is no error
if(!is.null(target) && length(VarInt)>=1) ChTM = CheckTargetModel(input,target,labeled,CT)$Error
if(!is.null(target) && taxo!="..." && is.null(ChTM) && length(VarInt)>=1) res = actionButton("RunDESeq",strong("Run analysis"),icon = icon("caret-right"))
return(res)
})
......@@ -2017,36 +2096,29 @@ shinyServer(function(input, output,session) {
observe({
input$TaxoSelect
testRank = FALSE
counts = dataMergeCounts()$counts
InterVar = input$InterestVar
if(length(InterVar)>0)
{
design = GetDesign(input)
target = dataInputTarget()$target
modelMatrix <- model.matrix(design, data = as.data.frame(target))
rk = qr(modelMatrix)$rank
testRank = (ncol(modelMatrix)==rk)
}
ChTM = ""
CT = dataInput()$data$counts
target = dataInputTarget()$target
labeled= dataInputTarget()$labeled
ChTM = CheckTargetModel(input,target,labeled,CT)$Error
if(input$AddFilter && !is.null(input$SliderThSamp) && !is.null(input$SliderThAb))
if(input$AddFilter && !is.null(input$SliderThSamp) && !is.null(input$SliderThAb) && is.null(ChTM))
{
ind.filter =Filtered_feature(counts,input$SliderThSamp,input$SliderThAb)$ind
counts = counts[-ind.filter,]
}
if (!is.null(counts)){
if (nrow(counts)>=2){
shinyjs::enable("RunDESeq")
}
if (nrow(counts)<2 || !testRank) {
if (nrow(counts)<2) {
shinyjs::disable("RunDESeq")
}
}
if (is.null(counts) || !testRank) {
if (is.null(counts) ) {
shinyjs::disable("RunDESeq")
}
......
......@@ -6,7 +6,7 @@ if (!require(devtools)) {
library(devtools)
}
if(!require(treeWeightD3)){
devtools::install_git('https://gitlab.pasteur.fr/plechat/treeWeightD3')
devtools::install_git('https://gitlab.pasteur.fr/plechat/treeWeightD3.git')
library(treeWeightD3)
}
if (!require(scatterD3)) {
......@@ -25,7 +25,11 @@ if(!require(shinyjs)){
install.packages("shinyjs")
library(shinyjs)
}
library(plotly)
if(!require(plotly)){
install.packages("plotly")
library(plotly)
}
function(request) {
sidebar <- dashboardSidebar(
useShinyjs(),
......@@ -309,10 +313,13 @@ body <- dashboardBody(
),
fluidRow(
column(width=6,uiOutput("SelectInterestVar")),
column(width=6,uiOutput("SelectInteraction2")),
column(width=6,uiOutput("SelectInteraction2"))
),
fluidRow(
column(width=6,uiOutput("RunButton"))
)
),
fluidRow(uiOutput("InfoModel"),uiOutput("InfoModelHowTo")),
uiOutput("BoxTarget"),
uiOutput("BoxCountsMerge")
),
......@@ -372,11 +379,11 @@ body <- dashboardBody(
fluidRow(
conditionalPanel(condition="input.AddFilter",
column(width=6,
plotlyOutput("Plot_ThSamp"),
plotOutput("Plot_ThSamp"),
column(width=10,uiOutput("ThSamp"),offset = 1)
),
column(width=6,
plotlyOutput("Plot_ThAb"),
plotOutput("Plot_ThAb"),
column(width=10,uiOutput("ThAb"),offset = 1)
),
column(width=12,
......@@ -589,7 +596,8 @@ body <- dashboardBody(
column(width=3,
box(title = "Select your plot", width = NULL, status = "primary", solidHeader = TRUE,collapsible = FALSE,collapsed= FALSE,
selectizeInput("PlotVisuSelect","",c("Barplot"="Barplot","Heatmap"="Heatmap","Boxplot"="Boxplot","Tree"="Tree","Scatterplot"="Scatterplot","Diversity"="Diversity","Rarefaction"="Rarefaction"),selected = "Barplot")
# selectizeInput("PlotVisuSelect","",c("Barplot"="Barplot","Heatmap"="Heatmap","Boxplot"="Boxplot","Tree"="Tree","Scatterplot"="Scatterplot","Diversity"="Diversity","Rarefaction"="Rarefaction"),selected = "Tree")
selectizeInput("PlotVisuSelect","",c("Barplot"="Barplot","Heatmap"="Heatmap","Boxplot"="Boxplot","Scatterplot"="Scatterplot","Diversity"="Diversity","Rarefaction"="Rarefaction"),selected = "Barplot")
),
......
Supports Markdown
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