{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"#1. Linkage Disequilibrium in an Infinite Population\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Gametic disequilibrium, which is more commonly referred to as linkage\n",
"disequilirium (LD), is the statistical dependence between alleles in a\n",
"haplotype. \n",
"\n",
"Under gametic equilibrium, alleles in a haplotype are\n",
"independent. This is also called linkage equilibrium. \n",
"\n",
"Note that it is\n",
"possible for two loci that are linked to be in gametic equilibrium;\n",
"also, loci that are unlinked can be in gametic disequilibrium.\n",
"\n",
"Suppose that starting from generations 0, all individuals are produced\n",
"by random mating. \n",
"\n",
"Then, the probability of haplotype ($A_{i},B_{j}$) in\n",
"generation 1 can be written as\n",
"\n",
"$${\\textstyle \\Pr_{1}(A_{i},B_{j})=(1-r)\\Pr_{0}(A_{i},B_{j})+r\\Pr(A_{i})\\Pr(B_{j})}$$\n",
"\n",
"where $r$ is the recombination rate between loci $A$ and $B$, and\n",
"$\\Pr_{0}(A_{i},B_{j})$ is the probability of ($A_{i},B_{j}$) in\n",
"generation 0. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"The disequilibrium in generation 1 is\n",
"\n",
"$$\\begin{split}\n",
"\\Delta_{1} & ={\\Pr}_{1}(A_{i},B_{j})-\\Pr(A_{i})\\Pr(B_{j})\\\\\n",
" & =(1-r){Pr}_{0}(A_{i},B_{j})+r\\Pr(A_{i})\\Pr(B_{j})-\\Pr(A_{i})\\Pr(B_{j})\\\\\n",
" & =(1-r){Pr}_{0}(A_{i},B_{j})-(1-r)\\Pr(A_{i})\\Pr(B_{j})\\\\\n",
" & =(1-r)\\Delta_{0}\n",
"\\end{split}$$ \n",
"\n",
"where $\\Delta_{0}$ is the disequilibrium in generation 0."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Similarly, the probability of haplotype ($A_{i},B_{j}$) in generation 2\n",
"is\n",
"\n",
"$${\\textstyle \\Pr_{2}(A_{i},B_{j})=(1-r)\\Pr_{1}(A_{i},B_{j})+r\\Pr(A_{i})\\Pr(B_{j})}$$\n",
"\n",
"and the disequilibrium in generation 2 is\n",
"\n",
"$$\\begin{split}\n",
"\\Delta_{2} & ={\\Pr}_{2}(A_{i},B_{j})-\\Pr(A_{i})\\Pr(B_{j})\\\\\n",
" & =(1-r){\\mbox{Pr}}_{1}(A_{i},B_{j})+r\\Pr(A_{i})\\Pr(B_{j})-\\Pr(A_{i})\\Pr(B_{j})\\\\\n",
" & =(1-r){\\mbox{Pr}}_{1}(A_{i},B_{j})-(1-r)\\Pr(A_{i})\\Pr(B_{j})\\\\\n",
" & =(1-r)\\Delta_{1}\\\\\n",
" & =(1-r)^{2}\\Delta_{0}\n",
"\\end{split}$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"It follows that in generation $n$, the disequilibrium is\n",
"\n",
"
\n",
"$$\\Delta_{n}=(1-r)^{n}\\Delta_{0}$$\n",
"\n",
"Thus, with each generation of random\n",
"mating the haplotype distribution moves closer to equilibrium\n",
"(statistical independence). \n",
"\n",
"For loci that are unlinked, $(1-r)=1/2$ and\n",
"equilibrium is reached quickly; for example, $(1/2)^{10}=1/1024$. \n",
"\n",
"On the\n",
"other hand, loci that are tightly linked will take much longer to reach\n",
"equilibrium; for example, $(1-r)^{10}>1/3$ for $r=0.1$. \n",
"\n",
"In the limit,\n",
"however, equilibrium is reached in an infinite population."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"#2. Linkage Disequilibrium in a Finite Population\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"In a closed finite population, in the absence of mutation, after a\n",
"sufficient number of generations of random mating all alleles will\n",
"become identical by descent; \n",
"\n",
"thus all alleles also will become identical\n",
"in state. \n",
"\n",
"In other words, all loci will become “fixed”.\n",
"\n",
"In such apopulation genetic variability is absent and LD is not defined."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Loci that are moving toward fixation will pass through a phase where\n",
"alleles that are identical by state (IBS) will also be identical by\n",
"descent (IBD). \n",
"\n",
"In other words, all alleles that are IBS will trace back\n",
"to a common ancestral mutant allele. \n",
"\n",
"Thus, at a biallelic locus A with\n",
"alleles $A_{1}$ and $A_{2}$, all the $A_{1}$ alleles will have a common\n",
"ancestor and all the $A_{2}$ alleles will have a different common\n",
"ancestor. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"In such loci, Sved (1971) has shown algebraically\n",
"that the expected value of LD between loci A and B as measured by the\n",
"squared correlation ($\\rho^{2}$) between allele states is related to the\n",
"probability of joint identity-by-descent at loci A and B for two\n",
"randomly sampled gametes. \n",
"\n",
"At first, this relationship may not be\n",
"obvious. \n",
"\n",
"To get an intuitive feel for this relationship, consider a\n",
"population where alleles $A_{1}$ and $A_{2}$ are segregating at the A\n",
"locus and alleles $B_{1}$ and $B_{2}$ at the B locus. \n",
"\n",
"Suppose all\n",
"haplotypes with $A_{1}$ also have $B_{1}$ and those with $A_{2}$ have\n",
"$B_{2}$. Then, $\\rho$ would be 1. \n",
"\n",
"If the allelic associations were\n",
"reversed, i.e., $A_{1}$ goes with $B_{2}$ and $A_{2}$ goes with $B_{1}$,\n",
"then $\\rho$ would be -1. \n",
"\n",
"In both cases $\\rho^{2}$ is 1. In such a\n",
"population for two randomly sampled gametes, the conditional probability\n",
"would be one that alleles at the B locus are identical-by-state (IBS)\n",
"given they are IBS at the A locus. \n",
"\n",
"On the other hand, if there are a few\n",
"haplotypes where the association between alleles is different from the\n",
"other haplotypes, then $\\rho^{2}$ would be less than 1. \n",
"\n",
"Further, in this\n",
"population for two randomly sampled gametes, the conditional probability\n",
"would be less than 1 that alleles at the B locus are identical-by-state\n",
"(IBS) given they are IBS at the A locus.\n",
"\n",
"Hopefully, this discussion\n",
"helps to see that LD as measured by $\\rho^{2}$ should be related to the\n",
"probability of joint IBS.\n",
"\n",
"Recall, however, that we assumed that all\n",
"$A_{1}$ alleles have descended from a common ancestor and similarly all\n",
"$A_{2}$ alleles have descended from a different common ancestor. \n",
"\n",
"Thus,\n",
"given alleles at the A locus are IBS, they have descended from a common\n",
"ancestor (IBD), and provided that no recombination between the two loci\n",
"has happened in the two paths descending from the common ancestor to the\n",
"two randomly sampled haplotypes, the B alleles will also be IBD. \n",
"\n",
"Sved\n",
"(1971) denoted this conditional probability by $Q$, and\n",
"reasoned that in the pool of gametes where alleles at the B locus are\n",
"IBD given they are IBD at the A locus, LD as measured by the squared\n",
"correlation ($\\rho^{2}$) between allele states will be 1.0\n",
"(Sved, 2009). \n",
"\n",
"On the other hand, in the pool of gametes where\n",
"recombination has taken place between loci A and B, given random mating,\n",
"$\\rho^{2}$ is expected to be null (Sved, 2009).\n",
"\n",
"Let $C=1$ denote the condition that for a randomly sampled pair of\n",
"gametes the alleles at locus B are IBD given that alleles are IBD at\n",
"locus A, and $C=0$ denote that this condition is not met. \n",
"\n",
"Then, the\n",
"expected value of $\\rho^{2}$ can be written as\n",
"\n",
"(1)\n",
"$$\\begin{aligned}\n",
"E(\\rho^{2}) \n",
"& = \\underset{C}{E}[E(\\rho^{2}|C)]\\\\\n",
"& = E(\\rho^{2}|C=1)\\Pr(C=1)+E(\\rho^{2}|C=0)\\Pr(C=0)\\\\\n",
"& = 1Q+0(1-Q)\\\\\n",
"& = Q.\n",
"\\end{aligned}$$\n",
"\n",
"Let $Q_{t}$ denote the probability of $C=1$ in generation $t$. This\n",
"probability can be recursively written as\n",
"\n",
"(3)\n",
"$$Q_{t}=[\\frac{1}{2N}+(1-\\frac{1}{2N})Q_{t-1}](1-r)^{2}$$\n",
"\n",
"where $N$ is the effective population size, \n",
"\n",
"$\\frac{1}{2N}$ is the\n",
"probability that two randomly sampled gametes in generation $t$ are both\n",
"inherited from the same gamete in the previous generation,\n",
"\n",
"(1-$\\frac{1}{2N}$) is the probability that they are inherited from\n",
"different gametes in the previous generation, \n",
"\n",
"and $(1-r)^{2}$ and the\n",
"probability that loci A and B do not recombine in these gametes in the\n",
"last generation. \n",
"\n",
"At equilibrium, $Q_{t}=Q_{t-1}=Q_{E}$. \n",
"\n",
"So, setting\n",
"$Q_{t}$ and $Q_{t-1}$ in to $Q_{E}$ gives (3)\n",
"\n",
"$$\\begin{split}\n",
"Q_{E} & =\\frac{(1-r)^{2}}{2N-(2N-1)(1-r)^{2}}\\\\\n",
" & \\approx\\frac{1}{4Nr+1}\n",
"\\end{split}$$\n",
"\n",
"In the derivation of (3) we only considered pairs of loci where\n",
"alleles are segregating at each locus. \n",
"\n",
"Further, it was assumed that all\n",
"haplotypes in the current generation descend from two ancestral\n",
"haplotypes. \n",
"\n",
"So, if we code the four possible haplotypes at a pair of\n",
"biallelic loci as: 00, 01, 10, an 11, the ancestral pair of haplotypes\n",
"must be either (00,11) or (01,10). \n",
"\n",
"Any other pair would lead to one\n",
"locus being fixed. \n",
"\n",
"For example, the pair (00,01) has locus one fixed for\n",
"the allele coded as 0. \n",
"\n",
"So, if $4Nr$ is close to zero, most haplotypes in\n",
"the current generations will be non-recombinants of the ancestral type\n",
"and $\\rho^{2}$ is close to 1. \n",
"\n",
"The consequences of a few recombinants are\n",
"examined below using the following Julia function."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"RSqr (generic function with 1 method)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function RSqr(nij) \n",
" N = sum(nij)\n",
" Exy = nij[4]/N\n",
" Ex = (nij[3] + nij[4])/N\n",
" Ey = (nij[2] + nij[4])/N\n",
" Vx = Ex * (1 - Ex)\n",
" Vy = Ey * (1 - Ey)\n",
" Cxy = Exy - Ex * Ey\n",
" res = Cxy^2/(Vx * Vy)\n",
" return(res)\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Here we only have non-recombinants:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"r^2 = 1.000 \n"
]
}
],
"source": [
"nij = [80, # ancestral haplotype 00\n",
" 0, # recombinant 01\n",
" 0, # recombinant 10\n",
" 20] # ancestral haplotype 11\n",
"@printf \"r^2 = %5.3f \\n\" RSqr(nij)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Now we introduce two recombinants:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"r^2 = 0.874 \n"
]
}
],
"source": [
"nij = [80, # ancestral haplotype 00\n",
" 1, # recombinant 01\n",
" 1, # recombinant 10\n",
" 18] # ancestral haplotype 11\n",
"@printf \"r^2 = %5.3f \\n\" RSqr(nij)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Here is another example, where one of the ancestral haplotypes has a\n",
"very low frequency:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"r^2 = 0.495 \n"
]
}
],
"source": [
"nij = [98, # ancestral haplotype 00\n",
" 1, # recombinant 01\n",
" 0, # recombinant 10\n",
" 1] # ancestral haplotype 11\n",
"@printf \"r^2 = %5.3f \\n\" RSqr(nij)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"In a finite population, however, most loci are fixed. \n",
"\n",
"Then, LD is not\n",
"defined for these loci. \n",
"\n",
"When mutation introduces variability into such a\n",
"locus, LD is defined but will be low. \n",
"\n",
"This is demonstrated in the\n",
"following example:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"r^2 = 0.002 \n"
]
}
],
"source": [
"nij = [80, # ancestral haplotype 00\n",
" 19, # recombinant 01\n",
" 1, # recombinant 10\n",
" 0] # ancestral haplotype 11\n",
"@printf \"r^2 = %5.3f \\n\" RSqr(nij)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"At mutation-drift equilibrium, most loci will be of this type, where\n",
"mutation has recently introduced variability. \n",
"\n",
"Thus, $E(\\rho^{2})$ can be\n",
"much lower in a population that has reached mutation-drift equilibrium\n",
"than indicated by (3). \n",
"\n",
"To examine this further, the exact\n",
"distribution of is $\\rho^{2}$ is recursively computed next."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## 3. Distribution of $\\rho^{2}$ in the Presence of Mutation\n",
"\n",
"\n",
"Computing the distribution of $\\rho^{2}$ involves computing the joint\n",
"distribution for allele frequencies at two loci. \n",
"\n",
"Thus, we will review\n",
"first how to compute the distribution for allele frequency at a single\n",
"locus."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### Computing the distribution of allele frequency at a single locus"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Consider a population of $2N$ gametes. \n",
"\n",
"Let $Y$ be the number of $A_{1}$\n",
"alleles at locus $A$. \n",
"\n",
"The value of $Y$ can take one of $2N+1$ values\n",
"ranging from 0 to $2N$. \n",
"\n",
"Suppose the distribution of allele frequency in\n",
"generation $t$ is given by the vector\n",
"$\\mathbf{p}_t$ with $2N+1$ probabilities\n",
"corresponding to each of the $2N+1$ possible values of $Y$. \n",
"\n",
"To model\n",
"random mating, assume $2N$ gametes are sampled with replacement from the\n",
"gametes of generation $t$. \n",
"\n",
"Then, ignoring mutation, migration and\n",
"selection, the distribution of allele frequency in generation $t+1$ can\n",
"be calculated as\n",
"\n",
"$$\\mathbf{p}_{t+1}=\\mathbf{Bp}_{t},$$\n",
"\n",
"where $\\mathbf{B}$ is a $(2N+1)\\times(2N+1)$ matrix with element $i,j$ containing the\n",
"probability that a random variable from a Binomial($2N,\\frac{j}{2N}$)\n",
"distribution would be equal to $i$ for $i,j=0,1,2,\\ldots,2N$. \n",
"\n",
"If this is\n",
"not obvious, section 3.10.2 of the notes given\n",
"[here](PopQuantGen.pdf) may\n",
"be useful.\n",
"\n",
"To model mutation, assume that an $A_{1}$ allele mutates to an $A_{2}$\n",
"with probability $u$ and an $A_{2}$ mutates to an $A_{1}$ with\n",
"probability $v$. \n",
"\n",
"Now, to accommodate mutation in computing the\n",
"distribution of allele frequency, $\\mathbf{B}$ is modified such that column $j$\n",
"contains probabilities from the binomial distribution\n",
"$$\\text{Binomial}[2N,\\frac{j}{2N}(1-u)+(1-\\frac{j}{2N})v].$$ \n",
"\n",
"Selection\n",
"can be similarly accommodated by modifying the binomial probabilities\n",
"for each $j$. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"# # Julia Implementation"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"using Distributions\n",
"Ne = 25\n",
"freq0 = 0.5\n",
"nGen = 10\n",
"u = 0.001 # mutation rate\n",
"coeffSel = 0.5\n",
"N = 2*Ne\n",
"p = Array(Float64,N+1,1)\n",
"B = Array(Float64,N+1,N+1)\n",
"x = [0:N];\n",
"for (j in 0:N) \n",
"\tprob = j/N\n",
"\tpMutation = prob*(1-u) + (1-prob)*u\n",
"\tq = 1 - pMutation\n",
"\tq1 = (q - 0.5*coeffSel*q - 0.5*coeffSel*q*q)/(1 - coeffSel*q);\n",
" d = Binomial(N, 1-q1) \n",
" B[:,j+1] = pdf(d,x)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"res = zeros(length(p),nGen)\n",
"d = Binomial(N, 0.5) \n",
"p = pdf(d,x) # initial distribution of gene frequency\n",
"res[:,1] = p\n",
"#repeat for the number of generations\n",
"for (i in 2:nGen)\n",
"\tp = B*p\n",
" res[:,i]=p\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Plot(...)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Gadfly\n",
"plot(layer (x=x/N,y=res[:,1], Geom.point), \n",
"layer (x=x/N,y=res[:,2], Geom.line),\n",
"layer (x=x/N,y=res[:,3], Geom.point),\n",
"layer (x=x/N,y=res[:,4], Geom.line), \n",
"layer (x=x/N,y=res[:,5], Geom.point),\n",
"layer (x=x/N,y=res[:,6], Geom.line), \n",
"Guide.title(\"Distribution of Gene Frequency\"),\n",
"Guide.ylabel(\"Frequency\"),\n",
"Guide.xlabel(\"Gene Frequency\"))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### Computing the joint distribution of allele frequencies at a two linked loci"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Consider a locus $A$ with alleles $A_{1}$ and $A_{2}$ and a linked locus\n",
"$B$ with alleles $B_{1}$ and $B_{2}$. \n",
"\n",
"A population of $2N$ gametes is\n",
"now characterized by a vector $\\mathbf{Y}$ with\n",
"four elements containing the numbers of gametes with haplotypes:\n",
"$A_{1}B_{1}$, $A_{1}B_{2}$, $A_{2}B_{1}$, and $A_{2}B_{2}$. \n",
"\n",
"Note that\n",
"these four numbers must sum to $2N$. \n",
"\n",
"Let $\\mathbf{X}$ be a $k\\times4$ matrix with\n",
"each row representing a possible value of\n",
"${Y}$, where the number $k$ of rows in $\\mathbf{X}$ is\n",
"equal to \n",
"\n",
"$$k=\\frac{(2N+3)!}{3!(2N)!}.$$ \n",
"\n",
"As before let\n",
"${p}_{t}$ denote the distribution of\n",
"haplotype frequencies in generation $t$. \n",
"\n",
"Then, ignoring recombination,\n",
"mutation, migration and selection, the distribution of haplotype\n",
"frequencies in the next generation are given by\n",
"\n",
"$${p}_{t+1}={Mp}_{t},$$\n",
"\n",
"where $\\mathbf{M}$ is a $k\\times k$ matrix with element $i,j$ containing the\n",
"probability that a random variable from a\n",
"Multinomial$(2N,\\frac{{x}_{j}'}{2N})$\n",
"distribution would be equal to ${x}_{i}'$,\n",
"for $i,j=1,2,\\ldots,k$.\n",
"\n",
"To model recombination, consider a population with frequency for\n",
"haplotype $A_{i}B_{j}$ given by $f_{ij}$. \n",
"\n",
"In gametes produced by this\n",
"population, the probability of a non-recombinant $A_{1}B_{1}$ is\n",
"$(1-r)\\frac{f_{11}}{2N}$. A recombinant $A_{1}B_{1}$ gamete can be\n",
"produced in one of four ways. They and their associated probabilities\n",
"are:\n",
"\n",
"1. alleles $A_{1}$ and $B_{1}$ originate from two different\n",
" $A_{1}B_{1}$ haplotypes with associated probability\n",
" $r\\frac{f_{11}}{2N}\\times\\frac{(f_{11}-1)}{2N-1}$;\n",
"\n",
"2. allele $A_{1}$ originates from an $A_{1}B_{1}$ haplotype and $B_{1}$\n",
" originates from an $A_{2}B_{1}$ with associated probability\n",
" $r\\frac{f_{11}}{2N}\\times\\frac{f_{21}}{2N-1}$;\n",
"\n",
"3. allele $A_{1}$ originates from an $A_{1}B_{2}$ haplotype and $B_{1}$\n",
" originates from an $A_{1}B_{1}$ with associated probability\n",
" $r\\frac{f_{12}}{2N}\\times\\frac{f_{11}}{2N-1}$; and\n",
"\n",
"4. allele $A_{1}$ originates from an $A_{1}B_{2}$ haplotype and $B_{1}$\n",
" originates from an $A_{2}B_{1}$ with associated probability\n",
" $r\\frac{f_{12}}{2N}\\times\\frac{f_{21}}{2N-1}$.\n",
"\n",
"Combining these probabilities gives\n",
"
\n",
"To accommodate recombination and mutation in computing\n",
"the distribution of haplotype frequencies, $\\mathbf{M}$ is modified such that column\n",
"$j$ contains probabilities from the\n",
"Multinomial$(2N,{\\beta}_{j}')$\n",
"distribution.\n",
"\n",
"Starting with an allele frequency of 0.5 at each locus and gametic\n",
"equilibrium between the two loci, the expected value of $\\rho^{2}$ was\n",
"computed for 2000 generations given a mutation rate of $u=v=1e^{-9}$, a\n",
"recombination rate of $r=0.002$ between the loci and an effective\n",
"population size of $N_{e}=5,10,25,$ or 50.\n",
"\n",
"The results are shown in the\n",
"figures [fig:r2N5] through [fig:r2N50]. \n",
"\n",
"In addition to $\\rho^{2}$, the\n",
"figures also plot the frequencies of three groups of populations.\n",
"\n",
"In\n",
"populations that belong to group 1, $\\rho^{2}<1$. \n",
"\n",
"In populations that\n",
"belong to group 2, two of the haplotypes are lost such that $\\rho^{2}=1$\n",
"(for example, when haplotypes $A_{1}B_{2}$ and $A_{2}B_{1}$ are lost and\n",
"only haplotypes $A_{1}B_{1}$ and $A_{2}B_{2}$ are segregating,\n",
"$\\rho=1$). \n",
"\n",
"In populations of group 3, one of the loci is fixed, and\n",
"therefore $\\rho$ is not defined.\n",
"\n",
"In all these cases, group 1 starts out having frequency close to 1.0.\n",
"\n",
"Due to drift, however, the frequency in group 1 drops rapidly and\n",
"frequencies in groups 2 and 3 rise.\n",
"\n",
"The expected value of $\\rho^{2}$\n",
"depends to a large extent on the relative magnitudes of the frequencies\n",
"of groups 1 and 2. \n",
"\n",
"Also, drift seems to reduce the frequency of group 1\n",
"faster than that of group 2. \n",
"\n",
"Therefore, $\\rho^{2}$ rises rapidly and\n",
"then stays high for a period. \n",
"\n",
"Once frequencies in groups 1 and 2 drop\n",
"sufficiently low, changes in group frequencies due to mutation become\n",
"significant. \n",
"\n",
"Mutation in group 3, which has the highest frequency, adds\n",
"to group 1 faster than to group 2. \n",
"\n",
"Further, recombination in group 2\n",
"also contributes to group 1. \n",
"\n",
"The balance between these forces and drift\n",
"back into group 3 from groups 1 and 2 determines the equilibrium value\n",
"for $\\rho^{2}$, which is reached by generation 2000 in all four cases."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"[fig:r2N5]\n",
"\n",
"\n",
"\n",
"[fig:r2N10]\n",
"\n",
"\n",
"\n",
"[fig:r2N25]\n",
"\n",
"\n",
"\n",
"[fig:r2N50]\n",
"\n",
"The relationship between the equilibrium value of $\\rho^{2}$, the\n",
"recombination rate, mutation rate, and effective population size can be\n",
"seen from figures [fig:N5] through [fig:N25], where for comparison\n",
"deterministic formulas by Sved @sved1971 and Hill @Hill75 are also\n",
"plotted. When mutation rate is zero, there is good agreement with Sved’s\n",
"formula! When mutation is present, agreement is better with Hill’s\n",
"formula.\n",
"\n",
"\n",
"\n",
"[fig:N5]\n",
"\n",
"\n",
"\n",
"[fig:N10]\n",
"\n",
"\n",
"\n",
"[fig:N25]\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"# ## Julia Implementation"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"400"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Input parameters\n",
"Ne = 5\t\t\t\t\t\t\t\t #effective population size\n",
"u = 1e-9\t\t\t\t\t\t\t\t#mutation rate\n",
"r = 0.002\t\t\t\t\t\t\t#recombination rate\n",
"nGen = 400\t\t\t\t\t\t#number of generations of random mating"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"# Setting up the matrices\n",
"n = 2*Ne\t\t\t\t\t\t\t\t#there are 2*Ne gamets in the population\n",
"sizeX = integer((n+3)*(n+2)*(n+1)/factorial(3))\n",
"x = zeros(sizeX,4)\n",
"row = 1\n",
"for (i in 0:n)\n",
"\tfor (j in 0:(n-i))\n",
"\t\tfor (k in 0:(n-i-j))\n",
"\t\t\tl = n-i-j-k\n",
" x[row,:] = [i,j,k,l]\n",
" row += 1\n",
" end\n",
" end\n",
"end\n",
"# Possibilities of mutations\t\n",
"a = (1-u)*(1-u)\t\t#no mutation at locus 1 and no mutation at locus 2\n",
"b = (1-u)*u\t\t\t#mutation at one locus \n",
"c = u*u\t\t\t\t#mutation at both loci\n",
"\t\n",
"Mu = [a b b c\n",
" b a c b \n",
" b c a b \n",
" c b b a]\n",
"A = zeros(sizeX,sizeX)\n",
"nMinus1Inv = 1/(n-1)\n",
"rsqr = zeros(sizeX) \n",
"for (i in 1:sizeX)\n",
"\t \n",
" x00 = x[i,1]\t#count for 00 haplotype from x matrix\n",
" x01 = x[i,2]\t#count for 01 haplotype from x matrix\n",
" x10 = x[i,3]\t#count for 10 haplotype from x matrix\n",
" x11 = x[i,4]\t#count for 11 haplotype from x matrix\n",
" p00 = x00/n\t\t#probability for 00 haplotype from x matrix\n",
" p01 = x01/n\t\t#probability for 01 haplotype from x matrix\n",
" p10 = x10/n\t\t#probability for 10 haplotype from x matrix\n",
" p11 = x11/n\t\t#probability for 11 haplotype from x matrix\n",
" \n",
" p1 = p10 + p11\t\t\t#total probability of 1 at locus 1\n",
" p2 = p01 + p11\t\t\t#total probability of 1 at locus 2\n",
" cov = p11 - p1*p2\t\t#covariance between the locci\n",
" rsqr[i] = cov^2/(p1*(1-p1)*p2*(1-p2)) #squared correlation between the locci\n",
" \n",
"# Adding recombinations \n",
" pRecomb = [p00*(1-r) + r*nMinus1Inv*( p00*(x00-1 + x10) + p01*(x00 + x10) ),\n",
" p01*(1-r) + r*nMinus1Inv*( p00*(x01 + x11) + p01*(x01-1 + x11) ),\n",
" p10*(1-r) + r*nMinus1Inv*( p10*(x10-1 + x00) + p11*(x10 + x00) ),\n",
" p11*(1-r) + r*nMinus1Inv*( p10*(x11 + x01) + p11*(x11-1 + x01) )]\n",
"\n",
"#putting together mutations and recombinations\t\t\t\t\n",
" pMut = Mu*pRecomb\n",
"#getting multinomial probabilities for all possible combinations of haplotypes \n",
" for (j in 1:sizeX)\n",
" d = Multinomial(n,pMut) \n",
" A[j,i] = pdf(d,vec(x[j,:]))\n",
" end\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"# only look at locus pairs where r^2 is defined--- alleles segregating at both loci\n",
"sel = rsqr .> 0; "
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false,
"scrolled": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"# initial probabilities for haplotypes\n",
"p = zeros(sizeX)\n",
"d = Multinomial(n,[0.25,0.25,0.25,0.25])\n",
"for (i in 1:sizeX)\n",
" p[i] = pdf(d,vec(x[i,:]))\n",
"end\n",
"res = zeros(nGen,2)\n",
"# compute probabilities of haplotypes for i = 1:nGen \n",
"for (i in 1:nGen)\n",
" pSel = p[sel]\n",
" pSel /= sum(pSel)\n",
" meanRSqr = sum(pSel.*rsqr[sel])\n",
" res[i,:] = [i,meanRSqr] \n",
" p = A*p \n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"