【ROSALIND】【练Python，学生信】28 随机序列的概率

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

题目：

随机字符串（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，以及一个包含20个以内在0和1之间数字的数组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.

需得：一个和A长度相同的数组B，B[k]是以A[k]为GC含量的一个序列与s完全相同的概率的常用对数。

测试数据

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

生物学背景

        基因组DNA并非核苷酸的随机排列，存在着很多反复出现的模序（motif）。如果一段序列在很多物种的基因组中都反复出现，我们就会怀疑它可能承担着某种生物学功能。但同时也要注意，如果一条序列足够长，理论上它可以包含所有可能的短序列，而人类的基因组包含多于三十亿个碱基对。因此如果一段序列反复出现，我们需要先确定这种现象并非仅仅由概率导致。

        理论上来说，越短的序列出现的概率越高，反之，长序列很难随机出现。怎样确定这个长短的临界值就成了关键点。为解决这个问题，本题先从用各碱基出现的概率随机构建序列入手，并计算此序列出现的概率。

数学背景

        数组是有序的元素序列，可以看做仅含有一行的矩阵。A[k]代表数组A中的第k个数。

        当我们假设一个DNA序列的GC含量为x时，G或C出现的频率即为x/2，A或T出现的频率为(1-x)/2,。

        概率相乘时往往得到很小的数字，在计算机存储时导致溢出，取对数是常用的方法。对任意正数x和y，有log10(x⋅y)=log10(x)+log10(y)成立。

思路

        本题思路较简单，可分为如下两步依次实现：

一．读入题目给出的GC含量，计算出G、C、A、T的概率；

二．依次读入序列，将在给定概率下每个位置出现该碱基的概率的对数相加。

代码

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): # 计算各碱基的频率
    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 = 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()