The minimal observable clade size of exchangeable coalescents
Fabian Freund Arno Siri-Jégousse


Observable clade sizes]The minimal observable clade size of exchangeable coalescents


 

  Crop Plant Biodiversity and Breeding Informatics Group (350b), Institute of Plant Breeding, Seed Science and Population Genetics,   University of Hohenheim, Fruwirthstrasse 21, 70599 Stuttgart, Germany
 fabian.freund@uni-hohenheim.de
 
  
  Departamento de Probabilidad y Estadística,
IIMAS, Universidad Nacional Autónoma de México, Mexico City, Mexico. 
 arno@sigma.iimas.unam.mx


For Λ-n-coalescents with mutation, we analyse the size O_n of the partition block of i∈{1,…,n} at the time where the first mutation appears on the tree that affects i and is shared with any other j∈{1,…,n}. We provide asymptotics of O_n for n→∞ and a recursion for all moments of O_n for finite n.
This variable gives an upper bound for the minimal clade size <cit.>, which is not observable in real data.
In applications to genetics, it has been shown to be useful to lower classification errors in genealogical model selection <cit.>. 



[2010]Primary 60C05; Secondary 92D20, 60F15, 60G09
[
    Arno Siri-Jégousse
    December 30, 2023
======================





§ INTRODUCTION

The potential for adaption of organisms to diverse environments is based on their genetic diversity. Moreover, the specific historic pattern of adaptation and demography leaves distinct marks in the genetic diversity of a sample taken from a population of said organisms. When observing the genetic diversity of a single non-recombining part of the genome, the diversity can be described by the inheritance pattern of the mutations on the genealogical tree of the sample, usually given by a Poisson point process on the genealogy. Modelling the genealogy is thus an important aspect of modelling genetic diversity. Usually, the exact genealogy is not known and cannot be reconstructed perfectly from observed genetic data (an example is provided later). Thus, genealogy models are usually defined as random variables on the set of possible genealogical trees.
  
Here, we are concerned with the genealogical tree of n alleles, i.e. the genetic information of a sample of size n from a genomic region.  Coalescent theory provides a rich class of genealogical tree models for a sample of n alleles as much as elegant tools and a convenient setting up for statistical inference. In particular, Kingman's coalescent <cit.> and the larger family of coalescents with multiple collisions <cit.> were widely studied in the past decades.
This class of Markov processes on the set of partitions of [n]:={1,…,n} is characterized by a finite measure Λ on [0,1], justifying their name of Λ-n-coalescents.
If a Λ-n-coalescent has b blocks, any given k-tuple of them will merge at rate

    λ_b,k=∫_0^1x^k-2(1-x)^b-kΛ(dx)

and the rate for the next coalescence event is

    λ_b=∑_k=2^bbkλ_b,k.

The starting partition is {1},…,{n}. The genealogical tree is recovered from the partition-valued process by first starting n branches from the leaves 1,…,n. Then, any merger of partition blocks corresponds to a joining of branches in a node (ancestor), where a single new branch starts. Partition blocks correspond to branches in the genealogical tree, the time a partition block is not merged  gives the length of this branch.  We refer to <cit.> for a survey.

Genealogical trees of the alleles in a genetic region come with an interpretation of relatedness and genetic similarity: An allele i is more closely related to allele j than to allele l if the common ancestor of (i,j) appears more recently than the common ancestor of (i,l), while the path lengths between leaves (alleles) measure the time available to accumulate mutations that decrease genetic similarity. Several statistics aim to capture these aspects and their biological meaning. For instance, the minimal clade size of an allele i gives the number of closest relatives of i, see <cit.>. Another example is the length of the external branch of i, i.e. the waiting time for the first merger of i, which gives a measure of the genetic uniqueness of an individual <cit.>. 
The mathematical properties of the minimal clade size <cit.>, the length of an external branch <cit.>, as much as the family (partition block) sizes at this time <cit.>, have been analysed recently.
In these works, asymptotic and exact behaviours are obtained for various examples of exchangeable coalescents.

However, these statistics cannot be observed directly from the genetic data. We will illustrate this for the minimal clade. By a clade we denote the set of all alleles that share a specific ancestor, and the minimal clade of i is the smallest clade including i. Assume the infinite-sites model of mutation, each mutation causes a change at a different position in the genomic region. Further assume that we know the ancestral state at the genomic region, i.e. we can identify mutations as changes compared to the ancestral state. Any clade can only be observed if there is at least one mutation that all its members share. This mutation is inherited from the common ancestor, thus has to be placed on the branch that connects this ancestor further towards the root of the genealogy (the most recent common ancestor of the whole sample). Thus, we can only observe a clade if there is a mutation on the branch directly above of the ancestor defining this clade.

Instead of looking at the minimal clade of an allele i, one could consider the smallest clade which includes i that can be observed from the data. We considered the sizes of these clades for all alleles sampled, the minimal observable clade sizes, in <cit.>. There we could show that they provide an additional set of statistics that faciliates the inference of a well-fitting genealogy model when coupled with standard statistics of genetic diversity as the site frequency spectrum.

In this article, we study some mathematical properties of the minimal observable clade size for an individual i. Its asymptotic behaviour for any Λ-n-coalescent for n→∞ as well as a recursion for all moments for finite n are established. For the Bolthausen-Sznitman coalescent, which provides a somewhat universal genealogical model for populations under strong selection, see e.g. <cit.>, <cit.>, <cit.>, we can show that the minimal observable clade size is asymptotically Beta-distributed.




§ A FORMAL DEFINITION OF THE MINIMAL OBSERVABLE CLADE SIZE

Let [n]:={1,…,n} for n∈ℕ, [n]_0=[n]∪{0}. For any Λ-n-coalescent and a sampled allele i∈[n], define 


  * 𝒞_n,i(t) as the partition block i is in at time t (a size-biased pick of a block of the n-coalescent at time t)

  * κ_n as the total number of jumps of the Λ-n-coalescent, κ_n,i as the total number of jumps (the block of) i participates in

  * K_n,i(0)(=0),K_n,i(1),…,K_n,i(κ_n,i)(=κ_n)∈[κ_n]_0 as the successive indices of jumps in the Λ-n-coalescent in which the block of i is involved   

  * 𝒞_n,i[k] as the partition block i is in at the time of its kth jump K_n,i(k), k∈[κ_n,i]_0. 𝒞_n,i[0]={i}, 𝒞_n,i[1] is the minimal clade of i, 𝒞_n,i[κ_n,i]=[n]. 

Given the Λ-n-coalescent tree, we set mutations on its branches via a homogeneous Poisson point process with rate . 
Mutations are interpreted under the infinite sites model, each mutation hits a site not hit by any other mutation, producing a new type. The new type is called derived type in contrast to the ancestral type of the most recent ancestor of the sample. Mutations on external branches are affecting only one individual, we will call these private mutations; they can also be referred to as singleton mutations. All other mutations are called non-private mutations. Since we are interested in the mutations carried by individual i, we have to record the mutations from t=0 to the time back to the most recent common ancestor of the sample (the root of the genealogy) on the path of i.  Let T^(i)_n be the waiting time until the first (youngest) mutation on the path of i that is non-private, i.e. does not fall on the external branch which ends in i (which has length E^(i)_n). If we continue the path of i after reaching the most recent common ancestor as a single ancestral line indefinitely (which we will do from now on), we have 

    T^(i)_nd=E^(i)_n+M^(i)

for an independent exponential random variable M^(i) with rate . Let 

    L^(i)_n:=max{k,}∈[κ_n,i],

i.e. the L^(i)_nth jump that i participates in is the last jump of it before T^(i)_n. The minimal observable clade of i is then given by 

    𝒞_n,i[L^(i)_n] ={j∈[n]:}

The definitions are equivalent since all non-private mutations of i are inherited from the youngest ancestor of i that bears at least one mutation on the branch connecting it to the next older ancestor. If i has no non-private mutations, 𝒞_n,i[L^(i)_n]=[n] almost surely, since in this case T^(i)_n is larger than the time back to the most recent common ancestor.

The statistic we are interested in is the size of the minimal observable clade of an allele i 

    O_n(i):=|𝒞_n,i[L^(i)_n]|.

See Figure <ref> for an example. Due to exchangeability, the distribution of O_n(i) does not depend on i, we can even choose i randomly without changing the distribution. For ease of notation, we fix the allele we are interested in to allele 1 and abbreviate O_n:=O_n(1).

Since the partition block including i can only increase in size over time, the minimal clade 𝒞_n,i[1] is a subset of the minimal observable clade 𝒞_n,i[L^(i)_n] for i∈[n]. Thus, the size M_n(i) of the minimal clade of i∈[n] satisfies M_n(i)≤ O_n(i). See Figure <ref> for an example.





§ ASYMPTOTICS OF THE OBSERVABLE CLADE SIZE

Asymptotically for sample size n→∞, the probabilistic structure of O_n simplifies considerably. 
First, we focus on coalescents without dust, a class which includes Beta(2-α,α)-n-coalescents for α∈(1,2) <cit.>, Kingman's n-coalescent (Λ=δ_0) and the Bolthausen-Sznitman n-coalescent (Λ uniform on [0,1]). A Λ-n-coalescent has dust if and only if μ_-1:=∫_[0,1] x^-1Λ(dx)<∞, see <cit.>.
 


Let O_n be defined for Λ-n-coalescents such that μ_-1=∞ (without dust) and with mutation rate . We have 
    O_n/na.s.→ S
 for n→∞ with S>0 a.s.. S is distributed as the size of the block of individual 1 (alternatively a size-biased pick of a block size) at a random time Md=Exp(). The distribution of S is uniquely determined by its moments

    E(S^k)=1-∑_r=2^k+1 a_k+1,r/λ_r+,
  
where λ_r is the total rate of the Λ-coalescent in a state with r blocks and a_k+1,r is a rational function of λ_2,…,λ_k+1, defined as in <cit.>. In particular,

    E(S)=Λ([0,1])/Λ([0,1])+, E(S^2)=1-3/2/Λ([0,1])++1/2/λ_3+.
 
  
 

Let E^(1)_n be the waiting time for the first collision of individual 1. By the consistency of the n-coalescents, we have
 E^(1)_2≥ E^(1)_3≥…. By a slight adaptation of <cit.>, we see that E^(1)_nd→0 for n→∞. Since (E^(1)_n)_n∈ℕ is monotonically decreasing, this convergence also holds almost surely.


All mutations of individual 1 on any n-coalescent lie on the path of leaf 1 to the root of the coalescent tree and are consistent for different values of n (any mutation on the path to the root in the m-coalescent is also a mutation on this path in every n-coalescent with n>m), since for m<n, the m-coalescent (seen as a tree) is the subtree of the n-coalescent which is spanned by the leaves [m], including mutations. Thus, we can represent the mutations of individual 1 on the n-coalescent by one common homogeneous Poisson process for all n, independent of the n-coalescents, on [0,∞) with rate  
This gives a new representation for T_n^(1), it is the smallest Poisson point T with T≥ E^(1)_n. Let M be the smallest Poisson point overall. Since E^(1)_n→ 0 a.s. for n→∞, for any realisation of the coalescent there exists a n_0∈ℕ s.t. E^(1)_m<M a.s. for all n≥ n_0. This shows that T_n^(1)=M a.s. for m≥ n_0, which implies 

    lim_n→∞O_n/n=lim_n→∞𝒞_n,1(M)/n=f_1(M)(=f_1(M^(1))) ,

where f_1(t) is the (asymptotic) frequency of the block individual 1 is in at time t≥ 0 in the Λ-coalescent (with values in ℕ, see <cit.>). The existence of the limit follows from Kingman's correspondence, since the coalescent stopped at the random time M (independent of the coalescent) gives an exchangeable partition of ℕ. Since we have a coalescent without dust, we have no singleton blocks a.s. at any time t>0 and a (potentially infinite) number of blocks with a.s. positive frequencies, again due to Kingman's correspondence. This shows f_1(M)>0 a.s..
 Due to exchangeability, the distribution of f_1(M) is the same as if we would make a size-biased pick from all blocks present.

Denoting f_1(M) by S, consider the moments 

    E(S^k)=∫_0^∞ E((f_1(t))^k)  e^- t dt.

From <cit.> we see that E((f_1(t))^k)=1-∑_r=2^k+1 a_k+1,r e^-λ_r t due to a connection to the exchangeable partition function of the coalescent at time t. Thus, we have

    E(S^k)=1- ∑_r=2^k+1 a_k+1,r∫_0^∞ e^-λ_r t e^- t dt=1-∑_r=2^k+1 a_k+1,r/λ_r+,

which is Eq. (<ref>). Using the explicit values a_2,2=1,a_3,2=3/2,a_3,3=-1/2 (essentially from <cit.>) and λ_2=Λ([0,1]) yields the first two moments. Since S takes values in [0,1], its distribution is uniquely determined by its moments.    




For the special case of the Bolthausen-Sznitman coalescent, the law of S can be identified

For the Bolthausen-Sznitman n-coalescent, 

    O_n/na.s.→ S,

for n→∞ where Sd=Beta(1/1+,/1+).

 

<cit.> shows that for the Bolthausen-Sznitman coalescent, (f_1(t))_t≥ 0 jumps at independent standard exponential times with ranked jump sizes given by a Poisson-Dirichlet distribution with parameters (0,1). The set of jump times is independent of the set of jump sizes. Comparing this with Eq. (<ref>), we see that to compute Sa.s.=f_1(M^(1)), we need to sum the sizes of all jumps of f_1 that happen before or at M^(1). Consider the jump times (T_k)_k∈ℕ of f_1 ordered according to the rank of their jump sizes (J_k)_k∈ℕ. Define B_k:=1_{T_k≤ M^(1)}. Hence, we have
P(B_k=1)=(1+)^-1. We can now express

    Sa.s.=∑_k∈ℕ B_kJ_k.

In other words, S can be seen as summing up a random thinning of a standard Poisson-Dirichlet distributed random variable.

The random variable S is Beta distributed. To see this, we will use the construction of the PD(0,θ') distribution from <cit.>, which is also summarised in <cit.>. Consider the points 𝒫:=(P_k)_k∈ℕ of a Poisson point process on [0,∞) with mean measure ν_θ':=θ' e^-x/xdx. Then, the size-ordered and normalized points  (P_[k]/P)_k∈ℕ with P=∑_i∈ℕP_k have the Poisson-Dirichlet distribution and are independent of 

    Pd=Gamma(θ',1),

where Gamma(α,β) is the Gamma distribution with shape parameter α and rate β.

We choose θ'=1 and make the correspondence between the ranked and normalised points (P_k)_k∈ℕ and the jump sizes (J_k)_k∈ℕ. To express Eq. (<ref>), we give each point P_k a mark B_k∈{0,1}. Marks are independent from (P_k)_k∈ℕ and from one another. We set the probability to be marked to m:=P(B_k=1)=(1+)^-1 for all k∈ℕ. (P_k,B_k)_k∈ℕ is a marked Poisson process. The Colouring Theorem <cit.> shows that all points P_k with marks 1 form a Poisson point process 𝒫_1 with mean measure mν_1, while all points with mark 0 form a Poisson point process 𝒫_0 with (1-m)ν_1.

We can now alternatively express (<ref>) as

    Sd=∑_p∈𝒫_1p/∑_p∈𝒫_1p+∑_p∈𝒫_0p=:X/X+Y.
         
where X and Y are independent due to the independence of 𝒫_0 and 𝒫_1. Since the mean measures of 𝒫_0 and 𝒫_1 are of the form θ' e^-x/xdx with θ' equal to m and 1-m, Eq. (<ref>) yields Xd=Gamma(m,1) and Yd=Gamma(1-m,1). Thus, X/X+Y is Beta-distributed with parameters m and 1-m.



Theorem <ref> can be generalised for some time-changed Λ-n-coalescents without dust, which appear when modelling genealogies in Cannings models with moderate fluctuations in population size, see <cit.>, <cit.>, <cit.> and <cit.>. 
Let (Π^(n)_g(t))_t≥ 0 be a time-changed Λ-n-coalescent, where g(t):=∫_0^tμ(s)ds  with continuous μ:[0,+∞)↦ [0,+∞), which includes some time changes proposed for Λ-n-coalescents in the references above. Observe that g is continuous, monotone and invertible with differentiable inverse. The time-changed Λ-n-coalescent is still exchangeable. The almost sure convergence of n^-1O_n for the time-changed Λ-n-coalescent works analogously as in Theorem <ref>. The time of the first merger of individual 1 is g(E^(1)_n), thus also converges to 0 almost surely. The waiting time M'_1 for the first mutation on the path of 1 to the root is an Exp()-distributed random variable, but on the time-changed path of (Π^(n)_g(t))_t≥ 0. Thus, the limit of n^-1O_n for the time-changed Λ-n-coalescent is the frequency of the block containing 1 at time M'_1 in (Π^(n)_g(t))_t≥ 0. This can also be expressed as f_1(g(M'_1)), where f_1 is said frequency in the Λ-n-coalescent (Π^(n)_t)_t≥ 0. The distribution of g(M'_1) is given by

    P(g(M'_1)≤ t)=P(M'_1≤ g^-1(t))=1-e^- g^-1(t),

which has density t↦μ(g^-1(t))^-1 e^- g^-1(t). Analogously to (<ref>) we can thus express, in terms of the a_k,r from Theorem <ref>, the kth moment of n^-1O_n for the n-coalescent with exponential growth as  

    E(S^k)=    1- ∑_r=2^k+1 a_k+1,r∫_0^∞ e^-λ_r tμ(g^-1(t))^-1 e^- g^-1(t)dt
    
    =    1- ∑_r=2^k+1 a_k+1,r∫_0^∞μ(g^-1(t))^-1 e^- g^-1(t)-λ_r tdt.

As an example, consider Kingman's n-coalescent with exponential growth with rate ρ=0. From <cit.>, we see that μ(t)=e^ρ t and thus g^-1(t)=ρ^-1log(1+ρ t). This leads to moments 

    E(S^k)=1-∑_r=2^k+1 a_k+1,r∫_0^∞ (1+ρ t)^-θ/2ρ-1e^-r2t.

 




Now, consider coalescents (Π^(n)_t)_t≥ 0 with dust which stay infinite. An example for this are Dirac n-coalescents with Λ=δ_p, p∈(0,1) <cit.>. For Λ-coalescents with dust, staying infinite is equivalent to Λ({1})=0.   

Let O_n be defined for Λ-n-coalescents with  μ_-1<∞ (with dust), P(lim sup_n→∞|Π^(n)_t|=∞ ∀ t>0)=1 and with mutation rate . We have 
    O_n/na.s.→ S
 for n→∞ with S>0 a.s.. We have
E(S)=1-θ/2μ_-1a/1-a with a=(1-Λ([0,1])/μ_-1)(/+μ_-1). 
    

From <cit.>, we see that the asymptotic frequency of the block of 1 forms an increasing jump-hold process f_1 with f_1a.s.=lim_n→∞𝒞_n,1(t)/n with values in [0,1], positive jumps and with i.i.d. Exp(μ_-1) waiting times between jumps. It fulfills E(f_1[k])=1-(1-Λ([0,1])/μ_-1)^k, where f_1[k] is the value of f_1 at its kth jump. We record between which indices of jumps K, K+1∈ℕ of f_1 the waiting time T^(1)_n for the first non-private mutation falls. From the proof of <cit.>, we know that there exists n_0∈ℕ so that E^(1)_n equals the time of the first jump of f_1 for all n≥ n_0 almost surely. Similarly to the proof of Theorem <ref>, we just have to trace back the first mutation after this first jump whose time of appearance does not depend on n. This implies that T^(1)_n falls between the same K,K+1 for all n≥ n_0. Thus we have lim_n→∞ n^-1O_n=f_1[K] a.s., where f_1[K] is the state of f_1 at the Kth jump. We only need to find the distribution of K. For n≥ n_0, T^(1)_n is the waiting time for the first jump of f_1 plus an independent Exp() random variable M^(1). Using that the waiting times (T_1,k)_k∈ℕ between the jumps of f_1 are i.i.d., we have K=1+Y, where Y is defined by ∑_i=1^Y T_1,i+1≤ M < ∑_i=1^Y+1 T_1,i+1 and thus Yd=Geo(/μ_-1+) on ℕ_0. This yields Kd=Geo(/μ_-1+) on ℕ. We compute

    E(S)   =∑_k∈ℕE(f_1[k])P(K=k)
       =1-∑_k∈ℕ(1-Λ([0,1])/μ_-1)^k(μ_-1/μ_-1+)^k-1/μ_-1+
       = 1- θ/2μ_-1∑_k∈ℕ((1-Λ([0,1])/μ_-1)(/+μ_-1)_=:a)^k
       = 1-θ/2μ_-1a/1-a
      



§ RECURSIONS FOR THE MOMENTS

To obtain recursive formulae for the moments of O_n for a Λ-n-coalescent, we first need to introduce
X_n, the size of the block of 1 at the exponential clock M of rate θ/2 in the n-coalescent. 






















Let j≥1.
The jth moments of X_n and O_n satisfy the following recursions: E(X_1^j)=1, E(O_2^j)=2^j and

    E(X_n^j)=/+λ_n+1/+λ_n∑_k=2^nnkλ_n,kE(X_n-k+1^j+∑_i=1^ja_i-1,j,kX_n-k+1^i/n-k+1).

and

    E(O_n^j)
    =∑_k=2^nn-1k-1λ_n,k/λ_n∑_i=0^j(a_i,j.kE(X_n-k+1^i)
    +(n-k+1) 1_i=j+b_i,j,k/n-kE(O_n-k+1^i)),

where a_-1,j,k=0,
a_i,j,k=ji(k-1)^j-i,
and
b_i,j,k=a_i-1,j,k-a_i,j,k.



Our proofs rely on tracking the number of blocks involved in the first jump of the Λ-n-coalescent, with some additional condition(s). 
Let us first prove (<ref>).
Let T_n be the waiting time for the first coalescence in the n-coalescent which is Exp(λ_n)-distributed. 
Let A_n  be the event that the first block merged is part of the block of 1 stopped at M.

    E(X_n^j)
    =   P(M^(1)≤ T^(1)_n)+P(M^(1) > T^(1)_n)∑_k=2^nnkλ_n,k/λ_n(E(( k-1+X_n-k+1)^j1_A_n)+E(X_n-k+1^j(1-1_A_n)))
    
    =   /+λ_n+1/+λ_n∑_k=2^nnkλ_n,k(E((k-1+X_n-k+1)^jX_n-k+1/n-k+1)+E(X_n-k+1^j(1-X_n-k+1/n-k+1)))
    
    =   /+λ_n+1/+λ_n∑_k=2^nnkλ_n,kE(X_n-k+1^j+∑_i=1^jji-1(k-1)^j-i+1X_n-k+1^i/n-k+1)


Now let us turn to the proof of (<ref>).
Let B_n={K_n,1(1)>1} be the event that 1 does participate in the first coalescence event.
Also let C_n  be the event that the first block merged is part of the observed clade.

    E(O_n^j)   =E(O_n^j1_B_n)+E(O_n^j(1-1_B_n))
    
       =∑_k=2^nn-1k-1λ_n,k/λ_nE((k-1+X_n-k+1)^j)
    
       +∑_k=2^nn-1kλ_n,k/λ_n
    E(O_n-k+1^j(1-1_C_n))
    
       +∑_k=2^nn-1kλ_n,k/λ_n
    E(( k-1+O_n-k+1)^j1_C_n)
       =∑_k=2^nn-1k-1λ_n,k/λ_nE((k-1+X_n-k+1)^j)
    
       +∑_k=2^nn-1kλ_n,k/λ_n
    E(O_n-k+1^j(1-O_n-k+1-1/n-k))
    
       +∑_k=2^nn-1kλ_n,k/λ_n
    E((k-1+O_n-k+1)^jO_n-k+1-1/n-k).

Expanding, we obtain the result.


For Kingman's n-coalescent, (<ref>) and (<ref>) considerably simplify.
In particular the two first moments of X_n are

    E(X_n)=/+n2+n2/+n2n/n-1E(X_n-1)

and

    E(X_n^2)=/+n2+n2/+n2(n+1/n-1E(X_n-1^2)+1/n-1E(X_n-1)).

and the two first moments of O_n are

    E(O_n)=2/n(1+E(X_n-1))-1/n+n-1/nE(O_n-1)

and

    E(O_n^2)=2/n(1+2E(X_n-1)+E(X_n-1^2))-1/n-1/nE(O_n-1)+E(O_n-1^2)





FF was funded by DFG grant FR 3633/2-1
through Priority Program 1590: Probabilistic Structures in Evolution. 



plain


