高通量DNA測(cè)序數(shù)據(jù)的生物信息學(xué)方法

詹曉娟，姚登舉，朱懷球

1. 黑龍江工程學(xué)院計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，黑龍江 哈爾濱 150050；

2. 哈爾濱理工大學(xué)軟件學(xué)院，黑龍江 哈爾濱 150040；3. 北京大學(xué)生物醫(yī)學(xué)工程系，北京 100871

◎研究

高通量DNA測(cè)序數(shù)據(jù)的生物信息學(xué)方法

詹曉娟1，姚登舉2，朱懷球3

1. 黑龍江工程學(xué)院計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，黑龍江 哈爾濱 150050；

2. 哈爾濱理工大學(xué)軟件學(xué)院，黑龍江 哈爾濱 150040；3. 北京大學(xué)生物醫(yī)學(xué)工程系，北京 100871

高通量測(cè)序技術(shù)產(chǎn)生的DNA序列數(shù)據(jù)長(zhǎng)度較短，而且數(shù)據(jù)量非常巨大。分析了高通量測(cè)序環(huán)境下大數(shù)據(jù)的挑戰(zhàn)和機(jī)遇，總結(jié)并討論了數(shù)據(jù)壓縮、宏基因組數(shù)據(jù)序列拼接、宏基因組數(shù)據(jù)序列分析方面的算法和工具等研究成果。最后，展望了高通量測(cè)序下DNA短讀序列數(shù)據(jù)研究的發(fā)展趨勢(shì)。

高通量DNA測(cè)序；生物信息學(xué)；短讀序列數(shù)據(jù)壓縮；短讀序列數(shù)據(jù)拼接；短讀序列數(shù)據(jù)分析

1 引言

高通量測(cè)序技術(shù)又稱“下一代”測(cè)序（next-generation sequencing, NGS）技術(shù)[1]，可以一次性測(cè)定幾十萬(wàn)甚至幾百萬(wàn)條序列，是現(xiàn)今應(yīng)用最廣泛的測(cè)序技術(shù)。相對(duì)于傳統(tǒng)的Sanger測(cè)序技術(shù)[2]，NGS具有高速、高通量、低價(jià)格等優(yōu)點(diǎn)[3]。高通量測(cè)序數(shù)據(jù)廣泛應(yīng)用于生物學(xué)、醫(yī)學(xué)、遺傳科學(xué)等諸多領(lǐng)域，具有重要研究?jī)r(jià)值。許多大型的科學(xué)研究項(xiàng)目，如千人基因組計(jì)劃（1 000 genome project)、DNA元件百科全書(shū)（encyclopedia of DNA elements）計(jì)劃、國(guó)際癌癥基因組計(jì)劃（international cancer genome project）等，正以前所未有的速度產(chǎn)生海量DNA序列。截至2014年2月，僅登錄在美國(guó)GenBank數(shù)據(jù)庫(kù)中的DNA序列數(shù)據(jù)就有十萬(wàn)億堿基對(duì)，所有高通量測(cè)序下的DNA短讀序列數(shù)據(jù)大小達(dá)到上千PB。隨著測(cè)序技術(shù)的不斷改善和測(cè)序成本的持續(xù)降低，每天都會(huì)有海量的DNA序列產(chǎn)生，使得生物數(shù)據(jù)量呈指數(shù)規(guī)模增長(zhǎng)，平均約每14個(gè)月增加一倍。圖1對(duì)高通量測(cè)序平臺(tái)下的短讀（short reads，以下簡(jiǎn)稱reads）序列數(shù)據(jù)和其他大數(shù)據(jù)領(lǐng)域的原始數(shù)據(jù)增長(zhǎng)方式進(jìn)行了比較，陰影區(qū)預(yù)報(bào)了未來(lái)的增長(zhǎng)趨勢(shì)，從圖11可以看出，高通量測(cè)序下的基因組序列數(shù)據(jù)即短讀序列數(shù)據(jù)的增長(zhǎng)遠(yuǎn)大于摩爾定律的增長(zhǎng)速度。計(jì)算機(jī)是存儲(chǔ)和處理DNA數(shù)據(jù)的主要工具，其微處理器性能和存儲(chǔ)設(shè)備容量平均18～24個(gè)月翻一番，而DNA測(cè)序數(shù)據(jù)平均4～5個(gè)月就翻一番，DNA測(cè)序數(shù)據(jù)的增長(zhǎng)速度已經(jīng)遠(yuǎn)遠(yuǎn)超過(guò)了計(jì)算機(jī)微處理器和存儲(chǔ)設(shè)備的增長(zhǎng)速度。面對(duì)如此迅速增長(zhǎng)的龐大的短讀序列數(shù)據(jù)集，如何有效管理、分析、充分利用這些信息，已成為生物信息學(xué)發(fā)展亟需解決的問(wèn)題[4]。

圖1 不同種類數(shù)據(jù)的近似增長(zhǎng)趨勢(shì)

2 生物大數(shù)據(jù)帶來(lái)的新挑戰(zhàn)

隨著高通量測(cè)序技術(shù)的發(fā)展，各種生物學(xué)數(shù)據(jù)呈現(xiàn)爆炸式增長(zhǎng)，并且這一趨勢(shì)將隨著生物測(cè)序技術(shù)的發(fā)展而進(jìn)一步增強(qiáng)。面對(duì)生命科學(xué)領(lǐng)域的大數(shù)據(jù)分析任務(wù)，多種不同維度的數(shù)據(jù)整合、多學(xué)科交叉的數(shù)據(jù)分析以及經(jīng)典的數(shù)據(jù)挖掘算法都面臨新的挑戰(zhàn)。

2.1 多學(xué)科交叉的挑戰(zhàn)

自從1990年人類基因組計(jì)劃正式啟動(dòng)以來(lái)，20余年間，各種基因組、蛋白質(zhì)組、轉(zhuǎn)錄組、宏基因組等國(guó)際生物學(xué)研究合作計(jì)劃開(kāi)始啟動(dòng)或已完成，目前國(guó)際上已經(jīng)成立了多個(gè)大的跨國(guó)科研合作機(jī)構(gòu)，生物信息領(lǐng)域的國(guó)際合作與交流也不斷加強(qiáng)（見(jiàn)表11）。各種組學(xué)和生物信息學(xué)領(lǐng)域的國(guó)際化和跨學(xué)科間的專家合作使得團(tuán)隊(duì)成員在該領(lǐng)域取得了突出的成果，不僅發(fā)表了很多有影響力的文章，而且開(kāi)發(fā)出許多新的數(shù)據(jù)集成和分析工具，以便資源和信息共享[5]。然而，面對(duì)飛速增長(zhǎng)的生物學(xué)大數(shù)據(jù)和日漸增多的生物信息學(xué)研究任務(wù)，跨學(xué)科的國(guó)際合作仍面臨巨大的挑戰(zhàn)，例如不同的實(shí)驗(yàn)室和平臺(tái)產(chǎn)生的大數(shù)據(jù)如何實(shí)現(xiàn)無(wú)障礙的共享和協(xié)作分析，不同組學(xué)產(chǎn)生的數(shù)據(jù)如何有效地進(jìn)行集成、管理、維護(hù)和更新，如何開(kāi)發(fā)新型的面向生物學(xué)大數(shù)據(jù)分析的算法和工具等。

2.2 數(shù)據(jù)和工具的整合問(wèn)題

目前主流的高通量測(cè)序平臺(tái)主要有Roche/454焦磷酸測(cè)序、Solexa/Illumina邊合成邊測(cè)序和ABI SOLiD連接測(cè)序。高通量測(cè)序技術(shù)的讀長(zhǎng)較短，但測(cè)序深度可以在一定程度上彌補(bǔ)讀長(zhǎng)較短帶來(lái)的問(wèn)題。其中，454測(cè)序平臺(tái)讀長(zhǎng)最長(zhǎng)有450～800 bp，適合對(duì)未知基因組從頭測(cè)序；Solexa/Illumina測(cè)序讀長(zhǎng)比454測(cè)序平臺(tái)短，但測(cè)序通量高、價(jià)位低，適合基因組重測(cè)序；SOLiD讀長(zhǎng)也較短，但測(cè)序精度高，特別適合SNP檢測(cè)等。目前應(yīng)用較普遍的是Illumina測(cè)序平臺(tái)，約占現(xiàn)有測(cè)序工具數(shù)量的一半。

表1 生物大數(shù)據(jù)項(xiàng)目合作計(jì)劃

不同的測(cè)序平臺(tái)產(chǎn)生的數(shù)據(jù)格式各不相同，常用的文件格式有.bam、.csfasta、.fasta、.fastq、.gvf、.sam、.tar、.tiff、.var、.vcf等?，F(xiàn)有的數(shù)據(jù)分析工具大多只能分析特定格式的數(shù)據(jù)，在實(shí)際的數(shù)據(jù)分析過(guò)程中往往需要把不同格式的數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化并重新整合，因此會(huì)浪費(fèi)很多時(shí)間進(jìn)行數(shù)據(jù)的預(yù)處理。例如，不同測(cè)序平臺(tái)會(huì)產(chǎn)生不同品質(zhì)和長(zhǎng)度的高通量短讀數(shù)據(jù)，由于沒(méi)有統(tǒng)一的行業(yè)標(biāo)準(zhǔn)來(lái)描述高通量測(cè)序下的核苷酸序列和質(zhì)量分?jǐn)?shù)值，導(dǎo)致需要跨平臺(tái)進(jìn)行序列分析。因此，開(kāi)發(fā)一組可以運(yùn)行在不同計(jì)算平臺(tái)下的互操作數(shù)據(jù)分析工具是一個(gè)具有挑戰(zhàn)性的課題。

表22列出了目前高通量測(cè)序下各種組學(xué)所使用的工具和方法。隨著這些多樣的組學(xué)數(shù)據(jù)的整合，數(shù)據(jù)分析和解釋的規(guī)模大大增加，這樣就對(duì)基因組學(xué)和生命科學(xué)領(lǐng)域的大數(shù)據(jù)工具和基礎(chǔ)設(shè)施提出更高的要求。對(duì)不同來(lái)源、不同形式的數(shù)據(jù)進(jìn)行挖掘、評(píng)估、整合和應(yīng)用還亟待加強(qiáng)。未來(lái)，多種組學(xué)數(shù)據(jù)的整合分析將會(huì)挑戰(zhàn)傳統(tǒng)的思維模式，發(fā)揮其至關(guān)重要的作用。

2.3 構(gòu)建新型學(xué)術(shù)交流平臺(tái)日益迫切

隨著高通量測(cè)序成本的降低，生物大數(shù)據(jù)對(duì)于傳統(tǒng)的數(shù)據(jù)存儲(chǔ)、分析和解釋提出了新的挑戰(zhàn)，而將這些數(shù)據(jù)和成果進(jìn)行系統(tǒng)整合并應(yīng)用于醫(yī)療實(shí)踐才剛剛開(kāi)始。當(dāng)前，一些小的實(shí)驗(yàn)室顯然不具備存儲(chǔ)和處理大數(shù)據(jù)的基礎(chǔ)設(shè)施和能力。隨著互聯(lián)網(wǎng)技術(shù)的快速發(fā)展，眾多的科學(xué)合作網(wǎng)絡(luò)平臺(tái)提供了實(shí)時(shí)的數(shù)據(jù)交換，使得人們可以通過(guò)互聯(lián)網(wǎng)方便地進(jìn)行數(shù)據(jù)分享和成果交流。例如，Illumina公司的新一代測(cè)序云計(jì)算平臺(tái)BaseSpace（www.basepace. c o m）、開(kāi)放科學(xué)框架平臺(tái)（h t t p:// openscienceframework.org）和Figshare（http://epic.org/privacy/medical）等。全球三大IT公司Amazon、Rackspace和Google都提供了云存儲(chǔ)和計(jì)算解決方案，通過(guò)云計(jì)算平臺(tái)可以實(shí)現(xiàn)大型數(shù)據(jù)中心的資源共享。然而，云計(jì)算基因組學(xué)也面臨著數(shù)據(jù)隱私和病人數(shù)據(jù)的合法性問(wèn)題，拓展新型的學(xué)術(shù)交流平臺(tái)成為生物大數(shù)據(jù)研究的一個(gè)重要任務(wù)。

表2 高通量測(cè)序下各種組學(xué)所使用的技術(shù)

2.4 數(shù)據(jù)挖掘技術(shù)在生物大數(shù)據(jù)處理中的挑戰(zhàn)

面對(duì)高通量測(cè)序數(shù)據(jù)的爆發(fā)式增長(zhǎng)，傳統(tǒng)的數(shù)據(jù)挖掘算法和工具遭遇巨大的挑戰(zhàn)：如何建立智能學(xué)習(xí)數(shù)據(jù)庫(kù)系統(tǒng)；如何對(duì)生物大數(shù)據(jù)存儲(chǔ)訪問(wèn)和計(jì)算；如何進(jìn)行隱私保護(hù)；如何結(jié)合領(lǐng)域知識(shí)設(shè)計(jì)新的適用于生物大數(shù)據(jù)挖掘分析的算法和工具。具體來(lái)說(shuō)，面向生物學(xué)數(shù)據(jù)挖掘的數(shù)據(jù)挖掘技術(shù)主要有3個(gè)層次的挑戰(zhàn)。第一個(gè)挑戰(zhàn)是數(shù)據(jù)的訪問(wèn)和程序的運(yùn)算。因?yàn)榇髷?shù)據(jù)都是分布式存儲(chǔ)的，隨著數(shù)據(jù)量的增長(zhǎng)，如何建立一個(gè)有效的平臺(tái)，使分散存儲(chǔ)的數(shù)據(jù)能夠擺脫計(jì)算機(jī)內(nèi)存的限制和大數(shù)據(jù)處理的障礙，進(jìn)行分布式計(jì)算。第二個(gè)挑戰(zhàn)是不同的大數(shù)據(jù)有不同的語(yǔ)義和領(lǐng)域知識(shí)，如何能夠更好地挖掘語(yǔ)義和領(lǐng)域知識(shí)，為數(shù)據(jù)所有者和消費(fèi)者服務(wù)。第三個(gè)挑戰(zhàn)集中在算法設(shè)計(jì)方面，生物大數(shù)據(jù)稀疏且具有各種各樣的混合數(shù)據(jù)，數(shù)據(jù)有不確定性、不完整性和多源性等特點(diǎn)，如何用數(shù)據(jù)融合技術(shù)進(jìn)行處理，并且挖掘出蘊(yùn)含其中的復(fù)雜和動(dòng)態(tài)信息；如何通過(guò)局部學(xué)習(xí)，得到一個(gè)反映全局問(wèn)題的融合模型[17]。

3 高通量DNA測(cè)序數(shù)據(jù)的生物信息學(xué)方法

隨著生物信息技術(shù)突飛猛進(jìn)地發(fā)展，越來(lái)越多的計(jì)算機(jī)和數(shù)學(xué)領(lǐng)域的專家加入生物信息學(xué)研究的隊(duì)伍，開(kāi)發(fā)出許多好用的生物信息學(xué)工具，使得生物學(xué)、醫(yī)學(xué)領(lǐng)域的專家可以利用這些先進(jìn)工具對(duì)生物大數(shù)據(jù)進(jìn)行分析，更準(zhǔn)確地揭示生物進(jìn)化的內(nèi)部規(guī)律，更好地解釋遺傳變異，為基礎(chǔ)醫(yī)學(xué)研究向醫(yī)學(xué)臨床應(yīng)用轉(zhuǎn)化提供新思路和新方法，取得了非常有意義的成果。但是NGS測(cè)序的樣本制備過(guò)程非常復(fù)雜，并且生成的序列難以處理，這給生物信息學(xué)專家?guī)?lái)了很大的挑戰(zhàn)。

3.1 高通量DNA測(cè)序數(shù)據(jù)的壓縮算法

NGS測(cè)序下的短讀序列的數(shù)據(jù)量呈爆炸性增長(zhǎng)，如果不對(duì)其進(jìn)行壓縮而直接存儲(chǔ)或傳輸會(huì)消耗巨大的硬件存儲(chǔ)設(shè)備，同時(shí)也會(huì)給網(wǎng)絡(luò)傳輸帶來(lái)很大的負(fù)擔(dān)。NGS測(cè)序數(shù)據(jù)有其自身的特點(diǎn)和規(guī)律，存在大量的信息冗余，傳統(tǒng)的數(shù)據(jù)壓縮算法并不能夠很好地壓縮DNA序列，這就需要開(kāi)發(fā)專門(mén)針對(duì)DNA序列的數(shù)據(jù)壓縮算法和工具。

近幾年，已經(jīng)研發(fā)了許多專門(mén)針對(duì)NGS數(shù)據(jù)的壓縮算法和工具，大多數(shù)是針對(duì)FASTQ格式的數(shù)據(jù)。根據(jù)DNA序列是否有參考基因組，壓縮方法分為有參考基因組的壓縮和無(wú)參考基因組的壓縮。有參考基因組的數(shù)據(jù)壓縮是利用參考基因組和短讀序列的差異信息來(lái)進(jìn)行壓縮。這種方法第一步先把短讀映射到參考基因組，記錄每條短讀在參考基因組上的位置以及與參考基因組的差異信息，然后再采用高效編碼方式存儲(chǔ)這些記錄，實(shí)現(xiàn)數(shù)據(jù)壓縮。其代表算法有DNAzip[18]、BWB[19]、SlimGene[20]、GRS[21]、mZIP[22]、NGC[23]、samcomp[24]等。由于同源物種基因組之間具有高度相似性，這種壓縮通常能達(dá)到很高的壓縮比，但這種方法有明顯的局限性，有些測(cè)序數(shù)據(jù)（如宏基因數(shù)據(jù)、從頭測(cè)序數(shù)據(jù)）并不存在現(xiàn)成的參考基因組，因此無(wú)法使用此算法；另外，該方法對(duì)于參考基因組依賴性太強(qiáng)，壓縮和解壓縮都需要相同參考基因組，這樣參考基因組必須事先保存在本地，如果參考基因組缺失將直接影響壓縮數(shù)據(jù)的使用。

無(wú)參考基因組的數(shù)據(jù)壓縮方法通常采用兩步法，首先最大限度地識(shí)別冗余DNA序列，然后再利用通用的壓縮方法（如gzip、bzip2）進(jìn)行處理。其代表算法工具有Beetl[25]、SCALCE[26]、SRComp[27]和ORCOM[28]。Beetl采用Burrows Wheeler變換算法，識(shí)別冗余；SCALCE采用局部一致性技術(shù)方法排序短讀序列，識(shí)別關(guān)鍵子串；SRComp采用burstsort排序的方法，使相同的字符串聚集在一起，然后再采用不同的編碼方式對(duì)其進(jìn)行編碼。ORCOM采用并行的Minimizers算法壓縮reads中的重疊區(qū)域（overlap）。另一種新穎的無(wú)參考基因組的數(shù)據(jù)壓縮方法是基于拼接的方法，代表算法有Quip[29]。Quip方法采用拼接的方式，用一小部分短讀拼接成疊連群作為臨時(shí)參考基因組，然后利用基于參考基因組的壓縮方法進(jìn)行壓縮。

盡管高通量測(cè)序數(shù)據(jù)的壓縮研究已取得一定成果，但其在計(jì)算資源、壓縮算法方面仍面臨巨大挑戰(zhàn)。隨著DNA測(cè)序數(shù)據(jù)量的增大，對(duì)計(jì)算資源的要求也越來(lái)越大，處理時(shí)間過(guò)長(zhǎng)是DNA測(cè)序數(shù)據(jù)分析最重要的問(wèn)題。另外，如何利用高通量測(cè)序技術(shù)產(chǎn)生有意義的冗余信息、采用并行化策略和基于索引的壓縮方法、建立統(tǒng)一的數(shù)據(jù)質(zhì)量評(píng)價(jià)標(biāo)準(zhǔn)等，都是重要的研究方向。

3.2 高通量DNA測(cè)序的序列拼接

由于測(cè)序技術(shù)的限制，新一代測(cè)序的讀長(zhǎng)較短（30～500 bp）[30]，測(cè)序所得序列無(wú)法滿足大多數(shù)序列分析的需要[31]，因此序列拼接成為基因組學(xué)研究中一個(gè)重要的環(huán)節(jié)。所謂序列拼接，是指將測(cè)序得到的短序列片段利用計(jì)算的方法拼接成較長(zhǎng)的連續(xù)序列片段（contig）或者中間帶有空隙的長(zhǎng)序列片段（scaffold）乃至整段基因組序列的方法。

序列拼接包括兩種不同的策略：從頭（D e N o v o）拼接的方法和對(duì)照（comparative）拼接的方法[32]。從頭拼接是指在沒(méi)有任何基因組序列參照的前提下，構(gòu)建全新基因組序列的策略，而對(duì)照拼接是指在參照基因組序列的指導(dǎo)下進(jìn)行的基因組序列的拼接。對(duì)照拼接適用于存在參照基因組序列的拼接，比如重測(cè)序項(xiàng)目中的序列拼接，而對(duì)于全新物種的大規(guī)模全基因組測(cè)序以及宏基因組測(cè)序項(xiàng)目主要使用從頭拼接。

拼接算法的主要挑戰(zhàn)來(lái)源于基因組中的重復(fù)序列片段。在不同區(qū)域的兩個(gè)完全一致的重復(fù)片段無(wú)法通過(guò)計(jì)算的方式來(lái)辨別。對(duì)于相似但不完全一致的重復(fù)片段，可以通過(guò)提高序列比對(duì)的相似度閾值區(qū)分不同的復(fù)本，這種方法一般還涉及對(duì)reads中測(cè)序錯(cuò)誤的估計(jì)[33]。重復(fù)片段的區(qū)分一般需要借助于reads或是mate-pair的跨越。所謂的mate-pair是指測(cè)序時(shí)從一段長(zhǎng)度已知的片段兩端測(cè)得的一對(duì)reads。對(duì)于reads來(lái)說(shuō)，如果reads的中間是重復(fù)序列，而兩端都有足夠長(zhǎng)的唯一片段，則可以區(qū)分中間的重復(fù)片段，這種方法針對(duì)短的重復(fù)片段有效，一般在k-mer圖算法中使用。對(duì)于mate-pair來(lái)說(shuō)，如果matepair分別處于重復(fù)序列的兩端，也可以指導(dǎo)正確的拼接路徑，而且mate-pair比reads更長(zhǎng)，因此可以區(qū)分更長(zhǎng)的重復(fù)片段。高的測(cè)序深度有利于重復(fù)片段的區(qū)分，因?yàn)楦叩臏y(cè)序深度可能提供更多的reads或者mate-pair跨越重復(fù)片段。對(duì)于新一代測(cè)序中短序列的拼接，重復(fù)片段的區(qū)分更加困難，因?yàn)閞eads更短，更多的重復(fù)片段無(wú)法通過(guò)reads來(lái)區(qū)分，因此提高測(cè)序深度和使用mate-pair尤為重要。

測(cè)序錯(cuò)誤也給重復(fù)片段問(wèn)題的解決增加了難度。因?yàn)槠唇铀惴ū仨氁驗(yàn)闇y(cè)序錯(cuò)誤而接受不完全一致的重疊，以免錯(cuò)漏了真實(shí)的重疊。然而對(duì)測(cè)序錯(cuò)誤的容忍又增加了拼接的假陽(yáng)性。更多不完全一致的重復(fù)片段會(huì)對(duì)算法造成麻煩。另外，序列拼接需要考慮的一個(gè)問(wèn)題是計(jì)算時(shí)間上的復(fù)雜度問(wèn)題，尤其對(duì)于reads數(shù)量越來(lái)越多的大規(guī)模測(cè)序數(shù)據(jù)。例如，為了提高拼接效率，所有的拼接軟件都在不同程度地以不同方式使用k-mer的概念。很直觀的一個(gè)結(jié)論是，reads之間的重疊區(qū)域必然共同享有k-mer。而對(duì)共享k-mer的搜索顯然要比計(jì)算序列比對(duì)簡(jiǎn)單得多。因此，幾乎所有的拼接算法都涉及對(duì)k-mer的計(jì)算。

理論上，序列拼接屬于一個(gè)NP難的問(wèn)題，尚無(wú)一個(gè)蓋棺定論的解答方法?，F(xiàn)有的拼接算法只能通過(guò)一系列復(fù)雜的推斷性質(zhì)的步驟來(lái)獲得近似的“解答”。這些算法仍有局限性，例如拼接結(jié)果錯(cuò)誤、拼接序列連續(xù)性差、計(jì)算時(shí)間長(zhǎng)、內(nèi)存消耗量大等。因此，序列拼接算法仍有很大的改進(jìn)空間。另外，測(cè)序技術(shù)的不斷變化和改進(jìn)，使得新數(shù)據(jù)對(duì)序列拼接不斷提出新的要求，以更好地適應(yīng)新數(shù)據(jù)的特點(diǎn)。

3.3 高通量測(cè)序下宏基因組的基因預(yù)測(cè)方法

基于高通量測(cè)序的宏基因組學(xué)研究給環(huán)境相關(guān)微生物的研究帶來(lái)了新的機(jī)遇。隨著越來(lái)越多的各種生態(tài)環(huán)境中宏基因組序列被測(cè)定并公開(kāi)，有效的宏基因組數(shù)據(jù)分析和功能預(yù)測(cè)軟件被開(kāi)發(fā)與應(yīng)用，這些都大大推動(dòng)了宏基因組學(xué)的發(fā)展。目前研究基因預(yù)測(cè)的方法主要有兩類：一類是基于序列相似性的預(yù)測(cè)方法，基于已知的基因序列通過(guò)搜索相似度較高的序列進(jìn)行預(yù)測(cè)；另一類是基于統(tǒng)計(jì)學(xué)模型的預(yù)測(cè)方法，即利用數(shù)學(xué)統(tǒng)計(jì)模型進(jìn)行基因預(yù)測(cè)，從已知的DNA序列中訓(xùn)練出統(tǒng)計(jì)學(xué)模型，應(yīng)用到宏基因組的測(cè)序結(jié)果上進(jìn)行預(yù)測(cè)。

（1）基于序列相似性比較的方法

序列比對(duì)是生物信息學(xué)的基礎(chǔ)，其基本問(wèn)題是比較兩個(gè)或兩個(gè)以上序列之間的相似性。兩個(gè)序列比對(duì)已有發(fā)展成熟的動(dòng)態(tài)規(guī)劃（dynamic programming）算法和在此基礎(chǔ)上發(fā)展起來(lái)的工具包BLAST[34]和FASTA[35]。事實(shí)上，在基于比對(duì)的方法中，高通量測(cè)序所得的序列較短，而這種短序列直接進(jìn)行比對(duì)的效果往往不理想，并且大量的原始數(shù)據(jù)進(jìn)行比對(duì)會(huì)耗費(fèi)很多時(shí)間，因此需要在比對(duì)前進(jìn)行序列拼接，將其拼接成較長(zhǎng)的序列，提高分析效率和分析效果[36]。由于必須與已知基因序列進(jìn)行相似性比較，故這種方法很難發(fā)現(xiàn)新基因。

基于序列相似性比較的高通量測(cè)序的宏基因組數(shù)據(jù)的應(yīng)用非常多。2010年，華大基因在Nature發(fā)表文章，對(duì)人體腸道微生物基因組研究計(jì)劃（MetaHIT）進(jìn)行了總結(jié)[37]。該計(jì)劃為研究人體腸道微生物群落與人類健康之間的關(guān)系，采集了124個(gè)歐洲人的糞便樣本，其中包括25個(gè)炎癥性腸?。╥nflammatory bowel disease，IBD）患者和99個(gè)健康志愿者的樣本，并用Illumina測(cè)序平臺(tái)進(jìn)行測(cè)序，產(chǎn)生了567.7 GB的測(cè)序數(shù)據(jù)，并對(duì)序列拼接、注釋、功能基因的分類、多態(tài)性分析等進(jìn)行了研究。2012年，華大基因在Nature發(fā)表了一篇研究人體腸道微生物與Ⅱ型糖尿病之間關(guān)系的文章[38]。該研究收集了345個(gè)中國(guó)人的腸道微生物樣本，用Illumina測(cè)序平臺(tái)對(duì)其進(jìn)行了深度測(cè)序，并在基因組關(guān)聯(lián)研究（genome wide association studies，GWAS）的基礎(chǔ)上開(kāi)發(fā)了一種全基因組相關(guān)聯(lián)研究（metagenome wide association studies，MGWAS）的方法，對(duì)Ⅱ型糖尿病與腸道微生物失調(diào)之間的關(guān)系進(jìn)行了深入研究。人體腸道中絕大多數(shù)種類的微生物是難以培養(yǎng)的，只有運(yùn)用宏基因組學(xué)技術(shù)才能研究人類腸道中的所有微生物群落，進(jìn)而了解人類腸道中細(xì)菌的物種分布。

（2）基于序列內(nèi)容統(tǒng)計(jì)特征的方法

基于序列內(nèi)容統(tǒng)計(jì)特征的基因預(yù)測(cè)方法一般是建立在密碼子的編碼區(qū)和非編碼區(qū)有不同相對(duì)出現(xiàn)頻率的基礎(chǔ)上的。除了一個(gè)區(qū)域堿基組成的特征外，基因長(zhǎng)度分布、CG含量、基因重疊區(qū)域的特征等因素也常被用于基因預(yù)測(cè)中。根據(jù)DNA序列中編碼蛋白質(zhì)區(qū)域和非編碼區(qū)域內(nèi)容統(tǒng)計(jì)特征的差別，建立其學(xué)習(xí)模型，可以有效地進(jìn)行基因預(yù)測(cè)。在單個(gè)基因組上具有代表性的方法包括采用馬爾科夫模型的GeneMark[39-41]系列、Glimmer[42,43]系列、FGENESB[44]和MED[45,46]系列。GeneMark對(duì)原核生物、真核生物和病毒均能進(jìn)行基因預(yù)測(cè)。Glimmer被廣泛應(yīng)用于微生物的基因預(yù)測(cè)。FGENESB主要用于細(xì)菌基因組的基因自動(dòng)預(yù)測(cè)和注釋。MED是筆者所在課題組開(kāi)發(fā)的一款基于多元熵距離法的原核生物基因預(yù)測(cè)算法，該算法的基礎(chǔ)為開(kāi)放閱讀框（ORF）和翻譯起始位點(diǎn)（TIS）的綜合統(tǒng)計(jì)模型。MED2.0在對(duì)DNA的GC核苷酸含量高的細(xì)菌基因組和古細(xì)菌基因組的基因預(yù)測(cè)上具有明顯優(yōu)勢(shì)，之后又推出了MED2.1，提高了預(yù)測(cè)精度，達(dá)到了國(guó)際水平。

針對(duì)宏基因組序列的研究，研究人員開(kāi)發(fā)了一系列宏基因組預(yù)測(cè)算法（見(jiàn)表33）。宏基因組預(yù)測(cè)算法借鑒了傳統(tǒng)的基于單基因組的基因預(yù)測(cè)方法，只是對(duì)原始數(shù)據(jù)增加了預(yù)處理的步驟。例如，MetaGUN算法基于序列組成的統(tǒng)計(jì)特征對(duì)輸入序列進(jìn)行分類，對(duì)同一類中的序列使用相同的統(tǒng)計(jì)模型刻畫(huà)，然后分別獨(dú)立地進(jìn)行基因預(yù)測(cè)，在模擬宏基因序列測(cè)試集和在兩個(gè)人體腸道微生物的真實(shí)數(shù)據(jù)上的測(cè)試表明，MetaGUN在發(fā)現(xiàn)新基因方面更具潛力。MetaGeneMark同時(shí)使用細(xì)菌—古細(xì)菌和嗜溫細(xì)菌—嗜熱細(xì)菌兩套模型進(jìn)行預(yù)測(cè)。FragGeneScan適用于有測(cè)序錯(cuò)誤的宏基因組序列。

近年來(lái)，專門(mén)針對(duì)宏基因組序列的基因預(yù)測(cè)方法目前面臨著新的挑戰(zhàn)，基于序列相似性比較的方法，使用BLAST系統(tǒng)工具對(duì)已知數(shù)據(jù)庫(kù)進(jìn)行相似性搜索，依賴性強(qiáng)，無(wú)法發(fā)現(xiàn)新基因?；诮y(tǒng)計(jì)建模的預(yù)測(cè)算法運(yùn)行速度快，在保證高特異性的條件下能獲得更高的敏感性。宏基因組序列來(lái)源于繁雜且大多為未知的物種，微生物中已知的細(xì)菌和古細(xì)菌只占全世界存在量的10%；同時(shí)高通量測(cè)序的宏基因組DNA序列很短，存在大量不完整基因，無(wú)法在單個(gè)序列片斷上完成自學(xué)習(xí)，為統(tǒng)計(jì)建模所能提供的信息有限；另外，如何把分析結(jié)果和已知的數(shù)據(jù)庫(kù)（Greengenes[55]、SILVA[56]等）結(jié)合起來(lái)、如何進(jìn)一步研究生物體之間以及生物體和環(huán)境之間的相互作用等，都成為亟待解決的問(wèn)題。

表3 宏基因組基因預(yù)測(cè)算法

4 結(jié)束語(yǔ)

高通量測(cè)序技術(shù)奠定了生物信息學(xué)的“大數(shù)據(jù)”基礎(chǔ)，面對(duì)如潮水般的基因序列數(shù)據(jù)，給后續(xù)基因組分析方法的研究和工具的發(fā)展帶來(lái)了巨大挑戰(zhàn)。本文總結(jié)討論了高通量測(cè)序數(shù)據(jù)的基因組分析及生物信息學(xué)方法。目前，基因組生物信息學(xué)研究正面臨從傳統(tǒng)的全基因組序列分析到當(dāng)前基于短讀的序列片段（含contigs）分析；從傳統(tǒng)的單個(gè)物種的全基因組序列分析到當(dāng)前多個(gè)物種混雜的序列片段數(shù)據(jù)集的分析；從本地計(jì)算機(jī)運(yùn)算分析到未來(lái)適應(yīng)“云計(jì)算”模式的遠(yuǎn)程、快速運(yùn)算分析這幾方面發(fā)展。面對(duì)如此快速的發(fā)展，現(xiàn)有的生物信息學(xué)方法和工具已經(jīng)不能滿足如此大量的數(shù)據(jù)資料的需求，只有進(jìn)一步發(fā)展出優(yōu)秀的生物信息學(xué)方法和工具，才能更好地利用高通量測(cè)序技術(shù)的優(yōu)勢(shì)和應(yīng)用價(jià)值。

[1] S C H US T ER S C. N e x t-g e n e r a t i o n sequencing transforms today’s biology[J]. Nature Methods, 2008, 5(1): 16-18.

[2] SANGER F, NICKLEN S, COULSON A R. DNA sequencing with chain-terminating inhibitors[J]. Proceeding of the National Academy of Sciences, 1977, B7(12): 5463-5467.

[3] SHENDURE J, JI H. Next-generation DNA sequencing[J]. Nature Biotechnology, 2008, 26(10): 1135-1145.

[4] HIGGINS G. Human Genomes and Big Data Challenges[R]. Mason: AssureRx Health Inc, 2013.

[5] WARD R M, SCHMIEDER R, HIGHNAM G, et al. Big data challenges and opportunities in highthrough-put sequencing[J]. Systems Biomedicine, 2013, 1(1): 29-34.

[6] DUNHAM I, BIRNEY E, LA JOIE B R, et al. An integrated encyclopedia of DNA elements in the human genome[J]. Nature, 2012, 489(7414): 57-74.

[7] COLLINS F S, BARKER A D. Mapping the cancer genome[J]. Scientific American, 2007, 296(3): 50-57.

[8] HAYDEN E C. International genome project launched[J]. Nature, 2008, 451(7177): 378-389.

[9] GEVERS D, KNIGHT R, PETROSINO J F, et al. The human microbiome project: a community resource for the healthy human microbiome[J]. PLoS Biology, 2012, 10(8): e1001377.

[10] HAUSSLER D, O’BRIEN S J, RYDER O A, et al. Genome 10K: a proposal to obtain whole-genome sequence for 10 000 vertebrate species[J]. The Journal of Heredity, 2008, 100(6): 659-674.

[11] O’ROAK B J, VIVES L, GIRIRAJAN S, et al. Sporadic autism exomes reveal a highly interconnected protein network of de novo mutations[J]. Nature, 2012, 485(7397): 246-250.

[12] EHRLICH S D. MetaHIT: the European union project on metagenomics of the human intestinal tract[M]// Metagenomics of the Human Body. New York: Springer, 2011: 307-316.

[13] L E G R A I N P, A E B E R S O L D R, A R C H A KO V A, e t a l. T h e h u m a n proteome project: current state and future direction[J]. Molecular & Cellular Proteomics, 2011, 10(7): M111. 009993.

[14] G I L B E R T J A, M E Y E R F, ANTONOPOULOS D, et al. Meeting report: the terabase metagenomics workshop and the vision of an earth microbiome project[J]. Standards in Genomic Sciences, 2010, 3(3): 243.

[15] R O B I N S O N G E, H A C K E T T K J, PURCELL M M, et al. Creating a buzz about insect genomes[J]. Science, 2011, 331(6023): 1386.

[16] JOLY Y, DOVE E S, KNOPPERS B M, et al. Data sharing in the post-genomic world: the experience of the international cancer genome consortium (ICGC) data access compliance office (DACO)[J]. PLoS Comput Biol, 2012, 8(7): e1002549.

[17] WU X D, ZHU X Q. Data mining with big data[J]. IEEE Transactions on Knowledge and Data Engineering, 2014, 26(1): 97-108.

[18] C H R I S T L E Y S, L U Y, L I C, e t a l. Human genomes as email attachments[J]. Bioinformatics, 2009, 25(2): 274-275.

[19] BRADON M C, WALLACE D C, BALDI P. Data structures and compression algorithms for genomic sequence data[J]. Bioinformatics, 2009, 25(14): 1731-1738.

[20] K O Z A N I T I S C, S A U N D E R S C, K RUGLYAK S, e t al. Co m p r e ssin g genomic sequence fragments using SlimGene[J]. Journal of Computational Biology, 2011, 18(3): 401-413.

[21] WANG C, ZHANG D. A novel compression tool for efficient storage of genome resequencing data[J]. Nucleic Acids Research, 2011, 39(7): e45.

[22] F R I T Z M H Y, L E I N O N E N R, COCHRANE G, et al. Efficient storage of high throughput DNA sequencing data using reference-based compression[J]. Genome Research, 2011, 21(5): 734-740.

[23] MILLER J R, KOREN S, SUT TON G. Assembly algorithms for next-generation sequencing data[J]. Genomics, 2010, 95(6): 315-327.

[24] B O N F I E L D J K, M A H O N E Y M V. Compression of FASTQ and SAM format sequencing data[J]. Plos One, 2013, 8(3): 1453-1456.

[25] COX A J, BAUER M J, JAKOBI T, et al. Large-scale compression of genomic sequence databases with the Burrows-Wheeler transform[J]. Bioinformatics, 2012, 28(11): 1415-1419.

[26] H A C H F, N U M A N A G I ? I, A L K A N C, et al. SCALCE: boosting sequence compression algorithms using locally consistent encoding[J]. Bioinformatics, 2012, 28(23): 3051-3057.

[27] SELVA J J, CHEN X. SRComp: short read sequence compression using burstsort and Elias omega coding[J]. PloS One, 2013, 8(12): e81414.

[28] P AT R O R, K I N G S F O R D C. D a t adependent bucketing improves referencefree compression of sequencing reads[J]. Bioinformatics, 2015: btv248.

[29] J O N E S D C, RUZ ZO W L, PEN G X, et al. Compression of next-generation sequencing reads aided by highly efficient de novo assembly[J]. Nucleic Acids Research, 2012, 40(22): e171.

[30] METZKER M L. Applications of nextgeneration sequencing technologies the next generation[J]. Nature Reviews Genetics, 2010, 11(1): 31-46.

[31] WOOLEY C, GODZIK A, FRIEDBERG I. A primer on metagenomics[J]. PLoS Comput Biol, 2010, 6(2): e1000667.

[32] POP M, PHILLIPPY A, DELCHER A L, et al. Comparative genome assembly[J]. Briefings in Bioinformatics, 2004, 5(3): 237-248.

[33] KECECIOGLU J, JU J. Separating repeats in DNA sequence assembly[C]// The 5th Annual International Conference on Computational Biology, April 22-25,2001, Montreal, Canada. [S.l.:s.n.], 2001: 176-183.

[34] P R I D E D T, M E I N E R S M A N N R J, WASSENAAR T M, et al. Evolutionary implications of microbial genome tetranucleotide frequency biases[J]. Genome Research, 2003, 13(2): 145-158.

[35] WU Y W, YE Y. A novel abundance-based algorithm for binning metagenomic sequences using l-tuples[J]. Journal of Computational Biology, 2011, 18(3): 523-534.

[36] PRAKASH T, TAYLOR T D. Functional assignment of metagenomic data: challenges and applications[J]. Briefings in Bioinformatics, 2012, 13(6): 711-727.

[37] QIN J, LI R, RAES J, et al. A human gut microbial gene catalogue established by metagenomic sequencing[J]. Nature, 2010, 464(7285): 59-65.

[38] QIN J, LI Y, CAI Z, et al. A metagenomewide association study of gut microbiota in type 2 diabetes[J]. Nature, 2012, 490(7418): 55-60.

[39] B O R O D O V S K Y M, M C I N I N C H J. GENMARK: parallel gene recognition for both DNA strands[J]. Computers & Chemistry, 1993, 17(2): 123-133.

[40] L U K A S H I N A, B O R O D O V S K Y M. GeneMark.hmm: new solutions for gene finding[J]. Nucleic Acids Research, 1998, 26(4): 1107-1115.

[41] B E S E M E R J, L O M S A D Z E A, BORODOVSKY M. GeneMarks: a selftraining method for prediction of gene starts in microbial genomes. Implications for finding sequence motifs in regulatory regions[J]. Nucleic Acids Research, 2001, 29(12): 2607-2618.

[42] SALZBERG S L, DELCHER A L, KASIF S, et al. Microbial gene identification using interpolated Markov models[J]. Nucleic Acids Research, 1998, 26(2): 544-548.

[43] DELCHER A L, BRATKE K A, POWERS E C, et al. Identifying bacterial genes and endosymbiont DNA with Glimmer[J]. Bioinformatics, 2007, 23(6): 673-679.

[44] FRIGAARD N U, MARTIMEZ A, MINCER T J, et al. Proteorhodopsin lateral gene transfer between marine planktonic bacteria and archaea[J]. Nature, 2006, 439(7078): 847-850.

[45] OUYANG Z, ZHU H, WANG J, et al. Multivariate entropy distance method for prokaryotic gene identification[J]. Journal of Bioinformatics and Computational Biology, 2004, 2(2): 353-373.

[46] ZHU H Q, HU G Q, YANG Y F, et al. MED: a new non-super vised gene prediction algorithm for bacterial and archaeal genomes[J]. BMC Bioinformatics, 2007, 8(1): 97.

[47] NOGUCHI H, TANIGUCHI T, ITOH T. MetaGeneAnnotator: detecting speciesspecific patterns of ribosomal binding site for precise gene prediction in anonymous prokaryotic and phage genomes[J]. DNA Research, 2008, 15(6): 387-396.

[48] H O F F K J, L I N G N E R T, M E I N I C K E P, et al. Orphelia: predicting genes in metagenomic sequencing reads[J]. Nucleic Acids Research, 2009, 37(suppl 2): W101-W105.

[49] ZHU W, LOMSADZE A, BORODOVSKY M. Ab initio gene identification in metagenomic sequences[J]. Nucleic Acids Research, 2010, 38(12): e132.

[50] RHO M, TANG H, YE Y. FragGeneScan: predicting genes in short and error-prone reads[J]. Nucleic Acids Research, 2010, 38(20): e191.

[51] KELLE Y D R, LIU B, DELCHER A L, et al. Gene prediction with Glimmer for metagenomic sequences augmented by classification and clustering[J]. Nucleic Acids Research, 2012, 40(1): e9.

[52] HYATT D, LOCASCIO P F, HAUSER L J,et al. Gene and translation initiation site prediction in metagenomic sequences[J]. Bioinformatics, 2012, 28(17): 2223-2230.

[53] WANG Y, LEUNG H C M, YIU S M, et al. MetaCluster 5.0: a two-round binning approach for metagenomic data for lowabundance species in a noisy sample[J]. Bioinformatics, 2012, 28(18): i356-i362.

[54] L I U Y, G U O J, H U G, e t a l. G e n e prediction in metagenomic fragments based on the SVM algorithm[J]. BMC Bioinformatics, 2013, 14(suppl 5): S12.

[55] D E S A N T I S T Z, H U G E N H O LT Z P, L A R S E N N, e t a l. G r e e n g e n e s, a chimera-checked 16S rRNA gene database and workbench compatible with ARB[J]. Applied and Environmental Microbiology, 2006, 72(7): 5069-5072.

[56] PRUESSE E, QUA ST C, K NIT TEL K, et al. SILVA: a comprehensive online resource for quality checked and aligned ribosomal RNA sequence data compatible with ARB[J]. Nucleic Acids Research, 2007, 35(21): 7188-7196.

Bioinformatics methods for high-throughput DNA sequencing data

ZHAN Xiaojuan1, YAO Dengju2, ZHU Huaiqiu3
1. College of Computer Science and Technology, Heilongjiang Institute of Technology, Harbin 150050, China
2. School of Software, Harbin University of Science and Technology, Harbin 150040, China 3. Department of Biomedical Engineering, Peking University, Beijing 100871, China

DNA sequence data generated by high-throughput sequencing technology is short in length, and the amount of data is enormous. The challenges and opportunities of the big data in high-throughput sequencing environment were analyzed. The data compression, the assembly of metagenomic sequence data, and algorithms and tools of metagenomic sequence data analysis also were summarized and discussed. Finally, the future of the study on short read DNA sequence data in high-throughput sequencing environment was discussed.

high-throughput DNA sequencing, bioinformatics, short read sequence data compression, short read sequence data splicing, short read sequence data analysis

TP399

A

10.11959/j.issn.2096-0271.2016021

詹曉娟（1978-），女，黑龍江工程學(xué)院講師，主要研究方向?yàn)閿?shù)據(jù)挖掘、機(jī)器學(xué)習(xí)、生物信息。

姚登舉（1980-），男，哈爾濱理工大學(xué)副教授，主要研究方向?yàn)閿?shù)據(jù)挖掘、機(jī)器學(xué)習(xí)、生物信息。

朱懷球（1970-），男，北京大學(xué)教授，主要研究方向?yàn)樯镝t(yī)學(xué)信息學(xué)和計(jì)算系統(tǒng)生物學(xué)。

2015-09-30

黑龍江省自然科學(xué)基金資助項(xiàng)目（No.F201313）；黑龍江省教育廳科學(xué)技術(shù)研究資助項(xiàng)目（No.12541124）；哈爾濱市科技創(chuàng)新人才資助項(xiàng)目（No.2013RFQXJ114）

Foundation Items: The Natural Science Foundation of Heilongjiang Province (No.F201313), The Foundation of Heilongjiang Province Educational Committee (No.12541124), The Harbin Special Funds for Technological Innovation Research of Heilongjiang Province of China (No.2013RFQXJ114)