如果第一次閱讀本系列文檔請先移步閱讀【ROSALIND】【練Python，學(xué)生信】00 寫在前面 ?謝謝配合~

題目：

酶切位點出現(xiàn)的期望

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..

所給：一個不超過一百萬的正整數(shù)n，一條長度不超過10的DNA序列s，以及一個所含元素數(shù)不超過20的數(shù)組A，A中的元素是介于0和1之間的數(shù)字。

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相同的數(shù)組B，B[i]代表長度為n的DNA序列t中出現(xiàn)s的次數(shù)期望，t的GC含量由A[i]確定。

?

測試數(shù)據(jù)

10

AG

0.25 0.5 0.75

測試輸出

0.422 0.563 0.422

?

生物學(xué)背景

? ? ? ? 在21 確定限制酶酶切位點中，我們已經(jīng)了解了限制酶的定義和特點。本題關(guān)注的是一個基因組上出現(xiàn)特定限制酶位點的概率有多大。

? ? ? ??從概率角度上來說，一段短的序列可以確保在基因組上反復(fù)出現(xiàn)。比如一段長6個堿基的序列，在隨機(jī)序列上平均每隔4^6=4096個堿基就會出現(xiàn)一次。

? ? ? ??把本題把問題又一般化了一步：計算GC含量已確定的序列某個序列出現(xiàn)的概率。

?

數(shù)學(xué)背景

? ? ? ??期望是概率論中的概念，是試驗中每次可能結(jié)果的概率乘以其結(jié)果的總和。舉例來說，你喜愛的球隊要參加三次比賽，每次獲勝的概率分別為0.3、0.8和0.6，那三次比賽總共獲勝的期望是0.3 + 0.8 + 0.6 = 1.7（不過當(dāng)然，球隊永遠(yuǎn)不可能贏1.7場）。

?

Python背景

? ? ? ??本題有一個我沒解決的問題，即四舍五入問題。Python中保留小數(shù)時常用的round()函數(shù)遵循“四舍六入五平分”原則，“五平分”指根據(jù)取舍的位數(shù)前的小數(shù)奇偶判斷是否進(jìn)位。因此示例結(jié)果本應(yīng)得到0.563時，我只能得到0.562。這其中似乎還涉及浮點數(shù)的存儲問題，我對相關(guān)細(xì)節(jié)力有不逮，歡迎熱心讀者補(bǔ)充解決方法。

?

思路

? ? ? ??我把本題看成兩部分。

? ? ? ??第一部分：想明白整個序列t中有多少條長度與s相同的序列。以示例數(shù)據(jù)為例，序列t長為10，酶切位點長為2。0-9可以數(shù)出01、12、23……78、89共9個長度為2的序列，即10-2+1。所以對任意長度的t和s，符合要求的序列數(shù)量為t的長度-n的長度+1。

? ? ? ??第二部分：計算期望。每個序列與s相等的概率在給定條件下相等，可以輕松求出，期望就是把它們連加起來，即用單個序列的概率乘序列總數(shù)即為所求結(jié)果。

?

?

代碼

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: # 按規(guī)定格式輸出
??? f.write(str(i))
??? f.write(' ')

f.close()