Real Data Tutorial
This final tutorial will show you how to take some existing real data from an Ebola outbreak in Sierra Leone, and use EpiFusion. The materials for this tutorial are can be found in Tutorial_3 in the examples folder, along with an Rmd script that makes it easy to follow along with the tutorial. Unlike Tutorial 2, where we use the prepare_epifusion_data function of EpiFusionUtilities to create an XML file with our data, we're going to manually populate an XML template in the Tutorial_3 folder, for a few reasons. First of all it is a good final challenge to fully understand EpiFusion XML. Secondly, while we are working every day to improve EpiFusion's versatility for use on real datasets, anecdotally we have found that finding the right parameter combination takes some tinkering - so the example XML file we include has a reasonable prior specification already set out.
For this example we will use an existing MCC tree and case incidence that are publicly available from other studies of this outbreak.
The tree can be obtained from the ebov/space-time repository on GitHub, associated with Dellicour et. al's Nature Comms 'Phylodynamic assessment of intervention strategies for the West African Ebola virus outbreak': https://github.com/ebov/space-time/blob/master/Analyses/Phylogenetic/Makona_1610_cds_ig.MCC.tree
Then we can get a line-list of cases from Fang et al's 'Transmission dynamics of Ebola virus disease and intervention effectiveness in Sierra Leone': https://www.pnas.org/doi/suppl/10.1073/pnas.1518587113/suppl_file/pnas.1518587113.sd02.xlsx
Download the tree from the above link and place it into your working directory (which should be the Tutorial_3 folder). The tree is a large MCC tree from a BEAST analysis, with sequences from Guinea, Sierra Leone, and Liberia. Let's load it into our workspace, and parse the metadata from the tip labels into a data frame using the ape package.
library(ape)
library(stringr)
raw_tree <- read.nexus("Makona_1610_cds_ig.MCC.tree")
#Turn tree labels into metadata table
alllabels <- raw_tree$tip.label
sequence_metadata <- data.frame(matrix(0, ncol = 7, nrow = length(alllabels)))
for (i in 1:length(alllabels)) {
label <- str_remove_all(alllabels[i], "'")
row <- c(label, unlist(str_split(label, "\\|")))
sequence_metadata[i,] <- row
}
colnames(sequence_metadata) <- c("tip", "Pathogen", "ID", "Accession", "Country", "Region", "Date")
raw_tree$tip.label <- str_remove_all(raw_tree$tip.label, "'")
Let's plot the tree with tips coloured by country, to see if there's any geographical clustering.
library(ggtree)
colourednational <- ggtree(raw_tree) %<+% sequence_metadata +
geom_tippoint(aes(color = Country)) +
theme_tree2()
colourednational
We can see from this plot that there is a nicely structured clade within SLE (Sierra Leone) which we can subset for our analysis. The clade has some GIN singletons and lineages which we will also prune. To subset the clade we just need to find the ancestral node we'd like to pick from. We can do this by labelling the internal nodes and highlighting our chosen clad to make sure we've picked the correct internal node.
raw_tree$node.label <- seq(1, raw_tree$Nnode)
ggtree(raw_tree) %<+% sequence_metadata +
geom_tippoint(aes(color = Country)) +
geom_nodelab(size = 3) +
theme_tree2()
When we take a look here we can see that node 219 seems to be the ancestor of the big SLE clade. To confirm, we can highlight the clade descended from that node, and check that the clade is correctly highlighted.
highlightedclade <- ggtree(raw_tree) %<+% sequence_metadata +
geom_hilight(node = "219", fill = "yellow", alpha = 0.2, type = "rect") +
geom_tippoint(aes(color = Country)) +
theme_tree2()
highlightedclade
Now we know what node we want to subset from, we can do this (again using a function in ape) and prune the GIN singletons, yielding a tree with just Sierra Leone sequences we can use for our analysis.
library(tidyverse)
library(dplyr)
subset_sle_tree <- extract.clade(raw_tree, "219", root.edge = T)
sierra_leone_sequence_metadata <- sequence_metadata %>%
dplyr::filter(tip %in% subset_sle_tree$tip.label) %>%
dplyr::filter(Country == "SLE")
subset_sle_tree <- keep.tip(subset_sle_tree, sierra_leone_sequence_metadata$tip)
ggtree(subset_sle_tree, mrsd = max(sierra_leone_sequence_metadata$Date)) %<+% sierra_leone_sequence_metadata +
geom_tippoint() +
theme_tree2()
Now that we have our tree, we just need to adjust our tree labels so that the internal node and tip labels have the 'time during outbreak' inside square brackets. This can also be done with the EpiFusionUtilities functions, but here we will do it manually. We use the get_all_distances_to_root from the castor package. We're going to start modelling from the beginning of the root edge of our tree (which is about 7 days prior to our ancestral node).
library(castor)
SLE_epifusion_tree <- subset_sle_tree
distances <- get_all_distances_to_root(SLE_epifusion_tree) + rnorm(length(get_all_distances_to_root(SLE_epifusion_tree)), 0.0001, 0.0001)
node_distances <- distances[(length(SLE_epifusion_tree$tip.label)+1):(length(SLE_epifusion_tree$tip.label)+SLE_epifusion_tree$Nnode)]
tip_distances <- distances[1:length(SLE_epifusion_tree$tip.label)]
SLE_epifusion_tree$node.label <- paste0("X[", (node_distances + subset_sle_tree$root.edge)*365, "]")
SLE_epifusion_tree$tip.label <- paste0("X[", (tip_distances + subset_sle_tree$root.edge)*365, "]")
treestring <- write.tree(SLE_epifusion_tree)
In the above code, we copy our subsetted tree into a new tree object, SLE_epifusion_tree, and use the get_all_distances_to_root function to get the distances for all tree nodes to the root node (in years, because the branch lengths of the tree are in terms of years). We use the rnorm package to stagger the node times ever so slightly (by just under one hour), just to make sure there are no concurrent nodes which EpiFusion is not the biggest fan of. Then we assign new node labels with these distances, plus the root edge which is where we are beginning modelling from, and multiplying by 363 to get the times in terms of days. Finally we use write.tree to save the newick string of this tree to treestring.
We can also calculate the 'start date' from which we are modelling by subtracting the tree length from the latest tip date:
start_date <- as.Date(max(sierra_leone_sequence_metadata$Date)) - ((max(distances) + subset_sle_tree$root.edge)*365)
In Fang et. al (2016) a line list of confirmed and suspected cases is provided in the supplementary in .xlsx. Download this and open in excel, where you can see the confirmed and suspected cases are in different tabs within the spreadsheet. We will want to save the 'confirmed' tab in .csv format to proceed with the tutorial. The resulting csv is provided in the tutorial folder, but you can also do this preprocessing yourself by getting the data from the provided link.
Now that we have it in CSV format, we can read it into our R session and process it for EpiFusion. We'll aggregate it into weekly counts starting from our start_date.
caselinelist <- read.csv("pnas.1518587113.sd02.csv") %>%
mutate(Date = as.Date(Date.of.symptom.onset, format = "%d-%b-%y")) %>%
filter(Date < max(sierra_leone_sequence_metadata$Date)) %>%
mutate(Week = ceiling(as.numeric(Date - start_date)/7)) %>%
group_by(Week) %>%
mutate(Cases = n()) %>%
dplyr::select(Week, Cases) %>%
distinct() %>%
ungroup()
Let's take a quick look at what that looks like by plotting it with ggplot:
ggplot(caselinelist, aes(x = Week, y = Cases)) +
geom_bar(stat = "identity")
This looks reasonable. For our own convenience, lets print out the values and times in a format that will go into the <incidence> section of the EpiFusion XML:
paste0(caselinelist$Cases, collapse = " ")
paste0(caselinelist$Week*7, collapse = " ")
It's time to open the XML template provided with this tutorial - ebola_xml_template.xml. Lines 8-11 house the 'incidence' block where we can put our incidence information. The first line we printed out in the above code block should go in incidenceValues, and the second line can go in incidenceTimes. The resulting block should now look like this:
<incidence>
<incidenceVals>13 18 20 71 54 55 83 85 46 81 76 75 32 90 181 185 214 383 365 369 451 369 433 369 485 448 468 359 254 304 268 305 219 190 108 91 90 71 87 64 64 79 48 42 18 18 8 11 13 6 4 14 1 15 15 15 8 10 8 16 3 2 2 2 1 3 3</incidenceVals>
<incidenceTimes type="exact">63 70 77 84 91 98 105 112 119 126 133 140 147 154 161 168 175 182 189 196 203 210 217 224 231 238 245 252 259 266 273 280 287 294 301 308 315 322 329 336 343 350 357 364 371 378 385 392 399 406 413 420 427 434 441 448 455 462 469 476 483 490 497 504 525 532 539</incidenceTimes>
</incidence>
[1] "13 18 20 71 54 55 83 85 46 81 76 75 32 90 181 185 214 383 365 369 451 369 433 369 485 448 468 359 254 304 268 305 219 190 108 91 90 71 87 64 64 79 48 42 18 18 8 11 13 6 4 14 1 15 15 15 8 10 8 16 3 2 2 2 1 3 3"
[1] "63 70 77 84 91 98 105 112 119 126 133 140 147 154 161 168 175 182 189 196 203 210 217 224 231 238 245 252 259 266 273 280 287 294 301 308 315 322 329 336 343 350 357 364 371 378 385 392 399 406 413 420 427 434 441 448 455 462 469 476 483 490 497 504 525 532 539"
We can do the same with the tree string by placing it in the tree block.
print(treestring)
With our XML now created, we can run EpiFusion. As before, you can run it in your R session (if you have EpiFusionUtilities installed and loaded) using the run_epifusion function, or from the command line in a terminal window with a command in this format:
java -jar path/to/EpiFusion/executable.jar path/to/XML/input/file.xml
As aforementioned, this example may take some time to run as it is necessary to run for many MCMC steps to get reasonable ESS figures for the parameters. We are making progress every day to make EpiFusion more user friendly and robust on real datasets. When your analysis is complete, you can parse and plot the results with the same instructions from Tutorial 1 and Tutorial 2.