Phase-type models for tract length distributions#
In this tutorial, we illustrate how to use the functions in tracts to compute and depict tract length distributions under the different models of admixture introduced in the paper. These are, in increasing order of complexity:
The Monoecious (M) model (only defined for autosomal admixture),
The Dioecious-Coarse (DC) model,
The Dioecious-Fine (DF) model,
The hybrid-pedidree refinements of the DC (H-DC) and the DF (H-DF) models.
We start by importing the required python libraries:
[ ]:
from tracts.phase_type import PhTMonoecious, PhTDioecious
from tracts.phase_type import hybrid_pedigree as HP
from tracts.phase_type.plot import plot_tractlength_density, plot_tractlength_histogram
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams.update(mpl.rcParamsDefault)
Note
This notebook does not illustrate the implementation of the usual inference strategy in tracts (see the user guide for that). Instead, it illustrates the use of the phase_type module to compute tract length distributions from a given pair of migration matrices. This allows to explore the effect of different migration scenarios on the tract length distribution, and to eventually compare the results with those obtained from simulations. This notebook is intended for users who want to explore the inner workings of the tracts package (see the phase_type module documentation for more details).
The notebook will illustrate the implementation of the following diagram, for different admixture models and transmission scenarios:
Throughout this notebook, we consider two source populations and a continuous pulse from the first to the second. This can be represented with the following migration matrix. The whichpop parameter sets the population of interest (0 or 1), that is, from which the tract length distribution is to be computed.
[2]:
whichpop = 1
mig_matrix = np.array([[0, 0],[0.2, 0],[0.2, 0],[0.2, 0],[0.2, 0],[0.2, 0],[0, 1]])
In what follows, we compute Phase-Type densities and histograms under autosomal and X chromosome admixture, using all the models presented in tracts. We consider different levels of sex-bias, in both migration and recombination rates.
Note
The previous migration matrix represents a demographic model of 6 generations. This is a toy framework where all the presented models are computationally efficient and can be quickly computed and compared.
For more realistic demographic models going up to 15–18 generations, the following models should be used depending on the context:
For autosomal admixture: Dioecious-Coarse (DC) or Monoecious (M) model.
For X chromosome admixture: Dioecious-Coarse (DC) model or its Hybrid-Pedigree refinement (H-DC). Typically, the number of pedigree generations should be set to
TP=2.
An example under both settings is presented at the end of the notebook.
1. Choose transmission#
[1]:
X_chr = False # Set to True for admixture on the X chromosome
X_chr_male = False # Set to True if admixture is considered on the X chromosome of a male individual (only maternally inherited alleles).
[ ]:
which_L = 1.7928357829 # Second chromosome
#which_L = 1.96 # For the X chromosome
# Set a point grid on the finite chromosome
bins = np.linspace(0, which_L+0.1, num = 50)
2. Set sex-bias#
The following parameters set the sex-bias for migration and recombination rates. The user can modify these parameters to observe how differences appear between all the considered models.
The migration sex-bias parameter sex_bias ranges from -1 (exclusively male migration) to 1 (exclusively female migration), with 0 corresponding to unbiased migration – matching the sex-bias parametrization used throughout tracts (see the :doc:user guide </user_guide/index>).
[ ]:
sex_bias = 0 # Sex-bias parameter, from -1 (exclusively male migration) to 1 (exclusively female migration). 0 corresponds to unbiased migration.
mig_f = mig_matrix + sex_bias*np.minimum(mig_matrix, 1-mig_matrix) # Female-specific migration matrix
mig_m = mig_matrix - sex_bias*np.minimum(mig_matrix, 1-mig_matrix) # Male-specific migration matrix
mig_m[-1,:] = mig_matrix[-1,:]
mig_f[-1,:] = mig_matrix[-1,:]
rho_f = 1 # Female-specific recombination rate
rho_m = 1 # Male-specific recombination rate
3. Build Markov model#
For the Monoecious (M), Dioecious Fine (DF), and Dioecious Coarse (DC) models, we first compute a PhTMonoecious or a PhTDioecious object, and then compute the corresponding phase-type density or histogram using the tractlength_histogram_windowed() method. For the Pedigree-Dioecious models, densities and histograms are computed directly from the model parameters using the hybrid_pedigree_distribution() function.
3.1 Monoecious model (only for autosomal admixture)#
The Monoecious model is build using the PhTMonoecious class. For comparison with the other settings, the Monoecious model is always computed using the mean migration matrix and the mean recombination rate.
[ ]:
PTmodel_mono = PhTMonoecious(migration_matrix=mig_matrix,
rho = 0.5*(rho_f+rho_m)) # PhTMonoecious object
3.2 Dioecious Fine and Coarse models#
The Dioecious Fine and Coarse models are build using the PhTDioecious class. The Dioecious models are defined for both autosomal and X chromosome admixture, so we compute them using the sex-specific migration matrices and recombination rates.
[ ]:
# PhTDioecious objects (DF and DC models)
PTmodel_DF = PhTDioecious(migration_matrix_f=mig_f,
migration_matrix_m=mig_m,
rho_f=rho_f,
rho_m=rho_m,
sex_model='DF',
X_chromosome=X_chr,
X_chromosome_male=X_chr_male)
PTmodel_DC = PhTDioecious(migration_matrix_f=mig_f,
migration_matrix_m=mig_m,
rho_f=rho_f,
rho_m=rho_m,
sex_model='DC',
X_chromosome=X_chr,
X_chromosome_male=X_chr_male)
4. Compute tract length distributions#
4.1 Compute tract length density functions#
Density functions are computed using the tractlength_histogram_windowed() (resp. tractlength_histogram_windowed()) function of the PhTMonoecious (resp. PhTDioecious) class, setting density=True. To set a frequency scale, set freq=True.
[ ]:
# Phase-Type density for the Monoecious (M) model
newbins, density_m, E = PTmodel_mono.tractlength_histogram_windowed(population_number=whichpop,
bins=bins,
L=which_L,
density=True,
freq=False)
# Phase-Type density for the Dioecious Fine (DF) model
newbins, density_df, E = PTmodel_DF.tractlength_histogram_windowed(population_number=whichpop,
bins=bins,
L=which_L,
density=True,
freq=False)
# Phase-Type density for the Dioecious Coarse (DC) model
newbins, density_dc, E = PTmodel_DC.tractlength_histogram_windowed(population_number=whichpop,
bins=bins,
L=which_L,
density=True,
freq=False)
4.2 Compute tract length histograms#
Tract length histograms are computed using the tractlength_histogram_windowed() (resp. tractlength_histogram_windowed()) function of the PhTMonoecious (resp. PhTDioecious) class, setting density=False.
[ ]:
# Tract length histogram for the Monoecious (M) model
bins, hist_m, E = PTmodel_mono.tractlength_histogram_windowed(population_number=whichpop,
bins=bins,
L=which_L,
density=False)
# Tract length histogram for the Dioecious Fine (DF) model
bins, hist_df, E = PTmodel_DF.tractlength_histogram_windowed(population_number=whichpop,
bins=bins,
L=which_L,
density=False)
# Tract length histogram for the Dioecious Coarse (DC) model
bins, hist_dc, E = PTmodel_DC.tractlength_histogram_windowed(population_number=whichpop,
bins=bins,
L=which_L,
density=False)
5. Hybrid-pedigree refinements#
For the hybrid-pedigree refinements of the Dioecious models, we compute the tract length distributions directly using the hybrid_pedigree_distribution() function. As the corresponding distribution is a Phase-Type mixture, the function internally defines one PhTDioecious object for each component, and then appropriately combine them using the mixture weights.
5.1 Phase-type mixture density functions#
[ ]:
# Phase-Type mixture density for the H-DF model (set freq = True for frequency scale)
newbins, density_hp_df = HP.hybrid_pedigree_distribution(mig_matrix_f=mig_f,
mig_matrix_m=mig_m,
TP=2,
Dioecious_model='DF',
L=which_L,
bingrid=bins,
whichpop=whichpop,
rho_f=rho_f,
rho_m=rho_m,
X_chr=X_chr,
X_chr_male=X_chr_male,
N_cores=5,
density=True,
freq=False)
# Phase-Type mixture density for the H-DC model (set freq = True for frequency scale)
newbins, density_hp_dc = HP.hybrid_pedigree_distribution(mig_matrix_f=mig_f,
mig_matrix_m=mig_m,
TP=2,
Dioecious_model='DC',
L=which_L,
bingrid=bins,
whichpop=whichpop,
rho_f=rho_f,
rho_m=rho_m,
X_chr=X_chr,
X_chr_male=X_chr_male,
N_cores=5,
density=True,
freq=False)
5.2 Tract length histograms for phase-type mixtures#
[ ]:
# Phase-Type mixture histogram for the H-DF model
bins, hist_hp_df = HP.hybrid_pedigree_distribution(mig_matrix_f=mig_f,
mig_matrix_m=mig_m,
TP=2,
Dioecious_model='DF',
L=which_L,
bingrid=bins,
whichpop=whichpop,
rho_f=rho_f,
rho_m=rho_m,
X_chr=X_chr,
X_chr_male=X_chr_male,
N_cores=5,
density=False)
# Phase-Type mixture histogram for the H-DC model
bins, hist_hp_dc = HP.hybrid_pedigree_distribution(mig_matrix_f=mig_f,
mig_matrix_m=mig_m,
TP=2,
Dioecious_model='DC',
L=which_L,
bingrid=bins,
whichpop=whichpop,
rho_f=rho_f,
rho_m=rho_m,
X_chr=X_chr,
X_chr_male=X_chr_male,
N_cores=5,
density=False)
6. Plot tract length distributions#
6.1 Plot densities#
6.1.1 Plot for autosomal admixture#
The following code produces a figure depicting all the computed Phase-Type densities for autosomal admixture. Run if X_chr = False.
[9]:
fig, ax = plot_tractlength_density(bins=newbins,
curves={'DF': density_df, 'DC': density_dc, 'M': density_m, 'H_DF': density_hp_df, 'H_DC': density_hp_dc},
L=which_L,
xlabel='Tract length on the second chromosome')
plt.show()
6.1.2 Plot for X chromosome admixture#
The following code produces a figure depicting all the computed Phase-Type densities for admixture on the X chromosome. Run if X_chr = True.
[10]:
fig, ax = plot_tractlength_density(bins=newbins,
curves={'DF': density_df, 'DC': density_dc, 'H_DF': density_hp_df, 'H_DC': density_hp_dc},
L=which_L,
xlabel='Tract length on the X chromosome')
plt.show()
6.2 Plot histograms#
6.2.1 Plot for autosomal admixture#
The following code produces a figure depicting all the computed Phase-Type histograms for autosomal admixture. Run if X_chr = False.
[11]:
fig, ax = plot_tractlength_histogram(bins=bins,
curves={'DF': hist_df, 'DC': hist_dc, 'M': hist_m, 'H_DF': hist_hp_df, 'H_DC': hist_hp_dc},
xlabel='Tract length on the second chromosome',
L=which_L)
plt.show()
6.2.2 Plot for X chromosome admixture#
The following code produces a figure depicting all the computed Phase-Type histograms for admixture on the X chromosome. Run if X_chr = True.
[12]:
fig, ax = plot_tractlength_histogram(bins=bins,
curves={'DF': hist_df, 'DC': hist_dc, 'H_DF': hist_hp_df, 'H_DC': hist_hp_dc},
xlabel='Tract length on the X chromosome',
L=which_L)
plt.show()
7. Continuous pulse during 15 generations#
To conclude, we consider a more realistic scenario with the following migration matrix, letting the user specify the sex-bias parameters as in the previous example. Once again, whichpop sets the population of interest (0 or 1).
[13]:
whichpop = 1
mig_matrix = np.array([[0, 0],[0.2, 0],[0.2, 0],[0.2, 0],[0.2, 0],[0.2, 0],[0.2, 0],[0.2, 0],[0.2, 0],[0.2, 0],[0.2, 0],[0.2, 0],[0.2, 0],[0.2, 0],[0.2, 0],[0, 1]])
print(np.shape(mig_matrix)[0]-1) # Number of generations
[14]:
sex_bias = 0 # Sex-bias parameter, from -1 (exclusively male migration) to 1 (exclusively female migration). 0 corresponds to unbiased migration.
mig_f = mig_matrix + sex_bias*np.minimum(mig_matrix, 1-mig_matrix)
mig_m = mig_matrix - sex_bias*np.minimum(mig_matrix, 1-mig_matrix)
mig_m[-1,:] = mig_matrix[-1,:]
mig_f[-1,:] = mig_matrix[-1,:]
rho_f = 1 # Female-specific recombination rate
rho_m = 1 # Male-specific recombination rate
7.1 Monoecious model for autosomal admixture#
[15]:
PTmodel_mono = PhTMonoecious(migration_matrix=mig_matrix,
rho = 0.5*(rho_f+rho_m)) # PhTMonoecious object
# Phase-Type density (set freq = True for frequency scale)
newbins, density_m, E = PTmodel_mono.tractlength_histogram_windowed(population_number=whichpop,
bins=bins,
L=which_L,
density=True,
freq=False)
# Phase-type histogram
bins, hist_m, E = PTmodel_mono.tractlength_histogram_windowed(population_number=whichpop,
bins=bins,
L=which_L,
density=False,
freq=False)
7.1.1 Plot density#
[16]:
fig, ax = plot_tractlength_density(bins=newbins,
curves={'M': density_m},
L=which_L,
xlabel='Tract length on the second chromosome')
plt.show()
7.1.2 Plot histogram#
[17]:
fig, ax = plot_tractlength_histogram(bins=bins,
curves={'M': hist_m},
xlabel='Tract length on the second chromosome',
L=which_L)
plt.show()
7.2 H-DC model for X chromosome admixture#
[18]:
# The following lines can take a while to run, depending on the number of cores used and the size of the grid.
# Phase-Type mixture density for the H-DC model (set freq = True for frequency scale)
newbins, density_hp_dc = HP.hybrid_pedigree_distribution(mig_matrix_f=mig_f,
mig_matrix_m=mig_m,
TP=2,
Dioecious_model='DC',
L=which_L,
bingrid=bins,
whichpop=whichpop,
rho_f=rho_f,
rho_m=rho_m,
X_chr=True,
X_chr_male=False,
N_cores=5,
density=True,
freq=True)
# Phase-Type histogram for the H-DC model
bins, hist_hp_dc = HP.hybrid_pedigree_distribution(mig_matrix_f=mig_f,
mig_matrix_m=mig_m,
TP=2,
Dioecious_model='DC',
L=which_L,
bingrid=bins,
whichpop=whichpop,
rho_f=rho_f,
rho_m=rho_m,
X_chr=True,
X_chr_male=False,
N_cores=5,
density=False,
freq=False)
7.2.1 Plot density#
[19]:
fig, ax = plot_tractlength_density(bins=newbins,
curves={'H_DC': density_hp_dc},
L=which_L,
xlabel='Tract length on the X chromosome')
plt.show()
7.2.2 Plot histogram#
[20]:
fig, ax = plot_tractlength_histogram(bins=bins,
curves={'H_DC': hist_hp_dc},
xlabel='Tract length on the X chromosome',
L=which_L)
plt.show()
