In microsud/microbiomeutilities-shiny: microbiomeutilities: An R package with utility functions for the microbiome R package
library(flexdashboard)
library(microbiome)
library(microbiomeutilities)
# install.packages("picante",repos="http://R-Forge.R-project.org")
# library(picante)
library(picante)
library(data.table)
library(DT)
library(RColorBrewer)
library(shiny)
library(shinythemes)
library(plotly)
library(tibble)
options(shiny.maxRequestSize = 30*1024^2) # 30mb limit
av = function(x){
if( isTRUE(all.equal(x, "")) | isTRUE(all.equal(x, "NULL")) ){
return(NULL)
}
return(x)
}
Data Anlaysis
Sidebar {.sidebar}
p("Data Upload")
# OTU table
h4('Biom/OTU table')
fileInput(inputId = 'upload_biom',
label = helpText("*.biom files must be either as JSON (BIOM-v1) or HDF5 (BIOM-v2)."),
accept = c('.biom'))
tags$hr()
# Mapping file
h4('Mapping file')
fileInput(inputId = 'upload_mapping',
label = helpText("Must be a valid mapping file in *.csv"),
accept = c('.csv'))
tags$hr()
# Tree
h4('Phylo tree file')
fileInput(inputId = 'upload_tree',
label = helpText("Tree must be in newick format with *.tre extension"),
accept = c('.tre'))
tags$hr()
# Example data
actionButton("load_example", label = 'Load Example Data', width = "225px")
h4(' ')
tags$strong()
h4('Notes: ')
tags$strong()
h5('Click on load example to test this app.
Remember this app is for pre-limanry investigation only!')
tags$strong()
h5('contact: sudarshanshetty9@gmail.com')
tags$strong()
h5('LICENSE: GNU General Public License v3.0')
tags$strong()
# Reactive values to store imported data
pseq <- reactive({
# Observe for example data
if(input$load_example > 0){
showNotification("Using example dataset...", type = "message")
data("biogeogut")
ps0 <- biogeogut
return(ps0)
}
# Err if null
if(is.null(input$upload_biom) || is.null(input$upload_mapping)){
return(NULL)
}
showNotification("Importing data and filtering out samples with less than 1000 reads...",
duration = NULL, type = "message")
# Run import function once all files are uploaded
ps0 <- read_phyloseq(otu.file = input$upload_biom$datapath,
metadata.file = input$upload_mapping$datapath,
type = "biom")
if(!is.null(input$upload_tree)){
treefil = read_tree(input$upload_tree$datapath)
ps0 <- merge_phyloseq(ps0, treefil)
} else{
warning("No tree provided. Phylogenetic investigation cannot be done.")
ps0 <- ps0
}
# Remove any OTUs with 0 abundance across all samples
ps0 <- prune_taxa(taxa_sums(ps0) > 0, ps0)
ps0 <- prune_samples(sample_sums(ps0) > 50, ps0)
# Return phyloseq object
return(ps0)
})
# Reactive values to store imported data
pseq_meta <- reactive({
# Err if null
if(is.null(input$upload_mapping)){
return(NULL)
}
# Show message
showNotification("Importing data and filtering out samples with less than 1000 reads...", type = "message")
# Run import function once all files are uploaded
pseq_metadf <- read.csv(input$upload_mapping$datapath, row.names=1,check.names=FALSE)
# Return meta object
return(pseq_metadf)
})
Column {data-width=1000} {.tabset .tabset-fade}
Library size distribution
output$LibSizeHist <- renderPlot({
validate(need(expr = !is.null(pseq()), message = "Please upload data..."))
ps0<- pseq()
# Plot sample depths
SeqDepth = colSums(otu_table(ps0))
sample_data(ps0)$SeqDepth = SeqDepth
meta_df <- meta(ps0)
p <- ggpubr::gghistogram(meta_df, "SeqDepth") + ggpubr::rotate_x_text()
return(p)
}, width = 400, height = 400)
fillPage(plotOutput("LibSizeHist"))
OTU counts
fluidPage(
titlePanel(""),
sidebarLayout(
sidebarPanel(width = 3,
numericInput(
"PrevalenceOTU",
label = "Prevalence: Number of samples OTU must be present 0.05 means 5%",0,
min = 0,
max = 0.9),
numericInput(
"DetectionOTU",
label = "Detection: Minimum count of OTUs",0,
min = 0,
max = 1000)
),
output$OTUHist_1 <- renderPlot({
validate(need(expr = !is.null(pseq()), message = "Please upload data..."))
ps0 <- pseq()
ps.fil <- core(ps0, detection = input$DetectionOTU, prevalence = input$PrevalenceOTU)
no_taxa <- phyloseq::ntaxa(ps.fil)
Text <- paste0("Numer of OTU/Taxa after filtering: ", no_taxa)
# Plot sample depths
px.log <- microbiome::transform(ps.fil, "log10")
px.clr <- microbiome::transform(ps.fil, "clr")
px.hel <- microbiome::transform(ps.fil, "hellinger")
p1 <- ggplot2::qplot(rowSums(otu_table(ps0))) + ggplot2::theme_bw() + ggplot2::xlab("Raw counts") + ggtitle("A", subtitle =Text )
p2 <- ggplot2::qplot(rowSums(otu_table(px.log))) + ggplot2::theme_bw() + ggplot2::xlab("log10(1+x) counts") + ggtitle("B", subtitle =Text )
p3 <- ggplot2::qplot(rowSums(otu_table(px.clr))) + ggplot2::theme_bw() + ggplot2::xlab("clr counts") + ggtitle("C", subtitle =Text )
p4 <- ggplot2::qplot(rowSums(otu_table(px.hel))) + ggplot2::theme_bw() + ggplot2::xlab("Hellinger counts") + ggtitle("D", subtitle =Text )
#t <- gridExtra::textGrob(paste0("Numer of OTU/Taxa after filtering: ", no_taxa) )
return(gridExtra::grid.arrange(p1,p2,p3,p4, nrow= 2, ncol= 2))
mainPanel(plotOutput("OTUHist_1"))
}, width = 600, height = 600)
))
Check different transformations and distribution of OTU counts in overall data.
Note down the minimum counts and prevalance used.
These parameters can be used to filter potentailly spurious OTUs/taxa before doing differential abundance analysis using:
1. DESeq2.
2. edgeR.
3. limma.
Note: The filtered data is not used for next analysis. For each of the tabs, raw data is used and individual filtering as specifed.
OTU counts varied taxa
# Select taxonomic level
selectInput(
inputId = "OtuTaxaDist",
label = "Choose taxa level",
choices = list(
"Domain" = "Domain",
"Phylum" = "Phylum",
"Class" = "Class",
"Order" = "Order",
"Family" = "Family",
"Genus" = "Genus"
),
selected = "Phylum")
output$OTUtaxa <- renderPlotly({
validate(need(expr = !is.null(pseq()), message = "Loading data..."))
ps0 <- pseq()
taxasums = rowSums(otu_table(ps0))
taxatable <- as.data.frame.matrix(tax_table(ps0))
tax_plot2 <- ggplotly(tax_plot1 <-
ggplot(taxatable, aes(x = taxasums, color = taxatable[, input$OtuTaxaDist])) +
geom_line(size = 1.5, stat = "density") +
xlab("Log10 OTU Counts") + theme_bw() + scale_x_log10() + guides(color=guide_legend(title=input$OtuTaxaDist, ncol = 3))
+ theme(legend.position="bottom"), width = 800, inline = T)
return(tax_plot2)
})
fluidPage(plotlyOutput("OTUtaxa"))
Summarise phyloseq
fillPage(renderPrint({
validate(need(expr = !is.null(pseq()), message = "Loading data..."))
ps0 <- pseq()
#ps_meta <- as.data.table(meta(ps0))
# Plot sample depths
smp <- summarize_phyloseq(ps0)
return(smp)
}))
Taxa Prevalence
# Select taxonomic level
selectInput(
inputId = "PrevalencePlot",
label = "Choose taxa level",
choices = list(
"Domain" = "Domain",
"Phylum" = "Phylum",
"Class" = "Class",
"Order" = "Order",
"Family" = "Family"
),
selected = "Phylum"
)
#Prevalence plot across all samples
output$PrevalencePlot <- renderPlot({
# Err if null
validate(need(expr = !is.null(pseq()), message = "Loading data..."))
ps0 <- pseq()
tax_prev_plot <- plot_taxa_prevalence(ps0, input$PrevalencePlot)
message("Generating prevalence plot")
tax_prev_plot <- tax_prev_plot + theme(legend.position = "bottom")
return(tax_prev_plot)
}, width = 1200, height = 800)
fluidPage(plotOutput("PrevalencePlot"))
Taxa barplot
fluidPage(titlePanel(""),
sidebarLayout(
sidebarPanel(
width = 2,
column = 5,
## plot variable
h5('Please select the plotting parameters'),
selectInput(
inputId = "TaxBarplot",
label = "Choose taxa level",
choices = list(
"Domain" = "Domain",
"Phylum" = "Phylum",
"Class" = "Class",
"Order" = "Order",
"Family" = "Family",
"Genus" = "Genus"
),
selected = "Phylum"
),
selectInput(
"taxbarcolors",
label = "Choice of palette default Paired",
choices = c("Set2", "Paired", "Set1", "Accent", "Dark2", "Set3"),
selected = "Paired"
),
numericInput(
"topOTUBar",
label = "Top taxa to be viewed",
5,
min = 3,
max = 100
)
),
output$boxplot <- renderPlot({
validate(need(expr = !is.null(pseq()), message = "Loading data..."))
ps.bar <- pseq()
taxic <- as.data.frame(ps.bar@tax_table)
if (input$taxbarcolors == "Set3" |
input$taxbarcolors == "Paired") {
getPalette = colorRampPalette(brewer.pal(12, input$taxbarcolors))
}
else if (input$taxbarcolors == "Set2" |
input$taxbarcolors == "Dark2" | input$taxbarcolors == "Accent")
{
getPalette = colorRampPalette(brewer.pal(8, input$taxbarcolors))
} else if (input$taxbarcolors == "Set1")
{
getPalette = colorRampPalette(brewer.pal(9, input$taxbarcolors))
}
otu.df <- as.data.frame(otu_table(ps.bar))
taxic$OTU <- rownames(taxic)
colnames(taxic)
taxmat <- as.matrix(taxic)
new.tax <- tax_table(taxmat)
tax_table(ps.bar) <- new.tax
guide_italics <-
guides(fill = guide_legend(
label.theme = element_text(
size = 15,
face = "italic",
colour = "Black",
angle = 0
)
))
if (!is.null(ps.bar@phy_tree)) {
ps.bar@phy_tree <- NULL
}
ps.bar.tlev <- aggregate_taxa(ps.bar,
level = input$TaxBarplot,
top = input$topOTUBar)
tax.lev <- as.data.frame.matrix(tax_table(ps.bar.tlev))
colourCount = length(unique(tax.lev$OTU))
ps.bar.rel <- microbiome::transform(ps.bar.tlev, "compositional")
plot.composition.relAbun <-
plot_composition(ps.bar.rel) +
scale_fill_manual(input$TaxBarplot, values = getPalette(colourCount)) + theme_bw() +
theme(axis.text.x = element_text(angle = 90)) +
ggtitle(paste0(
"Relative abundance of top ",
input$topOTUBar,
sep = " ",
input$TaxBarplot
)) + guide_italics + theme(legend.title = element_text(size =
18)) + theme(legend.position = "top")
return(plot.composition.relAbun)
mainPanel(plotOutput("boxplot"))
}, width = 1000, height = 600)
))
Taxa boxplot
fluidPage(
titlePanel(""),
sidebarLayout(
sidebarPanel(width = 2, column = 4,
## plot variable
h5('Please type the variable of choice from mapping file'),
p('No numeric values, only categorical'),
textInput('Sam_var2', 'Variable (eg. SampleType)'),
tags$div(id = 'placeholder1'),
hr(),
selectInput(
inputId = "Taxplot",
label = "Choose taxa level",
choices = list(
"Domain" = "Domain",
"Phylum" = "Phylum",
"Class" = "Class",
"Order" = "Order",
"Family" = "Family",
"Genus" = "Genus",
"Species" = "Species"),
selected = "Phylum"),
selectInput(
"taxcolors",
label = "Choice of palette default Paired",
choices = c("Set2", "Paired", "Set1", "Accent", "Dark2", "Set3"),
selected = "Paired"),
numericInput(
"topOTU",
label = "Top taxa to be viewed",5,
min = 0,
max = 50)
),
output$boxplot <- renderPlot({
# validate(
# need(expr = !is.null(pseq()), message = "Loading data...")
# )
validate(
need(expr = input$Sam_var2, message = "Type the variable!!")
)
ps.box <- pseq()
#req(input$Sam_var2)
ps.box.rel <- microbiome::transform(ps.box, "compositional")
pn1 <- plot_taxa_boxplot(ps.box.rel,
taxonomic.level = input$Taxplot,
top.otu = input$topOTU,
VariableA = input$Sam_var2,
title = "Relative abundance ",
color = input$taxcolors)
return(pn1 + ggpubr::rremove("grid"))
mainPanel(plotOutput("boxplot"))
}, width = 1000, height = 600)))
Taxa Heatmap
fluidPage(
titlePanel(""),
sidebarLayout(
sidebarPanel(width = 2, column = 4,
## link the JS file
tags$head(tags$script(src="script.js")),
## link the CSS file
tags$head(tags$link(rel="stylesheet",
type="text/css",
href="style.css")),
## plot variable
h5('Please type the variable of choice from mapping file'),
p('No numeric values, only categorical'),
textInput('txt1', 'Variable (eg. SampleType)'),
tags$div(id = 'placeholder1'),
hr(),
selectInput(
"transformation",
label = "Choice of transformation default log10",
choices = c("log10", "compositional", "clr", "Z-OTU"),
selected = "log10"),
selectInput(
"heatcolors",
label = "Choice of palette default GnBu",
choices = c("RdYlBu", "Blues", "Greys", "GnBu", "RdGy", "Spectral"),
selected = "GnBu"),
numericInput(
"SubsetOTU",
label = "Number of top OTUs to visualize (max.150)",10,
min = 5,
max = 150)
),
output$HeatMap <- renderPlot({
validate(
need(expr = !is.null(pseq()), message = "Please upload data...")
)
validate(
need(expr = input$txt1, message = "Type the variable!!")
)
ps0<- pseq()
#ps0.f <- microbiomeutilities::format_to_besthit(pseq())
mainPanel(plot_taxa_heatmap(ps0, subset.top = input$SubsetOTU,
VariableA = input$txt1,
heatcolors = brewer.pal(100, input$heatcolors),
transformation = input$transformation))
}, width = 1000, height = 700)))
Diversity (Non-Phylogenetic)
fluidPage(
titlePanel("Non-Phylogenetic diversity"),
sidebarLayout(
sidebarPanel(width = 2, column = 4,
## link the JS file
tags$head(tags$script(src="script.js")),
## link the CSS file
tags$head(tags$link(rel="stylesheet",
type="text/css",
href="style.css")),
## plot variable
h5('Please type the variable of choice from mapping file'),
p('No numeric values, only categorical'),
textInput('Sam_var3', 'Variable (eg. SampleType)'),
tags$div(id = 'placeholder1'),
hr(),
selectInput(
"typeDiv",
label = "Choice of diversities measure default: diversities",
choices = c("diversities", "dominance", "evenness"),
selected = "diversities"),
selectInput(
"divcolors",
label = "Choice of palette default: Set1",
choices = c("Set2", "Paired", "Set1", "Accent", "Dark2", "Set3"),
selected = "Set1"),
selectInput(
"figtype",
label = "Choice of plot default: stripchart",
choices = c("stripchart", "boxplot", "violin"),
selected = "stripchart")
),
output$diveplot <- renderPlot({
validate(
need(expr = !is.null(pseq()), message = "Loading data...")
)
validate(
need(expr = input$Sam_var3, message = "Type the variable!!")
)
ps0<- pseq()
p <- plot_alpha_diversities(ps0,
type = input$typeDiv,
index.val = "all",
plot.type = input$figtype,
variableA = input$Sam_var3,
palette = input$divcolors) + ggpubr::rremove("grid")
return(p)
mainPanel(plotOutput("diveplot"))
}, width = 900, height = 600)
))
Phylogenetic diversity
fluidPage(
titlePanel("Phylogenetic diversity"),
sidebarLayout(
sidebarPanel(width = 2, column = 3,
## link the JS file
tags$head(tags$script(src="script.js")),
## link the CSS file
tags$head(tags$link(rel="stylesheet",
type="text/css",
href="style.css")),
## plot variable
h5('Please type the variable of choice from mapping file'),
p('No numeric values, only categorical'),
textInput('Sam_var4', 'Variable (eg. SampleType)'),
tags$div(id = 'placeholder1'),
hr(),
selectInput(
"pdivcolors",
label = "Choice of palette default Set2",
choices = c("Set2", "Paired", "Set1", "Accent", "Dark2", "Set3"),
selected = "Set2"),
selectInput(
"figPtype",
label = "Choice of plot default: stripchart",
choices = c("stripchart", "boxplot", "violin"),
selected = "stripchart")
),
output$diveplot <- renderPlot({
av = function(x){
if( isTRUE(all.equal(x, "")) | isTRUE(all.equal(x, "NULL")) ){
return(NULL)
}
return(x)
}
validate(
need(expr = !is.null(pseq()), message = "Please upload data...")
)
validate(
need(expr = !is.null(pseq()@phy_tree), message = "No phylogenetic tree!")
)
ps0<- pseq()
if(!is.null(ps0@phy_tree)){
otu_table_ps3 <- as.data.frame(ps0@otu_table)
treefil2 = read_tree(input$upload_tree$datapath)
df.pd <- pd(t(otu_table_ps3), treefil2,include.root=F)
metadata_table_ps3 <- meta(ps0)
sample_data(ps0)$Phyogenetic_diversity <- df.pd$PD
#metadata_table_ps3$Phyogenetic_diversity <- df.pd$PD
mps0 <- phy_to_ldf(ps0, transform.counts = NULL)
colnames(mps0)[colnames(mps0) == input$Sam_var4] <- 'Variable_alias'
# mps0 <- reshape2::melt(metadata_table_ps3)
if (input$figPtype == "boxplot")
{
plot.pd <- ggpubr::ggboxplot(mps0, "Variable_alias", "Phyogenetic_diversity",
palette = input$pdivcolors,
fill = "Variable_alias",
legend = "right", add = "mean_sd"
)}
else if (input$figPtype == "stripchart")
{
plot.pd <- ggpubr::ggstripchart(mps0, "Variable_alias", "Phyogenetic_diversity",
palette = input$pdivcolors,
color = "Variable_alias",
legend = "right", add = "mean_sd", shape = 21)
}
else if (input$figPtype == "violin")
{
plot.pd <- ggpubr::ggviolin(mps0, "Variable_alias", "Phyogenetic_diversity",
palette = input$pdivcolors,
fill = "Variable_alias",
legend = "right", add = "mean_sd")
}
} else{
# Show message
showNotification("No tree provided. Phylogenetic diversity, unifrac distances cannot be calculated.", duration = NULL, type = "error")
plot.pd <- "No tree provided. Phylogenetic diversity, unifrac distances cannot be calculated"
}
return(plot.pd + ggpubr::rremove("grid"))
mainPanel(plotOutput("diveplot"))
}, width = 900, height = 600)))
Ordination
fluidPage(
titlePanel("Beta diversity"),
sidebarLayout(
sidebarPanel(width = 2, column = 3,
## link the JS file
tags$head(tags$script(src="script.js")),
## link the CSS file
tags$head(tags$link(rel="stylesheet",
type="text/css",
href="style.css")),
## plot variable
h5('Please type the variable of choice from mapping file'),
textInput('txt3', 'Variable (eg. SampleType)'),
tags$div(id = 'placeholder1'),
hr(),
selectInput(
"Distance",
label = "Choice of distance method default: Bray-Curtis",
choices = c("wunifrac", "unifrac", "bray"),
selected = "bray"),
selectInput(
"Varcolors",
label = "Choice of palette default Set1",
choices = c("Set2", "Paired", "Set1", "Accent", "Dark2", "Set3"),
selected = "Set1"),
numericInput(
"Prevalence",
label = "Prevalence: Number of samples OTU must be present 0.05 means 5%",0,
min = 0,
max = 0.9),
numericInput(
"Detection",
label = "Detection: Minimum relative abundance of OTUs",0,
min = 0,
max = 0.1)
),
output$Ordination <- renderPlot({
av = function(x) {
if (isTRUE(all.equal(x, "")) | isTRUE(all.equal(x, "NULL"))) {
return(NULL)
}
return(x)
}
validate(need(expr = !is.null(pseq()), message = "Please upload data..."))
ps0 <- pseq()
ps0a.rel <- microbiome::transform(ps0, "compositional")
ps0a.rel <-
core(ps0a.rel,
detection = input$Detection,
prevalence = input$Prevalence)
#ps0.f <- microbiomeutilities::format_to_besthit(pseq())
if (input$Distance == "wunifrac") {
ord <- ordinate(ps0a.rel, method = "PCoA", input$Distance)
ordiplot <-
plot_ordination(ps0a.rel, ord, color = av(input$txt3)) + ggtitle("Weighted Unifrac Compositional") + theme_bw() + geom_point(size = 3) + scale_color_brewer(palette = input$Varcolors) + ggpubr::rremove("grid")
return(ordiplot)
} else if (input$Distance == "unifrac") {
ord <-
ordinate(ps0a.rel,
method = "PCoA",
input$Distance,
weighted = F)
ordiplot <-
plot_ordination(ps0a.rel, ord, color = av(input$txt3)) + ggtitle("UnWeighted Unifrac Compositional") + theme_bw() + geom_point(size = 3) + scale_color_brewer(palette = input$Varcolors) + ggpubr::rremove("grid")
return(ordiplot)
} else if (input$Distance == "bray") {
ord <- ordinate(ps0a.rel, method = "PCoA", input$Distance)
ordiplot <-
plot_ordination(ps0a.rel, ord, color = av(input$txt3)) +
ggtitle("Bray-Curtis Compositional") + theme_bw() + geom_point(size = 3) + scale_color_brewer(palette = input$Varcolors) + ggpubr::rremove("grid")
return(ordiplot)
}
mainPanel(plotOutput("diveplot"))
}, width = 800, height = 600)
))
Core microbiota
fluidPage(
titlePanel(""),
sidebarLayout(
sidebarPanel(width = 2,column = 3,
numericInput(
"PrevalenceCore",
label = "Prevalence: Number of samples OTU must be present 0.05 means 5%",0.5,
min = 0.25,
max = 0.99),
h5('Minimum detection of taxa is set
at 0.0001% see x-axis for range')
),
output$core_plots <- renderPlot({
av = function(x){
if( isTRUE(all.equal(x, "")) | isTRUE(all.equal(x, "NULL")) ){
return(NULL)
}
return(x)
}
validate(
need(expr = !is.null(pseq()), message = "Please upload data...")
)
ps0<- pseq()
#ps0.sub <- phyloseq::subset_samples(ps0, av(input$txt5) == input$txt6)
ps0a.rel <- microbiome::transform(ps0, "compositional")
ps0a.rel <- format_to_besthit(ps0a.rel)
prevalences <- seq(.05, 1, .05)
detections <- 10^seq(log10(1e-4), log10(.2), length = 10)
core_plot <- plot_core(ps0a.rel, plot.type = "heatmap",
prevalences = prevalences,
detections = detections,
colours = rev(brewer.pal(5, "Spectral")),
min.prevalence = input$PrevalenceCore, horizontal = F)
return(core_plot + ggtitle(paste0("Core taxa at ", (input$PrevalenceCore*100), " % prevalence and min. abundance 0.0001%")))
mainPanel(plotOutput("core_plots"))
}, width = 800, height = 600)
))
Data tables
Tables {data-width=400} {.tabset .tabset-fade}
Metadata
get_meta_data = reactive({
ps0 <- pseq()
metadfx <- meta(ps0)
return(metadfx)
})
renderDataTable({
get_meta_data()
})
Diversities table
output$DiversitiesTab = DT::renderDataTable({
# Err if null
validate(
need(expr = !is.null(pseq()), message = "Please upload data...")
)
ps1<- pseq()
ps.diver <- diversities(ps1)
return(ps.diver)
},
server = F,
extensions = 'Buttons',
options = list(lengthChange = TRUE,
dom = 'Bfrtip',
buttons = c('copy', 'print'),
autoWidth = F,
pageLength = 15,
searchHighlight = TRUE,
scrollX = F)
)
fillPage(DT::dataTableOutput('DiversitiesTab'))
Dominance table
output$DominanceTab = DT::renderDataTable({
# Err if null
validate(
need(expr = !is.null(pseq()), message = "Please upload data...")
)
ps1<- pseq()
ps.diver <- dominance(ps1)
return(ps.diver)
},
server = F,
extensions = 'Buttons',
options = list(lengthChange = TRUE,
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
autoWidth = F,
pageLength = 15,
searchHighlight = TRUE,
scrollX = F)
)
fillPage(DT::dataTableOutput('DominanceTab'))
Evenness table
output$EvennessTab = DT::renderDataTable({
# Err if null
validate(
need(expr = !is.null(pseq()), message = "This will take time depending on your data size...")
)
showNotification("This will take time depending on your data size...", type = "message")
ps1<- pseq()
ps.diver <- evenness(ps1)
return(ps.diver)
},
server = F,
extensions = 'Buttons',
options = list(lengthChange = TRUE,
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
autoWidth = F,
pageLength = 15,
searchHighlight = TRUE,
scrollX = F)
)
fillPage(DT::dataTableOutput('EvennessTab'))
Phylogenetic diversity table
output$PhylogenyTab = DT::renderDataTable({
validate(
need(expr = !is.null(pseq()), message = "Please upload data...")
)
ps0<- pseq()
if(!is.null(ps0@phy_tree)){
otu_table_ps3 <- as.data.frame(ps0@otu_table)
treefil2 = read_tree(input$upload_tree$datapath)
df.pd <- pd(t(otu_table_ps3), treefil2,include.root=F)
return(df.pd)
}
else{
# Show message
showNotification("No tree provided. Phylogenetic diversity, unifrac distances cannot be calculated.", duration = NULL, type = "error")
plot.pd <- "No tree provided. Phylogenetic diversity, unifrac distances cannot be calculated"
}
},
server = F,
extensions = 'Buttons',
options = list(lengthChange = TRUE,
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
autoWidth = F,
pageLength = 15,
searchHighlight = TRUE,
scrollX = F))
fluidPage(DT::dataTableOutput('PhylogenyTab'))
Sample-Features {data-width=150}
Number of samples
output$sampleN <- renderValueBox({
# Err if NULL
if(is.null(pseq())){
return(valueBox(0, caption = "Number of samples", icon = "fa-cubes", color = "black"))
}
valueBox(nsamples(pseq()), caption = "Number of samples", icon = "fa-cubes", color = "black")
})
valueBoxOutput(outputId = "sampleN", width = 3)
Number of OTU's
output$taxaN <- renderValueBox({
# Err if NULL
if(is.null(pseq())){
return(valueBox(0, caption = "Number of taxa", icon = "fa-users", color = "black"))
}
valueBox(ntaxa(pseq()), caption = "Number of taxa", icon = "fa-cubes", color = "black")
})
valueBoxOutput(outputId = "taxaN", width = "auto")
Number of samples
output$totalDepth <- renderValueBox({
# Err if NULL
if(is.null(pseq())){
return(valueBox(0, caption = "Total sampling depth", icon = "fa-bar-chart", color = "black"))
}
valueBox(sum(sample_sums(pseq())), caption = "Total sampling depth", icon = "fa-info-circle", color = "black")
})
valueBoxOutput(outputId = "totalDepth")
Mean sampling depth
output$meanDepth <- renderValueBox({
# Err if NULL
if(is.null(pseq())){
return(valueBox(0, caption = "Mean sampling depth", icon = "fa-info-circle", color = "black"))
}
valueBox(round(mean(sample_sums(pseq())), 0), caption = "Mean sampling depth", icon = "fa-info-circle",color = "black" )
})
valueBoxOutput(outputId = "meanDepth")
Median sampling depth
output$medianDepth <- renderValueBox({
# Err if NULL
if(is.null(pseq())){
return(valueBox(0, caption = "Median sampling depth", icon = "fa-info-circle", color = "black"))
}
valueBox(round(median(sample_sums(pseq())), 0), caption = "Median sampling depth", icon = "fa-info-circle", color = "black")
})
valueBoxOutput(outputId = "medianDepth")
Range of sampling depth
output$rangeDepth <- renderValueBox({
# Err if NULL
if(is.null(pseq())){
return(valueBox(0, caption = "Range of sampling depth", icon = "fa-info-circle", color = "black"))
}
range <- paste0(min(sample_sums(pseq())), " - ", max(sample_sums(pseq())))
valueBox(range, caption = "Range of sampling depth", icon = "fa-info-circle", color = "black")
})
valueBoxOutput(outputId = "rangeDepth")
Standard deviation of sampling depth
output$sdDepth <- renderValueBox({
# Err if NULL
if(is.null(pseq())){
return(valueBox(0, caption = "Standard deviation of sampling depth", icon = "fa-info-circle", color = "black"))
}
sd <- round(sd(sample_sums(pseq())), digits = 1)
valueBox(sd, caption = "Standard deviation of sampling depth", icon = "fa-info-circle", color = "black")
})
valueBoxOutput(outputId = "sdDepth")
Singletons
output$sington <- renderValueBox({
# Err if NULL
if(is.null(pseq())){
return(valueBox(0, caption = "Singletons in data", icon = "fa-info-circle", color = "black"))
}
singletons <- round(length(taxa_sums(pseq())[taxa_sums(pseq()) <= 1]), digits = 1)
valueBox(singletons, caption = "Singletons in data", icon = "fa-info-circle", color = "black")
})
valueBoxOutput(outputId = "sington")
Sparsity
output$spar <- renderValueBox({
# Err if NULL
if(is.null(pseq())){
return(valueBox(0, caption = "Sparsity", icon = "fa-info-circle", color = "black"))
}
sparsit <- round(length(which(abundances(pseq()) == 0))/length(abundances(pseq())), digits = 1)
valueBox(sparsit, caption = "Sparsity", icon = "fa-info-circle", color = "black")
})
valueBoxOutput(outputId = "spar")
LICENSE
output$LICENSE <- renderText({
print("GNU General Public License v3.0")
})
textOutput(outputId = "LICENSE")
Exploratory
Exploratory-tabs {data-width=400} {.tabset .tabset-fade}
Taxa features
fluidPage(titlePanel(""),
sidebarLayout(
sidebarPanel(
width = 2,
column = 4,
## plot variable
h5('Please type the variable of choice from mapping file'),
p('No numeric values, only categorical'),
textInput('Sam_varFeat', 'Variable (eg. SampleType)'),
tags$div(id = 'placeholder1'),
hr(),
sliderInput(
"DetectionFeat",
label = "Detection: Minimum relative abundance of OTUs",
min = 0.0001,
max = 0.1,
value = 0.001,
step = 0.005
),
sliderInput(
"bins",
"Bins:",
min = 0,
max = 500,
value = 100
)
),
output$tableplot <- renderPlot({
validate(need(expr = !is.null(pseq()), message = "Please upload data..."))
ps1 <- pseq()
ps1.rel <- microbiome::transform(ps1, "compositional")
ps1.rel.filter <-
microbiome::core(ps1.rel,
detection = input$DetectionFeat,
prevalence = 0)
ps1.df <- phy_to_ldf(ps1.rel.filter, transform.counts = NULL)
ps1.df2 <-
ps1.df[, c("Phylum",
"Order",
"Class",
"Family" ,
"Abundance",
input$Sam_varFeat)]
return(
tabplot::tableplot(
ps1.df2,
fontsize = 10,
legend.lines = 20,
nBins = input$bins,
sampleBinSize = 1e2,
sortCol = Abundance,
relative = TRUE,
scales = "lin"
)
)
mainPanel(plotOutput("tableplot"))
}, width = 1200, height = 800)
))
Insights/Notes
Microbiome analytics tool
Pre-processing is an essential step and in many instance is not given importance, as we biologists do not take these steps as seriously as it should be taken. To avoid this there are some common approaches to look at the data. On purpose this webapp does not allow for downloading high-quality images. It is just to have a quick look at your data and then play around in R using Rmarkdown to make you codes easily traceble and results reproducible.
This shiny app is inspired from from other shiny tools such as shiny-phyloseq and btools-shiny.
This first-hand look at your data will aid in making different choices for more thorough and in-dept investagion.
You will need to install microbiomeutilities R package from github repo.
The main tools used here are Phyloseq and microbiome.
Check library size distribution
Check the library sizes for each of the samples within the main variable. Do you think it is okay or there is difference in library sizes and will this affect you downstream statistical analysis?
Check OTU/Taxa count distribution
If the OTU counts are higly skewed in distribution try normalising them.
Check OTU/Taxa count distribution
Which taxa are distributied differently. Check if they are of any interest in the ecosystem being investigated.
Summary of Phyloseq
Look if there are singletons, sparsity in your data.
Taxa prevalence
Which taxa are prevalant and what abundances. Are there OTUs with low abundace and low prevalance?
Taxa barplot and Taxa boxplot
Gives a firsthand view of the microbial composition.
Taxa heatmap
At OTU level, see if there is any difference in top N otus between variables of interest.
Alpha diversity
This includes the non-phylogenetic metrices. For more on this check Microbiome:Diversities
Phylogenetic diversity
For more on this check Picante. Is calculated only if tree is provided in the upload tab.
Beta Diversity (Ordinations)
Here, commonly used metrics for beta diversity are depicted. The default transformation is compositional (relative abundance). Choose filtering option Detections and Prevalance to remove potentially spurious OTUs/Taxa.
For more information:
Waste Not, Want Not: Why Rarefying Microbiome Data Is Inadmissible.
Normalisation and data transformation.
What is Constrained and Unconstrained Ordination.
R packages that work at the backend of this webapp are:
- picante
- data.table
- RColorBrewer
- shinythemes
- DT
- shiny
- reshape
- ggpubr
- tibble
- pheatmap
- Phyloseq
- microbiome
- plyr
- tidyverse R packages
- pheatmap
- Vegan
- flexdashboard
- Biobase
- plotly
- ape
Please let me know if I have missed any of the packages.
Terms and conditions
Microbiome analytics tool
Terms of Use
This web application is hosted on a Shinyapps.io server https://www.shinyapps.io/.
If you have concerns about the terms of use for this web hosted application, please run the app locally on your computer.
By using this app you are agreeing to the terms of use as described by Shinnyaps.io: https://www.rstudio.com/about/shinyapps-terms-use/
The authors and maintainers of microbiomeutilities-shiny app will not save your data on our servers. However, as the Shinyapps server is not HIPAA compliant, you must refrain from uploading protected health information or confidential data with this app. You may instead download the code and run the app locally on your private computer and network (see above). We the authors cannot be held responsible for the confidentiality, availability, security, loss, misuse or misappropriation of any data you submit to this application.
One of the Important information from Shinyapps terms of use (https://www.rstudio.com/about/shinyapps-terms-use/): "If you choose to upload data to an application you are using via the RStudio Service, you acknowledge and agree you are giving certain legal rights to the licensor of the application to process and otherwise use your data. Please carefully review any license terms accompanying any application to which you submit your data for the legal rights which you are giving the application licensor. Further, RStudio does not claim ownership of your data; however, you hereby grant RStudio a worldwide, perpetual, irrevocable, royalty-free, fully paid up, transferable and non-exclusive license, as applicable, to use and copy your data in connection with making your data available to the application (and the licensor of such application) to which you have submitted your data. You acknowledge and agree that even if you remove your data from the RStudio Service, your data may have been downloaded by, and remain accessible to, the licensors of those applications to which you submitted your data. Accordingly, do not submit data which you desire to remain confidential or which you wish to limit the right to access or use. You should never submit to the RStudio Service any data which consists of personally identifiable information, credit card information, or protected health information, as such terms are defined by relevant laws, rules, and regulations. RSTUDIO IS NOT RESPONSIBLE FOR THE CONFIDENTIALITY, AVAILABILITY, SECURITY, LOSS, MISUSE OR MISAPPROPRIATION OF ANY DATA YOU SUBMIT TO THE RSTUDIO SERVICE OR ANY APPLICATION MADE AVAILABLE VIA THE RSTUDIO SERVICE."
microsud/microbiomeutilities-shiny documentation built on May 7, 2019, 9:38 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(flexdashboard) library(microbiome) library(microbiomeutilities) # install.packages("picante",repos="http://R-Forge.R-project.org") # library(picante) library(picante) library(data.table) library(DT) library(RColorBrewer) library(shiny) library(shinythemes) library(plotly) library(tibble) options(shiny.maxRequestSize = 30*1024^2) # 30mb limit av = function(x){ if( isTRUE(all.equal(x, "")) | isTRUE(all.equal(x, "NULL")) ){ return(NULL) } return(x) }
Data Anlaysis
Sidebar {.sidebar}
p("Data Upload") # OTU table h4('Biom/OTU table') fileInput(inputId = 'upload_biom', label = helpText("*.biom files must be either as JSON (BIOM-v1) or HDF5 (BIOM-v2)."), accept = c('.biom')) tags$hr() # Mapping file h4('Mapping file') fileInput(inputId = 'upload_mapping', label = helpText("Must be a valid mapping file in *.csv"), accept = c('.csv')) tags$hr() # Tree h4('Phylo tree file') fileInput(inputId = 'upload_tree', label = helpText("Tree must be in newick format with *.tre extension"), accept = c('.tre')) tags$hr() # Example data actionButton("load_example", label = 'Load Example Data', width = "225px") h4(' ') tags$strong() h4('Notes: ') tags$strong() h5('Click on load example to test this app. Remember this app is for pre-limanry investigation only!') tags$strong() h5('contact: sudarshanshetty9@gmail.com') tags$strong() h5('LICENSE: GNU General Public License v3.0') tags$strong()
# Reactive values to store imported data pseq <- reactive({ # Observe for example data if(input$load_example > 0){ showNotification("Using example dataset...", type = "message") data("biogeogut") ps0 <- biogeogut return(ps0) } # Err if null if(is.null(input$upload_biom) || is.null(input$upload_mapping)){ return(NULL) } showNotification("Importing data and filtering out samples with less than 1000 reads...", duration = NULL, type = "message") # Run import function once all files are uploaded ps0 <- read_phyloseq(otu.file = input$upload_biom$datapath, metadata.file = input$upload_mapping$datapath, type = "biom") if(!is.null(input$upload_tree)){ treefil = read_tree(input$upload_tree$datapath) ps0 <- merge_phyloseq(ps0, treefil) } else{ warning("No tree provided. Phylogenetic investigation cannot be done.") ps0 <- ps0 } # Remove any OTUs with 0 abundance across all samples ps0 <- prune_taxa(taxa_sums(ps0) > 0, ps0) ps0 <- prune_samples(sample_sums(ps0) > 50, ps0) # Return phyloseq object return(ps0) })
# Reactive values to store imported data pseq_meta <- reactive({ # Err if null if(is.null(input$upload_mapping)){ return(NULL) } # Show message showNotification("Importing data and filtering out samples with less than 1000 reads...", type = "message") # Run import function once all files are uploaded pseq_metadf <- read.csv(input$upload_mapping$datapath, row.names=1,check.names=FALSE) # Return meta object return(pseq_metadf) })
Column {data-width=1000} {.tabset .tabset-fade}
Library size distribution
output$LibSizeHist <- renderPlot({ validate(need(expr = !is.null(pseq()), message = "Please upload data...")) ps0<- pseq() # Plot sample depths SeqDepth = colSums(otu_table(ps0)) sample_data(ps0)$SeqDepth = SeqDepth meta_df <- meta(ps0) p <- ggpubr::gghistogram(meta_df, "SeqDepth") + ggpubr::rotate_x_text() return(p) }, width = 400, height = 400) fillPage(plotOutput("LibSizeHist"))
OTU counts
fluidPage( titlePanel(""), sidebarLayout( sidebarPanel(width = 3, numericInput( "PrevalenceOTU", label = "Prevalence: Number of samples OTU must be present 0.05 means 5%",0, min = 0, max = 0.9), numericInput( "DetectionOTU", label = "Detection: Minimum count of OTUs",0, min = 0, max = 1000) ), output$OTUHist_1 <- renderPlot({ validate(need(expr = !is.null(pseq()), message = "Please upload data...")) ps0 <- pseq() ps.fil <- core(ps0, detection = input$DetectionOTU, prevalence = input$PrevalenceOTU) no_taxa <- phyloseq::ntaxa(ps.fil) Text <- paste0("Numer of OTU/Taxa after filtering: ", no_taxa) # Plot sample depths px.log <- microbiome::transform(ps.fil, "log10") px.clr <- microbiome::transform(ps.fil, "clr") px.hel <- microbiome::transform(ps.fil, "hellinger") p1 <- ggplot2::qplot(rowSums(otu_table(ps0))) + ggplot2::theme_bw() + ggplot2::xlab("Raw counts") + ggtitle("A", subtitle =Text ) p2 <- ggplot2::qplot(rowSums(otu_table(px.log))) + ggplot2::theme_bw() + ggplot2::xlab("log10(1+x) counts") + ggtitle("B", subtitle =Text ) p3 <- ggplot2::qplot(rowSums(otu_table(px.clr))) + ggplot2::theme_bw() + ggplot2::xlab("clr counts") + ggtitle("C", subtitle =Text ) p4 <- ggplot2::qplot(rowSums(otu_table(px.hel))) + ggplot2::theme_bw() + ggplot2::xlab("Hellinger counts") + ggtitle("D", subtitle =Text ) #t <- gridExtra::textGrob(paste0("Numer of OTU/Taxa after filtering: ", no_taxa) ) return(gridExtra::grid.arrange(p1,p2,p3,p4, nrow= 2, ncol= 2)) mainPanel(plotOutput("OTUHist_1")) }, width = 600, height = 600) ))
Check different transformations and distribution of OTU counts in overall data.
Note down the minimum counts and prevalance used.
These parameters can be used to filter potentailly spurious OTUs/taxa before doing differential abundance analysis using:
1. DESeq2.
2. edgeR.
3. limma.
Note: The filtered data is not used for next analysis. For each of the tabs, raw data is used and individual filtering as specifed.
OTU counts varied taxa
# Select taxonomic level selectInput( inputId = "OtuTaxaDist", label = "Choose taxa level", choices = list( "Domain" = "Domain", "Phylum" = "Phylum", "Class" = "Class", "Order" = "Order", "Family" = "Family", "Genus" = "Genus" ), selected = "Phylum") output$OTUtaxa <- renderPlotly({ validate(need(expr = !is.null(pseq()), message = "Loading data...")) ps0 <- pseq() taxasums = rowSums(otu_table(ps0)) taxatable <- as.data.frame.matrix(tax_table(ps0)) tax_plot2 <- ggplotly(tax_plot1 <- ggplot(taxatable, aes(x = taxasums, color = taxatable[, input$OtuTaxaDist])) + geom_line(size = 1.5, stat = "density") + xlab("Log10 OTU Counts") + theme_bw() + scale_x_log10() + guides(color=guide_legend(title=input$OtuTaxaDist, ncol = 3)) + theme(legend.position="bottom"), width = 800, inline = T) return(tax_plot2) }) fluidPage(plotlyOutput("OTUtaxa"))
Summarise phyloseq
fillPage(renderPrint({ validate(need(expr = !is.null(pseq()), message = "Loading data...")) ps0 <- pseq() #ps_meta <- as.data.table(meta(ps0)) # Plot sample depths smp <- summarize_phyloseq(ps0) return(smp) }))
Taxa Prevalence
# Select taxonomic level selectInput( inputId = "PrevalencePlot", label = "Choose taxa level", choices = list( "Domain" = "Domain", "Phylum" = "Phylum", "Class" = "Class", "Order" = "Order", "Family" = "Family" ), selected = "Phylum" ) #Prevalence plot across all samples output$PrevalencePlot <- renderPlot({ # Err if null validate(need(expr = !is.null(pseq()), message = "Loading data...")) ps0 <- pseq() tax_prev_plot <- plot_taxa_prevalence(ps0, input$PrevalencePlot) message("Generating prevalence plot") tax_prev_plot <- tax_prev_plot + theme(legend.position = "bottom") return(tax_prev_plot) }, width = 1200, height = 800) fluidPage(plotOutput("PrevalencePlot"))
Taxa barplot
fluidPage(titlePanel(""), sidebarLayout( sidebarPanel( width = 2, column = 5, ## plot variable h5('Please select the plotting parameters'), selectInput( inputId = "TaxBarplot", label = "Choose taxa level", choices = list( "Domain" = "Domain", "Phylum" = "Phylum", "Class" = "Class", "Order" = "Order", "Family" = "Family", "Genus" = "Genus" ), selected = "Phylum" ), selectInput( "taxbarcolors", label = "Choice of palette default Paired", choices = c("Set2", "Paired", "Set1", "Accent", "Dark2", "Set3"), selected = "Paired" ), numericInput( "topOTUBar", label = "Top taxa to be viewed", 5, min = 3, max = 100 ) ), output$boxplot <- renderPlot({ validate(need(expr = !is.null(pseq()), message = "Loading data...")) ps.bar <- pseq() taxic <- as.data.frame(ps.bar@tax_table) if (input$taxbarcolors == "Set3" | input$taxbarcolors == "Paired") { getPalette = colorRampPalette(brewer.pal(12, input$taxbarcolors)) } else if (input$taxbarcolors == "Set2" | input$taxbarcolors == "Dark2" | input$taxbarcolors == "Accent") { getPalette = colorRampPalette(brewer.pal(8, input$taxbarcolors)) } else if (input$taxbarcolors == "Set1") { getPalette = colorRampPalette(brewer.pal(9, input$taxbarcolors)) } otu.df <- as.data.frame(otu_table(ps.bar)) taxic$OTU <- rownames(taxic) colnames(taxic) taxmat <- as.matrix(taxic) new.tax <- tax_table(taxmat) tax_table(ps.bar) <- new.tax guide_italics <- guides(fill = guide_legend( label.theme = element_text( size = 15, face = "italic", colour = "Black", angle = 0 ) )) if (!is.null(ps.bar@phy_tree)) { ps.bar@phy_tree <- NULL } ps.bar.tlev <- aggregate_taxa(ps.bar, level = input$TaxBarplot, top = input$topOTUBar) tax.lev <- as.data.frame.matrix(tax_table(ps.bar.tlev)) colourCount = length(unique(tax.lev$OTU)) ps.bar.rel <- microbiome::transform(ps.bar.tlev, "compositional") plot.composition.relAbun <- plot_composition(ps.bar.rel) + scale_fill_manual(input$TaxBarplot, values = getPalette(colourCount)) + theme_bw() + theme(axis.text.x = element_text(angle = 90)) + ggtitle(paste0( "Relative abundance of top ", input$topOTUBar, sep = " ", input$TaxBarplot )) + guide_italics + theme(legend.title = element_text(size = 18)) + theme(legend.position = "top") return(plot.composition.relAbun) mainPanel(plotOutput("boxplot")) }, width = 1000, height = 600) ))
Taxa boxplot
fluidPage( titlePanel(""), sidebarLayout( sidebarPanel(width = 2, column = 4, ## plot variable h5('Please type the variable of choice from mapping file'), p('No numeric values, only categorical'), textInput('Sam_var2', 'Variable (eg. SampleType)'), tags$div(id = 'placeholder1'), hr(), selectInput( inputId = "Taxplot", label = "Choose taxa level", choices = list( "Domain" = "Domain", "Phylum" = "Phylum", "Class" = "Class", "Order" = "Order", "Family" = "Family", "Genus" = "Genus", "Species" = "Species"), selected = "Phylum"), selectInput( "taxcolors", label = "Choice of palette default Paired", choices = c("Set2", "Paired", "Set1", "Accent", "Dark2", "Set3"), selected = "Paired"), numericInput( "topOTU", label = "Top taxa to be viewed",5, min = 0, max = 50) ), output$boxplot <- renderPlot({ # validate( # need(expr = !is.null(pseq()), message = "Loading data...") # ) validate( need(expr = input$Sam_var2, message = "Type the variable!!") ) ps.box <- pseq() #req(input$Sam_var2) ps.box.rel <- microbiome::transform(ps.box, "compositional") pn1 <- plot_taxa_boxplot(ps.box.rel, taxonomic.level = input$Taxplot, top.otu = input$topOTU, VariableA = input$Sam_var2, title = "Relative abundance ", color = input$taxcolors) return(pn1 + ggpubr::rremove("grid")) mainPanel(plotOutput("boxplot")) }, width = 1000, height = 600)))
Taxa Heatmap
fluidPage( titlePanel(""), sidebarLayout( sidebarPanel(width = 2, column = 4, ## link the JS file tags$head(tags$script(src="script.js")), ## link the CSS file tags$head(tags$link(rel="stylesheet", type="text/css", href="style.css")), ## plot variable h5('Please type the variable of choice from mapping file'), p('No numeric values, only categorical'), textInput('txt1', 'Variable (eg. SampleType)'), tags$div(id = 'placeholder1'), hr(), selectInput( "transformation", label = "Choice of transformation default log10", choices = c("log10", "compositional", "clr", "Z-OTU"), selected = "log10"), selectInput( "heatcolors", label = "Choice of palette default GnBu", choices = c("RdYlBu", "Blues", "Greys", "GnBu", "RdGy", "Spectral"), selected = "GnBu"), numericInput( "SubsetOTU", label = "Number of top OTUs to visualize (max.150)",10, min = 5, max = 150) ), output$HeatMap <- renderPlot({ validate( need(expr = !is.null(pseq()), message = "Please upload data...") ) validate( need(expr = input$txt1, message = "Type the variable!!") ) ps0<- pseq() #ps0.f <- microbiomeutilities::format_to_besthit(pseq()) mainPanel(plot_taxa_heatmap(ps0, subset.top = input$SubsetOTU, VariableA = input$txt1, heatcolors = brewer.pal(100, input$heatcolors), transformation = input$transformation)) }, width = 1000, height = 700)))
Diversity (Non-Phylogenetic)
fluidPage( titlePanel("Non-Phylogenetic diversity"), sidebarLayout( sidebarPanel(width = 2, column = 4, ## link the JS file tags$head(tags$script(src="script.js")), ## link the CSS file tags$head(tags$link(rel="stylesheet", type="text/css", href="style.css")), ## plot variable h5('Please type the variable of choice from mapping file'), p('No numeric values, only categorical'), textInput('Sam_var3', 'Variable (eg. SampleType)'), tags$div(id = 'placeholder1'), hr(), selectInput( "typeDiv", label = "Choice of diversities measure default: diversities", choices = c("diversities", "dominance", "evenness"), selected = "diversities"), selectInput( "divcolors", label = "Choice of palette default: Set1", choices = c("Set2", "Paired", "Set1", "Accent", "Dark2", "Set3"), selected = "Set1"), selectInput( "figtype", label = "Choice of plot default: stripchart", choices = c("stripchart", "boxplot", "violin"), selected = "stripchart") ), output$diveplot <- renderPlot({ validate( need(expr = !is.null(pseq()), message = "Loading data...") ) validate( need(expr = input$Sam_var3, message = "Type the variable!!") ) ps0<- pseq() p <- plot_alpha_diversities(ps0, type = input$typeDiv, index.val = "all", plot.type = input$figtype, variableA = input$Sam_var3, palette = input$divcolors) + ggpubr::rremove("grid") return(p) mainPanel(plotOutput("diveplot")) }, width = 900, height = 600) ))
Phylogenetic diversity
fluidPage( titlePanel("Phylogenetic diversity"), sidebarLayout( sidebarPanel(width = 2, column = 3, ## link the JS file tags$head(tags$script(src="script.js")), ## link the CSS file tags$head(tags$link(rel="stylesheet", type="text/css", href="style.css")), ## plot variable h5('Please type the variable of choice from mapping file'), p('No numeric values, only categorical'), textInput('Sam_var4', 'Variable (eg. SampleType)'), tags$div(id = 'placeholder1'), hr(), selectInput( "pdivcolors", label = "Choice of palette default Set2", choices = c("Set2", "Paired", "Set1", "Accent", "Dark2", "Set3"), selected = "Set2"), selectInput( "figPtype", label = "Choice of plot default: stripchart", choices = c("stripchart", "boxplot", "violin"), selected = "stripchart") ), output$diveplot <- renderPlot({ av = function(x){ if( isTRUE(all.equal(x, "")) | isTRUE(all.equal(x, "NULL")) ){ return(NULL) } return(x) } validate( need(expr = !is.null(pseq()), message = "Please upload data...") ) validate( need(expr = !is.null(pseq()@phy_tree), message = "No phylogenetic tree!") ) ps0<- pseq() if(!is.null(ps0@phy_tree)){ otu_table_ps3 <- as.data.frame(ps0@otu_table) treefil2 = read_tree(input$upload_tree$datapath) df.pd <- pd(t(otu_table_ps3), treefil2,include.root=F) metadata_table_ps3 <- meta(ps0) sample_data(ps0)$Phyogenetic_diversity <- df.pd$PD #metadata_table_ps3$Phyogenetic_diversity <- df.pd$PD mps0 <- phy_to_ldf(ps0, transform.counts = NULL) colnames(mps0)[colnames(mps0) == input$Sam_var4] <- 'Variable_alias' # mps0 <- reshape2::melt(metadata_table_ps3) if (input$figPtype == "boxplot") { plot.pd <- ggpubr::ggboxplot(mps0, "Variable_alias", "Phyogenetic_diversity", palette = input$pdivcolors, fill = "Variable_alias", legend = "right", add = "mean_sd" )} else if (input$figPtype == "stripchart") { plot.pd <- ggpubr::ggstripchart(mps0, "Variable_alias", "Phyogenetic_diversity", palette = input$pdivcolors, color = "Variable_alias", legend = "right", add = "mean_sd", shape = 21) } else if (input$figPtype == "violin") { plot.pd <- ggpubr::ggviolin(mps0, "Variable_alias", "Phyogenetic_diversity", palette = input$pdivcolors, fill = "Variable_alias", legend = "right", add = "mean_sd") } } else{ # Show message showNotification("No tree provided. Phylogenetic diversity, unifrac distances cannot be calculated.", duration = NULL, type = "error") plot.pd <- "No tree provided. Phylogenetic diversity, unifrac distances cannot be calculated" } return(plot.pd + ggpubr::rremove("grid")) mainPanel(plotOutput("diveplot")) }, width = 900, height = 600)))
Ordination
fluidPage( titlePanel("Beta diversity"), sidebarLayout( sidebarPanel(width = 2, column = 3, ## link the JS file tags$head(tags$script(src="script.js")), ## link the CSS file tags$head(tags$link(rel="stylesheet", type="text/css", href="style.css")), ## plot variable h5('Please type the variable of choice from mapping file'), textInput('txt3', 'Variable (eg. SampleType)'), tags$div(id = 'placeholder1'), hr(), selectInput( "Distance", label = "Choice of distance method default: Bray-Curtis", choices = c("wunifrac", "unifrac", "bray"), selected = "bray"), selectInput( "Varcolors", label = "Choice of palette default Set1", choices = c("Set2", "Paired", "Set1", "Accent", "Dark2", "Set3"), selected = "Set1"), numericInput( "Prevalence", label = "Prevalence: Number of samples OTU must be present 0.05 means 5%",0, min = 0, max = 0.9), numericInput( "Detection", label = "Detection: Minimum relative abundance of OTUs",0, min = 0, max = 0.1) ), output$Ordination <- renderPlot({ av = function(x) { if (isTRUE(all.equal(x, "")) | isTRUE(all.equal(x, "NULL"))) { return(NULL) } return(x) } validate(need(expr = !is.null(pseq()), message = "Please upload data...")) ps0 <- pseq() ps0a.rel <- microbiome::transform(ps0, "compositional") ps0a.rel <- core(ps0a.rel, detection = input$Detection, prevalence = input$Prevalence) #ps0.f <- microbiomeutilities::format_to_besthit(pseq()) if (input$Distance == "wunifrac") { ord <- ordinate(ps0a.rel, method = "PCoA", input$Distance) ordiplot <- plot_ordination(ps0a.rel, ord, color = av(input$txt3)) + ggtitle("Weighted Unifrac Compositional") + theme_bw() + geom_point(size = 3) + scale_color_brewer(palette = input$Varcolors) + ggpubr::rremove("grid") return(ordiplot) } else if (input$Distance == "unifrac") { ord <- ordinate(ps0a.rel, method = "PCoA", input$Distance, weighted = F) ordiplot <- plot_ordination(ps0a.rel, ord, color = av(input$txt3)) + ggtitle("UnWeighted Unifrac Compositional") + theme_bw() + geom_point(size = 3) + scale_color_brewer(palette = input$Varcolors) + ggpubr::rremove("grid") return(ordiplot) } else if (input$Distance == "bray") { ord <- ordinate(ps0a.rel, method = "PCoA", input$Distance) ordiplot <- plot_ordination(ps0a.rel, ord, color = av(input$txt3)) + ggtitle("Bray-Curtis Compositional") + theme_bw() + geom_point(size = 3) + scale_color_brewer(palette = input$Varcolors) + ggpubr::rremove("grid") return(ordiplot) } mainPanel(plotOutput("diveplot")) }, width = 800, height = 600) ))
Core microbiota
fluidPage( titlePanel(""), sidebarLayout( sidebarPanel(width = 2,column = 3, numericInput( "PrevalenceCore", label = "Prevalence: Number of samples OTU must be present 0.05 means 5%",0.5, min = 0.25, max = 0.99), h5('Minimum detection of taxa is set at 0.0001% see x-axis for range') ), output$core_plots <- renderPlot({ av = function(x){ if( isTRUE(all.equal(x, "")) | isTRUE(all.equal(x, "NULL")) ){ return(NULL) } return(x) } validate( need(expr = !is.null(pseq()), message = "Please upload data...") ) ps0<- pseq() #ps0.sub <- phyloseq::subset_samples(ps0, av(input$txt5) == input$txt6) ps0a.rel <- microbiome::transform(ps0, "compositional") ps0a.rel <- format_to_besthit(ps0a.rel) prevalences <- seq(.05, 1, .05) detections <- 10^seq(log10(1e-4), log10(.2), length = 10) core_plot <- plot_core(ps0a.rel, plot.type = "heatmap", prevalences = prevalences, detections = detections, colours = rev(brewer.pal(5, "Spectral")), min.prevalence = input$PrevalenceCore, horizontal = F) return(core_plot + ggtitle(paste0("Core taxa at ", (input$PrevalenceCore*100), " % prevalence and min. abundance 0.0001%"))) mainPanel(plotOutput("core_plots")) }, width = 800, height = 600) ))
Data tables
Tables {data-width=400} {.tabset .tabset-fade}
Metadata
get_meta_data = reactive({ ps0 <- pseq() metadfx <- meta(ps0) return(metadfx) }) renderDataTable({ get_meta_data() })
Diversities table
output$DiversitiesTab = DT::renderDataTable({ # Err if null validate( need(expr = !is.null(pseq()), message = "Please upload data...") ) ps1<- pseq() ps.diver <- diversities(ps1) return(ps.diver) }, server = F, extensions = 'Buttons', options = list(lengthChange = TRUE, dom = 'Bfrtip', buttons = c('copy', 'print'), autoWidth = F, pageLength = 15, searchHighlight = TRUE, scrollX = F) ) fillPage(DT::dataTableOutput('DiversitiesTab'))
Dominance table
output$DominanceTab = DT::renderDataTable({ # Err if null validate( need(expr = !is.null(pseq()), message = "Please upload data...") ) ps1<- pseq() ps.diver <- dominance(ps1) return(ps.diver) }, server = F, extensions = 'Buttons', options = list(lengthChange = TRUE, dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print'), autoWidth = F, pageLength = 15, searchHighlight = TRUE, scrollX = F) ) fillPage(DT::dataTableOutput('DominanceTab'))
Evenness table
output$EvennessTab = DT::renderDataTable({ # Err if null validate( need(expr = !is.null(pseq()), message = "This will take time depending on your data size...") ) showNotification("This will take time depending on your data size...", type = "message") ps1<- pseq() ps.diver <- evenness(ps1) return(ps.diver) }, server = F, extensions = 'Buttons', options = list(lengthChange = TRUE, dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print'), autoWidth = F, pageLength = 15, searchHighlight = TRUE, scrollX = F) ) fillPage(DT::dataTableOutput('EvennessTab'))
Phylogenetic diversity table
output$PhylogenyTab = DT::renderDataTable({ validate( need(expr = !is.null(pseq()), message = "Please upload data...") ) ps0<- pseq() if(!is.null(ps0@phy_tree)){ otu_table_ps3 <- as.data.frame(ps0@otu_table) treefil2 = read_tree(input$upload_tree$datapath) df.pd <- pd(t(otu_table_ps3), treefil2,include.root=F) return(df.pd) } else{ # Show message showNotification("No tree provided. Phylogenetic diversity, unifrac distances cannot be calculated.", duration = NULL, type = "error") plot.pd <- "No tree provided. Phylogenetic diversity, unifrac distances cannot be calculated" } }, server = F, extensions = 'Buttons', options = list(lengthChange = TRUE, dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print'), autoWidth = F, pageLength = 15, searchHighlight = TRUE, scrollX = F)) fluidPage(DT::dataTableOutput('PhylogenyTab'))
Sample-Features {data-width=150}
Number of samples
output$sampleN <- renderValueBox({ # Err if NULL if(is.null(pseq())){ return(valueBox(0, caption = "Number of samples", icon = "fa-cubes", color = "black")) } valueBox(nsamples(pseq()), caption = "Number of samples", icon = "fa-cubes", color = "black") }) valueBoxOutput(outputId = "sampleN", width = 3)
Number of OTU's
output$taxaN <- renderValueBox({ # Err if NULL if(is.null(pseq())){ return(valueBox(0, caption = "Number of taxa", icon = "fa-users", color = "black")) } valueBox(ntaxa(pseq()), caption = "Number of taxa", icon = "fa-cubes", color = "black") }) valueBoxOutput(outputId = "taxaN", width = "auto")
Number of samples
output$totalDepth <- renderValueBox({ # Err if NULL if(is.null(pseq())){ return(valueBox(0, caption = "Total sampling depth", icon = "fa-bar-chart", color = "black")) } valueBox(sum(sample_sums(pseq())), caption = "Total sampling depth", icon = "fa-info-circle", color = "black") }) valueBoxOutput(outputId = "totalDepth")
Mean sampling depth
output$meanDepth <- renderValueBox({ # Err if NULL if(is.null(pseq())){ return(valueBox(0, caption = "Mean sampling depth", icon = "fa-info-circle", color = "black")) } valueBox(round(mean(sample_sums(pseq())), 0), caption = "Mean sampling depth", icon = "fa-info-circle",color = "black" ) }) valueBoxOutput(outputId = "meanDepth")
Median sampling depth
output$medianDepth <- renderValueBox({ # Err if NULL if(is.null(pseq())){ return(valueBox(0, caption = "Median sampling depth", icon = "fa-info-circle", color = "black")) } valueBox(round(median(sample_sums(pseq())), 0), caption = "Median sampling depth", icon = "fa-info-circle", color = "black") }) valueBoxOutput(outputId = "medianDepth")
Range of sampling depth
output$rangeDepth <- renderValueBox({ # Err if NULL if(is.null(pseq())){ return(valueBox(0, caption = "Range of sampling depth", icon = "fa-info-circle", color = "black")) } range <- paste0(min(sample_sums(pseq())), " - ", max(sample_sums(pseq()))) valueBox(range, caption = "Range of sampling depth", icon = "fa-info-circle", color = "black") }) valueBoxOutput(outputId = "rangeDepth")
Standard deviation of sampling depth
output$sdDepth <- renderValueBox({ # Err if NULL if(is.null(pseq())){ return(valueBox(0, caption = "Standard deviation of sampling depth", icon = "fa-info-circle", color = "black")) } sd <- round(sd(sample_sums(pseq())), digits = 1) valueBox(sd, caption = "Standard deviation of sampling depth", icon = "fa-info-circle", color = "black") }) valueBoxOutput(outputId = "sdDepth")
Singletons
output$sington <- renderValueBox({ # Err if NULL if(is.null(pseq())){ return(valueBox(0, caption = "Singletons in data", icon = "fa-info-circle", color = "black")) } singletons <- round(length(taxa_sums(pseq())[taxa_sums(pseq()) <= 1]), digits = 1) valueBox(singletons, caption = "Singletons in data", icon = "fa-info-circle", color = "black") }) valueBoxOutput(outputId = "sington")
Sparsity
output$spar <- renderValueBox({ # Err if NULL if(is.null(pseq())){ return(valueBox(0, caption = "Sparsity", icon = "fa-info-circle", color = "black")) } sparsit <- round(length(which(abundances(pseq()) == 0))/length(abundances(pseq())), digits = 1) valueBox(sparsit, caption = "Sparsity", icon = "fa-info-circle", color = "black") }) valueBoxOutput(outputId = "spar")
LICENSE
output$LICENSE <- renderText({ print("GNU General Public License v3.0") }) textOutput(outputId = "LICENSE")
Exploratory
Exploratory-tabs {data-width=400} {.tabset .tabset-fade}
Taxa features
fluidPage(titlePanel(""), sidebarLayout( sidebarPanel( width = 2, column = 4, ## plot variable h5('Please type the variable of choice from mapping file'), p('No numeric values, only categorical'), textInput('Sam_varFeat', 'Variable (eg. SampleType)'), tags$div(id = 'placeholder1'), hr(), sliderInput( "DetectionFeat", label = "Detection: Minimum relative abundance of OTUs", min = 0.0001, max = 0.1, value = 0.001, step = 0.005 ), sliderInput( "bins", "Bins:", min = 0, max = 500, value = 100 ) ), output$tableplot <- renderPlot({ validate(need(expr = !is.null(pseq()), message = "Please upload data...")) ps1 <- pseq() ps1.rel <- microbiome::transform(ps1, "compositional") ps1.rel.filter <- microbiome::core(ps1.rel, detection = input$DetectionFeat, prevalence = 0) ps1.df <- phy_to_ldf(ps1.rel.filter, transform.counts = NULL) ps1.df2 <- ps1.df[, c("Phylum", "Order", "Class", "Family" , "Abundance", input$Sam_varFeat)] return( tabplot::tableplot( ps1.df2, fontsize = 10, legend.lines = 20, nBins = input$bins, sampleBinSize = 1e2, sortCol = Abundance, relative = TRUE, scales = "lin" ) ) mainPanel(plotOutput("tableplot")) }, width = 1200, height = 800) ))
Insights/Notes
Microbiome analytics tool
Pre-processing is an essential step and in many instance is not given importance, as we biologists do not take these steps as seriously as it should be taken. To avoid this there are some common approaches to look at the data. On purpose this webapp does not allow for downloading high-quality images. It is just to have a quick look at your data and then play around in R using Rmarkdown to make you codes easily traceble and results reproducible.
This shiny app is inspired from from other shiny tools such as shiny-phyloseq and btools-shiny.
This first-hand look at your data will aid in making different choices for more thorough and in-dept investagion.
You will need to install microbiomeutilities R package from github repo.
The main tools used here are Phyloseq and microbiome.
Check library size distribution
Check the library sizes for each of the samples within the main variable. Do you think it is okay or there is difference in library sizes and will this affect you downstream statistical analysis?
Check OTU/Taxa count distribution
If the OTU counts are higly skewed in distribution try normalising them.
Check OTU/Taxa count distribution
Which taxa are distributied differently. Check if they are of any interest in the ecosystem being investigated.
Summary of Phyloseq
Look if there are singletons, sparsity in your data.
Taxa prevalence
Which taxa are prevalant and what abundances. Are there OTUs with low abundace and low prevalance?
Taxa barplot and Taxa boxplot
Gives a firsthand view of the microbial composition.
Taxa heatmap
At OTU level, see if there is any difference in top N otus between variables of interest.
Alpha diversity
This includes the non-phylogenetic metrices. For more on this check Microbiome:Diversities
Phylogenetic diversity
For more on this check Picante. Is calculated only if tree is provided in the upload tab.
Beta Diversity (Ordinations)
Here, commonly used metrics for beta diversity are depicted. The default transformation is compositional (relative abundance). Choose filtering option Detections and Prevalance to remove potentially spurious OTUs/Taxa.
For more information:
Waste Not, Want Not: Why Rarefying Microbiome Data Is Inadmissible.
Normalisation and data transformation.
What is Constrained and Unconstrained Ordination.
R packages that work at the backend of this webapp are:
- picante
- data.table
- RColorBrewer
- shinythemes
- DT
- shiny
- reshape
- ggpubr
- tibble
- pheatmap
- Phyloseq
- microbiome
- plyr
- tidyverse R packages
- pheatmap
- Vegan
- flexdashboard
- Biobase
- plotly
- ape
Please let me know if I have missed any of the packages.
Terms and conditions
Microbiome analytics tool
Terms of Use This web application is hosted on a Shinyapps.io server https://www.shinyapps.io/.
If you have concerns about the terms of use for this web hosted application, please run the app locally on your computer.
By using this app you are agreeing to the terms of use as described by Shinnyaps.io: https://www.rstudio.com/about/shinyapps-terms-use/
The authors and maintainers of microbiomeutilities-shiny app will not save your data on our servers. However, as the Shinyapps server is not HIPAA compliant, you must refrain from uploading protected health information or confidential data with this app. You may instead download the code and run the app locally on your private computer and network (see above). We the authors cannot be held responsible for the confidentiality, availability, security, loss, misuse or misappropriation of any data you submit to this application.
One of the Important information from Shinyapps terms of use (https://www.rstudio.com/about/shinyapps-terms-use/): "If you choose to upload data to an application you are using via the RStudio Service, you acknowledge and agree you are giving certain legal rights to the licensor of the application to process and otherwise use your data. Please carefully review any license terms accompanying any application to which you submit your data for the legal rights which you are giving the application licensor. Further, RStudio does not claim ownership of your data; however, you hereby grant RStudio a worldwide, perpetual, irrevocable, royalty-free, fully paid up, transferable and non-exclusive license, as applicable, to use and copy your data in connection with making your data available to the application (and the licensor of such application) to which you have submitted your data. You acknowledge and agree that even if you remove your data from the RStudio Service, your data may have been downloaded by, and remain accessible to, the licensors of those applications to which you submitted your data. Accordingly, do not submit data which you desire to remain confidential or which you wish to limit the right to access or use. You should never submit to the RStudio Service any data which consists of personally identifiable information, credit card information, or protected health information, as such terms are defined by relevant laws, rules, and regulations. RSTUDIO IS NOT RESPONSIBLE FOR THE CONFIDENTIALITY, AVAILABILITY, SECURITY, LOSS, MISUSE OR MISAPPROPRIATION OF ANY DATA YOU SUBMIT TO THE RSTUDIO SERVICE OR ANY APPLICATION MADE AVAILABLE VIA THE RSTUDIO SERVICE."
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.