macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

De novoアセンブリ関連のツールまとめ(主にショートリード)

2020 3/8  整頓、help追記

随時更新

 

これまでNGSデータをアセンブル様々なアセンブルやツールが発表されてきた。例えばOMICs toolでヒットするツールは100を超える(https://omictools.com/genome-assembly-category)。全てを整理しきるのは難しいが、代表的と思われるいくつかのツールをインストールしてランした結果をまとめる。

 

 

追加分

主にショートリードを使うアセンブラ

1、ラージゲノムのアセンブラ

2、スモールゲノムのアセンブラ

3、ウィルスゲノムのアセンブラ

4、メタゲノムのアセンブラ

5、オルガネラゲノムのアセンブラ

6、プラスミドのアセンブラ

 

主にロングリードを使うアセンブラ

1、ラージ〜スモールゲノム

2、メタゲノムのアセンブラ

 

エラーコレクションツール

1、ロングリード

2、ショートリード

 

その他

scaffoldingツール

ガイドアセンブリ

 前処理やフィルタリング

 Polishing

評価ツール

シミュレータ

Repetitive sequence

Primary contigs

リード情報は使わないscaffoldingツール

Eukaryotic genome Annotation (including gene predictions)

 

 

wiki

f:id:kazumaxneo:20180714114551p:plain

GAGE-Aのレシピ (論文 http://genome.cshlp.org/content/22/3/557.long)

http://gage.cbcb.umd.edu/recipes/index.html

 

AbySS

de Bruijn graph based assembler

マニュアル

http://crusade1096.web.fc2.com/abyss.html

Github

https://github.com/bcgsc/abyss/#install-abyss-on-debian-or-ubuntu

 

インストールは brewやconda で行える。

#bioconda
conda create -n abyss -y
conda activate abyss
conda install -c bioconda -y abyss

#homebrew
breq install homebrew/science/abyss

brewでインストールするとmax kは64だが、手動でインストールするなら、例えば ./configure --enable-maxk=96 としてビルドすれば、k=96まで扱えるように修正できる。

ただしソースからビルドするならboostが必要になる。

> abyss --help

u$ abyss --help

Usage: ABYSS -k<kmer> -o<output.fa> [OPTION]... FILE...

Assemble the input files, FILE, which may be in FASTA, FASTQ,

qseq, export, SAM or BAM format and compressed with gz, bz2 or xz.

 

 Options:

 

      --chastity        discard unchaste reads [default]

      --no-chastity     do not discard unchaste reads

      --trim-masked     trim masked bases from the ends of reads

                        [default]

      --no-trim-masked  do not trim masked bases from the ends of

                        reads

  -q, --trim-quality=N  trim bases from the ends of reads whose

                        quality is less than the threshold

  -Q, --mask-quality=N  mask all low quality bases as `N'

  --standard-quality    zero quality is `!' (33)

                        default for FASTQ and SAM files

  --illumina-quality    zero quality is `@' (64)

                        default for qseq and export files

      --SS              assemble in strand-specific mode

      --no-SS           do not assemble in strand-specific mode

  -o, --out=FILE        write the contigs to FILE

  -k, --kmer=N          the length of a k-mer (when -K is not set) [<=128]

                        or the span of a k-mer pair (when -K is set)

  -K, --single-kmer=N   the length of a single k-mer in a k-mer pair

  -t, --trim-length=N   maximum length of blunt contigs to trim [k]

  -c, --coverage=FLOAT  remove contigs with mean k-mer coverage

                        less than this threshold

      --kc=N            remove all k-mers with multiplicity < N [0]

  -b, --bubbles=N       pop bubbles shorter than N bp [3*k]

  -b0, --no-bubbles     do not pop bubbles

  -e, --erode=N         erode bases at the ends of blunt contigs with coverage

                        less than this threshold [round(sqrt(median))]

  -E, --erode-strand=N  erode bases at the ends of blunt contigs

                        with coverage less than this threshold on

                        either strand [1 if sqrt(median) > 2 else 0]

  --coverage-hist=FILE  write the k-mer coverage histogram to FILE

  -m, --mask-cov        do not include kmers containing masked bases in

                        coverage calculations [experimental]

  -s, --snp=FILE        record popped bubbles in FILE

  -v, --verbose         display verbose output

      --help            display this help and exit

      --version         output version information and exit

      --db=FILE         specify path of database repository in FILE

      --library=NAME    specify library NAME for database

      --strain=NAME     specify strain NAME for database

      --species=NAME    specify species NAME for database

 

 ABYSS Options: (won't work with ABYSS-P)

 

  -g, --graph=FILE      generate a graph in dot format

 

Report bugs to <abyss-users@bcgsc.ca>.

 

ラン

abyss-pe name=bacteria j=16 k=64 n=5 s=100 in='R1.fq R2.fq'
  •  number of threads [2]
  •  size of k-mer (bp)
  •  minimum unitig size required for building contigs (bp) [200]
  •  minimum number of pairs required for building contigs [10]
  • in=' 'でinputするfastqを指定。

 forループを使い、複数のk-merを自動検討する。

for k in {65..96}; do ¥
ABYSS -k$k R1.fq -o contigs-k$k.fa ¥
done

 

 SOAPdenovo2

de Bruijn graph based assembler

公式マニュアル

SOAPdenovo: short-read assembly

公式マニュアル2

SOAP :: Short Oligonucleotide Analysis Package

 

conda create -n soapdenovo2 -y 
conda activate soapdenovo2
conda install -c bioconda -y soapdenovo2

 > SOAPdenovo-127mer -h

$ SOAPdenovo-127mer -h

 

Version 2.04: released on July 13th, 2012

Compile Dec 23 2013 10:46:03

 

Usage: SOAPdenovo <command> [option]

    pregraph        construct kmer-graph

    sparse_pregraph construct sparse kmer-graph

    contig          eliminate errors and output contigs

    map             map reads to contigs

    scaff           construct scaffolds

    all             do pregraph-contig-map-scaff in turn

> SOAPdenovo-63mer -h

$ SOAPdenovo-63mer -h

 

Version 2.04: released on July 13th, 2012

Compile Dec 23 2013 10:45:55

 

Usage: SOAPdenovo <command> [option]

    pregraph        construct kmer-graph

    sparse_pregraph construct sparse kmer-graph

    contig          eliminate errors and output contigs

    map             map reads to contigs

    scaff           construct scaffolds

    all             do pregraph-contig-map-scaff in turn

 

ラン

ランはconfig_fileファイルを作成し、それを元にして実行する。mate-pairデータも扱うことができ、 contigをNでつなぎscaffoldsを作ってくれる。paired-endデータを使った場合もこの作業は行われる。ただしメモリ喰いで、1GB以上のfastqだとサーバークラスのメモリが必要となる。64GB程度だと途中でストップする。khmerなどのツールであらかじめシーケンスデータをダイエットしておく必要があるかもしれない。

SOAPdenovo-127mer all -s config_file -k 127 -R -o output

 

 

 SPAdes

de Bruijn graph based assembler

現在もバージョンアップが続けられているbacteriaむけのアセンブラ。paired-endデータを扱えるが、それだけでなくロングリードを使いhybrid-assemblyを行うことも可能である。

spadesはハプロタイプのような部位もできる限り繋げようとする傾向が強い。また、contigをpaired-ned情報を使いscaffoldingするため、ほかと比べてかなり長いcontigができる。ただしこれは精度と性能の天秤にになりやすく、微妙な部位をむりやり繋げる恐れも高いかもしれない。例えば、monoploidの生き物でも、ゲノムがマルチコピーだとmelody ploideの状態のindelが発生している可能性はある。そのようなシーケンスデータでも、spadesはarbitraryな判断で繋ぐ傾向が強い。これは例えWT:MT=50:50でも起こる(そもそもspadesは純化したバクテリアゲノム向け)。

spadesはbrewでインストールできる。

ラン

spades.py -k auto -t 8 -1 R1.fq -2 R2.fq -o output_dir

 

 

Edena

string graph-based assembler 

公式サイト

http://computing.bio.cam.ac.uk/local/doc/edena.pdf

conda create -n edena -y 
conda activate edena
conda install -c bioconda -y edena

edena

$ edena

Edena v3.131028

Copyright (C) 2008,2011,2012,2013

David Hernandez, Patrice Francois, Jacques Schrenzel

Genomic Research Laboratory, Geneva University Hospitals, Switzerland

All rights reserved.

 

PROGRAM OPTIONS:

  1) Overlapping mode:

    -nThreads <int>         Number of threads to use. Default 2

    -r

    -singleEnd <file1> <file2> ...

                            Reads file(s) in FASTA or FASTQ format.

                            Several files can be specified

    -DRpairs

    -paired <file1_1> <file1_2> <file2_1> <file_2_2> ...

                            Direct-Reverse paired reads files. Several

                            pairs of files can be specified.

    -RDpairs

    -matePairs <file1_1> <file1_2> <file2_1> <file_2_2> ...

                            Reverse-Direct paired reads files. Several

                            pairs of files can be specified.

    -p

    -prefix <name>          Prefix for the output files. Default is "out".

    -M

    -minOverlap <int>       Minimum size of the overlaps to compute.

                            Default is half of the reads length.

    -t

    -truncate <int>         Truncate the 3' end of the reads TO the specified

                            length

  2) Assembler mode:

    -e

    -edenaFile <file.ovl>   Edena overlap (.ovl) file. Required.

    -p

    -prefix <name>          Prefix for the output files.

    -m

    -overlapCutoff <int>    Only consider overlaps >= than the specified size.

                            The optimal setting of this parameter depends on the

                            coverage that was achieved by the sequencing run.

                            You should therefore try different values in order

                            to get the optimal one.

    -c

    -minContigSize <int>    Minimum size of the contigs to output.

                            Default is 1.5*readLength.

    -cc

    -contextualCleaning

                 <yes/no>   Contextual cleaning of spurious edges.

    -minCoverage <int>      Minimum required coverage for the contigs. This

                            value is automatically determined if not specified.

    -sph

    -shortPeHorizon <int>   Maximum search distance for short >< paired-end

                            sampling. Default: 1000

    -lph

    -longPeHorizon <int>    Maximum search distance for long <> paired-end

                            sampling. Default: 15000

    -peHorizon <int>        obsolete: Maximum search distance for both short

                            and long paired-end reads sampling.

    -trimRed <yes/no>       By default, possible redundant sequences, caused by

                            unresolved repeats, are trimmed from contigs ends.

                            Setting this flag to 'no' allows keeping such 

                            redundancies up to the length of the largest insert

                            size (maxJump). !! setting this setting to 'no'

                            produces an artificially increased assembly size !!

    -maxRed <int>           Max ending redundancy length. Default: 0 (equivalent

                            to '-trimRed yes'. Overrides -trimRed.

    -d

    -deadEnds <int>         Maximum length for dead-end paths removal.

                            Default value is set to 2*readLength-1.

    -discardNonUsable

                 <yes/no>   Reads that are involved in orphan nodes smaller than

                            1.5*readLength are considered as "non-usable".

                            This filter discards such nodes. Default: enabled

    -trim <int>             Coverage cutoff for contigs ends:

                            Contig ends ending in a gap may contain errors due

                            to low coverage. This option trim a few bases from

                            these ends until a minimum coverage is reached.

                            Default is 4. Ends are not trimmed if set to 1.

    -shell                  Interactive shell. Allows querying the assembly

                            graph.

REPORT BUGS:

   david.hernandez@genomic.ch

 

ラン

fastq配列の長さが同じでないと動かないので、短いリードを除去し、さらにスクリプトを書いて整形した。

edena -paired R1.fq R2.fq

 out.ovlが出力される。

edena -o out.ovl

mオプション次第でかなり変わる。-mの値は100-249まで選べるが、Miseq 301bpだと150前後がベストな気がする。

 

 

A5-miseq

quality check、data trimming、エラーコレクション、contig生成、scaffold構築、データビジュル化を全自動でやってくれるのが特徴のアセンブラである。A5-miseqは手法としての新しさよりも、パイプラインを自動化することに重きを置いているため、コマンドラインに不慣れなユーザーでも比較的手軽に扱うことができる。

ダウンロード

https://sourceforge.net/projects/ngopt/

conda create -n a5-miseq -y 
conda activate a5-miseq
conda install -c bioconda -y a5-miseq

a5_pipeline.pl -h

$ a5_pipeline.pl -h

 

A5-miseq version 20160825

Usage: a5_pipeline.pl [--begin=1-5] [--end=1-5] [--preprocessed] [--threads=4] [--debug] [--metagenome] <lib_file> <out_base>

 

Or: a5_pipeline.pl <Read 1 FastQ> <Read 2 FastQ> <out_base>

 

Or: a5_pipeline.pl <Read 1,2 Interleaved FastQ> <out_base>

 

<out_base> is the base file name for all output files. When assembling from 

a single library, the fastq files may be given directly on the command line.

If using more than one library, a library file must be given as <lib_file>.

The library file must contain the filenames of all read files.

 

If --preprocessed is used, <lib_file> is expected to be the library file

created during step 2 of the pipeline, named <out_base>.preproc.libs. Note 

that this flag only applies if beginning pipeline after step 2.

 

 docker版があり、仮想環境で気軽にテストもできる(docker-A5-Miseq)。

参考動画

 

 Velvet

インストール

#bioconda
conda create -n velvet -y
conda activate velvet
conda install -c bioconda -y velvet

#from source
git clone https://github.com/dzerbino/velvet.git
cd velvet/
vi MakeFile #開く。31mer以上でもアセンブル用にする。別のエディタで開いてもOK
6行目のMAXKMERLENGTH?=31を => MAXKMERLENGTH?=250 とかに修正。
#保存して閉じる。

make #makeはデフォルトでMakefileを自動認識する
パスを通しておく。

velvetg

$ velvetg

velvetg - de Bruijn graph construction, error removal and repeat resolution

Version 1.2.10

Copyright 2007, 2008 Daniel Zerbino (zerbino@ebi.ac.uk)

This is free software; see the source for copying conditions.  There is NO

warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Compilation settings:

CATEGORIES = 4

MAXKMERLENGTH = 191

OPENMP

LONGSEQUENCES

 

Usage:

./velvetg directory [options]

 

directory : working directory name

 

Standard options:

-cov_cutoff <floating-point|auto> : removal of low coverage nodes AFTER tour bus or allow the system to infer it

(default: no removal)

-ins_length <integer> : expected distance between two paired end reads (default: no read pairing)

-read_trkg <yes|no> : tracking of short read positions in assembly (default: no tracking)

-min_contig_lgth <integer> : minimum contig length exported to contigs.fa file (default: hash length * 2)

-amos_file <yes|no> : export assembly to AMOS file (default: no export)

-exp_cov <floating point|auto> : expected coverage of unique regions or allow the system to infer it

(default: no long or paired-end read resolution)

-long_cov_cutoff <floating-point>: removal of nodes with low long-read coverage AFTER tour bus

(default: no removal)

 

Advanced options:

-ins_length* <integer> : expected distance between two paired-end reads in the respective short-read dataset (default: no read pairing)

-ins_length_long <integer> : expected distance between two long paired-end reads (default: no read pairing)

-ins_length*_sd <integer> : est. standard deviation of respective dataset (default: 10% of corresponding length)

[replace '*' by nothing, '2' or '_long' as necessary]

-scaffolding <yes|no> : scaffolding of contigs used paired end information (default: on)

-max_branch_length <integer> : maximum length in base pair of bubble (default: 100)

-max_divergence <floating-point>: maximum divergence rate between two branches in a bubble (default: 0.2)

-max_gap_count <integer> : maximum number of gaps allowed in the alignment of the two branches of a bubble (default: 3)

-min_pair_count <integer> : minimum number of paired end connections to justify the scaffolding of two long contigs (default: 5)

-max_coverage <floating point> : removal of high coverage nodes AFTER tour bus (default: no removal)

-coverage_mask <int> : minimum coverage required for confident regions of contigs (default: 1)

-long_mult_cutoff <int> : minimum number of long reads required to merge contigs (default: 2)

-unused_reads <yes|no> : export unused reads in UnusedReads.fa file (default: no)

-alignments <yes|no> : export a summary of contig alignment to the reference sequences (default: no)

-exportFiltered <yes|no> : export the long nodes which were eliminated by the coverage filters (default: no)

-clean <yes|no> : remove all the intermediary files which are useless for recalculation (default : no)

-very_clean <yes|no> : remove all the intermediary files (no recalculation possible) (default: no)

-paired_exp_fraction <double> : remove all the paired end connections which less than the specified fraction of the expected count (default: 0.1)

-shortMatePaired* <yes|no> : for mate-pair libraries, indicate that the library might be contaminated with paired-end reads (default no)

-conserveLong <yes|no> : preserve sequences with long reads in them (default no)

 

Output:

directory/contigs.fa : fasta file of contigs longer than twice hash length

directory/stats.txt : stats file (tab-spaced) useful for determining appropriate coverage cutoff

directory/LastGraph : special formatted file with all the information on the final graph

directory/velvet_asm.afg : (if requested) AMOS compatible assembly file

 

ラン

1、velvethで処理するためのデータを準備。

velveth folder-name 31,45,2 -fastq -shortPaired1 input1.fastq -shortPaired2 input2.fastq #31,45,2はハッシュ値が31から45まで2ずつ検討するということ

シングルなら

velveth directory 31,45,2 -fastq single.fq 

 

2、velvetgでコンティグを作成 (例えばk=31に対して)

velvetg folder-name_31 -ins_length1 300 -ins_length2 3000
  •   -exp_cov <floating point|auto> expected coverage of unique regions or allow the system to infer it -exp_cov <floating point|auto> : expected coverage of unique regions or allow the system to infer it (default: no long or paired-end read resolution)
  •  -ins_length* <integer> expected distance between two paired-end reads in the respective short-read dataset (default: no read pairing) -ins_length* <integer> : expected distance between two paired-end reads in the respective short-read dataset (default: no read pairing)

 

追記

 

 

 以下のツールはデータの関係で使えなかった。

 

 

既に昔のデータになってしまっているが、他のアセンブルツールのコマンドについてはGAGEやAssemblathonのペーパー参照。

GAGE recipes

GAGE recipes

GAGE-B recipes

https://ccb.jhu.edu/gage_b/recipes/recipes.pdf

Assemblathon2

https://gigascience.biomedcentral.com/articles/10.1186/2047-217X-2-10

13742_2013_29_MOESM3_ESM.docx

 

 各ツールの特徴

Next Generation Sequencing (NGS)/De novo assembly - Wikibooks, open books for an open world

 

エラーコレクションツールのベンチマーク

https://biorxiv.org/cgi/content/short/2020.03.06.977975v1