深度解析：基因组De Novo组装技术原理与实战应用

时间：11-08 名人轶事提交错误

老铁们，大家好，相信还有很多朋友对于深度解析：基因组De Novo组装技术原理与实战应用和的相关问题不太懂，没关系，今天就由我来为大家分享分享深度解析：基因组De Novo组装技术原理与实战应用以及的问题，文章篇幅可能偏长，希望可以帮助到大家，下面一起来看看吧！

1. 基因组组装

2. 基于De-Bruijn Graph的组装算法

3. SOAPdenovo安装及使用说明：安装、描述、配置、操作

4. SOAPdenovo案例实践：数据下载、配置、操作、输出

写在前面

大家好，这是我们六点介绍生物信息学大数据分析和基因组数据分析系列文章的第一篇。我们将持续与大家分享生物医学大数据处理相关的知识和案例，希望能够帮助大家在自己的项目中更好的进行生物医学健康大数据处理。本文主要由四个部分向您介绍基因组从头的知识以及详细的应用案例。基因组组装，基于De-Bruijn Graph的组装算法，SOAPdenovo安装及使用说明：安装、说明、配置、操作，SOAPdenovo案例实践：数据下载、配置、操作、输出。

1. 基因组组装

基因组组装是生物信息学领域的核心问题。如果你想深入研究一个生物体，获得参考基因组是第一步也是必要的一步。基因组组装是将原始序列恢复为DNA序列片段，甚至整个物种的整个基因组序列的过程。

基因组组装是基因组分析的关键，对于物种起源和进化、挖掘功能基因、研究疾病的发生和发展具有重要意义。

然而，由于目前市场上广泛使用的二代测序技术获得的测序序列普遍较短，如何从短片段组装出完整的基因组成为迫切需要解决的问题。

基因组组装可分为基于参考基因组的组装（Mapping assembly）和从头组装（denovo assembly）。两者的主要区别在于是否有已知的基因组参考序列作为参考。在这篇文章中，我们主要介绍从头组装，这是一种不依赖于任何基因组参考序列相关信息的序列组装。目前，主流基因组从头组装使用的算法主要有两种：OLC方法（Overlap-Layout-Consensus）和DBG方法（De-Bruijn Graph）[2]。

DBG方法的核心思想是将序列拼接问题转化为众所周知的欧拉图问题[3]。

DBG方法内存消耗相对较低，运算速度较快，准确度较高。

目前主流的基因组装算法都是基于DBG方法进行设计和改进的。

2. 基于De-Bruijn Graph的组装算法

前面我们讨论了基因组从头组装的两种方法。接下来我们主要讨论基于De-Bruijn Graph的组装算法的基本原理。这里以目前广泛使用并由BGI团队开发的SOAPdenovo[4]为例。软件参考[5] 有兴趣的话可以看一下参考资料，阅读一下。

答：基因组DNA 被打碎成小片段，用于文库构建和双端测序。对于150~500bp的片段，直接进行双端测序，对于2-10kb的长片段，先进行环化，然后进行双端测序。

B：组装的核心部分，构建De-Bruijn Graph。构建De-Bruijn图的第一步是对测序reads进行k-merize，所谓k-mer是指将reads分成包含k个碱基的字符串，即在整个序列上取一个k长度的窗口read 1 通过逐个碱基滑动，每个滑动窗口内会生成一个k大小的序列，即一个k聚体。因此，一般长度为m的reads可以分为m-k+1 k-mers。其中，k必须为奇数。如果是偶数，则遇到回文序列时可能会产生完全相同的k-mers。我们使用k-mers 作为图的节点。如果两个节点有K-1个公共重叠子集，则将这两个节点连接在一起。这将形成一个De-Bruijn 图。您可以看到该图表非常好。显示序列的顺序信息。

C：简化图结构。虽然前面的步骤已经初步构建了图，但实际上由于测序错误、重复、杂交等原因，图上会出现很多提示、气泡等问题，因此需要进一步简化。这里的简化主要包括四个方面：1）去掉tips（可能是测序错误造成的）； 2) 移除低覆盖路径； 3）解开微小的重复区域（可以通过read pass来解决）4）合并气泡气泡区域（可能是由测序错误、重复或杂合性引起的）。

D: 分离出重叠群。切断重复的节点并输出重叠群。

E: 建造脚手架。读取和重叠群被重新比对，并且配对末端信息用于将单个重叠群连接到支架中。 1) 配对reads与contigs比对，以连接相邻的contigs； 3）使用不同的双端信息插入，由短到长逐步构建支架。

F:最终将多个支架组装成无GAP的基因组序列。

3.SOAPdenovo安装和使用说明

3.1 安装

SOAPdenovo 已更新为SOAPdenovo2，github[6] 链接：https://github.com/aquaskyline/SOAPdenovo2。

直接下载二进制文件[7] (https://sourceforge.net/projects/soapdenovo2/files/SOAPdenovo2/)

Mac 链接：linux 链接：

源码安装：

wgetunzip r242.zipcd SOAPdenovo2-r242make

安装完成后，可以看到两个可执行文件，SOAPdenovo-127mer和SOAPdenovo-63mer。 63mer表示支持的最大kmer长度为63，127mer表示支持的kmer最大长度为127。除了支持的kmer长度不同外，其他用法完全相同。

3.2 使用说明

由于计算量比较大，SOAPdenovo对计算机配置有一定的要求。官网对运行配置的解释是：SOAPdenovo适用于大型植物和动物基因组，尽管它也适用于细菌和真菌基因组。它运行在具有至少5G物理内存的64位Linux系统上。对于像人类这样的大型基因组，大约需要150 GB 的内存。运行SOAPdenovo-63mer，可以看到SOAPdenovo主要包含以下6个子命令：

1. pregraph 构造kmer-graph # 2.稀疏_pregraph 构造稀疏kmer-graph 3. contig 消除错误并输出contigs 4. map map 读取到contigs 5. scaff 构造脚手架6. 全部在tur 中执行pregraph-contig-map-scaff

其中1-5分别代表组装的4个步骤（1和2为两种组合方式，二选一），all用于一次性执行以上4个步骤。实际应用中，可以使用SOAPdenovo全部一步运行，也可以分4步单独运行。

3.3.软件运行前的准备工作——配置文件

soapdenovo 需要一个配置文件config_file，它提供输入文件和一些参数设置。

下面是配置文件的示例和说明：

max_rd_len=100#全局配置参数：如果序列大于这个长度，就会被切割成这个长度，然后分析[LIB]#每个库的配置以[LIB]开头avg_ins=200#平均长度文库插入，在实际设置时，可以参考文库大小分布图，取峰值。 verse_seq=0#序列是否需要反向互补。对于pair-end数据，不需要反向互补，设置为0；对于配对数据，这是必需的。反向互补，设置为1asm_flags=3#1表示仅组装contig。 2表示只组装scaffold，3表示同时组装contig和scaffold，4表示只互补gaprd_len_cutoff=100 #序列长度阈值，效果和max_rd_len一样，大于这个长度的序列会被切割到这个长度rank=1#设置不同库数据的优先顺序。取值范围为整数。组装scaffold时会同时使用多个具有相同rank值的库。 pair_num_cutoff=3# contig或scaffold之前的最小重叠数，对于pair-end数据，默认值为3；对于mate-paird数据，默认值为5map_len=32#比较长度的最小阈值，对于pair-end数据，默认值为32；对于配对数据，默认值为35q1=fastq1_read_1.fq# read1 序列文件q2=fastq1_read_2.fq# read2 序列文件

3.4.软件操作

分体式：

#step1:SOAPdenovo-63mer pregraph -s config_file -K 63 -R -o graph_prefix 1pregraph.log 2pregraph.err #ORSOAPdenovo-63mer稀疏_pregraph -s config_file -K 63 -z 5000000000 -R -o graph_prefix 1pregraph.log 2pregraph.err#step2:SOAP denovo - 63mer contig -g graph_prefix -R 1contig.log 2contig.err#step3:SOAPdenovo-63mer map -s config_file -g graph_prefix 1map.log 2map.err#step4:SOAPdenovo-63mer scaff -g graph_prefix -F 1scaff.log 2scaff.err

第一步：

SOAPdenovo-63mer all -s config_file -K 63 -R -o graph_prefix 1ass.log 2ass.err

输出文件：运行后会生成很多文件。后缀contig和scafSeq是对应的组装结果，分别对应contig和scaffold结果。

4. SOAPdenovo 实践

接下来我们找到NA12878样本的测序数据，详细实践一下。

4.1 数据下载

下载测序数据：

wget http://www.sixoclock.net/resources/data/NGS/Homo_sapiens/WGS/NA12878.WGS-100K_1.fastq.gzwget http://www.sixoclock.net/resources/data/NGS/Homo_sapiens/WGS/NA12878.WGS-100K_2.fastq。广州

4.2 准备配置文件

vi config_file，填写以下内容：

max_rd_len=100[LIB]avg_ins=200reverse_seq=0asm_flags=3rd_len_cutoff=100rank=1q1=NA12878.WGS-100K_1.fastq.gzq2=NA12878.WGS-100K_2.fastq.gz

4.3 运行命令

这里我们将程序的标准输出和标准错误分别重定向到相应的log和err文件中。一步操作：

SOAPdenovo-63mer all -s config_file -K 63 -R -o graph_prefix 1ass.log 2ass.err

单独运行的四个步骤：

SOAPdenovo-63mer pregraph -s config_file -K 63 -R -o graph_prefix 1pregraph.log 2pregraph.errSOAPdenovo-63mer contig -g graph_prefix -R 1contig.log 2contig.errSOAPdenovo2-r242/SOAPdenovo-63mer map -s config_file -g graph_prefix 1map。日志2map.errSOAPdenovo2-r242/SOAPdenovo-63mer scaff -g graph_prefix -F 1scaff.log 2scaff.err

4.4 输出结果

这里我们的测试数据已经被截取了，所以可以运行得很快。具体结果如下图所示。可以看到生成了很多中间结果文件。组装的重叠群和支架结果在图中圈出。出来的两个文件：

*.contig：contig序列文件，fasta格式；

*.scafSeq：脚手架序列文件，contig之间的间隙用N填充。

*.log和*.err是运行日志，里面包含很多统计信息，比如N50、N90、contig/Scaffold等信息。

*.scaf：包含scaffold中contig的详细信息；包括支架名称、重叠群长度和支架系中支架的长度。重叠群线包括重叠群名称、重叠群在支架上的起始位置、正向和反向链、长度以及重叠群之间的链接信息。

*.links：重叠群之间的对端连接信息

*.readOnContig：重叠群上的读取位置。

4.5 从Sixoclock 下载soapdenovo2

另外，六点老官网封装了基于CWL（通用工作流语言）的SOAPdenovo2软件，该软件可以通过我们开发的‘sixbox’软件快速运行。如果您不了解Sixbox，可以通过六点官方网站了解更多信息。具体运行步骤如下：

1）下载cwl源码

Sixbox拉取cad377c5-1a22-4a60-b761-d6e95e0d806b或从Sixbox官网下载soapdenovo2.cwl

2）下载数据

wgethttp://www.sixoclock.net/resources/data/NGS/Homo_sapiens/WGS/NA12878.WGS-100K_1.fastq.gzwgethttp://www.sixoclock.net/resources/data/NGS/Homo_sapiens/WGS/NA12878.WGS-100K_2.fastq.gz

3）使用sixbox生成参数模板文件（YAML）并配置yaml文件

Sixboxrun --make-template ./soapdenovo2.cwlsoapdenovo2.job.yamlvimsoapdenovo2.job.yaml#编辑参数配置文件，替换或设置参数，实现个性化分析

如果您不熟悉，可以直接将下面的示例内容粘贴到soapdenovo2.job.yaml中

q1: # "File" 类型数组（可选）-class: Filepath: NA12878.WGS-100K_1.fastq.gzq2: # "File" 类型数组-class: Filepath: NA12878.WGS-100K_2.fastq.gzo: graph_prefix #defaultvalue oftype "string".K:63# type "int" 的默认值。 (可选)R:true#"boolean"类型的默认值.rank:-1#"int"类型的数组(可选)rd_len_cutoff:-100#"int"类型的数组(可选)asm_flags:-3#"int"类型的数组(可选)reverse_seq:-0 大批oftype"int"（可选）avg_ins:-200# 类型"int"的数组（可选）max_rd_len:200#类型"int"的默认值。（可选）all:true#type"boolean"

4）使用sixbox运行

Sixboxrun ./soapdenovo2.cwl ./soapdenovo2.job.yaml#或者指定输出目录Sixboxrun --outdir /home/path ./soapdenovo2.cwl ./soapdenovo2.job.yaml

运行后，可以在当前目录或指定的输出目录中看到相应的SOAPdenovo组装结果文件输出。

至此，SOAPdenovo的实践体验就基本结束了。

以上就是我们为大家带来的基因组denovo的基本原理，以及在平台上运行经典的SOAPdenovo的详细操作流程。也欢迎您到我们六点官网看看我们发布的SOAPdenovo2的CWL处理工具。

如果您对生物医学健康大数据相关内容感兴趣，可以继续关注我们。如果想探索更多软件流程或者知识文档，可以查看六点官网。

参考

[1]六点钟官方网站：http://www.sixoclock.net

[2] OLC方法（Overlap-Layout-Consensus）和DBG方法（De-Bruijn Graph）:https://zh.wikipedia.org/wiki/%E5%BA%8F%E5%88%97%E7%B5%84%E8% A3%9D

[3]欧拉图问题：https://baike.baidu.com/item/Euler Graph/2587300

[4]SOAPdenovo:https://github.com/aquaskyline/SOAPdenovo2

[5] 参考文献：http://www.genome.org/cgi/doi/10.1101/gr.097261.109