基于因子模型对G蛋白偶联受体序列进行聚类

Clustering of G Protein-Coupled Receptor Sequences Based on Factor Model

摘要: G蛋白偶联受体(G Protein-Coupled Receptors, GPCRs)是一肽类膜蛋白家族，对GPCRs序列进行聚类分析有着重要的理论意义和应用价值。本文根据氨基酸的分类及其物化性质给出了蛋白质序列的特征向量表示，在此基础上用因子分析法对蛋白质序列的特征向量进行降维得到了因子模型，进而利用因子模型分析了40个GPCRs序列的相似性，并进行聚类分析，得到了较好的结果，为分析比较GPCRs序列提供新的手段。

Abstract: G protein-coupled receptors (GPCRs) is a family of peptide proteins, and it is of great theoretical and practical value to clustering analysis of GPCRs. In this paper, the eigenvector representations of protein sequences are given by the classification and physicochemical properties of amino acids. On the basis of this, dimensions of characteristic vectors of the protein sequences are reduced by factor analysis and obtain factor model. The factor model is used to analyze the similarity of 40 G protein-coupled receptor sequences, simultaneously carrying out the clustering analysis. Better results provide a new approach for analyzing and comparing GPCRs.

1. 引言

每个细胞信号的传递都是借助细胞膜的不同种类的受体，将细胞外的信号传递到细胞内。G蛋白偶联受体(G Protein-Coupled Receptors, GPCRs)就是一个因能结合和调节G蛋白活性而得名的超级膜蛋白家族。作用于GPCRs的信号物质，通过影响细胞的GPCRs而对G蛋白质起作用。因此GPCRs被认为是相似的分子机制而起作用。GPCRs在信号传导中的重要作用，不仅有助于了解细胞信号的传导机制、阐明疾病的致病机理，而且对药物的研究提供新的思路，GPCRs的功能失调会引发许多疾病，如疼痛、色盲症、哮喘等。通过调节有关GPCRs介导的信号传导，可以治疗高血压、紧张和消化道溃疡等病症。大部分药物可通过靶向作用于GPCRs而达到治疗的效果，所以GPCRs在制药领域成为重要的药物作用靶标。根据GPCRs的序列差异，准确地聚类GPCRs序列有着很重要的理论意义和应用价值 [1] [2] 。

蛋白质序列相似性分析是蛋白质序列聚类的关键所在。经典的相似性算法有Needleman-Wunsch算法、Smith-Waterman算法、接触度量矩阵法、矩阵法、T-coffee算法、SIM算法、基于氨基酸物化性的拓扑指数方法和基于LZ复杂度等方法 [3] [4] [5] [6] 。尽管这些算法都有它们各自的优点，但是计算量都比较大。在序列比对算法中，仅用相同或不同来说明两个残基之间的关系，无法描述残基替换对结构和功能的影响程度；矩阵法是构造一个适当的数值矩阵来描述蛋白质序列，然后选取矩阵的不变量把蛋白质序列之间的比较转化成矩阵不变量之间的比较，并且把初始的蛋白质字符串序列转换为特征向量，这些特征向量的维数可以按照自己的愿望进行选择。通过对这些不变量的比较来确定蛋白质序列的相似与不相似 [7] - [12] 。虽然这个方法简洁快速，但是在蛋白质序列转化成特征向量的过程中，总会丢失了一些生物信息；同样在LZ复杂度法中也只考虑了结构信息而忽略了氨基酸的物化性质对蛋白质的结构和功能的影响 [13] ，为此促使人们试图寻找其它更有效的方法来比较蛋白质序列的相似性。

本文中，首先在氨基酸的物化性质表征蛋白质序列的基础上，把蛋白质序列转化成11维的特征向量；其次，根据20种氨基酸的极性、非极性、疏水性、亲水性将其分为四类：极性且亲水性 $\left(pq\right)$ 、极性且疏水性 $\left(pr\right)$ 、非极性且亲水性 $\left(sq\right)$ 和非极性且疏水性 $\left(sr\right)$ ，将这四类氨基酸两两连接得到16个特征子列，并计算了16个特征子列在蛋白质序列中出现的频率，利用此频率将蛋白质序列转化成16维特征向量；最后，用因子分析法把蛋白质序列的特征向量进行降维得到因子模型，进而利用因子模型分析40个G蛋白偶联受体序列的相似性，并对其进行聚类分析。

2. 因子模型

$X={\left(X_1,X_2,\cdot \cdot \cdot ,X_p\right)}^{\text{T}}$ 是可观测的随机向量， $E\left(X\right)=\mu$ ， $Var\left(X\right)=\sum$ ， $F={\left(F_1,F_2,\cdot \cdot \cdot ,F_m\right)}^{\text{T}}\left(m<p\right)$ 是不可观测的随机向量， $E\left(F\right)=0$ ， $Var\left(F\right)=I_m$ ，又设 $\epsilon ={\left({\epsilon }_1,{\epsilon }_2,\cdot \cdot \cdot ,{\epsilon }_p\right)}^{\text{T}}$ 与 $F$ 不相关，且 $E\left(\epsilon \right)=0$ ， $D\left(\epsilon \right)=\text{diag}\left(d_1^2,d_2^2,\cdot \cdot \cdot ,d_p^2\right)=D$ 。设 $X$ 满足：

$\{\begin{array}{l}{X}_1=a_{11}{F}_1+a_{12}{F}_2+\cdot \cdot \cdot +a_{1m}{F}_m+{\epsilon }_1\\ X_2=a_{21}{F}_1+a_{22}{F}_2+\cdot \cdot \cdot +a_{2m}{F}_m+{\epsilon }_2\\ ⋮\\ X_p=a_{p1}{F}_1+a_{p2}{F}_2+\cdot \cdot \cdot +a_{pm}{F}_m+{\epsilon }_p\end{array}$ (2.1)

矩阵方程 $X_{p\times 1}=A_{p\times m}{F}_{m\times 1}+{\epsilon }_{p\times 1}$ 称正交因子模型，其中 $F_1,F_2,\cdots ,F_m$ 称 $X$ 的公共因子， ${\epsilon }_1,{\epsilon }_2,\cdots ,{\epsilon }_p$ 称 $X$ 的特殊因子。设 $F_1,F_2,\cdots ,F_m$ 分别是均值为0，方差为1的随机变量，即 $D\left(F\right)=I_m$ ；特殊因子 ${\epsilon }_1,{\epsilon }_2,\cdots ,{\epsilon }_p$ 分别是均值为0，方差为 $d_1^2,d_2^2,\cdots ,d_p^2$ 的随机变量，即 $D\left(\epsilon \right)=diag\left(d_1^2,d_2^2,\cdots ,d_p^2\right)=D$ ；各特殊因子之间及特殊因子与公共因子之间都是相互独立的，即 $Cov\left({\epsilon }_i,{\epsilon }_j\right)=0,i\ne j$ 及 $Cov\left(\epsilon ,F\right)=0$ 。 $a_{ij}$ 是第 $j$ 个变量在第 $i$ 个公共因子上的负荷， $A={\left(a_{ij}\right)}_{p\times m}$ 是待估系数矩阵(因子载荷矩阵) [14] 。

因子分析的目标是找出公共因素及特有的因素，即公共因子与特殊因子。在公因子分析中，特殊因子起到残差的作用，但被定义为彼此不相关且和公因子也不相关。而且每个公因子假定至少对两个变量有贡献，否则它将是一个特殊因子。在开始提取公因子时，为了简便，还假定公因子彼此不相关且具有单位方差。在这种情况下，向量 $X$ 的协方差矩阵 $\sum$ 可以表示为：

$\sum =D\left(X\right)=D\left(AF+\epsilon \right)=A{A}^{\prime }+D$ (2.2)

这里 $D=\text{diag}\left(d_1^2,d_2^2,\cdot \cdot \cdot ,d_p^2\right)$ ， $\text{diag}$ 表示对角矩阵。

如果已知 $X$ 协方差矩阵 $\sum$ 和 $D$ ，可以很容易地求出 $A$ 。根据式(2.2)有：

记 ${\Sigma }^{*}=\Sigma -D$ ，则 ${\Sigma }^{*}$ 是非负定矩阵。若记矩阵 ${\Sigma }^{*}$ 的 $p$ 个特征值 ${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_m,\cdots ,{\lambda }_p$ 且 ${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_m>{\lambda }_{m+1}=\cdots ={\lambda }_p=0$ ，及 $e_1,e_2,\cdots ,e_m$ 分别表示m个非零特征值所对应的标准化的特征向量，则 ${\Sigma }^{*}$ 的谱分解式为：

${\Sigma }^{*}={\lambda }_1{e}_1{e^{\prime }}_1+{\lambda }_2{e}_2{e^{\prime }}_2+\cdot \cdot \cdot +{\lambda }_m{e}_m{e^{\prime }}_m=\left(\sqrt{{\lambda }_1}{e}_1,\sqrt{{\lambda }_2}{e}_2,\cdot \cdot \cdot ,\sqrt{{\lambda }_m}{e}_m\right){\left(\sqrt{{\lambda }_1}{e}_1,\sqrt{{\lambda }_2}{e}_2,\cdot \cdot \cdot ,\sqrt{{\lambda }_m}{e}_m\right)}^{\prime }$ (2.3)

只要：

$A=\left(\sqrt{{\lambda }_1}{e}_1,\sqrt{{\lambda }_2}{e}_2,\cdot \cdot \cdot ,\sqrt{{\lambda }_m}{e}_m\right)$ (2.4)

就可以求出因子载荷矩阵 $A$ 。从而求得 $d_i^2=1-{\displaystyle \underset{t=1}{\overset{m}{\sum }}{a}_{ti}^2},i=1,2,\cdots ,p$ ，这时 $A$ 和 $D=diag\left(d_1^2,d_2^2,\cdot \cdot \cdot ,d_p^2\right)$ 为因子模型的一个解，这个解就称为主因子解。

因子模型被估计后，还必须对得到的公因子 $F$ 给出一种明确的解释，它用来反映在预测每个可观察变量中这个公因子的重要性，这个公因子的重要程度就是在因子模型矩阵中相应于这个因子的系数，显然这个因子的系数绝对值越大越重要，而接近0则表示对可观察变量没有什么影响。因子解释是一种主观的方法，有时侯，通过旋转公因子可以减少这种主观性，也就是要使用非奇异的线性变换。

设 $p$ 维可观察变量 $X$ 满足因子模型 $X=AF+\epsilon$ ，设 $\Gamma$ 是任一正交阵，则因子模型可改写为：

$X=A\Gamma {\Gamma }^{\prime }F+\epsilon =A^{*}{F}^{*}+\epsilon$ (2.5)

其中， $A^{*}=A\Gamma$ ， $F^{*}={\Gamma }^{\prime }F$ 。

3. 蛋白质序列的特征向量表示

3.1. 蛋白质序列的11维、2维特征向量表示

根据大量的实验结果，我们提取了对G偶联受体序列相似性分析有影响的氨基酸的11种物化属性 [15] ，它们是： ${\alpha }_c$ ( $\alpha$ 螺旋的 $C$ 端动力)， $C_{\alpha }$ ( $\alpha$ 螺旋接触面积)， $K^0$ (可压缩性)， $P_{\beta }$ ( $\beta$ 折叠趋势)， $R_a$ (在溶剂中的收缩率)， $\Delta ASA$ (溶剂可及表面积)， $PI$ (氨基酸的等电点)， $\Delta G_{hD}$ (吉布斯自由能变性蛋白水化的变化)， $N_m$ (平均中程接触)， $Mu$ (折射率)， $EL$ (长距离的非键能)，并对其进行了标准化和平均化 [15] ，见表1。

根据20种氨基酸的标准化和平均化后的物化属性，对40个G蛋白受体序(https://www.ncbi.nlm.nih. gov/protein/1LMB_4上下载)中所含的20个氨基酸进行统计并计算它们的算数平均数，由于数据篇幅比较大，表2只给出了部分G蛋白受体序列的算术平均数。

Table 1. 20 properties of amino acids

表1. 20种氨基酸的属性取值表

Table 2. The arithmetic mean of 40 proteins

表2. 40种蛋白质的算数平均值

根据得到蛋白质序列对应的11维特征向量，并计算11维特征向量的相关矩阵及其特征值，根据特征值累计贡献率提取相应的主因子，以主因子的方差贡献率作为权重得到因子模型如下：

$x_1=-\text{0}\text{.7378}{f}_1+0.2902f_2$ ， $x_2=\text{0}\text{.6938}{f}_1-\text{0}\text{.5964}{f}_2$ ， $x_3=-\text{0}\text{.4949}{f}_1+\text{0}\text{.8335}{f}_2$ ， $x_4=\text{0}\text{.8966}{f}_1-0.2191f_2$ ， $x_5=\text{0}\text{.9301}{f}_1+\text{0}\text{.0991}{f}_2$ ， $x_6=\text{0}\text{.8680}{f}_1-\text{0}\text{.3875}{f}_2$ ， $x_7=\text{0}.0469f_1+\text{0}\text{.0173}{f}_2$ ， $x_8=-0.0463f_1+\text{0}\text{.1128}{f}_2$ ， $x_9=\text{0}\text{.0220}{f}_1-\text{0}\text{.0383}{f}_2$ ， $x_{10}=\text{0}\text{.2801}{f}_1+\text{0}\text{.3878}{f}_2$ ， $x_{11}=\text{0}\text{.2922}{f}_1+\text{0}\text{.3878}{f}_2$

其中 $x_1$ 表示 ${\alpha }_c$ ， $x_2$ 表示 $C_{\alpha }$ ， $x_3$ 表示 $K^0$ ， $x_4$ 表示 $P_{\beta }$ ， $x_5$ 表示 $R_a$ ， $x_6$ 表示 $\Delta ASA$ ， $x_7$ 表示 $PI$ ， $x_8$ 表示 $\Delta G_{hD}$ ， $x_9$ 表示 $N_m$ ， $x_{10}$ 表示 $Mu$ ， $x_{11}$ 表示 $EL$ ，

$\begin{array}{l}{f}_1=-\text{0}\text{.10958}{x}_1^{*}-\text{0}\text{.01489}{x}_2^{*}\text{-0}\text{.0463}{x}_3^{*}+\text{0}\text{.15327}{x}_4^{*}+\text{0}\text{.31804}{x}_5^{*}-\text{0}\text{.04243}{x}_6^{*}\\ +\text{0}\text{.04687}{x}_7^{*}-\text{0}\text{.04615}{x}_8^{*}+\text{0}\text{.02197}{x}_9^{*}+\text{0}\text{.28006}{x}_{10}^{*}+\text{0}\text{.29216}{x}_{11}^{*}\\ f_2=-\text{0}\text{.04775}{x}_1^{*}-\text{0}\text{.15665}{x}_2^{*}+\text{0}\text{.57288}{x}_3^{*}+\text{0}\text{.03451}{x}_4^{*}+\text{0}\text{.25414}{x}_5^{*}-\text{0}\text{.23404}{x}_6^{*}\\ +\text{0}\text{.01727}{x}_7^{*}+\text{0}\text{.11279}{x}_8^{*}-\text{0}\text{.03826}{x}_9^{*}-\text{0}\text{.00144}{x}_{10}^{*}+\text{0}\text{.3878}{x}_{11}^{\text{*}}\end{array}$

$x_i^{*}\left(i=1,2,\cdots ,11\right)$ 是 $x_i\left(i=1,2,\cdots ,11\right)$ 的标准化，利用因子模型可得蛋白质序列对应的2维特征向量 $\left(f_1,f_2\right)$ 。例如，蛋白质序列Q8MXU2对应的二维特征向量 $\left(f_1,f_2\right)=\left(\text{0}\text{.5598},-\text{0}\text{.0536}\right)$ 。

3.2. 蛋白质序列的16维、6维特征向量表示

根据20种氨基酸的极性、非极性、疏水性、亲水性将其分为四类 [16] ：极性且亲水性 $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\}$ ，将这四类氨基酸两两连接得到16个特征子列： $pqpq$ ， $pqpr$ ， $prpq$ ， $pqsq$ ， $sqpq$ ， $pqsr$ ， $srpq$ ， $prpr$ ， $prsq$ ， $sqpr$ ， $prsr$ ， $srpr$ ， $sqsq$ ， $sqsr$ ， $srsq$ ， $srsr$ ，并计算这些特征子列在蛋白质序列中出现的频率得到蛋白质序列的16维特征向量表示。在此基础上，根据得到蛋白质序列对应的16维特征向量，并计算16维特征向量的相关矩阵及其特征值，根据特征值累计贡献率提取相应的主因子，以主因子的方差贡献率作为权重得到因子模型。利用因子模型可得蛋白质序列对应的6维特征向量 $(f_1,f_2,f_3,f_4,f_5,f_6)$ 。例如，蛋白质序列Q8MXU2对应的二维特征向量 $\left(f_1,f_2,f_3,f_4,f_5,f_6\right)=\left(0.0314,0.0750,-0.0189,-0.0911,0.0319,0.2247\right)$ 。

4. 蛋白质序列相似性及其聚类分析

从https://www.ncbi.nlm.nih.gov/protein/1LMB_4上下载了40个G蛋白受体序列 [17] 。下面采用向量端点之间的平方和距离来构造相似性矩阵。由于数据过多，40个G蛋白偶联受体序列对应的相似性矩阵不列在本文里。把相似性矩阵输入到SAS软件得到了40个G蛋白受体序列的聚类图，如图1~图4所示。

观察聚类图1~图4易看出，40个G蛋白偶联受体序列中：P97772和Q9UGT0，O00222和Q3MIV9，Q93564和Q622H2，Q6ZMQ2和Q14833，P31421和Q14416，Q5RAL3、Q9QYS2和P31422，Q863I4、O15303和P35349，Q93564和Q622H2最相似。根据文献 [13] ，如Q68EF4和Q14833，P47743和P70579，Q9V4U4和Q70GQ8，Q5TZ45和P31424是比较相似的。而通过观察图我们发现，只有16维向量和6维向量所得到的聚类图满足，11维向量和2维向量的到的聚类图只满足一部分。由此可见，根据氨基酸的极性和亲水性得到的蛋白质相似性比较更加合理。而图2和图4是利用因子分析法，把40个G蛋白偶联受体序列对应的11维和16维特征向量降维到2维和6维得到的聚类图。在图2中Q75QW7，Q90ZF3，Q4RJZ9，Q68EF4，Q9V4U4，Q8NHA9，Q4R3P0的位置和图1的有差别。但是通过参考文献 [13] ，我们知道Q75QW7应该和Q4REZ5相近，Q90ZF3、Q4RJZ9和Q4REI8相近，Q14833和Q68EF4相近，Q8NHA9和O15303相近。在图4中Q8MXU2，Q75QW7，Q5RDQ8，Q62916的位置和图3的有差别。而在参考文献 [13] 中，Q75QW7和P91685相近，Q5RDQ8

和Q62916、Q14833和Q68EF4相近。由此可见，基于2维和6维特征向量聚类比11维和16维特征向量聚类效果较好。再把图2与图4与文献 [13] 比较可知，图4与文献 [13] 的结果更一致，这说明，我们选取的20种氨基酸的极性和亲水性对G蛋白偶联受体序列相似性比较中有一定的影响，总之，用低维特征向量聚类过程简单，时间复杂度低等优势，而且实例验证，基于2维特征向量聚类效果也较好。

Figure 1. Cluster analysis of 40 G protein-coupled receptor sequences based on 11-dimensional factor vector

图1. 基于11维特征向量的40个G蛋白偶联受体序列聚类图

Figure 2. Cluster analysis of 40 G protein-coupled receptor sequences based on 2-dimensional factor vector

图2. 基于2维特征向量的40个G蛋白偶联受体序列聚类图

Figure 3. Cluster analysis of 40 G protein-coupled receptor sequences based on 16-dimensional factor vector

图3. 基于16维特征向量的40个G蛋白偶联受体序列聚类图

Figure 4. Cluster analysis of 40 G protein-coupled receptor sequences based on 6-dimensional factor vector

图4. 基于6维特征向量的40个G蛋白偶联受体序列聚类图

5. 总结

因子分析是主成分分析的推广，它是一种降维方法，其目的是用有限个不可观测的隐变量来解释原始变量之间的相关性。在因子模型中，由于因子载荷矩阵不是唯一的，我们利用这一特点可以通过因子的旋转，使得旋转后的因子有更鲜明的实际意义。本文利用因子分析法，对40个G蛋白偶联受体序列进行相似性分析，并将其聚类，得到的结果验证了本方法的可行性。而且，此方法对生物序列进化的研究具有操作简单、时间复杂度低、对序列的长度没有限制等优点。对于本文中我们选取的氨基酸的11种物化性质在不同功能的蛋白质中有没有规律以及分析这些物化性质与蛋白质功能和结构之间的影响都是有待于进一步研究的课题。

基金项目

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