-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathFigure_1C.Rmd
166 lines (144 loc) · 13.2 KB
/
Figure_1C.Rmd
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
---
title: "Figure_1C"
author: "QD"
date: "2024-02-13"
output: html_document
---
## Load libaries
```{r}
library(readxl)
library(tidyverse)
library(ggpubr)
library(kableExtra)
library(ggthemes)
library(ggembl)
```
## Calculations
```{r}
Mesnage_2023_metadata_before_after <- Mesnage_2023_metadata %>% filter(Timepoint == "Before" | Timepoint == "After") %>% count(Study_Patient) %>% filter(n > 1)
Mesnage_2023_metadata_before_after <- Mesnage_2023_metadata %>% filter(Timepoint == "Before" | Timepoint == "After") %>% filter(Study_Patient %in% Mesnage_2023_metadata_before_after$Study_Patient)
Mesnage_2023_metadata_before_after <- Mesnage_2023_metadata_before_after %>% mutate(BMI = Body_Weight / (Height/100)^2 )
clinical_markers_of_interest <- c("Glucose", "HbA1c", "HDL", "LDL", "Triglycerides", "GGT", "Uric_Acid", "Urea", "Body_Weight", "BMI")
Mesnage_2023_metadata_before_after_clinical_markers <- Mesnage_2023_metadata_before_after %>% select(1:12, all_of(clinical_markers_of_interest))
signif.num <- function(x) {
symnum(x, corr = FALSE, na = FALSE, legend = FALSE,
cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", " "))
}
wilcox_results <- list()
for (marker in clinical_markers_of_interest) {
# Perform comparison of means for the current marker
comparison <- compare_means(as.formula(paste(marker, "~Timepoint")), method="wilcox.test", paired = TRUE, data = Mesnage_2023_metadata_before_after_clinical_markers)
# Store the result in the list
wilcox_results[[marker]] <- comparison
}
wilcox_results_df <- bind_rows(wilcox_results)
wilcox_results_df <- wilcox_results_df %>% mutate(pval.adj = p.adjust (p, method = "fdr")) ## Calculate adjusted p-values, but use raw p-values as this is what Robin et al. always do for their more clinical variables, which is obviously different from an omics-type of experiment.
wilcox_results_df <- wilcox_results_df %>% rename("Clinical_Marker" = 1)
Mesnage_2023_metadata_before_after_clinical_markers <- Mesnage_2023_metadata_before_after_clinical_markers %>%
arrange(Study_Patient, Timepoint) ## you need to reorder the input dataframe to match things https://github.com/tidyverse/ggplot2/issues/3535
```
## Now make all the subpanels for Figure 1C
```{r}
clinical_markers_n <- Mesnage_2023_metadata_before_after_clinical_markers %>% select(sample.id, Timepoint, Glucose:BMI) %>% group_by(Timepoint) %>% summarise(across(Glucose:BMI, ~ sum(!is.na(.x)))) ## These numbers need to be displayed in the plots to have an n for each comparison
Glucose_Figure <- ggplot() +
theme_publication() +
geom_boxplot(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = Glucose, fill = Timepoint), outlier.shape = NA) +
geom_point(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = Glucose), size = 0.2, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
geom_line(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = Glucose, group = Study_Patient), linetype = "11", linewidth = 0.3, alpha = 0.4, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
stat_pvalue_manual(wilcox_results_df %>% filter(Clinical_Marker == "Glucose"), label = "{p.adj}", color = "black", size = 2) +
theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), legend.position = "none") +
scale_fill_manual(values=c('#366eb2','#bc3f60')) +
ylab("Glucose (mg/dL)") +
geom_text(data = clinical_markers_n %>% select(Timepoint, Glucose), aes(x = Timepoint, y = 52, label = paste0("n=", Glucose)), size = 2)
HbA1c_Figure <- ggplot() +
theme_publication() +
geom_boxplot(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = HbA1c, fill = Timepoint), outlier.shape = NA) +
geom_point(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = HbA1c), size = 0.2, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
geom_line(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = HbA1c, group = Study_Patient), linetype = "11", linewidth = 0.3, alpha = 0.4, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
stat_pvalue_manual(wilcox_results_df %>% filter(Clinical_Marker == "HbA1c"), label = "{p.adj}", color = "black", size = 2) +
theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), legend.position = "none") +
scale_fill_manual(values=c('#366eb2','#bc3f60')) +
ylab("HbA1c (%)") +
geom_text(data = clinical_markers_n %>% select(Timepoint, HbA1c), aes(x = Timepoint, y = 4.5, label = paste0("n=", HbA1c)), size = 2)
HDL_Figure <- ggplot() +
theme_publication() +
geom_boxplot(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = HDL, fill = Timepoint), outlier.shape = NA) +
geom_point(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = HDL), size = 0.2, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
geom_line(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = HDL, group = Study_Patient), linetype = "11", linewidth = 0.3, alpha = 0.4, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
stat_pvalue_manual(wilcox_results_df %>% filter(Clinical_Marker == "HDL"), label = "{p.adj}", color = "black", size = 2) +
theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), legend.position = "none") +
scale_fill_manual(values=c('#366eb2','#bc3f60')) +
ylab("HDL (mg/dL)") +
geom_text(data = clinical_markers_n %>% select(Timepoint, HDL), aes(x = Timepoint, y = 20, label = paste0("n=", HDL)), size = 2)
LDL_Figure <- ggplot() +
theme_publication() +
geom_boxplot(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = LDL, fill = Timepoint), outlier.shape = NA) +
geom_point(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = LDL), size = 0.2, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
geom_line(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = LDL, group = Study_Patient), linetype = "11", linewidth = 0.3, alpha = 0.4, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
stat_pvalue_manual(wilcox_results_df %>% filter(Clinical_Marker == "LDL"), label = "{p.adj}", color = "black", size = 2) +
theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), legend.position = "none") +
scale_fill_manual(values=c('#366eb2','#bc3f60')) +
ylab("LDL (mg/dL)") +
geom_text(data = clinical_markers_n %>% select(Timepoint, LDL), aes(x = Timepoint, y = 25, label = paste0("n=", LDL)), size = 2)
Triglycerides_Figure <- ggplot() +
theme_publication() +
geom_boxplot(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = Triglycerides, fill = Timepoint), outlier.shape = NA) +
geom_point(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = Triglycerides), size = 0.2, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
geom_line(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = Triglycerides, group = Study_Patient), linetype = "11", linewidth = 0.3, alpha = 0.4, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
stat_pvalue_manual(wilcox_results_df %>% filter(Clinical_Marker == "Triglycerides"), label = "{p.adj}", color = "black", size = 2) +
theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), legend.position = "none") +
scale_fill_manual(values=c('#366eb2','#bc3f60')) +
ylab("Triglycerides (mg/dl)") +
geom_text(data = clinical_markers_n %>% select(Timepoint, Triglycerides), aes(x = Timepoint, y = 39, label = paste0("n=", Triglycerides)), size = 2)
GGT_Figure <- ggplot() +
theme_publication() +
geom_boxplot(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = GGT, fill = Timepoint), outlier.shape = NA) +
geom_point(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = GGT), size = 0.2, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
geom_line(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = GGT, group = Study_Patient), linetype = "11", linewidth = 0.3, alpha = 0.4, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
stat_pvalue_manual(wilcox_results_df %>% filter(Clinical_Marker == "GGT"), label = "{p.adj}", color = "black", size = 2) +
theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), legend.position = "none") +
scale_fill_manual(values=c('#366eb2','#bc3f60')) +
ylab("Gamma-GT (U/L)") +
geom_text(data = clinical_markers_n %>% select(Timepoint, GGT), aes(x = Timepoint, y = 5, label = paste0("n=", GGT)), size = 2)
Uric_Acid_Figure <- ggplot() +
theme_publication() +
geom_boxplot(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = Uric_Acid, fill = Timepoint), outlier.shape = NA) +
geom_point(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = Uric_Acid), size = 0.2, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
geom_line(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = Uric_Acid, group = Study_Patient), linetype = "11", linewidth = 0.3, alpha = 0.4, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
stat_pvalue_manual(wilcox_results_df %>% filter(Clinical_Marker == "Uric_Acid"), label = "{p.adj}", color = "black", size = 2) +
theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), legend.position = "none") +
scale_fill_manual(values=c('#366eb2','#bc3f60')) +
ylab("Uric Acid (mg/dl)") +
geom_text(data = clinical_markers_n %>% select(Timepoint, Uric_Acid), aes(x = Timepoint, y = 2.9, label = paste0("n=", Uric_Acid)), size = 2)
Urea_Figure <- ggplot() +
theme_publication() +
geom_boxplot(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = Urea, fill = Timepoint), outlier.shape = NA) +
geom_point(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = Urea), size = 0.2, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
geom_line(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = Urea, group = Study_Patient), linetype = "11", linewidth = 0.3, alpha = 0.4, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
stat_pvalue_manual(wilcox_results_df %>% filter(Clinical_Marker == "Urea"), label = "{p.adj}", color = "black", size = 2) +
theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), legend.position = "none") +
scale_fill_manual(values=c('#366eb2','#bc3f60')) +
ylab("Urea (mg/dl)") +
geom_text(data = clinical_markers_n %>% select(Timepoint, Urea), aes(x = Timepoint, y = 5, label = paste0("n=", Urea)), size = 2)
Body_Weight_Figure <- ggplot() +
theme_publication() +
geom_boxplot(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = Body_Weight, fill = Timepoint), outlier.shape = NA) +
geom_point(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = Body_Weight), size = 0.2, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
geom_line(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = Body_Weight, group = Study_Patient), linetype = "11", linewidth = 0.3, alpha = 0.4, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
stat_pvalue_manual(wilcox_results_df %>% filter(Clinical_Marker == "Body_Weight"), label = "{p.adj}", color = "black", size = 2) +
theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), legend.position = "none") +
scale_fill_manual(values=c('#366eb2','#bc3f60')) +
ylab("Body Weight (kg)") +
geom_text(data = clinical_markers_n %>% select(Timepoint, Body_Weight), aes(x = Timepoint, y = 42, label = paste0("n=", Body_Weight)), size = 2)
BMI_Figure <- ggplot() +
theme_publication() +
geom_boxplot(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = BMI, fill = Timepoint), outlier.shape = NA) +
geom_point(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = BMI), size = 0.2, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
geom_line(data = Mesnage_2023_metadata_before_after_clinical_markers, aes(x = Timepoint, y = BMI, group = Study_Patient), linetype = "11", linewidth = 0.3, alpha = 0.4, position = position_jitter(height = 0, width = 0.2, seed = 100)) +
stat_pvalue_manual(wilcox_results_df %>% filter(Clinical_Marker == "Body_Weight"), label = "{p.adj}", color = "black", size = 2) +
theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), legend.position = "none") +
scale_fill_manual(values=c('#366eb2','#bc3f60')) +
ylab("BMI (kg/m2)") +
geom_text(data = clinical_markers_n %>% select(Timepoint, BMI), aes(x = Timepoint, y = 40, label = paste0("n=", BMI)), size = 2)
```