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

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  gganimate = list(
    nframes = 0
  ),
  out.width = "100%"
)
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)

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

Get variants data for India

gisaid_metadata <- qs::qread("~/data/epicov/metadata_tsv_2024_04_11.qs")
gisaid_india <- FilterGISAIDIndia(gisaid_metadata_all = gisaid_metadata)

vocs <- GetVOCs()
custom_voc_mapping <- list(
  `JN.1` = "JN.1",
  `JN.1.*` = "JN.1",
  `HV.1` = "HV.1",
  `HV.1.*` = "HV.1",
  `B.1` = "B.1",
  `B.1.1.306` = "B.1",
  `B.1.1.306.*` = "B.1",
  `B.1.1.326` = "B.1",
  `B.1.36.29` = "B.1",
  `B.1.560` = "B.1",
  `B.1.1` = "B.1",
  `B.1.210` = "B.1",
  `B.1.36.8` = "B.1",
  `B.1.36` = "B.1",
  `B.1.36.*` = "B.1"
)
gisaid_india <- gisaid_india %>%
  filter(pangolin_lineage != "None") %>%
  filter(pangolin_lineage != "Unassigned")


gisaid_india$District <- stringr::str_to_title(gisaid_india$District)
gisaid_india$City <- stringr::str_to_title(gisaid_india$City)

gisaid_india$custom_city <- gisaid_india$City
gisaid_india$custom_city[gisaid_india$custom_city == ""] <- gisaid_india$District[gisaid_india$custom_city == ""]

gisaid_india$custom_city <- stringr::str_to_title(gisaid_india$custom_city)

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

Distribution of variants

state_month_counts <- SummarizeVariantsMonthwise(gisaid_india)
state_month_counts$State <- "India"
state_month_prevalence <- CountsToPrevalence(state_month_counts)
state_month_prevalence <- CollapseLineageToVOCs(
  variant_df = state_month_prevalence,
  vocs = vocs,
  custom_voc_mapping = custom_voc_mapping, summarize = FALSE
)

p5 <- StackedBarPlotPrevalence(state_month_prevalence)
p5

Get weekly cases for India

GetIndiaCases <- function() {
  data <- read.csv("https://raw.githubusercontent.com/owid/covid-19-data/master/public/data/cases_deaths/new_cases.csv")
  confirmed <- data %>% select(date, India)
  colnames(confirmed)[2] <- c("cases")
  confirmed$MonthYear <- GetMonthYear(confirmed$date)
  confirmed$WeekYear <- tsibble::yearweek(confirmed$date)
  confirmed_subset_weekwise <- confirmed %>%
    group_by(WeekYear) %>%
    summarise(cases = mean(cases, na.rm = T)) %>%
    arrange(WeekYear)
  confirmed_subset_weekwise$cases <- ceiling(confirmed_subset_weekwise$cases)
  confirmed_subset_dateweekwise_long_india <- confirmed_subset_weekwise %>%
    rename(n = cases) %>%
    rename(WeekYearCollected = WeekYear)
}


confirmed_subset_dateweekwise_long_india <- GetIndiaCases()
gisaid_india_weekwise <- SummarizeVariantsWeekwise(gisaid_india)

Project weekly cases to variant prevalence data from GISAID

voc_to_keep <- gisaid_india_weekwise %>%
  group_by(lineage_collapsed) %>%
  summarise(n_sum = sum(n)) %>%
  filter(n_sum > 10) %>%
  pull(lineage_collapsed) %>%
  unique()
gisaid_india_weekwise <- gisaid_india_weekwise %>% filter(lineage_collapsed %in% voc_to_keep)

india_cases_pred_prob_sel_long <- FitMultinomWeekly(gisaid_india_weekwise, confirmed_subset_dateweekwise_long_india)
the_anim <- PlotVariantPrevalenceAnimated(india_cases_pred_prob_sel_long, title = "Estimated cases (weekly average) in India by variant", caption = "**Source: gisaid.org and ourworldindata.org/coronavirus**", date_breaks = "120 days")
gganimate::anim_save(filename = here::here("docs/articles/IN_animated.gif"), animation = the_anim)

Look at cases after January, 2022 only:

confirmed_subset_dateweekwise_long_india <- GetIndiaCases() %>%
  filter(WeekYearCollected >= tsibble::yearweek("2021 W35"))

gisaid_india <- gisaid_india %>% filter(MonthYearCollected > "Dec 2021")
gisaid_weekwise <- SummarizeVariantsWeekwise(gisaid_india)

voc_to_keep <- gisaid_weekwise %>%
  group_by(lineage_collapsed) %>%
  summarise(n_sum = sum(n)) %>%
  filter(n_sum > 10) %>%
  pull(lineage_collapsed) %>%
  unique()
gisaid_weekwise <- gisaid_weekwise %>% filter(lineage_collapsed %in% voc_to_keep)

cases_pred_prob_sel_long <- FitMultinomWeekly(gisaid_weekwise, confirmed_subset_dateweekwise_long_india)
the_anim <- PlotVariantPrevalenceAnimated(cases_pred_prob_sel_long, title = "Estimated cases (weekly average) in India by variant", caption = "**Source: gisaid.org and ourworldindata.org/coronavirus**<br>", date_breaks = "90 days") # , trans_y="log10")
gganimate::anim_save(filename = here::here("docs/articles/IN_animated_2021.gif"), animation = the_anim)

Look at cases in the last few months:

confirmed_subset_dateweekwise_long_india <- GetIndiaCases() %>%
  filter(WeekYearCollected >= tsibble::yearweek("2022 W12"))

gisaid_india <- gisaid_india %>% filter(MonthYearCollected >= "Mar 2022")
gisaid_weekwise <- SummarizeVariantsWeekwise(gisaid_india)

voc_to_keep <- gisaid_weekwise %>%
  group_by(lineage_collapsed) %>%
  summarise(n_sum = sum(n)) %>%
  filter(n_sum > 10) %>%
  pull(lineage_collapsed) %>%
  unique()
gisaid_weekwise <- gisaid_weekwise %>% filter(lineage_collapsed %in% voc_to_keep)

cases_pred_prob_sel_long <- FitMultinomWeekly(gisaid_weekwise, confirmed_subset_dateweekwise_long_india)
the_anim <- PlotVariantPrevalenceAnimated(cases_pred_prob_sel_long, title = "Estimated cases (weekly average) in India by variant", caption = "**Source: gisaid.org and ourworldindata.org/coronavirus**<br>", date_breaks = "90 days")
gganimate::anim_save(filename = here::here("docs/articles/IN_animated_2022.gif"), animation = the_anim)

Look at cases in the last few months:

confirmed_subset_dateweekwise_long_india <- GetIndiaCases() %>%
  filter(WeekYearCollected >= tsibble::yearweek("2022 W35"))

gisaid_india <- gisaid_india %>% filter(MonthYearCollected >= "Dec 2022")
gisaid_weekwise <- SummarizeVariantsWeekwise(gisaid_india)

voc_to_keep <- gisaid_weekwise %>%
  group_by(lineage_collapsed) %>%
  summarise(n_sum = sum(n)) %>%
  filter(n_sum > 10) %>%
  pull(lineage_collapsed) %>%
  unique()
gisaid_weekwise <- gisaid_weekwise %>% filter(lineage_collapsed %in% voc_to_keep)

cases_pred_prob_sel_long <- FitMultinomWeekly(gisaid_weekwise, confirmed_subset_dateweekwise_long_india)
the_anim <- PlotVariantPrevalenceAnimated(cases_pred_prob_sel_long, title = "Estimated cases (weekly average) in India by variant", caption = "**Source: gisaid.org and ourworldindata.org/coronavirus**<br>", date_breaks = "30 days")
gganimate::anim_save(filename = here::here("docs/articles/IN_animated_latest.gif"), animation = the_anim)

confirmed_subset_dateweekwise_long_india <- GetIndiaCases() %>%
  filter(WeekYearCollected >= tsibble::yearweek("2023 W23"))

gisaid_india <- gisaid_india %>% filter(MonthYearCollected >= "June 2023")
gisaid_weekwise <- SummarizeVariantsWeekwise(gisaid_india)

voc_to_keep <- gisaid_weekwise %>%
  group_by(lineage_collapsed) %>%
  summarise(n_sum = sum(n)) %>%
  filter(n_sum > 10) %>%
  pull(lineage_collapsed) %>%
  unique()
gisaid_weekwise <- gisaid_weekwise %>% filter(lineage_collapsed %in% voc_to_keep)

cases_pred_prob_sel_long <- FitMultinomWeekly(gisaid_weekwise, confirmed_subset_dateweekwise_long_india)
the_anim <- PlotVariantPrevalenceAnimated(cases_pred_prob_sel_long, title = "Estimated cases (weekly average) in India by variant", caption = "**Source: gisaid.org and ourworldindata.org/coronavirus**<br>", date_breaks = "30 days")
gganimate::anim_save(filename = here::here("docs/articles/IN_animated_2023.gif"), animation = the_anim)

JN.1 variant

How many cases of JN.1 variant have been deposited to GISAID across the states?

jn.1 <- gisaid_india_all %>% filter(lineage_collapsed %in% c("JN.1"))
jn.1.grouped <- jn.1 %>%
  group_by(State) %>%
  tally()
jn.1.grouped <- jn.1.grouped %>%
  filter(!State %in% c("Unknown", "Unassigned")) %>%
  arrange(desc(n))
jn.1.grouped$State <- factor(jn.1.grouped$State, levels = jn.1.grouped$State)
ggplot(jn.1.grouped, aes(State, n, label = n)) + # , color= "#4682b4"
  geom_col(position = "identity", fill = "#4682b4") +
  geom_bar_text(stat = "identity") +
  xlab("") +
  ylab("Number of JN.1 samples deposited") +
  scale_x_discrete(guide = guide_axis(angle = 90)) +
  labs(
    title = "Number of JN.1 samples deposited to GISAID from India",
    caption = paste0("**Source: gisaid.org** <br>", Sys.Date())
  )

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