彭寶成 張曉煒 劉 暘 樊國梁**

（1）內蒙古大學物理科學與技術學院，呼和浩特 010021；2）內蒙古醫(yī)科大學第一附屬醫(yī)院風濕免疫科，呼和浩特 010050）

啟動子是基因序列中一段可以調控基因表達的核苷酸序列（序列中含有起始位點），控制著基因的表達與否，因此啟動子在基因的轉錄和表達中具有重要的地位。

在大腸桿菌基因組中，啟動子是一段分布于起始位點上游60 bp及其下游20 bp、包含起始位點的長度為81 bp的DNA堿基序列。按照Sigma因子的類別，大腸桿菌共有7種啟動子，分別是Sigma19、Sigma24、 Sigma28、 Sigma32、 Sigma38、Sigma54、Sigma70。由于啟動子序列具有保守性，根據(jù)保守性片段區(qū)域的不同，啟動子又可以分為保守性區(qū)域在轉錄起始位點上游-10 至-35 位點附近的（以Sigma70 為代表）和在轉錄起始位點上游-12 至-24 位點附近的（以Sigma54 為代表）保守性啟動子［1-2］。

大腸桿菌的分布極廣，并且作為腸桿菌類的成員經常被作為細菌模式生物而被廣泛研究，因此人類對于大腸桿菌的研究非常深入（生物實驗方面）。大腸桿菌所有系列的基因序列已經在20 世紀被完全測量，但是基于生物實驗方法尋找啟動子的方式十分耗時、昂貴。雖然它可以較為準確地定位啟動子序列，但是在面對海量數(shù)據(jù)時，效率低的弊端開始凸顯，因此計算生物學的研究應運而生。以往的研究者提出了各種各樣的模型并且取得了良好的預測效果［3-7］。2015年，丁輝等構建的三聯(lián)體位置關聯(lián)矩陣的預測方法，預測Sigma54 的精度達到82.0%［8-9］；2015年閆妍等［10］利用位點特異性打分矩陣（positive-specific scoring matrices，PSSM）方法預測Sigma 啟動子，模型對Sigma54 的預測準確率為97.4%，對Sigma38 的預測準確率為96.0%，對Sigma70 的預測準確率為74.0%。此外還有基于柔性參數(shù)+二聯(lián)體位置關聯(lián)權重矩陣［11］、位置關聯(lián)打 分 特 征（position-correlation scoring feature，PCSF）+ 偽核苷酸特征（pseudo k-tuple nucleotide composition，PseKNC）［12］、啟動子元件的位置權重信息+k-聯(lián)體核苷酸頻率等方法，后來的趨勢表明，研究者趨向于組合多種特征來定義啟動子序列［13］，以至于分類特征維數(shù)急劇增加，因此不得不在分類時對數(shù)據(jù)進行降維處理。在算法方面，常用的機器學習算法有支持向量機（support vector machine，SⅤM）、隨機森林、K-近鄰、隱馬爾科夫、人工神經網絡、前向傳播算法等或其他算法如線性判別分析、二次判別分析等［14-19］。近些年深度學習算法也逐漸被研究者所關注，并且已有研究者將卷積神經網絡（convolutional neural networks，CNN）應用于大腸桿菌啟動子Sigma70和非啟動子的二分類預測中，對由A（1,0,0,0）、T（0,1,0,0）、G（0,0,1,0）和C（0,0,0,1）編碼的向量序列進行預測得到了敏感性Sn（sensitivity）=0.90、特 異 性Sp（specificity） =0.96、 準 確 率Acc（accuracy）=0.84 的結果，但準確率仍需進一步提高，編碼方式仍需優(yōu)化［20］。

1 數(shù)據(jù)集的構建和特征算法

1.1 數(shù)據(jù)的選取

論文中采用的Sigma38、Sigma54、Sigma70數(shù)據(jù)集可以從RegulonDB 10.8（http://regulondb.ccg.unam.mx/Downloadable）下載，為了便于后續(xù)數(shù)據(jù)處理，保證數(shù)據(jù)的一致性，對數(shù)據(jù)文件做如下處理：

a.剔除Sigma Factor that recognize the promoter標注為多種類型的行。

b.刪除文件信息（以“#”開頭的行）和序列信息列（如“acrDp2”等信息），只留下序列所在列。

c.刪除沒有序列內容的行。

經過處理后得到了Sigma38 序列146 條、Sigma54 序 列96 條、Sigma70 序 列810 條，共 計1 052條，作為正集。啟動子序列的長度均為81 bp（-60 bp～20 bp，設基因轉錄起始位點為0）。

分別選取了300 條編碼區(qū)序列和300 條基因間序列，長度均為81 bp，作為負集?；蜷g序列處于兩條基因之間的區(qū)域，選取中央區(qū)域可以最大限度的避免選入啟動子序列。

1.2 特征描述

大腸桿菌基因序列為字符序列，而CNN 算法（目前為止的大部分算法）能處理的只有數(shù)值序列，因此需要將字符序列轉換為數(shù)值序列，這一步也被稱之為特征提?。ɑ虿蓸樱?1-22］。在轉換過程中丟失的信息越少則對序列的描述也就越精確。首先，每一類序列的訓練集分別構建了位點特異性打分矩陣。其次，利用該PSSM對該種類的訓練集和測試集打分，從而將字母序列轉變?yōu)閿?shù)值序列。

1.2.1 構建位點特異性打分矩陣特征

a.構建頻數(shù)矩陣

將序列坐標化（從0 ～80共81個位點），統(tǒng)計每個坐標上A、G、C、T 4種元素出現(xiàn)的頻數(shù)，如果存在坐標點上面某種核苷酸出現(xiàn)的次數(shù)為0，則需要引入偽計數(shù)（即將該位點上所有類型核苷酸的頻數(shù)加一整數(shù)，本文中加1）。最終形成的頻數(shù)矩陣列名為序列坐標（s1～s81），行名為4 種核苷酸（A、G、C、T）。

b.生成頻率矩陣（偽計數(shù)頻率矩陣）P

計算每種核苷酸在該位點上出現(xiàn)的頻率，即用該位點上每一種元素的頻數(shù)除以該位點上所有元素的頻數(shù)和，其中xi,j是在i位點上第j種元素的頻數(shù)，Pi，j是在i位點上第j種元素的頻率：

c.生成對數(shù)幾率比矩陣oddratio：

其中P(x|M)為核苷酸出現(xiàn)的實際概率，P(x|R)為核苷酸出現(xiàn)的隨機概率（即0.25），因此將頻率矩陣P中的每個頻率值除以0.25 即可得到oddratio。隨后對oddratio求以2 為底的對數(shù)，該矩陣即為位點特異性打分矩陣，矩陣的大小為（4×81）維。

1.2.2 對特定序列打分

從PSSM 矩陣中查出特定位點上核苷酸的分值，對給定序列賦值。實現(xiàn)方式是矩陣點乘，最終可以得到1×81或者4×1×81的矩陣（即一維矩陣或者四維矩陣）（圖1）。

Fig.1 Sigma70 sequence matrix pixel matrix bitmap after scoring by PSSM

這里的維度指的是通道數(shù)：RGB圖片有3個通道（R、G、B）即通道數(shù)為3，每個維度上的每個位點的數(shù)值即為該通道在該位點上的通道顏色濃度，取值范圍為0～255，一張RGB圖片數(shù)值化后就是一個三維數(shù)值矩陣（圖2）。

類似的，可以把A、G、C、T 作為4 個通道，將序列數(shù)值化的方法就變得容易理解：一條序列在A通道（A通道長度為81）上每個位點的值可以是0（該位點不為A）或者1（該位點為A），對其他3個通道做同樣處理，結果如圖3c 所示。計算證實訓練上述方法得到的01 數(shù)值序列也可以獲得好的效果，但是01 離散數(shù)字序列攜帶的只有位置分布信息，信息量過小，因此過擬合的現(xiàn)象常常發(fā)生。

把A 通道中的1 替換成PSSM 中的打分值。在位點s1上，“g”的打分值為-0.497 5，那么G通道中位點s1 上的1 就可以被替換為-0.497 5，如果是0則不需要替換（表1）。這種方法顯然可以讓數(shù)值化序列帶有更多的信息，對其他通道做同樣處理，就構造出了4個維度的數(shù)值矩陣（圖3d）。

Table 1 Sigma70 position-specific scoring matrix（PSSM）for a certain training set

圖片上每個色塊的顏色是由3個通道顏色組合而成的，與此不同的是，啟動子序列的每個位點只有1個通道的值（AGCT中的一個）（圖3b），因此每個通道中依舊有很多的0 存在（圖3d），所以可以將4個通道對應位點的數(shù)值加和，得到一個沒有‘0’的長度為81 的向量（圖3e），將向量變換成9×9的矩陣繪圖以后形成的點陣圖片（圖1）。依舊是由于序列和圖片的不同點，序列的四維數(shù)值信息和一維數(shù)值信息沒有本質上的區(qū)別，反映在模型表現(xiàn)上就是預測的準確性沒有大的差別。

Fig.2 Picture of RGB three-channel schematic diagram

Fig.3 Overview of conversion method for certain Sigma70 sequence

與圖片信息相同的是，啟動子的保守性區(qū)域可以類比為被人類所理解的圖像信息（如數(shù)字“1”反映在圖像上是幾個白灰色像素塊組合，圖4），無論是保守性序列還是數(shù)字“1”都具有不變性，這為研究者們尋找啟動子的獨有特征奠定了理論基礎［23］。

Fig.4 Handwritten numbers “1” from http://deeplearning.net/data/mnist/

1.3 分類算法

CNN 是一個兼顧全連接和卷積取樣的前饋神經網絡算法，它的提出來源于對生物的視覺皮層研究。它在處理多層網格結構數(shù)據(jù)方面具有巨大優(yōu)勢，因此現(xiàn)在CNN 被廣泛應用于圖像識別、語音識別和自然語言處理等領域。CNN 的實現(xiàn)分為兩個主要步驟：a.數(shù)據(jù)的預處理：將圖像等信息數(shù)字化和去噪等操作；b.特征提取和分類：這一部分由CNN 網絡結構來實現(xiàn)，基礎的網絡結構包括卷積層、池化層和全連接層3 部分，深度CNN 是基礎網絡結構的多層疊加，這樣可以實現(xiàn)更多特征的提取從而提高識別精確度。

CNN 的主要特點有局部區(qū)域連接、權值共享和降采樣。a.局部區(qū)域連接（圖5 左），即前后兩層網絡的所有神經元并不是都互相連接的，目的是為了模擬視覺神經的選擇性聚焦，這有利于減少訓練參數(shù)；b.權重共享，即一個卷積核在提取特征時的權重在整張圖片上都相同，這意味著一個卷積核只能提取一種特征；c.降采樣，主要有最大值、平均值（圖5右）和方均值保留等方式，通過池化層（pooling layer）實現(xiàn)，目的是為了降低特征分辨率和減少過擬合風險。由于這3個特點，CNN在抽取特征時更能夠聚焦圖片的重要信息（圖4中黑色的背景顯然就不是有效信息，因此CNN 就不會著重關注）。由于啟動子序列的類圖片特點，將CNN作為特征提取和分類的算法是比較合適的選擇。

Fig.5 Fully connected layer and average pooling layer

2 檢驗方法

2.1 模型評價

2.1.1 自洽驗證

分別從各類序列中抽取96 條序列組成均衡樣本；構建均衡樣本中每種序列的位點特異性打分矩陣（1.2 節(jié)），并按照1.2 節(jié)方式分別使用各自的位點特異性打分矩陣將各自的所有序列（包括均衡樣本序列）轉變成數(shù)值序列。以均衡樣本序列作為訓練集，采用全部序列測試模型性能。

2.1.2 獨立檢驗

所有序列被分成兩部分（數(shù)量均等），一份用做訓練集，另一部分用作測試集，用訓練集的位點特異性打分矩陣將測試集和訓練集轉變?yōu)閿?shù)值化序列，使用訓練集訓練模型，測試集測試模型性能。

2.1.3 十交叉檢驗

該方法與前述獨立檢驗不同的是，構建位點特異性打分矩陣所使用的序列數(shù)目不同。十交叉檢驗隨機將數(shù)據(jù)集分成10份，取其中的9份構建位點特異性打分矩陣并被轉換成數(shù)值序列后用作訓練集，剩余1 份被轉換成數(shù)值序列后用作測試模型性能。如此循環(huán)，則共產生了10 份測試結果，對10 份結果取平均就得到了最終的驗證結果。

2.2 驗證參數(shù)

2.2.1 ROC曲線

ROC 曲線即受試者操作特性曲線（receiver operating characteristic curve，常稱ROC 曲線），描述的是在不同的標準下（主要是閾值），模型的擊中率和誤報率之間的函數(shù)關系。ROC 曲線常用在二分類模型當中，但是在多分類模型中也可以采取此參數(shù)（將被求ROC 參數(shù)的樣本設為正集，其余種類為負集）。

曲線下面積（area under curve，AUC）一般在0.5～1 之間，AUC 越大即代表模型的性能越好。ROC 曲線的優(yōu)勢在于當正負樣本數(shù)量對比發(fā)生變化時候曲線的形狀不會變化。

2.2.2Sn、Sp、Acc

Sn、Sp常用在二分類模型驗證當中，Sn表示正確預測正集樣本的概率（式3），Sp表示正確預測負集樣本的概率（式4）。

Acc無論是在多分類還是二分類都比較常用的參數(shù)，它表示樣本被正確預測的概率：

每個種類序列的準確率Acc公式為：

在Acc的計算公式中，Tn 和Tp 分別表示被預測成功的負集樣本數(shù)和正集樣本數(shù)，F(xiàn)p和Fn則表示被預測失敗的負集樣本數(shù)和正集樣本數(shù)。

在Acc的計算公式中，T代表該類序列中被預測正確的數(shù)目，P代表該序列中被預測失敗的數(shù)目。

3 結果評估

本文采用PSSM矩陣作為識別序列的特征，用打分的方式將字符序列轉變成數(shù)值序列。基于實際需要，本文做兩組四分類：Sigma 序列和編碼（Coding） 序 列、Sigma 序 列 和 非 編 碼（Noncoding）序列。

3.1 損失函數(shù)

模型訓練中使用的是交叉熵損失：

p(x)代表真實概率分布，q(x)代表預測概率分布，交叉熵損失是對兩個概率分布差距的評估。交叉熵的函數(shù)圖像（以Sigmod 函數(shù)作為激活函數(shù)）（圖6），可以看到交叉熵圖像具有單調性并且損失越大梯度越大，因此訓練時權重可以很好地進行更新。

3.2 自洽檢驗預測效果

自洽檢驗的目的是為了證明模型的合理性，由于樣本不均衡可能帶來收斂難度增大等問題，數(shù)量過多的樣本可以被人為地隨機丟棄一部分，使要進行分類的各種序列數(shù)量保持一致（即達到樣本均衡），并且將分布均勻的樣本作為訓練集，測試對全部樣本的預測能力，以此來檢驗均衡樣本對總體樣本的預測能力（圖7）。

Fig.6 Cross_Entropy curve

Fig.7 Self-consistent verification training curve

對樣本參數(shù)進行評估，可以看到整體的預測性能比較好，各項參數(shù)都有比較好的表現(xiàn)（表2）。

Table 2 Self-consistent verification overall evaluation Parameters

圖8為兩種預測模型的ROC 曲線，每個種類模型的AUC 值均在0.96 以上，模型表現(xiàn)出了良好的分類能力。

Fig.8 Self-consistent verification of two four-category ROC curves

需要注意的是，均衡樣本的優(yōu)勢在于可以快速收斂并且可以有效地防止過擬合現(xiàn)象的發(fā)生，它的訓練次數(shù)較少并且學習率較高，使用的分類器較少，可以很好地減少內存損耗，預測精度也十分可觀。人為制造均衡樣本的缺點在于訓練集可能不能完全包含樣本的所有特征從而影響預測準確率，為此通過對量大的樣本進行多次采樣，然后生成多個分類器，最后對所有分類器結果求平均可以解決此問題。

3.3 獨立檢驗

從圖9中可以看到獨立檢驗的訓練取得了良好的效果，損失函數(shù)曲線以及準確率曲線都比較平

Fig.9 Independent inspection training curve滑，這說明訓練的過程比較順利?？梢钥吹綔蚀_率曲線最終上升到了1.0，損失函數(shù)下降到了一個較低的位置（0.000 6）。

第999 次的準確率為tensor（1）；第999 次的損 失 函 數(shù) 為： tensor （0.000 6， grad_fn=＜NegBackward＞）。

下面的表中展示了獨立檢驗的驗證結果（平均結果），With-coding 表示Sigma 和Coding 序列的驗證 結 果（表3），With-Noncoding 表 示Sigma 和Noncoding序列的驗證結果（表4）。

Table 3 Independent inspection result（With-coding）

Table 4 Independent inspection result（With-noncoding）

表3、4 顯示Sigma38 序列相對于其他啟動子的準確率明顯較低，這可能是由于過擬合導致的，在訓練中Sigma70也較容易出現(xiàn)這種情況。

獨立檢驗繪出的ROC曲線（圖10、11），可以看到無論是哪種樣本組合，曲線的表現(xiàn)都十分讓人滿意。相對地來說，Sigma38的結果較差一些，這也是和上面的Acc結果有很好的對應。

Fig.10 Independent inspection ROC curve(With-coding)

Fig.11 Independent inspection ROC curve(With-noncoding)

獨立檢驗的結果證明（表5），訓練非均衡樣本依然是可行的途徑，這主要歸功于CNN 分類算法的高精度和PSSM 能夠較好地代表序列的特征。但是非均衡樣本會帶來收斂過慢甚至出現(xiàn)無法收斂的災難，因此常常需要對其進行過采樣處理（采用較低的學習率和較多的訓練次數(shù)、更多的分類器以提取更多特征，如本次卷積層使用了16×16甚至是20×20 的卷積核，而均衡樣本使用的卷積核為5×5），這樣就又增大了過擬合的風險，因此訓練成功的難度也相較于均衡樣本高，并且它對算力的損耗也增大了。

Table 5 Independent inspection overall evaluation parameters

3.4 十交叉檢驗

為了節(jié)省算力開銷，實驗結論通過十交叉檢驗法被驗證（表6～8）。

Table 6 Ten-fold inspection result(With-coding)

Table 7 Ten-fold inspection result(With-noncoding)

Table 8 Ten-fold inspection overall evaluation parameters

繪出的ROC曲線如下：

啟動子和Coding 區(qū)序列四分類的ROC 曲線（圖12）均在對角線上方，并且AUC值均比較接近1（Sigma38 的AUC 值為0.97，Sigma54 和Coding 序列的AUC 均為0.99，Sigma70 為0.96），這表明模型對于啟動子的分類效果比較理想。

Fig.12 Ten-fold inspection(With_coding)ROC curve

在啟動子和Non_coding 區(qū)序列四分類的ROC曲線（圖13）中，Sigma54序列和Non_coding區(qū)序列依舊達到了0.99，Sigma38 的AUC 值上升到了0.98。

Fig.13 Ten-fold inspection(With_noncoding)ROC curve

在十交叉檢驗的結果中，可以看到模型對于四分類的整體預測準確性都達到了0.97以上，并且對每一種序列的預測精確性也都達到了0.95以上。

3.5 對比分析

采用PSSM特征和采用二聯(lián)體+柔性參數(shù)［11］兩種方法得到的準確性（表9）結果的對比數(shù)據(jù)顯示：PSSM特征的預測效果更為理想，并且采用本論文中方法取得的效果遠遠好于二聯(lián)體+柔性參數(shù)的方法，這說明PSSM特征對于序列特征的描述更為精確。

單獨使用PSSM 特征分類和采用PSSM 特征+CNN 算法分類兩種方式對Sigma38、Sigma54 和Sigma70 的預測結果（表9）的對比數(shù)據(jù)顯示：CNN算法對3種啟動子的預測準確性都優(yōu)于僅僅使用PSSM 特征分類的方式，而且CNN 算法對每種啟動子的預測準確性都比較均衡，沒有出現(xiàn)對某個啟動子預測精度過小的現(xiàn)象.

為了更進一步探究算法對啟動子和非啟動子的分類效果，兩種方法（表10）對相同的數(shù)據(jù)集預測并且進行了十交叉驗證，對比結果顯示：本論文中的方法取得了更理想的結果。本論文方法在不同分類條件下得出的結果（表9，10）對比顯示：二分類的效果要好于多分類。

本論文取得的成果與同行的最新研究成果對比顯示（表11）：Grad-CAM 編碼方法（Feature by Grad-CAM）可以取得稍好的準確率，但其特異性尤其是AUC 值表現(xiàn)較為遜色，說明其模型穩(wěn)定性可能存在問題，本論文則兼顧了準確率和模型參數(shù)兩方面。獨熱編碼（one-hot encoding）方法取得的準確率為0.901，AUC值為0.957 2，本論文的結果相對較好［24-25］。

Table 9 Comparison of Acc（Ten-fold inspection）

Table 10 Comparison of comprehensive parameters（Ten-fold inspection）

Table 11 Comparison of comprehensive parameters(with the newest results)

4 討 論

CNN+PSSM 方法采用的特征簡單易用，并且多分類可以大幅提高預測效率。有研究者單獨采用PSSM打分進行分類，這種方式取得的效果稍遜于本文方法，主要原因可能是PSSM 特征稍顯簡單。在沒有有效算法輔助的情況下，這種分類方式相對利用更復雜、覆蓋差異更全面的特征描述方法表現(xiàn)確實遜色。但是這并不意味著利用簡單的特征描述方法得不到好的結果。Shujaat等［24］和Zhang等［25］都利用較簡單的特征描述（01 序列為基礎）得到了較好的預測結果（前者97%以上，后者80%左右）。其次是在分類算法上，Shujaat 等［24］采用了過于簡單的01序列造成了很大的泛化，本身01序列蘊含信息就比較少，再次提取特征只會讓CNN模型的過擬合風險增大。Zhang 等［25］雖然取得了較好的準確率，但是在AUC（僅在0.84 以上，且Sigma38的AUC為0.63）和Sp（0.78以下）等模型評估參數(shù)的表現(xiàn)上不如人意，其采用的序列轉換方法（由字符序列轉換為數(shù)字序列的方法）也是01編碼方式。

樣本不均衡容易造成訓練的泛化，目前還沒有好的調參辦法來完美解決這一問題，并且這一點也常被同行研究人員所忽視。目前大多數(shù)深度學習的研究者采用的方法是減少訓練樣本中數(shù)量過多的種類的數(shù)量，人為地調整樣本分布，或者采用過采樣的方式對數(shù)量較少的種類重復取樣。有研究者提出過采樣可能是在以ROC/AUC作為評價指標時最佳的處理方式［26-27］。為了減少模型擬合困難問題的發(fā)生，本文采用了人為調整樣本數(shù)量的方法，剔除了過多的樣本。

CNN 和PSSM 特征結合是采取了兩條腿走路的方法：選取良好的特征描述方法可以最大限度的覆蓋不同種類啟動子之間的差異，為CNN 提取差異進而進行有效的分類奠定基礎；CNN 自學習的優(yōu)良特性和對算法的優(yōu)化也是提升模型性能的關鍵。由于PSSM 特征對啟動子序列描述較為全面，因此本方法在序列轉換過程中丟失了更少的信息，同理如果可以選擇特征描述更為全面的轉換方法，模型的準確率會有進一步的提升。

5 結 論

本論文通過構建樣本的位點特異性打分矩陣，并且使用樣本的位點特異性打分矩陣將待預測的字符序列轉化成數(shù)值序列。利用PSSM特征訓練出來的CNN 模型對Sigma38、Sigma54 和Sigam70 3 種序列進行預測，分別得到了0.978 9、0.995、0.964 4的預測精確度。

在將序列數(shù)值化的過程中，PSSM特征能夠很好地表征每種序列的核苷酸分布信息，這使得每種序列之間的區(qū)分度比較明顯。由于PSSM特征的構建方法簡單，因此該特征對于序列的特征表述不會過于冗雜，這有效地降低了CNN 在訓練模型時發(fā)生過擬合的風險。本文提出了一種解釋序列特征的新思路，即利用類圖片特征構建序列的點陣圖像，這為下一步研究序列特征提取提供了一個新的方向：例如，如果可以基于PSSM再創(chuàng)造一套多通道的標準（類似于RGB標準），讓每個位點的數(shù)值由多個通道共同決定，那么將序列展開為多維矩陣的識別效果可能更好［28］。