R Notebook: Provides reproducible analysis for 16S rRNA data in the following manuscript:
Citation: Romanowicz KJ and Kling GW. (In Press) Summer thaw duration is a strong predictor of the soil microbiome and its response to permafrost thaw in arctic tundra. Environmental Microbiology. https://doi.org/10.1111/1462-2920.16218
GitHub Repository: https://github.com/kromanowicz/2022-Annual-Thaw-Microbes
NCBI BioProject: https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA794857
Accepted for Publication: 22 September 2022 Environmental Microbiology
This R Notebook provides complete reproducibility of the data analysis presented in “Summer thaw duration is a strong predictor of the soil microbiome and its response to permafrost thaw in arctic tundra” by Romanowicz and Kling.
This pipeline processes 16S rRNA gene sequences that were generated using the Illumina MiSeq platform using paired-end sequencing read amplicons.
# Make a vector of required packages
required.packages <- c("agricolae","corrr","data.table","devtools","dplyr","forcats","ggalluvial","ggdendro","ggplot2","ggpubr","grid","gridExtra","knitr","magrittr","microeco","patchwork","pheatmap","pvclust","qiime2R","RColorBrewer","tidyr","UpSetR","vegan")
# Load required packages
lapply(required.packages, library, character.only = TRUE)
Read in the Qiime2 Metadata file.
metadata <- read_q2metadata("QIIME/SILVA/qiime2R/TundraPro.SILVA.Metadata.txt")
rownames(metadata) <- as.character(metadata[, 1])
head(metadata) # show top lines of metadata
SampleID <fctr> | site_tundra <fctr> | site <fctr> | tundra <fctr> | increment <fctr> | depth <fctr> | pH <dbl> | conductivity <dbl> | ||
---|---|---|---|---|---|---|---|---|---|
ITT.0.10 | ITT.0.10 | ITT | Imnavait | Tussock | 0-10 | 0 | 4.61 | 73.7 | |
ITT.10.20 | ITT.10.20 | ITT | Imnavait | Tussock | 10-20 | 10 | 5.02 | 8.3 | |
ITT.20.30 | ITT.20.30 | ITT | Imnavait | Tussock | 20-30 | 20 | 5.28 | 5.9 | |
ITT.30.40 | ITT.30.40 | ITT | Imnavait | Tussock | 30-40 | 30 | 5.28 | 20.9 | |
ITT.40.50 | ITT.40.50 | ITT | Imnavait | Tussock | 40-50 | 40 | 5.48 | 17.4 | |
ITT.50.60 | ITT.50.60 | ITT | Imnavait | Tussock | 50-60 | 50 | 5.65 | 8.7 |
Read in the Qiime2 ASV data table.
ASV <- as.data.frame(read_qza("QIIME/SILVA/qiime2R/table-dada2.qza")$data)
Read in the Qiime2 Taxonomy table.
# Read taxonomy table
taxa_table <- read_qza("QIIME/SILVA/qiime2R/rep-seqs-dada2-taxonomy-SILVA138-FULL.qza")
taxa_table <- parse_taxonomy(taxa_table$data)
# Make the taxonomic table clean, this is very important.
taxa_table %<>% tidy_taxonomy
Create microtable of the imported Qiime2 data files.
dataset <- microtable$new(sample_table = metadata,
tax_table = taxa_table,
otu_table = ASV)
dataset
microtable-class object:
sample_table have 51 rows and 9 columns
otu_table have 25555 rows and 51 columns
tax_table have 25555 rows and 7 columns
Make sure that the data types of sample_table, otu_table and tax_table are all data.frame as the following part shows.
class(ASV)
[1] "data.frame"
ASV[1:5, 1:5]
ITT.0.10 <dbl> | ITT.10.20 <dbl> | ITT.20.30 <dbl> | ITT.30.40 <dbl> | ITT.40.50 <dbl> | |
---|---|---|---|---|---|
75ea2daf1a41ca52cecdd1e27b663173 | 0 | 0 | 3355 | 4534 | 7167 |
b4c1e218247b4ad7d670d7b305c43076 | 11 | 205 | 5010 | 6770 | 2470 |
f34cdcd73e535507d29c29b19d01f658 | 0 | 48 | 2947 | 1962 | 17937 |
f69f5292d08dd7ac727db58db850e379 | 0 | 0 | 868 | 1344 | 2102 |
2d51eee954d47bb2f71665bd2691bb94 | 39 | 8444 | 96 | 118 | 68 |
class(taxa_table)
[1] "data.frame"
taxa_table[1:5, 1:3]
Kingdom <chr> | Phylum <chr> | Class <chr> | |
---|---|---|---|
75ea2daf1a41ca52cecdd1e27b663173 | k__Bacteria | p__Bacteroidota | c__Bacteroidia |
b4c1e218247b4ad7d670d7b305c43076 | k__Bacteria | p__Actinobacteriota | c__Actinobacteria |
f34cdcd73e535507d29c29b19d01f658 | k__Bacteria | p__Proteobacteria | c__Gammaproteobacteria |
f69f5292d08dd7ac727db58db850e379 | k__Bacteria | p__Actinobacteriota | c__Thermoleophilia |
2d51eee954d47bb2f71665bd2691bb94 | k__Bacteria | p__Verrucomicrobiota | c__Verrucomicrobiae |
class(metadata)
[1] "data.frame"
metadata[1:5, ]
SampleID <fctr> | site_tundra <fctr> | site <fctr> | tundra <fctr> | increment <fctr> | depth <fctr> | pH <dbl> | conductivity <dbl> | ||
---|---|---|---|---|---|---|---|---|---|
ITT.0.10 | ITT.0.10 | ITT | Imnavait | Tussock | 0-10 | 0 | 4.61 | 73.7 | |
ITT.10.20 | ITT.10.20 | ITT | Imnavait | Tussock | 10-20 | 10 | 5.02 | 8.3 | |
ITT.20.30 | ITT.20.30 | ITT | Imnavait | Tussock | 20-30 | 20 | 5.28 | 5.9 | |
ITT.30.40 | ITT.30.40 | ITT | Imnavait | Tussock | 30-40 | 30 | 5.28 | 20.9 | |
ITT.40.50 | ITT.40.50 | ITT | Imnavait | Tussock | 40-50 | 40 | 5.48 | 17.4 |
class(dataset)
[1] "microtable" "R6"
print(dataset)
microtable-class object:
sample_table have 51 rows and 9 columns
otu_table have 25555 rows and 51 columns
tax_table have 25555 rows and 7 columns
To make the species and sample information consistent across different files in the dataset object, we can use function tidy_dataset() to trim the dataset.
dataset$tidy_dataset()
print(dataset)
microtable-class object:
sample_table have 51 rows and 9 columns
otu_table have 25555 rows and 51 columns
tax_table have 25555 rows and 7 columns
Remove ASVs which are not assigned in the Kingdom “Archaea” or “Bacteria”.
dataset$tax_table %<>% base::subset(Kingdom == "k__Archaea" | Kingdom == "k__Bacteria")
print(dataset)
microtable-class object:
sample_table have 51 rows and 9 columns
otu_table have 25555 rows and 51 columns
tax_table have 25514 rows and 7 columns
# Remove the lines containing the taxa word regardless of taxonomic ranks and ignoring word case in the tax_table.
dataset$filter_pollution(taxa = c("mitochondria", "chloroplast"))
Total 261 taxa are removed from tax_table ...
print(dataset)
microtable-class object:
sample_table have 51 rows and 9 columns
otu_table have 25555 rows and 51 columns
tax_table have 25253 rows and 7 columns
To make the ASVs same in otu_table and tax_table, we use tidy_dataset() again.
dataset$tidy_dataset()
print(dataset)
microtable-class object:
sample_table have 51 rows and 9 columns
otu_table have 25253 rows and 51 columns
tax_table have 25253 rows and 7 columns
Then we use sample_sums() to check the sequence numbers in each sample.
dataset$sample_sums() %>% range
[1] 50487 166753
dataset$sample_sums() %>% mean
[1] 116273.4
Sometimes, in order to reduce the effects of species number on the diversity measurements, we need to perform the resampling to make the sequence number equal for each sample. The function rarefy_samples can invoke the function tidy_dataset automatically before and after the rarefying.
# As an example, we use 10000 sequences in each sample
dataset$rarefy_samples(sample.size = 50487)
1455 OTUs were removed because they are no longer present in any sample after random subsampling ...
1455 taxa are removed from the otu_table, as the abundance is 0 ...
dataset$sample_sums() %>% range
[1] 50487 50487
Use tidy_dataset() again
dataset$tidy_dataset()
print(dataset)
microtable-class object:
sample_table have 51 rows and 9 columns
otu_table have 23798 rows and 51 columns
tax_table have 23798 rows and 7 columns
Then, we calculate the taxa abundance at each taxonomic rank using cal_abund(). This function return a list called taxa_abund containing several data frames of the abundance information at each taxonomic rank. The list is stored in the microtable object automatically.
dataset$cal_abund()
The result is stored in object$taxa_abund ...
# return dataset$taxa_abund
class(dataset$taxa_abund)
[1] "list"
Then, we calculate the alpha diversity. The result is also stored in the object microtable automatically. As an example, we do not calculate phylogenetic diversity.
# If you want to add Faith's phylogenetic diversity, use PD = TRUE, this will be a little slow
dataset$cal_alphadiv(PD = FALSE)
The result is stored in object$alpha_diversity ...
# return dataset$alpha_diversity
class(dataset$alpha_diversity)
[1] "data.frame"
We also calculate the distance matrix of beta diversity using function cal_betadiv(). We provide four most frequently used indexes: Bray-curtis, Jaccard, weighted Unifrac and unweighted unifrac.
# If you do not want to calculate unifrac metrics, use unifrac = FALSE
# Requires GUniFrac package
dataset$cal_betadiv(unifrac = FALSE)
The result is stored in object$beta_diversity ...
# return dataset$beta_diversity
class(dataset$beta_diversity)
[1] "list"
Clone a copy of the dataset before manipulating.
# Clone the full dataset
dataset.full <- clone(dataset)
print(dataset.full)
microtable-class object:
sample_table have 51 rows and 9 columns
otu_table have 23798 rows and 51 columns
tax_table have 23798 rows and 7 columns
Taxa abundance: calculated for Kingdom,Phylum,Class,Order,Family,Genus,Species
Alpha diversity: calculated for Observed,Chao1,se.chao1,ACE,se.ACE,Shannon,Simpson,InvSimpson,Fisher
Beta diversity: calculated for bray,jaccard
The trans_venn class is used for venn analysis. To analyze the unique and shared OTUs of groups, we first merge samples according to the “Group” column of sample_table.
# merge samples as one community for each group
dataset1 <- dataset.full$merge_samples(use_group = "site")
# dataset1 is a new microtable object
# create trans_venn object
t1 <- trans_venn$new(dataset1, ratio = "seqratio")
The details of each venn part is stored in object$data_details ...
The venn summary table used for plot is stored in object$data_summary ...
t1.plot.venn<-t1$plot_venn()
# The integer data is OTU number
# The percentage data is the sequence number/total sequence number
t1.plot.venn
#Export as .png (width:650, height:650; "venn.site.ASV.png")
# merge samples as one community for each group
dataset2 <- dataset.full$merge_samples(use_group = "tundra")
# dataset1 is a new microtable object
# create trans_venn object
t2 <- trans_venn$new(dataset2, ratio = "seqratio")
The details of each venn part is stored in object$data_details ...
The venn summary table used for plot is stored in object$data_summary ...
t2.plot.venn<-t2$plot_venn()
# The integer data is OTU number
# The percentage data is the sequence number/total sequence number
t2.plot.venn
#Export as .png (width:650, height:650; "venn.site.ASV.png")
#When the groups are too many to show with venn plot, we can use petal plot.
# Use "Type" column in sample_table
dataset3 <- dataset.full$merge_samples(use_group = "site_tundra")
t3 <- trans_venn$new(dataset3)
The details of each venn part is stored in object$data_details ...
The venn summary table used for plot is stored in object$data_summary ...
t3.plot.venn.petal<-t3$plot_venn(petal_plot = TRUE)
t3.plot.venn.petal
#Export as .png (width:1000, height:1000; "venn.petal.site.tundra.ASV.png")
Re-run the Beta Diversity metrics using NMDS ordination on the Bray-Curtis distance.
# we first create an object and select NMDS for ordination
t1 <- trans_beta$new(dataset = dataset.full, group = "site_tundra", measure = "bray")
# Use NMDS as an example, PCA is also available
t1$cal_ordination(ordination = "NMDS")
Run 0 stress 0.1198376
Run 1 stress 0.1198802
... Procrustes: rmse 0.009937913 max resid 0.04644796
Run 2 stress 0.1200816
... Procrustes: rmse 0.007228646 max resid 0.04866366
Run 3 stress 0.1584234
Run 4 stress 0.1198802
... Procrustes: rmse 0.009951994 max resid 0.04646245
Run 5 stress 0.1195786
... New best solution
... Procrustes: rmse 0.007204546 max resid 0.04643815
Run 6 stress 0.1195786
... New best solution
... Procrustes: rmse 3.916508e-06 max resid 1.073779e-05
... Similar to previous best
Run 7 stress 0.1555767
Run 8 stress 0.1195786
... New best solution
... Procrustes: rmse 2.517078e-06 max resid 1.26938e-05
... Similar to previous best
Run 9 stress 0.1195786
... Procrustes: rmse 5.207424e-06 max resid 2.890845e-05
... Similar to previous best
Run 10 stress 0.1195786
... Procrustes: rmse 7.424317e-06 max resid 4.128792e-05
... Similar to previous best
Run 11 stress 0.1198376
... Procrustes: rmse 0.007204105 max resid 0.0464588
Run 12 stress 0.1568478
Run 13 stress 0.1198802
... Procrustes: rmse 0.007040114 max resid 0.04743517
Run 14 stress 0.1195786
... Procrustes: rmse 4.295764e-06 max resid 2.387066e-05
... Similar to previous best
Run 15 stress 0.1195786
... Procrustes: rmse 6.202418e-06 max resid 3.347449e-05
... Similar to previous best
Run 16 stress 0.1195786
... Procrustes: rmse 4.144721e-06 max resid 2.278771e-05
... Similar to previous best
Run 17 stress 0.1195786
... Procrustes: rmse 9.206344e-06 max resid 3.218291e-05
... Similar to previous best
Run 18 stress 0.1195786
... Procrustes: rmse 4.58643e-06 max resid 1.349492e-05
... Similar to previous best
Run 19 stress 0.1195786
... Procrustes: rmse 1.371101e-05 max resid 4.871131e-05
... Similar to previous best
Run 20 stress 0.1198802
... Procrustes: rmse 0.007037372 max resid 0.04741693
*** Solution reached
The ordination result is stored in object$res_ordination ...
# t1$res_ordination is the ordination result list
class(t1$res_ordination)
[1] "list"
# plot the NMDS result
t1.plot.bray.nmds.site.tundra<-t1$plot_ordination(plot_color = "site_tundra", plot_shape = "site_tundra") + theme_classic() #+ theme(legend.position="bottom") + theme(legend.title = element_blank())
#, plot_group_ellipse = FALSE
t1.plot.bray.nmds.site.tundra
#Export as .png (width:800, height:700; "beta.nmds.site.tundra.png")
#For x- and y-coordinates, use:
#head(t1.plot.bray.nmds.site.tundra)
Then we plot and compare the group distances.
# calculate and plot sample distances within groups
t1$cal_group_distance()
The result is stored in object$res_group_distance ...
# return t1$res_group_distance
t1.plot.bray.anova <- t1$plot_group_distance(distance_pair_stat = TRUE)
The ordered groups are ITT IWS STT SWS TTT TWS ...
t1.plot.bray.anova
#Export as .eps (width:800, height:900; "boxplot.beta.bray.anova.silva.eps")
This class is used to transform taxonomic abundance data for plotting the taxa abundance with the ggplot2 package. We first use this class for the bar plot.
# create trans_abund object using 12 Phyla with the highest abundance in the dataset
t1 <- trans_abund$new(dataset = dataset.full, taxrank = "Phylum", ntaxa = 12)
The transformed abundance data is stored in object$data_abund ...
We remove the sample names in x axis and add the facet to show abundance according to groups.
# Place sites in order for plotting
t1$sample_table$site_tundra <- factor(t1$sample_table$site_tundra, levels = c("TTT","ITT","STT","TWS","IWS","SWS"))
# return a ggplot2 object
t1.plot.bar.v1.depths<-t1$plot_bar(facet = "site_tundra", xtext_type_hor = FALSE, xtext_size = 6) + facet_wrap(~site_tundra, scale="free_x", ncol=3) + theme_classic() + theme(axis.text.x = element_text(size = 6, angle = 45, hjust = 1)) #+ theme(legend.position="bottom")
Scale for 'y' is already present. Adding another scale for 'y', which will
replace the existing scale.
t1.plot.bar.v1.depths
#Export as .png (width:1000, height:800; "barplot.site.tundra.depth.png")
Then alluvial plot is implemented in the plot_bar function.
t1 <- trans_abund$new(dataset = dataset.full, taxrank = "Phylum", ntaxa = 12)
The transformed abundance data is stored in object$data_abund ...
# Place sites in order for plotting
t1$sample_table$site_tundra <- factor(t1$sample_table$site_tundra, levels = c("TTT","ITT","STT","TWS","IWS","SWS"))
# use_alluvium = TRUE make the alluvial plot, clustering = TRUE can be used to reorder the samples by clustering
t1.plot.bar.v3.alluvium<-t1$plot_bar(facet = "site_tundra", use_alluvium = TRUE, clustering = FALSE, xtext_type_hor = FALSE, xtext_size = 6) + facet_wrap(~site_tundra, scale="free_x", ncol=3) + theme_classic() + theme_classic() + theme(axis.text.x = element_text(size = 6, angle = 45, hjust = 1))
Scale for 'y' is already present. Adding another scale for 'y', which will
replace the existing scale.
t1.plot.bar.v3.alluvium
#Export as .png (width:1000, height:800; "barplot.site.tundra.depth.alluvium.png")
Alpha diversity can be transformed and plotted using trans_alpha class. Creating trans_alpha object can return two data frame: alpha_data and alpha_stat.
t1 <- trans_alpha$new(dataset = dataset.full, group = "site_tundra")
The transformed diversity data is stored in object$data_alpha ...
The group statistics are stored in object$data_stat ...
# return t1$alpha_stat
t1$alpha_stat[1:5, ]
NULL
Then, we test the differences among groups using anova with multiple comparisons.
t1$cal_diff(method = "anova")
The result is stored in object$res_diff ...
# return t1$res_alpha_diff
t1$res_alpha_diff
NULL
Now, let us plot the mean and se of alpha diversity for each group, and add the duncan.test (agricolae package) result.
t1.plot.alpha.chao1 <- t1$plot_alpha(add_letter = TRUE, measure = "Chao1")
t1.plot.alpha.chao1
#Export as .png (width:800, height:900; "dotplot.alpha.chao1.silva.png")
We can also use the boxplot to show the paired comparisons directly.
t1.plot.alpha.chao1.compare <- t1$plot_alpha(pair_compare = TRUE, measure = "Chao1")
t1.plot.alpha.chao1.compare
#Export as .png (width:800, height:600; "boxplot.alpha.chao1.site.tundra.png")
t1.plot.alpha.shannon <- t1$plot_alpha(add_letter = TRUE, measure = "Shannon")
t1.plot.alpha.shannon
#Export as .png (width:800, height:900; "dotplot.alpha.shannon.silva.png")
The distance matrix of beta diversity can be transformed and plotted using trans_beta class. The analysis referred to the beta diversity in this class mainly include ordination, group distance, clustering and manova. We first show the ordination using PCoA.
# we first create an object and select PCoA for ordination
t1 <- trans_beta$new(dataset = dataset.full, group = "site_tundra", measure = "bray")
# Use PCoA
t1$cal_ordination(ordination = "PCoA")
The ordination result is stored in object$res_ordination ...
# t1$res_ordination is the ordination result list
class(t1$res_ordination)
[1] "list"
# plot the PCoA result by site
t1.plot.bray.pcoa.site <- t1$plot_ordination(plot_color = "site",
plot_shape = "site") +
theme_classic()
#, plot_group_ellipse = TRUE
t1.plot.bray.pcoa.site
#Export as .png (width:800, height:700; "beta.bray.pcoa.site.png")
# plot the PCoA result by site_tundra
t1.plot.bray.pcoa.site.tundra <- t1$plot_ordination(plot_color = "site_tundra",
plot_shape = "site_tundra") + theme_classic()
#, plot_group_ellipse = TRUE
t1.plot.bray.pcoa.site.tundra
#Export as .png (width:800, height:700; "beta.pcoa.site.tundra.png")
Then we plot and compare the group distances.
# calculate and plot sample distances within groups
t1$cal_group_distance()
The result is stored in object$res_group_distance ...
# return t1$res_group_distance
t1.plot.bray.anova <- t1$plot_group_distance(distance_pair_stat = TRUE)
The ordered groups are ITT IWS STT SWS TTT TWS ...
t1.plot.bray.anova
#Export as .png (width:800, height:700; "boxplot.beta.bray.anova.site.tundra.png")
Clustering plot is also a frequently used method.
# use replace_name to set the label name, group parameter used to set the color
t1.plot.bray.clustering <- t1$plot_clustering(group = "tundra")
Scale for 'x' is already present. Adding another scale for 'x', which will
replace the existing scale.
t1.plot.bray.clustering
#Export as .png (width:600, height:800; "clustering.beta.bray.site.tundra.png")
Import the rarefied ASV table
# Transpose data
asv.rare <- t(asv.rare)
# Convert first row into column names
colnames(asv.rare) <- asv.rare[1,]
asv.rare <- asv.rare[-1, ]
# Make dataframe
asv.rare <- as.data.frame(asv.rare)
# Convert character columns to numeric
asv.rare[] <- lapply(asv.rare, function(x) as.numeric(as.character(x)))
Import the alpha diversity metrics
# Convert first column to row names
rownames(asv.alpha) <- asv.alpha[,1]
asv.alpha <- asv.alpha[,-1]
# Make dataframe
asv.alpha <- as.data.frame(asv.alpha)
Import the metadata metrics
# Convert first column to row names
rownames(asv.env)<-asv.env[,1]
asv.env<-asv.env[,-1]
# Make dataframe
asv.env <- as.data.frame(asv.env)
ASV Chao1 Diversity
# Alpha diversity for Site x Tundra Interactions - Chao1 Diversity
chao.asv <- aov(Chao1 ~ Site*Tundra, data=asv.alpha)
summary.aov(chao.asv)
Df Sum Sq Mean Sq F value Pr(>F)
Site 2 840172 420086 2.312 0.111
Tundra 1 102932 102932 0.567 0.456
Site:Tundra 2 147918 73959 0.407 0.668
Residuals 45 8175035 181667
# Chao1 post-hoc analysis
TukeyHSD(chao.asv)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Chao1 ~ Site * Tundra, data = asv.alpha)
$Site
diff lwr upr p adj
Sagwon-Imnavait 226.34428 -133.4667 586.15529 0.2891963
Toolik-Imnavait -79.67794 -429.0396 269.68372 0.8456006
Toolik-Sagwon -306.02222 -660.9542 48.90977 0.1033931
$Tundra
diff lwr upr p adj
WS-MAT 89.19561 -151.2676 329.6589 0.4588915
$`Site:Tundra`
diff lwr upr p adj
Sagwon:MAT-Imnavait:MAT 336.34602 -302.8853 975.5773 0.6245422
Toolik:MAT-Imnavait:MAT -92.10043 -731.3318 547.1309 0.9980334
Imnavait:WS-Imnavait:MAT 137.02136 -488.0704 762.1131 0.9861275
Sagwon:WS-Imnavait:MAT 269.14387 -408.8633 947.1511 0.8433217
Toolik:WS-Imnavait:MAT 93.94616 -545.2852 733.1775 0.9978381
Toolik:MAT-Sagwon:MAT -428.44646 -1026.3926 169.4997 0.2897154
Imnavait:WS-Sagwon:MAT -199.32466 -782.1305 383.4812 0.9095425
Sagwon:WS-Sagwon:MAT -67.20215 -706.4335 572.0292 0.9995714
Toolik:WS-Sagwon:MAT -242.39986 -840.3460 355.5463 0.8314365
Imnavait:WS-Toolik:MAT 229.12179 -353.6840 811.9276 0.8485767
Sagwon:WS-Toolik:MAT 361.24431 -277.9870 1000.4756 0.5505566
Toolik:WS-Toolik:MAT 186.04659 -411.8996 783.9927 0.9377073
Sagwon:WS-Imnavait:WS 132.12251 -492.9692 757.2143 0.9882301
Toolik:WS-Imnavait:WS -43.07520 -625.8810 539.7306 0.9999241
Toolik:WS-Sagwon:WS -175.19771 -814.4290 464.0336 0.9631388
ASV Shannon Diversity
# Alpha diversity for Site x Tundra Interactions - Shannon Diversity
shannon.asv <- aov(Shannon ~ Site*Tundra, data=asv.alpha)
summary.aov(shannon.asv)
Df Sum Sq Mean Sq F value Pr(>F)
Site 2 0.683 0.3417 0.570 0.569
Tundra 1 0.016 0.0165 0.027 0.869
Site:Tundra 2 1.383 0.6916 1.154 0.324
Residuals 45 26.964 0.5992
# Shannon post-hoc analysis
TukeyHSD(shannon.asv)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Shannon ~ Site * Tundra, data = asv.alpha)
$Site
diff lwr upr p adj
Sagwon-Imnavait 0.1716064 -0.4818620 0.8250747 0.8008846
Toolik-Imnavait -0.1112901 -0.7457810 0.5232007 0.9054049
Toolik-Sagwon -0.2828965 -0.9275039 0.3617108 0.5412920
$Tundra
diff lwr upr p adj
WS-MAT -0.03565892 -0.4723746 0.4010568 0.8701083
$`Site:Tundra`
diff lwr upr p adj
Sagwon:MAT-Imnavait:MAT 5.821860e-01 -0.5787497 1.7431216 0.6707333
Toolik:MAT-Imnavait:MAT 7.717905e-02 -1.0837566 1.2381147 0.9999551
Imnavait:WS-Imnavait:MAT 3.357135e-01 -0.7995427 1.4709696 0.9493673
Sagwon:WS-Imnavait:MAT 9.509779e-02 -1.1362604 1.3264560 0.9999057
Toolik:WS-Imnavait:MAT 9.519768e-02 -1.0657380 1.2561333 0.9998733
Toolik:MAT-Sagwon:MAT -5.050069e-01 -1.5909628 0.5809489 0.7362060
Imnavait:WS-Sagwon:MAT -2.464725e-01 -1.3049314 0.8119863 0.9818272
Sagwon:WS-Sagwon:MAT -4.870882e-01 -1.6480239 0.6738475 0.8104820
Toolik:WS-Sagwon:MAT -4.869883e-01 -1.5729442 0.5989676 0.7644371
Imnavait:WS-Toolik:MAT 2.585344e-01 -0.7999244 1.3169933 0.9775651
Sagwon:WS-Toolik:MAT 1.791874e-02 -1.1430169 1.1788544 1.0000000
Toolik:WS-Toolik:MAT 1.801863e-02 -1.0679372 1.1039745 1.0000000
Sagwon:WS-Imnavait:WS -2.406157e-01 -1.3758719 0.8946405 0.9880823
Toolik:WS-Imnavait:WS -2.405158e-01 -1.2989746 0.8179431 0.9837014
Toolik:WS-Sagwon:WS 9.989532e-05 -1.1608358 1.1610356 1.0000000
ASV Shannon Diversity – Site x Tundra x Layer
# Alpha diversity for Site x Tundra x Layer Interactions - Chao1 Diversity
chao1.layer.asv <- aov(Chao1 ~ site_tundra*Layer, data=asv.alpha)
summary.aov(chao1.layer.asv)
Df Sum Sq Mean Sq F value Pr(>F)
site_tundra 5 1091023 218205 1.907 0.115
Layer 1 3054165 3054165 26.690 7.4e-06 ***
site_tundra:Layer 5 658109 131622 1.150 0.351
Residuals 39 4462760 114430
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# Alpha diversity for Site x Tundra x Layer Interactions - Shannon Diversity
shannon.layer.asv <- aov(Shannon ~ site_tundra*Layer, data=asv.alpha)
summary.aov(shannon.layer.asv)
Df Sum Sq Mean Sq F value Pr(>F)
site_tundra 5 2.083 0.417 1.232 0.313
Layer 1 11.193 11.193 33.104 1.14e-06 ***
site_tundra:Layer 5 2.586 0.517 1.530 0.203
Residuals 39 13.186 0.338
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# Chao1 post-hoc analysis
TukeyHSD(chao1.layer.asv)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Chao1 ~ site_tundra * Layer, data = asv.alpha)
$site_tundra
diff lwr upr p adj
IWS-ITT 137.02136 -362.4206 636.46327 0.9616685
STT-ITT 336.34602 -174.3933 847.08530 0.3756576
SWS-ITT 269.14387 -272.5769 810.86468 0.6733734
TTT-ITT -92.10043 -602.8397 418.63885 0.9940923
TWS-ITT 93.94616 -416.7931 604.68544 0.9935227
STT-IWS 199.32466 -266.3312 664.98054 0.7925984
SWS-IWS 132.12251 -367.3194 631.56442 0.9671507
TTT-IWS -229.12179 -694.7777 236.53408 0.6821945
TWS-IWS -43.07520 -508.7311 422.58068 0.9997605
SWS-STT -67.20215 -577.9414 443.53713 0.9986718
TTT-STT -428.44646 -906.1993 49.30639 0.1010005
TWS-STT -242.39986 -720.1527 235.35299 0.6537812
TTT-SWS -361.24431 -871.9836 149.49498 0.2989892
TWS-SWS -175.19771 -685.9370 335.54157 0.9057232
TWS-TTT 186.04659 -291.7063 663.79944 0.8496780
$Layer
diff lwr upr p adj
permafrost-active -486.6383 -678.2965 -294.9801 8.2e-06
$`site_tundra:Layer`
diff lwr upr p adj
IWS:active-ITT:active 189.43747 -598.6237 977.49863 0.9993668
STT:active-ITT:active 408.21454 -422.4749 1238.90394 0.8544188
SWS:active-ITT:active 335.65252 -495.0369 1166.34192 0.9559226
TTT:active-ITT:active 123.56248 -707.1269 954.25188 0.9999944
TWS:active-ITT:active 391.94018 -396.1210 1180.00134 0.8447924
ITT:permafrost-ITT:active -272.28581 -1169.5329 624.96129 0.9950316
IWS:permafrost-ITT:active -148.78258 -936.8437 639.27858 0.9999377
STT:permafrost-ITT:active 68.80216 -719.2590 856.86332 1.0000000
SWS:permafrost-ITT:active -91.82013 -989.0672 805.42697 0.9999999
TTT:permafrost-ITT:active -474.67982 -1262.7410 313.38134 0.6310606
TWS:permafrost-ITT:active -541.10769 -1371.7971 289.58171 0.5181171
STT:active-IWS:active 218.77707 -569.2841 1006.83823 0.9976818
SWS:active-IWS:active 146.21505 -641.8461 934.27621 0.9999475
TTT:active-IWS:active -65.87499 -853.9361 722.18617 1.0000000
TWS:active-IWS:active 202.50272 -540.4885 945.49390 0.9980279
ITT:permafrost-IWS:active -461.72327 -1319.6556 396.20905 0.7701097
IWS:permafrost-IWS:active -338.22005 -1081.2112 404.77114 0.9057001
STT:permafrost-IWS:active -120.63531 -863.6265 622.35588 0.9999863
SWS:permafrost-IWS:active -281.25760 -1139.1899 576.67472 0.9906513
TTT:permafrost-IWS:active -664.11728 -1407.1085 78.87390 0.1178231
TWS:permafrost-IWS:active -730.54515 -1518.6063 57.51601 0.0911406
SWS:active-STT:active -72.56202 -903.2514 758.12738 1.0000000
TTT:active-STT:active -284.65206 -1115.3415 546.03734 0.9867477
TWS:active-STT:active -16.27436 -804.3355 771.78680 1.0000000
ITT:permafrost-STT:active -680.50034 -1577.7474 216.74676 0.2964835
IWS:permafrost-STT:active -556.99712 -1345.0583 231.06404 0.3963349
STT:permafrost-STT:active -339.41238 -1127.4735 448.64878 0.9329510
SWS:permafrost-STT:active -500.03467 -1397.2818 397.21243 0.7308222
TTT:permafrost-STT:active -882.89435 -1670.9555 -94.83319 0.0169998
TWS:permafrost-STT:active -949.32223 -1780.0116 -118.63283 0.0137564
TTT:active-SWS:active -212.09004 -1042.7794 618.59936 0.9989002
TWS:active-SWS:active 56.28766 -731.7735 844.34882 1.0000000
ITT:permafrost-SWS:active -607.93832 -1505.1854 289.30878 0.4591977
IWS:permafrost-SWS:active -484.43510 -1272.4963 303.62606 0.6026209
STT:permafrost-SWS:active -266.85036 -1054.9115 521.21080 0.9879221
SWS:permafrost-SWS:active -427.47265 -1324.7198 469.77445 0.8772357
TTT:permafrost-SWS:active -810.33233 -1598.3935 -22.27117 0.0391652
TWS:permafrost-SWS:active -876.76020 -1707.4496 -46.07081 0.0307885
TWS:active-TTT:active 268.37770 -519.6835 1056.43886 0.9873712
ITT:permafrost-TTT:active -395.84828 -1293.0954 501.39882 0.9220319
IWS:permafrost-TTT:active -272.34506 -1060.4062 515.71610 0.9858481
STT:permafrost-TTT:active -54.76032 -842.8215 733.30084 1.0000000
SWS:permafrost-TTT:active -215.38261 -1112.6297 681.86449 0.9993750
TTT:permafrost-TTT:active -598.24229 -1386.3035 189.81887 0.2952435
TWS:permafrost-TTT:active -664.67016 -1495.3596 166.01924 0.2284201
ITT:permafrost-TWS:active -664.22599 -1522.1583 193.70633 0.2693321
IWS:permafrost-TWS:active -540.72276 -1283.7139 202.26842 0.3538641
STT:permafrost-TWS:active -323.13802 -1066.1292 419.85316 0.9286862
SWS:permafrost-TWS:active -483.76032 -1341.6926 374.17201 0.7168151
TTT:permafrost-TWS:active -866.62000 -1609.6112 -123.62881 0.0109836
TWS:permafrost-TWS:active -933.04787 -1721.1090 -144.98671 0.0092689
IWS:permafrost-ITT:permafrost 123.50323 -734.4291 981.43555 0.9999960
STT:permafrost-ITT:permafrost 341.08796 -516.8444 1199.02029 0.9605462
SWS:permafrost-ITT:permafrost 180.46567 -778.7318 1139.66317 0.9999398
TTT:permafrost-ITT:permafrost -202.39401 -1060.3263 655.53831 0.9994683
TWS:permafrost-ITT:permafrost -268.82188 -1166.0690 628.42522 0.9955380
STT:permafrost-IWS:permafrost 217.58474 -525.4064 960.57592 0.9963225
SWS:permafrost-IWS:permafrost 56.96245 -800.9699 914.89477 1.0000000
TTT:permafrost-IWS:permafrost -325.89723 -1068.8884 417.09395 0.9247895
TWS:permafrost-IWS:permafrost -392.32511 -1180.3863 395.73605 0.8439848
SWS:permafrost-STT:permafrost -160.62229 -1018.5546 697.31003 0.9999426
TTT:permafrost-STT:permafrost -543.48197 -1286.4732 199.50921 0.3466045
TWS:permafrost-STT:permafrost -609.90984 -1397.9710 178.15131 0.2697991
TTT:permafrost-SWS:permafrost -382.85968 -1240.7920 475.07264 0.9163516
TWS:permafrost-SWS:permafrost -449.28755 -1346.5347 447.95955 0.8391373
TWS:permafrost-TTT:permafrost -66.42787 -854.4890 721.63329 1.0000000
# Shannon post-hoc analysis
TukeyHSD(shannon.layer.asv)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Shannon ~ site_tundra * Layer, data = asv.alpha)
$site_tundra
diff lwr upr p adj
IWS-ITT 3.357135e-01 -0.5227881 1.1942150 0.8474561
STT-ITT 5.821860e-01 -0.2957348 1.4601068 0.3680234
SWS-ITT 9.509779e-02 -0.8360779 1.0262734 0.9996113
TTT-ITT 7.717905e-02 -0.8007418 0.9550999 0.9998135
TWS-ITT 9.519768e-02 -0.7827231 0.9731185 0.9994793
STT-IWS 2.464725e-01 -0.5539535 1.0468985 0.9382947
SWS-IWS -2.406157e-01 -1.0991172 0.6178859 0.9580716
TTT-IWS -2.585344e-01 -1.0589604 0.5418916 0.9253558
TWS-IWS -2.405158e-01 -1.0409418 0.5599102 0.9441140
SWS-STT -4.870882e-01 -1.3650090 0.3908326 0.5639412
TTT-STT -5.050069e-01 -1.3262267 0.3162128 0.4516931
TWS-STT -4.869883e-01 -1.3082080 0.3342314 0.4919634
TTT-SWS -1.791874e-02 -0.8958396 0.8600021 0.9999999
TWS-SWS 9.989532e-05 -0.8778209 0.8780207 1.0000000
TWS-TTT 1.801863e-02 -0.8032011 0.8392384 0.9999998
$Layer
diff lwr upr p adj
permafrost-active -0.93159 -1.261035 -0.6021446 1.3e-06
$`site_tundra:Layer`
diff lwr upr p adj
IWS:active-ITT:active 0.341805387 -1.0128100 1.69642080 0.9990118
STT:active-ITT:active 0.772854716 -0.6550353 2.20074474 0.7639526
SWS:active-ITT:active 0.165515096 -1.2623749 1.59340512 0.9999996
TTT:active-ITT:active 0.569270216 -0.8586198 1.99716024 0.9597822
TWS:active-ITT:active 0.623228671 -0.7313867 1.97784409 0.8995751
ITT:permafrost-ITT:active -0.555502310 -2.0977998 0.98679519 0.9806176
IWS:permafrost-ITT:active -0.146523297 -1.5011387 1.20809212 0.9999998
STT:permafrost-ITT:active 0.001120652 -1.3534948 1.35573607 1.0000000
SWS:permafrost-ITT:active -0.554294270 -2.0965918 0.98800323 0.9809313
TTT:permafrost-ITT:active -0.745024243 -2.0996397 0.60959118 0.7461180
TWS:permafrost-ITT:active -1.100503998 -2.5283940 0.32738603 0.2752019
STT:active-IWS:active 0.431049330 -0.9235661 1.78566475 0.9926497
SWS:active-IWS:active -0.176290291 -1.5309057 1.17832513 0.9999986
TTT:active-IWS:active 0.227464829 -1.1271506 1.58208025 0.9999807
TWS:active-IWS:active 0.281423284 -0.9957204 1.55856695 0.9997203
ITT:permafrost-IWS:active -0.897307696 -2.3720262 0.57741078 0.6170242
IWS:permafrost-IWS:active -0.488328684 -1.7654723 0.78881498 0.9699901
STT:permafrost-IWS:active -0.340684735 -1.6178284 0.93645893 0.9983699
SWS:permafrost-IWS:active -0.896099657 -2.3708181 0.57861882 0.6189060
TTT:permafrost-IWS:active -1.086829629 -2.3639733 0.19031403 0.1614567
TWS:permafrost-IWS:active -1.442309384 -2.7969248 -0.08769397 0.0283291
SWS:active-STT:active -0.607339621 -2.0352296 0.82055040 0.9381069
TTT:active-STT:active -0.203584500 -1.6314745 1.22430552 0.9999964
TWS:active-STT:active -0.149626045 -1.5042415 1.20498937 0.9999998
ITT:permafrost-STT:active -1.328357026 -2.8706545 0.21394047 0.1499535
IWS:permafrost-STT:active -0.919378013 -2.2739934 0.43523740 0.4566773
STT:permafrost-STT:active -0.771734064 -2.1263495 0.58288135 0.7041754
SWS:permafrost-STT:active -1.327148987 -2.8694465 0.21514851 0.1508050
TTT:permafrost-STT:active -1.517878959 -2.8724944 -0.16326354 0.0169702
TWS:permafrost-STT:active -1.873358714 -3.3012487 -0.44546869 0.0025854
TTT:active-SWS:active 0.403755120 -1.0241349 1.83164515 0.9972794
TWS:active-SWS:active 0.457713575 -0.8969018 1.81232899 0.9881231
ITT:permafrost-SWS:active -0.721017405 -2.2633149 0.82128009 0.8897730
IWS:permafrost-SWS:active -0.312038393 -1.6666538 1.04257703 0.9995745
STT:permafrost-SWS:active -0.164394444 -1.5190099 1.19022097 0.9999993
SWS:permafrost-SWS:active -0.719809366 -2.2621069 0.82248813 0.8908342
TTT:permafrost-SWS:active -0.910539338 -2.2651548 0.44407608 0.4711549
TWS:permafrost-SWS:active -1.266019093 -2.6939091 0.16187093 0.1244017
TWS:active-TTT:active 0.053958455 -1.3006570 1.40857387 1.0000000
ITT:permafrost-TTT:active -1.124772526 -2.6670700 0.41752497 0.3508847
IWS:permafrost-TTT:active -0.715793513 -2.0704089 0.63882191 0.7892176
STT:permafrost-TTT:active -0.568149564 -1.9227650 0.78646585 0.9434910
SWS:permafrost-TTT:active -1.123564486 -2.6658620 0.41873301 0.3524193
TTT:permafrost-TTT:active -1.314294459 -2.6689099 0.04032096 0.0642265
TWS:permafrost-TTT:active -1.669774214 -3.0976642 -0.24188419 0.0106707
ITT:permafrost-TWS:active -1.178730981 -2.6534495 0.29598750 0.2296892
IWS:permafrost-TWS:active -0.769751968 -2.0468956 0.50739170 0.6302029
STT:permafrost-TWS:active -0.622108019 -1.8992517 0.65503565 0.8612003
SWS:permafrost-TWS:active -1.177522941 -2.6522414 0.29719554 0.2309161
TTT:permafrost-TWS:active -1.368252914 -2.6453966 -0.09110925 0.0266900
TWS:permafrost-TWS:active -1.723732669 -3.0783481 -0.36911725 0.0038604
IWS:permafrost-ITT:permafrost 0.408979013 -1.0657395 1.88369749 0.9977028
STT:permafrost-ITT:permafrost 0.556622962 -0.9180955 2.03134144 0.9726513
SWS:permafrost-ITT:permafrost 0.001208040 -1.6475773 1.64999342 1.0000000
TTT:permafrost-ITT:permafrost -0.189521933 -1.6642404 1.28519654 0.9999988
TWS:permafrost-ITT:permafrost -0.545001688 -2.0872992 0.99729581 0.9832165
STT:permafrost-IWS:permafrost 0.147643949 -1.1294997 1.42478761 0.9999996
SWS:permafrost-IWS:permafrost -0.407770973 -1.8824894 1.06694750 0.9977615
TTT:permafrost-IWS:permafrost -0.598500946 -1.8756446 0.67864272 0.8882324
TWS:permafrost-IWS:permafrost -0.953980701 -2.3085961 0.40063472 0.4016596
SWS:permafrost-STT:permafrost -0.555414922 -2.0301334 0.91930355 0.9730770
TTT:permafrost-STT:permafrost -0.746144895 -2.0232886 0.53099877 0.6720819
TWS:permafrost-STT:permafrost -1.101624650 -2.4562401 0.25299077 0.2095087
TTT:permafrost-SWS:permafrost -0.190729973 -1.6654484 1.28398850 0.9999987
TWS:permafrost-SWS:permafrost -0.546209727 -2.0885072 0.99608777 0.9829320
TWS:permafrost-TTT:permafrost -0.355479755 -1.7100952 0.99913566 0.9985917
ASV Beta diversity (Bray-Curtis)
#Create Bray-Curtis distance matrix for ASV community
asv.bc.dist <- vegdist(asv.rare, method = "bray", binary = FALSE)
PERMANOVA - Site effects on ASV beta diversity
adonis.asv.bc.dist <- adonis2(asv.bc.dist ~ Site*Tundra, data = asv.env)
adonis.asv.bc.dist
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
adonis2(formula = asv.bc.dist ~ Site * Tundra, data = asv.env)
Df SumOfSqs R2 F Pr(>F)
Site 2 2.7813 0.14468 4.7371 0.001 ***
Tundra 1 1.7665 0.09189 6.0172 0.001 ***
Site:Tundra 2 1.4654 0.07623 2.4959 0.001 ***
Residual 45 13.2107 0.68720
Total 50 19.2239 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
PERMANOVA - Site effects on ASV beta diversity
adonis.asv.layer.bc.dist <- adonis2(asv.bc.dist ~ site_tundra*Layer, data = asv.env)
adonis.asv.layer.bc.dist
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
adonis2(formula = asv.bc.dist ~ site_tundra * Layer, data = asv.env)
Df SumOfSqs R2 F Pr(>F)
site_tundra 5 6.0132 0.31280 5.7304 0.001 ***
Layer 1 1.9988 0.10398 9.5241 0.001 ***
site_tundra:Layer 5 3.0269 0.15745 2.8845 0.001 ***
Residual 39 8.1850 0.42577
Total 50 19.2239 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
NMDS Ordination
Plotting NMDS Ordination from microeco results
# Plotting bacteria NMDS in ggplot2
asv.ggplot2 <- ggplot(microeco.nmds, aes(x=NMDS1, y=NMDS2)) +
aes(shape=factor(Tundra)) +
geom_point(size=5, aes(fill=factor(Site))) +
scale_shape_manual("Tundra", values=c(21,24)) +
labs(fill="Site") +
theme_bw() +
scale_x_continuous(limits=c(-2.5,2)) +
scale_y_continuous(limits=c(-2.5, 2)) +
theme(legend.key=element_blank()) +
theme(axis.ticks.y = element_blank(),
axis.ticks.x = element_blank(),
axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14),
panel.background = element_blank(),
axis.title.x = element_text(size = 16),
axis.title.y = element_text(size = 16),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.text = element_text(size=14),
legend.title = element_text(size=16),
panel.border = element_rect(colour = "black", fill=NA, size=2))
asv.ggplot2
#Export as .eps (width:800, height:700; "NMDS.ASV.eps")
Toolik MAT Plotting
# Place taxa in order for plotting
ttt.phylum$Phylum <- factor(ttt.phylum$Phylum,levels = c("Acidobacteriota", "Actinobacteriota", "Bacteroidota", "Caldisericota", "Chloroflexi", "Desulfobacterota", "Firmicutes", "Gemmatimonadota", "Myxococcota", "Patescibacteria", "Planctomycetota", "Verrucomicrobiota", "Alphaproteobacteria", "Gammaproteobacteria", "Bacteria_Other", "Archaea"))
ttt.phylum$Phylum <- fct_rev(ttt.phylum$Phylum)
ttt.phylum$Depth <- factor(ttt.phylum$Depth,levels = c("80-90", "70-80", "60-70", "50-60", "40-50", "30-40", "20-30", "10-20", "0-10"))
colourCount = length(unique(ttt.phylum$Phylum))
getPalette = colorRampPalette(brewer.pal(12, "Paired"))
ttt.phylum.plot <- ggplot(ttt.phylum, aes(fill=Phylum, y=Value, x=Depth)) +
geom_bar(position = "stack", stat="identity", color=NA) +
coord_flip() +
xlab(expression(atop("Toolik MAT", paste("Soil Depth (cm)")))) +
ylab("Relative Abundance (%)") +
theme_minimal() +
theme(axis.line = element_line(colour = 'black', size = 1),
axis.ticks = element_line(colour = "black", size = 2),
plot.title = element_text(size = 18, hjust = 0.5),
axis.text.x = element_text(size = 16),
axis.text.y = element_text(size = 16),
panel.background = element_blank(),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
theme(legend.position = "none") +
scale_fill_manual(values = rev(getPalette(colourCount))) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 102)) +
guides(fill = guide_legend(reverse = TRUE))
ttt.phylum.plot
# Save as .png width=650, height=650 ("ttt.phyla.silva.png")
Toolik WS Plotting
# Place taxa in order for plotting
tws.phylum$Phylum <- factor(tws.phylum$Phylum,levels = c("Acidobacteriota", "Actinobacteriota", "Bacteroidota", "Caldisericota", "Chloroflexi", "Desulfobacterota", "Firmicutes", "Gemmatimonadota", "Myxococcota", "Patescibacteria", "Planctomycetota", "Verrucomicrobiota", "Alphaproteobacteria", "Gammaproteobacteria", "Bacteria_Other", "Archaea"))
tws.phylum$Phylum<-fct_rev(tws.phylum$Phylum)
tws.phylum$Depth<-factor(tws.phylum$Depth,levels = c("80-90", "70-80", "60-70", "50-60", "40-50", "30-40", "20-30", "10-20", "0-10"))
colourCount = length(unique(tws.phylum$Phylum))
getPalette = colorRampPalette(brewer.pal(12, "Paired"))
tws.phylum.plot<-ggplot(tws.phylum, aes(fill=Phylum, y=Value, x=Depth)) +
geom_bar(position = "stack", stat="identity", color=NA) +
coord_flip() +
xlab(expression(atop("Toolik WS", paste("Soil Depth (cm)")))) +
ylab("Relative Abundance (%)") +
theme_minimal() +
theme(axis.line = element_line(colour = 'black', size = 1),
axis.ticks = element_line(colour = "black", size = 2),
plot.title = element_text(size = 18, hjust = 0.5),
axis.text.x = element_text(size = 16),
axis.text.y = element_text(size = 16),
panel.background = element_blank(),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
theme(legend.position = "none") +
scale_fill_manual(values = rev(getPalette(colourCount))) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 102)) +
guides(fill = guide_legend(reverse = TRUE))
tws.phylum.plot
# Save as .png width=650, height=650 ("tws.phyla.silva.png")
Sagwon MAT Plotting
# Place taxa in order for plotting
stt.phylum$Phylum <- factor(stt.phylum$Phylum,levels = c("Acidobacteriota", "Actinobacteriota", "Bacteroidota", "Caldisericota", "Chloroflexi", "Desulfobacterota", "Firmicutes", "Gemmatimonadota", "Myxococcota", "Patescibacteria", "Planctomycetota", "Verrucomicrobiota", "Alphaproteobacteria", "Gammaproteobacteria", "Bacteria_Other", "Archaea"))
stt.phylum$Phylum <- fct_rev(stt.phylum$Phylum)
stt.phylum$Depth <- factor(stt.phylum$Depth,levels = c("80-90", "70-80", "60-70", "50-60", "40-50", "30-40", "20-30", "10-20", "0-10"))
colourCount = length(unique(stt.phylum$Phylum))
getPalette = colorRampPalette(brewer.pal(12, "Paired"))
stt.phylum.plot <- ggplot(stt.phylum, aes(fill=Phylum, y=Value, x=Depth)) +
geom_bar(position = "stack", stat="identity", color=NA) +
coord_flip() +
xlab(expression(atop("Sagwon MAT", paste("Soil Depth (cm)")))) +
ylab("Relative Abundance (%)") +
theme_minimal() +
theme(axis.line = element_line(colour = 'black', size = 1),
axis.ticks = element_line(colour = "black", size = 2),
plot.title = element_text(size = 18, hjust = 0.5),
axis.text.x = element_text(size = 16),
axis.text.y = element_text(size = 16),
panel.background = element_blank(),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
theme(legend.position = "none") +
scale_fill_manual(values = rev(getPalette(colourCount))) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 102)) +
guides(fill = guide_legend(reverse = TRUE))
stt.phylum.plot
# Save as .png width=650, height=650 ("stt.phyla.silva.png")
Sagwon WS Plotting
# Place taxa in order for plotting
sws.phylum$Phylum <- factor(sws.phylum$Phylum,levels = c("Acidobacteriota", "Actinobacteriota", "Bacteroidota", "Caldisericota", "Chloroflexi", "Desulfobacterota", "Firmicutes", "Gemmatimonadota", "Myxococcota", "Patescibacteria", "Planctomycetota", "Verrucomicrobiota", "Alphaproteobacteria", "Gammaproteobacteria", "Bacteria_Other", "Archaea"))
sws.phylum$Phylum <- fct_rev(sws.phylum$Phylum)
sws.phylum$Depth <- factor(sws.phylum$Depth,levels = c("60-70", "50-60", "40-50", "30-40", "20-30", "10-20", "0-10"))
colourCount = length(unique(sws.phylum$Phylum))
getPalette = colorRampPalette(brewer.pal(12, "Paired"))
sws.phylum.plot <- ggplot(sws.phylum, aes(fill=Phylum, y=Value, x=Depth)) +
geom_bar(position = "stack", stat="identity", color=NA) +
coord_flip() +
xlab(expression(atop("Sagwon WS", paste("Soil Depth (cm)")))) +
ylab("Relative Abundance (%)") +
theme_minimal() +
theme(axis.line = element_line(colour = 'black', size = 1),
axis.ticks = element_line(colour = "black", size = 2),
plot.title = element_text(size = 18, hjust = 0.5),
axis.text.x = element_text(size = 16),
axis.text.y = element_text(size = 16),
panel.background = element_blank(),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.text = element_text(size=18)) +
theme(legend.position = "none") +
scale_fill_manual(values = rev(getPalette(colourCount))) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 102)) +
guides(fill = guide_legend(reverse = TRUE))
sws.phylum.plot
# Save as .png width=900, height=650 ("sws.phyla.silva.png")
Save the taxonomy stackplots together
Import the taxonomy relative abundance at the phylum-class level
# Convert first column to row names
rownames(asv.taxa)<-asv.taxa[,1]
asv.taxa<-asv.taxa[,-1]
# Make dataframe
asv.taxa <- as.data.frame(asv.taxa)
Dominant Taxa (Order-Level) by Soil Layer – TTT – UPDATED
# Taxa ANOVA with Post-Hoc for Soil Layer Differences
#Acidobacteriales
ttt.taxa1<-aov(Acidobacteriales~Layer, data=aov.ttt.data)
#summary(ttt.taxa1)
TukeyHSD(ttt.taxa1)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Acidobacteriales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL -7.456969 -10.6635 -4.250439 0.0009074
#Solibacterales
ttt.taxa2<-aov(Solibacterales~Layer, data=aov.ttt.data)
#summary(ttt.taxa2)
TukeyHSD(ttt.taxa2)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Solibacterales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL -1.058688 -1.810467 -0.3069095 0.0125916
#Vicinamibacterales
ttt.taxa3<-aov(Vicinamibacterales~Layer, data=aov.ttt.data)
#summary(ttt.taxa3)
TukeyHSD(ttt.taxa3)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Vicinamibacterales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL -0.6191693 -2.343538 1.105199 0.4239269
#Rhizobiales
ttt.taxa4<-aov(Rhizobiales~Layer, data=aov.ttt.data)
#summary(ttt.taxa4)
TukeyHSD(ttt.taxa4)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Rhizobiales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL -3.665894 -4.956166 -2.375623 0.0002728
#Burkholderiales
ttt.taxa5<-aov(Burkholderiales~Layer, data=aov.ttt.data)
#summary(ttt.taxa5)
TukeyHSD(ttt.taxa5)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Burkholderiales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL -1.240914 -7.843932 5.362105 0.6701771
#Geobacterales
ttt.taxa6<-aov(Geobacterales~Layer, data=aov.ttt.data)
#summary(ttt.taxa6)
TukeyHSD(ttt.taxa6)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Geobacterales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL 0.6525442 -0.2396751 1.544764 0.1273576
#Syntrophales
ttt.taxa7<-aov(Syntrophales~Layer, data=aov.ttt.data)
#summary(ttt.taxa7)
TukeyHSD(ttt.taxa7)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Syntrophales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL 3.020381 2.480154 3.560609 3.3e-06
#Gemmatimonadales
ttt.taxa8<-aov(Gemmatimonadales~Layer, data=aov.ttt.data)
#summary(ttt.taxa8)
TukeyHSD(ttt.taxa8)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Gemmatimonadales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL -1.236358 -3.372046 0.8993301 0.2133397
#Myxococcales
ttt.taxa9<-aov(Myxococcales~Layer, data=aov.ttt.data)
#summary(ttt.taxa9)
TukeyHSD(ttt.taxa9)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Myxococcales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL -0.9077584 -1.89774 0.08222323 0.0667914
#Tepidisphaerales
ttt.taxa10<-aov(Tepidisphaerales~Layer, data=aov.ttt.data)
#summary(ttt.taxa10)
TukeyHSD(ttt.taxa10)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Tepidisphaerales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL -1.318062 -2.582066 -0.05405779 0.0431004
#Chthoniobacterales
ttt.taxa11<-aov(Chthoniobacterales~Layer, data=aov.ttt.data)
#summary(ttt.taxa11)
TukeyHSD(ttt.taxa11)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Chthoniobacterales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL -8.048904 -12.29545 -3.802355 0.0028597
#Pedosphaerales
ttt.taxa12<-aov(Pedosphaerales~Layer, data=aov.ttt.data)
#summary(ttt.taxa12)
TukeyHSD(ttt.taxa12)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Pedosphaerales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL -1.76075 -4.082526 0.5610259 0.1160257
#Gaiellales
ttt.taxa13<-aov(Gaiellales~Layer, data=aov.ttt.data)
#summary(ttt.taxa13)
TukeyHSD(ttt.taxa13)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Gaiellales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL 11.06037 5.252611 16.86813 0.0027876
#Micrococcales
ttt.taxa14<-aov(Micrococcales~Layer, data=aov.ttt.data)
#summary(ttt.taxa14)
TukeyHSD(ttt.taxa14)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Micrococcales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL 5.89221 1.739068 10.04535 0.012171
#Solirubrobacterales
ttt.taxa15<-aov(Solirubrobacterales~Layer, data=aov.ttt.data)
#summary(ttt.taxa15)
TukeyHSD(ttt.taxa15)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Solirubrobacterales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL -1.280231 -2.008927 -0.5515344 0.0042718
#RBG.16.55.12
ttt.taxa16<-aov(RBG.16.55.12~Layer, data=aov.ttt.data)
#summary(ttt.taxa16)
TukeyHSD(ttt.taxa16)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = RBG.16.55.12 ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL 5.141522 1.446805 8.836239 0.0132919
#WCHB1.81
ttt.taxa17<-aov(WCHB1.81~Layer, data=aov.ttt.data)
#summary(ttt.taxa17)
TukeyHSD(ttt.taxa17)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = WCHB1.81 ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL 0.7922832 0.2933988 1.291167 0.0071187
#Bacteroidales
ttt.taxa18<-aov(Bacteroidales~Layer, data=aov.ttt.data)
#summary(ttt.taxa18)
TukeyHSD(ttt.taxa18)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Bacteroidales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL 15.47398 9.237419 21.71055 0.0006198
#Chitinophagales
ttt.taxa19<-aov(Chitinophagales~Layer, data=aov.ttt.data)
#summary(ttt.taxa19)
TukeyHSD(ttt.taxa19)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Chitinophagales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL -1.55238 -3.42275 0.3179902 0.0904719
#Sphingobacteriales
ttt.taxa20<-aov(Sphingobacteriales~Layer, data=aov.ttt.data)
#summary(ttt.taxa20)
TukeyHSD(ttt.taxa20)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Sphingobacteriales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL -2.209084 -4.684698 0.2665307 0.0727838
#Caldisericales
ttt.taxa21<-aov(Caldisericales~Layer, data=aov.ttt.data)
#summary(ttt.taxa21)
TukeyHSD(ttt.taxa21)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Caldisericales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL 5.000495 -2.886344 12.88733 0.1774896
#Clostridiales
ttt.taxa22<-aov(Clostridiales~Layer, data=aov.ttt.data)
#summary(ttt.taxa22)
TukeyHSD(ttt.taxa22)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Clostridiales ~ Layer, data = aov.ttt.data)
$Layer
diff lwr upr p adj
PF-AL 7.461129 3.779817 11.14244 0.0019833
ANOVA – Dominant Taxa (Order-Level) by Soil Layer – TWS – UPDATED
# Taxa ANOVA with Post-Hoc for Soil Layer Differences
#Acidobacteriales
tws.taxa1<-aov(Acidobacteriales~Layer, data=aov.tws.data)
#summary(tws.taxa1)
TukeyHSD(tws.taxa1)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Acidobacteriales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL -1.027393 -3.295293 1.240507 0.3196124
#Solibacterales
tws.taxa2<-aov(Solibacterales~Layer, data=aov.tws.data)
#summary(tws.taxa2)
TukeyHSD(tws.taxa2)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Solibacterales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL -0.5178561 -0.8615218 -0.1741903 0.0091806
#Vicinamibacterales
tws.taxa3<-aov(Vicinamibacterales~Layer, data=aov.tws.data)
#summary(tws.taxa3)
TukeyHSD(tws.taxa3)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Vicinamibacterales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL 0.001683602 -0.5353787 0.5387459 0.9942924
#Rhizobiales
tws.taxa4<-aov(Rhizobiales~Layer, data=aov.tws.data)
#summary(tws.taxa4)
TukeyHSD(tws.taxa4)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Rhizobiales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL -1.122071 -3.516004 1.271861 0.3043403
#Burkholderiales
tws.taxa5<-aov(Burkholderiales~Layer, data=aov.tws.data)
#summary(tws.taxa5)
TukeyHSD(tws.taxa5)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Burkholderiales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL -1.767683 -6.975423 3.440057 0.4485754
#Geobacterales
tws.taxa6<-aov(Geobacterales~Layer, data=aov.tws.data)
#summary(tws.taxa6)
TukeyHSD(tws.taxa6)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Geobacterales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL -1.666667 -3.814981 0.4816481 0.1092256
#Syntrophales
tws.taxa7<-aov(Syntrophales~Layer, data=aov.tws.data)
#summary(tws.taxa7)
TukeyHSD(tws.taxa7)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Syntrophales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL -2.466278 -7.160364 2.227807 0.2541015
#Gemmatimonadales
tws.taxa8<-aov(Gemmatimonadales~Layer, data=aov.tws.data)
#summary(tws.taxa8)
TukeyHSD(tws.taxa8)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Gemmatimonadales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL -0.1250817 -0.2536737 0.00351026 0.054985
#Myxococcales
tws.taxa9<-aov(Myxococcales~Layer, data=aov.tws.data)
#summary(tws.taxa9)
TukeyHSD(tws.taxa9)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Myxococcales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL -0.131519 -0.4664813 0.2034433 0.3840759
#Tepidisphaerales
tws.taxa10<-aov(Tepidisphaerales~Layer, data=aov.tws.data)
#summary(tws.taxa10)
TukeyHSD(tws.taxa10)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Tepidisphaerales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL -0.004159487 -0.2233849 0.2150659 0.9654677
#Chthoniobacterales
tws.taxa11<-aov(Chthoniobacterales~Layer, data=aov.tws.data)
#summary(tws.taxa11)
TukeyHSD(tws.taxa11)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Chthoniobacterales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL -1.9815 -4.627208 0.6642079 0.1198641
#Pedosphaerales
tws.taxa12<-aov(Pedosphaerales~Layer, data=aov.tws.data)
#summary(tws.taxa12)
TukeyHSD(tws.taxa12)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Pedosphaerales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL -4.99277 -8.977154 -1.008387 0.0210132
#Gaiellales
tws.taxa13<-aov(Gaiellales~Layer, data=aov.tws.data)
#summary(tws.taxa13)
TukeyHSD(tws.taxa13)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Gaiellales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL 15.08656 8.442423 21.73069 0.0010423
#Micrococcales
tws.taxa14<-aov(Micrococcales~Layer, data=aov.tws.data)
#summary(tws.taxa14)
TukeyHSD(tws.taxa14)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Micrococcales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL 3.420286 1.043356 5.797217 0.0114024
#Solirubrobacterales
tws.taxa15<-aov(Solirubrobacterales~Layer, data=aov.tws.data)
#summary(tws.taxa15)
TukeyHSD(tws.taxa15)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Solirubrobacterales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL -0.1608335 -0.6005176 0.2788507 0.4157144
#RBG.16.55.12
tws.taxa16<-aov(RBG.16.55.12~Layer, data=aov.tws.data)
#summary(tws.taxa16)
TukeyHSD(tws.taxa16)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = RBG.16.55.12 ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL 4.058768 1.808661 6.308874 0.003722
#WCHB1.81
tws.taxa17<-aov(WCHB1.81~Layer, data=aov.tws.data)
#summary(tws.taxa17)
TukeyHSD(tws.taxa17)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = WCHB1.81 ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL 0.109038 -1.085328 1.303403 0.8352407
#Bacteroidales
tws.taxa18<-aov(Bacteroidales~Layer, data=aov.tws.data)
#summary(tws.taxa18)
TukeyHSD(tws.taxa18)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Bacteroidales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL 19.39054 11.48557 27.2955 0.0006633
#Chitinophagales
tws.taxa19<-aov(Chitinophagales~Layer, data=aov.tws.data)
#summary(tws.taxa19)
TukeyHSD(tws.taxa19)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Chitinophagales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL -0.1558817 -0.5367541 0.2249907 0.3653929
#Sphingobacteriales
tws.taxa20<-aov(Sphingobacteriales~Layer, data=aov.tws.data)
#summary(tws.taxa20)
TukeyHSD(tws.taxa20)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Sphingobacteriales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL -1.53178 -4.245748 1.182188 0.2237786
#Caldisericales
tws.taxa21<-aov(Caldisericales~Layer, data=aov.tws.data)
#summary(tws.taxa21)
TukeyHSD(tws.taxa21)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Caldisericales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL 2.275338 -0.6861649 5.236841 0.1121005
#Clostridiales
tws.taxa22<-aov(Clostridiales~Layer, data=aov.tws.data)
#summary(tws.taxa22)
TukeyHSD(tws.taxa22)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Clostridiales ~ Layer, data = aov.tws.data)
$Layer
diff lwr upr p adj
PF-AL 0.5067641 -2.967018 3.980546 0.7402592
ANOVA – Dominant Taxa (Order-Level) by Soil Layer – STT – UPDATED
# Taxa ANOVA with Post-Hoc for Soil Layer Differences
#Acidobacteriales
stt.taxa1<-aov(Acidobacteriales~Layer, data=aov.stt.data)
#summary(stt.taxa1)
TukeyHSD(stt.taxa1)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Acidobacteriales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL 0.002673956 -0.02616951 0.03151742 0.832738
#Solibacterales
stt.taxa2<-aov(Solibacterales~Layer, data=aov.stt.data)
#summary(stt.taxa2)
TukeyHSD(stt.taxa2)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Solibacterales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL -0.1013132 -0.2164469 0.01382052 0.0759979
#Vicinamibacterales
stt.taxa3<-aov(Vicinamibacterales~Layer, data=aov.stt.data)
#summary(stt.taxa3)
TukeyHSD(stt.taxa3)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Vicinamibacterales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL -4.058074 -7.279447 -0.8367017 0.0205498
#Rhizobiales
stt.taxa4<-aov(Rhizobiales~Layer, data=aov.stt.data)
#summary(stt.taxa4)
TukeyHSD(stt.taxa4)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Rhizobiales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL -2.778636 -6.755491 1.198219 0.1424813
#Burkholderiales
stt.taxa5<-aov(Burkholderiales~Layer, data=aov.stt.data)
#summary(stt.taxa5)
TukeyHSD(stt.taxa5)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Burkholderiales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL 2.355161 -2.122843 6.833164 0.2536599
#Geobacterales
stt.taxa6<-aov(Geobacterales~Layer, data=aov.stt.data)
#summary(stt.taxa6)
TukeyHSD(stt.taxa6)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Geobacterales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL -1.661022 -6.336707 3.014664 0.4286492
#Syntrophales
stt.taxa7<-aov(Syntrophales~Layer, data=aov.stt.data)
#summary(stt.taxa7)
TukeyHSD(stt.taxa7)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Syntrophales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL 1.43839 0.8797081 1.997072 0.0004969
#Gemmatimonadales
stt.taxa8<-aov(Gemmatimonadales~Layer, data=aov.stt.data)
#summary(stt.taxa8)
TukeyHSD(stt.taxa8)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Gemmatimonadales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL -1.221007 -2.622071 0.180056 0.0782799
#Myxococcales
stt.taxa9<-aov(Myxococcales~Layer, data=aov.stt.data)
#summary(stt.taxa9)
TukeyHSD(stt.taxa9)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Myxococcales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL -0.04466496 -0.1685632 0.07923324 0.422174
#Tepidisphaerales
stt.taxa10<-aov(Tepidisphaerales~Layer, data=aov.stt.data)
#summary(stt.taxa10)
TukeyHSD(stt.taxa10)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Tepidisphaerales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL -0.2329312 -0.4301706 -0.03569187 0.026811
#Chthoniobacterales
stt.taxa11<-aov(Chthoniobacterales~Layer, data=aov.stt.data)
#summary(stt.taxa11)
TukeyHSD(stt.taxa11)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Chthoniobacterales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL -9.776279 -14.28449 -5.268065 0.001357
#Pedosphaerales
stt.taxa12<-aov(Pedosphaerales~Layer, data=aov.stt.data)
#summary(stt.taxa12)
TukeyHSD(stt.taxa12)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Pedosphaerales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL -0.624022 -2.751118 1.503074 0.5102239
#Gaiellales
stt.taxa13<-aov(Gaiellales~Layer, data=aov.stt.data)
#summary(stt.taxa13)
TukeyHSD(stt.taxa13)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Gaiellales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL 3.842969 0.8675101 6.818429 0.0184777
#Micrococcales
stt.taxa14<-aov(Micrococcales~Layer, data=aov.stt.data)
#summary(stt.taxa14)
TukeyHSD(stt.taxa14)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Micrococcales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL 6.293897 3.096595 9.4912 0.0023288
#Solirubrobacterales
stt.taxa15<-aov(Solirubrobacterales~Layer, data=aov.stt.data)
#summary(stt.taxa15)
TukeyHSD(stt.taxa15)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Solirubrobacterales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL -1.200804 -2.721644 0.3200359 0.1041314
#RBG.16.55.12
stt.taxa16<-aov(RBG.16.55.12~Layer, data=aov.stt.data)
#summary(stt.taxa16)
TukeyHSD(stt.taxa16)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = RBG.16.55.12 ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL 1.522174 -0.2999356 3.344284 0.0887836
#WCHB1.81
stt.taxa17<-aov(WCHB1.81~Layer, data=aov.stt.data)
#summary(stt.taxa17)
TukeyHSD(stt.taxa17)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = WCHB1.81 ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL 2.472815 -0.169687 5.115317 0.06254
#Bacteroidales
stt.taxa18<-aov(Bacteroidales~Layer, data=aov.stt.data)
#summary(stt.taxa18)
TukeyHSD(stt.taxa18)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Bacteroidales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL 4.470557 3.05596 5.885154 0.0001405
#Chitinophagales
stt.taxa19<-aov(Chitinophagales~Layer, data=aov.stt.data)
#summary(stt.taxa19)
TukeyHSD(stt.taxa19)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Chitinophagales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL -1.711035 -2.431396 -0.9906727 0.0008019
#Sphingobacteriales
stt.taxa20<-aov(Sphingobacteriales~Layer, data=aov.stt.data)
#summary(stt.taxa20)
TukeyHSD(stt.taxa20)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Sphingobacteriales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL 0.439222 -0.5406204 1.419064 0.32436
#Caldisericales
stt.taxa21<-aov(Caldisericales~Layer, data=aov.stt.data)
#summary(stt.taxa21)
TukeyHSD(stt.taxa21)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Caldisericales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL 0.5441995 -0.04327552 1.131675 0.064637
#Clostridiales
stt.taxa22<-aov(Clostridiales~Layer, data=aov.stt.data)
#summary(stt.taxa22)
TukeyHSD(stt.taxa22)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Clostridiales ~ Layer, data = aov.stt.data)
$Layer
diff lwr upr p adj
PF-AL 0.9943154 0.601899 1.386732 0.0005468
ANOVA – Dominant Taxa (Order-Level) by Soil Layer – SWS – UPDATED
# Taxa ANOVA with Post-Hoc for Soil Layer Differences
#Acidobacteriales
sws.taxa1<-aov(Acidobacteriales~Layer, data=aov.sws.data)
#summary(sws.taxa1)
TukeyHSD(sws.taxa1)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Acidobacteriales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL -2.031711 -8.137826 4.074403 0.4314411
#Solibacterales
sws.taxa2<-aov(Solibacterales~Layer, data=aov.sws.data)
#summary(sws.taxa2)
TukeyHSD(sws.taxa2)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Solibacterales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL -0.3621394 -0.8571507 0.1328719 0.118788
#Vicinamibacterales
sws.taxa3<-aov(Vicinamibacterales~Layer, data=aov.sws.data)
#summary(sws.taxa3)
TukeyHSD(sws.taxa3)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Vicinamibacterales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL -1.962882 -3.970525 0.04476171 0.0536135
#Rhizobiales
sws.taxa4<-aov(Rhizobiales~Layer, data=aov.sws.data)
#summary(sws.taxa4)
TukeyHSD(sws.taxa4)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Rhizobiales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL 0.4552327 -2.98056 3.891025 0.7472528
#Burkholderiales
sws.taxa5<-aov(Burkholderiales~Layer, data=aov.sws.data)
#summary(sws.taxa5)
TukeyHSD(sws.taxa5)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Burkholderiales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL -8.041344 -15.40156 -0.6811273 0.0376161
#Geobacterales
sws.taxa6<-aov(Geobacterales~Layer, data=aov.sws.data)
#summary(sws.taxa6)
TukeyHSD(sws.taxa6)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Geobacterales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL 0.2114406 -0.0006069774 0.4234881 0.0504487
#Syntrophales
sws.taxa7<-aov(Syntrophales~Layer, data=aov.sws.data)
#summary(sws.taxa7)
TukeyHSD(sws.taxa7)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Syntrophales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL 3.480764 2.922865 4.038663 1.76e-05
#Gemmatimonadales
sws.taxa8<-aov(Gemmatimonadales~Layer, data=aov.sws.data)
#summary(sws.taxa8)
TukeyHSD(sws.taxa8)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Gemmatimonadales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL -0.8177022 -1.876331 0.2409267 0.1038376
#Myxococcales
sws.taxa9<-aov(Myxococcales~Layer, data=aov.sws.data)
#summary(sws.taxa9)
TukeyHSD(sws.taxa9)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Myxococcales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL -0.4235414 -1.012575 0.1654927 0.1238093
#Tepidisphaerales
sws.taxa10<-aov(Tepidisphaerales~Layer, data=aov.sws.data)
#summary(sws.taxa10)
TukeyHSD(sws.taxa10)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Tepidisphaerales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL -1.030133 -3.535851 1.475585 0.338966
#Chthoniobacterales
sws.taxa11<-aov(Chthoniobacterales~Layer, data=aov.sws.data)
#summary(sws.taxa11)
TukeyHSD(sws.taxa11)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Chthoniobacterales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL -12.04369 -30.8618 6.774414 0.1608507
#Pedosphaerales
sws.taxa12<-aov(Pedosphaerales~Layer, data=aov.sws.data)
#summary(sws.taxa12)
TukeyHSD(sws.taxa12)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Pedosphaerales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL -2.806828 -5.291689 -0.3219678 0.0336471
#Gaiellales
sws.taxa13<-aov(Gaiellales~Layer, data=aov.sws.data)
#summary(sws.taxa13)
TukeyHSD(sws.taxa13)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Gaiellales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL 6.228831 2.883229 9.574433 0.0049448
#Micrococcales
sws.taxa14<-aov(Micrococcales~Layer, data=aov.sws.data)
#summary(sws.taxa14)
TukeyHSD(sws.taxa14)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Micrococcales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL 3.199008 0.8956514 5.502365 0.0160413
#Solirubrobacterales
sws.taxa15<-aov(Solirubrobacterales~Layer, data=aov.sws.data)
#summary(sws.taxa15)
TukeyHSD(sws.taxa15)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Solirubrobacterales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL -1.45549 -1.865113 -1.045867 0.0002635
#RBG.16.55.12
sws.taxa16<-aov(RBG.16.55.12~Layer, data=aov.sws.data)
#summary(sws.taxa16)
TukeyHSD(sws.taxa16)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = RBG.16.55.12 ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL 2.860307 2.64874 3.071874 1e-06
#WCHB1.81
sws.taxa17<-aov(WCHB1.81~Layer, data=aov.sws.data)
#summary(sws.taxa17)
TukeyHSD(sws.taxa17)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = WCHB1.81 ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL 10.37495 8.658705 12.09119 2.04e-05
#Bacteroidales
sws.taxa18<-aov(Bacteroidales~Layer, data=aov.sws.data)
#summary(sws.taxa18)
TukeyHSD(sws.taxa18)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Bacteroidales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL 14.68365 11.90644 17.46086 3.87e-05
#Chitinophagales
sws.taxa19<-aov(Chitinophagales~Layer, data=aov.sws.data)
#summary(sws.taxa19)
TukeyHSD(sws.taxa19)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Chitinophagales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL -2.781739 -4.926194 -0.6372842 0.0206767
#Sphingobacteriales
sws.taxa20<-aov(Sphingobacteriales~Layer, data=aov.sws.data)
#summary(sws.taxa20)
TukeyHSD(sws.taxa20)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Sphingobacteriales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL -1.458626 -3.673014 0.7557612 0.1511861
#Caldisericales
sws.taxa21<-aov(Caldisericales~Layer, data=aov.sws.data)
#summary(sws.taxa21)
TukeyHSD(sws.taxa21)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Caldisericales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL 8.54923 4.477254 12.62121 0.0029494
#Clostridiales
sws.taxa22<-aov(Clostridiales~Layer, data=aov.sws.data)
#summary(sws.taxa22)
TukeyHSD(sws.taxa22)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Clostridiales ~ Layer, data = aov.sws.data)
$Layer
diff lwr upr p adj
PF-AL 2.748397 2.277266 3.219529 2.41e-05
Plot Archaea abundance for each site x tundra location
archaea.plot<-archaea.abund
archaea.plot$Increment<-factor(archaea.plot$Increment, levels = c("90-100","80-90","70-80","60-70","50-60","40-50","30-40","20-30","10-20","0-10"))
arch.plot<-ggplot(data=archaea.plot, aes(x=Increment, y=Abundance, group=site_tundra)) +
geom_line(aes(linetype=site_tundra)) +
geom_point(size=5, aes(shape=Tundra, fill=Site)) +
scale_shape_manual("Tundra", values=c(21,24)) +
coord_flip() +
xlab("Soil Depth (cm)") +
ylab("Relative Abundance (%)") +
theme_minimal() +
theme(axis.line = element_line(colour = 'black', size = 1),
axis.ticks = element_line(colour = "black", size = 2),
axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14),
panel.background = element_blank(),
axis.title.x = element_text(size = 16),
axis.title.y = element_text(size = 16),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none") +
scale_y_continuous(expand = c(0, 0), limits = c(0, 3), position = "right") +
guides(linetype = FALSE)
Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
arch.plot
# Save as .eps width=400, height=650 ("archaea.eps")
The session information is provided for full reproducibility.
devtools::session_info()
─ Session info ───────────────────────────────────────────────────────────────────
setting value
version R version 4.2.1 (2022-06-23)
os macOS Monterey 12.6
system x86_64, darwin17.0
ui RStudio
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz America/Los_Angeles
date 2022-09-22
rstudio 2022.07.1+554 Spotted Wakerobin (desktop)
pandoc 2.18 @ /Applications/RStudio.app/Contents/MacOS/quarto/bin/tools/ (via rmarkdown)
─ Packages ───────────────────────────────────────────────────────────────────────
package * version date (UTC) lib source
abind 1.4-5 2016-07-21 [1] CRAN (R 4.2.0)
ade4 1.7-19 2022-04-19 [1] CRAN (R 4.2.0)
agricolae * 1.3-5 2021-06-06 [1] CRAN (R 4.2.0)
AlgDesign 1.2.1 2022-05-25 [1] CRAN (R 4.2.0)
ape 5.6-2 2022-03-02 [1] CRAN (R 4.2.0)
backports 1.4.1 2021-12-13 [1] CRAN (R 4.2.0)
base64enc 0.1-3 2015-07-28 [1] CRAN (R 4.2.0)
Biobase 2.56.0 2022-04-26 [1] Bioconductor
BiocGenerics 0.42.0 2022-04-26 [1] Bioconductor
biomformat 1.24.0 2022-04-26 [1] Bioconductor
Biostrings 2.64.1 2022-08-18 [1] Bioconductor
bitops 1.0-7 2021-04-24 [1] CRAN (R 4.2.0)
broom 1.0.0 2022-07-01 [1] CRAN (R 4.2.0)
cachem 1.0.6 2021-08-19 [1] CRAN (R 4.2.0)
callr 3.7.2 2022-08-22 [1] CRAN (R 4.2.0)
car 3.1-0 2022-06-15 [1] CRAN (R 4.2.0)
carData 3.0-5 2022-01-06 [1] CRAN (R 4.2.0)
checkmate 2.1.0 2022-04-21 [1] CRAN (R 4.2.0)
cli 3.3.0 2022-04-25 [1] CRAN (R 4.2.0)
cluster 2.1.3 2022-03-28 [1] CRAN (R 4.2.1)
codetools 0.2-18 2020-11-04 [1] CRAN (R 4.2.1)
colorspace 2.0-3 2022-02-21 [1] CRAN (R 4.2.0)
combinat 0.0-8 2012-10-29 [1] CRAN (R 4.2.0)
corrr * 0.4.4 2022-08-16 [1] CRAN (R 4.2.0)
crayon 1.5.1 2022-03-26 [1] CRAN (R 4.2.0)
data.table * 1.14.2 2021-09-27 [1] CRAN (R 4.2.0)
deldir 1.0-6 2021-10-23 [1] CRAN (R 4.2.0)
devtools * 2.4.4 2022-07-20 [1] CRAN (R 4.2.0)
digest 0.6.29 2021-12-01 [1] CRAN (R 4.2.0)
dplyr * 1.0.9 2022-04-28 [1] CRAN (R 4.2.0)
DT 0.24 2022-08-09 [1] CRAN (R 4.2.0)
ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.2.0)
evaluate 0.16 2022-08-09 [1] CRAN (R 4.2.0)
fansi 1.0.3 2022-03-24 [1] CRAN (R 4.2.0)
farver 2.1.1 2022-07-06 [1] CRAN (R 4.2.0)
fastmap 1.1.0 2021-01-25 [1] CRAN (R 4.2.0)
forcats * 0.5.2 2022-08-19 [1] CRAN (R 4.2.0)
foreach 1.5.2 2022-02-02 [1] CRAN (R 4.2.0)
foreign 0.8-82 2022-01-16 [1] CRAN (R 4.2.1)
Formula 1.2-4 2020-10-16 [1] CRAN (R 4.2.0)
fs 1.5.2 2021-12-08 [1] CRAN (R 4.2.0)
generics 0.1.3 2022-07-05 [1] CRAN (R 4.2.0)
GenomeInfoDb 1.32.3 2022-08-09 [1] Bioconductor
GenomeInfoDbData 1.2.8 2022-08-29 [1] Bioconductor
ggalluvial * 0.12.3 2020-12-05 [1] CRAN (R 4.2.0)
ggdendro * 0.1.23 2022-02-16 [1] CRAN (R 4.2.0)
ggplot2 * 3.3.6 2022-05-03 [1] CRAN (R 4.2.0)
ggpubr * 0.4.0 2020-06-27 [1] CRAN (R 4.2.0)
ggsignif 0.6.3 2021-09-09 [1] CRAN (R 4.2.0)
glue 1.6.2 2022-02-24 [1] CRAN (R 4.2.0)
gridExtra * 2.3 2017-09-09 [1] CRAN (R 4.2.0)
gtable 0.3.0 2019-03-25 [1] CRAN (R 4.2.0)
haven 2.5.1 2022-08-22 [1] CRAN (R 4.2.0)
highr 0.9 2021-04-16 [1] CRAN (R 4.2.0)
Hmisc 4.7-1 2022-08-15 [1] CRAN (R 4.2.0)
hms 1.1.2 2022-08-19 [1] CRAN (R 4.2.0)
htmlTable 2.4.1 2022-07-07 [1] CRAN (R 4.2.0)
htmltools 0.5.3 2022-07-18 [1] CRAN (R 4.2.0)
htmlwidgets 1.5.4 2021-09-08 [1] CRAN (R 4.2.0)
httpuv 1.6.5 2022-01-05 [1] CRAN (R 4.2.0)
igraph 1.3.4 2022-07-19 [1] CRAN (R 4.2.0)
interp 1.1-3 2022-07-13 [1] CRAN (R 4.2.0)
IRanges 2.30.1 2022-08-18 [1] Bioconductor
iterators 1.0.14 2022-02-05 [1] CRAN (R 4.2.0)
jpeg 0.1-9 2021-07-24 [1] CRAN (R 4.2.0)
jsonlite 1.8.0 2022-02-22 [1] CRAN (R 4.2.0)
klaR 1.7-1 2022-06-27 [1] CRAN (R 4.2.0)
knitr * 1.40 2022-08-24 [1] CRAN (R 4.2.0)
labeling 0.4.2 2020-10-20 [1] CRAN (R 4.2.0)
labelled 2.9.1 2022-05-05 [1] CRAN (R 4.2.0)
later 1.3.0 2021-08-18 [1] CRAN (R 4.2.0)
lattice * 0.20-45 2021-09-22 [1] CRAN (R 4.2.1)
latticeExtra 0.6-30 2022-07-04 [1] CRAN (R 4.2.0)
lifecycle 1.0.1 2021-09-24 [1] CRAN (R 4.2.0)
magrittr * 2.0.3 2022-03-30 [1] CRAN (R 4.2.0)
MASS 7.3-57 2022-04-22 [1] CRAN (R 4.2.1)
Matrix 1.4-1 2022-03-23 [1] CRAN (R 4.2.1)
memoise 2.0.1 2021-11-26 [1] CRAN (R 4.2.0)
mgcv 1.8-40 2022-03-29 [1] CRAN (R 4.2.1)
microeco * 0.11.0 2022-06-22 [1] CRAN (R 4.2.0)
mime 0.12 2021-09-28 [1] CRAN (R 4.2.0)
miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.2.0)
multtest 2.52.0 2022-04-26 [1] Bioconductor
munsell 0.5.0 2018-06-12 [1] CRAN (R 4.2.0)
NADA 1.6-1.1 2020-03-22 [1] CRAN (R 4.2.0)
nlme 3.1-157 2022-03-25 [1] CRAN (R 4.2.1)
nnet 7.3-17 2022-01-16 [1] CRAN (R 4.2.1)
patchwork * 1.1.2 2022-08-19 [1] CRAN (R 4.2.0)
permute * 0.9-7 2022-01-27 [1] CRAN (R 4.2.0)
pheatmap * 1.0.12 2019-01-04 [1] CRAN (R 4.2.0)
phyloseq 1.40.0 2022-04-26 [1] Bioconductor
pillar 1.8.1 2022-08-19 [1] CRAN (R 4.2.0)
pkgbuild 1.3.1 2021-12-20 [1] CRAN (R 4.2.0)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.2.0)
pkgload 1.3.0 2022-06-27 [1] CRAN (R 4.2.0)
plyr 1.8.7 2022-03-24 [1] CRAN (R 4.2.0)
png 0.1-7 2013-12-03 [1] CRAN (R 4.2.0)
prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.2.0)
processx 3.7.0 2022-07-07 [1] CRAN (R 4.2.0)
profvis 0.3.7 2020-11-02 [1] CRAN (R 4.2.0)
promises 1.2.0.1 2021-02-11 [1] CRAN (R 4.2.0)
ps 1.7.1 2022-06-18 [1] CRAN (R 4.2.0)
purrr 0.3.4 2020-04-17 [1] CRAN (R 4.2.0)
pvclust * 2.2-0 2019-11-19 [1] CRAN (R 4.2.0)
qiime2R * 0.99.6 2022-08-29 [1] Github (jbisanz/qiime2R@2a3cee1)
questionr 0.7.7 2022-01-31 [1] CRAN (R 4.2.0)
R6 2.5.1 2021-08-19 [1] CRAN (R 4.2.0)
RColorBrewer * 1.1-3 2022-04-03 [1] CRAN (R 4.2.0)
Rcpp 1.0.9 2022-07-08 [1] CRAN (R 4.2.0)
RCurl 1.98-1.8 2022-07-30 [1] CRAN (R 4.2.0)
remotes 2.4.2 2021-11-30 [1] CRAN (R 4.2.0)
reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.2.0)
rhdf5 2.40.0 2022-04-26 [1] Bioconductor
rhdf5filters 1.8.0 2022-04-26 [1] Bioconductor
Rhdf5lib 1.18.2 2022-05-17 [1] Bioconductor
rlang 1.0.4 2022-07-12 [1] CRAN (R 4.2.0)
rmarkdown 2.16 2022-08-24 [1] CRAN (R 4.2.0)
rpart 4.1.16 2022-01-24 [1] CRAN (R 4.2.1)
rstatix 0.7.0 2021-02-13 [1] CRAN (R 4.2.0)
rstudioapi 0.14 2022-08-22 [1] CRAN (R 4.2.0)
S4Vectors 0.34.0 2022-04-26 [1] Bioconductor
scales 1.2.1 2022-08-20 [1] CRAN (R 4.2.0)
sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.2.0)
shiny 1.7.2 2022-07-19 [1] CRAN (R 4.2.0)
stringi 1.7.8 2022-07-11 [1] CRAN (R 4.2.0)
stringr 1.4.1 2022-08-20 [1] CRAN (R 4.2.0)
survival 3.3-1 2022-03-03 [1] CRAN (R 4.2.1)
tibble 3.1.8 2022-07-22 [1] CRAN (R 4.2.0)
tidyr * 1.2.0 2022-02-01 [1] CRAN (R 4.2.0)
tidyselect 1.1.2 2022-02-21 [1] CRAN (R 4.2.0)
truncnorm 1.0-8 2018-02-27 [1] CRAN (R 4.2.0)
UpSetR * 1.4.0 2019-05-22 [1] CRAN (R 4.2.0)
urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.2.0)
usethis * 2.1.6 2022-05-25 [1] CRAN (R 4.2.0)
utf8 1.2.2 2021-07-24 [1] CRAN (R 4.2.0)
vctrs 0.4.1 2022-04-13 [1] CRAN (R 4.2.0)
vegan * 2.6-2 2022-04-17 [1] CRAN (R 4.2.0)
withr 2.5.0 2022-03-03 [1] CRAN (R 4.2.0)
xfun 0.32 2022-08-10 [1] CRAN (R 4.2.0)
xtable 1.8-4 2019-04-21 [1] CRAN (R 4.2.0)
XVector 0.36.0 2022-04-26 [1] Bioconductor
yaml 2.3.5 2022-02-21 [1] CRAN (R 4.2.0)
zCompositions 1.4.0-1 2022-03-26 [1] CRAN (R 4.2.0)
zlibbioc 1.42.0 2022-04-26 [1] Bioconductor
[1] /Library/Frameworks/R.framework/Versions/4.2/Resources/library
──────────────────────────────────────────────────────────────────────────────────