AMF
M. Maltz
20180919
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] "2018Sep19"** 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("mapSal.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': 61 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]
unique(ssu$ssuspecies)## [1] "NES27"
## [2] "MH3_B"
## [3] NA
## [4] "sp"
## [5] "Acau16"
## [6] "MO_G47"
## [7] "laccatum"
## [8] "acnaGlo2"
## [9] "MO_Ar1"
## [10] "leptoticha"
## [11] "Whitfield_type_7"
## [12] "Torrecillas12b_Glo_G5"
## [13] "Sanchez_Castro12b_GLO12"
## [14] "Winther07_D"
## [15] "MO_G18"
## [16] "intraradices"
## [17] "Alguacil12a_Para_1"
## [18] "Douhan9"
## [19] "Alguacil12b_GLO_G3"
## [20] "Glo39"
## [21] "Alguacil10_Glo1"
## [22] "Ligrone07_sp"
## [23] "PT6"
## [24] "MO_G41"
## [25] "MO_GB1"
## [26] "A1"
## [27] "Liu2012b_Phylo_5"
## [28] "Yamato09_C1"
## [29] "acnaGlo7"
## [30] "MO_G7"
## [31] "Alguacil10_Glo6"
## [32] "Alguacil12b_PARA1"
## [33] "Glo58"
## [34] "lamellosum"
## [35] "VeGlo18"
## [36] "Liu2012a_Ar_1"
## [37] "Div"
## [38] "fennica"
## [39] "Voyriella_parviflora_symbiont_type_1"
## [40] "Para1_OTU2"
## [41] "Whitfield_type_3"
## [42] "Aca"
## [43] "MO_G20"
## [44] "pyriformis"
## [45] "Glo7"
## [46] "MO_GC1"
## [47] "Winther07_B"
## [48] "mosseae"
## [49] "caledonium"
## [50] "brasilianum"
## [51] "Alguacil09b_Glo_G8"
## [52] "Schechter08_Arch1"
## [53] "Wirsel_OTU21"
## [54] "Glo71"
## [55] "spurca"
## [56] "Glo59"
## [57] "Wirsel_OTU16"
## [58] "Glo32"
## [59] "MO_A8"
## [60] "MO_G27"
## [61] "MO_G5"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
# View(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), ]
unique(ssu$ssuspecies)## [1] "NES27"
## [2] "MH3_B"
## [3] "sp"
## [4] "Acau16"
## [5] "MO_G47"
## [6] "laccatum"
## [7] "acnaGlo2"
## [8] "MO_Ar1"
## [9] "leptoticha"
## [10] "Whitfield_type_7"
## [11] "Torrecillas12b_Glo_G5"
## [12] "Sanchez_Castro12b_GLO12"
## [13] "Winther07_D"
## [14] "MO_G18"
## [15] "intraradices"
## [16] "Alguacil12a_Para_1"
## [17] "Douhan9"
## [18] "Alguacil12b_GLO_G3"
## [19] "Glo39"
## [20] "Alguacil10_Glo1"
## [21] "Ligrone07_sp"
## [22] "PT6"
## [23] "MO_G41"
## [24] "MO_GB1"
## [25] "A1"
## [26] "Liu2012b_Phylo_5"
## [27] "Yamato09_C1"
## [28] "acnaGlo7"
## [29] "MO_G7"
## [30] "Alguacil10_Glo6"
## [31] "Alguacil12b_PARA1"
## [32] "Glo58"
## [33] "lamellosum"
## [34] "VeGlo18"
## [35] "Liu2012a_Ar_1"
## [36] "Div"
## [37] "fennica"
## [38] "Voyriella_parviflora_symbiont_type_1"
## [39] "Para1_OTU2"
## [40] "Whitfield_type_3"
## [41] "Aca"
## [42] "MO_G20"
## [43] "Glo7"
## [44] "MO_GC1"
## [45] "Winther07_B"
## [46] "mosseae"
## [47] "caledonium"
## [48] "brasilianum"
## [49] "Alguacil09b_Glo_G8"
## [50] "Schechter08_Arch1"
## [51] "Wirsel_OTU21"
## [52] "Glo71"
## [53] "spurca"
## [54] "Glo59"
## [55] "Wirsel_OTU16"
## [56] "Glo32"
## [57] "MO_A8"
## [58] "MO_G27"
## [59] "MO_G5"Add in functional groups
ssul <- melt(ssu, variable.name = "ID", value.name = "ssureads")
str(ssul)## 'data.frame': 17568 obs. of 11 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 "NES27" "MH3_B" "NES27" "sp" ...
## $ 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 ...OTU richnesses and read abundance for each taxonomic group
Below is in format for graphing
ssuspeciesREADRICH <- data.frame(ssul %>% group_by(ID, ssuspecies) %>% summarise(OTU_Richness_Sample = length(unique(ssuotu[ssureads >
0])), Read_Abundance_Sample = sum(ssureads)))Make SSU data frame with species level OTU richness and Read abundance for each species
ssuSpecific <- data.frame(ssul %>% group_by(ID, ssuspecies) %>% summarise(species_OTU_Richness = length(unique(ssuotu[ssureads >
0])), species_Read_Abundance = sum(ssureads)))Add metadata into ssu species
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"
ssuSpecific <- ssuSpecific %>% left_join(metassu, by = "ID")
str(ssuSpecific)## 'data.frame': 3599 obs. of 15 variables:
## $ ID : chr "MM933M" "MM933M" "MM933M" "MM933M" ...
## $ ssuspecies : chr "A1" "Aca" "Acau16" "acnaGlo2" ...
## $ species_OTU_Richness : int 1 0 1 11 0 1 1 1 4 1 ...
## $ species_Read_Abundance: num 2.88 0 3.78 63.57 0 ...
## $ #SampleID : int 5 5 5 5 5 5 5 5 5 5 ...
## $ 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 5 5 5 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...Adds meta data to data for box plot
ssuspeciesREADRICH <- ssuspeciesREADRICH %>% left_join(metassu, by = "ID")
str(ssuspeciesREADRICH)## 'data.frame': 3599 obs. of 15 variables:
## $ ID : chr "MM933M" "MM933M" "MM933M" "MM933M" ...
## $ ssuspecies : chr "A1" "Aca" "Acau16" "acnaGlo2" ...
## $ OTU_Richness_Sample : int 1 0 1 11 0 1 1 1 4 1 ...
## $ Read_Abundance_Sample: num 2.88 0 3.78 63.57 0 ...
## $ #SampleID : int 5 5 5 5 5 5 5 5 5 5 ...
## $ 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 5 5 5 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...filename <- paste0(date, "_ssu_genera_read_rich.csv")
write.csv(ssuspeciesREADRICH, file = filename, row.names = FALSE)move ID to single line and species to the length
### For ssufamily, moves ID to single line and species to the length
LssuspeciesOTURICH <- dcast(ssuSpecific, ID ~ ssuspecies, value.var = "species_OTU_Richness",
fun.aggregate = sum)
LssuspeciesOTUREAD <- dcast(ssuSpecific, ID ~ ssuspecies, value.var = "species_Read_Abundance",
fun.aggregate = sum)
head(LssuspeciesOTUREAD)## ID A1 Aca Acau16 acnaGlo2 acnaGlo7 Alguacil09b_Glo_G8
## 1 MM933M 2.8816 0.000 3.7802 63.5668 0.0000 3.7802
## 2 MM9362CU 0.0000 0.000 7.9631 116.0020 0.0000 0.0000
## 3 MM936M 0.0000 0.000 3.3956 52.6969 0.0000 4.8861
## 4 MM941M 9.8250 2.947 5.1109 63.1901 0.0000 0.0000
## 5 MM943M 0.0000 0.000 3.1296 20.7585 0.0000 0.0000
## 6 MMS700 0.0000 0.000 2.5041 48.2501 3.9087 0.0000
## Alguacil10_Glo1 Alguacil10_Glo6 Alguacil12a_Para_1 Alguacil12b_GLO_G3
## 1 5.0377 3.7802 28.4424 2.8816
## 2 6.6992 2.7709 30.9235 1.5567
## 3 5.8615 4.3254 34.4490 0.0000
## 4 4.4015 3.8503 20.5668 0.0000
## 5 5.5702 4.0447 28.4957 0.0000
## 6 2.5041 0.0000 19.2025 0.0000
## Alguacil12b_PARA1 brasilianum caledonium Div Douhan9 fennica Glo32
## 1 5.7632 0.0000 5.2934 10.1640 10.2350 7.8474 0.0000
## 2 18.5031 6.2617 4.3517 0.0000 5.1680 16.6448 6.2977
## 3 0.0000 0.0000 4.8861 0.0000 3.3956 6.7912 3.3956
## 4 0.0000 0.0000 2.9470 0.0000 8.3058 7.2795 15.2070
## 5 3.1296 0.0000 3.1296 13.2846 5.5702 9.4324 4.0447
## 6 11.0195 2.5041 2.5041 0.0000 3.3710 6.7420 5.4281
## Glo39 Glo58 Glo59 Glo7 Glo71 intraradices laccatum lamellosum
## 1 18.0629 6.1510 0.0000 0.000 0 0.0000 77.6862 12.8366
## 2 32.1172 8.4292 7.4733 0.000 0 4.7138 89.7367 10.7354
## 3 20.6828 0.0000 0.0000 0.000 0 7.7178 26.3343 15.4420
## 4 23.3941 0.0000 0.0000 2.947 0 4.7994 44.3578 14.6970
## 5 6.2592 3.1296 0.0000 0.000 0 0.0000 59.1740 7.1743
## 6 7.9777 0.0000 0.0000 0.000 0 2.5041 67.2955 27.7956
## leptoticha Ligrone07_sp Liu2012a_Ar_1 Liu2012b_Phylo_5 MH3_B MO_A8
## 1 13.2127 0.000 2.8816 0.0000 13.8731 0.0000
## 2 15.0458 0.000 1.5567 1.5567 9.6662 0.0000
## 3 15.7734 0.000 3.3956 0.0000 31.9729 0.0000
## 4 4.4015 9.996 2.9470 0.0000 12.6913 3.8503
## 5 14.0330 0.000 0.0000 0.0000 13.0897 7.0529
## 6 8.5154 0.000 4.8598 0.0000 22.2010 0.0000
## MO_Ar1 MO_G18 MO_G20 MO_G27 MO_G41 MO_G47 MO_G5 MO_G7 MO_GB1 MO_GC1
## 1 2.8816 0 3.7802 4.7267 5.0377 0 4.7267 0 0 0
## 2 13.1906 0 1.5567 6.8766 0.0000 0 5.2439 0 0 0
## 3 10.1868 0 0.0000 4.3254 0.0000 0 5.6034 0 0 0
## 4 0.0000 0 2.9470 4.7994 9.6553 0 5.7732 0 0 0
## 5 5.0003 0 0.0000 5.0003 0.0000 0 5.7883 0 0 0
## 6 16.8826 0 5.5768 5.4281 6.8221 0 4.8598 0 0 0
## mosseae NES27 Para1_OTU2 PT6 Sanchez_Castro12b_GLO12
## 1 8.3195 44.9255 2.8816 5.7632 76.8618
## 2 10.9900 48.3770 8.3194 3.1134 73.1450
## 3 8.0643 39.6186 7.5810 6.7912 43.8708
## 4 9.4211 20.3713 4.4015 11.0910 114.9387
## 5 8.1296 42.5804 0.0000 30.6556 55.5911
## 6 8.9687 45.6336 0.0000 10.8650 91.9122
## Schechter08_Arch1 sp spurca Torrecillas12b_Glo_G5 VeGlo18
## 1 0.0000 280.6634 0 67.5200 0.0000
## 2 0.0000 161.5516 0 85.3482 5.0879
## 3 3.3956 196.8241 0 72.7028 0.0000
## 4 0.0000 209.2776 0 70.6345 2.9470
## 5 0.0000 216.9472 0 58.8557 0.0000
## 6 0.0000 238.9615 0 71.1061 0.0000
## Voyriella_parviflora_symbiont_type_1 Whitfield_type_3 Whitfield_type_7
## 1 0.0000 0.0000 13.8353
## 2 4.7138 6.6992 10.9385
## 3 0.0000 3.3956 24.9410
## 4 9.4629 6.3494 17.7806
## 5 0.0000 0.0000 19.1941
## 6 8.5137 3.3710 8.9478
## Winther07_B Winther07_D Wirsel_OTU16 Wirsel_OTU21 Yamato09_C1
## 1 0.0000 6.4946 4.3297 0.0000 2.8816
## 2 0.0000 12.8325 8.7551 1.5567 13.5632
## 3 10.2190 6.7912 5.2889 0.0000 0.0000
## 4 4.4015 13.4651 2.9470 0.0000 0.0000
## 5 6.1450 13.9366 0.0000 0.0000 0.0000
## 6 0.0000 10.5125 5.4281 0.0000 12.5768Genera reads and richness together ###WORKED ON THIS 20180826
REPLACED Functional.Group with ssuspecies
Work on changing to species in the next iteration
ssuspeciesREADRICHlong<- data.frame(ssul %>%
group_by(ID) %>%
summarise(OTU_Richness_Sample = length(unique(ssuotu[ssureads>0])),
NES27_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="NES27"])),
MH3_B_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="MH3_B"])),
Acau16_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Acau16"])),
MO_G47_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="MO_G47"])),
laccatum_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="laccatum"])),
acnaGlo2_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="acnaGlo2"])),
MO_Ar1_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="MO_Ar1"])),
leptoticha_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="leptoticha"])),
Whitfield_type_7_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Whitfield_type_7"])),
Torrecillas12b_Glo_G5_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Torrecillas12b_Glo_G5"])),
Sanchez_Castro12b_GLO12_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Sanchez_Castro12b_GLO12"])),
Winther07_D_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Winther07_D"])),
MO_G18_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="MO_G18"])),
intraradices_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="intraradices"])),
Alguacil12a_Para_1_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Alguacil12a_Para_1"])),
Douhan9_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Douhan9"])),
Alguacil12b_GLO_G3_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Alguacil12b_GLO_G3"])),
Glo39_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Glo39"])),
Alguacil10_Glo1_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Alguacil10_Glo1"])),
Ligrone07_sp_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Ligrone07_sp"])),
PT6_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="PT6"])),
MO_G41_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="MO_G41"])),
MO_GB1_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="MO_GB1"])),
A1_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="A1"])),
Liu2012b_Phylo_5_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Liu2012b_Phylo_5"])),
Yamato09_C1_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Yamato09_C1"])),
acnaGlo7_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="acnaGlo7"])),
MO_G7_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="MO_G7"])),
Alguacil10_Glo6_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Alguacil10_Glo6"])),
Alguacil12b_PARA1_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Alguacil12b_PARA1"])),
Glo58_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Glo58"])),
lamellosum_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="lamellosum"])),
VeGlo18_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="VeGlo18"])),
Liu2012a_Ar_1_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Liu2012a_Ar_1"])),
Div_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Div"])),
fennica_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="fennica"])),
Voyriella_parviflora_symbiont_type_1_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Voyriella_parviflora_symbiont_type_1"])),
Para1_OTU2_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Para1_OTU2"])),
Whitfield_type_3_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Whitfield_type_3"])),
Aca_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Aca"])),
MO_G20_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="MO_G20"])),
Glo7_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Glo7"])),
MO_GC1_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="MO_GC1"])),
Winther07_B_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Winther07_B"])),
mosseae_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="mosseae"])),
caledonium_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="caledonium"])),
brasilianum_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="brasilianum"])),
Alguacil09b_Glo_G8_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Alguacil09b_Glo_G8"])),
Schechter08_Arch1_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Schechter08_Arch1"])),
Wirsel_OTU21_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Wirsel_OTU21"])),
Glo71_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Glo71"])),
spurca_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="spurca"])),
Glo59_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Glo59"])),
Wirsel_OTU16_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Wirsel_OTU16"])),
Glo32_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="Glo32"])),
MO_A8_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="MO_A8"])),
MO_G27_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="MO_G27"])),
MO_G5_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="MO_G5"])),
Glomussp_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="sp" & ssugenus == "Glomus"])),
Archaeosporasp_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="sp" & ssugenus == "Archaeospora"])),
Diversisporasp_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="sp" & ssugenus == "Diversispora"])),
Paraglomussp_Richness = length(unique(ssuotu[ssureads>0 & ssuspecies =="sp" & ssugenus == "Paraglomus"])),
NES27_Read_Abundance = sum(ssureads[ssuspecies == "NES27"]),
MH3_B_Read_Abundance = sum(ssureads[ssuspecies == "MH3_B"]),
Acau16_Read_Abundance = sum(ssureads[ssuspecies == "Acau16"]),
MO_G47_Read_Abundance = sum(ssureads[ssuspecies == "MO_G47"]),
laccatum_Read_Abundance = sum(ssureads[ssuspecies == "laccatum"]),
acnaGlo2_Read_Abundance = sum(ssureads[ssuspecies == "acnaGlo2"]),
MO_Ar1 = sum(ssureads[ssuspecies == "MO_Ar1"]),
leptoticha_Read_Abundance = sum(ssureads[ssuspecies == "leptoticha"]),
Whitfield_type_7_Read_Abundance = sum(ssureads[ssuspecies == "Whitfield_type_7"]),
Torrecillas12b_Glo_G5_Read_Abundance = sum(ssureads[ssuspecies == "Torrecillas12b_Glo_G5"]),
Sanchez_Castro12b_GLO12_Read_Abundance = sum(ssureads[ssuspecies == "Sanchez_Castro12b_GLO12"]),
Winther07_D_Read_Abundance = sum(ssureads[ssuspecies == "Winther07_D"]),
MO_G18_Read_Abundance = sum(ssureads[ssuspecies == "MO_G18"]),
intraradices_Read_Abundance = sum(ssureads[ssuspecies == "intraradices"]),
Alguacil12a_Para_1_Read_Abundance = sum(ssureads[ssuspecies == "Alguacil12a_Para_1"]),
Douhan9_Read_Abundance = sum(ssureads[ssuspecies == "Douhan9"]),
Alguacil12b_GLO_G3_Read_Abundance = sum(ssureads[ssuspecies == "Alguacil12b_GLO_G3"]),
Glo39_Read_Abundance = sum(ssureads[ssuspecies == "Glo39"]),
Alguacil10_Glo1_Read_Abundance = sum(ssureads[ssuspecies == "Alguacil10_Glo1"]),
Ligrone07_sp_Read_Abundance = sum(ssureads[ssuspecies == "Ligrone07_sp"]),
PT6_Read_Abundance = sum(ssureads[ssuspecies == "PT6"]),
MO_G41_Read_Abundance = sum(ssureads[ssuspecies == "MO_G41"]),
NES27_Read_Abundance = sum(ssureads[ssuspecies == "MO_GB1"]),
A1_Read_Abundance = sum(ssureads[ssuspecies == "A1"]),
Liu2012b_Phylo_5_Read_Abundance = sum(ssureads[ssuspecies == "Liu2012b_Phylo_5"]),
Yamato09_C1_Read_Abundance = sum(ssureads[ssuspecies == "Yamato09_C1"]),
NES27_Read_Abundance = sum(ssureads[ssuspecies == "acnaGlo7"]),
MO_G7_Read_Abundance = sum(ssureads[ssuspecies == "MO_G7"]),
Alguacil10_Glo6_Read_Abundance = sum(ssureads[ssuspecies == "Alguacil10_Glo6"]),
Alguacil12b_PARA1_Read_Abundance = sum(ssureads[ssuspecies == "Alguacil12b_PARA1"]),
Glo58_Read_Abundance = sum(ssureads[ssuspecies == "Glo58"]),
lamellosum_Read_Abundance = sum(ssureads[ssuspecies == "lamellosum"]),
VeGlo18_Read_Abundance = sum(ssureads[ssuspecies == "VeGlo18"]),
Liu2012a_Ar_1_Read_Abundance = sum(ssureads[ssuspecies == "Liu2012a_Ar_1"]),
Div_Read_Abundance = sum(ssureads[ssuspecies == "Div"]),
fennica_Read_Abundance = sum(ssureads[ssuspecies == "fennica"]),
Voyriella_parviflora_symbiont_type_1_Read_Abundance = sum(ssureads[ssuspecies == "Voyriella_parviflora_symbiont_type_1"]),
Para1_OTU2_Read_Abundance = sum(ssureads[ssuspecies == "Para1_OTU2"]),
Whitfield_type_3_Read_Abundance = sum(ssureads[ssuspecies == "Whitfield_type_3"]),
Aca_Read_Abundance = sum(ssureads[ssuspecies == "Aca"]),
MO_G20_Read_Abundance = sum(ssureads[ssuspecies == "MO_G20"]),
Glo7_Read_Abundance = sum(ssureads[ssuspecies == "Glo7"]),
MO_GC1_Read_Abundance = sum(ssureads[ssuspecies == "MO_GC1"]),
Winther07_B7_Read_Abundance = sum(ssureads[ssuspecies == "Winther07_B"]),
mosseae_Read_Abundance = sum(ssureads[ssuspecies == "mosseae"]),
caledonium_Read_Abundance = sum(ssureads[ssuspecies == "caledonium"]),
brasilianum_Read_Abundance = sum(ssureads[ssuspecies == "brasilianum"]),
Alguacil09b_Glo_G8_Read_Abundance = sum(ssureads[ssuspecies == "Alguacil09b_Glo_G8"]),
Schechter08_Arch1_Read_Abundance = sum(ssureads[ssuspecies == "Schechter08_Arch1"]),
Wirsel_OTU21_Read_Abundance = sum(ssureads[ssuspecies == "Wirsel_OTU21"]),
Glo71_Read_Abundance = sum(ssureads[ssuspecies == "Glo71"]),
spurca_Read_Abundance = sum(ssureads[ssuspecies == "spurca"]),
Glo59_Read_Abundance = sum(ssureads[ssuspecies == "Glo59"]),
Wirsel_OTU16_Read_Abundance = sum(ssureads[ssuspecies == "Wirsel_OTU16"]),
Glo32_Read_Abundance = sum(ssureads[ssuspecies == "Glo32"]),
MO_A8_Read_Abundance = sum(ssureads[ssuspecies == "MO_A8"]),
MO_G27_Read_Abundance = sum(ssureads[ssuspecies == "MO_G27"]),
MO_G5_Read_Abundance = sum(ssureads[ssuspecies == "MO_G5"]),
Glomussp_Read_Abundance = sum(ssureads[ssuspecies == "sp"& ssugenus == "Glomus"]),
Archaeosporasp_Read_Abundance = sum(ssureads[ssuspecies == "sp"& ssugenus == "Archaeospora"]),
Diversisporasp_Read_Abundance = sum(ssureads[ssuspecies == "sp"& ssugenus == "Diversispora"]),
Paraglomussp_Read_Abundance = sum(ssureads[ssuspecies == "sp"& ssugenus == "Paraglomus"]),
Read_Abundance_Sample = sum(ssureads)))
##View(ssuspeciesREADRICHlong)
ssuspeciesREADRICHlong <- merge(ssuspeciesREADRICHlong, metassu, by = "ID")
filename<-paste0(date, '_ssu_species_read_rich.csv')
write.csv(ssuspeciesREADRICHlong, file = filename, row.names=FALSE)
str(ssuspeciesREADRICHlong)## 'data.frame': 61 obs. of 136 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 ...
## $ NES27_Richness : int 9 11 6 5 6 11 13 10 15 14 ...
## $ MH3_B_Richness : int 4 4 4 4 3 4 4 4 2 2 ...
## $ Acau16_Richness : int 1 1 1 1 1 1 2 2 0 1 ...
## $ MO_G47_Richness : int 0 0 0 0 0 0 0 0 1 0 ...
## $ laccatum_Richness : int 12 17 6 7 8 12 18 20 11 17 ...
## $ acnaGlo2_Richness : int 11 15 7 10 3 11 10 11 10 12 ...
## $ MO_Ar1_Richness : int 1 2 3 0 1 2 2 3 2 2 ...
## $ leptoticha_Richness : int 3 3 2 1 2 2 2 2 4 4 ...
## $ Whitfield_type_7_Richness : int 3 2 2 3 2 2 2 2 2 3 ...
## $ Torrecillas12b_Glo_G5_Richness : int 9 12 8 9 7 11 13 14 10 11 ...
## $ Sanchez_Castro12b_GLO12_Richness : int 9 14 7 12 9 16 19 17 14 14 ...
## $ Winther07_D_Richness : int 1 3 2 3 2 2 3 1 2 3 ...
## $ MO_G18_Richness : int 0 0 0 0 0 0 0 0 1 1 ...
## $ intraradices_Richness : int 0 1 1 1 0 1 1 1 1 0 ...
## $ Alguacil12a_Para_1_Richness : int 4 5 3 3 3 3 5 4 2 5 ...
## $ Douhan9_Richness : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Alguacil12b_GLO_G3_Richness : int 1 1 0 0 0 0 0 2 3 1 ...
## $ Glo39_Richness : int 3 4 4 4 2 2 2 1 4 4 ...
## $ Alguacil10_Glo1_Richness : int 1 1 1 1 1 1 1 0 1 1 ...
## $ Ligrone07_sp_Richness : int 0 0 0 1 0 0 0 1 0 0 ...
## $ PT6_Richness : int 2 2 2 2 3 1 0 2 1 1 ...
## $ MO_G41_Richness : int 1 0 0 1 0 1 1 0 0 0 ...
## $ MO_GB1_Richness : int 0 0 0 0 0 0 0 0 1 1 ...
## $ A1_Richness : int 1 0 0 1 0 0 0 0 0 0 ...
## $ Liu2012b_Phylo_5_Richness : int 0 1 0 0 0 0 0 0 0 0 ...
## $ Yamato09_C1_Richness : int 1 3 0 0 0 3 2 2 3 2 ...
## $ acnaGlo7_Richness : int 0 0 0 0 0 1 0 0 0 0 ...
## $ MO_G7_Richness : int 0 0 0 0 0 0 0 1 1 0 ...
## $ Alguacil10_Glo6_Richness : int 1 1 1 1 1 0 0 1 2 0 ...
## $ Alguacil12b_PARA1_Richness : int 2 3 0 0 1 3 3 3 3 2 ...
## $ Glo58_Richness : int 1 1 0 0 1 0 0 0 0 0 ...
## $ lamellosum_Richness : int 3 4 4 4 2 6 3 4 5 4 ...
## $ VeGlo18_Richness : int 0 1 0 1 0 0 0 1 1 0 ...
## $ Liu2012a_Ar_1_Richness : int 1 1 1 1 0 1 0 2 1 0 ...
## $ Div_Richness : int 1 0 0 0 2 0 0 2 1 0 ...
## $ fennica_Richness : int 1 3 2 1 1 2 3 3 3 2 ...
## $ Voyriella_parviflora_symbiont_type_1_Richness : int 0 1 0 1 0 1 1 1 0 0 ...
## $ Para1_OTU2_Richness : int 1 1 1 1 0 0 0 0 0 0 ...
## $ Whitfield_type_3_Richness : int 0 1 1 1 0 1 0 0 1 0 ...
## $ Aca_Richness : int 0 0 0 1 0 0 0 0 0 0 ...
## $ MO_G20_Richness : int 1 1 0 1 0 1 0 0 0 0 ...
## $ Glo7_Richness : int 0 0 0 1 0 0 0 0 0 0 ...
## $ MO_GC1_Richness : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Winther07_B_Richness : int 0 0 1 1 1 0 0 0 0 1 ...
## $ mosseae_Richness : int 1 1 1 1 1 1 1 1 1 1 ...
## $ caledonium_Richness : int 1 1 1 1 1 1 1 1 0 2 ...
## $ brasilianum_Richness : int 0 1 0 0 0 1 1 1 0 1 ...
## $ Alguacil09b_Glo_G8_Richness : int 1 0 1 0 0 0 1 1 1 0 ...
## $ Schechter08_Arch1_Richness : int 0 0 1 0 0 0 0 0 0 0 ...
## $ Wirsel_OTU21_Richness : int 0 1 0 0 0 0 0 0 0 1 ...
## $ Glo71_Richness : int 0 0 0 0 0 0 1 0 0 0 ...
## $ spurca_Richness : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Glo59_Richness : int 0 1 0 0 0 0 0 0 0 1 ...
## $ Wirsel_OTU16_Richness : int 1 1 1 1 0 1 1 1 1 0 ...
## $ Glo32_Richness : int 0 1 1 1 1 1 0 1 1 1 ...
## $ MO_A8_Richness : int 0 0 0 1 1 0 1 0 0 0 ...
## $ MO_G27_Richness : int 1 1 1 1 1 1 1 1 1 1 ...
## $ MO_G5_Richness : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Glomussp_Richness : int 13 15 15 15 16 20 20 19 12 14 ...
## $ Archaeosporasp_Richness : int 16 13 8 9 10 13 11 10 6 10 ...
## $ Diversisporasp_Richness : int 1 3 1 1 1 3 2 2 1 1 ...
## $ Paraglomussp_Richness : int 3 4 2 2 2 4 4 4 2 3 ...
## $ NES27_Read_Abundance : num 0 0 0 0 0 ...
## $ MH3_B_Read_Abundance : num 13.87 9.67 31.97 12.69 13.09 ...
## $ Acau16_Read_Abundance : num 3.78 7.96 3.4 5.11 3.13 ...
## $ MO_G47_Read_Abundance : num 0 0 0 0 0 ...
## $ laccatum_Read_Abundance : num 77.7 89.7 26.3 44.4 59.2 ...
## $ acnaGlo2_Read_Abundance : num 63.6 116 52.7 63.2 20.8 ...
## $ MO_Ar1 : num 2.88 13.19 10.19 0 5 ...
## $ leptoticha_Read_Abundance : num 13.2 15 15.8 4.4 14 ...
## $ Whitfield_type_7_Read_Abundance : num 13.8 10.9 24.9 17.8 19.2 ...
## $ Torrecillas12b_Glo_G5_Read_Abundance : num 67.5 85.3 72.7 70.6 58.9 ...
## $ Sanchez_Castro12b_GLO12_Read_Abundance : num 76.9 73.1 43.9 114.9 55.6 ...
## $ Winther07_D_Read_Abundance : num 6.49 12.83 6.79 13.47 13.94 ...
## $ MO_G18_Read_Abundance : num 0 0 0 0 0 ...
## $ intraradices_Read_Abundance : num 0 4.71 7.72 4.8 0 ...
## $ Alguacil12a_Para_1_Read_Abundance : num 28.4 30.9 34.4 20.6 28.5 ...
## $ Douhan9_Read_Abundance : num 10.23 5.17 3.4 8.31 5.57 ...
## $ Alguacil12b_GLO_G3_Read_Abundance : num 2.88 1.56 0 0 0 ...
## $ Glo39_Read_Abundance : num 18.06 32.12 20.68 23.39 6.26 ...
## $ Alguacil10_Glo1_Read_Abundance : num 5.04 6.7 5.86 4.4 5.57 ...
## $ Ligrone07_sp_Read_Abundance : num 0 0 0 10 0 ...
## $ PT6_Read_Abundance : num 5.76 3.11 6.79 11.09 30.66 ...
## $ MO_G41_Read_Abundance : num 5.04 0 0 9.66 0 ...
## $ A1_Read_Abundance : num 2.88 0 0 9.82 0 ...
## $ Liu2012b_Phylo_5_Read_Abundance : num 0 1.56 0 0 0 ...
## $ Yamato09_C1_Read_Abundance : num 2.88 13.56 0 0 0 ...
## $ MO_G7_Read_Abundance : num 0 0 0 0 0 ...
## $ Alguacil10_Glo6_Read_Abundance : num 3.78 2.77 4.33 3.85 4.04 ...
## $ Alguacil12b_PARA1_Read_Abundance : num 5.76 18.5 0 0 3.13 ...
## $ Glo58_Read_Abundance : num 6.15 8.43 0 0 3.13 ...
## $ lamellosum_Read_Abundance : num 12.84 10.74 15.44 14.7 7.17 ...
## $ VeGlo18_Read_Abundance : num 0 5.09 0 2.95 0 ...
## $ Liu2012a_Ar_1_Read_Abundance : num 2.88 1.56 3.4 2.95 0 ...
## $ Div_Read_Abundance : num 10.2 0 0 0 13.3 ...
## $ fennica_Read_Abundance : num 7.85 16.64 6.79 7.28 9.43 ...
## $ Voyriella_parviflora_symbiont_type_1_Read_Abundance: num 0 4.71 0 9.46 0 ...
## [list output truncated]##View(ssuspeciesREADRICHlong)Non Unifrac distances, try with bray
WORK ON THIS
Permanova
Permanova tests for differences in composition among groups
Reminder - permanova is always based on pairwise distances/dissimilarities.
Can either nest Treat in Year with Year/Treat Or, can restrict permutations within year, using strata=‘Year’ Or, can remove Year, and focus on Treat and other grouping varuables, like Site
Also, dist.j<-vegdist(ssuspeciesREADRICHlong[,-1], method=‘jaccard’, na.rm=TRUE)
USE BRAY #dist.bc<-vegdist(ssuspeciesREADRICHlong[,-1], method=‘bray’, na.rm=TRUE)
Part 3 - multivariate community analyses and data visualization
library(dplyr)
library(reshape2)
library(tidyr)
library(vegan)
library(ggplot2) # ggplot resource -http://rpubs.com/collnell/ggplot2
library(tidyverse) # useful packages in 1 - dplyr, ggplot2, tidyr +community matrix at order level with mapping data
otu_map<-read.csv(‘/Users/maltz/Desktop/RdataBreathe/otus_by_class.csv’)
head(ssuspeciesREADRICHlong)## ID OTU_Richness_Sample NES27_Richness MH3_B_Richness
## 1 MM933M 130 9 4
## 2 MM9362CU 167 11 4
## 3 MM936M 105 6 4
## 4 MM941M 121 5 4
## 5 MM943M 99 6 3
## 6 MMS700 151 11 4
## Acau16_Richness MO_G47_Richness laccatum_Richness acnaGlo2_Richness
## 1 1 0 12 11
## 2 1 0 17 15
## 3 1 0 6 7
## 4 1 0 7 10
## 5 1 0 8 3
## 6 1 0 12 11
## MO_Ar1_Richness leptoticha_Richness Whitfield_type_7_Richness
## 1 1 3 3
## 2 2 3 2
## 3 3 2 2
## 4 0 1 3
## 5 1 2 2
## 6 2 2 2
## Torrecillas12b_Glo_G5_Richness Sanchez_Castro12b_GLO12_Richness
## 1 9 9
## 2 12 14
## 3 8 7
## 4 9 12
## 5 7 9
## 6 11 16
## Winther07_D_Richness MO_G18_Richness intraradices_Richness
## 1 1 0 0
## 2 3 0 1
## 3 2 0 1
## 4 3 0 1
## 5 2 0 0
## 6 2 0 1
## Alguacil12a_Para_1_Richness Douhan9_Richness Alguacil12b_GLO_G3_Richness
## 1 4 1 1
## 2 5 1 1
## 3 3 1 0
## 4 3 1 0
## 5 3 1 0
## 6 3 1 0
## Glo39_Richness Alguacil10_Glo1_Richness Ligrone07_sp_Richness
## 1 3 1 0
## 2 4 1 0
## 3 4 1 0
## 4 4 1 1
## 5 2 1 0
## 6 2 1 0
## PT6_Richness MO_G41_Richness MO_GB1_Richness A1_Richness
## 1 2 1 0 1
## 2 2 0 0 0
## 3 2 0 0 0
## 4 2 1 0 1
## 5 3 0 0 0
## 6 1 1 0 0
## Liu2012b_Phylo_5_Richness Yamato09_C1_Richness acnaGlo7_Richness
## 1 0 1 0
## 2 1 3 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 3 1
## MO_G7_Richness Alguacil10_Glo6_Richness Alguacil12b_PARA1_Richness
## 1 0 1 2
## 2 0 1 3
## 3 0 1 0
## 4 0 1 0
## 5 0 1 1
## 6 0 0 3
## Glo58_Richness lamellosum_Richness VeGlo18_Richness
## 1 1 3 0
## 2 1 4 1
## 3 0 4 0
## 4 0 4 1
## 5 1 2 0
## 6 0 6 0
## Liu2012a_Ar_1_Richness Div_Richness fennica_Richness
## 1 1 1 1
## 2 1 0 3
## 3 1 0 2
## 4 1 0 1
## 5 0 2 1
## 6 1 0 2
## Voyriella_parviflora_symbiont_type_1_Richness Para1_OTU2_Richness
## 1 0 1
## 2 1 1
## 3 0 1
## 4 1 1
## 5 0 0
## 6 1 0
## Whitfield_type_3_Richness Aca_Richness MO_G20_Richness Glo7_Richness
## 1 0 0 1 0
## 2 1 0 1 0
## 3 1 0 0 0
## 4 1 1 1 1
## 5 0 0 0 0
## 6 1 0 1 0
## MO_GC1_Richness Winther07_B_Richness mosseae_Richness
## 1 0 0 1
## 2 0 0 1
## 3 0 1 1
## 4 0 1 1
## 5 0 1 1
## 6 0 0 1
## caledonium_Richness brasilianum_Richness Alguacil09b_Glo_G8_Richness
## 1 1 0 1
## 2 1 1 0
## 3 1 0 1
## 4 1 0 0
## 5 1 0 0
## 6 1 1 0
## Schechter08_Arch1_Richness Wirsel_OTU21_Richness Glo71_Richness
## 1 0 0 0
## 2 0 1 0
## 3 1 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 0 0
## spurca_Richness Glo59_Richness Wirsel_OTU16_Richness Glo32_Richness
## 1 0 0 1 0
## 2 0 1 1 1
## 3 0 0 1 1
## 4 0 0 1 1
## 5 0 0 0 1
## 6 0 0 1 1
## MO_A8_Richness MO_G27_Richness MO_G5_Richness Glomussp_Richness
## 1 0 1 1 13
## 2 0 1 1 15
## 3 0 1 1 15
## 4 1 1 1 15
## 5 1 1 1 16
## 6 0 1 1 20
## Archaeosporasp_Richness Diversisporasp_Richness Paraglomussp_Richness
## 1 16 1 3
## 2 13 3 4
## 3 8 1 2
## 4 9 1 2
## 5 10 1 2
## 6 13 3 4
## NES27_Read_Abundance MH3_B_Read_Abundance Acau16_Read_Abundance
## 1 0.0000 13.8731 3.7802
## 2 0.0000 9.6662 7.9631
## 3 0.0000 31.9729 3.3956
## 4 0.0000 12.6913 5.1109
## 5 0.0000 13.0897 3.1296
## 6 3.9087 22.2010 2.5041
## MO_G47_Read_Abundance laccatum_Read_Abundance acnaGlo2_Read_Abundance
## 1 0 77.6862 63.5668
## 2 0 89.7367 116.0020
## 3 0 26.3343 52.6969
## 4 0 44.3578 63.1901
## 5 0 59.1740 20.7585
## 6 0 67.2955 48.2501
## MO_Ar1 leptoticha_Read_Abundance Whitfield_type_7_Read_Abundance
## 1 2.8816 13.2127 13.8353
## 2 13.1906 15.0458 10.9385
## 3 10.1868 15.7734 24.9410
## 4 0.0000 4.4015 17.7806
## 5 5.0003 14.0330 19.1941
## 6 16.8826 8.5154 8.9478
## Torrecillas12b_Glo_G5_Read_Abundance
## 1 67.5200
## 2 85.3482
## 3 72.7028
## 4 70.6345
## 5 58.8557
## 6 71.1061
## Sanchez_Castro12b_GLO12_Read_Abundance Winther07_D_Read_Abundance
## 1 76.8618 6.4946
## 2 73.1450 12.8325
## 3 43.8708 6.7912
## 4 114.9387 13.4651
## 5 55.5911 13.9366
## 6 91.9122 10.5125
## MO_G18_Read_Abundance intraradices_Read_Abundance
## 1 0 0.0000
## 2 0 4.7138
## 3 0 7.7178
## 4 0 4.7994
## 5 0 0.0000
## 6 0 2.5041
## Alguacil12a_Para_1_Read_Abundance Douhan9_Read_Abundance
## 1 28.4424 10.2350
## 2 30.9235 5.1680
## 3 34.4490 3.3956
## 4 20.5668 8.3058
## 5 28.4957 5.5702
## 6 19.2025 3.3710
## Alguacil12b_GLO_G3_Read_Abundance Glo39_Read_Abundance
## 1 2.8816 18.0629
## 2 1.5567 32.1172
## 3 0.0000 20.6828
## 4 0.0000 23.3941
## 5 0.0000 6.2592
## 6 0.0000 7.9777
## Alguacil10_Glo1_Read_Abundance Ligrone07_sp_Read_Abundance
## 1 5.0377 0.000
## 2 6.6992 0.000
## 3 5.8615 0.000
## 4 4.4015 9.996
## 5 5.5702 0.000
## 6 2.5041 0.000
## PT6_Read_Abundance MO_G41_Read_Abundance A1_Read_Abundance
## 1 5.7632 5.0377 2.8816
## 2 3.1134 0.0000 0.0000
## 3 6.7912 0.0000 0.0000
## 4 11.0910 9.6553 9.8250
## 5 30.6556 0.0000 0.0000
## 6 10.8650 6.8221 0.0000
## Liu2012b_Phylo_5_Read_Abundance Yamato09_C1_Read_Abundance
## 1 0.0000 2.8816
## 2 1.5567 13.5632
## 3 0.0000 0.0000
## 4 0.0000 0.0000
## 5 0.0000 0.0000
## 6 0.0000 12.5768
## MO_G7_Read_Abundance Alguacil10_Glo6_Read_Abundance
## 1 0 3.7802
## 2 0 2.7709
## 3 0 4.3254
## 4 0 3.8503
## 5 0 4.0447
## 6 0 0.0000
## Alguacil12b_PARA1_Read_Abundance Glo58_Read_Abundance
## 1 5.7632 6.1510
## 2 18.5031 8.4292
## 3 0.0000 0.0000
## 4 0.0000 0.0000
## 5 3.1296 3.1296
## 6 11.0195 0.0000
## lamellosum_Read_Abundance VeGlo18_Read_Abundance
## 1 12.8366 0.0000
## 2 10.7354 5.0879
## 3 15.4420 0.0000
## 4 14.6970 2.9470
## 5 7.1743 0.0000
## 6 27.7956 0.0000
## Liu2012a_Ar_1_Read_Abundance Div_Read_Abundance fennica_Read_Abundance
## 1 2.8816 10.1640 7.8474
## 2 1.5567 0.0000 16.6448
## 3 3.3956 0.0000 6.7912
## 4 2.9470 0.0000 7.2795
## 5 0.0000 13.2846 9.4324
## 6 4.8598 0.0000 6.7420
## Voyriella_parviflora_symbiont_type_1_Read_Abundance
## 1 0.0000
## 2 4.7138
## 3 0.0000
## 4 9.4629
## 5 0.0000
## 6 8.5137
## Para1_OTU2_Read_Abundance Whitfield_type_3_Read_Abundance
## 1 2.8816 0.0000
## 2 8.3194 6.6992
## 3 7.5810 3.3956
## 4 4.4015 6.3494
## 5 0.0000 0.0000
## 6 0.0000 3.3710
## Aca_Read_Abundance MO_G20_Read_Abundance Glo7_Read_Abundance
## 1 0.000 3.7802 0.000
## 2 0.000 1.5567 0.000
## 3 0.000 0.0000 0.000
## 4 2.947 2.9470 2.947
## 5 0.000 0.0000 0.000
## 6 0.000 5.5768 0.000
## MO_GC1_Read_Abundance Winther07_B7_Read_Abundance mosseae_Read_Abundance
## 1 0 0.0000 8.3195
## 2 0 0.0000 10.9900
## 3 0 10.2190 8.0643
## 4 0 4.4015 9.4211
## 5 0 6.1450 8.1296
## 6 0 0.0000 8.9687
## caledonium_Read_Abundance brasilianum_Read_Abundance
## 1 5.2934 0.0000
## 2 4.3517 6.2617
## 3 4.8861 0.0000
## 4 2.9470 0.0000
## 5 3.1296 0.0000
## 6 2.5041 2.5041
## Alguacil09b_Glo_G8_Read_Abundance Schechter08_Arch1_Read_Abundance
## 1 3.7802 0.0000
## 2 0.0000 0.0000
## 3 4.8861 3.3956
## 4 0.0000 0.0000
## 5 0.0000 0.0000
## 6 0.0000 0.0000
## Wirsel_OTU21_Read_Abundance Glo71_Read_Abundance spurca_Read_Abundance
## 1 0.0000 0 0
## 2 1.5567 0 0
## 3 0.0000 0 0
## 4 0.0000 0 0
## 5 0.0000 0 0
## 6 0.0000 0 0
## Glo59_Read_Abundance Wirsel_OTU16_Read_Abundance Glo32_Read_Abundance
## 1 0.0000 4.3297 0.0000
## 2 7.4733 8.7551 6.2977
## 3 0.0000 5.2889 3.3956
## 4 0.0000 2.9470 15.2070
## 5 0.0000 0.0000 4.0447
## 6 0.0000 5.4281 5.4281
## MO_A8_Read_Abundance MO_G27_Read_Abundance MO_G5_Read_Abundance
## 1 0.0000 4.7267 4.7267
## 2 0.0000 6.8766 5.2439
## 3 0.0000 4.3254 5.6034
## 4 3.8503 4.7994 5.7732
## 5 7.0529 5.0003 5.7883
## 6 0.0000 5.4281 4.8598
## Glomussp_Read_Abundance Archaeosporasp_Read_Abundance
## 1 95.1801 154.7074
## 2 71.2978 48.5210
## 3 117.6183 54.1959
## 4 125.3506 44.9167
## 5 103.8362 87.1494
## 6 119.1184 73.4597
## Diversisporasp_Read_Abundance Paraglomussp_Read_Abundance
## 1 13.0340 14.8603
## 2 18.3924 19.0128
## 3 12.2380 12.7719
## 4 13.9450 12.0393
## 5 9.6978 12.2191
## 6 18.6176 18.3299
## Read_Abundance_Sample #SampleID BarcodeSequence LinkerPrimerSequence
## 1 833.7869 5 ATCCCGTATCGATTGG NA
## 2 886.0327 4 GCAACCTTTCGATTGG NA
## 3 695.0015 3 GCAACCTTAGAGTGTG NA
## 4 788.3764 1 GCAACCTTTGGGTGAT NA
## 5 678.3267 2 GCAACCTTTACTGTGC NA
## 6 805.4538 47 GAAGATCCCTCTCAAG NA
## Location Project Year Site Treat Rep TSDel Description
## 1 E1 MM-Salvage 2015 WL O 5 <NA> RecipientPre
## 2 D1 MM-Salvage 2015 WL O 4 <NA> RecipientPre
## 3 C1 MM-Salvage 2015 WL O 3 <NA> RecipientPre
## 4 A1 MM-Salvage 2015 WL O 1 <NA> RecipientPre
## 5 B1 MM-Salvage 2015 WL O 2 <NA> RecipientPre
## 6 G6 MM-Salvage 2017 WL D 1 B1A2C RecipientPostcolnames(ssuspeciesREADRICHlong)## [1] "ID"
## [2] "OTU_Richness_Sample"
## [3] "NES27_Richness"
## [4] "MH3_B_Richness"
## [5] "Acau16_Richness"
## [6] "MO_G47_Richness"
## [7] "laccatum_Richness"
## [8] "acnaGlo2_Richness"
## [9] "MO_Ar1_Richness"
## [10] "leptoticha_Richness"
## [11] "Whitfield_type_7_Richness"
## [12] "Torrecillas12b_Glo_G5_Richness"
## [13] "Sanchez_Castro12b_GLO12_Richness"
## [14] "Winther07_D_Richness"
## [15] "MO_G18_Richness"
## [16] "intraradices_Richness"
## [17] "Alguacil12a_Para_1_Richness"
## [18] "Douhan9_Richness"
## [19] "Alguacil12b_GLO_G3_Richness"
## [20] "Glo39_Richness"
## [21] "Alguacil10_Glo1_Richness"
## [22] "Ligrone07_sp_Richness"
## [23] "PT6_Richness"
## [24] "MO_G41_Richness"
## [25] "MO_GB1_Richness"
## [26] "A1_Richness"
## [27] "Liu2012b_Phylo_5_Richness"
## [28] "Yamato09_C1_Richness"
## [29] "acnaGlo7_Richness"
## [30] "MO_G7_Richness"
## [31] "Alguacil10_Glo6_Richness"
## [32] "Alguacil12b_PARA1_Richness"
## [33] "Glo58_Richness"
## [34] "lamellosum_Richness"
## [35] "VeGlo18_Richness"
## [36] "Liu2012a_Ar_1_Richness"
## [37] "Div_Richness"
## [38] "fennica_Richness"
## [39] "Voyriella_parviflora_symbiont_type_1_Richness"
## [40] "Para1_OTU2_Richness"
## [41] "Whitfield_type_3_Richness"
## [42] "Aca_Richness"
## [43] "MO_G20_Richness"
## [44] "Glo7_Richness"
## [45] "MO_GC1_Richness"
## [46] "Winther07_B_Richness"
## [47] "mosseae_Richness"
## [48] "caledonium_Richness"
## [49] "brasilianum_Richness"
## [50] "Alguacil09b_Glo_G8_Richness"
## [51] "Schechter08_Arch1_Richness"
## [52] "Wirsel_OTU21_Richness"
## [53] "Glo71_Richness"
## [54] "spurca_Richness"
## [55] "Glo59_Richness"
## [56] "Wirsel_OTU16_Richness"
## [57] "Glo32_Richness"
## [58] "MO_A8_Richness"
## [59] "MO_G27_Richness"
## [60] "MO_G5_Richness"
## [61] "Glomussp_Richness"
## [62] "Archaeosporasp_Richness"
## [63] "Diversisporasp_Richness"
## [64] "Paraglomussp_Richness"
## [65] "NES27_Read_Abundance"
## [66] "MH3_B_Read_Abundance"
## [67] "Acau16_Read_Abundance"
## [68] "MO_G47_Read_Abundance"
## [69] "laccatum_Read_Abundance"
## [70] "acnaGlo2_Read_Abundance"
## [71] "MO_Ar1"
## [72] "leptoticha_Read_Abundance"
## [73] "Whitfield_type_7_Read_Abundance"
## [74] "Torrecillas12b_Glo_G5_Read_Abundance"
## [75] "Sanchez_Castro12b_GLO12_Read_Abundance"
## [76] "Winther07_D_Read_Abundance"
## [77] "MO_G18_Read_Abundance"
## [78] "intraradices_Read_Abundance"
## [79] "Alguacil12a_Para_1_Read_Abundance"
## [80] "Douhan9_Read_Abundance"
## [81] "Alguacil12b_GLO_G3_Read_Abundance"
## [82] "Glo39_Read_Abundance"
## [83] "Alguacil10_Glo1_Read_Abundance"
## [84] "Ligrone07_sp_Read_Abundance"
## [85] "PT6_Read_Abundance"
## [86] "MO_G41_Read_Abundance"
## [87] "A1_Read_Abundance"
## [88] "Liu2012b_Phylo_5_Read_Abundance"
## [89] "Yamato09_C1_Read_Abundance"
## [90] "MO_G7_Read_Abundance"
## [91] "Alguacil10_Glo6_Read_Abundance"
## [92] "Alguacil12b_PARA1_Read_Abundance"
## [93] "Glo58_Read_Abundance"
## [94] "lamellosum_Read_Abundance"
## [95] "VeGlo18_Read_Abundance"
## [96] "Liu2012a_Ar_1_Read_Abundance"
## [97] "Div_Read_Abundance"
## [98] "fennica_Read_Abundance"
## [99] "Voyriella_parviflora_symbiont_type_1_Read_Abundance"
## [100] "Para1_OTU2_Read_Abundance"
## [101] "Whitfield_type_3_Read_Abundance"
## [102] "Aca_Read_Abundance"
## [103] "MO_G20_Read_Abundance"
## [104] "Glo7_Read_Abundance"
## [105] "MO_GC1_Read_Abundance"
## [106] "Winther07_B7_Read_Abundance"
## [107] "mosseae_Read_Abundance"
## [108] "caledonium_Read_Abundance"
## [109] "brasilianum_Read_Abundance"
## [110] "Alguacil09b_Glo_G8_Read_Abundance"
## [111] "Schechter08_Arch1_Read_Abundance"
## [112] "Wirsel_OTU21_Read_Abundance"
## [113] "Glo71_Read_Abundance"
## [114] "spurca_Read_Abundance"
## [115] "Glo59_Read_Abundance"
## [116] "Wirsel_OTU16_Read_Abundance"
## [117] "Glo32_Read_Abundance"
## [118] "MO_A8_Read_Abundance"
## [119] "MO_G27_Read_Abundance"
## [120] "MO_G5_Read_Abundance"
## [121] "Glomussp_Read_Abundance"
## [122] "Archaeosporasp_Read_Abundance"
## [123] "Diversisporasp_Read_Abundance"
## [124] "Paraglomussp_Read_Abundance"
## [125] "Read_Abundance_Sample"
## [126] "#SampleID"
## [127] "BarcodeSequence"
## [128] "LinkerPrimerSequence"
## [129] "Location"
## [130] "Project"
## [131] "Year"
## [132] "Site"
## [133] "Treat"
## [134] "Rep"
## [135] "TSDel"
## [136] "Description"Create the comm.grps from ssuspeciesREADRICHlong
To make the comm.grps I want the first column (col 1), then I don’t want (want to de-select) cols 2 through 22, and then I want the rest of the data = cols 23-30
This didn’t work:
comm.grps<-ssuspeciesREADRICHlong%>%dplyr::select(‘#SampleID’:Description) #mapping data #heads[1,-c(12:17)] ##This will not work because it uses numeric syntax and the text, only use one or the other (only select; or only numeric -c) comm.grps<-ssuspeciesREADRICHlong%>%dplyr::select[‘ID’:Description,-c(2:22)] #mapping data
#So I used the longhand way of typing out all the column names that I wanted to keep
`?`(dplyr::select)
comm.grps <- ssuspeciesREADRICHlong %>% dplyr::select("ID", "Location", "Year",
"Site", "Treat", "Rep", "TSDel", "Description") #mapping data
colnames(comm.grps)## [1] "ID" "Location" "Year" "Site" "Treat"
## [6] "Rep" "TSDel" "Description"WORK ON THIS!!!
Make the comm.mat ##Make the changes for this for species for the next iteration
comm.mat <- ssuspeciesREADRICHlong %>% dplyr::select(NES27_Richness:Paraglomussp_Richness) # community matrix - all but mapping datacomparing ecological communities
diversity vs composition
abundance and richness are univariate response variables used to quantify communities
in multivariate analyses we have these variables for multiple entities
similarly, multivariate analyses have counterparts in univariate stats - t-test, ANOVA, mutliple regression
univariate analyses of diversity
list.files() #shows what is in folder## [1] "2018Sep19_ssu_funcguild_read_rich.csv"
## [2] "2018Sep19_ssu_genera_read_rich.csv"
## [3] "2018Sep19_ssu_genus_read_rich.csv"
## [4] "2018Sep19_ssu_species_read_rich.csv"
## [5] "AMF.Rmd"
## [6] "AMF2.Rmd"
## [7] "AMF4.Rmd"
## [8] "AMF5.1.Rmd"
## [9] "AMFModeling.Rmd"
## [10] "AMFModelingRevise180913.Rmd"
## [11] "AMFModelingSpecies_1.Rmd"
## [12] "AMFModelingTaxa_3.Rmd"
## [13] "AMFModelingTaxa1 2.Rmd"
## [14] "AMFModelingTaxa1.Rmd"
## [15] "AMFModelingTaxa2.Rmd"
## [16] "AMFSubsetRevise180913.Rmd"
## [17] "bray_curtis_nmds_converted.txt"
## [18] "CN.csv"
## [19] "CSS_table_sorted.txt"
## [20] "HL.xlsx"
## [21] "HyphalSlidesDataInput_LD.xlsx"
## [22] "IdentifyingIssues"
## [23] "map_SSU_Salvage_qiimeformat.txt"
## [24] "mapSal.txt"
## [25] "SalvageHyphalExtractionSlidesMM.calc_in process.xlsx"
## [26] "SalvageProjectDescription"
## [27] "SSU_species.rtf"
## [28] "Test_cache"
## [29] "Test_files"
## [30] "Test.html"
## [31] "Test.Rmd"
## [32] "TestSpecies_cache"
## [33] "TestSpecies.Rmd"
## [34] "TypesOfAMF"# View(comm.mat)
str(comm.grps)## 'data.frame': 61 obs. of 8 variables:
## $ ID : Factor w/ 61 levels "MM933M","MM936M",..: 1 5 2 3 4 6 7 8 9 10 ...
## $ Location : chr "E1" "D1" "C1" "A1" ...
## $ 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" ...head(comm.grps)## ID Location Year Site Treat Rep TSDel Description
## 1 MM933M E1 2015 WL O 5 <NA> RecipientPre
## 2 MM9362CU D1 2015 WL O 4 <NA> RecipientPre
## 3 MM936M C1 2015 WL O 3 <NA> RecipientPre
## 4 MM941M A1 2015 WL O 1 <NA> RecipientPre
## 5 MM943M B1 2015 WL O 2 <NA> RecipientPre
## 6 MMS700 G6 2017 WL D 1 B1A2C RecipientPoststr(comm.mat)## 'data.frame': 61 obs. of 62 variables:
## $ NES27_Richness : int 9 11 6 5 6 11 13 10 15 14 ...
## $ MH3_B_Richness : int 4 4 4 4 3 4 4 4 2 2 ...
## $ Acau16_Richness : int 1 1 1 1 1 1 2 2 0 1 ...
## $ MO_G47_Richness : int 0 0 0 0 0 0 0 0 1 0 ...
## $ laccatum_Richness : int 12 17 6 7 8 12 18 20 11 17 ...
## $ acnaGlo2_Richness : int 11 15 7 10 3 11 10 11 10 12 ...
## $ MO_Ar1_Richness : int 1 2 3 0 1 2 2 3 2 2 ...
## $ leptoticha_Richness : int 3 3 2 1 2 2 2 2 4 4 ...
## $ Whitfield_type_7_Richness : int 3 2 2 3 2 2 2 2 2 3 ...
## $ Torrecillas12b_Glo_G5_Richness : int 9 12 8 9 7 11 13 14 10 11 ...
## $ Sanchez_Castro12b_GLO12_Richness : int 9 14 7 12 9 16 19 17 14 14 ...
## $ Winther07_D_Richness : int 1 3 2 3 2 2 3 1 2 3 ...
## $ MO_G18_Richness : int 0 0 0 0 0 0 0 0 1 1 ...
## $ intraradices_Richness : int 0 1 1 1 0 1 1 1 1 0 ...
## $ Alguacil12a_Para_1_Richness : int 4 5 3 3 3 3 5 4 2 5 ...
## $ Douhan9_Richness : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Alguacil12b_GLO_G3_Richness : int 1 1 0 0 0 0 0 2 3 1 ...
## $ Glo39_Richness : int 3 4 4 4 2 2 2 1 4 4 ...
## $ Alguacil10_Glo1_Richness : int 1 1 1 1 1 1 1 0 1 1 ...
## $ Ligrone07_sp_Richness : int 0 0 0 1 0 0 0 1 0 0 ...
## $ PT6_Richness : int 2 2 2 2 3 1 0 2 1 1 ...
## $ MO_G41_Richness : int 1 0 0 1 0 1 1 0 0 0 ...
## $ MO_GB1_Richness : int 0 0 0 0 0 0 0 0 1 1 ...
## $ A1_Richness : int 1 0 0 1 0 0 0 0 0 0 ...
## $ Liu2012b_Phylo_5_Richness : int 0 1 0 0 0 0 0 0 0 0 ...
## $ Yamato09_C1_Richness : int 1 3 0 0 0 3 2 2 3 2 ...
## $ acnaGlo7_Richness : int 0 0 0 0 0 1 0 0 0 0 ...
## $ MO_G7_Richness : int 0 0 0 0 0 0 0 1 1 0 ...
## $ Alguacil10_Glo6_Richness : int 1 1 1 1 1 0 0 1 2 0 ...
## $ Alguacil12b_PARA1_Richness : int 2 3 0 0 1 3 3 3 3 2 ...
## $ Glo58_Richness : int 1 1 0 0 1 0 0 0 0 0 ...
## $ lamellosum_Richness : int 3 4 4 4 2 6 3 4 5 4 ...
## $ VeGlo18_Richness : int 0 1 0 1 0 0 0 1 1 0 ...
## $ Liu2012a_Ar_1_Richness : int 1 1 1 1 0 1 0 2 1 0 ...
## $ Div_Richness : int 1 0 0 0 2 0 0 2 1 0 ...
## $ fennica_Richness : int 1 3 2 1 1 2 3 3 3 2 ...
## $ Voyriella_parviflora_symbiont_type_1_Richness: int 0 1 0 1 0 1 1 1 0 0 ...
## $ Para1_OTU2_Richness : int 1 1 1 1 0 0 0 0 0 0 ...
## $ Whitfield_type_3_Richness : int 0 1 1 1 0 1 0 0 1 0 ...
## $ Aca_Richness : int 0 0 0 1 0 0 0 0 0 0 ...
## $ MO_G20_Richness : int 1 1 0 1 0 1 0 0 0 0 ...
## $ Glo7_Richness : int 0 0 0 1 0 0 0 0 0 0 ...
## $ MO_GC1_Richness : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Winther07_B_Richness : int 0 0 1 1 1 0 0 0 0 1 ...
## $ mosseae_Richness : int 1 1 1 1 1 1 1 1 1 1 ...
## $ caledonium_Richness : int 1 1 1 1 1 1 1 1 0 2 ...
## $ brasilianum_Richness : int 0 1 0 0 0 1 1 1 0 1 ...
## $ Alguacil09b_Glo_G8_Richness : int 1 0 1 0 0 0 1 1 1 0 ...
## $ Schechter08_Arch1_Richness : int 0 0 1 0 0 0 0 0 0 0 ...
## $ Wirsel_OTU21_Richness : int 0 1 0 0 0 0 0 0 0 1 ...
## $ Glo71_Richness : int 0 0 0 0 0 0 1 0 0 0 ...
## $ spurca_Richness : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Glo59_Richness : int 0 1 0 0 0 0 0 0 0 1 ...
## $ Wirsel_OTU16_Richness : int 1 1 1 1 0 1 1 1 1 0 ...
## $ Glo32_Richness : int 0 1 1 1 1 1 0 1 1 1 ...
## $ MO_A8_Richness : int 0 0 0 1 1 0 1 0 0 0 ...
## $ MO_G27_Richness : int 1 1 1 1 1 1 1 1 1 1 ...
## $ MO_G5_Richness : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Glomussp_Richness : int 13 15 15 15 16 20 20 19 12 14 ...
## $ Archaeosporasp_Richness : int 16 13 8 9 10 13 11 10 6 10 ...
## $ Diversisporasp_Richness : int 1 3 1 1 1 3 2 2 1 1 ...
## $ Paraglomussp_Richness : int 3 4 2 2 2 4 4 4 2 3 ...head(comm.mat)## NES27_Richness MH3_B_Richness Acau16_Richness MO_G47_Richness
## 1 9 4 1 0
## 2 11 4 1 0
## 3 6 4 1 0
## 4 5 4 1 0
## 5 6 3 1 0
## 6 11 4 1 0
## laccatum_Richness acnaGlo2_Richness MO_Ar1_Richness leptoticha_Richness
## 1 12 11 1 3
## 2 17 15 2 3
## 3 6 7 3 2
## 4 7 10 0 1
## 5 8 3 1 2
## 6 12 11 2 2
## Whitfield_type_7_Richness Torrecillas12b_Glo_G5_Richness
## 1 3 9
## 2 2 12
## 3 2 8
## 4 3 9
## 5 2 7
## 6 2 11
## Sanchez_Castro12b_GLO12_Richness Winther07_D_Richness MO_G18_Richness
## 1 9 1 0
## 2 14 3 0
## 3 7 2 0
## 4 12 3 0
## 5 9 2 0
## 6 16 2 0
## intraradices_Richness Alguacil12a_Para_1_Richness Douhan9_Richness
## 1 0 4 1
## 2 1 5 1
## 3 1 3 1
## 4 1 3 1
## 5 0 3 1
## 6 1 3 1
## Alguacil12b_GLO_G3_Richness Glo39_Richness Alguacil10_Glo1_Richness
## 1 1 3 1
## 2 1 4 1
## 3 0 4 1
## 4 0 4 1
## 5 0 2 1
## 6 0 2 1
## Ligrone07_sp_Richness PT6_Richness MO_G41_Richness MO_GB1_Richness
## 1 0 2 1 0
## 2 0 2 0 0
## 3 0 2 0 0
## 4 1 2 1 0
## 5 0 3 0 0
## 6 0 1 1 0
## A1_Richness Liu2012b_Phylo_5_Richness Yamato09_C1_Richness
## 1 1 0 1
## 2 0 1 3
## 3 0 0 0
## 4 1 0 0
## 5 0 0 0
## 6 0 0 3
## acnaGlo7_Richness MO_G7_Richness Alguacil10_Glo6_Richness
## 1 0 0 1
## 2 0 0 1
## 3 0 0 1
## 4 0 0 1
## 5 0 0 1
## 6 1 0 0
## Alguacil12b_PARA1_Richness Glo58_Richness lamellosum_Richness
## 1 2 1 3
## 2 3 1 4
## 3 0 0 4
## 4 0 0 4
## 5 1 1 2
## 6 3 0 6
## VeGlo18_Richness Liu2012a_Ar_1_Richness Div_Richness fennica_Richness
## 1 0 1 1 1
## 2 1 1 0 3
## 3 0 1 0 2
## 4 1 1 0 1
## 5 0 0 2 1
## 6 0 1 0 2
## Voyriella_parviflora_symbiont_type_1_Richness Para1_OTU2_Richness
## 1 0 1
## 2 1 1
## 3 0 1
## 4 1 1
## 5 0 0
## 6 1 0
## Whitfield_type_3_Richness Aca_Richness MO_G20_Richness Glo7_Richness
## 1 0 0 1 0
## 2 1 0 1 0
## 3 1 0 0 0
## 4 1 1 1 1
## 5 0 0 0 0
## 6 1 0 1 0
## MO_GC1_Richness Winther07_B_Richness mosseae_Richness
## 1 0 0 1
## 2 0 0 1
## 3 0 1 1
## 4 0 1 1
## 5 0 1 1
## 6 0 0 1
## caledonium_Richness brasilianum_Richness Alguacil09b_Glo_G8_Richness
## 1 1 0 1
## 2 1 1 0
## 3 1 0 1
## 4 1 0 0
## 5 1 0 0
## 6 1 1 0
## Schechter08_Arch1_Richness Wirsel_OTU21_Richness Glo71_Richness
## 1 0 0 0
## 2 0 1 0
## 3 1 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 0 0
## spurca_Richness Glo59_Richness Wirsel_OTU16_Richness Glo32_Richness
## 1 0 0 1 0
## 2 0 1 1 1
## 3 0 0 1 1
## 4 0 0 1 1
## 5 0 0 0 1
## 6 0 0 1 1
## MO_A8_Richness MO_G27_Richness MO_G5_Richness Glomussp_Richness
## 1 0 1 1 13
## 2 0 1 1 15
## 3 0 1 1 15
## 4 1 1 1 15
## 5 1 1 1 16
## 6 0 1 1 20
## Archaeosporasp_Richness Diversisporasp_Richness Paraglomussp_Richness
## 1 16 1 3
## 2 13 3 4
## 3 8 1 2
## 4 9 1 2
## 5 10 1 2
## 6 13 3 4does diversity vary across groups?
compute diversity indices
indices <- comm.grps
indices\(richness <- rowSums(comm.mat >0) #indices\)shannon <- diversity(comm.mat, index=‘shannon’)
indices$rarified <- c(rarefy(comm.mat, sample=50)) # rarefied diversity for a given sample size
indices <- comm.grps
indices$richness <- rowSums(comm.mat > 0)
indices$shannon <- diversity(comm.mat, index = "shannon")
# View(indices) indices$rarified <- c(rarefy(comm.mat, sample=18103)) #
# rarefied diversity for a given sample sizeWORK ON THIS!!!
Depending on what I put first in the model (Treat or Description) The terms are either significant (p-value) or not-significant.
Subset the data
prevpost.mat = comm.mat[comm.grps$Site == "WL", ]
prevpost.grps = comm.grps[comm.grps$Site == "WL", ]
distWL.bc <- vegdist(prevpost.mat, method = "bray", na.rm = TRUE)
WLTreat.div <- adonis2(distWL.bc ~ Description, strata = Treat, data = prevpost.grps,
permutations = 999, method = "bray")
WLTreat.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = distWL.bc ~ Description, data = prevpost.grps, permutations = 999, method = "bray", strata = Treat)
## Df SumOfSqs R2 F Pr(>F)
## Description 1 0.084955 0.32949 4.914 0.009 **
## Residual 10 0.172884 0.67051
## Total 11 0.257839 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1dvpostWL.mat = comm.mat[(comm.grps$Site == "WL" & comm.grps$Description == "RecipientPost") |
comm.grps$Site == "DS", ]
dvpostWL.grps = comm.grps[(comm.grps$Site == "WL" & comm.grps$Description ==
"RecipientPost") | comm.grps$Site == "DS", ]
distDWL.bc <- vegdist(dvpostWL.mat, method = "bray", na.rm = TRUE)
DWLTreat.div <- adonis2(distDWL.bc ~ Description, strata = Treat, data = dvpostWL.grps,
permutations = 999, method = "bray")
DWLTreat.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = distDWL.bc ~ Description, data = dvpostWL.grps, permutations = 999, method = "bray", strata = Treat)
## Df SumOfSqs R2 F Pr(>F)
## Description 1 0.031320 0.27535 3.4197 0.005 **
## Residual 9 0.082428 0.72465
## Total 10 0.113748 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1prevpostPS.mat = comm.mat[comm.grps$Site == "PS", ]
prevpostPS.grps = comm.grps[comm.grps$Site == "PS", ]
distPS.bc <- vegdist(prevpostPS.mat, method = "bray", na.rm = TRUE)
PSTreat.div <- adonis2(distPS.bc ~ Description, strata = Treat, data = prevpostPS.grps,
permutations = 999, method = "bray")
PSTreat.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = distPS.bc ~ Description, data = prevpostPS.grps, permutations = 999, method = "bray", strata = Treat)
## Df SumOfSqs R2 F Pr(>F)
## Description 1 0.02512 0.10627 2.0215 0.025 *
## Residual 17 0.21125 0.89373
## Total 18 0.23637 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1dim(prevpostPS.grps)## [1] 19 8View(prevpostPS.grps)
dvpostPS.mat = comm.mat[(comm.grps$Site == "PS" & comm.grps$Description == "RecipientPost") |
comm.grps$Site == "DS", ]
dvpostPS.grps = comm.grps[(comm.grps$Site == "PS" & comm.grps$Description ==
"RecipientPost") | comm.grps$Site == "DS", ]
distDPS.bc <- vegdist(dvpostPS.mat, method = "bray", na.rm = TRUE)
DPSTreat.div <- adonis2(distDPS.bc ~ Description, strata = Treat, data = dvpostPS.grps,
permutations = 999, method = "bray")
DPSTreat.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = distDPS.bc ~ Description, data = dvpostPS.grps, permutations = 999, method = "bray", strata = Treat)
## Df SumOfSqs R2 F Pr(>F)
## Description 1 0.01833 0.07542 1.4682 0.147
## Residual 18 0.22472 0.92458
## Total 19 0.24305 1.00000dim(dvpostPS.grps)## [1] 20 8postTreatPS.mat = comm.mat[comm.grps$Site == "PS", ]
postTreatPS.grps = comm.grps[comm.grps$Site == "PS", ]
distTPS.bc <- vegdist(postTreatPS.mat, method = "bray", na.rm = TRUE)
# View(postTreatPS.grps)
contr.mineCD <- function(...) cbind(c(-1, 1, 0, 0, 0, 0))
PSTreatCD <- adonis2(distTPS.bc ~ Treat, data = postTreatPS.grps, permutations = 999,
method = "bray", contr.unordered = "contr.mineCD")
PSTreatCD## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = distTPS.bc ~ Treat, data = postTreatPS.grps, permutations = 999, method = "bray", contr.unordered = "contr.mineCD")
## Df SumOfSqs R2 F Pr(>F)
## Treat 5 0.084771 0.35864 1.4539 0.036 *
## Residual 13 0.151595 0.64136
## Total 18 0.236366 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1contr.mineCT <- function(...) cbind(c(-1, 0, 1, 0, 0, 0))
PSTreatCT <- adonis2(distTPS.bc ~ Treat, data = postTreatPS.grps, permutations = 999,
method = "bray", contr.unordered = "contr.mineCT")
PSTreatCT## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = distTPS.bc ~ Treat, data = postTreatPS.grps, permutations = 999, method = "bray", contr.unordered = "contr.mineCT")
## Df SumOfSqs R2 F Pr(>F)
## Treat 5 0.084771 0.35864 1.4539 0.051 .
## Residual 13 0.151595 0.64136
## Total 18 0.236366 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1contr.mineOC <- function(...) cbind(c(-1, 0, 0, 0, 0, 1))
PSTreatOC <- adonis2(distTPS.bc ~ Treat, data = postTreatPS.grps, permutations = 999,
method = "bray", contr.unordered = "contr.mineOC")
PSTreatOC## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = distTPS.bc ~ Treat, data = postTreatPS.grps, permutations = 999, method = "bray", contr.unordered = "contr.mineOC")
## Df SumOfSqs R2 F Pr(>F)
## Treat 5 0.084771 0.35864 1.4539 0.037 *
## Residual 13 0.151595 0.64136
## Total 18 0.236366 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1postTreatPS.div <- adonis2(distPS.bc ~ Description, strata = Treat, data = postTreatPS.grps,
permutations = 999, method = "bray")
postTreatPS.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = distPS.bc ~ Description, data = postTreatPS.grps, permutations = 999, method = "bray", strata = Treat)
## Df SumOfSqs R2 F Pr(>F)
## Description 1 0.02512 0.10627 2.0215 0.032 *
## Residual 17 0.21125 0.89373
## Total 18 0.23637 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1dim(postTreatPS.grps)## [1] 19 8# View(postTreatPS.grps)
dvpostPS.mat = comm.mat[(comm.grps$Site == "PS" & comm.grps$Description == "RecipientPost") |
comm.grps$Site == "DS", ]
dvpostPS.grps = comm.grps[(comm.grps$Site == "PS" & comm.grps$Description ==
"RecipientPost") | comm.grps$Site == "DS", ]
distDPS.bc <- vegdist(dvpostPS.mat, method = "bray", na.rm = TRUE)
DPSTreat.div <- adonis2(distDPS.bc ~ Description, strata = Treat, data = dvpostPS.grps,
permutations = 999, method = "bray")
DPSTreat.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = distDPS.bc ~ Description, data = dvpostPS.grps, permutations = 999, method = "bray", strata = Treat)
## Df SumOfSqs R2 F Pr(>F)
## Description 1 0.01833 0.07542 1.4682 0.13
## Residual 18 0.22472 0.92458
## Total 19 0.24305 1.00000dim(dvpostPS.grps)## [1] 20 8set.seed(304)
`?`(vegdist)
## are the results the same with other (non evolutionary) dissimiarlity
## indices? dist.j<-vegdist(ssuspeciesREADRICHlong[,-1], method='jaccard')
## dist.bc<-vegdist(ssuspeciesREADRICHlong[,-1], method='bray', na.rm=TRUE)
## dist.j<-vegdist(comm.mat, method='jaccard')
dist.bc <- vegdist(comm.mat, method = "bray", na.rm = TRUE)
AMTreat.div <- adonis2(dist.bc ~ Treat, data = comm.grps, permutations = 999,
method = "bray")
AMTreat.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.bc ~ Treat, data = comm.grps, permutations = 999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Treat 5 0.15520 0.15878 2.0762 0.003 **
## Residual 55 0.82226 0.84122
## Total 60 0.97746 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AMTD.div <- adonis2(dist.bc ~ Treat + Description, data = comm.grps, permutations = 999,
method = "bray")
AMTD.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.bc ~ Treat + Description, data = comm.grps, permutations = 999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Treat 5 0.15520 0.15878 2.0814 0.002 **
## Description 1 0.01696 0.01735 1.1374 0.305
## Residual 54 0.80530 0.82387
## Total 60 0.97746 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AMTSD.div <- adonis2(dist.bc ~ Treat + Site * Description, data = comm.grps,
permutations = 999, method = "bray")
AMTSD.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.bc ~ Treat + Site * Description, data = comm.grps, permutations = 999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Treat 5 0.15520 0.15878 2.3932 0.001 ***
## Site 3 0.11468 0.11732 2.9473 0.001 ***
## Site:Description 2 0.05910 0.06046 2.2784 0.010 **
## Residual 50 0.64849 0.66344
## Total 60 0.97746 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# AMD_ST_J.div<-adonis2(dist.bc~Description+Site*Treat, data=comm.grps,
# permutations = 999, method='bray') AMD_ST_J.div
AMD_ST_BC.div <- adonis2(dist.bc ~ Description + Site + Treat, data = comm.grps,
permutations = 999, method = "bray")
AMD_ST_BC.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.bc ~ Description + Site + Treat, data = comm.grps, permutations = 999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Description 2 0.11225 0.11484 4.1245 0.001 ***
## Site 2 0.09636 0.09858 3.5408 0.001 ***
## Treat 4 0.06127 0.06268 1.1256 0.290
## Residual 52 0.70759 0.72390
## Total 60 0.97746 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Depending on what comes first in the model, we see different significance of these factors
AMTreat.div <- adonis2(dist.bc ~ Treat * Description, data = comm.grps, permutations = 999,
method = "bray")
AMTreat.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.bc ~ Treat * Description, data = comm.grps, permutations = 999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Treat 5 0.15520 0.15878 2.0814 0.002 **
## Description 1 0.01696 0.01735 1.1374 0.297
## Residual 54 0.80530 0.82387
## Total 60 0.97746 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1amfMDS<-metaMDS(comm.mat, distance="bray", k=2, trymax=35, autotransform=TRUE) ##k is the number of dimensions## Wisconsin double standardization
## Run 0 stress 0.2319388
## Run 1 stress 0.2308257
## ... New best solution
## ... Procrustes: rmse 0.05042144 max resid 0.3599927
## Run 2 stress 0.2367229
## Run 3 stress 0.2434353
## Run 4 stress 0.2471315
## Run 5 stress 0.243741
## Run 6 stress 0.2433168
## Run 7 stress 0.2409587
## Run 8 stress 0.2371608
## Run 9 stress 0.2441726
## Run 10 stress 0.2367508
## Run 11 stress 0.2422019
## Run 12 stress 0.2428477
## Run 13 stress 0.2310441
## ... Procrustes: rmse 0.06038808 max resid 0.3893195
## Run 14 stress 0.2374624
## Run 15 stress 0.248226
## Run 16 stress 0.2389516
## Run 17 stress 0.2405907
## Run 18 stress 0.2395373
## Run 19 stress 0.2463107
## Run 20 stress 0.2435754
## Run 21 stress 0.2386304
## Run 22 stress 0.2453639
## Run 23 stress 0.2386625
## Run 24 stress 0.2390593
## Run 25 stress 0.2304344
## ... New best solution
## ... Procrustes: rmse 0.02334888 max resid 0.1714081
## Run 26 stress 0.2372084
## Run 27 stress 0.2361776
## Run 28 stress 0.2401444
## Run 29 stress 0.2367545
## Run 30 stress 0.2354733
## Run 31 stress 0.2334984
## Run 32 stress 0.2353978
## Run 33 stress 0.2431827
## Run 34 stress 0.239765
## Run 35 stress 0.2462569
## *** No convergence -- monoMDS stopping criteria:
## 35: stress ratio > sratmaxamfMDS ##metaMDS takes eaither a distance matrix or your community matrix (then requires method for 'distance=')##
## Call:
## metaMDS(comm = comm.mat, distance = "bray", k = 2, trymax = 35, autotransform = TRUE)
##
## global Multidimensional Scaling using monoMDS
##
## Data: wisconsin(comm.mat)
## Distance: bray
##
## Dimensions: 2
## Stress: 0.2304344
## Stress type 1, weak ties
## No convergent solutions - best solution after 35 tries
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'wisconsin(comm.mat)'stressplot(amfMDS)
#install.packages('ggplot2') ##plotting package
library(ggplot2)
NMDS1 <- amfMDS$points[,1] ##also found using scores(amfMDS)
NMDS2 <- amfMDS$points[,2]
?cbind
amf.plot<-cbind(comm.grps, NMDS1, NMDS2, comm.mat)
p<-ggplot(amf.plot, aes(NMDS1, NMDS2, color=Year))+
geom_point(position=position_jitter(.1), shape=3)+##separates overlapping points
stat_ellipse(type='t',size =1)+ ##draws 95% confidence interval ellipses
theme_minimal()
p
plot<-ggplot(amf.plot, aes(NMDS1, NMDS2, color=Treat))+
stat_ellipse(type='t',size =1)+
theme_minimal()+geom_text(data=amf.plot,aes(NMDS1, NMDS2, label=Site), position=position_jitter(.35))+
annotate("text", x=min(NMDS1), y=min(NMDS2), label=paste('Stress =',round(amfMDS$stress,3))) #add stress to plot
plot
Subset the community data matrix Number of samples that are left
summary(comm.grps)## ID Location Year Site
## MM933M : 1 Length:61 Length:61 Length:61
## MM936M : 1 Class :character Class :character Class :character
## MM941M : 1 Mode :character Mode :character Mode :character
## MM943M : 1
## MM9362CU: 1
## MMS700 : 1
## (Other) :55
## Treat Rep TSDel Description
## C: 9 1:20 Length:61 Length:61
## D: 9 2:19 Class :character Class :character
## F: 9 3:16 Mode :character Mode :character
## O:17 4: 3
## S: 8 5: 2
## T: 9 6: 1
## class(comm.mat)## [1] "data.frame"# subset the comm.mat and comm.grps x=comm.mat[comm.grps$Year == '2015',]
# y=comm.grps[comm.grps$Year == '2015',]commented out some difficult lines
# what is x?
# SSMDS<-metaMDS(x, distance='bray', k = 2, trymax = 35, autotransform =
# TRUE) ##k is the number of dimensions SSMDS ##metaMDS takes eaither a
# distance matrix or your community matrix (then requires method for
# 'distance=')
# stressplot(SSMDS)
# install.packages('ggplot2') ##plotting package
library(ggplot2)
# NMDS1 <- SSMDS$points[,1] ##also found using scores(amfMDS) NMDS2 <-
# SSMDS$points[,2]
`?`(cbind)The code was stopping at this point, and I created a new chunk below
SSamf.plot<-cbind(y, NMDS1, NMDS2, x)
Sp<-ggplot(SSamf.plot, aes(NMDS1, NMDS2, color=Site))+ geom_point(position=position_jitter(.1), shape=7)+##separates overlapping points # stat_ellipse(type=‘t’,size =1)+ ##draws 95% confidence interval ellipses theme_minimal() Sp
SSp<-ggplot(SSamf.plot, aes(NMDS1, NMDS2, color=Site))+ geom_point(position=position_jitter(.1), shape=3)+##separates overlapping points stat_ellipse(type=‘t’,size =1)+ ##draws 95% confidence interval ellipses theme_minimal() SSp
plot<-ggplot(SSamf.plot, aes(NMDS1, NMDS2, color=Treat))+ stat_ellipse(type=‘t’,size =1)+ theme_minimal()+geom_text(data=amf.plot,aes(NMDS1, NMDS2, label=Site), position=position_jitter(.35))+ annotate(“text”, x=min(NMDS1), y=min(NMDS2), label=paste(‘Stress =’,round(amfMDS$stress,3))) #add stress to plot plot
ad.bc<-adonis2(dist.bc~Treat+Metadata1+Metadata, data=grps, permutations=1000)
ad.bc
Examining the statistic support for increased rhizophilic richness in recipientpost plots
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"))
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 "NES27" "MH3_B" "NES27" "sp" ...
## $ 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 ...ssuguildREADRICH <- data.frame(ssul %>% group_by(ID, Functional.Group) %>% summarise(OTU_Richness_Sample = length(unique(ssuotu[ssureads >
0])), Read_Abundance_Sample = sum(ssureads)))
ssufamily <- data.frame(ssul %>% group_by(ID, ssufamily) %>% summarise(Family_OTU_Richness = length(unique(ssuotu[ssureads >
0])), Family_Read_Abundance = sum(ssureads)))
colnames(metassu) # need to edit for appropriate headers## [1] "#SampleID" "BarcodeSequence" "LinkerPrimerSequence"
## [4] "Location" "Project" "Year"
## [7] "Site" "Treat" "Rep"
## [10] "TSDel" "Description" "ID"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)
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.3515ssuguildREADRICHlong <- 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" ...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_Sample Ancestral_Richness Edaphophilic_Richness
## <fct> <int> <int> <int>
## 1 MM933M 130 26 2
## 2 MM936M 105 20 1
## 3 MM941M 121 19 1
## 4 MM943M 99 20 3
## 5 MM9362CU 167 28 3
## 6 MMS700 151 23 3
## # ... with 1 more variable: Rhizophilic_Richness <int>ssuguildRICH <- merge(ssuguildRICH, metassu, by = "ID")
ssuguildREAD <- merge(ssuguildREAD, metassu, by = "ID")
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" ...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")
# LOGNORMAL
qqp(ssuguildREADRICHlong$OTU_Richness_Sample, "lnorm")
# 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]])
# POISSON
poisson <- fitdistr(ssuguildREADRICHlong$OTU_Richness_Sample.t, "Poisson")
qqp(ssuguildREADRICHlong$OTU_Richness_Sample.t, "pois", lambda = poisson$estimate)
# GAMMA
gamma <- fitdistr(ssuguildREADRICHlong$OTU_Richness_Sample.t, "gamma")
qqp(ssuguildREADRICHlong$OTU_Richness_Sample.t, "gamma", shape = gamma$estimate[[1]],
rate = gamma$estimate[[2]])
# 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.4modelingAMF_MP <- glm(OTU_Richness_Sample.t ~ Site + Treat + Description + Rep +
Year, data = ssuguildREADRICHlong, family = gaussian(link = "identity"))
plot(modelingAMF_MP)
stepAIC(modelingAMF_MP)## Start: AIC=533.64
## OTU_Richness_Sample.t ~ Site + Treat + Description + Rep + Year
##
##
## Step: AIC=533.64
## OTU_Richness_Sample.t ~ Site + Treat + Description + Rep
##
##
## Step: AIC=533.64
## OTU_Richness_Sample.t ~ Site + Treat + Rep
##
## Df Deviance AIC
## - Site 3 15052 533.12
## <none> 13759 533.64
## - Rep 5 16673 535.36
## - Treat 5 21667 551.34
##
## Step: AIC=533.12
## OTU_Richness_Sample.t ~ Treat + Rep
##
## Df Deviance AIC
## - Rep 5 17693 532.98
## <none> 15052 533.12
## - Treat 5 24372 552.52
##
## 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: 533modelingAMF_MP##
## Call: glm(formula = OTU_Richness_Sample.t ~ Site + Treat + Description +
## Rep + Year, family = gaussian(link = "identity"), data = ssuguildREADRICHlong)
##
## Coefficients:
## (Intercept) SiteHH
## 173.606 -4.303
## SitePS SiteWL
## -11.545 -15.011
## TreatD TreatF
## -5.124 17.677
## TreatO TreatS
## -22.530 10.321
## TreatT DescriptionRecipientPost
## 3.024 NA
## DescriptionRecipientPre Rep2
## NA -1.802
## Rep3 Rep4
## -10.783 16.776
## Rep5 Rep6
## -7.570 36.936
## Year2017
## NA
##
## Degrees of Freedom: 60 Total (i.e. Null); 47 Residual
## Null Deviance: 26340
## Residual Deviance: 13760 AIC: 533.6stepAIC(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"))Plotting
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 #Tried with negative binomial but it didn’t work
# salvage.div<-adonis2(bird.dist~DIVERSITY, data=birds, permutations = 999,
# method='bray', strata='PLOT') bird.div