-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathCOBS_ITS_Ascomycota.R
341 lines (309 loc) · 16.7 KB
/
COBS_ITS_Ascomycota.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
#Elizabeth Bach
#COBS ITS data 2012
#Analysis of Ascomycota only
# 15 July 2014
#Use Asco.data generated in COBS_ITS_PhylaSubset.R
#Graphs can befound in COBS_ITS_PhylaSubset.R as well
#will need these libraries for mixed models, data manipulation
library(reshape)
library(lme4)
library(lmerTest)
library(bbmle)
#Some orders show distinct SoilFrac and Agg, not much date effect
head(Asco.data)
Asco.null<-lmer(value~1+(1|Block), data=Asco.data, REML=FALSE)
Asco.model.full<-lmer(value~Date*Crop*SoilFrac+(1|Block), data=Asco.data, REML=FALSE)
Asco.model.main<-lmer(value~Date+Crop+SoilFrac+(1|Block), data=Asco.data, REML=FALSE)
AICtab(Asco.null,Asco.model.full,Asco.model.main)
anova(Asco.null,Asco.model.full,Asco.model.main)
#Main model best fit by 5.6
anova(Asco.model.main)
anova(Asco.model.full)
#No interactions in full model, proceed with main
difflsmeans(Asco.model.main, ddf="Satterthwaite",type=3,method.grad="simple")
#Significant effect of Date (Oct>July, P=0.001) and Crop (P=PF>CC, P=0.007)
#unidentified=Orbiliales
Asco.unid<-subset(Asco.data, Asco.data$Order=="o__unidentified")
head(Asco.unid)
str(Asco.unid)
Asco.orbi<-subset(Asco.unid, Asco.unid$Species=="s__Orbiliomycetessp")
str(Asco.orbi)
#All "unidentified order are identified as species Orbiliomycetes sp., without identification of higher taxanomic status
#Orbiliomycetes is class, Orbiliales is order, Orbiliaceae is family
Asco.null<-lmer(value~1+(1|Block), data=Asco.orbi, REML=FALSE)
Asco.model.full<-lmer(value~Date*Crop*SoilFrac+(1|Block), data=Asco.orbi, REML=FALSE)
Asco.model.main<-lmer(value~Date+Crop+SoilFrac+(1|Block), data=Asco.orbi, REML=FALSE)
AICtab(Asco.null,Asco.model.full,Asco.model.main)
anova(Asco.null,Asco.model.full,Asco.model.main)
#Main model best fit by 4.5
anova(Asco.model.main)
anova(Asco.model.full)
#No significant interactions, proceed with main effects model
difflsmeans(Asco.model.main, ddf="Satterthwaite",type=3,method.grad="simple")
#Significant effect of Date: Oct>July (P=0.02)
#Significant effect of Crop: PF>P=CC (P=0.0005)
Asco.orbiC<-ddply(Asco.orbi, .(Crop), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value),N=length(value),SE=(sd(value)/sqrt(N-1))
)
head(Asco.orbiC)
ggplot(Asco.orbiC)+geom_pointrange(aes(x=Species,y=mean,ymax=high95,ymin=low95, color=Crop),position=position_dodge(width=1))+coord_flip()+scale_y_log10())
Asco.orbiD<-ddply(Asco.orbi, .(Date, Species), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value)
)
head(Asco.orbiD)
ggplot(Asco.orbiD)+geom_pointrange(aes(x=Species,y=mean,ymax=high95,ymin=low95, color=Date),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
#Sordariales
Asco.sord<-subset(Asco.data, Asco.data$Order=="o__Sordariales")
head(Asco.sord)
str(Asco.sord)
#Only detected in 2 samples, P46.MM Oct. 2012 and PF41.micro Oct. 2012, very low abundance
#single species, Humicola nigescens
#Note, this order is Andy Miller (new INHS boss) area of expertise
#Pleosporales
Asco.pleo<-subset(Asco.data, Asco.data$Order=="o__Pleosporales")
head(Asco.pleo)
str(Asco.pleo)
Asco.null<-lmer(value~1+(1|Block), data=Asco.pleo, REML=FALSE)
Asco.model.full<-lmer(value~Date*Crop*SoilFrac+(1|Block), data=Asco.pleo, REML=FALSE)
Asco.model.main<-lmer(value~Date+Crop+SoilFrac+(1|Block), data=Asco.pleo, REML=FALSE)
AICtab(Asco.null,Asco.model.full,Asco.model.main)
anova(Asco.null,Asco.model.full,Asco.model.main)
#Null model best fit, by 5, main effects next best fit
anova(Asco.model.main)
anova(Asco.model.full)
#No significant effects
Asco.pleoC<-ddply(Asco.pleo, .(Crop, Genus), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value)
)
head(Asco.pleoC)
ggplot(Asco.pleoC)+geom_pointrange(aes(x=Genus,y=mean,ymax=high95,ymin=low95, color=Crop),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
Asco.pleoD<-ddply(Asco.pleo, .(Date, Genus), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value)
)
head(Asco.pleoD)
ggplot(Asco.pleoD)+geom_pointrange(aes(x=Genus,y=mean,ymax=high95,ymin=low95, color=Date),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
Asco.pleoSF<-ddply(Asco.pleo, .(SoilFrac, Genus), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value)
)
head(Asco.pleoSF)
ggplot(Asco.pleoSF)+geom_pointrange(aes(x=Genus,y=mean,ymax=high95,ymin=low95, color=SoilFrac),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
#Genus Leptosphaerulina is most abundant, but doesn't differ by crop, date, or SoilFrac
#Genus Drechslera only found in 3 CC samples, abundance ranges 16-90, all Drechslera sp. BAFC3419
#Genus Alternaria very low abundance in only found in 2 samples, all Alternaria sp. 3MU_2012
Asco.drec<-subset(Asco.data, Asco.data$Genus=="g__Drechslera")
head(Asco.drec)
Asco.alt<-subset(Asco.data, Asco.data$Genus=="g__Alternaria")
head(Asco.alt)
#Pezizales
Asco.pez<-subset(Asco.data, Asco.data$Order=="o__Pezizales")
head(Asco.pez)
str(Asco.pez)
Asco.null<-lmer(value~1+(1|Block), data=Asco.pez, REML=FALSE)
Asco.model.full<-lmer(value~Date*Crop*SoilFrac+(1|Block), data=Asco.pez, REML=FALSE)
Asco.model.main<-lmer(value~Date+Crop+SoilFrac+(1|Block), data=Asco.pez, REML=FALSE)
AICtab(Asco.null,Asco.model.full,Asco.model.main)
anova(Asco.null,Asco.model.full,Asco.model.main)
#Null model lowest AIC by 1, main effects next best fit
anova(Asco.model.main)
anova(Asco.model.full)
#No significant interactions, so using main effects model
#Significant effect of date, Oct>July P=0.047
#SoilFrac P=0.08
#All samples either Peziza varia or unkn, both show date effect, unkn species in genus Peziza
Asco.pezD<-ddply(Asco.pez, .(Date, Species), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value)
)
head(Asco.pezD)
ggplot(Asco.pezD)+geom_pointrange(aes(x=Species,y=mean,ymax=high95,ymin=low95, color=Date),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
Asco.pezSF<-ddply(Asco.pez, .(SoilFrac, Species), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value)
)
head(Asco.pezSF)
ggplot(Asco.pezSF)+geom_pointrange(aes(x=Species,y=mean,ymax=high95,ymin=low95, color=SoilFrac),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
#Peziza varia is common wood-decomposition cup fungus
#Only found in CC and PF plots, again no P indicating fertilizer effect or fluke? No diff in CC and PF
Asco.pezC<-ddply(Asco.pez, .(Crop), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value),N=length(value),SE=(sd(value)/sqrt(N-1))
)
head(Asco.pezC)
ggplot(Asco.pezC)+geom_pointrange(aes(x=Species,y=mean,ymax=high95,ymin=low95, color=Crop),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
#Incertaesedis
Asco.inc<-subset(Asco.data, Asco.data$Order=="o__Incertaesedis")
head(Asco.inc)
str(Asco.inc)
#15 observations
Asco.null<-lmer(value~1+(1|Block), data=Asco.inc, REML=FALSE)
Asco.model.full<-lmer(value~Date*Crop*SoilFrac+(1|Block), data=Asco.inc, REML=FALSE)
Asco.model.main<-lmer(value~Date+Crop+SoilFrac+(1|Block), data=Asco.inc, REML=FALSE)
AICtab(Asco.null,Asco.model.full,Asco.model.main)
anova(Asco.null,Asco.model.full,Asco.model.main)
#Data suffeciently unbalanced to make contrast models not possible
Asco.incD<-ddply(Asco.inc, .(Date, Species), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value)
)
head(Asco.incD)
ggplot(Asco.incD)+geom_pointrange(aes(x=Species,y=mean,ymax=high95,ymin=low95, color=Date),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
Asco.incSF<-ddply(Asco.inc, .(SoilFrac, Species), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value)
)
head(Asco.incSF)
ggplot(Asco.incSF)+geom_pointrange(aes(x=Species,y=mean,ymax=high95,ymin=low95, color=SoilFrac),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
Asco.incC<-ddply(Asco.inc, .(Crop, Species), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value)
)
head(Asco.incC)
ggplot(Asco.incC)+geom_pointrange(aes(x=Species,y=mean,ymax=high95,ymin=low95, color=Crop),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
#Three species: Vermispora spAS60291 (appears to be plant pathogen?), Scolecobasidium humicola (endophyte, promotes growth in organic tomato production), Polyscytalum pustulans (a potato plant pathogen)
#Also a couple of unkw species
#Very sporadic distribution among samples, no clear pattern, except all species consistently appear in PF, but not necessarily every plot
#Hypocreales
Asco.hypo<-subset(Asco.data, Asco.data$Order=="o__Hypocreales")
head(Asco.hypo)
str(Asco.hypo)
#150 observations!
Asco.null<-lmer(value~1+(1|Block), data=Asco.hypo, REML=FALSE)
Asco.model.full<-lmer(value~Date*Crop*SoilFrac+(1|Block), data=Asco.hypo, REML=FALSE)
Asco.model.main<-lmer(value~Date+Crop+SoilFrac+(1|Block), data=Asco.hypo, REML=FALSE)
AICtab(Asco.null,Asco.model.full,Asco.model.main)
anova(Asco.null,Asco.model.full,Asco.model.main)
#Null model and main model essentially same AIC
anova(Asco.model.main)
anova(Asco.model.full)
#Date*Crop interaction is significant in full model, P=0.03, Crop is significant in main effects model (P=0.01)
#Truncated model to include main factors + Date*Crop interaction
Hypo.model<-lmer(value~Date+Crop+SoilFrac+Date*Crop+(1|Block), data=Asco.hypo, REML=FALSE)
AICtab(Asco.null,Asco.model.full,Asco.model.main, Hypo.model)
#Best AIC fit
anova(Hypo.model)
#Date*Crop interaction significant P=0.04
Asco.hypoCD<-ddply(Asco.hypo, .(Date, Crop, Order), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value)
)
head(Asco.hypoCD)
ggplot(Asco.hypoCD)+geom_pointrange(aes(x=Order,y=mean,ymax=high95,ymin=low95, group=Crop, color=Date),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
#Interaction is from similar abundance in P and PF at both sampling dates, but greater presence in CC in Oct
#Includes three families, including unk, and there are some different patterns within families, so will look more closely at that
Asco.hypoCD<-ddply(Asco.hypo, .(Date, Crop, Family), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value)
)
head(Asco.hypoCD)
ggplot(Asco.hypoCD)+geom_pointrange(aes(x=Family,y=mean,ymax=high95,ymin=low95, group=Crop, color=Date),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
#Family Hypocreaceae
Asco.hypoc<-subset(Asco.data, Asco.data$Family=="f__Hypocreaceae")
head(Asco.hypoc)
str(Asco.hypoc)
#64 observations
Asco.null<-lmer(value~1+(1|Block), data=Asco.hypoc, REML=FALSE)
Asco.model.full<-lmer(value~Date*Crop*SoilFrac+(1|Block), data=Asco.hypoc, REML=FALSE)
Asco.model.main<-lmer(value~Date+Crop+SoilFrac+(1|Block), data=Asco.hypoc, REML=FALSE)
AICtab(Asco.null,Asco.model.full,Asco.model.main)
anova(Asco.null,Asco.model.full,Asco.model.main)
#Null model lowest AIC by 4
anova(Asco.model.main)
anova(Asco.model.full)
#No interactions, Crop significant (P=0.04), PF=CC>P, all observations are Trichoderma citrinoviride, except 1 unk
difflsmeans(Asco.model.main, ddf="Satterthwaite",type=3,method.grad="simple")
Asco.hypoC<-ddply(Asco.hypoc, .(Crop, Species), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value),N=length(value),SE=(sd(value)/sqrt(N-1))
)
head(Asco.hypoC)
ggplot(Asco.hypoC)+geom_pointrange(aes(x=Species,y=mean,ymax=high95,ymin=low95, color=Crop),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
#Family Bionectriaceae
Asco.bion<-subset(Asco.data, Asco.data$Family=="f__Bionectriaceae")
head(Asco.bion)
str(Asco.bion)
#33 observations
Asco.null<-lmer(value~1+(1|Block), data=Asco.bion, REML=FALSE)
Asco.model.full<-lmer(value~Date*Crop*SoilFrac+(1|Block), data=Asco.bion, REML=FALSE)
Asco.model.main<-lmer(value~Date+Crop+SoilFrac+(1|Block), data=Asco.bion, REML=FALSE)
AICtab(Asco.null,Asco.model.full,Asco.model.main)
anova(Asco.null,Asco.model.full,Asco.model.main)
#Main effects model has lowest AIC
anova(Asco.model.main)
anova(Asco.model.full)
#No interactions, SoilFrac highly significant: P=0.006, highly abundant in LM and micro, small presence in MM, only 1 in SM and WS
difflsmeans(Asco.model.main, ddf="Satterthwaite",type=3,method.grad="simple")
Asco.bionSF<-ddply(Asco.bion, .(SoilFrac, Family), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value)
)
head(Asco.bionSF)
ggplot(Asco.bionSF)+geom_pointrange(aes(x=Family,y=mean,ymax=high95,ymin=low95, color=SoilFrac),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
#Species: Clonostachysroseaf.catenulata + unkw, family-level relationship is same for both species
Asco.bionSF<-ddply(Asco.bion, .(SoilFrac, Species), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value)
)
head(Asco.bionSF)
ggplot(Asco.bionSF)+geom_pointrange(aes(x=Species,y=mean,ymax=high95,ymin=low95, color=SoilFrac),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
#unk Family
Asco.unk<-subset(Asco.hypo, Asco.hypo$Family=="")
head(Asco.unk)
str(Asco.unk)
#53 observations, abundances ~10-50
Asco.null<-lmer(value~1+(1|Block), data=Asco.unk, REML=FALSE)
Asco.model.full<-lmer(value~Date*Crop*SoilFrac+(1|Block), data=Asco.unk, REML=FALSE)
Asco.model.main<-lmer(value~Date+Crop+SoilFrac+(1|Block), data=Asco.unk, REML=FALSE)
AICtab(Asco.null,Asco.model.full,Asco.model.main)
anova(Asco.null,Asco.model.full,Asco.model.main)
#Null model is lowest AIC by 3, main model is next best fit
anova(Asco.model.full)
anova(Asco.model.main)
#Full model shows significant interaction of Date*Crop (P=0.004), main effects of Date and Crop both marginal at P=0.06
#Probably not worth interpreting this, just leave at order level
Asco.unkC<-ddply(Asco.unk, .(Crop, Family), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value)
)
head(Asco.unkC)
ggplot(Asco.unkC)+geom_pointrange(aes(x=Family,y=mean,ymax=high95,ymin=low95, color=Crop),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
#Helotiales
Asco.helo<-subset(Asco.data, Asco.data$Order=="o__Helotiales")
head(Asco.helo)
str(Asco.helo)
#13 observations, low abundances, probably not contributing much
Asco.heloD<-ddply(Asco.helo, .(Date, Species), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value)
)
head(Asco.heloD)
ggplot(Asco.heloD)+geom_pointrange(aes(x=Species,y=mean,ymax=high95,ymin=low95, color=Date),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
Asco.heloSF<-ddply(Asco.helo, .(SoilFrac, Species), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value)
)
head(Asco.heloSF)
ggplot(Asco.heloSF)+geom_pointrange(aes(x=Species,y=mean,ymax=high95,ymin=low95, color=SoilFrac),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
Asco.heloC<-ddply(Asco.helo, .(Crop, Species), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value)
)
head(Asco.heloC)
ggplot(Asco.heloC)+geom_pointrange(aes(x=Species,y=mean,ymax=high95,ymin=low95, color=Crop),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
#Generally low abundance, consistently present in Oct., but not July. Not consistent across crop or aggregate fraction
#Capnodiales
Asco.cap<-subset(Asco.data, Asco.data$Order=="o__Capnodiales")
head(Asco.cap)
str(Asco.cap)
#Only 6 observations, all Cladosporium sp4MU_2012, value=1, no consistency in samples, super common air-borne fungus
Asco.cap
#unk
Asco.unk<-subset(Asco.data, Asco.data$Order=="")
head(Asco.unk)
str(Asco.unk)
#113 observations, can't do much interpretation, but looking a patterns just incase something pops out
Asco.null<-lmer(value~1+(1|Block), data=Asco.unk, REML=FALSE)
Asco.model.full<-lmer(value~Date*Crop*SoilFrac+(1|Block), data=Asco.unk, REML=FALSE)
Asco.model.main<-lmer(value~Date+Crop+SoilFrac+(1|Block), data=Asco.unk, REML=FALSE)
AICtab(Asco.null,Asco.model.full,Asco.model.main)
anova(Asco.null,Asco.model.full,Asco.model.main)
#Full model best fit by 5
anova(Asco.model.full)
#Significant Crop*SoilFrac interaction, truncated model
Asco.model.unk<-lmer(value~Date+Crop+SoilFrac+Crop*SoilFrac+(1|Block), data=Asco.unk, REML=FALSE)
AICtab(Asco.null,Asco.model.full,Asco.model.main,Asco.model.unk)
#Much better fit
anova(Asco.model.unk)
difflsmeans(Asco.model.unk, ddf="Satterthwaite",type=3,method.grad="simple")
#Highly significnat Crop*SoilFrac interaction, P<0.0001,
Asco.unkCS<-ddply(Asco.unk, .(Crop, SoilFrac, Family), summarise,.progress="text",mean=mean(value),high95=boot.high(value),
low95=boot.low(value)
)
head(Asco.unkCS)
ggplot(Asco.unkCS)+geom_pointrange(aes(x=Crop,y=mean,ymax=high95,ymin=low95, group=Crop, color=SoilFrac),position=position_dodge(width=1))+coord_flip()+scale_y_log10()
#Each cropping system has unique rank of aggregate fractions, lots of overlap
#Again, no power to pull this apart further, but interesting to note amount of unidentified Ascomycota