AMF
C. Nell
20180821
Logic
Salvage topsoil was removed from a donor site and delivered to three recipient sites in 2015. We sampled all sites (Donor and recipients) prior to topsoil delivery and sequenced AM Fungi from all sites in 2015 and 2017 from manipulative plots, set up in a randomized block design at each site with: control treatments (no topsoil added), a dusting of topsoil, and three levels of topsoil thicknesses (2“, 4” and 6" thick layers of this topsoil, originating from the donor site)), delivered to all recipient sites).
knitr::opts_chunk$set(cache=TRUE, tidy=TRUE, error = FALSE, eval = TRUE, message = FALSE, warning = FALSE, rows.print=5, cols.min.print=4, fig.width=6, fig.height=4.5)Questions:
**Q1: Propagule determination: Do AM fungal communities resemble the donor site, regardless of where soil was delivered? Is provenance a driver of AM fungal community composition?
**Q2: Environmental filtering: Do AM fungal communities resemble the recipient sites, and are more dissimilar to the AM fungal communities from the donor site?
**Q3: Propagule pressure: Do AM fungal communities resemble the donor site more in higher level topsoil treatments, than they do in the control or dusted sites?
Outputs:
Analysis of similarity:
Determine the similarity of AMF communities at each recipient site to the donor site, from where topsoils originated from. Generate heatmaps with distances, to visually determine - which sites/samples are more similar to each other, and which are more dissimilar. Generate NMDS / Pcoa, plotted with Year as shapes and Sites as colors, and Description as …? Run permanova Do AM fungal communities at each site resemble the recipient sites, and are they more dissimilar to the donor sites?
Functional group richness:
Determine the richness of AMF functional groups in the donor site, from where topsoils originated from, and in the recipient sites pre-topsoil delivery (pre-treatment), and post-topsoil delivery (post-treatment)
Taxonomic diversity:
Determine the OTU taxa richness (alpha diversity) of AMF communities in the donor site, from where topsoils originated from, and in the recipient sites pre-topsoil delivery (pre-treatment), and post-topsoil delivery (post-treatment). Determine the beta diversity of AMF communities in the donor site, from where topsoils originated from, and in the recipient sites pre-topsoil delivery (pre-treatment), and post-topsoil delivery (post-treatment).
Topsoil level treatments
Propagule pressure: Do AM fungal communities resemble the donor site more in thicker topsoil layer treatments than they do in the control or dusted sites? Are the sites with Treat groups = S, more similar to the donor site than F, and are treatment groups S and F more similar to the donor site than T? Is this relationship clinal, with S being the thickest and most similar to donor site, followed by F, and then by T. TO FIGURE OUT - how to address this question statistically??multiple regression with Treat on the X and similarity to Donor Site on the Y?? A permanova with multiple comparisons?
knitr::opts_chunk$set(cache = TRUE, tidy = TRUE, error = FALSE, eval = TRUE,
message = FALSE, warning = FALSE, rows.print = 5, cols.min.print = 4, fig.width = 6,
fig.height = 4.5)Load Required Packages
library(tidyverse)
library(dplyr) ## for data wrangling - %>% function
library(reshape2) ##melt and cast data
library(ggplot2) # plotting
library(data.table)
library(stringr)
library(tidyr) # 'separate' function
library(readxl) #read xlsx files into r on mac computer
library(vegan) # dissimilarity matrix, permanova functions
library(magrittr)
library(cowplot)
library(formatR)
date <- format(Sys.time(), "%Y%b%d")
date## [1] "2018Aug21"** the following is adapted from 2017-07-07_Chapter_2_SSU_IntialClean **
Import data files
# read in OTU table
ssu <- fread("CSS_table_sorted.txt")
str(ssu)## Classes 'data.table' and 'data.frame': 386 obs. of 63 variables:
## $ #OTU ID : chr "denovo538" "denovo1813859" "denovo1031897" "denovo34" ...
## $ MM933M : num 5.29 4.33 7.32 0 0 ...
## $ MM936M : num 6.27 11.24 3.4 0 0 ...
## $ MM941M : num 4.4 2.95 0 0 0 ...
## $ MM943M : num 5.31 4.04 0 0 0 ...
## $ MM9362CU: num 3.13 3.13 8.77 0 1.56 ...
## $ MMS700 : num 2.5 6.82 7.82 0 2.5 ...
## $ MMS701 : num 1.87 2.66 8.27 2.66 0 ...
## $ MMS702 : num 0 2.79 8.24 2.79 0 ...
## $ MMS703 : num 0 2.58 7.18 3.45 0 ...
## $ MMS704 : num 2.75 0 7.43 9.75 0 ...
## $ MMS705 : num 3.28 0 6.99 0 4.76 ...
## $ MMS706 : num 1.83 2.61 8.32 0 0 ...
## $ MMS707 : num 3.03 7.49 7.43 12.04 0 ...
## $ MMS708 : num 0 3.88 7.64 2.97 0 ...
## $ MMS709 : num 5.69 8.79 7.95 0 1.56 ...
## $ MMS710 : num 3.98 3.07 8.5 2.23 0 ...
## $ MMS711 : num 4.44 1.99 8.14 1.99 0 ...
## $ MMS712 : num 2.23 3.07 5.91 7.39 4.75 ...
## $ MMS713 : num 2.87 2.05 8.23 2.05 0 ...
## $ MMS714 : num 4.03 4.34 6.54 2.28 8.78 ...
## $ MMS715 : num 3.03 6.25 4.71 9.76 3.03 ...
## $ MMS716 : num 4.22 2.43 8.93 10.72 0 ...
## $ MMS717 : num 3.12 5.56 8.23 2.28 0 ...
## $ MMS890 : num 4.73 2.89 9.24 0 0 ...
## $ MMS891 : num 4.35 2.55 8.66 0 0 ...
## $ MMS900 : num 0 2.95 7.12 2.95 0 ...
## $ MMS901 : num 2.56 6.47 6.11 2.56 0 ...
## $ MMS902 : num 4.22 7.82 6.94 7.64 0 ...
## $ MMS903 : num 4.52 6.31 9.14 2.86 0 ...
## $ MMS904 : num 3.7 2.32 9.18 2.32 0 ...
## $ MMS905 : num 2.95 6.57 7.42 0 0 ...
## $ MMS906 : num 3.5 3.5 9.15 0 0 ...
## $ MMS907 : num 3.54 4.88 7.56 4.69 0 ...
## $ MMS908 : num 5.26 6.77 9 3.22 0 ...
## $ MMS909 : num 3.76 0 6.25 0 0 ...
## $ MMS910 : num 3.4 0 7.02 6.86 0 ...
## $ MMS911 : num 0 0 7.01 4.39 5.73 ...
## $ MMS912 : num 3.01 2.18 8.48 0 0 ...
## $ MMS913 : num 3.94 1.95 9.29 1.95 0 ...
## $ MMS914 : num 3.83 0 7.66 2.44 4.53 ...
## $ MMS915 : num 0 3.72 6.21 2.02 0 ...
## $ MMS916 : num 4.04 2.03 7.68 4.69 0 ...
## $ MMS917 : num 2.58 0 0 0 7.61 ...
## $ MMS918 : num 2.78 8.18 0 0 0 ...
## $ MMS919 : num 3.79 2.4 8.47 0 0 ...
## $ MMS920 : num 4.31 2.26 6.97 3.1 0 ...
## $ MMS921 : num 4.16 2.6 7.69 3.68 1.43 ...
## $ MMS922 : num 3.6 1.92 7.99 2.71 1.92 ...
## $ MMS923 : num 3.55 3.03 8.85 3.55 0 ...
## $ MMS924 : num 3.68 0 8.09 2.79 0 ...
## $ MMS925 : num 3.69 2.31 7.78 3.16 0 ...
## $ MMS928 : num 3.4 3.4 7.87 2.53 0 ...
## $ MMS929 : num 4.19 4.89 9.17 1.95 0 ...
## $ MMS930 : num 4.75 3.81 9.89 0 0 ...
## $ MMS931 : num 0 4.61 8.61 0 0 ...
## $ MMS932 : num 4.63 0 9.28 0 0 ...
## $ MMS938 : num 2.5 2.5 8.11 0 0 ...
## $ MMS939 : num 6.1 3.73 8.08 0 0 ...
## $ MMS942 : num 0 0 8.1 2.97 0 ...
## $ MMS943 : num 0 2.91 7.07 4.92 0 ...
## $ MMS944 : num 5.13 3.87 8.23 2.96 0 ...
## $ taxonomy: chr "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MH3_B" "No blast hit" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27" ...
## - attr(*, ".internal.selfref")=<externalptr>## mapping data
metassu <- fread("map_SSU_Salvage_qiimeformat.txt", header = TRUE)
metassu$Year <- as.character(metassu$Year)
metassu$Rep <- as.factor(metassu$Rep)
metassu$Treat <- as.factor(metassu$Treat)
str(metassu)## Classes 'data.table' and 'data.frame': 62 obs. of 12 variables:
## $ #SampleID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ BarcodeSequence : chr "GCAACCTTTGGGTGAT" "GCAACCTTTACTGTGC" "GCAACCTTAGAGTGTG" "GCAACCTTTCGATTGG" ...
## $ LinkerPrimerSequence: logi NA NA NA NA NA NA ...
## $ Location : chr "A1" "B1" "C1" "D1" ...
## $ Project : chr "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "WL" "WL" "WL" "WL" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 4 1 6 5 1 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 1 2 3 4 5 1 1 2 3 1 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
## $ SampleID2 : chr "MM941M" "MM943M" "MM936M" "MM9362CU" ...
## - attr(*, ".internal.selfref")=<externalptr>Clean SSU headers
colnames(ssu)## [1] "#OTU ID" "MM933M" "MM936M" "MM941M" "MM943M" "MM9362CU"
## [7] "MMS700" "MMS701" "MMS702" "MMS703" "MMS704" "MMS705"
## [13] "MMS706" "MMS707" "MMS708" "MMS709" "MMS710" "MMS711"
## [19] "MMS712" "MMS713" "MMS714" "MMS715" "MMS716" "MMS717"
## [25] "MMS890" "MMS891" "MMS900" "MMS901" "MMS902" "MMS903"
## [31] "MMS904" "MMS905" "MMS906" "MMS907" "MMS908" "MMS909"
## [37] "MMS910" "MMS911" "MMS912" "MMS913" "MMS914" "MMS915"
## [43] "MMS916" "MMS917" "MMS918" "MMS919" "MMS920" "MMS921"
## [49] "MMS922" "MMS923" "MMS924" "MMS925" "MMS928" "MMS929"
## [55] "MMS930" "MMS931" "MMS932" "MMS938" "MMS939" "MMS942"
## [61] "MMS943" "MMS944" "taxonomy"head(ssu$taxonomy)## [1] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27"
## [2] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MH3_B"
## [3] "No blast hit"
## [4] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27"
## [5] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__sp"
## [6] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Acaulosporaceae; g__Acaulospora; s__Acau16"# rename columns
names(ssu)[1] <- "ssuotu" #rename first column
names(ssu)[names(ssu) == "taxonomy"] <- "ssutaxonomy" # rename column that is currently called taxonomy
# split taxonomy column
`?`(str_match)
ssu$ssukingdom <- str_match(ssu$ssutaxonomy, "k__(.*?);")[, 2]
ssu$ssuphylum <- str_match(ssu$ssutaxonomy, "p__(.*?);")[, 2]
ssu$ssuclass <- str_match(ssu$ssutaxonomy, "c__(.*?);")[, 2]
ssu$ssuorder <- str_match(ssu$ssutaxonomy, "o__(.*?);")[, 2]
ssu$ssufamily <- str_match(ssu$ssutaxonomy, "f__(.*?);")[, 2]
ssu$ssugenus <- str_match(ssu$ssutaxonomy, "g__(.*?);")[, 2]
ssu$ssuspecies <- str_match(ssu$ssutaxonomy, "s__(.*?)")[, 2]
colnames(ssu)## [1] "ssuotu" "MM933M" "MM936M" "MM941M" "MM943M"
## [6] "MM9362CU" "MMS700" "MMS701" "MMS702" "MMS703"
## [11] "MMS704" "MMS705" "MMS706" "MMS707" "MMS708"
## [16] "MMS709" "MMS710" "MMS711" "MMS712" "MMS713"
## [21] "MMS714" "MMS715" "MMS716" "MMS717" "MMS890"
## [26] "MMS891" "MMS900" "MMS901" "MMS902" "MMS903"
## [31] "MMS904" "MMS905" "MMS906" "MMS907" "MMS908"
## [36] "MMS909" "MMS910" "MMS911" "MMS912" "MMS913"
## [41] "MMS914" "MMS915" "MMS916" "MMS917" "MMS918"
## [46] "MMS919" "MMS920" "MMS921" "MMS922" "MMS923"
## [51] "MMS924" "MMS925" "MMS928" "MMS929" "MMS930"
## [56] "MMS931" "MMS932" "MMS938" "MMS939" "MMS942"
## [61] "MMS943" "MMS944" "ssutaxonomy" "ssukingdom" "ssuphylum"
## [66] "ssuclass" "ssuorder" "ssufamily" "ssugenus" "ssuspecies"ssu.save <- ssu
# remove Geosiphonaceae
unique(ssu$ssufamily)## [1] "Glomeraceae" NA "Acaulosporaceae"
## [4] "Paraglomeraceae" "Archaeosporaceae" "Ambisporaceae"
## [7] "Claroideoglomeraceae" "Diversisporaceae" "Geosiphonaceae"drop.family <- c("Geosiphonaceae")
ssu <- ssu[-which(ssu$ssufamily %in% drop.family), ]
# remove no blast hits
unique(ssu$ssutaxonomy)## [1] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27"
## [2] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MH3_B"
## [3] "No blast hit"
## [4] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__sp"
## [5] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Acaulosporaceae; g__Acaulospora; s__Acau16"
## [6] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G47"
## [7] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Paraglomerales; f__Paraglomeraceae; g__Paraglomus; s__laccatum"
## [8] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__acnaGlo2"
## [9] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Archaeosporaceae; g__Archaeospora; s__MO_Ar1"
## [10] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Archaeosporaceae; g__Archaeospora; s__sp"
## [11] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Ambisporaceae; g__Ambispora; s__leptoticha"
## [12] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Whitfield_type_7"
## [13] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__Torrecillas12b_Glo_G5"
## [14] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Sanchez_Castro12b_GLO12"
## [15] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Winther07_D"
## [16] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G18"
## [17] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__intraradices"
## [18] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Paraglomerales; f__Paraglomeraceae; g__Paraglomus; s__Alguacil12a_Para_1"
## [19] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__Douhan9"
## [20] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__Alguacil12b_GLO_G3"
## [21] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Glo39"
## [22] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Alguacil10_Glo1"
## [23] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Ligrone07_sp"
## [24] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Acaulosporaceae; g__Kuklospora; s__PT6"
## [25] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Diversisporaceae; g__Diversispora; s__sp"
## [26] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G41"
## [27] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__MO_GB1"
## [28] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__A1"
## [29] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Liu2012b_Phylo_5"
## [30] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Yamato09_C1"
## [31] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__acnaGlo7"
## [32] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G7"
## [33] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Alguacil10_Glo6"
## [34] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Paraglomerales; f__Paraglomeraceae; g__Paraglomus; s__Alguacil12b_PARA1"
## [35] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__Glo58"
## [36] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__lamellosum"
## [37] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__VeGlo18"
## [38] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Ambisporaceae; g__Ambispora; s__Liu2012a_Ar_1"
## [39] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Diversisporaceae; g__Diversispora; s__Div"
## [40] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Ambisporaceae; g__Ambispora; s__fennica"
## [41] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Voyriella_parviflora_symbiont_type_1"
## [42] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Paraglomerales; f__Paraglomeraceae; g__Paraglomus; s__sp"
## [43] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Paraglomerales; f__Paraglomeraceae; g__Paraglomus; s__Para1_OTU2"
## [44] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Whitfield_type_3"
## [45] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Archaeosporaceae; g__Archaeospora; s__Aca"
## [46] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G20"
## [47] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Acaulosporaceae; g__Acaulospora; s__sp"
## [48] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Glo7"
## [49] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Diversisporaceae; g__Diversispora; s__MO_GC1"
## [50] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Winther07_B"
## [51] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__mosseae"
## [52] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__caledonium"
## [53] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Paraglomerales; f__Paraglomeraceae; g__Paraglomus; s__brasilianum"
## [54] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Alguacil09b_Glo_G8"
## [55] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Archaeosporaceae; g__Archaeospora; s__Schechter08_Arch1"
## [56] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Archaeosporaceae; g__Archaeospora; s__Wirsel_OTU21"
## [57] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__Glo71"
## [58] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Diversisporaceae; g__Diversispora; s__spurca"
## [59] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__Glo59"
## [60] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Wirsel_OTU16"
## [61] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Glo32"
## [62] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Acaulosporaceae; g__Acaulospora; s__MO_A8"
## [63] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G27"
## [64] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G5"drop <- c("No blast hit")
ssu <- ssu[-which(ssu$ssutaxonomy %in% drop), ]
## alternatively, can use filtering to drop these and avoid rewriting things
ssu <- ssu.save %>% filter(ssufamily != "Geosiphonaceae" & ssutaxonomy != "No blast hit")
ssu <- ssu[complete.cases(ssu), ]Add in functional groups
# this does the same as above but avoids hard coding names, more generalized
ssu.fg <- ssu %>% mutate(Functional.Group = case_when(ssufamily %in% c("Glomeraceae",
"Claroideoglomeraceae", "Paraglomeraceae") ~ "Rhizophilic", ssufamily %in%
c("Gigasporaceae", "Diversisporaceae") ~ "Edaphophilic", ssufamily %in%
c("Ambisporaceae", "Archaeosporaceae", "Acaulosporaceae") ~ "Ancestral"))Take data from wide to long and renaming variables
ssul <- melt(ssu.fg, variable.name = "ID", value.name = "ssureads")
str(ssul)## 'data.frame': 17568 obs. of 12 variables:
## $ ssuotu : chr "denovo538" "denovo1813859" "denovo34" "denovo11952" ...
## $ ssutaxonomy : chr "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MH3_B" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__sp" ...
## $ ssukingdom : chr "Fungi" "Fungi" "Fungi" "Fungi" ...
## $ ssuphylum : chr "Glomeromycota" "Glomeromycota" "Glomeromycota" "Glomeromycota" ...
## $ ssuclass : chr "Glomeromycetes" "Glomeromycetes" "Glomeromycetes" "Glomeromycetes" ...
## $ ssuorder : chr "Glomerales" "Glomerales" "Glomerales" "Glomerales" ...
## $ ssufamily : chr "Glomeraceae" "Glomeraceae" "Glomeraceae" "Glomeraceae" ...
## $ ssugenus : chr "Glomus" "Glomus" "Glomus" "Glomus" ...
## $ ssuspecies : chr "" "" "" "" ...
## $ Functional.Group: chr "Rhizophilic" "Rhizophilic" "Rhizophilic" "Rhizophilic" ...
## $ ID : Factor w/ 61 levels "MM933M","MM936M",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ ssureads : num 5.29 4.33 0 0 3.78 ...Functional group OTU richnesses and read abundance for each functional group
Below is in format for graphing
ssuguildREADRICH <- data.frame(ssul %>% group_by(ID, Functional.Group) %>% summarise(OTU_Richness_Sample = length(unique(ssuotu[ssureads >
0])), Read_Abundance_Sample = sum(ssureads)))Make SSU data frame with family level OTU richness and Read abundance for each family
ssufamily <- data.frame(ssul %>% group_by(ID, ssufamily) %>% summarise(Family_OTU_Richness = length(unique(ssuotu[ssureads >
0])), Family_Read_Abundance = sum(ssureads)))Add metadata into ssu family
colnames(metassu) # need to edit for appropriate headers## [1] "#SampleID" "BarcodeSequence" "LinkerPrimerSequence"
## [4] "Location" "Project" "Year"
## [7] "Site" "Treat" "Rep"
## [10] "TSDel" "Description" "SampleID2"names(metassu)[names(metassu) == "SampleID2"] <- "ID"
ssufamily <- ssufamily %>% left_join(metassu, by = "ID")
str(ssufamily)## 'data.frame': 427 obs. of 15 variables:
## $ ID : chr "MM933M" "MM933M" "MM933M" "MM933M" ...
## $ ssufamily : chr "Acaulosporaceae" "Ambisporaceae" "Archaeosporaceae" "Claroideoglomeraceae" ...
## $ Family_OTU_Richness : int 4 5 17 15 2 65 22 3 5 12 ...
## $ Family_Read_Abundance: num 12.4 23.9 157.6 99.6 23.2 ...
## $ #SampleID : int 5 5 5 5 5 5 5 3 3 3 ...
## $ BarcodeSequence : chr "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" ...
## $ LinkerPrimerSequence : logi NA NA NA NA NA NA ...
## $ Location : chr "E1" "E1" "E1" "E1" ...
## $ Project : chr "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "WL" "WL" "WL" "WL" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 5 5 5 5 5 5 5 3 3 3 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...ssuguildREADRICH <- ssuguildREADRICH %>% left_join(metassu, by = "ID")
str(ssuguildREADRICH)## 'data.frame': 183 obs. of 15 variables:
## $ ID : chr "MM933M" "MM933M" "MM933M" "MM936M" ...
## $ Functional.Group : chr "Ancestral" "Edaphophilic" "Rhizophilic" "Ancestral" ...
## $ OTU_Richness_Sample : int 26 2 102 20 1 84 19 1 101 20 ...
## $ Read_Abundance_Sample: num 194 23.2 616.6 103.9 12.2 ...
## $ #SampleID : int 5 5 5 3 3 3 1 1 1 2 ...
## $ BarcodeSequence : chr "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "GCAACCTTAGAGTGTG" ...
## $ LinkerPrimerSequence : logi NA NA NA NA NA NA ...
## $ Location : chr "E1" "E1" "E1" "C1" ...
## $ Project : chr "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "WL" "WL" "WL" "WL" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 5 5 5 3 3 3 1 1 1 2 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...filename <- paste0(date, "_ssu_funcguild_read_rich.csv")
write.csv(ssuguildREADRICH, file = filename, row.names = FALSE)Adds meta data to data for box plot
### For ssufamily, moves ID to single line and family to the length
LssufamilyOTURICH <- dcast(ssufamily, ID ~ ssufamily, value.var = "Family_OTU_Richness",
fun.aggregate = sum)
LssufamilyOTUREAD <- dcast(ssufamily, ID ~ ssufamily, value.var = "Family_Read_Abundance",
fun.aggregate = sum)
head(LssufamilyOTUREAD)## ID Acaulosporaceae Ambisporaceae Archaeosporaceae
## 1 MM933M 12.4250 23.9417 157.5890
## 2 MM9362CU 15.4041 33.2473 63.2683
## 3 MM936M 10.1868 25.9602 67.7783
## 4 MM941M 33.0782 14.6280 47.8637
## 5 MM943M 44.8828 23.4654 92.1497
## 6 MMS700 22.8050 20.1172 90.3423
## Claroideoglomeraceae Diversisporaceae Glomeraceae Paraglomeraceae
## 1 99.6242 23.1980 387.3753 129.6337
## 2 118.7108 18.3924 464.2526 172.7572
## 3 91.5404 12.2380 406.1616 81.1362
## 4 93.6373 13.9450 503.8588 81.3654
## 5 74.7298 22.9824 317.0982 103.0184
## 6 106.1814 18.6176 429.0388 118.3515Functional group read and richness together
ssuguildREADRICHlong <- data.frame(ssul %>% group_by(ID) %>% summarise(OTU_Richness_Sample = length(unique(ssuotu[ssureads >
0])), Ancestral_Richness = length(unique(ssuotu[ssureads > 0 & Functional.Group ==
"Ancestral"])), Edaphophilic_Richness = length(unique(ssuotu[ssureads >
0 & Functional.Group == "Edaphophilic"])), Rhizophilic_Richness = length(unique(ssuotu[ssureads >
0 & Functional.Group == "Rhizophilic"])), Ancestral_Read_Abundance = sum(ssureads[Functional.Group ==
"Ancestral"]), Edaphophilic_Read_Abundance = sum(ssureads[Functional.Group ==
"Edaphophilic"]), Rhizophilic_Read_Abundance = sum(ssureads[Functional.Group ==
"Rhizophilic"]), Read_Abundance_Sample = sum(ssureads)))
View(ssuguildREADRICHlong)
ssuguildREADRICHlong <- merge(ssuguildREADRICHlong, metassu, by = "ID")
filename <- paste0(date, "_ssu_funcguild_read_rich.csv")
write.csv(ssuguildREADRICHlong, file = filename, row.names = FALSE)
str(ssuguildREADRICHlong)## 'data.frame': 61 obs. of 20 variables:
## $ ID : Factor w/ 61 levels "MM933M","MM936M",..: 1 5 2 3 4 6 7 8 9 10 ...
## $ OTU_Richness_Sample : int 130 167 105 121 99 151 158 161 138 148 ...
## $ Ancestral_Richness : int 26 28 20 19 20 23 22 24 19 23 ...
## $ Edaphophilic_Richness : int 2 3 1 1 3 3 2 4 2 1 ...
## $ Rhizophilic_Richness : int 102 136 84 101 76 125 134 133 117 124 ...
## $ Ancestral_Read_Abundance : num 194 111.9 103.9 95.6 160.5 ...
## $ Edaphophilic_Read_Abundance: num 23.2 18.4 12.2 13.9 23 ...
## $ Rhizophilic_Read_Abundance : num 617 756 579 679 495 ...
## $ Read_Abundance_Sample : num 834 886 695 788 678 ...
## $ #SampleID : int 5 4 3 1 2 47 48 49 50 51 ...
## $ BarcodeSequence : chr "ATCCCGTATCGATTGG" "GCAACCTTTCGATTGG" "GCAACCTTAGAGTGTG" "GCAACCTTTGGGTGAT" ...
## $ LinkerPrimerSequence : logi NA NA NA NA NA NA ...
## $ Location : chr "E1" "D1" "C1" "A1" ...
## $ Project : chr "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "WL" "WL" "WL" "WL" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 2 2 2 2 2 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 5 4 3 1 2 1 2 3 1 2 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...Functional group read and richness separate
ssuguildREAD <- ssul %>% group_by(ID) %>% summarise(Ancestral_Read_Abundance = sum(ssureads[Functional.Group ==
"Ancestral"]), Edaphophilic_Read_Abundance = sum(ssureads[Functional.Group ==
"Edaphophilic"]), Rhizophilic_Read_Abundance = sum(ssureads[Functional.Group ==
"Rhizophilic"]))
ssuguildRICH <- ssul %>% group_by(ID) %>% summarise(OTU_Richness_Sample = length(unique(ssuotu[ssureads >
0])), Ancestral_Richness = length(unique(ssuotu[ssureads > 0 & Functional.Group ==
"Ancestral"])), Edaphophilic_Richness = length(unique(ssuotu[ssureads >
0 & Functional.Group == "Edaphophilic"])), Rhizophilic_Richness = length(unique(ssuotu[ssureads >
0 & Functional.Group == "Rhizophilic"])))
head(ssuguildRICH)## # A tibble: 6 x 5
## ID OTU_Richness_Sa… Ancestral_Richn… Edaphophilic_Ri… Rhizophilic_Ri…
## <fct> <int> <int> <int> <int>
## 1 MM93… 130 26 2 102
## 2 MM93… 105 20 1 84
## 3 MM94… 121 19 1 101
## 4 MM94… 99 20 3 76
## 5 MM93… 167 28 3 136
## 6 MMS7… 151 23 3 125Add metadata to above #Everytime I run this the output is different, sometime it merges the data with the metadata with two copies of the metadata (e.g., Site.x and Site.y) #To address this, DON’T RUN NEXT CHUNK WITHOUT RERUNNING THE PREVIOUS CHUNK GENERATING THE ORIGINAL ssuguildRICH FILE, WITHOUT THE METADATA, tried to use “all.y = FALSE” after merge(ssuguildRICH, metassu, by=“ID”,“”
ssuguildRICH <- merge(ssuguildRICH, metassu, by = "ID", all.y = FALSE)
ssuguildREAD <- merge(ssuguildREAD, metassu, by = "ID", all.y = FALSE)
str(ssuguildRICH)## 'data.frame': 61 obs. of 16 variables:
## $ ID : Factor w/ 61 levels "MM933M","MM936M",..: 1 5 2 3 4 6 7 8 9 10 ...
## $ OTU_Richness_Sample : int 130 167 105 121 99 151 158 161 138 148 ...
## $ Ancestral_Richness : int 26 28 20 19 20 23 22 24 19 23 ...
## $ Edaphophilic_Richness: int 2 3 1 1 3 3 2 4 2 1 ...
## $ Rhizophilic_Richness : int 102 136 84 101 76 125 134 133 117 124 ...
## $ #SampleID : int 5 4 3 1 2 47 48 49 50 51 ...
## $ BarcodeSequence : chr "ATCCCGTATCGATTGG" "GCAACCTTTCGATTGG" "GCAACCTTAGAGTGTG" "GCAACCTTTGGGTGAT" ...
## $ LinkerPrimerSequence : logi NA NA NA NA NA NA ...
## $ Location : chr "E1" "D1" "C1" "A1" ...
## $ Project : chr "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "WL" "WL" "WL" "WL" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 2 2 2 2 2 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 5 4 3 1 2 1 2 3 1 2 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...str(ssuguildREAD)## 'data.frame': 61 obs. of 15 variables:
## $ ID : Factor w/ 61 levels "MM933M","MM936M",..: 1 5 2 3 4 6 7 8 9 10 ...
## $ Ancestral_Read_Abundance : num 194 111.9 103.9 95.6 160.5 ...
## $ Edaphophilic_Read_Abundance: num 23.2 18.4 12.2 13.9 23 ...
## $ Rhizophilic_Read_Abundance : num 617 756 579 679 495 ...
## $ #SampleID : int 5 4 3 1 2 47 48 49 50 51 ...
## $ BarcodeSequence : chr "ATCCCGTATCGATTGG" "GCAACCTTTCGATTGG" "GCAACCTTAGAGTGTG" "GCAACCTTTGGGTGAT" ...
## $ LinkerPrimerSequence : logi NA NA NA NA NA NA ...
## $ Location : chr "E1" "D1" "C1" "A1" ...
## $ Project : chr "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "WL" "WL" "WL" "WL" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 2 2 2 2 2 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 5 4 3 1 2 1 2 3 1 2 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...View(ssuguildREAD)Examining the statistic support for increased rhizophilic richness in recipientpost plots
y <- ssuguildRICH$Rhizophilic_Richness
str(y)## int [1:61] 102 136 84 101 76 125 134 133 117 124 ...D <- ssuguildRICH$Description
S <- ssuguildRICH$Site
length(S)## [1] 61length(D)## [1] 61`?`(lm)
# lm(y~D)
SRhizRich <- lm(y ~ S)
str(SRhizRich)## List of 13
## $ coefficients : Named num [1:4] 125.4 12.56 6.18 -2.98
## ..- attr(*, "names")= chr [1:4] "(Intercept)" "SHH" "SPS" "SWL"
## $ residuals : Named num [1:61] -20.4 13.6 -38.4 -21.4 -46.4 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ effects : Named num [1:61] -1030.1 39.5 24.8 5.6 -39.2 ...
## ..- attr(*, "names")= chr [1:61] "(Intercept)" "SHH" "SPS" "SWL" ...
## $ rank : int 4
## $ fitted.values: Named num [1:61] 122 122 122 122 122 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ assign : int [1:4] 0 1 1 1
## $ qr :List of 5
## ..$ qr : num [1:61, 1:4] -7.81 0.128 0.128 0.128 0.128 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:61] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:4] "(Intercept)" "SHH" "SPS" "SWL"
## .. ..- attr(*, "assign")= int [1:4] 0 1 1 1
## .. ..- attr(*, "contrasts")=List of 1
## .. .. ..$ S: chr "contr.treatment"
## ..$ qraux: num [1:4] 1.13 1.09 1.14 1.11
## ..$ pivot: int [1:4] 1 2 3 4
## ..$ tol : num 1e-07
## ..$ rank : int 4
## ..- attr(*, "class")= chr "qr"
## $ df.residual : int 57
## $ contrasts :List of 1
## ..$ S: chr "contr.treatment"
## $ xlevels :List of 1
## ..$ S: chr [1:4] "DS" "HH" "PS" "WL"
## $ call : language lm(formula = y ~ S)
## $ terms :Classes 'terms', 'formula' language y ~ S
## .. ..- attr(*, "variables")= language list(y, S)
## .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:2] "y" "S"
## .. .. .. ..$ : chr "S"
## .. ..- attr(*, "term.labels")= chr "S"
## .. ..- attr(*, "order")= int 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(y, S)
## .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
## .. .. ..- attr(*, "names")= chr [1:2] "y" "S"
## $ model :'data.frame': 61 obs. of 2 variables:
## ..$ y: int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
## ..$ S: chr [1:61] "WL" "WL" "WL" "WL" ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' language y ~ S
## .. .. ..- attr(*, "variables")= language list(y, S)
## .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:2] "y" "S"
## .. .. .. .. ..$ : chr "S"
## .. .. ..- attr(*, "term.labels")= chr "S"
## .. .. ..- attr(*, "order")= int 1
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(y, S)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
## .. .. .. ..- attr(*, "names")= chr [1:2] "y" "S"
## - attr(*, "class")= chr "lm"DRhizRich <- lm(y ~ D)
str(DRhizRich)## List of 13
## $ coefficients : Named num [1:3] 125.4 11.44 -8.98
## ..- attr(*, "names")= chr [1:3] "(Intercept)" "DRecipientPost" "DRecipientPre"
## $ residuals : Named num [1:61] -14.4 19.6 -32.4 -15.4 -40.4 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ effects : Named num [1:61] -1030.1 62.3 16.9 -13.5 -38.5 ...
## ..- attr(*, "names")= chr [1:61] "(Intercept)" "DRecipientPost" "DRecipientPre" "" ...
## $ rank : int 3
## $ fitted.values: Named num [1:61] 116 116 116 116 116 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ assign : int [1:3] 0 1 1
## $ qr :List of 5
## ..$ qr : num [1:61, 1:3] -7.81 0.128 0.128 0.128 0.128 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:61] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:3] "(Intercept)" "DRecipientPost" "DRecipientPre"
## .. ..- attr(*, "assign")= int [1:3] 0 1 1
## .. ..- attr(*, "contrasts")=List of 1
## .. .. ..$ D: chr "contr.treatment"
## ..$ qraux: num [1:3] 1.13 1.18 1.12
## ..$ pivot: int [1:3] 1 2 3
## ..$ tol : num 1e-07
## ..$ rank : int 3
## ..- attr(*, "class")= chr "qr"
## $ df.residual : int 58
## $ contrasts :List of 1
## ..$ D: chr "contr.treatment"
## $ xlevels :List of 1
## ..$ D: chr [1:3] "Donor" "RecipientPost" "RecipientPre"
## $ call : language lm(formula = y ~ D)
## $ terms :Classes 'terms', 'formula' language y ~ D
## .. ..- attr(*, "variables")= language list(y, D)
## .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:2] "y" "D"
## .. .. .. ..$ : chr "D"
## .. ..- attr(*, "term.labels")= chr "D"
## .. ..- attr(*, "order")= int 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(y, D)
## .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
## .. .. ..- attr(*, "names")= chr [1:2] "y" "D"
## $ model :'data.frame': 61 obs. of 2 variables:
## ..$ y: int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
## ..$ D: chr [1:61] "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' language y ~ D
## .. .. ..- attr(*, "variables")= language list(y, D)
## .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:2] "y" "D"
## .. .. .. .. ..$ : chr "D"
## .. .. ..- attr(*, "term.labels")= chr "D"
## .. .. ..- attr(*, "order")= int 1
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(y, D)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
## .. .. .. ..- attr(*, "names")= chr [1:2] "y" "D"
## - attr(*, "class")= chr "lm"DSRhizRich <- lm(y ~ S * D)
str(DSRhizRich)## List of 13
## $ coefficients : Named num [1:12] 125.4 8.6 -3.4 -18.6 31.2 ...
## ..- attr(*, "names")= chr [1:12] "(Intercept)" "SHH" "SPS" "SWL" ...
## $ residuals : Named num [1:61] -4.83 29.17 -22.83 -5.83 -30.83 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ effects : Named num [1:61] -1030.1 39.5 24.8 5.6 49.4 ...
## ..- attr(*, "names")= chr [1:61] "(Intercept)" "SHH" "SPS" "SWL" ...
## $ rank : int 7
## $ fitted.values: Named num [1:61] 107 107 107 107 107 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ assign : int [1:12] 0 1 1 1 2 2 3 3 3 3 ...
## $ qr :List of 5
## ..$ qr : num [1:61, 1:12] -7.81 0.128 0.128 0.128 0.128 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:61] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:12] "(Intercept)" "SHH" "SPS" "SWL" ...
## .. ..- attr(*, "assign")= int [1:12] 0 1 1 1 2 2 3 3 3 3 ...
## .. ..- attr(*, "contrasts")=List of 2
## .. .. ..$ S: chr "contr.treatment"
## .. .. ..$ D: chr "contr.treatment"
## ..$ qraux: num [1:12] 1.13 1.09 1.14 1.11 1.11 ...
## ..$ pivot: int [1:12] 1 2 3 4 5 7 8 6 9 10 ...
## ..$ tol : num 1e-07
## ..$ rank : int 7
## ..- attr(*, "class")= chr "qr"
## $ df.residual : int 54
## $ contrasts :List of 2
## ..$ S: chr "contr.treatment"
## ..$ D: chr "contr.treatment"
## $ xlevels :List of 2
## ..$ S: chr [1:4] "DS" "HH" "PS" "WL"
## ..$ D: chr [1:3] "Donor" "RecipientPost" "RecipientPre"
## $ call : language lm(formula = y ~ S * D)
## $ terms :Classes 'terms', 'formula' language y ~ S * D
## .. ..- attr(*, "variables")= language list(y, S, D)
## .. ..- attr(*, "factors")= int [1:3, 1:3] 0 1 0 0 0 1 0 1 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:3] "y" "S" "D"
## .. .. .. ..$ : chr [1:3] "S" "D" "S:D"
## .. ..- attr(*, "term.labels")= chr [1:3] "S" "D" "S:D"
## .. ..- attr(*, "order")= int [1:3] 1 1 2
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(y, S, D)
## .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "character" "character"
## .. .. ..- attr(*, "names")= chr [1:3] "y" "S" "D"
## $ model :'data.frame': 61 obs. of 3 variables:
## ..$ y: int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
## ..$ S: chr [1:61] "WL" "WL" "WL" "WL" ...
## ..$ D: chr [1:61] "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' language y ~ S * D
## .. .. ..- attr(*, "variables")= language list(y, S, D)
## .. .. ..- attr(*, "factors")= int [1:3, 1:3] 0 1 0 0 0 1 0 1 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:3] "y" "S" "D"
## .. .. .. .. ..$ : chr [1:3] "S" "D" "S:D"
## .. .. ..- attr(*, "term.labels")= chr [1:3] "S" "D" "S:D"
## .. .. ..- attr(*, "order")= int [1:3] 1 1 2
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(y, S, D)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "character" "character"
## .. .. .. ..- attr(*, "names")= chr [1:3] "y" "S" "D"
## - attr(*, "class")= chr "lm"summary(DSRhizRich)##
## Call:
## lm(formula = y ~ S * D)
##
## Residuals:
## Min 1Q Median 3Q Max
## -42.304 -7.304 2.000 7.000 35.167
##
## Coefficients: (5 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 125.400 6.942 18.065 < 2e-16 ***
## SHH 8.600 12.987 0.662 0.51065
## SPS -3.400 10.413 -0.327 0.74529
## SWL -18.567 9.399 -1.975 0.05335 .
## DRecipientPost 31.167 8.962 3.478 0.00101 **
## DRecipientPre NA NA NA NA
## SHH:DRecipientPost -26.862 14.535 -1.848 0.07006 .
## SPS:DRecipientPost -19.033 12.514 -1.521 0.13411
## SWL:DRecipientPost NA NA NA NA
## SHH:DRecipientPre NA NA NA NA
## SPS:DRecipientPre NA NA NA NA
## SWL:DRecipientPre NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.52 on 54 degrees of freedom
## Multiple R-squared: 0.3018, Adjusted R-squared: 0.2242
## F-statistic: 3.89 on 6 and 54 DF, p-value: 0.002672y <- ssuguildRICH$Rhizophilic_Richness
str(y)## int [1:61] 102 136 84 101 76 125 134 133 117 124 ...D <- ssuguildRICH$Description
S <- ssuguildRICH$Site
length(S)## [1] 61length(D)## [1] 61`?`(lm)
glm(y ~ D)##
## Call: glm(formula = y ~ D)
##
## Coefficients:
## (Intercept) DRecipientPost DRecipientPre
## 125.400 11.441 -8.983
##
## Degrees of Freedom: 60 Total (i.e. Null); 58 Residual
## Null Deviance: 18630
## Residual Deviance: 14470 AIC: 514.7G_SRhizRich <- glm(y ~ S)
str(G_SRhizRich)## List of 30
## $ coefficients : Named num [1:4] 125.4 12.56 6.18 -2.98
## ..- attr(*, "names")= chr [1:4] "(Intercept)" "SHH" "SPS" "SWL"
## $ residuals : Named num [1:61] -20.4 13.6 -38.4 -21.4 -46.4 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ fitted.values : Named num [1:61] 122 122 122 122 122 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ effects : Named num [1:61] -1030.1 39.5 24.8 5.6 -39.2 ...
## ..- attr(*, "names")= chr [1:61] "(Intercept)" "SHH" "SPS" "SWL" ...
## $ R : num [1:4, 1:4] -7.81 0 0 0 -3.2 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:4] "(Intercept)" "SHH" "SPS" "SWL"
## .. ..$ : chr [1:4] "(Intercept)" "SHH" "SPS" "SWL"
## $ rank : int 4
## $ qr :List of 5
## ..$ qr : num [1:61, 1:4] -7.81 0.128 0.128 0.128 0.128 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:61] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:4] "(Intercept)" "SHH" "SPS" "SWL"
## ..$ rank : int 4
## ..$ qraux: num [1:4] 1.13 1.09 1.14 1.11
## ..$ pivot: int [1:4] 1 2 3 4
## ..$ tol : num 1e-11
## ..- attr(*, "class")= chr "qr"
## $ family :List of 11
## ..$ family : chr "gaussian"
## ..$ link : chr "identity"
## ..$ linkfun :function (mu)
## ..$ linkinv :function (eta)
## ..$ variance :function (mu)
## ..$ dev.resids:function (y, mu, wt)
## ..$ aic :function (y, n, mu, wt, dev)
## ..$ mu.eta :function (eta)
## ..$ initialize: expression({ n <- rep.int(1, nobs) if (is.null(etastart) && is.null(start) && is.null(mustart) && ((family$lin| __truncated__
## ..$ validmu :function (mu)
## ..$ valideta :function (eta)
## ..- attr(*, "class")= chr "family"
## $ linear.predictors: Named num [1:61] 122 122 122 122 122 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ deviance : num 16424
## $ aic : num 524
## $ null.deviance : num 18634
## $ iter : int 2
## $ weights : Named num [1:61] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ prior.weights : Named num [1:61] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ df.residual : int 57
## $ df.null : int 60
## $ y : Named int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ converged : logi TRUE
## $ boundary : logi FALSE
## $ model :'data.frame': 61 obs. of 2 variables:
## ..$ y: int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
## ..$ S: chr [1:61] "WL" "WL" "WL" "WL" ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' language y ~ S
## .. .. ..- attr(*, "variables")= language list(y, S)
## .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:2] "y" "S"
## .. .. .. .. ..$ : chr "S"
## .. .. ..- attr(*, "term.labels")= chr "S"
## .. .. ..- attr(*, "order")= int 1
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(y, S)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
## .. .. .. ..- attr(*, "names")= chr [1:2] "y" "S"
## $ call : language glm(formula = y ~ S)
## $ formula :Class 'formula' language y ~ S
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## $ terms :Classes 'terms', 'formula' language y ~ S
## .. ..- attr(*, "variables")= language list(y, S)
## .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:2] "y" "S"
## .. .. .. ..$ : chr "S"
## .. ..- attr(*, "term.labels")= chr "S"
## .. ..- attr(*, "order")= int 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(y, S)
## .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
## .. .. ..- attr(*, "names")= chr [1:2] "y" "S"
## $ data :<environment: R_GlobalEnv>
## $ offset : NULL
## $ control :List of 3
## ..$ epsilon: num 1e-08
## ..$ maxit : num 25
## ..$ trace : logi FALSE
## $ method : chr "glm.fit"
## $ contrasts :List of 1
## ..$ S: chr "contr.treatment"
## $ xlevels :List of 1
## ..$ S: chr [1:4] "DS" "HH" "PS" "WL"
## - attr(*, "class")= chr [1:2] "glm" "lm"G_DRhizRich <- glm(y ~ D)
str(G_DRhizRich)## List of 30
## $ coefficients : Named num [1:3] 125.4 11.44 -8.98
## ..- attr(*, "names")= chr [1:3] "(Intercept)" "DRecipientPost" "DRecipientPre"
## $ residuals : Named num [1:61] -14.4 19.6 -32.4 -15.4 -40.4 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ fitted.values : Named num [1:61] 116 116 116 116 116 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ effects : Named num [1:61] -1030.1 62.3 16.9 -13.5 -38.5 ...
## ..- attr(*, "names")= chr [1:61] "(Intercept)" "DRecipientPost" "DRecipientPre" "" ...
## $ R : num [1:3, 1:3] -7.81 0 0 -5.63 3.5 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:3] "(Intercept)" "DRecipientPost" "DRecipientPre"
## .. ..$ : chr [1:3] "(Intercept)" "DRecipientPost" "DRecipientPre"
## $ rank : int 3
## $ qr :List of 5
## ..$ qr : num [1:61, 1:3] -7.81 0.128 0.128 0.128 0.128 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:61] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:3] "(Intercept)" "DRecipientPost" "DRecipientPre"
## ..$ rank : int 3
## ..$ qraux: num [1:3] 1.13 1.18 1.12
## ..$ pivot: int [1:3] 1 2 3
## ..$ tol : num 1e-11
## ..- attr(*, "class")= chr "qr"
## $ family :List of 11
## ..$ family : chr "gaussian"
## ..$ link : chr "identity"
## ..$ linkfun :function (mu)
## ..$ linkinv :function (eta)
## ..$ variance :function (mu)
## ..$ dev.resids:function (y, mu, wt)
## ..$ aic :function (y, n, mu, wt, dev)
## ..$ mu.eta :function (eta)
## ..$ initialize: expression({ n <- rep.int(1, nobs) if (is.null(etastart) && is.null(start) && is.null(mustart) && ((family$lin| __truncated__
## ..$ validmu :function (mu)
## ..$ valideta :function (eta)
## ..- attr(*, "class")= chr "family"
## $ linear.predictors: Named num [1:61] 116 116 116 116 116 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ deviance : num 14472
## $ aic : num 515
## $ null.deviance : num 18634
## $ iter : int 2
## $ weights : Named num [1:61] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ prior.weights : Named num [1:61] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ df.residual : int 58
## $ df.null : int 60
## $ y : Named int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ converged : logi TRUE
## $ boundary : logi FALSE
## $ model :'data.frame': 61 obs. of 2 variables:
## ..$ y: int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
## ..$ D: chr [1:61] "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' language y ~ D
## .. .. ..- attr(*, "variables")= language list(y, D)
## .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:2] "y" "D"
## .. .. .. .. ..$ : chr "D"
## .. .. ..- attr(*, "term.labels")= chr "D"
## .. .. ..- attr(*, "order")= int 1
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(y, D)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
## .. .. .. ..- attr(*, "names")= chr [1:2] "y" "D"
## $ call : language glm(formula = y ~ D)
## $ formula :Class 'formula' language y ~ D
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## $ terms :Classes 'terms', 'formula' language y ~ D
## .. ..- attr(*, "variables")= language list(y, D)
## .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:2] "y" "D"
## .. .. .. ..$ : chr "D"
## .. ..- attr(*, "term.labels")= chr "D"
## .. ..- attr(*, "order")= int 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(y, D)
## .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
## .. .. ..- attr(*, "names")= chr [1:2] "y" "D"
## $ data :<environment: R_GlobalEnv>
## $ offset : NULL
## $ control :List of 3
## ..$ epsilon: num 1e-08
## ..$ maxit : num 25
## ..$ trace : logi FALSE
## $ method : chr "glm.fit"
## $ contrasts :List of 1
## ..$ D: chr "contr.treatment"
## $ xlevels :List of 1
## ..$ D: chr [1:3] "Donor" "RecipientPost" "RecipientPre"
## - attr(*, "class")= chr [1:2] "glm" "lm"summary(G_DRhizRich)##
## Call:
## glm(formula = y ~ D)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -40.841 -8.400 3.159 10.159 26.159
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 125.400 7.064 17.751 <2e-16 ***
## DRecipientPost 11.441 7.455 1.535 0.13
## DRecipientPre -8.983 8.408 -1.068 0.29
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 249.5173)
##
## Null deviance: 18634 on 60 degrees of freedom
## Residual deviance: 14472 on 58 degrees of freedom
## AIC: 514.73
##
## Number of Fisher Scoring iterations: 2DSRhizRich <- lm(y ~ S * D)
str(DSRhizRich)## List of 13
## $ coefficients : Named num [1:12] 125.4 8.6 -3.4 -18.6 31.2 ...
## ..- attr(*, "names")= chr [1:12] "(Intercept)" "SHH" "SPS" "SWL" ...
## $ residuals : Named num [1:61] -4.83 29.17 -22.83 -5.83 -30.83 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ effects : Named num [1:61] -1030.1 39.5 24.8 5.6 49.4 ...
## ..- attr(*, "names")= chr [1:61] "(Intercept)" "SHH" "SPS" "SWL" ...
## $ rank : int 7
## $ fitted.values: Named num [1:61] 107 107 107 107 107 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ assign : int [1:12] 0 1 1 1 2 2 3 3 3 3 ...
## $ qr :List of 5
## ..$ qr : num [1:61, 1:12] -7.81 0.128 0.128 0.128 0.128 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:61] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:12] "(Intercept)" "SHH" "SPS" "SWL" ...
## .. ..- attr(*, "assign")= int [1:12] 0 1 1 1 2 2 3 3 3 3 ...
## .. ..- attr(*, "contrasts")=List of 2
## .. .. ..$ S: chr "contr.treatment"
## .. .. ..$ D: chr "contr.treatment"
## ..$ qraux: num [1:12] 1.13 1.09 1.14 1.11 1.11 ...
## ..$ pivot: int [1:12] 1 2 3 4 5 7 8 6 9 10 ...
## ..$ tol : num 1e-07
## ..$ rank : int 7
## ..- attr(*, "class")= chr "qr"
## $ df.residual : int 54
## $ contrasts :List of 2
## ..$ S: chr "contr.treatment"
## ..$ D: chr "contr.treatment"
## $ xlevels :List of 2
## ..$ S: chr [1:4] "DS" "HH" "PS" "WL"
## ..$ D: chr [1:3] "Donor" "RecipientPost" "RecipientPre"
## $ call : language lm(formula = y ~ S * D)
## $ terms :Classes 'terms', 'formula' language y ~ S * D
## .. ..- attr(*, "variables")= language list(y, S, D)
## .. ..- attr(*, "factors")= int [1:3, 1:3] 0 1 0 0 0 1 0 1 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:3] "y" "S" "D"
## .. .. .. ..$ : chr [1:3] "S" "D" "S:D"
## .. ..- attr(*, "term.labels")= chr [1:3] "S" "D" "S:D"
## .. ..- attr(*, "order")= int [1:3] 1 1 2
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(y, S, D)
## .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "character" "character"
## .. .. ..- attr(*, "names")= chr [1:3] "y" "S" "D"
## $ model :'data.frame': 61 obs. of 3 variables:
## ..$ y: int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
## ..$ S: chr [1:61] "WL" "WL" "WL" "WL" ...
## ..$ D: chr [1:61] "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' language y ~ S * D
## .. .. ..- attr(*, "variables")= language list(y, S, D)
## .. .. ..- attr(*, "factors")= int [1:3, 1:3] 0 1 0 0 0 1 0 1 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:3] "y" "S" "D"
## .. .. .. .. ..$ : chr [1:3] "S" "D" "S:D"
## .. .. ..- attr(*, "term.labels")= chr [1:3] "S" "D" "S:D"
## .. .. ..- attr(*, "order")= int [1:3] 1 1 2
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(y, S, D)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "character" "character"
## .. .. .. ..- attr(*, "names")= chr [1:3] "y" "S" "D"
## - attr(*, "class")= chr "lm"summary(DSRhizRich)##
## Call:
## lm(formula = y ~ S * D)
##
## Residuals:
## Min 1Q Median 3Q Max
## -42.304 -7.304 2.000 7.000 35.167
##
## Coefficients: (5 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 125.400 6.942 18.065 < 2e-16 ***
## SHH 8.600 12.987 0.662 0.51065
## SPS -3.400 10.413 -0.327 0.74529
## SWL -18.567 9.399 -1.975 0.05335 .
## DRecipientPost 31.167 8.962 3.478 0.00101 **
## DRecipientPre NA NA NA NA
## SHH:DRecipientPost -26.862 14.535 -1.848 0.07006 .
## SPS:DRecipientPost -19.033 12.514 -1.521 0.13411
## SWL:DRecipientPost NA NA NA NA
## SHH:DRecipientPre NA NA NA NA
## SPS:DRecipientPre NA NA NA NA
## SWL:DRecipientPre NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.52 on 54 degrees of freedom
## Multiple R-squared: 0.3018, Adjusted R-squared: 0.2242
## F-statistic: 3.89 on 6 and 54 DF, p-value: 0.002672lm(y ~ D * S)##
## Call:
## lm(formula = y ~ D * S)
##
## Coefficients:
## (Intercept) DRecipientPost DRecipientPre
## 125.40 12.60 -18.57
## SHH SPS SWL
## 27.17 15.17 NA
## DRecipientPost:SHH DRecipientPre:SHH DRecipientPost:SPS
## -26.86 NA -19.03
## DRecipientPre:SPS DRecipientPost:SWL DRecipientPre:SWL
## NA NA NAlm(y ~ D + S)##
## Call:
## lm(formula = y ~ D + S)
##
## Coefficients:
## (Intercept) DRecipientPost DRecipientPre SHH
## 125.400 5.752 -11.719 8.205
## SPS SWL
## 4.105 NA`?`(glm)
glm(y ~ D)##
## Call: glm(formula = y ~ D)
##
## Coefficients:
## (Intercept) DRecipientPost DRecipientPre
## 125.400 11.441 -8.983
##
## Degrees of Freedom: 60 Total (i.e. Null); 58 Residual
## Null Deviance: 18630
## Residual Deviance: 14470 AIC: 514.7# Description+Sitestart ‘2017_16_07_SSU_Boxplots.R’ Plotting
# richness
plot.df <- ssuguildREADRICH
str(plot.df)## 'data.frame': 183 obs. of 15 variables:
## $ ID : chr "MM933M" "MM933M" "MM933M" "MM936M" ...
## $ Functional.Group : chr "Ancestral" "Edaphophilic" "Rhizophilic" "Ancestral" ...
## $ OTU_Richness_Sample : int 26 2 102 20 1 84 19 1 101 20 ...
## $ Read_Abundance_Sample: num 194 23.2 616.6 103.9 12.2 ...
## $ #SampleID : int 5 5 5 3 3 3 1 1 1 2 ...
## $ BarcodeSequence : chr "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "GCAACCTTAGAGTGTG" ...
## $ LinkerPrimerSequence : logi NA NA NA NA NA NA ...
## $ Location : chr "E1" "E1" "E1" "C1" ...
## $ Project : chr "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "WL" "WL" "WL" "WL" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 5 5 5 3 3 3 1 1 1 2 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...richbox <- ggplot(plot.df, aes(x = Functional.Group, y = OTU_Richness_Sample)) +
geom_boxplot(aes(fill = Description))
richbox
richbox <- richbox + theme_bw(base_size = 15) + xlab("Functional Group") + ylab("AMF Taxa Richness")
richbox
richbox <- richbox + scale_fill_manual(values = c("red", "darkslateblue", "gold3"))
richbox
## reads
plot.read <- ssuguildREADRICH
str(plot.read)
reachbox <- ggplot(plot.read, aes(x = Functional.Group, y = OTU_Richness_Sample)) +
geom_boxplot(aes(fill = Description))
reachbox
reachbox <- reachbox + theme_bw(base_size = 15) + xlab("Functional Group") +
ylab("AMF Taxa Reads")
reachbox
reachbox <- reachbox + scale_fill_manual(values = c("red", "darkslateblue",
"gold3"))
reachboxADAPTED FROM #GLMM for roots and soil OTU Richness SSU by functional group-modified from “GLM_NAU_6_2017” #—different models for each functional group #Author: Michala Phillips #MODIFY FOR THIS DATA
library(car)
require(MASS)
library(lme4)Choose distribution based on data
find the right distribution for data
ssuguildREADRICHlong$OTU_Richness_Sample.t <- ssuguildREADRICHlong$OTU_Richness_Sample +
1
qqp(ssuguildREADRICHlong$OTU_Richness_Sample.t, "norm")
## [1] 5 3# LOGNORMAL
qqp(ssuguildREADRICHlong$OTU_Richness_Sample, "lnorm")
## [1] 5 3# NEG BINOMIAL
nbinom <- fitdistr(ssuguildREADRICHlong$OTU_Richness_Sample.t, "Negative Binomial")
qqp(ssuguildREADRICHlong$OTU_Richness_Sample.t, "nbinom", size = nbinom$estimate[[1]],
mu = nbinom$estimate[[2]])
## [1] 5 3# POISSON
poisson <- fitdistr(ssuguildREADRICHlong$OTU_Richness_Sample.t, "Poisson")
qqp(ssuguildREADRICHlong$OTU_Richness_Sample.t, "pois", lambda = poisson$estimate)
## [1] 5 3# GAMMA
gamma <- fitdistr(ssuguildREADRICHlong$OTU_Richness_Sample.t, "gamma")
qqp(ssuguildREADRICHlong$OTU_Richness_Sample.t, "gamma", shape = gamma$estimate[[1]],
rate = gamma$estimate[[2]])
## [1] 5 3# modelingAMF <- glm(OTU_Richness_Sample.t ~ Site + Treat, data = ssuguildREADRICHlong,
family = gaussian(link = "identity"))
modelingAMF##
## Call: glm(formula = OTU_Richness_Sample.t ~ Site + Treat, family = gaussian(link = "identity"),
## data = ssuguildREADRICHlong)
##
## Coefficients:
## (Intercept) SiteHH SitePS SiteWL TreatD
## 170.6650 -1.2149 -8.8985 -10.4627 -9.8063
## TreatF TreatO TreatS TreatT
## 14.1937 -20.2650 -0.4437 1.2222
##
## Degrees of Freedom: 60 Total (i.e. Null); 52 Residual
## Null Deviance: 26340
## Residual Deviance: 16670 AIC: 535.4stepAIC(modelingAMF)## Start: AIC=535.36
## OTU_Richness_Sample.t ~ Site + Treat
##
## Df Deviance AIC
## - Site 3 17693 532.98
## <none> 16673 535.36
## - Treat 5 23069 545.17
##
## Step: AIC=532.98
## OTU_Richness_Sample.t ~ Treat
##
## Df Deviance AIC
## <none> 17693 532.98
## - Treat 5 26339 547.25##
## Call: glm(formula = OTU_Richness_Sample.t ~ Treat, family = gaussian(link = "identity"),
## data = ssuguildREADRICHlong)
##
## Coefficients:
## (Intercept) TreatD TreatF TreatO TreatS
## 166.8889 -12.8889 11.1111 -22.4183 -0.7639
## TreatT
## 1.2222
##
## Degrees of Freedom: 60 Total (i.e. Null); 55 Residual
## Null Deviance: 26340
## Residual Deviance: 17690 AIC: 533modelingAMFfinal <- glm(OTU_Richness_Sample.t ~ Site + Treat, family = gaussian(link = "identity"),
data = ssuguildREADRICHlong)
summary(modelingAMFfinal)##
## Call:
## glm(formula = OTU_Richness_Sample.t ~ Site + Treat, family = gaussian(link = "identity"),
## data = ssuguildREADRICHlong)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -41.644 -8.501 3.328 9.604 33.063
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 170.6650 11.8523 14.399 <2e-16 ***
## SiteHH -1.2149 10.6902 -0.114 0.9100
## SitePS -8.8985 10.3362 -0.861 0.3932
## SiteWL -10.4627 10.1178 -1.034 0.3059
## TreatD -9.8063 8.8021 -1.114 0.2704
## TreatF 14.1937 8.8021 1.613 0.1129
## TreatO -20.2650 8.7378 -2.319 0.0243 *
## TreatS -0.4437 8.7041 -0.051 0.9595
## TreatT 1.2222 8.4412 0.145 0.8854
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 320.6386)
##
## Null deviance: 26339 on 60 degrees of freedom
## Residual deviance: 16673 on 52 degrees of freedom
## AIC: 535.36
##
## Number of Fisher Scoring iterations: 2plot(modelingAMFfinal)
modelingAMF <- glm(OTU_Richness_Sample.t ~ Site + Treat, data = ssuguildREADRICHlong,
family = gaussian(link = "identity"))
modelingAMF##
## Call: glm(formula = OTU_Richness_Sample.t ~ Site + Treat, family = gaussian(link = "identity"),
## data = ssuguildREADRICHlong)
##
## Coefficients:
## (Intercept) SiteHH SitePS SiteWL TreatD
## 170.6650 -1.2149 -8.8985 -10.4627 -9.8063
## TreatF TreatO TreatS TreatT
## 14.1937 -20.2650 -0.4437 1.2222
##
## Degrees of Freedom: 60 Total (i.e. Null); 52 Residual
## Null Deviance: 26340
## Residual Deviance: 16670 AIC: 535.4stepAIC(modelingAMF)## Start: AIC=535.36
## OTU_Richness_Sample.t ~ Site + Treat
##
## Df Deviance AIC
## - Site 3 17693 532.98
## <none> 16673 535.36
## - Treat 5 23069 545.17
##
## Step: AIC=532.98
## OTU_Richness_Sample.t ~ Treat
##
## Df Deviance AIC
## <none> 17693 532.98
## - Treat 5 26339 547.25##
## Call: glm(formula = OTU_Richness_Sample.t ~ Treat, family = gaussian(link = "identity"),
## data = ssuguildREADRICHlong)
##
## Coefficients:
## (Intercept) TreatD TreatF TreatO TreatS
## 166.8889 -12.8889 11.1111 -22.4183 -0.7639
## TreatT
## 1.2222
##
## Degrees of Freedom: 60 Total (i.e. Null); 55 Residual
## Null Deviance: 26340
## Residual Deviance: 17690 AIC: 533modelingAMFfinal <- glm(OTU_Richness_Sample.t ~ Site + Treat, family = gaussian(link = "identity"),
data = ssuguildREADRICHlong)
summary(modelingAMFfinal)##
## Call:
## glm(formula = OTU_Richness_Sample.t ~ Site + Treat, family = gaussian(link = "identity"),
## data = ssuguildREADRICHlong)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -41.644 -8.501 3.328 9.604 33.063
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 170.6650 11.8523 14.399 <2e-16 ***
## SiteHH -1.2149 10.6902 -0.114 0.9100
## SitePS -8.8985 10.3362 -0.861 0.3932
## SiteWL -10.4627 10.1178 -1.034 0.3059
## TreatD -9.8063 8.8021 -1.114 0.2704
## TreatF 14.1937 8.8021 1.613 0.1129
## TreatO -20.2650 8.7378 -2.319 0.0243 *
## TreatS -0.4437 8.7041 -0.051 0.9595
## TreatT 1.2222 8.4412 0.145 0.8854
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 320.6386)
##
## Null deviance: 26339 on 60 degrees of freedom
## Residual deviance: 16673 on 52 degrees of freedom
## AIC: 535.36
##
## Number of Fisher Scoring iterations: 2plot(modelingAMFfinal)
``` #Working on examining other variables
modelingAMF <- glm(OTU_Richness_Sample.t ~ Treat + Site, data = ssuguildREADRICHlong,
family = gaussian(link = "identity"))plot.df <- ssuguildREADRICH
str(plot.df)## 'data.frame': 183 obs. of 15 variables:
## $ ID : chr "MM933M" "MM933M" "MM933M" "MM936M" ...
## $ Functional.Group : chr "Ancestral" "Edaphophilic" "Rhizophilic" "Ancestral" ...
## $ OTU_Richness_Sample : int 26 2 102 20 1 84 19 1 101 20 ...
## $ Read_Abundance_Sample: num 194 23.2 616.6 103.9 12.2 ...
## $ #SampleID : int 5 5 5 3 3 3 1 1 1 2 ...
## $ BarcodeSequence : chr "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "GCAACCTTAGAGTGTG" ...
## $ LinkerPrimerSequence : logi NA NA NA NA NA NA ...
## $ Location : chr "E1" "E1" "E1" "C1" ...
## $ Project : chr "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "WL" "WL" "WL" "WL" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 5 5 5 3 3 3 1 1 1 2 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...richbox <- ggplot(plot.df, aes(x = Functional.Group, y = OTU_Richness_Sample)) +
geom_boxplot(aes(fill = Description))
richbox
richbox <- richbox + theme_bw(base_size = 15) + xlab("Functional Group") + ylab("AMF Taxa Richness")
richbox
richbox <- richbox + scale_fill_manual(values = c("red", "darkslateblue", "gold3"))
richbox
ADAPTED FROM #GLMM for roots and soil OTU Richness SSU by functional group-modified from “GLM_NAU_6_2017” #—different models for each functional group #Author: Michala Phillips #MODIFY FOR THIS DATA