forked from kunstner/2022_canine_atopic_dermatitis_paper
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path03_Modelling.R
422 lines (339 loc) · 19.5 KB
/
03_Modelling.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
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
# General Information -----------------------------------------------------
# Authors: Axel Künstner
# Date: 2022-05-31
# Libraries ---------------------------------------------------------------
library(phyloseq)
library(tidyverse)
library(patchwork)
library(ggpubr)
# devtools::install_github("davidgohel/flextable")
library(flextable)
theme_user <-
theme(axis.line = element_line(colour = "black"),
legend.position = "right",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.text.x = element_text(angle = 0, size=12),
strip.text.x = element_text(angle = 0, size = 12),
axis.text.y = element_text(size=12, face = "italic"),
axis.ticks.x=element_blank(),
strip.background = element_rect(colour="white", fill="white", size=1.5, linetype="solid"))
layout <- "
AABB
DDCC
"
layout3 <- "
ABC
"
layout2 <- "
AB
"
# Staphylococcus ----------------------------------------------------------
ps <- readRDS(file = "data/phyloseq.OTU.RDS")
sample_data(ps)$group <- gsub(pattern = "healthy",
replacement = "Healthy",
x = sample_data(ps)$group)
ps_clr <- microbiome::transform(ps, "clr")
ps_ra <- transform_sample_counts(ps, function(x){x / sum(x)})
ps.sta <- phyloseq::tax_glom(physeq = ps_ra %>% subset_samples(., group != "ADpost"),
taxrank = "Genus")
table(sample_data(ps.sta)$group)
table(sample_data(ps.sta)$location)
genus_table <- otu_table(ps.sta) %>% data.frame() %>% t()
colnames(genus_table) <- data.frame(tax_table(ps.sta))$Genus
cova <- sample_data(ps.sta) %>% data.frame()
sum( rownames(cova) == rownames(genus_table) ) == nrow(cova)
sort( colnames(genus_table) )
cova$Staphylococcus <- genus_table[, 'Staphylococcus']
cova$Fusobacterium <- genus_table[, 'Fusobacterium']
cova$Test <- genus_table[, 'Megamonas']
cova$group <- factor(cova$group, levels = c("ADpre", "ADtreatment", "Healthy" ))
cova$CADESI[is.na(cova$CADESI)] <- 0
cova$PVAS[is.na(cova$PVAS)] <- 0
m1 <- lm(formula = Staphylococcus ~ CADESI + group, data = cova %>% dplyr::filter(location == "abdomen"))
m2 <- lm(formula = Staphylococcus ~ CADESI + group, data = cova %>% dplyr::filter(location == "axilla"))
m3 <- lm(formula = Staphylococcus ~ CADESI + group, data = cova %>% dplyr::filter(location == "ear"))
m4 <- lm(formula = Staphylococcus ~ CADESI + group, data = cova %>% dplyr::filter(location == "elbow"))
m5 <- lm(formula = Staphylococcus ~ CADESI + group, data = cova %>% dplyr::filter(location == "flank"))
m6 <- lm(formula = Staphylococcus ~ CADESI + group, data = cova %>% dplyr::filter(location == "frontPaw"))
m7 <- lm(formula = Staphylococcus ~ CADESI + group, data = cova %>% dplyr::filter(location == "groin"))
m8 <- lm(formula = Staphylococcus ~ CADESI + group, data = cova %>% dplyr::filter(location == "hindPaw"))
m9 <- lm(formula = Staphylococcus ~ CADESI + group, data = cova %>% dplyr::filter(location == "lips"))
m10<- lm(formula = Staphylococcus ~ CADESI + group, data = cova %>% dplyr::filter(location == "palm"))
m11<- lm(formula = Staphylococcus ~ CADESI + group, data = cova %>% dplyr::filter(location == "perineum"))
m12<- lm(formula = Staphylococcus ~ CADESI + group, data = cova %>% dplyr::filter(location == "tail"))
m13<- lm(formula = Test ~ PVAS + group, data = cova %>% dplyr::filter(location == "stool"))
m14<- lm(formula = Fusobacterium ~ PVAS + group, data = cova %>% dplyr::filter(location == "stool"))
summary(m1) #
summary(m2) #
summary(m3)
summary(m4) #
summary(m5) #
summary(m6) #
summary(m7)
summary(m8) #
summary(m9) #
summary(m10) #
summary(m11)
summary(m12)
summary(m13)
summary(m1)
summary(m2)
summary(m4)
summary(m5)
summary(m6)
summary(m8)
summary(m9)
summary(m10)
sjPlot::plot_model(m10, type = "pred", terms = "CADESI") + theme_user + ggtitle('') + ylab('Staphylococcus')
# abdomen
m1.group <- update(m1, Staphylococcus ~ group)
p1 <- sjPlot::plot_model(m1, show.values = T, show.p = T, sort.est = FALSE, vline.color = "grey45", colors = "grey25") + theme_user + ggtitle('')
p2 <- ggplot(cova %>% dplyr::filter(location == "abdomen"), aes(CADESI, Staphylococcus)) +
geom_point() + geom_smooth(method='lm') + theme_user + ggtitle('') + ylab('Staphylococcus') + ylim(-0.1,1.4)
p3 <- sjPlot::plot_model(m1.group, type = "pred", terms = "group", show.values = T, ) + theme_user + xlab('') + ylab('Staphylococcus') + ggtitle('')
emmeans::emmeans(m1.group, pairwise ~ group)$contrasts
em_df <- emmeans::emmeans(m1.group, pairwise ~ group)$contrasts %>% data.frame()
p4 <- flextable::flextable( emmeans::emmeans(m1.group, pairwise ~ group)$contrasts %>% as.data.frame() )
p1 + p2 + p3 +
grid::rasterGrob(flextable::as_raster(x = p4)) +
plot_layout(design = layout) +
plot_annotation(title = 'Abdomen')
ggsave(filename = "plots/model_Staph_abdomen.pdf", height = 9, width = 9)
effectsize::eta_squared(car::Anova(m1, type = 2), partial = FALSE, ci = 0.95)
plot_abdomen <-
p3 + geom_point(size = 4) + theme(axis.text.y = element_text(size = 16), axis.title.y = element_text(face = "italic", size = 16)) +
annotate("text", x = 2.5, y = 0.6, label = paste0("p = ", round(em_df$p.value[2],4)), size = 8) + ylim(-0.1,0.8) +
p2 + geom_point(size = 4) + theme(axis.text.y = element_text(size = 16), axis.title.y = element_text(face = "italic", size = 16)) +
plot_layout(design = layout2)
plot_abdomen
# Axilla
m2.group <- update(m2, Staphylococcus ~ group)
p1 <- sjPlot::plot_model(m2, show.values = T, show.p = T, sort.est = FALSE, vline.color = "grey45", colors = "grey25") + theme_user + ggtitle('')
p2 <- ggplot(cova %>% dplyr::filter(location == "axilla"), aes(CADESI, Staphylococcus)) +
geom_point() + geom_smooth(method='lm') + theme_user + ggtitle('') + ylab('Staphylococcus') + ylim(-0.1,1.4)
p3 <- sjPlot::plot_model(m2.group, type = "pred", terms = "group", show.values = T, ) + theme_user + xlab('') + ylab('Staphylococcus') + ggtitle('')
emmeans::emmeans(m2.group, pairwise ~ group)$contrasts
em_df <- emmeans::emmeans(m2.group, pairwise ~ group)$contrasts %>% data.frame()
p4 <- flextable::flextable( emmeans::emmeans(m2.group, pairwise ~ group)$contrasts %>% as.data.frame() )
p1 + p2 + p3 +
grid::rasterGrob(flextable::as_raster(x = p4)) +
plot_layout(design = layout) +
plot_annotation(title = 'Axilla')
ggsave(filename = "plots/model_Staph_axilla.pdf", height = 9, width = 9)
effectsize::eta_squared(car::Anova(m2, type = 2), partial = FALSE, ci = 0.95)
plot_axilla <-
p3 + geom_point(size = 4) + theme(axis.text.y = element_text(size = 16), axis.title.y = element_text(face = "italic", size = 16)) +
annotate("text", x = 2.5, y = 0.6, label = paste0("p = ", round(em_df$p.value[2],4)), size = 8) + ylim(-0.1,0.8) +
p2 + geom_point(size = 4) + theme(axis.text.y = element_text(size = 16), axis.title.y = element_text(face = "italic", size = 16)) +
plot_layout(design = layout2)
# Elbow
m4.group <- update(m4, Staphylococcus ~ group)
p1 <- sjPlot::plot_model(m4, show.values = T, sort.est = FALSE, transform = NULL, vline.color = "grey45", colors = "grey25") + theme_user + ggtitle('')
p2 <- ggplot(cova %>% dplyr::filter(location == "elbow"), aes(CADESI, Staphylococcus)) +
geom_point() + geom_smooth(method='lm') + theme_user + ggtitle('') + ylab('Staphylococcus') + ylim(-0.1,1.4)
p3 <- sjPlot::plot_model(m4.group, type = "pred", terms = "group", show.values = T, ) + theme_user + xlab('') + ylab('Staphylococcus') + ggtitle('')
emmeans::emmeans(m4.group, pairwise ~ group)$contrasts
em_df <- emmeans::emmeans(m4.group, pairwise ~ group)$contrasts %>% data.frame()
p4 <- flextable::flextable( emmeans::emmeans(m4.group, pairwise ~ group)$contrasts %>% as.data.frame() )
p1 + p2 + p3 +
grid::rasterGrob(flextable::as_raster(x = p4)) +
plot_layout(design = layout) +
plot_annotation(title = 'Elbow')
ggsave(filename = "plots/model_Staph_elbow.pdf", height = 9, width = 9)
effectsize::eta_squared(car::Anova(m4, type = 2), partial = FALSE, ci = 0.95)
plot_elbow <-
p3 + geom_point(size = 4) + theme(axis.text.y = element_text(size = 16), axis.title.y = element_text(face = "italic", size = 16)) +
annotate("text", x = 2.5, y = 0.6, label = paste0("p = ", round(em_df$p.value[2],4)), size = 8) + ylim(-0.1,0.8) +
p2 + geom_point(size = 4) + theme(axis.text.y = element_text(size = 16), axis.title.y = element_text(face = "italic", size = 16)) +
plot_layout(design = layout2)
# Flank
m5.group <- update(m5, Staphylococcus ~ group)
p1 <- sjPlot::plot_model(m5, show.values = T, show.p = T, sort.est = FALSE, transform = NULL, vline.color = "grey45", colors = "grey25") + theme_user + ggtitle('')
p2 <- ggplot(cova %>% dplyr::filter(location == "flank"), aes(CADESI, Staphylococcus)) +
geom_point() + geom_smooth(method='lm') + theme_user + ggtitle('') + ylab('Staphylococcus') + ylim(-0.1,1.4)
p3 <- sjPlot::plot_model(m5.group, type = "pred", terms = "group", show.values = T, ) + theme_user + xlab('') + ylab('Staphylococcus') + ggtitle('')
emmeans::emmeans(m5.group, pairwise ~ group)$contrasts
em_df <- emmeans::emmeans(m5.group, pairwise ~ group)$contrasts %>% data.frame()
p4 <- flextable::flextable( emmeans::emmeans(m5.group, pairwise ~ group)$contrasts %>% as.data.frame() )
p1 + p2 + p3 +
grid::rasterGrob(flextable::as_raster(x = p4)) +
plot_layout(design = layout) +
plot_annotation(title = 'Flank')
ggsave(filename = "plots/model_Staph_flank.pdf", height = 9, width = 9)
effectsize::eta_squared(car::Anova(m5, type = 2), partial = FALSE, ci = 0.95)
plot_flank <-
p3 + geom_point(size = 4) + theme(axis.text.y = element_text(size = 16), axis.title.y = element_text(face = "italic", size = 16)) +
annotate("text", x = 2.5, y = 0.6, label = paste0("p = ", round(em_df$p.value[2],4)), size = 8) + ylim(-0.1,0.8) +
p2 + geom_point(size = 4) + theme(axis.text.y = element_text(size = 16), axis.title.y = element_text(face = "italic", size = 16)) +
plot_layout(design = layout2)
# Front Paw
m6.group <- update(m6, Staphylococcus ~ group)
p1 <- sjPlot::plot_model(m6, show.values = T, show.p = T, sort.est = FALSE, transform = NULL, vline.color = "grey45", colors = "grey25") + theme_user + ggtitle('')
p2 <- ggplot(cova %>% dplyr::filter(location == "frontPaw"), aes(CADESI, Staphylococcus)) +
geom_point() + geom_smooth(method='lm') + theme_user + ggtitle('') + ylab('Staphylococcus') + ylim(-0.1,1.4)
p3 <- sjPlot::plot_model(m6.group, type = "pred", terms = "group", show.values = T, ) + theme_user + xlab('') + ylab('Staphylococcus') + ggtitle('')
emmeans::emmeans(m6.group, pairwise ~ group)$contrasts
em_df <- emmeans::emmeans(m6.group, pairwise ~ group)$contrasts %>% data.frame()
p4 <- flextable::flextable( emmeans::emmeans(m6.group, pairwise ~ group)$contrasts %>% as.data.frame() )
p1 + p2 + p3 +
grid::rasterGrob(flextable::as_raster(x = p4)) +
plot_layout(design = layout) +
plot_annotation(title = 'Front Paw')
ggsave(filename = "plots/model_Staph_frontPaw.pdf", height = 9, width = 9)
effectsize::eta_squared(car::Anova(m6, type = 2), partial = FALSE, ci = 0.95)
plot_frontpaw <-
p3 + geom_point(size = 4) + theme(axis.text.y = element_text(size = 16), axis.title.y = element_text(face = "italic", size = 16)) +
annotate("text", x = 2.5, y = 0.6, label = paste0("p = ", round(em_df$p.value[2],4)), size = 8) + ylim(-0.1,0.8) +
p2 + geom_point(size = 4) + theme(axis.text.y = element_text(size = 16), axis.title.y = element_text(face = "italic", size = 16)) +
plot_layout(design = layout2)
# hindPaw
m8.group <- update(m8, Staphylococcus ~ group)
p1 <- sjPlot::plot_model(m8, show.values = T, sort.est = FALSE, transform = NULL, vline.color = "grey45", colors = "grey25") + theme_user + ggtitle('')
p2 <- ggplot(cova %>% dplyr::filter(location == "hindPaw"), aes(CADESI, Staphylococcus)) +
geom_point() + geom_smooth(method='lm') + theme_user + ggtitle('') + ylab('Staphylococcus') + ylim(-0.1,1.4)
p3 <- sjPlot::plot_model(m8.group, type = "pred", terms = "group", show.values = T, ) + theme_user + xlab('') + ylab('Staphylococcus') + ggtitle('')
emmeans::emmeans(m8.group, pairwise ~ group)$contrasts
em_df <- emmeans::emmeans(m8.group, pairwise ~ group)$contrasts %>% data.frame()
p4 <- flextable::flextable( emmeans::emmeans(m8.group, pairwise ~ group)$contrasts %>% as.data.frame() )
p1 + p2 + p3 +
grid::rasterGrob(flextable::as_raster(x = p4)) +
plot_layout(design = layout) +
plot_annotation(title = 'Hind Paw')
ggsave(filename = "plots/model_Staph_hindPaw.pdf", height = 9, width = 9)
effectsize::eta_squared(car::Anova(m8, type = 2), partial = FALSE, ci = 0.95)
plot_hindpaw <-
p3 + geom_point(size = 4) + theme(axis.text.y = element_text(size = 16), axis.title.y = element_text(face = "italic", size = 16)) +
annotate("text", x = 2.5, y = 0.6, label = paste0("p = ", round(em_df$p.value[2],4)), size = 8) + ylim(-0.1,0.8) +
p2 + geom_point(size = 4) + theme(axis.text.y = element_text(size = 16), axis.title.y = element_text(face = "italic", size = 16)) +
plot_layout(design = layout2)
# Lips
m9.group <- update(m9, Staphylococcus ~ group)
p1 <- sjPlot::plot_model(m9, show.values = T, sort.est = FALSE, transform = NULL, vline.color = "grey45", colors = "grey25") + theme_user + ggtitle('')
p2 <- ggplot(cova %>% dplyr::filter(location == "lips"), aes(CADESI, Staphylococcus)) +
geom_point() + geom_smooth(method='lm') + theme_user + ggtitle('') + ylab('Staphylococcus') + ylim(-0.1,1.4)
p3 <- sjPlot::plot_model(m9.group, type = "pred", terms = "group", show.values = T, ) + theme_user + xlab('') + ylab('Staphylococcus') + ggtitle('')
emmeans::emmeans(m9.group, pairwise ~ group)$contrasts
em_df <- emmeans::emmeans(m9.group, pairwise ~ group)$contrasts %>% data.frame()
p4 <- flextable::flextable( emmeans::emmeans(m9.group, pairwise ~ group)$contrasts %>% as.data.frame() )
p1 + p2 + p3 +
grid::rasterGrob(flextable::as_raster(x = p4)) +
plot_layout(design = layout) +
plot_annotation(title = 'Lips')
ggsave(filename = "plots/model_Staph_lips.pdf", height = 9, width = 9)
effectsize::eta_squared(car::Anova(m9, type = 2), partial = FALSE, ci = 0.95)
plot_lips <-
p3 + geom_point(size = 4) + theme(axis.text.y = element_text(size = 16), axis.title.y = element_text(face = "italic", size = 16)) +
annotate("text", x = 2.5, y = 0.6, label = paste0("p = ", round(em_df$p.value[2],4)), size = 8) + ylim(-0.1,0.8) +
p2 + geom_point(size = 4) + theme(axis.text.y = element_text(size = 16), axis.title.y = element_text(face = "italic", size = 16)) +
plot_layout(design = layout2)
# palm
m10.group <- update(m10, Staphylococcus ~ group)
p1 <- sjPlot::plot_model(m10, show.values = T, show.p = T, sort.est = FALSE, transform = NULL, vline.color = "grey45", colors = "grey25") + theme_user + ggtitle('')
p2 <- ggplot(cova %>% dplyr::filter(location == "palm"), aes(CADESI, Staphylococcus)) +
geom_point() + geom_smooth(method='lm') + theme_user + ggtitle('') + ylab('Staphylococcus') + ylim(-0.1,1.4)
p3 <- sjPlot::plot_model(m10.group, type = "pred", terms = "group", show.values = T, ) + theme_user + xlab('') + ylab('Staphylococcus') + ggtitle('')
emmeans::emmeans(m10.group, pairwise ~ group)$contrasts
em_df <- emmeans::emmeans(m10.group, pairwise ~ group)$contrasts %>% data.frame()
p4 <- flextable::flextable( emmeans::emmeans(m10.group, pairwise ~ group)$contrasts %>% as.data.frame() )
p1 + p2 + p3 +
grid::rasterGrob(flextable::as_raster(x = p4)) +
plot_layout(design = layout) +
plot_annotation(title = 'Palm')
ggsave(filename = "plots/model_Staph_palm.pdf", height = 9, width = 9)
effectsize::eta_squared(car::Anova(m10, type = 2), partial = FALSE, ci = 0.95)
plot_palm <-
p3 + geom_point(size = 4) + theme(axis.text.y = element_text(size = 16), axis.title.y = element_text(face = "italic", size = 16)) +
annotate("text", x = 2.5, y = 0.6, label = paste0("p = ", round(em_df$p.value[2],4)), size = 8) + ylim(-0.1,0.8) +
p2 + geom_point(size = 4) + theme(axis.text.y = element_text(size = 16), axis.title.y = element_text(face = "italic", size = 16)) +
plot_layout(design = layout2)
# Plot interaction effect
# sjPlot::plot_model(m5, type = "int")
# Fusobacterium
summary(m14)
m14.group <- update(m14, Fusobacterium ~ group)
summary(m14.group)
emmeans::emmeans(m14.group, pairwise ~ group)$contrasts
# Pretty plotting ---------------------------------------------------------
layout3 <- "
AA
BB
CC
DD
EE
FF
GG
HH
"
p_model <- plot_abdomen / # sig.
plot_axilla / # sig.
plot_elbow / # sig.
plot_palm / # sig.
plot_lips / # sig.
plot_hindpaw / # not sig.
plot_flank / # not sig.
plot_frontpaw + # not sig.
plot_layout(design = layout3) +
plot_annotation(tag_levels = 'a')
ggsave(filename = "plots/FigS6.pdf", height = 27, width = 12, plot = p_model)
# Staphylococcus vs Neisseria ---------------------------------------------
ps <- readRDS(file = "phyloseq.OTU.RDS")
ps_clr <- microbiome::transform(ps, "clr")
ps_ra <- transform_sample_counts(ps, function(x){x / sum(x)})
# ps.genus <- phyloseq::tax_glom(physeq = ps_ra, taxrank = "Genus")
ps.genus <- microbiome::aggregate_taxa(x = ps_ra, level = "Genus")
table(sample_data(ps.genus)$group)
table(sample_data(ps.genus)$location)
genus_table <- otu_table(ps.genus) %>% data.frame() %>% t()
colnames(genus_table) <- data.frame(tax_table(ps.genus))$unique
cova <- sample_data(ps.genus) %>% data.frame()
sum( rownames(cova) == rownames(genus_table) ) == nrow(cova)
sort(colnames(genus_table))
cova$Staphylococcus <- genus_table[, 'Staphylococcus']
cova$Fusobacterium <- genus_table[, 'Fusobacterium']
cova$Neisseria <- genus_table[, 'Neisseria']
cova$Neisseria_unk <- genus_table[, 'Bacteria_Proteobacteria_Betaproteobacteria_Neisseriales_Neisseriaceae_']
summary(cova$Staphylococcus)
summary(cova$Neisseria)
summary(cova$Neisseria_unk)
cova$log_Staph <- log1p(cova$Staphylococcus)
cova$log_Neiss <- log1p(cova$Neisseria)
cova$log_Neiss_un <- log1p(cova$Neisseria_unk)
cova$log_Neiss_total <- log1p(cova$Neisseria + cova$Neisseria_unk)
cova$Staph_Neis <- cova$log_Staph - cova$log_Neiss
cova$group <- factor(cova$group, levels = c("ADpre", "ADpost", "ADtreatment", "Healthy" ))
comp <- list(
c('ADpre', 'ADpost'),
c('ADpre', 'ADtreatment'),
c('ADpre', 'Healthy'),
c('ADpost', 'ADtreatment'),
c('ADpost', 'Healthy'),
c('ADtreatment', 'Healthy')
)
cova %>%
dplyr::filter(location != "stool") %>%
ggviolin(data = ., x = 'group', y = 'Staph_Neis',
trim = T, facet.by = 'location',
draw_quantiles = c(0.1, 0.5, 0.9), add = 'jitter') +
# facet_wrap(~location, ncol = 3) +
stat_compare_means() +
stat_compare_means(comparisons = comp, method = "wilcox.test") +
ylim(-1.,2) + ylab('log ratio Staphylococcus vs Neisseria') +
xlab('') +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey45") +
theme(axis.line = element_line(colour = "black"),
legend.position = "right",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1, size=12),
strip.text.x = element_text(angle = 0, size = 12),
axis.text.y = element_text(size=12, face = "italic"),
axis.ticks.x=element_blank(),
strip.background = element_rect(colour="white", fill="white", size=1.5, linetype="solid"))
ggsave(filename = "plots/Staph_Neisseria.pdf", height = 15, width = 15)
# Session Info ------------------------------------------------------------
sessioninfo::session_info()