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.genomesize.com/

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.com/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.genomics.org.cn/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.genomics.org.cn/soapdenovo.html#intro2

 

基因組組裝結果評估：

 

*****

 

 

 

 

 

參考：

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

http://www.life.umd.edu/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