Supplementary material for:

Effective number of breeders from sibship reconstruction: empirical evaluations using hatchery steelhead

This supplementary notebook was written by Ryan Waples (ryan.waples@gmail.com)

This document originated as a Jupyter notebook running Python 2.7.

In [1]:
import numpy as np
from __future__ import division
import itertools
import collections

Parentage analysis without parents (PWOP)

Waples et. al. (2011) equation 2a, also Supplemental equation 2 in the current manuscipt.

Produces an Ne estimate from a vector of k_i values using the equation:

\begin{equation} N_e = \frac{\sum_{}^{} k_i - 1}{\frac{\sum_{}^{} k_i^2}{\sum_{}^{} k_i} -1 } \end{equation}

Where k_i is the number of offspring produced by parent i.

In [2]:
def pwop(k_i):
    assert type(k_i) is np.ndarray
    print '{} contributing parents'.format(len(k_i))
    print '{} gametes'.format(np.sum(k_i))
    numerator = np.sum(k_i)-1
    denominator = np.sum(k_i**2)/np.sum(k_i) - 1
    Ne = numerator/denominator 
    return(Ne)

example

In [3]:
K_i = np.array([2,2,3,7,2,5,1])
print('The PWOP Ne estimate is: {:1.3f}'.format(pwop(K_i)))

7 contributing parents
22 gametes
The PWOP Ne estimate is: 6.243

p_same

calculates the probability that two gametes sampled at random from a k_i vector share the same parent.

\begin{equation} p\_same = \frac{1}{N(N-1)} \sum_{i=1}^{N} (k_i(k_i-1)) \end{equation}

Where N is the number of parent (or offspring) gametes.

In [4]:
def p_same(k_i):
    """chance that two random gametes present in the offspring generation are from the same parent"""
    N = np.int(np.sum(k_i)) # 2x number of offspring or parents
    print '{} contributing parents'.format(len(k_i))
    print '{} gametes'.format(N)
    p = 1./(N*(N-1)) * np.sum(k_i*(k_i-1))
    return(p)
In [5]:
p_same(np.array([4,2]))

2 contributing parents
6 gametes
Out[5]:
0.4666666666666667

example

In [6]:
# use the same k_i vector
print(K_i)
print('The chance two gametes share a parent: {:1.3f}'.format(p_same(K_i)))

[2 2 3 7 2 5 1]
7 contributing parents
22 gametes
The chance two gametes share a parent: 0.160

COLONY method of Ne estimation.

Wang (2009): "Equations for the effective size ($N_e$) of a population were derived in terms of the frequencies of a pair of offspring taken at random from the population being sibs sharing the same one or two parents"

Wang (2009) eq. 10:

\begin{equation} \frac{1}{N_e} = \frac{1+3\alpha}{4}(Q_1+Q_2+2Q_3)-\frac{\alpha}{2}(\frac{1}{N_1}+\frac{1}{N_2}) \end{equation}

Where $\alpha$ is a measure of departure from Hardy-Weinberg equilibrium, $Q_1$, $Q_2$, and $Q_3$ are the probabilities of a pair of offspring being paternal half-siblings, maternal half-siblings, and full-siblings, respectively. $N_1$ and $N_2$ are the number of male and female parents.

We will assume alpha = 0 (i.e. complete HWE). The equation simplifies to:

\begin{equation} \frac{1}{N_e} = \frac{1}{4}(Q_1+Q_2+2Q_3) \end{equation}

(Q1+Q2+2*Q3)/4 is eqivalent to p_same (Wang 2009, equation 8), so we simplify to: \begin{equation} \frac{1}{N_e} = p_s \end{equation}

But notice the impact of the simplifications, we lose all consideration of separate sexes and assuming alpha = 0, may not be reasonable, especially for small populations.

In [7]:
def COLONY(p_same_parent):
    """Given the chance that two random gametes share a parent, return the Ne"""
    Ne = 1./p_same_parent
    return(Ne)
In [8]:
print('The COLONY Ne estimate is: {:1.3f}'.format(COLONY(p_same(K_i))))

7 contributing parents
22 gametes
The COLONY Ne estimate is: 6.243

Simulate an ideal mating population

Generates a K_i vector from one generation of random mating of N parents producing N offspring (total of 2N gametes). The result excludes parents that do not contribute offspring.

In [9]:
def sim_ideal(N):
    """Generates a K_i vector from one generation of random mating of N parents"""
    a = collections.defaultdict(int)
    for xx in np.random.randint(0, N, 2*N): # For each of 2N gametes select a parent
        a[xx]+=1 # assign the gamete to that parent
    k_i = np.array(a.values())
    print '{} contributing parents'.format(len(k_i))
    print '{} gametes'.format(np.sum(k_i))
    return(k_i)
In [10]:
K_i = sim_ideal(200)
print (K_i)

168 contributing parents
400 gametes
[1 1 2 3 1 1 3 4 1 3 3 1 2 1 2 2 3 2 5 1 2 1 1 3 2 1 4 3 2 2 1 1 2 1 5 2 3
 2 4 1 1 1 1 3 3 3 1 1 1 2 2 1 3 3 1 1 2 1 2 3 5 3 2 3 3 2 2 2 2 4 3 7 2 4
 3 2 4 4 1 1 3 1 2 3 5 2 3 1 2 2 3 2 2 1 1 4 1 6 1 2 3 2 3 5 2 2 3 1 4 2 1
 5 3 1 3 3 6 3 2 5 2 2 2 3 1 2 3 4 1 2 6 1 4 4 1 5 2 4 1 4 1 2 1 2 3 4 2 2
 1 2 1 2 3 4 3 2 3 2 2 1 4 1 1 3 2 1 3 2]

Apply both methods to the simulated k_i vector

In [11]:
p_same(K_i)

168 contributing parents
400 gametes
Out[11]:
0.0051879699248120305
In [12]:
pwop(K_i)

168 contributing parents
400 gametes
Out[12]:
192.75362318840581
In [13]:
COLONY(p_same(K_i))

168 contributing parents
400 gametes
Out[13]:
192.75362318840578

Also compare to Crow (1954) & Crow and Denniston (1988) eq 1:

\begin{equation} \frac{1}{N_e} = \frac{ \bar{k_i} - 1 + \frac{Var(k_i)}{\bar{k_i}} } {N \bar{k_i} -1} \end{equation}
In [14]:
def Crow_Denniston(k_i):
    top = np.mean(k_i) - 1 + np.var(k_i)/np.mean(k_i)
    bottom = len(k_i)*np.mean(k_i)-1
    Ne = bottom/top
    return(Ne)
In [15]:
Crow_Denniston(K_i)
Out[15]:
192.75362318840581

References

Crow, James F., and Carter Denniston. "Inbreeding and variance effective population numbers." Evolution (1988): 482-495.

Jones, Owen R., and Jinliang Wang. "COLONY: a program for parentage and sibship inference from multilocus genotype data." Molecular ecology resources 10.3 (2010): 551-555.

Wang, Jinliang. "A new method for estimating effective population sizes from a single sample of multilocus genotypes." Molecular Ecology 18.10 (2009): 2148-2164.

Waples, Robin S., and Ryan K. Waples. "Inbreeding effective population size and parentage analysis without parents." Molecular ecology resource 11.s1 (2011): 162-171.

In [ ]: