## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Calculate confidence intervals for relative risk by simulation
# https://api.rpubs.com/juanhklopper/RelativeRiskConfidenceIntervals
simulate_group <- function(n, p){
xs = runif(n)
k = sum(xs < p)
return(k / n)
}
simulate_trial <- function(n1, p1, n2, p2){
risk1 = simulate_group(n1, p1)
risk2 = simulate_group(n2, p2)
efficacy = 1 - risk2 / risk1
return(efficacy)
}
#Group 1-3 DST1
DST1_endpoints = DST1_tot_def %>% filter(t==12 & is.element(Group,c(1,2,3))) %>% group_by(Treatment) %>% summarise(N=sum(S+R+C),I=sum(R+C))
print(DST1_endpoints)
## # A tibble: 2 × 3
## Treatment N I
## <chr> <dbl> <dbl>
## 1 C 48 43
## 2 V 55 30
DST1_eff <- comprehenr::to_vec(for (i in 1:10000) simulate_trial(DST1_endpoints$N[1],(DST1_endpoints$I[1])/(DST1_endpoints$N[1]),
DST1_endpoints$N[2],(DST1_endpoints$I[2])/(DST1_endpoints$N[2])))
DST1_eff <- signif(100*quantile(DST1_eff,c(0.025,0.5,0.975)),2)
DST 1 endpoint efficacy is 39%; 95% CI: 22-54%
DST10_endpoints = DST10_tot_def %>% filter(t==12 & is.element(Group,c(1,2,3))) %>% group_by(Treatment) %>% summarise(N=sum(S+R+C),I=sum(R+C))
print(DST10_endpoints)
## # A tibble: 2 × 3
## Treatment N I
## <chr> <dbl> <dbl>
## 1 C 52 41
## 2 V 59 28
DST10_eff <- comprehenr::to_vec(for (i in 1:10000) simulate_trial(DST10_endpoints$N[1],(DST10_endpoints$I[1])/(DST10_endpoints$N[1]),
DST10_endpoints$N[2],(DST10_endpoints$I[2])/(DST10_endpoints$N[2])))
DST10_eff <- signif(100*quantile(DST10_eff,c(0.025,0.5,0.975)),2)
DST 10 endpoint efficacy is 40%; 95% CI: 19-57%
vl_endpoints <- vl_df %>% filter(is.element(Group,c(3,4)) & OrganPM=='Total') %>% group_by(Treatment) %>% summarise(N=length(VL),I=sum(VL=='Yes'))
print(vl_endpoints)
## # A tibble: 2 × 3
## Treatment N I
## <chr> <int> <int>
## 1 Control 34 29
## 2 Vaccinate 33 21
vl_eff <- comprehenr::to_vec(for (i in 1:10000) simulate_trial(vl_endpoints$N[1],(vl_endpoints$I[1])/(vl_endpoints$N[1]),
vl_endpoints$N[2],(vl_endpoints$I[2])/(vl_endpoints$N[2])))
vl_eff <- signif(100*quantile(vl_eff,c(0.025,0.5,0.975)),2)
VL endpoint efficacy is 25%; 95% CI: 1.1-47%
# This script generates two box and whisker plots for Fig S2 to visually represent post-mortem (PM) data focused on visible lesions (VL) scores, complementing the data presented in Fig2.
# v122423 - fixed color schema for plot 2 to match Fig2 panels 3 and 4
# Install and load required libraries if not already installed
required_libraries <- c("tidyverse", "ggplot2", "ggsignif", "ggthemes", "scales", "patchwork")
new_libraries <- required_libraries[!required_libraries %in% installed.packages()[,"Package"]]
if(length(new_libraries)) install.packages(new_libraries)
lapply(required_libraries, library, character.only = TRUE)
##
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
##
## discard
## The following object is masked from 'package:readr':
##
## col_factor
## [[1]]
## [1] "lubridate" "forcats" "stringr" "dplyr" "purrr" "readr"
## [7] "tidyr" "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [13] "grDevices" "utils" "datasets" "methods" "base"
##
## [[2]]
## [1] "lubridate" "forcats" "stringr" "dplyr" "purrr" "readr"
## [7] "tidyr" "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [13] "grDevices" "utils" "datasets" "methods" "base"
##
## [[3]]
## [1] "ggsignif" "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [7] "readr" "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [13] "graphics" "grDevices" "utils" "datasets" "methods" "base"
##
## [[4]]
## [1] "ggthemes" "ggsignif" "lubridate" "forcats" "stringr" "dplyr"
## [7] "purrr" "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [13] "stats" "graphics" "grDevices" "utils" "datasets" "methods"
## [19] "base"
##
## [[5]]
## [1] "scales" "ggthemes" "ggsignif" "lubridate" "forcats" "stringr"
## [7] "dplyr" "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [13] "tidyverse" "stats" "graphics" "grDevices" "utils" "datasets"
## [19] "methods" "base"
##
## [[6]]
## [1] "patchwork" "scales" "ggthemes" "ggsignif" "lubridate" "forcats"
## [7] "stringr" "dplyr" "purrr" "readr" "tidyr" "tibble"
## [13] "ggplot2" "tidyverse" "stats" "graphics" "grDevices" "utils"
## [19] "datasets" "methods" "base"
# Import Total PM score data of Sentinels and seeders
TotalVLScoresVivek <- read.csv('./TotalVLScoresVivekcsv.csv', stringsAsFactors = FALSE)
# Define color vectors for each phase for Plot 1
TotalVLScoresVivek$CombinedGroup <- with(TotalVLScoresVivek,
ifelse(GroupCh == "Phase-I",
paste(GroupCh, Treatment, sep="_"),
paste(GroupCh, Treatment, sep="_")))
plot1_custom_colors <- c("Phase-I_Control" = "#F090D4", "Phase-I_Vaccinated" = "#893DF6",
"Phase-II_Control" = "#8ED5FB", "Phase-II_Vaccinated" = "#78CC6F")
# Plot 1 - Phase I and Phase-II
p1 = ggplot(data=TotalVLScoresVivek, aes(x = GroupCh, y = Score, fill = CombinedGroup)) +
geom_boxplot(width = 0.9, alpha=0.5) +
stat_boxplot(geom = "errorbar", width=0.4, size=0.8, position = position_dodge(0.9), alpha=0.7) +
geom_dotplot(binaxis = "y", stackdir = "center", binwidth = 1.6, position = position_dodge(0.9), dotsize = 1.4, alpha=0.8) +
scale_y_continuous(expand = c(0, 0)) +
expand_limits(y = c(-3, 125)) +
scale_y_continuous(breaks = seq(0, 125, 20)) +
labs(title='A',
subtitle = "Visible Lesion Score by Study Phase",
x = "Study Phase",
y = "Visible Lesion Score",
fill = "Combined Group") +
scale_fill_manual(values = plot1_custom_colors) +
theme_bw() +
theme(legend.position=c(0.7,0.8),
legend.text = element_text(size = 8),
text = element_text(size = 10),
plot.subtitle = element_text(hjust = 0.5)) +
annotate("text",x=1,y=125,label="p = 0.01") +
annotate("text",x=2,y=75,label="p = 0.01") +
annotate("segment",x=0.75,y=120,xend=1.25,yend=120)+
annotate("segment",x=1.75,y=70,xend=2.25,yend=70)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
# Save Plot 1
# ggsave(file = "Boxlot_of_VL_Score_Phase-I_and_Phase-II.tiff", units = "in", width = 5.1, height = 3.5, dpi = 1200)
# Import the sentinel PM data and filter the total score only for Plot 2
TotalScoreSentinels <- read.csv('./PMScoreSentinelsscsv.csv', stringsAsFactors = FALSE) %>%
select(ID1, Exposure, Treatment, Group, OrganPM, Score) %>%
filter(OrganPM == "Total")
# Ensure that Exposure and Treatment are factors with the correct levels for Plot 2
TotalScoreSentinels$Exposure <- factor(TotalScoreSentinels$Exposure, levels = c("Control", "Vaccinate"))
TotalScoreSentinels$Treatment <- factor(TotalScoreSentinels$Treatment, levels = c("Control", "Vaccinate"))
# Create a new column that combines Phase-II Exposure and Treatment
TotalScoreSentinels$ExposureTreatment <- with(TotalScoreSentinels, paste(Exposure, Treatment, sep = "_"))
# Define custom colors based on the combined Exposure and Treatment for Plot 2
custom_colors_plot2 <- c("Control_Control" = "#643F95",
"Control_Vaccinate" = "#AECDE1",
"Vaccinate_Control" = "#EF8632",
"Vaccinate_Vaccinate" = "#549E3F")
# Plot 2 - Exposure to C and V with custom colors
p2 = ggplot(TotalScoreSentinels, aes(x = Exposure, y = Score, fill = ExposureTreatment)) +
geom_boxplot(width = 0.9, alpha=0.5) +
stat_boxplot(geom = "errorbar", width=0.4, size=0.8, position = position_dodge(0.9), alpha=0.7) +
geom_dotplot(binaxis = "y", stackdir = "center", binwidth = 1.6, position = position_dodge(0.9), dotsize = 1.5, alpha=0.8) +
scale_y_continuous(expand = c(0, 0)) +
expand_limits(y = c(-3, 125)) +
scale_y_continuous(breaks = seq(0, 125, 20)) +
labs(title='B',subtitle = "Visible Lesion Score by Exposure",
x = "Exposure",
y = "Visible Lesion Score",
fill="Exposure Treatment") +
scale_fill_manual(values = custom_colors_plot2) +
theme_bw() +
theme(legend.position=c(0.7,0.8),
legend.text = element_text(size = 8),
text = element_text(size = 10),
plot.subtitle = element_text(hjust = 0.5)) +
annotate("text",x=1,y=75,label="p = 0.01") +
annotate("text",x=2,y=75,label="p = 0.25") +
annotate("segment",x=0.75,y=70,xend=1.25,yend=70)+
annotate("segment",x=1.75,y=70,xend=2.25,yend=70)
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
# Save Plot 2
#ggsave(file = "Boxplot_of_VL_Score_Sentinels_Phase-II_Exposure_Treatment.tiff", units = "in", width = 5.1, height = 3.5, dpi = 1200)
print(p1+p2)
ggsave(p1+p2,file = "../Manuscript_figures/Fig_S2.pdf", units = "in", width = 7.25, height = 5.0)