SOAPdenovo组装软件使用记录

背景：

 

 

1、为什么要从头测序组装基因组？

基因组是不同表型的遗传基础；获得参考基因组是深入研究一个生物体全基因组的第一步也是必须的一步；从头测序组装能够对新的测序物种构建参考基因组；

2、为什么要研究全基因组？

确定基因组中缺失了什么；确定难以生化研究的基因和pathways；研究感兴趣的pathway通路中的每一个基因；研究基因组的非编码区域（introns内含子、promoters启动子、telomeres端粒等）的调控机理和结构特征；基因组提供了一个可以进行各种统计的大型数据库（provide large databases that are amenable to statistical methods）;识别不同的可能有细微表型的序列；研究物种和基因组的进化过程。

 

测序前物种信息的准备：

 

 

　　搜集物种相关信息（染色体的倍型、基因组大小、杂合度、重复序列比例、是否有可用的遗传图谱、GC含量 和 GC分布）。提供已经发表的近源物种。根据近源物种分析以上信息，尤其是GC含量以及对应的GC分布，重复程度。

1、获取基因组大小

作用：

1）、  基因组太大（>10Gb）超出了目前denovo组装基因组软件对机器内存的要求，则无法实现组装。

2）、  对组装结果的大小的进行正确性与否判断。

途径：

1）、  动物基因组大小数据库（ANIMAL GENOME SIZE DATABASE）: http://www./

2）、  对于查不到的物种基因组大小的，可以通过一些方法去估计。

一个通过实验（流式细胞仪）估计基因组大小的例子：Yoshida, S., J. K. Ishida, et al. (2010). "A full-length enriched cDNA library and expressed sequence tag analysis of the parasitic weed, Striga hermonthica." BMC Plant Biol 10: 55.

基于福尔摩根染色估计基因组大小：

一本经典书，Gregory, T. (2005). The evolution of the genome《基因组进化》, Academic Press.

通过定量PCR估计基因组大小的例子：Wilhelm, J., A. Pingoud, et al. (2003). "Real-time PCR-based method for the estimation of genome sizes." Nucleic Acids Res 31(10): e56.

Jeyaprakash, A. and M. A. Hoy (2009). "The nuclear genome of the phytoseiid Metaseiulus occidentalis (Acari: Phytoseiidae) is among the smallest known in arthropods." Exp Appl Acarol 47(4): 263-273.

一个通过Kmer估计基因组大小的例子：Kim, E. B., X. Fang, et al. (2011). "Genome sequencing reveals insights into physiology and longevity of the naked mole rat." Nature 479(7372): 223-227.

2、杂合度评估

影响：

1）、  主要体现在不能合并姊妹染色体，杂合度高的区域，会把两条姊妹染色单体都组装出来，从而造成组装的基因组偏大于实际的基因组大小。

2）、  杂合度>0.5%,则组装有一定难度。杂合度>1%,则很难组装出来。

3）、  杂合度高，则组装出来的序列不适合用于后续生物学分析（eg:拷贝数、基因完整结构等）。

途径：

1）、  通过kmer分析的例子：http://www./nature/journal/vaop/ncurrent/full/nature11413.html

2）、  等

降低杂合度：可通过很多代近交来实现。

3、是否有遗传图谱可用

作用：随着测序对质量要求越来越高和相关技术的逐渐成熟，遗传图谱也快成了denovo基因组的必须组成。

遗传图谱构建相关概念，推荐参考书：The handbook of plant genome mapping: genetic and physical mapping 

4、生物学问题的调研

这一步也是很重要的。

 

测序：

 

 

1、测序——技术发展史：

 

二代测序NGS:next generation sequencing  or now generation sequencing

注：SOAPdenovo最初是为illumina测序平台设计的。

2、测序——策略选择：

　　一般都是用不同梯度的插入片段来测序，小片段（200,500,800）和大片段（1k, 2kb 5kb 10kb 20kb 40kb）。如果是杂合度高和重复序列较多的物种，可能要采取fosmid-by-fosmid或者fosmid pooling的策略。不言而喻，后者花费是相当高的。

 

3、基因组De Novo测序：

 

4、基因组重测序：

 

 

测序后的基因组组装原理：

 

 

 1、什么是基因组组装？

即测序序列组装，指通过aligning对齐和merging合并片段为一个更长的DNA序列，来重构建原始序列。

 

 2、测序和组装的两种策略：

 BAC-by-BAC:测序和组装每一个BAC,然后，合并BAC和移除BAC冗余部分，从而获得参考基因组序列。

 whole genome shotgun：全基因组鸟枪法，染色体DNA被随机打断成片段，然后依次测序和组装。

 

评估：

BAC-by-BAC：复杂，耗时长，劳动密集型，低复杂度计算 ，高成本，高质量，使用少。

whole genome shotgun：容易，实验步骤快速，计算步骤困难，性价比高，广泛应用

 3、二代测序数据从头组装的解决overlap的三种算法：

overlap-layout-consensus：重叠布局一致OLC法，【软件：PHRAP.NEWBLER.CABOG.CELERA.SHORTY.EDENA,popular for long reads】，1. Overlap discovery involves all-against-all, pair-wise read comparison. 2. Construction an approximate read layout according to the pair-wise alignment 3. Multiple sequence alignment determines the precise layout and the consensus.

De bruijn graph:DBG法，【软件：SOAPdenovo2.Velvet.EULER,popular for illumina ，for short reads】,1.所有的测序reads都被切割成某一固定Kmer长度的序列（21bp=<kmer<=127bp）.2.相邻kmers链接是来自read序列，所以它不需要成对序列比对(The links between neighboring Kmers are derived from read sequences,so it doesn’t need pair-wise reads alignment.）3.冗余的数据自动被压缩。

greedy method：贪婪法(use OLC or DBG)，【软件：SSAKE.SHARCGS.VCAKE】，从给定的reads和contigs开始，使用下一个得分最高的overlap去做下一个连接，这样一直做下去，直到不能进行下去为止。 

评估：组装short reads的挑战是，1.基因组的复杂性，重复序列、杂合的二倍体基因组heterozygous diploid genome、多倍性polyploidy.2.illumina reads 的数据特征，测序错误率~1%、short read 长度~100bp、~100X的高测序深度、不同级别的文库插入片段（200bp~40Kbp）。3.Complexity of computation

 

SOAPdenovo组装软件介绍：

 

 

官网：http://soap./soapdenovo.html#intro2

可下载地址：https://github.com/aquaskyline/SOAPdenovo2

论文：

《SOAPdenovo2:an empirically improved memory-efficient short-read de novo assembler》https://wenku.baidu.com/view/6fa2546069eae009581becd3.html?re=view###

Ruiqiang Li, et al. De novo assembly of human genomes with massively parallel short read sequencing. 2009, Genome Research.

1、说明：他是一种新型的short read组装软件，设计服务于大型的植物和动物基因组，尽管他对细菌和真菌的基因组也有效。利用De bruijn graph组装算法。是第一个利用short read的组装软件去组装哺乳动物基因组。他已经组装了数百种动植物的基因组，发表的文章有很多。

 

2、流程：

 

contiging:

a.  基因组DNA被随机打断，并且使用paired-end测序。长度在150-500bp的short clones扩增直接测序。然而，在2-10kb的长paired-end libraries 通过DNA环化、fragmentation破碎，然后净化400-600bp的碎片为了cluster 结构。

b.  raw reads 或者预修正reads被装入计算机内存中，并且，de Bruijn graph data structure 被用于表示reads间的overlap。

Kmer-graph构建：所有的测序reads都被切割成某一固定Kmer长度的序列（21bp=<kmer<=127bp），形成等长的Kmers。将Kmers连成图。相邻的kmers是通过K-1  overlaping连接的，所以它不需要成对序列比对（The neighboring kmers are K-1 overlaping which generated from read sequences, so it doesn’t need pair-wise reads alignment.）。重复序列在图中被压缩。

c.  会产生tips翼尖、bubbles气泡、low coverage links低覆盖率链接、tiny repeat微小重复等问题。

tips翼尖（a图）和bubbles气泡（c图）：由于测序错误或者杂合或者高重复序列相似性，将会导致翼尖和气泡出现。

low coverage links低覆盖率链接：（b图）（d图）。

tiny repeat微小重复(e图)：重复在graph中被压缩，并作为不同路径的共享边缘，但是能够通过reads 穿过他来解决。

 

 

 

                 e

 

移除错误链接和graph simplification图形简化，得到contigs or contig  graphs：tips翼尖移除；删除低覆盖链接；bubbles合并气泡；解决微小重复；

 

 

 

 

 

 d.  contig  graphs,在重复的节点处剪断，输出contigs

 

 scaffolding:

e.  重新用reads和contigs进行比对，使用paired-end信息来把单一的contigs连接成scaffolds。reads 比对到contigs上，临近的contig建立连接；repeat将会引来冲突矛盾信息；在组装成scaffold时，repeat contigs将会被屏蔽；paired-end信息的不同插入片段被用来一步步从短到长的建立scaffold graph(Scaffolding iteratively from short to long insert PEs./Various insert size of paired-end information is used to build contig graph step by step from short to long)。

 

 

 Gap Filling:

f.  使用paired-end reads来填补scaffolds内部可能是由重复序列所造成的Gap。contig N50 通常比较小（<3KB）,但是，gap filling之后能够显著提高N50值（i.e.,>20KB）;Most of the gaps are repeat relative sequences.;Reads locate at gaps can collected by their paired-end which uniquely map to the contig.

 

3.软件使用：

a.   SOAPdenovo可以一步跑完，也可以分成四步单独跑

一步跑完的脚本:

./SOAPdenovo all -s config_file -K 25 -R -D 1 -d  -o graph_prefix 1>ass-K25.log 2>ass-K25.err

四步单独跑的脚本:

./SOAPdenovo pregraph -s config_file  -K 25 -R -d 1 -p -o graph_prefix  >pregraph.log

./SOAPdenovo contig   -g graph_prefix -R -D 1 -M 1  >contig.log

./SOAPdenovo map      -s config_file  -g graph_prefix -p  -f >map.log

./SOAPdenovo scaff    -g graph_prefix -F -u -G -p >scaff.log

 

b.  参数说明

all:

     -s  <string>      solexa reads 的配置文件

    -p  <int>         程序运行时设定的cpu线程数，默认值[8]

pregraph: 

    -o  <string>      图形输出的文件名前缀 

    -K  <int>         输入的K-mer值大小，默认值[23]，取值范围 13-127 #K-mer值必须是奇数；组装杂合子基因组的K-mer值应该小一点；组装含有高repeats基因组且要求其有高的测序深度和长的reads,的K-mer应该大一点。

    -a  <int>         假设初始化内存， 为了避免进一步的再分配，默认[0]G

    -d  <int>         去除kmers频数不大于该值（kmerFreqCutoff）的k-mer，默认值[0] ##最小化错误测序带来的影响

    -R  (optional)    利用read鉴别短的重复序列，默认值[NO]

contigs：

      -g  <string>      输入graph file文件名前缀 

    -R  (optional)    移除repeats，使用pregraph步骤中产生的结果，如果参数-R在pregraph步骤中被设置的话，默认[NO]

    -D  <int>         去除频数不大于该值（edgeCovCutoff）的由k-mer连接的边，默认值[1]，即该边上每个点的频数都小于等于1时才去除。edges with coverage no larger than EdgeCovCutoff will be deleted, [1] #最小化错误测序带来的影响

    -M  <int>         在contiging操作时，合并相似序列的强度，默认值为[1]，最小值0，最大值3。#deal with heterozygosis

    -e  <int>         两边缘之间的弧的权重大于该值（arcWeight），将被线性化，默认值[0]

    -m  <int>         最大的kmer size（max 127）用于multi-kmer，默认[NO]

    -E  (optional)    在iterate迭代之前，合并clean bubble功能，仅在当使用multi-kmer时且设置-M参数，默认[NO]

map：

    -g  <string>      输入graph file文件名前缀 

    -k  <int>         该值（kmer_R2C）是kmer size用于mapping reads到contigs上时的值，默认值[K]

    -f  (optional)    在map那一步中，对于使用SRkgf去填充gap，输出与gap相关的reads,默认[NO]

scaffold：

      -F  (optional)    利用read对scaffold中的gap进行填补，默认[NO]

    -u  (optional)    构建scaffolding前不屏蔽高/低覆盖度的contigs，这里高频率覆盖度指平均contig覆盖深度的2倍。默认[mask]屏蔽

    -w  (optional)    在scaffold中，保持contigs弱连接于其他contigs，默认[NO]

    -G  <int>         估计gap的大小和实际补gap的大小的差异值，默认值[50]bp。

    -L  <int>         用于构建scaffold的contig的最短长度(minContigLen)，默认为：[Kmer参数值+2]

    -c  <float>       在scaffolding之前，contigs小于100bp，且覆盖率小于该最小contig覆盖率（c*avgCvg），将被屏蔽，除非参数-u已设定，默认值[0.1]

    -C  <float>       在scaffolding之前，contigs覆盖率大于该最大contigs覆盖率（C*avgCvg），或者contigs小于100bp且覆盖率大于0.8*（C*avgCvg），将被屏蔽，除非参数-u已设定，默认值[2]

    -b  <float>       当处理contigs间的pair-end连接时，如果参数-b>1，该插入片段的上限值（b*avg_ins）将被用来作为大的插入片段（>1000）的上限，默认值[1.5]

    -B  <float>       如果两个contigs的覆盖率都小于bubble覆盖率（bubbleCoverage）乘以contigs平均覆盖率（bubbleCoverage*avgCvg），则去除在bubble结构中的低覆盖率contig，默认值[0.6]

    -N  <int>         统计基因组大小，默认值[0]

    -V  (optional)    输出相关信息为了可视化组装，默认[NO]

    

c.   solexa reads配置文件config_file(需要手动配置)：

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

[LIB]                 #文库信息以此开头 #在其后，可以整多个文库，仍以[LIB]开头

avg_ins=300           #文库平均插入长度，一般取插入片段分布图中给出的文库大小，illumina测序数据平均插入片段长度为300bp

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

asm_flags=3           #该文库中的read序列在组装的哪些过程（contig/scaff/fill）中用到：              

                              设为1：只用于构建contig；                #454  or sanger or merging reads

 

                              设为2：只用于构建scaffold；              #长插入片段>=2k  , mate-pair libraries

 

                              设为3：同时用于构建contig和scaffold；    #短插入片段<2k 

 

                              设为4：只用于补洞gap filling             #454  single 长 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=./pathfastq_read_1.fq

#read 1的fastq格式的序列文件

q2=./pathfastq_read_2.fq

#read 2的fastq格式的序列文件(与read1对应的read2文件紧接在read1之后）

f1=./pathfasta_read_1.fa

#read 1的fasta格式的序列文件

f2=./pathfasta_read_2.fa

#read 2的fasta格式的序列文件

q=./pathfastq_read_single.fq

#单向测序得到的fastq格式的序列文件

f=./pathfasta_read_single.fa

#单向测序得到的fasta格式的序列文件

p=./pathpairs_in_one_file.fa

#双向测序得到的一个fasta格式的序列文件

#The assembler accepts read file in two formats: FASTA or FASTQ. Mate-pair relationship could be indicated in two ways: two sequence files with reads in the same order belonging to a pair, or two adjacent reads in a single file (FASTA only) belonging to a pair. In the configuration file single end files are indicated by “f=/path/filename” or “q=/pah/filename” for fasta or fastq formats separately. Paired reads in two fasta sequence files are indicated by “f1=” and “f2=”. While paired reads in two fastq sequences files are indicated by “q1=” and “q2=”. Paired reads in a single fasta sequence file is indicated by “p=” item. 

 

d.  输出文件

组装结果文件：

*.contig                #没有使用mate pair 信息的contig sequences 。

*.scafSeq              #SOAPdenovo软件最终的组装序列结果，可用于后续研究。

*.scafStatistics       #contigs和scaffolds的最终统计信息。

组装过程中产生的其他文件，详见官网：

*.kmerFreq             #每行显示一个数，这个数是kmer值出现的频率等于行号的kmer个数。

 

http://soap./soapdenovo.html#intro2

 

基因组组装结果评估：

 

*****

 

 

 

 

 

参考：

http://teacher.bmc./costuppsala2012/COSTUPPSALA2012/Lectures_files/SOAPdenovo-COST-XiaodongFANG-BGI.pdf

http://www.life./labs/delwiche/bsci348s/lec/Genomics.html

Jason R. Miller et al., Assembly algorithms for next‐generation sequencing data. Genomics

Li R, et al. De novo assembly of human genomes with massively parallel short read sequencing. Genome Research (2010).

http://blog.sina.com.cn/s/blog_5d1edf6a0100w56l.html

http://blog.sina.com.cn/s/blog_78c527410102w7ek.html

http://blog.sina.com.cn/s/blog_14ece68cc0102wagf.html