code for estimating split time and ancestral population size in two-island split model

Introduction

This repository contains code used in Huang et al.. It estimates genetic split time and the ancestral population size in a simple two-island split model by fitting the empirical IBD decay curve. Note that this simplistic model assumes no gene flow after the split, an assumption often violated in real populations. Gene flow after the initial split will cause underestimated split time.

Usage

The main function is "two_island_noGeneFlow_constNe_MLE_multiStart", located in ./code/analytic.py. It needs the following required arguments:

• ibd_segments: a list of IBD segment length shared between pop1 and pop2, e.g, [4.05, 4.678, 5.65]. Length should be in cM.

• minLen: minimum length cutoff (in cM) for IBD used in inference.

• maxlen: maximum length cutoff (in cM) for IBD used in inference.

• step: binning width (in cM).

• G: a list of chromosome length (in cM). For human genetics, one can use [286.279, 268.840, 223.361, 214.688, 204.089, 192.040, 187.221, 168.003, 166.359, 181.144, 158.219, 174.679, 125.706, 120.203, 141.860, 134.038, 128.491, 117.709, 107.734, 108.267, 62.786, 74.110].

• numPairs: number of pairs of haplotype pairs (that is, 4number of diploid sample pairs, or equivalently, 4(# of samples in pop1)*(# of samples in pop2)).