Re: Fixation rates for mutations by genetic drift

From: Malcolm (malcolm_at_55bank.freeserve.co.uk)
Date: 09/25/04 

Date: Sat, 25 Sep 2004 21:18:40 +0000 (UTC)

"phillip smith" <deletethis-phills@ihug.co.nz> wrote
>
> From Kimura the time for a neutral mutation to be fixed is the mutation
> rate.
> If I have ten genes each of which can have a neutral mutation. What is the
> mean time taken for there to be no individuals which are free of
mutation.
> That is the whole population carries at least one of the ten mutations
>

Here's a simulation for you, written in C.

(It is portable and should work on any computer with a C compiler. Post here
if you want me to send you a PC executable compiled by me).

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <assert.h>

#define WILD 0
#define MUTANT 1

typedef struct
{
  char *maternal; /* maternal genes */
  char *paternal; /* paternal genes */
} ORGANISM;

typedef struct
{
  int Nmales; /* number of males */
  int Nfemales; /* number of females */
  ORGANISM *males; /* list of males */
  ORGANISM *females; /* list of females */
  int Ngenes; /* no genes each organism has */
} GENERATION;

GENERATION *generation(int N, int Ngenes);
void killgeneration(GENERATION *g);
void addwild(GENERATION *g);
void nextgen(GENERATION *prev, GENERATION *next);
int countmutants(GENERATION *g);
void breed(ORGANISM *mother, ORGANISM *father, ORGANISM *child, int Ngenes);
int hasmutant(ORGANISM *org, int Ngenes);
int random(int N);

int main(int argc, char **argv)
{
  int N;
  int Ngenes;
  int gen = 0;
  GENERATION *g1;
  GENERATION *g2;
  GENERATION *gtemp;
  int muts;

  srand(time(0));

  printf("Traces fixation of mutant genes through a population\n, by genetic
drift.");
  printf("Each mutant is originally present in a single copy\n");
  printf("The population is diploid and sexual\n");

  printf("Population size (min 4)\n");
  while( scanf("%d", &N) != 1)
    printf("Bad input\n");
  if(N < 4)
  {
    printf("Bad input\n");
 exit(EXIT_FAILURE);
  }

  printf("Number of genes (min 1)\n");
  while(scanf("%d", &Ngenes) != 1)
    printf("Bad input\n");
  if(Ngenes < 1)
  {
    printf("Bad Input\n");
 exit(EXIT_FAILURE);
  }

  printf("N = %d, Ngenes = %d\n", N, Ngenes);

  g1 = generation(N, Ngenes);
  if(!g1)
  {
    printf("Out of memory\n");
 exit(EXIT_FAILURE);
  }

  g2 = generation(N, Ngenes);
  if(!g2)
  {
    killgeneration(g1);
 printf("Out of memory\n");
 exit(EXIT_FAILURE);
  }

  addwild(g1);

  printf("Gen 0 mutants %d\n", countmutants(g1));

  do
  {
 gen++;
    nextgen( g1, g2);
 muts = countmutants(g2);

 printf("Gen %d mutants %d\n", gen, muts);

    gtemp = g1;
 g1 = g2;
 g2 = gtemp;
  }
  while(muts > 0 && muts < N);

  return 0;
}

/*
  Creates the gereration structure
  Params: N - population size
          Ngenes - number of loci
*/
GENERATION *generation(int N, int Ngenes)
{
  GENERATION *answer;
  int i;
  int ii;

  assert(N >= 4);
  assert(Ngenes > 0);

  answer = malloc(sizeof(GENERATION));
  if(!answer)
    return 0;

  answer->Nmales = N/2;
  answer->Nfemales = N - answer->Nmales;

  answer->males = malloc( answer->Nmales * sizeof(ORGANISM));
  if(!answer->males)
  {
    killgeneration(answer);
 return 0;
  }

  for(i=0;i<answer->Nmales;i++)
  {
    answer->males[i].maternal = 0;
 answer->males[i].paternal = 0;
  }

  answer->females = malloc(answer->Nfemales * sizeof(ORGANISM));
  if(!answer->females)
  {
    killgeneration(answer);
 return 0;
  }

  for(i=0;i<answer->Nfemales;i++)
  {
    answer->females[i].maternal = 0;
    answer->females[i].paternal = 0;
  }

  for(i=0;i<answer->Nmales;i++)
  {
    answer->males[i].maternal = malloc(Ngenes);

    if(!answer->males[i].maternal)
 {
   killgeneration(answer);
   return 0;
 }

 for(ii=0;ii<Ngenes;ii++)
      answer->males[i].maternal[ii] = WILD;

 answer->males[i].paternal = malloc(Ngenes);

 if(!answer->males[i].paternal)
 {
   killgeneration(answer);
   return 0;
 }

 for(ii=0;ii<Ngenes;ii++)
   answer->males[i].paternal[ii] = WILD;
  }

  for(i=0;i<answer->Nfemales;i++)
  {
    answer->females[i].maternal = malloc(Ngenes);

 if(!answer->females[i].maternal)
 {
   killgeneration(answer);
   return 0;
 }

 for(ii=0;ii<Ngenes;ii++)
   answer->females[i].maternal[ii] = WILD;

 answer->females[i].paternal = malloc(Ngenes);

 if(!answer->females[i].paternal)
 {
   killgeneration(answer);
   return 0;
 }

 for(ii=0;ii<Ngenes;ii++)
   answer->females[i].paternal[ii] = WILD;
  }

  answer->Ngenes = Ngenes;

  return answer;
}

/*
  GENERATION structure destructor
*/
void killgeneration(GENERATION *g)
{
  int i;

  if(g->males)
  {
    for(i=0;i<g->Nmales;i++)
    {
      free(g->males[i].maternal);
   free(g->males[i].paternal);
    }
  }

  if(g->females)
  {
    for(i=0;i<g->Nfemales;i++)
    {
      free(g->females[i].maternal);
   free(g->females[i].paternal);
    }
  }

  free(g->males);
  free(g->females);

  free(g);
}

/*
  add the mutant alleles to the population
  adds one per gene.
*/
void addwild(GENERATION *g)
{
  int i;
  int target;

  for(i=0;i<g->Ngenes;i++)
  {
    if(random(2) == 0)
 {
   target = random(g->Nmales);
   if(random(2) == 0)
     g->males[target].maternal[i] = MUTANT;
   else
     g->males[target].paternal[i] = MUTANT;
 }
 else
 {
   target = random(g->Nfemales);
   if(random(2) == 0)
     g->females[target].maternal[i] = MUTANT;
   else
     g->females[target].paternal[i] = MUTANT;
 }
  }
}

/*
  run the simulation for one generation
*/
void nextgen(GENERATION *prev, GENERATION *next)
{
  int i;
  int mummy;
  int daddy;

  for(i=0;i<next->Nmales;i++)
  {
    mummy = random(prev->Nfemales);
 daddy = random(prev->Nmales);

 breed(&prev->females[mummy], &prev->males[daddy], &next->males[i],
prev->Ngenes);
  }

  for(i=0;i<next->Nfemales;i++)
  {
    mummy = random(prev->Nfemales);
 daddy = random(prev->Nmales);

 breed(&prev->females[mummy], &prev->males[daddy], &next->females[i],
prev->Ngenes);
  }

}

/*
  return the number of organisms with at least one mutant allele
*/
int countmutants(GENERATION *g)
{
  int answer = 0;
  int i;

  for(i=0;i<g->Nmales;i++)
    if(hasmutant(&g->males[i], g->Ngenes))
   answer++;

  for(i=0;i<g->Nfemales;i++)
    if(hasmutant(&g->females[i], g->Ngenes))
   answer++;

  return answer;

}

/*
  produce random offspring.
  Params: mother - parent 1
          father - parent 2
    child - child to fill
    Ngenes - number of genes
*/
void breed(ORGANISM *mother, ORGANISM *father, ORGANISM *child, int Ngenes)
{
  int i;

  for(i=0;i<Ngenes;i++)
  {
    if(random(2) == 0)
   child->maternal[i] = mother->maternal[i];
 else
   child->maternal[i] = mother->paternal[i];

 if(random(2) == 0)
   child->paternal[i] = father->maternal[i];
 else
   child->paternal[i] = father->paternal[i];
  }
}

/*
  test if the organism has any mutant alleles
  Returns 1 if has mutant, 0 if not.
*/
int hasmutant(ORGANISM *org, int Ngenes)
{
  int i;

  for(i=0;i<Ngenes;i++)
  {
    if(org->maternal[i] == MUTANT)
   return 1;
 if(org->paternal[i] == MUTANT)
   return 1;
  }

  return 0;
}

/*
  random - produce an random number 0 - N-1.
  Notes could be tweaked to make simulation better.
*/
int random(int N)
{
  if(N == 2)
   return rand() < RAND_MAX / 2 ? 0 : 1;

  return rand() % N;
}

Relevant Pages