Reanalysis Puhach et al, 2022

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"
  )