改进：利用拓扑熵进行基因预测

 本文将引入入拓扑熵来预测DNA序列中的编码片段。因为在利用上文思路进行预测时，发现结果很不理想。。。。。上文的代码逻辑是寻找到全部可能到orf，并设置一个长度的阈值，凭借此筛选出长度大于此阈值的orf，然后对筛选出的orf引入香农熵进行预测。但是根据返回的结果发现，最终预测出来的序列的长度总是接近于阈值的长度。也就是说，如果我设置阈值为100，那么它返回结果片段长度就是100附近，阈值是1000，它返回结果片段长度就是1000附近。据此现象，我们意识到了一个重要的问题：比较不同长度的序列的碱基分布随机度用香农熵可能是不准确的。
  那现在面临的问题是如何消除长度的影响？比如两个orf片段，一个长度为10，一个长度为1000。长度为1000的片段中可能有很多重复的密码子，而长度为10的片段有很大的可能是具有不同的密码子。据此，我们认为长度为10的片段是更随机的，更有可能是编码片段。显然这是不准确的。因为当序列长度是1000时候，我们三个碱基三个碱基进行统计，当长度为10的序列，我们仍是三个碱基三个碱基来统计。而共有64种不同的三联体（4的3次方），所以长度为1000 的片段是更容易出现重复的三联体的。所以我们需要引入一个新的量来比较不同长度片段的熵。
 经过文献查阅，我查到了Kloslicki先生的一篇论文《Topological entropy of DNA sequences》。里面提到了利用拓扑熵进行序列的分析

This definition will allow the comparison of entropies of sequences of differing length ,a property no other implementation of topological entropy has been able to incorporate.

  对于拓扑熵， 我们期望它有以下的性质：
 1，它的取值区间范围是0~1。

 2，它的值为0，当且仅当序列有高度重复的子字（subwords)。

 3，它的值接近于1时，表明它有着尽可能多的不同的subwords。

 4，对于不同序列，其对应的拓扑熵是具有可比性的。

并且对于不同的序列，其对应的subwords的长度也是不同的。这里的subwords(子字)可以理解为n联体。比如我们在香农熵预测基因时候是统计三联体（密码子）的出现的频率。这里的subwords就是三联体，即以三个碱基为一组进行统计。然而在拓扑熵中，不同长度的片段所对应的subwords的长度可能是不同的。举例来说，对于长度为10的片段，如果设置滑动窗口为1的话，它最多最多出现10种三联体。而对于长度为10000的片段来说。它很有可能出现全部的64种三联体。故随着序列长度的增长，它对应的subwords的长度也得随值增长。具体的对应关系，由以下的公式给出：
$4^n+n-1\le w<4^{n+1}+(n+1)-1$

n是subwords的长度，即n联体。w：是序列的长度。

 之后，我们引出拓扑熵的计算公式： $H_{\top (w)}=\frac{{\log }_4(p(n))}{n}$
p(n):该段序列中出现的不同的n联体的个数。具体的数学背景不再赘述，这篇博客主要是使用拓扑熵
以下是实现的代码：

import math
f=open('D:\se.txt')
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)
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
def topological_entropy(sequence):
    length=len(sequence)
    temp_n=math.log(length,4)
    up_n=math.ceil(temp_n)#向上取整
    down_n=math.floor(temp_n)#向下取整
    if length>=math.pow(4,down_n)+down_n-1:#得出对应的subwords的长度
        n=down_n
    else:
        n=up_n
    store=[]
    stt=[]
    cnt=0
    while cnt < length-n:
        i=0
        while i<n:
            i_i=cnt+i
            store.append(sequence[i_i])
            i=i+1
        st=''.join(store)#以n个碱基连接成一个字符串
        stt.append(st)
        store.clear()
        cnt=cnt+1
    stt=set(stt)#利用set进行去重
    kind_n=math.log(len(stt),4)#直接统计去重后的结果
    entropy=kind_n/n
    return entropy#返回计算得到的拓扑熵
    
find_start()
find_end()
ex=[]
to_entropy=[]
i=0
flag=0
temp_to=0
index=0
#ex=se_map[start[0]:end[0]+3]
while i<len(start):
    if(end[i]-start[i]>=67):
        ex=se_map[start[i]:end[i]+3]
        temp_to=topological_entropy(ex)
        to_entropy.append(temp_to)
        if flag<=temp_to:
            flag=temp_to;
            index=i
    i=i+1
t=0
index_x=[]
while t<len(to_entropy):
    if flag==to_entropy[t]:
        index_x.append(t)
    t+=1
print(flag)#输出最大的拓扑熵值
print(index)#存放熵值最大的片段的起始密码子的位置，可能有多个，如本例有两组
print(start[index])#输出熵值最大的第一组的密码子位置
print(end[index])#对应的终止密码子的位置
print(se_map[start[index]:end[index]+3])#输出对应的片段

 以上是运行的结果结果截图。