本文提出了一種基于對偶四元數配準的蛋白質局部螺旋結構參數擬合（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% 左右。上述結果說明本文提出的螺旋參數擬合方法可以成功地應用于蛋白質 α 螺旋檢測。將本文方法擴展應用于蛋白質其他類型二級結構的檢測是下一步要完成的工作。
