Counting Disease Carriers

Once again we are looking at probability, this time regarding the occurrence of alleles in a population following the Hardy-Weinberg principle. For a population in genetic equilibrium for some given alleles, we are given the frequency of homozygous recessive individuals for each one and are asked to return the probabilities of any randomly selected individual having at least one recessive allele.

We can divide the population into three groups:

A - Homozygous recessive
B - Heterozygous
C - Homozygous dominant

We know that A + B + C = 1. The probability of picking individuals who are not homozygous dominants is P = 1 - C = A + B, which is what we are after.

If we denote the probability of a chromosome having a recessive allele q, and the probability of a chromosome having a dominant allele p, then since each individual carries two chromosomes a homozygous recessive can be described as A = q^2, a heterozygous individual is B = 2pq, and a homozygous dominant is C = p^2.

If we put this together we get that P = q^2 + 2pq. We also know that p + q = 1. So p = 1-q, which gives us P=q^2 + 2(1-q)q which in turn can be simplified to P = 2A^0.5-A.

Now we have our formula, so all we need to do is write a program that makes the calculation for all the given alleles. You can find my version here or below.

A = []
with open('rosalind_afrq.txt','r') as f:
    for nr in f.readline().split(' '):
        A.append(float(nr))

for i in A:
    print(round(2*i**0.5-i,3),end=' ')

Expected Number of Restriction Sites

This problem is very similar to "Introduction to Random Strings". As in that problem, we are given a string, s, and an array, A, containing some GC contents. However, in this case we are also given an integer, n, representing the length of a second string, t. t is the random string formed with each given GC content and we are asked to find the probability of finding s as a substring of each t. What we need to realise to solve this problem is that the number of opportunities for finding s in t is equal to n-len(s)+1. To get the overall probability of finding s in t, we simply need to add all the individual probabilities of randomly forming s with the specified GC contents.

Sample Dataset
10
AG
0.25 0.5 0.75

Expected Output
0.422 0.563 0.422

The following is the code I wrote to solve the problem. When I wrote it I became aware that there is a difference in how Python3 and Python2 handles rounding of 0.5. In Python2 it is rounded up, but in Python3 it is rounded down. On an earlier problem I had to run my program using Python2 because the answer I got from Python3 wasn't accepted by Rosalind. In this case however, the sample data set I got yielded the same answer regardless of which version I used (suggesting it no cases of having to round 0.5 occurred in this data set). (Note: if the following code is to be run with Python2 the formating of the output needs to be rewritten for it to work).

data = []
with open('rosalind_eval.txt', 'r') as f:
    for line in f:
        data.append(line.strip('\n'))
n = int(data[0])
s = data[1]
A = [float(x) for x in data[2].split()]

AT, GC = 0, 0
for nt in s:
    if nt == 'A' or nt == 'T':
        AT += 1
    elif nt == 'G' or nt == 'C':
        GC += 1

B = [None]*len(A)
for i, j in enumerate(A):
    P = (((1 - j)/2)**AT)*((j/2)**GC)*(n - len(s)+1)
    B[i] = '%0.3f' % P
print(*B, sep=' ')

Matching Random Motifs

This problem is very similar to “Introduction to Random Strings”. In this case however, we are given a sting s, and are asked to calculate the probability of randomly forming s when generating N nr of stings of equal length of s and a GC content of x.

Sample Dataset
90000 0.6
ATAGCCGA

Sample Output
0.689

To calculate the probability of randomly getting a string that equals s when constructing it with GC content x, we can use the same equation as in introduction to random strings. To go from this probability (lets call it s_prob) to the probability of getting at least one string that is equal to s we use the following equation:

P(at least 1 match of s) = 1 − P(no matches out of N strings) = 1 − [1 - s_prob]^N

The following program makes the above calculations and outputs the answer with three significant figures:

N = 90000
x = 0.6
s = 'ATAGCCGA'
AT = 0
GC = 0
for nt in s:
    if nt == 'A' or nt == 'T':
        AT += 1
    elif nt == 'G' or nt == 'C':
        GC += 1

s_prob = (((1 - x) / 2)**AT) * (((x) / 2)**GC)
prob = 1 - (1 - s_prob)**N
print('%0.3f' % prob)

Introduction to Random Strings

We know that genomes are not just a random collection of A,T,C and G, and that across a genome there are loads of different motifs, some of which are similar in many different species. But since the human genome is huge, we need to account for the possibility that the subsequence we are looking at was created at random. In this problem we look at a simplified way to calculate the probability that a certain subsequence occurs randomly. We are given a substring of at most 80 bp and an array containing the GC-content of up to 20 random strings, and we are asked to find log10 of the probability that the substring will match the random strings exactly, and return the result in an array of equal length of the given array.

Sample dataset:

ACGATACAA

0.129 0.287 0.423 0.476 0.641 0.742 0.783

Expected output:

-5.737 -5.217 -5.263 -5.360 -5.958 -6.628 -7.009

The probability, P, of the subsequence occuring in a sequence of a GC content x, can be simplified and written as follows, where AC is the total nr of A and C in the subsequence and GC is the total number of G and C:

To solve this problem, we need to write a program that extracts the GC contents of the array and counts AC and GC of the subsequence. Then it's just a simple matter of iterating over the list of GC contents to calculate the probabilities using the above equation. Here is my final version:

import math

AT = 0

GC = 0

with open('sampledata.txt', 'r') as f:

    for line in f:

        if line[0] != 'A' and line[0] != 'T' and line[0] != 'G' and line[

                0] != 'C':

            numbers = line.split()

            GC_contents = [float(x) for x in numbers]

        for i in line:

            if i == 'A' or i == 'T':

                AT += 1

            elif i == 'G' or i == 'C':

                GC += 1

probabilities = []

for j in range(len(GC_contents)):

    prob = math.log10((((1 - GC_contents[j]) / 2)**AT) * (GC_contents[j] / 2)

                      **GC)

    probabilities.append('%0.3f' % prob)

print(*probabilities, sep=' ')

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)