高通量測序生物信息學(xué)分析

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

高通量測序分析

高通量測序，一次性對幾百萬到十億條DNA分子進(jìn)行并行測序，又稱為下一代測序技術(shù)，其使得可對一個(gè)物種的轉(zhuǎn)錄組和基因組進(jìn)行深入、細(xì)致、全貌的分析，所以又被稱為深度測序。主要包括：High-throughput Sequencing，Next Generation Sequencing，Deep Sequencing。

圖1 高通量測序流程

高通量測序應(yīng)用范圍廣泛：

1 DNA測序：全基因組de novo測序，基因組重測序，宏基因組測序，人類外顯子組捕獲測序。

2 RNA測序：轉(zhuǎn)錄組測序，小RNA測序，電子表達(dá)譜測序。

3 表觀基因組研究：ChIP-Seq，DNA甲基化測序。

基因組測序

基因組測序是對物種的基因組DNA打斷后進(jìn)行高通量測序，根據(jù)是否有已知基因組數(shù)據(jù)主要分為de novo全基因組測序和基因組重測序。De novo 基因組測序是對未知基因組序列的物種進(jìn)行基因組從頭測序，利用生物信息學(xué)分析手段對序列進(jìn)行拼接、組裝，從而獲得該物種的基因組圖譜。全基因組重測序是對已知基因組序列的物種進(jìn)行不同個(gè)體的基因組測序，并在此基礎(chǔ)上對個(gè)體或群體進(jìn)行差異性分析。

圖2 基因組測序策略

圖3 Paired-end原理

Paired-End方法，基因組打斷后，選擇一定長度（200－500bp）的序列連接兩端接頭進(jìn)行兩頭測序。Mate-end建庫較復(fù)雜，序列打斷后，選取一定長度序列（3-5kb），需先連接生物素，再環(huán)化，再打斷，生物素富集，連接兩端接頭進(jìn)行兩端測序。

基因組測序應(yīng)用生物信息學(xué)分析其結(jié)果，主要涵蓋以下幾方面。

1 數(shù)據(jù)產(chǎn)出處理：圖像識(shí)別與Base Calling\去除接頭序列、檢測與去除污染序列等；

2 基因組組裝：原始數(shù)據(jù)統(tǒng)計(jì)、測序深度分析、組裝結(jié)果統(tǒng)計(jì)等；

3 基因組注釋：Coding Gene注釋、RNA分類注釋、重復(fù)序列注釋等；

4 基因功能注釋：GO功能分類、Interpro功能分類等；

5 比較基因組及分子進(jìn)化分析：SNP/InDel/CNV檢測等。

宏基因組測序

宏基因組測序是對某一特定環(huán)境，如腸道、土壤、海水等中的所有微生物進(jìn)行基因組測序。通過此方法可對該環(huán)境中的微生物種類和優(yōu)勢物種進(jìn)行檢測，揭示微生物群落多樣性、種群結(jié)構(gòu)、進(jìn)化關(guān)系、功能活性、相互協(xié)作關(guān)系及與環(huán)境之間的關(guān)系 。自然環(huán)境中很多微生物無法分離培養(yǎng)，而此方法無需對微生物進(jìn)行分離培養(yǎng)。宏基因組測序方法現(xiàn)在有全基因組的宏基因組測序和16S/18S rRNA宏基因組測序。

1 全基因組的宏基因組測序

通過高通量測序技術(shù)，對環(huán)境樣品的總 DNA 直接進(jìn)行全基因組的宏基因組測序，能夠?qū)崿F(xiàn)微生物群落的物種分類研究、群落結(jié)構(gòu)、系統(tǒng)進(jìn)化、功能注釋以及物種間的代謝網(wǎng)絡(luò)研究，挖掘具有應(yīng)用價(jià)值的基因資源，開發(fā)新的微生物活性物質(zhì)。與傳統(tǒng)的 Sanger法相比，速度快，性價(jià)比高，周期短，單個(gè)樣品的測序量可以接近飽和。

宏基因組測序信息分析主要包括：拼接組裝，物種分類組成分析，基因預(yù)測和功能注釋，生成Profiling table，主成分分析（PCA），篩選與樣品分組顯著相關(guān)的因子，多樣品間比較分析等。

2 16S/18S rRNA宏基因組測序

16S/18S rRNA是微生物群落分析和細(xì)菌進(jìn)化研究以及分類研究最常用的靶分子，采用新一代測序技術(shù)，對16S/18S rDNA的可變區(qū)進(jìn)行測序分析，不需進(jìn)行克隆篩選，能全面的反映微生物群體的物種組成，真實(shí)的物種分布及豐度信息。

16S/18S rRNA測序信息分析主要包括：物種分類、物種豐度分析，OTU（Operational Taxonomic Units）分析，多樣性分析，系統(tǒng)進(jìn)化分析，多樣品間的比較分析等。

人類外顯子組捕獲測序

外顯子組是指全部外顯子區(qū)域的集合，該區(qū)域包含合成蛋白質(zhì)所需要的重要信息，涵蓋了與個(gè)體表型相關(guān)的大部分功能性變異。與全基因組重測序相比，外顯子組測序只需針對外顯子區(qū)域的DNA，覆蓋度更深、數(shù)據(jù)準(zhǔn)確性更高，更加簡便、經(jīng)濟(jì)、高效。

圖4 人類外顯子組捕獲測序原理

外顯子捕獲是指用外顯子芯片雜交，把基因組外顯子序列進(jìn)行捕獲，然后對所捕獲的序列進(jìn)行測序。現(xiàn)在常用外顯子芯片有Roche NimbleGen Sequence Capture 2.1M Human Exome Array和Agilent SureSelect Target Enrichment System（Human Exome）。

圖5 人類外顯子組捕獲測序分析流程

轉(zhuǎn)錄組測序

轉(zhuǎn)錄組即特定細(xì)胞在某一功能狀態(tài)下所能轉(zhuǎn)錄出來的所有RNA的總和，包括mRNA和非編碼RNA（Non-coding RNA）。

第二代測序系統(tǒng)可精確檢測單個(gè)堿基，并且不受到研究中先驗(yàn)信息的干擾，科研人員能夠快速地獲得某一物種特定器官或組織在某一狀態(tài)下幾乎所有mRNA轉(zhuǎn)錄本序列，從而能夠開展：UTRs區(qū)域界定、可變剪切研究、低豐度新轉(zhuǎn)錄本發(fā)現(xiàn)、融合基因鑒定、cSNP（編碼序列單核苷酸多態(tài)性）研究等。

圖6 轉(zhuǎn)錄組測序流程

圖7 無參考序列及有參考序列轉(zhuǎn)錄組測序流程

無參考序列轉(zhuǎn)錄組分析內(nèi)容包括：1 測序數(shù)據(jù)產(chǎn)量統(tǒng)計(jì)，數(shù)據(jù)成分和質(zhì)量評(píng)估；2 Contig及Scaffold長度分布；3 Unigene的長度分布和功能注釋，GO分類，Pathway分析，差異表達(dá)分析；4 蛋白功能預(yù)測與分類，差異表達(dá)基因GO富集和 Pathway富集分析。

有參考序列轉(zhuǎn)錄組分析內(nèi)容包括：1 基本數(shù)據(jù)統(tǒng)計(jì)，比對參考序列；2 序列在基因組上在分布；3 測序深度分析、隨機(jī)性評(píng)估和基因差異表達(dá)分析；4 新基因預(yù)測，基因可變剪接鑒定和基因融合鑒定等。

電子表達(dá)譜測序

電子表達(dá)譜測序（Digital Gene Expression， DGE）又稱為基因表達(dá)標(biāo)簽測序（mRNA tag profiling），又稱Tag-SAGE。其原理是通過兩種酶切作用對基因中一段長度為21nt的序列標(biāo)簽進(jìn)行測序。由于其測序只針對表達(dá)的基因進(jìn)行測序，產(chǎn)生的數(shù)據(jù)量相對較小，是研究基因表達(dá)譜的經(jīng)濟(jì)而快速的研究手段。是對特定處理?xiàng)l件下的全基因組基因表達(dá)譜進(jìn)行分析，已被廣泛用于功能基因組學(xué)和醫(yī)學(xué)等研究領(lǐng)域。

圖8 電子表達(dá)譜測序流程圖

電子表達(dá)譜分析內(nèi)容包括：圖像識(shí)別與原始?jí)A基數(shù)據(jù)讀取，去污染、去接頭，標(biāo)簽序列計(jì)數(shù)統(tǒng)計(jì)，基因組比對與統(tǒng)計(jì)，基因序列比對獲得所表達(dá)的基因列表，基因差異表達(dá)分析，聚類與表達(dá)類型分析，GO基因富集與分類分析，Pathway富集與分類分析，蛋白相互作用網(wǎng)絡(luò)分析，反義鏈轉(zhuǎn)錄本與新轉(zhuǎn)錄本檢測等。

小RNA測序

小 RNA是指長度在21-31nt的內(nèi)源性非蛋白質(zhì)編碼RNA，廣泛存在于高等和低等生物體內(nèi)，其對mRNA的轉(zhuǎn)錄及轉(zhuǎn)錄后水平等生命過程起到調(diào)節(jié)作用。現(xiàn)已知小RNA可歸納成三類：微RNA （miRNA），小干擾RNA（siRNA）和與piwi相互作用的RNA（piRNA）。

miRNA長度為21~24nt，產(chǎn)生于有典型莖環(huán)二級(jí)結(jié)構(gòu)的原轉(zhuǎn)錄本（pri-miRNA），在動(dòng)植物的目標(biāo)mRNA的降解與抑制方面發(fā)揮重要作用。siRNA，長度在19~25nt，產(chǎn)生于長雙鏈RNA，同樣在動(dòng)植物的目標(biāo)mRNA的降解與抑制方面發(fā)揮重要作用。piRNA，長度26~31nt，由與其相互作用的Piwi蛋白定義，目前研究表明其在配子形成的過程中起作用。

圖9 小RNA測序流程圖

小RNA測序分析內(nèi)容包括以下兩個(gè)主要方面：

1 基本分析：原始數(shù)據(jù)讀取，去接頭、去污染序列，長度分布統(tǒng)計(jì)，基因組比對等。

2 高級(jí)分析：Small RNA的分類注釋，miRNA / siRNA / piRNA的鑒定，新miRNA預(yù)測，差異表達(dá)miRNA聚類分析等。

ChIP-Seq

ChIP-Chromatin Immunoprecipitation染色質(zhì)免疫共沉淀，是指通過蛋白免疫相互作用，用抗體把和染色質(zhì)相互作用的蛋白，如組蛋白、轉(zhuǎn)錄因子等，沉淀下來，從而所獲取與其相結(jié)合的DNA序列。ChIP-Seq就是通過高通量測序?qū)hIP所得到的序列進(jìn)行測序，從而進(jìn)行蛋白和DNA相互作用相關(guān)研究。

ChIP-Seq分析內(nèi)容包括：

1 ChIP Sequencing結(jié)果與參考基因組序列進(jìn)行比對。

2 ChIP Sequencing reads 在全基因組的分布：唯一比對reads 在repeats 區(qū)域的分布，唯一比對reads 在各基因功能元件上的分布，唯一比對reads 的全基因組覆蓋深度。

3 全基因組peak 掃描：peak 掃描，peak 長度分布統(tǒng)計(jì)，peak 的全基因組覆蓋度，peak 在基因功能元件上的分布特征，

4 Peak相關(guān)基因分析篩選與GO功能富集分析。

5 多個(gè)樣品的差異分析：基于peak 相關(guān)基因的差異分析，基于peak 的差異分析。

圖10 ChIP-Seq分析流程

DNA甲基化測序

DNA甲基化對機(jī)體發(fā)育和基因表達(dá)有很重要的調(diào)控作用，和各種癌癥的發(fā)生和發(fā)展也有很大相關(guān)性，所以對基因組DNA甲基化進(jìn)行研究是一直來的熱門課題。通過高通量測序來研究DNA甲基化現(xiàn)在主要有兩種方法，一種是MeDIP，是通過與DNA甲基化位點(diǎn)相結(jié)合的抗體，進(jìn)行免疫共沉淀，然后對所得DNA序列進(jìn)行測序。另一種是Bisulfite Sequencing，是通過Bisulfite處理基因組來區(qū)分甲基化位點(diǎn)。

圖11 MeDIP 原理

MeDIP-Seq分析內(nèi)容包括：

1 MeDIP-seq 序列與參考序列的比對。

2 MeDIP-seq 序列數(shù)據(jù)在全基因組的分布趨勢： MeDIP-seq 測序reads 在全基因組上每條染色體上的分布，MeDIP-seq 測序reads 在全基因組上的覆蓋深度，MeDIP-Seq 測序reads 在CG、CHG和CHH位點(diǎn)上的覆蓋深度，MeDIP-Seq 測序reads 在不同基因功能元件上的分布，MeDIP-Seq 測序reads 在不同OE含量區(qū)域中的分布。

3 統(tǒng)計(jì)MeDIP-seq 序列富集區(qū)域（peak）的信息：Peak 掃描，Peak 長度數(shù)量及比例分布統(tǒng)計(jì)，單個(gè)樣品Peak 的OE含量分布統(tǒng)計(jì)，尋找Peak 相關(guān)基因，統(tǒng)計(jì)Peak 在不同基因功能元件上的分布。

4 基于Peak 的多樣品間差異分析：分析兩個(gè)樣品間的Peak 相關(guān)差異基因，對兩個(gè)樣品間的差異基因進(jìn)行GO功能富集分析及pathway 功能分析。

圖12 Bisulfite Sequencing原理

Bisulfite Sequencing分析內(nèi)容包括：

1 Bisulfite-seq序列與參考序列的比對。

2 深度和覆蓋度分析：C堿基有效測序深度的累積分布，不同reads 測序深度下的基因組覆蓋度。

3 計(jì)算C堿基的甲基化水平。

4 全基因組甲基化數(shù)據(jù)分布趨勢分析：甲基化C堿基中CG， CHG 與CHH的分布比例（H=A、C or T），CG、CHG和CHH中的所有C的甲基化水平，各條染色體中CG、CHG和CHH中C的甲基化水平（該項(xiàng)分析目前只用于“人”），統(tǒng)計(jì)不同基因區(qū)域內(nèi)CG、CHG和CHH中C的甲基化水平，不同基因元件區(qū)域中CG、CHG和CHH中C的甲基化水平，CHG，CHH中甲基化C附近的9bp序列的序列特征分析。

5 全基因組DNA 甲基化圖譜：染色體水平的甲基化C堿基的密度分布（該項(xiàng)分析目前只用于“人”），Scaffold的甲基化C堿基密度分布（該項(xiàng)分析針對物種：非人），不同基因組區(qū)域的甲基化分布特征，基因組不同轉(zhuǎn)錄元件中的DNA甲基化水平。

6 差異甲基化區(qū)域（DMR）分析。