http://www./article/26
這是一份關(guān)于基因組組裝軟件SOAPdenovo的使用說明，內(nèi)容包括了程序使用、參數(shù)的詳細(xì)說明、參數(shù)如何調(diào)整、各個(gè)主要輸出文件的格式說明等。

簡(jiǎn)介：
 
SOAPdenovo（目前最新版是SOAPdenovo2）是利用一種新的組裝短read的方法，它以kerm為節(jié)點(diǎn)單位，利用de Bruijn圖的方法實(shí)現(xiàn)全基因組的組裝，和其他短序列組裝軟件相比，它可以進(jìn)行大型基因組比如人類基因組的組裝，組裝結(jié)果更加準(zhǔn)確可靠，可以通過組裝的結(jié)果非常準(zhǔn)確地鑒別出基因組上的序列結(jié)構(gòu)性變異，為構(gòu)建全基因組參考序列和以低測(cè)序成本對(duì)未知基因組實(shí)施精確分析創(chuàng)造了可能。
程序的下載及安裝：
下載地址：http://soap./soapdenovo.html
安裝：
(a) 下載SOAPdenovo的壓縮包          
(b) 解壓縮     
(c)將得到可執(zhí)行文件SOAPdenovo和一個(gè)配置文件的模板e(cuò)xample.contig
 
1 使用程序及參數(shù)： 
SOAPdenovo可以一步跑完，也可以分成四步單獨(dú)跑
一步跑完的腳本:

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

四步單獨(dú)跑的腳本:

./ 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 參數(shù)說明

用法：/PathToProgram/SOAPdenovo all -s configFile [-K kmer -d KmerFreqCutOff -D EdgeCovCutoff -M mergeLevel -R -u -G gapLenDiff -L minContigLen -p n_cpu] -o Output

    -s    STR     配置文件

    -o    STR     輸出文件的文件名前綴

    -g    STR     輸入文件的文件名前綴

    -K    INT     輸入的K-mer值大小，默認(rèn)值23，取值范圍 13-63

    -p    INT     程序運(yùn)行時(shí)設(shè)定的線程數(shù)，默認(rèn)值8

    -R            利用read鑒別短的重復(fù)序列，默認(rèn)值不進(jìn)行此操作

    -d    INT     去除頻數(shù)不大于該值的k-mer，默認(rèn)值為0

    -D    INT     去除頻數(shù)不大于該值的由k-mer連接的邊，默認(rèn)值為1，即該邊上每個(gè)點(diǎn)的頻數(shù)都小于等于1時(shí)才去除

    -M    INT     連接contig時(shí)合并相似序列的等級(jí)，默認(rèn)值為1，最大值3。

    -F            利用read對(duì)scaffold中的gap進(jìn)行填補(bǔ)，默認(rèn)不執(zhí)行

    -u            構(gòu)建scaffold前不屏蔽高覆蓋度的contig，這里高頻率覆蓋度指平均contig覆蓋深度的2倍。默認(rèn)屏蔽

    -G    INT     估計(jì)gap的大小和實(shí)際補(bǔ)gap的大小的差異，默認(rèn)值為50bp。

    -L            用于構(gòu)建scaffold的contig的最短長(zhǎng)度，默認(rèn)為：Kmer參數(shù)值 ×2

3 使用方法及示例
（1）示例

SOAPdenovo all -s HCB.lib -K 25 -d -o test

（2） 輸入文件
configFile，配置文件內(nèi)容如下，非程序生成，需要軟件使用者自己配置。各個(gè)說明參考如下：

#maximal read length （read的最大長(zhǎng)度）

以“#”開頭的行是注釋內(nèi)容

max_rd_len=50    

#該值一般設(shè)置的比實(shí)際read讀長(zhǎng)稍微短一些，截去測(cè)序最后的部分，具體長(zhǎng)度看測(cè)序質(zhì)量

[LIB]

#文庫信息以此開頭

avg_ins=200

#文庫平均插入長(zhǎng)度，一般取插入片段分布圖中給出的文庫大小

reverse_seq=0

#序列是否需要被反轉(zhuǎn)，目前的測(cè)序技術(shù)，插入片段大于等于2k的采用了環(huán)化，所以對(duì)于插入長(zhǎng)度大于等于2k文庫，序列需要反轉(zhuǎn)，reverse_seq＝1，小片段設(shè)為0

asm_flags=3

#該文庫中的read序列在組裝的哪些過程（contig/scaff/fill）中用到

設(shè)為1：只用于構(gòu)建contig；

設(shè)為2：只用于構(gòu)建scaffold；

設(shè)為3：同時(shí)用于構(gòu)建contig和scaffold；

設(shè)為4：只用于補(bǔ)洞



【注意】短插入片段(<2K)的設(shè)為3，同時(shí)用于構(gòu)建contig和scaffold，長(zhǎng)插入片段(>=2k)設(shè)為2，不用于構(gòu)建contig，只用于構(gòu)建scaffold，454single 長(zhǎng)reads只用于補(bǔ)洞。

rank=1

#rank該值取整數(shù)，決定了reads用于構(gòu)建scaffold的次序，值越低，數(shù)據(jù)越優(yōu)先用于構(gòu)建scaffold。設(shè)置了同樣rank的文庫數(shù)據(jù)會(huì)同時(shí)用于組裝scaffold。一般將短插入片段設(shè)為1；2k設(shè)為2；5k設(shè)為3；10k設(shè)為4；當(dāng)某個(gè)檔的數(shù)據(jù)量較大時(shí)，也可以將其分為多個(gè)檔，同樣，當(dāng)某檔數(shù)據(jù)量不足夠時(shí)，可以將多個(gè)檔的數(shù)據(jù)合在一起構(gòu)建scaffold。這里說的數(shù)據(jù)量夠與不夠是從該檔的測(cè)序覆蓋度和物理覆蓋度兩個(gè)方面來考慮的。

pair_num_cutoff=3

#可選參數(shù)，pair_num_cutoff該參數(shù)規(guī)定了連接兩個(gè)contig 或者是pre-scaffold 的可信連接的閾值，即，當(dāng)連接數(shù)大于該值，連接才算有效。短插入片段(<2k)默認(rèn)值為3，長(zhǎng)插入長(zhǎng)度序列默認(rèn)值為5

map_len=32

#map_len該參數(shù)規(guī)定了在map過程中 reads和contig的比對(duì)長(zhǎng)度必須達(dá)到該值（比對(duì)不容mismacth和gap），該比對(duì)才能作為一個(gè)可信的比對(duì)?？蛇x參數(shù)，短插入片段(<2k)一般設(shè)置為32，長(zhǎng)插入片段設(shè)置為35，默認(rèn)值是K＋2。

q1=/path/**LIBNAMEA**/fastq_read_1.fq

#read 1的fastq格式的序列文件，“/path/**LIBNAMEA**/fastq_read_1.fq”為read的存儲(chǔ)路徑

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

#read 2的fastq格式的序列文件，與read1對(duì)應(yīng)的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

#單向測(cè)序得到的fastq格式的序列文件

f=/path/**LIBNAMEA**/fasta_read_single.fa

#單向測(cè)序得到的fasta格式的序列文件

p=/path/**LIBNAMEA**/pairs_in_one_file.fa

#雙向測(cè)序得到的一個(gè)fasta格式的序列文件

4 輸出文件及說明
SOAPdenovo 分四部分別對(duì)應(yīng)的輸出文件：

• pregraph  生成7個(gè)文件 *.kmerFreq  *.edge  *.preArc  *.markOnEdge  *.path  *.vertex  *.preGraphBasic
• contig       生成4個(gè)文件 *.contig  *.ContigIndex  *.updated.edge  *.Arc
• map          生成3個(gè)文件 *.readOnContig  *.peGrads  *.readInGap
• scaff          生成6個(gè)文件 *.newContigIndex  *.links  *.scaf  *.scaf_gap  *.scafSeq  *.gapSeq

 
*.contig：contig序列文件，fasta格式；
*.scafSeq：fasta格式的scaffold序列文件，contig之間的gap用N填充；
對(duì)于得到的*.scafSeq文件還需要用GapCloser去合并其中的gap，最后的contig文件則是對(duì)補(bǔ)洞之后的scaffold文件通過打斷N區(qū)的方法得到。
以上兩個(gè)文件是組裝結(jié)果中最主要的輸出。

*.scaf：包括scaffold中contig的詳細(xì)信息；
在scaffold行中包括scaffold名字、contig長(zhǎng)度和該scaffold長(zhǎng)度。
在contig行包括contig名字、contig在scaffold上的起始位置、正反鏈、長(zhǎng)度和contig間的鏈接信息
*.links：contig間的pair-end連接信息
*.readOnContig：reads在contig上的位置
*.peGrads： 主要可以通過調(diào)整本文件中的參數(shù)來顯示構(gòu)建scaffold所用到的插入片段庫的個(gè)數(shù)，總共要到的read數(shù)，最長(zhǎng)的read的長(zhǎng)度，每個(gè)庫對(duì)應(yīng)的哪些reads，rank設(shè)置，pair_num_cutoff設(shè)置。例如：

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個(gè)文庫，該文件則有n+1行，且按照文庫大小順序排列。
第1行中，第二三四列分別是 所用文庫，reads總數(shù)和組裝中用到的最長(zhǎng)的reads長(zhǎng)度。
第2行中，四列分別是文庫大小，文庫中的reads數(shù)目，該文庫reads用到的rank等級(jí)和該文庫中reads用到的pair_num_cutoff。
第3～n+1行，四列分別是文庫大小，文庫中的reads數(shù)目加上前面的文庫中的reads總數(shù)，該文庫reads用到的rank等級(jí)和該文庫中reads用到的pair_num_cutoff。
如果配置文件中沒有設(shè)置pair_num_cutoff，即使用默認(rèn)參數(shù)，則最后一列顯示為0。

對(duì)于SOAPdenovo的每個(gè)步驟都有日志文件輸出，要保存好日志文件，日志文件中包含有很多有用的信息。
 
SOAPdenovo日志輸出說明：
 
1）pregraph.log:  其中有很多的統(tǒng)計(jì)信息，包括構(gòu)建debruijn-graph時(shí)用到多少reads數(shù)，構(gòu)圖中生成了多少uniq的kmer以及設(shè)置-d參數(shù)后去除了多少kmer。
在pregraph中，可選參數(shù)有 –R –K –d  結(jié)果如：

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

3283081670 kmer removed

其中Kmer 數(shù)是取決于所設(shè)k值大小以及數(shù)據(jù)量，nodes數(shù)即特異性的kmer數(shù)目，當(dāng)nodes數(shù)目過高（一般和基因組大小差不多大小），可能是數(shù)據(jù)的錯(cuò)誤率比較高，也可能是存在雜合。若nodes數(shù)目偏小，并且kmer數(shù)目很多，則基因組本身可能存在一定的重復(fù)度。對(duì)于k值的選取，當(dāng)數(shù)據(jù)量充足時(shí)（>=40X），植物基因組一般采用大kmer會(huì)有比較好的效果，而對(duì)于動(dòng)物基因組，k值一般多取27和29則足夠。kmer removed表示的 –d 參數(shù)所去除的低頻的kmer。

2）contig.log: contig 中，可選參數(shù) –R –D –M，注，-R 參數(shù)的選定，必須pregraph和contig中同時(shí)選擇才有效。結(jié)果例子：

16430183 pairs found, 2334584 pairs of paths compared, 1674493 pairs merged

從merged的數(shù)量可以作為估計(jì)雜合以及測(cè)序錯(cuò)誤的程度。

sum up 1932549703bp, with average length 1170

the longest is 36165bp, contig N50 is 2871 bp,contig N90 is 553 bp

相關(guān)的統(tǒng)計(jì)信息，一般情況下，植物基因組的contig N90都在200bp左右，若N90高于200bp，則該基因組的scaffold構(gòu)建都不會(huì)有太大的問題。動(dòng)物基因組，scaffold構(gòu)建很少有出問題的。

3）map.log: 

Output 415219610 out of 1956217742 (21.2)% reads in gaps

1661094582 out of 1956217742 (84.9)% reads mapped to contigs

一般情況下，reads in gap的比例和map to contig 的比例總和大于1?？赡苁且?yàn)閞eads map到多個(gè)地方都被算在其中的原因。當(dāng)map to contig的比例很高（80%左右時(shí)），但是組裝效果并不很好，可能是重復(fù)序列比較多。reads in gap比例較高（大于40%），是因?yàn)榛蚪M組裝的較碎，gap區(qū)域較多。
map_len 默認(rèn)值=K+5，當(dāng)默認(rèn)值大于設(shè)置的map_len時(shí)，以默認(rèn)值為準(zhǔn)，當(dāng)默認(rèn)值小于map_len值時(shí)，設(shè)置的map_len為準(zhǔn)。

4）scaff.log:

average contig coverage is 23, 5832270 contig masked

構(gòu)建scaffold是對(duì)高頻覆蓋的contig進(jìn)行屏蔽（即頻率高于average contig coverage的兩倍的contig不用于構(gòu)建scaffold），從這里可以看出組裝的基因組一定的重復(fù)情況。

estimated PE size 162, by 40034765 pairs

on contigs longer than 173, 38257479 pairs found,SD=8, insert_size estimated: 163

173 是配置文件中該文庫的insertsize，163 是根據(jù)reads  map到contig上的距離的估計(jì)值，8是這個(gè)分布的標(biāo)準(zhǔn)偏差。一般考慮 比對(duì)上去的pair數(shù)目和SD值。若pair對(duì)數(shù)很多且SD值很小（小片段文庫數(shù)據(jù)不超過三位數(shù)，大片段文庫數(shù)據(jù)部超過500），那我們一般可以將配置文件中的文庫插入片段的值改對(duì)短插入片段文庫（<1k）的大小估計(jì)值，一般是比較準(zhǔn)確的，下次組裝以及補(bǔ)洞時(shí)應(yīng)根據(jù)這個(gè)值對(duì)原來配置文件中的insertsize信息做修正。對(duì)于大片段文庫（>=2K），因?yàn)槭前裷eads map到contig上，若最長(zhǎng)contig較短時(shí)，可能找不到成pair比對(duì)上去的reads，這時(shí)，無法估計(jì)文庫大小，需要自己將大片段一級(jí)一級(jí)的map到前一級(jí)的組裝結(jié)果上，然后再分析大片段文庫的插入片段大小。注，需要調(diào)整insertsize信息時(shí)，只需要修改* .peGrads文件中的第一列，然后刪除*.links文件，重新跑scaff這一步即可。即構(gòu)建scaffold時(shí)，主要是根據(jù)*.links文件的信息進(jìn)行連接。

Cutoff for number of pairs to make a reliable connection: 3

1124104 weak connects removed (there were 4773564 active cnnects))

Cutoff for number是在配置文件中設(shè)的pair_num_cutoff值，weak connects是低于這個(gè)值被認(rèn)定為無效的連接數(shù)，active connects是滿足cutoff的連接數(shù)，根據(jù)這個(gè)數(shù)值可對(duì)pair_num_cutoff做調(diào)整

Picked  25241 subgraphs,4 have conflicting connections

conflicting connections 是表示構(gòu)建scaffold時(shí)的矛盾數(shù)，矛盾數(shù)比較高（>100）時(shí)，可根據(jù)前面的有效連接數(shù)，適當(dāng)提高pair_num_cutoff值，即提高scaffold連接要求的最少關(guān)系數(shù)

182483 scaffolds&singleton sum up 1990259817bp, with average length 10906

the longest is 6561520bp,scaffold N50 is 836795 bp, scaffold N90 is 157667 bp

scaffold 統(tǒng)計(jì)信息，將是根據(jù)rank分梯度的統(tǒng)計(jì)

Done with 13301 scaffolds, 2161915 gaps finished, 2527441 gaps overall

-F 參數(shù)補(bǔ)洞的統(tǒng)計(jì)信息。

5 參數(shù)調(diào)整
一般組裝時(shí)需要調(diào)整的參數(shù)，主要分兩種：
一種是針對(duì)腳本中的參數(shù)改動(dòng)：如調(diào)整  -K  -R  -d  -D  -M
-K 值一般與基因組的特性和數(shù)據(jù)量相關(guān)，目前用到的SOAPdenovo軟件主要有兩個(gè)版本，grape1123和grape63mer，其中g(shù)rape1123是最新版的組裝軟件，K值范圍13-31，grape63mer是可以使用大kmer的組裝版本，K值范圍13-63。 
【經(jīng)驗(yàn)】：植物基因組的組裝采用大kmer效果會(huì)比較好（要求短片段reads長(zhǎng)度75bp），動(dòng)物基因組很少有用到大kmer后有明顯改進(jìn)效果的，且動(dòng)物基因組的組裝K值一般設(shè)置為27和29較多。
-R參數(shù)，對(duì)于動(dòng)物基因組，R參數(shù)一般不設(shè)置，植物基因組由于較多的repeat區(qū)，則設(shè)置R參數(shù)后，效果更好。注意，設(shè)置-R時(shí)，一般使用-M 的默認(rèn)值。（熊貓基因組組裝時(shí)得出的結(jié)論）
-M 參數(shù)，0-3,默認(rèn)值1。一般雜合率為千分之幾就設(shè)為幾。熊貓基因組組裝時(shí)-M 2 。
-d 參數(shù)，對(duì)于沒有糾錯(cuò)，沒有處理的質(zhì)量又較差的原始數(shù)據(jù)，kmer的頻數(shù)為1的很多的數(shù)據(jù)的組裝，一般設(shè)置為-d 1 則足夠。對(duì)于處理過，或者是測(cè)序質(zhì)量較好的數(shù)據(jù)，可以不用設(shè)置。數(shù)據(jù)量很多時(shí)，也可以以-d 參數(shù)去除部分質(zhì)量稍差的數(shù)據(jù)。
-D 參數(shù)，默認(rèn)為1，一般不用另行設(shè)置。
第二種，從map這一過程去調(diào)節(jié)參數(shù)。可以調(diào)整配置文件的map_len的值和調(diào)整文件*.peGrads。
當(dāng)文庫插入片段分布圖中文庫大小與實(shí)驗(yàn)給出的文庫大小差異很大時(shí)，調(diào)整*.peGrads文件中的插入片段大小。
根據(jù)每一檔數(shù)據(jù)的數(shù)據(jù)量去調(diào)整文庫的rank等級(jí)。當(dāng)該文庫的數(shù)據(jù)量很多或者是在構(gòu)建scaffold的過程中的沖突數(shù)很多時(shí)，可是適當(dāng)?shù)恼{(diào)大第四列 的pair_num_cutoff，把條件設(shè)置的更嚴(yán)一些。
6 內(nèi)存估計(jì)
SOAPdenovo的四個(gè)步驟消耗的內(nèi)存是不一樣的，其中第一步消耗的內(nèi)存最多，使用沒有糾錯(cuò)的的reads，(K<=31)第一步消耗的內(nèi)存在基因組大小的80－100倍左右，糾過錯(cuò)則在40－50倍左右，第二步相對(duì)消耗的內(nèi)存會(huì)少很多，第三步消耗的內(nèi)存是僅次于第一步的，在第一步的一半左右，第四步消耗的內(nèi)存也會(huì)比較少。對(duì)于CPU的使用，默認(rèn)是8個(gè)，如果申請(qǐng)內(nèi)存時(shí)申請(qǐng)一個(gè)計(jì)算節(jié)點(diǎn)的所有內(nèi)存測(cè)將CPU就設(shè)置為該計(jì)算節(jié)點(diǎn)的CPU個(gè)數(shù)充分利用計(jì)算資源，如果僅申請(qǐng)一個(gè)節(jié)點(diǎn)的部分內(nèi)存則根據(jù)實(shí)際情況考慮。對(duì)于大kemr(K>31)其內(nèi)存使用是(k<=31)的1.5倍左右，有時(shí)甚至更多，要充分估計(jì)內(nèi)存的使用，在第一次運(yùn)行的時(shí)候考慮不能太保守。
    
7 常見錯(cuò)誤
1）配置文件中read存儲(chǔ)路徑錯(cuò)誤
只輸出日志文件。
pregraph.log中的錯(cuò)誤信息：“Cannot open /path/**LIBNAMEA**/fastq_read_1.fq. Now exit to system...”
2）-g 后所跟參數(shù)與pregraph（第一步） -o  后所跟參數(shù)名不一致
contig  map  scaff 這三個(gè)步驟都只是輸出日志文件。
contig.log中的錯(cuò)誤信息：“Cannot open *.preGraphBasic. Now exit to system...”
map.log中的錯(cuò)誤信息：“Cannot open *.contig. Now exit to system...”
scaff.log中的錯(cuò)誤信息：“Cannot open *.preGraphBasic. Now exit to system...”
3）從map開始重新跑時(shí)，需要?jiǎng)h除*.links文件，否則會(huì)生成core文件，程序退出。