macでインフォマティクス

macでインフォマティクス

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

(vertebrateなどの) ラージゲノムをアセンブルするためのパイプライン CSA

 

 脊椎動物ラージゲノムの全ゲノムショットガン(WGS)アセンブリは、過去20年間のバイオインフォマティクス研究の重要なテーマだが、脊椎動物の大型ゲノムについては、単一のバイオインフォマティクスツールを用いて完全にアセンブリされた染色体を得ることはまだ実現されていない。De novoゲノムプロジェクトでは、ショートリードからロングリードへの置き換えが進んでいるが、脊椎動物の染色体レベルでのアセンブリには、特に最先端のゲノムマップやscaffoldsデータが利用できないプロジェクトでは、バイオインフォマティクスの専門知識が必要となる。

 今日、ほとんどの脊椎動物ゲノムは、ノイズの多いロングリードを用いてアセンブルすることができ、その結果、コンティグN50として測定されるアセンブルコンティグ性はショートリードシーケンスで得られた結果を100倍以上上回ることがある。今日のノイズの多いロングリードアセンブリのコンティグN50は、数年前に得られた高品質のショートリードゲノムアセンブリーのscaffolds N50に近い長さに達している。それでも、現在のアセンブリツールは、その祖先から利益を得ることができる。これまでのところ、ほとんどのツールはコンティグのみを生成し、これらのコンティグをscaffoldsに整列させるための付加的な情報を組み込んでいないが、付加的な情報によってさらなるギャップクローズと染色体レベルのアセンブリが可能になる。

 ゲノムプロジェクトの最終目標である染色体レベルのゲノムアセンブリには、やはり追加のスキャフォールディングやマッピングデータ(Hi-Cやオプティカルマッピング、高密度遺伝的連鎖地図)が必要であり、その結果、シーケンシングプロジェクトに多大な人的、時間的、資金的リソースを追加することになる。多くの、特に希少種の場合、de novoゲノムシークエンシングのためのDNAはアーカイブ組織(例えば、冷凍、エタノール固定、または他の保存媒体に保存されたもの)から得られるため、生きた細胞を必要とすることが多いHi-Cを適用することができず、その結果、マッピングパネルを確立することはほとんどできない。このような場合には、進化的に関連するゲノム間のシンテニーおよび遺伝子オーダー解析が、ゲノムアセンブリプロセスを改善するための唯一の選択肢であるかもしれない。  

 現在(2019年8月)、約71,000種が科学的に記述されている脊椎動物はすべて、2つの先祖代々の全ゲノム重複(WGD)を経験しており、約38,000種の四足類が現存している一方で、約33,000種のテレポストのほとんどは3回目のWGDを経験している[ref.21, 22]。脊椎動物のゲノムサイズは、欠失、サイレンシング、仮性化、サブ機能化、新機能化の進化に影響を与えた異なる祖先のWGD「基質」を超えて、大きく異なる(典型的なサイズは0.4~4Gbの範囲で、この点では両生類の巨大なゲノムサイズは強力な例外となる)。C値(ハプロイドゲノムサイズ)3.3~57pg [ref.23])。形態学的複雑さや遺伝子数とは無関係だが、ゲノムサイズの違いは様々な反復エレメントやその他のノンコーディング DNA の量に起因しており、 脊椎動物のゲノムの最大 98%を構成している[ref.24]。

 

(一部省略)  

 実際、脊椎動物の近縁種では染色体構造(シンテニーや遺伝子順序)が高度に保存されているような進化的関係があり[ref.42]、いくつかの分類群では複数の分類群間であっても[ref.43]、比較ゲノミクスによる染色体スケールのアセンブリを低コストで近似することが可能である[ref.44-47]。ショートリードアプリケーションの成功例としては、ゲノムアセンブリーパイプラインIMAP [ref.48]がある。もちろん、このようなアプローチを行うためには、1つ以上の適切な高品質リファレンスゲノムが必要であり問題となっている。

 ここでは、「染色体スケールアセンブラ」(CSA)と呼ぶ新しいバイオインフォマティクスパイプラインを紹介する。CSAは、オプティカルマッピング、Hi-C、または10X Genomicsから得られた発散したリファレンスゲノムやscaffolds間の比較をde novoアセンブリプロセスに統合することで、現在のロングリードアセンブラの限界を克服している。また、ロングリードゲノムアセンブリ、全ゲノムアラインメント、リファレンスアシスト染色体アセンブリのための計算効率の高いツールを反復的に実行している。脊椎動物の小型ゲノム(魚類、鳥類)の染色体レベルのアセンブリを、低コストの計算機(1,000~2,000ドル、Intel i7、128GB RAM)を用いて、ロングリードデータと発散したリファレンスゲノム(発散時間約6,500万年前[Mya])を入力として使用するだけで、12時間以内に作成できることを示した。ヒトのような大きな哺乳類のゲノムは、サーバーマシン(Intel Xeon、1 TB RAM)で16時間以内にアセンブリできる。入力データの種類やカバレッジにもよるが、コンティグN50の長さを初期から最終的なコンティグアセンブリまで最大4.5倍に改善することができる。

 CSAの最初のステップ(論文図1)は、ノイズの多いロングリードデータ(Pacific Biosciences [PacBio]またはOxford Nanoporeデータのいずれか)のde novoアセンブリである。これにはWTDBG2(バージョン2.2、2018年12月11日)アセンブラを使用しているが、それは、それが今日までの最も計算効率の高いde novoゲノムアセンブラであり、RuanとLi [ref.14]によると、その最も近い競合の2〜17倍の速さである。CSAは、パラメータを少し変えた2つのWTDBGアセンブリを実行し、両方のアセンブリー間の不一致部分でコンティグを分割することで、まれに発生するミスアセンブリを排除している。さらに、部分的にしかWTDBG2アセンブリマッピングできないロングリード(リード長の10%未満)を再アセンブルする。このステップにより、現在のWTDBG2プライマリーアセンブリでは欠落しているゲノム配列(大きなコンティグ)を最大1~2%回復できる。最終的なアセンブリキュレーションは、アセンブリへのロングリードの再マッピング(minimap2)[ref.51]で、カバレッジゼロの領域でコンティグを分割する。必要に応じてWTDBG2アセンブリを省略し、別のゲノムアセンブリーツールからのコンティグファイルをCSAで使用することができる;したがって、CSAは既存のアセンブリを更新するためにも使用することができる。

 ステップ2では、結果として得られたキュレーションされたコンティグを、LASTアライナー[ref.52]により、1つ以上のリファレンスにマッピングする。これらのリファレンスは、様々な方法(例えば、10X Genomics、オプティカルマッピング、またはHi-C)を用いて構築された同一種からのscaffoldsであってもよい。CSAの優れた特徴は、発散したリファレンスゲノムとドラフトゲノムであっても適切な入力となることである。RAGOUT [ref.53]では、LASTアラインメントを用いて、(ステップ1で得られた)キュレーションされたコンティグをscaffoldsにアラインさせる。

 ステップ3では、ノイズの多いロングリードはすべてMINIMAP2によってscaffoldsにマッピングされる。scaffoldsのギャップやコンティグエンドの周りに20kbのウィンドウでマップされたリードは、別個のfastaファイルに抽出される。これらのファイルはWTDBG2アセンブラに送られ、並列にローカル再アセンブルされる。結果として得られた各ギャップ/コンティグエンドのローカル再アセンブルは、その後、WTDBG2のリード長の制限(バージョン2.2では256kb、バージョン2.4では制限はないが、テストではパフォーマンスが低下している)を満たすように、オーバーラップのある疑似リードに分割された一次WTDBG2コンティグを用いてアセンブルされる。WTDBG2では、より高い精度(コンセンサス精度98~99%)であらかじめアセンブリされたリードをアセンブリするようになったため、より厳しいパラメータが設定されている。この反復的なアセンブリステップは、以下の異なる試験で示されるように、典型的にはN50コンティグサイズを2倍にすることができる。改良されたコンティグの先行スキャフォールドへのアラインメントは、少数のスキャフォールド内およびスキャフォールド間のミスアセンブリを除去するために使用される。

 ステップ4では、再び、ステップ2で使用されたリファレンスに対し、LASTアライナーで改良されたコンティグをマッピングし、コンティグをscaffoldsにアラインさせるためにRAGOUTを実行する。最後に、隣り合うコンティグエンドが重なっているいくつかのギャップは、LASTによって特定され、閉じられる。

 CSAパイプラインを構築するために選択されたツールは、主に発散したリファレンスゲノムを使用する際の性能と感度に基づいて選択されている。比較的単純な手順では、支持率の低い領域(例えば、2つのde novoアセンブリを比較する際に連続的にカバーされていない領域やリードカバレッジの低い領域)で一次de novoアセンブリを分割し、非常に大きなゲノムを扱う場合でも比較的高速に処理することができる。アセンブリのリコンサイルについては、SMSC [ref.54]やBIGMAC [ref.55]など、より洗練されたツールが開発されているが、これらのツールはラージゲノムではほとんどテストされておらず、スモールゲノム(細菌、酵母)での公表されているベンチマークによると、これらのモジュールはCSAパイプライン全体よりも時間がかかる可能性がある。それにもかかわらず、リードリマッピングによるミスアッセンブリーの検出速度が向上すれば、将来のCSAのバージョンではアセンブリの品質がさらに向上する可能性がある。

 MeDuSa [ref.56]やRACA [ref.44]のような同様のツールは、計算量が多い(RACAのLASTZ)か感度が低いアライナーを使用しており、発散したリファレンスゲノム上ではうまく機能しない(MeDuSaのMUMMER [ref.57])。RAGOUT1 [ref.53]は、RAGOUT2 [ref.45]よりも、わずかに良好な染色体アセンブリをもたらしたために好まれた。RAGOUTはギャップを閉じることもできるが、この機能はショートリードアセンブリからのコンティグの短い、エラーのないオーバーラップのために設計されたものであり、ロングリードアセンブリでは効率的には働かない。そこで、ギャップのローカル再アセンブリと最終的なコンティグのステッチのための独自のソリューションを実装した。原理的には、ギャップクローズはPBJelly [ref.58] や LR_gapcloser [ref.59] などのツールでも可能であるが、これらのツールではクローズ可能なギャップの中にクローズされないものがあることがわかった。

 現在のCSAのバージョンには、コンセンサスポリッシュの方法は含まれていないが、高品質なコンセンサス配列を得るためには、ロングリードを用いたコンセンサスポリッシュを1回、ショートリードを用いたコンセンサスポリッシュを2回繰り返してから、配列アノテーションを行う必要がある。この点では、MEDAKA(Oxford Nanopore Technologies [ONT] [ref.60])や PILON を使用した経験がある(ベンチマークシナリオ 3)。最近では、MEDAKA の代わりに FLYE-POLISH ref.[12] を使用しているが、これは SMRT (single molecule real time; PacBio) と ONT の両方のデータを使用できるからである。

 

 

f:id:kazumaxneo:20200607154631p:plain

CSA pipeline. Githubより転載。

 

 

CANU、FALCON、FLYE、そしてSHASTAのような広く使われているゲノムアセンブリツールはすべて、Q40以下のエラーレートを得るために必要なポリッシングアルゴリズムを含んでいない(これはフレームシフトエラーが過剰になることなくアノテーションを行うために必要とされる)。現在一般的に行われているように、コンセンサスポリッシュの選択はユーザに任せている。これは、異なるロングリード技術(SMRTやONT)に適したポリッシュパイプラインを選択することは複雑なテーマであり(例えば、異なるシーケンシングライブラリが関与し、また、シーケンシングする種の特性にも依存する)、著者らの意見ではヒトの意思決定を必要とするからである。著者らは、FLYEポリッシャとロングリードデータによるポリッシュの1反復を行い、続いてPILONとイルミナショートリードデータを使用して2回ポリッシュで良い経験をしてきた。ロングリードポリッシャ(MEDAKAやFLYEなど)の中には、"n "文字を嫌うものもあるので、ポリッシュ前にCSAのscaffoldsをコンティグに分割し、ポリッシュ後にscaffoldsを再構築する必要がある。

インストール

ubuntu18.04LTSでテストした。

依存

  • Miniconda, wtdbg2.2, RAGOUT v1.0, Minimap2, last aligner v941, bedtools v2.27.1, samtools v1.7, pigz, GNU parallel, seqtk, seqkit, hgFakeAgp

Github

#ここではpython3.6の仮想環境を作ってテストする。
conda create -n CSA python=3.6
conda activate CSA

#clone repository
git clone https://github.com/HMPNK/CSA2.6.git
cd CSA2.6/CSA2.6/INSTALL/
bash INSTALL.bash

./CSA2.6c.pl

# ./CSA2.6c.pl 

 

CSA version 2.6c (Chromosome Scale Assembler, employing synteny with diverged reference genomes)

Author: Heiner Kuhl, Phd (kuhl@igb-berlin.de)

 

THIS SCRIPT RUNS FOUR STEP CSA2 PIPELINE

 

INPUT DATA: pacbio or oxford nanopore long reads (provided as fa.gz)

            suitable references (scaffolded assembly, 10X Genomics assembly, diverged reference genomes)

 

      Options:

 

              -d directory to run CSA in (will be created)

 

              -r input long reads as fa.gz (required!)

 

              -C input genome assembly contigs (overrides initial WTDBG assembly step)

 

              -g reference sequences (suitable backbone assembly or reference genomes,

                                     ordered by evolutionary distance (most similar first), 

                                     as comma separated list. Input format fasta or fasta.gz)

 

              -o output-prefix

 

              -t number threads (default=2, use multiples of 2)

 

              -k kmer wtdbg assembler (default=21, use 21 for high coverage data)

 

              -s step wtdbg assembler (default=4, use 4 for high coverage data and/or to save RAM)

 

              -e minimum edge coverage wtdbg assembler (default=4, decrease for low coverage data) 

 

              -m Min similarity in WTDBG (default=0.05)

 

              -x off/on (default on). Set "-x on" to further improve chromosomal assembly on the cost of more structural missassemblies.

                 Use "-x on" if you a have a closely related high quality reference genome with similar karyotype (family or genus level, same haploid chr no.)

                 If "-x off" a second assembly using guessed edges in the last assembly step will be output anyway since version 2.5 ("xyz.final.B.fa")

 

              -p 0,1,2 (0 = wtdbg-cons (default), 1 = wtpoa-cons; consensus calculation in step 1, 2 = wtdbg-cons -S 0

                        wtpoa-cns is slower but a bit more accurate; Parameter -S 0 is more stable on very long reads (e.g. ONT Ultra Long Reads))

 

              -l special parameters for wtdbg2 (default= -l "-L 5000 -A"; e.g. use -l "-L 70000 --aln-min-length 25000 --keep-multiple-alignment-parts 1 -A" for ultra long reads)

 

This script generates a bash script for running the pipeline! Write script to file and run by: nohup bash <script> & !

 

 

テストラン

1、データのダウンロード

mkdir CSA-TEST
cd CSA-TEST/
wget https://nimbus.igb-berlin.de/index.php/s/njM2jqplusn17OZ/download
mv download SRR6476833.fa.gz
wget https://nimbus.igb-berlin.de/index.php/s/4VekUKms8tdL4V4/download
mv download sacCer.fa.gz

ls -lh

# ls -lh

total 821M

-rw-r--r-- 1 root root 1.1K Jun  7 09:14 RUN-CSA-TEST.bash

-rw-r--r-- 1 root root 3.7M Jun  7 09:13 sacCer.fa.gz

 

2、 configファイル作成

#CREATE CSA-PIPELINE SCRIPT
../CSA2.6c.pl -r SRR6476833.fa.gz -g sacCer.fa.gz -t 4 -o SC_CSA -d SC_CSA > RUN-CSA-TEST.bash

cat RUN-CSA-TEST.bash

# cat RUN-CSA-TEST.bash 

#!/usr/bash

set -e

set -o pipefail

 

 

#RUN CSA2.6 assembly pipeline:

 

 

 

export PATH=/root/CSA2.6/INSTALL/../bin/RAGOUT_V1.0/lib:$PATH

export PYTHONPATH=/root/CSA2.6/INSTALL/../bin/RAGOUT_V1.0/lib

 

mkdir SC_CSA

cd SC_CSA

 

mkdir 01_WTDBG

cd 01_WTDBG

/root/CSA2.6/INSTALL/../script/01_ASSEMBLY.pl -r ../../SRR6476833.fa.gz -o SC_CSA -t 4 -k 21 -s 4 -e 4 -m 0.05 -p 0 -l "-L 5000 -A" > STEP_01.bash

time bash STEP_01.bash

cd ..

 

 

mkdir 02_ORDERCTG

cd 02_ORDERCTG

/root/CSA2.6/INSTALL/../script/02_ORDERCONTIGS.pl -c ../01_WTDBG/SC_CSA.step1.fa -g ../../sacCer.fa.gz -o SC_CSA -t 4 -x on > STEP_02.bash

time bash STEP_02.bash

cd ..

 

 

mkdir 03_REASSEMBLE

cd 03_REASSEMBLE

/root/CSA2.6/INSTALL/../script/03_REASSEMBLE.pl -r ../../SRR6476833.fa.gz -a ../02_ORDERCTG/SC_CSA.step2.fa -t 4 -o SC_CSA > STEP_03.bash

time bash STEP_03.bash

cd ..

 

mkdir 04_ORDERCTG

cd 04_ORDERCTG

/root/CSA2.6/INSTALL/../script/04_ORDERCONTIGS.pl -c ../03_REASSEMBLE/SC_CSA.step3.fa -g ../../sacCer.fa.gz -o SC_CSA -t 4 -x on > STEP_04.bash

time bash STEP_04.bash

cd ..

 

 

#END CSA2.6 assembly pipeline

 

3、ラン

bash RUN-CSA-TEST.bash

出力

>ls -lh SC_CSA/

T# ls -lh SC_CSA/

total 13M

drwxr-xr-x 2 root root 4.0K Jun  7 09:44 01_WTDBG

drwxr-xr-x 3 root root 4.0K Jun  7 09:44 02_ORDERCTG

drwxr-xr-x 2 root root  24K Jun  7 09:50 03_REASSEMBLE

drwxr-xr-x 3 root root 4.0K Jun  7 09:50 04_ORDERCTG

-rw-r--r-- 1 root root  13M Jun  7 09:50 SC_CSA.final.A.fa

-rw-r--r-- 1 root root 2.3K Jun  7 09:50 SC_CSA.final.A.fa.fai

  

 論文にヒトを含むいくつかの動物ゲノムのアセンブリコマンドが掲載されています。確認して下さい。

引用

CSA: A high-throughput chromosome-scale assembly pipeline for vertebrate genomes
Heiner Kuhl, Ling Li, Sven Wuertz, Matthias Stöck, Xu-Fang Liang, Christophe Klopp

GigaScience, Volume 9, Issue 5, May 2020

 

関連