Working example giving rough outline of organism models.

example.py

@@ -0,0 +1,28 @@
+"""
+Example of using the stdpopsim library with msprime.
+"""
+import msprime
+import stdpopsim.h_sapiens as h_sap
+
+
+model = h_sap.models.GutenkunstThreePopOutOfAfrica()
+
+model.debug()
+
+# One sample each from YRI, CEU and CHB. There's no point in pushing
+# the sampling strategy into the model generation
+samples = [
+    msprime.Sample(population=0, time=0),
+    msprime.Sample(population=1, time=0),
+    msprime.Sample(population=2, time=0)]
+
+ts = msprime.simulate(
+    samples=samples,
+    length=h_sap.chr22.length,
+    recombination_rate=h_sap.chr22.mean_recombination_rate,
+    mutation_rate=h_sap.chr22.mean_mutation_rate,
+    **model.asdict())
+
+# print(ts.tables)
+print("simulated:", ts.num_trees, ts.num_sites)
+

stdpopsim/__init__.py

@@ -7,3 +7,5 @@
     __version__ = _version.version
 except ImportError:
     pass
+
+from . import h_sapiens

stdpopsim/h_sapiens/__init__.py

@@ -0,0 +1,31 @@
+"""
+Human models, recombination and mutation rates.
+"""
+
+from . import models
+
+
+# TODO this infrastructure should live somewhere else in the package
+# hierarchy so it can be reused across the different species.
+
+class Chromosome(object):
+
+    def __init__(self, length, mean_recombination_rate, mean_mutation_rate):
+        self.length = length
+        self.mean_recombination_rate = mean_recombination_rate
+        self.mean_mutation_rate = mean_mutation_rate
+
+    # TODO add methods to return recombination maps
+
+    # Add methods to print this out. __str__ should give a nice summary.
+
+
+# Define the chromosomes.
+
+chr22 = Chromosome(
+    length=50818468,  # Taken from wikipedia, but should really be based on GRCh38.
+    mean_mutation_rate=1e-8,  # WRONG!
+    mean_recombination_rate=1e-8)  # WRONG!
+
+
+# ETC

stdpopsim/h_sapiens/models.py

@@ -0,0 +1,99 @@
+"""
+Simulation models for Homo Sapiens.
+"""
+import math
+
+import msprime
+
+
+class Model(object):
+    """
+    Class representing a simulation model that can be run in msprime.
+    """
+    def __init__(self):
+        self.population_configurations = None
+        self.migration_matrix = None
+        self.demographic_events = None
+
+    def debug(self):
+        # Use the demography debugger to print out the demographic history
+        # that we have just described.
+        dd = msprime.DemographyDebugger(
+            population_configurations=self.population_configurations,
+            migration_matrix=self.migration_matrix,
+            demographic_events=self.demographic_events)
+        dd.print_history()
+
+    def asdict(self):
+        return {
+            "population_configurations": self.population_configurations,
+            "migration_matrix": self.migration_matrix,
+            "demographic_events": self.demographic_events}
+
+
+
+class GutenkunstThreePopOutOfAfrica(Model):
+    """
+    The three population Out-of-Africa model from Gutenkunst et al.
+
+    TODO:
+
+    Clearly document that the different population indexes are.
+    """
+
+    def __init__(self):
+
+        # First we set out the maximum likelihood values of the various parameters
+        # given in Table 1.
+        N_A = 7300
+        N_B = 2100
+        N_AF = 12300
+        N_EU0 = 1000
+        N_AS0 = 510
+        # Times are provided in years, so we convert into generations.
+        generation_time = 25
+        T_AF = 220e3 / generation_time
+        T_B = 140e3 / generation_time
+        T_EU_AS = 21.2e3 / generation_time
+        # We need to work out the starting (diploid) population sizes based on
+        # the growth rates provided for these two populations
+        r_EU = 0.004
+        r_AS = 0.0055
+        N_EU = N_EU0 / math.exp(-r_EU * T_EU_AS)
+        N_AS = N_AS0 / math.exp(-r_AS * T_EU_AS)
+        # Migration rates during the various epochs.
+        m_AF_B = 25e-5
+        m_AF_EU = 3e-5
+        m_AF_AS = 1.9e-5
+        m_EU_AS = 9.6e-5
+        # Population IDs correspond to their indexes in the population
+        # configuration array. Therefore, we have 0=YRI, 1=CEU and 2=CHB
+        # initially.
+        self.population_configurations = [
+            msprime.PopulationConfiguration(initial_size=N_AF),
+            msprime.PopulationConfiguration(initial_size=N_EU, growth_rate=r_EU),
+            msprime.PopulationConfiguration(initial_size=N_AS, growth_rate=r_AS)
+        ]
+        self.migration_matrix = [
+            [      0, m_AF_EU, m_AF_AS],  # noqa
+            [m_AF_EU,       0, m_EU_AS],  # noqa
+            [m_AF_AS, m_EU_AS,       0],  # noqa
+        ]
+        self.demographic_events = [
+            # CEU and CHB merge into B with rate changes at T_EU_AS
+            msprime.MassMigration(
+                time=T_EU_AS, source=2, destination=1, proportion=1.0),
+            msprime.MigrationRateChange(time=T_EU_AS, rate=0),
+            msprime.MigrationRateChange(
+                time=T_EU_AS, rate=m_AF_B, matrix_index=(0, 1)),
+            msprime.MigrationRateChange(
+                time=T_EU_AS, rate=m_AF_B, matrix_index=(1, 0)),
+            msprime.PopulationParametersChange(
+                time=T_EU_AS, initial_size=N_B, growth_rate=0, population_id=1),
+            # Population B merges into YRI at T_B
+            msprime.MassMigration(
+                time=T_B, source=1, destination=0, proportion=1.0),
+            # Size changes to N_A at T_AF
+            msprime.PopulationParametersChange(
+                time=T_AF, initial_size=N_A, population_id=0)
+        ]