real

real_samples.Rmd

@@ -38,6 +38,8 @@ library(tidyverse)
 library(ComplexHeatmap)
 library(phyloseqCompanion)
 library(circlize)
+library(ggVennDiagram)
+library(ggplot2)
 ```
 
 
@@ -49,89 +51,6 @@ options(DT.options = list(pageLength = 10,
                          buttons = c('copy', 'csv', 'excel', 'pdf', 'print')))
 ```
 
-# MEAT
-
-```{bash, eval=F}
-# Genologin
-cd /work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/scripts
-sh dada2_frogs.sh
-```
-
-```{bash, eval=F}
-#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/ITS1/frogs.biom MEAT/ITS1.biom
-#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/ITS2/frogs.biom MEAT/ITS2.biom
-#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/D1D2/frogs.biom MEAT/D1D2.biom
-scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/RPB2/frogs.biom MEAT/RPB2.biom
-
-#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/ITS1/frogs.tsv MEAT/ITS1.tsv
-#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/ITS2/frogs.tsv MEAT/ITS2.tsv
-#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/D1D2/frogs.tsv MEAT/D1D2.tsv
-scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/RPB2/frogs.tsv MEAT/RPB2.tsv
-
-#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/ITS1/multiaff.tsv MEAT/ITS1_multiaff.tsv
-#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/ITS2/multiaff.tsv MEAT/ITS2_multiaff.tsv
-#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/D1D2/multiaff.tsv MEAT/D1D2_multiaff.tsv
-scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/RPB2/multiaff.tsv MEAT/RPB2_multiaff.tsv
-```
-
-## Preparation of final BIOM files {.tabset}
-
-### ITS1
-
-### ITS2
-
-### D1D2
-
-### RPB2
-
-## Analysis
-
-```{r, eval=F}
-if (!file.exists("REAL/ITS1.rds")) {
-  biomfile <- "REAL/ITS1.biom"
-  physeq_its1 <- import_frogs(biomfile, taxMethod = "blast")
-  metadata <- read.table("REAL/metadata.tsv", row.names = 1, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
-  sample_data(physeq_its1) <- metadata
-  sample_data(physeq_its1)$Marker <- "ITS1"
-  saveRDS(physeq_its1,"REAL/ITS1.rds")
-} else {
-  physeq_its1 <- readRDS("REAL/ITS1.rds")
-}
-
-if (!file.exists("REAL/ITS2.rds")) {
-  biomfile <- "REAL/ITS2.biom"
-  physeq_its2 <- import_frogs(biomfile, taxMethod = "blast")
-  metadata <- read.table("REAL/metadata.tsv", row.names = 1, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
-  sample_data(physeq_its2) <- metadata
-  sample_data(physeq_its2)$Marker <- "ITS2"
-  saveRDS(physeq_its2,"REAL/ITS2.rds")
-} else {
-  physeq_its2 <- readRDS("REAL/ITS2.rds")
-}
-
-if (!file.exists("REAL/D1D2.rds")) {
-  biomfile <- "REAL/D1D2.biom"
-  physeq_d1d2 <- import_frogs(biomfile, taxMethod = "blast")
-  metadata <- read.table("REAL/metadata.tsv", row.names = 1, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
-  sample_data(physeq_d1d2) <- metadata
-  sample_data(physeq_d1d2)$Marker <- "D1D2"
-  saveRDS(physeq_d1d2,"REAL/D1D2.rds")
-} else {
-  physeq_d1d2 <- readRDS("REAL/D1D2.rds")
-}
-
-if (!file.exists("REAL/RPB2.rds")) {
-  biomfile <- "REAL/RPB2.biom"
-  physeq_rpb2 <- import_frogs(biomfile, taxMethod = "blast")
-  metadata <- read.table("REAL/metadata.tsv", row.names = 1, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
-  sample_data(physeq_rpb2) <- metadata
-  sample_data(physeq_rpb2)$Marker <- "RPB2"
-  saveRDS(physeq_rpb2,"REAL/RPB2.rds")
-} else {
-  physeq_rpb2 <- readRDS("REAL/RPB2.rds")
-}
-```
-
 
 ```{r, eval=T}
 create_frogs_object <- function(path, marker, type, file, multifile){
@@ -183,40 +102,125 @@ keep_only_interest_species <- function(physeq, true_species){
 }
 ```
 
+# MEAT
 
-```{r, eval=T}
-rpb2 <- create_frogs_object("REAL","RPB2", "DNA", "RPB2.biom", "RPB2_multiaff.tsv")
-tax_table(rpb2)[, colnames(tax_table(rpb2))] <- gsub(tax_table(rpb2)[, colnames(tax_table(rpb2))],     pattern = "[a-z]:", replacement = "")
-d1d2 <- create_frogs_object("REAL","D1D2", "DNA", "D1D2.biom", "D1D2_multiaff.tsv")
-tax_table(d1d2)[, colnames(tax_table(d1d2))] <- gsub(tax_table(d1d2)[, colnames(tax_table(d1d2))],     pattern = "[a-z]:", replacement = "")
-its1 <- create_frogs_object("REAL","ITS1", "DNA", "ITS1.biom", "ITS1_multiaff.tsv")
-its2 <- create_frogs_object("REAL","ITS2", "DNA", "ITS2.biom", "ITS2_multiaff.tsv")
+```{bash, eval=F}
+# Genologin
+cd /work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/scripts
+sh dada2_frogs.sh
+```
 
-all <- merge_phyloseq(its1, its2, d1d2, rpb2)
+```{bash, eval=F}
+scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/ITS1/affiliation.biom REAL_MEAT/ITS1.biom
+scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/ITS2/affiliation.biom REAL_MEAT/ITS2.biom
+scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/D1D2/affiliation.biom REAL_MEAT/D1D2.biom
+scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/RPB2/affiliation.biom REAL_MEAT/RPB2.biom
+
+scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/ITS1/affiliation.tsv REAL_MEAT/ITS1.tsv
+scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/ITS2/affiliation.tsv REAL_MEAT/ITS2.tsv
+scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/D1D2/affiliation.tsv REAL_MEAT/D1D2.tsv
+scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/RPB2/affiliation.tsv REAL_MEAT/RPB2.tsv
+
+#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/ITS1/multiaff.tsv REAL_MEAT/ITS1_multiaff.tsv
+#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/ITS2/multiaff.tsv REAL_MEAT/ITS2_multiaff.tsv
+#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/D1D2/multiaff.tsv REAL_MEAT/D1D2_multiaff.tsv
+#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/REAL/results/DADA2_FROGS/RPB2/multiaff.tsv REAL_MEAT/RPB2_multiaff.tsv
+```
 
-sample_data(all)$Marker <- factor(sample_data(all)$Marker, levels=c("ITS1","ITS2","D1D2","RPB2"))
-Method <- get_variable(all, "method")
+## Analysis of raw BIOM files
 
-saveRDS(all, "REAL/real_all.rds")
+```{r, eval=T}
+if (!file.exists("REAL_MEAT/ITS1.rds")) {
+  biomfile <- "REAL_MEAT/ITS1.biom"
+  physeq_its1 <- import_frogs(biomfile, taxMethod = "blast")
+  tax_table(physeq_its1) <- gsub(tax_table(physeq_its1),     pattern = "[a-z]__", replacement = "")
+  metadata <- read.table("REAL_MEAT/metadata.tsv", row.names = 1, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
+  sample_data(physeq_its1) <- metadata
+  sample_data(physeq_its1)$Marker <- "ITS1"
+  sample_names(physeq_its1) <- glue::glue(paste("{sample_names(physeq_its1)}","ITS1", sep="_"))
+  saveRDS(physeq_its1,"REAL_MEAT/ITS1.rds")
+} else {
+  physeq_its1 <- readRDS("REAL_MEAT/ITS1.rds")
+}
 
+if (!file.exists("REAL_MEAT/ITS2.rds")) {
+  biomfile <- "REAL_MEAT/ITS2.biom"
+  physeq_its2 <- import_frogs(biomfile, taxMethod = "blast")
+  tax_table(physeq_its2) <- gsub(tax_table(physeq_its2),     pattern = "[a-z]__", replacement = "")
+  metadata <- read.table("REAL_MEAT/metadata.tsv", row.names = 1, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
+  sample_data(physeq_its2) <- metadata
+  sample_data(physeq_its2)$Marker <- "ITS2"
+  sample_names(physeq_its2) <- glue::glue(paste("{sample_names(physeq_its2)}","ITS2", sep="_"))
+  saveRDS(physeq_its2,"REAL_MEAT/ITS2.rds")
+} else {
+  physeq_its2 <- readRDS("REAL_MEAT/ITS2.rds")
+}
 
+if (!file.exists("REAL_MEAT/D1D2.rds")) {
+  biomfile <- "REAL_MEAT/D1D2.biom"
+  physeq_d1d2 <- import_frogs(biomfile, taxMethod = "blast")
+  tax_table(physeq_d1d2) <- gsub(tax_table(physeq_d1d2),     pattern = "[a-z]__", replacement = "")
+  tax_table(physeq_d1d2) <- gsub(tax_table(physeq_d1d2),     pattern = "Eukaryota", replacement = "Fungi")
+  metadata <- read.table("REAL_MEAT/metadata.tsv", row.names = 1, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
+  sample_data(physeq_d1d2) <- metadata
+  sample_data(physeq_d1d2)$Marker <- "D1D2"
+  sample_names(physeq_d1d2) <- glue::glue(paste("{sample_names(physeq_d1d2)}","D1D2", sep="_"))
+  saveRDS(physeq_d1d2,"REAL_MEAT/D1D2.rds")
+} else {
+  physeq_d1d2 <- readRDS("REAL_MEAT/D1D2.rds")
+}
 
+if (!file.exists("REAL_MEAT/RPB2.rds")) {
+  biomfile <- "REAL_MEAT/RPB2.biom"
+  physeq_rpb2 <- import_frogs(biomfile, taxMethod = "blast")
+  tax_table(physeq_rpb2) <- gsub(tax_table(physeq_rpb2),     pattern = "[a-z]__", replacement = "")
+  metadata <- read.table("REAL_MEAT/metadata.tsv", row.names = 1, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
+  sample_data(physeq_rpb2) <- metadata
+  sample_data(physeq_rpb2)$Marker <- "RPB2"
+  sample_names(physeq_rpb2) <- glue::glue(paste("{sample_names(physeq_rpb2)}","RPB2", sep="_"))
+  saveRDS(physeq_rpb2,"REAL_MEAT/RPB2.rds")
+} else {
+  physeq_rpb2 <- readRDS("REAL_MEAT/RPB2.rds")
+}
 ```
 
-
 ```{r, eval=T}
-p <- plot_bar(all)
+all <- merge_phyloseq(physeq_its1, physeq_its2, physeq_d1d2, physeq_rpb2)
+
+sample_data(all)$Marker <- factor(sample_data(all)$Marker, levels=c("ITS1","ITS2","D1D2","RPB2"))
+
+p <- plot_composition(physeq = all, taxaRank1 = "Kingdom", taxaSet1 = "Fungi", taxaRank2 = "Species", numberOfTaxa = 22, x = "Sample")
 p + facet_grid(". ~ Marker", scales = "free_x", space = "free")
+
 ```
 
 ```{r, eval=T}
-p <- plot_composition(physeq = all, taxaRank1 = "Kingdom", taxaSet1 = "Fungi", taxaRank2 = "Genus", numberOfTaxa = 22, x = "Sample")
-p + facet_grid(". ~ Marker", scales = "free_x", space = "free")
+p <- plot_richness(physeq = all, x = "Marker", color = "Marker", shape = NULL, title = "Alpha diversity graphics", measures = c("Observed", "Chao1", "ACE", "Shannon", "Simpson", "InvSimpson", "Fisher"))
+p + geom_boxplot() + NULL
+```
+
+
+```{r, eval=T}
+beta.dist <- distance(all, method = "bray")
+ord <- ordinate(all, method = "MDS", distance = beta.dist)
+p <- plot_ordination(physeq = all, ordination = ord, type = "samples", axes = c(1, 2), color = "Sample", shape = "Marker", label = NULL, title = "Samples ordination graphic")
+p <- p + stat_ellipse(aes_string(group = "Marker"))
+p + theme_bw()
 ```
 
 ```{r, eval=T}
-p <- plot_composition(physeq = all, taxaRank1 = "Kingdom", taxaSet1 = "Fungi", taxaRank2 = "Species", numberOfTaxa = 22, x = "Sample")
-p + facet_grid(". ~ Marker", scales = "free_x", space = "free")
+x <- list(
+  ITS1 = as_tibble(tax_table(physeq_its1))$Genus %>% as.vector(), 
+  ITS2 = as_tibble(tax_table(physeq_its2))$Genus %>% as.vector(),
+  D1D2 = as_tibble(tax_table(physeq_d1d2))$Genus %>% as.vector(),
+  RPB2 = as_tibble(tax_table(physeq_rpb2))$Genus %>% as.vector()
+  )
+ggVennDiagram(x,label_alpha = 0) + 
+  scale_fill_gradient(low = "#F4FAFE", high = "#4981BF") 
+  
+
+
+
 ```
 
 # CHEESE
@@ -228,34 +232,113 @@ sh dada2_frogs.sh
 ```
 
 ```{bash, eval=F}
-#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/ITS1/frogs.biom CHEESE/ITS1.biom
-#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/ITS2/frogs.biom CHEESE/ITS2.biom
-#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/D1D2/frogs.biom CHEESE/D1D2.biom
-scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/RPB2/frogs.biom CHEESE/RPB2.biom
-
-#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/ITS1/frogs.tsv CHEESE/ITS1.tsv
-#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/ITS2/frogs.tsv CHEESE/ITS2.tsv
-#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/D1D2/frogs.tsv CHEESE/D1D2.tsv
-scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/RPB2/frogs.tsv CHEESE/RPB2.tsv
-
-#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/ITS1/multiaff.tsv CHEESE/ITS1_multiaff.tsv
-#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/ITS2/multiaff.tsv CHEESE/ITS2_multiaff.tsv
-#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/D1D2/multiaff.tsv CHEESE/D1D2_multiaff.tsv
-scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/RPB2/multiaff.tsv CHEESE/RPB2_multiaff.tsv
+scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/ITS1/affiliation.biom REAL_CHEESE/ITS1.biom
+scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/ITS2/frogs.biom REAL_CHEESE/ITS2.biom
+scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/D1D2/affiliation.biom REAL_CHEESE/D1D2.biom
+scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/RPB2/affiliation.biom REAL_CHEESE/RPB2.biom
+
+scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/ITS1/affiliation.tsv REAL_CHEESE/ITS1.tsv
+scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/ITS2/affiliation.tsv REAL_CHEESE/ITS2.tsv
+scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/D1D2/affiliation.tsv REAL_CHEESE/D1D2.tsv
+scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/RPB2/affiliation.tsv REAL_CHEESE/RPB2.tsv
+
+#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/ITS1/multiaff.tsv REAL_CHEESE/ITS1_multiaff.tsv
+#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/ITS2/multiaff.tsv REAL_CHEESE/ITS2_multiaff.tsv
+#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/D1D2/multiaff.tsv REAL_CHEESE/D1D2_multiaff.tsv
+#scp orue@genologin.toulouse.inra.fr:/work/frogsfungi/FungiPubli/Real_mock/METABARFOOD/CHEESE_REAL/results/DADA2_FROGS/RPB2/multiaff.tsv REAL_CHEESE/RPB2_multiaff.tsv
+```
+
+
+## Analysis of raw BIOM files {.tabset}
+
+```{r, eval=T}
+if (!file.exists("REAL_CHEESE/ITS1.rds")) {
+  biomfile <- "REAL_CHEESE/ITS1.biom"
+  physeq_its1 <- import_frogs(biomfile, taxMethod = "blast")
+  tax_table(physeq_its1) <- gsub(tax_table(physeq_its1),     pattern = "[a-z]__", replacement = "")
+  metadata <- read.table("REAL_CHEESE/metadata.tsv", row.names = 1, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
+  sample_data(physeq_its1) <- metadata
+  sample_data(physeq_its1)$Marker <- "ITS1"
+  sample_names(physeq_its1) <- glue::glue(paste("{sample_names(physeq_its1)}","ITS1", sep="_"))
+  saveRDS(physeq_its1,"REAL_CHEESE/ITS1.rds")
+} else {
+  physeq_its1 <- readRDS("REAL_CHEESE/ITS1.rds")
+}
+
+if (!file.exists("REAL_CHEESE/ITS2.rds")) {
+  biomfile <- "REAL_CHEESE/ITS2.biom"
+  physeq_its2 <- import_frogs(biomfile, taxMethod = "blast")
+  tax_table(physeq_its2) <- gsub(tax_table(physeq_its2),     pattern = "[a-z]__", replacement = "")
+  metadata <- read.table("REAL_CHEESE/metadata.tsv", row.names = 1, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
+  sample_data(physeq_its2) <- metadata
+  sample_data(physeq_its2)$Marker <- "ITS2"
+  sample_names(physeq_its2) <- glue::glue(paste("{sample_names(physeq_its2)}","ITS2", sep="_"))
+  saveRDS(physeq_its2,"REAL_CHEESE/ITS2.rds")
+} else {
+  physeq_its2 <- readRDS("REAL_CHEESE/ITS2.rds")
+}
+
+if (!file.exists("REAL_CHEESE/D1D2.rds")) {
+  biomfile <- "REAL_CHEESE/D1D2.biom"
+  physeq_d1d2 <- import_frogs(biomfile, taxMethod = "blast")
+  tax_table(physeq_d1d2) <- gsub(tax_table(physeq_d1d2),     pattern = "[a-z]__", replacement = "")
+  tax_table(physeq_d1d2) <- gsub(tax_table(physeq_d1d2),     pattern = "Eukaryota", replacement = "Fungi")
+  metadata <- read.table("REAL_CHEESE/metadata.tsv", row.names = 1, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
+  sample_data(physeq_d1d2) <- metadata
+  sample_data(physeq_d1d2)$Marker <- "D1D2"
+  sample_names(physeq_d1d2) <- glue::glue(paste("{sample_names(physeq_d1d2)}","D1D2", sep="_"))
+  saveRDS(physeq_d1d2,"REAL_CHEESE/D1D2.rds")
+} else {
+  physeq_d1d2 <- readRDS("REAL_CHEESE/D1D2.rds")
+}
+
+if (!file.exists("REAL_CHEESE/RPB2.rds")) {
+  biomfile <- "REAL_CHEESE/RPB2.biom"
+  physeq_rpb2 <- import_frogs(biomfile, taxMethod = "blast")
+  tax_table(physeq_rpb2) <- gsub(tax_table(physeq_rpb2),     pattern = "[a-z]__", replacement = "")
+  metadata <- read.table("REAL_CHEESE/metadata.tsv", row.names = 1, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
+  sample_data(physeq_rpb2) <- metadata
+  sample_data(physeq_rpb2)$Marker <- "RPB2"
+  sample_names(physeq_rpb2) <- glue::glue(paste("{sample_names(physeq_rpb2)}","RPB2", sep="_"))
+  saveRDS(physeq_rpb2,"REAL_CHEESE/RPB2.rds")
+} else {
+  physeq_rpb2 <- readRDS("REAL_CHEESE/RPB2.rds")
+}
 ```
 
+```{r, eval=T}
+all <- merge_phyloseq(physeq_its1, physeq_its2, physeq_d1d2, physeq_rpb2)
+
+sample_data(all)$Marker <- factor(sample_data(all)$Marker, levels=c("ITS1","ITS2","D1D2","RPB2"))
 
-## Preparation of final BIOM files {.tabset}
+p <- plot_composition(physeq = all, taxaRank1 = "Kingdom", taxaSet1 = "Fungi", taxaRank2 = "Species", numberOfTaxa = 22, x = "Sample")
+p + facet_grid(". ~ Marker", scales = "free_x", space = "free")
 
-### ITS1
+```
 
-### ITS2
+```{r, eval=T}
+p <- plot_richness(physeq = all, x = "Marker", color = "Marker", shape = NULL, title = "Alpha diversity graphics", measures = c("Observed", "Chao1", "ACE", "Shannon", "Simpson", "InvSimpson", "Fisher"))
+p + geom_boxplot() + NULL
+```
 
-### D1D2
+```{r, eval=T}
+beta.dist <- distance(all, method = "bray")
+ord <- ordinate(all, method = "MDS", distance = beta.dist)
+p <- plot_ordination(physeq = all, ordination = ord, type = "samples", axes = c(1, 2), color = "Sample", shape = "Marker", label = NULL, title = "Samples ordination graphic")
+p <- p + stat_ellipse(aes_string(group = "Marker"))
+p + theme_bw()
+```
 
-### RPB2
+```{r, eval=T}
+x <- list(
+  ITS1 = as_tibble(tax_table(physeq_its1))$Genus %>% as.vector(), 
+  ITS2 = as_tibble(tax_table(physeq_its2))$Genus %>% as.vector(),
+  D1D2 = as_tibble(tax_table(physeq_d1d2))$Genus %>% as.vector(),
+  RPB2 = as_tibble(tax_table(physeq_rpb2))$Genus %>% as.vector()
+  )
+ggVennDiagram(x, show_intersect = T)
 
-## Analysis
+```