2020 3/8 整頓、help追記, 2024/01/22 タイトルからショートリード削除
随時更新
これまでNGSデータをアセンブル 様々なアセンブル やツールが発表されてきた。例えばOMICs toolでヒットするツールは100を超える(https://omictools.com/genome-assembly-category )。全てを整理しきるのは難しいが、代表的と思われるいくつかのツールをインストールしてランした結果をまとめる。
追加分
主にショートリードを使うアセンブラ
1、ラージゲノムにも対応したアセンブラ
2、スモールゲノム向けのアセンブラ
3、ウィルスゲノムのアセンブラ
4、メタゲノムのアセンブラ
5、オルガネラゲノムのアセンブラ
6、プラスミドのアセンブラ
7、シード配列から拡張
8,その他
主にロングリードを使うアセンブラ
1、ラージ〜スモールゲノム
2、メタゲノムのアセンブラ
3,ターゲットアセンブラ
4、プラスミドのアセンブラ
エラーコレクションツール
1、ロングリード
2、ショートリード
その他
scaffoldingツール
ガイドアセンブリ
前処理やフィルタリング
Merge
Polishing
評価ツール
シミュレータ
Repetitive sequence
Primary contigs
リード情報は使わないscaffoldingツール
Ab initio gene prediction
Evidence-driven gene prediction
EVidenceModeler
Eukaryotic genome Annotation (including gene predictions)
wiki
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 )。
参考動画
VIDEO
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)
真核生物ゲノムの遺伝子予測
遺伝子予測の実践的な手順についてはこちらがよくまとまっています。
基礎的な知識を学ぶならこちら
https://www.science.smith.edu/cmbs/wp-content/uploads/sites/36/2015/09/euk_genome_annotation_review.pdf
PASAに記載されている真核生物アノテーション ワークフローの手順
以下のソフトウェアツールを用いたab initio遺伝子探索。GeneMarkHMM、FGENESH、Augustus、およびSNAP、GlimmerHMM
GeneWiseソフトウェアとuniref90非冗長タンパク質データベースを用いたタンパク質相同性検出とイントロン 分析
既知のEST、完全長cDNA、最近ではTrinity RNA -Seqアセンブリ のゲノムへのアラインメント
ステップ( C)で得られたオーバーラップする転写産物アラインメントに基づくPASAアラインメントアセンブリ
EVidenceModeler (EVM) を用いて、上記 (A, B, C, D) に基づく重み付きコンセンサス遺伝子構造アノテーション の計算
PASAを使用してEVMコンセンサス予測を更新し、UTRアノテーション と交互スプライシング アイソフォームのモデルを追加(DとEを利用)
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