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

題目：

隨機(jī)字符串（Introduction to Random Strings）

Given: A DNA string s of length at most 100 bp and an array A containing at most 20 numbers between 0 and 1.

所給：一段長度不大于100bp的DNA序列s，以及一個(gè)包含20個(gè)以內(nèi)在0和1之間數(shù)字的數(shù)組A。

Return: An array B having the same length as A in which B[k] represents the common logarithm of the probability that a random string constructed with the GC-content found in A[k] will match s exactly.

需得：一個(gè)和A長度相同的數(shù)組B，B[k]是以A[k]為GC含量的一個(gè)序列與s完全相同的概率的常用對數(shù)。

?

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

ACGATACAA

0.129 0.287 0.423 0.476 0.641 0.742 0.783

測試輸出

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

?

生物學(xué)背景

? ? ? ? 基因組DNA并非核苷酸的隨機(jī)排列，存在著很多反復(fù)出現(xiàn)的模序（motif）。如果一段序列在很多物種的基因組中都反復(fù)出現(xiàn)，我們就會(huì)懷疑它可能承擔(dān)著某種生物學(xué)功能。但同時(shí)也要注意，如果一條序列足夠長，理論上它可以包含所有可能的短序列，而人類的基因組包含多于三十億個(gè)堿基對。因此如果一段序列反復(fù)出現(xiàn)，我們需要先確定這種現(xiàn)象并非僅僅由概率導(dǎo)致。

????????理論上來說，越短的序列出現(xiàn)的概率越高，反之，長序列很難隨機(jī)出現(xiàn)。怎樣確定這個(gè)長短的臨界值就成了關(guān)鍵點(diǎn)。為解決這個(gè)問題，本題先從用各堿基出現(xiàn)的概率隨機(jī)構(gòu)建序列入手，并計(jì)算此序列出現(xiàn)的概率。

?

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

? ? ? ? 數(shù)組是有序的元素序列，可以看做僅含有一行的矩陣。A[k]代表數(shù)組A中的第k個(gè)數(shù)。

? ? ? ? 當(dāng)我們假設(shè)一個(gè)DNA序列的GC含量為x時(shí)，G或C出現(xiàn)的頻率即為x/2，A或T出現(xiàn)的頻率為(1-x)/2,。

? ? ? ? 概率相乘時(shí)往往得到很小的數(shù)字，在計(jì)算機(jī)存儲(chǔ)時(shí)導(dǎo)致溢出，取對數(shù)是常用的方法。對任意正數(shù)x和y，有l(wèi)og10(x?y)=log10(x)+log10(y)成立。

?

思路

????????本題思路較簡單，可分為如下兩步依次實(shí)現(xiàn)：

一．讀入題目給出的GC含量，計(jì)算出G、C、A、T的概率；

二．依次讀入序列，將在給定概率下每個(gè)位置出現(xiàn)該堿基的概率的對數(shù)相加。

?

代碼

import? math

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

seq = input[0].replace('\n','')
per = input[1].split(' ')
i = 0
GC = []
AT = []
while i < len(per): # 計(jì)算各堿基的頻率
??? per[i] = float(per[i])
??? GC.append(per[i] / 2)
??? AT.append(0.5 - per[i] / 2)
??? i += 1
i = 0
result = []

while i < len(per): # 計(jì)算序列出現(xiàn)的概率
??? j = 0
??? proGC = GC[i]
??? proAT = AT[i]
??? pro = 0
??? while j < len(seq):
?????? if seq[j] == 'G' or seq[j] == 'C':
?????????? pro = pro + math.log10(proGC)
?????? else:
?????????? pro = pro + math.log10(proAT)
?????? j += 1
??? result.append(round(pro,3 ))
??? i += 1
# print(result)

i = 0
f = open('output.txt', 'a')
while i < len(per):
??? f.write(str(result[i]))
??? f.write(' ')
??? i += 1
f.close()

?