A set of R functions that help faciliate a lot of tedious processing
install.packages('devtools')
devtools::install_github('twbattaglia/btools')
library(btools)
# Run in Terminal
biom convert \
-i otu_table.biom \
-o otu_table_json.biom \
--to-json \
--otu-table 'OTU table'
# 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")
For each PCR run, remove the blanks that correspond to each plate.
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")
Thanks to NanoStringNorm
genes <- import_rcc("cel_files/")
# Calculate BF ratio
phyloseq <- bf_ratio(phyloseq)
# View log2 BF ratio's
phyloseq$log2_bf_ratio
Thanks to DESeq2
plotPCA3D(deseq2, intgroup = "Treatment")