<> # Action performed on: Nov 5 2020 # Action performed by: Sara Kowalski & Elizabeth Vaz # Imported the data using the qiime import command qiime tools import \ --type 'SampleData[SequencesWithQuality]' \ --input-path /mnt/datasets/project_2/soil/soil_metadata.txt \ --output-path soil_demux.qza \ --input-format SingleEndFastqManifestPhred33V2 # Summarized data qiime demux summarize \ --i-data soil_demux.qza \ --o-visualization soil_demux.qzv #Created screen for quality control screen -s QC #ASV's truncated based on 150 cutoff qiime dada2 denoise-single --i-demultiplexed-seqs soil-demux.qza --p-trim-left 0 --p-trunc-len 150 --o-representative-sequences trunc-150-rep-seqs.qza --o-table trunc-150-table.qza --o-denoising-stats trunc-150-stats.qza #ASV's truncated based on 200 cutoff qiime dada2 denoise-single --i-demultiplexed-seqs soil-demux.qza --p-trim-left 0 --p-trunc-len 200 --o-representative-sequences trunc-200-rep-seqs.qza --o-table trunc-200-table.qza --o-denoising-stats trunc-200-stats.qza #Visualize table with 200 truncation qiime feature-table summarize --i-table trunc_200_table.qza --o-visualization trunc_200_table.qzv --m-sample-metadata-file /mnt/datasets/project_2/soil/soil_metadata.txt #Output files created in working directory (/data) from session #soil-demux.qza #soil_demux.qzv #trunc-150-stats.qza #trunc-150-table.qza #trunc-150-rep-seqs.qza #trunc-200-rep-seqs.qza #trunc-200-table.qza #trunc-200-stats.qza #trunc-200-rep-seqs.qza #trunc_200_table.qzv # Action performed on: Nov 6 2020 # Action performed by: Mark Pitblado # Removed observations that contained either an N/A or 0 for the Total Nitrogen, pH, Total Carbon, or Moisture Content variables using the 200 truncated table qiime feature-table filter-samples --i-table trunc_200_table.qza --m-metadata-file /mnt/datasets/project_2/soil/soil_metadata.txt --p-where "[Total Nitrogen]!='0' AND [Total Nitrogen] != 'NA' AND [pH] != '0' AND [pH] != 'NA' AND [Total Carbon] != 'NA' AND [Moisture Content] != 'NA' AND [Moisture Content] != '0'" --o-filtered-table All_filtered_200_table.qza # Generated visualization of filtered table and saved to the working directory qiime feature-table summarize --i-table All_filtered_200_table.qza --o-visualization all_filtered_200_table.qzv --m-sample-metadata-file /mnt/datasets/project_2/soil/soil_metadata.txt #Output files created in working directory (/data) from session # All_filtered_200_table.qza # all_filtered_200_table.qzv # Action performed on: Nov 7 2020 # Action performed by: Mark Pitblado # Screen started screen -s "Root Tree" # Generation of phylogentic tree files qiime phylogeny align-to-tree-mafft-fasttree --i-sequences trunc_200_rep-seqs.qza --o-alignment aligned-rep-seqs.qza --o-masked-alignment masked-aligned-rep-seqs.qza --o-tree unrooted-tree.qza --o-rooted-tree rooted-tree.qza #Output files created in working directory (/data) from session # aligned-rep-seqs.qza # masked-aligned-rep-seqs.qza # unrooted-tree.qza # rooted-tree.qza <> --- title: "Association of Abiotic Soil Factors with Bacterial Diversity and Abundance in North American Regions of Long-Term Reforestation" author: "Mark Pitblado" date: Dec 15, 2020 output: pdf_document: highlight: zenburn toc: true --- ```{r Setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) #install.packages(tidyverse) #install.packages(vegan) #install.packages(phyloseq) #install.packages(DESeq2) #install.packages(ggplot2) #install.packages(dplyr) #install.packages(ggthemes) #install.packages(ggpubr) #install.packages(viridis) #install.packages(ape) #install.packages(plyr) #install.packages(maps) #install.packages(mapdata) #install.packages(scales) #install.packages(ggsignif) #install.packages(Hmisc) #install.packages(aod) require(tidyverse) require(vegan) require(phyloseq) require(DESeq2) require(ggplot2) require(dplyr) require(ggthemes) require(ggpubr) require(viridis) require(ape) require(plyr) require(maps) require(mapdata) require(maptools) require(scales) require(ggsignif) require(Hmisc) require(aod) library(tidyverse) library(vegan) library(phyloseq) library(DESeq2) library(ggplot2) library(dplyr) library(ggthemes) library(ggpubr) library(viridis) library(ape) library(plyr) library(maps) library(mapdata) library(maptools) library(scales) library(ggsignif) library(Hmisc) library(aod) ``` ```{r Import the Data} #Importing Files from Qiime2 workflow setwd("~/Documents/University/Year_4/MICB_447/Project_2/MICB_447") Sys.setlocale("LC_ALL","C") biom_file <- import_biom("table-with-taxonomy.biom") metadata <- import_qiime_sample_data("soil_metadata.tsv") metadata$Sampling.Depth<-as.character(metadata$Sampling.Depth) tree <- read_tree_greengenes("tree.nwk") tree <- multi2di(tree) # assuming you read in your tree.nwk file as the variable tree ``` ```{r Cleaning Data} #Ensure that 0 and NA values in abiotic factors are removed from the metadata file metadata<-subset(metadata, Total.Carbon != 0 & Total.Carbon != "NA" & Total.Nitrogen != 0 & Total.Nitrogen != "NA" & pH != 0 & pH != "NA" & Moisture.Content != 0 & Moisture.Content != "NA") ``` ```{r Setting up Phyloseq} #Setup of Phyloseq #Combine all objects into phyloseq object physeq_untrimmed <- merge_phyloseq(biom_file, metadata, tree) #Filter out low abundance ASV's physeq = prune_samples(sample_sums(physeq_untrimmed)>=10,physeq_untrimmed) #Overview of phyloseq object physeq # Set set of random numbers (711 is default) set.seed(711) ``` ```{r Renaming Taxonomic Ranks} #Renaming #Rename column names for taxonomic ranks colnames(tax_table(physeq)) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species") #Check that renaming was successful rank_names(physeq) ``` ```{r Generating Region Level PCoA Plot} #PCoA Plots (ALL) #Rarefy Data, convert to relative abundance, and produce weighted unifrac PCoA coloured by region physeq_rar <- rarefy_even_depth(physeq, sample.size = 5000) physeq_rar_RA <- transform_sample_counts(physeq_rar, function(x) x/sum(x)) ord <- ordinate(physeq_rar_RA, method = "PCoA", distance = "wunifrac") PCoA_Region <- plot_ordination(physeq_rar_RA, ord, type = "sample", color = ("Region"), title = "PCoA (Weighted Unifrac) - Region") PCoA_Region + theme_classic(base_size = 20) + scale_color_manual(values = c("#f35e5a", "#6ba204", "#17b3b7", "#b95eff")) + stat_ellipse(type = "norm", linetype = 2 ) ``` ```{r Subsetting Data by Region} #Split Regions #Seperate by region, as distinct clustering observed in ALL PCoA plot region_BC <- subset_samples(physeq, Region == "British Columbia") region_ON <- subset_samples(physeq, Region == "Ontario") region_CA <- subset_samples(physeq, Region == "California") region_TX <- subset_samples(physeq, Region == "Texas") ``` ``` {r PCoA Plots: British Columbia} #PCoA Plots (British Columbia) #Rarefy Data, convert to relative abundance, produce PCoA plots coloured by abiotic factors physeq_rar_BC <- rarefy_even_depth(region_BC, sample.size = 5000) physeq_rar_RA_BC <- transform_sample_counts(physeq_rar_BC, function(x) x/sum(x)) ord_BC <- ordinate(physeq_rar_RA_BC, method = "PCoA", distance = "wunifrac") #Total Carbon BC_PCoA_Total.Carbon <- plot_ordination(physeq_rar_RA_BC, ord_BC, type = "sample", color = ("Total.Carbon"), title = "PCoA (Weighted Unifrac) British Columbia - Total Carbon") BC_PCoA_Total.Carbon + theme_classic() + scale_color_gradient(low="#F805CE", high="#0511E8") #Total Nitrogen BC_PCoA_Total.Nitrogen <- plot_ordination(physeq_rar_RA_BC, ord_BC, type = "sample", color = ("Total.Nitrogen"), title = "PCoA (Weighted Unifrac) British Columbia - Total Nitrogen") BC_PCoA_Total.Nitrogen + theme_classic() + scale_color_gradient(low="#F805CE", high="#0511E8") #pH BC_PCoA_pH <- plot_ordination(physeq_rar_RA_BC, ord_BC, type = "sample", color = ("pH"), title = "PCoA (Weighted Unifrac) British Columbia - pH") BC_PCoA_pH + theme_classic() + scale_color_gradient(low="#F805CE", high="#0511E8") #Moisture Content BC_PCoA_Moisture.Content <- plot_ordination(physeq_rar_RA_BC, ord_BC, type = "sample", color = ("Moisture.Content"), title = "PCoA (Weighted Unifrac) British Columbia - Moisture Content") BC_PCoA_Moisture.Content + theme_classic() + scale_color_gradient(low="#F805CE", high="#0511E8") #Sampling Depth BC_PCoA_Sampling.Depth <- plot_ordination(physeq_rar_RA_BC, ord_BC, type = "sample", color = ("Sampling.Depth"), title = "PCoA (Weighted Unifrac) British Columbia - Sampling Depth") BC_PCoA_Sampling.Depth + theme_classic() + scale_color_manual(values =c("#F805CE", "#0511E8")) ``` ``` {r PCoA Plots Ontario} #PCoA Plots (Ontario) #Rarefy Data, convert to relative abundance, and create ordination physeq_rar_ON <- rarefy_even_depth(region_ON, sample.size = 5000) physeq_rar_RA_ON <- transform_sample_counts(physeq_rar_ON, function(x) x/sum(x)) ord_ON <- ordinate(physeq_rar_RA_ON, method = "PCoA", distance = "wunifrac") #Total Carbon ON_PCoA_Total.Carbon <- plot_ordination(physeq_rar_RA_ON, ord_ON, type = "sample", color = ("Total.Carbon"), title = "PCoA (Weighted Unifrac) Ontario - Total Carbon") ON_PCoA_Total.Carbon + theme_classic() + scale_color_gradient(low="#F805CE", high="#0511E8") #Total Nitrogen ON_PCoA_Total.Nitrogen <- plot_ordination(physeq_rar_RA_ON, ord_ON, type = "sample", color = ("Total.Nitrogen"), title = "PCoA (Weighted Unifrac) Ontario - Total Nitrogen") ON_PCoA_Total.Nitrogen + theme_classic() + scale_color_gradient(low="#F805CE", high="#0511E8") #pH ON_PCoA_pH <- plot_ordination(physeq_rar_RA_ON, ord_ON, type = "sample", color = ("pH"), title = "PCoA (Weighted Unifrac) Ontario - pH") ON_PCoA_pH + theme_classic() + scale_color_gradient(low="#F805CE", high="#0511E8") #Moisture Content ON_PCoA_Moisture.Content <- plot_ordination(physeq_rar_RA_ON, ord_ON, type = "sample", color = ("Moisture.Content"), title = "PCoA (Weighted Unifrac) Ontario - Moisture Content") ON_PCoA_Moisture.Content + theme_classic() + scale_color_gradient(low="#F805CE", high="#0511E8") #Sampling Depth ON_PCoA_Sampling.Depth <- plot_ordination(physeq_rar_RA_ON, ord_ON, type = "sample", color = ("Sampling.Depth"), title = "PCoA (Weighted Unifrac) Ontario - Sampling Depth") ON_PCoA_Sampling.Depth + theme_classic() + scale_color_manual(values =c("#F805CE", "#0511E8")) ``` ``` {r PCoA Plots Texas} #Texas #Rarefy Data, convert to relative abundance, produce ordination and PCoA plots by abiotic factors physeq_rar_TX <- rarefy_even_depth(region_TX, sample.size = 5000) physeq_rar_RA_TX <- transform_sample_counts(physeq_rar_TX, function(x) x/sum(x)) ord_TX <- ordinate(physeq_rar_RA_TX, method = "PCoA", distance = "wunifrac") #Total Carbon TX_PCoA_Total.Carbon <- plot_ordination(physeq_rar_RA_TX, ord_TX, type = "sample", color = ("Total.Carbon"), title = "PCoA (Weighted Unifrac) Texas - Total Carbon") TX_PCoA_Total.Carbon + theme_classic() + scale_color_gradient(low="#F805CE", high="#0511E8") #Total Nitrogen TX_PCoA_Total.Nitrogen <- plot_ordination(physeq_rar_RA_TX, ord_TX, type = "sample", color = ("Total.Nitrogen"), title = "PCoA (Weighted Unifrac) Texas - Total Nitrogen") TX_PCoA_Total.Nitrogen + theme_classic() + scale_color_gradient(low="#F805CE", high="#0511E8") #pH TX_PCoA_pH <- plot_ordination(physeq_rar_RA_TX, ord_TX, type = "sample", color = ("pH"), title = "PCoA (Weighted Unifrac) Texas - pH") TX_PCoA_pH + theme_classic() + scale_color_gradient(low="#F805CE", high="#0511E8") #Moisture Content TX_PCoA_Moisture.Content <- plot_ordination(physeq_rar_RA_TX, ord_TX, type = "sample", color = ("Moisture.Content"), title = "PCoA (Weighted Unifrac) Texas - Moisture Content") TX_PCoA_Moisture.Content + theme_classic() + scale_color_gradient(low="#F805CE", high="#0511E8") #Sampling Depth TX_PCoA_Sampling.Depth <- plot_ordination(physeq_rar_RA_TX, ord_TX, type = "sample", color = ("Sampling.Depth"), title = "PCoA (Weighted Unifrac) Texas - Sampling Depth") TX_PCoA_Sampling.Depth + theme_classic() + scale_color_manual(values =c("#F805CE", "#0511E8")) ``` ```{r PCoA Plots California} #California #Rarefy Data, calculate relative abundance, create ordination and produce PCoA with abiotic factors physeq_rar_CA <- rarefy_even_depth(region_CA, sample.size = 5000) physeq_rar_RA_CA <- transform_sample_counts(physeq_rar_CA, function(x) x/sum(x)) ord_CA <- ordinate(physeq_rar_RA_CA, method = "PCoA", distance = "wunifrac") #Total Carbon CA_PCoA_Total.Carbon <- plot_ordination(physeq_rar_RA_CA, ord_CA, type = "sample", color = ("Total.Carbon"), title = "PCoA (Weighted Unifrac) California - Total Carbon") CA_PCoA_Total.Carbon + theme_classic() + scale_color_gradient(low="#F805CE", high="#0511E8") #Total Nitrogen CA_PCoA_Total.Nitrogen <- plot_ordination(physeq_rar_RA_CA, ord_CA, type = "sample", color = ("Total.Nitrogen"), title = "PCoA (Weighted Unifrac) California - Total Nitrogen") CA_PCoA_Total.Nitrogen + theme_classic() + scale_color_gradient(low="#F805CE", high="#0511E8") #pH CA_PCoA_pH <- plot_ordination(physeq_rar_RA_CA, ord_CA, type = "sample", color = ("pH"), title = "PCoA (Weighted Unifrac) California - pH") CA_PCoA_pH + theme_classic() + scale_color_gradient(low="#FFB336", high="#0511E8") #Moisture Content CA_PCoA_Moisture.Content <- plot_ordination(physeq_rar_RA_CA, ord_CA, type = "sample", color = ("Moisture.Content"), title = "PCoA (Weighted Unifrac) California - Moisture Content") CA_PCoA_Moisture.Content + theme_classic() + scale_color_gradient(low="#F805CE", high="#0511E8") #Sampling Depth CA_PCoA_Sampling.Depth <- plot_ordination(physeq_rar_RA_CA, ord_CA, type = "sample", color = ("Sampling.Depth"), title = "PCoA (Weighted Unifrac) California - Sampling Depth") CA_PCoA_Sampling.Depth + theme_classic() + scale_color_manual(values =c("#F805CE", "#0511E8")) ``` ``` {r Abundance Bar Plots} #Abundance Bar Plots (ALL) # physeq.phylum = tax_glom(physeq_rar_RA, taxrank ="Phylum", NArm=FALSE) physeq.phylum.BC = tax_glom(physeq_rar_RA_BC, taxrank ="Phylum", NArm=TRUE) plot_bar(physeq.phylum, fill="Phylum") + facet_wrap(~Region, scales= "free_x", nrow=4) + theme_classic() ``` ```{r Boxplots} #BOX PLOTs #ALL //////////////////////////////////////////////////////////////////// plot_richness(physeq_rar, x="Region", title = "Alpha Diversity by Region", measures=c("Chao1", "Observed","Simpson")) + theme_classic() + geom_boxplot(outlier.shape = NA, notch = TRUE) #Nitrogen by Sampling Depth Box Plot Total_Nitrogen_Boxplot <- ggplot(metadata, aes(x=Sampling.Depth, y=Total.Nitrogen, fill=Region)) + geom_boxplot() + facet_wrap(~Region) + theme(legend.position = "none") + stat_compare_means(aes(group = Sampling.Depth)) Total_Nitrogen_Boxplot #Carbon by Sampling Depth Box Plot Total_Carbon_Boxplot <- ggplot(metadata, aes(x=Sampling.Depth, y=Total.Carbon, fill=Region)) + geom_boxplot() + facet_wrap(~Region) + theme(legend.position = "none") + stat_compare_means(aes(group = Sampling.Depth)) Total_Carbon_Boxplot #pH by Sampling Depth Box Plot Total_pH_Boxplot <- ggplot(metadata, aes(x=Sampling.Depth, y=pH, fill=Region)) + geom_boxplot() + facet_wrap(~Region) + theme(legend.position = "none") + stat_compare_means(aes(group = Sampling.Depth)) Total_pH_Boxplot #pH by Sampling Depth Box Plot Total_Moisture_Boxplot <- ggplot(metadata, aes(x=Sampling.Depth, y=Moisture.Content, fill=Region)) + geom_boxplot() + facet_wrap(~Region) + theme(legend.position = "none") + stat_compare_means(aes(group = Sampling.Depth)) Total_Moisture_Boxplot #Nitrogen by OM group Box Plot Total_Nitrogen_Boxplot_OM <- ggplot(metadata, aes(x=LTSP.Treatment, y=Total.Nitrogen, fill=Region)) + geom_boxplot() + facet_wrap(~Region) + theme(legend.position = "none") + stat_compare_means(aes(group = LTSP.Treatment)) Total_Nitrogen_Boxplot_OM #Carbon by OM group Box Plot Total_Carbon_Boxplot_OM <- ggplot(metadata, aes(x=LTSP.Treatment, y=Total.Carbon, fill=Region)) + geom_boxplot() + facet_wrap(~Region) + theme(legend.position = "none") + stat_compare_means(aes(group = LTSP.Treatment)) Total_Carbon_Boxplot_OM #oH by OM group Box Plot Total_pH_Boxplot_OM <- ggplot(metadata, aes(x=LTSP.Treatment, y=pH, fill=Region)) + geom_boxplot() + facet_wrap(~Region) + theme(legend.position = "none") + stat_compare_means(aes(group = LTSP.Treatment)) Total_pH_Boxplot_OM #Moisture Content by OM group Box Plot Total_Moisture_Boxplot_OM <- ggplot(metadata, aes(x=LTSP.Treatment, y=Moisture.Content, fill=Region)) + geom_boxplot() + facet_wrap(~Region) + theme(legend.position = "none") + stat_compare_means(aes(group = LTSP.Treatment)) Total_Moisture_Boxplot_OM ``` ```{r Sample Site Maps} #MAP of Sampling Sites map("worldHires","Canada", xlim=c(-141,-53),ylim=c(40,85), col="gray90", fill=TRUE) points(metadata$Longitude, metadata$Latitude, pch=19, col="red", cex=1) map("worldHires","USA", xlim=c(-125.75583,-68.01197),ylim=c(19.50139,64.85694), col="gray90", fill=TRUE) points(metadata$Longitude, metadata$Latitude, pch=19, col="red", cex=0.5) ``` ```{r Creation of Ecozones} #Creation of Ecozones IDFBC <- subset_samples(physeq_untrimmed, Ecozone == "IDFBC") SBSBC <- subset_samples(physeq_untrimmed, Ecozone == "SBSBC") BSON <- subset_samples(physeq_untrimmed, Ecozone == "BSON") JPON <- subset_samples(physeq_untrimmed, Ecozone == "JPON") LPTX <- subset_samples(physeq_untrimmed, Ecozone == "LPTX") PPCA <- subset_samples(physeq_untrimmed, Ecozone == "PPCA") ``` ```{r Calculating Alpha Diversity for Ecozones} #Shannon Alpha Diversity estimates by ecozone alpha_IDFBC <- estimate_richness(IDFBC, measures = "Shannon") alpha_SBSBC <- estimate_richness(SBSBC, measures = "Shannon") alpha_BSON <- estimate_richness(BSON, measures = "Shannon") alpha_JPON <- estimate_richness(JPON, measures = "Shannon") alpha_LPTX <- estimate_richness(LPTX, measures = "Shannon") alpha_PPCA <- estimate_richness(PPCA, measures = "Shannon") ``` ```{r Linear Regression Models Nitrogen} #Linear Regressions x = Nitrogen y = Shannon Diversity by Ecozone #IDFBC df = cbind(sample_data(IDFBC)[, "Total.Nitrogen", drop=FALSE], alpha_IDFBC[1]) fit = lm(unlist(alpha_IDFBC) ~ Total.Nitrogen, data = df) summary(fit) #SBSBC df = cbind(sample_data(SBSBC)[, "Total.Nitrogen", drop=FALSE], alpha_SBSBC[1]) fit = lm(unlist(alpha_SBSBC) ~ Total.Nitrogen, data = df) summary(fit) #BSON df = cbind(sample_data(BSON)[, "Total.Nitrogen", drop=FALSE], alpha_BSON[1]) fit = lm(unlist(alpha_BSON) ~ Total.Nitrogen, data = df) summary(fit) #JPON df = cbind(sample_data(JPON)[, "Total.Nitrogen", drop=FALSE], alpha_JPON[1]) fit = lm(unlist(alpha_JPON) ~ Total.Nitrogen, data = df) summary(fit) #LPTX df = cbind(sample_data(LPTX)[, "Total.Nitrogen", drop=FALSE], alpha_LPTX[1]) fit = lm(unlist(alpha_LPTX) ~ Total.Nitrogen, data = df) summary(fit) #PPCA df = cbind(sample_data(PPCA)[, "Total.Nitrogen", drop=FALSE], alpha_PPCA[1]) fit = lm(unlist(alpha_PPCA) ~ Total.Nitrogen, data = df) summary(fit) ``` ```{r Linear Regression Models Carbon} #Linear Regressions x = Carbon y = Shannon Diversity by Ecozone #IDFBC df = cbind(sample_data(IDFBC)[, "Total.Carbon", drop=FALSE], alpha_IDFBC[1]) fit = lm(unlist(alpha_IDFBC) ~ Total.Carbon, data = df) summary(fit) #SBSBC df = cbind(sample_data(SBSBC)[, "Total.Carbon", drop=FALSE], alpha_SBSBC[1]) fit = lm(unlist(alpha_SBSBC) ~ Total.Carbon, data = df) summary(fit) #BSON df = cbind(sample_data(BSON)[, "Total.Carbon", drop=FALSE], alpha_BSON[1]) fit = lm(unlist(alpha_BSON) ~ Total.Carbon, data = df) summary(fit) #JPON df = cbind(sample_data(JPON)[, "Total.Carbon", drop=FALSE], alpha_JPON[1]) fit = lm(unlist(alpha_JPON) ~ Total.Carbon, data = df) summary(fit) #LPTX df = cbind(sample_data(LPTX)[, "Total.Carbon", drop=FALSE], alpha_LPTX[1]) fit = lm(unlist(alpha_LPTX) ~ Total.Carbon, data = df) summary(fit) #PPCA df = cbind(sample_data(PPCA)[, "Total.Carbon", drop=FALSE], alpha_PPCA[1]) fit = lm(unlist(alpha_PPCA) ~ Total.Carbon, data = df) summary(fit) #pH #IDFBC df = cbind(sample_data(IDFBC)[, "pH", drop=FALSE], alpha_IDFBC[1]) fit = lm(unlist(alpha_IDFBC) ~ pH, data = df) summary(fit) #SBSBC df = cbind(sample_data(SBSBC)[, "pH", drop=FALSE], alpha_SBSBC[1]) fit = lm(unlist(alpha_SBSBC) ~ pH, data = df) summary(fit) #BSON df = cbind(sample_data(BSON)[, "pH", drop=FALSE], alpha_BSON[1]) fit = lm(unlist(alpha_BSON) ~ pH, data = df) summary(fit) #JPON df = cbind(sample_data(JPON)[, "pH", drop=FALSE], alpha_JPON[1]) fit = lm(unlist(alpha_JPON) ~ pH, data = df) summary(fit) #LPTX df = cbind(sample_data(LPTX)[, "pH", drop=FALSE], alpha_LPTX[1]) fit = lm(unlist(alpha_LPTX) ~ pH, data = df) summary(fit) #PPCA df = cbind(sample_data(PPCA)[, "pH", drop=FALSE], alpha_PPCA[1]) fit = lm(unlist(alpha_PPCA) ~ pH, data = df) summary(fit) ``` ```{r Linear Regression Models Moisture Content} #Linear Regressions x = Moisture Content y = Shannon Diversity by Ecozone #IDFBC df = cbind(sample_data(IDFBC)[, "Moisture.Content", drop=FALSE], alpha_IDFBC[1]) fit = lm(unlist(alpha_IDFBC) ~ Moisture.Content, data = df) summary(fit) #SBSBC df = cbind(sample_data(SBSBC)[, "Moisture.Content", drop=FALSE], alpha_SBSBC[1]) fit = lm(unlist(alpha_SBSBC) ~ Moisture.Content, data = df) summary(fit) #BSON df = cbind(sample_data(BSON)[, "Moisture.Content", drop=FALSE], alpha_BSON[1]) fit = lm(unlist(alpha_BSON) ~ Moisture.Content, data = df) summary(fit) #JPON df = cbind(sample_data(JPON)[, "Moisture.Content", drop=FALSE], alpha_JPON[1]) fit = lm(unlist(alpha_JPON) ~ Moisture.Content, data = df) summary(fit) #LPTX df = cbind(sample_data(LPTX)[, "Moisture.Content", drop=FALSE], alpha_LPTX[1]) fit = lm(unlist(alpha_LPTX) ~ Moisture.Content, data = df) summary(fit) #PPCA df = cbind(sample_data(PPCA)[, "Moisture.Content", drop=FALSE], alpha_PPCA[1]) fit = lm(unlist(alpha_PPCA) ~ Moisture.Content, data = df) summary(fit) ``` ``` {r Citations} citation("tidyverse") citation("vegan") citation("phyloseq") citation("DESeq2") citation("ggplot2") #citation(dplyr) #Hadley Wickham, Romain François, Lionel Henry and Kirill Müller (2020). dplyr: A Grammar of Data Manipulation. R package version 1.0.2. https://CRAN.R-project.org/package=dplyr citation("ggthemes") citation("ggpubr") citation("viridis") citation("ape") citation("plyr") citation("maps") citation("mapdata") #citation("maptools") citation("scales") citation("ggsignif") citation("Hmisc") citation("aod") ```