譚生龍

摘要：測序技術(shù)的快速進(jìn)步產(chǎn)出了大量生物序列，DNA序列是生物大數(shù)據(jù)的重要組成部分，僅有極少部分DNA序列已通過實(shí)驗(yàn)驗(yàn)證了功能；通過機(jī)器學(xué)習(xí)方法快速預(yù)測DNA序列的功能是確實(shí)可行的途徑。本文探討了將DNA序列轉(zhuǎn)化為特征向量的方法，并使用機(jī)器學(xué)習(xí)方法對未知功能序列進(jìn)行功能標(biāo)注一般步驟。

關(guān)鍵詞：DNA序列；特征提??；功能預(yù)測；序列數(shù)據(jù)庫

中圖分類號：TP18 文獻(xiàn)標(biāo)識碼：A 文章編號：1009-3044（2016）25-0151-02

1 引言

隨著測序技術(shù)的迅速進(jìn)步，各類生物數(shù)據(jù)庫中的序列數(shù)據(jù)正在快速增長，生物大數(shù)據(jù)正在高速填充世界各生物公共平臺的后臺數(shù)據(jù)庫；僅以美國國家生物技術(shù)信息中心（National Center for Biotechnology Information）的DNA序列數(shù)據(jù)庫GenBank為例，截止2014年8月，GenBank數(shù)據(jù)庫中的堿基對數(shù)量已超過13630億對（Base Pair），較上一年增長了45%[1]，該數(shù)據(jù)庫的堿基對數(shù)量在2013年、2012年和2011年的年增長率分別是43%[2]、45.1%[3]和33.1%[4]。如此快速增長的序列數(shù)據(jù)，僅通過實(shí)驗(yàn)手段對這些序列數(shù)據(jù)進(jìn)行功能注釋顯然已不現(xiàn)實(shí)，基于計(jì)算技術(shù)的快速功能注釋已經(jīng)變得勢在必行。

DNA序列是由A，T，C和G四個字母組成的字符串，而目前的機(jī)器學(xué)習(xí)方法僅以特征向量作為輸入；因此，將DNA序列轉(zhuǎn)化為特征向量并盡可能保留序列內(nèi)部的信息是特征提取技術(shù)的關(guān)鍵。

對新測序或者未知功能的DNA序列，對其功能進(jìn)行驗(yàn)證的可靠方法是人工實(shí)驗(yàn)，但在數(shù)量龐大的DNA序列面前，全部由實(shí)驗(yàn)方法驗(yàn)證其功能顯然已不可行，借助計(jì)算機(jī)領(lǐng)域的機(jī)器學(xué)習(xí)方法快速注釋新序列的潛在功能便是一種可行的途徑。這種功能注釋方法的理論基礎(chǔ)是序列的相似性意味著功能上的相似性。機(jī)器學(xué)習(xí)方法首先要獲得一組DNA序列的訓(xùn)練集，該集合中的序列是已確定其功能的序列，由該訓(xùn)練集構(gòu)建學(xué)習(xí)模型，并在訓(xùn)練集上進(jìn)行交叉檢驗(yàn)來驗(yàn)證該學(xué)習(xí)模型的預(yù)測性能，然后應(yīng)用該模型對未知功能DNA序列進(jìn)行功能預(yù)測。當(dāng)然，并不是所有機(jī)器學(xué)習(xí)方法都適合對DNA序列的功能進(jìn)行預(yù)測，因此，本文對DNA序列的特征向量提取方法及構(gòu)建機(jī)器學(xué)習(xí)模型等問題進(jìn)行了探討。

2 DNA序列的特征提取策略

DNA序列由4種核苷酸堿基組成，分別是腺嘌呤（Adenine， A）、鳥嘌呤（Guanine，G）、胞嘧啶（Cytosine，C）和胸腺嘧啶（Thymine， T）。DNA序列的特征提取就是將由A、G、C和T四個字母組成的長串序列（字符串）轉(zhuǎn)化成用數(shù)值表示的特征向量的過程。

基于k-mer的特征提取方法是一種常用方法?？紤]由字母表∑={A，G，C，T}生成長度為k的序列片段（即k-mer），并統(tǒng)計(jì)這些片段在DNA序列中的出現(xiàn)頻率，由這些頻率值構(gòu)造特征向量。當(dāng)k=1時，即統(tǒng)計(jì)字母表∑中4個字母在序列中的出現(xiàn)頻率，生成一個有4個分量的特征向量。當(dāng)k=1時，一個特征向量僅有4個分量，一般沒有意義。當(dāng)k=2時，即計(jì)算集合∑2={AA， AG， AC， ...， TC， TT}中的16個雙核苷酸堿基在DNA序列中的出現(xiàn)頻率，由此構(gòu)成一個有16個分量的特征向量。例如，一條DNA序列為“ACGT”，則該序列包含三個2-mer分別為AC、CG和GT，這三個2-mer的出現(xiàn)頻率均為1/3=0.33；該序列生成一個有16個分量的向量，其中有三個分量為0.33，即為前面所提到的3個2-mer的出現(xiàn)頻率。當(dāng)k=3時，特征向量的長度為43，即64維。隨著k的增大，特征向量的維度迅速升高，例如，當(dāng)k=8時，表示這條DNA序列的特征向量長度為65536維（48=65536），如此高維的特征向量已引起維度災(zāi)難，機(jī)器學(xué)習(xí)算法在處理高維向量時，其性能會顯著下降，k值并不是越高越好。

基于k-mer的特征提取方案，衍生出一系列的特征提取方法。比如，將不同k值的k-mer組合，生成混合特征向量。例如將k=1、k=2和k=3三類特征向量進(jìn)行組合，生成具有84個分量（41+42+43=84）的特征向量?；趉-mer的編碼思想，王樹林[5]等人提出了基于k-mer的哈希編碼方案。在他們的論文[5]中，將字母表∑中4個字母進(jìn)行二進(jìn)制編碼：Code（A）=（00）2，Code（G）=（01）2，Code（C）=（10）2和Code（T）=（11）2，括號外的下標(biāo)2表示二進(jìn)制，編碼函數(shù)Code（si）表示對字母表∑中的單字符si進(jìn)行二進(jìn)制編碼，并將k-mer短序列通過哈希函數(shù)映射為離散的數(shù)值向量，其哈希函數(shù)f：∑k→N定義為：

s[1..k]表示長度為k的DNA短序列片段，即k-mer。例如，DNA序列為“ACGT”，當(dāng)k=2時，2-mer“AC”的哈希映射可以表示為：f（AC）=42-1·Code（A） + 42-2·Code（C），計(jì)算得f（AC）=2；2-mer“CG” 可表示為f（CG）=42-1·Code（C） + 42-2·Code（G），計(jì)算得f（CG）=9；同理2-mer“GT” 可表示為f（GT）=42-1·Code（G） + 42-2·Code（T），計(jì)算f（GT）=7。因此，DNA序列“ACGT”可表示為向量（2，9，7）。

由于DNA序列由雙鏈構(gòu)成，字母表∑中的4個字符在雙鏈上以互補(bǔ)配對方式出現(xiàn)，即A與T配對，C與G配對。為了消除雙鏈中的單鏈特異性，Noble[6]和劉濱[7]等人應(yīng)用反向互補(bǔ)k-mer對DNA序列進(jìn)行向量化。比如，當(dāng)k=2時，基本的2-mer有16個，即AA， AC， AG， AT， CA， CC， CG， CT， GA， GC， GG， GT， TA，TC， TG， TT，而考慮反向互補(bǔ)后的2-mer則縮減為10個，即為AA，AC， AG， AT， CA， CC， ‘CG， GA， GC， TA，生成的特征向量為10維。

與基于k-mer的特征提取方法不同，基于偽核苷酸組成成分[8]（Pseudo dinucleotide composition）的特征向量提取方法不僅考慮了DNA序列中的局部順序關(guān)系（即k-mer），也考慮了堿基的全局次序模式（即偽核苷酸組成）。例如，基于雙核苷酸的偽雙核苷酸組成的特征向量提取方法可以參考[9]。

3 常用DNA序列的分類算法

當(dāng)使用合適的特征提取方法將DNA序列轉(zhuǎn)化為特征向量后，便可以使用正負(fù)樣本集合生成訓(xùn)練集。生物信息學(xué)領(lǐng)域應(yīng)用比較廣泛的機(jī)器學(xué)習(xí)算法包括支持向量機(jī)（Support Vector Machine， SVM）算法、隨機(jī)森林（Random Forests， RF）算法和人工神經(jīng)網(wǎng)絡(luò)（Artificial Neutral Networks， ANN）等算法。

使用機(jī)器學(xué)習(xí)分類算法的流程一般包括，構(gòu)建已知功能DNA序列的訓(xùn)練集，訓(xùn)練集包括正樣本和負(fù)樣本。正樣本是指具有確定功能的DNA序列，負(fù)樣本指不具有該功能的一般DNA序列。訓(xùn)練集中正負(fù)樣本數(shù)量（DNA序列）應(yīng)該保存一致。為了提高學(xué)習(xí)模型的預(yù)測性能，需要對訓(xùn)練集中的DNA序列通過去重復(fù)來消除偏倚，可使用blast[10]或cd-hit[11]等軟件消除正負(fù)訓(xùn)練集中相似度較高的序列，正負(fù)訓(xùn)練集之間不能有重復(fù)序列，即交集為空，負(fù)訓(xùn)練集應(yīng)該具有一般序列的代表性，然后應(yīng)用特征提取方法將DNA正負(fù)樣本序列轉(zhuǎn)化為帶有標(biāo)號的特征向量集合，并輸入到指定的機(jī)器學(xué)習(xí)算法構(gòu)建預(yù)測模型。為了評價所構(gòu)建模型的分類性能，一般在訓(xùn)練集上使用交叉檢驗(yàn)來對模型進(jìn)行評測。留一法[12] 交叉檢驗(yàn)（leave one out cross validation）和5倍交叉檢驗(yàn)[13]（5-fold cross validation）都是比較常用的方法。5倍交叉檢驗(yàn)是將訓(xùn)練集均分成5份，依次選其中的4份用作模型訓(xùn)練，用剩余一份用作測試集，如此重復(fù)5次；留一法則是選取訓(xùn)練集中的一個樣本用作測試樣本，剩余樣本作為訓(xùn)練樣本，依次讓每個樣本用作測試集僅一次，并用幾項(xiàng)評測指標(biāo)對分類模型進(jìn)行評價。評測指標(biāo)包括：靈敏度（Sensitivity， Sn）、特異度（Specificity， Sp）、準(zhǔn)確度（Precision， Pr）、馬修相關(guān)系數(shù)（Mathews correlation coefficient， MCC）等。

公式中的TP表示真陽性（True Positive），表示訓(xùn)練集正樣本中被預(yù)測為正的樣本數(shù)；TN表示真陰性（True Negative），即負(fù)樣本中預(yù)測為負(fù)的樣本數(shù)；FP表示假陽性（False Positive），表示訓(xùn)練集中負(fù)樣本被預(yù)測為正的樣本數(shù)；FN表示假陰性（False Negative），即正樣本中被預(yù)測為負(fù)的樣本數(shù)。通過上述指標(biāo)，我們可以評價一個機(jī)器學(xué)習(xí)分類模型的性能，并通過修正機(jī)器學(xué)習(xí)算法中的參數(shù)使分類器的性能達(dá)到最優(yōu)。

4 結(jié)束語

本文簡單探討了應(yīng)用機(jī)器學(xué)習(xí)方法對DNA序列進(jìn)行功能預(yù)測的方法，其步驟包括構(gòu)建具有特定功能的DNA序列訓(xùn)練集，訓(xùn)練集中包括正樣本和負(fù)樣本；將正負(fù)樣本通過DNA特征提取方法轉(zhuǎn)化為特征向量集，然后應(yīng)用一種機(jī)器學(xué)習(xí)算法對特征向量集進(jìn)行訓(xùn)練，生成預(yù)測模型，使用交叉檢驗(yàn)方法對預(yù)測模型進(jìn)行參數(shù)調(diào)優(yōu)，應(yīng)用該模型即可對未知功能的DNA序列進(jìn)行功能預(yù)測，判斷未知功能DNA序列是否具有相應(yīng)的功能。在本文中，我們簡單介紹了使用機(jī)器學(xué)習(xí)算法對DNA序列進(jìn)行功能預(yù)測的一般過程，希望能對機(jī)器學(xué)習(xí)方法在生物信息學(xué)領(lǐng)域的應(yīng)用起到拋磚引玉的作用。

參考文獻(xiàn)：

[1] K. Clark， I. Karsch-Mizrachi， D. J. Lipman， et al.， "GenBank，" Nucleic Acids Res， 2016（44）：67-72.

[2] D. A. Benson， K. Clark， I. et al.， "GenBank，" Nucleic Acids Res， 2015（43）：30.

[3] D. A. Benson， K. Clark， I. Karsch-Mizrachi， et al.， "GenBank，" Nucleic Acids Res，2014（42）：32.

[4] D. A. Benson， M. Cavanaugh， K. Clark， I. Karsch-Mizrachi， D. J. Lipman， J. Ostell， et al.， "GenBank，" Nucleic Acids Res，2013（41）：36-42.

[5] 王樹林， 王戟， 陳火旺， 等.k-長DNA子序列計(jì)數(shù)算法研究[J].計(jì)算機(jī)工程，2007（33）：3.

[6] W. S. Noble， S. Kuehn， R. Thurman， et al.， "Predicting the in vivo signature of human gene regulatory sequences，" Bioinformatics， 2005，21（1）：338.

[7] B. Liu， R. Long， K. C. Chou， "iDHS-EL： identifying DNase I hypersensitive sites by fusing three different modes of pseudo nucleotide composition into an ensemble learning framework，" Bioinformatics， Apr 8 2016.

[8] K. C. Chou， "Prediction of protein cellular attributes using pseudo-amino acid composition，" Proteins，2001（43）：246-55.

[9] W. Chen， T. Y. Lei， D. C. Jin， et al.， "PseKNC： a flexible web server for generating pseudo K-tuple nucleotide composition，" Anal Biochem，2014（456）：53-60.

[10] S. F. Altschul， W. Gish， W. Miller， E. et al.， "Basic local alignment search tool，" J Mol Biol， 1990（215）：403-10.

[11] W. Li，A. Godzik， "Cd-hit： a fast program for clustering and comparing large sets of protein or nucleotide sequences，" Bioinformatics， 2006（22）：1658.

[12] H. Bhaskar， D. C. Hoyle， S. Singh， "Machine learning in bioinformatics： a brief survey and recommendations for practitioners，" Comput Biol Med， 2006（36）：1104-25.

[13] G. Liu， "Using weighted features to predict recombination hotspots in Saccharomyces cerevisiae，" Journal of Theoretical Biology， 2016.