本文提出了一种基于对偶四元数配准的蛋白质局部螺旋结构参数拟合（DQRFit）方法。该方法首先提取蛋白质结构数据中各残基的 C、N 原子坐标，然后用滑动窗口分别构造各段的待配准数据和参考数据。接着以配准前后数据点的距离平方和最小为寻优目标，利用对偶四元数配准算法求解出最优的旋转矩阵和平移向量并计算出该段二级结构的螺旋参数：每周残基数（τ）、螺旋半径（ρ）和螺距（p）。本文通过对偶四元数配准可实现 τ、ρ、p 三个螺旋参数的同时拟合，并且滑动窗口宽度可调以适应不同的误差等级。与传统螺旋拟合方法相比，具有计算复杂度低、抗干扰性好、拟合精度高的优点。将本文算法应用于蛋白质 α 螺旋结构检测，结果表明检测精度与蛋白质二级结构词典（DSSP）相当，而明显优于其它几种传统方法。本文研究结果或可在今后的蛋白质结构分类和功能预测、药物设计、蛋白质结构可视化等领域具有重大意义。
引用本文： 徐永红, 张少伟, 景军, 赵勇, 候飞翔. 基于对偶四元数配准的蛋白质局部螺旋参数拟合方法. 生物医学工程学杂志, 2018, 35(1): 131-138. doi: 10.7507/1001-5515.201610020 复制
引言
蛋白质二级结构检测是蛋白质结构信息学领域的一项重要课题[1]。蛋白质二级结构主要包括 α 螺旋、β 折叠和 β 转角等。α 螺旋是蛋白质二级结构中最为常见的类型。目前自动检测蛋白质二级结构的传统方法有基于氢键模型的蛋白质二级结构词典（dictionary of secondary structure of proteins, DSSP）[2]、结构识别（structural identification, STRIDE）[3]，以及利用各种几何约束条件和特征的方法，例如：蛋白质二级元素指定（protein secondary element assignment, P-SEA）[4]、蛋白质曲率分析（protein curve, P-CURVE）[5]、线性二级结构元素预测指定（predictive assignment of linear secondary structure elements, PALSSE）[6]、蛋白质结构线段编码（STICKS）[7]、蛋白质二级结构指定程序（XTLSSTR）[8]、蛋白质二级结构指定方法（KAKSI）[9]等。这些方法的检测结果存在一定差异，主要区别在于对二级结构片段末端的定位及更多细微差异之间的识别结果不同（例如：α 螺旋、310 螺旋或 π 螺旋）[1]。目前还没有统一的标准评价这些方法检测效果的好坏。基于氢键模型的方法具有较强的物理化学意义，但忽略了蛋白质原子坐标数据中隐含的几何信息。基于各种几何约束条件和特征的方法只运用 Cα 原子坐标估计参数，所用信息单一。近年来陆续提出了一些利用肽平面内多种原子坐标信息拟合局部螺旋参数的方法，主要用于检测蛋白质二级结构中的螺旋结构。Enkhbayar 等[10]提出了一种基于总体最小二乘法的螺旋拟合方法（helix fitting by a total least squares method, HELFIT），通过拟合三类螺旋参数：螺旋轴、螺旋半径和螺距检测蛋白质的螺旋结构。该方法不易受到噪声的干扰，可以对短序列的结构进行最佳拟合，但是需要给定初值并用迭代法求解，计算复杂度较高。Kneller 等[11]于 2006 年提出了基于螺旋运动的蛋白质二级结构高效表征方法，并将其命名为 ScrewFit 算法，该方法应用四元数对蛋白质的二级结构进行描述，利用相邻肽平面主链上的 C、O、N 原子坐标计算螺旋参数。该算法用到的原子坐标信息太少，只能表达相邻两个肽平面的旋转，无法表达相邻肽平面之间的平移，更不能对多个连续肽平面的空间局部结构进行最佳拟合，影响了螺旋参数计算的准确性。另外，文献[12]利用四元数描述蛋白质的局部螺旋结构，定义了四元数螺旋轴球面距离（quaternion helix axis spherical distance, QHASD）并在单位球上进行螺旋轴可视化，直观地展示了蛋白质的空间取向和分布，但是只根据螺旋轴的球面距离和空间分布无法确定拟合出 α 螺旋在序列上是否连续。
针对上述问题，考虑到蛋白质的原子坐标信息符合点云数据的特征，因此可将螺旋结构的拟合问题转化为点云数据配准问题。配准时需要对目标点进行旋转和平移，对偶四元数正好可以同时表达刚体的旋转和平移[13]，故本文拟采用对偶四元数进行参数拟合。基于以上原因，本文首先介绍了利用配准思想拟合螺旋参数的具体算法，然后验证了本文算法的抗干扰性，之后对所选蛋白质数据集的 α 螺旋结构进行检测并对检测结果进行讨论分析。研究结果显示，本文方法能够拟合出更加精确的螺旋参数，这不仅为蛋白质二级结构检测领域提供了行之有效的工具，也为今后预测蛋白质三级结构的功能奠定了基础。
1 局部螺旋参数拟合方法
本文方法首先将螺旋结构的拟合问题转化为点云配准问题，然后用对偶四元数表示蛋白质序列，最后应用拟合算法计算相应的螺旋参数：每周残基数（τ）、螺旋半径（ρ）和螺距（p）。
1.1 将拟合问题转化为配准问题
螺旋是蛋白质中规则的二级结构，其原子坐标分布符合点云数据的特征，本文将该配准思想应用到基于对偶四元数配准的蛋白质局部螺旋结构参数拟合（dual quaternions registration fitting, DQRFit）的过程中，利用拟合出的螺旋参数 τ、ρ、p 对局部螺旋结构进行检测，从而将螺旋拟合问题转化为特征点集配准问题。首先提取某段蛋白质结构数据中各残基的 C、N 原子坐标，构造数据点集，如式（1）、（2）所示，然后将两者对应坐标求差，记作 P，如式（3）所示。拟合时利用滑动窗口分别构造各段的待配准数据和参考数据，滑动窗口的长度为 n。式（3）中，p1 至 pn 为 1 × 3 的行向量，C、N、P 均为 n × 3 的矩阵。设所提取蛋白质结构的总残基个数为 m，则滑动窗口需满足 3 ≤ n ≤ m。为了保证较好的拟合效果，以 P 中前 n — 1 行的数据为参考数据，第 2 至 n 行的数据为待配准数据，分别记为
和
，
，拟合过程如图 1 所示。
 |
 |
 |
待配准数据
经配准后得到新的数据集
，如式（4）所示：
'/> |
以
和
中对应点的距离平方和为误差函数，如式（5）所示，利用对偶四元数配准算法对误差函数寻优求解出最优的旋转矩阵和平移向量。
 |
图1
利用 C、N 原子坐标拟合 α 螺旋示意图
Figure1.
α helix fitting using the coordinates of C and N atoms
1.2 蛋白质序列的对偶四元数表示
对偶四元数为四元数理论和对偶数理论的结合及扩展[14]，可以同时表示刚体的旋转和平移。它是由两个四元数组成，含有 8 个标量，如式（6）所示：
 |
其中，
是一个对偶四元数，
是对偶四元数的实部，
是对偶四元数的对偶部，两者皆为四元数，下角标 r 代表刚体的旋转，下角标 t 代表刚体的平移。另外，可用对偶向量
和对偶角
表示对偶四元数
，如式（7）所示：
 |
其中，
，为对偶向量；
，为对偶角度；ε 是对偶四元数的对偶标记，该式可进一步表示为：
 |
 |
其中，n 为单位旋转轴向量，
为刚体上某一点绕旋转轴 n 的旋转角度，d 为刚体沿旋转轴方向的平移距离。
在三维空间中描述刚体运动时只需要 6 个独立变量，上式需满足以下两个约束条件：
 |
 |
按照文献[13]中将点云数据转换成四元数的方法，将待配准数据和参考数据表示为如式（12）所示的对偶四元数形式，其中，p 代表由 C、N 原子坐标差构成的向量：
 |
1.3 拟合算法
对偶四元数的实部代表刚体的旋转，根据旋转矩阵与四元数的转换关系[15]，对应的旋转矩阵表示如下：
 |
其中，
，
，
为斜对称矩阵。
式（13）还可以表示成另一种矩阵形式：
 |
其中，
，
。
对偶四元数的对偶部代表刚体的平移，其四元数形式表示为：
 |
其中，
，t 即为所求的平移向量。
 |
因此，配准后得到的新四元数序列
和参考四元数序列
的最小距离平方和公式如式（17）所示：
'/> |
将式（17）展开得到：
'/> |
其中，
，
，
，
。
通过将约束方程（10）、（11）与误差方程（5）联立得到式（19）的最优解：
 |
其中，
和
为拉格朗日乘子，对误差函数求导可得：
 |
 |
方程（21）左右两边同时乘以
可得：
 |
由于 C2、C3 均为斜对称矩阵，所以：
 |
由式（21）可得：
 |
将式（24）代入（20）化简可得：
 |
此时求出的
使得 D 取得最小值 Dmin。
将式（20）的左右两边同时乘以
，化简得到：
 |
式（21）的左右两边同时乘以
，得到：
 |
将式（26）、（27）代入式（18）化简整理得：
 |
由上述推导可知，当 λ1 取最大值 λmax 时，D 取得最小值 Dmin，令
 |
因为 C1 为实对称矩阵，C2、C3 均为斜对称矩阵，所以式（29）可化简为：
 |
则式（25）可改写为如下形式：
 |
由于 C1，C2，C3 为已知量，故可结合式（30）、（31）求得
，再将
分别代入式（13）和（15）方可求出旋转矩阵 R 及平移向量 t。
具体拟合方法步骤如下：
（1）以肽平面中 C、N 原子的坐标差作为原始输入，构建待配准数据集
和参考数据集
，根据公式（12）并将其表示为相应的对偶四元数形式；
（2）根据式（31）计算 M 的最大特征值
以及最大特征值对应的特征向量
；
（3）根据步骤（2）中得到的特征向量
分别代入式（13）、（15），利用对偶四元数实部对应旋转矩阵，对偶部对应平移向量的关系，计算旋转矩阵 R 和平移向量 t；
（4）将步骤（3）中求得的
代入式（8），计算对应的旋转角
和螺旋轴 n；
（5）将步骤（4）中的旋转角
和螺旋轴 n 代入表达式（32）、（33）、（34），计算 τ、ρ、p 三个螺旋参数，其计算表达式如下[16]：
 |
 |
 |
其中
，
。
2 算法验证与应用
实验所用数据均选自蛋白质数据库（protein data bank, PDB, http://www.rcsb.org/pdb/）。本文方法在两个实验数据集 T 和 S 上进行了验证，数据集 T 含有 100 个由 X 射线晶体衍射法获得的蛋白质结构序列[17]，其中有 50 个蛋白质的结构分辨率在 1.0～2.0 Å之间，其余蛋白质的结构分辨率大于等于 2.5 Å。数据集 S 为 PDB Select25 数据集的一个子集[17]，该子集包含 2 144 种序列同源性小于 25% 的蛋白质结构序列，所有实验均在 Matlab 2009b 环境下进行。
2.1 算法验证
对于数据集 S，当 n 取不同值时，对应螺旋参数的误差棒图如图 2 所示，滑动窗口宽度可调以适应不同的误差等级。当 n ≤ 9 时，螺旋参数随着滑动窗口 n 的逐渐增大无明显变化，误差较小。当 n > 9 时，随着 n 的增大 τ、ρ、p 的值均偏离理想值，且方差越来越大。这是因为随着 n 的增大，非螺旋结构被包含到滑动窗口中，使得对某些短序列蛋白质的螺旋参数拟合误差逐渐增大（图中未画出 n > 14 的部分）。所以本文考虑到数据集中的某些蛋白质可能具有较短的螺旋序列，所有实验中 n 值取 4。
图2
n 取不同值时三种螺旋参数的误差棒图
Figure2.
The error bar charts of the three helix parameters with n of different values
分别应用本文算法和 ScrewFit 算法对数据集 S 加噪前和加噪后的特征参数 τ、ρ、p 进行求解并比较，三个特征参数的取值范围分别为
，
，
，连续有 3 个及 3 个以上的螺旋参数同时满足阈值条件，则确定为 α 螺旋结构。应用本文算法和 ScrewFit 算法对数据集 S 中的 α 螺旋结构进行拟合，比较 3 个特征参数的均值与方差大小。如表 1 所示，本文算法的 3 个螺旋参数的均值更接近理想值，除了 ρ 的方差比由 ScrewFit 算法得到的 ρ 的方差略大之外，其余两个参数的方差均缩小了一半左右。这是因为 ScrewFit 算法只利用相邻两个肽平面的原子坐标进行拟合，可用信息有限，本文算法利用滑动窗口同时结合 n 个肽平面的原子坐标信息进行拟合，从而降低了计算误差。
表1
DQRFit 和 ScrewFit 两种算法对应的螺旋参数均值方差比较
Table1.
Comparison of mean and square deviation values between DQRFit and ScrewFit algorithms
在相同噪声环境下，本文算法和 ScrewFit 算法的抗噪性对比图如图 3 所示，实验时在提取的 C、N、O 原子坐标中加入随机正态分布噪声，当噪声强度低于 0.06 nm 时，本文方法对应的 τ、ρ 值无明显变化，p 值在噪声强度高于 0.03 nm 时有所下降，所有螺旋参数变化范围均在误差允许范围之内；而 ScrewFit 算法对应的 τ、ρ 及 p 值在噪声强度高于 0.02 nm 时开始降低，严重偏离理想值。实验结果表明：本文算法比 ScrewFit 算法的稳定性更好。
图3
不同噪声强度下 DQRFit 和 ScrewFit 算法稳定性对比图
Figure3.
Comparison of the stability of DQRFit and ScrewFit algorithms under different noise intensities
另外在运行时间上，对于同一组标准 α 螺旋数据进行拟合，将本文 DQRFit 算法与 ScrewFit 算法的运行时间进行比较，如表 2 所示。
表2
DQRFit 和 ScrewFit 两种算法计算 α 螺旋参数的运行时 间比较
Table2.
Comparison of running time of calculating the α helix parameter between DQRFit and ScrewFit algorithms
由表 2 可知，本文算法的运行时间要低于 ScrewFit 算法的运行时间，这就体现了对偶四元数运行时间低，效率高的优点。
2.2 应用研究
为验证本文方法的有效性，将本文算法与其他 8 种方法相比较，对蛋白质结构数据集 T 的 α 螺旋序列进行检测。序列一致性指标计算方式为
[18]。其中 k 是进行比较的两种方法同时测定为 α 螺旋结构的蛋白质序列的长度，而 k1、k2 分别为两种方法单独测定时彼此不一致的 α 螺旋序列的长度。各方法检测出的 α 螺旋序列所含残基总数分别为：8 919、9 052、9 379、12 332、8 282、8 309、9 392、9 704、9 245。
如表 3 所示，结果一表明，应用本文算法检测出的 α 螺旋序列和 DSSP 的序列一致性达到了 93.5%，比和 DSSP 序列一致性较高的 STRIDE 还高出了约 2.9%。而 ScrewFit 和 DSSP 的序列一致性仅有 85.5%，和本文方法的序列一致性仅有 83.2%。这是由于 ScrewFit 倾向于短螺旋序列的检测，这些短螺旋序列一般含有 4～10 个残基。而 DQRFit 拟合螺旋序列时考虑了连续多个肽平面内的原子坐标信息，因此对含有任意残基个数螺旋序列拟合效果均较好。另外，DQRFit 和 PALSSE、STICKS 两种方法的一致性较低，分别为 75.6% 和 76.1%，原因是 PALSSE 和 STICKS 只能检测出线性度高的螺旋结构，对于含有特殊结构［例如：结节（kink）］的蛋白质的 α 螺旋结构，其检测能力大幅下降。另外，表 3 及表 4 中加粗部分的结果旨在便于观测本文方法与其他方法的结果差异，并无其他含义。
表3
不同方法精度比较结果一（数据集 T）
Table3.
Comparison of accuracy among different methods on set T
由于数据集 T 中没有包含由核磁共振技术（nuclear magnetic resonance, NMR）获得的蛋白质结构数据，为进一步验证本文算法的准确性，应用本文算法对数据集 S 进行检测，该数据集既含有由 X 射线晶体衍射法获得的蛋白质结构数据，又含有由 NMR 法获得的蛋白质结构数据。检测结果如表 4 所示。
表4
不同方法精度比较结果二（数据集 S）
Table4.
Comparison of accuracy among different methods on set S
如表 4 所示的结果二表明，应用本文算法检测出的 α 螺旋序列和 DSSP 的序列一致性达到了 93.8%，比和 DSSP 序列一致性较高的 STRIDE 还高出了约 2.7%。DSSP 和 STRIDE 均是基于氢键模式对蛋白质二级结构进行检测，故两者所测得的序列一致性较高。而 ScrewFit 和 DSSP 的序列一致性仅有 84.0%，和本文方法的序列一致性仅有 81.9%。本文方法和 ScrewFit 都是通过拟合螺旋参数 τ、ρ、p 对 α 螺旋进行描述，充分考虑了蛋白质序列的空间几何特征，但是 ScrewFit 只考虑了相邻两个肽平面的旋转，而 DQRFit 考虑了连续 n 个肽平面的旋转和平移，故拟合出的螺旋参数更加精确。DQRFit 和 PALSSE、STICKS 两种方法的一致性仍然较低，分别为 76.0% 和 78.7%。上述结果表明本文方法和目前被广泛采纳的 DSSP 方法检测能力相当。
3 结论
本文提出了一种基于对偶四元数配准的螺旋参数计算方法，利用点云配准思想拟合螺旋参数 τ、ρ、p。与 HELFIT 相比，DQRFit 无须赋初值，除了可以同时计算螺旋轴 n、螺旋半径 ρ 和螺距 p 外，还可计算新的参数 τ，拟合精度提高。与 ScrewFit 相比，DQRFit 利用对偶四元数理论可同时求解出旋转矩阵和平移向量，计算复杂度降低。由于滑动窗口的引入，融合了连续多个肽面内多种原子的坐标信息，使得拟合的螺旋参数更加精确，从而对蛋白质结构片段末端的定位更加准确。最后在蛋白质结构数据集 S 和 T 上的验证结果表明：本文方法比 ScrewFit 更稳定；同其他 7 种经典方法检测出的 α 螺旋序列一致性较高；和检测效果较好的 DSSP 方法序列一致性在 93% 以上，而 ScrewFit 与 DSSP 方法序列一致性仅达到 84% 左右。上述结果说明本文提出的螺旋参数拟合方法可以成功地应用于蛋白质 α 螺旋检测。将本文方法扩展应用于蛋白质其他类型二级结构的检测是下一步要完成的工作。
引言
蛋白质二级结构检测是蛋白质结构信息学领域的一项重要课题[1]。蛋白质二级结构主要包括 α 螺旋、β 折叠和 β 转角等。α 螺旋是蛋白质二级结构中最为常见的类型。目前自动检测蛋白质二级结构的传统方法有基于氢键模型的蛋白质二级结构词典（dictionary of secondary structure of proteins, DSSP）[2]、结构识别（structural identification, STRIDE）[3]，以及利用各种几何约束条件和特征的方法，例如：蛋白质二级元素指定（protein secondary element assignment, P-SEA）[4]、蛋白质曲率分析（protein curve, P-CURVE）[5]、线性二级结构元素预测指定（predictive assignment of linear secondary structure elements, PALSSE）[6]、蛋白质结构线段编码（STICKS）[7]、蛋白质二级结构指定程序（XTLSSTR）[8]、蛋白质二级结构指定方法（KAKSI）[9]等。这些方法的检测结果存在一定差异，主要区别在于对二级结构片段末端的定位及更多细微差异之间的识别结果不同（例如：α 螺旋、310 螺旋或 π 螺旋）[1]。目前还没有统一的标准评价这些方法检测效果的好坏。基于氢键模型的方法具有较强的物理化学意义，但忽略了蛋白质原子坐标数据中隐含的几何信息。基于各种几何约束条件和特征的方法只运用 Cα 原子坐标估计参数，所用信息单一。近年来陆续提出了一些利用肽平面内多种原子坐标信息拟合局部螺旋参数的方法，主要用于检测蛋白质二级结构中的螺旋结构。Enkhbayar 等[10]提出了一种基于总体最小二乘法的螺旋拟合方法（helix fitting by a total least squares method, HELFIT），通过拟合三类螺旋参数：螺旋轴、螺旋半径和螺距检测蛋白质的螺旋结构。该方法不易受到噪声的干扰，可以对短序列的结构进行最佳拟合，但是需要给定初值并用迭代法求解，计算复杂度较高。Kneller 等[11]于 2006 年提出了基于螺旋运动的蛋白质二级结构高效表征方法，并将其命名为 ScrewFit 算法，该方法应用四元数对蛋白质的二级结构进行描述，利用相邻肽平面主链上的 C、O、N 原子坐标计算螺旋参数。该算法用到的原子坐标信息太少，只能表达相邻两个肽平面的旋转，无法表达相邻肽平面之间的平移，更不能对多个连续肽平面的空间局部结构进行最佳拟合，影响了螺旋参数计算的准确性。另外，文献[12]利用四元数描述蛋白质的局部螺旋结构，定义了四元数螺旋轴球面距离（quaternion helix axis spherical distance, QHASD）并在单位球上进行螺旋轴可视化，直观地展示了蛋白质的空间取向和分布，但是只根据螺旋轴的球面距离和空间分布无法确定拟合出 α 螺旋在序列上是否连续。
针对上述问题，考虑到蛋白质的原子坐标信息符合点云数据的特征，因此可将螺旋结构的拟合问题转化为点云数据配准问题。配准时需要对目标点进行旋转和平移，对偶四元数正好可以同时表达刚体的旋转和平移[13]，故本文拟采用对偶四元数进行参数拟合。基于以上原因，本文首先介绍了利用配准思想拟合螺旋参数的具体算法，然后验证了本文算法的抗干扰性，之后对所选蛋白质数据集的 α 螺旋结构进行检测并对检测结果进行讨论分析。研究结果显示，本文方法能够拟合出更加精确的螺旋参数，这不仅为蛋白质二级结构检测领域提供了行之有效的工具，也为今后预测蛋白质三级结构的功能奠定了基础。
1 局部螺旋参数拟合方法
本文方法首先将螺旋结构的拟合问题转化为点云配准问题，然后用对偶四元数表示蛋白质序列，最后应用拟合算法计算相应的螺旋参数：每周残基数（τ）、螺旋半径（ρ）和螺距（p）。
1.1 将拟合问题转化为配准问题
螺旋是蛋白质中规则的二级结构，其原子坐标分布符合点云数据的特征，本文将该配准思想应用到基于对偶四元数配准的蛋白质局部螺旋结构参数拟合（dual quaternions registration fitting, DQRFit）的过程中，利用拟合出的螺旋参数 τ、ρ、p 对局部螺旋结构进行检测，从而将螺旋拟合问题转化为特征点集配准问题。首先提取某段蛋白质结构数据中各残基的 C、N 原子坐标，构造数据点集，如式（1）、（2）所示，然后将两者对应坐标求差，记作 P，如式（3）所示。拟合时利用滑动窗口分别构造各段的待配准数据和参考数据，滑动窗口的长度为 n。式（3）中，p1 至 pn 为 1 × 3 的行向量，C、N、P 均为 n × 3 的矩阵。设所提取蛋白质结构的总残基个数为 m，则滑动窗口需满足 3 ≤ n ≤ m。为了保证较好的拟合效果，以 P 中前 n — 1 行的数据为参考数据，第 2 至 n 行的数据为待配准数据，分别记为
和
，
，拟合过程如图 1 所示。
|
|
|
待配准数据
经配准后得到新的数据集
，如式（4）所示：
|
'/> |
以
和
中对应点的距离平方和为误差函数，如式（5）所示，利用对偶四元数配准算法对误差函数寻优求解出最优的旋转矩阵和平移向量。
|
图1
利用 C、N 原子坐标拟合 α 螺旋示意图
Figure1.
α helix fitting using the coordinates of C and N atoms
1.2 蛋白质序列的对偶四元数表示
对偶四元数为四元数理论和对偶数理论的结合及扩展[14]，可以同时表示刚体的旋转和平移。它是由两个四元数组成，含有 8 个标量，如式（6）所示：
|
其中，
是一个对偶四元数，
是对偶四元数的实部，
是对偶四元数的对偶部，两者皆为四元数，下角标 r 代表刚体的旋转，下角标 t 代表刚体的平移。另外，可用对偶向量
和对偶角
表示对偶四元数
，如式（7）所示：
|
其中，
，为对偶向量；
，为对偶角度；ε 是对偶四元数的对偶标记，该式可进一步表示为：
|
|
其中，n 为单位旋转轴向量，
为刚体上某一点绕旋转轴 n 的旋转角度，d 为刚体沿旋转轴方向的平移距离。
在三维空间中描述刚体运动时只需要 6 个独立变量，上式需满足以下两个约束条件：
|
|
按照文献[13]中将点云数据转换成四元数的方法，将待配准数据和参考数据表示为如式（12）所示的对偶四元数形式，其中，p 代表由 C、N 原子坐标差构成的向量：
|
1.3 拟合算法
对偶四元数的实部代表刚体的旋转，根据旋转矩阵与四元数的转换关系[15]，对应的旋转矩阵表示如下：
|
其中，
，
，
为斜对称矩阵。
式（13）还可以表示成另一种矩阵形式：
|
其中，
，
。
对偶四元数的对偶部代表刚体的平移，其四元数形式表示为：
|
其中，
，t 即为所求的平移向量。
|
因此，配准后得到的新四元数序列
和参考四元数序列
的最小距离平方和公式如式（17）所示：
|
'/> |
将式（17）展开得到：
|
'/> |
其中，
，
，
，
。
通过将约束方程（10）、（11）与误差方程（5）联立得到式（19）的最优解：
|
其中，
和
为拉格朗日乘子，对误差函数求导可得：
|
|
方程（21）左右两边同时乘以
可得：
|
由于 C2、C3 均为斜对称矩阵，所以：
|
由式（21）可得：
|
将式（24）代入（20）化简可得：
|
此时求出的
使得 D 取得最小值 Dmin。
将式（20）的左右两边同时乘以
，化简得到：
|
式（21）的左右两边同时乘以
，得到：
|
将式（26）、（27）代入式（18）化简整理得：
|
由上述推导可知，当 λ1 取最大值 λmax 时，D 取得最小值 Dmin，令
|
因为 C1 为实对称矩阵，C2、C3 均为斜对称矩阵，所以式（29）可化简为：
|
则式（25）可改写为如下形式：
|
由于 C1，C2，C3 为已知量，故可结合式（30）、（31）求得
，再将
分别代入式（13）和（15）方可求出旋转矩阵 R 及平移向量 t。
具体拟合方法步骤如下：
（1）以肽平面中 C、N 原子的坐标差作为原始输入，构建待配准数据集
和参考数据集
，根据公式（12）并将其表示为相应的对偶四元数形式；
（2）根据式（31）计算 M 的最大特征值
以及最大特征值对应的特征向量
；
（3）根据步骤（2）中得到的特征向量
分别代入式（13）、（15），利用对偶四元数实部对应旋转矩阵，对偶部对应平移向量的关系，计算旋转矩阵 R 和平移向量 t；
（4）将步骤（3）中求得的
代入式（8），计算对应的旋转角
和螺旋轴 n；
（5）将步骤（4）中的旋转角
和螺旋轴 n 代入表达式（32）、（33）、（34），计算 τ、ρ、p 三个螺旋参数，其计算表达式如下[16]：
|
|
|
其中
，
。
2 算法验证与应用
实验所用数据均选自蛋白质数据库（protein data bank, PDB, http://www.rcsb.org/pdb/）。本文方法在两个实验数据集 T 和 S 上进行了验证，数据集 T 含有 100 个由 X 射线晶体衍射法获得的蛋白质结构序列[17]，其中有 50 个蛋白质的结构分辨率在 1.0～2.0 Å之间，其余蛋白质的结构分辨率大于等于 2.5 Å。数据集 S 为 PDB Select25 数据集的一个子集[17]，该子集包含 2 144 种序列同源性小于 25% 的蛋白质结构序列，所有实验均在 Matlab 2009b 环境下进行。
2.1 算法验证
对于数据集 S，当 n 取不同值时，对应螺旋参数的误差棒图如图 2 所示，滑动窗口宽度可调以适应不同的误差等级。当 n ≤ 9 时，螺旋参数随着滑动窗口 n 的逐渐增大无明显变化，误差较小。当 n > 9 时，随着 n 的增大 τ、ρ、p 的值均偏离理想值，且方差越来越大。这是因为随着 n 的增大，非螺旋结构被包含到滑动窗口中，使得对某些短序列蛋白质的螺旋参数拟合误差逐渐增大（图中未画出 n > 14 的部分）。所以本文考虑到数据集中的某些蛋白质可能具有较短的螺旋序列，所有实验中 n 值取 4。
图2
n 取不同值时三种螺旋参数的误差棒图
Figure2.
The error bar charts of the three helix parameters with n of different values
分别应用本文算法和 ScrewFit 算法对数据集 S 加噪前和加噪后的特征参数 τ、ρ、p 进行求解并比较，三个特征参数的取值范围分别为
，
，
，连续有 3 个及 3 个以上的螺旋参数同时满足阈值条件，则确定为 α 螺旋结构。应用本文算法和 ScrewFit 算法对数据集 S 中的 α 螺旋结构进行拟合，比较 3 个特征参数的均值与方差大小。如表 1 所示，本文算法的 3 个螺旋参数的均值更接近理想值，除了 ρ 的方差比由 ScrewFit 算法得到的 ρ 的方差略大之外，其余两个参数的方差均缩小了一半左右。这是因为 ScrewFit 算法只利用相邻两个肽平面的原子坐标进行拟合，可用信息有限，本文算法利用滑动窗口同时结合 n 个肽平面的原子坐标信息进行拟合，从而降低了计算误差。
表1
DQRFit 和 ScrewFit 两种算法对应的螺旋参数均值方差比较
Table1.
Comparison of mean and square deviation values between DQRFit and ScrewFit algorithms
在相同噪声环境下，本文算法和 ScrewFit 算法的抗噪性对比图如图 3 所示，实验时在提取的 C、N、O 原子坐标中加入随机正态分布噪声，当噪声强度低于 0.06 nm 时，本文方法对应的 τ、ρ 值无明显变化，p 值在噪声强度高于 0.03 nm 时有所下降，所有螺旋参数变化范围均在误差允许范围之内；而 ScrewFit 算法对应的 τ、ρ 及 p 值在噪声强度高于 0.02 nm 时开始降低，严重偏离理想值。实验结果表明：本文算法比 ScrewFit 算法的稳定性更好。
图3
不同噪声强度下 DQRFit 和 ScrewFit 算法稳定性对比图
Figure3.
Comparison of the stability of DQRFit and ScrewFit algorithms under different noise intensities
另外在运行时间上，对于同一组标准 α 螺旋数据进行拟合，将本文 DQRFit 算法与 ScrewFit 算法的运行时间进行比较，如表 2 所示。
表2
DQRFit 和 ScrewFit 两种算法计算 α 螺旋参数的运行时 间比较
Table2.
Comparison of running time of calculating the α helix parameter between DQRFit and ScrewFit algorithms
由表 2 可知，本文算法的运行时间要低于 ScrewFit 算法的运行时间，这就体现了对偶四元数运行时间低，效率高的优点。
2.2 应用研究
为验证本文方法的有效性，将本文算法与其他 8 种方法相比较，对蛋白质结构数据集 T 的 α 螺旋序列进行检测。序列一致性指标计算方式为
[18]。其中 k 是进行比较的两种方法同时测定为 α 螺旋结构的蛋白质序列的长度，而 k1、k2 分别为两种方法单独测定时彼此不一致的 α 螺旋序列的长度。各方法检测出的 α 螺旋序列所含残基总数分别为：8 919、9 052、9 379、12 332、8 282、8 309、9 392、9 704、9 245。
如表 3 所示，结果一表明，应用本文算法检测出的 α 螺旋序列和 DSSP 的序列一致性达到了 93.5%，比和 DSSP 序列一致性较高的 STRIDE 还高出了约 2.9%。而 ScrewFit 和 DSSP 的序列一致性仅有 85.5%，和本文方法的序列一致性仅有 83.2%。这是由于 ScrewFit 倾向于短螺旋序列的检测，这些短螺旋序列一般含有 4～10 个残基。而 DQRFit 拟合螺旋序列时考虑了连续多个肽平面内的原子坐标信息，因此对含有任意残基个数螺旋序列拟合效果均较好。另外，DQRFit 和 PALSSE、STICKS 两种方法的一致性较低，分别为 75.6% 和 76.1%，原因是 PALSSE 和 STICKS 只能检测出线性度高的螺旋结构，对于含有特殊结构［例如：结节（kink）］的蛋白质的 α 螺旋结构，其检测能力大幅下降。另外，表 3 及表 4 中加粗部分的结果旨在便于观测本文方法与其他方法的结果差异，并无其他含义。
表3
不同方法精度比较结果一（数据集 T）
Table3.
Comparison of accuracy among different methods on set T
由于数据集 T 中没有包含由核磁共振技术（nuclear magnetic resonance, NMR）获得的蛋白质结构数据，为进一步验证本文算法的准确性，应用本文算法对数据集 S 进行检测，该数据集既含有由 X 射线晶体衍射法获得的蛋白质结构数据，又含有由 NMR 法获得的蛋白质结构数据。检测结果如表 4 所示。
表4
不同方法精度比较结果二（数据集 S）
Table4.
Comparison of accuracy among different methods on set S
如表 4 所示的结果二表明，应用本文算法检测出的 α 螺旋序列和 DSSP 的序列一致性达到了 93.8%，比和 DSSP 序列一致性较高的 STRIDE 还高出了约 2.7%。DSSP 和 STRIDE 均是基于氢键模式对蛋白质二级结构进行检测，故两者所测得的序列一致性较高。而 ScrewFit 和 DSSP 的序列一致性仅有 84.0%，和本文方法的序列一致性仅有 81.9%。本文方法和 ScrewFit 都是通过拟合螺旋参数 τ、ρ、p 对 α 螺旋进行描述，充分考虑了蛋白质序列的空间几何特征，但是 ScrewFit 只考虑了相邻两个肽平面的旋转，而 DQRFit 考虑了连续 n 个肽平面的旋转和平移，故拟合出的螺旋参数更加精确。DQRFit 和 PALSSE、STICKS 两种方法的一致性仍然较低，分别为 76.0% 和 78.7%。上述结果表明本文方法和目前被广泛采纳的 DSSP 方法检测能力相当。
3 结论
本文提出了一种基于对偶四元数配准的螺旋参数计算方法，利用点云配准思想拟合螺旋参数 τ、ρ、p。与 HELFIT 相比，DQRFit 无须赋初值，除了可以同时计算螺旋轴 n、螺旋半径 ρ 和螺距 p 外，还可计算新的参数 τ，拟合精度提高。与 ScrewFit 相比，DQRFit 利用对偶四元数理论可同时求解出旋转矩阵和平移向量，计算复杂度降低。由于滑动窗口的引入，融合了连续多个肽面内多种原子的坐标信息，使得拟合的螺旋参数更加精确，从而对蛋白质结构片段末端的定位更加准确。最后在蛋白质结构数据集 S 和 T 上的验证结果表明：本文方法比 ScrewFit 更稳定；同其他 7 种经典方法检测出的 α 螺旋序列一致性较高；和检测效果较好的 DSSP 方法序列一致性在 93% 以上，而 ScrewFit 与 DSSP 方法序列一致性仅达到 84% 左右。上述结果说明本文提出的螺旋参数拟合方法可以成功地应用于蛋白质 α 螺旋检测。将本文方法扩展应用于蛋白质其他类型二级结构的检测是下一步要完成的工作。
