Skip to content

Instantly share code, notes, and snippets.

@USMortality

USMortality/puhach_et_al_2022.r

Last active August 16, 2024 16:38
Show Gist options
  • Save USMortality/9e6b16fa25e9785dfd659c402c170bd2 to your computer and use it in GitHub Desktop.
Save USMortality/9e6b16fa25e9785dfd659c402c170bd2 to your computer and use it in GitHub Desktop.
Download ZIP
Reanalysis Puhach et al, 2022
Raw
puhach_et_al_2022.r
# Study: https://www.nature.com/articles/s41591-022-01816-0
library(readxl)
library(dplyr)
library(ggplot2)
sf <- 2
options(vsc.dev.args = list(width = 600 * sf, height = 335 * sf, res = 72 * sf))
download.file(
paste0(
"https://static-content.springer.com/esm/art%3A10.1038%2Fs41591-022",
"-01816-0/MediaObjects/41591_2022_1816_MOESM3_ESM.xlsx"
),
"/tmp/data.xlsx"
)
df <-
read_excel("/tmp/data.xlsx", sheet = "Original SARS-CoV-2") |>
select(DPOS, "Vaccination status", "RNA load/ml", "FFU/ml", "Variant") |>
setNames(c("dpos", "vaxx_status", "rna_load_ml", "ffu_ml", "variant")) |>
mutate(dpos = as.integer(dpos), vaxx_status = as.integer(vaxx_status)) |>
filter(variant == "Delta") |>
mutate(vaxx_status_label = case_when(
vaxx_status == 0 ~ "Unvaccinated",
vaxx_status == 1 ~ "Vaccinated",
TRUE ~ as.character(vaxx_status)
))
df2 <- df |>
group_by(vaxx_status_label, dpos) |>
summarize(
rna_load_ml_mean = mean(rna_load_ml),
rna_load_ml_mean_lower = rna_load_ml_mean - qt(0.975, df = n() - 1) * sd(rna_load_ml) / sqrt(n()),
rna_load_ml_mean_upper = rna_load_ml_mean + qt(0.975, df = n() - 1) * sd(rna_load_ml) / sqrt(n()),
rna_load_ml_median = median(rna_load_ml),
rna_load_ml_median_lower = rna_load_ml_median - qt(0.975, df = n() - 1) * sd(rna_load_ml) / sqrt(n()),
rna_load_ml_median_upper = rna_load_ml_median + qt(0.975, df = n() - 1) * sd(rna_load_ml) / sqrt(n()),
ffu_ml_mean = mean(ffu_ml),
ffu_ml_mean_lower = ffu_ml_mean - qt(0.975, df = n() - 1) * sd(ffu_ml) / sqrt(n()),
ffu_ml_mean_upper = ffu_ml_mean + qt(0.975, df = n() - 1) * sd(ffu_ml) / sqrt(n()),
ffu_ml_median = median(ffu_ml),
ffu_ml_median_lower = ffu_ml_median - qt(0.975, df = n() - 1) * sd(ffu_ml) / sqrt(n()),
ffu_ml_median_upper = ffu_ml_median + qt(0.975, df = n() - 1) * sd(ffu_ml) / sqrt(n())
) |>
arrange(dpos)
# RNA Viral loads
## LOESS
ggplot(
df, aes(x = dpos, color = vaxx_status_label, fill = vaxx_status_label)
) +
geom_smooth(aes(y = rna_load_ml), method = "loess", se = TRUE) +
geom_point(aes(y = rna_load_ml)) +
labs(
title = "RNA viral loads [LOESS]",
subtitle = paste0(
"Source: Puhach et al., 2022: ",
"https://www.nature.com/articles/s41591-022-01816-0"
),
x = "DPOS",
y = "genome copies per ml"
)
## Mean
ggplot(
df2,
aes(x = dpos, color = vaxx_status_label, fill = vaxx_status_label)
) +
geom_ribbon(
aes(
ymin = rna_load_ml_mean_lower,
ymax = rna_load_ml_mean_upper,
x = dpos
),
alpha = 0.2
) +
geom_line(aes(y = rna_load_ml_mean), linetype = "dashed") +
geom_point(aes(y = rna_load_ml), data = df) +
labs(
title = "RNA viral loads [MEAN]",
subtitle = paste0(
"Source: Puhach et al., 2022: ",
"https://www.nature.com/articles/s41591-022-01816-0"
),
x = "DPOS",
y = "genome copies per ml"
)
## Median
ggplot(
df2, aes(x = dpos, color = vaxx_status_label, fill = vaxx_status_label)
) +
geom_ribbon(
aes(
ymin = rna_load_ml_median_lower,
ymax = rna_load_ml_median_upper,
x = dpos
),
alpha = 0.2
) +
geom_line(aes(y = rna_load_ml_median), linetype = "dashed") +
geom_point(aes(y = rna_load_ml), data = df) +
labs(
title = "RNA viral loads [MEDIAN]",
subtitle = paste0(
"Source: Puhach et al., 2022: ",
"https://www.nature.com/articles/s41591-022-01816-0"
),
x = "DPOS",
y = "genome copies per ml"
)
# Infectious viral loads
## LOESS
ggplot(
df, aes(x = dpos, color = vaxx_status_label, fill = vaxx_status_label)
) +
geom_smooth(aes(y = rna_load_ml), method = "loess", se = TRUE) +
geom_point(aes(y = rna_load_ml)) +
scale_y_log10() +
labs(
title = "Infectious viral loads [LOESS]",
subtitle = paste0(
"Source: Puhach et al., 2022: ",
"https://www.nature.com/articles/s41591-022-01816-0"
),
x = "DPOS",
y = "log10 FFU/ml"
)
## MEAN
ggplot(
df2,
aes(x = dpos, color = vaxx_status_label, fill = vaxx_status_label)
) +
geom_ribbon(
aes(
ymin = ffu_ml_mean_lower,
ymax = ffu_ml_mean_upper,
x = dpos
),
alpha = 0.2
) +
geom_line(aes(y = ffu_ml_mean), linetype = "dashed") +
geom_point(aes(y = ffu_ml), data = df) +
labs(
title = "Infectious viral loads [MEAN]",
subtitle = paste0(
"Source: Puhach et al., 2022: ",
"https://www.nature.com/articles/s41591-022-01816-0"
),
x = "DPOS",
y = "FFU/ml"
)
## MEDIAN
ggplot(
df2,
aes(x = dpos, color = vaxx_status_label, fill = vaxx_status_label)
) +
geom_ribbon(
aes(
ymin = ffu_ml_median_lower,
ymax = ffu_ml_median_upper,
x = dpos
),
alpha = 0.2
) +
geom_line(aes(y = ffu_ml_median), linetype = "dashed") +
geom_point(aes(y = ffu_ml), data = df) +
labs(
title = "Infectious viral loads [MEAN]",
subtitle = paste0(
"Source: Puhach et al., 2022: ",
"https://www.nature.com/articles/s41591-022-01816-0"
),
x = "DPOS",
y = "FFU/ml"
)
# Plot by day and vaccination status
ggplot(df, aes(x = vaxx_status_label, y = rna_load_ml)) +
geom_boxplot() +
geom_dotplot(binaxis = "y", stackdir = "center", dotsize = 1) +
facet_wrap(vars(dpos), scales = "free") +
scale_y_log10() +
stat_summary(
fun.data = function(x) {
mean_cl_normal(x, conf.int = 0.999) # 99% confidence interval
},
geom = "errorbar",
width = 0.2,
color = "red"
) +
stat_summary(
fun = mean, # Add a point for the mean
geom = "point",
color = "blue",
size = 3
) +
labs(
title = "RNA viral loads",
subtitle = paste0(
"Median/Quartiles · Mean/99.9% CI · Source: Puhach et al., 2022: ",
"nature.com/articles/s41591-022-01816-0"
),
x = "DPOS",
y = "genome copies per ml"
)
ggplot(df, aes(x = vaxx_status_label, y = ffu_ml)) +
geom_boxplot() +
geom_dotplot(binaxis = "y", stackdir = "center", dotsize = 1) +
facet_wrap(vars(dpos), scales = "free") +
scale_y_log10() +
stat_summary(
fun.data = function(x) {
mean_cl_normal(x, conf.int = 0.999) # 99% confidence interval
},
geom = "errorbar",
width = 0.2,
color = "red"
) +
stat_summary(
fun = mean, # Add a point for the mean
geom = "point",
color = "blue",
size = 3
) +
labs(
title = "Infectious viral loads",
subtitle = paste0(
"Median/Quartiles · Mean/99.9% CI · Source: Puhach et al., 2022: ",
"nature.com/articles/s41591-022-01816-0"
),
x = "DPOS",
y = "log10 FFU/ml"
)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment