Regional Seasonality

In phylodyn: Statistical Tools for Phylodynamics

This vignette shows a case study of seasonality in influenza data from nine geographic regions, using phylodyn. We start by loading the phylodyn package.

library(phylodyn)

In preparation, we aligned the sequences by region using the software MUSCLE, and inferred a maximum clade credibility genealogy for each region using the software BEAST[^1]. We packaged the genealogies in a list of phylo objects regional_flu, and we load it now.

[^1]: We infer the genealogy branch lengths in units of years using a strict molecular clock, a constant effective population size prior, and an HKY substitution model with the first two nucleotides of a codon sharing the same estimated transition matrix, while the third nucleotide's transition matrix is estimated separately.

data(regional_flu)

We supply region names, plotting parameters, and other statistics unique to each region.

mains <- list(USACanada = "USA / Canada", Europe = "Europe", NorthChina = "NorthChina",
              JapanKorea = "Japan / Korea", India = "India", SouthChina = "South China",
              SouthAmerica = "South America", SoutheastAsia = "Southeast Asia", Oceania = "Oceania")

axlabs <- list(x = seq(1, 0, by=-1/12),
              labs = c("Sep", "Oct", "Nov", "Dec", "Jan", "Feb", "Mar",
                       "Apr", "May", "Jun", "Jul", "Aug", "Sep"))
start <- 8/12

zero_dates <- list(USACanada = 2012.301, Europe = 2011.044, NorthChina = 2011.285,
                   JapanKorea = 2012.29, India = 2010.814, SouthChina = 2011.282,
                   SouthAmerica = 2011.518, SoutheastAsia = 2011.995, Oceania = 2010.964)
years = list(USACanada = 12, Europe = 11, NorthChina = 10,
             JapanKorea = 12, India = 10, SouthChina = 11,
             SouthAmerica = 11, SoutheastAsia = 11, Oceania = 10)

We select a region (more than one can be chosen) and use BNPR and BNPR_PS.

regions <- list("USACanada")

results <- list()
for (region in regions)
{
  results[[region]] <- list(bnpr    = BNPR(data = regional_flu[[region]], lengthout = 100),
                            bnpr_ps = BNPR_PS(data = regional_flu[[region]], lengthout = 100))
}

We plot with years overlaid to show seasonal patterns.

par(mfrow=c(length(regions), 2), cex=1.0, cex.lab=1.4, cex.main=1.5,
    oma=c(2,3,0,0)+0.1, mar=c(3,2,2,1), xpd=NA)
for (region in regions)
{
  plot_seasonality(BNPR_out = results[[region]]$bnpr, zero_date = zero_dates[[region]],
                   start = start, years = years[[region]], ylim = c(0.5, 15), log_y = TRUE,
                   main=mains[[region]], axlabs = axlabs, xlab="Time", 
                   col_years = rgb(0.943, 0.893, 0.769), col_mean = rgb(0.829, 0.680, 0.306),
                   ylab="Effective Population Size", heatmap_labels = TRUE,
                   heatmap_width = 7, heatmap_labels_side = "left", legend = "BNPR")
  plot_seasonality(BNPR_out = results[[region]]$bnpr_ps, zero_date = zero_dates[[region]],
                   start = start, years = years[[region]], ylim = c(0.5, 15), log_y = TRUE,
                   main=mains[[region]], axlabs = axlabs, xlab="Time",
                   col_years = rgb(0.777, 0.828, 0.943), col_mean = rgb(0.330, 0.484, 0.828),
                   ylab="", heatmap_labels = FALSE, heatmap_width = 7,
                   heatmap_labels_side = "left", legend = "BNPR-PS")
}

References

1. R. C. Edgar. MUSCLE: Multiple sequence alignment with high accuracy and high through- put. Nucleic Acids Research, 32:1792–1797, 2004.

2. A. J. Drummond, M. A. Suchard, D. Xie, and A. Rambaut. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Molecular Biology and Evolution, 29:1969–1973, 2012.

3. D. Zinder, T. Bedford, E. B. Baskerville, R. J. Woods, M. Roy, and M. Pascual. Seasonality in the migration and establishment of H3N2 influenza lineages with epidemic growth and decline. BMC Evolutionary Biology, 14(1):272, 2014.

4. M. D. Karcher, J. A. Palacios, T. Bedford, M. A. Suchard, and V. N. Minin. Quantifying and mitigating the effect of preferential sampling on phylodynamic inference. arXiv preprint arXiv:1510.00775, 2015.

phylodyn documentation built on May 29, 2017, 1:28 p.m.