pedkin

pedkin's algorithm

recursive algorithm

The algorithm used by pedkin is known in the literature as the 'recursive algorithm'. Given individuals a and b, where a is not an ancestor of b, it defines their kinship coefficient k(a,b) as follows:

• If a or b is a founder, then k(a,b) = 0 (by definition)
• Otherwise, for a single individual a with parents fa and ma, k(a,a) = (1 + k(fa,ma))/2 (kinship with self is half of one plus inbreeding coefficient)
• Otherwise, k(a,b) = (k(fa,b) + k(ma,b))/2 (kinship with other is 'inherited' from parents, each contributing half)

Mind the condition "a is not an ancestor of b". We need this because if a were an ancestor of b, then recursing along the parents of a will never get us to b. Swapping them and starting from b solves this for any pair of individuals. Obviously then, the condition can always be met.

implementation

The algorithm described above is easily implemented in code, once you have ordered all the individuals in such a way that no individual can be an ancestor of an individual preceding it.

Below is the C++ code from kcoeff.cpp which implements the algorithm for pedkin. Parameters p1 and p2 are the sequence numbers of two individuals (0 meaning founder); p is the pedigree object, which can return the parents for any indivual.

double kcoef(unsigned p1, unsigned p2, const Pedigree &p) 
{
    if (p1 == 0 || p2 == 0)
        return 0;
    else if (p1 == p2)
        return (1 + kcoef(p.pa(p1), p.ma(p1), p)) / 2;
    else if (p1 < p2)
        return (kcoef(p1, p.pa(p2), p) + kcoef(p1, p.ma(p2), p)) / 2;
    else 
        return (kcoef(p2, p.pa(p1), p) + kcoef(p2, p.ma(p1), p)) / 2;
}

efficiency

The algorithm is pretty memory efficient since it does not store any intermediate results. The downside to this is that it is not very computation efficient; it calculates the same coefficients over and over. The implementation could be made a lot faster by caching coefficients once they have been computed.

The speed gain would go at the expense of memory consumption, which rises quadratically with the number of individuals in the pedigree. Assuming 8 bytes for a double, a pedigree of size 10,000 would require nearly 0.8 GB, which is still a considerable amount of RAM these day. Though there are ways around that, I didn't bother to dive into them. The program runs fast enough as it is, and I like the elegance and simplicity of the current solution.