Commit 2b5bd3f5 authored by Amine  GHOZLANE's avatar Amine GHOZLANE
Browse files

Debug relative abundance and add KL and JSD distance

parent aefbaef8
...@@ -149,6 +149,10 @@ if (!require(phytools)) { ...@@ -149,6 +149,10 @@ if (!require(phytools)) {
library(phytools) library(phytools)
} }
if(!require(philentropy)){
devtools::install_github("HajkD/philentropy", build_vignettes = TRUE, dependencies = TRUE)
library(philentropy)
}
# if(!require(plotly)){ # if(!require(plotly)){
......
...@@ -411,9 +411,19 @@ PCoAPlot_meta <-function (input, dds, group_init, CT,tree,col = c("SpringGreen", ...@@ -411,9 +411,19 @@ PCoAPlot_meta <-function (input, dds, group_init, CT,tree,col = c("SpringGreen",
counts.norm = counts.norm[,ind_kept] counts.norm = counts.norm[,ind_kept]
# print(head(counts.norm)) # print(head(counts.norm))
## Get the distance ## Get the distance
if(input$DistClust!="sere" && input$DistClust!="Unifrac") dist.counts.norm = vegdist(t(counts.norm), method = input$DistClust)
if(input$DistClust=="sere") dist.counts.norm = as.dist(SEREcoef(counts.norm)) if(input$DistClust=="sere") dist.counts.norm = as.dist(SEREcoef(counts.norm))
if(input$DistClust=="Unifrac") else if(input$DistClust=="jsd"){
#dist_jsd = JSD(sweep(counts.norm,2,colSums(counts.norm)))
dist_jsd = JSD(t(sweep(counts.norm,2,colSums(counts.norm),`/`)))
dist_jsd[is.na(dist_jsd)]=0
dist.counts.norm = as.dist(dist_jsd)
}
else if(input$DistClust=="kl"){
dist_kl = KL(t(sweep(counts.norm,2,colSums(counts.norm),`/`)))
dist_kl[is.na(dist_kl)]=0
dist.counts.norm = as.dist(dist_kl)
}
else if(input$DistClust=="Unifrac")
{ {
tmp = UniFracDist(CT,tree) tmp = UniFracDist(CT,tree)
if(is.null(tree) || is.null(tmp)) dist.counts.norm = NULL if(is.null(tree) || is.null(tmp)) dist.counts.norm = NULL
...@@ -427,7 +437,10 @@ PCoAPlot_meta <-function (input, dds, group_init, CT,tree,col = c("SpringGreen", ...@@ -427,7 +437,10 @@ PCoAPlot_meta <-function (input, dds, group_init, CT,tree,col = c("SpringGreen",
} }
} }
else{
dist.counts.norm = vegdist(t(counts.norm), method = input$DistClust)
}
if(!is.null(dist.counts.norm)) if(!is.null(dist.counts.norm))
{ {
## Do PCoA ## Do PCoA
...@@ -672,9 +685,18 @@ Get_pcoa_table <-function (input, dds, group_init,CT,tree) ...@@ -672,9 +685,18 @@ Get_pcoa_table <-function (input, dds, group_init,CT,tree)
counts.norm = counts.norm[,ind_kept] counts.norm = counts.norm[,ind_kept]
## Get the distance ## Get the distance
if(input$DistClust!="sere" && input$DistClust!="Unifrac") dist.counts.norm = vegdist(t(counts.norm), method = input$DistClust)
if(input$DistClust=="sere") dist.counts.norm = as.dist(SEREcoef(counts.norm)) if(input$DistClust=="sere") dist.counts.norm = as.dist(SEREcoef(counts.norm))
if(input$DistClust=="Unifrac") { else if(input$DistClust=="jsd"){
dist_jsd = JSD(t(sweep(counts.norm,2,colSums(counts.norm),`/`)))
dist_jsd[is.na(dist_jsd)]=0.0
dist.counts.norm = as.dist(dist_jsd)
}
else if(input$DistClust=="kl"){
dist_kl = KL(t(sweep(counts.norm,2,colSums(counts.norm),`/`)))
dist_kl[is.na(dist_kl)]=0.0
dist.counts.norm = as.dist(dist_kl)
}
else if(input$DistClust=="Unifrac") {
tmp = UniFracDist(CT,tree) tmp = UniFracDist(CT,tree)
if(is.null(tree) || is.null(tmp)) dist.counts.norm = NULL if(is.null(tree) || is.null(tmp)) dist.counts.norm = NULL
if(!is.null(tree)) {dist.counts.norm = switch(input$DistClustUnifrac, if(!is.null(tree)) {dist.counts.norm = switch(input$DistClustUnifrac,
...@@ -684,6 +706,9 @@ Get_pcoa_table <-function (input, dds, group_init,CT,tree) ...@@ -684,6 +706,9 @@ Get_pcoa_table <-function (input, dds, group_init,CT,tree)
) )
} }
} }
else {
dist.counts.norm = vegdist(t(counts.norm), method = input$DistClust)
}
if(!is.null(dist.counts.norm)){ if(!is.null(dist.counts.norm)){
## To get always the same result ## To get always the same result
......
...@@ -571,6 +571,7 @@ GetDataToPlot <- function(input,resDiff,VarInt,ind_taxo,aggregate=TRUE,rarefy=FA ...@@ -571,6 +571,7 @@ GetDataToPlot <- function(input,resDiff,VarInt,ind_taxo,aggregate=TRUE,rarefy=FA
counts_tmp = counts_tmp[,colnames(counts_tmp)%in%rownames(targetInt)] counts_tmp = counts_tmp[,colnames(counts_tmp)%in%rownames(targetInt)]
## Proportions over all the taxonomies ## Proportions over all the taxonomies
## Proportion verified
prop_all = t(counts)/rowSums(t(counts)) prop_all = t(counts)/rowSums(t(counts))
prop_all = as.data.frame(prop_all[,Taxonomy%in%ind_taxo]) prop_all = as.data.frame(prop_all[,Taxonomy%in%ind_taxo])
prop_all = as.matrix(prop_all[rownames(prop_all)%in%rownames(targetInt),]) prop_all = as.matrix(prop_all[rownames(prop_all)%in%rownames(targetInt),])
...@@ -586,6 +587,7 @@ GetDataToPlot <- function(input,resDiff,VarInt,ind_taxo,aggregate=TRUE,rarefy=FA ...@@ -586,6 +587,7 @@ GetDataToPlot <- function(input,resDiff,VarInt,ind_taxo,aggregate=TRUE,rarefy=FA
} }
if(!aggregate && nrow(counts_tmp)>0 && nrow(targetInt)>0) if(!aggregate && nrow(counts_tmp)>0 && nrow(targetInt)>0)
{ {
## Proportion verified
counts_tmp_combined = t(counts_tmp) counts_tmp_combined = t(counts_tmp)
prop_tmp_combined = counts_tmp_combined/colSums(counts_tmp) prop_tmp_combined = counts_tmp_combined/colSums(counts_tmp)
rownames(counts_tmp_combined) = targetInt$AllVar rownames(counts_tmp_combined) = targetInt$AllVar
......
...@@ -756,7 +756,7 @@ shinyServer(function(input, output,session) { ...@@ -756,7 +756,7 @@ shinyServer(function(input, output,session) {
## Export in .csv ## Export in .csv
output$ExportRelative <- downloadHandler( output$ExportRelative <- downloadHandler(
filename = function() { 'RelativeAb.csv' }, filename = function() { 'RelativeAb.csv' },
content = function(file){write.csv(dataMergeCounts()$counts/colSums(dataMergeCounts()$counts), file)} content = function(file){write.csv(sweep(dataMergeCounts()$counts,2,colSums(dataMergeCounts()$counts),`/`), file)}
) )
## Export size factors ## Export size factors
...@@ -1551,16 +1551,16 @@ shinyServer(function(input, output,session) { ...@@ -1551,16 +1551,16 @@ shinyServer(function(input, output,session) {
ErrorTree = dataInputTree()$Error ErrorTree = dataInputTree()$Error
TaxoSelect = input$TaxoSelect TaxoSelect = input$TaxoSelect
res = selectInput("DistClust","Distance",c("euclidean", "SERE"="sere", "canberra", "bray", "kulczynski", "jaccard", res = selectInput("DistClust","Distance", c("altGower", "binomial", "bray", "canberra", "cao", "chao", "euclidean","gower", "horn",
"gower", "altGower", "morisita", "horn","mountford","raup","binomial", "jaccard","jsd","kl", "kulczynski", "mahalanobis", "morisita", "mountford","raup",
"chao","cao","mahalanobis"),selected="canberra") "SERE"="sere"),selected="bray")
## Add the unifrac distance ## Add the unifrac distance
if(!is.null(tree) && !is.null(input$fileTree) && is.null(ErrorTree) && TaxoSelect == "OTU/Gene") if(!is.null(tree) && !is.null(input$fileTree) && is.null(ErrorTree) && TaxoSelect == "OTU/Gene")
{ {
res = selectInput("DistClust","Distance",c("euclidean", "SERE"="sere", "canberra", "bray", "kulczynski", "jaccard", res = selectInput("DistClust","Distance",c("altGower", "binomial", "bray", "canberra", "cao", "chao", "euclidean","gower", "horn",
"gower", "altGower", "morisita", "horn","mountford","raup","binomial", "jaccard","jsd","kl", "kulczynski", "mahalanobis", "morisita", "mountford","raup",
"chao","cao","mahalanobis","Unifrac"),selected="canberra") "SERE"="sere","Unifrac"),selected="bray")
} }
return(res) return(res)
......
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