组装组参考文档version1

附录四 基因组de novo组装

基因测序仍然是基因研究的重点。其应用包括鉴定一个新物种的基因组序列，鉴定一个种群中个体的基因组序列，测定一个特定的样品的RNA分子序列，还有把DNA序列作为一个解析出来的资料应用在分子生物学技术里。尽管人类、老鼠和很多其他的基因组已经被测定出来了，鉴定一个物种的全部基因组序列依然是测序的一个重要应用。它们只是生物界中百万种生物的一个很小的部分。

随着新一代高通量测序仪器的使用，测序的费用得到相当大的降低。但是，所测得到的数据读取的长度很短，这对 de novo 组装提出了非常严峻的挑战，下面将主要介绍SOAPdenovo组装软件的使用。

1 SOAPdenovo软件

背景

SOAPdenovo是利用一种新的组装短read的方法，它能为诸如人类基因组大小(3G)的基因组de novo组装。这个程序是为组装Illumina GA short reads特别设计的。它为构建全基因组参考序列序列和以低测序成本对未知基因组实施精确分析创造了可能。

程序的下载及安装：

下载地址：http://soap.genomics.org.cn/soapdenovo.html 安装：

a 下载SOAPdenovo的压缩包 b 解压缩

c 将得到可执行文件SOAPdenovo和一个配置文件的模板example.contig

使用程序及参数：

SOAPdenovo可以一步跑完，也可以分成四步单独跑 一步跑完的脚本:

./SOAPdenovo all -s lib.cfg -K 29 -D 1 -o ant >>ass.log

四步单独跑的脚本:

./SOAPdenovo pregraph -s lib.cfg -d 1 -K 29 -o ant >pregraph.log ./SOAPdenovo contig -g ant -D 1 -M 3 >contig.log ./SOAPdenovo map -s lib23.cfg -g ant >map.log ./SOAPdenovo scaff -g ant -F >scaff.log

2 参数说明

-s STR 配置文件

-o STR 输出文件的文件名前缀 -g STR 输入文件的文件名前缀

-K INT 输入的K-mer值大小，默认值23，取值范围 13-63 -p INT 程序运行时设定的线程数，默认值8

-R 利用read鉴别短的重复序列，默认值不进行此操作 -d INT 去除频数不大于该值的k-mer，默认值为0

-D INT 去除频数不大于该值的由k-mer连接的边，默认值为1，即该边上每个点的频数都小于等于1时才去除

-M INT 连接contig时合并相似序列的等级，默认值为1，最大值3。 -F 利用read对scaffold中的gap进行填补，默认不执行

-u 构建scaffold前不屏蔽高覆盖度的contig，这里高频率覆盖度指平均contig覆盖深度的2倍。默认屏蔽

-G INT 估计gap的大小和实际补gap的大小的差异，默认值为50bp。

-L 用于构建scaffold的contig的最短长度，默认为：Kmer参数值 ×2

3 使用方法及示例

1）示例

SOAPdenovo_Release1.0/SOAPdenovo all -s Data/HCB.lib -K 25 -d -o test 2) 输入文件

configFile (配置文件内容如下，非程序生成，需要软件使用者自己配置) #maximal read length （read的最大长度） 以“#”开头的行是注释内容 max_rd_len=50

#该值一般设置的比实际read读长稍微短一些，截去测序最后的部分，具体长度看测序质量 [LIB]

#文库信息以此开头 avg_ins=200

#文库平均插入长度，一般取插入片段分布图中给出的文库大小 reverse_seq=0

#序列是否需要被反转，目前的测序技术，插入片段大于等于2k的采用了环化，所以对于插入长度大于等于2k文库，序列需要反转，reverse_seq＝1，小片段设为0 asm_flags=3

#该文库中的read序列在组装的哪些过程（contig/scaff/fill）中用到 设为1：只用于构建contig； 设为2：只用于构建scaffold；

设为3：同时用于构建contig和scaffold； 设为4：只用于补洞

短插入片段(<2K)的设为3，同时用于构建contig和scaffold，长插入片段(>=2k)设为2，

不用于构建contig，只用于构建scaffold，454single 长reads只用于补洞。 rank=1

#rank该值取整数，决定了reads用于构建scaffold的次序，值越低，数据越优先用于构建scaffold。设置了同样rank的文库数据会同时用于组装scaffold。一般将短插入片段设为1；2k设为2；5k设为3；10k设为4；当某个档的数据量较大时，也可以将其分为多个档，同样，当某档数据量不足够时，可以将多个档的数据合在一起构建scaffold。这里说的数据量够与不够是从该档的测序覆盖度和物理覆盖度两个方面来考虑的。 pair_num_cutoff=3

#可选参数，pair_num_cutoff该参数规定了连接两个contig 或者是pre-scaffold 的可信连接的阈值，即，当连接数大于该值，连接才算有效。短插入片段(<2k)默认值为3，长插入长度序列默认值为5 map_len=32

#map_len该参数规定了在map过程中 reads和contig的比对长度必须达到该值（比对不容mismacth和gap），该比对才能作为一个可信的比对。可选参数，短插入片段(<2k)一般设置为32，长插入片段设置为35，默认值是K＋2。 q1=/path/**LIBNAMEA**/fastq_read_1.fq

#read 1的fastq格式的序列文件，“/path/**LIBNAMEA**/fastq_read_1.fq”为read的存储路径

q2=/path/**LIBNAMEA**/fastq_read_2.fq

#read 2的fastq格式的序列文件，与read1对应的read2文件紧接在read1之后） f1=/path/**LIBNAMEA**/fasta_read_1.fa #read 1的fasta格式的序列文件

f2=/path/**LIBNAMEA**/fasta_read_2.fa #read 2的fasta格式的序列文件

q=/path/**LIBNAMEA**/fastq_read_single.fq #单向测序得到的fastq格式的序列文件

f=/path/**LIBNAMEA**/fasta_read_single.fa #单向测序得到的fasta格式的序列文件

p=/path/**LIBNAMEA**/pairs_in_one_file.fa #双向测序得到的一个fasta格式的序列文件

4 输出文件及说明

SOAPdenovo分四部分别对应的输出文件：

pregraph生成7个文件 *.kmerFreq *.edge *.preArc *.markOnEdge *.path *.vertex *.preGraphBasic

contig生成4个文件 *.contig *.ContigIndex *.updated.edge *.Arc map生成3个文件 *.readOnContig *.peGrads *.readInGap

scaff生成6个文件 *.newContigIndex *.links *.scaf *.scaf_gap *.scafSeq *.gapSeq

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

*.scafSeq：fasta格式的scaffold序列文件，contig之间的gap用N填充

#对于得到的*.scafSeq文件还需要用GapCloser去合并其中的gap，最后的contig文件则是对补洞之后的scaffold文件通过打断N区的方法得到。 以上两个文件是组装结果中最主要的输出。 *.scaf：包括scaffold中contig的详细信息。

在scaffold行中包括scaffold名字、contig长度和该scaffold长度。

在contig行包括contig名字、contig在scaffold上的起始位置、正反链、长度和contig间的链接信息

*.links：contig间的pair-end连接信息 *.readOnContig：reads在contig上的位置

*.peGrads： 主要可以通过调整本文件中的参数来显示构建scaffold所用到的插入片段库的个数，总共要到的read数，最长的read的长度，每个库对应的哪些reads，rank设置，pair_num_cutoff设置。例如：

grads&num: 10 522083934 70 323 104577616 1 3 334 180770522 1 3 345 226070520 1 3 486 361955834 2 3 2200 392088076 3 5 2290 422272580 3 5 2400 445522690 3 5 4870 475666064 4 5 9000 511030930 5 8 9110 522083934 5 5

#文件中共分成4列。组装的配置文件中有n个文库，该文件则有n+1行，且按照文库大小顺序排列。

第1行中，第二三四列分别是 所用文库，reads总数和组装中用到的最长的reads长度。 第2行中，四列分别是文库大小，文库中的reads数目，该文库reads用到的rank等级和该文库中reads用到的pair_num_cutoff。 第3～n+1行，四列分别是文库大小，文库中的reads数目加上前面的文库中的reads总数，该文库reads用到的rank等级和该文库中reads用到的pair_num_cutoff。

如果配置文件中没有设置pair_num_cutoff，即使用默认参数，则最后一列显示为0。 对于SOAPdenovo的每个步骤都有日志文件输出，要保存好日志文件，日志文件中包含有很多有用的信息。

SOAPdenovo日志输出说明 pregraph.log:

其中有很多的统计信息，包括构建debruijn-graph时用到多少reads数，构图中生成了多少uniq的kmer以及设置-d参数后去除了多少kmer。 在pregraph中，可选参数有 –R –K –d

5467781332 nodes allocated, 70662750348 kmer in reads, 70662750348 kmer processed 3283081670 kmer removed

Kmer 数是取决于所设k值大小以及数据量，nodes数即特异性的kmer数目，当nodes数目过高（一般和基因组大小差不多大小），可能是数据的错误率比较高，也可能是存在杂合。若nodes数目偏小，并且kmer数目很多，则基因组本身可能存在一定的重复度。对于k值