macでインフォマティクス

macでインフォマティクス

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

De novoアセンブリ関連のツールまとめ

assembly summary

2020 3/8  整頓、help追記, 2024/01/22 タイトルからショートリード削除

随時更新

 

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

 

追加分

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

１、ラージゲノムにも対応したアセンブラ

２、スモールゲノム向けのアセンブラ

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

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

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

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

7、シード配列から拡張

８，その他

 

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

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

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

３，ターゲットアセンブラ

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

 

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

１、ロングリード

２、ショートリード

 

その他

scaffoldingツール

ガイドアセンブリ

 前処理やフィルタリング

Merge

Polishing

評価ツール

シミュレータ

Repetitive sequence

Primary contigs

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

Ab initio gene prediction

Evidence-driven gene prediction

EVidenceModeler

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'
  • j  number of threads [2]
  • k  size of k-mer (bp)
  • s  minimum unitig size required for building contigs (bp) [200]
  • n  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

 

mamba create -n soapdenovo2 -y 
conda activate soapdenovo2
mamba install -c conda-forge -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-127mer all

$ SOAPdenovo-127mer all

 

Version 2.04: released on July 13th, 2012

Compile Feb 14 2018    20:39:13

 

SOAPdenovo all -s configFile -o outputGraph [-R -F -u -w] [-K kmer -p

n_cpu -a initMemoryAssumption -d KmerFreqCutOff -D EdgeCovCutoff -M

mergeLevel -k kmer_R2C, -G gapLenDiff -L minContigLen -c minContigCvg

-C maxContigCvg -b insertSizeUpperBound -B bubbleCoverage -N

genomeSize]

  -s <string>    configFile: the config file of solexa reads

  -o <string>    outputGraph: prefix of output graph file name

  -K <int>       kmer(min 13, max 127): kmer size, [23]

  -p <int>       n_cpu: number of cpu for use, [8]

  -a <int>       initMemoryAssumption: memory assumption initialized

to avoid further reallocation, unit G, [0]

  -d <int>       kmerFreqCutoff: kmers with frequency no larger than

KmerFreqCutoff will be deleted, [0]

  -R (optional)  resolve repeats by reads, [NO]

  -D <int>       edgeCovCutoff: edges with coverage no larger than

EdgeCovCutoff will be deleted, [1]

  -M <int>       mergeLevel(min 0, max 3): the strength of merging

similar sequences during contiging, [1]

  -e <int>       arcWeight: two edges, between which the arc's weight

is larger than arcWeight, will be linerized, [0]

  -m <int>       maxKmer (max 127): maximum kmer size used for multi-kmer, [NO]

  -E (optional)  merge clean bubble before iterate, works only if -M

is set when using multi-kmer, [NO]

  -k <int>       kmer_R2C(min 13, max 127): kmer size used for mapping

reads to contigs, [K]

  -F (optional)  fill gaps in scaffolds, [NO]

  -u (optional)  un-mask contigs with high/low coverage before

scaffolding, [mask]

  -w (optional)  keep contigs weakly connected to other contigs in

scaffold, [NO]

  -G <int>       gapLenDiff: allowed length difference between

estimated and filled gap, [50]

  -L <int>       minContigLen: shortest contig for scaffolding, [K+2]

  -c <float>     minContigCvg: minimum contig coverage (c*avgCvg),

contigs shorter than 100bp with coverage smaller than c*avgCvg will be

masked before scaffolding unless -u is set, [0.1]

  -C <float>     maxContigCvg: maximum contig coverage (C*avgCvg),

contigs with coverage larger than C*avgCvg or contigs shorter than

100bp with coverage larger than 0.8*C*avgCvg will be masked before

scaffolding unless -u is set, [2]

  -b <float>     insertSizeUpperBound: (b*avg_ins) will be used as

upper bound of insert size for large insert size ( > 1000) when

handling pair-end connections between contigs if b is set to larger

than 1, [1.5]

  -B <float>     bubbleCoverage: remove contig with lower cvoerage in

bubble structure if both contigs' coverage are smaller than

bubbleCoverage*avgCvg, [0.6]

  -N <int>       genomeSize: genome size for statistics, [0]

  -V (optional)  output information for Hawkeye to visualize the assembly, [NO]

 

ラン

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

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

 

 

 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

mamba create -n edena -y 
conda activate edena
mamba 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
６行目の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

 

ラン

１、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

 

２、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)

 

 

真核生物ゲノムの遺伝子予測

遺伝子予測の実践的な手順についてはこちらがよくまとまっています。


基礎的な知識を学ぶならこちら

https://www.science.smith.edu/cmbs/wp-content/uploads/sites/36/2015/09/euk_genome_annotation_review.pdf

 

PASAに記載されている真核生物アノテーションワークフローの手順

  1. 以下のソフトウェアツールを用いたab initio遺伝子探索。GeneMarkHMM、FGENESH、Augustus、およびSNAP、GlimmerHMM
  2. GeneWiseソフトウェアとuniref90非冗長タンパク質データベースを用いたタンパク質相同性検出とイントロン分析
  3. 既知のEST、完全長cDNA、最近ではTrinity RNA-Seqアセンブリのゲノムへのアラインメント
  4. ステップ( C)で得られたオーバーラップする転写産物アラインメントに基づくPASAアラインメントアセンブリ
  5. EVidenceModeler (EVM) を用いて、上記 (A, B, C, D) に基づく重み付きコンセンサス遺伝子構造アノテーションの計算
  6. PASAを使用してEVMコンセンサス予測を更新し、UTRアノテーションと交互スプライシングアイソフォームのモデルを追加（DとEを利用）
  7. Argo または Apollo を用いたゲノムアノテーション (F) の限定的な手動改良

 

追記

 

 

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

 

 

既に昔のデータになってしまっているが、他のアセンブルツールのコマンドについては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

 

2020 10/31 

バクテリアゲノムアセンブリのガイド

 

2020 12/13

 

2020 12/27

 

2021 1/27

 

(Published: November 12, 2020)

 

2021 08/01


2022 06/25

バイオインフォマティクスの前提を崩すような奇妙な遺伝子アノテーションのリスト

 

2023

"様々なヘテロ接合度を持つ6つのゲノムを用いて、ゲノムサイズやヘテロ接合度などのゲノム特性の推定、de novoアセンブル、ポリッシング、対立遺伝子コンティグの除去などの一連のプロセスを評価した。5つのロングリードのみのアセンブラ（Canu、Flye、miniasm、NextDenovo、Redbean）と、ショートリードとロングリードを組み合わせた5つのハイブリッドアセンブラ（HASLR、MaSuRCA、Platanus-allee、SPAdes、WENGAN）を評価し、安定で高性能なアセンブラを用いて、ヘテロ接合の程度に応じたハプロタイプ表現の構築、それに続くハプロティグのポリッシングとパージについて、具体的なガイドラインを提案した。"

A practical assembly guideline for genomes with various levels of heterozygosity 
Takako Mochizuki, Mika Sakamoto, Yasuhiro Tanizawa, Takuro Nakayama, Goro Tanifuji, Ryoma Kamikawa, Yasukazu Nakamura
Briefings in Bioinformatics, Volume 24, Issue 6, November 2023, bbad337

 

 

https://academic.oup.com/bib/article/24/6/bbad337/7291993?login=false