The main locus for the blood type of mice is called Ag-B (B). Several alleles are associated to this locus but for some crossovers Mendel’s laws do not seem to hold. A mating AaBb x AaBb, originated a progeny, yielding

Genotype Frequency Probability
AABB 11
AABb 14
AAbb 1
AaBB 10
AaBb 27
Aabb 12
aaBB 3
aaBb 13
aabb 11

Estimate the recombination fraction, , from these data by the EM algorithm.

Step 1

• Read the data and state ao many recombinant gametes are there for each genotype.

``````nAABB<-11  # 0 recombinant gametes
nAABb<-14  # 1 recombinant gamete
nAAbb<-1   # 2 recombinant gametes
nAaBB<-10  # 1 recombinant gamete
nAaBb<-27  # 0 or 2 recombinant gametes
nAabb<-12  # 1 recombinant gamete
naaBB<-3   # 2 recombinant gametes
naaBb<-13  # 1 recombinant gamete
naabb<-11  # 0 recombinant gametes
``````

• Calculate , the number of individuals from 1 recombinant gametes ().

``````n1 <- nAABb + nAaBB + nAabb + naaBb
n1
``````

• Calculate , the number of individuals from 2 recombinant gametes ().

Note that .

``````n2.star <- NULL
n2 <- nAAbb + naaBB + n2.star
``````

• Calculate n, the total number of individuals ().

``````n <- n1 + nAAbb + naaBB + nAABB + nAaBb + naabb
n
``````

Step 2

• Initialize .

``````r <- 0.3
``````

Step 3 - E (Expectation)

• Create function in order to calculate the expected value for ,

: random variable representing the number of individuals from 2 recombinant gametes, among individuals.

with

then, .

``````expected <- function(r)
{
n2.star <- nAaBb*r^2/(r^2+(1-r)^2)
n2.star
}
``````

Step 4 - M (Maximization)

• Create function \texttt{update.theta} in order to update according to:

meaning that the proportion of recombinant gametes is calculated as the total number of recombinant gametes (0, 1 or 2 for each individual) over the total number of gametes for n individuals.

``````update.theta <- function(n2.star)
{
r <- (n1+2*(nAAbb+naaBB+n2.star))/(2*n)
r
}
``````

Step 5 - Iterative procedure

• Compute the cycle.

``````i<-0
er<-1
error<-10^(-5)
while(er>=error)
{
# Step E
n2.star<-expected(r)
# Step M
r.updated<-update.theta(n2.star)
# Stop criteria
er<-abs(r-r.updated)
i<-i+1
r<-r.updated
cat(i,r,"\n")
}
``````

Step 6 - Print the results

``````cat("\nThe final solution, after",i,"iterations, is r* =",r,"\n")