diff --git a/Dockerfile b/Dockerfile index d4b7986..1725dac 100755 --- a/Dockerfile +++ b/Dockerfile @@ -7,7 +7,7 @@ USER root RUN apt-get update -y && apt-get install -y bc # Install treesimulator -RUN cd /usr/local/ && pip3 install --no-cache-dir treesimulator==0.1.22 +RUN cd /usr/local/ && pip3 install --no-cache-dir treesimulator==0.2.0 # Switch to your new user in the docker image USER evolbioinfo diff --git a/README.md b/README.md index 6421fcc..3bdc211 100755 --- a/README.md +++ b/README.md @@ -1,11 +1,41 @@ # treesimulator -Simulation of rooted phylogenetic trees under a given Multitype Birth–Death (MTBD) model. +Simulation of rooted phylogenetic trees under a given Multitype Birth–Death (MTBD) model, with or without partner notification (PN). + +[//]: # ([![DOI:10.1093/sysbio/syad059](https://zenodo.org/badge/DOI/10.1093/sysbio/syad059.svg)](https://doi.org/10.1093/sysbio/syad059)) +[![GitHub release](https://img.shields.io/github/v/release/evolbioinfo/treesimulator.svg)](https://github.com/evolbioinfo/treesimulator/releases) +[![PyPI version](https://badge.fury.io/py/treesimulator.svg)](https://pypi.org/project/treesimulator/) +[![PyPI downloads](https://shields.io/pypi/dm/treesimulator)](https://pypi.org/project/treesimulator/) +[![Docker pulls](https://img.shields.io/docker/pulls/evolbioinfo/treesimulator)](https://hub.docker.com/r/evolbioinfo/treesimulator/tags) + +## MTBD The MTBD models were introduced by Stadler & Bonhoeffer [[_Philos. Trans. R. Soc. B_ 2013]](https://royalsocietypublishing.org/doi/10.1098/rstb.2012.0198). +An MTBD model with m states has + +m(m-1) state transition rate parameters: +* μij -- transition rate from state i to state j (1 ≤ i, j ≤ m; i ≠ j), where μij ≥ 0 + +m2 transmission rate parameters: +* λij -- transmission rate from state i (donor) to state j (recipient) (1 ≤ i, j ≤ m), where λij ≥ 0 + +m removal (becoming non-infectious) rate parameters: +* ψi -- removal rate of state i (1 ≤ i ≤ m), where ψi ≥ 0 + +m sampling probability upon removal parameters: +* pi -- probability to sample the pathogen of an individual in state i upon removal (1 ≤ i ≤ m), where 0 < pi ≤ 1 + +## Partner Notification (PN) +Partner notification adds two parameters to the initial MTBD model: + +* υ -- probability to notify partners upon sampling +* φ -- notified partner removal and sampling rate: φ >> ψi ∀i (1 ≤ i ≤ m). The pathogen of a notified partner is sampled automatically (with a probability of 1) upon removal. + + We pay particular interest to the classical BD model, the BD Exposed-Infectious (BDEI) model, and BD with superspreading (BDSS), -as they are described in [[Voznica _et al._ 2021]]((https://www.biorxiv.org/content/10.1101/2021.03.11.435006v1)). +as they are described in [[Voznica _et al._ 2021]]((https://www.biorxiv.org/content/10.1101/2021.03.11.435006v1)), and to their -PN versions. + ## BD 3 parameters: @@ -17,8 +47,23 @@ Epidemiological parameters: * R0=λ/ψ -- reproduction number * 1/ψ -- infectious time + +## BD-PN +5 parameters: +* λ -- transmission rate +* ψ -- removal rate +* p -- sampling probability upon removal +* υ -- probability to notify partners upon sampling +* φ -- notified partner removal and sampling rate: φ >> ψ + + +Epidemiological parameters: +* R0=λ/ψ -- reproduction number +* 1/ψ -- infectious time +* 1/φ -- notified partner removal time + ## BDEI -2 compartments: +2 states: * E, exposed, i.e. infected but not yet infectious * I, infectious @@ -33,6 +78,26 @@ Epidemiological parameters: * 1/ψ -- infectious time * 1/μ -- incubation period + +## BDEI-PN +2 states: +* E, exposed, i.e. infected but not yet infectious +* I, infectious + +6 parameters: +* μ -- transition rate from E to I (becoming infectious) +* λ -- transmission rate from I to E +* ψ -- removal rate of I +* p -- sampling probability upon removal +* υ -- probability to notify partners upon sampling +* φ -- notified partner removal and sampling rate: φ >> ψ + +Epidemiological parameters: +* R0=λ/ψ -- reproduction number +* 1/ψ -- infectious time +* 1/μ -- incubation period +* 1/φ -- notified partner removal time + ## BDSS 2 compartments: * N, standard infectious individual @@ -48,104 +113,242 @@ Epidemiological parameters: * ψ -- removal rate of S and of N (the same) * p -- sampling probability upon removal (the same for N and S) +Epidemiological parameters: +* R0=(λnn + λss)/ψ -- reproduction number +* 1/ψ -- infectious time +* X=λss/λns=λsn/λnn -- super-spreading transmission ratio +* f=λss/(λsn + λss) -- super-spreading fraction + +## BDSS-PN +2 states: +* N, standard infectious individual +* S, superspreader +8 parameters: +* λnn -- transmission rate from N to N +* λns -- transmission rate from N to S +* λsn -- transmission rate from S to N +* λss -- transmission rate from S to S + + (with a constraint that λss/λns=λsn/λnn) +* ψ -- removal rate of S and of N (the same) +* p -- sampling probability upon removal (the same for N and S) +* υ -- probability to notify partners upon sampling +* φ -- notified partner removal and sampling rate: φ >> ψ Epidemiological parameters: * R0=(λnn + λss)/ψ -- reproduction number * 1/ψ -- infectious time * X=λss/λns=λsn/λnn -- super-spreading transmission ratio * f=λss/(λsn + λss) -- super-spreading fraction +* 1/φ -- notified partner removal time + ## Installation -To install treesimulator: + +There are 4 alternative ways to run __treesimulator__ on your computer: +with [docker](https://www.docker.com/community-edition), +[apptainer](https://apptainer.org/), +in Python3, or via command line (requires installation with Python3). + + +### Installation in python3 or command-line + +You could either install python (version 3.6 or higher) system-wide and then install treesimulator via pip: ```bash +sudo apt install -y python3 python3-pip python3-setuptools python3-distutils pip3 install treesimulator ``` -## Usage -### Command line +or alternatively, you could install python (version 3.6 or higher) and treesimulator via [conda](https://conda.io/docs/) (make sure that conda is installed first). -## BD -The following command simulates a tree with 200-500 tips under BD model, with λ=0.5, ψ=0.25, p=0.5, -and saves it to a file tree.nwk, while saving the parameters to a comma-separated file params.csv: +(Optional) to install treesimulator in a new conda environment (e.g., called _phyloenv_ below), first create and activate the environment: ```bash -generate_bd --min_tips 200 --max_tips 500 --la 0.5 --psi 0.25 --p 0.5 --nwk tree.nwk --log params.csv +conda create --name phyloenv python=3.6 +conda activate phyloenv ``` -To seee detailed options, run: + +Install treesimulator with conda +```bash +conda install treesimulator +``` + + +### Basic usage in a command line +If you installed __treesimulator__ in a conda environment (here named _phyloenv_), do not forget to first activate it, e.g. + +```bash +conda activate phyloenv +``` + + +#### BD and BD-PN +The following command simulates a tree with 200-500 tips under the BD model, with λ=0.5, ψ=0.25, p=0.5, +and saves it to the file tree.nwk, while saving the parameters to the comma-separated file params.csv: +```bash +generate_bd --min_tips 200 --max_tips 500 \ +--la 0.5 --psi 0.25 --p 0.5 \ +--nwk tree.nwk --log params.csv +``` +The following command simulates a tree with 200-500 tips under the BD-PN model, with λ=0.5, ψ=0.25, p=0.5, φ=2.5, υ=0.2, +and allowing to notify only the most recent partner of each sampled index case. +The simulated tree is saved to the file tree.nwk, while the model parameters are saved to the comma-separated file params.csv: +```bash +generate_bd --min_tips 200 --max_tips 500 \ +--la 0.5 --psi 0.25 --p 0.5 \ +--phi 2.5 --upsilon 0.2 --max_notified_partners 1 \ +--nwk tree.nwk --log params.csv +``` +To see detailed options, run: ```bash generate_bd --help ``` -## BDEI -The following command simulates a tree with 200-500 tips under BDEI model, with μ=1, λ=0.5, ψ=0.25, p=0.5, -and saves it to a file tree.nwk, while saving the parameters to a comma-separated file params.csv: +#### BDEI and BDEI-PN +The following command simulates a tree with 200-500 tips under the BDEI model, with μ=1, λ=0.5, ψ=0.25, p=0.5, +and saves it to the file tree.nwk, while saving the parameters to the comma-separated file params.csv: +```bash +generate_bdei --min_tips 200 --max_tips 500 \ +--mu 1 --la 0.5 --psi 0.25 --p 0.5 \ +--nwk tree.nwk --log params.csv +``` +The following command simulates a tree with 200-500 tips under the BDEI-PN model, with μ=1, λ=0.5, ψ=0.25, p=0.5, φ=2.5, υ=0.2, +and allowing to notify only the most recent partner of each sampled index case. +The simulated tree is saved to the file tree.nwk, while the model parameters are saved to the comma-separated file params.csv: ```bash -generate_bdei --min_tips 200 --max_tips 500 --mu 1 --la 0.5 --psi 0.25 --p 0.5 --nwk tree.nwk --log params.csv +generate_bdei --min_tips 200 --max_tips 500 \ +--mu 1 --la 0.5 --psi 0.25 --p 0.5 \ +--phi 2.5 --upsilon 0.2 --max_notified_partners 1 \ +--nwk tree.nwk --log params.csv ``` -To seee detailed options, run: +To see detailed options, run: ```bash generate_bdei --help ``` -## BDSS -The following command simulates a tree with 200-500 tips under BDSS model, +#### BDSS and BDSS-PN +The following command simulates a tree with 200-500 tips under the BDSS model, with λnn=0.1, λns=0.3, λsn=0.5, λss=1.5, ψ=0.25, p=0.5, -and saves it to a file tree.nwk, while saving the parameters to a comma-separated file params.csv: +and saves it to the file tree.nwk, while saving the parameters to the comma-separated file params.csv: +```bash +generate_bdss --min_tips 200 --max_tips 500 \ +--la_nn 0.1 --la_ns 0.3 --la_sn 0.5 --la_ss 1.5 --psi 0.25 --p 0.5 \ +--nwk tree.nwk --log params.csv +``` +The following command simulates a tree with 200-500 tips under the BDSS-PN model, +with λnn=0.1, λns=0.3, λsn=0.5, λss=1.5, ψ=0.25, p=0.5, φ=2.5, υ=0.2, +and allowing to notify only the most recent partner of each sampled index case. +The simulated tree is saved to the file tree.nwk, while the model parameters are saved to the comma-separated file params.csv: ```bash -generate_bdss --min_tips 200 --max_tips 500 --la_nn 0.1 --la_ns 0.3 --la_sn 0.5 --la_ss 1.5 --psi 0.25 --p 0.5 --nwk tree.nwk --log params.csv +generate_bdss --min_tips 200 --max_tips 500 \ +--la_nn 0.1 --la_ns 0.3 --la_sn 0.5 --la_ss 1.5 --psi 0.25 --p 0.5 \ +--phi 2.5 --upsilon 0.2 --max_notified_partners 1 \ +--nwk tree.nwk --log params.csv ``` -To seee detailed options, run: +To see detailed options, run: ```bash generate_bdss --help ``` -## User-defined MTBD model +#### User-defined MTBD and MTBD-PN models The following command simulates a tree with 200-500 tips under a generic MTBD model, with two states A and B, with μaa=0.5, μab=0.6, μba=0.7, μbb=0.8, λaa=0.1, λab=0.2, λba=0.3, λbb=0.4, ψa=0.05, ψb=0.08, p=a0.15, p=b0.65, -and saves it to a file tree.nwk, while saving the parameters to a comma-separated file params.csv: +and saves it to the file tree.nwk, while saving the parameters to the comma-separated file params.csv: ```bash -generate_mtbd --min_tips 200 --max_tips 500 --nwk tree.nwk --log params.csv \ +generate_mtbd --min_tips 200 --max_tips 500 \ --states A B \ --transition_rates 0.5 0.6 0.7 0.8 \ --transmission_rates 0.1 0.2 0.3 0.4 \ --removal_rates 0.05 0.08 \ ---sampling_probabilities 0.15 0.65 +--sampling_probabilities 0.15 0.65 \ +--nwk tree.nwk --log params.csv ``` -To seee detailed options, run: +The following command simulates a tree with 200-500 tips under a generic MTBD model, with two states A and B, +with μaa=0.5, μab=0.6, μba=0.7, μbb=0.8, +λaa=0.1, λab=0.2, λba=0.3, λbb=0.4, +ψa=0.05, ψb=0.08, +p=a0.15, p=b0.65, +φ=2.5, υ=0.2, +and allowing to notify only the most recent partner of each sampled index case. +The simulated tree is saved to the file tree.nwk, while the model parameters are saved to the comma-separated file params.csv: +```bash +generate_mtbd --min_tips 200 --max_tips 500 \ +--states A B \ +--transition_rates 0.5 0.6 0.7 0.8 \ +--transmission_rates 0.1 0.2 0.3 0.4 \ +--removal_rates 0.05 0.08 \ +--sampling_probabilities 0.15 0.65 \ +--phi 2.5 --upsilon 0.2 --max_notified_partners 1 \ +--nwk tree.nwk --log params.csv +``` +To see detailed options, run: ```bash generate_mtbd --help ``` -### Python3 +#### Basic usage in Python3 To simulate trees with 200-500 tips under the above models and settings: -```python +```python3 from treesimulator.generator import generate from treesimulator import save_forest -from treesimulator.mtbd_models import Model, BirthDeathModel, BirthDeathExposedInfectiousModel, BirthDeathWithSuperSpreadingModel +from treesimulator.mtbd_models import Model, BirthDeathModel, BirthDeathExposedInfectiousModel, \ + BirthDeathWithSuperSpreadingModel, PNModel +# BD and BD-PN bd_model = BirthDeathModel(p=0.5, la=0.5, psi=0.25) print(bd_model.get_epidemiological_parameters()) -[bd_tree], (bd_total_tips, _, bd_total_time) = generate(bd_model, min_tips=200, max_tips=500) +[bd_tree], _, _ = generate(bd_model, min_tips=200, max_tips=500) save_forest([bd_tree], 'BD_tree.nwk') +bdpn_model = PNModel(model=bd_model, upsilon=0.2, partner_removal_rate=2.5) +[bdpn_tree], _, _ = generate(bdpn_model, min_tips=200, max_tips=500) +save_forest([bdpn_tree], 'BDPN_tree.nwk') +# BDEI and BDEI-PN bdei_model = BirthDeathExposedInfectiousModel(p=0.5, mu=1, la=0.5, psi=0.25) print(bdei_model.get_epidemiological_parameters()) -[bdei_tree], (bdei_total_tips, _, bdei_total_time) = generate(bdei_model, min_tips=200, max_tips=500) +[bdei_tree], _, _ = generate(bdei_model, min_tips=200, max_tips=500) save_forest([bdei_tree], 'BDEI_tree.nwk') +bdeipn_model = PNModel(model=bdei_model, upsilon=0.2, partner_removal_rate=2.5) +[bdeipn_tree], _, _ = generate(bdeipn_model, min_tips=200, max_tips=500) +save_forest([bdeipn_tree], 'BDEIPN_tree.nwk') +# BDSS and BDSS-PN bdss_model = BirthDeathWithSuperSpreadingModel(p=0.5, la_nn=0.1, la_ns=0.3, la_sn=0.5, la_ss=1.5, psi=0.25) print(bdss_model.get_epidemiological_parameters()) -[bdss_tree], (bdss_total_tips, _, bdss_total_time) = generate(bdss_model, min_tips=200, max_tips=500) +[bdss_tree], _, _ = generate(bdss_model, min_tips=200, max_tips=500) save_forest([bdss_tree], 'BDSS_tree.nwk') +bdsspn_model = PNModel(model=bdss_model, upsilon=0.2, partner_removal_rate=2.5) +[bdsspn_tree], _, _ = generate(bdsspn_model, min_tips=200, max_tips=500) +save_forest([bdsspn_tree], 'BDSSPN_tree.nwk') +# MTBD and MTBD-PN mtbd_model = Model(states=['A', 'B'], transition_rates=[[0.5, 0.6], [0.7, 0.8]], transmission_rates=[[0.1, 0.2], [0.3, 0.4]], removal_rates=[0.05, 0.08], ps=[0.15, 0.65]) -[mtbd_tree], (mtbd_total_tips, _, mtbd_total_time) = generate(mtbd_model, min_tips=200, max_tips=500) +[mtbd_tree], _, _ = generate(mtbd_model, min_tips=200, max_tips=500) save_forest([mtbd_tree], 'MTBD_tree.nwk') -``` \ No newline at end of file +mtbdpn_model = PNModel(model=mtbd_model, upsilon=0.2, partner_removal_rate=2.5) +[mtbdpn_tree], _, _ = generate(mtbdpn_model, min_tips=200, max_tips=500) +save_forest([mtbdpn_tree], 'MTBDPN_tree.nwk') +``` + +### Run with apptainer + +Once [apptainer](https://apptainer.org/docs/user/latest/quick_start.html#installation) is installed, +run the following command: + +```bash +apptainer run docker://evolbioinfo/treesimlator +``` + +This will launch a terminal session within the container, +in which you can run treesimulator following the instructions for the command line ("Basic usage in a command line") above. + + + + diff --git a/setup.py b/setup.py index e1e07f2..7dc7121 100755 --- a/setup.py +++ b/setup.py @@ -13,7 +13,7 @@ 'Topic :: Software Development', 'Topic :: Software Development :: Libraries :: Python Modules', ], - version='0.1.22', + version='0.2.0', description='Simulation of rooted phylogenetic trees under a given Multitype Birth–Death model.', author='Anna Zhukova', author_email='anna.zhukova@pasteur.fr', diff --git a/treesimulator/__init__.py b/treesimulator/__init__.py index 2d9ad31..c72e900 100644 --- a/treesimulator/__init__.py +++ b/treesimulator/__init__.py @@ -1,7 +1,7 @@ import os STATE = 'state' -DIST_TO_START ='D' +DIST_TO_START = 'D' TIME_TILL_NOW = 'T' diff --git a/treesimulator/generator.py b/treesimulator/generator.py index a968633..b75b208 100755 --- a/treesimulator/generator.py +++ b/treesimulator/generator.py @@ -1,6 +1,6 @@ import logging import random -from collections import Counter, defaultdict +from collections import Counter import numpy as np from ete3 import TreeNode