罗莎琳德“孟德尔第一定律”IPRB答案

【问题标题】：Rosalind "Mendel's First Law" IPRB罗莎琳德“孟德尔第一定律”IPRB
【发布时间】：2014-08-04 12:50:31
【问题描述】：

作为即将到来的生物信息学课程的准备，我正在做一些来自 rosalind.info 的作业。我目前被困在作业“Mendel's First Law”中。

我想我可以强行通过这个，但不知何故，我的想法一定太复杂了。我的方法是这样的：

构建一个包含三个级别的概率树。有两种生物交配，生物 A 和生物 B。第一级是，选择生物 A 纯合显性 (k)、杂合 (m) 或纯合隐性 (n) 的概率是多少。似乎以纯合显性为例，由于总共有(k+m+n)个生物，其中k个是纯合显性的，概率是k/(k+m+n)。

然后在这棵树中，假设我们知道 A 被选为什么生物，那么在每一个下面都会出现生物 B 为 k / m / n 的概率。例如，如果生物 A 被选为杂合子 (m)，那么生物 B 也是杂合子的概率为 (m-1)/(k+m+n-1)，因为现在剩下的杂合子少了一个。

这将给出两个级别的概率，并且会涉及大量代码以达到这一点，因为我实际上是在构建一个树结构，并且为每个分支手动编写了该部分的代码。

现在选择生物 A 和 B 后，他们每个人都有两条染色体。这些染色体之一可以随机挑选。所以对于 A 染色体 1 或 2 可以被挑选出来，对于 B 来说是一样的。所以有 4 种不同的选择：挑选 A 的 1，挑选 B 的 1。挑选 A 的 2，挑选 B 的 1。挑选 A 的 1，挑选 B 的 2。挑选A 中的 2 个，B 中的 2 个。这些中的每一个的概率都是 1/4。所以最后这棵树会有这些叶子概率。

然后我会以某种方式通过魔法将所有这些概率相加，看看两种生物产生具有显性等位基因的生物的概率是多少。

我怀疑这项任务是否设计为需要数小时才能完成。我在想什么？

更新：

以最荒谬的蛮力方式解决了这个问题。只需运行数千次模拟交配并找出最终具有显性等位基因的部分，直到有足够的精度通过分配。

import random
k = 26
m = 18
n = 25

trials = 0
dominants = 0

while True:
    s = ['AA'] * k + ['Aa'] * m + ['aa'] * n
    first = random.choice(s)
    s.remove(first)
    second = random.choice(s)
    has_dominant_allele = 'A' in [random.choice(first), random.choice(second)]
    trials += 1
    if has_dominant_allele:
        dominants += 1
    print "%.5f" % (dominants / float(trials))

【问题讨论】：

你应该问biostars.org
iprb 代表什么？

标签： python bioinformatics rosalind

【解决方案1】：

具有显性等位基因的物种是AA 或Aa。

您的总人口（@987654323@ 由 k (hom) 纯合优势生物与 AA、m (het) 杂合优势生物与 Aa 和 @9876543330@ (@987654 @) 带有aa 的纯合隐性生物。它们中的每一个都可以与任何其他的交配。

具有显性等位基因的生物的概率是：

P_dom = n_dominant/n_total or 1 - n_recessive/n_total

对这些组合中的每一个都做 Punnett 方格并不是一个坏主意：

  hom + het

  |  A | a
-----------
A | AA | Aa
a | Aa | aa


  het + rec

  |  a | a
-----------
A | Aa | Aa
a | aa | aa

显然，两种生物的交配会产生四个可能的孩子。 hom + het 产生 4 个具有隐性等位基因的生物体中的 1 个，het + rec 产生 4 个具有隐性等位基因的生物体中的 2 个。

您可能也想对其他组合执行此操作。

由于我们不只是将生物体一对一交配，而是将整个k + m + n 束组合在一起，因此很高兴知道后代的总数和具有特定等位基因的“孩子”的数量。

如果您不介意一点 Python，来自scipy.misc 的comb 在这里可能会有所帮助。在计算中，不要忘记 (a) 您从每个组合中得到 4 子代，并且 (b) 您需要一个因子（来自 Punnett 方格）来确定组合中的隐性（或显性）后代。

更新

    # total population
    pop_total = 4 * comb(hom + het + rec, 2)

    # use PUNNETT squares!

    # dominant organisms         
    dom_total = 4*comb(hom,2) + 4*hom*het + 4*hom*rec + 3*comb(het,2) + 2*het*rec

    # probability for dominant organisms
    phom = dom_total/pop_total
    print phom

    # probability for dominant organisms + 
    # probability for recessive organisms should be 1
    # let's check that:
    rec_total = 4 * comb(rec, 2) + 2*rec*het + comb(het, 2)
    prec = totalrec/totalpop
    print 1 - prec

【讨论】：

我不能让它工作（不能访问'from scipy.misc import comb'）来测试，但是，我认为这可能会导致一个不正确的答案在这个情况 因为这个问题代表了有限的人口，从混合中抽取一个人限制了以后的选择 - 即并非所有哈迪/温伯格条件都得到满足。在样本数据中，我们分别得到 2、2、2 个 AA、Aa、aa，并被要求计算获得显性表型后代的概率。 H/W 会给你 75%。 Rosalind 给出了 78333% 的正确答案，我也对这个问题感到茫然。
包改成：from scipy.special import comb

【解决方案2】：

这更像是一个概率/计数问题，而不是编码问题。首先计算只有隐性特征的后代的概率更容易。如果您在理解任何内容时遇到任何困难，请告诉我。我运行了以下代码，我的输出通过了 rosalind 分级机。

def mendel(x, y, z):
    #calculate the probability of recessive traits only
    total = x+y+z
    twoRecess = (z/total)*((z-1)/(total-1))
    twoHetero = (y/total)*((y-1)/(total-1))
    heteroRecess = (z/total)*(y/(total-1)) + (y/total)*(z/(total-1))

    recessProb = twoRecess + twoHetero*1/4 + heteroRecess*1/2
    print(1-recessProb) # take the complement

#mendel(2, 2, 2)

with open ("rosalind_iprb.txt", "r") as file:
    line =file.readline().split()
    x, y, z= [int(n) for n in line]
    print(x, y, z)
file.close()
print(mendel(x, y, z))

【讨论】：

【解决方案3】：

克劳斯的解决方案大部分是正确的；但是，在计算具有至少一个显性等位基因的组合数量时会出现错误。这部分是不正确的，因为虽然组合 2 个等位基因形成后代时有 4 种可能性，但实际上只执行了一种可能性。因此，克劳斯的解决方案计算的百分比明显高于应有的百分比。

计算具有至少一个显性等位基因的生物体组合数的正确方法如下：

# k  = number of homozygous dominant organisms
# n = number of heterozygous organisms
# m  = number of homozygous recessive organisms

dom_total = comb(k, 2) + k*m + k*n + .5*m*n + .75*comb(m, 2)

# Instead of:
# 4*comb(k,2) + 4*k*n + 4*k*m + 3*comb(n,2) + 2*n*m

上述代码段用于计算主导组合的总数因为它将每个部分乘以它将产生主导组合的百分比（100% 为 1）后代。您可以将每个部分视为每种类型组合（k&k、k&m、k&n、m&n、m&m）的小方块数。

所以整个正确的代码段应该是这样的：

# Import comb (combination operation) from the scipy library 
from scipy.special import comb

def calculateProbability(k, m, n):
    # Calculate total number of organisms in the population:
    totalPop = k + m + n 
    # Calculate the number of combos that could be made (valid or not):
    totalCombos = comb(totalPop, 2)
    # Calculate the number of combos that have a dominant allele therefore are valid:
    validCombos = comb(k, 2) + k*m + k*n + .5*m*n + .75*comb(m, 2)
    probability = validCombos/totalCombos
    return probability

# Example Call: 
calculateProbability(2, 2, 2)
# Example Output: 0.783333333333

【讨论】：

【解决方案4】：

您不需要在 while 循环中运行数千次模拟。您可以运行一次模拟，并根据其结果计算概率。

from itertools import product

k = 2  # AA  homozygous dominant
m = 2  # Aa  heterozygous
n = 2  # aa  homozygous recessive

population = (['AA'] * k) + (['Aa'] * m) + (['aa'] * n)

all_children = []
for parent1 in population:
    # remove selected parent from population.
    chosen = population[:]
    chosen.remove(parent1)
    for parent2 in chosen:
        # get all possible children from 2 parents. Punnet square
        children = product(parent1, parent2)
        all_children.extend([''.join(c) for c in children])

dominants = filter(lambda c: 'A' in c, all_children)
# float for python2
print float(len(list(dominants))) / len(all_children)
# 0.7833333

【讨论】：

【解决方案5】：

我刚刚找到了答案的公式。您有 8 种可能的交配相互作用，可以产生显性后代：

DDxDD, DDxDd, DdxDD, DdxDd, DDxdd, ddxDD, Ddxdd, ddxDd

分别产生以下优势后代的概率：

1.0, 1.0, 1.0, 0.75, 1.0, 1.0, 0.5, 0.5

一开始我觉得DDxdd 和ddxDD 是两个独立的交配事件似乎很奇怪，但如果你仔细想想，它们在概念上会略有不同。 DDxdd 的概率为k/(k+m+n) * n/((k-1)+m+n)，ddxDD 的概率为n/(k+m+n) * k/(k+m+(n-1))。从数学上讲，这些是相同的，但从概率的角度来看，这是两个独立的事件。所以你的总概率是这些不同交配事件的概率之和乘以该交配事件产生显性后代的概率。我不会在这里一步一步地简化它，但它会给你代码：

total_probability = ((k ** 2 - k) + (2 * k * m) + (3 / 4 * (m ** 2 - m)) + (2* k * n) + (m * n)) / (total_pop ** 2 - total_pop)

您所需要做的就是插入您的k、m 和n 值，您就会得到他们要求的概率。

【讨论】：

【解决方案6】：

我怀疑这项任务是否设计为需要数小时才能完成。我在想什么？

我也有同样的问题。在阅读了整个线程之后，我想出了代码。

希望代码本身能解释概率计算：

def get_prob_of_dominant(k, m, n):
    # A - dominant factor
    # a - recessive factor
    # k - amount of organisms with AA factors (homozygous dominant)
    # m - amount of organisms with Aa factors (heterozygous)
    # n - amount of organisms with aa factors (homozygous recessive)
    events = ['AA+Aa', 'AA+aa', 'Aa+aa', 'AA+AA', 'Aa+Aa', 'aa+aa']

    # get the probability of dominant traits (set up Punnett square)
    punnett_probabilities = {
        'AA+Aa': 1,
        'AA+aa': 1,
        'Aa+aa': 1 / 2,
        'AA+AA': 1,
        'Aa+Aa': 3 / 4,
        'aa+aa': 0,
    }
    event_probabilities = {}
    totals = k + m + n

    # Event: AA+Aa -> P(X=k, Y=m) + P(X=m, Y=k):
    P_km = k / totals * m / (totals - 1)
    P_mk = m / totals * k / (totals - 1)
    event_probabilities['AA+Aa'] = P_km + P_mk

    # Event: AA+aa -> P(X=k, Y=n) + P(X=n, Y=k):
    P_kn = k / totals * n / (totals - 1)
    P_nk = n / totals * k / (totals - 1)
    event_probabilities['AA+aa'] = P_kn + P_nk

    # Event: Aa+aa -> P(X=m, Y=n) +P(X=n, Y=m):
    P_mn = m / totals * n / (totals - 1)
    P_nm = n / totals * m / (totals - 1)
    event_probabilities['Aa+aa'] = P_mn + P_nm

    # Event: AA+AA -> P(X=k, Y=k):
    P_kk = k / totals * (k - 1) / (totals - 1)
    event_probabilities['AA+AA'] = P_kk

    # Event: Aa+Aa -> P(X=m, Y=m):
    P_mm = m / totals * (m - 1) / (totals - 1)
    event_probabilities['Aa+Aa'] = P_mm

    # Event: aa+aa -> P(X=n, Y=n) + P(X=n, Y=n) = 0 (will be * 0, so just don't use it)
    event_probabilities['aa+aa'] = 0

    # Total probability is the sum of (prob of dominant factor * prob of the event)
    total_probability = 0
    for event in events:
        total_probability += punnett_probabilities[event] * event_probabilities[event]
    return round(total_probability, 5)

【讨论】：