一種基于信息容量的模體比較非比對(duì)度量算法

郭麗娟，張少?gòu)?qiáng)，花季偉

（天津師范大學(xué)計(jì)算機(jī)與信息工程學(xué)院，天津300387）

破譯基因組中復(fù)雜的基因調(diào)控網(wǎng)絡(luò)是一項(xiàng)極具挑戰(zhàn)性的課題[1-2]，要實(shí)現(xiàn)這個(gè)目標(biāo)，首先要在基因組中鑒別所有轉(zhuǎn)錄因子的結(jié)合位點(diǎn)[1，3-4].轉(zhuǎn)錄因子結(jié)合位點(diǎn)是轉(zhuǎn)錄因子結(jié)合的一組短基因序列，長(zhǎng)度通常為6～25個(gè)堿基對(duì)（base pair，bp），具有觸發(fā)細(xì)胞轉(zhuǎn)錄調(diào)控的功能.屬于同一轉(zhuǎn)錄因子的結(jié)合位點(diǎn)通常具有特定的保守性和相同的長(zhǎng)度，但它們也可以顯示出一定程度的變異，而這些結(jié)合位點(diǎn)位于一段非常長(zhǎng)的非編碼序列中，這些都會(huì)導(dǎo)致它們的預(yù)測(cè)計(jì)算變得非常困難.同一轉(zhuǎn)錄因子的一組具有高保守性的結(jié)合位點(diǎn)通常被稱作模體，可以由實(shí)驗(yàn)驗(yàn)證，或者通過比較一組可能含有轉(zhuǎn)錄因子結(jié)合位點(diǎn)的方法預(yù)測(cè).由于轉(zhuǎn)錄因子結(jié)合位點(diǎn)比它們周圍的DNA片段更保守，所以，許多從頭測(cè)序模體查找算法被開發(fā)出來用于識(shí)別轉(zhuǎn)錄因子結(jié)合位點(diǎn).模體可以由位置賦權(quán)矩陣（position weight matrix，PWM） 和位置頻率矩陣（position frequency matrix，PFM）精確地表述出來[5-6].這2個(gè)矩陣是模體結(jié)合位點(diǎn)序列比對(duì)的變形，它們極大程度上反映了相應(yīng)的轉(zhuǎn)錄因子的位置結(jié)合傾向性.因此，通過在這2個(gè)矩陣中掃描可能包含TFBS的序列即可以發(fā)現(xiàn)模體.

利用模體查找工具獲得一些假定的模體后，通過在一個(gè)模體數(shù)據(jù)庫(kù)中找到與假定模體匹配的模體，從而推斷出轉(zhuǎn)錄因子附屬于這些假定的模體[7]；或聚類相同轉(zhuǎn)錄因子的相似子模體從而去除冗余模體、形成一個(gè)完整的模體.因此，在上述提到的應(yīng)用中，需要一種有效的度量法用于捕獲不相關(guān)模體之間的細(xì)微差別，強(qiáng)調(diào)同種群間模體的相似度.目前，計(jì)算模體相似度的比對(duì)度量方法包括兩大類：一類是列相似度度量法，即從2個(gè)模體的位置頻率矩陣（或者位置賦權(quán)矩陣）中各取一列計(jì)算相似性，如SSD（sum of squared distances）[8-9]、pCS （p-value of Chi-square）[10]、ALLR （average log-likelihood ratio）[11]、AKL（average Kullback-Leibler，AKL）[12]和 PCC（pearson’scorrelation coefficient，PCC）[13]等；另一類是雙序列比對(duì)算法，利用列相似度度量法和一個(gè)空位罰分函數(shù)作為分?jǐn)?shù)比對(duì)2個(gè)模體[14]，在假定具有空位罰分函數(shù)的情況下，Needleman-Wunsch 算法[15]和 Smith-Waterman（SW）算法[16]都可用于查找最優(yōu)比對(duì).文獻(xiàn)[7]和文獻(xiàn)[17]對(duì)這些度量和比對(duì)算法進(jìn)行評(píng)估后，建立了網(wǎng)絡(luò)服務(wù)器STAMP，用于集成這些帶有比對(duì)的度量法.除此之外，還有2個(gè)用于比較模體的非比對(duì)度量方法KFV（kmer frequency vector） 和 Mosta 包 的 AC（Asymptotic Covariance），它們分別由文獻(xiàn)[7]和文獻(xiàn)[15]提出.

上述度量方法僅用到了位置頻率矩陣，均沒有使用列信息容量（column information contents）和位置賦權(quán)矩陣.實(shí)際上，由于矩陣中所對(duì)應(yīng)的列具有很高的相關(guān)性，2個(gè)總體信息容量低的模體也可能具有高的相似度分?jǐn)?shù).因此，如果2個(gè)模體某些列的信息容量很低，在應(yīng)用這些度量方法比對(duì)前，低信息容量的列就要被刪除.上述帶有比對(duì)的度量方法在聚類相似模體方面具有較好的效果，但是它們基本不能從帶有低信息容量列的混雜模體中分離出真模體.此外，帶有比對(duì)的度量公式需要用到比對(duì)算法，而比對(duì)算法依賴的參數(shù)較多，運(yùn)行所需時(shí)間也較多，因此基于非比對(duì)的度量公式可以更快速精確地進(jìn)行模體比較.綜上所述，本研究提出一種帶有位置信息容量的相似度非比對(duì)度量法（information contents based similarity metric，ICBSM）.算法中不僅包含位置頻率矩陣和位置賦權(quán)矩陣，還加入了每個(gè)位置的信息內(nèi)容，并利用來自于 STAMP[17]、KFV[14]和 GLECLUBS[18-19]中的數(shù)據(jù)集對(duì)ICBSM算法進(jìn)行評(píng)估，將該算法與國(guó)際上已經(jīng)提出的算法進(jìn)行比較分析.

1 ICBSM方法

1.1 ICBSM度量法的提出

設(shè)模體Motif由n個(gè)長(zhǎng)度為L(zhǎng)的序列組成，定義其位置頻率矩陣為

式（3）用于表示PWM1生成PFM2的可能性，其中 Alignment（1，2）是通過固定 PWM1、滑動(dòng) PFM2得到的矩陣列的比對(duì).圖1為PFM在PWM上的滑動(dòng)示意圖.

圖1 PFM在PWM上的滑動(dòng)示意圖Fig.1 Representation of PFM sliding on PWM

由圖1可以看出，當(dāng)用PFM在PWM上逐列由左向右滑動(dòng)形成比對(duì)s時(shí)，在該比對(duì)s中，PWM的第i列與PFM 的第 s（i）列對(duì)應(yīng).

再用Motif2的PFM2和Motif1的PWM1進(jìn)行比對(duì)，計(jì)算相似性：

1.2 算法驗(yàn)證

為了驗(yàn)證算法的性能，利用經(jīng)過驗(yàn)證的3個(gè)數(shù)據(jù)集對(duì)ICBSM進(jìn)行測(cè)試和評(píng)估.數(shù)據(jù)集-1由Mahony等[3]從JASPAR庫(kù)中首次選出，該數(shù)據(jù)庫(kù)由96個(gè)真實(shí)的模體組成，這些模體屬于13個(gè)已知的不同結(jié)構(gòu)的TF類.文獻(xiàn)[6]創(chuàng)建了數(shù)據(jù)集-2，用以測(cè)試KFV度量法對(duì)于識(shí)別冗余的位置頻率矩陣的顯著性能.該數(shù)據(jù)集由124個(gè)JASPAR的核心模體及每個(gè)核心模體的3個(gè)子模體組成，這些子模體通過隨機(jī)選取每個(gè)模體的2/3序列得到.數(shù)據(jù)集-3可由http：//gleclubs.uncc.edu/pbs頁(yè)面下載，包含了大約105個(gè)假定的模體[18-19]，這些模體來自大腸桿菌2000多組全基因組的同源基因間序列以及其他54個(gè)γ-變形菌門的參照基因組.關(guān)于3個(gè)數(shù)據(jù)集的詳細(xì)參數(shù)參見表1.

表1 用于測(cè)試與評(píng)估的3個(gè)數(shù)據(jù)集的參數(shù)Tab.1 Parameters of three datasets for testing and assessing

將ICBSM算法、STAMP工具包中的5個(gè)算法（ALLR，AKL，SSD，pCS，PCC）、KFV 法以及 Mosta算法中的AC應(yīng)用到1個(gè)數(shù)據(jù)集上，針對(duì)聚類相關(guān)的真模體、過濾偽模體和找回模體等方面進(jìn)行性能比較.利用STAMP平臺(tái)計(jì)算5個(gè)依靠比對(duì)的度量法得分（http：//www.benoslab.pitt.edu/stamp/），Mosta包計(jì)算AC得分（http：//mosta.molgen.mpg.de），KFV的網(wǎng)絡(luò)服務(wù)器計(jì)算KFV得分（http：//bioinfo.uncc.edu/kfv/）.

1.2.1 模體找回

帶有比對(duì)算法的列相似度度量法和非比對(duì)相似度度量法可以用于將待查模體與數(shù)據(jù)庫(kù)中的每一個(gè)模體進(jìn)行比較，從而找回模體.如果在1個(gè)數(shù)據(jù)集中模體的相似度分?jǐn)?shù)超過閾值，則表明這些模體被待查模體“命中”；如果有多個(gè)“命中”[6]，則相似度分?jǐn)?shù)最高的“命中”稱為“最佳命中”.通過使用“最佳命中法”把在數(shù)據(jù)庫(kù)中搜索模體的正確找回率定義為度量法的“性能精確度”.

與其他3個(gè)帶有最優(yōu)比對(duì)[6]的列相似度度量法和非比對(duì)度量的AC法[14]相比，SSD、PCC和KFV度量法在查找模體時(shí)具有更高精確度，因此選出它們與ICBSM度量法進(jìn)行對(duì)比，比較它們?cè)跀?shù)據(jù)集-1中找回同一個(gè)轉(zhuǎn)錄因子家族模體的能力.在數(shù)據(jù)集-1上，STAMP包的5個(gè)帶有比對(duì)設(shè)置的列相似度度量法中，結(jié)合SW的非空位比對(duì)算法PCC（PCC/SWU）和結(jié)合SW、空位延伸為0.5、空位開放為1的比對(duì)算法SSD/SW是最好的2個(gè)度量法和比對(duì)設(shè)置[7].根據(jù)文獻(xiàn)[6]的描述，當(dāng)把4-mer和夾角余弦值用于向量構(gòu)建和比較時(shí)，KFV會(huì)獲得最優(yōu)結(jié)果.

本研究利用 ROC（receiver operating characteristic）曲線考察度量法在數(shù)據(jù)集-1和數(shù)據(jù)集-2中識(shí)別出相同轉(zhuǎn)錄因子的模體的性能.ROC曲線的繪制方法依據(jù)下述規(guī)則：給定1個(gè)由n個(gè)模體組成的數(shù)據(jù)集，其中這些模體的轉(zhuǎn)錄因子結(jié)構(gòu)類已知，n個(gè)模體具有n（n+1）/2個(gè)組對(duì)，應(yīng)用度量法分別計(jì)算出每一對(duì)的相似度分?jǐn)?shù).如果2個(gè)模體的相似度分?jǐn)?shù)小于1個(gè)閾值或大于閾值但沒有“最佳命中”，則設(shè)定這2個(gè)模體為錯(cuò)誤匹配，否則為正確匹配.如果由度量法計(jì)算出2個(gè)模體正確匹配，且這2個(gè)模體確實(shí)同屬于1個(gè)轉(zhuǎn)錄因子，則該正確匹配稱為“真陽(yáng)性（true positive，TP）”，否則這個(gè)正確匹配為“假陽(yáng)性（false positive，F(xiàn)P）”；如果2個(gè)模體由度量法計(jì)算出是錯(cuò)誤匹配，且這2個(gè)模體確實(shí)屬于不同的轉(zhuǎn)錄因子，則該匹配稱為“真陰性（true negative，TN）”，否則這個(gè)錯(cuò)誤匹配為“假陰性（false negative，F(xiàn)N）”.ROC曲線是在不同的模體相似度閾值下由真陽(yáng)性率對(duì)比假陽(yáng)性率的描述.

1.2.2 從混雜的模體中分離出真模體

一些基于遺傳系譜印技術(shù)的轉(zhuǎn)錄因子綁定位點(diǎn)的全基因組測(cè)序算法需要把任意轉(zhuǎn)錄因子的子模體和冗余模體合并成一個(gè)獨(dú)立的模體并剔除偽模體[8-9，13]，即聚類相似模體，區(qū)分出不相關(guān)的模體.因此，研究人員需要一個(gè)不僅能精確計(jì)算出一對(duì)模體的精確度，而且還能有效區(qū)分出無關(guān)模體的度量法，這個(gè)算法可以為相同轉(zhuǎn)錄因子模體的2個(gè)子模體賦予足夠高的相似度值，為沒有任何進(jìn)化關(guān)系的2個(gè)模體賦予足夠低的相似度值，從而在混雜的模體中分離出真模體.由GLECLUBS生成的數(shù)據(jù)集-3[8]由大量的混雜模體和一小部分的真模體構(gòu)成，為從數(shù)據(jù)集-3中發(fā)現(xiàn)真模體，在Regulon數(shù)據(jù)庫(kù)中選出一組真模體用于在數(shù)據(jù)集-3上進(jìn)行評(píng)估.該組真實(shí)模體是大腸桿菌的122個(gè)轉(zhuǎn)錄因子模體生成的大量的真的子模體.每個(gè)轉(zhuǎn)錄因子模體均是由n個(gè)結(jié)合位點(diǎn)構(gòu)成（n≥3），度量法把每個(gè)轉(zhuǎn)錄因子模體隨機(jī)分成1個(gè)大小為k的子模體和1個(gè)大小為（n-k）的子模體，其中 k∈{1，2，…，[n/2]}.因此，每個(gè)大小為n的模體都可以生成［n/2］對(duì)的子模體.度量法對(duì)每個(gè)大小為k的子模體重復(fù)前面的分離過程，生成［k/2］對(duì)子模體的子模體.當(dāng)每個(gè)子模體的大小為1時(shí)，過程停止.然后，利用這些度量法計(jì)算每對(duì)子模體間的相似度值[7，11]，并在數(shù)據(jù)集-3上計(jì)算每對(duì)模體的相似度值.通過計(jì)算數(shù)據(jù)庫(kù)-3中每對(duì)模體相似度分?jǐn)?shù)標(biāo)準(zhǔn)化后的分布和每對(duì)真的子模體的相似度分?jǐn)?shù)標(biāo)準(zhǔn)化后的分布，查看2個(gè)分布的重疊區(qū)域.

2 ICBSM算法性能分析與結(jié)果

2.1 模體找回

對(duì)于從一個(gè)數(shù)據(jù)集中找回模體，本研究將模體比較的閾值設(shè)置為0.6，然后將ICBSM、KFV、PCC/SWU和SSD/SW算法在數(shù)據(jù)集-1上計(jì)算精確度，結(jié)果如表2所示.

表2 在數(shù)據(jù)集-1上，ICBSM、KFV、PCC/SWU及SSD/SW模體找回的精確度Tab.2 Accuracy for searching motifs of ICBSM，KFV，PCC/SWU and SSD/SW on dataset-1

數(shù)據(jù)集-1可以分為包含25個(gè)真實(shí)模體的鋅指狀結(jié)構(gòu)蛋白質(zhì)家族（zinc-finger，ZF）和包含71個(gè)真實(shí)模體的非鋅指狀結(jié)構(gòu)蛋白質(zhì)家族（non-ZF）.由表2中結(jié)果可知，對(duì)于ZF蛋白質(zhì)家族、Non-ZF蛋白質(zhì)家族以及整個(gè)蛋白質(zhì)家族集合，ICBSM算法的模體找回精度最高，說明該算法在數(shù)據(jù)庫(kù)中能夠正確找回模體的能力最強(qiáng)，比其他3種度量法具有更卓越的策略.

為了將ICBSM與PCC/SWU、KFV（4-mer夾角余弦值）的最優(yōu)策略做進(jìn)一步比較，在模體比較閾值設(shè)置為0.6的情況下，在數(shù)據(jù)集-1和數(shù)據(jù)集-2上，對(duì)這3種策略的性能進(jìn)行ROC分析，結(jié)果如圖2所示.由圖2可知，假陽(yáng)性率相同的情況下，ICBSM度量法的真陽(yáng)性率最高，即對(duì)于同1個(gè)數(shù)據(jù)集，ICBSM度量法能夠正確找回模體的能力比其他2種方法更強(qiáng).

圖23 種度量法的ROC曲線圖Fig.2 ROC curves of three metrics

2.2 從混雜的模體中分離出真模體

用ICBSM度量法、STAMP工具包、AC度量法以及KFV度量法分別計(jì)算數(shù)據(jù)集-3的每對(duì)模體相似度分?jǐn)?shù)以及每對(duì)真的子模體的相似度分?jǐn)?shù)，并將這2個(gè)分?jǐn)?shù)標(biāo)準(zhǔn)化形成曲線分布圖，以ICBSM度量法與AKL度量法曲線分布效果為例，結(jié)果如圖3所示.

圖3中“數(shù)據(jù)集-3模體”的曲線是在數(shù)據(jù)集-3中計(jì)算每對(duì)模體相似度分?jǐn)?shù)標(biāo)準(zhǔn)化后的分布曲線，標(biāo)有“真的子模體”曲線是每對(duì)真的子模體的標(biāo)準(zhǔn)化相似度分?jǐn)?shù)的分布曲線.在數(shù)據(jù)集-3中，由于每對(duì)真的子模體具有相關(guān)性而大多數(shù)模體具有無關(guān)性，因此性能好的度量法應(yīng)該可以把前一個(gè)相似度分布區(qū)域與后一個(gè)相似度分布區(qū)域分離出來，即圖3中2個(gè)曲線所圍成的2個(gè)區(qū)域的重疊部分越小，分離效果越好.ICBSM在計(jì)算模體的相似度分?jǐn)?shù)時(shí)考慮了信息容量，因此可以從帶有低信息容量的混雜的模體中分離出真模體.

將ICBSM與其他度量法生成的相似度分布曲線的重疊區(qū)域比率進(jìn)行比較，結(jié)果如圖4所示.在ICBSM的分布曲線下，2塊區(qū)域具有最小的重疊部分，這說明與其他度量法相比，ICBSM能夠更加精確地從混亂模體中分離出真模體.

3 結(jié)論

在生物信息處理過程中，由于很多應(yīng)用都包含了模體比較的過程，因此提出一種基于列信息內(nèi)容的用于模體比較的非比對(duì)度量法ICBSM，通過對(duì)比分析，結(jié)果表明：

（1）ICBSM度量法采用了帶有信息容量的非比對(duì)策略計(jì)算模體間的相似度分?jǐn)?shù)，將信息容量添加到模體的位置賦權(quán)矩陣上，將一個(gè)模體的位置頻率矩陣在另一個(gè)模體的位置賦權(quán)矩陣上滑動(dòng)，計(jì)算2個(gè)模體間的相似度.該算法依賴參數(shù)少，提升了計(jì)算效率.

（2）在模體比較的閾值設(shè)置為0.6的情況下，在數(shù)據(jù)集-1上，ICBSM度量法與KFV、PCC/SWU及SSD/SW相比較，其模體找回的精確度最高；同時(shí)，與KFV、PCC/SWU相比較，ICBSM的ROC曲線的真陽(yáng)性率值也最高，這說明該方法在數(shù)據(jù)庫(kù)中找回模體的效果更好.

（3）由于ICBSM在計(jì)算模體相似度時(shí)考慮了模體的信息容量，因此它計(jì)算出的真的子模體的相似度分?jǐn)?shù)標(biāo)準(zhǔn)化后的分布曲線與數(shù)據(jù)集中所有模體的相似度分?jǐn)?shù)標(biāo)準(zhǔn)化后的分布曲線重疊率最低，說明該方法能夠精確地將真模體從混雜的模體中區(qū)分出來，為聚類相似模體、分組不相關(guān)模體提供了有效工具.

[1]CELNIKER S E，DILLON L A，GERSTEIN M B，et al.Unlocking the secrets of the genome[J].Nature，2009，459（7249）：927-930.

[2] RISTER J，DESPLAN C.Deciphering the genome's regulatory code：the many languages of DNA[J].Bioessays，2010，32（5）：381-384.

[3] REED J L，F(xiàn)AMILI I，THIELE I，et al.Towards multidimensional genome annotation[J].Nat Rev Genet，2006，7（2）：130-141.

[4]ALEXANDER RP，F(xiàn)ANG G，ROZOWSKY J，SNYDER M，et al.Annotating non-coding regions of the genome[J].Nat Rev Genet，2010，11（8）：559-571.

[5] GUHATHAKURTA D.Computational identification of transcriptional regulatory elements in DNA sequence[J].Nucleic Acids Res，2006，34（12）：3585-3598.

[6] STORMO G D.DNA binding sites：representation and discovery[J].Bioinformatics，2000，16（1）：16-23.

[7]MAHONY S，AURON PE，BENOS P V.DNA familial binding profiles made easy：comparison of various motif alignment and clustering strategies[J].PLoS Comput Biol，2007，3（3）：61.

[8]SANDELIN A，WASSERMAN W W.Constrained binding site diversity within families of transcription factors enhances pattern discovery bioinformatics[J].J Mol Biol，2004，338（2）：207-215.

[9] WANG T，STORMO G D.Identifying the conserved network of cisregulatory sites of a eukaryotic genome[J].Proc Natl Acad Sci USA，2005，102（48）：17400-17405.

[10]SCHONES D E，SUMAZIN P，ZHANG M Q.Similarity of position frequency matrices for transcription factor binding sites[J].Bioinformatics，2005，21（3）：307-313.

[11]WANG T，STORMO G D.Combining phylogenetic data with coregulated genes to identify regulatory motifs[J].Bioinformatics，2003，19（18）：2369-2380.

[12]KULLBACK S，LEIBLER R A.On information and sufficiency[J].Ann Math Statist，1951，22（1）：79-86.

[13]PIETROKOVSKI S.Searching databases of conserved sequence regions by aligning protein multiple-alignments[J].Nucleic Acids Res，1996，24（19）：3836-3845.

[14]XU M，SU Z.A novel alignment-free method for comparing transcription factor binding site motifs[J].PLoS One，2010，5（1）：87-97.

[15]NEEDLEMAN S B，WUNSCH C D.A general method applicable to the search for similarities in the amino acid sequence of two proteins[J].J Mol Biol，1970，48（3）：443-453.

[16]SMITH T F，WATERMAN M S.Identification of common molecular subsequences[J].J Mol Boil，1981，147（1）：195-197.

[17]MAHONY S，BENOS P V.STAMP：a web tool for exploring DNA-binding motif similarities[J].Nucleic Acids Res，2007，35：253-258.

[18]ZHANG S，XU M，LI S，et al.Genome-wide de novo prediction of cisregulatory binding sites in prokaryotes[J].Nucleic Acids Res，2009，37（10）：72.

[19]ZHANG S，LI S，PHAM P T，et al.Simultaneous prediction of transcri-ption factor binding sites in a group of prokaryotic genomes[J].BMC Bi-oinformatics，2010，11：397.