install.packages('devtools') devtools::install_github('twbattaglia/btools') library(btools)

# Install necessary packages source("https://bioconductor.org/biocLite.R") biocLite("phyloseq") biocLite("ggplot2") biocLite("vegan") # Load necessary packages library(phyloseq) library(ggplot2) library(btools) # Import OTU table + tree + map phylo <- create_phylo(biom_fp = "otu_table_json.biom", mappingfile_fp = "mapping_file.txt", tree_fp = "rep_set.tre")

phylo_noblanks <- remove_blanks(phylo = phylo, runID = "PCR_Plate", blankID = "Group", blankName = "blank", removeBlank = FALSE)

compare_alpha_diversity(phylo, x = "Time", group = "Treatment", diversity = "Observed", test_type = "nonparametric", write = T, filename = "adiv_results") 

compare_beta_diversity(phylo, x = "Time", group = "Treatment", bdiv = "unweighted", test = "adonis", write = T, fdr = TRUE, filename = "bdiv_results")

estimate_pd(phylo)

contributions <- analyze_contributions(contributions_fp = "metagenomic_contributions.tab", mappingfile_fp = "mapping_file.txt") # Plot contributions contributions %>% group_by(Gene, Treatment) %>% mutate(Contribution_perc = ContributionPercentOfAllSamples * 100) %>% filter(Contribution_perc >= 0) %>% select(Gene, family, Contribution_perc) %>% mutate(Contribution = Contribution_perc/sum(Contribution_perc) * 100) %>% ggplot(aes(x = Treatment, y = Contribution, fill = family)) + geom_bar(stat = "identity") + theme(axis.text.x = element_text(size = 6)) + scale_y_continuous(expand = c(0, 0), limits = c(0, 100)) + theme_light(base_size = 18) + scale_fill_brewer(palette = "Set1")

jaccard <- diversity_comparison(phylo, distance = "jaccard") jsd <- diversity_comparison(phylo, distance = "jsd") unweighted <- diversity_comparison(phylo, distance = "unifrac") weighted <- diversity_comparison(phylo, distance = "wunifrac")

genes <- import_rcc("cel_files/")

# Calculate BF ratio phyloseq <- bf_ratio(phyloseq) # View log2 BF ratio's phyloseq$log2_bf_ratio

plotPCA3D(deseq2, intgroup = "Treatment")