秦雪瑞, 劉雄恩

(福建農林大學計算機與信息學院，福建 福州 350002)

分子系統(tǒng)發(fā)育分析是生物信息計算的一個重要分支，推算分子系統(tǒng)發(fā)育樹可以重建祖先序列和估計分歧時間.通過分子系統(tǒng)發(fā)育研究可以探索生命的起源和物種間的進化歷史，開展分類與區(qū)系研究以及流行病學、微生物生態(tài)學等的研究[1].

分子系統(tǒng)進化研究的第一步是建立同源性假設.DNA多序列比對數據代表了最初的同源性假設[2].一般情況下，參與比對的序列長度不是完全相同的，為了對齊需要插入間隔.比對序列的1列為1個位點，至少含有1個間隔的位點稱為間隔位點，由共同祖先分歧后發(fā)生的插入或缺失事件引起.由于間隔起源于這種特殊的突變事件，包含適合于系統(tǒng)發(fā)育分析的歷史信息[3]，因此在分子系統(tǒng)發(fā)育分析中融合間隔位點的信息是有必要的.

常用的DNA進化馬爾可夫模型，如JC69、K80、F81、F84、HKY85、TN93、REV94等，都只描述了4種核苷酸的置換(substitution)過程(本文統(tǒng)稱這類DNA進化模型為4-狀態(tài)模型)，忽略了插入/缺失事件，在分子系統(tǒng)發(fā)育分析中應用這類模型勢必會低估同源序列間的進化距離.將比對間隔視為堿基的第5種狀態(tài)，Mcguire et al[4]首次提出了包含間隔位點信息的JC69+gap、F81+gap、F84+gap等模型，但這3種改進模型將插入/缺失事件與堿基置換或顛換(transversion)同等對待.2015年林碧嬌等[5]在上述改進模型基礎上引入新的參數，進一步區(qū)分了插入/缺失與堿基置換在性質上的差異，提出JC69+gap′、 F81+gap′、F84+gap′等模型，改進后的5種狀態(tài)模型的參數較多，計算復雜，且僅在系統(tǒng)發(fā)育重建方法中的最大似然法上評估了應用效果.對于以上融合間隔位點信息的DNA進化馬爾可夫模型(統(tǒng)稱為5-狀態(tài)模型)，未在距離計算偏差上進行過有效分析.

為了在分子系統(tǒng)發(fā)育分析中盡可能多地融合indel信息，本文將多序列比對后出現的間隔視為統(tǒng)計抽樣過程中產生的隨機缺失數據.嘗試以比對多序列的p距離矩陣表示序列間親緣關系，依據最近鄰原則選擇堿基插補于特定序列的特定間隔位點，并比較分析插補前與插補后序列基于4-狀態(tài)模型及插補前序列基于5-狀態(tài)模型的序列間進化距離的大小，進而評估核苷酸最近鄰插補法的有效性.

1 研究方法

1.1 最近鄰插補的一般方法

最近鄰插補根據研究對象在輔助變量上的接近程度來選擇賦值單元，即利用輔助變量定義一個衡量單元間距離的函數，在無回答單元臨近的回答單元中，選擇與無回答單元距離最接近的回答單元所對應的值插補無回答值[6].其中，距離函數可根據應用實際采用不同的距離測度.

對于離散化矩陣，最近鄰插補法一般采用匹配度來計算樣本單元i和j之間的距離[7].記m為樣本單元個數，n為屬性類數，則樣本間的距離為：

(1)

式中，i，j= 1,2,…,m.Ai為含缺失數據的樣本單元i的屬性向量;Aj為其他各樣本單元j的屬性向量;ail為樣本單元i在屬性l上的值.要求屬性值向量Ai和Aj在屬性l上無缺失值.

(2)

可見，dij為第i個和第j個樣本單元的屬性值向量Ai和Aj中屬性不同的屬性個數.屬性個數越少，兩樣本單元距離越近.

設樣本單元i的屬性t缺失，則插補函數為

(3)

樣本單元i與k距離最小且k在t上的屬性值akt存在，則樣本單元i在t上的屬性值ait插補為akt.

由于最近鄰插補算法計算簡便，效果明顯，在缺失數據處理中都有著廣泛的應用.

1.2 同源DNA序列中間隔位點的核苷酸插補

將同源DNA比對序列視為統(tǒng)計抽樣的多個樣本單元，每個位點獨立進化，位點視為樣本單元的屬性，比對序列中的間隔即為缺失數據.在分子系統(tǒng)發(fā)育分析中，由于針對比對后的同源序列間無論使用觀察距離(即p距離[8])還是基于核苷酸替代模型的進化距離，反映的物種間親緣關系遠近的順序是一致的，而p距離通過2個序列中非同一核苷酸位點的比例來測度分歧大小，即

(4)

式中，mij和nij分別為序列i與序列j中非同一核苷酸位點數和位點總數.

基于最小進化原理[9]，以序列間p距離中最短距離作為最近鄰的依據，間隔位點核苷酸插補函數為:

(5)

即，序列Si與Sk的p距離最短且Sk在位點t上存在核苷酸Skt，則Sit的間隔插補為Skt.

同源DNA序列中間隔位點核苷酸最近鄰插補算法描述如下：

Algorithm Nucleotide Interpolation by NNI

Begin

Input multi-aligned DNA sequencesS

Computingp-distance matrixP

Fort← 1 st To the last gap site Do

Begin

Fori← 1 st To the last sequence with gap attDo

Ifpik=minj(pij) andSkt∈{A,T,C,G} ThenSit←Skt

Forj← 2 nd To the last sequence Do

IfSjt≠SitThen break Else continue loop

Ifj> count of sequences Then delete sitetElse remaint

End

OutputSafter Nucleotide Interpolation at gap sites

End

核苷酸最近鄰插補算法：先計算p距離矩陣；然后針對多序列比對的核苷酸矩陣，對含有間隔的所有列中的每個存在間隔的序列，選擇與該序列距離最近且在該位點沒有間隔的核苷酸，將其在該位點的間隔進行替換，即插補.若插補后各序列在該位點的核苷酸完全相同，則刪除該位點(整列)，否則保留插補后位點.剔除插補后核苷酸相同的位點，是因為原來的間隔位點代表可能的indel突變事件，而插補后這種位點在分子進化分析中不提供進化信息，反而會減低序列間進化距離的估算.

假設在p距離上，與序列①最近的是序列②，與序列②最近的是序列①，與序列③最近的是序列②，與序列④最近的是序列⑤，與序列⑤最近的是序列④.綠色線框標注的是最近鄰插補后核苷酸不同的位點，紅色線框標注的是插補后核苷酸相同而刪除的位點.

圖1 核苷酸最近鄰法插補示意圖Fig.1 Schematic diagram of nucleotide interpolation by nearest neighbor method

1.3 幾種DNA進化模型下的進化距離

DNA進化的馬爾可夫模型以不同狀態(tài)(4種核苷酸，或再加上1個gap狀態(tài))間置換的速率進行矩陣描述.基于這類模型可以推導出序列間進化距離(平均每個位點核苷酸置換次數)的計算公式[8].為便于測試和比較本文提出的核苷酸最近鄰插補處理方法與傳統(tǒng)的直接忽略間隔位點的方法，以及將間隔視為第5種狀態(tài)的改進模型的方法，本文采用F81、F84、F81+gap、F84+gap和F81+gap′等模型下的距離.

1.3.1 F81模型 Felsenstein[10]在JC69模型置換速率矩陣中引入4個核苷酸的比例(平衡頻率)，提出F81模型.當同源序列中4種核苷酸的比例存在偏倚，且轉換和顛換位點比例均衡時F81模型較為有效.F81模型下導出的距離為：

(6)

式中，a=2(πTπC+πAπG+πYπR)，πR=πA+πG，πY=πT+πC，πT、πC、πA和πG分別為4種核苷酸T、C、A和G的平衡頻率，計算時以觀察值估算；p為距離.

1.3.2 F84模型 Felsenstein et al[11]將核苷酸置換區(qū)分為兩類事件，類型Ⅰ僅含有轉換，類型Ⅱ既有轉換又有顛換，同時兩種類型里都有1個核苷酸都可以被相同的核苷酸置換，即核苷酸不發(fā)生改變.F84模型能較好地擬合進化過程中核苷酸的變化情況[5].該模型導出的距離為：

(7)

式中，S是轉換位點的概率，V是顛換位點的概率.顯然有p=S+V.

1.3.3 F81+gap模型 Mcguire et al[4]在F81模型中引入第5種狀態(tài)，即比對間隔，其與4種核苷酸的置換采用核苷酸之間的置換速率，該模型導出的距離，表示如下：

(8)

式中，a=2[πTπC+πAπG+πYπR+π_(1-π_)]，π_是間隔的平衡頻率.

1.3.4 F84+gap模型 Mcguire et al[4]在F84模型中同樣引入間隔狀態(tài)，將核苷酸轉換用速率α表示，而嘧啶和嘌呤的之間的顛換、4種核苷酸與間隔的置換用另一速率β表示.F84+gap模型導出的距離為：

(9)

1.3.5 F81+gap′模型 在考慮核苷酸平衡頻率因素的同時，將核苷酸之間的置換與核苷酸和間隔之間的置換(即插入/缺失)區(qū)別對待，在F81+gap′模型中，引入參數γ表示核苷酸與間隔間的置換速率.該模型導出的距離為：

(10)

式中，a=2(πTπC+πAπG+πYπR)，b=1/[a+2π_(1-π_)]，S是核苷酸置換位點的概率，I是核苷酸與間隔間置換位點的概率.顯然有，p=S+I.

2 結果與分析

2.1 測試序列

分別選取3組同源DNA序列進行測試.第1組為7種猿類物種的線粒體DNA全序列，物種及其序列GenBank檢索號分別為Pantroglodytes(NC_001643.1)、Panpaniscus(NC_001644.1)、Homosapiens(NC_012920.1)、Pongopygmaeus(NC_001646.1)、Pongoabelii(NC_002083.1)、Gorillagorilla(NC_001645.1)、Hylobateslar(NC_002082.1).用ClustalX2默認的參數進行多比對、手工優(yōu)化后，序列長度為16 644 bp，其中間隔位點為419 bp，數據缺失率為2.5%.

第2組為6屬6種睡蓮科植物的核糖體DNA中的內轉錄間隔區(qū)(ITS)序列[12]，分別為Nelumbopentapetala(AY620419.1)、Nymphaeacaerulea(AY620420.1)、Victoriacruziana(AY620423.1)、Cabombafurcata(AY620425.1)、Braseniaschreberi(AY620426.1)、Nupharlutea(AY620427.1).比對后序列長度為673 bp，其中間隔位點181 bp，數據缺失率為26.9%.

第3組為真菌側耳屬8個種的25S rDNA序列，分別為Pleurotusabieticola(AF135176.1)、Pleurotusaustralis(AF261432.1)、Pleurotuscalyptratus(AF135177.1)、Pleurotuscornucopiae(U04146.1)、Pleurotusdryinus(AF135178.1)、Pleurotusfossulatus(U04136.1)、Pleurotuspopulinus(U04159.1)、Pleurotussmithii(U04150.1).比對后序列長度為903 bp，其中間隔位點53 bp，數據缺失率為5.9%.

2.2 測試結果

表1顯示第1組數據分別在刪除間隔位點后4-狀態(tài)模型、融合間隔位點5-狀態(tài)模型和間隔插補核苷酸后4-狀態(tài)模型下的成對序列間進化距離.

表1 猿類7個物種線粒體DNA序列在幾種處理和模型下成對進化距離1)Table 1 Evolutionary distances of mitochondrial DNA sequences of 7 apes under several processings and models

1)F81(D)為刪除所有間隔位點后采用F81模型的距離，F81(NNI)為最近鄰法核苷酸插補間隔后采用F81模型的距離，F84(D)為刪除所有間隔位點后采用F84模型的距離，F84(NNI)為最近鄰法核苷酸插補間隔后采用F84模型的距離.

表2顯示第2組數據分別在刪除間隔位點后4-狀態(tài)模型、融合間隔位點5-狀態(tài)模型和間隔插補核苷酸后4-狀態(tài)模型下的成對序列間進化距離.

表3顯示第3組數據分別在刪除間隔位點后4-狀態(tài)模型、融合間隔位點5-狀態(tài)模型和間隔插補核苷酸后4-狀態(tài)模型下的成對序列間進化距離.

表1～3中的序列間平均距離以及圖2～4中針對間隔位點的不同處理或模型下估算距離的對照直觀地表明：融合間隔位點信息的5-狀態(tài)模型中的F81+gap和F84+gap的距離估算明顯偏低，改進的5-狀態(tài)模型F81+gap′、傳統(tǒng)的刪除間隔位點的處理和本文提出核苷酸最近鄰插補處理后4-狀態(tài)模型估算的距離相對接近，而改進的F81+gap′模型和核苷酸最近鄰插補處理后在4-狀態(tài)模型下估算的距離略高于直接忽略間隔位點信息在4-狀態(tài)模型下的估算，且核苷酸最近鄰插補處理方法估算的距離又略高一些.其次，序列間間隔位點數越大，忽略間隔位點方法造成的進化距離偏低估計越加突出.

表2 睡蓮科6種植物核糖體DNA中ITS序列的成對進化距離Table 2 Evolutionary distances of ITS sequences in ribosomal DNA of 6 Nymphaeaceae plants

表3 側耳屬8種真菌25S rDNA序列的成對進化距離Table 3 Evolutionary distances of 25S rDNA sequences of 8 Pleurotus fungus

由于間隔位點代表DNA突變中的核苷酸插入/缺失事件，直接刪除同源多序列比對后的間隔位點的簡單處理方法勢必導致序列間進化距離的偏低估計，應用于分子系統(tǒng)發(fā)育分析和進化樹推斷時將低估序列間距離，造成枝長偏低估計.5-狀態(tài)模型中的F81+gap和F84+gap更加低估了序列間距離，本文認為這是由于這兩個模型均沒有區(qū)分核苷酸之間的替代與核苷酸與間隔之間的置換(插入/缺失)，簡單地處理為相同性質、同一置換速率的狀態(tài)轉換過程.

圖2 不同方法估算的7種猿類線粒體DNA序列間距離對照Fig.2 Comparison of estimated distances of mitochondrial DNA sequences of 7 apes under different methods

圖3 不同方法估算的6種睡蓮科植物核糖體DNA中ITS序列間距離對照Fig.3 Comparison of estimated distances of ITS sequences in ribosomal DNA of 6 Nymphaeaceae plants

圖4 不同方法估算的8種側耳屬真菌25S rDNA序列間距離對照Fig.4 Comparison of estimated distances of 25S rDNA sequences of 8 Pleurotus fungus under different methods

睡蓮科6種植物核糖體DNA中ITS序列的成對進化距離的測試結果(表2和圖2)表明，當序列間隔位點數較多，即DNA進化過程中核苷酸插入/缺失事件的比例較高時，本文提出的最近鄰核苷酸插補方法在進化距離和進化樹枝長估算上能更為有效地消除偏低估計，對間隔位點進行核苷酸插補的處理方法使得傳統(tǒng)的4-狀態(tài)模型在序列間分歧度的估算中能夠更有效地融合DNA進化的插入/缺失信息.

改進的F81+gap′模型和核苷酸最近鄰插補處理方法至少能夠減少同源序列間距離的偏低估計.如果核苷酸最近鄰插補的方法沒有導致進化距離的偏高估計，無疑是一種有效的融合InDel信息的方法.

3 小結

鑒于分子系統(tǒng)發(fā)育重建研究中忽略多序列比對出現的間隔位點而導致低估序列間進化距離或進化樹枝長的問題，本文借鑒統(tǒng)計學中處理缺失數據的最近鄰插補法，提出一種核苷酸最近鄰插補間隔位點的處理方法.通過對3組同源DNA序列在不同的處理方法下的距離估算對照測試和上述分析，本文發(fā)現將間隔視為4種核苷酸外的第5種狀態(tài)的F81+gap和F84+gap模型不能有效融合間隔所表示的indel進化信息，反而更加低估了序列間距離，改進的同類模型F81+gap′能夠在一定程度上融合間隔所攜帶的indel信息，而本文所提出的核苷酸最近鄰插補法能夠有效運用DNA進化的4-狀態(tài)馬爾可夫模型估算進化距離，至少它能減小序列間進化距離的偏低估計，至于是否出現偏高估計還需要通過對DNA模擬進化序列進行分子系統(tǒng)發(fā)育重建和分析做出進一步判斷.