-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy path.Rhistory
512 lines (512 loc) · 31.3 KB
/
.Rhistory
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
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
test = merge(basal_richness, data.frame(Subject_Id=row.names(test_combine),test_combine) )
#test = merge(basal_richness,basal_diet[which(basal_diet$visite=="before_study"),], by.x="Subject_Id", by.y="sujet")
#test = merge(basal_richness,basal_diet[which(basal_diet$visite=="after_study"),], by.x="Subject_Id", by.y="sujet")
#test = merge(basal_richness,basal_diet, by.x="Subject_Id", by.y="sujet")
for(i in 5:dim(test)[2]){
test[,i] = as.numeric(as.character(test[,i]))
}
test = test[, c(2,4+which(apply(test[-c(1:4)],2,sum)!=0))]
cor(test$richness, test[,-1], method="spearman")
apply(test[,-1], 2, function(x) cor.test(test$richness, x, method="spearman")$p.value)
diet_test_result = data.frame(rho=t(cor(test$richness, test[,-1], method="spearman")), p=apply(test[,-1], 2, function(x) wilcox.test(test$richness~x)$p.value))
diet_test_result = data.frame(rho=t(cor(test$richness, test[,-1], method="spearman")), p=apply(test[,-1], 2, function(x) cor.test(test$richness,x, method="spearman")$p.value))
diet_test_result[which(diet_test_result$p<0.05),]
ggplot(melt(test, id.vars=c("richness")), aes(fill=as.character(value), y=richness, x=variable)) + geom_boxplot() + coord_flip()
p_ind = ggplot(data.frame(Richness = test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li)) +
coord_tern() + geom_point(aes(
x=(Axis1-min(Axis1))/(max(Axis1)-min(Axis1)),
y=(Axis2-min(Axis2))/(max(Axis2)-min(Axis2)),
z=(Axis3-min(Axis3))/(max(Axis3)-min(Axis3)),
size=log2(Richness))) + theme_rgbw() +
xlab("PC1") + ylab("PC2") + zlab("PC3")
p_diet = ggplot(data.frame(dudi.pca(test[,-1], scannf=F, nf=10)$co)) +
coord_tern() + geom_text(aes(
x=(Comp1-min(Comp1))/(max(Comp1)-min(Comp1)),
y=(Comp2-min(Comp2))/(max(Comp2)-min(Comp2)),
z=(Comp3-min(Comp3))/(max(Comp3)-min(Comp3)),
label=row.names(dudi.pca(test[,-1], scannf=F, nf=10)$co)), size=3) + theme_rgbw() +
xlab("PC1 load.") + ylab("PC2 load.") + zlab("PC3 load.")
p_diet
p_ind_diet = grid.arrange(p_ind, p_diet, ncol=2)
cor.test(test$richness, dudi.pca(test[-1], scannf=F, nf=10)$li[,2], method="spearman")
cor.test(test$richness, dudi.pca(test[-1], scannf=F, nf=10)$li[,1], method="spearman")
cor.test(test$richness, dudi.pca(test[-1], scannf=F, nf=10)$li[,3], method="spearman")
p_ind
test_combine = apply(basal_diet[,-c(1:3)] , 2, function(x){ tapply(as.numeric(x), as.factor(basal_diet[,1]) ,sum) } )
test = merge(basal_richness, data.frame(Subject_Id=row.names(test_combine),test_combine) )
#test = merge(basal_richness,basal_diet[which(basal_diet$visite=="before_study"),], by.x="Subject_Id", by.y="sujet")
#test = merge(basal_richness,basal_diet[which(basal_diet$visite=="after_study"),], by.x="Subject_Id", by.y="sujet")
#test = merge(basal_richness,basal_diet, by.x="Subject_Id", by.y="sujet")
for(i in 5:dim(test)[2]){
test[,i] = as.numeric(as.character(test[,i]))
}
test = test[, c(2,4+which(apply(test[-c(1:4)],2,sum)!=0))]
cor(test$richness, test[,-1], method="spearman")
apply(test[,-1], 2, function(x) cor.test(test$richness, x, method="spearman")$p.value)
diet_test_result = data.frame(rho=t(cor(test$richness, test[,-1], method="spearman")), p=apply(test[,-1], 2, function(x) wilcox.test(test$richness~x)$p.value))
diet_test_result = data.frame(rho=t(cor(test$richness, test[,-1], method="spearman")), p=apply(test[,-1], 2, function(x) cor.test(test$richness,x, method="spearman")$p.value))
diet_test_result[which(diet_test_result$p<0.05),]
ggplot(melt(test, id.vars=c("richness")), aes(fill=as.character(value), y=richness, x=variable)) + geom_boxplot() + coord_flip()
p_ind = ggplot(data.frame(Richness = test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li)) +
coord_tern() + geom_point(aes(
x=(Axis1-min(Axis1))/(max(Axis1)-min(Axis1)),
y=(Axis2-min(Axis2))/(max(Axis2)-min(Axis2)),
z=(Axis3-min(Axis3))/(max(Axis3)-min(Axis3)),
size=log2(Richness))) + theme_rgbw() +
xlab("PC1") + ylab("PC2") + zlab("PC3")
p_diet = ggplot(data.frame(dudi.pca(test[,-1], scannf=F, nf=10)$co)) +
coord_tern() + geom_text(aes(
x=(Comp1-min(Comp1))/(max(Comp1)-min(Comp1)),
y=(Comp2-min(Comp2))/(max(Comp2)-min(Comp2)),
z=(Comp3-min(Comp3))/(max(Comp3)-min(Comp3)),
label=row.names(dudi.pca(test[,-1], scannf=F, nf=10)$co)), size=3) + theme_rgbw() +
xlab("PC1 load.") + ylab("PC2 load.") + zlab("PC3 load.")
p_diet
p_ind_diet = grid.arrange(p_ind, p_diet, ncol=2)
cor.test(test$richness, dudi.pca(test[-1], scannf=F, nf=10)$li[,2], method="spearman")
cor.test(test$richness, dudi.pca(test[-1], scannf=F, nf=10)$li[,1], method="spearman")
cor.test(test$richness, dudi.pca(test[-1], scannf=F, nf=10)$li[,3], method="spearman")
p_ind
p_diet
q()
basal_diet = read.csv2("data-raw/diet_basal_alimintest.csv")
basal_diet$sujet = as.character(basal_diet$sujet)
basal_diet$jour = as.character(basal_diet$jour)
for(i in 4:dim(basal_diet)[2]){
basal_diet[,i] = as.factor(basal_diet[,i])
}
save(basal_diet,file="data/basal_diet.RData")
library(ade4)
library(ggplot2)
basal_diet_tmp = basal_diet[which(basal_diet$visite=="before_study"),]
for(i in 4:dim(basal_diet_tmp)[2]){
basal_diet_tmp[,i] = as.factor(as.character(basal_diet_tmp[,i]))
}
diet.acm = dudi.acm(basal_diet_tmp[,-c(1:3)], nf=3, scannf=FALSE)
diet.bca = bca(diet.acm, fac=as.factor(basal_diet_tmp[,1]), nf=3, scannf=F)
randtest(diet.bca)
basal_richness = metadata76[which(metadata76$Time.point=="2"),c("Subject_Id","richness")]
diet.richness.bca = merge(data.frame(Subject_Id=row.names(diet.bca$li), diet.bca$li), basal_richness, by="Subject_Id")
ggplot(diet.richness.bca, aes(x=Axis1, y=Axis2, col=richness, size=richness)) + geom_point()
diet.richness.bca
basal_diet_tmp[,1]
test_combine = apply(basal_diet[which(basal_diet$visite=="after_study"),-c(1:3)] , 2, function(x){ tapply(as.numeric(x), as.factor(basal_diet[which(basal_diet$visite=="after_study"),1]) ,sum) } )
test_combine = apply(basal_diet[,-c(1:3)] , 2, function(x){ tapply(as.numeric(x), as.factor(basal_diet[,1]) ,sum) } )
test = merge(basal_richness, data.frame(Subject_Id=row.names(test_combine),test_combine) )
#test = merge(basal_richness,basal_diet[which(basal_diet$visite=="before_study"),], by.x="Subject_Id", by.y="sujet")
#test = merge(basal_richness,basal_diet[which(basal_diet$visite=="after_study"),], by.x="Subject_Id", by.y="sujet")
#test = merge(basal_richness,basal_diet, by.x="Subject_Id", by.y="sujet")
for(i in 5:dim(test)[2]){
test[,i] = as.numeric(as.character(test[,i]))
}
test = test[, c(2,4+which(apply(test[-c(1:4)],2,sum)!=0))]
cor(test$richness, test[,-1], method="spearman")
test[,-1]
apply(test[,-1], 2, function(x) cor.test(test$richness, x, method="spearman")$p.value)
diet_test_result = data.frame(rho=t(cor(test$richness, test[,-1], method="spearman")), p=apply(test[,-1], 2, function(x) wilcox.test(test$richness~x)$p.value))
diet_test_result = data.frame(rho=t(cor(test$richness, test[,-1], method="spearman")), p=apply(test[,-1], 2, function(x) cor.test(test$richness,x, method="spearman")$p.value))
diet_test_result[which(diet_test_result$p<0.05),]
ggplot(melt(test, id.vars=c("richness")), aes(fill=as.character(value), y=richness, x=variable)) + geom_boxplot() + coord_flip()
library(reshape2)
ggplot(melt(test, id.vars=c("richness")), aes(fill=as.character(value), y=richness, x=variable)) + geom_boxplot() + coord_flip()
library(ggtern)
install.packages("ggtern")
p_ind = ggplot(data.frame(Richness = test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li)) +
coord_tern() + geom_point(aes(
x=(Axis1-min(Axis1))/(max(Axis1)-min(Axis1)),
y=(Axis2-min(Axis2))/(max(Axis2)-min(Axis2)),
z=(Axis3-min(Axis3))/(max(Axis3)-min(Axis3)),
size=log2(Richness))) + theme_rgbw() +
xlab("PC1") + ylab("PC2") + zlab("PC3")
p_diet = ggplot(data.frame(dudi.pca(test[,-1], scannf=F, nf=10)$co)) +
coord_tern() + geom_text(aes(
x=(Comp1-min(Comp1))/(max(Comp1)-min(Comp1)),
y=(Comp2-min(Comp2))/(max(Comp2)-min(Comp2)),
z=(Comp3-min(Comp3))/(max(Comp3)-min(Comp3)),
label=row.names(dudi.pca(test[,-1], scannf=F, nf=10)$co)), size=3) + theme_rgbw() +
xlab("PC1 load.") + ylab("PC2 load.") + zlab("PC3 load.")
p_diet
p_ind_diet = grid.arrange(p_ind, p_diet, ncol=2)
library(gridExtra)
p_ind = ggplot(data.frame(Richness = test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li)) +
coord_tern() + geom_point(aes(
x=(Axis1-min(Axis1))/(max(Axis1)-min(Axis1)),
y=(Axis2-min(Axis2))/(max(Axis2)-min(Axis2)),
z=(Axis3-min(Axis3))/(max(Axis3)-min(Axis3)),
size=log2(Richness))) + theme_rgbw() +
xlab("PC1") + ylab("PC2") + zlab("PC3")
p_diet = ggplot(data.frame(dudi.pca(test[,-1], scannf=F, nf=10)$co)) +
coord_tern() + geom_text(aes(
x=(Comp1-min(Comp1))/(max(Comp1)-min(Comp1)),
y=(Comp2-min(Comp2))/(max(Comp2)-min(Comp2)),
z=(Comp3-min(Comp3))/(max(Comp3)-min(Comp3)),
label=row.names(dudi.pca(test[,-1], scannf=F, nf=10)$co)), size=3) + theme_rgbw() +
xlab("PC1 load.") + ylab("PC2 load.") + zlab("PC3 load.")
p_diet
p_ind_diet = grid.arrange(p_ind, p_diet, ncol=2)
library(ggtern)
p_ind_diet = grid.arrange(p_ind, p_diet, ncol=2)
library(ade4)
library(ggplot2)
library(reshape2)
library(ggtern)
library(gridExtra)
basal_diet_tmp = basal_diet[which(basal_diet$visite=="before_study"),]
for(i in 4:dim(basal_diet_tmp)[2]){
basal_diet_tmp[,i] = as.factor(as.character(basal_diet_tmp[,i]))
}
diet.acm = dudi.acm(basal_diet_tmp[,-c(1:3)], nf=3, scannf=FALSE)
diet.bca = bca(diet.acm, fac=as.factor(basal_diet_tmp[,1]), nf=3, scannf=F)
randtest(diet.bca)
basal_richness = metadata76[which(metadata76$Time.point=="2"),c("Subject_Id","richness")]
diet.richness.bca = merge(data.frame(Subject_Id=row.names(diet.bca$li), diet.bca$li), basal_richness, by="Subject_Id")
ggplot(diet.richness.bca, aes(x=Axis1, y=Axis2, col=richness, size=richness)) + geom_point()
#test = merge(basal_richness,basal_diet[which(basal_diet$visite=="before_study" & basal_diet$jour=="-1"),], by.x="Subject_Id", by.y="sujet")
#addition des jours precedent
test_combine = apply(basal_diet[which(basal_diet$visite=="after_study"),-c(1:3)] , 2, function(x){ tapply(as.numeric(x), as.factor(basal_diet[which(basal_diet$visite=="after_study"),1]) ,sum) } )
test_combine = apply(basal_diet[,-c(1:3)] , 2, function(x){ tapply(as.numeric(x), as.factor(basal_diet[,1]) ,sum) } )
test = merge(basal_richness, data.frame(Subject_Id=row.names(test_combine),test_combine) )
#test = merge(basal_richness,basal_diet[which(basal_diet$visite=="before_study"),], by.x="Subject_Id", by.y="sujet")
#test = merge(basal_richness,basal_diet[which(basal_diet$visite=="after_study"),], by.x="Subject_Id", by.y="sujet")
#test = merge(basal_richness,basal_diet, by.x="Subject_Id", by.y="sujet")
for(i in 5:dim(test)[2]){
test[,i] = as.numeric(as.character(test[,i]))
}
test = test[, c(2,4+which(apply(test[-c(1:4)],2,sum)!=0))]
cor(test$richness, test[,-1], method="spearman")
apply(test[,-1], 2, function(x) cor.test(test$richness, x, method="spearman")$p.value)
diet_test_result = data.frame(rho=t(cor(test$richness, test[,-1], method="spearman")), p=apply(test[,-1], 2, function(x) wilcox.test(test$richness~x)$p.value))
diet_test_result = data.frame(rho=t(cor(test$richness, test[,-1], method="spearman")), p=apply(test[,-1], 2, function(x) cor.test(test$richness,x, method="spearman")$p.value))
diet_test_result[which(diet_test_result$p<0.05),]
ggplot(melt(test, id.vars=c("richness")), aes(fill=as.character(value), y=richness, x=variable)) + geom_boxplot() + coord_flip()
p_ind = ggplot(data.frame(Richness = test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li)) +
coord_tern() + geom_point(aes(
x=(Axis1-min(Axis1))/(max(Axis1)-min(Axis1)),
y=(Axis2-min(Axis2))/(max(Axis2)-min(Axis2)),
z=(Axis3-min(Axis3))/(max(Axis3)-min(Axis3)),
size=log2(Richness))) + theme_rgbw() +
xlab("PC1") + ylab("PC2") + zlab("PC3")
p_diet = ggplot(data.frame(dudi.pca(test[,-1], scannf=F, nf=10)$co)) +
coord_tern() + geom_text(aes(
x=(Comp1-min(Comp1))/(max(Comp1)-min(Comp1)),
y=(Comp2-min(Comp2))/(max(Comp2)-min(Comp2)),
z=(Comp3-min(Comp3))/(max(Comp3)-min(Comp3)),
label=row.names(dudi.pca(test[,-1], scannf=F, nf=10)$co)), size=3) + theme_rgbw() +
xlab("PC1 load.") + ylab("PC2 load.") + zlab("PC3 load.")
p_diet
p_ind_diet = grid.arrange(p_ind, p_diet, ncol=2)
p_ind_diet
p_diet
p_ind
p_ind = ggplot(data.frame(Richness = test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li)) +
coord_tern() + geom_point(aes(
x=(Axis1-min(Axis1))/(max(Axis1)-min(Axis1)),
y=(Axis2-min(Axis2))/(max(Axis2)-min(Axis2)),
z=(Axis3-min(Axis3))/(max(Axis3)-min(Axis3)),
size=Richness)) + theme_rgbw() +
xlab("PC1") + ylab("PC2") + zlab("PC3")
p_diet = ggplot(data.frame(dudi.pca(test[,-1], scannf=F, nf=10)$co)) +
coord_tern() + geom_text(aes(
x=(Comp1-min(Comp1))/(max(Comp1)-min(Comp1)),
y=(Comp2-min(Comp2))/(max(Comp2)-min(Comp2)),
z=(Comp3-min(Comp3))/(max(Comp3)-min(Comp3)),
label=row.names(dudi.pca(test[,-1], scannf=F, nf=10)$co)), size=3) + theme_rgbw() +
xlab("PC1 load.") + ylab("PC2 load.") + zlab("PC3 load.")
grid.arrange(p_ind, p_diet, ncol=2)
warnings()
p_ind
data.frame(Richness = test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li)
cor(data.frame(Richness = test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li), method="spearman")
plot(test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li[,3], cex=test$richness)
plot(test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li[,3])
plot(log2(test$richness), dudi.pca(test[,-1], scannf=F, nf=10)$li[,3])
plot(log(test$richness), dudi.pca(test[,-1], scannf=F, nf=10)$li[,3])
plot(test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li[,3])
plot(10^test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li[,3])
plot(1/log10(test$richness), dudi.pca(test[,-1], scannf=F, nf=10)$li[,3])
plot(log10(test$richness), log10(dudi.pca(test[,-1], scannf=F, nf=10)$li[,3]))
plot(test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li[,3])
plot(1/test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li[,3])
plot(sqrt(test$richness), dudi.pca(test[,-1], scannf=F, nf=10)$li[,3])
plot(1/test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li[,3])
plot(log10(1/test$richness), dudi.pca(test[,-1], scannf=F, nf=10)$li[,3])
plot(1/test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li[,3])
plot(log10(1/test$richness), dudi.pca(test[,-1], scannf=F, nf=10)$li[,3])
hist(test$richness)
plot(test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li[,3])
p_ind
p_diet
basal_diet[,1]
#addition des jours precedent
test_combine = apply(basal_diet[which(basal_diet$visite=="before_study_study"),-c(1:3)] , 2, function(x){ tapply(as.numeric(x), as.factor(basal_diet[which(basal_diet$visite=="after_study"),1]) ,sum) } )
#test_combine = apply(basal_diet[,-c(1:3)] , 2, function(x){ tapply(as.numeric(x), as.factor(basal_diet[,1]) ,sum) } )
test = merge(basal_richness, data.frame(Subject_Id=row.names(test_combine),test_combine) )
#addition des jours precedent
test_combine = apply(basal_diet[which(basal_diet$visite=="before_study"),-c(1:3)] , 2, function(x){ tapply(as.numeric(x), as.factor(basal_diet[which(basal_diet$visite=="after_study"),1]) ,sum) } )
#test_combine = apply(basal_diet[,-c(1:3)] , 2, function(x){ tapply(as.numeric(x), as.factor(basal_diet[,1]) ,sum) } )
test = merge(basal_richness, data.frame(Subject_Id=row.names(test_combine),test_combine) )
#test = merge(basal_richness,basal_diet[which(basal_diet$visite=="before_study"),], by.x="Subject_Id", by.y="sujet")
#test = merge(basal_richness,basal_diet[which(basal_diet$visite=="after_study"),], by.x="Subject_Id", by.y="sujet")
#test = merge(basal_richness,basal_diet, by.x="Subject_Id", by.y="sujet")
for(i in 5:dim(test)[2]){
test[,i] = as.numeric(as.character(test[,i]))
}
test = test[, c(2,4+which(apply(test[-c(1:4)],2,sum)!=0))]
cor(test$richness, test[,-1], method="spearman")
apply(test[,-1], 2, function(x) cor.test(test$richness, x, method="spearman")$p.value)
diet_test_result = data.frame(rho=t(cor(test$richness, test[,-1], method="spearman")), p=apply(test[,-1], 2, function(x) wilcox.test(test$richness~x)$p.value))
diet_test_result = data.frame(rho=t(cor(test$richness, test[,-1], method="spearman")), p=apply(test[,-1], 2, function(x) cor.test(test$richness,x, method="spearman")$p.value))
diet_test_result[which(diet_test_result$p<0.05),]
ggplot(melt(test, id.vars=c("richness")), aes(fill=as.character(value), y=richness, x=variable)) + geom_boxplot() + coord_flip()
p_ind = ggplot(data.frame(Richness = test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li)) +
coord_tern() + geom_point(aes(
x=(Axis1-min(Axis1))/(max(Axis1)-min(Axis1)),
y=(Axis2-min(Axis2))/(max(Axis2)-min(Axis2)),
z=(Axis3-min(Axis3))/(max(Axis3)-min(Axis3)),
size=Richness)) + theme_rgbw() +
xlab("PC1") + ylab("PC2") + zlab("PC3")
p_diet = ggplot(data.frame(dudi.pca(test[,-1], scannf=F, nf=10)$co)) +
coord_tern() + geom_text(aes(
x=(Comp1-min(Comp1))/(max(Comp1)-min(Comp1)),
y=(Comp2-min(Comp2))/(max(Comp2)-min(Comp2)),
z=(Comp3-min(Comp3))/(max(Comp3)-min(Comp3)),
label=row.names(dudi.pca(test[,-1], scannf=F, nf=10)$co)), size=3) + theme_rgbw() +
xlab("PC1 load.") + ylab("PC2 load.") + zlab("PC3 load.")
p_diet
cor(test$richness, test[,-1], method="spearman")
cor(data.frame(Richness = test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li), method="spearman")
plot(test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li[,3])
plot(test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li[,2])
plot(test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li[,1])
plot(dudi.pca(test[,-1], scannf=F, nf=10)$li[,2], dudi.pca(test[,-1], scannf=F, nf=10)$li[,3])
plot(dudi.pca(test[,-1], scannf=F, nf=10)$li[,2], dudi.pca(test[,-1], scannf=F, nf=10)$li[,3], cex=log10(test$richness)/2)
plot(dudi.pca(test[,-1], scannf=F, nf=10)$li[,2], dudi.pca(test[,-1], scannf=F, nf=10)$li[,3], cex=log10(test$richness))
plot(dudi.pca(test[,-1], scannf=F, nf=10)$li[,2], dudi.pca(test[,-1], scannf=F, nf=10)$li[,3], cex=test$richness)
plot(dudi.pca(test[,-1], scannf=F, nf=10)$li[,2], dudi.pca(test[,-1], scannf=F, nf=10)$li[,3], cex=test$richness/10)
plot(dudi.pca(test[,-1], scannf=F, nf=10)$li[,2], dudi.pca(test[,-1], scannf=F, nf=10)$li[,3], cex=test$richness/100)
plot(dudi.pca(test[,-1], scannf=F, nf=10)$li[,2], dudi.pca(test[,-1], scannf=F, nf=10)$li[,3], cex=sqrt(test$richness))
plot(dudi.pca(test[,-1], scannf=F, nf=10)$li[,2], dudi.pca(test[,-1], scannf=F, nf=10)$li[,3], cex=1/test$richness)
plot(dudi.pca(test[,-1], scannf=F, nf=10)$li[,2], dudi.pca(test[,-1], scannf=F, nf=10)$li[,3], cex=test$richness/10)
plot(dudi.pca(test[,-1], scannf=F, nf=10)$li[,2], dudi.pca(test[,-1], scannf=F, nf=10)$li[,3], cex=test$richness/100)
plot(test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li[,3])
plot(test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li[,3], cex=test$richness/100)
cor.test(test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li[,3], method="spearman")
cor.test(test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li[,3]+dudi.pca(test[,-1], scannf=F, nf=10)$li[,2], method="spearman")
cor.test(test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li[,3]-dudi.pca(test[,-1], scannf=F, nf=10)$li[,2], method="spearman")
plot(test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li[,3]-dudi.pca(test[,-1], scannf=F, nf=10)$li[,2], method="spearman")
plot(test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li[,3]-dudi.pca(test[,-1], scannf=F, nf=10)$li[,2])
p_ind = ggplot(data.frame(Richness = test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li)) +
geom_point(aes(
y=(Axis3-min(Axis3))/(max(Axis3)-min(Axis3)),
x=(Axis2-min(Axis2))/(max(Axis2)-min(Axis2)),
size=Richness)) +
xlab("PC2") + ylab("PC3")
p_ind
p_ind = ggplot(data.frame(Richness = test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li)) +
geom_point(aes(
y=Axis3-Axis2,
x=Richness,
size=Richness)) +
xlab("PC2") + ylab("PC3")
p_ind
p_ind = ggplot(data.frame(Richness = test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li)) +
geom_point(aes(
y=Axis3-Axis2,
x=Richness,
size=Richness))
p_ind
p_ind = ggplot(data.frame(Richness = test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li)) +
geom_point(aes(
y=sqrt(Axis3^2+Axis2^2),
x=Richness,
size=Richness))
p_ind
tmp = data.frame(Richness = test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li)
cor.test(data=tmp, x=Richness, y=sqrt(Axis3^2+Axis2^2))
cor.test(x=tmp$Richness, y=sqrt(tmp$Axis3^2+tmp$Axis2^2))
cor.test(x=tmp$Richness, y=sqrt(tmp$Axis3^2+tmp$Axis2^2), method="spearman")
cor.test(x=tmp$Richness, y=(tmp$Axis3-tmp$Axis2), method="spearman")
cor.test(x=tmp$Richness, y=(tmp$Axis3-tmp$Axis2-tmp$Axis1), method="spearman")
lm(tmp$Richness~ tmp$Axis3-tmp$Axis2-tmp$Axis1)
anova(lm(tmp$Richness~ tmp$Axis3-tmp$Axis2-tmp$Axis1))
summary(anova(lm(tmp$Richness~ tmp$Axis3-tmp$Axis2-tmp$Axis1)))
summary(lm(tmp$Richness~ tmp$Axis3 + tmp$Axis2 + tmp$Axis1))
plot(lm(tmp$Richness~ tmp$Axis3 + tmp$Axis2 + tmp$Axis1))
plot(lm(tmp$Richness~ tmp$Axis3 - tmp$Axis2 - tmp$Axis1))
summary(lm(tmp$Richness~ tmp$Axis3 - tmp$Axis2 - tmp$Axis1))
summary(lm(tmp$Richness~ tmp$Axis3 + (-tmp$Axis2) + (-tmp$Axis1)))
cor(data.frame(Richness = test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li), method="pearson")
round(cor(data.frame(Richness = test$richness, dudi.pca(test[,-1], scannf=F, nf=10)$li), method="pearson"),3)
basal_diet_tmp = basal_diet[which(basal_diet$visite=="before_study"),]
for(i in 4:dim(basal_diet_tmp)[2]){
basal_diet_tmp[,i] = as.factor(as.character(basal_diet_tmp[,i]))
}
diet.acm = dudi.acm(basal_diet_tmp[,-c(1:3)], nf=3, scannf=FALSE)
diet.bca = bca(diet.acm, fac=as.factor(basal_diet_tmp[,1]), nf=3, scannf=F)
randtest(diet.bca)
basal_diet_tmp
basal_richness = metadata76[which(metadata76$Time.point=="2"),c("Subject_Id","richness")]
basal_richness
diet.richness.bca = merge(data.frame(Subject_Id=row.names(diet.bca$li), diet.bca$li), basal_richness, by="Subject_Id")
diet.richness.bca
plot(diet.bca)
ggplot(diet.richness.bca, aes(x=Axis1, y=Axis2, col=richness, size=richness)) + geom_point()
ggplot(diet.richness.bca, aes(x=Axis1, y=Axis3, col=richness, size=richness)) + geom_point()
cor(diet.richness.bca[2:5])
round(cor(diet.richness.bca[2:5]),2)
round(cor(diet.richness.bca[2:5]),3)
ggplot(diet.richness.bca, aes(x=Axis2+Axis3, y=richness, col=richness, size=richness)) + geom_point()
ggplot(diet.richness.bca, aes(x=Axis2+Axis3, y=richness, col=richness, size=richness)) + geom_point() + geom_smooth()
ggplot(diet.richness.bca, aes(x=Axis2+Axis3, y=richness, col=richness, size=richness)) + geom_point()
diet.bca$co
diet.bca$co[,2]+diet.bca$co[,3]
diet.bca$co[2]+diet.bca$co[3]
diet.bca$co[2]+diet.bca$co[3][order(diet.bca$co[,3]),]
diet.bca$co[3][order(diet.bca$co[,3]),]
diet.bca$co[3][order(diet.bca$co[,3])]
diet.bca$co[3][order(diet.bca$co[,3]),1]
diet.bca$co[order(diet.bca$co[,3]),]
diet.bca$co[order(diet.bca$co[,3]+diet.bca$co[,2]),]
diet.bca
basal_diet_tmp
diet.bca$co[order(diet.bca$co[,3]+diet.bca$co[,2]),]
diet.bca$co[order(diet.bca$co[,3]+diet.bca$co[,2]),]ggplot(diet.richness.bca, aes(x=Axis2+Axis3, y=richness, col=richness, size=richness)) + geom_point()
ggplot(diet.richness.bca, aes(x=Axis2+Axis3, y=richness, col=richness, size=richness)) + geom_point()
plot(diet.bca)
plot(diet.bca, yax=3)
plot(diet.bca, xax=2, yax=3)
dev.new()
ggplot(diet.richness.bca, aes(x=Axis2+Axis3, y=richness, col=richness, size=richness)) + geom_point()
ggplot(diet.richness.bca, aes(x=Axis2, y=Axis3, col=richness, size=richness)) + geom_point()
diet.bca
diet.bca$ls
as.factor(basal_diet_tmp[,1])
diet.richness.bca
merge(data.frame(diet.bca$ls, Subject_Id="basal_diet_tmp[,1]"),diet.richness.bca, by= "Subject_Id")
data.frame(diet.bca$ls, Subject_Id="basal_diet_tmp[,1]")
merge(data.frame(diet.bca$ls, Subject_Id=basal_diet_tmp[,1]),diet.richness.bca, by= "Subject_Id")
diet.richness.bca2 = merge(data.frame(diet.bca$ls, Subject_Id=basal_diet_tmp[,1]),diet.richness.bca, by= "Subject_Id")
ggplot(diet.richness.bca2) + geom_segment(aes(x=CS2,y=CS3, xend=Axis2, yend=Axis3)) + geom_text(aes(x=Axis2, y=Axis3, label=Subject_Id))
ggplot(diet.richness.bca2) + geom_point(aes(x=CS2,y=CS3)) + geom_segment(aes(x=CS2,y=CS3, xend=Axis2, yend=Axis3)) + geom_text(aes(x=Axis2, y=Axis3, label=Subject_Id))
ggplot(diet.richness.bca2) + geom_point(aes(x=CS2,y=CS3)) + geom_segment(aes(x=CS2,y=CS3, xend=Axis2, yend=Axis3, col=Subject_Id)) + geom_text(aes(x=Axis2, y=Axis3, label=Subject_Id))
p1_diet = ggplot(diet.richness.bca2) +
geom_point(aes(x=CS2,y=CS3)) +
geom_segment(aes(x=CS2,y=CS3, xend=Axis2, yend=Axis3, col=Subject_Id)) +
geom_text(aes(x=Axis2, y=Axis3, label=Subject_Id))
p2_diet = ggplot(diet.richness.bca, aes(x=Axis2+Axis3, y=richness, col=richness, size=richness)) + geom_point()
grid.arrange(p1_diet,p2_diet, ncol=2)
p1_diet = ggplot(diet.richness.bca2) +
geom_point(aes(x=CS2+CS3,y=richness)) +
geom_segment(aes(x=CS2+CS3,y=richness, xend=Axis2+Axis3, yend=richness, col=Subject_Id)) +
geom_text(aes(x=Axis2+Axis3, y=richness, label=Subject_Id))
p2_diet = ggplot(diet.richness.bca, aes(x=Axis2+Axis3, y=richness, col=richness, size=richness)) + geom_point()
grid.arrange(p1_diet,p2_diet, ncol=2)
p1_diet = ggplot(diet.richness.bca2) +
geom_point(aes(x=CS2,y=CS3)) +
geom_segment(aes(x=CS2,y=CS3, xend=Axis2, yend=Axis3, col=Subject_Id)) +
geom_text(aes(x=Axis2, y=Axis3, label=Subject_Id)) + xlab("Axis2") +ylab("Axis3")
# p1_diet = ggplot(diet.richness.bca2) +
# geom_point(aes(x=CS2+CS3,y=richness)) +
# geom_segment(aes(x=CS2+CS3,y=richness, xend=Axis2+Axis3, yend=richness, col=Subject_Id)) +
# geom_text(aes(x=Axis2+Axis3, y=richness, label=Subject_Id))
p2_diet = ggplot(diet.richness.bca, aes(x=Axis2+Axis3, y=richness, col=richness, size=richness)) + geom_point()
grid.arrange(p1_diet,p2_diet, ncol=2)
basal_diet = read.csv2("data-raw/diet_basal_alimintest.csv")
basal_diet$sujet = as.character(basal_diet$sujet)
basal_diet$jour = as.character(basal_diet$jour)
for(i in 4:dim(basal_diet)[2]){
basal_diet[,i] = as.factor(basal_diet[,i])
}
save(basal_diet,file="data/basal_diet.RData")
#analysis
library(ade4)
library(ggplot2)
library(reshape2)
library(ggtern)
library(gridExtra)
basal_diet_tmp = basal_diet[which(basal_diet$visite=="before_study"),]
for(i in 4:dim(basal_diet_tmp)[2]){
basal_diet_tmp[,i] = as.factor(as.character(basal_diet_tmp[,i]))
}
diet.acm = dudi.acm(basal_diet_tmp[,-c(1:3)], nf=3, scannf=FALSE)
diet.bca = bca(diet.acm, fac=as.factor(basal_diet_tmp[,1]), nf=3, scannf=F)
randtest(diet.bca)
#conclusion: individual can be clustered according their diet.
basal_richness = metadata76[which(metadata76$Time.point=="2"),c("Subject_Id","richness")]
diet.richness.bca = merge(data.frame(Subject_Id=row.names(diet.bca$li), diet.bca$li), basal_richness, by="Subject_Id")
#ggplot(diet.richness.bca, aes(x=Axis1, y=Axis2, col=richness, size=richness)) + geom_point()
diet.bca$co[order(diet.bca$co[,3]+diet.bca$co[,2]),]
diet.richness.bca2 = merge(data.frame(diet.bca$ls, Subject_Id=basal_diet_tmp[,1]),diet.richness.bca, by= "Subject_Id")
p1_diet = ggplot(diet.richness.bca2) +
geom_point(aes(x=CS2,y=CS3)) +
geom_segment(aes(x=CS2,y=CS3, xend=Axis2, yend=Axis3, col=Subject_Id)) +
geom_text(aes(x=Axis2, y=Axis3, label=Subject_Id)) + xlab("Axis2") +ylab("Axis3")
# p1_diet = ggplot(diet.richness.bca2) +
# geom_point(aes(x=CS2+CS3,y=richness)) +
# geom_segment(aes(x=CS2+CS3,y=richness, xend=Axis2+Axis3, yend=richness, col=Subject_Id)) +
# geom_text(aes(x=Axis2+Axis3, y=richness, label=Subject_Id))
p2_diet = ggplot(diet.richness.bca, aes(x=Axis2+Axis3, y=richness, col=richness, size=richness)) + geom_point()
grid.arrange(p1_diet,p2_diet, ncol=2)
p1_diet = ggplot(diet.richness.bca2) +
geom_point(aes(x=CS2,y=CS3)) +
geom_segment(aes(x=CS2,y=CS3, xend=Axis2, yend=Axis3, col=Subject_Id)) +
geom_text(aes(x=Axis2, y=Axis3, label=Subject_Id)) + xlab("Axis2") +ylab("Axis3") + guides(col=FALSE)
# p1_diet = ggplot(diet.richness.bca2) +
# geom_point(aes(x=CS2+CS3,y=richness)) +
# geom_segment(aes(x=CS2+CS3,y=richness, xend=Axis2+Axis3, yend=richness, col=Subject_Id)) +
# geom_text(aes(x=Axis2+Axis3, y=richness, label=Subject_Id))
p2_diet = ggplot(diet.richness.bca, aes(x=Axis2+Axis3, y=richness, col=richness)) + geom_point()
grid.arrange(p1_diet,p2_diet, ncol=2)
p1_diet = ggplot(diet.richness.bca2) +
geom_point(aes(x=CS2,y=CS3)) +
geom_segment(aes(x=CS2,y=CS3, xend=Axis2, yend=Axis3, col=Subject_Id)) +
geom_text(aes(x=Axis2, y=Axis3, label=Subject_Id)) + xlab("PC2") +ylab("PC3") + guides(col=FALSE)
# p1_diet = ggplot(diet.richness.bca2) +
# geom_point(aes(x=CS2+CS3,y=richness)) +
# geom_segment(aes(x=CS2+CS3,y=richness, xend=Axis2+Axis3, yend=richness, col=Subject_Id)) +
# geom_text(aes(x=Axis2+Axis3, y=richness, label=Subject_Id))
p2_diet = ggplot(diet.richness.bca, aes(x=Axis2+Axis3, y=richness, col=richness)) + geom_point() +
xlab("PC2 + PC3") +ylab("Microbiota Richness")
grid.arrange(p1_diet,p2_diet, ncol=2)
p1_diet = ggplot(diet.richness.bca2) +
geom_point(aes(x=CS2,y=CS3)) +
geom_segment(aes(x=CS2,y=CS3, xend=Axis2, yend=Axis3, col=Subject_Id)) +
geom_text(aes(x=Axis2, y=Axis3, label=Subject_Id)) + xlab("PC2") +ylab("PC3") + guides(col=FALSE)
# p1_diet = ggplot(diet.richness.bca2) +
# geom_point(aes(x=CS2+CS3,y=richness)) +
# geom_segment(aes(x=CS2+CS3,y=richness, xend=Axis2+Axis3, yend=richness, col=Subject_Id)) +
# geom_text(aes(x=Axis2+Axis3, y=richness, label=Subject_Id))
p2_diet = ggplot(diet.richness.bca, aes(x=Axis2+Axis3, y=richness, col=richness), size=5) + geom_point() +
xlab("PC2 + PC3") +ylab("Microbiota Richness")
grid.arrange(p1_diet,p2_diet, ncol=2)
p1_diet = ggplot(diet.richness.bca2) +
geom_point(aes(x=CS2,y=CS3)) +
geom_segment(aes(x=CS2,y=CS3, xend=Axis2, yend=Axis3, col=Subject_Id)) +
geom_text(aes(x=Axis2, y=Axis3, label=Subject_Id)) + xlab("PC2") +ylab("PC3") + guides(col=FALSE)
# p1_diet = ggplot(diet.richness.bca2) +
# geom_point(aes(x=CS2+CS3,y=richness)) +
# geom_segment(aes(x=CS2+CS3,y=richness, xend=Axis2+Axis3, yend=richness, col=Subject_Id)) +
# geom_text(aes(x=Axis2+Axis3, y=richness, label=Subject_Id))
p2_diet = ggplot(diet.richness.bca, aes(x=Axis2+Axis3, y=richness, col=richness)) + geom_point(size=5) +
xlab("PC2 + PC3") +ylab("Microbiota Richness")
grid.arrange(p1_diet,p2_diet, ncol=2)
dev.size()
dim(basal_diet)
head(basal_diet)
library(rmarkdown)
getwd()
load_all(Alimintest)
load_all()
library(devtools)
load_all()
data(basal_diet)
render("vignettes/basal_diet.Rmd")
render("vignettes/basal_diet.Rmd")
p1_diet = ggplot(diet.richness.bca2) +
geom_point(aes(x=CS2,y=CS3)) +
geom_segment(aes(x=CS2,y=CS3, xend=Axis2, yend=Axis3, col=Subject_Id)) +
geom_text(aes(x=Axis2, y=Axis3, label=Subject_Id)) + xlab("PC2") +ylab("PC3") + guides(col=FALSE)
p2_diet = ggplot(diet.richness.bca, aes(x=Axis2+Axis3, y=richness, col=richness)) +
geom_point(size=4) + xlab("PC2 + PC3") + ylab("Microbiota richness")
grid.arrange(p1_diet,p2_diet, ncol=2)
q()