基因组二代测序数据的自动化分析流程 南京廖华

  本文关键词：基因组二代测序数据的自动化分析流程，由笔耕文化传播整理发布。

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

技术与方法

基因组二代测序数据的自动化分析流程

李文轲1, 李丰余1,2, 张思瑶1, 蔡斌1, 郑娜1, 聂宇1, 周到2, 赵倩1

1. 中国医学科学院, 北京协和医学院, 国家心血管病中心, 阜外心血管病医院, 心血管疾病国家重点实验室, 北京 100037; 2. 中南民族大学生物医学工程学院, 武汉430074

摘要： 二代测序技术的发展对测序数据的处理分析提出了很高的要求。目前二代测序数据分析软件很多, 但是

绝大多数软件仅能完成单一的分析功能(例如：仅进行序列比对或变异读取或功能注释等), 如何能正确高效地选择整合这些软件已成为迫切需求。文章设计了一套基于perl语言和SGE资源管理的自动化处理流程来分析Illumina平台基因组测序数据。该流程以测序原始序列数据作为输入, 调用业界标准的数据处理软件(如：BWA, Samtools, GATK, ANNOVAR等), 最终生成带有相应功能注释、便于研究者进一步分析的变异位点列表。该流程通过自动化并行脚本控制流程的高效运行, 一站式输出分析结果和报告, 简化了数据分析过程中的人工操作, 大大提高了运行效率。用户只需填写配置文件或使用图形界面输入即可完成全部操作。该工作为广大研究者分析二代测序数据提供了便利的途径。

关键词：　二代测序; 自动化数据分析; 流程; 变异检测

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

基金项目:国家重点基础研究发展计划(973计划)项目(编号：2010CB529505)和中央高校基本科研业务费专项资金(编号：2012-XHGX02)资助

作者简介:李文轲, 硕士, 助理研究员, 研究方向：生物信息学。Tel: 010-88396071; E-mail: wksofia@gmail.com

通讯作者:赵倩, 博士, 副研究员, 研究方向：遗传学, 生物信息学。E-mail: zhaoqian82@gmail.com DOI: 10.3724/SP.J.1005.2014.0618 网络出版时间: 2014-3-25 17:22:30

URL:

第6期

李文轲等: 基因组二代测序数据分析流程的自动化实现

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

二代测序技术(Next-generation sequencing)大幅度降低了测序的时间和成本, 使得大规模测序逐渐成为常规的实验室研究和临床检测手段。测序产生的数据量急剧增加, 如何高效地分析这些数据, 已成为迫切需要解决的问题。目前, 分析序列信息的生物信息学软件纷繁复杂, 但基本上每个软件只能完成单一的分析功能, 实现一个完整的分析流程则需要对众多软件进行整合, 而手动串联的效率往往不尽人意; 同时, 这些软件需要在Linux工作环境下以命令行运行, 要求用户具备较好的计算机背景; 另外, 即便一些实验室完成了分析流程的构建, 他们往往不会公开许多细节, 新用户仍然要从头建起。本研究致力于构建经典的二代测序数据分析流程, 并实现各个环节的高效自动化管理和分析, 减轻研究者前期的工作负担, 促进相关领域进一步对基因组测序研究项目的顺利开展。

(Hiseq2500), 甚至更高的250个碱基(Miseq)。测序读长不断增加, 测序通量也在不断上升。Illumina Hiseq2500是目前世界上通量最高的测序平台, 最多可以在大约10 d的时间内测定3000亿个碱基——即6～7个人类全基因组或60～80个人类全外显子组的序列测定。

Illumina平台以FASTQ格式[2]存储测序结果, 这也是本流程的输入文件。FASTQ文件记录内容包括所测的碱基读段和质量, 其数据格式如图1所示。每条读段(reads)占四行：第一行和第三行为读段识别码, 包含测序仪SN号、产生读段的巷道(lane)、该读段的编号等信息; 第二行为读段测到的碱基序列; 第四行为所测到碱基的质量分数, 每一个碱基都会对应一个质量分数。

1.2 数据处理流程及软件简介

目前测序数据处理软件很多, 我们综合考虑了适用性和效率, 整合出了一套标准的数据处理流程。具体来说, 获得FASTQ格式的原始测序数据后, 需要对数据进行以下处理：(1)使用BWA软件把这些短序列和参考基因组进行对比, 确定短序列在基因组上的位置, 把短序列组装成完整的人类参考基因组; (2)使用Samtools软件把这些短序列调整成按一定顺序(1~22, X, Y, 其他)排列的序列, 并进行数据格式的转换; (3)使用Picard软件把测序产生的冗

1 数据的获取和分析流程的构建

1.1 Illumina测序数据

本流程适用于Illumina测序平台产出的双端(Paired ends)测序数据。Illumina测序技术采用边合成边测序(Sequencing by synthesis, SBS)的方法, 早期的GA测序仪测序读长只有100个碱基, 随着技术的改进, 目前的读长已经增加到150个碱基

图1 FASTQ格式示例

620

Hereditas (Beijing) 2014

第36卷

余信息和噪声去掉; (4)使用GATK寻找样本测序数据与参考基因组的差异, 列出这些差异点; (5)使用Annovar对这些变异位点进行功能注释, 得到一个易于理解的变异位点列表。处理流程图如图 2 所示。

的序列比对软件, 能高效地比对短序列和参考基因组, 找到短序列在参考基因组上的位置, 该软件最长支持至1 Mb的短序列比对。BWT方法通过B-W转换将基因组序列按一定规则压缩并建立索引, 再通过查找和回溯来定位读段, 在查找时可通过碱基替代来实现允许的错配。采用Burrows-Wheeler转换的代表软件是Bowtie和BWA。比对结果如图3所示：界面上方是测到的短序列, 下方是短序列所比对到的参考基因组。

1.2.2 SAM文件处理软件Samtools

读段定位到基因组后推荐采用SAM(Sequence Alignment/Map)格式或其二进制版本BAM格式来存储。二进制版本可大大节省存储空间, 但不能直接用普通文本编辑工具显示。

SAM文件处理软件Samtools可以很好的对SAM/BAM格式数据进行操作, 因此, 本文使用它来进行数据格式转换和排序。

1.2.3 测序噪声去除和测序数据评价软件Picard

对组装好的全基因组数据, 需要将过度重复测

到的数据进行剔除, 并且需要对数据质量进行评价, Picard软件可以很好地完成这两项工作。 1.2.4 变异检测软件GATK

GATK主要用于在测序数据中寻找变异, 包括单碱基变异(SNV)、短插入缺失(INDEL), 是当前业界用来寻找变异的主流软件。

图2 处理流程图

1.2.1 读句比对软件BWA

BWA(Burrows-Wheeler Alignment tool)是基于Burrows-Wheeler变换(Burrows-Wheeler Transform)

图3 比对结果示例

利用Broad Institute的IGV(Integrated Genomics Viewer)对数据进行可视化, 图4同。

第6期

李文轲等: 基因组二代测序数据分析流程的自动化实现

621

变异指测序序列和参考序列的差异。如图4所示, 参考序列上的碱基是胸腺嘧啶(T), 而测序数据上的碱基是鸟嘌呤(G), 说明此处有一个T→G 的 突变。

1.2.5 变异注释软件ANNOVAR

ANNOVAR是一个用于高效注释变异的工具。注释信息包括变异所在的染色体, 开始位置, 结束位置, 参考序列信息和观察到的序列信息的列表。一个变异经过ANNOVAR注释之后, 其功能一目了然, 便于进一步的生物学分析。

务正在进行时, Perl对它进行监控; 当计算完成, Perl去查看它的输出的计算结果, 并把结果作为下一个计算任务的输入, 往计算节点上投放新的计算任务。如此循环, 直到流程运行完毕。

同时, 由于每次运行的样本不同, 数据的输入输出位置也有差异。如果每处理一个新的样本, 就要对流程源码进行大量修改, 将不利于流程的使用。为此, 本流程定义了一个配置文件(config file)。通过配置文件可以指定：流程处理的样品名、数据输入输出路径、参考序列文件, 甚至流程中涉及到的软件的位置、软件的运行方式; 另外, 我们还提供了对流程中主要软件参数的修改, 以满足高级用户需求。每次进行一个新样本的分析, 不需要修改主程序代码, 只要为其创建一个配置文件, 主程序会自动读取配置文件, 生成相应的执行代码。

流程文件构成如图5所示。

2 自动化实现

2.1 基于Perl语言的流程设计

本数据处理流程主要使用Perl编程语言实现对各个软件的高效串接和自动化操作。一项计算任

图4 单碱基突变示例

2.2 基于资源管理软件(SGE)的并行设计

流程的运行环境是计算机集群, 其有别于普通PC机, 一般是由一台管理主机来协调许多计算主机

来完成大型的计算任务。根据这样的硬件特点来设计流程, 需要考虑以下两个问题：(1)如何让众多计

图5 分析流程结构

  本文关键词：基因组二代测序数据的自动化分析流程，由笔耕文化传播整理发布。