本文提出一種基于四元數表示的蛋白質α螺旋檢測方法。本方法將蛋白質Cα坐標系螺旋軸映射到單位球面上，在此基礎上定義了四元數螺旋軸球面距離(Quaternion Helix Axis Spherical Distance, QHASD)并用該指標對二級結構中的α螺旋進行檢測。基于PDBselect數據庫應用本文方法進行驗證，對α螺旋檢測精度達到91.7%。本文方法具有檢測精度高、計算復雜度低和幾何意義明確的顯著優點。
引用本文： 徐永紅, 褚澤斐. 蛋白質α螺旋檢測的四元數方法. 生物醫學工程學雜志, 2016, 33(1): 155-160. doi: 10.7507/1001-5515.20160028 復制
引言
蛋白質結構與功能是當前生命科學中一項重要研究內容，準確地識別蛋白質的二級結構，對其三維結構的理解至關重要[1]。膜蛋白占蛋白質分子總數的30%，其跨膜區結構大部分為α螺旋，α螺旋檢測精度的高低直接影響膜蛋白跨膜螺旋區段的定位精度，對于藥物設計中作用靶點的確定至關重要。α螺旋是蛋白質中最典型、含量最豐富的二級結構類型。對α螺旋結構的深入研究和精確描述能夠在很大程度上提高蛋白質結構檢測的準確率，并對蛋白質折疊和蛋白質設計研究有重要的推動作用。本文將重點對α螺旋結構進行研究。
蛋白質的α螺旋、β折疊等二級結構是對蛋白質局部結構的一種定性描述。檢測蛋白質的二級結構有多種方法，其中最準確的是利用核磁共振測定出每個原子的三維坐標，從而根據不同二級結構的幾何特點計算出不同的二級結構。應用最廣泛的DSSP(definition of secondary structure of proteins)算法是以主鏈中氨基之間的氫鍵模式來定義蛋白質的二級結構。對蛋白質二級結構的含量還可以用光譜來初步估計，最常用的方法是圓二色譜。利用波長范圍為170~250 nm長紫外線，在所獲得光譜吸收曲線上，α螺旋結構會在208 nm及222 nm兩處同時出現極小值，而無規卷曲和β折疊結構則分別會在204 nm和207 nm處出現單個極小值。另一個較常用的方法是紅外光譜，它可以偵測因氫鍵所造成的胺基震蕩。近幾年來，四元數在分子建模中的應用越來越廣泛[2]，涉及到的方面主要有分子的最小二乘擬合[3]、分子的平均取向和隨機取向分析，以及分子的四元數可視化[4]。
1993年，Srinivasan等[5]提出一種應用局部螺旋軸來測量螺旋規律性的方法。這種方法將兩個不同的蛋白質坐標系在四元數的基礎上進行最小二乘擬合疊加，來進行兩個不同蛋白質空間結構的對比。然而他方法中的參數Φ并不能通過簡單計算得到，而且不能直接與傳統參數Ramachandran相聯系。1999年，Quine[6]在理想蛋白質幾何形狀的螺旋結構定義中引入了四元數，但其方法只能用于理想蛋白質結構中，局限性較大，不能廣泛應用。2004年，Coutsias等[7]應用四元數來尋找蛋白質的最佳旋轉，用于兩蛋白質空間結構的配準和結構比較。2006年，Kneller等[8]利用局部螺旋運動對蛋白質二級結構的幾何情況進行描述。他們同時定義了取向距離σ、螺半徑ρ和直線度參數ξ，通過這些參數值對二級結構進行檢測。他們首次提出了straightness的概念，但是對于相鄰殘基間取向的標量描述不能與Ramachandran角直接聯系。這一年，Hanson等[9]在Kneller和Calligari的基礎上完善了straightness的定義，對蛋白質定義的坐標系可以直接與Ramachandran角相關聯，可以廣泛應用于實際蛋白質數據中。2012年，Hanson等[10]根據Frenet坐標系對蛋白質建造唯一的坐標系，轉換為四元數后在三維球體上對蛋白質的四元數數據進行可視化，構造出四元數地圖，進而區分出不同的二級結構，但并未涉及進行二級結構檢測的指標或方法。
本文在文獻[10]所提出Cα坐標系的基礎上，根據蛋白質局部螺旋結構的幾何意義計算Cα坐標系并從中提取螺旋軸向量，將螺旋軸向量進行球面可視化，提出參數四元數螺旋軸球面距離(Quaternion Helix Axis Spherical Distance，QHASD)，作為蛋白質α螺旋的檢測的指標。
1 理論基礎
1.1 四元數
由于四元數在旋轉的描述方面有著顯著的優越性，現已有許多研究者將四元數應用到生物信息學領域的蛋白質結構分析中。文獻[11]將四元數用于RNA分子研究中的分子取向分析的問題。
四元數q包括4個分量，即一個實部和三個虛部，
| $q\left( \theta ,n \right)=\left( {{q}_{0}},{{q}_{x}}i,{{q}_{y}}j,{{q}_{z}}k \right)$ |
其中，虛部單位i,j,k需滿足：
| ${{i}^{2}}={{j}^{2}}={{k}^{2}}=-1,ij=-ji=k,jk=-kj=i,ki=-ik=j$ |
四元數在表示旋轉時，n為螺旋軸向量，θ為旋轉角，則四元數又表示為
| $q\left( \theta ,n \right)=\left( \cos \left( \frac{\theta }{2} \right),n\sin \left( \frac{\theta }{2} \right) \right)$ |
1.2 四元數與旋轉矩陣的映射關系
傳統方法中通常采用旋轉矩陣描述旋轉。
1.2.1 四元數到旋轉矩陣
若已知單位四元數q(θ,n)=(q0,qxi,qyj,qzk)，可根據關系式(4)。
| $\left\lfloor \begin{matrix} {{q}_{0}}^{2}+{{q}_{x}}^{2}-{{q}_{y}}^{2}-{{q}_{z}}^{2} & 2{{q}_{x}}{{q}_{y}}-2{{q}_{0}}{{q}_{z}} & 2{{q}_{z}}{{q}_{x}}+2{{q}_{0}}{{q}_{y}} \\ 2{{q}_{x}}{{q}_{y}}+2{{q}_{0}}{{q}_{z}} & {{q}_{0}}^{2}-{{q}_{x}}^{2}+{{q}_{y}}^{2}-{{q}_{z}}^{2} & 2{{q}_{y}}{{q}_{z}}-2{{q}_{0}}{{q}_{z}} \\ 2{{q}_{x}}{{q}_{z}}-2{{q}_{0}}{{q}_{y}} & 2{{q}_{y}}{{q}_{z}}+2{{q}_{0}}{{q}_{x}} & {{q}_{0}}^{2}-{{q}_{x}}^{2}-{{q}_{y}}^{2}+{{q}_{z}}^{2} \\ \end{matrix} \right\rfloor $ |
R就是我們所要找的旋轉矩陣，其中，R為正交矩陣。
1.2.2 旋轉矩陣到四元數
對于3×3的旋轉矩陣R，可以通過一個角度θ和一個螺旋軸n表示為R(θ,n)[12]，軸角表示為
| $\left\lfloor \begin{matrix} c+{{n}_{x}}^{2}\left( 1-c \right) & {{n}_{x}}{{n}_{y}}\left( 1-c \right)-s{{n}_{z}} & {{n}_{z}}{{n}_{x}}\left( 1-c \right)+s{{n}_{y}} \\ {{n}_{x}}{{n}_{y}}\left( 1-c \right)+s{{n}_{z}} & c+{{n}_{y}}^{2}\left( 1-c \right) & {{n}_{z}}{{n}_{y}}\left( 1-c \right)-s{{n}_{x}} \\ {{n}_{x}}{{n}_{z}}\left( 1-c \right)-s{{n}_{y}} & {{n}_{y}}{{n}_{z}}\left( 1-c \right)+s{{n}_{x}} & c+{{n}_{z}}^{2}\left( 1-c \right) \\ \end{matrix} \right\rfloor $ |
其中c=cosθ,s=sinθ,n·n=1。
已知R后，可以通過如下的方法得到旋轉角θ和螺旋軸n：
| $trR=1+2\cos \theta $ |
| ${{\cos }^{2}}\left( \frac{\theta }{2} \right)=\frac{1}{4}\left( trR+1 \right)$ |
| $R-{{R}^{'}}=\left\lfloor \begin{matrix} 0 & 2{{n}_{z}}\sin \theta & 2{{n}_{y}}\sin \theta \\ 2{{n}_{z}}\sin \theta & 0 & -2{{n}_{x}}\sin \theta \\ -2{{n}_{y}}\sin \theta & 2{{n}_{x}}\sin \theta & 0 \\ \end{matrix} \right\rfloor $ |
| $R-{{R}^{'}}=\left\lfloor \begin{matrix} {{r}_{11}} & {{r}_{12}} & {{r}_{13}} \\ {{r}_{21}} & {{r}_{22}} & {{r}_{23}} \\ {{r}_{31}} & {{r}_{32}} & {{r}_{33}} \\ \end{matrix} \right\rfloor =\left\lfloor \begin{matrix} 0 & -2{{n}_{z}}\sin \theta & 2{{n}_{y}}\sin \theta \\ 2{{n}_{z}}\sin \theta & 0 & -2{{n}_{x}}\sin \theta \\ -2{{n}_{y}}\sin \theta & 2{{n}_{x}}\sin \theta & 0 \\ \end{matrix} \right\rfloor $ |
根據式(9)可以求得nx,ny,nz的值，進而得到螺旋軸n。
根據旋轉角θ和螺旋軸n，由公式(3)得到旋轉矩陣R(θ,n)對應的四元數q(θ,n)=(q0,qxi,qyj,qzk)。
2 本文算法
標準α螺旋繞固定軸旋進，該旋進可分解為繞一定點的旋轉和沿固定方向的平移，故同一α螺旋的旋轉軸為一條直線。本文利用四元數來表達蛋白質Cα坐標系間的相對旋轉，定義了反映螺旋軸在單位球面上聚類結果的參數——QHASD，根據QHASD值對α螺旋結構進行檢測。
QHASD的計算步驟：
(1) 從蛋白質空間結構數據中提取出Cα坐標系的主要原子Cα、C、N的空間坐標。 忽略掉H原子，氨基酸的骨架是一個頂端為Cα原子的四面體，另外三個頂點分別為氨基、羧基和R基。如圖 1所示為蛋白質1CRN中兩個氨基酸的基本原子框架，其中CB-SG為R基，可以看到每個氨基酸的骨架原子為N-Cα-C的組合。
圖1
氨基酸的基本結構
Figure1.
Amino acid neighboring structure
(2)在弗萊納公式(Frenet-Serret)的基礎上根據Cα、C、N原子的空間坐標構建蛋白質的Cα坐標系：
| ${{X}_{i}}=\frac{{{C}_{i}}-{{C}_{\alpha i}}}{{{C}_{i}}-{{C}_{\alpha i}}},{{U}_{i}}=\frac{{{N}_{i}}-{{C}_{\alpha i}}}{{{N}_{i}}-{{C}_{\alpha i}}},{{Z}_{i}}=\frac{{{X}_{i}}\times {{U}_{i}}}{{{X}_{i}}\times {{U}_{i}}},{{Y}_{i}}={{Z}_{i}}\times {{X}_{i}}$ |
其中，XYZ分別為列向量，因此得到3×3的單位正交矩陣構成的Cα坐標系序列[F1,F2,…,Fi,…,Fn]，其中Fi=(Xi,Yi,Zi)。弗萊納公式常用來描述粒子在連續可微的曲線上運動時曲線的切向、法向、副法方向之間的關系。每個Cα坐標系Fi均表征了對應氨基酸骨架的空間取向。
(3) 計算表示相鄰Cα坐標系之間的相對旋轉矩陣：
| ${{R}_{i}}={{F}_{i+1}}\centerdot F_{i}^{-1}$ |
其中Fi、Fi+1分別為第i個和第i+1個氨基酸的Cα坐標系，Ri表示相鄰氨基酸對應Cα坐標系之間的相對旋轉矩陣，最終得到旋轉矩陣序列[R1，R2，…，Ri，…，Rn]。
(4) 根據1.2中介紹的旋轉矩陣和四元數的對應關系，將蛋白質的旋轉矩陣序列[R1，R2，…，Ri，…，Rn]轉換為四元數序列[q1，q2，…，qn]。一個3×3的旋轉矩陣R可能被表示為q或-q兩種形式，必須選取其中之一，因此計算相鄰四元數的內積qi·qi+1，當內積為負時，將qi+1改寫為-qi+1。 (5)根據公式(3)從每個四元數中提取螺旋軸，從而可以得到對應的螺旋軸序列[n1,n2,…,ni,…,nn]，因為ni為單位向量，故將螺旋軸序列在單位球面上進行可視化，同時計算相鄰螺旋軸ni與ni+1之間的球面距離QHASD。經過多次計算驗證后設定QHASD的閾值，對蛋白質進行α螺旋檢測。
本文根據同一α螺旋中相鄰氨基酸的相對旋轉是圍繞同一螺旋軸這一幾何特征，基于四元數對蛋白質局部螺旋結構進行描述，在文獻[10]提出的Cα坐標系的基礎上進一步實現四元數螺旋軸序列的球面可視化，利用四元數螺旋軸球面距離來描述蛋白質局部結構的三維空間取向，進而檢測α螺旋。與傳統蛋白質α螺旋檢測算法相比，本文將蛋白質可視化與結構檢測相結合，在描述蛋白質中氨基酸空間取向與分布的同時，利用四元數表達局部螺旋結構進而檢測到α螺旋，該算法效率高，大大降低了計算復雜度。
3 實驗分析
本文利用實際數據對該方法進行驗證，并將本文方法得到的真實蛋白質數據的檢測結果與DSSP、KASKI等方法進行比較。
3.1 實驗及結果分析
從蛋白質結構分類數據庫SCOP(Structural Classification of Proteins)中選取以α螺旋為主的蛋白質數據1A6G，利用MATLAB與R語言混合編程對蛋白質PDB文件進行文本處理后，進行α螺旋結構檢測并與其他算法的檢測結果進行對比。SCOP數據庫是提供關于已知結構蛋白質之間的結構和進化關系的信息，所涉及的蛋白質包括結構數據庫PDB中的所有條目。PDB數據庫(Protein Data Bank)是美國Brookhaven實驗室于1971年建立的大分子蛋白質晶體結構資料數據庫，是互連網上唯一有生物高分子三維結構的數據庫，它的內容主要是根據X射線結晶和核磁共振得到的實驗數據。截止2014年2月，PDB數據庫已含有97 980個結構數據，其中約93%是蛋白質的結構數據。SCOP從總體上將蛋白質進行分類，例如全α型、全β型、以平行折疊為主的α/β型和以反平行折疊為主的α+β型等。本文選取的蛋白質1A6G為全α型。
如圖 2所示為蛋白質1A6G的三維空間結構，可以看到其中包含有6條α螺旋鏈。如圖 3所示為蛋白質1A6G的α螺旋檢測，圖 3(a)為蛋白質1A6G對應四元數序列螺旋軸的球面可視化，可以看到螺旋軸在球面有6個聚集區域，分別對應蛋白質1A6G的6條α螺旋。圖 3(b)為計算得到的相鄰螺旋軸的球面距離，有6段的球面距離明顯減少。經計算驗證，取QHASD的值為0.18時，對α螺旋的檢測精度最高，設定閾值為0.18。
圖2
蛋白質1A6G的空間結構
Figure2.
rotein structure of 1A6G
圖3
蛋白質1A6G的α螺旋檢測
(a)蛋白質1A6G對應四元數序列螺旋軸向量的球面可視化；(b)蛋白質1A6G相鄰螺旋軸球面距離
Figure3. α helix characterization of protein 1A6G(a) protein 1A6G corresponding helix axis vectors of quaternion sequences′ spherical visualization; (b) protein 1A6G′s QHASD
將蛋白質1A6G的Cα坐標系矩陣序列轉換為四元數表示，并對四元數序列中21-37段四元數的(ω,y,z)和(x,y,z)坐標進行球面可視化，根據檢測結果，該段對應的為α螺旋結構。如圖 4所示，(ω,y,z)和(x,y,z)均表示為單位球內的閉合圓環。
圖4
蛋白質1A6G的Cα坐標系對應四元數序列可視化
Figure4.
Visualization of the quaternions corresponding to the Cα frames of protein 1A6G
3.2 檢測結果與精度比較
應用標準算法DSSP和KAKSI等多種廣泛應用的二級結構檢測算法對蛋白質1A6G的α螺旋進行檢測，并與本文的四元數檢測算法結果進行比較。
蛋白質α螺旋的檢測算法中，目前應用最為廣泛的是Ramachandran(1963)[13]經典著作中提出的拉氏圖方法與DSSP算法。拉氏圖算法是提取了局部方位角信息作為參數，來描述蛋白質結構的空間構象。DSSP算法[14]是用于對蛋白質結構中的氨基酸殘基進行二級結構構象分類的標準化算法，由Wolfgang Kabsch和Chris Sander設計，該算法使用PDB格式的原子級分辨率的蛋白質三維結構坐標集數據，依靠以靜電學定義進行的氫鍵識別，以及對主鏈和側鏈二面角的計算，從而得到每個氨基酸殘基的二級結構構象參數。DSSP數據庫是由此算法生成的一個存放蛋白質二級結構分類數據的數據庫，其中包括了PDB數據庫中的所有條目。
如圖 5所示，為上述5種算法對蛋白質1A6G的α螺旋的檢測結果比較，其中α螺旋標識為“H”。
圖5
蛋白質1A6G的α螺旋檢測結果比較
Figure5.
Comparison of α helix characterization result of protein 1A6G
本文算法與其他多種經典算法對蛋白質1A6G的α螺旋的檢測精度進行了比較，具體結果如表 1所示。
表1
蛋白質1A6G的α螺旋檢測精度比較
Table1.
Comparison of α helix characterization accuracy of protein 1A6G
將本文算法與DSSP、PALSSE、STRIDE和STICKS等多種經典二級結構檢測算法對PDBselect數據庫的一個子集中蛋白質數據的α螺旋結構進行檢測，并將檢測精度進行對比。本文所使用的PDBselect數據庫是一套經實驗測定的，其序列之間的同源性小于25%的非冗余PDB數據庫的子集，共2 144種蛋白質結構序列。結果如表 2所示，與其他蛋白質二級結構檢測算法相比，本文算法比四種方法中效果最好的DSSP方法的精度略低，與效果較好STRIDE方法的相比，檢測精度高出約2%，比剩余算法的檢測精度均高出3%以上，體現出明顯的優越性。
表2
PDBselect數據庫中蛋白質數據的α螺旋檢測精度比較
Table2.
Comparison of α helix characterization accuracy of protein data in PDBselect database
4 結論
本文基于四元數描述蛋白質的局部螺旋結構，首先，利用四元數坐標系表示蛋白質相鄰氨基酸的相對旋轉，從中提取螺旋軸，并在單位球面上進行螺旋軸可視化，其次，自定義參數——QHASD，并利用該參數值對α螺旋結構進行檢測，最后，在PDBselect蛋白質數據庫的基礎上對本文算法進行驗證，最終α螺旋檢測精度為91.7%。本文在完成蛋白質α螺旋結構檢測的同時通過可視化更直觀地了解蛋白質的空間取向和分布，從而進一步掌握蛋白質的空間結構，并通過與其他蛋白質二級結構檢測算法比較，證明本文算法具有計算復雜度低，檢測精度高的優點。
深入理解蛋白質的二級結構，是理解蛋白質復雜結構和功能的先決條件。完成蛋白質α螺旋結構的檢測只是其中的一小部分，如何將本文算法應用到蛋白質二級結構檢測中是下一步需要完成的工作。
引言
蛋白質結構與功能是當前生命科學中一項重要研究內容，準確地識別蛋白質的二級結構，對其三維結構的理解至關重要[1]。膜蛋白占蛋白質分子總數的30%，其跨膜區結構大部分為α螺旋，α螺旋檢測精度的高低直接影響膜蛋白跨膜螺旋區段的定位精度，對于藥物設計中作用靶點的確定至關重要。α螺旋是蛋白質中最典型、含量最豐富的二級結構類型。對α螺旋結構的深入研究和精確描述能夠在很大程度上提高蛋白質結構檢測的準確率，并對蛋白質折疊和蛋白質設計研究有重要的推動作用。本文將重點對α螺旋結構進行研究。
蛋白質的α螺旋、β折疊等二級結構是對蛋白質局部結構的一種定性描述。檢測蛋白質的二級結構有多種方法，其中最準確的是利用核磁共振測定出每個原子的三維坐標，從而根據不同二級結構的幾何特點計算出不同的二級結構。應用最廣泛的DSSP(definition of secondary structure of proteins)算法是以主鏈中氨基之間的氫鍵模式來定義蛋白質的二級結構。對蛋白質二級結構的含量還可以用光譜來初步估計，最常用的方法是圓二色譜。利用波長范圍為170~250 nm長紫外線，在所獲得光譜吸收曲線上，α螺旋結構會在208 nm及222 nm兩處同時出現極小值，而無規卷曲和β折疊結構則分別會在204 nm和207 nm處出現單個極小值。另一個較常用的方法是紅外光譜，它可以偵測因氫鍵所造成的胺基震蕩。近幾年來，四元數在分子建模中的應用越來越廣泛[2]，涉及到的方面主要有分子的最小二乘擬合[3]、分子的平均取向和隨機取向分析，以及分子的四元數可視化[4]。
1993年，Srinivasan等[5]提出一種應用局部螺旋軸來測量螺旋規律性的方法。這種方法將兩個不同的蛋白質坐標系在四元數的基礎上進行最小二乘擬合疊加，來進行兩個不同蛋白質空間結構的對比。然而他方法中的參數Φ并不能通過簡單計算得到，而且不能直接與傳統參數Ramachandran相聯系。1999年，Quine[6]在理想蛋白質幾何形狀的螺旋結構定義中引入了四元數，但其方法只能用于理想蛋白質結構中，局限性較大，不能廣泛應用。2004年，Coutsias等[7]應用四元數來尋找蛋白質的最佳旋轉，用于兩蛋白質空間結構的配準和結構比較。2006年，Kneller等[8]利用局部螺旋運動對蛋白質二級結構的幾何情況進行描述。他們同時定義了取向距離σ、螺半徑ρ和直線度參數ξ，通過這些參數值對二級結構進行檢測。他們首次提出了straightness的概念，但是對于相鄰殘基間取向的標量描述不能與Ramachandran角直接聯系。這一年，Hanson等[9]在Kneller和Calligari的基礎上完善了straightness的定義，對蛋白質定義的坐標系可以直接與Ramachandran角相關聯，可以廣泛應用于實際蛋白質數據中。2012年，Hanson等[10]根據Frenet坐標系對蛋白質建造唯一的坐標系，轉換為四元數后在三維球體上對蛋白質的四元數數據進行可視化，構造出四元數地圖，進而區分出不同的二級結構，但并未涉及進行二級結構檢測的指標或方法。
本文在文獻[10]所提出Cα坐標系的基礎上，根據蛋白質局部螺旋結構的幾何意義計算Cα坐標系并從中提取螺旋軸向量，將螺旋軸向量進行球面可視化，提出參數四元數螺旋軸球面距離(Quaternion Helix Axis Spherical Distance，QHASD)，作為蛋白質α螺旋的檢測的指標。
1 理論基礎
1.1 四元數
由于四元數在旋轉的描述方面有著顯著的優越性，現已有許多研究者將四元數應用到生物信息學領域的蛋白質結構分析中。文獻[11]將四元數用于RNA分子研究中的分子取向分析的問題。
四元數q包括4個分量，即一個實部和三個虛部，
| $q\left( \theta ,n \right)=\left( {{q}_{0}},{{q}_{x}}i,{{q}_{y}}j,{{q}_{z}}k \right)$ |
其中，虛部單位i,j,k需滿足：
| ${{i}^{2}}={{j}^{2}}={{k}^{2}}=-1,ij=-ji=k,jk=-kj=i,ki=-ik=j$ |
四元數在表示旋轉時，n為螺旋軸向量，θ為旋轉角，則四元數又表示為
| $q\left( \theta ,n \right)=\left( \cos \left( \frac{\theta }{2} \right),n\sin \left( \frac{\theta }{2} \right) \right)$ |
1.2 四元數與旋轉矩陣的映射關系
傳統方法中通常采用旋轉矩陣描述旋轉。
1.2.1 四元數到旋轉矩陣
若已知單位四元數q(θ,n)=(q0,qxi,qyj,qzk)，可根據關系式(4)。
| $\left\lfloor \begin{matrix} {{q}_{0}}^{2}+{{q}_{x}}^{2}-{{q}_{y}}^{2}-{{q}_{z}}^{2} & 2{{q}_{x}}{{q}_{y}}-2{{q}_{0}}{{q}_{z}} & 2{{q}_{z}}{{q}_{x}}+2{{q}_{0}}{{q}_{y}} \\ 2{{q}_{x}}{{q}_{y}}+2{{q}_{0}}{{q}_{z}} & {{q}_{0}}^{2}-{{q}_{x}}^{2}+{{q}_{y}}^{2}-{{q}_{z}}^{2} & 2{{q}_{y}}{{q}_{z}}-2{{q}_{0}}{{q}_{z}} \\ 2{{q}_{x}}{{q}_{z}}-2{{q}_{0}}{{q}_{y}} & 2{{q}_{y}}{{q}_{z}}+2{{q}_{0}}{{q}_{x}} & {{q}_{0}}^{2}-{{q}_{x}}^{2}-{{q}_{y}}^{2}+{{q}_{z}}^{2} \\ \end{matrix} \right\rfloor $ |
R就是我們所要找的旋轉矩陣，其中，R為正交矩陣。
1.2.2 旋轉矩陣到四元數
對于3×3的旋轉矩陣R，可以通過一個角度θ和一個螺旋軸n表示為R(θ,n)[12]，軸角表示為
| $\left\lfloor \begin{matrix} c+{{n}_{x}}^{2}\left( 1-c \right) & {{n}_{x}}{{n}_{y}}\left( 1-c \right)-s{{n}_{z}} & {{n}_{z}}{{n}_{x}}\left( 1-c \right)+s{{n}_{y}} \\ {{n}_{x}}{{n}_{y}}\left( 1-c \right)+s{{n}_{z}} & c+{{n}_{y}}^{2}\left( 1-c \right) & {{n}_{z}}{{n}_{y}}\left( 1-c \right)-s{{n}_{x}} \\ {{n}_{x}}{{n}_{z}}\left( 1-c \right)-s{{n}_{y}} & {{n}_{y}}{{n}_{z}}\left( 1-c \right)+s{{n}_{x}} & c+{{n}_{z}}^{2}\left( 1-c \right) \\ \end{matrix} \right\rfloor $ |
其中c=cosθ,s=sinθ,n·n=1。
已知R后，可以通過如下的方法得到旋轉角θ和螺旋軸n：
| $trR=1+2\cos \theta $ |
| ${{\cos }^{2}}\left( \frac{\theta }{2} \right)=\frac{1}{4}\left( trR+1 \right)$ |
| $R-{{R}^{'}}=\left\lfloor \begin{matrix} 0 & 2{{n}_{z}}\sin \theta & 2{{n}_{y}}\sin \theta \\ 2{{n}_{z}}\sin \theta & 0 & -2{{n}_{x}}\sin \theta \\ -2{{n}_{y}}\sin \theta & 2{{n}_{x}}\sin \theta & 0 \\ \end{matrix} \right\rfloor $ |
| $R-{{R}^{'}}=\left\lfloor \begin{matrix} {{r}_{11}} & {{r}_{12}} & {{r}_{13}} \\ {{r}_{21}} & {{r}_{22}} & {{r}_{23}} \\ {{r}_{31}} & {{r}_{32}} & {{r}_{33}} \\ \end{matrix} \right\rfloor =\left\lfloor \begin{matrix} 0 & -2{{n}_{z}}\sin \theta & 2{{n}_{y}}\sin \theta \\ 2{{n}_{z}}\sin \theta & 0 & -2{{n}_{x}}\sin \theta \\ -2{{n}_{y}}\sin \theta & 2{{n}_{x}}\sin \theta & 0 \\ \end{matrix} \right\rfloor $ |
根據式(9)可以求得nx,ny,nz的值，進而得到螺旋軸n。
根據旋轉角θ和螺旋軸n，由公式(3)得到旋轉矩陣R(θ,n)對應的四元數q(θ,n)=(q0,qxi,qyj,qzk)。
2 本文算法
標準α螺旋繞固定軸旋進，該旋進可分解為繞一定點的旋轉和沿固定方向的平移，故同一α螺旋的旋轉軸為一條直線。本文利用四元數來表達蛋白質Cα坐標系間的相對旋轉，定義了反映螺旋軸在單位球面上聚類結果的參數——QHASD，根據QHASD值對α螺旋結構進行檢測。
QHASD的計算步驟：
(1) 從蛋白質空間結構數據中提取出Cα坐標系的主要原子Cα、C、N的空間坐標。 忽略掉H原子，氨基酸的骨架是一個頂端為Cα原子的四面體，另外三個頂點分別為氨基、羧基和R基。如圖 1所示為蛋白質1CRN中兩個氨基酸的基本原子框架，其中CB-SG為R基，可以看到每個氨基酸的骨架原子為N-Cα-C的組合。
圖1
氨基酸的基本結構
Figure1.
Amino acid neighboring structure
(2)在弗萊納公式(Frenet-Serret)的基礎上根據Cα、C、N原子的空間坐標構建蛋白質的Cα坐標系：
| ${{X}_{i}}=\frac{{{C}_{i}}-{{C}_{\alpha i}}}{{{C}_{i}}-{{C}_{\alpha i}}},{{U}_{i}}=\frac{{{N}_{i}}-{{C}_{\alpha i}}}{{{N}_{i}}-{{C}_{\alpha i}}},{{Z}_{i}}=\frac{{{X}_{i}}\times {{U}_{i}}}{{{X}_{i}}\times {{U}_{i}}},{{Y}_{i}}={{Z}_{i}}\times {{X}_{i}}$ |
其中，XYZ分別為列向量，因此得到3×3的單位正交矩陣構成的Cα坐標系序列[F1,F2,…,Fi,…,Fn]，其中Fi=(Xi,Yi,Zi)。弗萊納公式常用來描述粒子在連續可微的曲線上運動時曲線的切向、法向、副法方向之間的關系。每個Cα坐標系Fi均表征了對應氨基酸骨架的空間取向。
(3) 計算表示相鄰Cα坐標系之間的相對旋轉矩陣：
| ${{R}_{i}}={{F}_{i+1}}\centerdot F_{i}^{-1}$ |
其中Fi、Fi+1分別為第i個和第i+1個氨基酸的Cα坐標系，Ri表示相鄰氨基酸對應Cα坐標系之間的相對旋轉矩陣，最終得到旋轉矩陣序列[R1，R2，…，Ri，…，Rn]。
(4) 根據1.2中介紹的旋轉矩陣和四元數的對應關系，將蛋白質的旋轉矩陣序列[R1，R2，…，Ri，…，Rn]轉換為四元數序列[q1，q2，…，qn]。一個3×3的旋轉矩陣R可能被表示為q或-q兩種形式，必須選取其中之一，因此計算相鄰四元數的內積qi·qi+1，當內積為負時，將qi+1改寫為-qi+1。 (5)根據公式(3)從每個四元數中提取螺旋軸，從而可以得到對應的螺旋軸序列[n1,n2,…,ni,…,nn]，因為ni為單位向量，故將螺旋軸序列在單位球面上進行可視化，同時計算相鄰螺旋軸ni與ni+1之間的球面距離QHASD。經過多次計算驗證后設定QHASD的閾值，對蛋白質進行α螺旋檢測。
本文根據同一α螺旋中相鄰氨基酸的相對旋轉是圍繞同一螺旋軸這一幾何特征，基于四元數對蛋白質局部螺旋結構進行描述，在文獻[10]提出的Cα坐標系的基礎上進一步實現四元數螺旋軸序列的球面可視化，利用四元數螺旋軸球面距離來描述蛋白質局部結構的三維空間取向，進而檢測α螺旋。與傳統蛋白質α螺旋檢測算法相比，本文將蛋白質可視化與結構檢測相結合，在描述蛋白質中氨基酸空間取向與分布的同時，利用四元數表達局部螺旋結構進而檢測到α螺旋，該算法效率高，大大降低了計算復雜度。
3 實驗分析
本文利用實際數據對該方法進行驗證，并將本文方法得到的真實蛋白質數據的檢測結果與DSSP、KASKI等方法進行比較。
3.1 實驗及結果分析
從蛋白質結構分類數據庫SCOP(Structural Classification of Proteins)中選取以α螺旋為主的蛋白質數據1A6G，利用MATLAB與R語言混合編程對蛋白質PDB文件進行文本處理后，進行α螺旋結構檢測并與其他算法的檢測結果進行對比。SCOP數據庫是提供關于已知結構蛋白質之間的結構和進化關系的信息，所涉及的蛋白質包括結構數據庫PDB中的所有條目。PDB數據庫(Protein Data Bank)是美國Brookhaven實驗室于1971年建立的大分子蛋白質晶體結構資料數據庫，是互連網上唯一有生物高分子三維結構的數據庫，它的內容主要是根據X射線結晶和核磁共振得到的實驗數據。截止2014年2月，PDB數據庫已含有97 980個結構數據，其中約93%是蛋白質的結構數據。SCOP從總體上將蛋白質進行分類，例如全α型、全β型、以平行折疊為主的α/β型和以反平行折疊為主的α+β型等。本文選取的蛋白質1A6G為全α型。
如圖 2所示為蛋白質1A6G的三維空間結構，可以看到其中包含有6條α螺旋鏈。如圖 3所示為蛋白質1A6G的α螺旋檢測，圖 3(a)為蛋白質1A6G對應四元數序列螺旋軸的球面可視化，可以看到螺旋軸在球面有6個聚集區域，分別對應蛋白質1A6G的6條α螺旋。圖 3(b)為計算得到的相鄰螺旋軸的球面距離，有6段的球面距離明顯減少。經計算驗證，取QHASD的值為0.18時，對α螺旋的檢測精度最高，設定閾值為0.18。
圖2
蛋白質1A6G的空間結構
Figure2.
rotein structure of 1A6G
圖3
蛋白質1A6G的α螺旋檢測
(a)蛋白質1A6G對應四元數序列螺旋軸向量的球面可視化；(b)蛋白質1A6G相鄰螺旋軸球面距離
Figure3. α helix characterization of protein 1A6G(a) protein 1A6G corresponding helix axis vectors of quaternion sequences′ spherical visualization; (b) protein 1A6G′s QHASD
將蛋白質1A6G的Cα坐標系矩陣序列轉換為四元數表示，并對四元數序列中21-37段四元數的(ω,y,z)和(x,y,z)坐標進行球面可視化，根據檢測結果，該段對應的為α螺旋結構。如圖 4所示，(ω,y,z)和(x,y,z)均表示為單位球內的閉合圓環。
圖4
蛋白質1A6G的Cα坐標系對應四元數序列可視化
Figure4.
Visualization of the quaternions corresponding to the Cα frames of protein 1A6G
3.2 檢測結果與精度比較
應用標準算法DSSP和KAKSI等多種廣泛應用的二級結構檢測算法對蛋白質1A6G的α螺旋進行檢測，并與本文的四元數檢測算法結果進行比較。
蛋白質α螺旋的檢測算法中，目前應用最為廣泛的是Ramachandran(1963)[13]經典著作中提出的拉氏圖方法與DSSP算法。拉氏圖算法是提取了局部方位角信息作為參數，來描述蛋白質結構的空間構象。DSSP算法[14]是用于對蛋白質結構中的氨基酸殘基進行二級結構構象分類的標準化算法，由Wolfgang Kabsch和Chris Sander設計，該算法使用PDB格式的原子級分辨率的蛋白質三維結構坐標集數據，依靠以靜電學定義進行的氫鍵識別，以及對主鏈和側鏈二面角的計算，從而得到每個氨基酸殘基的二級結構構象參數。DSSP數據庫是由此算法生成的一個存放蛋白質二級結構分類數據的數據庫，其中包括了PDB數據庫中的所有條目。
如圖 5所示，為上述5種算法對蛋白質1A6G的α螺旋的檢測結果比較，其中α螺旋標識為“H”。
圖5
蛋白質1A6G的α螺旋檢測結果比較
Figure5.
Comparison of α helix characterization result of protein 1A6G
本文算法與其他多種經典算法對蛋白質1A6G的α螺旋的檢測精度進行了比較，具體結果如表 1所示。
表1
蛋白質1A6G的α螺旋檢測精度比較
Table1.
Comparison of α helix characterization accuracy of protein 1A6G
將本文算法與DSSP、PALSSE、STRIDE和STICKS等多種經典二級結構檢測算法對PDBselect數據庫的一個子集中蛋白質數據的α螺旋結構進行檢測，并將檢測精度進行對比。本文所使用的PDBselect數據庫是一套經實驗測定的，其序列之間的同源性小于25%的非冗余PDB數據庫的子集，共2 144種蛋白質結構序列。結果如表 2所示，與其他蛋白質二級結構檢測算法相比，本文算法比四種方法中效果最好的DSSP方法的精度略低，與效果較好STRIDE方法的相比，檢測精度高出約2%，比剩余算法的檢測精度均高出3%以上，體現出明顯的優越性。
表2
PDBselect數據庫中蛋白質數據的α螺旋檢測精度比較
Table2.
Comparison of α helix characterization accuracy of protein data in PDBselect database
4 結論
本文基于四元數描述蛋白質的局部螺旋結構，首先，利用四元數坐標系表示蛋白質相鄰氨基酸的相對旋轉，從中提取螺旋軸，并在單位球面上進行螺旋軸可視化，其次，自定義參數——QHASD，并利用該參數值對α螺旋結構進行檢測，最后，在PDBselect蛋白質數據庫的基礎上對本文算法進行驗證，最終α螺旋檢測精度為91.7%。本文在完成蛋白質α螺旋結構檢測的同時通過可視化更直觀地了解蛋白質的空間取向和分布，從而進一步掌握蛋白質的空間結構，并通過與其他蛋白質二級結構檢測算法比較，證明本文算法具有計算復雜度低，檢測精度高的優點。
深入理解蛋白質的二級結構，是理解蛋白質復雜結構和功能的先決條件。完成蛋白質α螺旋結構的檢測只是其中的一小部分，如何將本文算法應用到蛋白質二級結構檢測中是下一步需要完成的工作。
