Background

Several papers have developed algorithms for calculating identity coefficients, given a pedigree (Karigl, 1981; Abney, 2009;  Garcia-Cortes, 2015; Kirkpatrick et al., 2018). These assume a single genetic locus, and primarily consider the 9 condensed identity coefficients that describe the relationship between two diploid individuals; however, Karigl (1982) considers identities amongst more than two individuals.  This work has developed algorithms that can efficiently calculate identity coefficients involving two individuals, across large pedigrees .

Here, we define and implement a (fairly) simple algorithm that deals with multiple sets of genes across multiple individuals.  This is unlikely to be as efficient as existing algorithms for identities amongst one set of genes across two individuals; it is limited by the need to calculate and store identities amongst very many sets of ancestral genes, corresponding to the very many routes by which genes may descend through the pedigree.

Notation

Genes are defined by individual and gene: 𝕚={i,1} or {i,2}, for two genes in individual i; for compactness, we may write 𝕚={i,1}=(Subscript[i, 1]).  We define the probability F[{(Subscript[S, 1]),(Subscript[S, 2]),…, (Subscript[S, n])}] that the sets of genes (Subscript[S, 1]), (Subscript[S, 2]), …,(Subscript[S, n]) are identical by descent, tracing back to n distinct founders in the ancestral population.  For example,   F[{{𝕚},{𝕛,𝕜},{𝕝,𝕞}}] is the probability that these 3 sets of genes each trace back to 3 distinct founders.  Necessarily, F[{{𝕚}}]=1, and F[{{𝕚},{𝕛}}]=1-F[{{𝕚,𝕛}}] .  Identities in generation t are denoted (Subscript[F, t]).

Given the pedigree, the identity (Subscript[F, t]) is a linear combination of identities in the previous generation,  (Subscript[F, t-1]).  The algorithm identifies all sets of genes in the parents that can contribute to the identities in the specified set in the offspring.  We first outline the algorithm, and then explain it in more detail, in terms of the Mathematica code.

Outline

In generation t=0, (Subscript[F, 0])[S] is, set to 1 if the set S contains zero or more sets each with a single distinct gene (S={{𝕚},{𝕛},…}, 𝕚≠𝕛 ….  Otherwise, it is zero.

A single gene traces back to one or other parent with equal probability; two genes in the same individual must trace back to the two  parents, though those may be the same individual if there is selfing.  This generates a sum over possible sets of parents.

For each parent, the descendant genes must come from one or other of the two genes in that parent. Descendant genes are allocated at random to two sets, generating a list of sets of genes that contribute to the identities.

This list contains many permutations that are equivalent, differing only by order; these are tallied to reduce the number of configurations that need to be stored.

This gives a recursion back to the founder generation. All intermediate identities are stored, which speeds computation. The number of generations and size of pedigree is limited by the amount of memory needed to store intermediate coefficients.

Mathematica code

Setting up a pedigree

Evaluate initialisation sets up the function definitions, defined below.

The pedigree is labelled 𝒫,  and corresponding identities are calculated using Fx[𝒫, t][{(Subscript[S, 1]),…}].

This sets up lists of the parents, over 50 generations, for 10 individuals:

pedigreeList[𝒫]=Table[Parents[randomPedigree[{10,10}]],{50}];

This shows the pedigree in the first 3 generations. For example, the first individual in generation 3 has parents {1,7} in the previous generation.

TableForm[pedigreeList[𝒫][[1;;3]],TableDepth->2]

Simple examples

This is the identity between the two genes in individual 1, at t=3:

Fx[𝒫,3][{{{1,1},{1,2}}}]

This is the probability that the two genes in individual 1 are ibd, and that the two genes in individual 2 are ibd, but that these two sets trace back to different ancestors:

Fx[𝒫,3][{{{1,1},{1,2}},{{2,1},{2,2}}}]

This is a check that for all pairs of individuals {j,k}, F[{{(Subscript[k, 1]),(Subscript[k, 2])}}] = F[{{(Subscript[k, 1]),(Subscript[k, 2])},{(Subscript[j, 1])}}]+F[{{(Subscript[k, 1]),(Subscript[k, 2]),(Subscript[j, 1])}}]:

Table[{Fx[𝒫,5][{{{k,1},{k,2}}}],Fx[𝒫,5][{{{k,1},{k,2}},{{j,1}}}]+Fx[𝒫,5][{{{k,1},{k,2},{j,1}}}]},{j,2,10},{k,1,j-1}]

How the algorithm works

The function Fx[𝒫,t][{(Subscript[S, 1]),(Subscript[S, 2]),…}] calculates and stores identity coefficients amongst multiple one or more sets of genes, (Subscript[S, i])={𝕚,𝕛,…}.

Fx[𝒫,0][{(Subscript[S, 1]),(Subscript[S, 2]),…}] is set to zero if any of the sets are contain more than 1 gene, and zero otherwise. It checks that all the genes are distinct.

For t>0, Fx[𝒫,t][{(Subscript[S, 1]),(Subscript[S, 2]),…}] checks whether all the genes are distinct. If not, the (Subscript[S, i]) are replaced by Union[(Subscript[S, i])], since repeated genes are necessarily identical. If the genes are still not distinct, that implies that the same gene appears in more than one set (Subscript[S, i]), and Fx is set to zero.  If all the genes are now distinct, Fx[𝒫,t] is called.

When Fx[𝒫,t][{(Subscript[S, 1]),(Subscript[S, 2]),…}] is called with entirely distinct genes, there is a check whether the sets {(Subscript[S, 1]),(Subscript[S, 2]),…} are each ordered, and the set of sets is {(Subscript[S, 1]),(Subscript[S, 2]),…} is ordered. If not, the sets and the set of sets is sorted, and Fx[𝒫,t] is called.

When Fx[𝒫,t][{(Subscript[S, 1]),(Subscript[S, 2]),…}] is called with entirely distinct genes, and with ordered {(Subscript[S, 1]),(Subscript[S, 2]),…}, iterF[Fx[𝒫,t-1],{(Subscript[S, 1]),(Subscript[S, 2]),…},{{6,9},…}] is called, to construct the linear recursion back to the identities in the previous generation. The last argument lists the two parents of each individual; {i,i} denotes an individual born by selfing of parent i.

Within iterF, genesToParents[{(Subscript[S, 1]),(Subscript[S, 2]),…},parents] is used to find where each of the genes might have come from. For example, consider the set of genes {(Subscript[1, 1]),(Subscript[1, 2]),(Subscript[2, 1])}.  The two genes (Subscript[1, 1]),(Subscript[1, 2]) in individual 1 must have come one from each parent - i.e., either {6,7} or {7,6}.  The gene (Subscript[2, 1]) could have come from either of the two parents of #2, i.e. from 6 or from 8.  This gives 4 possible sets of parental origins:

gtp=genesToParents[{{{1,1},{1,2},{2,1}}},{{6,7},{6,8}}]

For each of these parental origins, we use parentsToGenes to find which genes within those parents actually contributed.  {6,7,6} and {7,6,6} could have contributed in two ways - either by parent 6 passing on the same allele, or different alleles, to the two offspring.  However, {6,7,8} and {7,6,8} only contribute in one way each. (Note that we use the symmetry between the two alleles, denoting them by label 1 by default).

parentsToGenes/@gtp

Flattening, we have a weighted sum over 6 identities:

fgtp=Flatten[parentsToGenes/@gtp,1]

Weighting each appropriately, we obtain a sum of identities:

Total[(((#[[2]])/(Length[gtp]))f[#[[1]]])&/@fgtp]

iterf[f,{(Subscript[S, 1]),…},{{6,7},…}] achieves the same.  Replacing f by Fx[𝒫,t-1] gives he required recursion:

iterF[Fx[𝒫,t-1],{{{1,1},{1,2},{2,1}}},{{6,7},{6,8}}]

Collapsing equivalent sets

These 4 are in fact equivalent to two distinct sets , so we only need to consider these two sets, each with multiplicity 4:

tgtp=Tally[Map[Sort,gtp,2]]

Now, we apply parentsToGenes in the same way, but adjust the weights by the multiplicity of each distinct parental origin:

ftgtp=Flatten[(ptg=parentsToGenes[#[[1]]];Transpose[{ptg[[All,1]],ptg[[All ,2]]#[[2]]}])&/@tgtp,1]

This gives the same result, but slightly more efficiently:

Total[(((#[[2]])/(Length[gtp]))f[#[[1]]])&/@ftgtp]

References

Abney M., 2009 A graphical algorithm for fast computation of identity coefficients and generalized kinship coefficients. Bioinformatics 25: 1561–1563.

García-Cortés L. A., 2015 A novel recursive algorithm for the calculation of the detailed identity coefficients. Genetics selection evolution 47: 1319–4.

Karigl G., 1981 A recursive algorithm for the calculation of identity coefficients. Annals Of Human Genetics 45: 290–305.

Karigl G., 1982 A mathematical approach to multiple genetic relationships. Theor. Popul. Biol.: 1–15.

Kirkpatrick B., Ge S., Wang L., 2018 Efficient computation of the kinship coefficients (O Stegle, Ed.). Bioinformatics 35: 1002–1008.

Definitions

pedigreeList

pedigreeList::usage="pedigreeList[𝒫] stores a list of pedigrees associated with the symbol 𝒫.  This allows functions such as Fx to refer to the list without storing it.  It needs to be defined before use.  ** Here, a pedigree is represented as {{(Subscript[par, 1]),(Subscript[par, 2])},…}, using Parents[]";

randomPedigree

randomPedigree::usage="randomPedigree[{n,p}] makes a random WF pedigree (but with no selfing), with n offspring and p parents. randomPedigree[{{(Subscript[n, 1]),(Subscript[n, 2]),…,(Subscript[n, k])},{(Subscript[p, 1]),…,(Subscript[p, k])}},{{(Subscript[m, 1, 1]),(Subscript[m, 1, 2]),…},…}] makes a random pedigree for the k-island model, with (Subscript[n, i]) offspring, (Subscript[p, i]) parents, and a fraction (Subscript[m, i, j]) in deme i coming from deme j. randomPedigree[{n,p},{{(Subscript[m, 1, 1]),(Subscript[m, 1, 2]),…},…}] assumes that all k demes have n offspring and p parents. There must be at least two parents per deme to avoid selfing. Returns a SparseArray. Selfing->True allows selfing, giving the WF model - only defined for one deme as yet";

Selfing::usage="Selfing is an option for randomPedigree that allows selfing";

Options[randomPedigree]={Selfing->False};

randomPedigree[{n_Integer,p_Integer},opts___Rule]:=Module[{c=ConstantArray[0,p],rp},
If[Selfing/.{opts}/.Options[randomPedigree],
SparseArray[Table[rp=RandomInteger[{1,p},{2}];
If[rp[[1]]==rp[[2]],
ReplacePart[c,{{rp[[1]]}}->2],
ReplacePart[c,(List/@rp)->1]],{n}]],
SparseArray[Table[ReplacePart[c,(List/@RandomSample[Range[p],2])->1],{n}]]]];

randomPedigree[{n_Integer,p_Integer},m_List,opts___Rule]:=Module[{k=Length[m]},
randomPedigree[{ConstantArray[n,k],ConstantArray[p,k]},m,opts]];

randomPedigree[{n:{___Integer},p:{___Integer}},m_List]:=Module[{demes=Range[Length[m]],ds,pars,c=ConstantArray[0,Total[p]],rp},
rp=partitionList[p];
ds=Flatten[MapThread[RandomChoice[#1->demes,{#2,2}]&,{m,n}],1];
pars=If[#[[1]]==#[[2]],RandomSample[rp[[#[[1]]]],2],{RandomChoice[rp[[#[[1]]]]],RandomChoice[rp[[#[[2]]]]]}]&/@ds;
SparseArray[ReplacePart[c,(List/@#)->1]&/@pars]];

genesToParents

genesToParents::usage="genesToParents[{{{1,1},{3,2}},…},{{3,4},…}] gives sets of IBD genes that point to the diploid parents. The first argument lists sets of genes that are IBD within, but not between, sets. The second argument lists parents, numbered consecutively. The result is expressed as terms of the form {{4,5},…}, which are weighted equally. ";

genesToParents[g:{___List},pars:{{_Integer,_Integer}..}]:=
If[OrderedQ[g]∧(And@@(OrderedQ/@g)),
Module[{inds,inds1,inds2,h,f,n1,n2},
inds=First/@Sort[Flatten[g,1]];
inds1=Select[Union[inds],Count[inds,#]==1&];
inds2=Select[Union[inds],Count[inds,#]==2&];
n1=Length[inds1];n2=Length[inds2];
h=Flatten[Outer[List,
Sequence@@pars[[inds1]],Sequence@@ConstantArray[{1,2},n2]],n1+n2-1];
f[s:{___Integer}]:=Join[
MapThread[{#1,_}->#2&,{inds1,Take[s,n1]}],
Flatten[MapThread[{{#1,1}->pars[[#1,#2]],{#1,2}->pars[[#1,3-#2]]}&,{inds2,Take[s,-n2]}],1]];Replace[g,f[#],{2}]&/@h],
genesToParents[Sort[Sort/@g],pars](* put into ordered form *)];

parentsToGenes

parentsToGenes::usage="parentsToGenes[{{1,3,5},{2,5}}] gives the joint probability of identity between genes drawn randomly from the specified individuals, as a list {{{…},w},…}.  parentsToGenes[{{{1,3,5},{2,5}},w}] returns a weighted list. For each distinct individual, genes are partitioned randomly into two sets, and the corresponding sets are merged if genes come from the same set. ";

parentsToGenes[g:{___List}]:=Module[{inds,sets,spl},
inds=Union[Flatten[g]];
sets=(First/@Position[g,#])&/@inds;
spl=makeAllSplits[sets];
{Replace[gatherSets[#[[1]]],{j_Integer,k:(1|2)}:>{inds[[j]],k},{2}],#[[2]]}&/@spl];

distinctQ

distinctQ::usage="distinctQ[{{a,b},{c,d}…}] tests whether the elements {a,b,c,d …} are distinct";

distinctQ[s:{___List}]:=Module[{sl},
sl=Flatten[s,1];Length[sl]==Length[Union[sl]]];

Fx

Fx::usage="Fx[𝒫,t][{{{1,1},{3,2}},…}] stores the probability that the listed sets of genes are all IBD within sets, but distinct from each other. This assumes that pedigreeList[𝒫] lists of parents for each individual, numbered consecutively, for each generation. The result is expressed in terms of a sum over Fx[𝒫,t][{{{4,1},{5,2}},{…},…},…}], where the sets point to genes in diploid individuals. Fx[𝓅,0][] is set to 0 or 1 as appropriate";

Fx[_,0][g:{___List}]:=
If[distinctQ[g],
If[Or@@(Length[#]>1&/@g),0.,1.],
Message[Fx::multipleGenes,g];$Failed];
(* If the sets are all distinct individuals, and all of length 1 (or 0), then F is 1. The initial F is set to a numerical value to avoid Mathematica getting bogged down in exact arithmetic *)

Fx[𝒫_,t_Integer][g:{___List}]:=Fx[𝒫,t][g]=Module[{g1},
If[distinctQ[g],
If[g=={},1,
If[OrderedQ[g]∧(And@@(OrderedQ/@g)),
Fx[𝒫,t][g]=
iterF[Fx[𝒫,t-1],g,pedigreeList[𝒫][[t]]],
Fx[𝒫,t][Sort[Sort/@g]](* put into ordered form *)]],
(* if genes appear more than once within sets, replace by the union *)
g1=Union/@g;
(* if the same gene appears in multiple sets, then F=0 *)
If[distinctQ[g1],Fx[𝒫,t][g1],0]
]];

Fx::multipleGenes="Genes appear more than once in `1`";

iterF

iterF::usage="iterF[f,{{{1,1},{3,2}},…},{{1,3},…}] iterates IBD amongst sets of genes, given a list of parents of diploid individuals, returning the result as an expression in the form ∑ w*f[{{{1,1},…},…}]. Fx sets f to Fx[𝒫,t-1]";

iterF[f_,g:{___List},pars:{{_Integer,_Integer}..}]:=Module[{gtp,tgtp,ptg},
gtp=genesToParents[g,pars];
(* The original version (below) is inefficient, because gtp contains many equivalent sets.
Total[(((#[[2]])/(Length[gtp]))f[#[[1]]])&/@Flatten[parentsToGenes/@gtp,1]*)
tgtp=Tally[Map[Sort,gtp,2]];(* tally equivalent permutations *)
Total[(((#[[2]])/(Length[gtp]))f[#[[1]]])&/@Flatten[(ptg=parentsToGenes[#[[1]]];Transpose[{ptg[[All,1]],ptg[[All ,2]]#[[2]]}])&/@tgtp,1]]
(* for each permutation generated by genesToParents, produce the list of IBD sets, and weight by the multiplicity of that permutation *) ];

makeAllSplits

makeAllSplits::usage="makeAllSplits[{{1,2},…,{3,4}}] takes a list of the sets that descend from each individual, and lists all the ways that these can be validly assigned, together with the weights, (Subscript[w, i]).  The result is a list in the form {{{{{1},{2}},{{},{3}},…},(Subscript[w, 1])},…}, where the first entry lists the sets which descend from the two genes in each individual. ";

makeAllSplits[sets_List]:=Module[{aa,bb,ts,wts},
aa=(ts=Tally[Sort/@findSplit[#]];ts/.{a_List,j_Integer}:>{Union/@a,j (2^(-Length[#]))})&/@sets;
bb=Flatten[Outer[List,Sequence@@aa,1],Length[aa]-1];
wts=(Times@@#[[All,2]])&/@bb;
Transpose[{bb[[All,All,1]],wts}]];

For example, suppose that we have four sets of IBD genes, two of which descend from one individual, and two from the other.  This can only happen if 1 & 2 descend from the different genes in one individual, and 3 & 4 from the different genes in the other, which has probability 1/4:

makeAllSplits[{{1,2},{3,4}}]

gatherSets

gatherSets::usage="gatherSets[{{{1},{2}},…,{{},{3}}}] takes a list of the sets of genes that descend from each individual (Subscript[i, 1]),…,(Subscript[i, n]), and gathers together sets of genes that are IBD. An entry {{x},{y}} for individual a corresponds to genes {a,1} in set x and {a,2} in set y, whilst {{},{x}} corresponds to a gene {a,1} in set x.";

gatherSets[g:{{{___Integer},{_Integer}}...}]:=Module[{n=Length[g],xx,ss},
xx=Flatten[Table[g[[j]]/.{{{},{s_Integer}}:>{{s,{j,1}}},{{s1_Integer},{s2_Integer}}:>{{s1,{j,1}},{s2,{j,2}}}},{j,n}],1];
ss=Union[First/@xx];(Cases[xx,{#,_}]&/@ss)[[All,All,2]]];

findSplit

findSplit::usage="findSplit[{1,2,2}] finds valid partitions into two sets, out of the total (2^n)";

findSplit[x:{___Integer}]:=Select[{#,TakeAway[x,#]}&/@Subsets[x],(Length[Union[#[[1]]]]≤1)∧(Length[Union[#[[2]]]]≤1)&];

TakeAway

TakeAway::usage="TakeAway[U,V]  gives the set U, minus elements V. Note that this is not the same as Complement[U,V] if U contains repeated elements.  ";

TakeAway[U_List,{}]:=U;
TakeAway[U_List,{v_}]:=Drop[U,Position[U,v][[1]]]/;Complement[{v},U]=={};
TakeAway[U_List,V_List]:=TakeAway[TakeAway[U,Take[V,1]],Drop[V,1]]/;Complement[V,U]=={};
TakeAway[U_List,W__List,V_List]:=TakeAway[TakeAway[U,W],V];