本文關(guān)鍵詞：基因組二代測(cè)序數(shù)據(jù)的自動(dòng)化分析流程，由筆耕文化傳播整理發(fā)布。

Hereditas (Beijing) 2014年6月, 36(6): 618―624

技術(shù)與方法

基因組二代測(cè)序數(shù)據(jù)的自動(dòng)化分析流程

李文軻1, 李豐余1,2, 張思瑤1, 蔡斌1, 鄭娜1, 聶宇1, 周到2, 趙倩1

1. 中國(guó)醫(yī)學(xué)科學(xué)院, 北京協(xié)和醫(yī)學(xué)院, 國(guó)家心血管病中心, 阜外心血管病醫(yī)院, 心血管疾病國(guó)家重點(diǎn)實(shí)驗(yàn)室, 北京 100037; 2. 中南民族大學(xué)生物醫(yī)學(xué)工程學(xué)院, 武漢430074

摘要： 二代測(cè)序技術(shù)的發(fā)展對(duì)測(cè)序數(shù)據(jù)的處理分析提出了很高的要求。目前二代測(cè)序數(shù)據(jù)分析軟件很多, 但是

絕大多數(shù)軟件僅能完成單一的分析功能(例如：僅進(jìn)行序列比對(duì)或變異讀取或功能注釋等), 如何能正確高效地選擇整合這些軟件已成為迫切需求。文章設(shè)計(jì)了一套基于perl語(yǔ)言和SGE資源管理的自動(dòng)化處理流程來(lái)分析Illumina平臺(tái)基因組測(cè)序數(shù)據(jù)。該流程以測(cè)序原始序列數(shù)據(jù)作為輸入, 調(diào)用業(yè)界標(biāo)準(zhǔn)的數(shù)據(jù)處理軟件(如：BWA, Samtools, GATK, ANNOVAR等), 最終生成帶有相應(yīng)功能注釋、便于研究者進(jìn)一步分析的變異位點(diǎn)列表。該流程通過(guò)自動(dòng)化并行腳本控制流程的高效運(yùn)行, 一站式輸出分析結(jié)果和報(bào)告, 簡(jiǎn)化了數(shù)據(jù)分析過(guò)程中的人工操作, 大大提高了運(yùn)行效率。用戶只需填寫配置文件或使用圖形界面輸入即可完成全部操作。該工作為廣大研究者分析二代測(cè)序數(shù)據(jù)提供了便利的途徑。

關(guān)鍵詞：　二代測(cè)序; 自動(dòng)化數(shù)據(jù)分析; 流程; 變異檢測(cè)

Automatic analysis pipeline of next-generation sequencing data

Wenke Li1, Fengyu Li1, 2, Siyao Zhang1, Bin Cai1, Na Zheng1, Yu Nie1, Dao Zhou2, Qian Zhao1

1. State Key Laboratory of Cardiovascular Disease, Fuwai Hospital, National Center for Cardiovascular Disease, Chinese Academy of Medical Sciences and Peking Union Medical College, Beijing 100037, China;

2. College of Biomedical Engineering, South-Central University for Nationalities, Wuhan 430074, China

Abstract: The development of next-generation sequencing has generated high demand for data processing and analysis. Although there are a lot of software for analyzing next-generation sequencing data, most of them are designed for one specific function (e.g., alignment, variant calling or annotation). Therefore, it is necessary to combine them together for data analysis and to generate interpretable results for biologists. This study designed a pipeline to process Illumina sequencing data based on Perl programming language and SGE system. The pipeline takes original sequence data (fastq format) as input, calls the standard data processing software (e.g., BWA, Samtools, GATK, and Annovar), and finally outputs a list of annotated vari-ants that researchers can further analyze. The pipeline simplifies the manual operation and improves the efficiency by

收稿日期: 2013?09?07; 修回日期: 2014?01?20

基金項(xiàng)目:國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973計(jì)劃)項(xiàng)目(編號(hào)：2010CB529505)和中央高�；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金(編號(hào)：2012-XHGX02)資助

作者簡(jiǎn)介:李文軻, 碩士, 助理研究員, 研究方向：生物信息學(xué)。Tel: 010-88396071; E-mail: wksofia@gmail.com

通訊作者:趙倩, 博士, 副研究員, 研究方向：遺傳學(xué), 生物信息學(xué)。E-mail: zhaoqian82@gmail.com DOI: 10.3724/SP.J.1005.2014.0618 網(wǎng)絡(luò)出版時(shí)間: 2014-3-25 17:22:30

URL:

第6期

李文軻等: 基因組二代測(cè)序數(shù)據(jù)分析流程的自動(dòng)化實(shí)現(xiàn)

619

automatization and parallel computation. Users can easily run the pipeline by editing the configuration file or clicking the graphical interface. Our work will facilitate the research projects using the sequencing technology.

Keywords: next generation sequencing; automatic data analysis; pipeline; variantion detection

二代測(cè)序技術(shù)(Next-generation sequencing)大幅度降低了測(cè)序的時(shí)間和成本, 使得大規(guī)模測(cè)序逐漸成為常規(guī)的實(shí)驗(yàn)室研究和臨床檢測(cè)手段。測(cè)序產(chǎn)生的數(shù)據(jù)量急劇增加, 如何高效地分析這些數(shù)據(jù), 已成為迫切需要解決的問(wèn)題。目前, 分析序列信息的生物信息學(xué)軟件紛繁復(fù)雜, 但基本上每個(gè)軟件只能完成單一的分析功能, 實(shí)現(xiàn)一個(gè)完整的分析流程則需要對(duì)眾多軟件進(jìn)行整合, 而手動(dòng)串聯(lián)的效率往往不盡人意; 同時(shí), 這些軟件需要在Linux工作環(huán)境下以命令行運(yùn)行, 要求用戶具備較好的計(jì)算機(jī)背景; 另外, 即便一些實(shí)驗(yàn)室完成了分析流程的構(gòu)建, 他們往往不會(huì)公開(kāi)許多細(xì)節(jié), 新用戶仍然要從頭建起。本研究致力于構(gòu)建經(jīng)典的二代測(cè)序數(shù)據(jù)分析流程, 并實(shí)現(xiàn)各個(gè)環(huán)節(jié)的高效自動(dòng)化管理和分析, 減輕研究者前期的工作負(fù)擔(dān), 促進(jìn)相關(guān)領(lǐng)域進(jìn)一步對(duì)基因組測(cè)序研究項(xiàng)目的順利開(kāi)展。

(Hiseq2500), 甚至更高的250個(gè)堿基(Miseq)。測(cè)序讀長(zhǎng)不斷增加, 測(cè)序通量也在不斷上升。Illumina Hiseq2500是目前世界上通量最高的測(cè)序平臺(tái), 最多可以在大約10 d的時(shí)間內(nèi)測(cè)定3000億個(gè)堿基——即6～7個(gè)人類全基因組或60～80個(gè)人類全外顯子組的序列測(cè)定。

Illumina平臺(tái)以FASTQ格式[2]存儲(chǔ)測(cè)序結(jié)果, 這也是本流程的輸入文件。FASTQ文件記錄內(nèi)容包括所測(cè)的堿基讀段和質(zhì)量, 其數(shù)據(jù)格式如圖1所示。每條讀段(reads)占四行：第一行和第三行為讀段識(shí)別碼, 包含測(cè)序儀SN號(hào)、產(chǎn)生讀段的巷道(lane)、該讀段的編號(hào)等信息; 第二行為讀段測(cè)到的堿基序列; 第四行為所測(cè)到堿基的質(zhì)量分?jǐn)?shù), 每一個(gè)堿基都會(huì)對(duì)應(yīng)一個(gè)質(zhì)量分?jǐn)?shù)。

1.2 數(shù)據(jù)處理流程及軟件簡(jiǎn)介

目前測(cè)序數(shù)據(jù)處理軟件很多, 我們綜合考慮了適用性和效率, 整合出了一套標(biāo)準(zhǔn)的數(shù)據(jù)處理流程。具體來(lái)說(shuō), 獲得FASTQ格式的原始測(cè)序數(shù)據(jù)后, 需要對(duì)數(shù)據(jù)進(jìn)行以下處理：(1)使用BWA軟件把這些短序列和參考基因組進(jìn)行對(duì)比, 確定短序列在基因組上的位置, 把短序列組裝成完整的人類參考基因組; (2)使用Samtools軟件把這些短序列調(diào)整成按一定順序(1~22, X, Y, 其他)排列的序列, 并進(jìn)行數(shù)據(jù)格式的轉(zhuǎn)換; (3)使用Picard軟件把測(cè)序產(chǎn)生的冗

1 數(shù)據(jù)的獲取和分析流程的構(gòu)建

1.1 Illumina測(cè)序數(shù)據(jù)

本流程適用于Illumina測(cè)序平臺(tái)產(chǎn)出的雙端(Paired ends)測(cè)序數(shù)據(jù)。Illumina測(cè)序技術(shù)采用邊合成邊測(cè)序(Sequencing by synthesis, SBS)的方法, 早期的GA測(cè)序儀測(cè)序讀長(zhǎng)只有100個(gè)堿基, 隨著技術(shù)的改進(jìn), 目前的讀長(zhǎng)已經(jīng)增加到150個(gè)堿基

圖1 FASTQ格式示例

620

Hereditas (Beijing) 2014

第36卷

余信息和噪聲去掉; (4)使用GATK尋找樣本測(cè)序數(shù)據(jù)與參考基因組的差異, 列出這些差異點(diǎn); (5)使用Annovar對(duì)這些變異位點(diǎn)進(jìn)行功能注釋, 得到一個(gè)易于理解的變異位點(diǎn)列表。處理流程圖如圖 2 所示。

的序列比對(duì)軟件, 能高效地比對(duì)短序列和參考基因組, 找到短序列在參考基因組上的位置, 該軟件最長(zhǎng)支持至1 Mb的短序列比對(duì)。BWT方法通過(guò)B-W轉(zhuǎn)換將基因組序列按一定規(guī)則壓縮并建立索引, 再通過(guò)查找和回溯來(lái)定位讀段, 在查找時(shí)可通過(guò)堿基替代來(lái)實(shí)現(xiàn)允許的錯(cuò)配。采用Burrows-Wheeler轉(zhuǎn)換的代表軟件是Bowtie和BWA。比對(duì)結(jié)果如圖3所示：界面上方是測(cè)到的短序列, 下方是短序列所比對(duì)到的參考基因組。

1.2.2 SAM文件處理軟件Samtools

讀段定位到基因組后推薦采用SAM(Sequence Alignment/Map)格式或其二進(jìn)制版本BAM格式來(lái)存儲(chǔ)。二進(jìn)制版本可大大節(jié)省存儲(chǔ)空間, 但不能直接用普通文本編輯工具顯示。

SAM文件處理軟件Samtools可以很好的對(duì)SAM/BAM格式數(shù)據(jù)進(jìn)行操作, 因此, 本文使用它來(lái)進(jìn)行數(shù)據(jù)格式轉(zhuǎn)換和排序。

1.2.3 測(cè)序噪聲去除和測(cè)序數(shù)據(jù)評(píng)價(jià)軟件Picard

對(duì)組裝好的全基因組數(shù)據(jù), 需要將過(guò)度重復(fù)測(cè)

到的數(shù)據(jù)進(jìn)行剔除, 并且需要對(duì)數(shù)據(jù)質(zhì)量進(jìn)行評(píng)價(jià), Picard軟件可以很好地完成這兩項(xiàng)工作。 1.2.4 變異檢測(cè)軟件GATK

GATK主要用于在測(cè)序數(shù)據(jù)中尋找變異, 包括單堿基變異(SNV)、短插入缺失(INDEL), 是當(dāng)前業(yè)界用來(lái)尋找變異的主流軟件。

圖2 處理流程圖

1.2.1 讀句比對(duì)軟件BWA

BWA(Burrows-Wheeler Alignment tool)是基于Burrows-Wheeler變換(Burrows-Wheeler Transform)

圖3 比對(duì)結(jié)果示例

利用Broad Institute的IGV(Integrated Genomics Viewer)對(duì)數(shù)據(jù)進(jìn)行可視化, 圖4同。

第6期

李文軻等: 基因組二代測(cè)序數(shù)據(jù)分析流程的自動(dòng)化實(shí)現(xiàn)

621

變異指測(cè)序序列和參考序列的差異。如圖4所示, 參考序列上的堿基是胸腺嘧啶(T), 而測(cè)序數(shù)據(jù)上的堿基是鳥(niǎo)嘌呤(G), 說(shuō)明此處有一個(gè)T→G 的 突變。

1.2.5 變異注釋軟件ANNOVAR

ANNOVAR是一個(gè)用于高效注釋變異的工具。注釋信息包括變異所在的染色體, 開(kāi)始位置, 結(jié)束位置, 參考序列信息和觀察到的序列信息的列表。一個(gè)變異經(jīng)過(guò)ANNOVAR注釋之后, 其功能一目了然, 便于進(jìn)一步的生物學(xué)分析。

務(wù)正在進(jìn)行時(shí), Perl對(duì)它進(jìn)行監(jiān)控; 當(dāng)計(jì)算完成, Perl去查看它的輸出的計(jì)算結(jié)果, 并把結(jié)果作為下一個(gè)計(jì)算任務(wù)的輸入, 往計(jì)算節(jié)點(diǎn)上投放新的計(jì)算任務(wù)。如此循環(huán), 直到流程運(yùn)行完畢。

同時(shí), 由于每次運(yùn)行的樣本不同, 數(shù)據(jù)的輸入輸出位置也有差異。如果每處理一個(gè)新的樣本, 就要對(duì)流程源碼進(jìn)行大量修改, 將不利于流程的使用。為此, 本流程定義了一個(gè)配置文件(config file)。通過(guò)配置文件可以指定：流程處理的樣品名、數(shù)據(jù)輸入輸出路徑、參考序列文件, 甚至流程中涉及到的軟件的位置、軟件的運(yùn)行方式; 另外, 我們還提供了對(duì)流程中主要軟件參數(shù)的修改, 以滿足高級(jí)用戶需求。每次進(jìn)行一個(gè)新樣本的分析, 不需要修改主程序代碼, 只要為其創(chuàng)建一個(gè)配置文件, 主程序會(huì)自動(dòng)讀取配置文件, 生成相應(yīng)的執(zhí)行代碼。

流程文件構(gòu)成如圖5所示。

2 自動(dòng)化實(shí)現(xiàn)

2.1 基于Perl語(yǔ)言的流程設(shè)計(jì)

本數(shù)據(jù)處理流程主要使用Perl編程語(yǔ)言實(shí)現(xiàn)對(duì)各個(gè)軟件的高效串接和自動(dòng)化操作。一項(xiàng)計(jì)算任

圖4 單堿基突變示例

2.2 基于資源管理軟件(SGE)的并行設(shè)計(jì)

流程的運(yùn)行環(huán)境是計(jì)算機(jī)集群, 其有別于普通PC機(jī), 一般是由一臺(tái)管理主機(jī)來(lái)協(xié)調(diào)許多計(jì)算主機(jī)

來(lái)完成大型的計(jì)算任務(wù)。根據(jù)這樣的硬件特點(diǎn)來(lái)設(shè)計(jì)流程, 需要考慮以下兩個(gè)問(wèn)題：(1)如何讓眾多計(jì)

圖5 分析流程結(jié)構(gòu)

  本文關(guān)鍵詞：基因組二代測(cè)序數(shù)據(jù)的自動(dòng)化分析流程，由筆耕文化傳播整理發(fā)布。