CENH3 lethality probabilities

We have CENH3+/CENH3- heterozygote that is selfed. On some ears we see 3:1 viable:unviable kernels, suggestive of recessive lethal. But when we genotype unviable kernels, some are not CENH3-/CENH3- and when we genotype viable kernles (seedlings) we find that there are CENH3-/CENH3- among those. Suspicion is a second, possibly linked locus that is recessive lethal.

In the below, “n” is negative and “p” is positive, so “nn” = CENH3-/CENH3- and so on. “D” and “L” refer to dead and live kernels. Kernel/plant phenotype is defined as CENH3 genotype and live/dead status.

Check: all sixteen genotypes sum to 1:

R = 0
0.5^2 * (1 - R)^2 * 2 + 0.5^2 * (1 - R)^2 + 0.5^2 * (1 - R) * R * 2 + 2 * R * 
    (1 - R) * 0.5^2 + 0.5^2 * (1 - R) * R * 2 + 0.5^2 * R^2 + 0.5^2 * R^2 + 
    0.5^2 * R^2 * 2 + 0.5^2 * (1 - R)^2 + (1 - R) * 0.5^2 * R * 2

## [1] 1

Check: all six phenotypes sum to 1:

R = 0.37
nnD = 0.5^2 * (1 - R)^2
nnL = 0.5^2 * R * (2 - R)
pnL = (1 - R + R^2)/2
pnD = 0.5 * (1 - R) * R
ppD = 0.5^2 * R^2
ppL = 0.5^2 * (1 - R) * (1 + R)
nnD + nnL + pnD + pnL + ppD + ppL

## [1] 1

Function to report total phenotypic proportions given recombination fraction

cross = function(R) {

    nnD = 0.5^2 * (1 - R)^2
    nnL = 0.5^2 * R * (2 - R)
    pnL = (1 - R + R^2)/2
    pnD = 0.5 * (1 - R) * R
    ppD = 0.5^2 * R^2
    ppL = 0.5^2 * (1 - R) * (1 + R)
    Dtotal = pnD + nnD + ppD
    Ltotal = pnL + nnL + ppL

    return(paste("total: +-D", pnD, "++D", ppD, "--D", nnD, "+-L", pnL, "++L", 
        ppL, "--L", nnL))
}

Function to report phenotypic proportions conditional on being dead or alive given recombination fraction

conditional_cross = function(R) {

    nnD = 0.5^2 * (1 - R)^2
    nnL = 0.5^2 * R * (2 - R)
    pnL = (1 - R + R^2)/2
    pnD = 0.5 * (1 - R) * R
    ppD = 0.5^2 * R^2
    ppL = 0.5^2 * (1 - R) * (1 + R)
    Dtotal = pnD + nnD + ppD
    Ltotal = pnL + nnL + ppL

    return(paste("conditional: +-D", pnD/Dtotal, "++D", ppD/Dtotal, "--D", nnD/Dtotal, 
        "+-L", pnL/Ltotal, "++L", ppL/Ltotal, "--L", nnL/Ltotal))
}

We see e.g. out of 20 live individuals, 4 are CENH3-/CENH3-. That corresponds to R~0.5:

conditional_cross(0.5)

## [1] "conditional: +-D 0.5 ++D 0.25 --D 0.25 +-L 0.5 ++L 0.25 --L 0.25"

For likelihood, this function:

like_live <- function(R, gens) {
    nnL = 0.5^2 * R * (2 - R)
    pnL = (1 - R + R^2)/2
    ppL = 0.5^2 * (1 - R) * (1 + R)
    Ltotal = pnL + nnL + ppL
    return(dmultinom(gens, prob = c(nnL/Ltotal, pnL/Ltotal, ppL/Ltotal)))
}

Sayuri's observed, from 16 live offspring: 4 -/-, 10 +/-, 2 +/+. Thus:

# data
gens = c(4, 10, 2)

# range of possible recombination rates
rvec = c(0:50/100)

# likelihoods
likeLEE = sapply(rvec, function(X) like_live(X, gens))

# plot
plot(likeLEE ~ rvec, type = "b", pch = 16, cex = 0.5, lwd = 2, ylab = "likelihood", 
    xlab = "recombination")