一次探索：基于香农熵预测DNA中编码序列，python实现。

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

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

很佩服大刘的文笔，能把硬核的知识描写的如此有诗意。
在前些日子，我读到了一篇有些年头的论文《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)=−∑iP(xi)log⁡P(xi)H ( x ) = - \sum _ { i } P \left( x _ { i } \right) \log P \left( x _ { i } \right)H(x)=−∑iP(xi)logP(xi)
其中p(x_i)是第i种密码子出现的频率。
以下是具体的程序，用py实现的。恩。。。我只考虑了orf正向的三种读取方式。（后期会有改进的版本）

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=0for 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=0t=0while 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+3i=0t+=1i=t
def find_end():t=0t_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)breakif se_map[i] =="T" and se_map[i+1]=="A" and se_map[i+2]=="A":end.append(i)breakif se_map[i] =="T" and se_map[i+1]=="G" and se_map[i+2]=="A":end.append(i)break  i+=3t+=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=0while i_i<64:fre_singel[i_i]=0i_i+=1s=(end[cnt]-start[cnt])fre_all.append(s)st=start[cnt]en=end[cnt]while st <en:i=0while 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。breaki+=1st+=1#步长为1i=0singel_score=0final_score=0while i<64:ss=fre_singel[i]/s#计算密码子的频率。因为步长为1，所以理论密码子总数为该段序列的碱基数目s。if ss!=0:singel_score=-ss*(math.log(ss,2))#计算香农熵else:singel_score=0final_score+=singel_scorei+=1cnt+=1entropy.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=0length_noncode=len(noncode)final=0i_i=0while i_i<64:singal[i_i]=0i_i+=1while t<length_noncode-2:t_t=0while 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]+=1breakt_t+=1t+=1t_t_t=0while 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+=1noncode_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=ii+=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])#终止密码子的位置。

运行的结果：

一次探索：基于香农熵预测DNA中编码序列，python实现。相关推荐

BERT6mA：使用基于深度学习的方法预测DNA N6甲基腺嘌呤位点
<BERT6mA: prediction of DNA N6-methyladenine site using deep learning-based approaches> Sho Ts ...
Nat. Mach. Intell. | 基于深度学习预测DNA甲基化位点
研究人员开发了一种预测DNA甲基化位点的机器学习算法可以帮助识别致病机制.该论文2020年8月3日发表在"Nature Machine Intelligence"上. 研究人员通过 ...
Kaggle比赛—预测 DNA、RNA 和蛋白质测量如何在单细胞中共同变化
Kaggle比赛-预测 DNA.RNA 和蛋白质测量如何在单细胞中共同变化本次比赛的目标是预测随着骨髓干细胞发育成更成熟的血细胞,DNA.RNA 和蛋白质测量值如何在单个细胞中共同变化.您将开发一个 ...
车联网中基于轨迹预测的无人机动态协同优化覆盖算法
本文由吴壮,唐伦,蒲昊,汪智平,陈前斌联合创作摘要针对城市车联网中出现的基站覆盖空洞及局部流量过载等问题,该文提出了一种基于车辆轨迹预测信息的动态预部署方案.首先,为了训练得到统一的 Seq2Se ...
前馈神经网络中的前馈_前馈神经网络在基于趋势的交易中的有效性（1）
前馈神经网络中的前馈 This is a preliminary showcase of a collaborative research by Seouk Jun Kim (Daniel) and ...
干货！探索单目车辆估计中的中间几何表示
点击蓝字关注我们 AI TIME欢迎每一位AI爱好者的加入! 我们在这项工作中提出一种从单张RGB图片估计车辆在相机坐标系中姿态的方法.与传统方法不同的是,我们不采用先估计观测角再进行转换的二步方法 ...
基于链接预测和卷积学习的Web服务网络嵌入
Web Service Network Embedding based on Link Prediction and Convolutional Learning 这是我读研的第一篇论文,也是花了好几 ...
Tensorflow 使用Bidirectional()包装器构建双向LSTM模型，预测DNA序列功能
循环神经网络(RNN) 循环神经网络RNN能处理时间序列,过去几年中,应用 RNN 在语音识别,语言建模,翻译,图片描述等问题上已经取得一定成功,并且这个列表还在增长.RNN模型的一个代表是LSTM ...
将整本《绿野仙踪》存入纳米级DNA中，高效准确，读取无压力
萧箫发自凹非寺量子位报道 | 公众号 QbitAI 如何将一整本<绿野仙踪>,存进纳米级的DNA里? 现在,德克萨斯大学奥斯汀分校的科学家们做到了. 他们开创了一套新的DNA数据编 ...

一次探索：基于香农熵预测DNA中编码序列，python实现。

一次探索：基于香农熵预测DNA中编码序列，python实现。相关推荐

最新文章

热门文章