利用計算生物學方法識別原核啟動子的研究進展

蘇偉，孫自杰，岳鵬，林昊*

(1.電子科技大學生命科學與技術學院 成都 611731；2.成都東軟學院健康醫療科技學院 成都 611844)

啟動子通常位于基因上游，能與RNA 聚合酶特異性結合并起始轉錄的一段DNA 序列，作為轉錄起始過程的關鍵元件，激活RNA 聚合酶與模板DNA 結合，是基因表達和轉錄調節的起始步驟[1]。

原核生物RNA 聚合酶中的σ 因子可以特異性識別并結合啟動子。在大腸桿菌中，存在多種σ 因子，根據分子量可以分為7 類，σ70、σ54、σ38、σ32、σ28、σ24、σ19，在已知的7 類σ 因子中前6 類保守性極強，而σ19 在大多數基因組中是缺失的[2]。每一類σ 因子具有特定的生物學功能[3-6]，σ70 主要負責持家基因的轉錄；σ54 被認為是參與氮代謝的調控因子以及控制一些輔助進程；σ38 參與穩定期基因的調節；σ32 是熱休克σ 因子(熱激因子)；σ28 參與鞭毛的合成；σ24 與極端熱應激反應有關；σ19 則參與對鐵離子轉運系統的調控。根據σ 因子的同源性，可將其大致分為兩類：一類是σ70 家族，包括σ70、σ38、σ32、σ28、σ24、σ19；另一類是σ54 家族。大腸桿菌基因組內的啟動子類型依據與之結合的σ 因子種類也可分為相應的類型。不同類型的啟動子共有序列也有所差異。因此，啟動子也依據被識別的片段分為σ70 家族和σ54 家族。如σ70 啟動子具有兩個重要的基序區域，?10 區和?35 區，分別位于轉錄起始位點上游 約10bp 和35bp 處。?10 區含有保守序列“TATAAT”，又被稱為Pribnow box 或TATA box，富含腺嘌呤(adenine,A)和胸腺嘧啶(thymine,T)，有助于DNA 雙鏈解螺旋分離；?35 區則由6 個保守的核苷酸“TTGACA”組成[7]。除了σ70 因子，?10 區和?35 區也是被σ70 家族其他因子識別的重要片段。相比之下，σ54 啟動子的共有序列及其位置與σ70 啟動子具有明顯差異，在σ54 啟動子的?24 區和?12 區存在保守區域，其保守序列分別是“TGGCA[CT][GA]”和“TGC[AT][TA]”[8]。

啟動子序列的鑒定對于研究基因表達、分析基因調控機制、研究基因結構以及注釋基因信息至關重要。準確識別啟動子的方法一般是依靠昂貴且耗時費力的實驗檢測方法，然而，在全基因組范圍內進行檢測是一項艱巨的任務。隨著測序技術以及計算機技術的發展，越來越多生物的全基因組被測序出來，尤其是原核生物，因此出現了基于計算生物學的啟動子預測方法，這些預測方法在不斷地改進，有助于鑒別啟動子序列。

1 原核啟動子識別方法

原核生物RNA 聚合酶中的σ 因子可以特異性識別并結合啟動子，如圖1 所示。

圖1 大腸桿菌σ70 啟動子與RNA 聚合酶結合

2005 年至今已經開發了30 多種計算方法來預測原核生物啟動子，大致流程如圖2 所示。這些方法在許多方面有所不同，包括使用的基準數據集、特征提取方法、特征選擇技術以及分類方法等。本文總結了39 種原核啟動子預測方法，從基準數據集信息、特征表示、特征選擇、性能評估策略等多方面進行了比較和分析，如表1 所示。

圖2 基于計算方法的原核啟動子預測流程

表1 39 個原核啟動子預測工具比較

續表

39 個預測工具根據其功能可分為以下3 類。

1)普通啟動子的識別。工具1～9[9-17]屬于這一類，這些工具收集各種原核生物的啟動子作為基準數據集，包含大腸桿菌、枯草芽孢桿菌、結核桿菌、乳酸乳球菌、天藍色鏈霉菌、分枝桿菌以及假單胞菌等。并沒有指出這一類啟動子具體的類型，因此這些方法只是簡單地對啟動子序列進行預測。

2)特殊類型啟動子的預測。這一類方法包含工具10～30[18-38]。這些工具以具體類型的啟動子作為基準數據集，如大腸桿菌的6 類啟動子，原核生物的σ54 啟動子，藍細菌的5 類啟動子等。不同類型的啟動子在基因表達調控過程中起著不同且重要的作用，如目前已知的σ54 啟動子僅有數百條，而原核生物有3 萬多種，還有大量σ54 啟動子未被發現。σ54啟動子參與了氮代謝的調控，因此σ54 啟動子的預測對于了解原核生物氮代謝過程具有重要意義。

3)啟動子的預測與分類。剩余的9 個方法[39-47]均屬于這一類，以大腸桿菌啟動子作為數據集。這類方法具有一個典型的特征，即模型具有兩層結構，第一層均是對啟動子的預測，第二層是對啟動子屬性分類。工具31～36 除了預測啟動子和非啟動子，第二層還判斷啟動子的具體類型(σ70,σ54,σ38,σ32,σ28,σ24)。實際上，啟動子還有強弱之分。強啟動子能增加轉錄頻率從而提高基因的表達水平，所以預測啟動子的強度也很重要。基于此，模型37～39 的第二層鑒定啟動子的強弱(Strong,Weak)。

隨著后基因組時代的到來以及計算機的發展，對于原核啟動子的預測方法也不局限于初步的分類，還增加了對啟動子類型和強度的鑒定，為了解基因調控過程提供新信息。

2 數據集構建

建立原核啟動子預測模型的第一步需要構建一個高質量的基準數據集。大腸桿菌(E.coli)作為原核生物中被廣泛使用、研究的模式生物，其經過實驗驗證的轉錄調控信息已被系統地收錄在RegulonDB數據庫[48]中。DBTBS 數據庫[49]則收集整理了關于枯草芽孢桿菌(B.subtilis)的啟動子數據。因此，RegulonDB 和DBTBS 數據庫為預測方法提供了數據基礎。39 個工具中共有35 個工具的數據集包含大腸桿菌和枯草芽孢桿菌啟動子。

另外，為了減少由序列同源性引起的潛在誤差，通常會使用CD-HIT[50]工具以75%～85%的序列相似性閾值來去除掉數據集中序列冗余。原核啟動子相較真核啟動子，其結構相對較為簡單、功能元件也相對較少，因此一般選擇轉錄起始位點(transcriptional start site,TSS)上游60bp 以及下游20bp 作為原核啟動子序列，不僅包含了重要的共有序列，如?35 區、?10 區、起始位點等，也避免了序列過長導致引入不必要的信息，具體數據可見原核啟動子數據庫(prokaryotic promoter database,PPD)[51]。

3 特征提取

幾乎所有的機器學習方法是以數值向量作為輸入，因此需要一個合適的特征描述方法將數據集中的每一個樣本轉換為能夠反映序列信息的數值向量。在原核啟動子識別工作中，這些特征大致可以分為5 類：核苷酸組成、核苷酸理化性質、偽核苷酸組成、二進制編碼以及位置權重矩陣，以下對這5 類特征進行簡單的介紹。

3.1 核苷酸組成

核苷酸組成，也叫k-mer，統計了DNA 序列片段的所有可能組合的k長度子串出現頻率，其計算公式為：

式中，i代表某一k聯體，有4k種可能性；N(t)表示DNA 序列中某一k聯體出現的次數；L表示DNA 序列的長度。隨著k值的增加，DNA 序列的局部或短程信息也會逐漸增加。

此外，核苷酸組成還包括了g-gapped k-mer，GC 含量，累積核苷酸頻率(accumulated nucleotide frequency,ANF)等。ANF 表示了每一個堿基在序列中的分布密度，表達式為：

式中，|si|代表第i個堿基的位置；N(si)表示某一堿基出現頻數；q∈{A,C,G,T}。

3.2 理化性質

DNA 序列中堿基的理化性質也可作為啟動子預測的重要特征，包括核苷酸的化學性質、雙鏈的穩定性、自由能、應激誘導的DNA 雙鏈不穩定性等。

根據表2 中對不同核苷酸的分類，DNA 序列中第i個核苷酸可以表示為：

表2 核苷酸化學性質

式中，xi,yi,zi分別表示指環結構(ring structure)，功能組別(function group)，以及氫鍵(hydrogen bond)，如：

因此4 種堿基(A,C,G,T)可以分別表示為(1,1,1)，(0,1,0)，(1,0,0)和(0,0,1)。

3.3 偽核苷酸組成

偽核苷酸組成(pseudo k-tuple nucleotide composition,PseKNC)最初是由文獻[52]提出，分為I 型和II 型。這兩種方法基于核苷酸的物化性質引入了DNA 序列的全局或長程順序信息。

I 型PseKNC，也叫平行相關偽核苷酸組成，將每一條DNA 序列轉化為4k+λ 維的向量，具體表示為：

II 型PseKNC，也叫串聯相關偽核苷酸組成，可產生4k+λ Λ維向量：

式(5)和式(6)中的fu與式(1)相同；前4k個元素是核苷酸組成特征，后面的元素是偽核苷酸組成特征；λ是一個正整數，反映序列順序關聯階數；ω是權重因子，用于權衡核苷酸組分和DNA 序列局部結構性質的影響；τj代表的是m階關聯因子，反映了每條DNA 序列所有二核苷酸的m階順序關聯性。

3.4 二進制編碼

二進制編碼通過將4 種核苷酸轉換成包含4 個元素的向量作為特征，其中一個元素為1，其余為0，既A、C、G 和T 分別表示為(1,0,0,0)，(0,1,0,0)，(0,0,1,0)以及(0,0,0,1)。因此，一段長為L的DNA序列可以用L×4 的二維矩陣表示。

3.5 位置權重矩陣

位置權重矩陣(position weight matrix,PWM)可用來表示序列的保守片段，以序列每一位置的堿基保守程度為參量，分別計算每種堿基的保守指數，以此作為特征，具體表示為：

式中，Si,j表示堿基i在第j個位置的保守指數；qi,j是指在背景序列中堿基i出現在第j個位置的頻率；bi是背景概率。

因此，PWM 是一個4×L的二維矩陣：

4 特征選擇

從式(1)以及式(5)、式(6)可以看出，隨著k值的增加，特征維度呈指數級增長，會導致“維度災難”以及過擬合問題，而且由不同特征提取方法整合形成的融合特征集合往往會夾雜一些冗余或不相關的信息，所以為了避免出現上述問題并且提高計算效率，篩選有用的特征也是必不可少的步驟。

4.1 最小冗余最大相關

最小冗余最大相關(minimum redundancy maximum relevance,mRMR)[53]是一種通過篩選相關性最大的特征來減少信息冗余的方法。mRMR 的應用大大減少了特征維數和模型訓練的時間，幾乎不丟失有效信息。

對于兩個隨機變量x和y，其互信息為：

式中，p()表示概率密度函數。

最大相關性為：

式中，c為類別變量；S為特征子集。

最小冗余度則表示為：

最后的評選標準如式(12)所示：

mRMR 會將所有特征的最大相關最小冗余打分按從大從小排序，值越大表明該特征越重要。

4.2 最大相關最大距離

當兩個特征高度依賴時，它們對模型的貢獻不能疊加，文獻[54]基于距離函數提出了最大相關最大距離(max-relevance-max-distance,MRMD)來衡量每個特征的獨立性。

MRMD 包含兩個方面的特征排序度量：1)特征子集與目標類別的相關性；2)特征子集的冗余度。采用皮爾遜相關系數來衡量相關性、多種距離函數來計算冗余度。皮爾遜相關系數越大，特征與目標類別之間的相關性越高；特征距離越大，特征子集的冗余度越低；相關性與距離之和大的特征被選入最終的特征子集。因此，MRMD 生成的特征子集冗余度最低，與目標類別的相關性最強。

4.3 F-score

F-score 是一種基于filter 的特征選擇方法，對每一個特征進行重要性打分，其具體計算方法為：

式中，n+、n?分別表示正負樣本的數量；分別指第i個特征在正樣本、負樣本以及所有樣本中的平均值；分別指的是正負樣本中第k條序列的第i個特征的數值。

F-score 通常與增量特征選擇技術相結合來確定最優特征子集。

4.4 增量特征選擇

增量特征選擇(incremental feature selection,IFS)方法適用于確定最優特征子集。該方法的核心思想是將按重要性評分降序的特征依次加入到特征子集中，形成新的子集，將每一個子集輸入至模型中，從而根據結果決策出最優特征子集。

5 分類方法

選擇合適的算法可以使最終的模型具有良好的性能和泛化能力，各種監督學習方法已經被廣泛應用于預測原核啟動子，大致有以下4 類。

5.1 支持向量機

支持向量機(support vector machine,SVM)[55]是基于監督學習方式對數據進行二元分類，在樣本空間中尋找最優分類超平面使得兩類的間隔最大。

對于線性可分的情況，存在一個分類超平面能將訓練樣本正確分類。而對于線性不可分的情況，需要使用核函數將低維不可分樣本映射到更高維的特征空間，使得樣本在高維空間中線性可分。

5.2 神經網絡

神經網絡(neural networks,NN)學習是一種模擬生物大腦神經網絡的自適應計算模型。隨著近年來人工智能的快速發展，人工神經網絡(artificial neural network,ANN)及其卷積神經網絡(convolutional neural network,CNN)已成為研究生物信息學問題的重要方法。

基本的ANN 結構包括輸入層、隱藏層和輸出層，主要特點是信號正向傳播，誤差反向傳播。通過最小化誤差函數，修正神經元間的連接權重，當其誤差小于一定閾值的時候，即停止訓練。

CNN 目前在很多研究領域都取得了巨大的成功，如語音識別、圖像識別、自然語言處理等，是深度學習的代表算法之一。CNN 通常由輸入層、卷積層、激活函數、池化層、全連接層和輸出層組成。與傳統的神經網絡不同的是CNN 采用局部連接和權值共享，使得網絡易于優化并且降低了模型的復雜度，減小過擬合風險。

5.3 集成學習

集成學習(ensemble learning,EL)通過構建并結合多個學習器來完成學習任務。在預測原核啟動子的方法中，集成學習也是被廣泛應用的，如隨機森林(random forest,RF)。

RF 是一種基于決策樹的集成學習方法，在決策樹的訓練過程中引入了隨機屬性選擇。對于基決策樹的每個結點，隨機選擇該結點屬性集合中的一個子集，再從這個子集中選擇一個最優屬性用于劃分。RF 的每一個決策樹都會產生一個分類結果，通過投票決定最終輸出。與單一的決策樹相比，RF 具有較強的魯棒性，并且對大數據具有較好的處理效果。

5.4 線性判別分析

線性判別分析(linear discriminant analysis,LDA)在二分類問題上最初是由文獻[56]提出的，亦稱為“Fisher 判別分析”。

LDA 的核心思想相對簡單：首先將訓練集中的樣本投影到一條直線上，使得同一類樣本盡可能靠近，不同類樣本盡可能遠離；當新樣本進來時，將其投影到同一直線上，從而根據投影點的位置判斷其類別。

6 性能評估

在統計分析中，獨立測試集和K 折疊交叉驗證已經被廣泛地應用于驗證分類器性能。當樣本數量足夠多時，會將基準數據集劃分為訓練集和獨立測試集。獨立測試集由于未參與模型的訓練，可以更好地評價模型性能。在原核啟動子識別模型中，K 折疊交叉驗證的應用最為廣泛，其基本思想是重復利用數據，每一個樣本既可以作為訓練集參與模型訓練，也會作為測試集參與模型評估。方法是將數據平均分成K份，K?1 個子集用作訓練，剩余一份用作測試，重復K次，最后返回K次結果的平均值。K 折疊交叉驗證最大程度上利用了每一個數據，能更好地反應模型的預測性能。

另外，受試者工作特征曲線(receiver operating characteristic curve,ROC)下面積AUC 值也可以反應模型性能，其值越接近于1，表明模型性能越好。

7 結束語

近年來，基于生物信息學的原核啟動子預測方法備受學者關注，已有多種方法被提出。為了充分了解這個領域的發展現狀，本文收集并系統地分析了2005 年至今共計39 個原核啟動子預測方法，詳細闡述了這些方法的數據集構建、特征選擇、特征提取、分類算法以及性能評估，詳細信息如表1所示。

目前，對原核啟動子預測的研究取得了令人滿意的結果。隨著更多原核生物的基因組被測序出來，被研究的物種也不局限于少數幾個模式生物，使用這些預測算法有助于了解原核生物基因調控機制。本文系統地比較了原核啟動子預測方法，為研究此問題提供新思路、新角度。