popsim-consortium / stdpopsim

A library of standard population genetic models
GNU General Public License v3.0
122 stars 86 forks source link

Handling of sex chromosomes #383

Open jeromekelleher opened 4 years ago

jeromekelleher commented 4 years ago

In the meeting 2020-01-14 we discuss the current (mis)handling of sex chromosomes. We're not handling these properly at the moment in the shipped software, so we need to do something about it. Options discussed:

  1. Remove the X chromosome from the species definition so people can't use it
  2. Print out a big warning when people use the X chromosome in stdpopsim
  3. Update the code so that we do a reasonable job of simulating the X by changing Ne according to the expected sex ratio in the species

We decided on 2 initially, followed by 3.

petrelharp commented 4 years ago

Proposal was, for coalescent simulations:

(Since: 1/(2Ne) is the chance of coalescence; for two X chromosomes to coalesce they must pick the same sex parent (probability 1/4); then, if female, they must pick the same female (prob 1/Nf) and then both pick the same X chromosome (1/2); if male, they must only pick the same male (prob 1/Nm); putting this together,

(prob of coal) = (1/4) * (1/(2Nf) + 1/Nm) = (1/8) * (2 Nf + Nm) / (Nf Nm)

and putting in Nm = p N and Nf = (1-p)N we get that

(prob of coal) = (1/4) * ( (2-p) / (p (1-p)) ) * (1/(2N))

so we need to set Ne so this is equal to 1/(2Ne); which is obtained if Ne = (4 p (1-p) /(2 - p)) * N.)

However: if we're really adding this new parameter, p, we should rescale the usual coalescence rate as well! I kinda think we should just stick with p=1/2.

reedacartwright commented 4 years ago

I would add two more things to implement when needed.

reedacartwright commented 4 years ago

This is what I got for p = 1/2.

X-chromosome coalescent rate is 2 / (3 Ne).

(2 / 3 Ne) / (1 / (2 Ne) is 4/3. So Ne_X = 3/4 Ne_A.

Explanation

Looking into the past, you have three types of X-chromosomes.

The frequencies of these are p/(2-p), (1-p)/(2-p), and (1-p)/(2-p). Therefore,

Putting this together to total rate of coalescence is

[(1-p)/(2-p)]^2 (1/Nm) + [1/(2-p)]^2 (1 / (2 Nf))

The inverse of this is 2N*p(1-p)(2-p) / (1- 2 p(1-p)) and thus the correction factor for Ne is

p(1-p)(2-p) / (1-2 p(1-p))

petrelharp commented 4 years ago

Whoops - good point - I forgot to reweight by sex.

petrelharp commented 1 year ago

Adding notes from @mathbionerd's stdpopsim talk just now: Screenshot from 2023-03-08 12-35-36 Screenshot from 2023-03-08 12-34-43 Screenshot from 2023-03-08 12-36-50