real samples

real_samples.Rmd

@@ -126,6 +126,8 @@ for(eco in ecos){
       physeq <- import_frogs(biomfile, taxMethod = "blast")
       tax_table(physeq) <- gsub(tax_table(physeq),pattern = "[a-z]__", replacement = "")
       metadata <- read.table(file.path(dir,"metadata.tsv"), row.names = 1, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
+      tax_table(physeq)[,colnames(tax_table(physeq))] <- gsub(tax_table(physeq)[,colnames(tax_table(physeq))],pattern="Debaromycetaceae",replacement="Debaryomycetaceae")
+      tax_table(physeq)[,colnames(tax_table(physeq))] <- gsub(tax_table(physeq)[,colnames(tax_table(physeq))],pattern="Eukaryota",replacement="Fungi")
       sample_data(physeq) <- metadata
       sample_data(physeq)$Marker <- marker
       sample_data(physeq)$Type <- "Raw"
@@ -163,6 +165,9 @@ for(eco in ecos){
       sample_data(physeq_final)$Type <- "Curated"
       #sample_data(physeq_final)$Reads <- sample_sums(physeq_final)
       sample_names(physeq_final) <- glue::glue(paste("{sample_data(physeq)$SampleID}",marker,"CURATED", sep="_"))
+      tax_table(physeq_final)[,colnames(tax_table(physeq_final))] <- gsub(tax_table(physeq_final)[,colnames(tax_table(physeq_final))],pattern="Debaromycetaceae",replacement="Debaryomycetaceae")
+      saveRDS(physeq_final,file.path(dir,paste0(marker,"_final.rds")))
+      tax_table(physeq_final)[,colnames(tax_table(physeq_final))] <- gsub(tax_table(physeq_final)[,colnames(tax_table(physeq_final))],pattern="Eukaryota",replacement="Fungi")
       saveRDS(physeq_final,file.path(dir,paste0(marker,"_final.rds")))
     #}
   }
@@ -172,7 +177,7 @@ for(eco in ecos){
 
 
 ```{r, eval=T, results="asis"}
-ecos <- c("MEAT","CHEESE","WINE")#,"WINE","CHEESE")
+ecos <- c("MEAT","CHEESE","WINE")
 markers <- c("ITS1","ITS2","D1D2","RPB2")
 for(eco in ecos){
   dir <- paste("REAL",eco, sep="_")
@@ -221,6 +226,68 @@ for(eco in ecos){
 }
 ```
 
+# Comparison
+
+
+```{r, eval=T, results="asis"}
+ecos <- c("MEAT","CHEESE","WINE")
+markers <- c("ITS1","ITS2","D1D2","RPB2")
+ranks <- c("Phylum","Class","Order","Family","Genus", "Species")
+for(eco in ecos){
+  cat("\n## ", eco, "{.tabset}\n\n")
+  dir <- paste("REAL",eco, sep="_")
+  its1 <- readRDS(file.path(dir,paste0("ITS1_final",".rds")))
+  its2 <- readRDS(file.path(dir,paste0("ITS2_final",".rds")))
+  d1d2 <- readRDS(file.path(dir,paste0("D1D2_final",".rds")))
+  rpb2 <- readRDS(file.path(dir,paste0("RPB2_final",".rds")))
+  all <- merge_phyloseq(its1, its2, d1d2, rpb2)
+  all_rel  = transform_sample_counts(all, function(x) x / sum(x) )
+  all_rel <- prune_samples(sample_sums(all_rel) > 0, all_rel)
+  
+  for(rank in ranks){
+    cat("\n### ",rank, " \n")
+    all_rel_glom <- tax_glom(all_rel, taxrank=rank)
+    p <- plot_heatmap(all_rel_glom, distance="bray", sample.label="SampleID", taxa.label=rank, low="#1a9850", high="#FF3300", na.value = "white")
+    
+    p <- p + facet_grid(~Marker, scales = "free_x", space = "free_x")
+    print(p)
+    cat("\n")
+    x <- list(
+      ITS1 = tax_table(its1) %>% as("matrix") %>% as_tibble() %>% select(rank) %>% as.vector() %>% unlist() %>% unname(),
+      ITS2 = tax_table(its2) %>% as("matrix") %>% as_tibble() %>% select(rank) %>% as.vector() %>% unlist() %>% unname(),
+      D1D2 = tax_table(d1d2) %>% as("matrix") %>% as_tibble() %>% select(rank) %>% as.vector() %>% unlist() %>% unname(),
+      RPB2 = tax_table(rpb2) %>% as("matrix") %>% as_tibble() %>% select(rank) %>% as.vector() %>% unlist() %>% unname()
+      )
+    print(ggVennDiagram(x,label_alpha = 0) + 
+      scale_fill_gradient(low = "#F4FAFE", high = "#4981BF") )
+    cat("\n")
+    x <- list(
+      ITS1 = tax_table(its1) %>% as("matrix") %>% as_tibble() %>% select(rank) %>% as.vector() %>% unlist() %>% unname() %>% unique(),
+      ITS2 = tax_table(its2) %>% as("matrix") %>% as_tibble() %>% select(rank) %>% as.vector() %>% unlist() %>% unname() %>% unique(),
+      D1D2 = tax_table(d1d2) %>% as("matrix") %>% as_tibble() %>% select(rank) %>% as.vector() %>% unlist() %>% unname() %>% unique(),
+      RPB2 = tax_table(rpb2) %>% as("matrix") %>% as_tibble() %>% select(rank) %>% as.vector() %>% unlist() %>% unname() %>% unique()
+      )
+    print(x)
+    cat("\n")
+  }
+}
+
+
+```
+
+
+```{r}
+x <- list(
+  ITS1 = as_tibble(tax_table(physeq_its1_final))$Family %>% as.vector(), 
+  ITS2 = as_tibble(tax_table(physeq_its2_final))$Family %>% as.vector(),
+  D1D2 = as_tibble(tax_table(physeq_d1d2_final))$Family %>% as.vector(),
+  RPB2 = as_tibble(tax_table(physeq_rpb2_final))$Family %>% as.vector()
+  )
+ggVennDiagram(x,label_alpha = 0) + 
+  scale_fill_gradient(low = "#F4FAFE", high = "#4981BF") 
+```
+
+