nspope
commented
6 months ago I think the big API question is whether we want singletons to be present in the input tree sequence or not. If the former, they could be rephased (e.g. the output might swap node labels for some mutations). If the latter, the question is how to pass them into tsdate and how to return them (e.g. should they just be used for dating? Should they be added to the mutation table?)
In general we'll have to assume singletons from a particular individual are all phased or all unphased.
hyanwong
jeromekelleher
petrelharp
duncanMR
It'll frequently be the case that singletons are not phased. Singleton variants aren't used by tsinfer, but are crucial for correct dating. Here's how we'll handle unphased singletons:
Split edges leading to the same individual into "phase blocks" -- e.g. the intersection of intervals for these edges. For example, over some region of the genome the pair of edges leading to the same individual might have coordinates $[l_i, r_i]$ and $[l_j, r_j]$. The phase block for these two edges would have coordinates $[\max(l_i, l_j), \min(r_i, r_j)]$.
Use an expectation propagation update by phase block rather than separately for the two constituent edges. For example, say we have a phase block where the parents of the constituent edges have ages $t_i$ and $tj$. The children of the edges are both present-day haplotypes (of age 0). If the phase block has (mutation-weighted) span $\mu$ and contains $y$ singleton mutations, then the total mutational target is $\lambda{ij} = \mu (t_i + tj)$ and the mutation likelihood is $y \sim \mathrm{Poisson}(\lambda{ij})$.
So, if we have variational posteriors $t_i \sim \mathrm{Gamma}(\alpha_i, \beta_i), t_j \sim \mathrm{Gamma}(\alpha_j, \beta_j)$, the partition function used in the expectation propagation update for the phase block is, $$\int_0^{\infty} \int_0^{\infty} t_i^{\alpha_i - 1} e^{-\beta_i t_i} t_j^{\alpha_j - 1} e^{-\beta_j t_j} (t_i + t_j)^{y} e^{-\mu (t_i + t_j)} dt_i dt_j = \mathcal{B}(\alpha_i, \alpha_j) \frac{\Gamma(\alpha_i + \alpha_j + y)}{(\beta_i + \mu)^{\alpha_i + \alpha_j + y}} {_2} F_1 (\alpha_j, \alpha_i + \alpha_j + y; \alpha_i + \alpha_j; 1 - \frac{\beta_j + \mu}{\beta_i + \mu})$$ and we get the updated moments of $t_i, t_j$ by differentiating this wrt $\alpha, \beta$ as usual.
This gets substantially more complicated if the individual in question is not in the present day (especially if the age of the individual is unknown). So for now we'll have to assume that mutations on ancient samples are phased (I think these are often pseudohaploid sequences anyways).