Skyline-plots | past population dynamics using Bayesian Coalescent Skyline | Machine Learning library

While the coalescent-based Bayesian Skyline Plot is integrated in the BEAST2 core, we need to install the BDSKY package, which contains the Birth-Death Skyline model. Installation of packages is done using the package manager, which is integrated into BEAUti. Open the BEAST2 Package Manager by navigating to File > Manage Packages. Install the BDSKY package by selecting it and clicking the Install/Upgrade button (Figure 2). After the installation of a package, the program is on your computer, but BEAUti is unable to load the template files for the newly installed model unless it is restarted. So, let's restart BEAUti to make sure we have the BDSKY model at hand. Close the BEAST2 Package Manager and restart BEAUti to fully load the BDSKY package.
To start we have to import the alignment into BEAUti. In the Partitions panel, import the nexus file with the alignment by navigating to File > Import Alignment in the menu and then finding the hcv.nexus file on your computer or simply drag and drop the file into the BEAUti window. BEAUti will recognize the sequences from the *.nexus file as nucleotide data. It will do so for sequence files with the character set of A C G T N, where N indicates an unknown nucleotide. As soon as other non-gap characters are included (e.g. using R or Y to indicate purines and pyramidines) BEAUti will not recognize the data as nucleotides anymore (unless the type of data is specified in the *.nexus file) and open a dialogue box to confirm the data type. The sequences were all sampled in 1993 so we are dealing with a homochronous alignment and do not need to specify tip dates. Skip the Tip Dates panel and navigate to the Site Model panel. The next step is to specify the model of nucleotide evolution (the site model). We will be using the GTR model, which is the most general reversible model and estimates transition probabilities between individual nucleotides separately. That means that the transition probabilities between e.g. A and T will be inferred separately to the ones between A and C, however transition probabilities from A to C will be the same as C to A etc. Additionally, we allow for rate heterogeneity among sites. We do this by changing the Gamma Category Count to 4 (normally between 4 and 6). Change the Gamma Category Count to 4, make sure that the estimate box next to the Shape parameter of the Gamma distribution is ticked and set Subst Model to GTR. Make sure that the estimate box is ticked for all but one of the 6 rates (there should be 5 ticked boxes) and that Frequencies are estimated (Figure 3). Topic for discussion: Why are only 5 of the 6 rates of the GTR model estimated?. Because our sequences are contemporaneous (homochronous data) there is no information in our dataset to estimate the clock rate (for more information on this refer to the prior-selection tutorial) and we have to use external information to calibrate the clock. We will use an estimate inferred in {% cite Pybus2001 --file Skyline-plots/master-refs %} to fix the clock rate. In this case all the samples were contemporaneous (sampled at the same time) and the clock rate is simply a scaling of the estimated tree branch lengths (in substitutions/site) into calendar time. Navigate to the Clock Model panel. Leave the clock model as a Strict Clock and set Clock.rate to 0.00079 s/s/y (Figure 4). (Note that BEAUti is smart enough to know that the clock rate cannot be estimated on this dataset and grays out the estimate checkbox). Now we are ready to set up the Coalescent Bayesian Skyline as a tree-prior. Navigate to the Priors panel and select Coalescent Bayesian Skyline as the tree prior (Figure 5). The Coalescent Bayesian Skyline divides the time between the present and the root of the tree (the tMRCA) into segments, and estimates a different effective population size ({% eqinline N_e %}) for each segment. The endpoints of segments are tied to the branching times (also called coalescent events) in the tree (Figure 6), and the size of segments is measured in the number of coalescent events included in each segment. The Coalescent Bayesian Skyline groups coalescent events into segments and jointly estimates the {% eqinline N_e %} (bPopSizes parameter in BEAST) and the size of each segment (bGroupSizes parameter). To set the number of segments we have to change the dimension of bPopSizes and bGroupSizes (note that the dimension of both parameters always has to be the same). Note that the length of a segment is not fixed, but dependent on the timing of coalescent events in the tree (Figure 6), as well as the number of events contained within a segment (bGroupSizes). To change the number of segments we have to navigate to the Initialialization panel, which is by default not visible. Navigate to View > Show Initialization Panel to make it visible and navigate to it (Figure 7). Set the dimension of bPopSizes and bGroupSizes to 4 (the default value is 5) after expanding the boxes for the two parameters (Figure 8). This sets the number of segments equal to 4 (the parameter dimension), which means {% eqinline N_e %} will be allowed to change 3 times between the tMRCA and the present (if we have {% eqinline d %} segments, {% eqinline N_e %} is allowed to change {% eqinline d-1 %} times). We can leave the rest of the priors as they are and save the XML file. We want to shorten the chain length and decrease the sampling frequency so the analysis completes in a reasonable time and the output files stay small. (Keep in mind that it will be necessary to run a longer chain for parameters to mix properly). Navigate to the MCMC panel. Change the Chain Length from 10'000'000 to 3'000'000. Click on the arrow next to the tracelog and change the File Name to $(filebase).log and set the Log Every to 3'000. Click on the arrow next to the treelog and change the File Name to $(filebase)-$(tree).log and set the Log Every to 3'000. Leave all other settings at their default values and save the file as hcv_coal.xml. (Note that since BEAST 2.7 the filenames used here are the default filenames and should not need to be changed!). When we run the analysis $(filebase) in the name of the *.log and *.trees files will be replaced by the name of the XML file. This is a good idea, since it makes it easy to keep track of which XML files produced which output files. Now we are ready to run the analysis. Start BEAST2 and choose the file hcv_coal.xml. If you have BEAGLE installed tick the box to Use BEAGLE library if available, which will make the analysis run faster. Hit Run to start the analysis. The analysis will take about 10 minutes to complete. Read through the next section while waiting for your results or start preparing the XML file for the birth-death skyline analysis.
In the first analysis, we used the coalescent approach to estimate population dynamics. We now want to repeat the analysis using the Birth-Death Skyline model. We will use the same model setup as in the previous analysis and only change the tree prior. Restart BEAUti, load hcv.nexus as before and set up the same site and clock model as in the Coalescent Bayesian Skyline analysis. We will need to set the prior to Birth Death Skyline Contemporary, since the sequences were all sampled at the same point in time. For heterochronous data (sequences sampled at different times), we would use Birth Death Skyline Serial. As with the Coalescent Bayesian Skyline, we need to set the number of dimensions. Here we set the dimension for {% eqinline R_e %}, the effective reproduction number, which denotes the average number of secondary infections caused by an infected person at a given time during the epidemic, i.e. an {% eqinline R_e %} of 2 would mean that every infected person causes two new infections on average. In other words, an {% eqinline R_e %} above 1 means that the number of cases are increasing, therefore the disease will cause an exponentially growing epidemic, and an {% eqinline R_e %} below 1 means that the epidemic will die out. Navigate to the Priors panel and select Birth Death Skyline Contemporary as the tree prior (Figure 15). Then, click on the button where it says initial = [2.0] [0.0, Infinity] next to reproductiveNumber. A pop-up window will open which allows us to change the dimension of the parameter (Figure 16). In this case we will keep the default dimension of 10. Press OK to close the pop-up window. This means that {% eqinline R_e %} will be allowed to change at 9 equally spaced times between the origin of the epidemic and the present time. Choosing this dimension can again be arbitrary and may require the testing of a few different values. Too few intervals and not all rate shifts are captured. Too many intervals and the intervals may not contain enough information to infer parameters. (As with setting the dimension of the Coalescent Bayesian Skyline the dimension of {% eqinline R_e %} can also be set in the initialization panel). Besides {% eqinline R_e %} (reproductiveNumber), the Birth Death Skyline Contemporary model has 3 more parameters, becomeUninfectiousRate (the rate at which infected patients become uninfectious, {% eqinline \delta %}, through recovery, death or isolation), rho (the proportion of lineages sampled in the present, {% eqinline \rho %}) and the origin (the time at which the index case became infected, which is always earlier than the tMRCA of the tree). We may know some of these parameters from literature or be able to estimate them from external sources. For example, the average time that patients are able to transmit a disease is informative about the becomeUninfectiousRate. This prior knowledge we can incorporate in our analysis by setting appropriate priors for these parameters. We will use a lognormal prior for {% eqinline R_e %}. This is a good prior distribution to use for rates since it is always positive (a rate cannot be negative) and has a long tail defined over all positive numbers. The long tail allows arbitrarily high estimates of {% eqinline R_e %}, but does not place much weight on very high rates. This agrees with our prior knowledge about {% eqinline R_e %} (most diseases have an {% eqinline R_e %} between 1.2 and 5. Measles is one of the most infectious diseases we know about and has {% eqinline R_e \approx 18 %}). If an epidemic is neither growing or declining, it has an {% eqinline R_e %} of 1, which we will use as a null hypothesis, by setting a prior on {% eqinline R_e %} centered around 1 (we assume that if there isn't a strong signal in an interval for an epidemic to grow or decline that {% eqinline R_e = 1 %}, i.e. the epidemic size stays constant). Note that this prior is used for each of the {% eqinline R_e %} intervals (the Birth-Death Skyline assumes that {% eqinline R_e %} is independent in each of the intervals). Select a Log Normal distribution for the reproductiveNumber prior. Click on the arrow to the left of reproductiveNumber to open all the options for {% eqinline R_e %} settings Set M to 0, which results in a median of 1. We set the variance to 1.25, which places most weight below 7.82 (95% quantile). (Figure 17). For the becoming uninfectious rate we will again use a log normal prior. The inverse of the becoming uninfectious rate is the average infectious period. In some patients an HCV infection only lasts a few weeks, while in others it is a chronic infection lasting for many years. Setting {% eqinline M=0 %} and {% eqinline S=1.25 %} results in the same prior we used for the {% eqinline R_e %}. In terms of the becoming uninfectious rate, this translates to the 95% quantiles for the infectious period falling between 0.0862 years (31.5 days) and 11.59 years, with a median of 1 year. We will see later that there is a strong signal in the data for a longer becoming uninfectious period. Set the same prior for becomeUninfectiousRate as for reproductiveNumber (Log Normal, with M=0.0, S=1.25) (Figure 18). The sampling proportion, {% eqinline \rho %}, represents the proportion of HCV cases in Egypt in 1993 that are included in the analysis. In 1993 Egypt had a population of roughly 60 million people, and with a prevalence of at least 15% this translates into millions of cases, while we only have 63 sequences. We will use a beta distribution for the prior on {% eqinline \rho %}. Beta distributions are a very flexible class of distributions that are only defined between 0 and 1, making them ideal to use for proportions. Select a Beta distribution for the rho prior. Click on the arrow to the left of rho to open all the options for the prior settings. Alpha to 1 and Beta to 9999, reflecting our prior knowledge that our dataset represents only a miniscule fraction of cases (Figure 19). Finally, we need to set a prior for the origin of the epidemic. We will once again use a log normal distribution for this parameter. Note that the origin also has to be positive and needs to be bigger than the MRCA of the tree. We know that HCV has been circulating in Egypt for at least a hundred years, so we set a prior with a median value greater than 100. Set a Log Normal prior for origin with M = 5 and S = 0.5 (Figure 20), resulting in a median prior estimate for the origin of 148 years. The rest of the priors pertain to the site model parameters and we can leave them as they are. Navigate to the MCMC panel. Change the Chain Length from 10'000'000 to 3'000'000. Click on the arrow next to the tracelog and change the File Name to $(filebase).log and set the Log Every to 3'000. Click on the arrow next to the treelog and change the File Name to $(filebase)-$(tree).log and set the Log Every to 3'000. Leave all other settings at their default values and save the file as hcv_bdsky.xml. (Note that since BEAST 2.7 the filenames used here are the default filenames and should not need to be changed!). Now we are ready to run the analysis. Start BEAST2 and choose the file hcv_bdsky.xml. If you have BEAGLE installed tick the box to Use BEAGLE library if available, which will make the analysis run faster. Hit Run to start the analysis. Look at the topics for discussion below and read through the next section while waiting for the analysis to finish. Topics for discussion: We set a prior on {% eqinline R_e %} in the Birth-Death Skyline analysis, but did not set any prior for {% eqinline N_e %} in the Coalescent Bayesian Skyline analysis. Is there a prior on {% eqinline N_e %}? If so, what is it? We fixed the clock rate to an independent estimate and set a strict clock. If we had strong prior knowledge that there is substitution rate variation over time in the Egyptian HCV epidemic, could we use a relaxed clock here?.
We set a prior on {% eqinline R_e %} in the Birth-Death Skyline analysis, but did not set any prior for {% eqinline N_e %} in the Coalescent Bayesian Skyline analysis. Is there a prior on {% eqinline N_e %}? If so, what is it?
We fixed the clock rate to an independent estimate and set a strict clock. If we had strong prior knowledge that there is substitution rate variation over time in the Egyptian HCV epidemic, could we use a relaxed clock here?