Various Maxima and R code for the main calculations

Calculations/AILprob.R

@@ -0,0 +1,113 @@
+######################################################################
+# AILprob.R
+#
+# Karl Broman, 13 July 2011
+#
+# calculations of haplotype probabilities in AIL
+######################################################################
+
+######################################################################
+# hapAAail:
+#
+# Probability of AA haplotype in AIL
+# 
+# For X chr, we assume balance of A and B (that is, start with male
+# F1s being 1/2 from A x B and 1/2 from B x A)
+#           
+######################################################################
+hapAAail <-
+function(r, k, what=c("autosome", "femaleX", "maleX"))
+{
+  what <- match.arg(what)
+  if(any(k < 2)) stop("must have k > 2")
+  if(any(r < 0 | r > 1/2)) stop("r must be between 0 and 1/2")
+
+  if(what != "autosome") {
+    z <- sqrt((9-r)*(1-r))
+    wk <- ((1-r-z)/4)^(k-2)
+    yk <- ((1-r+z)/4)^(k-2)
+  }
+
+  switch(what,
+         "autosome"= 1/4 + (1-2*r)*(1-r)^(k-2)/4,
+
+         "femaleX" = 1/4 + (yk + wk)*(1-r)/8   + (3-6*r+r^2)*(yk-wk)/(8*z),
+
+         "maleX" =   1/4 + (yk + wk)*(1-2*r)/8 + (3-5*r+2*r^2)*(yk-wk)/(8*z))
+}
+
+######################################################################
+# ailRprob:
+#
+# Probability of recombinant haplotype in AIL
+#
+# For X chr, we assume balance of A and B (that is, start with male
+# F1s being 1/2 from A x B and 1/2 from B x A)
+#           
+######################################################################
+ailRprob <-
+function(r, k, what=c("autosome", "femaleX", "maleX"))
+1 - 2*hapAAail(r, k, what)
+
+
+hapAAailXub <-
+function(r, maxk)
+{
+  if(maxk < 2) stop("maxk >= 2")
+  res <- matrix(nrow=maxk, ncol=2)
+  dimnames(res) <- list(1:maxk, c("femaleAA", "maleAA"))
+  res[1,] <- c(1/2, 1)
+  res[2,] <- c((1-r)/4 + 1/2, (1-r)/2)
+
+  freqAfem <- 2/3 + 1/3 * (-1/2)^(1:maxk)
+  freqAmal <- 2/3 + 1/3 * (-1/2)^(1:maxk-1)
+
+  for(k in 3:maxk) {
+    res[k,2] <- res[k-1,1] * (1-r) + r*freqAfem[k-2]*freqAmal[k-2]
+    res[k,1] <- res[k-1,2]/2 + res[k-1,1]* (1-r)/2 + r/2 * freqAfem[k-2]*freqAmal[k-2]
+  }
+  res
+}
+
+MEailXub <-
+function(maxk)
+{  
+  if(maxk < 2) stop("maxk >= 2")
+  derivhap <- hap <- matrix(nrow=maxk, ncol=2)
+  dimnames(derivhap) <- list(1:maxk, c("female", "male"))
+  hap[1,] <- c(1/2, 1)
+  hap[2,] <- c(1/4 + 1/2, 1/2)
+  derivhap[1,] <- c(0, 0)
+  derivhap[2,] <- c(-1/4, -1/2)
+
+  freqAfem <- 2/3 + 1/3 * (-1/2)^(1:maxk)
+  freqAmal <- 2/3 + 1/3 * (-1/2)^(1:maxk-1)
+
+  for(k in 3:maxk) {
+    hap[k,2] <- hap[k-1,1]
+    hap[k,1] <- hap[k-1,2]/2 + hap[k-1,1]/2
+
+    derivhap[k,2] <- -hap[k-1,1] + derivhap[k-1,1] + freqAfem[k-2]*freqAmal[k-2]
+    derivhap[k,1] <- derivhap[k-1,2]/2 - hap[k-1,1]/2 + derivhap[k-1,1]/2 + freqAfem[k-2]*freqAmal[k-2]/2
+  }
+
+  cbind(-2*derivhap, "overall"= -4/3*derivhap[,1] - 2/3*derivhap[,2])
+}
+
+MEailXub2 <-
+function(maxk)
+{  
+  if(maxk < 2) stop("maxk >= 2")
+  me <- rep(NA, maxk)
+  me[1] <- 0
+  me[2] <- 2/3
+
+  freqAfem <- 2/3 + 1/3 * (-1/2)^(1:maxk)
+  freqAmal <- 2/3 + 1/3 * (-1/2)^(1:maxk-1)
+
+  for(k in 3:maxk) 
+    me[k] <- me[k-1] + 4/3*(freqAfem[k-1] - freqAfem[k-2]*freqAmal[k-2])
+  me
+}
+
+# end of AILprob.R

Calculations/DOprob.R

@@ -0,0 +1,135 @@
+######################################################################
+# DOprob.R
+#
+# Karl Broman, 13 July 2011
+#
+# calculations of haplotype probabilities in the diversity outcross
+# population
+######################################################################
+
+######################################################################
+# hapAA8way:
+#
+# Probability of AA haplotype in pre-CC at G2:Fk
+# 
+# This assumes that the order of the 8 strains in the cross is random.
+#           
+######################################################################
+hapAA8way <-
+function(r, k, what=c("autosome", "femaleX", "maleX"))
+{
+  what <- match.arg(what)
+  if(length(r) > 1) return(sapply(r, hapAA8way, k, what))
+  if(r < 0 || r > 1/2) stop("r must be between 0 and 1/2")
+
+  if(what=="autosome" && r==0.5) {
+    if(k==0) return(1/16)
+    else if(k==1) return(1/32)
+    else return(1/64)
+  }
+
+  switch(what,
+         autosome=
+  (1-r)/2*(((1)/(4*(1+6*r))) - (((6*r^2-7*r-3*r*sqrt(4*r^2-12*r+5))/(4*(1+6*r)*sqrt(4*r^2-12*r+5))))*(((1 - 2*r + sqrt(4*r^2-12*r+5))/(4)))^k + (((6*r^2-7*r+3*r*sqrt(4*r^2-12*r+5))/(4*(1+6*r)*sqrt(4*r^2-12*r+5))))*(((1 - 2*r - sqrt(4*r^2-12*r+5))/(4)))^k)
+
+         ,
+
+         femaleX =
+( (2-r) * (((1)/(3*(1+4*r))) + ((1)/(6*(1+r)))*(-((1)/(2)))^k - (((4*r^3 - (4*r^2+3*r)*sqrt(r^2-10*r+5)+3*r^2-5*r)/(4*(4*r^2+5*r+1)*sqrt(r^2-10*r+5))))*(((1 - r + sqrt(r^2-10*r+5))/(4)))^k + (((4*r^3 + (4*r^2+3*r)*sqrt(r^2-10*r+5)+3*r^2-5*r)/(4*(4*r^2+5*r+1)*sqrt(r^2-10*r+5))))*(((1 - r - sqrt(r^2-10*r+5))/(4)))^k) +
+ (1-r) * (((1)/(3*(1+4*r))) - ((1)/(3*(1+r)))*(-((1)/(2)))^k + (((9*r^2 +5*r + r*sqrt(r^2-10*r+5))/(2*(4*r^2+5*r+1)*sqrt(r^2-10*r+5))))*(((1 - r + sqrt(r^2-10*r+5))/(4)))^k - (((9*r^2 +5*r - r*sqrt(r^2-10*r+5))/(2*(4*r^2+5*r+1)*sqrt(r^2-10*r+5))))*(((1 - r - sqrt(r^2-10*r+5))/(4)))^k) ) / 8
+
+         ,
+
+         maleX=
+( (2-r) * (((1)/(3*(1+4*r))) - ((1)/(3*(1+r)))*(-((1)/(2)))^k + (((r^3 - (8*r^3+r^2-3*r)*sqrt(r^2-10*r+5)-10*r^2+5*r)/(2*(4*r^4-35*r^3-29*r^2+15*r+5)))) *(((1 - r + sqrt(r^2-10*r+5))/(4)))^k + (((r^3 + (8*r^3+r^2-3*r)*sqrt(r^2-10*r+5)-10*r^2+5*r)/(2*(4*r^4-35*r^3-29*r^2+15*r+5)))) *(((1 - r - sqrt(r^2-10*r+5))/(4)))^k) +
+ (1-r) * (((1)/(3*(1+4*r))) + ((2)/(3*(1+r)))*(-((1)/(2)))^k + (((2*r^4 + (2*r^3-r^2+r)*sqrt(r^2-10*r+5)-19*r^3+5*r)/(4*r^4-35*r^3-29*r^2+15*r+5))) *(((1 - r + sqrt(r^2-10*r+5))/(4)))^k + (((2*r^4 - (2*r^3-r^2+r)*sqrt(r^2-10*r+5)-19*r^3+5*r)/(4*r^4-35*r^3-29*r^2+15*r+5))) *(((1 - r - sqrt(r^2-10*r+5))/(4)))^k) ) / 8       
+
+         )
+}
+
+######################################################################
+# rikRprob:
+#
+# Probability of recombinant haplotype in pre-CC at G2:Fk
+#
+# Again, assuming that the order of the 8 strains in the initial cross 
+# is random.
+######################################################################
+rikRprob <-
+function(r, k, what=c("autosome", "femaleX", "maleX"))
+{
+  if(k < 0) stop("k must be >= 0")
+  if(any(r < 0 | r > 1/2)) stop("r must be between 0 and 1/2")
+
+  1-8*hapAA8way(r, k, what)
+}
+
+######################################################################
+# doRprob:
+#
+# Probability of recombinant haplotype at generation s in the 
+# diversity outcross, assuming a large initial population all of
+# pre-CC mice at generation k, and with the order of the 8 strains
+# in the initial crosses to form the pre-CC were random and independent 
+#
+# s = generation of DO
+#
+# At s=1 (first DO generation), each individual is cross between
+# two random pre-CC individuals at generation G2:Fk
+######################################################################
+doRprob <-
+function(r, k, s, what=c("autosome", "femaleX", "maleX"))
+{
+  if(any(k < 0)) stop("k must be >= 0")
+  if(s < 1) stop("s must be >= 1")
+  if(any(r < 0 | r > 1/2)) stop("r must be between 0 and 1/2")
+
+  what <- match.arg(what)
+  if(what=="femaleX" || what=="maleX") {
+    z <- sqrt((1-r)*(9-r))
+    ws <- ((1-r+z)/4)^(s-1)
+    ys <- ((1-r-z)/4)^(s-1)
+
+    malehapprob <- hapAA8way(r, k[1]+1, "male")
+    femalehapprob <- hapAA8way(r, k[1]+1, "female")
+    if(length(k) > 1) {
+      for(j in 2:length(k)) {
+        malehapprob <- malehapprob + hapAA8way(r, k[j]+1, "male")
+        femalehapprob <- femalehapprob + hapAA8way(r, k[j]+1, "female")
+      }
+      malehapprob <- malehapprob/length(k)
+      femalehapprob <- femalehapprob/length(k)
+    }
+  }
+  else {
+    happrob <- hapAA8way(r, k[1]+1, what)
+    if(length(k) > 1) {
+      for(j in 2:length(k)) {
+        happrob <- happrob + hapAA8way(r, k[j]+1, what)
+      }
+      happrob <- happrob/length(k)
+    }
+  }
+
+  switch(what,
+         autosome= 1-8*((1-r)^(s-1)*(happrob - 1/64) + 1/64),
+
+         femaleX=  1-8*(1/64 + 1/128 * ((128*malehapprob + 64*femalehapprob*(1-r) - (3-r))*(ws-ys)/z -
+                                        (1-64*femalehapprob)*(ws+ys))),
+
+         maleX =   1-8*(1/64 + 1/128 * ((64*malehapprob - 256*femalehapprob + 3) * (ys-ws)*(1-r)/z -
+                                        (1-64*malehapprob)*(ws+ys)))
+         )
+}
+
+# HS probabilities, Pr(a1|a2) along single haplotype, from Mott et al (2000)
+hsprobMott <- 
+function(r, g)
+{
+  m <- exp(-g*imf.h(r)/100)
+  cbind(m + (1-m)/8, (1-m)/8)
+}
+
+# end of DOprob.R
+
+

Calculations/Fk_distribution.txt

@@ -0,0 +1,10 @@
+gen number
+F4 21  
+F5 64  
+F6 24  
+F7 10  
+F8 5  
+F9 9  
+F10 5  
+F11 3  
+F12 3  

Calculations/ReadMe.txt

@@ -0,0 +1,35 @@
+
+Various files for calculations related to the AIL probability paper
+
+Fk_distribution.txt         Distribution of generation numbers for founders
+                            of the DO population at the Jackson Lab
+                            [from Karen Svenson]
+
+ail_map_expansion.txt       Maxima code for ail map expansion for autosomes
+ail_xchr.txt                Same for the X chr
+ail_xchr_unbalanced.txt     Same for the X chr in the unbalanced case
+
+do_map_expansion.txt        Maxima code related to map expansion in the DO
+do_xchr.txt                 The same, for the X chromosome
+
+hs.txt                      Optima code for calculations related to
+                            heterogeneous stock
+
+AILprob.R                   Calculations of haplotype probabilities for AIL
+DOprob.R                    Calculations of haplotype probabilities for DO
+
+---------
+Simulation code, not really used in paper but makes me more
+comfortable with the theoretical results
+
+cross_func.R                R code for simulating meiosis and RI lines
+
+diversity_outcross.R        Calculations regarding the DO
+diversity_outcross_sims.R   Simulate the DO  [SUB gets replaced by an integer]
+sims.R                      More DO simulations
+
+dosim_comparisons.R         R code to compare the simulation results
+                            to the theoretical results
+
+sim_ail.R                   simulate AIL
+sim_ail_cluster.R           simulate AIL on the cluster

Calculations/ail_map_expansion.txt

@@ -0,0 +1,16 @@
+/* autosome */
+ailAAhap : 1/4*( 1 + (1-2*r)*(1-r)^(k-2))$
+ailR : radcan(1 - 2*ailAAhap)$
+g : radcan(diff(ailR, r))$
+me : radcan(ev(g, r=0));
+
+/* X chr (balanced case) */
+z : sqrt(r-9)*sqrt(r-1)$
+wk : ((1-r-z)/4)^(k-2)$
+yk : ((1-r+z)/4)^(k-2)$
+malkalt : 1/4 + (yk+wk)*(1-2*r)/8 + (3-5*r+2*r^2)*(yk-wk)/(8*z)$
+femkalt : 1/4 + (yk+wk)*(1-r)/8 + (3-6*r+r^2)*(yk-wk)/(z*8)$
+
+ailRX : 2/3*(1-2*femkalt) + 1/3*(1-2*malkalt)$
+g : radcan(diff(ailRX, r))$
+meX : radcan(ev(g, r=0));

Calculations/ail_xchr.txt

@@ -0,0 +1,59 @@
+A : matrix([0, 1/2, 0], [1-r, (1-r)/2, 0], [r/4, r/8, 1])$
+pi0 : matrix([(1-r)/2, 1/4 + (1-r)/4, 1])$
+b1 : matrix([1], [0], [0])$
+b2 : matrix([0], [1], [0])$
+
+e : eigenvectors(A)$
+D : matrix([e[1][1][1], 0, 0], [0, e[1][1][2], 0], [0, 0, 1])$
+Dk : matrix([e[1][1][1]^k, 0, 0], [0, e[1][1][2]^k, 0], [0, 0, 1])$
+v : transpose(matrix(e[2][1][1], e[2][2][1], [0, 0, 1]))$
+vinv : invert(v)$
+
+
+/* male */
+malk : -2^(-2*k-3)*((2*(-r-sqrt(r-9)*sqrt(r-1)+1)^k-2*(-r+sqrt(r-9)*sqrt(r-1)+1)^k)
+            *r^2
+            +sqrt(r-9)*sqrt(r-1)
+                      *((2*(-r+sqrt(r-9)*sqrt(r-1)+1)^k
+                       +2*(-r-sqrt(r-9)*sqrt(r-1)+1)^k)
+                       *r
+                       -(-r+sqrt(r-9)*sqrt(r-1)+1)^k
+                       -(-r-sqrt(r-9)*sqrt(r-1)+1)^k-2^(2*k+1))
+            +(5*(-r+sqrt(r-9)*sqrt(r-1)+1)^k-5*(-r-sqrt(r-9)*sqrt(r-1)+1)^k)*r
+            -3*(-r+sqrt(r-9)*sqrt(r-1)+1)^k+3*(-r-sqrt(r-9)*sqrt(r-1)+1)^k)/(sqrt(r-9)*sqrt(r-1))$
+
+z : sqrt(r-9)*sqrt(r-1)$
+wk : ((1-r-z)/4)^k$
+yk : ((1-r+z)/4)^k$
+w : (1-r-z)$
+y : (1-r+z)$
+malkalt : 1/4 + (yk+wk)*(1-2*r)/8 + (3-5*r+2*r^2)*(yk-wk)/(8*z)$
+
+
+
+/* female */
+femk : -2^(-2*k-3)*(((-r-sqrt(r-9)*sqrt(r-1)+1)^k-(-r+sqrt(r-9)*sqrt(r-1)+1)^k)*r^2
+            +sqrt(r-9)*sqrt(r-1)
+                      *(((-r+sqrt(r-9)*sqrt(r-1)+1)^k
+                       +(-r-sqrt(r-9)*sqrt(r-1)+1)^k)
+                       *r
+                       -(-r+sqrt(r-9)*sqrt(r-1)+1)^k
+                       -(-r-sqrt(r-9)*sqrt(r-1)+1)^k-2^(2*k+1))
+            +(6*(-r+sqrt(r-9)*sqrt(r-1)+1)^k-6*(-r-sqrt(r-9)*sqrt(r-1)+1)^k)*r
+            -3*(-r+sqrt(r-9)*sqrt(r-1)+1)^k+3*(-r-sqrt(r-9)*sqrt(r-1)+1)^k)/(sqrt(r-9)*sqrt(r-1))$
+
+femkalt : 1/4 + (yk+wk)*(1-r)/8 + (3-6*r+r^2)*(yk-wk)/(z*8)$
+
+
+/* unbalanced case */
+A : transpose(matrix([(1-r)/2, 1/2, r/2, 0], [1-r, 0, r, 0], [0, 0, 0, 1], [1/4, 0, 1/4, 1/2]))$
+pi0 : matrix([1/2, 1, 0, 1/2])$
+pi0p : matrix([1/2, 1/2, 0, 1/2])$
+b1 : matrix([1], [0], [0], [0])$
+b2 : matrix([0], [1], [0], [0])$
+
+/* e : eigenvectors(A);  /* <- this chokes */
+
+/* check with the other results */
+malkalt : 1/4 + (yk+wk)*(1-2*r)/8 + (3-5*r+2*r^2)*(yk-wk)/(8*z)$
+femkalt : 1/4 + (yk+wk)*(1-r)/8 + (3-6*r+r^2)*(yk-wk)/(z*8)$