我正试图在http://rosalind.info/problems/iprb/
给定:三个正整数
k
、m
和n
,表示总体含有k+m+n
生物体:k
个体对a因子,m
为杂合子,n
为纯合隐性。返回:两个随机选择的交配生物的概率会产生一个拥有显性等位基因的个体(因此显示显性表型)。假设任何两个生物体都可以朋友
我的解决方案适用于示例,但不适用于生成的任何问题。经过进一步的研究,我似乎应该找到随机选择任何一个生物体的概率,找到选择第二个生物体的可能性,然后找到配对产生具有显性等位基因的后代的概率。
我的问题是:我下面的代码找到的概率是多少?它能找到所有可能交配中具有显性等位基因的后代的百分比吗?所以,我的代码不是一次随机交配的概率,而是求解如果所有配对都进行了测试,则具有显性等位点的后代的比例?
f = open('rosalind_iprb.txt', 'r')
r = f.read()
s = r.split()
############# k = # homozygotes dominant, m = #heterozygotes, n = # homozygotes recessive
k = float(s[0])
m = float(s[1])
n = float(s[2])
############# Counts for pairing between each group and within groups
k_k = 0
k_m = 0
k_n = 0
m_m = 0
m_n = 0
n_n = 0
##############
if k > 1:
k_k = 1.0 + (k-2) * 2.0
k_m = k * m
k_n = k * n
if m > 1:
m_m = 1.0 + (m-2) * 2.0
m_n = m * n
if n> 1:
n_n = 1.0 + (n-2) * 2.0
#################
dom = k_k + k_m + k_n + 0.75*m_m + 0.5*m_n
total = k_k + k_m + k_n + m_m + m_n + n_n
chance = dom/total
print chance
查看您的代码,我很难弄清楚它应该做什么。我将在这里解决这个问题。
让我们简化措辞。有n1个类型1、n2个类型2和n3个类型3项目。
有多少种方法可以从所有物品中选择一套2号的?(n1+n2+n3)选择2。
每一对项目都有对应于以下六个无序多集之一的项目类型:{1,1}、{2,2}、{3,3}、{1,2}、{1,3}、{2,3}
形式为{i,i}的多重集有多少?ni选择2。
有多少形式为{i,j}的多重集,其中i!=jni*nj。
因此,六个多集的概率如下:
- P({1,1})=[n1选择2]/[(n1+n2+n3)选择2]
- P({2,2})=[n2选择2]/[(n1+n2+n3)选择2]
- P({3,3})=[n3选择2]/[(n1+n2+n3)选择2]
- P({1,2})=[n1*n2]/[(n1+n2+n3)选择2]
- P({1,3})=[n1*n3]/[(n1+n2+n3)选择2]
- P({2,3})=[n2*n3]/[(n1+n2+n3)选择2]
这些总和为1。请注意,对于X>1,[X选择2]仅为[X*(X-1)/2],对于X=0或1,则为0。
回报:两个随机选择的交配生物体产生具有显性等位基因的个体(从而表现出显性表型)的概率。
要回答这个问题,您只需要确定六个多集中的哪一个对应于此事件。由于缺乏遗传学知识来回答这个问题,我将把这个问题留给你。
例如,假设两个父母中的任何一个是1型,就会产生显性等位基因。那么感兴趣的事件是{1,1},{1,2},{1,3},并且事件的概率是P({1,1})+P({1,2})+P({1,3})。
我在这个问题上花了一些时间,所以,为了在python中澄清:
lst = ['2', '2', '2']
k, m, n = map(float, lst)
t = sum(map(float, lst))
# organize a list with allele one * allele two (possibles) * dominant probability
# multiplications by one were ignored
# remember to substract the haplotype from the total when they're the same for the second haplotype choosed
couples = [
k*(k-1), # AA x AA
k*m, # AA x Aa
k*n, # AA x aa
m*k, # Aa x AA
m*(m-1)*0.75, # Aa x Aa
m*n*0.5, # Aa x aa
n*k, # aa x AA
n*m*0.5, # aa x Aa
n*(n-1)*0 # aa x aa
]
# (t-1) indicate that the first haplotype was select
print(round(sum(couples)/t/(t-1), 5))
如果您感兴趣,我刚刚找到了一个解决方案并将其放入C#中。
public double mendel(double k, double m, double n)
{
double prob;
prob = ((k*k - k) + 2*(k*m) + 2*(k*n) + (.75*(m*m - m)) + 2*(.5*m*n))/((k + m + n)*(k + m + n -1));
return prob;
}
我们的参数是k(显性),m(杂),&n(隐性)。首先,我找到了每个可能的繁殖配对选择的概率,以种群的百分比表示。因此,k的第一轮选择看起来像k/(k+m+n),并且在k的第一回合选择之后的k的第二回合选择看起来像(k-1)/(k+m+n)。然后将这两者相乘得到结果。由于有三个已确定的人群,因此有九个可能的结果。
然后我把每个结果乘以它的优势概率——对于任何有k的结果都是100%,对于m&m、 50%用于m&n、 n&m、 并且对于n&n.现在把结果加在一起,你就有了解决方案。
http://rosalind.info/problems/iprb/
以下是我在python中所做的代码:我们不希望后代是完全隐性的,所以我们应该制作概率树,看看事件可能发生的情况和概率。那么我们想要的概率是1-p_reccesive。以下代码的注释部分提供了更多解释。
"""
Let d: dominant, h: hetero, r: recessive
Let a = k+m+n
Let X = the r.v. associated with the first person randomly selected
Let Y = the r.v. associated with the second person randomly selected without replacement
Then:
k = f_d => p(X=d) = k/a => p(Y=d| X=d) = (k-1)/(a-1) ,
p(Y=h| X=d) = (m)/(a-1) ,
p(Y=r| X=d) = (n)/(a-1)
m = f_h => p(X=h) = m/a => p(Y=d| X=h) = (k)/(a-1) ,
p(Y=h| X=h) = (m-1)/(a-1)
p(Y=r| X=h) = (n)/(a-1)
n = f_r => p(X=r) = n/a => p(Y=d| X=r) = (k)/(a-1) ,
p(Y=h| X=r) = (m)/(a-1) ,
p(Y=r| X=r) = (n-1)/(a-1)
Now the joint would be:
| offspring possibilites given X and Y choice
-------------------------------------------------------------------------
X Y | P(X,Y) | d(dominant) h(hetero) r(recessive)
-------------------------------------------------------------------------
d d k/a*(k-1)/(a-1) | 1 0 0
d h k/a*(m)/(a-1) | 1/2 1/2 0
d r k/a*(n)/(a-1) | 0 1 0
|
h d m/a*(k)/(a-1) | 1/2 1/2 0
h h m/a*(m-1)/(a-1) | 1/4 1/2 1/4
h r m/a*(n)/(a-1) | 0 1/2 1/2
|
r d n/a*(k)/(a-1) | 0 0 0
r h n/a*(m)/(a-1) | 0 1/2 1/2
r r n/a*(n-1)/(a-1) | 0 0 1
Here what we don't want is the element in the very last column where the offspring is completely recessive.
so P = 1 - those situations as follow
"""
path = 'rosalind_iprb.txt'
with open(path, 'r') as file:
lines = file.readlines()
k, m, n = [int(i) for i in lines[0].split(' ')]
a = k + m + n
p_recessive = (1/4*m*(m-1) + 1/2*m*n + 1/2*m*n + n*(n-1))/(a*(a-1))
p_wanted = 1 - p_recessive
p_wanted = round(p_wanted, 5)
print(p_wanted)