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

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

1. 中國醫(yī)學科學院, 北京協(xié)和醫(yī)學院, 國家心血管病中心, 阜外心血管病醫(yī)院, 心血管疾病國家重點實驗室, 北京 100037;

2. 中南民族大學生物醫(yī)學工程學院, 武漢430074

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

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

1. 中國醫(yī)學科學院, 北京協(xié)和醫(yī)學院, 國家心血管病中心, 阜外心血管病醫(yī)院, 心血管疾病國家重點實驗室, 北京 100037;

2. 中南民族大學生物醫(yī)學工程學院, 武漢430074

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

二代測序; 自動化數(shù)據(jù)分析; 流程; 變異檢測

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

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

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

本流程適用于 Illumina測序平臺產(chǎn)出的雙端(Paired ends)測序數(shù)據(jù)。Illumina測序技術(shù)采用邊合成邊測序(Sequencing by synthesis, SBS)的方法[1],早期的GA測序儀測序讀長只有100個堿基, 隨著技術(shù)的改進, 目前的讀長已經(jīng)增加到 150個堿基(Hiseq2500), 甚至更高的 250個堿基(Miseq)。測序讀長不斷增加, 測序通量也在不斷上升。Illumina Hiseq2500是目前世界上通量最高的測序平臺, 最多可以在大約10 d的時間內(nèi)測定3000億個堿基——即6～7個人類全基因組或60～80個人類全外顯子組的序列測定。

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

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

目前測序數(shù)據(jù)處理軟件很多, 我們綜合考慮了適用性和效率, 整合出了一套標準的數(shù)據(jù)處理流程。具體來說, 獲得FASTQ格式的原始測序數(shù)據(jù)后,需要對數(shù)據(jù)進行以下處理：(1)使用 BWA軟件把這些短序列和參考基因組進行對比, 確定短序列在基因組上的位置, 把短序列組裝成完整的人類參考基因組; (2)使用Samtools軟件把這些短序列調(diào)整成按一定順序(1～22, X, Y, 其他)排列的序列, 并進行數(shù)據(jù)格式的轉(zhuǎn)換; (3)使用 Picard軟件把測序產(chǎn)生的冗余信息和噪聲去掉; (4)使用GATK尋找樣本測序數(shù)據(jù)與參考基因組的差異, 列出這些差異點; (5)使用Annovar對這些變異位點進行功能注釋, 得到一個易于理解的變異位點列表。處理流程圖如圖 2 所示。

圖1 FASTQ格式示例

圖2 處理流程圖

1.2.1 讀句比對軟件BWA

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

1.2.2 SAM文件處理軟件Samtools

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

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

1.2.3 測序噪聲去除和測序數(shù)據(jù)評價軟件Picard

對組裝好的全基因組數(shù)據(jù), 需要將過度重復測到的數(shù)據(jù)進行剔除, 并且需要對數(shù)據(jù)質(zhì)量進行評價, Picard軟件可以很好地完成這兩項工作[5]。

1.2.4 變異檢測軟件GATK

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

圖3 比對結(jié)果示例利用Broad Institute的IGV(Integrated Genomics Viewer)對數(shù)據(jù)進行可視化, 圖4同。

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

1.2.5 變異注釋軟件ANNOVAR

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

2 自動化實現(xiàn)

2.1 基于Perl語言的流程設(shè)計

本數(shù)據(jù)處理流程主要使用Perl編程語言實現(xiàn)對各個軟件的高效串接和自動化操作[8]。一項計算任務(wù)正在進行時, Perl對它進行監(jiān)控; 當計算完成, Perl去查看它的輸出的計算結(jié)果, 并把結(jié)果作為下一個計算任務(wù)的輸入, 往計算節(jié)點上投放新的計算任務(wù)。如此循環(huán), 直到流程運行完畢。

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

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

圖4 單堿基突變示例

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

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

流程的運行環(huán)境是計算機集群, 其有別于普通PC機, 一般是由一臺管理主機來協(xié)調(diào)許多計算主機來完成大型的計算任務(wù)。根據(jù)這樣的硬件特點來設(shè)計流程, 需要考慮以下兩個問題：(1)如何讓眾多計算機協(xié)同工作; (2)程序設(shè)計盡可能讓計算任務(wù)并行,充分利用計算資源, 縮短計算時間。

SGE(Sun Grid Engine)是使用最廣泛的分布式資源管理器(DRM)。SGE軟件為用戶提供了SGE系統(tǒng)提交要求計算的任務(wù)的方法, 動態(tài)分配工作負荷[9]。主節(jié)點接受用戶提交的計算任務(wù), 根據(jù)計算節(jié)點的負載情況, 動態(tài)決定把計算任務(wù)分配到哪個計算節(jié)點上進行, 使眾多計算機協(xié)同工作。

通過分析各個軟件的工作方式, 我們對中間多個步驟進行了并行設(shè)計, 配合 SGE, 對計算機資源進行高效調(diào)用, 從而大大縮短流程運行時間。流程的并行設(shè)計如圖6所示。

圖6 流程并行設(shè)計

2.3 基于Java的圖形界面設(shè)計

按照預定格式填寫配置文件后, 本流程即可在終端直接運行, 不過為了進一步改善用戶體驗, 使操作更加更加簡潔直觀, 我們還提供了一個基于Java開發(fā)的圖形界面：WGS_Pipeline_Runner。使用時, 用戶可以直接調(diào)用已有的配置文件, 并進行修改, 也可以直接在表單界面進行填寫, 完成后可以保存至本地。完成配置文件后, 點擊 Run即可自動化完成分析流程, 一步輸出分析報告。圖形界面如圖7所示。

本流程所涉及的所有軟件說明、自動化代碼及使用說明、配置文件說明等均可在 https://github.com/ wksofia/wgs_pipeline中下載使用。

3 流程測試

3.1 運行效率統(tǒng)計

一個135 GB的人類全基因組測序數(shù)據(jù), 在計算機集群上使用該流程來處理大約耗時 50 h, 各階段運行耗時如表 1所示。與該自動化流程相比, 在不考慮中間銜接耗時、不采用并行的情況下, 執(zhí)行同樣流程用時在一周以上?？梢? 本流程不僅簡化了分析操作, 更極大地節(jié)約了時間, 從而加速科研進展。

3.2 結(jié)果展示

自動化運行完全部流程, 得到一系列結(jié)果, 包括：BWA align讀句定位生成的sai文件, BWA sampe整合pair-end信息得到的sam文件, Samtools convert轉(zhuǎn)換sam得到的bam文件, Samtools sort對bam文件排序得到的sorted.bam文件, Picard rmdup去除重復得到的sample_duprmed.bam文件, GATK UG和GATK VQSR得到的一系列raw.vcf文件, Filter過濾后得到的 filtered.vcf文件, 以及 Annotation注釋后的 csv變異文件。此外還給出了一個包含對實驗數(shù)據(jù)質(zhì)量評價的summary文件。綜合以上結(jié)果, 用戶能夠從中挖掘出感興趣的變異信息。

在這些結(jié)果中, 用戶最值得關(guān)注的主要有兩個文件：(1)經(jīng)過功能注釋的變異列表(見Annotation文件夾); (2)對實驗數(shù)據(jù)質(zhì)量的評價表(見sample.Summary)。

變異列表(部分)如表2所示。

每個個體大約會攜帶大約3百萬個所謂的“變異”, 其中一些跟某些疾病的患病風險有關(guān), 科研人員正是希望找到這種致病變異。表中每一行代表一個變異, 這個列表包含的信息主要有：這個變異在所在的基因, 變異的功能, 是否處于重復序列, 是否被前人報道過。

測序數(shù)據(jù)評價如表3所示。該表主要關(guān)注兩個方面：平均測序深度和參考序列的覆蓋度。

圖7 圖形界面

表1 流程詳細運行時間

該評價表會在分析完成后, 通過電子郵件自動發(fā)送到用戶郵箱, 既便于用戶第一時間知道自己的數(shù)據(jù)質(zhì)量, 也方便了監(jiān)控。

表2 變異列表(部分)

4 結(jié) 語

本項目成功整合了一系列二代測序數(shù)據(jù)分析軟件, 形成了一套經(jīng)典的數(shù)據(jù)分析流程。本流程通過并行化設(shè)計和自動化處理, 一方面簡化了操作成本、縮短了數(shù)據(jù)分析周期, 另一方面也使本流程可以引入更完善的數(shù)據(jù)校驗步驟, 增強結(jié)果的可信度。本流程針對Illumina平臺雙端測序數(shù)據(jù)開發(fā), 滿足了大部分處理需求, 并對其他用戶提供了一個很好的參考, 后續(xù)我們將根據(jù)用戶需求對該自動化流程進行持續(xù)維護。

表3 測序數(shù)據(jù)評價表(部分)

隨著二代測序技術(shù)的逐步發(fā)展, 二代測序已經(jīng)廣泛應(yīng)用于科研和臨床研究。本流程提高了二代測序數(shù)據(jù)分析的入門和運轉(zhuǎn)效率, 其必將在二代測序相關(guān)基因組學研究中, 促進廣大科研人員工作的高效進行。

[1] Illumina Inc. Illumina Sequencing Technology. http://www. illumina.com/documents/products/techspotlights/techspotl ight_sequencing.pdf.

[2] Cock PJA, Fields CJ, Goto N, Heuer ML, Rice PM. The Sanger FASTQ file format for sequences with quality scores, and the Solexa/Illumina FASTQ variants. Nucleic Acids Res, 2010, 38(6): 1767–1771.

[3] Li H, Durbin R. Fast and accurate long-read alignment with Burrows-Wheeler Transform. Bioinformatics, 2010, 26(5): 589–595.

[4] Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, 1000 Genome Project Data Processing Subgroup. The sequence alignment/map (SAM) format and SAMtools. Bioinformatics, 2009, 25(16): 2078–2079.

[5] Picard. http://picard.sourceforge.net

[6] McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M, DePristo MA. The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res, 2010, 20(9): 1297–1303.

[7] Wang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res, 2010, 38(16): e164.

[8] Scbwartz RL, Pboenix T, Foy BD著. 盛春, 蔣永清, 王暉譯. Perl語言入門 (第五版). 南京: 東南大學出版社, 2009, 200.

[9] ORACLE INC. N1 Grid Engine 6 用 戶 指 南 . http://docs.oracle.com/cd/E19080-01/n1.grid.eng6/817-7681/esqcr/index.html.

(責任編委：胡松年)

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

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 variants that researchers can further analyze. The pipeline simplifies the manual operation and improves the efficiency byautomatization 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.

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

2013-09-07;

2014-01-20

國家重點基礎(chǔ)研究發(fā)展計劃(973計劃)項目(編號：2010CB529505)和中央高?；究蒲袠I(yè)務(wù)費專項資金(編號：2012-XHGX02)資助

李文軻, 碩士, 助理研究員, 研究方向：生物信息學。Tel: 010-88396071; E-mail: wksofia@gmail.com

趙倩, 博士, 副研究員, 研究方向：遺傳學, 生物信息學。E-mail: zhaoqian82@gmail.com

10.3724/SP.J.1005.2014.0618

時間: 2014-3-25 17:22:30

URL: http://www.cnki.net/kcms/detail/11.1913.R.20140325.1722.002.html