Skip to content

GitLab

Commit a9956d3a authored by svolant's avatar svolant
Browse files

bug biom+tree abundance

parent 23fa289e
Hide whitespace changes
Inline Side-by-side
Rfunctions/Data_Management.R
View file @ a9956d3a
......@@ -47,8 +47,6 @@ read_rdp <- function(filename, threshold_annot)
## Check the format of the counts table
CheckCountsTable <- function(counts)
{
......@@ -165,13 +163,15 @@ GetCountsMerge <- function(input,dataInput,taxoSelect,target,design)
CT_noNorm = NULL
normFactors = NULL
FeatureSize = NULL
CT_Norm = NULL
## Counts and taxo tables
CT = dataInput$counts
taxo = dataInput$taxo
# save(CT,target,taxo,file="testMerge.RData")
## Select cols in the target
labels = target[,1]
labels = rownames(target)
ind = which(colnames(CT)%in%labels)
## Get the normalization variable (normalization can be done according to this variable)
......@@ -200,10 +200,10 @@ GetCountsMerge <- function(input,dataInput,taxoSelect,target,design)
CT = OrderCounts(counts=CT,labels=labels)$CountsOrder
CT_noNorm = CT
RowProd = sum(apply(CT_noNorm,1,prod))
## Create the dds object
dds <- DESeqDataSetFromMatrix(countData=CT, colData=target, design=design)
dds <- DESeqDataSetFromMatrix(countData=CT, colData=target, design=design,ignoreRank=TRUE)
save(dds,file="testdds.RData")
if(is.null(VarNorm)){
## Counts normalisation
## Normalisation with or without 0
......@@ -225,7 +225,7 @@ GetCountsMerge <- function(input,dataInput,taxoSelect,target,design)
CT_tmp = CT[,indgrp]
CT_tmp = removeNulCounts(CT_tmp)
target_tmp = data.frame(labels = rownames(target)[indgrp])
dds_tmp <- DESeqDataSetFromMatrix(countData=CT_tmp, colData=target_tmp, design=~labels)
dds_tmp <- DESeqDataSetFromMatrix(countData=CT_tmp, colData=target_tmp, design=~labels,ignoreRank=TRUE)
if(input$AccountForNA=="NonNull") {dds_tmp = estimateSizeFactors(dds_tmp,locfunc=eval(as.name(input$locfunc)),geoMeans=GeoMeansCT(CT_tmp)); normFactors[indgrp] = sizeFactors(dds_tmp)}
if(input$AccountForNA=="All") {dds_tmp = estimateSizeFactors(dds_tmp,locfunc=eval(as.name(input$locfunc))); normFactors[indgrp] = sizeFactors(dds_tmp)}
if(input$AccountForNA=="Weighted" && input$AccountForNA!="NonNull" ) {dds_tmp = estimateSizeFactors(dds_tmp,locfunc=eval(as.name(input$locfunc)),geoMeans=GeoMeansCT(CT_tmp)); normFactors[indgrp] = w.sizefactor(CT_tmp)}
......@@ -247,9 +247,13 @@ 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]
#ordOTU = order(rownames(taxo))
......@@ -258,9 +262,9 @@ GetCountsMerge <- function(input,dataInput,taxoSelect,target,design)
# indOTU_annot = which(rownames(CT)%in%rownames(taxo))
# counts_annot = CT[indOTU_annot[ordOTU],]
## Aggregate matrix
if(taxoSelect=="OTU/Gene" || input$FileFormat=="fileBiom") counts = counts_annot
if(taxoSelect=="OTU/Gene") counts = counts_annot
else{
if(input$TypeTable == "MGS"){
if(input$TypeTable == "MGS" && input$FileFormat!="fileBiom"){
taxoS = taxo[,input$TypeTable]
counts = aggregate(counts_annot,by=list(Taxonomy = taxoS),mean)
rownames(counts)=counts[,1]
......@@ -270,7 +274,7 @@ GetCountsMerge <- function(input,dataInput,taxoSelect,target,design)
colnames(counts_int)=colnames(counts)
counts=counts_int
}
if(taxoSelect != "MGS"){
if(taxoSelect != "MGS" || input$FileFormat=="fileBiom"){
#taxoS = taxo[ordOTU,taxoSelect]
taxoS = taxo[,taxoSelect]
counts = aggregate(counts_annot,by=list(Taxonomy = taxoS),sum)
......@@ -528,6 +532,7 @@ 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
......
Rfunctions/Stat_Model.R
View file @ a9956d3a
......@@ -33,7 +33,9 @@ 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)
dds <- DESeqDataSetFromMatrix(countData=counts, colData=target, design=design,ignoreRank = TRUE)
sizeFactors(dds) = normFactorsOTU
dds <- estimateDispersions(dds, fitType=input$fitType)
if(as.numeric(R.Version()$major)>=3 && as.numeric(R.Version()$minor) >=1.3){
......
Rfunctions/VisuPlot.R
View file @ a9956d3a
......@@ -529,7 +529,7 @@ GetDataToPlot <- function(input,resDiff,VarInt,ind_taxo,aggregate=TRUE,rarefy=FA
## Create the variable to plot
targetInt = as.data.frame(target[,VarInt])
rownames(targetInt)=target[,1]
rownames(targetInt)=rownames(target)
## Combining the Varint
if(length(VarInt)>1){targetInt$AllVar = apply(targetInt,1,paste, collapse = "-"); targetInt$AllVar = factor(targetInt$AllVar,levels = expand.grid2.list(list.val))}
if(length(VarInt)<=1){targetInt$AllVar = target[,VarInt]; targetInt$AllVar = factor(targetInt$AllVar,levels = val)}
......@@ -586,4 +586,126 @@ GetDataToPlot <- function(input,resDiff,VarInt,ind_taxo,aggregate=TRUE,rarefy=FA
}
## Create the data table for the tree representation
CreateTableTree <- function(input,resDiff,CT_Norm_OTU,taxo_table,VarInt,ind_taxo=rownames(CT_Norm_OTU))
{
dds = resDiff$dds
val = c()
list.val = list()
counts = CT_Norm_OTU
target = resDiff$target
counts_tmp_combined = NULL
prop_tmp_combined = NULL
targetInt = NULL
namesCounts = NULL
levelsMod = NULL
## Select a subset within the taxonomy level (default is the 12 most abundant)
nbKept = length(ind_taxo)
Taxonomy = rownames(counts)
if (length(VarInt)>0 && nbKept>0)
{
## Get the modalities to keep
for(i in 1:length(VarInt))
{
## Replace "-" by "."
target[,VarInt[i]] = gsub("-",".",target[,VarInt[i]])
Tinput = paste("input$","ModVisu",VarInt[i],sep="")
expr=parse(text=Tinput)
## All the modalities for all the var of interest
val = c(val,eval(expr))
list.val[[i]] = eval(expr)
}
if (!is.null(val) && !is.null(list.val))
{
## Create the variable to plot
targetInt = as.data.frame(target[,VarInt])
rownames(targetInt)=rownames(target)
## Combining the Varint
if(length(VarInt)>1){targetInt$AllVar = apply(targetInt,1,paste, collapse = "-"); targetInt$AllVar = factor(targetInt$AllVar,levels = expand.grid2.list(list.val))}
if(length(VarInt)<=1){targetInt$AllVar = target[,VarInt]; targetInt$AllVar = factor(targetInt$AllVar,levels = val)}
colnames(targetInt) = c(VarInt,"AllVar")
## Keep only the selected modalities
ind_kept = which(!is.na(targetInt$AllVar))
targetInt = targetInt[ind_kept,]
levelsMod = levels(targetInt$AllVar)
## Create the counts matrix only for the selected subset
counts_tmp = counts[Taxonomy%in%ind_taxo,]
counts_tmp = counts_tmp[,colnames(counts_tmp)%in%rownames(targetInt)]
## Be careful transposition !
if(nrow(counts_tmp)>0 && nrow(targetInt)>0)
{
counts_tmp_combined = aggregate(t(counts_tmp),by=list(targetInt$AllVar),mean)
rownames(counts_tmp_combined) = counts_tmp_combined$Group.1
namesCounts = counts_tmp_combined$Group.1
counts_tmp_combined = as.matrix(counts_tmp_combined[,-1])
}
## Ordering the counts
if(!is.null(counts_tmp_combined))
{
MeanCounts = apply(counts_tmp_combined,2,mean)
ord = order(MeanCounts,decreasing=TRUE)
counts_tmp_combined = as.matrix(counts_tmp_combined[,ord])
}
}
}
return(list(counts = counts_tmp_combined,targetInt=targetInt,namesCounts=namesCounts,levelsMod=levelsMod))
}
###########################
## Tree
###########################
## The count matrix must be given at the leaf level.
Plot_Visu_Tree <- function(input,resDiff,CT_Norm_OTU,taxo_table)
{
res = NULL
## Get Input for BarPlot
VarInt = input$VisuVarInt
ind_taxo = input$selectTaxoPlot
## Removed column with only 1 modality
ind = which(apply(taxo_table,2,FUN = function(x) length(unique(x[!is.na(x)])))==1)
if(length(ind)>0) taxo_table = taxo_table[,-ind]
# tmp_combined = GetDataToPlot(input,resDiff,VarInt,ind_taxo,CT_Norm_OTU=CT_Norm_OTU)
if(nrow(CT_Norm_OTU)>0 && !is.null(CT_Norm_OTU) && nrow(taxo_table)>0 && !is.null(taxo_table))
{
tmp = CreateTableTree(input,resDiff,CT_Norm_OTU,taxo_table,VarInt)
if(nrow(tmp$counts)>0 && !is.null(tmp$counts))
{
merge_dat = merge(taxo_table,round(t(tmp$counts)),by="row.names")
colnames(merge_dat)[1] = "OTU"
levels <- c("OTU", colnames(taxo_table))
conditions <- rownames(tmp$counts)
nodeFind = input$TaxoTree
if(input$TaxoTree=="...") nodeFind = NULL
res = treeWeightD3(merge_dat,conditions,levels,nodeFind=nodeFind, height =input$heightVisu+10)
}
}
return(res)
}
server.R
View file @ a9956d3a
source('LoadPackages.R')
library(plotly)
library(treeWeightD3)
shinyServer(function(input, output,session) {
hide(id = "loading-content", anim = TRUE, animType = "fade",time=1.5)
......@@ -471,18 +471,26 @@ shinyServer(function(input, output,session) {
data = read.csv(inFile$datapath,sep=input$septarget,header=TRUE)
data = as.data.frame(data)
names = colnames(data)
## Change the rownames
rownames(data) <- as.character(data[, 1])
## Keep only the row which are in the count table
ind = which(rownames(data)%in%colnames(counts))
data = as.data.frame(data[ind,])
colnames(data) = names
## Replace "-" by "."
ind_num = which(sapply(data,is.numeric))
ind_num = which(sapply(as.data.frame(data[,-1]),is.numeric)) + 1
if(length(ind_num)>0){
data_tmp =cbind( as.data.frame(apply(data[,-ind_num],2,gsub,pattern = "-",replacement = ".")),data[,ind_num])
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 = "."))}
rownames(data) <- as.character(data[, 1])
ind = which(rownames(data)%in%colnames(counts))
target = data[ind,]
target = data
# target = as.data.frame(apply(target,2,gsub,pattern = "-",replacement = "."))
......@@ -1071,6 +1079,7 @@ shinyServer(function(input, output,session) {
counts = dMC$counts
CT_noNorm = dMC$CT_noNorm
CT_Norm = dMC$CT_Norm
## If no file, size factors are estimated
normFactors = SizeFactors_fromFile()$normFactors
......@@ -1508,6 +1517,16 @@ shinyServer(function(input, output,session) {
##
#####################################################
output$PlotVisuTree <- renderTreeWeightD3({
resDiff = ResDiffAnal()
taxo_table = dataInput()$data$taxo
CT_Norm_OTU = dataMergeCounts()$CT_Norm
res = NULL
if(!is.null(resDiff$dds) && length(input$VisuVarInt)>=1) res = Plot_Visu_Tree(input,resDiff,CT_Norm_OTU,taxo_table)
return(res)
})
output$PlotVisuBar <- renderChart({
resDiff = ResDiffAnal()
......@@ -1699,8 +1718,10 @@ shinyServer(function(input, output,session) {
if(input$PlotVisuSelect=="Barplot") res = showOutput("PlotVisuBar")
if(input$PlotVisuSelect=="Heatmap") res = d3heatmapOutput("heatmap", height = input$heightVisu+10)
if(input$PlotVisuSelect=="Boxplot") res = plotOutput("Boxplot", height = input$heightVisu+10)
if(input$PlotVisuSelect=="Tree") res = treeWeightD3Output('PlotVisuTree', height = input$heightVisu+10,width="100%")
if(input$PlotVisuSelect=="Scatterplot" && !input$AddRegScatter) res = scatterD3Output("ScatterplotD3", height = input$heightVisu+10)
if(input$PlotVisuSelect=="Scatterplot" && input$AddRegScatter) res = plotOutput("Scatterplotgg", height = input$heightVisu+10)
if(input$PlotVisuSelect=="Diversity") res = plotOutput("DiversityPlot", height = input$heightVisu+10)
if(input$PlotVisuSelect=="Rarefaction") res = plotOutput("RarefactionPlot",dblclick = "RarefactionPlot_dblclick",brush = brushOpts(id = "RarefactionPlot_brush",resetOnNew = TRUE), height = input$heightVisu+10)
return(res)
......@@ -1941,6 +1962,28 @@ shinyServer(function(input, output,session) {
})
output$VarIntVisuTree <- renderUI({
target=dataInputTarget()$target
data = dataInput()$data
taxo = input$TaxoSelect
resDiff = ResDiffAnal()
res = NULL
if(!is.null(data$counts) && !is.null(data$taxo) && nrow(data$counts)>0 && nrow(data$taxo)>0 && !is.null(taxo) && taxo!="..." && !is.null(target))
{
counts = dataMergeCounts()$counts
Available_x = sort(rownames(counts))
res = selectizeInput("TaxoTree",h6(strong(paste("Select a specific",taxo,sep=" "))),c("...",Available_x),multiple = FALSE)
}
return(res)
})
#####################################################
##
## KRONA
......@@ -1972,8 +2015,23 @@ shinyServer(function(input, output,session) {
## Disable the actionbutton if the number of feature is lower than 2
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)
}
if(input$AddFilter && !is.null(input$SliderThSamp) && !is.null(input$SliderThAb))
{
ind.filter =Filtered_feature(counts,input$SliderThSamp,input$SliderThAb)$ind
......@@ -1984,13 +2042,14 @@ shinyServer(function(input, output,session) {
if (nrow(counts)>=2){
shinyjs::enable("RunDESeq")
}
if (nrow(counts)<2) {
if (nrow(counts)<2 || !testRank) {
shinyjs::disable("RunDESeq")
}
}
if (is.null(counts)) {
if (is.null(counts) || !testRank) {
shinyjs::disable("RunDESeq")
}
})
......
ui.R
View file @ a9956d3a
options(shiny.sanitize.errors = FALSE)
source("css/owncss.R")
source("Rfunctions/Data_Management.R")
library(treeWeightD3)
if (!require(scatterD3)) {
devtools::install_github('aghozlane/scatterD3')
library(scatterD3)
}
if(!require(shinydashboard)){
installed.packages("shinydashboard")
library(shinydashboard)
......@@ -577,7 +582,7 @@ 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","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 = "Barplot")
),
......@@ -592,6 +597,8 @@ body <- dashboardBody(
h5(strong("Select the modalities")),
uiOutput("ModVisu")
),
conditionalPanel(condition="input.PlotVisuSelect=='Tree' ",
uiOutput("VarIntVisuTree")),
conditionalPanel(condition="input.PlotVisuSelect=='Scatterplot' ",
uiOutput("VarIntVisuScatter"),
radioButtons("TransDataScatter","Data transformation",c("Log2 +1" = "log2","None" = "none"),inline=TRUE),
......@@ -599,14 +606,14 @@ body <- dashboardBody(
radioButtons("CorMeth","Correlation method",c("Pearson" = "pearson","Spearman" = "spearman"),inline=TRUE),
checkboxInput("AddRegScatter","Add regression line",FALSE)
),
conditionalPanel(condition="input.PlotVisuSelect!='Rarefaction' && input.PlotVisuSelect!='Diversity' && input.PlotVisuSelect!='Scatterplot'",
conditionalPanel(condition="input.PlotVisuSelect!='Rarefaction' && input.PlotVisuSelect!='Diversity' && input.PlotVisuSelect!='Scatterplot' && input.PlotVisuSelect!='Tree'",
radioButtons("SelectSpecifTaxo","Select the features",c("Most abundant"="Most","All"="All", "Differential features" = "Diff", "Non differential features" = "NoDiff"))
),
conditionalPanel(condition="input.PlotVisuSelect!='Rarefaction' && input.PlotVisuSelect!='Diversity' && input.PlotVisuSelect!='Scatterplot' && (input.SelectSpecifTaxo=='Diff' || input.SelectSpecifTaxo=='NoDiff') ",
conditionalPanel(condition="input.PlotVisuSelect!='Rarefaction' && input.PlotVisuSelect!='Diversity' && input.PlotVisuSelect!='Scatterplot' && input.PlotVisuSelect!='Tree' && (input.SelectSpecifTaxo=='Diff' || input.SelectSpecifTaxo=='NoDiff') ",
selectizeInput("ContrastList_table_Visu","",choices = "", multiple = TRUE),
radioButtons("UnionInterContrasts","Union or intersection ?",c("Union"="Union","Intersection"="Inter"),inline = TRUE)
),
conditionalPanel(condition="input.PlotVisuSelect!='Rarefaction' && input.PlotVisuSelect!='Diversity' && input.PlotVisuSelect!='Scatterplot'",
conditionalPanel(condition="input.PlotVisuSelect!='Rarefaction' && input.PlotVisuSelect!='Diversity' && input.PlotVisuSelect!='Scatterplot' && input.PlotVisuSelect!='Tree'",
uiOutput("TaxoToPlotVisu")
),
......@@ -700,7 +707,7 @@ body <- dashboardBody(
##################
## ALL
##################
conditionalPanel(condition="input.PlotVisuSelect!='Rarefaction' && input.PlotVisuSelect!='Scatterplot'",
conditionalPanel(condition="input.PlotVisuSelect!='Rarefaction' && input.PlotVisuSelect!='Scatterplot' && input.PlotVisuSelect!='Tree'",
radioButtons(inputId = "SensPlotVisu",label = h6(strong("Orientation")),choices = c("Vertical" = "Vertical", "Horizontal" = "Horizontal"),selected = "Vertical",inline = TRUE)
)
),
......
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