关于更多生物医疗大数据分析工具和软件的介绍和使用请看六点了官网[1]。
目录
1、基因组组装
2、基于De-Bruijn Graph的组装算法
3、SOAPdenovo的安装和使用说明：安装、说明、配置、运行
4、SOAPdenovo案例实战：数据下载、配置、运行、输出
写在前面
大家好，这是我们六点了给大家介绍生物信息大数据分析 基因组数据分析系列文章第一篇。我们会持续为大家分享关于生物医疗大数据处理相关的知识和案例，希望帮助大家更好地进行自己项目中生物医疗健康大数据处理工作。本篇文章主要四部分来为大家介绍基因组的denovo的知识和以及详细应用案例。①基因组组装、 ②基于De-Bruijn Graph的组装算法、 ③SOAPdenovo的安装和使用说明：安装、说明、配置、运行，以及 ④SOAPdenovo案例实战：数据下载、配置、运行、输出。
1. 基因组组装
基因组组装 (Genome assembly)是生物信息学领域的核心问题, 想要深入研究一个生物体，获得参考基因组是第一步也是必须的一步。基因组组装是将原始的下机序列还原成DNA序列片段、以至于整个物种全基因组序列的过程。
基因组组装是基因组分析的关键，对物种起源与进化，挖掘功能基因进而研究疾病发生和发展具有重大意义。
然而由于目前市面上广为应用的二代测序技术获得的测序序列一般都较短，因此如何通过短片段组装成完整的基因组成了亟待解决的问题。
基因组组装可分为基于参考基因组的组装（Mapping assembly）和从头组装（denovo assembly）。两者主要的区别在于是否存在已知的基因组参考序列作为参照。本文我们主要介绍的是denovo组装，即不依赖任何基因组参考序列相关信息而进行的序列组装。目前，应用于主流的基因组denovo组装的算法主要有两个[1]：OLC方法 (Overlap-Layout-Consensus)和 DBG方法 (De-Bruijn Graph)[2]。
而DBG方法的核心思想是将序列拼接问题转化为人们所熟知的欧拉图（Euler Graph）问题[3]。
DBG方法内存消耗相对较低，运算速度快，且准确率高。
目前主流的基因组装算法都是基于DBG方法改进设计的。
2.基于De-Bruijn Graph的组装算法
前面我们说到基因组denovo组装两种方法，下面主要展开说说基于De-Bruijn Graph的组装算法的基本原理。此处，就以目前使用比较广泛，由华大基因团队开发的SOAPdenovo[4]为例。软件的参考文献[5]有兴趣可以在参考资料看一下读读。
A：基因组DNA打断成小的片段，进行建库和双端测序。150~500bp的进行直接双端测序，长的片段2-10kb的则先进行环化再进行双端测序。
B：组装的核心部分，进行De-Bruijn Graph的构建。构建De-Bruijn图的第一步是将测序read k-mer化，而所谓的k-mer是指将reads分成包含k个碱基的字符串，即拿一个k长度的窗口在整个read上1个碱基一个碱基的滑动，每次滑动窗口内部都会产生一个k大小的序列,即为一个k-mer，因此一般长短为m的reads可以分成m-k+1个k-mers。其中k一定是奇数，如果是偶数遇到回文序列可能会产生完全相同的k-mers。我们将k-mers作为图的节点，如果两个节点有 K-1个共同重叠子集，就把两个节点连接在一起，这样就会形成De-Bruijn Graph，可以看到该图可以很好地展现出序列的顺序信息。
C：进行图结构的精简。尽管前面步骤已经初步构建出图形，但是实际上由于测序错误，重复，杂合等原因，图上会出现很多类似翼尖（tips）、气泡（bubbles）等问题，因此还需要进一步简化。此处简化主要包含四个方面：1）去除tips（可能为测序错误导致的）；2）去除低覆盖度的路径；3）解开微小重复的区域（可以通过read穿过来解决）4）合并bubbles气泡区（可能为测序错误，重复或者杂合导致的）。
D: 拆分出contig。在重复的节点处剪断，输出contigs。
E: 构建scaffolds。重新用reads和contigs进行比对，使用paired-end信息来把单一的contigs连接成scaffolds。1）paired reads 比对到contigs上，使临近的contig建立连接；3）paired-end信息的不同插入片段被用来一步步从短到长的建立scaffold.
F: 最终是把多个scaffold组装成无GAP的基因组序列。
3.SOAPdenovo的安装和使用说明
3.1安装
SOAPdenovo目前已更新到SOAPdenovo2,github[6]链接:https://github.com/aquaskyline/SOAPdenovo2。
直接下载二进制[7]（https://sourceforge.net/projects/soapdenovo2/files/SOAPdenovo2/）
mac链接: <https://sourceforge.net/projects/soapdenovo2/files/SOAPdenovo2/bin/r240/SOAPdenovo2-bin-r240-mac.tgz/download>linux链接: <https://sourceforge.net/projects/soapdenovo2/files/SOAPdenovo2/bin/r240/SOAPdenovo2-bin-LINUX-generic-r240.tgz/download>
源代码安装:
wget <https://github.com/aquaskyline/SOAPdenovo2/archive/refs/tags/r242.zip>unzip 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 construct kmer-graph # 2. sparse_pregraph construct sparse kmer-graph 3. contig eliminate errors and output contigs 4. map map reads to contigs 5. scaff construct scaffolds 6. all do pregraph-contig-map-scaff in tur
其中，1-5分别表示组装的4个步骤（1,2是两种构图方式，二选一），all则用于一次执行以上的4个步骤。实际应用中，可以使用SOAPdenovo all 一步式跑完，也可以分成4步单独去跑。
3.3. 软件运行前准备工作-配置文件
soapdenovo需要一个配置文件config_file,里面给定输入文件和一些参数设置。
下面是配置文件的示例和说明：
max_rd_len=100# 全局配置参数：如果序列大于该长度，会被切成该长度，然后在分析[LIB]#每个文库的配置以[LIB]开头avg_ins=200#文库插入片段的平均长度，在实际设置时，可以参考文库size分布图，取峰值即可reverse_seq=0#是否需要将序列反向互补，对于pair-end数据，不需要反向互补，设置为0；对于mate-pair数据，需要反向互补，设置为1asm_flags=3#1表示只组装contig. 2表示只组装scaffold,3表示同时组装contig和scaffold,4表示只补gaprd_len_cutoff=100#序列长度阈值，作用和max_rd_len相同，大于该长度的序列会被切除到该长度rank=1#设置不同文库数据的优先级顺序，取值范围为整数，rank值相同的多个文库，在组装scaffold时，会同时使用。pair_num_cutoff=3# contig或者scaffold之前的最小overlap个数，对于pair-end数据，默认值为3；对于mate-paird数据，默认值为5map_len=32#比对长度的最小阈值，对于pair-end数据，默认值为32；对于mate-pair数据，默认值为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 1>pregraph.log 2>pregraph.err#ORSOAPdenovo-63mer sparse_pregraph -s config_file -K 63 -z 5000000000 -R -o graph_prefix 1>pregraph.log 2>pregraph.err#step2:SOAPdenovo-63mer contig -g graph_prefix -R 1>contig.log 2>contig.err#step3:SOAPdenovo-63mer map -s config_file -g graph_prefix 1>map.log 2>map.err#step4:SOAPdenovo-63mer scaff -g graph_prefix -F 1>scaff.log 2>scaff.err
一步式：
SOAPdenovo-63mer all -s config_file -K 63 -R -o graph_prefix 1>ass.log 2>ass.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.gz
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 1>ass.log 2>ass.err
四步单独运行：
SOAPdenovo-63mer pregraph -s config_file -K 63 -R -o graph_prefix 1>pregraph.log 2>pregraph.errSOAPdenovo-63mer contig -g graph_prefix -R 1>contig.log 2>contig.errSOAPdenovo2-r242/SOAPdenovo-63mer map -s config_file -g graph_prefix 1>map.log 2>map.errSOAPdenovo2-r242/SOAPdenovo-63mer scaff -g graph_prefix -F 1>scaff.log 2>scaff.err
4.4 输出结果
此处我们的测试数据做了截取，因此可以非常快速的跑完，具体的结果如下图所示，可以看到生成了不少的中间结果文件，其中组装出来的contig和scafford结果即图上圈出来的两个文件：
*.contig：contig序列文件，fasta格式；
*.scafSeq：scaffold序列文件，contig之间的gap用N填充 。
*.log和*.err是运行的日志，里面包含很多的统计信息，如N50，N90，contig/Scaffold等信息。
*.scaf：包括scaffold中contig的详细信息；在scaffold行中包括scaffold名字、contig长度和该scaffold长度。在contig行包括contig名字、contig在scaffold上的起始位置、正反链、长度和contig间的链接信息
*.links：contig间的pair-end连接信息
*.readOnContig：reads在contig上的位置。
4.5 从sixoclock下载soapdenovo2
此外，六点了官网基于CWL (common workflow language) 对SOAPdenovo2软件进行了封装，通过我们开发的`sixbox` 软件可以快速进行软件的运行。对sixbox不了解可以通过六点了官网了解下。下面是具体的运行步骤如下：
1)下载cwl 源码
sixbox pull cad377c5-1a22-4a60-b761-d6e95e0d806b 或 在六点了官网上下载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.cwl > soapdenovo2.job.yamlvim soapdenovo2.job.yaml# 编辑参数配置文件，替换或设置参数以实现个性化分析
不熟悉的，可以直接粘贴下方示例内容到soapdenovo2.job.yaml
q1: # array oftype"File"(optional)-class: Filepath: NA12878.WGS-100K_1.fastq.gzq2: # array oftype"File"-class: Filepath: NA12878.WGS-100K_2.fastq.gzo: graph_prefix #defaultvalue oftype"string".K:63#defaultvalue oftype"int". (optional)R:true#defaultvalue oftype"boolean".rank:-1# array oftype"int"(optional)rd_len_cutoff:-100# array oftype"int"(optional)asm_flags:-3# array oftype"int"(optional)reverse_seq:-0# array oftype"int"(optional)avg_ins:-200# array oftype"int"(optional)max_rd_len:200#defaultvalue oftype"int". (optional)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流程工具。
如果对生物医疗健康大数据相关内容感兴趣也可以持续关注我们。想要探索更多的软件流程或者知识文档，可以到六点了官网查看。
References
[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]欧拉图（Euler Graph）问题:https://baike.baidu.com/item/欧拉图/2587300
[4]SOAPdenovo:https://github.com/aquaskyline/SOAPdenovo2
[5]参考文献:http://www.genome.org/cgi/doi/10.1101/gr.097261.109
[6]github:https://github.com/aquaskyline/SOAPdenovo2
[7]二进制:https://sourceforge.net/projects/soapdenovo2/files/SOAPdenovo2/
网友评论