本发明属于系统发育树构建方法，特别涉及基于DFT的RNA二级结构距离计算构建系统发育树的方法。
背景技术：
核糖核酸(缩写为RNA，即RibonucleicAcid)，近几年，越来越多的研究表明RNA在抵御细菌侵扰和治疗肿瘤方面作用重大，研究RNA分子结构的组成和特点也随之成为了一个热点问题。事实上，RNA的二级结构相比原始序列更加具有保守性，如何计算RNA二级结构间的距离用于RNA二级结构的功能预测变得尤为重要。基于比对计算RNA序列间的距离方法，消耗了大量的计算资源，耗时较长，也得不到比较好的比对结果；同时基于比对计算RNA序列间的距离方法，只考虑了RNA的一级序列顺序结构，忽略了RNA二级结构固有的自身的回折后依据特定碱基配对A-U，C-G，G-U之间的氢键形成稳定的二级结构，计算RNA二级结构间的距离不准确。无需比对计算RNA二级结构在不丢失精度的情况下，能够节约计算资源，降低计算复杂度(算法整体的时间复杂度达到O(n2)，其中n为RNA二级结构序列中最长的序列长度)，相比基于比对计算RNA序列间的距离方法(时间复杂度达到了O(m*n2)，其中m表示基于比对的RNA二级结构序列的数量，n表示所有RNA二级结构序列中最长序列的碱基数目)具有明显的优势。
技术实现要素：
针对上述技术问题，本发明提供基于DFT的RNA二级结构距离计算构建系统发育树的方法。
本发明解决上述技术问题的技术方案如下：
基于DFT的RNA二级结构距离计算构建系统发育树的方法，具体步骤如下：
1)RNA二级结构的映射规则：
用A，G，C，U符号分别表示RNA二级结构中未配对的碱基，用A’，G’，C’和U’分别表示RNA二级结构中的配对碱基，则得到RNA二级结构5’端到3’端的特征序列，通过如下公式(1)将RNA二级结构中的自由基和配对基表示为一个(3×N)的二维矩阵：
式(1)中：
i＝1,2,3....N；xi,yi,zi∈[-1,1]；
Ai、Ui、Gi、Ci、A'i、U'i、G'i、C'i分别对应表示从上述特征序列中第一个碱基到第i个碱基中A、U、G、C和A’、U’、G’、C’中的碱基数量，N为RNA二级结构的所有碱基数量；
2)RNA二级结构x、y和z轴序列的离散傅里叶变换：
将步骤1)得到的(3×N)的二维矩阵拆分成只包含[x1,x2,x3......xN]的x轴序列Lx(N)，只包含[y1,y2,y3......yN]的y轴序列Ly(N)，只包含[z1,z2,z3......zN]的z轴序列Lz(N)，分别将Lx(N)、Ly(N)和Lz(N)进行DFT变换，对应得到序列Fx(k)、Fy(k)和Fz(k)；所述Lx(N)的DFT变换公式如下：
式(2)中，k＝0,1,2,...,N-1，N表示RNA二级结构序列的长度；同理Ly(N)和Lz(N)分别进行DFT变换对应得到Fy(k)和Fz(k)序列；Fx(k)在特定的k频率下的DFT功率谱PSx(k)表示如式(3)：
PSx(k)＝|Fx(k)|2,k＝0,1,...,N-1 (3)
同理分别得到Fy(k)和Fz(k)对应的功率谱Psy(k)和PSz(k)；
3)提取功率谱的特征值：
根据步骤2)所得的功率谱PSx(k)、Psy(k)和PSz(k)，按公式(4)提取RNA二级结构x轴上功率谱的特征值：
同理分别得到y轴和z轴上的特征值；则每一个RNA二级结构表示成唯一的三维向量[Px，Py，Pz]’；
4)物种系统发育树的建立：
具体步骤如下：
4.1)使用欧式距离计算不同物种RNA二级结构间的距离dis(i,j)，具体公式如下：
式(5)中，i，j分别表示不同物种，ik和jk表示RNA二级结构的特征向量，用式(5)分别计算出t个物种任意两个物种i和j之间的距离，则得到一个t×t的距离矩阵D；D是一个对角线全部为0的对称矩阵，每一个元素D[i,j]表示了物种i与物种j之间的RNA二级结构的距离；
4.2)将步骤4.1)得到的距离矩阵D导入到分子进化遗传分析软件MEGA6.0中，构建系统发育树。
本发明的有益效果是：无需通过RNA序列比对计算RNA序列间的距离，通过本发明方法的实施，同时可以兼顾RNA序列特征和RNA的二级结构的特征，相比RNA序列比对求解距离的方法来说，节约了计算资源，缩短了计算时间；利用DFT在不丢失信息的情况下挖掘RNA二级结构x、y和z轴坐标系隐藏的信息，从而获得RNA二级结构x、y和z轴序列功率谱，通过分析功率谱的固有特征，提取功率谱的有效特征值作为RNA二级结构的特征值，从而建立物种RNA二级结构间的距离矩阵，根据距离矩阵能够快速且准确生成物种系统发育树，分析物种间的进化关系。
附图说明
图1是CVV-3的部分RNA二级结构的示意图。
图2是用本发明方法对9种RNA病毒的二级结构构造的系统发育树。
图3是用现有技术(参见Li Y,Duan M,Liang Y.Multi-scale RNA comparison based on RNA triple vector curve representation[J].BMC Bioinformatics,2012,13(1):1-12.)的方法对9种RNA病毒的二级结构构造的系统发育树。
图4是用现有技术(参见Zhang Y,Huang H,Dong X,et al.A Dynamic 3D Graphical Representation for RNA Structure Analysis and Its Application in Non-Coding RNA Classification[J].Plos One,2016,11(5).)中的方法对9种RNA病毒的二级结构构造的系统发育树。
具体实施方式
基于DFT的RNA二级结构距离计算构建系统发育树的方法，具体步骤如下：
1)RNA二级结构的映射规则：
近年来的研究显示RNA序列是由核苷酸单元{A，G，C，U}组成的一长链，并能通过自身的回折后依据特定碱基配对A-U，C-G，G-U之间的氢键形成稳定的二级结构，分别用A，G，C，U符号表示RNA二级结构中未配对的碱基，分别用A’，G’，C’和U’表示RNA二级结构中的配对碱，则可得到RNA二级结构的特征序列。如图1所示的CVV-3的部分RNA二级结构，将图1中的RNA二级结构从5’端到3’端进行序列化表示的特征序列为：GCC’U’C’C’GAAG’G’A’G’AU。使用如下公式(1)表示RNA二级结构的可视化。
式(1)中，i＝1,2,3....N,xi,yi,zi∈[-1,1]；
Ai、Ui、Gi、Ci、A'i、U'i、G'i、C'i分别对应表示从上述的RNA二级结构特征序列中第一个碱基到第i个碱基中A，U，G，C和A’，U’，G’，C’中的碱基数量，N为RNA二级结构的所有碱基数量。通过式(1)可以将RNA二级结构中的自由基和配对基表示为一个(3×N)的二维矩阵。
2)RNA二级结构x、y和z轴序列的离散傅里叶变换：
将步骤1)得到的(3×N)的二维矩阵拆分成只包含[x1,x2,x3......xN]的x轴序列Lx(N)，只包含[y1,y2,y3......yN]的y轴序列Ly(N)，只包含[z1,z2,z3......zN]的z轴序列Lz(N)，分别将Lx(N)、Ly(N)和Lz(N)进行DFT变换(离散傅里叶变换，Discrete Fourier Transform，缩写为DFT)，对应得到序列Fx(k)、Fy(k)和Fz(k)。
上述Lx(N)变换的DFT公式如下：
式(2)中的N表示RNA二级结构序列的长度。
同理分别将Ly(N)和Lz(N)进行DFT变换对应得到Fy(k)和Fz(k)序列。Fx(k)在特定的k频率下的DFT功率谱PSx(k)表示为：
PSx(k)＝|Fx(k)|2，k＝0,1,...,N-1 (3)
同理分别得到Fy(k)和Fz(k)对应的功率谱Psy(k)和PSz(k)。不同的RNA二级结构序列长度N不相同，生成的DFT功率谱长度也不一样，不能直接根据功率谱计算RNA二级结构之间的距离。
3)提取功率谱的特征值：
RNA二级结构表示的坐标系对应的功率谱表示为[PSx(k)Psy(k)PSz(k)]’，其中0≤k≤(N-1)。综合分析离散傅里叶功率谱的特征，可以得到：在x轴上，PSx(0)值比较大，对应于1～N-1位置的各功率谱值关于轴(N-1)/2对称，因此累加对应于1～(N-1)/2位置上的各功率谱值表示RNA二级结构的特征值，具体公式如下：
同理分别得到y轴和z轴上的特征值Py和Pz。则每一个RNA二级结构表示成唯一的三维向量[Px，Py，Pz]’。
4)物种系统发育树的建立：根据步骤3)得到的RNA二级结构的特征向量可以计算出RNA二级结构间的距离矩阵，构建系统发育树。
具体步骤如下：
4.1)使用欧式距离计算不同物种RNA二级结构间的距离，具体公式如下：
式(5)i，j分别表示不同物种，ik和jk表示RNA二级结构的特征向量；用式(5)计算可以得到t个物种任意两个物种RNA二级结构之间的距离，得到一个t×t的距离矩阵D。D是一个对角线全部为0的对称矩阵，每一个元素D[i,j]表示了物种i与物种j之间的RNA二级结构的距离；距离越小，表明物种i与物种j之间的亲缘关系越紧密，反之亦然。
4.2)将得到的距离矩阵D导入到分子进化遗传分析软件MEGA6.0中，构建UPGMA系统发育树。
按本发明方法选取常用的9种RNA病毒的二级结构(分别是：AIMV-3,CiLRV-3,TSV-3,CVV-3,APMV-3,PDV-3,LRMV-3,EMV-3和AVII)，构建的系统发育树，如图2所示。由图2可知，9种病毒中LRMV-3和EMV-3，CiRV-3和EMV-3，APMV-3和AlMV-3是比较相似的，而TSV-3与其他病毒的RNA二级结构的距离较远。这表明TSV-3与其他病毒RNA二级结构之间差异性较为明显。用现有技术(参见Li Y,Duan M,Liang Y.Multi-scale RNA comparison based on RNA triple vector curve representation[J].BMC Bioinformatics,2012,13(1):1-12.)的方法对9种RNA病毒的二级结构构造的系统发育树如图3所示，从图3构建的系统发育树可以看到，LRMV-3和EMV-3的相似性关系上距离较远。用现有技术(参见Zhang Y,Huang H,Dong X,et al.A Dynamic 3D Graphical Representation for RNA Structure Analysis and Its Application in Non-Coding RNA Classification[J].Plos One,2016,11(5).)中的方法对9种RNA病毒的二级结构构造的系统发育树，如图4所示，图4的系统发育树，无法反映出CVV-3和EMV-3较近的相似性关系。这说明，采用本发明的方法与现有技术的RNA二级结构距离计算方法相比较，本发明的方法能快速获得较准确的物种发育树。