soapdenovo [生物資訊實驗室]
 

Introduction

SOAPdenovo為SOAP(Short Oligonucleotide Analysis Package)相關應用程式中,用來進行de novo assembly的工具,

目前只開放編譯好的執行檔,最新版本為1.04版(2009/12/21 updated),

相關資訊可以在SOAPdenovo網頁查閱 http://soap.genomics.org.cn/soapdenovo.html

SOAPdenovo有三個相關tools

Correction tool

Correction tool是利用quality score和k-mer frequency對pair-end reads做修正的動作,

將Correction tool解開後會包含4支程式,

主要會用到的有KmerFreq、Corrector和merge_pair.pl這三支程式.

KmerFreq

KmerFreq就是在統計所有K-mer出現的次數。

Usage:
       KmerFreq [options]
       -i      <string>        input a file list
       -o      <string>        output filename prefix
       -q      <int>   quality cutoff [default 5]
       -s      <int>   seed length [default 17]
       -n      <int>   output kmer index along with frequency ? 0: no, 1: yes.[default 0]
       -f      <int>   file format: 1: fq, 2: fa.[default 1]
       -h      -?      help

通常會用到的參數是-i、-o和-s,其他我都會使用default值,

特別提一下-s這個參數,-s後要接一個K-mer length,這邊不建議設太大,

因為程式會先跟系統要一塊可以儲存所有K-mer組合的記憶體,

簡單來說,17-mer需要4的17次方大小的記憶體空間(也就是16G),

設越大越容易將記憶體吃垮,所以這邊要先了解系統的記憶體大小後再謹慎的選擇K-mer size;

-i是要input一個”file list“,

舉例來說,假設你現在有2個存放read的檔案要進行correction的動作,

其檔名分別為file1.fq和file2.fq,路徑存放在/home/user/下,

這時候你需要新開一個file_list檔,內容放你要correct的路徑+檔名:

/home/user/file1.fq
/home/user/file2.fq

最後在執行程式時將file_list加到-i參數後面就大功告成了。

  • Eample:
./KmerFreq -i file_list -o test -s 17
  • Output:

假設output的prefix為”test”,會產生兩個檔,

一個為test.freq另一個為test.stat,

test.freq是一個binary檔,在Corrector時會用到,

test.stat統計所有K-mer frequency的個數。

Corrector

Corrector利用KmerFreq統計出來的K-mer frequency將read上不正確的核甘酸做修正。

Usage:
       Corrector [options]
       -i      <string>        input a read file list
       -r      <string>        input kmer frequence file name
       -n      <int>   kmer frequency along with index?: 0: no, 1: yes. [default 0]
       -k      <int>   start of kmer frequence cutoff [default 5]
       -e      <int>   end of kmer frequence cutoff [default 5]
       -d      <int>   maximum error bases allowed [default 2]
       -s      <int>   seed length [default 17]
       -t      <int>   thread number [default 4]
       -f      <int>   file format: 1: fq, 2: fa. [default 1]
       -h      -?      help 

主要用到-i、-r、-s和-t這四的參數,

-i跟KmerFreq一樣,需要輸入一個file list,這裡就不再多做說明;

-r需要KmerFreq產生的.freq檔當做input;

-s設跟KmerFreq所使用的值即可;

-t為thread數目,可以看系統架構適度平行化。

  • Eample:
./Corrector -i file_list -r test.freq -s 17 -t 8
  • Output:

會根據file_list中的路徑和檔名,將修正過後的reads存放在同一路徑下,並加上.corr以示區別,

假設file_list中欲修正兩個檔(路徑為/home/user/):

/home/user/file1.fq
/home/user/file2.fq

Corrector完成後會產生”新”的兩個檔file1.fq.corr和file2.fq.corr(路徑也為/home/user/)。

merge_pair.pl

merge_pair是要將兩個分散的fasta檔案做合併,

通常pair-end兩端的reads會存在兩個不同的檔案,

而在做完Corrector後,有些reads可能會因為quality過低太多N而被砍掉,

所以兩端的reads個數可能會不一樣多(有些reads因為被刪除而變成single-end),

但是SOAPdenovo又無法處理這樣的情形(會有Floating Point Exception),

這時候merge_pair會將兩端都完整的paired-end取出來之外,

還會把剩下的single-end reads存到另一個檔案中。

Usage: <pairFile1> <pairFile2> <outputFile>

<pairFile1>會放左端修正過後的read檔,同理<pairFile2>則是放右端修正過後的read檔,

如果follow上面的例子,以下為執行的example:

  • Example:
./merge_pair.pl /home/user/file1.fq.corr /home/user/file2.fq.corr file_merge
  • Output:

正確執行結束後,會產生file_merge.pair、file_merge.single和file_merge.readsum這三個檔,

file_merge.pair會記錄完整的paired-end reads,格式如下:

>Read1左端header
Read1左端sequence
>Read1右端header
Read1右端sequence
>Read2左端header
Read2左端sequence
>Read2右端header
Read2右端sequence
...

而file_merge.single則是記錄那些另一端已經被砍掉的reads;

file_merge.readsum則是記錄paired-end和single-end reads各幾條。

SOAPdenovo

Config

Read data config(example.config):
#maximal read length
max_rd_len=36
[LIB]
#average insert size
avg_ins=0
#if sequence needs to be reversed
reverse_seq=0
#in which part(s) the reads are used
asm_flags=3
#in which order the reads are used while scaffolding
rank=1
#fastq file for read
q=/home/w40/DNA-seq/k562.fastq

Execution

SOAPdenovo configuration:
./SOAPdenovo pregraph -s example.config -o k562_21 -p 16 -d -K 21
./SOAPdenovo contig -g k562_21
./SOAPdenovo map -s example.config -g k562_21 -p 16 -K 21
./SOAPdenovo all -s config -o k45_d -K 45 -d -p 8
# SOAPdenovo可分好幾階段,下 all 可以全部一起跑
# -s config file
# -o output prefix
# -K k-mer
# -d remove k-mer freq <= 1
# -p mutithread
# -a 80 限制記憶體使用量(以G為單位),超過限制值,aborted

GapCloser

Usage:
       GapCloser [options]
       -o      <string>        Output file
       -b      <string>        Input configure file
       -a      <string>        Input scaffold file
       -p      <int>   Overlap length/K value(default 25, max 31)
       -t      <int>   Thread number
       -h      -?      Help
 
soapdenovo.txt · 上一次變更: 2014/06/05 16:33 (external edit)
 
Recent changes RSS feed Creative Commons License Donate Powered by PHP Valid XHTML 1.0 Valid CSS Driven by DokuWiki