Adding traits to stdpopsim

[andrewkern]

the goal here is to add quantitative traits to stdpopsim!

There is a whole bunch of leg work to do. I'm putting up this issue as a place to:

• hash out the API. What do we want in terms of a user-facing interface? What parameters will be input? What options do we want to make available?
• list possible papers that we might implement from. thoughts include human GWAS (e.g. height), quantitative traits from agricultural systems (e.g. milk production in cows).
• thoughts on how to align whatever we do with future improvements to SLiM
• other stuff!

feel free to start posting ideas here

[jeffspence]

Thanks for this, @andrewkern! I think our goal should be trying to balance biological realism with matching what is already done in practice. In keeping with the stdpopsim model, it would be nice to have estimates of whatever parameters we end up using to simulate phenotypes that we can pull from the literature and QC in some meaningful way.

Some thoughts:

• I think a forward-looking way to set it up would be to have distributions of effect sizes (possibly in windows or annotations) across the genome. That would allow us to simulate neutral traits by dropping effect sizes on mutations, and also would eventually (along with some specification of a fitness function) allow us to simulate traits under selection. For multiple traits, this would be a joint distribution of effect sizes. A noncomprehensive list of relevant papers in terms of estimating what these look like (in humans, each with a different parameterization): paper 1, paper 2, paper 3.
• A sort of moonshot would be to use existing GWAS data to parameterize a model (e.g., estimate an effect size distribution -- maybe even in windows or annotations across the genome -- directly from GWAS results and then use that to simulate). That would open us up to being able to take in GWAS results, lots of which are available through the GWAS catalog. An alternative would be to use the weights from polygenic scores (maintained in the PGS catalog), which might be closer to point estimates of effect sizes than the marginal association test results reported by GWAS. In either case, we would still need some way of turning these into simulation parameters.
• The other end of the spectrum would be implementing the alpha model. It's phenomenological in the sense that it doesn't correspond to an evolutionary model, and would make it more difficult once we start thinking about selection on traits. The upside would be that there are lots of published estimates in the literature (e.g., this paper)
• Following up on a point Roshni made in the call, it is probably worth thinking about how all of this interacts with what we already have implemented. E.g., if effect sizes depend on frequency (like if we parameterize a model based on contributions to heritability instead of effect sizes) that interacts with demography (particularly multi-population models) in a not entirely straightforward way. With multiple populations we may want to eventually consider things like GxG or GxE to generate heterogeneous effect sizes (relevant papers exploring this in humans: paper 1, paper 2, paper 3). Relatedly we may want to (eventually) consider different fitness functions or different environments (relevant conceptual paper).

[roshnipatel]

Re: @jeffspence's point about the alpha model, I agree with all of the possible advantages: it's something lots of people already use, it plugs in nicely to fast coalescent simulations, and most of all, we have many many estimates across human traits. I think the trickiest thing (that I vaguely alluded to on the call) is that there's no principled way for the alpha model to intersect with the demographic models implemented in stdpopsim.

The alpha model parameterizes the relationship between effect size and frequency, which means we expect this relationship (and the corresponding estimates of alpha) to differ across populations with very different demographic histories. A couple potential options and their limitations:

• We have lots of papers that have estimated alpha from the UK Biobank and other predominantly European ancestry GWAS, which is an argument for simulating human traits such that effect sizes are conditional on frequencies in European-like populations (by which I mean populations with European-like demographies). But under most demographic models, this approach will misspecify the relationship between effect size and frequency in any non-European populations by treating P(effect | p_Eur) as equivalent to P(effect | p_Eur, p_pop2, p_pop3, ...). In the process we lose a lot of biological realism - traits will look like they're under selection in European-like populations and vaguely neutral in others.
• There's another argument for simulating traits such that effect sizes are conditional on the average frequency across populations in your demographic model -- this takes into account the joint frequency across populations, albeit in a rather crude way (though note there is some precedent, e.g. this paper). The downfall is that it's not obvious that the published estimates of alpha (from predominantly European GWAS) are a sensible way to parameterize the relationship between effect size and global allele frequency.
• A possible third option is to parameterize the relationship between effect size and the joint frequency spectrum across populations in the demographic model, i.e. fitting a parameter to describe P(effect | p_Eur, p_pop2, p_pop3, ...). This is getting closer to Jeff's moonshot suggestion, but potentially worth thinking about if we really want the ability to simulate certain traits under a coalescent framework.

Also, a couple more papers with estimates of selection coefficients from human GWAS: Zeng et al 2021 and Xue et al 2024

[gregorgorjanc]

Given that tstrait implements the "alpha" model for a post-festum trait simulation after tree sequence has already been simulated, should we be considering stdpopsim implementation at all?

Or is our goal to create a call to 1) msprime&tstrait for the backward-in-time simulation (where having ability to simulate traits according to the alpha model would be useful) and 2) SLiM for the forward-in-time (where alpha mode makes less sense, but could be used to show how using the approach 1) could skew results in some cases?!)