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: ./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]
output_filename_format: "MXL_test_output_{label}"
log_filename: 'ASW_continuous_pulse.log'
output_directory: ./output_continuous_pulse/
verbose_log: 1
verbose_screen: 30
log_scale : True
model_filename: ../models/ccc.yaml
start_params:
t1: 13.5
REUR: 0.2
RAFR: 0.02
RNAT: 0.2
t2: 6.8
REUR_sex_bias: -0.99 # more males
RNAT_sex_bias: 0.99 # more females
RAFR_sex_bias: -0.1
repetitions: 3
seed: 100
maximum_iterations: 1000
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.
exclude_tracts_below_cm: 2
npts : 50
#fix_parameters_from_ancestry_proportions: ['REUR', 'RAFR','REUR_sex_bias', 'RAFR_sex_bias']
ad_model_autosomes : M
ad_model_allosomes : DC
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#
parameter |
value |
|---|---|
REUR |
0.17171759557643532 |
REUR_sex_bias |
-0.5139335360921073 |
RNAT |
0.18230823146034544 |
RNAT_sex_bias |
-0.3299330409958683 |
RAFR |
0.02269813340105542 |
RAFR_sex_bias |
0.08710315047572315 |
t1 |
14.142002078753531 |
t2 |
6.568268268907068 |
likelihood -1130.17 |
Tract length histograms#
Autosomal admixture#
X chromosome admixture in females#
X chromosome admixture in males#
------------------------------------------------------------------------------------------------
Running tracts 2.0 with driver file: MXL_continuous.yaml
------------------------------------------------------------------------------------------------
Results will be written to: output_continuous_pulse.
Using log file: output_continuous_pulse/ASW_continuous_pulse.log.
excluding_tracts_below set to 2.0 cM.
Ancestries: EUR, NAT, AFR
Data autosome proportions: [0.468066 0.49277278 0.03920868]
Data allosome proportions: [0.33731709 0.62309703 0.03958588]
Model parameters: REUR, REUR_sex_bias, RNAT, RNAT_sex_bias, RAFR, RAFR_sex_bias, t1, t2
Admixture is modelled with the M model for autosomes and with the DC model for allosomes.
Multiple starting parameters will be generated and used for multiple optimization runs.
----------------------------------------------------------------------------------------------
Step 1 : Optimizing autosomal likelihood over parameters ['REUR', 'RNAT', 'RAFR', 't1', 't2'].
----------------------------------------------------------------------------------------------
Starting parameters for step 1 optimization
------------------------------------------------------------------------------
Run | REUR | RNAT | RAFR | t1 | t2
------------------------------------------------------------------------------
1 | 0.2 | 0.2 | 0.02 | 13.58 | 7.821
2 | 0.2 | 0.2 | 0.02 | 13.2 | 7.07
3 | 0.2 | 0.2 | 0.02 | 12.79 | 6.644
------------------------------------------------------------------------------
Optimization run #1
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
30 , -891.942 , array([ 0.201816 , 0 , 0.20138 , 0 , 0.0205409 , 0 , 13.5728 , 7.14392 ]), Autosomes
60 , -870.746 , array([ 0.197281 , 0 , 0.200909 , 0 , 0.0209357 , 0 , 13.5733 , 6.92884 ]), Autosomes
90 , -864.537 , array([ 0.195854 , 0 , 0.200363 , 0 , 0.0211367 , 0 , 13.5143 , 6.86632 ]), Autosomes
120 , -863.424 , array([ 0.195482 , 0 , 0.200276 , 0 , 0.0211756 , 0 , 13.5191 , 6.8601 ]), Autosomes
150 , -862.512 , array([ 0.195171 , 0 , 0.200206 , 0 , 0.0212078 , 0 , 13.5235 , 6.85463 ]), Autosomes
180 , -861.594 , array([ 0.194857 , 0 , 0.200122 , 0 , 0.0212402 , 0 , 13.5296 , 6.84909 ]), Autosomes
Optimization completed
----------------------
Optimization run #2
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
30 , -856.817 , array([ 0.192746 , 0 , 0.200003 , 0 , 0.021417 , 0 , 13.5088 , 6.82596 ]), Autosomes
60 , -849.169 , array([ 0.189168 , 0 , 0.19848 , 0 , 0.0216601 , 0 , 13.5648 , 6.80239 ]), Autosomes
90 , -847.076 , array([ 0.188371 , 0 , 0.197955 , 0 , 0.021765 , 0 , 13.5855 , 6.79094 ]), Autosomes
120 , -846.155 , array([ 0.188004 , 0 , 0.197657 , 0 , 0.0218046 , 0 , 13.5975 , 6.7864 ]), Autosomes
150 , -845.529 , array([ 0.187765 , 0 , 0.197481 , 0 , 0.021838 , 0 , 13.6056 , 6.78331 ]), Autosomes
180 , -844.958 , array([ 0.187529 , 0 , 0.197315 , 0 , 0.0218643 , 0 , 13.6152 , 6.78033 ]), Autosomes
Optimization completed
----------------------
Optimization run #3
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
30 , -839.57 , array([ 0.183856 , 0 , 0.195493 , 0 , 0.0221893 , 0 , 13.4096 , 6.75129 ]), Autosomes
60 , -827.174 , array([ 0.176409 , 0 , 0.187956 , 0 , 0.0227785 , 0 , 13.7219 , 6.66409 ]), Autosomes
90 , -821.985 , array([ 0.173436 , 0 , 0.184106 , 0 , 0.0228318 , 0 , 14.0421 , 6.58555 ]), Autosomes
120 , -821.129 , array([ 0.172708 , 0 , 0.183349 , 0 , 0.0227807 , 0 , 14.0742 , 6.58266 ]), Autosomes
150 , -820.108 , array([ 0.171919 , 0 , 0.182519 , 0 , 0.0227173 , 0 , 14.1276 , 6.57129 ]), Autosomes
Optimization completed
----------------------
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 | -861.593 | 0.1949 | 0 | 0.2001 | 0 | 0.02124 | 0 | 13.53 | 6.849
2 | -844.747 | 0.1874 | 0 | 0.1972 | 0 | 0.02187 | 0 | 13.62 | 6.779
3 | -819.84 | 0.1717 | 0 | 0.1823 | 0 | 0.0227 | 0 | 14.14 | 6.568
---------------------------------------------------------------------------------------------------------------------------------------------
Selecting best parameters from step 1 and proceeding to step 2 optimization.
---------------------------------------------------------------------------------------------------------------
Step 2 : Optimizing allosomal likelihood over parameters : ['REUR_sex_bias', 'RNAT_sex_bias', 'RAFR_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.1717 | 0.1545 | 0.1823 | -0.7313 | 0.0227 | 0.09647 | 14.14 | 6.568
2 | 0.1717 | -0.496 | 0.1823 | -0.3444 | 0.0227 | 0.08104 | 14.14 | 6.568
3 | 0.1717 | 0.4789 | 0.1823 | -0.1464 | 0.0227 | -0.3684 | 14.14 | 6.568
------------------------------------------------------------------------------------------------------------------------------
Optimization run #1
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
30 , 2.22045e+16 , array([ 0.171718 , 0.0717114 , 0.182308 , -0.716576 , 0.0226981 , 0.0991307 , 14.142 , 6.56827 ]), OOB (oob=-2.220446049250313e-16)
60 , -193.843 , array([ 0.171718 , 0.0408261 , 0.182308 , -0.713516 , 0.0226981 , 0.0983642 , 14.142 , 6.56827 ]), Female allosomes
60 , -138.608 , array([ 0.171718 , 0.0408261 , 0.182308 , -0.713516 , 0.0226981 , 0.0983642 , 14.142 , 6.56827 ]), Male allosomes
90 , -193.396 , array([ 0.171718 , 0.0212418 , 0.182308 , -0.710113 , 0.0226981 , 0.0980743 , 14.142 , 6.56827 ]), Female allosomes
90 , -138.414 , array([ 0.171718 , 0.0212418 , 0.182308 , -0.710113 , 0.0226981 , 0.0980743 , 14.142 , 6.56827 ]), Male allosomes
120 , -193.261 , array([ 0.171718 , 0.015284 , 0.182308 , -0.709179 , 0.0226981 , 0.0980429 , 14.142 , 6.56827 ]), Female allosomes
120 , -138.356 , array([ 0.171718 , 0.015284 , 0.182308 , -0.709179 , 0.0226981 , 0.0980429 , 14.142 , 6.56827 ]), Male allosomes
150 , -193.156 , array([ 0.171718 , 0.010592 , 0.182308 , -0.708501 , 0.0226981 , 0.0981328 , 14.142 , 6.56827 ]), Female allosomes
150 , -138.312 , array([ 0.171718 , 0.010592 , 0.182308 , -0.708501 , 0.0226981 , 0.0981328 , 14.142 , 6.56827 ]), Male allosomes
180 , -193.056 , array([ 0.171718 , 0.00612233 , 0.182308 , -0.707805 , 0.0226981 , 0.0980413 , 14.142 , 6.56827 ]), Female allosomes
180 , -138.269 , array([ 0.171718 , 0.00612233 , 0.182308 , -0.707805 , 0.0226981 , 0.0980413 , 14.142 , 6.56827 ]), Male allosomes
Optimization completed.
-----------------------
Optimization run #2
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
30 , -182.116 , array([ 0.171718 , -0.508668 , 0.182308 , -0.332799 , 0.0226981 , 0.0874006 , 14.142 , 6.56827 ]), Female allosomes
30 , -131.399 , array([ 0.171718 , -0.508668 , 0.182308 , -0.332799 , 0.0226981 , 0.0874006 , 14.142 , 6.56827 ]), Male allosomes
60 , 2.22045e+16 , array([ 0.171718 , -0.513851 , 0.182308 , -0.329976 , 0.0226981 , 0.087103 , 14.142 , 6.56827 ]), OOB (oob=-2.220446049250313e-16)
Optimization completed.
-----------------------
Optimization run #3
Iter. Log-likelihood Model parameters Transmission
-------------------------------------------------------------
30 , -202.118 , array([ 0.171718 , 0.443036 , 0.182308 , -0.132377 , 0.0226981 , -0.358753 , 14.142 , 6.56827 ]), Female allosomes
30 , -137.313 , array([ 0.171718 , 0.443036 , 0.182308 , -0.132377 , 0.0226981 , -0.358753 , 14.142 , 6.56827 ]), Male allosomes
60 , -201.99 , array([ 0.171718 , 0.437971 , 0.182308 , -0.131184 , 0.0226981 , -0.357928 , 14.142 , 6.56827 ]), Female allosomes
60 , -137.256 , array([ 0.171718 , 0.437971 , 0.182308 , -0.131184 , 0.0226981 , -0.357928 , 14.142 , 6.56827 ]), Male allosomes
Optimization completed.
-----------------------
Step 2: 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 | -331.272 | 0.1717 | 0.004472 | 0.1823 | -0.7076 | 0.0227 | 0.09811 | 14.14 | 6.568
2 | -313.382 | 0.1717 | -0.5139 | 0.1823 | -0.3299 | 0.0227 | 0.0871 | 14.14 | 6.568
3 | -339.237 | 0.1717 | 0.4377 | 0.1823 | -0.1312 | 0.0227 | -0.3579 | 14.14 | 6.568
---------------------------------------------------------------------------------------------------------------------------------------------
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.
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
---------------------------------------------------------------------------------------------------------------------------------------
-1130.17 | 0.1717 | -0.5139 | 0.1823 | -0.3299 | 0.0227 | 0.0871 | 14.14 | 6.568
---------------------------------------------------------------------------------------------------------------------------------------
Results saved to : output_continuous_pulse
{'destination_dir': PosixPath('/home/runner/work/tracts/tracts/docs/source/auto_examples/MXL/output_continuous_pulse'), 'table_file': PosixPath('/home/runner/work/tracts/tracts/docs/source/auto_examples/MXL/output_continuous_pulse/MXL_test_output_optimal_parameters.txt')}
import sys
from pathlib import Path
sys.path.append('.')
from tracts.driver import run_tracts
script_path = Path(sys.argv[0]).resolve()
driver_filename = "MXL_continuous.yaml"
run_tracts(driver_filename = driver_filename, script_dir = script_path.parent)
# Don't run the code below: for documentation purposes only.
from tracts.doc_utils import prepare_example_outputs_for_docs
prepare_example_outputs_for_docs("output_continuous_pulse")
Total running time of the script: (168 minutes 11.878 seconds)