% Created 2011-08-15 Mon 17:05
\documentclass[11pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{fixltx2e}
\usepackage{graphicx}
\usepackage{longtable}
\usepackage{float}
\usepackage{wrapfig}
\usepackage{soul}
\usepackage{textcomp}
\usepackage{marvosym}
\usepackage{wasysym}
\usepackage{latexsym}
\usepackage{amssymb}
\usepackage{hyperref}
\tolerance=1000
\usepackage{amsmath}
\usepackage[left=2cm,top=2cm,right=3cm,head=2cm,foot=2cm]{geometry}
\newcommand{\Pr}{\text{Pr}}
\newcommand{\pipe}{\arrowvert}
\providecommand{\alert}[1]{\textbf{#1}}

\title{Genetic drift}
%\author{Stephen Eglen}
\date{\today}

\begin{document}

\maketitle


\section*{A simple model of evolution}
\label{sec-1}

  Evolution is the result of changes in the genetic composition of
  populations over time. One of the simplest models of evolution is as
  follows. There is a population of N individuals, among which there
  are two genetic types: red and blue\footnote{Every individual is chacterised by a single type; no sex,
recombination, mutation, selection, etc. }. Here is the initial
  generation of the population (N=10).



\includegraphics[width=30em]{drift-1-gen.png}

  There is no mutation, no selection and no sex; the next generation
  is made up by randomly choosing 10 individuals from the previous
  generation\footnote{All members of the previous generation die as the next
generation is formed. }. A single individual can be chosen more than once,
  or not at all; the number of times an individual is chosen
  corresponds to the number of progeny it has in the next
  generation. Even without mutation or natural selection the
  proportions of red and blue types will change, because different
  individuals will have different numbers of offspring, by chance.

  So the first two generations might look like this.



\includegraphics[width=30em]{drift-2-gen.png}

  This is a form of evolution called ``genetic drift''. It is inevitable,
  although if the population is very large it will have less effect.

  Let X$_t$ be the number of red individuals in generation t, and let p$_t$
  be the proportion of red individuals. So X$_1$ is 6 and p$_1$ is 0.6. To
  choose the number of red individuals in generation 2 we make 10
  choices, each time having probability 6/10 of getting a red
  individual. So X$_2$ is a \emph{binomial} random variable, with 10 trials and
  success probability 0.6. In general, the random process is described
  by the following transition probabilities.



\includegraphics[width=30em]{transprob.png}


  We can simulate the evolution over many generations in R. This code
  simulates the change in frequency in a single population over 100
  generations. We'll make the population larger (N=1000) but still
  start off with 60\% red individuals.


\begin{verbatim}
p <- numeric(ngens)
p[1] <- X1/N
for(g in 2:ngens)
    p[g] <- rbinom(1, size=N, prob=p[g-1]) / N
plot(p, type="l", ylim=c(0,1), xlab="Generation", ylab="Proportion red")
\end{verbatim}




\includegraphics[width=30em]{simpledrift.png}


  But how variable is this process? To answer this we need to repeat
  the simulation many times (i.e. simulate many identical but
  independent populations). We could do that as follows


\begin{verbatim}
drift.slow <- function(N, X1, ngens, nreps) {
    p <- matrix(NA, nrow=ngens, ncol=nreps)
    p[1,] <- X1/N
    for(rep in 1:nreps) {
        for(g in 2:ngens)
            p[g,rep] <- rbinom(1, size=N, prob=p[g-1,rep]) / N
    }
    p
}
\end{verbatim}




  But this is not a good implementation. One should make use of
  ``vectorisation'', which makes the simulation much more efficient when
  there are many replicates\footnote{Note that we can't vectorise the entire simulation because
drift is a Markov process. }. Note the way that rbinom simulates
  all replicates at once, but still one generation at a time.


\begin{verbatim}
drift.faster <- function(N, X1, ngens, nreps) {
    p <- matrix(NA, nrow=ngens, ncol=nreps)
    p[1,] <- X1/N
    for(gen in 2:ngens)
        p[gen,] <- rbinom(n=nreps, size=N, prob=p[gen-1,]) / N
    p
}
\end{verbatim}




  To run the simulation:


\begin{verbatim}
p <- drift.faster(N, X1, ngens, nreps)
matplot(p, type="l", ylim=c(0,1), lty=1)
\end{verbatim}




\includegraphics[width=30em]{repdrift.png}

  And let's quickly see how much of a speed difference the vectorisation
  makes.


\begin{verbatim}
functions <- c(drift.slow=drift.slow, drift.faster=drift.faster)
times <- sapply(functions, function(f) as.numeric(system.time(f(N, X1, ngens, nreps))[1]))
print(times)
cat(sprintf("\nFactor speed-up = %.1f\n", times[1] / times[2]))
\end{verbatim}




\begin{verbatim}
   drift.slow drift.faster 
        1.383        0.030
 
 Factor speed-up = 46.1
\end{verbatim}

\end{document}