In saketkc/covmuller: A tool to model COVID19 variant prevalence using data from GISAID

suppressPackageStartupMessages({
  library(covmuller)
  library(tidyverse)
  library(magick)
})
theme_set(CovmullerTheme())

Get data

current_date <- "2024_04_11"
fpath.tar <- paste0("~/data/epicov/metadata_tsv_", current_date, ".tar.xz")
fpath.qs <- paste0("~/data/epicov/metadata_tsv_", current_date, ".qs")
gisaid_metadata <- qs::qread(file = fpath.qs)
gisaid_metadata <- FormatGISAIDMetadata(gisaid_metadata) %>% filter(pangolin_lineage != "Unassigned")

vocs <- GetVOCs()
omicron <- vocs[["omicron"]]
vocs[["omicron"]] <- NULL
vocs[["omicron-others"]] <- omicron


custom_voc_mapping <- list(
  `BA\\.1` = "BA.1+",
  `BA\\.1\\.*` = "BA.1+",
  `BA\\.2\\.10` = "BA.2.10+",
  `BA\\.2\\.10\\.*` = "BA.2.10+",
  `BA\\.2\\.12` = "BA.2.12+",
  `BA\\.2\\.12\\.*` = "BA.2.12+",
  `BA\\.2` = "BA.2",
  `BA\\.3` = "BA.3",
  `BA\\.4` = "BA.4",
  `BA\\.5` = "BA.5"
)

month_prevalence <- CollapseLineageToVOCs(
  variant_df = gisaid_metadata,
  vocs = vocs,
  custom_voc_mapping = custom_voc_mapping, summarize = FALSE
)

confirmed_subset_dateweekwise_long <- GetIndiaConfirmedCasesWeeklyLong()

confirmed_subset_dateweekwise_long_india <- confirmed_subset_dateweekwise_long %>%
  filter(State == "India") %>%
  rename(n = value) %>%
  rename(WeekYearCollected = WeekYear) %>%
  dplyr::select(-contains("type"))
confirmed_subset_dateweekwise_long_india$State <- NULL


confirmed.tmp <- COVID19::covid19(country = "South Africa", verbose = FALSE, level = 2) %>%
  select(date, administrative_area_level_2, confirmed) %>%
  rename(State = administrative_area_level_2)

confirmed <- COVID19::covid19(country = "South Africa", verbose = FALSE) %>%
  select(date, confirmed) %>%
  filter(!is.na(confirmed))

confirmed$daily_cases <- c(confirmed$confirmed[1], diff(confirmed$confirmed))


confirmed$WeekYear <- tsibble::yearweek(confirmed$date)

confirmed_subset_dateweekwise_long_southafrica <- confirmed %>%
  group_by(WeekYear) %>%
  summarise(n = ceiling(mean(daily_cases, na.rm = T))) %>%
  arrange(WeekYear) %>%
  rename(WeekYearCollected = WeekYear)


gisaid_metadata_2022 <- gisaid_metadata %>%
  filter(MonthYearCollected >= "Dec 2021") %>%
  filter(pangolin_lineage != "") %>%
  filter(Country %in% c("South Africa", "India"))


gisaid_southafrica <- gisaid_metadata_2022 %>% filter(Country == "South Africa")
gisaid_southafrica$State <- "South Africa"
gisaid_southafrica$Country <- "X"


gisaid_india <- gisaid_metadata_2022 %>% filter(Country == "India")
gisaid_india$State <- "India"
gisaid_india$Country <- "X"

gisaid_sel <- rbind(gisaid_southafrica, gisaid_india) %>% arrange(State, MonthYearCollected)


gisaid_sel_collapsed <- CollapseLineageToVOCs(
  variant_df = gisaid_sel,
  vocs = vocs,
  custom_voc_mapping = custom_voc_mapping,
  summarize = FALSE
)


gisaid_southafrica_collapsed <- CollapseLineageToVOCs(
  variant_df = gisaid_southafrica,
  vocs = vocs,
  custom_voc_mapping = custom_voc_mapping,
  summarize = FALSE
)

gisaid_india_collapsed <- CollapseLineageToVOCs(
  variant_df = gisaid_india,
  vocs = vocs,
  custom_voc_mapping = custom_voc_mapping,
  summarize = FALSE
)

gisaid_sel_dateweek <- SummarizeVariantsDatewise(gisaid_sel_collapsed, by_state = TRUE)


fit <- FitMultinomStatewiseDaily(gisaid_sel_dateweek)
head(fit)

muller <- PlotMullerDailyPrevalence(fit)
muller

gisaid_southafrica <- gisaid_metadata_2022 %>% filter(Country == "South Africa")

gisaid_india <- gisaid_metadata_2022 %>% filter(Country == "India")

gisaid_southafrica_collapsed <- CollapseLineageToVOCs(
  variant_df = gisaid_southafrica,
  vocs = vocs,
  custom_voc_mapping = custom_voc_mapping,
  summarize = FALSE
)

gisaid_india_collapsed <- CollapseLineageToVOCs(
  variant_df = gisaid_india,
  vocs = vocs,
  custom_voc_mapping = custom_voc_mapping,
  summarize = FALSE
)

collapsed_unique_variants <- union(unique(gisaid_india_collapsed$lineage_collapsed), unique(gisaid_southafrica_collapsed$lineage_collapsed))

colors_to_use <- c("#ebac23", "#b80058", "#008cf9", "#006e00", "#00bbad", "#d163e6", "#b24502", "#ff9287", "#5954d6", "#00c6f8", "#878500", "#00a76c")
names(colors_to_use) <- collapsed_unique_variants

gisaid_southafrica_weekwise <- SummarizeVariantsWeekwise(gisaid_southafrica_collapsed)


gisaid_india_weekwise <- SummarizeVariantsWeekwise(gisaid_india_collapsed)

cases_southafrica <- confirmed_subset_dateweekwise_long_southafrica %>% filter(WeekYearCollected >= tsibble::yearweek("2021 W44"))

cases_india <- confirmed_subset_dateweekwise_long_india %>% filter(WeekYearCollected >= tsibble::yearweek("2021 W44"))

preds_southafrica <- FitMultinomWeekly(gisaid_southafrica_weekwise, cases_southafrica)

preds_india <- FitMultinomWeekly(gisaid_india_weekwise, cases_india)


sa_anim <- PlotVariantPrevalenceAnimated(preds_southafrica, title = "Estimated cases (weekly average) in South Africa by variant", caption = "**Source: gisaid.org and covid19datahub.io**<br>", colors = colors_to_use)
gganimate::anim_save(filename = here::here("docs/articles/SA_animated_compare.gif"), animation = sa_anim)


india_anim <- PlotVariantPrevalenceAnimated(preds_india, title = "Estimated cases (weekly average) in India by variant", caption = "**Source: gisaid.org and covid19bharat.org**<br>", colors = colors_to_use)
gganimate::anim_save(filename = here::here("docs/articles/IN_animated_compare.gif"), animation = india_anim)

india_mgif <- image_read(india_anim)
southafrica_mgif <- image_read(sa_anim)



final_gif <- image_append(c(india_mgif[1], southafrica_mgif[1]))
rm(india_mgif)
rm(southafrica_mgif)
for (i in 2:length(india_mgif)) {
  combined <- image_append(c(india_mgif[i], southafrica_mgif[i]))
  final_gif <- c(final_gif, combined)
}

gganimate::anim_save(filename = here::here("docs/articles/SA_IN_animated.gif"), animation = final_gif)

saketkc/covmuller documentation built on April 19, 2024, 10:14 a.m.