Active microorganisms in forest soils differ from the total community yet are shaped by the same environmental factors: the influence of pH and soil moisture
GitHub Repository: http://zaklab-soils.github.io/MI_Gradient_Romanowicz_2016/
In Review for: FEMS Microbiology Ecology
Submitted for Review: March 22, 2016
Set the Global Knitr Options
# Set the global options
knitr::opts_chunk$set(echo = TRUE, eval = TRUE, fig.width = 6, fig.height = 5, fig.align = 'center', warning = FALSE, message = FALSE)
# Set the seed
set.seed(183409)Load in R-Packages and Set the Working Environment
# Load packages
necessary_packages <- c("data.table","vegan","ggplot2","grid","gridExtra","devtools","plyr","pander","reshape2","scales")
# Install the necessary packages
packages <- lapply(necessary_packages, library, character.only = TRUE)Set Table Options
#Set pander (table) output functions
panderOptions('table.split.table', Inf)
panderOptions('table.style', 'rmarkdown')
# Do not show scientific notation
options(scipen = 999, digits = 7)Import MOTHUR data
MOTHUR groups.summary Files
# Call in the MOTHUR "groups.summary" file for Bacterial Alpha Diversity metrics
bact.alpha <- read.table("c.g.DNA.16S.combined.trim.unique.good.filter.unique.precluster.pick.pick.an.unique_list.0.03.subsample.0.03.groups.summary", header=TRUE)
# Convert to data.table format
bact.alpha2 <- data.table(bact.alpha)
# Add new columns for "Community" and "Site" to the end of the data table
bact.alpha2$Site <- c("A","A","A","B","B","B","C","C","C","D","D","D","A","A","A","B","B","B","C","C","C","D","D","D")
bact.alpha2$Community <- c("Active","Active","Active","Active","Active","Active","Active","Active","Active","Active","Active","Active","Total","Total","Total","Total","Total","Total","Total","Total","Total","Total","Total","Total")
# Separate Total and Active Communities
bact.total.alpha<-bact.alpha2[c(13:24),]
bact.active.alpha<-bact.alpha2[c(1:12),]
# Call in the MOTHUR "groups.summary" file for Fungal Alpha Diversity metrics
fungi.alpha <- read.table("c.g.DNA.28S.combined.trim.unique.good.filter.unique.precluster.pick.pick.an.unique_list.0.01.subsample.0.01.groups.summary", header=TRUE)
# Convert to data.table format
fungi.alpha2 <- data.table(fungi.alpha)
# Add new columns for "Community" and "Site" to the end of the data table
fungi.alpha2$Site <- c("A","A","A","B","B","B","C","C","C","D","D","D","A","A","A","B","B","B","C","C","C","D","D","D")
fungi.alpha2$Community <- c("Active","Active","Active","Active","Active","Active","Active","Active","Active","Active","Active","Active","Total","Total","Total","Total","Total","Total","Total","Total","Total","Total","Total","Total")
# Separate Total and Active Communities
fungi.total.alpha<-fungi.alpha2[c(13:24),]
fungi.active.alpha<-fungi.alpha2[c(1:12),]Import Environmental Metadata
# Call in "env.txt" file with meta-data
env<-read.table("env.txt",header=TRUE)
new.env<-env
# Add new columns for "Community" and "Site" to the end of the data table
new.env$Site<-c("A","A","A","B","B","B","C","C","C","D","D","D","A","A","A","B","B","B","C","C","C","D","D","D")
new.env$Community<-c("Total","Total","Total","Total","Total","Total","Total","Total","Total","Total","Total","Total","Active","Active","Active","Active","Active","Active","Active","Active","Active","Active","Active","Active")
new.env.total<-new.env[1:12,]
new.env.active<-new.env[13:24,]Environmental Metadata Statistics
# Remove Community column from env file
new.env2<-new.env[,-11]
# Remove rows 13:24 (redundant env between active and total)
new.env2<-new.env2[-c(13:24), ]
# Calculate mean of each column by Site
table1.mean<-ddply(new.env2, .(Site), colwise(mean))
# Calculate sd of each column by Site
table1.sd<-ddply(new.env2, .(Site), colwise(sd))Recreating Table 1: Environmental Metadata
### Transpose data frame
# Remember Sites as non-numeric values
table1.mean.sites<-table1.mean$Site
table1.sd.sites<-table1.sd$Site
# Transpose all but first column (Site)
table1.mean<-as.data.frame(t(table1.mean[,-1]))
colnames(table1.mean)<-table1.mean.sites
table1.mean<-format(round(table1.mean,1))
table1.sd<-as.data.frame(t(table1.sd[,-1]))
colnames(table1.sd)<-table1.sd.sites
table1.sd<-format(round(table1.sd,1))
# Combine mean and sd into single column with ± divider
table1.combine<-as.data.frame(do.call(cbind, lapply(1:ncol(table1.mean), function(i) paste0(table1.mean[ , i], " ±", table1.sd[ , i]))))| Site A | Site B | Site C | Site D | |
|---|---|---|---|---|
| Moisture | 24.2 ± 3.9 | 66.3 ± 1.1 | 69.9 ± 1.9 | 63.6 ±1.1 |
| pH | 5.1 ± 0.1 | 5.7 ± 0.3 | 5.5 ± 0.1 | 5.5 ±0.1 |
| Aryl | 15.9 ± 0.5 | 15.6 ± 0.4 | 14.5 ± 0.4 | 14.4 ±0.6 |
| Alkyl | 20.0 ± 1.2 | 20.8 ± 1.6 | 22.2 ± 0.6 | 22.3 ±1.8 |
| O/N-alkyl | 56.7 ± 0.9 | 56.8 ± 1.5 | 57.1 ± 1.0 | 57.1 ±2.8 |
| Carboxyl | 7.5 ± 0.3 | 6.7 ± 0.5 | 6.2 ± 0.3 | 6.3 ±0.5 |
| C | 430.3 ±24.5 | 392.3 ±13.7 | 410.4 ±15.8 | 403.6 ±8.6 |
| N | 13.6 ± 0.7 | 17.3 ± 1.4 | 16.6 ± 0.5 | 14.9 ±1.4 |
| C:N | 31.6 ± 1.8 | 22.8 ± 2.1 | 24.7 ± 1.0 | 27.3 ±3.1 |
Statistics on Environmental Differences Across Sites
#Create matrix of numeric columns
env.bind.mat<-as.matrix(cbind(new.env2[,1:9]))
# MANOVA
env.manova<-manova(env.bind.mat~Site, data=new.env2)
env.man.sum<-summary.aov(env.manova)
env.man.sum## Response Moisture :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 4101.2 1367.07 260.62 0.00000002579 ***
## Residuals 8 42.0 5.25
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response pH :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 0.64737 0.215789 8.4485 0.007327 **
## Residuals 8 0.20433 0.025542
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Aryl :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 5.4016 1.80053 8.6916 0.006737 **
## Residuals 8 1.6573 0.20716
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Alkyl :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 11.396 3.7985 1.9545 0.1995
## Residuals 8 15.548 1.9435
##
## Response O_alkyl :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 0.3788 0.12627 0.0428 0.9873
## Residuals 8 23.5862 2.94828
##
## Response Carboxyl :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 3.0566 1.0189 6.5104 0.01536 *
## Residuals 8 1.2520 0.1565
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response C :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 2282.5 760.83 2.7394 0.1131
## Residuals 8 2221.9 277.74
##
## Response N :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 24.5907 8.1969 6.7322 0.01401 *
## Residuals 8 9.7406 1.2176
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response CN :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 129.820 43.273 9.5535 0.005069 **
## Residuals 8 36.237 4.530
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Run individual ANOVA on those Environmental factors with significant Site differences to get at Post-hoc analysis
# Moisture ANOVA with Post-Hoc for Site Differences
env.moisture<-aov(Moisture~Site, data=new.env2)
TukeyHSD(env.moisture)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Moisture ~ Site, data = new.env2)
##
## $Site
## diff lwr upr p adj
## B-A 42.027165 36.038743 48.0155878 0.0000001
## C-A 45.720123 39.731701 51.7085457 0.0000000
## D-A 39.380817 33.392395 45.3692396 0.0000001
## C-B 3.692958 -2.295464 9.6813802 0.2727584
## D-B -2.646348 -8.634770 3.3420742 0.5247361
## D-C -6.339306 -12.327728 -0.3508837 0.0384333# pH ANOVA with Post-Hoc for Site Differences
env.ph<-aov(pH~Site, data=new.env2)
TukeyHSD(env.ph)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = pH ~ Site, data = new.env2)
##
## $Site
## diff lwr upr p adj
## B-A 0.63666667 0.218790064 1.0545433 0.0053643
## C-A 0.44333333 0.025456731 0.8612099 0.0380357
## D-A 0.42000000 0.002123397 0.8378766 0.0488685
## C-B -0.19333333 -0.611209936 0.2245433 0.4895945
## D-B -0.21666667 -0.634543269 0.2012099 0.4011996
## D-C -0.02333333 -0.441209936 0.3945433 0.9977972# Aryl ANOVA with Post-Hoc for Site Differences
env.aryl<-aov(Aryl~Site, data=new.env2)
TukeyHSD(env.aryl)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Aryl ~ Site, data = new.env2)
##
## $Site
## diff lwr upr p adj
## B-A -0.2222217 -1.412297 0.9678537 0.9298607
## C-A -1.3956699 -2.585745 -0.2055946 0.0231855
## D-A -1.4885274 -2.678603 -0.2984521 0.0165323
## C-B -1.1734482 -2.363524 0.0166271 0.0532509
## D-B -1.2663057 -2.456381 -0.0762304 0.0375062
## D-C -0.0928575 -1.282933 1.0972178 0.9940820# Carboxyl ANOVA with Post-Hoc for Site Differences
env.carboxyl<-aov(Carboxyl~Site, data=new.env2)
TukeyHSD(env.carboxyl)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Carboxyl ~ Site, data = new.env2)
##
## $Site
## diff lwr upr p adj
## B-A -0.7126191 -1.7469863 0.3217481 0.2010875
## C-A -1.2790177 -2.3133849 -0.2446504 0.0175779
## D-A -1.1768443 -2.2112115 -0.1424770 0.0270401
## C-B -0.5663986 -1.6007658 0.4679686 0.3593749
## D-B -0.4642252 -1.4985924 0.5701420 0.5129485
## D-C 0.1021734 -0.9321938 1.1365406 0.9882196# N ANOVA with Post-Hoc for Site Differences
env.n<-aov(N~Site, data=new.env2)
TukeyHSD(env.n)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = N ~ Site, data = new.env2)
##
## $Site
## diff lwr upr p adj
## B-A 3.6421667 0.75700160 6.5273317 0.0157310
## C-A 2.9818333 0.09666827 5.8669984 0.0430058
## D-A 1.2656667 -1.61949840 4.1508317 0.5302980
## C-B -0.6603333 -3.54549840 2.2248317 0.8812625
## D-B -2.3765000 -5.26166506 0.5086651 0.1108683
## D-C -1.7161667 -4.60133173 1.1689984 0.2981955# CN ANOVA with Post-Hoc for Site Differences
env.cn<-aov(CN~Site, data=new.env2)
TukeyHSD(env.cn)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = CN ~ Site, data = new.env2)
##
## $Site
## diff lwr upr p adj
## B-A -8.776146 -14.340997 -3.211296 0.0043500
## C-A -6.876124 -12.440974 -1.311273 0.0176452
## D-A -4.296414 -9.861264 1.268437 0.1396704
## C-B 1.900023 -3.664828 7.464873 0.7031021
## D-B 4.479733 -1.085117 10.044583 0.1205691
## D-C 2.579710 -2.985140 8.144561 0.4880686Alpha Diversity
Bacterial Chao1 Richness
# Bacterial alpha diversity for Site x Community Interactions - Chao1 Richness
chao.bact<-aov(chao~Site*Community, data=bact.alpha2)
summary.aov(chao.bact)## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 13742424388 4580808129 8.759 0.00114 **
## Community 1 2719585260 2719585260 5.200 0.03663 *
## Site:Community 3 5124152936 1708050979 3.266 0.04885 *
## Residuals 16 8367603967 522975248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Tukey's HSD Post-hoc analysis
TukeyHSD(chao.bact)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = chao ~ Site * Community, data = bact.alpha2)
##
## $Site
## diff lwr upr p adj
## B-A 60087.162 22312.48 97861.84 0.0016743
## C-A 27400.111 -10374.57 65174.79 0.2030673
## D-A 3850.716 -33923.96 41625.40 0.9910238
## C-B -32687.051 -70461.73 5087.63 0.1024371
## D-B -56236.446 -94011.13 -18461.76 0.0030225
## D-C -23549.395 -61324.08 14225.29 0.3164802
##
## $Community
## diff lwr upr p adj
## Total-Active 21290 1498.38 41081.62 0.0366284
##
## $`Site:Community`
## diff lwr upr p adj
## B:Active-A:Active 28626.4935 -36019.394 93272.38 0.7799784
## C:Active-A:Active 31397.8902 -33247.998 96043.78 0.6981843
## D:Active-A:Active 6224.5320 -58421.356 70870.42 0.9999658
## A:Total-A:Active 8745.4656 -55900.422 73391.35 0.9996688
## B:Total-A:Active 100293.2961 35647.408 164939.18 0.0012547
## C:Total-A:Active 32147.7979 -32498.090 96793.69 0.6749062
## D:Total-A:Active 10222.3664 -54423.521 74868.25 0.9990887
## C:Active-B:Active 2771.3966 -61874.491 67417.28 0.9999999
## D:Active-B:Active -22401.9616 -87047.849 42243.93 0.9206056
## A:Total-B:Active -19881.0280 -84526.916 44764.86 0.9555236
## B:Total-B:Active 71666.8025 7020.915 136312.69 0.0244538
## C:Total-B:Active 3521.3044 -61124.583 68167.19 0.9999993
## D:Total-B:Active -18404.1272 -83050.015 46241.76 0.9701644
## D:Active-C:Active -25173.3582 -89819.246 39472.53 0.8670168
## A:Total-C:Active -22652.4246 -87298.312 41993.46 0.9164168
## B:Total-C:Active 68895.4059 4249.518 133541.29 0.0325076
## C:Total-C:Active 749.9078 -63895.980 65395.80 1.0000000
## D:Total-C:Active -21175.5238 -85821.412 43470.36 0.9392220
## A:Total-D:Active 2520.9336 -62124.954 67166.82 0.9999999
## B:Total-D:Active 94068.7641 29422.876 158714.65 0.0023763
## C:Total-D:Active 25923.2660 -38722.622 90569.15 0.8498766
## D:Total-D:Active 3997.8344 -60648.053 68643.72 0.9999984
## B:Total-A:Total 91547.8305 26901.943 156193.72 0.0030845
## C:Total-A:Total 23402.3324 -41243.556 88048.22 0.9030890
## D:Total-A:Total 1476.9008 -63168.987 66122.79 1.0000000
## C:Total-B:Total -68145.4982 -132791.386 -3499.61 0.0350942
## D:Total-B:Total -90070.9297 -154716.818 -25425.04 0.0035956
## D:Total-C:Total -21925.4315 -86571.319 42720.46 0.9282119Recreating Table 2. Alpha Diversity Bacterial Chao1 Richness
| Site A | Site B | Site C | Site D | |
|---|---|---|---|---|
| Total | 44395 ± 21017 | 135943 ± 42170 | 67797 ± 20669 | 45872 ± 5332 |
| Active | 35649 ± 11111 | 64276 ± 34446 | 67047 ± 5591 | 41874 ± 12922 |
Fungal Chao1 Richness
# Fungal alpha diversity - Chao1 Richness
chao.fungi<-aov(chao~Site*Community, data=fungi.alpha2)
summary.aov(chao.fungi)## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 1077854702 359284901 5.403 0.00924 **
## Community 1 308018855 308018855 4.632 0.04700 *
## Site:Community 3 1323902561 441300854 6.636 0.00403 **
## Residuals 16 1064007626 66500477
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Tukey's HSD post-hoc analysis
TukeyHSD(chao.fungi)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = chao ~ Site * Community, data = fungi.alpha2)
##
## $Site
## diff lwr upr p adj
## B-A -1367.870 -14838.03 12102.285 0.9911241
## C-A -10333.345 -23803.50 3136.811 0.1668916
## D-A -16370.788 -29840.94 -2900.633 0.0148108
## C-B -8965.474 -22435.63 4504.682 0.2651359
## D-B -15002.918 -28473.07 -1532.762 0.0264850
## D-C -6037.444 -19507.60 7432.712 0.5865062
##
## $Community
## diff lwr upr p adj
## Total-Active 7164.948 107.4101 14222.49 0.0469956
##
## $`Site:Community`
## diff lwr upr p adj
## B:Active-A:Active -2732.046594 -25784.264514 20320.171 0.9998619
## C:Active-A:Active 5584.968422 -17467.249498 28637.186 0.9877208
## D:Active-A:Active -4590.783599 -27643.001520 18461.434 0.9961412
## A:Total-A:Active 20332.018569 -2720.199351 43384.236 0.1053141
## B:Total-A:Active 20328.324390 -2723.893530 43380.542 0.1054173
## C:Total-A:Active -5919.639133 -28971.857054 17132.579 0.9829561
## D:Total-A:Active -7818.774711 -30870.992632 15233.443 0.9281973
## C:Active-B:Active 8317.015016 -14735.202904 31369.233 0.9045434
## D:Active-B:Active -1858.737005 -24910.954926 21193.481 0.9999898
## A:Total-B:Active 23064.065163 11.847243 46116.283 0.0498337
## B:Total-B:Active 23060.370984 8.153064 46112.589 0.0498855
## C:Total-B:Active -3187.592539 -26239.810460 19864.625 0.9996177
## D:Total-B:Active -5086.728117 -28138.946038 17965.490 0.9928612
## D:Active-C:Active -10175.752022 -33227.969942 12876.466 0.7824973
## A:Total-C:Active 14747.050147 -8305.167774 37799.268 0.3925915
## B:Total-C:Active 14743.355968 -8308.861952 37795.574 0.3928773
## C:Total-C:Active -11504.607556 -34556.825476 11547.610 0.6713124
## D:Total-C:Active -13403.743134 -36455.961054 9648.475 0.5031589
## A:Total-D:Active 24922.802168 1870.584248 47975.020 0.0293541
## B:Total-D:Active 24919.107990 1866.890069 47971.326 0.0293854
## C:Total-D:Active -1328.855534 -24381.073454 21723.362 0.9999990
## D:Total-D:Active -3227.991112 -26280.209032 19824.227 0.9995849
## B:Total-A:Total -3.694179 -23055.912099 23048.524 1.0000000
## C:Total-A:Total -26251.657702 -49303.875623 -3199.440 0.0199854
## D:Total-A:Total -28150.793280 -51203.011201 -5098.575 0.0114811
## C:Total-B:Total -26247.963524 -49300.181444 -3195.746 0.0200069
## D:Total-B:Total -28147.099102 -51199.317022 -5094.881 0.0114935
## D:Total-C:Total -1899.135578 -24951.353498 21153.082 0.9999882Recreating Table 2. Alpha Diversity Fungal Chao1 Richness
| Site A | Site B | Site C | Site D | |
|---|---|---|---|---|
| Total | 34051 ± 17222 | 34047 ± 6896 | 7799 ± 3426 | 5900 ± 1627 |
| Active | 13719 ± 885 | 10987 ± 623 | 19304 ± 12705 | 9128 ± 3300 |
Beta Diversity - PERMANOVA for Site and Community Influences
BACTERIA
# Bacteria - Bray-Curtis Dissimilarity (OTU relative abundance)
bact.bc.dist <- vegdist(bacteria, method = "bray", binary = FALSE)
# PERMANOVA - Bacteria - Bray-Curtis
adonis.bc.bact<-adonis(bact.bc.dist~Site*Community, data=new.env)
adonis.bc.bact##
## Call:
## adonis(formula = bact.bc.dist ~ Site * Community, data = new.env)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Site 3 1.6198 0.53993 1.1715 0.14863 0.001 ***
## Community 1 0.5102 0.51021 1.1070 0.04682 0.013 *
## Site:Community 3 1.3937 0.46457 1.0080 0.12789 0.289
## Residuals 16 7.3741 0.46088 0.67666
## Total 23 10.8979 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1FUNGI
# Fungi - Bray-Curtis Dissimilarity (OTU relative abundance)
fungi.bc.dist <- vegdist(fungi, method = "bray", binary = FALSE)
# PERMANOVA - Fungi - Bray-Curtis
adonis.bc.fungi<-adonis(fungi.bc.dist~Site*Community, data=new.env)
adonis.bc.fungi##
## Call:
## adonis(formula = fungi.bc.dist ~ Site * Community, data = new.env)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Site 3 2.0060 0.66865 2.0041 0.22815 0.001 ***
## Community 1 0.4807 0.48075 1.4409 0.05468 0.009 **
## Site:Community 3 0.9673 0.32243 0.9664 0.11002 0.630
## Residuals 16 5.3382 0.33364 0.60715
## Total 23 8.7922 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Site Effects within Total and Active Communities
# Separate total and active bacteria
bact.total <- bacteria[1:12,]
bact.active <- bacteria[13:24,]
#Create Bray-Curtis distance matrix for both communities
bact.total.bc.dist <- vegdist(bact.total, method = "bray", binary = FALSE)
bact.active.bc.dist <- vegdist(bact.active, method = "bray", binary = FALSE)
# Separate total and active fungi
fungi.total <- fungi[1:12,]
fungi.active <- fungi[13:24,]
#Create Bray-Curtis distance matrix for both communities
fungi.total.bc.dist <- vegdist(fungi.total, method = "bray", binary = FALSE)
fungi.active.bc.dist <- vegdist(fungi.active, method = "bray", binary = FALSE)PERMANOVA - Site effects on TOTAL BACTERIA
adonis.bc.bact.total <- adonis(bact.total.bc.dist~Site, data = new.env.total)
adonis.bc.bact.total##
## Call:
## adonis(formula = bact.total.bc.dist ~ Site, data = new.env.total)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Site 3 1.5268 0.50892 1.1081 0.29356 0.001 ***
## Residuals 8 3.6740 0.45926 0.70644
## Total 11 5.2008 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PERMANOVA - Site effects on ACTIVE BACTERIA
adonis.bc.bact.active <- adonis(bact.active.bc.dist~Site, data = new.env.active)
adonis.bc.bact.active##
## Call:
## adonis(formula = bact.active.bc.dist ~ Site, data = new.env.active)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Site 3 1.4868 0.49558 1.0715 0.28664 0.002 **
## Residuals 8 3.7001 0.46251 0.71336
## Total 11 5.1868 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PERMANOVA - Site effects on TOTAL FUNGI
adonis.bc.fungi.total <- adonis(fungi.total.bc.dist~Site, data = new.env.total)
adonis.bc.fungi.total##
## Call:
## adonis(formula = fungi.total.bc.dist ~ Site, data = new.env.total)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Site 3 1.4989 0.49964 1.4873 0.35804 0.001 ***
## Residuals 8 2.6875 0.33594 0.64196
## Total 11 4.1864 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PERMANOVA - Site effects on ACTIVE FUNGI
adonis.bc.fungi.active <- adonis(fungi.active.bc.dist~Site, data = new.env.active)
adonis.bc.fungi.active##
## Call:
## adonis(formula = fungi.active.bc.dist ~ Site, data = new.env.active)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Site 3 1.4743 0.49144 1.4832 0.35741 0.002 **
## Residuals 8 2.6507 0.33134 0.64259
## Total 11 4.1250 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Visualizing Beta Diversity Through Unconstrained Ordination - NMDS
# NMDS Ordination for Bacteria
ord1<-metaMDS(bacteria, trymax=1000, k=3)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1514181
## Run 1 stress 0.1448613
## ... New best solution
## ... procrustes: rmse 0.09897082 max resid 0.3561536
## Run 2 stress 0.1514174
## Run 3 stress 0.1444905
## ... New best solution
## ... procrustes: rmse 0.05249811 max resid 0.121341
## Run 4 stress 0.1438198
## ... New best solution
## ... procrustes: rmse 0.072648 max resid 0.1706673
## Run 5 stress 0.1438159
## ... New best solution
## ... procrustes: rmse 0.00420372 max resid 0.01549874
## Run 6 stress 0.1447424
## Run 7 stress 0.1447323
## Run 8 stress 0.1544887
## Run 9 stress 0.143821
## ... procrustes: rmse 0.001006251 max resid 0.003029177
## *** Solution reachedord1##
## Call:
## metaMDS(comm = bacteria, k = 3, trymax = 1000)
##
## global Multidimensional Scaling using monoMDS
##
## Data: wisconsin(sqrt(bacteria))
## Distance: bray
##
## Dimensions: 3
## Stress: 0.1438159
## Stress type 1, weak ties
## Two convergent solutions found after 9 tries
## Scaling: centring, PC rotation
## Species: expanded scores based on 'wisconsin(sqrt(bacteria))'scores(ord1)## NMDS1 NMDS2 NMDS3
## Total_A_1 -0.26868394 -0.246538197 -0.110917893
## Total_A_2 -0.12469932 -0.078296151 0.004690917
## Total_A_3 -0.37894268 -0.177221379 0.055682327
## Total_B_1 0.41719338 -0.430407934 0.050003664
## Total_B_2 0.42458874 -0.052267609 -0.294630921
## Total_B_3 0.22444248 0.300684934 0.542212773
## Total_C_1 -0.17595019 -0.287994326 0.231450279
## Total_C_2 0.23514415 -0.160334824 0.257881598
## Total_C_3 -0.04816849 -0.073054754 0.125697454
## Total_D_1 0.09125526 -0.290097734 -0.192765927
## Total_D_2 0.23636646 -0.181728864 -0.031655117
## Total_D_3 0.10662068 -0.027849955 0.104405739
## Active_A_1 -0.17375875 -0.036728955 -0.093598896
## Active_A_2 -0.20847768 0.034668086 -0.063015639
## Active_A_3 -0.37605822 -0.002443412 0.010180899
## Active_B_1 0.03171583 0.195684123 -0.105137971
## Active_B_2 0.22729108 0.399515672 -0.360082693
## Active_B_3 -0.05276909 0.327641494 0.029150638
## Active_C_1 -0.17563736 0.055520538 -0.211245603
## Active_C_2 0.20884903 0.271470198 0.128206002
## Active_C_3 -0.30814482 0.255044189 0.015931263
## Active_D_1 0.07704519 0.043810282 -0.204896120
## Active_D_2 0.05032799 0.020621871 -0.034351053
## Active_D_3 -0.03954973 0.140302708 0.146804278stressplot(ord1)
# Data for plotting
bact.NMDS.data<-new.env
bact.NMDS.data$NMDS1<-ord1$points[ ,1]
bact.NMDS.data$NMDS2<-ord1$points[ ,2]
# R2 Superscript for plot
bact.ss<-0.859
bact.ss<-paste("R^2 == ", round(bact.ss,2))# NMDS Ordination for Fungi
ord2<-metaMDS(fungi, trymax=1000, k=3)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1306279
## Run 1 stress 0.1304346
## ... New best solution
## ... procrustes: rmse 0.121362 max resid 0.3652675
## Run 2 stress 0.1326594
## Run 3 stress 0.1305449
## ... procrustes: rmse 0.01328777 max resid 0.04678179
## Run 4 stress 0.1303758
## ... New best solution
## ... procrustes: rmse 0.08682962 max resid 0.2865087
## Run 5 stress 0.1320392
## Run 6 stress 0.1341844
## Run 7 stress 0.1320326
## Run 8 stress 0.1320177
## Run 9 stress 0.130441
## ... procrustes: rmse 0.08519884 max resid 0.2816304
## Run 10 stress 0.1313978
## Run 11 stress 0.1391071
## Run 12 stress 0.1304984
## ... procrustes: rmse 0.024495 max resid 0.06126887
## Run 13 stress 0.13044
## ... procrustes: rmse 0.08728216 max resid 0.2853151
## Run 14 stress 0.1320258
## Run 15 stress 0.1320142
## Run 16 stress 0.1293662
## ... New best solution
## ... procrustes: rmse 0.09194912 max resid 0.3051822
## Run 17 stress 0.1325901
## Run 18 stress 0.1293693
## ... procrustes: rmse 0.002940371 max resid 0.008278037
## *** Solution reachedord2##
## Call:
## metaMDS(comm = fungi, k = 3, trymax = 1000)
##
## global Multidimensional Scaling using monoMDS
##
## Data: wisconsin(sqrt(fungi))
## Distance: bray
##
## Dimensions: 3
## Stress: 0.1293662
## Stress type 1, weak ties
## Two convergent solutions found after 18 tries
## Scaling: centring, PC rotation
## Species: expanded scores based on 'wisconsin(sqrt(fungi))'scores(ord2)## NMDS1 NMDS2 NMDS3
## Total_A_1 -0.349235916 0.189011319 0.134381024
## Total_A_2 -0.288571424 0.073579519 0.006254044
## Total_A_3 -0.412166374 0.389216740 -0.159824733
## Total_B_1 0.203517749 0.394541779 0.387434244
## Total_B_2 0.524239551 0.122718190 0.052251865
## Total_B_3 0.356679439 0.353931359 -0.362827947
## Total_C_1 -0.051440714 0.192130253 0.088283003
## Total_C_2 0.058303303 0.005608434 0.035573923
## Total_C_3 0.012456414 0.104941548 0.177724165
## Total_D_1 0.027704457 -0.041883724 -0.078563679
## Total_D_2 0.016421225 -0.035246171 0.010370622
## Total_D_3 -0.003288533 0.011068391 -0.011954368
## Active_A_1 -0.268817101 -0.086624691 -0.109755430
## Active_A_2 -0.326687895 -0.200222737 -0.051292122
## Active_A_3 -0.210290046 0.019504011 -0.233241265
## Active_B_1 0.158224598 -0.221635406 -0.030687077
## Active_B_2 0.304794010 -0.167109619 -0.152576947
## Active_B_3 0.216188228 -0.046876049 -0.238498466
## Active_C_1 -0.113192520 -0.107409027 0.089452393
## Active_C_2 0.140478381 -0.103523681 0.148640328
## Active_C_3 0.011072457 -0.205169098 0.342184125
## Active_D_1 -0.013916713 -0.192291961 -0.114725060
## Active_D_2 -0.010612084 -0.294668523 0.054670981
## Active_D_3 0.018139507 -0.153590855 0.016726378stressplot(ord2)
# Data for plotting
fungi.NMDS.data<-new.env
fungi.NMDS.data$NMDS1<-ord2$points[ ,1]
fungi.NMDS.data$NMDS2<-ord2$points[ ,2]
# R2 Superscript
fung.ss<-0.904
fung.ss<-paste("R^2 == ", round(fung.ss,2))Reproducing Figure 1 - NMDS
# Set geom_point shape to new default prior to plotting
update_geom_defaults("point", list(shape=21))
# Plotting bacteria NMDS in ggplot2
bact.ggplot<-ggplot(bact.NMDS.data, aes(x=NMDS1, y=NMDS2)) + aes(shape=factor(Site)) + geom_point(colour="black", size=3.5, aes(fill=factor(Community))) + scale_fill_manual(values=c("white", "black")) + scale_shape_manual("Site", values=c(21,22,23,24))+ scale_colour_manual(values=c("black", "black")) + labs(fill="Community") + theme_bw() + scale_x_continuous(limits=c(-0.5, 0.5)) + scale_y_continuous(limits=c(-0.55, 0.45)) + annotate("text", x=-0.48, y=0.45, label = "(a)", size=6) + theme(legend.key=element_blank()) + annotate("text", x=-0.5, y=-0.50, label=bact.ss, parse=TRUE, hjust=0, size=4) + annotate("text", x=-0.5, y=-0.55, label="Stress=0.14", hjust=0, size=4) + theme(axis.ticks.y = element_blank(), axis.ticks.x = element_blank(), panel.background = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank())
# Plotting fungi NMDS in ggplot2
fungi.ggplot<-ggplot(fungi.NMDS.data, aes(x=NMDS1, y=NMDS2)) + aes(shape=factor(Site)) + geom_point(colour="black", size=3.5, aes(fill=factor(Community))) + scale_fill_manual(values=c("white", "black")) + scale_shape_manual("Site", values=c(21,22,23,24)) + scale_colour_manual(values=c("black", "black")) + labs(fill="Community") + theme_bw() + scale_x_continuous(limits=c(-0.6, 0.6)) + scale_y_continuous(limits=c(-0.55, 0.45)) + annotate("text", x=-0.58, y=0.45, label = "(b)", size=6) + theme(legend.key=element_blank()) + annotate("text", x=-0.6, y=-0.5, label=fung.ss, parse=TRUE, hjust=0, size=4) + annotate("text", x=-0.6, y=-0.55, label="Stress=0.13", hjust=0, size=4) + theme(axis.ticks.y = element_blank(), axis.ticks.x = element_blank(), panel.background = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank())# Arrange ggplot2 graphs with a specific width
grid.arrange(bact.ggplot2, fungi.ggplot2, legend, ncol=3, widths=c(5,5,2))
Figure 1. Nonmetric multidimensional scaling (NMDS) ordinations based on Bray-Curtis dissimilarity of OTU abundance between (a) bacterial and (b) fungal communities reveals spatial variation by site and between total and active communities. Each marker represents an individual plot in the NMDS charts.
Taxonomy - Bacterial Community Relative Abundance
# Call in bacteria taxonomy table (absolute abundance)
bact.tax<-read.table("c.g.DNA.16S.combined.trim.unique.good.filter.unique.precluster.pick.pick.an.unique_list.0.03.subsample.0.03.cons.tax.summary", header = TRUE, sep = "\t")
# Subset taxonomy file to pull out only order level taxa ("taxlevel" = 4)
bact.order<-subset(bact.tax, taxlevel==4)
# Remove N-deposition plots
bact.order<-bact.order[, -c(4:5,9:11,15:17,21:23,27:29,33:35,39:41,45:47,51:54)]
# Create new column with rankID and taxon
bact.order$taxID<-paste(bact.order$taxon,bact.order$rankID,sep="_")
# Preserve taxID as row names
row.names(bact.order)<-bact.order$taxID
# Remove non-numeric columns
bact.order.new<-bact.order[,-c(1:3,28)]
# Convert from absolute to relative abundance
bact.order.sums<-colSums(bact.order.new)
bact.order.test1<-as.matrix(bact.order.new)
bact.order.test2<-1/bact.order.sums
bact.order.test3<-t(bact.order.test1)
bact.order.test4<-(bact.order.test3*bact.order.test2)*100
bact.order.test5<-t(bact.order.test4)
bact.order.rel<-as.data.frame(bact.order.test5)
bact.order.rel$taxID<-row.names(bact.order.rel)
# Add factor columns back in
bact.order.factor<-bact.order[,c(1:3,28)]
# Merge factors back into Rel Abund table
bact.order.final<-merge(bact.order.factor,bact.order.rel,by="taxID",all=TRUE)
# Remove certain factors for downstream processing
bact.order.final<-bact.order.final[,-c(2:3)]Separate Active Bacteria
# Active data frame
bact.order.active<-bact.order.final[,1:14]
# Calculate taxon means and add new means column
bact.order.active$MEANS_Active<-rowMeans(bact.order.active[,3:14])
# Calculate taxon sd and add new sd column
bact.order.active$SD_Active<-apply(as.matrix(bact.order.active[,3:14]),1,sd)
# Preserve taxID as row names
row.names(bact.order.active)<-bact.order.active$taxID
# Remove non-numeric columns
bact.order.active.new<-bact.order.active[,-1]
# Collect all taxon greater than 1% rel abund
bact.order.active.relevant<-subset(bact.order.active.new, MEANS_Active>=1 | taxon=="Flavobacteriales")
bact.order.active.relevant<-subset(bact.order.active.relevant, taxon!="unclassified")
#Remove taxon column
bact.order.active.relevant<-bact.order.active.relevant[,-1]
# Collect all taxon less than 1% rel abund
bact.order.active.other<-subset(bact.order.active.new, MEANS_Active<1 | taxon=="unclassified")
bact.order.active.other<-subset(bact.order.active.other, taxon!="Flavobacteriales")
# Sum "unclassfied" into single row
bact.order.active.other<-ddply(bact.order.active.other, "taxon", numcolwise(sum))
# Pull "unclassified" from "other" table
bact.order.active.other.noclass<-subset(bact.order.active.other, taxon!="unclassified")
# Create "unclassified" factor to pull out these taxa from "other" table
bact.active.unclass<-"unclassified"
bact.order.active.other.unclass<-bact.order.active.other[bact.order.active.other$taxon %in% bact.active.unclass,]
# Preserve taxon as row names
row.names(bact.order.active.other.noclass)<-bact.order.active.other.noclass$taxon
row.names(bact.order.active.other.unclass)<-bact.order.active.other.unclass$taxon
# Save names of all taxa placed in "Other" category
bact.other.names<-row.names(bact.order.active.other.noclass)
#Remove taxon column
bact.order.active.other.noclass<-bact.order.active.other.noclass[,-1]
bact.order.active.other.unclass<-bact.order.active.other.unclass[,-1]
# Column sums
bact.order.active.other.vec<-colSums(bact.order.active.other.noclass)
# Bind "Other" to relevant data frame
bact.order.active.final<-rbind(bact.order.active.relevant,bact.order.active.other.vec,bact.order.active.other.unclass)
# Rename row names with "Other" column
row.names(bact.order.active.final)<-c("Acidobacteria_Gp1","Actinomycetales","Burkholderiales","Caulobacterales","Flavobacteriales","Planctomycetales","Pseudomonadales","Rhizobiales","Rhodospirillales","Sphingobacteriales","Sphingomonadales","Other","Unclassified")
# Add back in Taxon column
bact.order.active.final$Taxon<-row.names(bact.order.active.final)Separate Total Bacteria
# Total data frame
bact.order.total<-bact.order.final[,c(1:2,15:26)]
# Calculate taxon means and add new means column
bact.order.total$MEANS_Total<-rowMeans(bact.order.total[,3:14])
# Calculate taxon sd and add new sd column
bact.order.total$SD_Total<-apply(as.matrix(bact.order.total[,3:14]),1,sd)
# Preserve taxID as row names
row.names(bact.order.total)<-bact.order.total$taxID
# Remove non-numeric columns
bact.order.total.new<-bact.order.total[,-1]
# Collect all taxon greater than 1% rel abund
bact.order.total.relevant<-subset(bact.order.total.new, MEANS_Total>=1)
# Remove unclassified taxa
bact.order.total.relevant<-subset(bact.order.total.relevant, taxon!="unclassified")
# Remove TM7 taxa
bact.order.total.relevant<-subset(bact.order.total.relevant, taxon!="TM7_order_incertae_sedis")
# Remove taxon column
bact.order.total.relevant<-bact.order.total.relevant[,-1]
# Collect all taxon less than 1% rel abund including all unclassified and TM7
bact.order.total.other<-subset(bact.order.total.new, MEANS_Total<1 | taxon=="unclassified" | taxon=="TM7_order_incertae_sedis")
# Sum "unclassfied" into single row
bact.order.total.other<-ddply(bact.order.total.other, "taxon", numcolwise(sum))
# Pull "unclassified" from "other" table
bact.order.total.other.noclass<-subset(bact.order.total.other, taxon!="unclassified")
# Create "unclassified" object to pull these taxa out of "other" table
bact.total.unclass<-"unclassified"
bact.order.total.other.unclass<-bact.order.total.other[bact.order.total.other$taxon %in% bact.total.unclass,]
# Preserve taxon as row names
row.names(bact.order.total.other.noclass)<-bact.order.total.other.noclass$taxon
row.names(bact.order.total.other.unclass)<-bact.order.total.other.unclass$taxon
#Remove taxon column
bact.order.total.other.noclass<-bact.order.total.other.noclass[,-1]
bact.order.total.other.unclass<-bact.order.total.other.unclass[,-1]
# Column sums
bact.order.total.other.vec<-colSums(bact.order.total.other.noclass)
# Bind "Other" to relevant data frame
bact.order.total.final<-rbind(bact.order.total.relevant,bact.order.total.other.vec,bact.order.total.other.unclass)
# Rename row names with "Other" column
row.names(bact.order.total.final)<-c("Acidobacteria_Gp1","Actinomycetales","Burkholderiales","Caulobacterales","Flavobacteriales","Planctomycetales","Pseudomonadales","Rhizobiales","Rhodospirillales","Sphingobacteriales","Sphingomonadales","Other","Unclassified")
# Add back in Taxon column
bact.order.total.final$Taxon<-row.names(bact.order.total.final)
# Merge Active and Total data frames
bact.order.merge<-merge(bact.order.total.final,bact.order.active.final,by="Taxon",all=TRUE)MANOVA on change in relative abundance between bacterial communities at the order-level
# Final Bacteria Taxonomy table for MANOVA
bact.tax.final<-bact.order.merge[, -c(14,15,28,29)]
### Transpose taxonomy table
# Remember taxon as non-numeric values
bact.tax.final.taxon<-bact.tax.final$Taxon
# Transpose all but first column (Taxon)
bact.tax.final<-as.data.frame(t(bact.tax.final[,-1]))
colnames(bact.tax.final)<-bact.tax.final.taxon
# Add new columns for "Site" and "Community" to the end of the data table
bact.tax.final$Site <- c("A","A","A","B","B","B","C","C","C","D","D","D","A","A","A","B","B","B","C","C","C","D","D","D")
bact.tax.final$Community <- c("Total","Total","Total","Total","Total","Total","Total","Total","Total","Total","Total","Total","Active","Active","Active","Active","Active","Active","Active","Active","Active","Active","Active","Active")
#Rename "Acidobacteria Gp1" without spaces
names(bact.tax.final)[names(bact.tax.final)=='Acidobacteria Gp1']<-'Acidobacteria_Gp1'
#Create matrix of numeric columns
bact.bind.mat<-as.matrix(cbind(bact.tax.final[,1:13]))
# MANOVA
bact.manova<-manova(bact.bind.mat~Community, data=bact.tax.final)
bact.man.sum<-summary.aov(bact.manova)
bact.man.sum## Response Acidobacteria_Gp1 :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 58.45 58.448 0.9802 0.3329
## Residuals 22 1311.79 59.627
##
## Response Actinomycetales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 12.57 12.568 0.6624 0.4244
## Residuals 22 417.43 18.974
##
## Response Burkholderiales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 16.655 16.6551 3.26 0.08469 .
## Residuals 22 112.395 5.1089
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Caulobacterales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 17.301 17.3009 4.7434 0.04041 *
## Residuals 22 80.242 3.6474
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Flavobacteriales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 0.0572 0.05716 0.0582 0.8116
## Residuals 22 21.6129 0.98241
##
## Response Other :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 0.625 0.62503 0.2161 0.6466
## Residuals 22 63.623 2.89193
##
## Response Planctomycetales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 23.836 23.8360 4.4975 0.04545 *
## Residuals 22 116.596 5.2998
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Pseudomonadales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 37.569 37.569 4.9176 0.03722 *
## Residuals 22 168.072 7.640
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Rhizobiales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 41.362 41.362 5.5291 0.02807 *
## Residuals 22 164.578 7.481
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Rhodospirillales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 16.454 16.4542 3.0497 0.0947 .
## Residuals 22 118.699 5.3954
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Sphingobacteriales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 18.858 18.8584 5.3614 0.0303 *
## Residuals 22 77.383 3.5174
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Sphingomonadales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 0.102 0.10179 0.066 0.7997
## Residuals 22 33.956 1.54345
##
## Response Unclassified :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 3.775 3.7755 0.4254 0.521
## Residuals 22 195.255 8.8752# Extract p-values
bact.man.p<-lapply(bact.man.sum, function(i) {
i[1,5]
})
# Make into df
bact.man.p.df<-as.data.frame(do.call(rbind,bact.man.p))
# Clean up the format of df
colnames(bact.man.p.df)<-"P-Value"
# Get rid of "Response"
bact.man.p.df$ID<-row.names(bact.man.p.df)
bact.man.p.df$Taxon<-gsub(" Response ","",bact.man.p.df$ID)
# Make Taxon row names
row.names(bact.man.p.df)<-bact.man.p.df$Taxon
# Remove ID and Taxon columns
bact.man.p.df<-subset(bact.man.p.df,select="P-Value")Recreating Table 3. Relative abundance of total and active bacteria with % difference between communities.
# Final Bacteria Means table (Table 3 recreated)
bact.tax.final.means<-subset(bact.order.merge, select=c("Taxon","MEANS_Total","SD_Total","MEANS_Active","SD_Active"))
# Sort in decending order from Total community
bact.tax.final.means<-bact.tax.final.means[order(bact.tax.final.means$MEANS_Total,decreasing = TRUE),]
# Round numbers in each column
bact.tax.final.means$MEANS_Total<-round(bact.tax.final.means$MEANS_Total,0)
bact.tax.final.means$MEANS_Active<-round(bact.tax.final.means$MEANS_Active,0)
bact.tax.final.means$SD_Total<-round(bact.tax.final.means$SD_Total,1)
bact.tax.final.means$SD_Active<-round(bact.tax.final.means$SD_Active,1)
# Calculate difference between Total and Active Means
bact.tax.final.means$Difference<-bact.tax.final.means$MEANS_Active - bact.tax.final.means$MEANS_Total
bact.tax.final.means$Difference<-round(bact.tax.final.means$Difference,0)
# Preserve taxon as row names
row.names(bact.tax.final.means)<-bact.tax.final.means$Taxon
bact.tax.final.means<-bact.tax.final.means[,-1]
# Paste Total means with sd into single column
bact.tax.final.means$Total<-paste(bact.tax.final.means$MEANS_Total,"±",bact.tax.final.means$SD_Total,sep=" ")
# Paste Active means with sd into single column
bact.tax.final.means$Active<-paste(bact.tax.final.means$MEANS_Active,"±",bact.tax.final.means$SD_Active,sep=" ")
# Subset out columns of interest
bact.table2<-subset(bact.tax.final.means,select=c("Total","Active","Difference"))
# Merging everything into Table 3
bact.table2$Taxon<-row.names(bact.table2)
bact.man.p.df$Taxon<-row.names(bact.man.p.df)
#Final table
bact.table2.final<-merge(bact.table2,bact.man.p.df,by="Taxon",all=TRUE)
# Round p-values to 2 decimals
bact.table2.final$`P-Value`<-round(bact.table2.final$`P-Value`,3)
# Put P-values in parantheses
bact.table2.final$P2<-paste("(",bact.table2.final$`P-Value`,")",sep="")
# Add + signs to Difference factors
for(i in 1:length(bact.table2.final$Taxon)) {
if (bact.table2.final[i,4]>0) {
bact.table2.final[i,4]<-paste("+",bact.table2.final[i,4],sep="")}
}
# Combine Difference with P-values in single column
bact.table2.final$`Difference % (P-Value)`<-paste(bact.table2.final$Difference,bact.table2.final$P2,sep=" ")
# Remove unnecessary columns from final table
bact.table2.final<-bact.table2.final[,-c(4,5,6)]
bact.table2.final<-bact.table2.final[c("Taxon","Total","Active","Difference % (P-Value)")]
# Reorder Taxon by Total Rel Abund
bact.table2.final$Taxon<-as.factor(bact.table2.final$Taxon)
bact.table2.final$Taxon<-reorder(bact.table2.final$Taxon,c(3,1,5,6,10,13,7,11,2,8,4,9,12))
bact.table2.final<-bact.table2.final[order(bact.table2.final$Taxon),]
# Name Taxon row name
row.names(bact.table2.final)<-bact.table2.final$Taxon
# Remove Taxon column
bact.table2.final<-bact.table2.final[,-1]| % Total | % Active | % Difference (P-value) | |
|---|---|---|---|
| Actinomycetales | 22 ± 5.5 | 20 ± 2.8 | -2 (0.424) |
| Rhizobiales | 14 ± 3.6 | 11 ± 1.5 | -3 (0.028) |
| Acidobacteria Gp1 | 14 ± 8.8 | 11 ± 6.5 | -3 (0.333) |
| Sphingobacteriales | 8 ± 1.2 | 10 ± 2.4 | +2 (0.03) |
| Burkholderiales | 5 ± 1.6 | 7 ± 2.7 | +2 (0.085) |
| Caulobacterales | 5 ± 2.4 | 3 ± 1.3 | -2 (0.04) |
| Planctomycetales | 4 ± 2.5 | 6 ± 2.1 | +2 (0.045) |
| Rhodospirillales | 4 ± 2.1 | 6 ± 2.5 | +2 (0.095) |
| Sphingomonadales | 4 ± 1.2 | 4 ± 1.3 | 0 (0.8) |
| Flavobacteriales | 1 ± 0.9 | 1 ± 1.1 | 0 (0.812) |
| Pseudomonadales | 1 ± 1 | 4 ± 3.8 | +3 (0.037) |
| Unclassified | 13 ± 5.1 | 12 ± 4.4 | -1 (0.521) |
| Other | 5 ± 4 | 6 ± 4 | +1 (0.647) |
Supplemental Table S3. List of “Other” bacteria (<1% relative abundance)
bact.other.names## [1] "Acidimicrobiales"
## [2] "Acidobacteria_Gp10_order_incertae_sedis"
## [3] "Acidobacteria_Gp16_order_incertae_sedis"
## [4] "Acidobacteria_Gp2_order_incertae_sedis"
## [5] "Acidobacteria_Gp3_order_incertae_sedis"
## [6] "Acidobacteria_Gp3_order_incetae_sedis"
## [7] "Acidobacteria_Gp4_order_incertae_sedis"
## [8] "Acidobacteria_Gp5_order_incertae_sedis"
## [9] "Acidobacteria_Gp6_order_incertae_sedis"
## [10] "Acidobacteria_Gp7_order_incertae_sedis"
## [11] "Alphaproteobacteria_order_incertae_sedis"
## [12] "Alteromonadales"
## [13] "Anaerolineales"
## [14] "Armatimonadales"
## [15] "Armatimonadetes_gp4_order_incetae_sedis"
## [16] "Armatimonadetes_gp5_order_incetae_sedis"
## [17] "Bacillales"
## [18] "Bacteroidales"
## [19] "Bacteroidetes_incertae_sedis_order_incertae_sedis"
## [20] "Bdellovibrionales"
## [21] "Caldilineales"
## [22] "Chloroflexales"
## [23] "Chthonomonadales"
## [24] "Clostridiales"
## [25] "Coriobacteriales"
## [26] "Deinococcales"
## [27] "Enterobacteriales"
## [28] "Entomoplasmatales"
## [29] "Erysipelotrichales"
## [30] "Fusobacteriales"
## [31] "Gammaproteobacteria_order_incertae_sedis"
## [32] "Gemmatimonadales"
## [33] "Ktedonobacterales"
## [34] "Lactobacillales"
## [35] "Legionellales"
## [36] "Methylophilales"
## [37] "Myxococcales"
## [38] "Neisseriales"
## [39] "Nitrosomonadales"
## [40] "Oceanospirillales"
## [41] "OD1_order_incertae_sedis"
## [42] "Opitutales"
## [43] "Pasteurellales"
## [44] "Phycisphaerales"
## [45] "Rhodobacterales"
## [46] "Rhodocyclales"
## [47] "Rickettsiales"
## [48] "Selenomonadales"
## [49] "Solirubrobacterales"
## [50] "Spartobacteria_order_incertae_sedis"
## [51] "Spirochaetales"
## [52] "SR1_order_incertae_sedis"
## [53] "Subdivision3_order_incertae_sedis"
## [54] "Thiotrichales"
## [55] "TM7_order_incertae_sedis"
## [56] "Verrucomicrobiales"
## [57] "Xanthomonadales"Taxonomy - Fungal Community Relative Abundance
# Call in fungal taxonomy table (absolute abundance)
fungi.tax<-read.table("c.g.DNA.28S.combined.trim.unique.good.filter.unique.precluster.pick.pick.an.unique_list.0.01.subsample.0.01.cons.tax.summary", header = TRUE, sep = "\t")
# Subset taxonomy file to pull out only class level taxa ("taxlevel" = 4)
fungi.phyla<-subset(fungi.tax, taxlevel==4)
# Remove N-deposition plot columns
fungi.phyla<-fungi.phyla[, -c(4:5,9:11,15:17,21:23,27:29,33:35,39:41,45:47,51:54)]
# Create new column with rankID and taxon
fungi.phyla$taxID<-paste(fungi.phyla$taxon,fungi.phyla$rankID,sep="_")
# Preserve taxID as row names
row.names(fungi.phyla)<-fungi.phyla$taxID
# Remove non-numeric columns
fungi.phyla.new<-fungi.phyla[,-c(1:3,28)]
# Convert from absolute to relative abundance
fungi.sums<-colSums(fungi.phyla.new)
fungi.test1<-as.matrix(fungi.phyla.new)
fungi.test2<-1/fungi.sums
fungi.test3<-t(fungi.test1)
fungi.test4<-(fungi.test3*fungi.test2)*100
fungi.test5<-t(fungi.test4)
fungi.phyla.rel<-as.data.frame(fungi.test5)
fungi.phyla.rel$taxID<-row.names(fungi.phyla.rel)
# Add factor columns back in
fungi.phyla.factor<-fungi.phyla[,c(1:3,28)]
# Merge factors back into Rel Abund table
fungi.phyla.final<-merge(fungi.phyla.factor,fungi.phyla.rel,by="taxID",all=TRUE)
# Remove certain factors for downstream processing
fungi.phyla.final<-fungi.phyla.final[,-c(2:3)]Separate Active Fungi
# Active data frame
fungi.active<-fungi.phyla.final[,1:14]
# Calculate taxon means and add new means column
fungi.active$MEANS_Active<-rowMeans(fungi.active[,3:14])
# Calculate taxon sd and add new sd column
fungi.active$SD_Active<-apply(as.matrix(fungi.active[,3:14]),1,sd)
# Preserve taxID as row names
row.names(fungi.active)<-fungi.active$taxID
# Remove non-numeric columns
fungi.active.new<-fungi.active[,-1]
# Collect all taxon greater than 1% rel abund
fungi.active.relevant<-subset(fungi.active.new, MEANS_Active>=1)
fungi.active.relevant<-subset(fungi.active.relevant, taxon!="unclassified")
#Remove taxon column
fungi.active.relevant<-fungi.active.relevant[,-1]
# Collect all taxon less than 1% rel abund
fungi.active.other<-subset(fungi.active.new, MEANS_Active<1 | taxon=="unclassified")
# Sum "unclassfied" into single row
fungi.active.other<-ddply(fungi.active.other, "taxon", numcolwise(sum))
# Pull "unclassified" from "other" table
fungi.active.other.noclass<-subset(fungi.active.other, taxon!="unclassified")
# Create "unclassified" factor to pull out these taxa from "other" table
fungi.active.unclass<-"unclassified"
fungi.active.other.unclass<-fungi.active.other[fungi.active.other$taxon %in% fungi.active.unclass,]
# Preserve taxon as row names
row.names(fungi.active.other.noclass)<-fungi.active.other.noclass$taxon
row.names(fungi.active.other.unclass)<-fungi.active.other.unclass$taxon
# Save names of all taxa placed in "Other" category
fungi.other.names<-row.names(fungi.active.other.noclass)
#Remove taxon column
fungi.active.other.noclass<-fungi.active.other.noclass[,-1]
fungi.active.other.unclass<-fungi.active.other.unclass[,-1]
# Column sums
fungi.active.other.vec<-colSums(fungi.active.other.noclass)
# Bind "Other" and "unclassified" to relevant data frame
fungi.active.final<-rbind(fungi.active.relevant,fungi.active.other.vec,fungi.active.other.unclass)
# Rename row names with "Other" column
row.names(fungi.active.final)<-c("Agaricomycetes","Dothideomycetes","Eurotiomycetes","Leotiomycetes","Microbotryomycetes","Sordariomycetes","Tremellomycetes","Other","Unclassified")
# Add back in Taxon column
fungi.active.final$Taxon<-row.names(fungi.active.final)Separate Total Fungi
# Total data frame
fungi.total<-fungi.phyla.final[,c(1:2,15:26)]
# Calculate taxon means and add new means column
fungi.total$MEANS_Total<-rowMeans(fungi.total[,3:14])
# Calculate taxon sd and add new sd column
fungi.total$SD_Total<-apply(as.matrix(fungi.total[,3:14]),1,sd)
# Preserve taxID as row names
row.names(fungi.total)<-fungi.total$taxID
# Remove non-numeric columns
fungi.total.new<-fungi.total[,-1]
# Collect all taxon greater than 1% rel abund
fungi.total.relevant<-subset(fungi.total.new, MEANS_Total>=1)
fungi.total.relevant<-subset(fungi.total.relevant, taxon!="unclassified")
#Remove taxon column
fungi.total.relevant<-fungi.total.relevant[,-1]
# Collect all taxon less than 1% rel abund
fungi.total.other<-subset(fungi.total.new, MEANS_Total<1 | taxon=="unclassified")
# Sum "unclassfied" into single row
fungi.total.other<-ddply(fungi.total.other, "taxon", numcolwise(sum))
# Pull "unclassified" from "other" table
fungi.total.other.noclass<-subset(fungi.total.other, taxon!="unclassified")
# Create "unclassified" factor to pull out these taxa from "other" table
fungi.total.unclass<-"unclassified"
fungi.total.other.unclass<-fungi.total.other[fungi.total.other$taxon %in% fungi.total.unclass,]
# Preserve taxon as row names
row.names(fungi.total.other.noclass)<-fungi.total.other.noclass$taxon
row.names(fungi.total.other.unclass)<-fungi.total.other.unclass$taxon
#Remove taxon column
fungi.total.other.noclass<-fungi.total.other.noclass[,-1]
fungi.total.other.unclass<-fungi.total.other.unclass[,-1]
# Column sums
fungi.total.other.vec<-colSums(fungi.total.other.noclass)
# Bind "Other" and "unclassified" to relevant data frame
fungi.total.final<-rbind(fungi.total.relevant,fungi.total.other.vec,fungi.total.other.unclass)
# Rename row names with "Other" column
row.names(fungi.total.final)<-c("Agaricomycetes","Dothideomycetes","Eurotiomycetes","Leotiomycetes","Microbotryomycetes","Sordariomycetes","Tremellomycetes","Other","Unclassified")
# Add back in Taxon column
fungi.total.final$Taxon<-row.names(fungi.total.final)
# Merge Active and Total data frames
fungi.merge<-merge(fungi.total.final,fungi.active.final,by="Taxon",all=TRUE)MANOVA on change in relative abundance between fungal communities at the class-level
# Final Fungal taxonomy table for MANOVA
fungi.tax.final<-fungi.merge[, -c(14,15,28,29)]
### Transpose taxonomy table
# Remember taxon as non-numeric values
fungi.tax.final.taxon<-fungi.tax.final$Taxon
# Transpose all but first column (Taxon)
fungi.tax.final<-as.data.frame(t(fungi.tax.final[,-1]))
colnames(fungi.tax.final)<-fungi.tax.final.taxon
# Add new columns for "Site" and "Community" to the end of the data table
fungi.tax.final$Site <- c("A","A","A","B","B","B","C","C","C","D","D","D","A","A","A","B","B","B","C","C","C","D","D","D")
fungi.tax.final$Community <- c("Total","Total","Total","Total","Total","Total","Total","Total","Total","Total","Total","Total","Active","Active","Active","Active","Active","Active","Active","Active","Active","Active","Active","Active")
#Create matrix of numeric columns
fungi.bind.mat<-as.matrix(cbind(fungi.tax.final[,1:9]))
# MANOVA
fungi.manova<-manova(fungi.bind.mat~Community, data=fungi.tax.final)
fungi.man.sum<-summary.aov(fungi.manova)
fungi.man.sum## Response Agaricomycetes :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 295.51 295.5 11.068 0.003061 **
## Residuals 22 587.41 26.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Dothideomycetes :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 81.935 81.935 8.7148 0.007367 **
## Residuals 22 206.841 9.402
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Eurotiomycetes :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 0.5667 0.56673 0.5564 0.4636
## Residuals 22 22.4099 1.01863
##
## Response Leotiomycetes :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 1134.8 1134.80 6.3918 0.01915 *
## Residuals 22 3905.8 177.54
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Microbotryomycetes :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 62.323 62.323 17.527 0.0003826 ***
## Residuals 22 78.229 3.556
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Other :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 3.8686 3.8686 5.6927 0.02607 *
## Residuals 22 14.9507 0.6796
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Sordariomycetes :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 142.65 142.646 1.6572 0.2114
## Residuals 22 1893.66 86.075
##
## Response Tremellomycetes :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 65.470 65.470 23.814 0.00007048 ***
## Residuals 22 60.482 2.749
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Unclassified :
## Df Sum Sq Mean Sq F value Pr(>F)
## Community 1 0.47 0.4659 0.03 0.864
## Residuals 22 341.16 15.5072# Extract p-values
fungi.man.p<-lapply(fungi.man.sum, function(i) {
i[1,5]
})
# Make into df
fungi.man.p.df<-as.data.frame(do.call(rbind,fungi.man.p))
# Clean up the format of df
colnames(fungi.man.p.df)<-"P-Value"
# Get rid of "Response"
fungi.man.p.df$ID<-row.names(fungi.man.p.df)
fungi.man.p.df$Taxon<-gsub(" Response ","",fungi.man.p.df$ID)
# Make Taxon row names
row.names(fungi.man.p.df)<-fungi.man.p.df$Taxon
# Remove ID and Taxon columns
fungi.man.p.df<-subset(fungi.man.p.df,select="P-Value")Recreating Table 3. Relative abundance of total and active fungi with % difference between communities.
# Final Fungi Means table (Table 3 recreated)
fungi.tax.final.means<-subset(fungi.merge, select=c("Taxon","MEANS_Total","SD_Total","MEANS_Active","SD_Active"))
# Sort in decending order from Total community
fungi.tax.final.means<-fungi.tax.final.means[order(fungi.tax.final.means$MEANS_Total,decreasing = TRUE),]
# Round numbers in each column
fungi.tax.final.means$MEANS_Total<-round(fungi.tax.final.means$MEANS_Total,0)
fungi.tax.final.means$MEANS_Active<-round(fungi.tax.final.means$MEANS_Active,0)
fungi.tax.final.means$SD_Total<-round(fungi.tax.final.means$SD_Total,1)
fungi.tax.final.means$SD_Active<-round(fungi.tax.final.means$SD_Active,1)
# Calculate difference between Total and Active Means
fungi.tax.final.means$Difference<-fungi.tax.final.means$MEANS_Active - fungi.tax.final.means$MEANS_Total
fungi.tax.final.means$Difference<-round(fungi.tax.final.means$Difference,0)
# Preserve taxon as row names
row.names(fungi.tax.final.means)<-fungi.tax.final.means$Taxon
fungi.tax.final.means<-fungi.tax.final.means[,-1]
#Paste Total means with sd into single column
fungi.tax.final.means$Total<-paste(fungi.tax.final.means$MEANS_Total,"±",fungi.tax.final.means$SD_Total,sep=" ")
#Paste Active means with sd into single column
fungi.tax.final.means$Active<-paste(fungi.tax.final.means$MEANS_Active,"±",fungi.tax.final.means$SD_Active,sep=" ")
# Subset out columns of interest
fungi.table2<-subset(fungi.tax.final.means,select=c("Total","Active","Difference"))
# Merging everything into Table 2
fungi.table2$Taxon<-row.names(fungi.table2)
fungi.man.p.df$Taxon<-row.names(fungi.man.p.df)
#Final table
fungi.table2.final<-merge(fungi.table2,fungi.man.p.df,by="Taxon",all=TRUE)
# Round p-values to 2 decimals
fungi.table2.final$`P-Value`<-round(fungi.table2.final$`P-Value`,3)
# Put P-values in parantheses
fungi.table2.final$P2<-paste("(",fungi.table2.final$`P-Value`,")",sep="")
# Add + signs to Difference factors
for(i in 1:length(fungi.table2.final$Taxon)) {
if (fungi.table2.final[i,4]>0) {
fungi.table2.final[i,4]<-paste("+",fungi.table2.final[i,4],sep="")}
}
# Combine Difference with P-values in single column
fungi.table2.final$`Difference % (P-Value)`<-paste(fungi.table2.final$Difference,fungi.table2.final$P2,sep=" ")
# Remove unnecessary columns from final table
fungi.table2.final<-fungi.table2.final[,-c(4,5,6)]
fungi.table2.final<-fungi.table2.final[c("Taxon","Total","Active","Difference % (P-Value)")]
# Reorder Taxon by Total Rel Abund
fungi.table2.final$Taxon<-as.factor(fungi.table2.final$Taxon)
fungi.table2.final$Taxon<-reorder(fungi.table2.final$Taxon,c(4,3,7,1,6,9,2,5,8))
fungi.table2.final<-fungi.table2.final[order(fungi.table2.final$Taxon),]
# Name Taxon row name
row.names(fungi.table2.final)<-fungi.table2.final$Taxon
# Remove Taxon column
fungi.table2.final<-fungi.table2.final[,-1]| % Total | % Active | % Difference (P-value) | |
|---|---|---|---|
| Leotiomycetes | 45 ± 12.8 | 32 ± 13.8 | -13 (0.019) |
| Sordariomycetes | 20 ± 11 | 15 ± 7.1 | -5 (0.211) |
| Dothideomycetes | 8 ± 2.4 | 12 ± 3.6 | +4 (0.007) |
| Agaricomycetes | 5 ± 3.4 | 12 ± 6.5 | +7 (0.003) |
| Tremellomycetes | 3 ± 1.1 | 6 ± 2.1 | +3 (0) |
| Microbotryomycetes | 2 ± 1.1 | 5 ± 2.4 | +3 (0) |
| Eurotiomycetes | 1 ± 0.9 | 2 ± 1.1 | +1 (0.464) |
| Unclassified | 15 ± 4.8 | 15 ± 4.4 | 0 (0.864) |
| Other | 1 ± 1.3 | 2 ± 2.1 | +1 (0.026) |
Supplemental Table S4. List of “Other” fungi (<1% relative abundance)
fungi.other.names## [1] "Agaricostilbomycetes" "Arthoniomycetes"
## [3] "Ascomycota_incertae_sedis" "Chytridiomycetes"
## [5] "Cystobasidiomycetes" "Exobasidiomycetes"
## [7] "Fungi_incertae_sedis" "Lecanoromycetes"
## [9] "Pezizomycetes" "Pucciniomycetes"
## [11] "Saccharomycetes" "Taphrinomycetes"
## [13] "Ustilaginomycetes"Supplemental Figure S2
Initial BACTERIAL Site*Community setup
# Create new dataframe
bact.tax.sites<-bact.order.merge[, -c(14,15,28,29)]
# Assign row names as Taxon
row.names(bact.tax.sites)<-bact.tax.sites$Taxon
# Remove taxon column
bact.tax.sites<-bact.tax.sites[2:25]
# Convert to matrix
bact.tax.sites.final<-as.data.frame(t(as.matrix(bact.tax.sites)))
# Add in columns for Site and Community
bact.tax.sites.final$Site<-c("A","A","A","B","B","B","C","C","C","D","D","D","A","A","A","B","B","B","C","C","C","D","D","D")
bact.tax.sites.final$Community<-c("Total","Total","Total","Total","Total","Total","Total","Total","Total","Total","Total","Total","Active","Active","Active","Active","Active","Active","Active","Active","Active","Active","Active","Active")
#Create matrix of numeric columns
bact.tax.sites.bind<-as.matrix(cbind(bact.tax.sites.final[1:13]))
bact.site.manova<-manova(bact.tax.sites.bind~Site*Community, data=bact.tax.sites.final)
bact.site.man.sum<-summary.aov(bact.site.manova)
bact.site.man.sum ## Response Acidobacteria_Gp1 :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 1014.28 338.09 22.9278 0.000004899 ***
## Community 1 58.45 58.45 3.9636 0.06386 .
## Site:Community 3 61.58 20.53 1.3919 0.28145
## Residuals 16 235.94 14.75
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Actinomycetales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 257.158 85.719 12.6797 0.0001692 ***
## Community 1 12.568 12.568 1.8591 0.1916089
## Site:Community 3 52.104 17.368 2.5691 0.0906014 .
## Residuals 16 108.166 6.760
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Burkholderiales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 50.796 16.9319 5.2106 0.01060 *
## Community 1 16.655 16.6551 5.1254 0.03783 *
## Site:Community 3 9.608 3.2025 0.9855 0.42441
## Residuals 16 51.992 3.2495
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Caulobacterales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 46.590 15.5301 8.4049 0.001393 **
## Community 1 17.301 17.3009 9.3632 0.007482 **
## Site:Community 3 4.088 1.3626 0.7375 0.544939
## Residuals 16 29.564 1.8478
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Flavobacteriales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 15.6818 5.2273 14.4379 0.0000816 ***
## Community 1 0.0572 0.0572 0.1579 0.6964
## Site:Community 3 0.1383 0.0461 0.1273 0.9425
## Residuals 16 5.7928 0.3621
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Other :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 17.024 5.6747 1.9964 0.1552
## Community 1 0.625 0.6250 0.2199 0.6454
## Site:Community 3 1.120 0.3733 0.1313 0.9400
## Residuals 16 45.479 2.8424
##
## Response Planctomycetales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 55.219 18.4062 4.8858 0.01344 *
## Community 1 23.836 23.8360 6.3271 0.02295 *
## Site:Community 3 1.102 0.3672 0.0975 0.96030
## Residuals 16 60.276 3.7673
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Pseudomonadales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 23.304 7.768 1.0085 0.41467
## Community 1 37.569 37.569 4.8772 0.04215 *
## Site:Community 3 21.519 7.173 0.9312 0.44840
## Residuals 16 123.248 7.703
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Rhizobiales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 36.534 12.178 5.0041 0.0123176 *
## Community 1 41.362 41.362 16.9962 0.0007978 ***
## Site:Community 3 89.106 29.702 12.2049 0.0002085 ***
## Residuals 16 38.938 2.434
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Rhodospirillales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 107.642 35.881 57.8047 0.000000008323 ***
## Community 1 16.454 16.454 26.5081 0.000097100063 ***
## Site:Community 3 1.125 0.375 0.6042 0.6217
## Residuals 16 9.932 0.621
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Sphingobacteriales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 16.827 5.6091 1.8872 0.17245
## Community 1 18.858 18.8584 6.3450 0.02279 *
## Site:Community 3 13.001 4.3337 1.4581 0.26338
## Residuals 16 47.555 2.9722
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Sphingomonadales :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 9.3279 3.10930 2.2127 0.1262
## Community 1 0.1018 0.10179 0.0724 0.7913
## Site:Community 3 2.1452 0.71507 0.5089 0.6818
## Residuals 16 22.4828 1.40518
##
## Response Unclassified :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 93.615 31.2049 5.7873 0.007071 **
## Community 1 3.775 3.7755 0.7002 0.415038
## Site:Community 3 15.370 5.1232 0.9502 0.439887
## Residuals 16 86.271 5.3919
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Setting up plotting data
# Separate Total and Active Communities into dataframes
bact.tax.sites.total<-bact.tax.sites.final[1:12,]
bact.tax.sites.active<-bact.tax.sites.final[13:24,]
# Calculate Site Means for Total
bact.tax.sites.total.means<-aggregate(cbind(Acidobacteria_Gp1,Actinomycetales,Burkholderiales,Caulobacterales,Flavobacteriales,Other,Planctomycetales,Pseudomonadales,Rhizobiales,Rhodospirillales,Sphingobacteriales,Sphingomonadales,Unclassified)~Site+Community,bact.tax.sites.total,mean)
# Create new ID column with Site and Community combined for Active Bacteria
bact.tax.sites.active$ID<-paste(bact.tax.sites.active$Site,bact.tax.sites.active$Community,sep="_")
# Create new ID column with Site and Community combined for Total Means
bact.tax.sites.total.means$ID<-paste(bact.tax.sites.total.means$Site,bact.tax.sites.total.means$Community,sep="_")
# Adding "Total" prefix to column names
colnames(bact.tax.sites.total.means)<-paste("Total", colnames(bact.tax.sites.total.means),sep="_")
bact.tax.sites.total.means$Site<-bact.tax.sites.total.means$Total_Site
bact.tax.sites.total.means<-bact.tax.sites.total.means[c(1:15,17)]
# Merge Total and Active df's by ID
bact.tax.sites.merge<-merge(bact.tax.sites.active,bact.tax.sites.total.means,by="Site", all=TRUE)
# Make new column with Site_Plot_Community
bact.tax.sites.merge$ID2<-c("Active_A_1","Active_A_2","Active_A_3","Active_B_1","Active_B_2","Active_B_3","Active_C_1","Active_C_2","Active_C_3","Active_D_1","Active_D_2","Active_D_3")
# Make ID as rownames
row.names(bact.tax.sites.merge)<-bact.tax.sites.merge$ID2
# Split into two matrices of numeric values
# Active plot abundances
bact.tax.sites.active.abund<-as.matrix(bact.tax.sites.merge[2:14])
# Total mean abundace
bact.tax.sites.total.means.abund<-as.matrix(bact.tax.sites.merge[19:31])
# Calculate site differences
bact.tax.sites.diff<-as.data.frame(bact.tax.sites.active.abund-bact.tax.sites.total.means.abund)
# Add column for Site and Community
bact.tax.sites.diff$Site<-c("A","A","A","B","B","B","C","C","C","D","D","D")
# Aggregate by site mean
bact.tax.sites.diff.means<-aggregate(cbind(Acidobacteria_Gp1,Actinomycetales,Burkholderiales,Caulobacterales,Flavobacteriales,Other,Planctomycetales,Pseudomonadales,Rhizobiales,Rhodospirillales,Sphingobacteriales,Sphingomonadales,Unclassified)~Site,bact.tax.sites.diff,mean)
# Calculate sd
bact.tax.sites.diff.sd<-aggregate(cbind(Acidobacteria_Gp1,Actinomycetales,Burkholderiales,Caulobacterales,Flavobacteriales,Other,Planctomycetales,Pseudomonadales,Rhizobiales,Rhodospirillales,Sphingobacteriales,Sphingomonadales,Unclassified)~Site,bact.tax.sites.diff,sd)
# Melt Mean df
bact.mean.melt<-melt(bact.tax.sites.diff.means,id="Site")
# Rename columns
colnames(bact.mean.melt)<-c("Site","Taxon","Mean_Diff")
# Create Taxon vector
tax.order<-c(13,13,13,13,11,11,11,11,4,4,4,4,10,10,10,10,7,7,7,7,6,6,6,6,3,3,3,3,1,1,1,1,12,12,12,12,2,2,2,2,5,5,5,5,8,8,8,8,9,9,9,9)
bact.mean.melt$Taxon<-reorder(bact.mean.melt$Taxon,tax.order)Build figure
# Initial ggplot figure
bact.tax.sites.diff.fig<-ggplot(bact.mean.melt,aes(x=Taxon,y=Mean_Diff,fill=Site))+geom_bar(width=0.5,colour="black",position=position_dodge(width = 0.5),stat="identity") + geom_hline(yintercept=0) + labs(x="",y="% Change from Total Bacterial Community") + theme(axis.ticks.y = element_blank(), axis.ticks.x = element_blank(), panel.background = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.key = element_blank()) + scale_fill_manual("Site",labels=c("A","B","C","D"),values=c("white","light gray","dark gray","black")) + theme(axis.text.x=element_text(angle=90,vjust=1), panel.border=element_rect(colour = "black",fill=NA)) + theme(axis.text.x=element_text(size=12, color = "black"),axis.text.y=element_text(size = 12, color = "black"),axis.title.y=element_text(size = 12,color = "black"))+ annotate("text",x=12,y=5,label="*",size=12)Initial FUNGAL Site*Community setup
# Create new dataframe
fungi.tax.sites<-fungi.merge[, -c(14,15,28,29)]
# Assign row names as Taxon
row.names(fungi.tax.sites)<-fungi.tax.sites$Taxon
# Remove taxon column
fungi.tax.sites<-fungi.tax.sites[2:25]
# Convert to matrix
fungi.tax.sites.final<-as.data.frame(t(as.matrix(fungi.tax.sites)))
# Add in columns for Site and Community
fungi.tax.sites.final$Site<-c("A","A","A","B","B","B","C","C","C","D","D","D","A","A","A","B","B","B","C","C","C","D","D","D")
fungi.tax.sites.final$Community<-c("Total","Total","Total","Total","Total","Total","Total","Total","Total","Total","Total","Total","Active","Active","Active","Active","Active","Active","Active","Active","Active","Active","Active","Active")
#Create matrix of numeric columns
fungi.tax.sites.bind<-as.matrix(cbind(fungi.tax.sites.final[1:9]))
# MANOVA
fungi.site.manova<-manova(fungi.tax.sites.bind~Site*Community, data=fungi.tax.sites.final)
fungi.site.man.sum<-summary.aov(fungi.site.manova)
fungi.site.man.sum## Response Agaricomycetes :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 388.79 129.596 16.2293 0.00004130 ***
## Community 1 295.51 295.505 37.0063 0.00001583 ***
## Site:Community 3 70.85 23.618 2.9577 0.06388 .
## Residuals 16 127.76 7.985
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Dothideomycetes :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 65.272 21.757 2.6525 0.08397 .
## Community 1 81.935 81.935 9.9889 0.00606 **
## Site:Community 3 10.327 3.442 0.4197 0.74137
## Residuals 16 131.242 8.203
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Eurotiomycetes :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 5.0494 1.68312 1.7268 0.2017
## Community 1 0.5667 0.56673 0.5814 0.4568
## Site:Community 3 1.7655 0.58849 0.6038 0.6220
## Residuals 16 15.5951 0.97469
##
## Response Leotiomycetes :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 3247.5 1082.50 56.508 0.000000009817 ***
## Community 1 1134.8 1134.80 59.238 0.000000912355 ***
## Site:Community 3 351.8 117.28 6.122 0.005638 **
## Residuals 16 306.5 19.16
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Microbotryomycetes :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 34.592 11.531 6.7811 0.003673 **
## Community 1 62.323 62.323 36.6514 0.00001673 ***
## Site:Community 3 16.430 5.477 3.2208 0.050782 .
## Residuals 16 27.207 1.700
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Other :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 7.0159 2.3386 5.8922 0.006582 **
## Community 1 3.8686 3.8686 9.7470 0.006570 **
## Site:Community 3 1.5843 0.5281 1.3306 0.299357
## Residuals 16 6.3504 0.3969
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Sordariomycetes :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 1048.70 349.57 8.2025 0.001563 **
## Community 1 142.65 142.65 3.3472 0.086018 .
## Site:Community 3 163.09 54.36 1.2756 0.316400
## Residuals 16 681.87 42.62
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Tremellomycetes :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 15.726 5.242 3.0927 0.05675 .
## Community 1 65.470 65.470 38.6269 0.00001237 ***
## Site:Community 3 17.637 5.879 3.4686 0.04112 *
## Residuals 16 27.119 1.695
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Unclassified :
## Df Sum Sq Mean Sq F value Pr(>F)
## Site 3 138.601 46.200 3.6744 0.03464 *
## Community 1 0.466 0.466 0.0371 0.84977
## Site:Community 3 1.383 0.461 0.0367 0.99024
## Residuals 16 201.174 12.573
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Setting up plotting data
# Separate Total and Active Communities into dataframes
fungi.tax.sites.total<-fungi.tax.sites.final[1:12,]
fungi.tax.sites.active<-fungi.tax.sites.final[13:24,]
# Calculate Site Means for Total
fungi.tax.sites.total.means<-aggregate(cbind(Agaricomycetes,Dothideomycetes,Eurotiomycetes,Leotiomycetes,Microbotryomycetes,Other,Sordariomycetes,Tremellomycetes,Unclassified)~Site+Community,fungi.tax.sites.total,mean)
# Create new ID column with Site and Community combined for Active Fungi
fungi.tax.sites.active$ID<-paste(fungi.tax.sites.active$Site,fungi.tax.sites.active$Community,sep="_")
# Create new ID column with Site and Community combined for Total Means
fungi.tax.sites.total.means$ID<-paste(fungi.tax.sites.total.means$Site,fungi.tax.sites.total.means$Community,sep="_")
# Adding "Total" prefix to column names
colnames(fungi.tax.sites.total.means)<-paste("Total", colnames(fungi.tax.sites.total.means),sep="_")
fungi.tax.sites.total.means$Site<-fungi.tax.sites.total.means$Total_Site
fungi.tax.sites.total.means<-fungi.tax.sites.total.means[c(1:11,13)]
# Merge Total and Active df's by ID
fungi.tax.sites.merge<-merge(fungi.tax.sites.active,fungi.tax.sites.total.means,by="Site", all=TRUE)
# Make new column with Site_Plot_Community
fungi.tax.sites.merge$ID2<-c("Active_A_1","Active_A_2","Active_A_3","Active_B_1","Active_B_2","Active_B_3","Active_C_1","Active_C_2","Active_C_3","Active_D_1","Active_D_2","Active_D_3")
# Make ID as rownames
row.names(fungi.tax.sites.merge)<-fungi.tax.sites.merge$ID2
# Split into two matrices of numeric values
# Active plot abundances
fungi.tax.sites.active.abund<-as.matrix(fungi.tax.sites.merge[2:10])
# Total mean abundace
fungi.tax.sites.total.means.abund<-as.matrix(fungi.tax.sites.merge[15:23])
# Calculate site differences
fungi.tax.sites.diff<-as.data.frame(fungi.tax.sites.active.abund-fungi.tax.sites.total.means.abund)
# Add column for Site
fungi.tax.sites.diff$Site<-c("A","A","A","B","B","B","C","C","C","D","D","D")
# Aggregate by site mean
fungi.tax.sites.diff.means<-aggregate(cbind(Agaricomycetes,Dothideomycetes,Eurotiomycetes,Leotiomycetes,Microbotryomycetes,Other,Sordariomycetes,Tremellomycetes,Unclassified)~Site,fungi.tax.sites.diff,mean)
# Calculate sd
fungi.tax.sites.diff.sd<-aggregate(cbind(Agaricomycetes,Dothideomycetes,Eurotiomycetes,Leotiomycetes,Microbotryomycetes,Other,Sordariomycetes,Tremellomycetes,Unclassified)~Site,fungi.tax.sites.diff,sd)
# Melt Mean df
fungi.mean.melt<-melt(fungi.tax.sites.diff.means,id="Site")
# Rename columns
colnames(fungi.mean.melt)<-c("Site","Taxon","Mean_Diff")
# Create Taxon vector
tax.fungi.order<-c(1,1,1,1,4,4,4,4,7,7,7,7,9,9,9,9,2,2,2,2,5,5,5,5,8,8,8,8,3,3,3,3,6,6,6,6)
fungi.mean.melt$Taxon<-reorder(fungi.mean.melt$Taxon,tax.fungi.order)Build figure
# Initial ggplot figure
fungi.tax.sites.diff.fig<-ggplot(fungi.mean.melt,aes(x=Taxon,y=Mean_Diff,fill=Site))+geom_bar(width=0.5,colour="black",position=position_dodge(width = 0.5),stat="identity") + geom_hline(yintercept=0) + labs(x="",y="% Change from Total Fungal Community") + theme(axis.ticks.y = element_blank(), axis.ticks.x = element_blank(), panel.background = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.key = element_blank()) + scale_fill_manual("Site",labels=c("A","B","C","D"),values=c("white","light gray","dark gray","black")) + theme(axis.text.x=element_text(angle=90,vjust=1), panel.border=element_rect(colour = "black",fill=NA)) + theme(axis.text.x=element_text(size=12, color = "black"),axis.text.y=element_text(size = 12, color = "black"),axis.title.y=element_text(size = 12,color = "black"))+ annotate("text",x=3,y=8,label="*",size=12) + annotate("text",x=9,y=8,label="*",size=12)Recreating Supplemental Figure S2.
grid.arrange(bact.tax.sites.diff.fig,fungi.tax.sites.diff.fig,ncol=1)
Supplemental Figure S2. Percent change in relative abundance from the Total community by Sites.
Environmental Influence
Transform environmental data
#Z-score transform environmental data
new.env.qual<-new.env
new.env.qual<-new.env.qual[1:9]
new.env.qual.z<-scale(new.env.qual)
new.env.qual.z<-as.data.frame(new.env.qual.z)Test for collinearity of environmental data (|r| > 0.80)
env.corr<-cor(new.env.qual.z)
# Round to 2 decimal places
env.corr<-round(env.corr,2)| Moisture | pH | Aryl | Alkyl | O_alkyl | Carboxyl | C | N | CN | |
|---|---|---|---|---|---|---|---|---|---|
| Moisture | 1 | 0.78 | -0.58 | 0.52 | 0.09 | -0.78 | -0.57 | 0.72 | -0.8 |
| pH | 0.78 | 1 | -0.22 | 0.24 | 0.03 | -0.4 | -0.59 | 0.76 | -0.77 |
| Aryl | -0.58 | -0.22 | 1 | -0.51 | -0.38 | 0.91 | 0.1 | -0.03 | 0.1 |
| Alkyl | 0.52 | 0.24 | -0.51 | 1 | -0.58 | -0.47 | -0.27 | 0.34 | -0.41 |
| O_alkyl | 0.09 | 0.03 | -0.38 | -0.58 | 1 | -0.41 | 0.12 | -0.24 | 0.24 |
| Carboxyl | -0.78 | -0.4 | 0.91 | -0.47 | -0.41 | 1 | 0.25 | -0.25 | 0.33 |
| C | -0.57 | -0.59 | 0.1 | -0.27 | 0.12 | 0.25 | 1 | -0.49 | 0.73 |
| N | 0.72 | 0.76 | -0.03 | 0.34 | -0.24 | -0.25 | -0.49 | 1 | -0.95 |
| CN | -0.8 | -0.77 | 0.1 | -0.41 | 0.24 | 0.33 | 0.73 | -0.95 | 1 |
Taxonomy Correlations with Environmental Factors
# Remove Mean and SD Columns in Bacteria Total Community
bact.order.total.final.corr<-bact.order.total.final[1:12]
# Remove Mean and SD Columns in Bacteria Active Community
bact.order.active.final.corr<-bact.order.active.final[1:12]
# Remove non-significant taxa from further downstream correlations
bact.order.total.final.corr<-bact.order.total.final.corr[-c(1:3,5,9,11:13),]
# Remove non-significant taxa from further downstream correlations
bact.order.active.final.corr<-bact.order.active.final.corr[-c(1:3,5,9,11:13),]
# Remove Mean and SD Columns in Fungal Total Community
fungi.total.final.corr<-fungi.total.final[1:12]
# Remove Mean and SD Columns in Fungal Active Community
fungi.active.final.corr<-fungi.active.final[1:12]
# Remove non-significant taxa from further downstream correlations
fungi.total.final.corr<-fungi.total.final.corr[-c(3,6,9),]
# Remove non-significant taxa from further downstream correlations
fungi.active.final.corr<-fungi.active.final.corr[-c(3,6,9),]
# Combine taxa tables into a list
taxa.final.corr.list1<-list(bact.order.total.final.corr,bact.order.active.final.corr,fungi.total.final.corr,fungi.active.final.corr)
# Rename list
names(taxa.final.corr.list1)<- c("bact.total","bact.active","fungi.total","fungi.active")
# Convert to matrices
for (i in 1:length(taxa.final.corr.list1)) {
taxa.final.corr.list1[[i]] <- as.data.frame(t(as.matrix(taxa.final.corr.list1[[i]])))
}
# Rename row names to match env file
row.names(taxa.final.corr.list1[[1]])<-c("Total_A_1","Total_A_2","Total_A_3","Total_B_1","Total_B_2","Total_B_3","Total_C_1","Total_C_2","Total_C_3","Total_D_1","Total_D_2","Total_D_3")
row.names(taxa.final.corr.list1[[2]])<-c("Active_A_1","Active_A_2","Active_A_3","Active_B_1","Active_B_2","Active_B_3","Active_C_1","Active_C_2","Active_C_3","Active_D_1","Active_D_2","Active_D_3")
row.names(taxa.final.corr.list1[[3]])<-c("Total_A_1","Total_A_2","Total_A_3","Total_B_1","Total_B_2","Total_B_3","Total_C_1","Total_C_2","Total_C_3","Total_D_1","Total_D_2","Total_D_3")
row.names(taxa.final.corr.list1[[4]])<-c("Active_A_1","Active_A_2","Active_A_3","Active_B_1","Active_B_2","Active_B_3","Active_C_1","Active_C_2","Active_C_3","Active_D_1","Active_D_2","Active_D_3")
# Convert taxa dataframes to matrices
for (i in 1:length(taxa.final.corr.list1)) {
taxa.final.corr.list1[[i]] <- as.matrix(taxa.final.corr.list1[[i]])
}
# Separate out total env factors
env.total<-as.matrix(env[1:12,])
# Separate out active env factors
env.active<-as.matrix(env[13:24,])
# Z-Score Env Factors
env.total.z<-scale(env.total)
# Z-Score Env Factors
env.active.z<-scale(env.active)
# Final Bacteria Total Correlation
bact.total.corr.final<-cor(taxa.final.corr.list1[[1]],env.total.z,method = "spearman")
# Final Bacteria Active Correlation
bact.active.corr.final<-cor(taxa.final.corr.list1[[2]],env.active.z,method = "spearman")
# Final Fungi Total Correlation
fungi.total.corr.final<-cor(taxa.final.corr.list1[[3]],env.total.z,method = "spearman")
# Final Fungi Active Correlation
fungi.active.corr.final<-cor(taxa.final.corr.list1[[4]],env.active.z,method = "spearman")Recreating Table 5 - BACTERIAL Taxa Correlations to Environmental factors
#TOTAL BACTERIA
bact.total.corr.final<-round(bact.total.corr.final,2)
#ACTIVE BACTERIA
bact.active.corr.final<-round(bact.active.corr.final,2)
# Rename rows in the data frame to environmental variables
row.names(bact.total.corr.final)<-c("Caulobacterales","Planctomycetales","Pseudomonadales","Rhizobiales","Sphingobacteriales")
# Rename rows in the data frame to environmental variables
row.names(bact.active.corr.final)<-c("Caulobacterales","Planctomycetales","Pseudomonadales","Rhizobiales","Sphingobacteriales")
# Remove collinear environmental factors (Aryl and C:N) and non-significant factors from marginal DistLM (Alkyl,O/N-alkyl, C)
bact.total.corr.final<-bact.total.corr.final[,-c(3,4,5,7,9)]
# Remove collinear environmental factors (Aryl and C:N) and non-significant factors from marginal DistLM (Alkyl and O/N-alkyl)
bact.active.corr.final<-bact.active.corr.final[,-c(3,4,5,9)]| Moisture | pH | Carboxyl | N | |
|---|---|---|---|---|
| Caulobacterales | -0.5 | -0.71 | 0.5 | -0.5 |
| Planctomycetales | 0.26 | 0.57 | -0.29 | 0.34 |
| Pseudomonadales | -0.38 | -0.13 | 0.61 | -0.1 |
| Rhizobiales | 0.46 | 0.55 | -0.57 | 0.45 |
| Sphingobacteriales | 0.12 | -0.13 | 0.4 | 0.21 |
| Moisture | pH | Carboxyl | C | N | |
|---|---|---|---|---|---|
| Caulobacterales | -0.82 | -0.66 | 0.41 | 0.59 | -0.92 |
| Planctomycetales | 0.32 | 0.37 | -0.49 | -0.43 | 0.27 |
| Pseudomonadales | 0.22 | 0.28 | -0.05 | -0.31 | 0.37 |
| Rhizobiales | -0.7 | -0.66 | 0.32 | 0.33 | -0.73 |
| Sphingobacteriales | 0.29 | 0.21 | -0.71 | 0.03 | 0.13 |
Recreating Table 5 - FUNGAL Taxa Correlations to Environmental factors
#TOTAL FUNGI
fungi.total.corr.final<-round(fungi.total.corr.final,2)
#ACTIVE FUNGI
fungi.active.corr.final<-round(fungi.active.corr.final,2)
# Rename rows in the data frame to environmental variables
row.names(fungi.total.corr.final)<-c("Agaricomycetes","Dothideomycetes","Leotiomycetes","Microbotryomycetes","Tremellomycetes","Other")
# Rename rows in the data frame to environmental variables
row.names(fungi.active.corr.final)<-c("Agaricomycetes","Dothideomycetes","Leotiomycetes","Microbotryomycetes","Tremellomycetes","Other")
# Remove collinear environmental factors (Aryl and C:N) and non-significant factors from marginal DistLM (Alkyl,O/N-alkyl, C)
fungi.total.corr.final<-fungi.total.corr.final[,-c(3,4,5,7,9)]
# Remove collinear environmental factors (Aryl and C:N) and non-significant factors from marginal DistLM (Alkyl and O/N-alkyl)
fungi.active.corr.final<-fungi.active.corr.final[,-c(3,4,5,9)]| Moisture | pH | Carboxyl | N | |
|---|---|---|---|---|
| Agaricomycetes | 0.57 | 0.71 | -0.58 | 0.69 |
| Dothideomycetes | 0.34 | 0.34 | -0.55 | 0.31 |
| Leotiomycetes | -0.41 | -0.36 | 0.71 | -0.13 |
| Microbotryomycetes | 0.33 | 0.37 | -0.54 | 0.38 |
| Tremellomycetes | -0.2 | -0.23 | -0.05 | -0.24 |
| Other | 0.38 | 0.55 | -0.56 | 0.29 |
| Moisture | pH | Carboxyl | C | N | |
|---|---|---|---|---|---|
| Agaricomycetes | 0.71 | 0.64 | -0.52 | -0.44 | 0.83 |
| Dothideomycetes | 0.53 | 0.21 | -0.28 | -0.09 | 0.51 |
| Leotiomycetes | -0.29 | -0.46 | 0.56 | 0.34 | -0.13 |
| Microbotryomycetes | 0.24 | 0.56 | -0.3 | -0.57 | 0.44 |
| Tremellomycetes | 0.01 | 0.38 | -0.27 | -0.55 | 0.18 |
| Other | 0.34 | 0.73 | -0.48 | -0.46 | 0.34 |
Reproducibility
# Provide the Session Information for reproducibility
devtools::session_info()## setting value
## version R version 3.2.4 (2016-03-10)
## system x86_64, darwin13.4.0
## ui X11
## language (EN)
## collate en_US.UTF-8
## tz America/Detroit
## date 2016-06-21
##
## package * version date source
## chron 2.3-47 2015-06-24 CRAN (R 3.2.0)
## cluster 2.0.4 2016-04-18 CRAN (R 3.2.5)
## colorspace 1.2-6 2015-03-11 CRAN (R 3.2.0)
## data.table * 1.9.6 2015-09-19 CRAN (R 3.2.0)
## devtools * 1.11.1 2016-04-21 CRAN (R 3.2.5)
## digest 0.6.9 2016-01-08 CRAN (R 3.2.3)
## evaluate 0.8.3 2016-03-05 CRAN (R 3.2.4)
## formatR 1.3 2016-03-05 CRAN (R 3.2.4)
## ggplot2 * 2.1.0 2016-03-01 CRAN (R 3.2.4)
## gridExtra * 2.2.1 2016-02-29 CRAN (R 3.2.4)
## gtable 0.2.0 2016-02-26 CRAN (R 3.2.3)
## htmltools 0.3.5 2016-03-21 CRAN (R 3.2.4)
## knitr 1.12.3 2016-01-22 CRAN (R 3.2.3)
## labeling 0.3 2014-08-23 CRAN (R 3.2.0)
## lattice * 0.20-33 2015-07-14 CRAN (R 3.2.4)
## magrittr 1.5 2014-11-22 CRAN (R 3.2.0)
## MASS 7.3-45 2015-11-10 CRAN (R 3.2.4)
## Matrix 1.2-5 2016-04-17 CRAN (R 3.2.5)
## memoise 1.0.0 2016-01-29 CRAN (R 3.2.3)
## mgcv 1.8-12 2016-03-03 CRAN (R 3.2.4)
## munsell 0.4.3 2016-02-13 CRAN (R 3.2.3)
## nlme 3.1-127 2016-04-16 CRAN (R 3.2.5)
## pander * 0.6.0 2015-11-23 CRAN (R 3.2.2)
## permute * 0.9-0 2016-01-24 CRAN (R 3.2.3)
## plyr * 1.8.3 2015-06-12 CRAN (R 3.2.0)
## Rcpp 0.12.4 2016-03-26 CRAN (R 3.2.4)
## reshape2 * 1.4.1 2014-12-06 CRAN (R 3.2.0)
## rmarkdown 0.9.5 2016-02-22 CRAN (R 3.2.3)
## scales * 0.4.0 2016-02-26 CRAN (R 3.2.3)
## stringi 1.0-1 2015-10-22 CRAN (R 3.2.0)
## stringr 1.0.0 2015-04-30 CRAN (R 3.2.0)
## vegan * 2.3-5 2016-04-09 CRAN (R 3.2.4)
## withr 1.0.1 2016-02-04 CRAN (R 3.2.3)
## yaml 2.1.13 2014-06-12 CRAN (R 3.2.0)