由前面的博客介绍可知，我们用熵来衡量系统的混沌度。写到这里，我突然想起了《三体》对低熵体和高熵体的精彩描述:

宇宙的熵在升高，有序度在降低，像平衡鹏那无边无际的黑翅膀，向存在的一切压下来，压下来。可是低熵体不一样，低熵体的熵还在降低。有序度还在上升，像漆黑海面上升起的磷火，这就是意义，最高层的意义，比乐趣的意义层次要高 。

 很佩服大刘的文笔，能把硬核的知识描写的如此有诗意。
 在前些日子，我读到了一篇有些年头的论文《THE USE OF ENTROPY TO CHARACTERIZE CODING GENES ​FOR TRANSLATION OF PROTEINS FROM DNA》,其中的一句话引起了我的兴趣即：

This reflects the expectation that sequences containing a definite biological information should have a lower entropy than those which lack any that information.
 这提示我们可以用信息熵来找承担生命信息的序列片段。进一步的查阅相关资料可知：DNA中，能编码蛋白质的编码区域，由于此区域的核苷酸的改变会导致生物性状发生改变，所以在长期进化过程，该区域会非常保守，碱基分布的随机度较低，故而其熵较小。并且在真核生物的DNA序列中编码区的信息熵是小于非编码区的信息熵值，而且两者的差值还比较大。

 在实验室里，我们得到一段DNA序列后，如果想知道编码出的蛋白质序列等，我们还得对它的mRNA进行测序，但是mRNA降解很快，半衰期极短。所以我们能否绕过测序mRNA，直接由DNA序列预测编码序列。
  要想做到这一点，我们需要有一点点的生物学知识。以真核生物为例子，DNA转录成rna,而后rna进行修饰（其中的过程很精细，mrna得进行剪接和加polyA尾巴、加“帽子”），最后成熟的mrna从5‘端起始密码子到3’端终止密码子构成一个开放阅读框（orf)，然后由携带反密码子的trna和细胞机器——核糖体进行多肽链的合成。即生命信息从dna流向rna再流向氨基酸序列。而我们现在想要预测dna中哪些片段指导了蛋白质的合成。就相当于预测mrna中的开放阅读框是由dna中哪部分片段转录而来。

 结合以上的图片我们知道开放阅读框（ORF）是结构基因的正常核苷酸序列，它开始于起始密码子截至到终止密码子可编码完整的氨基酸序列，其间不存在使翻译中断的终止密码子。由于一段RNA序列有多种不同读取方式，因此可能存在许多不同的开放阅读框架。比如：
AGTTTTGTCC
正向方向有三种方式：
第一种读取方式：AGT TTT GTC

第二种读取方式：GTT TTG TCC

第三种读取方式：TTT TGT CC
反向也有三种，只不过是向左读取，这里不再举例。
所以一段序列中因为有多种起始密码子和终止密码子和多种读取方式，所以存在多种开放阅读框。并不是所有读码框都能表达出蛋白产物，承担一定的生物功能，但是表达蛋白质的区域一定是在某一个开放阅读框中
 至此，我们形成了一个初步的思路。如果有一段DNA序列，我们先找出其中的所有的可能的orf。那如何找呢？我们要知道真核生物中mRNA的起始密码子是AUG，所以在对应的dna编码链中是ATG。终止密码子是UAG、UAA、UGA，在DNA模板链中是TAG、TAA 、TGA。所以我们要在DNA序列中按正向的三种读取方式和反向的三种读取方式找到以ATG开头的，同时以TAG或TAA或TGA结束的序列。满足以上条件的序列片段是一个可能的编码序列。 所以我们可能找到许多条候选的orf。
 下面我们来考虑如何找到其中最有可能承担遗传信息的序列。 结合本文一开始所写：编码蛋白质的片段的熵是低于非编码片段的熵。那如何用这一结论来进行预测呢？思路如下：如果现在有10个碱基的序列，其中第4个碱基到第8个碱基是一个可能的orf，认为它是编码序列，命名为序列A,那么第1个碱基到第3个碱基、第9个到第10个碱基是非编码的,分别命名为序列B、C。则可以计算出一个信息熵差值为[entropy_B]+[entropy_C]-[entropy_A]。所以当找到多个orf时，同时也有多个熵差值，最后我们可找到其中最大的那组熵差值所对应的orf，就可以认为它是最有可能承担遗传信息的序列片段，即预测它为编码序列。
  香农熵函数：
$H ( x ) = - sum _ { i } P left( x _ { i } right) log P left( x _ { i } right)$

import math
f=open('D:se.txt')#从ncbi中下载的基因序列文件。并把第一行序列信息删去。
se_map=[]
start=[]#记录起始密码子的位置
end=[]#记录终止密码子的位置
for line in f:
    s=line.strip()
    ss=list(s)
    for i in ss:
        se_map.append(i)
length=len(se_map)
base=['A','G','C','T']
k_tuple=[[]for i in range(64)]
def creat_tuple():  
    t=0
    for i in range(4):
        for ii in range(4):
            for iii in range(4):
                k_tuple[t].append(base[i])
                k_tuple[t].append(base[ii])
                k_tuple[t].append(base[iii])
                t+=1
def find_start():
    i=0
    t=0
    while t<3:
        while i<length-2:
            if se_map[i]=='A' and se_map[i+1]=='T' and se_map[i+2]=='G':   
                start.append(i)
               
            i=i+3
        i=0
        t+=1
        i=t
def find_end():
    t=0
    t_t=len(start)
    while t<t_t:
        i=start[t]
        while i<length:
            if se_map[i] =="T" and se_map[i+1]=="A" and se_map[i+2]=="G":
                end.append(i)
                break
            if se_map[i] =="T" and se_map[i+1]=="A" and se_map[i+2]=="A":
                end.append(i)
                break
            if se_map[i] =="T" and se_map[i+1]=="G" and se_map[i+2]=="A":
                end.append(i)
                break  
            i+=3
        t+=1
find_start()#记录起始密码子位置
find_end()#记录对应的终止密码子位置
creat_tuple()#生成4*4*4种密码子
entropy=[]#存储香农熵
fre_all=[]#存储每个orf中包含的碱基个数
fre_singel=[]#存储每种密码子的频数
cnt=0
s=0
i_i=0
while i_i<64:
    fre_singel.append(0)
    i_i+=1
while cnt<len(start):
    i_i=0
    while i_i<64:
        fre_singel[i_i]=0
        i_i+=1
    s=(end[cnt]-start[cnt])
    fre_all.append(s)
    st=start[cnt]
    en=end[cnt]
    while st <en:
        i=0
        while i < 64:
            if se_map[st]==k_tuple[i][0] and se_map[st+1]==k_tuple[i][1] and se_map[st+2]==k_tuple[i][2]:
                fre_singel[i]+=1#记录此个orf中每个密码子出现的频数。需注意：我们是窗口滑动式遍历。即，窗口大小为3，步长为1。
                break
            i+=1
        st+=1#步长为1
    i=0
    singel_score=0
    final_score=0
    while i<64:
        ss=fre_singel[i]/s#计算密码子的频率。因为步长为1，所以理论密码子总数为该段序列的碱基数目s。
        if ss!=0:
            singel_score=-ss*(math.log(ss,2))#计算香农熵
        else:
            singel_score=0
        final_score+=singel_score
        i+=1
    cnt+=1
    entropy.append(final_score)
cnt=0
length_start=len(start)
end_d=[]
start_t=[]
entropy_y=[]
while cnt < length_start:
    if end[cnt]-start[cnt]>=100:#我们只要长度大于100的
        end_d.append(end[cnt])
        start_t.append(start[cnt])
        entropy_y.append(entropy[cnt])
    cnt+=1
#computes the entropies of the noncode区域.主体逻辑与上面一致。
length_sequence=len(se_map)
i=0
noncode_entropy=[]
noncode=[]
singal=[]
i=0
while i<64:
    singal.append(0)
    i+=1
length_start=len(start_t)
i=0
while i< length_start:
    noncode=se_map[0:start_t[i]]+se_map[end_d[i]:length_sequence]
    t=0
    length_noncode=len(noncode)
    final=0
    i_i=0
    while i_i<64:
        singal[i_i]=0
        i_i+=1
    while t<length_noncode-2:
        t_t=0
        while t_t<64:
            if noncode[t]==k_tuple[t_t][0] and noncode[t+1]==k_tuple[t_t][1] and noncode[t+2]==k_tuple[t_t][2]:
                singal[t_t]+=1
                break
            t_t+=1
        t+=1
    t_t_t=0
    while t_t_t<64:
        if singal[t_t_t]!=0:
            final+=-singal[t_t_t]/length_noncode*math.log(singal[t_t_t]/length_noncode,2)
        t_t_t+=1
    noncode_entropy.append(final)
    i+=1
i=0
max=0
index=0
while i<length_start:
    if noncode_entropy[i]-entropy_y[i]>0:

        if max<(noncode_entropy[i]-entropy_y[i]):
            max=noncode_entropy[i]-entropy_y[i]
            index=i
    i+=1  
f.close()  
print(max)#  max(非编码区的香农熵减去编码区的香农熵的值)
print(index)#最大熵差值组对应的下标
print(noncode_entropy[index])#非编码区的熵
print(entropy_y[index])#编码区的熵
print(start_t[index])#预测出来的，最有可能编码蛋白的orf的起始密码子的位置。（从零开始计数）
print(end_d[index])#终止密码子的位置。

运行的结果：