【ROSALIND】【练Python，学生信】46 计算酶切位点出现的期望

如果第一次阅读本系列文档请先移步阅读【ROSALIND】【练Python，学生信】00 写在前面  谢谢配合~

题目：

酶切位点出现的期望

Given: A positive integer n (n≤1,000,000), a DNA string s of even length at most 10, and an array A of length at most 20, containing numbers between 0 and 1..

所给：一个不超过一百万的正整数n，一条长度不超过10的DNA序列s，以及一个所含元素数不超过20的数组A，A中的元素是介于0和1之间的数字。

Return: An array B having the same length as A in which B[i] represents the expected number of times that s will appear as a substring of a random DNA string t of length n, where t is formed with GC-content A[i].

需得：长度与A相同的数组B，B[i]代表长度为n的DNA序列t中出现s的次数期望，t的GC含量由A[i]确定。

测试数据

10

AG

0.25 0.5 0.75

测试输出

0.422 0.563 0.422

生物学背景

        在21 确定限制酶酶切位点中，我们已经了解了限制酶的定义和特点。本题关注的是一个基因组上出现特定限制酶位点的概率有多大。

        从概率角度上来说，一段短的序列可以确保在基因组上反复出现。比如一段长6个碱基的序列，在随机序列上平均每隔4^6=4096个碱基就会出现一次。

        把本题把问题又一般化了一步：计算GC含量已确定的序列某个序列出现的概率。

数学背景

        期望是概率论中的概念，是试验中每次可能结果的概率乘以其结果的总和。举例来说，你喜爱的球队要参加三次比赛，每次获胜的概率分别为0.3、0.8和0.6，那三次比赛总共获胜的期望是0.3 + 0.8 + 0.6 = 1.7（不过当然，球队永远不可能赢1.7场）。

Python背景

        本题有一个我没解决的问题，即四舍五入问题。Python中保留小数时常用的round()函数遵循“四舍六入五平分”原则，“五平分”指根据取舍的位数前的小数奇偶判断是否进位。因此示例结果本应得到0.563时，我只能得到0.562。这其中似乎还涉及浮点数的存储问题，我对相关细节力有不逮，欢迎热心读者补充解决方法。

思路

        我把本题看成两部分。

        第一部分：想明白整个序列t中有多少条长度与s相同的序列。以示例数据为例，序列t长为10，酶切位点长为2。0-9可以数出01、12、23……78、89共9个长度为2的序列，即10-2+1。所以对任意长度的t和s，符合要求的序列数量为t的长度-n的长度+1。

        第二部分：计算期望。每个序列与s相等的概率在给定条件下相等，可以轻松求出，期望就是把它们连加起来，即用单个序列的概率乘序列总数即为所求结果。

代码

f = open('input.txt', 'r')
lines = f.readlines()
f.close()

l = lines[0].replace('\n', '')
l = int(l) # 读入是序列的长度
sites = lines[1].replace('\n', '') # 读入酶切位点的序列
siteslen = len(sites) # 计算酶切位点的长度
GC = lines[2].split(' ')
for i in range(0,len(GC)): # 读入GC含量
    GC[i] = float(GC[i])

expect = []
for i in GC:
    G = i/2
    A = 0.5- i/2
    tep = 1
    for j in sites:
        if j == 'G' or j == 'C':
            tep *= G
        if j == 'A' or j == 'T':
            tep *= A
    tep = tep * (l - siteslen + 1) # 期望计算
    expect.append(round(tep, 3))


print(expect)
f = open('output.txt', 'a')
for i in expect: # 按规定格式输出
    f.write(str(i))
    f.write(' ')

f.close()