In pratikunterwegs/patho-move-evol: Analysis code for the _pathomove_ model.

Plot social-spatial structure

library(data.table)
library(glue)
library(tidygraph)
library(ggraph)
library(colorspace)
library(patchwork)

Load networks from default scenario

params <- list.files(
  "data/parameters",
  pattern = "default", full.names = T
)
params <- lapply(params, fread) |>
  rbindlist(use.names = TRUE)

params <- params[scenario == 1 & ((!infect_percent) &
  (costInfect == 0.25) &
  (regen_time == 50) &
  (dispersal == 2))]

data_files <- glue("data/results/networks/data_networks_{params$seed}.Rds")

# read in networks
data_ntwk <- lapply(data_files, function(file) {
  nt <- readRDS(file)
  nt
})

Plot degree distribution

Prepare degree data

# prepare generations for degree distributions
gen_patho_intro <- as.character(3000)
gen_patho_adapt <- as.character(3500)

# calculate centrality degree
degree_data <- lapply(
  data_ntwk, function(le) {
    le_pre <- le[[gen_patho_intro]]
    le_post <- le[[gen_patho_adapt]]

    # pre pathogen degree
    le_pre <- mutate(
      le_pre,
      degree = tidygraph::centrality_degree(
        normalized = F
      ),
      type = "pre"
    ) |>
      as_tibble()

    # post pathogen degree
    le_post <- mutate(
      le_post,
      degree = tidygraph::centrality_degree(
        normalized = F
      ),
      type = "post"
    ) |>
      as_tibble()

    rbindlist(
      list(
        le_pre,
        le_post
      )
    )
  }
) |>
  rbindlist()

# save degree data
fwrite(
  degree_data,
  "data/results/data_default_degree_distribution.csv"
)

# set the factor levels
degree_data <- fread("data/results/data_default_degree_distribution.csv")
degree_data$type <- factor(degree_data$type, levels = c("pre", "post"))

Prepare degree distribution data

# popsize per sim
popsize <- 500
# how many replicates
replicates <- 10

# plot degree distribution data
plot_degree <-
  ggplot(degree_data) +
  geom_histogram(
    aes(
      degree,
      fill = type,
      y = after_stat(count) / (popsize * replicates)
    ),
    position = "identity",
    bins = 25, alpha = 0.9,
    col = NA,
    show.legend = F
  ) +
  facet_grid(
    cols = vars(type),
    labeller = labeller(
      type = c(
        "pre" = "Pre-introduction",
        "post" = "Post-introduction"
      )
    )
  ) +
  scale_fill_discrete_diverging(
    palette = "Blue-Red",
    rev = F,
    l1 = 50,
    name = NULL,
    labels = c(
      "Pathogen-naive (G = 500)",
      "Pathogen-adapted (G = 700)"
    )
  ) +
  scale_x_continuous(
    breaks = c(0, 10, 50, 500),
    labels = function(x) {
      scales::percent(
        accuracy = 1,
        as.numeric(x) / popsize
      )
    },
    name = "% Pop. encountered"
  ) +
  scale_y_continuous(
    labels = scales::percent_format(accuracy = 1)
  ) +
  theme_test(
    base_family = "Arial",
    base_size = 8
  ) +
  theme(
    legend.position = "top",
    legend.key.height = unit(0.5, "mm"),
    legend.key.width = unit(2, "mm"),
    axis.text.y = element_text(
      angle = 90,
      hjust = 0.5
    ),
    plot.background = element_blank(),
    panel.background = element_blank()
  ) +
  labs(
    y = "% Indiv."
  ) +
  guides(
    fill = guide_legend()
  )

Prepare networks data

# select a replicate
repl <- 8
# select nice network
ntwks_example <- data_ntwk[[repl]]

# get the landscape
landscape <- list.files(
  path = "data/output", pattern = as.character(params$seed[repl]),
  full.names = TRUE
) |> readRDS()
landscape <- landscape$output@landscape

# select before and after disease
ntwk_pre <- ntwks_example[[gen_patho_intro]] %>%
  activate(edges) %>%
  filter(weight > quantile(weight, probs = 0.33)) %>%
  activate(nodes)

# select a nice network
# networks are plotted for illustration only
ntwks_example <- data_ntwk[[repl]]
ntwk_post <- ntwks_example[[gen_patho_adapt]] %>%
  activate(edges) %>%
  filter(weight > quantile(weight, probs = 0.33)) %>%
  activate(nodes)

# sanity check
ggraph(ntwk_pre, x = xn, y = yn) +
  geom_edge_fan(
    edge_width = 0.1
  ) +
  geom_node_point(
    aes(
      fill = t_infec,
      size = assoc
    ),
    shape = 21,
    show.legend = F
  ) +
  coord_equal(
    xlim = c(0, 60),
    ylim = c(0, 60)
  )

Plot network data

# make network figures
networkplots <- lapply(
  list(ntwk_pre, ntwk_post), function(n) {
    ggraph(n, x = xn, y = yn) +
      geom_point(
        data = landscape,
        aes(x, y, col = "food"),
        size = 0.05, alpha = 0.5
      ) +
      geom_edge_fan(
        edge_width = 0.1,
        edge_color = "grey",
        show.legend = FALSE
      ) +
      geom_node_point(
        aes(
          fill = t_infec,
          size = assoc
        ),
        shape = 21,
        colour = "grey40",
        show.legend = T
      ) +
      scale_radius(
        range = c(0.5, 5)
      ) +
      scale_fill_continuous_sequential(
        palette = "Rocket",
        limit = c(1, 101),
        breaks = c(1, 50, 100),
        na.value = "lightblue"
      ) +
      scale_colour_manual(
        values = c(
          food = "forestgreen"
        ),
        name = NULL,
        labels = "Food item locations"
      ) +
      coord_cartesian(
        expand = TRUE,
        xlim = c(0, 60),
        ylim = c(0, 60)
      ) +
      theme_graph(
        base_family = "Arial",
        background = "white",
        border = T,
        base_size = 8,
        plot_margin = margin(rep(0, 3))
      ) +
      theme(
        legend.margin = margin(rep(0, 4)),
        legend.position = "top",
        legend.title = element_text(size = 6),
        legend.key.height = unit(1, units = "mm"),
        legend.key.width = unit(3, units = "mm"),
        plot.background = element_blank()
      ) +
      labs(
        fill = "Time infected"
      ) +
      guides(
        size = "none",
        edge_alpha = "none",
        colour = guide_legend(
          override.aes = list(
            size = 1,
            shape = 16,
            colour = "forestgreen"
          )
        )
      )
  }
)

Prepare intermediate networks plot

# wrap plots
plot_networks <-
  wrap_plots(networkplots, guides = "collect", ncol = 2) &
    plot_annotation(
      tag_levels = c("A")
    ) &
    theme(
      plot.tag = element_text(
        face = "bold"
      ),
      legend.position = "bottom"
    )

Load and plot transmission chains

Load chain data

# read data
repl <- 9 # pick nice replicate
chains <- list.files(
  "data/results/transmission_chains",
  pattern = as.character(params$seed[[repl]]),
  full.names = TRUE
) |> readRDS()

# set factor levels
chains <- lapply(chains, function(g) {
  g |>
    activate(nodes) |>
    mutate(
      social_strat = factor(
        social_strat,
        levels = c("agent avoiding", "handler tracking", "agent tracking")
      )
    )
})

Plot and save transmission chains

plot_chains <- (lapply(chains, function(g) {
  ggraph(g, layout = "circlepack") +
    geom_node_circle(
      aes(fill = social_strat),
      colour = "grey40"
    ) +
    scale_fill_discrete_sequential(
      palette = "Viridis",
      l1 = 15, l2 = 80,
      rev = F,
      limits = c("agent avoiding", "handler tracking", "agent tracking"),
      order = c(1, 3, 2),
      name = NULL,
      na.value = "grey",
      labels = stringr::str_to_sentence
    ) +
    theme_graph(
      base_family = "Arial",
      background = "white",
      border = T,
      base_size = 8,
      plot_margin = margin(rep(0, 3))
    ) +
    theme(
      legend.margin = margin(rep(0, 4)),
      legend.position = "top",
      legend.title = element_text(size = 6),
      legend.key.height = unit(1, units = "mm"),
      legend.key.width = unit(3, units = "mm"),
      plot.background = element_blank()
    ) +
    labs(
      fill = "Time infected"
    ) +
    coord_equal(
      xlim = c(-21, 21),
      ylim = c(-21, 21)
    )
}) |> wrap_plots(guides = "collect", ncol = 2) &
  plot_annotation(
    tag_levels = c("A")
  ) &
  theme(
    plot.tag = element_text(
      face = "bold"
    ),
    legend.position = "bottom"
  ))

ggsave(
  plot_chains,
  filename = "supplement/figures/fig_default_chains.png",
  width = 160, height = 80, units = "mm"
)

Plot chain size

chain_size_data <- fread("data/results/transmission_data.csv")
chain_size_data <- chain_size_data[scenario_tag == "default", ]

# set factor levels
chain_size_data[, stage := factor(stage, levels = c("pre", "post"))]

Fit distributions on data

# prepare for fitdistrplus
chain_size_summary <- chain_size_data[, c(
  "stage", "n_infected", "mean_freq",
  "replicate"
)]
chain_size_summary <- split(chain_size_summary, by = "stage")

chain_size_summary <- lapply(
  chain_size_summary, function(df) {
    v <- rep(df$n_infected, times = df$mean_freq)

    fitdistrplus::fitdist(v, distr = "nbinom", discrete = TRUE)
  }
)

plot_chain_size <-
  ggplot(chain_size_data) +
  stat_summary(
    aes(
      x = as.integer(n_infected), y = mean_freq,
      fill = as.factor(stage)
    ),
    geom = "col",
    show.legend = FALSE
  ) +
  scale_fill_viridis_d(
    option = "F",
    begin = 0.5, end = 0.2,
    # l1 = 50, l2 = 70,
    name = NULL,
    labels = c(
      "Pathogen-naive (G = 3000)",
      "Pathogen-adapted (G = 3500)"
    )
  ) +
  facet_grid(
    cols = vars(stage),
    labeller = labeller(
      stage = c(
        "pre" = "3,000 ≤ G ≤ 3,100",
        "post" = "3,500 ≤ G ≤ 3,600"
      )
    )
  ) +
  theme_test(
    base_family = "Arial",
    base_size = 8
  ) +
  theme(
    legend.position = "top",
    legend.key.height = unit(0.5, "mm"),
    legend.key.width = unit(2, "mm"),
    axis.text.y = element_text(
      angle = 90,
      hjust = 0.5
    ),
    plot.background = element_blank(),
    panel.background = element_blank()
  ) +
  labs(
    x = "Secondary infections",
    y = "Mean count"
  ) +
  guides(
    fill = guide_legend()
  )

Save figure: networks and sociality metrics

plot_distributions <- wrap_plots(
  plot_degree, plot_chain_size,
  guides = "collect"
)

plot_sociality <-
  wrap_plots(
    plot_networks,
    plot_distributions,
    design = "AA\nAA\nAA\nBB"
  ) &
    plot_annotation(
      tag_levels = c("A")
    ) &
    theme(
      plot.tag = element_text(
        face = "bold"
      ),
      legend.position = "bottom"
    )

# save plot
ggsave(
  plot = plot_sociality,
  filename = "figures/fig_networks.png",
  height = 120,
  width = 150,
  units = "mm"
)

pratikunterwegs/patho-move-evol documentation built on April 15, 2023, 5:54 p.m.