Independent Alleles

We're back in probability again. This time we are looking at independent events, in our case two alleles that are inherited independently (Mendel's 2nd law). In the problem we are to start with an organism with the genotype Aa Bb for traits A and B. He mates with an organism of genotype Aa Bb and gets two offspring. These in turn both mate with organisms of genotype Aa Bb and get two offspring each. This continues for k generations, always mating with organisms of genotype Aa Bb and getting two offspring. We are given the task to calculate the probability that at least N organisms in generation k (excluding the mates) have the genotype Aa Bb.

The key to this problem is realizing  that no matter what genotype the organism that mates with someone of genotype Aa Bb has, the probability that the offspring is Aa Bb is always 0.25 (if you want to you can draw the 9 different 4x4 Punnett squares to convince yourself). When you have realized this you just need to figure out how to calculate the probability that N or more of the population in generation k have the correct genotype, which leads us to binomial distribution.

For the final population P, we want to know the probability that N to P organisms have the correct allele. This means that we want the sum of all the separate probabilities [N, N+1, N+2, ..., P]. We can use the formula for the binomial distribution to calculate all the probabilities separately, and then sum them up to get the overall probability that we are after. For a population of P = k², where k is the number of generations, and when N is the least number of the population with genotype Aa Bb we are looking for, we get the following formula for the overall probability (derived from the general formula of the binomial distribution):

Turning this formula into Python code we get:

import math                                                                    

k = 5                                                                          

N = 8                                                                          

P = 2**k                                                                       

probability = 0                                                                

for i in range(N, P + 1):                                                      

    prob = (math.factorial(P) /                                                

            (math.factorial(i) * math.factorial(P - i))) * (0.25**i) * (0.75**(

                P - i))                                                        

    probability += prob                                                        

print(probability)                                                             

k and N are the variables given by Rosalind.

Calculating Expected Offspring

Once again we are back in the beautiful world of probability! This time we are asked to find the expected number of offspring with a dominant phenotype (i.e. having at least one dominant allele), for a given population of couples, in which every couple receives exactly 2 children. The population consists of the following couples:

1. AA-AA
2. AA-Aa
3. AA-aa
4. Aa-Aa
5. Aa-aa
6. aa-aa

The problem mentions the formula for calculating the expected value of a uniform random variable. However, in the case above, the random variable is not uniform since the probability that couple 1 receives a child with a dominant phenotype is not the same as for couple 4, 5 and 6. We must instead take into account the probability of each type of couple when calculating the expected values. The expected values of the couple types can then be multiplied with the amount of couples of each type and finally added together to give the overall expected value for the population. Note that since there is a 0 probability that a couple 6 receives any child with a dominant allele so their expectation vallue also becomes 0 and we can exclude them from the calculations. The following is the code I wrote to perform the calculations. P1-5 are the number of couples in each type, E1-5 is the expected value for the couples (as you can see E6 would have been 2*0).

P1, P2, P3, P4, P5 = 19100, 18788, 17003, 18906, 16994

E1, E2, E3 = 2, 2, 2                                  

E4 = 2 * 0.75                                         

E5 = 2 * 0.5                                          

E = E1 * P1 + E2 * P2 + E3 * P3 + E4 * P4 + E5 * P5   

print(E)                                              

Mendel's First Law

This time we are asked to calculate the probability that the offspring of two people in a given population receives a dominant allele for a trait. The people in the population are either homozygous dominant (k), heterozygous (m), or homozygous recessive (n) for the trait.

To solve this problem I felt that the easiest way would be to derive an equation for the probability that the offspring gets a dominant allele and then make a program that makes the calculation based on this equation.

The equation must take two things in consideration. First, it needs to include the probability that each parents can have a set of alleles that is either k, m or n (that is AA, Aa or aa, where A is dominant and a is recessive). Then it also needs to consider the probability of each possible couple receiving a child with one or two dominant alleles. For example, if either of the parents is homozygous dominant, then the probability that the child will have at least one dominant allele is 1, but if both parents are homozygous recessive the probability is 0.

To set up the equation I started drawing the possible combinations of parents and the probability for each of them. Below is a sketch of this (pop = total population = k + m + n).

To receive the probability of a specific couple being randomly selected you simply multiply the probability for the first selection event with that of the second selection event. For exmple, the probability that both parents are homozygous dominant is (k/pop)((k-1)/(pop-1)).

After this step I also needed to consider the probability that the offspring actually gets the dominant allele. This is quite easily figured out using a Punnett square. Then you just multiply the probability for the couple with the probability of them being randomly selected. To get the overall probability of the offspring of a randomly selected couple having a dominant allele you just add them all together.

After simplifying the equation, this is what I ended up with:

Then it was just a matter of writing a simple program to make the calculation given k, m and n. The following is what I ended up with, including the values of k, m and n that I received from Rosalind.

k = 23                                                        

m = 26                                                        

n = 22                                                        

pop = k + m + n                                               

prob = (4*(k*(k-1)+2*k*m+2*k*n+m*n)+3*m*(m-1))/(4*pop*(pop-1))

print(prob)