#Inputting phyloseq Object

## 
## microbiome R package (microbiome.github.com)
##     
## 
## 
##  Copyright (C) 2011-2022 Leo Lahti, 
##     Sudarshan Shetty et al. <microbiome.github.io>
## 
## Attaching package: 'microbiome'
## The following object is masked from 'package:ggplot2':
## 
##     alpha
## The following object is masked from 'package:base':
## 
##     transform
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
rm(list = ls())
ps1 <- readRDS("../00_FilesFromWynton/BoomBustPhylo.rds")
pseq0 <- ps1 %>%
  subset_samples(!Category %in% c( "Control_Slice", "Control_Slice", "Control_Lawn", "Grazed_Slice","Worm") & !SampleID %in% c("ACME29", "PCME64",
"PCME71", 
"PCME79", 
"BCME82", 
"BCME83", 
"BCME84", 
"PCME85",  
"ACME70", 
"ACME89", 
"ACME90", 
"ACME91", 
"ACME92", 
"ACME93", 
"ACME94", 
"ACME95") & !CME %in% c("BCME", "PCME"))
pseq <- subset_samples(pseq0,  Day %in% c( "Day22", "control"))
taxa_names(pseq) <- paste0("SV", seq(ntaxa(pseq)))
metadata <- sample_data(pseq)
#View(metadata)
otus = as(otu_table(pseq), "matrix")
#head(otus)

extract only those samples in common between the two tables

common.sample.ids <- intersect(rownames(metadata), rownames(otus))
otus <- otus[common.sample.ids,]
metadata <- metadata[common.sample.ids,]
# double-check that the mapping file and otu table
# had overlapping samples
if(length(common.sample.ids) <= 1) {
    message <- paste(sprintf('Error: there are %d sample ids in common '),
                    'between the metadata file and data table')
    stop(message)
}

extract the source environments and source/sink indices

train.ix <- which(metadata$SourceSink=='Source')
test.ix <- which(metadata$SourceSink=='Sink')
envs <- metadata$Category
if(is.element('Description',colnames(metadata))) desc <- metadata$Category

load SourceTracker package

source('SourceTracker.r')

#NOTE skipping tuning

alpha1 <- alpha2 <- 0.001

train SourceTracker object on training data

st <- sourcetracker(otus[train.ix,], envs[train.ix])
## Rarefying training data at 1000

Estimate source proportions in test data

results <- predict(st,otus[test.ix,], alpha1=alpha1, alpha2=alpha2)
##                                           Control        Unknown
## ..........   1 of    3, depth= 1000:  0.01 (0.00)    0.99 (0.00)
## ..........   2 of    3, depth= 1000:  0.01 (0.00)    0.99 (0.00)
## ..........   3 of    3, depth= 1000:  0.01 (0.00)    0.99 (0.00)

Estimate leave-one-out source proportions in training data

results.train <- predict(st, alpha1=alpha1, alpha2=alpha2)
## ndraws=10, V=2, T=39552, N=3
## Rarefying training data at 1000
## ..........   1 of    3, depth= 1000:  0.93 (0.00)    0.07 (0.00)
## Rarefying training data at 1000
## ..........   2 of    3, depth= 1000:  0.67 (0.00)    0.33 (0.00)
## Rarefying training data at 1000
## ..........   3 of    3, depth= 1000:  0.83 (0.00)    0.17 (0.00)

plot results

desc <- metadata$Category

labels <- sprintf('%s %s', envs,desc)
plot(results, labels[test.ix], type='pie')
plot(results, labels[test.ix], type='bar')

plot(results, labels[test.ix], type='dist')

plot(results.train, labels[train.ix], type='pie')

plot(results.train, labels[train.ix], type='bar')

plot(results.train, labels[train.ix], type='dist')

# plot results with legend
plot(results, labels[test.ix], type='pie', include.legend=TRUE, env.colors=c("#33ff33",  "#00cccc","#006666","#ff66ff", "darkblue","#800000","#4B0082","#fb9a99","#FF00FF","dodgerblue3","darkgoldenrod1"))

downstream <- data.frame(results$proportions)
downstream$id <- row.names(downstream)
meltdown <- reshape2::melt(downstream, id.vars = "id")
std.env.colors <- c("Control" = "red", "Unknown" = "darkgreen")

# Re-create the plot
p <- ggplot(meltdown, aes(x = id, y = value)) + 
  geom_bar(aes(fill = variable), stat = "identity") + 
  scale_fill_manual(values = std.env.colors) + 
  scale_x_discrete(drop = TRUE) + 
  guides(fill = guide_legend(keywidth = 1, keyheight = 1)) +
  theme_minimal(base_size = 16) + 
  theme(
    panel.grid.major = element_blank(),
    axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, face = "bold"),
    axis.text.y = element_text(hjust = 1, face = "bold"),
    axis.text = element_text(face = "bold"),
    strip.text.x = element_text(size = 12),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.8),
    axis.line = element_line(color = "black"), 
    axis.ticks = element_line(color = "black") ) 

p

ggsave("BustPhaseSamples.pdf", height = 4, width = 6)