Note
Go to the end to download the full example code.
MXL inference - continuous pulse model#
This example implements inference for the MXL population under a continuous pulse model of admixture, using the tracts package. Inference is performed using autosomal and X chromosome data, allowing for the specification of sex-biased admixture.
To implement this example, we use the following driver file:
samples:
directory: ./MXL_TrioPhased
individual_names: [
"NA19648","NA19649","NA19651","NA19652","NA19654","NA19655","NA19657","NA19658","NA19661","NA19663",
"NA19664","NA19669","NA19670","NA19676","NA19678","NA19679","NA19681","NA19682","NA19684","NA19716",
"NA19717","NA19719","NA19720","NA19722","NA19723","NA19725","NA19726","NA19728","NA19729","NA19731",
"NA19732","NA19734","NA19735","NA19740","NA19741","NA19746","NA19747","NA19749","NA19750","NA19752",
"NA19755","NA19756","NA19758","NA19759","NA19761","NA19762","NA19764","NA19770","NA19771","NA19773",
"NA19774","NA19776","NA19777","NA19779","NA19780","NA19782","NA19783","NA19785","NA19786","NA19788",
"NA19789","NA19792","NA19794","NA19795"]
male_names : [
"NA19649","NA19652","NA19655","NA19658","NA19661","NA19664","NA19670","NA19676","NA19679","NA19682",
"NA19717","NA19720","NA19723","NA19726","NA19729","NA19732","NA19735","NA19741","NA19747","NA19750",
"NA19756","NA19759","NA19762","NA19771","NA19774","NA19777","NA19780","NA19783","NA19786","NA19789",
"NA19792","NA19795"] #see Readme_dataprocessing.md for how this was generated
filename_format: "{name}_{label}_final.bed"
labels: [A, B] #If this field is omitted, 'A' and 'B' will be used by default
chromosomes: 1-22 #The chromosomes to use for analysis. Can be specified as a list or a range
allosomes: [X]
models:
model_filename: ../models/ccp.yaml
ad_model_autosomes: M
ad_model_allosomes: DC
start_params:
t1: 10:15
REUR: 0.07
RAFR: 0.08
RNAT: 0.095
t2: 3:5
REUR_sex_bias: -0.99 # more males
RNAT_sex_bias: 0.99 # more females
RAFR_sex_bias: -0.1
optim:
repetitions: 10
seed: 100
maximum_iterations: 100
npts: 50
exclude_tracts_below_cm: 2
unknown_labels_for_smoothing: ["UNK", "centromere","miscall"] # segments with these labels will be smoother over, that is, will be filled with neighbouring ancestries up to their midpoints.
fix_parameters_from_ancestry_proportions: ['REUR', 'RAFR','REUR_sex_bias', 'RAFR_sex_bias']
output:
output_directory: "./output_ccp/{date}/"
output_filename_format: "MXL_output_{label}"
log_filename: 'MXL_continuous_pulse.log'
verbose_log: 1
verbose_screen: 30
log_scale: True
Complete results from this analysis are saved in the output directory specified in the driver file. Below, we display the optimal parameters estimated from this analysis, as well as the plots illustrating the inferred tract length distributions, compared to the observed histograms, for every source population and chromosome type (autosomes and X chromosome).
Optimal parameters#
Optimal migration matrices#
Tract length histograms#
Autosomal admixture#
X chromosome admixture in females#
X chromosome admixture in males#
------------------------------------------------------------------------------------------------
Running tracts 2.0 with driver file: /home/runner/work/tracts/tracts/example/documentation_examples/MXL/MXL_continuous.yaml
------------------------------------------------------------------------------------------------
Results will be written to: output_ccp.
Using log file: output_ccp/MXL_continuous_pulse.log.
excluding_tracts_below set to 2.0 cM.
Re-optimization will be performed until convergence or maximum 5 times.
Ancestries: EUR, NAT, AFR
Data autosome proportions: [0.468066 0.49277278 0.03920868]
Data allosome proportions: [0.33731709 0.62309703 0.03958588]
The following parameters have been fixed from ancestry proportions: REUR, RAFR, REUR_sex_bias, RAFR_sex_bias
Model parameters and bounds:
-------------------------------------------
Parameter | Lower bound | Upper bound
-------------------------------------------
REUR | 1e-09 | 1
REUR_sex_bias | -1 | 1
RNAT | 1e-09 | 1
RNAT_sex_bias | -1 | 1
RAFR | 1e-09 | 1
RAFR_sex_bias | -1 | 1
t1 | 1 | inf
t2 | 1 | inf
-------------------------------------------
Multiple starting parameters will be generated and used for multiple optimization runs.
------------------------------------------------------------------------------
Step 1 : Optimizing autosomal likelihood over parameters ['RNAT', 't1', 't2'].
------------------------------------------------------------------------------
Starting parameters for step 1 optimization
------------------------------------------------
Run | RNAT | t1 | t2
------------------------------------------------
Warning: starting parameters are out of bounds.
1 | 0.095 | 13.95 | 4.821
2 | 0.095 | 13 | 4.07
3 | 0.095 | 11.97 | 3.644
Warning: starting parameters are out of bounds.
4 | 0.095 | 13.13 | 3.165
Warning: starting parameters are out of bounds.
5 | 0.095 | 12.73 | 3.454
6 | 0.095 | 10.48 | 4.912
7 | 0.095 | 12.9 | 3.37
8 | 0.095 | 12.33 | 3.147
9 | 0.095 | 11.78 | 4.34
10 | 0.095 | 12.88 | 4.499
------------------------------------------------
Optimization run #1
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.046. This can happen when no sex bias parameter allows for the observed ancestry proportions.
7 , -1e+23 , array([ 0.0953154 , -0.999999 , 0.095 , 0 , 8.03472e-103, -1 , 13.9513 , 4.82068 ]), Autosomes
Optimization completed.
-----------------------
Optimization run #2
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.045. This can happen when no sex bias parameter allows for the observed ancestry proportions.
10 , -1e+23 , array([ 0.095885 , -0.999999 , 0.095 , 0 , 0 , -1 , 12.9996 , 4.0705 ]), Autosomes
Optimization completed.
-----------------------
Optimization run #3
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.045. This can happen when no sex bias parameter allows for the observed ancestry proportions.
10 , -1e+23 , array([ 0.0974928 , -0.99993 , 0.095 , 0 , 0 , -1 , 11.9678 , 3.64385 ]), Autosomes
Optimization completed.
-----------------------
Optimization run #4
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.05. This can happen when no sex bias parameter allows for the observed ancestry proportions.
21 , -879.968 , array([ 0.0924973 , -1 , 0.0949899 , 0 , 0.086392 , -1 , 13.1331 , 3.19656 ]), Autosomes
Optimization completed.
-----------------------
Optimization run #5
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.047. This can happen when no sex bias parameter allows for the observed ancestry proportions.
10 , -1e+23 , array([ 0.0946854 , -0.999999 , 0.095 , 0 , 0 , -1 , 12.7328 , 3.45449 ]), Autosomes
Optimization completed.
-----------------------
Optimization run #6
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.019. This can happen when no sex bias parameter allows for the observed ancestry proportions.
4 , -722.234 , array([ 0.109281 , -1 , 0.095 , 0 , 0.0404696 , -1 , 10.4814 , 4.91161 ]), Autosomes
Optimization completed.
-----------------------
Optimization run #7
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.046. This can happen when no sex bias parameter allows for the observed ancestry proportions.
7 , -1e+23 , array([ 0.0943261 , -0.999999 , 0.095 , 0 , 0 , -1 , 12.9049 , 3.3699 ]), Autosomes
Optimization completed.
-----------------------
Optimization run #8
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.048. This can happen when no sex bias parameter allows for the observed ancestry proportions.
6 , -1e+23 , array([ 0.094375 , -1 , 0.095 , 0 , 0 , 1 , 12.3315 , 3.14739 ]), Autosomes
Optimization completed.
-----------------------
Optimization run #9
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.044. This can happen when no sex bias parameter allows for the observed ancestry proportions.
18 , -9.99999e+22, array([ 0.100065 , -1 , 0.095 , 0 , 1.35974e-15, -1 , 11.7842 , 4.34026 ]), Autosomes
Optimization completed.
-----------------------
Optimization run #10
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.045. This can happen when no sex bias parameter allows for the observed ancestry proportions.
10 , -1e+23 , array([ 0.0972069 , -1 , 0.095 , 0 , 0 , 1 , 12.8758 , 4.49881 ]), Autosomes
Optimization completed.
-----------------------
In Step 1: Results from multiple optimization runs with different starting parameters:
---------------------------------------------------------------------------------------------------------------------------------------------
Run | LogLik | REUR | REUR_sex_bias | RNAT | RNAT_sex_bias | RAFR | RAFR_sex_bias | t1 | t2
---------------------------------------------------------------------------------------------------------------------------------------------
1 | -1e+23 | 0.09532 | -1 | 0.095 | 0 | 8.035e-103 | -1 | 13.95 | 4.821
2 | -1e+23 | 0.09588 | -1 | 0.095 | 0 | 0 | -1 | 13 | 4.07
3 | -1e+23 | 0.09749 | -0.9999 | 0.095 | 0 | 0 | -1 | 11.97 | 3.644
4 | -879.968 | 0.0925 | -1 | 0.09499 | 0 | 0.08639 | -1 | 13.13 | 3.197
5 | -1e+23 | 0.09469 | -1 | 0.095 | 0 | 0 | -1 | 12.73 | 3.454
6 | -722.234 | 0.1093 | -1 | 0.095 | 0 | 0.04047 | -1 | 10.48 | 4.912
7 | -1e+23 | 0.09433 | -1 | 0.095 | 0 | 0 | -1 | 12.9 | 3.37
8 | -1e+23 | 0.09437 | -1 | 0.095 | 0 | 0 | 1 | 12.33 | 3.147
9 | -9.99999e+22 | 0.1001 | -1 | 0.095 | 0 | 1.36e-15 | -1 | 11.78 | 4.34
10 | -1e+23 | 0.09721 | -1 | 0.095 | 0 | 0 | 1 | 12.88 | 4.499
---------------------------------------------------------------------------------------------------------------------------------------------
Selecting best parameters from step 1 and proceeding to step 2 optimization.
-----------------------------------------------------------------------------
Step 2 : Optimizing allosomal likelihood over parameters : ['RNAT_sex_bias'].
-----------------------------------------------------------------------------
Starting parameters for step 2 optimization (non-sex-bias parameters are fixed to the best step 1 estimates).
------------------------------------------------------------------------------------------------------------------------------
Run | REUR | REUR_sex_bias | RNAT | RNAT_sex_bias | RAFR | RAFR_sex_bias | t1 | t2
------------------------------------------------------------------------------------------------------------------------------
1 | 0.1093 | -0.6242 | 0.095 | 0.99 | 0.04047 | 0.5738 | 10.48 | 4.912
------------------------------------------------------------------------------------------------------------------------------
Optimization run #1
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.01. This can happen when no sex bias parameter allows for the observed ancestry proportions.
10 , -189.319 , array([ 0.109281 , -0.624203 , 0.095 , 0.99 , 0.0404696 , 0.57377 , 10.4814 , 4.91161 ]), Female allosomes
10 , -124.663 , array([ 0.109281 , -0.624203 , 0.095 , 0.99 , 0.0404696 , 0.57377 , 10.4814 , 4.91161 ]), Male allosomes
Optimization completed.
-----------------------
Selecting best parameters from step 2.
Step 2 used allosomal data only. Final likelihood is evaluated on autosomal + allosomal data at the selected optimal parameters.
--------------------------------------------------------------------------------------------------
Launching re-optimization until convergence is achieved or 5 re-optimizations have been performed.
--------------------------------------------------------------------------------------------------
Re-optimization 1/5: re-optimizing starting from the current optimal parameters (likelihood = -1088.224553).
------------------------------------------------------------------------------
Step 1 : Optimizing autosomal likelihood over parameters ['RNAT', 't1', 't2'].
------------------------------------------------------------------------------
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.01. This can happen when no sex bias parameter allows for the observed ancestry proportions.
16 , -621.762 , array([ 0.084462 , -0.617181 , 0.0951435 , 0.99 , 0.0195364 , -1 , 10.4882 , 4.9254 ]), Autosomes
Optimization completed.
-----------------------
-----------------------------------------------------------------------------
Step 2 : Optimizing allosomal likelihood over parameters : ['RNAT_sex_bias'].
-----------------------------------------------------------------------------
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.003. This can happen when no sex bias parameter allows for the observed ancestry proportions.
16 , -191.67 , array([ 0.084462 , -0.619666 , 0.0951435 , 0.989801 , 0.0195364 , 1 , 10.4882 , 4.9254 ]), Female allosomes
16 , -124.571 , array([ 0.084462 , -0.619666 , 0.0951435 , 0.989801 , 0.0195364 , 1 , 10.4882 , 4.9254 ]), Male allosomes
Optimization completed.
-----------------------
Change in likelihood from -1088.224553 to -1080.535650 after re-optimizing.
Re-optimization 2/5: re-optimizing starting from the current optimal parameters (likelihood = -1080.535650).
------------------------------------------------------------------------------
Step 1 : Optimizing autosomal likelihood over parameters ['RNAT', 't1', 't2'].
------------------------------------------------------------------------------
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.039. This can happen when no sex bias parameter allows for the observed ancestry proportions.
4 , -689.335 , array([ 0.0730464 , -0.611348 , 0.0960079 , 0.989801 , 0.0200552 , 1 , 10.4882 , 4.9254 ]), Autosomes
Optimization completed.
-----------------------
-----------------------------------------------------------------------------
Step 2 : Optimizing allosomal likelihood over parameters : ['RNAT_sex_bias'].
-----------------------------------------------------------------------------
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.003. This can happen when no sex bias parameter allows for the observed ancestry proportions.
12 , -190.769 , array([ 0.0730464 , -0.627402 , 0.0960079 , 0.989978 , 0.0200552 , 0.999872 , 10.4882 , 4.9254 ]), Female allosomes
12 , -123.626 , array([ 0.0730464 , -0.627402 , 0.0960079 , 0.989978 , 0.0200552 , 0.999872 , 10.4882 , 4.9254 ]), Male allosomes
Optimization completed.
-----------------------
Change in likelihood from -1080.535650 to -1001.981667 after re-optimizing.
Re-optimization 3/5: re-optimizing starting from the current optimal parameters (likelihood = -1001.981667).
------------------------------------------------------------------------------
Step 1 : Optimizing autosomal likelihood over parameters ['RNAT', 't1', 't2'].
------------------------------------------------------------------------------
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
21 , -651.6 , array([ 0.0777449 , -0.627199 , 0.0968491 , 0.989978 , 0.0195463 , 0.575594 , 10.4975 , 4.92622 ]), Autosomes
Optimization completed.
-----------------------
-----------------------------------------------------------------------------
Step 2 : Optimizing allosomal likelihood over parameters : ['RNAT_sex_bias'].
-----------------------------------------------------------------------------
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.003. This can happen when no sex bias parameter allows for the observed ancestry proportions.
12 , -192.042 , array([ 0.0777449 , -0.626587 , 0.0968491 , 0.990074 , 0.0195463 , 1 , 10.4975 , 4.92622 ]), Female allosomes
12 , -124.147 , array([ 0.0777449 , -0.626587 , 0.0968491 , 0.990074 , 0.0195463 , 1 , 10.4975 , 4.92622 ]), Male allosomes
Optimization completed.
-----------------------
Change in likelihood from -1001.981667 to -1043.480427 after re-optimizing.
Re-optimization 4/5: re-optimizing starting from the current optimal parameters (likelihood = -1043.480427).
------------------------------------------------------------------------------
Step 1 : Optimizing autosomal likelihood over parameters ['RNAT', 't1', 't2'].
------------------------------------------------------------------------------
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.003. This can happen when no sex bias parameter allows for the observed ancestry proportions.
5 , -678.17 , array([ 0.0760479 , -0.644618 , 0.0977274 , 0.990074 , 0.0214625 , 1 , 10.603 , 4.92622 ]), Autosomes
Optimization completed.
-----------------------
-----------------------------------------------------------------------------
Step 2 : Optimizing allosomal likelihood over parameters : ['RNAT_sex_bias'].
-----------------------------------------------------------------------------
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.003. This can happen when no sex bias parameter allows for the observed ancestry proportions.
13 , -190.636 , array([ 0.0760479 , -0.645414 , 0.0977274 , 0.99007 , 0.0214625 , 1 , 10.603 , 4.92622 ]), Female allosomes
13 , -123.457 , array([ 0.0760479 , -0.645414 , 0.0977274 , 0.99007 , 0.0214625 , 1 , 10.603 , 4.92622 ]), Male allosomes
Optimization completed.
-----------------------
Change in likelihood from -1043.480427 to -992.175850 after re-optimizing.
Re-optimization 5/5: re-optimizing starting from the current optimal parameters (likelihood = -992.175850).
------------------------------------------------------------------------------
Step 1 : Optimizing autosomal likelihood over parameters ['RNAT', 't1', 't2'].
------------------------------------------------------------------------------
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.01. This can happen when no sex bias parameter allows for the observed ancestry proportions.
7 , -607.854 , array([ 0.0896673 , -0.634655 , 0.0986127 , 0.99007 , 0.0233015 , -1 , 10.7096 , 4.97573 ]), Autosomes
Optimization completed.
-----------------------
-----------------------------------------------------------------------------
Step 2 : Optimizing allosomal likelihood over parameters : ['RNAT_sex_bias'].
-----------------------------------------------------------------------------
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.039. This can happen when no sex bias parameter allows for the observed ancestry proportions.
13 , -190.665 , array([ 0.0896673 , -0.594004 , 0.0986127 , 0.99016 , 0.0233015 , 1 , 10.7096 , 4.97573 ]), Female allosomes
13 , -124.204 , array([ 0.0896673 , -0.594004 , 0.0986127 , 0.99016 , 0.0233015 , 1 , 10.7096 , 4.97573 ]), Male allosomes
Optimization completed.
-----------------------
Change in likelihood from -992.175850 to -1046.365993 after re-optimizing.
Convergence not achieved after 5 repetitions. Stopping re-optimization.
Final parameters and corresponding likelihood computed on autosomal + allosomal data:
---------------------------------------------------------------------------------------------------------------------------------------
LogLik | REUR | REUR_sex_bias | RNAT | RNAT_sex_bias | RAFR | RAFR_sex_bias | t1 | t2
---------------------------------------------------------------------------------------------------------------------------------------
-1046.37 | 0.08967 | -0.594 | 0.09861 | 0.9902 | 0.0233 | 1 | 10.71 | 4.976
---------------------------------------------------------------------------------------------------------------------------------------
The optimal solution has sex-bias parameter(s) RNAT_sex_bias near their ±1 boundary.
Re-optimizing with sex-bias parameter(s) RNAT_sex_bias fixed at their boundary value.
------------------------------------------------------------------------------
Step 1 : Optimizing autosomal likelihood over parameters ['RNAT', 't1', 't2'].
------------------------------------------------------------------------------
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
At the end of the current optimization step: Could not solve for parameters fixed by ancestry proportions. Final error: 0.003. This can happen when no sex bias parameter allows for the observed ancestry proportions.
21 , -607.392 , array([ 0.092209 , -0.663993 , 0.0995051 , 0.999 , 0.0259941 , -1 , 10.7096 , 4.97573 ]), Autosomes
Optimization completed.
-----------------------
All sex-bias parameters are fixed (REUR_sex_bias, RAFR_sex_bias by ancestry proportions; RNAT_sex_bias by user-provided values). Step 2 has no free parameters to optimize and will be skipped.
Final parameters and corresponding likelihood computed on autosomal + allosomal data:
---------------------------------------------------------------------------------------------------------------------------------------
LogLik | REUR | REUR_sex_bias | RNAT | RNAT_sex_bias | RAFR | RAFR_sex_bias | t1 | t2
---------------------------------------------------------------------------------------------------------------------------------------
-932.109 | 0.09221 | -0.664 | 0.09951 | 0.999 | 0.02599 | -1 | 10.71 | 4.976
---------------------------------------------------------------------------------------------------------------------------------------
All sex-bias parameters have been fixed by value. Boundary re-optimization completed.
Predicted autosome proportions: [0.4647373 0.48917425 0.04608846]
Predicted allosome proportions: [0.34285309 0.6289362 0.02821072]
Results saved to : output_ccp
{'destination_dir': PosixPath('/home/runner/work/tracts/tracts/docs/source/auto_examples/MXL/output_ccp'), 'table_file': PosixPath('/home/runner/work/tracts/tracts/docs/source/auto_examples/MXL/output_ccp/MXL_output_optimal_parameters.txt')}
import sys
from pathlib import Path
from tracts.driver import run_tracts
# Read files automatically for online documentation
sys.path.append('.')
script_dir = Path.cwd()
driver_filename = script_dir / "MXL_continuous.yaml"
run_tracts(
driver_filename=str(driver_filename),
script_dir=str(script_dir),
)
Total running time of the script: (4 minutes 37.176 seconds)