transition probabilities for 19-way Arabidopsis MAGIC

arab_magic.R

# calculating transition probabilities for 19-way Arabidopsis MAGIC
#
# full diallel both directions (F1)
# 3 generations random mating (F4) - 384 families
# selfed for 6 generations; could be "cousin" lines from an F4 family
# 1026 MAGIC lines; 527 in paper
#
# Kover et al (2009) PLOS Genet 5: e1000551
# https://doi.org/10.1371/journal.pgen.1000551
#
# Scarcelli et al (2007) PNAS 104:16986-16991
# https://doi.org/10.1073/pnas.0708209104



# these are the two-locus probabilities
arab19 <-
    function(r)
{
    lastgen <- c(AA=1/19, AB=0)
    result <- rbind(lastgen)


    for(i in 2:4) {
        lastgen <- c(AA=lastgen[1]*(1-r) + r/19^2,
                     AB=lastgen[2]*(1-r) + r/19^2)

        result <- rbind(result, lastgen)
    }

    rownames(result) <- paste0("f", 1:4)

    result <- rbind(result,
                    ril=c(AA=lastgen[1]/(1+2*r) + 2*r/(1+2*r)/19^2,
                          AB=lastgen[2]/(1+2*r) + 2*r/(1+2*r)/19^2))


    result
}


# the same thing, algebraically
arab19B <-
    function(r)
{
    result <- rbind(f1=c(1/19, 0),
                    f2=c((19-18*r)/19^2,        r/19^2),
                    f3=c( (19-36*r+18*r^2)/19^2,
                         (1 - (19-36*r+18*r^2)/19)/19/18),
                    f4=c( ((1-r)^3*19 + (3*r-3*r^2+r^3))/19^2,
                         (1 - ((1-r)^3*19 + (3*r-3*r^2+r^3))/19)/19/18),
                    ril=c( (19 - 52*r + 54*r^2 - 18*r^3)/(1+2*r)/19^2,
                           r*(90 - 54*r + 18*r^2)/(1+2*r)/19^2/18 ))

    dimnames(result) <- list(c("f1", "f2", "f3", "f4", "ril"),
                             c("AA", "AB"))
    result

}

# meiosis recombination fraction -> breakpoint probability
R <- function(r)
    r*(90-54*r+18*r^2)/19/(1+2*r)

# inverse: breakpoint probability -> meiosis recombination fraction
r <- function(R) {

    A <- 3^(-9/2) * sqrt(2475-304*R) * (18-19*R)/4
    B <- (18-19*R)/12
    C <- -(A+B)^(1/3)

    D <- 18-19*R
    E <- sqrt(2475 - 304*R)*(18-19*R)/4/3^(9/2)
    F <- (18-19*R)/12
    G <- (E+F)^(1/3) * 27

    C + D/G + 1
}