基于傅里叶功率谱的H1N1病毒血凝素蛋白质序列的比较分析

Comparison and Analysis of H1N1 Hemagglutinin Protein Sequences Based on Fourier Power Spectrum

摘要: 基于经典的HP模型，将不同特征下的H1N1病毒血凝素蛋白质序列转换为数字序列并且用离散傅里叶变换求出相应序列的功率谱。根据这些功率谱建立数学矩函数，并将数字序列转换为多维的矩向量，得到蛋白质序列对应的特征向量。再利用特征向量之间的中间距离对蛋白质序列进行聚类比较分析，得到了较好的结果。这一方法将不同长度的蛋白质序列通过功率谱和力矩将其转化为相同维数的向量，使我们更加容易比较分析生物序列。

Abstract: Based on the classical HP model, the H1N1 hemagglutinin protein sequence under different characteristics was converted into a digital sequence and the power spectrum of the corresponding sequence was calculated using a discrete Fourier transform. According to these power spectra, a mathematical moment function is established, and the digital sequence is converted into a multi-dimensional moment vector to obtain the corresponding feature vector of the protein sequence. Then using the middle distances between the feature vectors to compare and analyze the protein sequences, a good result was obtained. This method converts protein sequences of different lengths through power spectrum and moments into vectors of the same dimension, which makes it easier for us to compare and analyze biological sequences.

1. 引言

在过去几十年中，研究者已经提出了几种分类生物序列的方法。这些方法中的大多数是基于对齐的，其中通过使用选定的评分系统获得最佳比对。这些方法提供了生物序列的准确分类，并且已经开发并成功应用了几种算法 [1] [2] [3] 。然而，它的主要缺点是消耗时间长，这在需要进行快速聚类时(例如新的致命病毒)是不合适的 [4] 。以后，无对齐技术是一种趋势方法，它通常能在同一数据集上给出更快的分类 [5] [6] [7] [8] 。例如，k-mer方法是最流行的无对准方法。为了测量两个序列的不同，收集两个生物序列中k个集合或长度为k的子序列，然后计算它们之间的进化距离 [5] [9] 。k-mer方法给出与基于对准的方法相当的结果，同时计算速度更快 [10] 。离散傅里叶变换(DFT)是信号和图像处理中的强大工具。近年来，DFT越来越多地用于信息处理的各个领域，如基因预测，蛋白质编码区和周期性分析 [11] [12] 。DNA序列的DFT功率谱反映了该序列的核苷酸分布和周期性模式，并且已经应用于鉴定基因组序列中的蛋白质编码区 [13] [14] [15] 。

目前，聚类方法在处理大数据的各个方面扮演着越来越重要的角色，如分析蛋白质之间的相似性、提取蛋白质结构信息等 [16] [17] 。对于H1N1禽流感病毒，谢佳新等人 [18] 采用蛋白质序列进化距离进行聚类并构建进化树，在此基础上对病毒的变异性进化进行了研究；Zhao等 [19] 在傅里叶变换的基础上应用不同的聚类法，对生物序列构建进化树；赵剑等人 [20] 在蛋白质的二位数字表达的基础上结合使用向量的傅里叶变换理论，提出高维共鸣识别法来判别蛋白质序列的相似性。李巍巍等人 [21] 在不同特征描述下对多条H1N1病毒血素蛋白质序列进行比较分析，在不影响表征蛋白质序列的前提下，用16维的特征向量代替已有的表征蛋白质序列的40维特征向量，大大减少了计算的复杂度。但是，此方法只考虑了物化性质而忽略了蛋白质序列的内部结构，得到的结果有一定的局限性。在本文中，在20种氨基酸以及氨基酸的四类理化性质(极性且亲水pq，极性且疏水pr，非极性且亲水sq和非极性且疏水sr)两两连接所得的特征下用傅里叶变换将蛋白质的符号序列转换为数字序列，基于本文的方法得到蛋白质序列对应的特征向量，通过特征向量之间的中间距离对蛋白质序列进行相似性分析并聚类。

2. 材料

自流感病毒H1N1出现以来，世界各地的人们对其进行了研究，通过大量研究表明这种病毒是由禽流感、猪流感和人流感混合而成的。这种病毒的基因由8个长短不一的可编码的10个病毒蛋白的线状负链RNA片段组成。这10个病毒蛋白分别是PB2、PB1、PA、HA、NP、NA、M1、M2、NS1、NS2，其中NS1和NS2为非结构蛋白外，其他均是结构蛋白 [21] 。

本文从NCBI网站中Molecular Databases的Protein Sequence下载了在1902~2013年全球22,455条H1N1型流感病毒中，选取了31条含有血凝素蛋白的蛋白质序列进行研究，如表1。

3. 方法

在信号处理中，时域中的序列通常被转换成频域，使一些重要特征直观化。通过这种转换，没有信息丢失，而且一些隐藏的属性可以被揭示。

离散傅里叶变换是较最常见的转换方法之一。对于长度为N的信号 $f\left(n\right)$ ， $n=1,2,\cdots ,N$ 。在频率k的信号的DFT为 $F\left(k\right)={\displaystyle \underset{n=1}{\overset{N}{\sum }}f\left(n\right){\text{e}}^{-i\left(\text{2π}/N\right)kn}}$ ， $k=1,2,\cdots ,N$ 。频率k处的信号的功率谱被定义为 $PS\left(k\right)={\left|F\left(k\right)\right|}^2$ ， $k=1,2,\cdots ,N$ 。

通过DFT功率谱将蛋白质序列转换成相应的数字序列之后，不同长度的数字序列之间进行相似性比较仍然很困难，解决这个问题的一个常用方法是矩向量，将不同长度的数字序列转换为相同维数的距向量，求出向量之间的中间距离矩阵，利用SAS软件建立基于距离矩阵的系统聚类树。我们将PS-M方法建立在蛋白质序列的不同属性上进行了比较。

Table 1. 31 influenza viruses and their corresponding serial numbers

表1. 31条流感病毒及其对应的序号

3.1. 符号序列的数字表达HP模型

3.1.1. 基于20种氨基酸的功率谱

对于一个长度为N的蛋白质序列，序列中 $s_n$ ( $n=1,2,\cdots ,N$ )属于一个有限的符号集合 $T=\left\{A_1,A_2,\cdots ,A_{20}\right\}$ ，其中 $A_t$ 为20种氨基酸中的一种， $s_n$ 是 $A_1,A_2,\cdots ,A_{20}$ 中的某个字母，符号序列 $s_1,s_2,\cdots ,s_N$ 中符号 $A_t$ 的指示函数为 [22] ：

$u_{A_t}\left(n\right)=\{\begin{array}{l}1s_n=A_t\\ 0其他\end{array}n=1,2,\cdots ,N$ (1)

例如，蛋白质序列EVLVLWGVHHPPTGTDQQS，核苷酸V的相应指示剂序列是 $u_{A_{15}}=01010001000000000000$ 。

通过指示函数得到20个长度为N的二进制数列设为 $w_1,w_2,\cdots ,w_{20}$ ，那么符号序列 $s_1,s_2,\cdots ,s_N$ 可以表示为 ${\displaystyle \underset{t=1}{\overset{20}{\sum }}{u}_{A_t}\left(1\right)w_t},{\displaystyle \underset{t=1}{\overset{20}{\sum }}{u}_{A_t}\left(2\right)w_t},\cdots ,{\displaystyle \underset{t=1}{\overset{20}{\sum }}{u}_{A_t}\left(N\right)w_t}$ ，记 $s\left(n\right)={\displaystyle \underset{t=1}{\overset{20}{\sum }}{u}_{A_t}\left(n\right)w_t}$ ， $n=1,2,\cdots ,N$ 。

因此蛋白质序列对应的离散傅里叶变换为:

$\begin{array}{c}S\left(k\right)={\displaystyle \underset{n=1}{\overset{N}{\sum }}s\left(n\right){\text{e}}^{-\frac{i2\text{π}}{N}kn}}\\ ={\displaystyle \underset{n=1}{\overset{N}{\sum }}{u}_{A_1}\left(n\right)w_1{\text{e}}^{-\frac{i2\text{π}}{N}kn}}+{\displaystyle \underset{n=1}{\overset{N}{\sum }}{u}_{A_2}\left(n\right)w_2{\text{e}}^{-\frac{i2\text{π}}{N}kn}}+\cdots +{\displaystyle \underset{n=1}{\overset{N}{\sum }}{u}_{A_{20}}\left(n\right)w_{20}{\text{e}}^{-\frac{i2\text{π}}{N}kn}}\\ =U_{A_1}\left(k\right)w_1+U_{A_2}\left(k\right)w_2+\cdots +U_{A_{20}}\left(k\right)w_{20}\end{array}$ (2)

其中 $U_{A_t}\left(n\right)$ 为 $u_{A_t}\left(n\right)$ 的离散傅里叶变换，即

$U_{A_t}\left(k\right)={\displaystyle \underset{n=1}{\overset{N}{\sum }}{u}_{A_t}\left(n\right){\text{e}}^{-\frac{i2\text{π}}{N}kn}}={\displaystyle \underset{n=1}{\overset{N}{\sum }}{u}_{A_t}\left(n\right)}\left(\cos \frac{2\text{π}nk}{N}-i\sin \frac{2\text{π}nk}{N}\right)$ , $t=1,2,\cdots ,20$ , $k=1,2,\cdots ,N$ (3)

因此，得到数列的离散傅里叶变换为

$S\left(k\right)={\displaystyle \underset{t=1}{\overset{20}{\sum }}{U}_{A_t}}\left(k\right)w_t$ , $k=1,2,\cdots ,N$ (4)

定义 $u_{A_t}$ 的功率谱为 $P_{A_t}\left(k\right)={\left|U_{A_t}\left(k\right)\right|}^2$ ， $k=1,2,\cdots ,N$ ，原蛋白质序列的功率谱函数为

$P\left(k\right)={\displaystyle \underset{t=1}{\overset{20}{\sum }}{P}_{A_t}}\left(k\right)$ , $k=1,2,\cdots ,N$ (5)

3.1.2. 基于电荷和极性性质的功率谱

蛋白质序列的经典HP模型是以构成蛋白质序列的氨基酸的结构分类到物化特征间的对应关系为基础，将20种氨基酸分为4大类，分别是极性且亲水性(pq)极性且疏水性(pr)、非极性且亲水性(sq)和非极性且疏水性(sr)， $pq=\left\{G\right\}$ ， $pr=\left\{A,V,L,I,F,P\right\}$ ， $sq=\left\{S,C,N,E,T,Q,K,R,H\right\}$ 和 $sr=\left\{W,Y,M\right\}$ 。这也为蛋白质序列的结构与功能的研究提供了新思路，蛋白质序列的组成相似，进而推测出它们的结构和功能也相似，这就是经典HP模型的意义所在 [23] 。

经过分类之后，对任意一个长度为N的蛋白质序列 $S=s_1{s}_2\cdots s_N$ ，其中 $s_i$ ， $i=1,2,\cdots ,N$ 为20种氨基酸中的某一种，进行数据化定义，以非极性氨基酸(NP)为例说明：

$u_{NP}=\{\begin{array}{l}1,s_i\in pq\\ 2,s_i\in pr\\ -1,s_i\in sq\\ -2,s_i\in sr\end{array}$ (6)

显然， $u_{pq}$ 是一长度为N的二进制的数列，将20个氨基酸一一对应于4个不同的向量 $w_1,w_2,w_3,w_4$ 。利用离散的傅里叶变换，可将指示函数得到的蛋白质序列数据离散化：

$U_{pq}\left(n\right)={\displaystyle \underset{n=1}{\overset{N}{\sum }}{u}_{pq}\left(n\right){\text{e}}^{-i\frac{\text{2π}}{N}kn}}$ , $n=1,2,\cdots ,N$ , $k=1,2,\cdots ,N$ (7)

序列的功率谱： $P_{pq}\left(k\right)={\left|U_{pq}\left(k\right)\right|}^2$ ， $k=1,2,\cdots ,N$ ，同样可以得到 $P_{pr}\left(k\right)$ ， $P_{UP}\left(k\right)$ 和 $P_{PP}\left(k\right)$ 。原蛋白质序列的功率谱为

$$ , $k=1,2,\cdots ,N$ (8)

定义j阶距 [21]

$M_j^{pq}=\frac{1}{N_{pq}^{j-1}{\left(N-N_{pq}\right)}^{j-1}}{\displaystyle \underset{k=1}{\overset{N/2}{\sum }}{\left(P_{pq}\left(k\right)\right)}^j}$ (9)

同样可以求得 $M_{pr}^j$ ， $M_{sq}^j$ 和 $M_{sr}^j$ 。我们的实验结果表明 $j=1,2,3,4,5$ 对于精确聚类来说是足够的。因此，每个蛋白质序列可以在20维欧氏空间中作为几何点来实现，即

$\left(M_1^{pq},M_1^{pr},M_1^{sq},M_1^{sr},M_2^{pq},M_2^{pr},M_2^{sq},M_2^{sr},M_3^{pq},M_3^{pr},M_3^{sq},M_3^{sr},M_4^{pq},M_4^{pr},M_4^{sq},M_4^{sr},M_5^{pq},M_5^{pr},M_5^{sq},M_5^{sr}\right)$

3.1.3. 聚类

聚类分析在数据分析领域应用甚广，如在数据挖掘、生物信息学和统计学等领域中扮演这非常重要的角色。聚类分析不仅可以达到物以类聚的效果，还可以探索和提取数据中隐含的新规律和新知识。本

文将基于Q型系统聚类法，对所获得的数据进行聚类分析。设n个样本构成的有限集为 $X=\left\{x_1,x_2,\cdots ,x_n\right\}$ ， $d=d\left(x_i,x_j\right)$ $\left(x_i,x_j\in X\right)$ 是任意两个样本之间的中间距离，记 $D=\left\{d\left(x_i,x_j\right)|x_i,x_j\in X\right\}=\left\{d_0,d_1,d_2,\cdots ,d_m\right\}$ ，其中 $d_0=0<d_1<\cdots <d_m$ 。

4. 结果分析与讨论

根据公式(1)和公式(6)将31条含有血凝素蛋白的蛋白质序列转换为二进制序列和四元序列，利用离散的傅里叶变换及上述的二进制序列和四元序列，可将蛋白质序列数据离散化。由于不同长度的蛋白质序列通过傅里叶变换转换得到的数字序列的长度依然不同，使得分析蛋白质序列之间的相似性仍然很困难，为了解决这一难题，依据公式(9)将不同维数的特征向量转换维相同维数的特征向量，以此来达到蛋白质序列相似性分析的目的。如基于20种氨基酸的功率谱得到血凝素蛋白质的20维特征向量，表2是6种血凝素蛋白质序列的20种氨基酸的部分氨基酸数据(由于篇幅的问题这里不一一列举)。

应用SAS软件对31条血凝素蛋白质序列进行Q型系统聚类，根据上述的特征向量矩阵，先将各研究样本看成单独的一类，确定样本之间的‘距离’公式，再计算新样本与其他类之间的距离(本文采用中间距离法)，重复此过程，直到将所有的变量都找到各自的类别，最后通过SAS软件得到相应的聚类图定义，见图1、图2。

在流感病毒编码的10种病毒蛋白质中，本章选取了有血凝素蛋白质的病毒进行了研究。图1是依据20种氨基酸构造成20维特征向量得到的流感病毒蛋白质序列的聚类图，图2是依据氨基酸的四种理化性质并通过数学力矩函数的思想构造了20维特征向量通过聚类得到流感病毒蛋白质序列的聚类图。例如，两者将31条H1N1病毒血凝素蛋白质序列分为不同类，图1是基于20种氨基酸对血凝素蛋白质序列进行聚类，分类结果为(1)，(2)，(3)，(5)，(13)，(15)，(25)，(26)，(28)，(29)，(8，16、20、22、30)，(4、6、7、10、12)，(17、21、23、24)，(11、14、19)，(9、18、27、31)；图2是基于氨基酸的四类理化性质对血凝素蛋白质序列进行聚类，分类结果为(1)，(2)，(3)，(5)，(13)，(15)，(25)，(26)，(29)，(8、17、21、23、30)，(16、20、22、24、28)，(4、6、10、11、9、12、14、19)，(7、18、27、31)。由图1和图2可知，基于蛋白质序列的不同特征属性，应用本文的方法对血凝素蛋白质序列进行分类的差异较小。将这两种分类结果与文献 [24] 进行比较发现，基于氨基酸的四种理化性质对血凝素蛋白质进行的分类结果更加接近文献 [24] 。第9条病毒蛋白质与第4条病毒蛋白质、第6条病毒蛋白质、第10条病毒蛋白质、

Table 2. Eigenvector data based on 20 amino acids

表2. 基于20种氨基酸的特征向量数据

Figure 1. A clustering diagram based on the power spectrum of 20 amino acids

图1. 基于20种氨基酸的功率谱的聚类图

Figure 2. Clustering of power spectrum based on charge and polarity properties

图2. 基于电荷和极性性质的功率谱的聚类图

第11条病毒蛋白质、第12条病毒蛋白质、第14条病毒蛋白质和第19条病毒蛋白质属于同一类，第28条病毒蛋白质与第24条病毒蛋白质属于同一类，以氨基酸的四种理化性质所进行的分类结果与文献 [24] 高度相似，而基于20种氨基酸所得到的结果并没有做到这一点。第8条病毒蛋白质与第16条病毒蛋白质、第20条病毒蛋白质和第22条病毒蛋白质属于同一类，以氨基酸的四种理化性质所进行的分类结果没有做到这一点，而第7条和第30条在不同特征下的聚类不一样。在本章中，我们只选取了蛋白质的20种氨基酸和四种理化性质对蛋白质序列进行研究，尽管得到的结论与文献 [24] 中的结果很相似，但是在大数据的处理过程中综合应用蛋白质的性质越多，对蛋白质的相似性比较越准确，这是我们在今后的研究中重点进行的工作。此外，傅里叶功率谱构造特征向量来表征蛋白质序列，并结合其数字编码的蛋白质可以完全包含序列的所有信息，可自动提取蛋白质序列特征信息，这正是本章内容研究的重点。

5. 总结

本章基于蛋白质二维数字表达结合高维共鸣识别法判别双序列蛋白质的相似性和在频率域上表示DNA序列的基础上，提出了应用傅里叶功率谱分析多个蛋白质序列的相似性。将DNA序列上传统的研究方法转换到研究蛋白质序列上，主要包括：在经典的HP模型之上，以20种氨基酸和氨基酸的四种理化性质为基础上将蛋白质序列数值化。在此基础上，通过离散的傅里叶变换将数字序列离散化，为了统一离散化序列的维数，再根据定义计算序列的功率谱，并构造向量矩阵，计算中间距离。在上述的基础上，采用系统聚类算法获取分层结构，构造聚类树讨论蛋白质序列的相似性。

本章选取了31条H1N1病毒血凝素蛋白质序列对提出的方法进行验证，在不同属性的基础上，经过反复的验证，将多维的数字序列进行降维，最终我们采用20维的特征向量表征整条蛋白质序列，利用系统聚类算法，对31条蛋白质序列进行分类，实验结果与文献 [24] 高度吻合。因此将基于蛋白质二维数字表达结合高维共鸣识别法判别双序列蛋白质的相似性和在频率域上表示DNA序列方法的结合扩展到对多个蛋白质序列在频率域上的特征的提取，这为研究蛋白质序列提供了更为严谨的方法。这些对基于大数据和结构分析的研究具有积极的意义，将大大降低计算的复杂度。

致 谢

感谢基金项目：辽宁省教育厅科学研究一般项目(No. L2015093)对本论文的支持。同时，也要衷心的感谢本文中引用文章的作者。

参考文献

NOTES

^*通讯作者。