This repository hold the main scripts and SEIR model configurations for all analysis done in the paper:.
"Full genome viral sequences inform patterns of SARS-CoV-2 spear into and within Israel".

The phylodynamic analyses presented in this report depend on having BEAST2.6.2 (https://github.com/CompEvol/beast2) and PhyDyn 1.3.6 (https://github.com/mrc-ide/PhyDyn) installed. It also depends on having already generated a ML starting tree (data/tree.nwk) and an alignment (data/masked.fasta) of your sequences. In the future we plan to extend these scripts to allow more flexibility in the generation of the Beast XML file, particularly the use of random starting trees. The file scripts/process_fasta.py takes in the starting tree and alignment and labels sequence names in the tree and alignment with date of sampling and whether the sequences are considered part of the global reservoir of sequences ("exog") or whether they are sequences from Israel ("Il" or "Ih"). Sequences from Israel are randomly assigned to "Ih" with probability p_h, else they are assigned to "Il." This script then calls the scripts/generate_xml.py file which inserts your model parameters into the XML template (config/SEIR_TEMPLATE.xml). This XML template was originally generated using BEAUti (https://github.com/CompEvol/beast2) but has been extensively edited by hand. As it is set up right now, the script generates 7 XML files across a range of p_h values. These XML files are all saved in subdirectories labelled by the p_h value which are all within a directory called beast. The XML files which we used are not currently provided in this repository, however, final XML files upon journal submission will be available, less the alignments due to GISAID rules.

Once you have your BEAST2 XML files generated, they are ready to be run with beast (beast beast/{ph}/ph{ph}.xml). By default they are set to run for 10,000,000 MCMC steps. Once you have your Beast log and PhyDyn trajectory files, these can be analyzed using the scripts/combine_plot_traj.py file which will generate an MCC tree (using TreeAnnotator ((https://github.com/CompEvol/beast2))) and plot both the tree as well as PhyDyn trajectories for each p_h value. This script also contains a function which calls LogCombiner (https://github.com/CompEvol/beast2), however, this is not currently used as we only have one chain for each p_h value. Plotting of the trees relies on Baltic (https://github.com/evogytis/baltic).