2019 4/12 dockerリンク追加
NGSデータを分析するためのシーケンスアライメント、アセンブリ、変異検出、またはそれらのいくつかの組み合わせは、通常、バイオインフォマティクスパイプラインの主要なモジュールである[論文より ref.1,2,3,4,5,6]。微生物ゲノムシーケンシングの重要な用途は、食物サプライチェーン[ref.7,8,9]および病院[ref.10,11,12,13]での病原菌の大発生を検出することである。NGSを使った食品媒介性病原菌のサーベイランスおよびアウトブレイク調査の利点およびバイオインフォマティクスの課題を、Listeria Monocytogenesを例として、またリアルタイムアウトブレイク[ref.15]例を引用して検討した。両方のレビューでは、情報を使用する上での重大な挑戦であるNGSの新しいアセンブリが確認された。
米国の州、連邦機関、および国際パートナー間の協力で、食品由来病原菌のシーケンシングデータをsubmitするNCBIのPathogen Detection Project （PDP）によりNGSベースのアウトブレイクの調査が加速している。サルモネラ、リステリア、エシェリヒアおよびシゲラ、カンピロバクターの4つの主要な食品媒介性病原菌でsubmitされたシーケンシングデータセットの数は、2016年の初めにPDPによって発表された44017から、2018年の85823に急拡大している。
(一部略）
シーケンシングリードのde-novoアセンブラがいくつか公開されている[ref.21,22,23,24,25,26,27,28]。ploidity[ref.29]やメタゲノム[ref.30,31,32,33,34,35]、シングルセル[ref.36]、シーケンシング技術[ref.37]、またはいくつかのアセンブリを1つに結合することに特化してたりと[ref.38]、それぞれ特徴がある。アセンブラは、一倍体ゲノムであっても「エラーのない」アセンブリを保証しない。微生物ゲノムに加えて、ハプロイドアセンブリは、ヒトの特別なケース: hydatidiform mole (胞状奇胎) のような場合に重要である[39]。 PDPのようなユースケースでは、微生物ゲノムのアプリケーションは、母系から娘細胞へのゲノムデータの垂直継承に依存し、4Mbゲノムではほとんど変化のない、典型的には変異が10未満ということに基づいてリードセットクローン性パターンを解決する[ref.15]。このようなアプリケーションでは、非常に高いシーケンス品質のアセンブリが必要となるが、真のバリエーションを確実に検出することができる。
シーケンシングマシンによって生成されたリードセットを使用してアセンブリを評価するには、同じサンプルのリードとほぼ完全な高品質のドラフトアセンブリの両方を含む公的に利用可能なベンチマークセットが必要になる。 FDA-ARGOS（HP link）は、Food and Drug Administration（US FDA）[ref.19]によって開発された、微生物ゲノムのアセンブリ品質評価のベンチマーク要件を満たすデータベースであり、バクテリアのregulatory grade のシーケンスから構成されている。
ここでは、イルミナシーケンシング技術を使用して生成された微生物ゲノムのリードを迅速かつ高品質にデノボアセンブリするSKESA [skee-sa]（strategic k-mer extension for scrupulous assemblies)）を紹介する。 SKESAで使用されるHeuristicsは、Illuminaシーケンシングの低レベルのコンタミネーションやストランド固有のエラーによるアセンブリ品質への影響が低減されるように設計されている。エラー率の高い他のシーケンシングテクノロジでは、SKESAで使用される控えめなヒューリスティックにより、他のアセンブラで生成されたものよりアセンブリのcontiguityは低くなる。 SKESAは微生物ゲノムより大きなゲノムもアセンブリできるが、他のアセンブラと比較して、そのようなゲノムのためのプロファイリングはされていない。例えば、SKILSAは約40Mbの長さであるMonilinia fructigenaのSRR7262862データ（49.8 million reads with total length 9.2 Gb）とMonilinia laxaのSRR6748693データ（73.2 million reads with total length 18.3 Gb）を、10コア使用で3時間以下でアセンブリし、100コア使用で30分でアセンブリする。
多くのデノボアセンブラ、例えばSPAdes [ref.24]やMegaHit [ref.35]は、アセンブリ時に複数のk-merサイズを使用する。 ALLPATHS_LG [ref.40]は、100塩基対が40塩基オーバーラップしている特別短いインサートライブラリープロトコルを使用する。オーバーラップを利用して160 bpのマージされたリードが生成されたが、最大k-merサイズは96が使用された。 SKESAの顕著な特徴は、リードのサブセットのミニアセンブリを使ってペアエンドのインサートサイズまでの長いk-merを生成することである。リードより長いk-merを使用することで、SKESAは、インサートサイズよりも短いが、リードよりも長いリピート領域を正確にアセンブリできる。これに対して、現在のアセンブラはすべて、ペアエンドのmateリードサイズまでのk-mersだけ使用している。
本稿では、SKESAとSPAdes、MegaHitを、5つのタイプの微生物検査セットに対して適用した。（一部略） MegaHitは非常に高速なアセンブラであるため、SPAdesは、さまざまなテクノロジやサンプルのオプションを提供する多目的で広く使用されているアセンブラのため、この２つのアセンブラを選択した。ラージゲノムとスモールゲノムの両方をアセンブリできるよう設計されたABySS [ref.28]（バージョン2.0.2）と不均一なカバレッジを扱うように設計されたIDBA [36]（バージョン1.1.1）もテストされた（Additional file 1）。
（一部略）
著者らは、微生物ゲノムのアセンブリでは、SKESAとMegaHitは同じくらい速く、SPAdesよりもはるかに速いことを示す。 SKESAはSRAから直接アクセスすることもできる。これにより、ファイルからの入力を読み取るよりも高速に処理できる。 QUAST [ref.42]によって計算された100Kbあたりのミスマッチ数、N50統計量によって測定されたアセンブリの連続性、およびリファレンスアセンブリの長さからのずれによって測定されたアセンブリ品質は、SKESAアセンブリの品質がSPAdesとMegaHitの品質より良好であることを示している。同じ入力では、スレッド数、メモリ、または実行回数に関係なく、同じ結果が得られる。大量のデータを処理し、回帰テストを必要とする本番システムでは、これは非常に重要な要件である。著者らのテストでは、SPAdesとMegaHitの両方は、同じ数のスレッドとメモリの設定であっても、繰り返しにより異なるアセンブリを生成していた。したがって、SKESAは、高い基本レベルの品質と下流の分析に十分なcontiguityを備え、低レベルの汚染を扱い、迅速かつ高精度なアセンブリを作成するすべての要件を満たしている。すなわち生産環境でrobustである。 SKESAは現在、SRAの微生物ゲノムを組み立てるためにNCBIのproductionに使用されており、PDPのワークフローに組み込まれている。 SKESAのソフトウェアは自由に利用可能であり[ref.43、44]（see “Availability and requirements”）、クラウドでも利用できるようにする予定である。
New Docker image available: SKESA
— Curtis Kapsak (@CurtisKapsak) April 11, 2019
'docker pull staphb/skesa'
'docker run -v ${PWD}:/data staphb/skesa:latest skesa --fastq /data/Ecoli-PE-interleaved.fastq.gz --use_paired_ends --contigs_out /data/Ecoli.contigs.fasta'https://t.co/WfbHGUfTcr#staphb #docker #bioinformatics
インストール
mac os10.12のanaconda2-4.2.0 環境でテストした。
git clone https://github.com/ncbi/SKESA
cd SKESA/
#If you would like to build NGS library for accessing reads from SRA, then do
make
#Otherwise, if reading inputs only from files, do
make -f Makefile.nongs
#またはcondaで導入(sraの直接アセンブリ機能はなし）
conda install -c bioconda skesa
> skesa
$ skesa
skesa
Provide some input reads
General options:
-h [ --help ] Produce help message
-v [ --version ] Print version
--cores arg (=0) Number of cores to use (default all) [integer]
--memory arg (=32) Memory available (GB, only for sorted counter)
[integer]
--hash_count Use hash counter [flag]
--estimated_kmers arg (=100) Estimated number of unique kmers for bloom
filter (M, only for hash counter) [integer]
--skip_bloom_filter Don't do bloom filter; use --estimated_kmers as
the hash table size (only for hash counter)
[flag]
Input/output options : at least one input providing reads for assembly must be specified:
--fasta arg Input fasta file(s) (could be used multiple
times for different runs) [string]
--fastq arg Input fastq file(s) (could be used multiple
times for different runs) [string]
--use_paired_ends Indicates that a single (not comma separated)
fasta/fastq file contains paired reads [flag]
--contigs_out arg Output file for contigs (stdout if not
specified) [string]
Assembly options:
--kmer arg (=21) Minimal kmer length for assembly [integer]
--min_count arg Minimal count for kmers retained for comparing
alternate choices [integer]
--max_kmer_count arg Minimum acceptable average count for estimating
the maximal kmer length in reads [integer]
--vector_percent arg (=0.05) Count for vectors as a fraction of the read
number (1. disables) [float (0,1]]
--insert_size arg Expected insert size for paired reads (if not
provided, it will be estimated) [integer]
--steps arg (=11) Number of assembly iterations from minimal to
maximal kmer length in reads [integer]
--fraction arg (=0.1) Maximum noise to signal ratio acceptable for
extension [float [0,1)]
--max_snp_len arg (=150) Maximal snp length [integer]
--min_contig arg (=200) Minimal contig length reported in output
[integer]
--allow_snps Allow additional step for snp discovery [flag]
Debugging options:
--force_single_ends Don't use paired-end information [flag]
--seeds arg Input file with seeds [string]
--all arg Output fasta for each iteration [string]
--dbg_out arg Output kmer file [string]
--hist arg File for histogram [string]
--connected_reads arg File for connected paired reads [string]
実行方法
ペアエンドのfastqを指定する。
skesa --fastq pair1.fq,pair2.fq --cores 8 --memory 48 > skesa_assembly.fa
SRAのシーケンシングデータを直接アセンブリする。
skesa --sra_run SRR1960353 --cores 4 --memory 48 > SRR1960353.skesa.fa
* --sra_runのフラグはビルド時に -f Makefile.nongsを指定していると存在しない。
２つ以上のSRAを指定
skesa --sra_run SRR1695624 --sra_run SRR1745628 --cores 4 --memory 48 > SAMN03218571.skesa.fa
他の例はGIthubのREADMEを参考にして下さい。
簡単なベンチマーク
他のスモールゲノム用アセンブラとごくごく簡単に比較してみる。シーケンシングデータはbacteriaのアセンブリコンペティション GAGE-B（リンク）で使われたRhodobacter sphaeroidesのMiSeqシーケンシングデータを使う( Miseqのペアエンドデータ: SRA link）。MiSeqのリードならオーバーラップが十分あり、SKESAのアルゴリズムが活かせると考えた。SKESAやMEGAHITの特徴であるランタイムの短さはここでは調べない（*1）。
1、fastqの準備。
#1 GATGE-Bよりダウンロード
mkdir test && cd test/
wget https://ccb.jhu.edu/gage_b/datasets/R_sphaeroides_MiSeq.tar.gz
tar -zxvf R_sphaeroides_MiSeq.tar.gz
mv R_sphaeroides_MiSeq/raw/insert_540* .
#2 reenameと圧縮
mv insert_540_1__cov100x.fastq pair1.fq && gzip pair1.fq
mv insert_540_2__cov100x.fastq pair2.fq && gzip pair2.fq
#リファレンス
wget ftp://ftp.ncbi.nih.gov/genomes/Bacteria/Rhodobacter_sphaeroides_2_4_1_uid57653
#ftpのlinkが消えていたのでNCBI Refseqからダウンロードした
#===================================================================
#GAGE-Bでなく、例えばSRAをデータ "SRR1703350" をダウンロードするなら
fastq-dump --split-spot --skip-technical --split-files SRR1703350
#SRAからダウンロードした場合、BBsketchを使いベストマッチするゲノムを推定（リンク）
sendsketch.sh in=SRR1703350_1.fq reads=100k
#ANI近似値が最大のゲノムをダウンロードする
pair1.fq.gzとpair2.fq.gzができる。
２、de novo assembly
#1 skesa
skesa --fastq pair1.fq.gz,pair2.fq.gz --cores 12 --memory 32 > skesa.fa
#2 megahit(紹介)
megahit -1 pair1.fq.gz -2 pair2.fq.gz -o megahit_out --k-min 21 --k-max 141 --k-step 12 -t 12
#3 spades1(紹介)
spades.py -k 21,31,41,51,61,71,81,91,101,111,127 \
-1 pair1.fq.gz -2 pair2.fq.gz \
-t 12 --careful -o spades1
#4 spades2
#spades v3.12で追加された--mergedのリードを組み込む。
#merged.fqの作成
bbmerge.sh in1=pair1.fq.gz in2=pair2.fq.gz \
out=merged.fq outu1=unmerged1.fq outu2=unmerged2.fq \
ihist=ihist.txt
#圧縮
gzip merged.fq
#assembly
spades.py -k 21,31,41,51,61,71,81,91,101,111,127 \
-1 pair1.fq.gz -2 pair2.fq.gz --merged merged.fq.gz \
-t 12 --careful -o spades2
#5 unicycler(紹介)
unicycler -1 pair1.fq.gz -2 pair2.fq.gz -l merged.fq.gz \
-o unicycler -t 12
#6 stride(紹介)
stride all pair1.fq.gz pair2.fq.gz -r 250 -i 550 -t 12
#a7 A5miseq(紹介)
a5_pipeline.pl pair1.fq.gz pair2.fq.gz a5miseq_out_dir
spadesは、v3.12からmerged.fqを使いパフォーマンスを上げるオプションがついたので、"--merged merged.fq"の有無による差も比較した。
３、Quastによる評価（リンク）
#gene情報もあるが、ここでは無しで実行する。
quast skesa.fa megahit.fasta spades1.fasta spades2.fasta unicycler.fasta StriDe.fasta a5output.fasta\
-R GCF_000012905.2_ASM1290v2_genomic.fna \
-1 pair1.fq.gz -2 pair2.fq.gz \
-o quast_output \
-t 12
report.html
SKESAは比較したアセンブラの中で最もキメラアセンブリ、ミスマッチ、indelsが少なかった。
連続性に関しては、NG50がspades やunicyclerの約半分だった。
icarus.html
全体を眺めていくと、下の図のように、SPAdesでアセンブリできていない領域がSKESAではアセンブリできていることもあった。
引用
SKESA: strategic k-mer extension for scrupulous assemblies
Souvorov A, Agarwala R, Lipman DJ
Genome Biol. 2018 Oct 4;19(1):153
*1
主観では、spadesより4倍~5倍以上短い時間で終わる。
