8. Estimating recombination rate

First we simulate for 10,000 loci and r=0.001

  start.p<-c(0.25,0.25,0.25,0.25)
  recombination.rate<-0.001 #simulating with r=0.001

  trans.p<-GetTransitionProb(recombination.rate)
  child1.chromatid.mom<-SimulateChromatid(10000,recombination.rate) 
  child1.chromatid.dad<-SimulateChromatid(10000,recombination.rate)
  simulated<-paste("mom.",child1.chromatid.mom$parent.path,"/","dad.",child1.chromatid.dad$parent.path,sep="")
  obs<-child1.chromatid.mom$child+child1.chromatid.dad$child

And now we try to estimate using an arbitrary r=0.4 in our transition matrix

  trans.p<-GetTransitionProb(0.4); #we will estimate with r=0.4
  v.2<-Viterbi3 (obs,start.p, trans.p,child1.chromatid.mom,child1.chromatid.dad)
  vitrowmax.2 <- apply(v.2, 1, function(x) which.max(x))
  estimated<-rownames(trans.p)[vitrowmax.2]
  CountRecombEvents(simulated)

## [1] 21

  CountRecombEvents(estimated)

## [1] 26

And our estimated recombination rate could be calculated like this:

  CountRecombEvents(estimated)/(10000*2)

## [1] 0.0013

Now testing different initial recombination rates

  rates<-c(0.5,0.4,0.3,0.3,0.2,0.1,0.01,0.001)
  estimated.rate<-c()
  recomb.count<-c()
  index<-1;
  for(rate in rates){
    trans.p<-GetTransitionProb(rate); 
    v.2<-Viterbi3 (obs,start.p, trans.p,child1.chromatid.mom,child1.chromatid.dad)
    vitrowmax.2 <- apply(v.2, 1, function(x) which.max(x))
    estimated<-rownames(trans.p)[vitrowmax.2]
    recomb.count[index]<-CountRecombEvents(estimated)
    estimated.rate[index]<- recomb.count[index]/(10000*2);
    index<-index+1
  }
  
  rates

## [1] 0.500 0.400 0.300 0.300 0.200 0.100 0.010 0.001

  estimated.rate

## [1] 0.22125 0.00130 0.00130 0.00130 0.00130 0.00130 0.00130 0.00130

  recomb.count

## [1] 4425   26   26   26   26   26   26   26

  plot(rates,estimated.rate,xlab="arbitrary recomb. rate given to transition matrix",ylab="estimated recomb. rate",main="r.rate estimations with different starting r.rate\n on one simulation with 10,000 loci and r=0.001")

Now we simulate with 10,000 loci and r=0.01

  recombination.rate<-0.01 #simulating with r=0.01
  trans.p<-GetTransitionProb(recombination.rate)
  child1.chromatid.mom<-SimulateChromatid(10000,recombination.rate) 
  child1.chromatid.dad<-SimulateChromatid(10000,recombination.rate)
  simulated<-paste("mom.",child1.chromatid.mom$parent.path,"/","dad.",child1.chromatid.dad$parent.path,sep="")
  obs<-child1.chromatid.mom$child+child1.chromatid.dad$child

And now we try to estimate using an arbitrary r=0.4 in our transition matrix

  trans.p<-GetTransitionProb(0.4); #we will estimate with r=0.4
  v.2<-Viterbi3 (obs,start.p, trans.p,child1.chromatid.mom,child1.chromatid.dad)
  vitrowmax.2 <- apply(v.2, 1, function(x) which.max(x))
  estimated<-rownames(trans.p)[vitrowmax.2]
  CountRecombEvents(simulated)

## [1] 207

  CountRecombEvents(estimated)

## [1] 233

And our estimated recombination rate could be calculated like this:

  CountRecombEvents(estimated)/(10000*2)

## [1] 0.01165

Now testing different initial recombination rates

  rates<-c(0.5,0.4,0.3,0.3,0.2,0.1,0.01,0.001)
  estimated.rate<-c()
  recomb.count<-c()
  index<-1;
  for(rate in rates){
    trans.p<-GetTransitionProb(rate); 
    v.2<-Viterbi3 (obs,start.p, trans.p,child1.chromatid.mom,child1.chromatid.dad)
    vitrowmax.2 <- apply(v.2, 1, function(x) which.max(x))
    estimated<-rownames(trans.p)[vitrowmax.2]
    recomb.count[index]<-CountRecombEvents(estimated)
    estimated.rate[index]<- recomb.count[index]/(10000*2);
    index<-index+1
  }
  
  rates

## [1] 0.500 0.400 0.300 0.300 0.200 0.100 0.010 0.001

  estimated.rate

## [1] 0.21910 0.01165 0.01165 0.01165 0.01165 0.01165 0.01165 0.01155

  recomb.count

## [1] 4382  233  233  233  233  233  233  231

  plot(rates,estimated.rate,xlab="arbitrary recomb. rate given to transition matrix",ylab="estimated recomb. rate",main="r.rate estimations with different starting r.rate\n on one simulation with 10,000 loci and r=0.01")

Now we simulate wih 10,000 loci and r=0.05

  recombination.rate<-0.05 #simulating with r=0.05
  trans.p<-GetTransitionProb(recombination.rate)
  child1.chromatid.mom<-SimulateChromatid(10000,recombination.rate) 
  child1.chromatid.dad<-SimulateChromatid(10000,recombination.rate)
  simulated<-paste("mom.",child1.chromatid.mom$parent.path,"/","dad.",child1.chromatid.dad$parent.path,sep="")
  obs<-child1.chromatid.mom$child+child1.chromatid.dad$child

Now testing different initial recombination rates

  rates<-c(0.5,0.4,0.3,0.3,0.2,0.1,0.01,0.001)
  estimated.rate<-c()
  recomb.count<-c()
  index<-1;
  for(rate in rates){
    trans.p<-GetTransitionProb(rate);
    v.2<-Viterbi3 (obs,start.p, trans.p,child1.chromatid.mom,child1.chromatid.dad)
    vitrowmax.2 <- apply(v.2, 1, function(x) which.max(x))
    estimated<-rownames(trans.p)[vitrowmax.2]
    recomb.count[index]<-CountRecombEvents(estimated)
    estimated.rate[index]<- recomb.count[index]/(10000*2);
    index<-index+1
  }
  
  rates

## [1] 0.500 0.400 0.300 0.300 0.200 0.100 0.010 0.001

  estimated.rate

## [1] 0.23760 0.04995 0.04985 0.04985 0.05010 0.05005 0.04980 0.05010

  recomb.count

## [1] 4752  999  997  997 1002 1001  996 1002

  plot(rates,estimated.rate,xlab="arbitrary recomb. rate given to transition matrix",ylab="estimated recomb. rate",main="r.rate estimations with different starting r.rate\n on one simulation with 10,000 loci and r=0.05")

Now we simulate wih 10,000 loci and r=0.1

  recombination.rate<-0.1 #simulating with r=0.1
  trans.p<-GetTransitionProb(recombination.rate)
  child1.chromatid.mom<-SimulateChromatid(10000,recombination.rate) 
  child1.chromatid.dad<-SimulateChromatid(10000,recombination.rate)
  simulated<-paste("mom.",child1.chromatid.mom$parent.path,"/","dad.",child1.chromatid.dad$parent.path,sep="")
  obs<-child1.chromatid.mom$child+child1.chromatid.dad$child

Now testing different initial recombination rates

  rates<-c(0.5,0.4,0.3,0.3,0.2,0.1,0.01,0.001)
  estimated.rate<-c()
  recomb.count<-c()
  index<-1;
  for(rate in rates){
    trans.p<-GetTransitionProb(rate);
    v.2<-Viterbi3 (obs,start.p, trans.p,child1.chromatid.mom,child1.chromatid.dad)
    vitrowmax.2 <- apply(v.2, 1, function(x) which.max(x))
    estimated<-rownames(trans.p)[vitrowmax.2]
    recomb.count[index]<-CountRecombEvents(estimated)
    estimated.rate[index]<- recomb.count[index]/(10000*2);
    index<-index+1
  }
  
  rates

## [1] 0.500 0.400 0.300 0.300 0.200 0.100 0.010 0.001

  estimated.rate

## [1] 0.24170 0.08740 0.08755 0.08755 0.08755 0.08750 0.08740 0.08745

  recomb.count

## [1] 4834 1748 1751 1751 1751 1750 1748 1749

  plot(rates,estimated.rate,xlab="arbitrary recomb. rate given to transition matrix",ylab="estimated recomb. rate",main="r.rate estimations with different starting r.rate\n on one simulation with 10,000 loci and r=0.1")

Author

Luis Avila (lmavila@gmail.com)