From ebda8d81b59c899f02d345acae965010760c873d Mon Sep 17 00:00:00 2001 From: Karl Broman Date: Wed, 13 Mar 2013 18:16:25 -0400 Subject: [PATCH] Various Maxima and R code for the main calculations --- Calculations/AILprob.R | 113 +++++++ Calculations/DOprob.R | 135 +++++++++ Calculations/Fk_distribution.txt | 10 + Calculations/ReadMe.txt | 35 +++ Calculations/ail_map_expansion.txt | 16 + Calculations/ail_xchr.txt | 59 ++++ Calculations/ail_xchr_unbalanced.txt | 94 ++++++ Calculations/cross_func.R | 404 +++++++++++++++++++++++++ Calculations/diversity_outcross.R | 18 ++ Calculations/diversity_outcross_sims.R | 26 ++ Calculations/do_map_expansion.txt | 50 +++ Calculations/do_xchr.txt | 28 ++ Calculations/dosim_comparisons.R | 73 +++++ Calculations/hs.txt | 29 ++ Calculations/sim_ail.R | 172 +++++++++++ Calculations/sim_ail_cluster.R | 44 +++ Calculations/sims.R | 127 ++++++++ 17 files changed, 1433 insertions(+) create mode 100644 Calculations/AILprob.R create mode 100644 Calculations/DOprob.R create mode 100644 Calculations/Fk_distribution.txt create mode 100644 Calculations/ReadMe.txt create mode 100644 Calculations/ail_map_expansion.txt create mode 100644 Calculations/ail_xchr.txt create mode 100644 Calculations/ail_xchr_unbalanced.txt create mode 100644 Calculations/cross_func.R create mode 100644 Calculations/diversity_outcross.R create mode 100644 Calculations/diversity_outcross_sims.R create mode 100644 Calculations/do_map_expansion.txt create mode 100644 Calculations/do_xchr.txt create mode 100644 Calculations/dosim_comparisons.R create mode 100644 Calculations/hs.txt create mode 100644 Calculations/sim_ail.R create mode 100644 Calculations/sim_ail_cluster.R create mode 100644 Calculations/sims.R diff --git a/Calculations/AILprob.R b/Calculations/AILprob.R new file mode 100644 index 0000000..34fb182 --- /dev/null +++ b/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 diff --git a/Calculations/DOprob.R b/Calculations/DOprob.R new file mode 100644 index 0000000..da40802 --- /dev/null +++ b/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 + + diff --git a/Calculations/Fk_distribution.txt b/Calculations/Fk_distribution.txt new file mode 100644 index 0000000..cfc285c --- /dev/null +++ b/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 diff --git a/Calculations/ReadMe.txt b/Calculations/ReadMe.txt new file mode 100644 index 0000000..2558e79 --- /dev/null +++ b/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 diff --git a/Calculations/ail_map_expansion.txt b/Calculations/ail_map_expansion.txt new file mode 100644 index 0000000..b79ee72 --- /dev/null +++ b/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)); diff --git a/Calculations/ail_xchr.txt b/Calculations/ail_xchr.txt new file mode 100644 index 0000000..f440535 --- /dev/null +++ b/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)$ diff --git a/Calculations/ail_xchr_unbalanced.txt b/Calculations/ail_xchr_unbalanced.txt new file mode 100644 index 0000000..ae5bcdb --- /dev/null +++ b/Calculations/ail_xchr_unbalanced.txt @@ -0,0 +1,94 @@ +fAfem : radcan((2/3) + (1/3)*(-1/2)^k)$ +fAmal : radcan((2/3) + (1/3)*(-1/2)^(k-1))$ + +f1 : 1/2$ +m1 : 1$ + +f2 : (1-r)/4 + 1/2$ +m2 : (1-r)/2$ +gf : diff((ev(fAfem,k=2) - f2)*2, r)$ +mef2 : radcan(ev(gf, r=0)); +gm : diff((ev(fAmal,k=2) - m2)*2, r)$ +mem2 : radcan(ev(gm, r=0)); +me2 : radcan((2/3)*mef2+(1/3)*mem2); + +m3 : radcan((1-r)*f2 + r*ev(fAfem, k=1)*ev(fAmal, k=1))$ +f3 : radcan(m2/2 + (1-r)/2*f2 + r/2*ev(fAfem, k=1)*ev(fAmal, k=1))$ +gf : diff((ev(fAfem,k=3) - f3)*2, r)$ +mef3 : radcan(ev(gf, r=0)); +gm : diff((ev(fAmal,k=3) - m3)*2, r)$ +mem3 : radcan(ev(gm, r=0)); +me3 : radcan((2/3)*mef3+(1/3)*mem3); + +m4 : radcan((1-r)*f3 + r*ev(fAfem, k=2)*ev(fAmal, k=2))$ +f4 : radcan(m3/2 + (1-r)/2*f3 + r/2*ev(fAfem, k=2)*ev(fAmal, k=2))$ +gf : diff((ev(fAfem,k=4) - f4)*2, r)$ +mef4 : radcan(ev(gf, r=0)); +gm : diff((ev(fAmal,k=4) - m4)*2, r)$ +mem4 : radcan(ev(gm, r=0)); +me4 : radcan((2/3)*mef4+(1/3)*mem4); + +m5 : radcan((1-r)*f4 + r*ev(fAfem, k=3)*ev(fAmal, k=3))$ +f5 : radcan(m4/2 + (1-r)/2*f4 + r/2*ev(fAfem, k=3)*ev(fAmal, k=3))$ +gf : diff((ev(fAfem,k=5) - f5)*2, r)$ +mef5 : radcan(ev(gf, r=0)); +gm : diff((ev(fAmal,k=5) - m5)*2, r)$ +mem5 : radcan(ev(gm, r=0)); +me5 : radcan((2/3)*mef5+(1/3)*mem5); + +m6 : radcan((1-r)*f5 + r*ev(fAfem, k=4)*ev(fAmal, k=4))$ +f6 : radcan(m5/2 + (1-r)/2*f5 + r/2*ev(fAfem, k=4)*ev(fAmal, k=4))$ +gf : diff((ev(fAfem,k=6) - f6)*2, r)$ +mef6 : radcan(ev(gf, r=0)); +gm : diff((ev(fAmal,k=6) - m6)*2, r)$ +mem6 : radcan(ev(gm, r=0)); +me6 : radcan((2/3)*mef6+(1/3)*mem6); + +m7 : radcan((1-r)*f6 + r*ev(fAfem, k=5)*ev(fAmal, k=5))$ +f7 : radcan(m6/2 + (1-r)/2*f6 + r/2*ev(fAfem, k=5)*ev(fAmal, k=5))$ +gf : diff((ev(fAfem,k=7) - f7)*2, r)$ +mef7 : radcan(ev(gf, r=0)); +gm : diff((ev(fAmal,k=7) - m7)*2, r)$ +mem7 : radcan(ev(gm, r=0)); +me7 : radcan((2/3)*mef7+(1/3)*mem7); + +m8 : radcan((1-r)*f7 + r*ev(fAfem, k=6)*ev(fAmal, k=6))$ +f8 : radcan(m7/2 + (1-r)/2*f7 + r/2*ev(fAfem, k=6)*ev(fAmal, k=6))$ +gf : diff((ev(fAfem,k=8) - f8)*2, r)$ +mef8 : radcan(ev(gf, r=0))*1.0; +gm : diff((ev(fAmal,k=8) - m8)*2, r)$ +mem8 : radcan(ev(gm, r=0))*1.0; +me8 : radcan((2/3)*mef8+(1/3)*mem8)*1.0; + +m9 : radcan((1-r)*f8 + r*ev(fAfem, k=7)*ev(fAmal, k=7))$ +f9 : radcan(m8/2 + (1-r)/2*f8 + r/2*ev(fAfem, k=7)*ev(fAmal, k=7))$ +gf : diff((ev(fAfem,k=9) - f9)*2, r)$ +mef9 : radcan(ev(gf, r=0))*1.0; +gm : diff((ev(fAmal,k=9) - m9)*2, r)$ +mem9 : radcan(ev(gm, r=0))*1.0; +me9 : radcan((2/3)*mef9+(1/3)*mem9)*1.0; + +m10 : radcan((1-r)*f9 + r*ev(fAfem, k=8)*ev(fAmal, k=8))$ +f10 : radcan(m9/2 + (1-r)/2*f9 + r/2*ev(fAfem, k=8)*ev(fAmal, k=8))$ +gf : diff((ev(fAfem,k=10) - f10)*2, r)$ +mef10 : radcan(ev(gf, r=0))*1.0; +gm : diff((ev(fAmal,k=10) - m10)*2, r)$ +mem10 : radcan(ev(gm, r=0))*1.0; +me10 : radcan((2/3)*mef10+(1/3)*mem10)*1.0; + +m11 : radcan((1-r)*f10 + r*ev(fAfem, k=9)*ev(fAmal, k=9))$ +f11 : radcan(m10/2 + (1-r)/2*f10 + r/2*ev(fAfem, k=9)*ev(fAmal, k=9))$ +gf : diff((ev(fAfem,k=11) - f11)*2, r)$ +mef11 : radcan(ev(gf, r=0))*1.0; +gm : diff((ev(fAmal,k=11) - m11)*2, r)$ +mem11 : radcan(ev(gm, r=0))*1.0; +me11 : radcan((2/3)*mef11+(1/3)*mem11)*1.0; + +m12 : radcan((1-r)*f11 + r*ev(fAfem, k=10)*ev(fAmal, k=10))$ +f12 : radcan(m11/2 + (1-r)/2*f11 + r/2*ev(fAfem, k=10)*ev(fAmal, k=10))$ +gf : diff((ev(fAfem,k=12) - f12)*2, r)$ +mef12 : radcan(ev(gf, r=0))*1.0; +gm : diff((ev(fAmal,k=12) - m12)*2, r)$ +mem12 : radcan(ev(gm, r=0))*1.0; +me12 : radcan((2/3)*mef12+(1/3)*mem12)*1.0; + diff --git a/Calculations/cross_func.R b/Calculations/cross_func.R new file mode 100644 index 0000000..3f9c90c --- /dev/null +++ b/Calculations/cross_func.R @@ -0,0 +1,404 @@ +###################################################################### +# +# func.R +# +# Karl W Broman, Johns Hopkins University +# 18 Nov 2003 +# +# functions for simulating meiotic products and RI lines +# +###################################################################### + +############################## +# meiosis.sub <- +# +# Simulate the locations of crossovers on a meiotic product +# via the chi-square model (m=0 corresponds to no interference) +# and potentially with an obligate chiasma. +############################## + +meiosis.sub <- +function(L, m=10, obligate.chiasma=TRUE) +{ + if(obligate.chiasma) { # adjust mean no. chiasmata + if(L <= 50) stop("L must be > 50 cM") + if(m==0) f <- function(Lstar,f.L,f.m) f.L-Lstar/(1-exp(-Lstar/50)) + else { + f <- function(Lstar,f.L,f.m=0) + { + lambdastar <- Lstar/50*(f.m+1) + temp <- lambdastar + for(i in 1:length(temp)) + temp[i] <- sum(exp(-lambdastar[i] + (0:f.m)*log(lambdastar[i])- + lgamma((0:f.m)+1)) * (f.m+1-(0:f.m))/(f.m+1)) + f.L - Lstar/(1-temp) + } + } + + Lstar <- uniroot(f,c(1e-5,L+1),f.m=m,f.L=L)$root + } + else Lstar <- L + + if(m==0) { # no interference + if(!obligate.chiasma) # no obligate chiasma + n.xo <- rpois(1,Lstar/100) + else { + up <- qpois(1e-14,Lstar/50,lower=FALSE) + p <- dpois(1:up,Lstar/50)/ppois(0,Lstar/50) + n.chi <- sample(1:up,1,prob=p) + n.xo <- rbinom(1,n.chi,0.5) + } + if(n.xo==0) xo <- NULL + else xo <- sort(runif(n.xo,0,L)) + } + else { # chi-square model + n.chi <- 0 + while(n.chi == 0) { + n.pts <- rpois(1,Lstar/50*(m+1)) + first <- sample(1:(m+1),1) + if(first <= n.pts || !obligate.chiasma) n.chi <- 1 + } + if(first > n.pts) + xo <- NULL + else { + pt.loc <- sort(runif(n.pts,0,L)) + chi.loc <- pt.loc[seq(first,length(pt.loc),by=m+1)] + n.xo <- rbinom(1,length(chi.loc),0.5) + if(n.xo==0) xo <- NULL + else if(length(chi.loc)==1) xo <- chi.loc + else xo <- sort(sample(chi.loc,n.xo,repl=FALSE)) + } + } + + if(length(xo) == 0) xo <- NULL + xo +} + +############################## +# create.par +# +# create a parental individual +############################## + +create.par <- +function(L, allele=1) +{ + if(length(allele) == 1) allele <- rep(allele,2) + if(length(allele) != 2) + stop("allele should be of length 1 or 2") + + list(mat=rbind(c(0,L),allele[1]), + pat=rbind(c(0,L),allele[2])) +} + + +############################## +# meiosis +# +# Output a random meiotic product from an +# input individual. +############################## + +meiosis <- +function(parent, m=10, obligate.chiasma=TRUE) +{ + L <- parent$mat[1,ncol(parent$mat)] + if(abs(parent$pat[1,ncol(parent$pat)] - L) > 1e-13) + stop("There is a problem with the parent's data structure.") + + product <- meiosis.sub(L, m, obligate.chiasma) + a <- sample(1:2,1) +# print(product) +# print(a) + + if(length(product)==0) return(parent[[a]]) + + else { + for(i in 1:length(product)) { +# cat("\t",i,a,"\n") + if(i == 1) + result <- parent[[a]][,parent[[a]][1,]=product[i-1] & parent[[a]][1,]=product[i]] + result <- cbind(result,c(product[i],u[1])) + a <- 3-a + } + temp <- parent[[a]][1,]>=product[length(product)] + result <- cbind(result,parent[[a]][,temp]) + } + + # clean out excess stuff in the result + if(ncol(result)>2) { + keep <- rep(TRUE,ncol(result)) + for(i in 2:(ncol(result)-1)) + if(result[2,i] == result[2,i+1]) + keep[i] <- FALSE + } +# print(result) + result[,keep,drop=FALSE] +} + + +############################## +# cross +# +# cross two individuals to create a +# single progeny +############################## + +cross <- +function(mom, dad, m=10, obligate.chiasma=TRUE, xchr=FALSE, male=FALSE) +{ + if(!xchr) { + return(list(mat=meiosis(mom,m,obligate.chiasma), + pat=meiosis(dad,m,obligate.chiasma))) + } + else { + if(male) + return(list(mat=meiosis(mom,m,obligate.chiasma), + pat=dad$pat)) + else + return(list(mat=meiosis(mom,m,obligate.chiasma), + pat=dad$mat)) + } +} + +ri2 <- +function(L,n.gen=20,m=10,obligate.chiasma=TRUE) +{ + f1 <- create.par(L,c(1,2)) + par1 <- cross(f1,f1,m,obligate.chiasma) + par2 <- cross(f1,f1,m,obligate.chiasma) + for(i in 1:n.gen) { + c1 <- cross(par1,par2,m,obligate.chiasma) + c2 <- cross(par1,par2,m,obligate.chiasma) + par1 <- c1 + par2 <- c2 + } + par1 +} + +ri8 <- +function(L,n.gen=20,m=10,obligate.chiasma=TRUE) +{ + f1a <- create.par(L,c(1,2)) + f1b <- create.par(L,c(3,4)) + f1c <- create.par(L,c(5,6)) + f1d <- create.par(L,c(7,8)) + par1 <- cross(f1a,f1b,m,obligate.chiasma) + par2 <- cross(f1c,f1d,m,obligate.chiasma) + if(length(n.gen)==1) { + for(i in 1:(n.gen+1)) { + c1 <- cross(par1,par2,m,obligate.chiasma) + c2 <- cross(par1,par2,m,obligate.chiasma) + par1 <- c1 + par2 <- c2 + } + return(par1) + } + else { + result <- vector("list",length(n.gen)) + names(result) <- n.gen + n.gen <- c(-1,n.gen) + for(j in 2:length(n.gen)) { + for(i in (n.gen[j-1]+2):(n.gen[j]+1)) { + c1 <- cross(par1,par2,m,obligate.chiasma) + c2 <- cross(par1,par2,m,obligate.chiasma) + par1 <- c1 + par2 <- c2 + } + result[[j-1]] <- par1 + } + return(result) + } + +} + + +############################## +# where.het +# +# find regions of heterozygosity +# in an individual +############################## +where.het <- +function(ind) +{ + if(ncol(ind$mat)==ncol(ind$pat) && all(ind$mat == ind$pat)) { +# cat(" --No regions of heterozygosity\n") + return(NULL) + } + u <- sort(unique(c(ind$mat[1,],ind$pat[1,]))) + het <- NULL + for(i in 2:length(u)) { + mat <- ind$mat[,ind$mat[1,] >= u[i],drop=FALSE] + mat <- mat[2,1] + + pat <- ind$pat[,ind$pat[1,] >= u[i],drop=FALSE] + pat <- pat[2,1] + + if(mat!=pat) { # heterozygous + if(is.null(het)) het <- cbind(u[i-1],u[i]) + else het <- rbind(het,c(u[i-1],u[i])) + } + } + + # clean up + if(nrow(het) > 1) { + keep <- rep(TRUE,nrow(het)) + for(j in 2:nrow(het)) { + if(het[j,1] == het[j-1,2]) { + het[j,1] <- het[j-1,1] + keep[j-1] <- FALSE + } + } + het <- het[keep,,drop=FALSE] + } + het +} + +simDO <- +function(n.pairs=500, k=6, maxs=10, L=200, xchr=FALSE, m=0, obligate.chiasma=FALSE) +{ + if(k < 1) stop("k must be >= 1") + if(n.pairs < 2) stop("n.pairs must be >= 2") + if(maxs < 1) stop("maxs must be >= 1") + if(L < 0) stop("L must be > 0") + + nalle <- 8 + res <- vector("list", 1+maxs) + for(i in 1:(maxs+1)) + res[[i]] <- vector("list", n.pairs) + + for(i in 1:n.pairs) { + # simulate F0 of collaborate cross + crossdir <- sample(1:nalle) + v <- vector("list", nalle/2) + mom <- cross(create.par(L, crossdir[1:2]), + create.par(L, crossdir[3:4]), xchr=xchr, male=FALSE, + m=m, obligate.chiasma=obligate.chiasma) + dad <- cross(create.par(L, crossdir[5:6]), + create.par(L, crossdir[7:8]), xchr=xchr, male=TRUE, + m=m, obligate.chiasma=obligate.chiasma) + + # simulate k generations of inbreeding + for(j in 1:k) { + fkid <- cross(mom, dad, xchr=xchr, male=FALSE, + m=m, obligate.chiasma=obligate.chiasma) + mkid <- cross(mom, dad, xchr=xchr, male=TRUE, + m=m, obligate.chiasma=obligate.chiasma) + mom <- fkid + dad <- mkid + } + res[[1]][[i]] <- list(fkid, mkid) + } + + # simulate DO + for(i in 2:(maxs+1)) { + # randomly swap dads + z <- sample(n.pairs) + while(any(z==1:n.pairs)) + z <- sample(n.pairs) + + # generate two offspring per pair + for(j in 1:n.pairs) + res[[i]][[j]] <- list(cross(res[[i-1]][[j]][[1]], res[[i-1]][[z[j]]][[2]], xchr=xchr, male=FALSE, + m=m, obligate.chiasma=obligate.chiasma), + cross(res[[i-1]][[j]][[1]], res[[i-1]][[z[j]]][[2]], xchr=xchr, male=TRUE, + m=m, obligate.chiasma=obligate.chiasma)) + } + res +} + +convert2discrete <- +function(individual, locations) { + f <- + function(x,y) { + if(!any(y[1,]= 2") + if(n.pairs < 2) stop("n.pairs must be >= 2") + if(L < 0) stop("L must be > 0") + + res <- vector("list", maxk-1) + for(i in 1:(maxk-1)) + res[[i]] <- vector("list", n.pairs) + + # create F1 + mom <- create.par(L, 1:2) + dad <- create.par(L, 1:2) + dad2 <- create.par(L, 2:1) + + # generate F2 + for(i in 1:n.pairs) { + if(!balanced || i <=n.pairs/2) { + fkid <- cross(mom, dad, xchr=xchr, male=FALSE, + m=m, obligate.chiasma=obligate.chiasma) + mkid <- cross(mom, dad, xchr=xchr, male=TRUE, + m=m, obligate.chiasma=obligate.chiasma) + } + else { + fkid <- cross(mom, dad2, xchr=xchr, male=FALSE, + m=m, obligate.chiasma=obligate.chiasma) + mkid <- cross(mom, dad2, xchr=xchr, male=TRUE, + m=m, obligate.chiasma=obligate.chiasma) + } + res[[1]][[i]] <- list(fkid, mkid) + } + if(maxk==2) return(res[[1]]) + + for(i in 2:(maxk-1)) { + # randomly swap dads + z <- sample(n.pairs) + while(any(z==1:n.pairs)) + z <- sample(n.pairs) + + # generate two offspring per pair + for(j in 1:n.pairs) + res[[i]][[j]] <- list(cross(res[[i-1]][[j]][[1]], res[[i-1]][[z[j]]][[2]], xchr=xchr, male=FALSE, + m=m, obligate.chiasma=obligate.chiasma), + cross(res[[i-1]][[j]][[1]], res[[i-1]][[z[j]]][[2]], xchr=xchr, male=TRUE, + m=m, obligate.chiasma=obligate.chiasma)) + } + res +} diff --git a/Calculations/diversity_outcross.R b/Calculations/diversity_outcross.R new file mode 100644 index 0000000..9442689 --- /dev/null +++ b/Calculations/diversity_outcross.R @@ -0,0 +1,18 @@ +A <- cbind(c(0, 1-0.05, 0.05/64), + c(1/2, (1-0.05)/2, 0.05/128), + c(0, 0, 1)) +pi0 <- c(0.03, 0.08, 1) +b <- cbind(c(1, 0, 0), c(0, 1, 0)) +e <- eigen(A) +D <- diag(e$values) +v <- e$vectors +vi <- solve(v) +pi0 %*% v %*% D %*% vi %*% b + +mat <- matrix(ncol=2, nrow=51) +mat[1,] <- pi0[1:2] +for(i in 1:50) + mat[i+1,] <- pi0 %*% v %*% D^i %*% vi %*% b + +plot(0:50, mat[,2], type="l", col="red") +lines(0:50, mat[,1], col="blue") diff --git a/Calculations/diversity_outcross_sims.R b/Calculations/diversity_outcross_sims.R new file mode 100644 index 0000000..88e2ecb --- /dev/null +++ b/Calculations/diversity_outcross_sims.R @@ -0,0 +1,26 @@ +set.seed(27729779+SUB) +source("cross_func.R") + +doA <- simDO(n.pairs=40000, k=SUB) # 22 sec +doX <- simDO(n.pairs=40000, xchr=TRUE, k=SUB) # 10 sec + +r <- seq(0, 0.49, by=0.01) +loc <- imf.h(r) +doA <- convertall2discrete(doA, loc) # 15 sec +doX <- convertall2discrete(doX, loc) # 15 sec + +gA <- vector("list", length(doA)) +for(i in 1:length(doA)) + g[[i]] <- gtable(lapply(doA[[i]], function(a) a[[1]])) + + gtable(lapply(doA[[i]], function(a) a[[2]])) + +gXf <- vector("list", length(doX)) +for(i in 1:length(doX)) + gXf[[i]] <- gtable(lapply(doX[[i]], function(a) a[[1]])) + +gXm <- vector("list", length(doX)) +for(i in 1:length(doX)) + gXm[[i]] <- gtable(lapply(doX[[i]], function(a) a[[1]]), "first") + + +save(gA, gXf, gXm, file="diversity_outcross_resultsSUB.RData") diff --git a/Calculations/do_map_expansion.txt b/Calculations/do_map_expansion.txt new file mode 100644 index 0000000..ce240a8 --- /dev/null +++ b/Calculations/do_map_expansion.txt @@ -0,0 +1,50 @@ + +hapAAauto : (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)$ + +hapAAfemaleX : ( (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$ + +hapAAmaleX : ( (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$ + +hapAAauto : radcan(ev(hapAAauto, k=k+1))$ +hapAAfemaleX : radcan(ev(hapAAfemaleX, k=k+1))$ +hapAAmaleX : radcan(ev(hapAAmaleX, k=k+1))$ + +z : sqrt((1-r)*(9-r))$ +ws : ((1-r+z)/4)^(s-1)$ +ys : ((1-r-z)/4)^(s-1)$ + + +Rauto : radcan(1-8*((1-r)^(s-1)*(hapAAauto- 1/64) + 1/64))$ + +RfemaleX : radcan(1-8*(1/64 + 1/128 * ((128*hapAAmaleX + 64*hapAAfemaleX*(1-r) - (3-r))*(ws-ys)/z - + (1-64*hapAAfemaleX)*(ws+ys))))$ + +RmaleX : radcan(1-8*(1/64 + 1/128 * ((64*hapAAmaleX - 256*hapAAfemaleX + 3) * (ys-ws)*(1-r)/z - + (1-64*hapAAmaleX)*(ws+ys))))$ + +gauto : radcan(diff(Rauto, r))$ +meauto : radcan(ev(gauto, r=0))$ + +meautoalt : 7*(7+s)/8+((7*sqrt(5)-15)/5)*((1-sqrt(5))/4)^(k+1) + -((7*sqrt(5)+15)/5)*((1+sqrt(5))/4)^(k+1)$ + + + +gfemaleX : radcan(diff(RfemaleX, r))$ +mefemaleX : radcan(ev(gfemaleX, r=0))$ + + +gmaleX : radcan(diff(RmaleX, r))$ +memaleX : radcan(ev(gmaleX, r=0))$ + +memaleXalt : (21*s+154)/36 + -((30+14*sqrt(5))/15)*((1+sqrt(5))/4)^(k+1) + -((30-14*sqrt(5))/15)*((1-sqrt(5))/4)^(k+1) + +(-1/2)^s*7/18 -2*(-1/2)^(k+1+s)/3 + -(-1/2)^(s)*(4*sqrt(5)/15)*( ((1+sqrt(5))/4)^(k+1) - ((1-sqrt(5))/4)^(k+1) )$ + + + +radcan((mefemaleX * 2/3 + memaleX * 1/3) / meauto); /* X chr 2/3 that of autosome */ diff --git a/Calculations/do_xchr.txt b/Calculations/do_xchr.txt new file mode 100644 index 0000000..4df1b18 --- /dev/null +++ b/Calculations/do_xchr.txt @@ -0,0 +1,28 @@ +A : matrix([0, 1/2, 0], + [1-r, (1-r)/2, 0], + [r/64, r/128, 1])$ + +pi0 : matrix([m0, f0, 1])$ + +b : transpose(matrix([1, 0, 0], [0, 1, 0]))$ + +ev : eigenvectors(A)$ +D : matrix([ev[1][1][1], 0, 0], + [0, ev[1][1][2], 0], + [0, 0, 1])$ +Dk : matrix([ev[1][1][1]^k, 0, 0], + [0, ev[1][1][2]^k, 0], + [0, 0, 1])$ +v : transpose(matrix(ev[2][1][1], ev[2][2][1], ev[2][3][1]))$ +vinv : invert(v)$ +radcan(A - v . D . vinv); + +z : sqrt((1-r)*(9-r)); +w : ((1-r + z)/4)^k; +y : ((1-r - z)/4)^k; + + +/* use grind() */ +maleX : 1/64 + ((64*m0-256*f0+3)*(y-w)*(1-r)/z - (1-64*m0)*(w+y)) / 128$ + +femaleX : 1/64 + ( (128*m0 + 64*f0*(1-r) - (3-r))*(w-y)/z - (1-64*f0)*(w+y)) / 128$ diff --git a/Calculations/dosim_comparisons.R b/Calculations/dosim_comparisons.R new file mode 100644 index 0000000..8bfa144 --- /dev/null +++ b/Calculations/dosim_comparisons.R @@ -0,0 +1,73 @@ +source("DOprob.R") + +# load data +r <- seq(0, 0.49, by=0.01) +oRA <- array(dim=c(12,11,50)) +dimnames(oRA) <- list(1:12, 0:10, myround(r, 2)) +oRXf <- oRXm <- oRA +for(i in 1:12) { + cat(i,"\n") + attachfile(i, "diversity_outcross_results") + oRA[i,,] <- cbind(0, t(sapply(gA, apply, 3, function(a) 1-sum(diag(a))/sum(a)))) + oRXf[i,,] <- cbind(0, t(sapply(gXf, apply, 3, function(a) 1-sum(diag(a))/sum(a)))) + oRXm[i,,] <- cbind(0, t(sapply(gXm, apply, 3, function(a) 1-sum(diag(a))/sum(a)))) + detach(2) +} + +RA <- RXf <- RXm <- oRA +for(k in 1:12) { + RA[k,1,] <- rikRprob(r, k, "autosome") + RXf[k,1,] <- rikRprob(r, k, "femaleX") + RXm[k,1,] <- rikRprob(r, k, "maleX") + for(s in 1:10) { + RA[k,s+1,] <- doRprob(r, k, s, "autosome") + RXf[k,s+1,] <- doRprob(r, k, s, "femaleX") + RXm[k,s+1,] <- doRprob(r, k, s, "maleX") + } +} + +max((abs(oRA - RA)/sqrt(RA*(1-RA)/160000)), na.rm=TRUE) +max((abs(oRXf - RXf)/sqrt(RXf*(1-RXf)/80000)), na.rm=TRUE) +max((abs(oRXm - RXm)/sqrt(RXm*(1-RXm)/40000)), na.rm=TRUE) +median((abs(oRA - RA)/sqrt(RA*(1-RA)/160000)), na.rm=TRUE) +median((abs(oRXf - RXf)/sqrt(RXf*(1-RXf)/80000)), na.rm=TRUE) +median((abs(oRXm - RXm)/sqrt(RXm*(1-RXm)/40000)), na.rm=TRUE) + +hist(((oRA - RA)/sqrt(RA*(1-RA)/160000)), breaks=150) +hist(((oRXf - RXf)/sqrt(RXf*(1-RXf)/80000)), breaks=150) +hist(((oRXm - RXm)/sqrt(RXm*(1-RXm)/40000)), breaks=150) + +pdf("results.pdf", height=7.5, width=10) +par(las=1) +for(k in 1:12) { + par(mfrow=c(2,2)) + for(s in 0:10) { + plot(r, oRA[k, s+1,], ylim=c(0, 7/8), xlab="r", ylab="Pr(rec't haplotype)", + main=paste("k = ", k, ", s = ", s, sep=""), cex=0.5) + lines(r, RA[k,s+1,]) + points(r, oRXf[k, s+1,], col="red", cex=0.5) + lines(r, RXf[k, s+1,], col="red") + points(r, oRXm[k, s+1,], col="blue", cex=0.5) + lines(r, RXm[k, s+1,], col="blue") + } +} +dev.off() + +pdf("simulations.pdf", height=7.5, width=10) +par(las=1) +par(mfrow=c(2,2)) +for(k in c(5,12)) { + for(s in c(3,10)) { + plot(r, oRA[k, s+1,], ylim=c(0, 7/8), xlab="recombination fraction (at meiosis)", + ylab="Pr(recombinant haplotype)", + main=paste("k = ", k, ", s = ", s, sep=""), cex=0.5, type="n") + abline(lty=3, h=(0:7)/8, col="gray80", v=seq(0, 0.5, by=0.05)) + points(r, oRA[k, s+1,], cex=0.5) + lines(r, RA[k,s+1,]) + points(r, oRXf[k, s+1,], col="red", cex=0.5) + lines(r, RXf[k, s+1,], col="red") + points(r, oRXm[k, s+1,], col="blue", cex=0.5) + lines(r, RXm[k, s+1,], col="blue") + } +} +dev.off() diff --git a/Calculations/hs.txt b/Calculations/hs.txt new file mode 100644 index 0000000..a6e8fa1 --- /dev/null +++ b/Calculations/hs.txt @@ -0,0 +1,29 @@ + +/* autosome */ +hapAA4way : ((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$ + +hapAA8way : radcan((1-r)/2 * hapAA4way)$ + +hapAAhs : radcan((1-r)^(s-1) * (ev(hapAA8way, k=2) - 1/64) + 1/64)$ +hapAAR : radcan(1 - 8*hapAAhs); +g : radcan(diff(hapAAR, r))$ +me : radcan(ev(g, r=0)); +radcan(me - (7*(s+2)+3)/8); + + +/* X chr */ +femaleX : radcan(ev(( (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, k=1))$ + +maleX : radcan(ev(( (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 , k=1))$ diff --git a/Calculations/sim_ail.R b/Calculations/sim_ail.R new file mode 100644 index 0000000..025e0e6 --- /dev/null +++ b/Calculations/sim_ail.R @@ -0,0 +1,172 @@ +ailA <- simAIL(n.pairs=1000) +ailX <- simAIL(n.pairs=1000, xchr=TRUE) + +r <- seq(0, 0.49, by=0.01) +loc <- imf.h(r) +ailA <- convertall2discrete(ailA, loc) +ailX <- convertall2discrete(ailX, loc) + +g.ailA <- vector("list", length(ailA)) +for(i in 1:length(ailA)) + g.ailA[[i]] <- gtable(lapply(ailA[[i]], function(a) a[[1]])) + + gtable(lapply(ailA[[i]], function(a) a[[2]])) + +g.ailXf <- vector("list", length(ailX)) +for(i in 1:length(ailX)) + g.ailXf[[i]] <- gtable(lapply(ailX[[i]], function(a) a[[1]])) + +g.ailXm <- vector("list", length(ailX)) +for(i in 1:length(ailX)) + g.ailXm[[i]] <- gtable(lapply(ailX[[i]], function(a) a[[1]]), "first") + +###################################################################### + +system.time(ailXbal <- simAIL(n.pairs=1000, balanced=TRUE, xchr=TRUE)) +system.time(ailXbal <- convertall2discrete(ailXbal, loc)) + +g.ailXfbal <- vector("list", length(ailXbal)) +for(i in 1:length(ailXbal)) + g.ailXfbal[[i]] <- gtable(lapply(ailXbal[[i]], function(a) a[[1]])) + +g.ailXmbal <- vector("list", length(ailXbal)) +for(i in 1:length(ailXbal)) + g.ailXmbal[[i]] <- gtable(lapply(ailXbal[[i]], function(a) a[[1]]), "first") + +###################################################################### + +R.ailA <- cbind(0, t(sapply(g.ailA, apply, 3, function(a) 1-sum(diag(a))/sum(a)))) +R.ailXfbal <- cbind(0, t(sapply(g.ailXfbal, apply, 3, function(a) 1-sum(diag(a))/sum(a)))) +R.ailXmbal <- cbind(0, t(sapply(g.ailXmbal, apply, 3, function(a) 1-sum(diag(a))/sum(a)))) + +f <- function(r, k) +{ + if(length(r) > 1) return(sapply(r, f, k)) + A <- cbind(c(0, 1-r, r/4), c(1/2, (1-r)/2, r/8), c(0, 0, 1)) + pi0 <- cbind((1-r)/2, 1/4+(1-r)/4, 1) + if(k==2) return(pi0[1:2]) + cur <- pi0 + for(j in 3:k) + cur <- cur %*% A + cur[1:2] +} + + +k <- 4 +plot(r, R.ailA[k-1,], lwd=2) +lines(r, 1-2*(1/4 + (1-2*r)*(1-r)^(k-2)/4), lwd=2) +points(r, R.ailXfbal[k-1,], lwd=2, col="red") +points(r, R.ailXmbal[k-1,], lwd=2, col="blue") +temp <- f(r, k) +lines(r, 1-2*temp[2,], lwd=2, col="red") +lines(r, 1-2*temp[1,], lwd=2, col="blue") + +###################################################################### + +par(ask=TRUE) +for(k in 2:20) { + plot(r, ailRprob(r, k=k, "autosome"), type="l", lwd=2, main=paste("k =", k), + xlab="recombination fraction (at meiosis)", + ylab="Pr(Recombinant haplotype)", las=1) + lines(r, ailRprob(r, k=k, "femaleX"), lwd=2, col="red") + lines(r, ailRprob(r, k=k, "maleX"), lwd=2, col="blue") +} + +###################################################################### +# comparisons to simulations +###################################################################### +rm(g.ailA, g.ailXf, g.ailXfbal, g.ailXm, g.ailXmbal) + +for(i in 13:20) { + attachfile(i, "ail_results") + if(i==13) { + gailA <- g.ailA + gailXf <- g.ailXf + gailXfbal <- g.ailXfbal + gailXm <- g.ailXm + gailXmbal <- g.ailXmbal + } + else { + for(j in seq(along=g.ailA)) { + gailA[[j]] <- gailA[[j]] + g.ailA[[j]] + gailXf[[j]] <- gailXf[[j]] + g.ailXf[[j]] + gailXfbal[[j]] <- gailXfbal[[j]] + g.ailXfbal[[j]] + gailXm[[j]] <- gailXm[[j]] + g.ailXm[[j]] + gailXmbal[[j]] <- gailXmbal[[j]] + g.ailXmbal[[j]] + } + } + detach(2) +} +oR.ailA <- cbind(0, t(sapply(gailA, apply, 3, function(a) 1-sum(diag(a))/sum(a)))) +oR.ailXfbal <- cbind(0, t(sapply(gailXfbal, apply, 3, function(a) 1-sum(diag(a))/sum(a)))) +oR.ailXmbal <- cbind(0, t(sapply(gailXmbal, apply, 3, function(a) 1-sum(diag(a))/sum(a)))) +oR.ailXf <- cbind(0, t(sapply(gailXf, apply, 3, function(a) 1-sum(diag(a))/sum(a)))) +oR.ailXm <- cbind(0, t(sapply(gailXm, apply, 3, function(a) 1-sum(diag(a))/sum(a)))) + +r <- seq(0, 0.49, by=0.01) + +pdf("ail_results.pdf", height=7.5, width=10) +for(k in 2:12) { + plot(r, oR.ailA[k-1,], lwd=2, main=paste("k =", k), + xlab="recombination fraction (at meiosis)", + ylab="Pr(recombinant haplotype)") + lines(c(r, 0.5), ailRprob(c(r, 0.5), k), lwd=2) + points(r, oR.ailXfbal[k-1,], lwd=2, col="red") + lines(c(r, 0.5), ailRprob(c(r, 0.5), k, "femaleX"), lwd=2, col="red") + points(r, oR.ailXmbal[k-1,], lwd=2, col="blue") + lines(c(r, 0.5), ailRprob(c(r, 0.5), k, "maleX"), lwd=2, col="blue") + + points(r, oR.ailXf[k-1,], lwd=2, col="red", pch=4) + points(r, oR.ailXm[k-1,], lwd=2, col="blue", pch=4) +} +abline(h=c(1/2,4/9), lty=2, col="gray") +dev.off() + +###################################################################### +# haplotype probabilities: unbalanced case +###################################################################### +hapAAXf <- t(sapply(gailXf, apply, 3, function(a) a[1,1]/sum(a))) +hapBBXf <- t(sapply(gailXf, apply, 3, function(a) a[2,2]/sum(a))) +hapABXf <- t(sapply(gailXf, apply, 3, function(a) a[1,2]/sum(a))) +hapBAXf <- t(sapply(gailXf, apply, 3, function(a) a[2,1]/sum(a))) + +hapAAXm <- t(sapply(gailXm, apply, 3, function(a) a[1,1]/sum(a))) +hapBBXm <- t(sapply(gailXm, apply, 3, function(a) a[2,2]/sum(a))) +hapABXm <- t(sapply(gailXm, apply, 3, function(a) a[1,2]/sum(a))) +hapBAXm <- t(sapply(gailXm, apply, 3, function(a) a[2,1]/sum(a))) + +# single locus frequencies +plot(2:12, (hapAAXf + hapABXf)[,1], col="red", lwd=2, ylim=c(0.5, 0.75), + xlab="Generation", ylab="Pr(A)", las=1) +lines(2:12, 2/3 + (1/3)*(-1/2)^(2:12), col="red", lwd=2) +points(2:12, (hapAAXm + hapABXm)[,1], col="blue", lwd=2) +lines(2:12, 2/3 + (1/3)*(-1/2)^(2:12-1), col="blue", lwd=2) + +# haplotype frequencies +plot(2:12, hapAAXf[,1], lwd=2, ylim=c(0.25, 0.75), col="red") +lines(2:12, hapAAailub(0.01, 12)[-1,1], col="red", lwd=2,lty=1) +points(2:12, hapAAXm[,1], lwd=2, pch=1, col="blue") +lines(2:12, hapAAailub(0.01, 12)[-1,2], col="blue", lwd=2, lty=1) +points(2:12, hapAAXf[,25], lwd=2, pch=4, col="red") +lines(2:12, hapAAailub(0.25, 12)[-1,1], col="red", lwd=2, lty=2) +points(2:12, hapAAXm[,25], lwd=2, pch=4, col="blue") +lines(2:12, hapAAailub(0.25, 12)[-1,2], col="blue", lwd=2, lty=2) +points(2:12, hapAAXf[,49], lwd=2, pch=0, col="red") +lines(2:12, hapAAailub(0.49, 12)[-1,1], col="red", lwd=2, lty=3) +points(2:12, hapAAXm[,49], lwd=2, pch=0, col="blue") +lines(2:12, hapAAailub(0.49, 12)[-1,2], col="blue", lwd=2, lty=3) + + +###################################################################### +d <- imf.h(0.25) +f1 <- create.par(d, c(1,2)) +n.pairs <- 1000 +f2m <- f2f <- vector("list", n.pairs) +for(i in 1:n.pairs) { + f2f[[i]] <- cross(f1, f1, xchr=TRUE, male=FALSE, obl=FALSE, m=0) + f2m[[i]] <- cross(f1, f1, xchr=TRUE, male=TRUE, obl=FALSE, m=0) +} +f3m <- f3f <- vector("list", n.pairs) +for(i in 1:n.pairs) { + f3f[[i]] <- cross(f2f[[i]], f2m[[i]], xchr=TRUE, male=FALSE, obl=FALSE, m=0) + f3m[[i]] <- cross(f2f[[i]], f2m[[i]], xchr=TRUE, male=TRUE, obl=FALSE, m=0) +} diff --git a/Calculations/sim_ail_cluster.R b/Calculations/sim_ail_cluster.R new file mode 100644 index 0000000..ed735fb --- /dev/null +++ b/Calculations/sim_ail_cluster.R @@ -0,0 +1,44 @@ +source("cross_func.R") +set.seed(18439413+SUB) +n.sim <- 4000 +ailA <- simAIL(n.pairs=n.sim) + +r <- seq(0, 0.49, by=0.01) +loc <- imf.h(r) +ailA <- convertall2discrete(ailA, loc) + +g.ailA <- vector("list", length(ailA)) +for(i in 1:length(ailA)) + g.ailA[[i]] <- gtable(lapply(ailA[[i]], function(a) a[[1]])) + + gtable(lapply(ailA[[i]], function(a) a[[2]])) + +save(g.ailA, file="ail_resultsSUB.RData") + +ailX <- simAIL(n.pairs=n.sim, xchr=TRUE) +ailX <- convertall2discrete(ailX, loc) + +g.ailXf <- vector("list", length(ailX)) +for(i in 1:length(ailX)) + g.ailXf[[i]] <- gtable(lapply(ailX[[i]], function(a) a[[1]])) + +g.ailXm <- vector("list", length(ailX)) +for(i in 1:length(ailX)) + g.ailXm[[i]] <- gtable(lapply(ailX[[i]], function(a) a[[1]]), "first") + +save(g.ailA, g.ailXf, g.ailXm, file="ail_resultsSUB.RData") + +###################################################################### + +ailXbal <- simAIL(n.pairs=n.sim, balanced=TRUE, xchr=TRUE) +system.time(ailXbal <- convertall2discrete(ailXbal, loc)) + +g.ailXfbal <- vector("list", length(ailXbal)) +for(i in 1:length(ailXbal)) + g.ailXfbal[[i]] <- gtable(lapply(ailXbal[[i]], function(a) a[[1]])) + +g.ailXmbal <- vector("list", length(ailXbal)) +for(i in 1:length(ailXbal)) + g.ailXmbal[[i]] <- gtable(lapply(ailXbal[[i]], function(a) a[[1]]), "first") + +save(g.ailA, g.ailXf, g.ailXm, g.ailXfbal, g.ailXmbal, file="ail_resultsSUB.RData") + diff --git a/Calculations/sims.R b/Calculations/sims.R new file mode 100644 index 0000000..c4a2769 --- /dev/null +++ b/Calculations/sims.R @@ -0,0 +1,127 @@ +source("cross_func.R") +source("DOprob.R") + +doA <- simDO(n.pairs=500) # 22 sec +doX <- simDO(n.pairs=500, xchr=TRUE) # 10 sec + +r <- seq(0, 0.49, by=0.01) +loc <- imf.h(r) +doA <- convertall2discrete(doA, loc) # 15 sec +doX <- convertall2discrete(doX, loc) # 15 sec + +gA <- vector("list", length(doA)) +for(i in 1:length(doA)) + g[[i]] <- gtable(lapply(doA[[i]], function(a) a[[1]])) + + gtable(lapply(doA[[i]], function(a) a[[2]])) + +gXf <- vector("list", length(doX)) +for(i in 1:length(doX)) + gXf[[i]] <- gtable(lapply(doX[[i]], function(a) a[[1]])) + +gXm <- vector("list", length(doX)) +for(i in 1:length(doX)) + gXm[[i]] <- gtable(lapply(doX[[i]], function(a) a[[1]]), "first") + + +RA <- cbind(0, t(sapply(gA, apply, 3, function(a) 1-sum(diag(a))/sum(a)))) +RXf <- cbind(0, t(sapply(gXf, apply, 3, function(a) 1-sum(diag(a))/sum(a)))) +RXm <- cbind(0, t(sapply(gXm, apply, 3, function(a) 1-sum(diag(a))/sum(a)))) + + +plot(r, RA[1,], type="l", ylim=c(0,1)) +for(i in 2:length(x)) lines(r, RA[i,]) + +plot(r, RA[1,], type="l", ylim=c(0,1)) +for(i in 2:length(x)) lines(r, RA[i,]) +for(i in 1:length(x)) lines(r, RXf[i,], col="red") +for(i in 1:length(x)) lines(r, RXm[i,], col="blue") + + +###################################################################### + +n.sim <- 1000 +v <- vector("list", n.sim) +for(i in 1:n.sim) + v[[i]] <- ri8(L=200, n.gen=50, m=0, obligate.chiasma=FALSE) + +r <- seq(0, 0.49, by=0.01) +loc <- imf.h(r) +vd <- convertall2discrete(v, loc) + +g <- gtable(vd) +R <- apply(g, 3, function(a) 1-sum(diag(a))/sum(a)) + +fsnp <- simFounderSnps(list(a=loc), n.str="8") +x <- sim.cross(list(a=loc), n.ind=n.sim, type="ri8sib", foun=fsnp) +x <- x$geno[[1]]$truegeno +v <- vector("list", n.sim) +for(i in 1:n.sim) v[[i]] <- rbind(x[i,], x[i,]) +g <- gtable(v) +R <- apply(g, 3, function(a) 1-sum(diag(a))/sum(a)) +plot(r[-1], R) +lines(r, 7*r/(1+6*r)) + +fsnp4 <- simFounderSnps(list(a=loc), n.str="4") +x <- sim.cross(list(a=loc), n.ind=n.sim, type="ri4sib", foun=fsnp4) +x <- x$geno[[1]]$truegeno +v <- vector("list", n.sim) +for(i in 1:n.sim) v[[i]] <- rbind(x[i,], x[i,]) +g <- gtable(v) +R <- apply(g, 3, function(a) 1-sum(diag(a))/sum(a)) +points(r[-1], R, col="green") +lines(r, 6*r/(1+6*r), col="green") + + +n.sim <- 1000 +v <- vector("list", n.sim) +for(i in 1:n.sim) + v[[i]] <- ri8(L=200, n.gen=50, m=0, obligate.chiasma=FALSE) +vd <- convertall2discrete(v, loc) +g <- gtable(vd) +R <- apply(g, 3, function(a) 1-sum(diag(a))/sum(a)) + +doA <- simDO() +doA <- convertall2discrete(doA, loc) +R <- matrix(ncol=length(loc)-1, nrow=length(doA)) +for(i in 1:length(doA)) { + cat(i,"\n") + g <- gtable(lapply(doA[[i]], function(a) a[[1]])) + + gtable(lapply(doA[[i]], function(a) a[[2]])) + R[i,] <- apply(g, 3, function(a) 1-sum(diag(a))/sum(a)) +} + +plot(r[-1], R[1,], type="l", ylim=c(0,1)) +for(i in 2:length(doA)) lines(r[-1], R[i,]) + +source("DOprob.R") + + +doA <- simDO(maxs=10, k=1) +doA <- convertall2discrete(doA, loc) +R <- matrix(ncol=length(loc)-1, nrow=length(doA)) +for(i in 1:length(doA)) { + cat(i,"\n") + g <- gtable(lapply(doA[[i]], function(a) a[[1]])) + + gtable(lapply(doA[[i]], function(a) a[[2]])) + R[i,] <- apply(g, 3, function(a) 1-sum(diag(a))/sum(a)) +} +plot(r[-1], R[1,]) +lines(r, rikRprob(r, 1)) + +par(mfrow=c(3,3)) +for(i in 1:9) { + plot(r[-1], R[i+1,]) + lines(r, doRprob(r, k=1, s=i-1)) +} + +par(mfrow=c(3,3)) +for(i in 1:9) { + plot(r[-1], R[i+1,]) + lines(r, doRprob(r, k=0, s=i+1)) +} + +par(mfrow=c(3,3)) +for(i in 1:9) { + plot(r[-1], R[i+1,]) + lines(r, doRprob(r, k=0, s=i)) +}