macでインフォマティクス

macでインフォマティクス

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

MaSuRCA アセンブラ

8/28,29 dockerコマンド等、分かりにくい部分を修正

 

  2001年にヒトゲノムのドラフトバージョンが作成された後、800bpを超えるリード長を有する第1世代(すなわちSanger)シーケンス技術を使用して、多くのスモールゲノムおよびラージゲノムのシーケンスが決定された。より最近では、様々な種類のsecond-generation sequencing(SGS)技術が、50〜400bpのリード長で出現している(論文が執筆された当時)。近年、着実に改良されてきた新しいアセンブリ方法が、ショートリードアセンブリの課題に対応して開発されてきた。しかし、この進歩にもかかわらず、ゲノムを決定する問題は、解決から遠い。事実上、現在公開されているすべてのアセンブリは、さまざまなレベルの品質を持つ「下書き」ゲノムであり、下流の分析にこれらのゲノムに頼っている科学者にとって重大な問題となる多くのギャップとアセンブリエラーが含まれている。この記事では、ゲノムアセンブリの新しいアプローチによって促進されるゲノムのアセンブリ方法を報告する。最初に、全ゲノムショットガンシーケンシングデータのアセンブリに使用された2つの一般的なアプローチについて簡単に説明する。

 Overlap–layout–consensus (OLC) assembly。簡単に述べると、OLCパラダイムは、最初に、配列の類似性を用いて重複を決定することによって、リード間のすべてのペアワイズ重複を計算しようと試みる。次に、OLCアルゴリズムはレイアウトを作成する。これはすべての重複するリードのアライメントである。このレイアウトから、アルゴリズムはマルチアライメントアラインメントを列ごとにスキャンすることによってコンセンサス配列を抽出する。 Celera Assembler(論文より Miller et al、2008; Myers et al、2000)、PCAP(Huang、2003)、Arachne(Batzoglou et al、2002)およびPhusion(Mullikin and Ning、2003)を含むSangerシーケンシングデータのためのほとんどのアセンブラ)は、OLCのアプローチに基づいている。OLCアプローチの主な利点の2つは、リード長に関する柔軟性とシーケンシングエラーに対する堅牢性である。見かけ上の重なりが現実である(そしてリピートが誘発されない)可能性を高めるために、OLCアルゴリズムは、典型的には、ある最小長を超えることを要求する。 Celera Assemblerは40 bp以上のオーバーラップを必要とするため、オーバーラップ領域でのエラー率(1〜2%)が小さくなる。

 SGSのデノボ・アセンブリ・プロジェクトでは、サンガーシークエンス・プロジェクトの100倍のリードを生成する。例えばオリジナルのヒトゲノムプロジェクト(Lander et al., 2001; Venter et al., 2001)とマウスゲノムプロジェクト(Mouse Genome Sequencing Consortium et al、 2002)は、それぞれ3500万のリード生成した。一方、最近のヒトゲノムシーケンス(Li et al., 2010) では、3-4億リード生成した。 de Bruijn graphアプローチは、ペアワイズ重複計算を完全に回避する。これがSGSアセンブリの主要な方法になった理由の1つである。

The de Bruijn graph approach。 de BruijnグラフアセンブリアルゴリズムはPevzner et alが先鞭をつけた手法で (Idury and Waterman, 1995; Pevzner, 1989)、オイラーアセンブラで初めて実装された(Pevzner et al、2001)。オイラーはサンガーシーケンスリードのために設計されているが、SGSデータをアセンブリするためのプログラムや、特にイルミナのデータ用のプログラムでも、同じ枠組みが採用されている。 de Bruijn戦略を使用する最近開発されたアセンブラにはAllpaths-LG(Gnerre et al、2010)、SOAPdenovo(Li et al、2008)、Velvet(Zerbino and Birney、2008)、EULER-SR(Chaisson and Pevzner、2008) )およびABySS(Simpsonら、2009)がある。

 このアプローチは、以下のように、シーケンスデータからde Bruijn graphを作成することから始まる。固定値kに対して、各リードから長さk(a k-mer)の各部分文字列は、ノードAとノードBとを接続するgraphの有向枝に割り当てられる。ノードAとBは、それぞれの最初と最後のk-1個のヌクレオチドに対応する。正式にはオイラーパスとして知られているすべてのエッジを正確に1回訪れるgraph内のパスは、リードのドラフトアセンブリを形成している。実際には、これらのグラフは多くの交差するサイクルと多くの代替のオイラーパスを持つ複雑なものであるため、graphを作成することは、優れたドラフトアセンブリを作成する最初の小さなステップに過ぎない。アセンブリを完了するには、graphにメイトペア情報を組み込み、リピートシーケンスによって作成された複雑なサイクルを解消する必要がある。 k-mersはリードよりも短いので、graphにはリードよりも少ない情報しか含まれていないため、後でgraphの曖昧さをなくすためにリードを保持する必要がある。

 このアプローチの主な利点は計算効率である。膨大な数のオーバーラップが計算されないという事実から利益を得ている。主な欠点は、グラフ中のk-mer隣接情報の損失およびデータのエラーによって引き起こされる偽のブランチングである。

このジャーナルでは、著者らが super-readsと呼ぶものを生成する、ショートリードデータアセンブリの第3のパラダイムを提案する。目的は、元のリードに存在するすべてのシーケンス情報を含む super-readsのセットを作成することである。理想的なエラーのない場合については、下記の定理を参照。

  Super-readsの基本的な概念は、extension(以後、伸長)が一意である限り、各元のリードを前後にベースごとに伸長することである。この概念は以下のように説明することができる。効率的なハッシュテーブルを使用してk-merカウントルックアップテーブルを作成し、各k-merが何回発生するかを素早く判断する。リードの終わりに見つかったk-merがある場合、ゲノムのシーケンスの次のk-merになる可能性のある4つのk-merがある。これらは、最後にA、C、GまたはTを付けることによって形成された文字列ですk-1個の塩基が読み取られる。著者らのアルゴリズムは、これらのk-mersのどれがテーブル内に現れるかを調べる。 4つの可能なk-merのうちの1つのみが出現する場合、その読みは固有のk-merを有し、その塩基をリードに付加する。リードが1通りに伸長されなくなるまで続ける。すなわち、複数の可能性が存在するか、または行き止まりに達するまで。この伸長は、リードの3 '末端および5'末端の両方で行う。新しい長い文字列は、super-readsと呼ばれる。多くのリードは、論文図1(リンク)に示すように同じsuper-readsに伸長される。たとえば、2つの同一でないリピートまたは2つの異なるハプロタイプに由来している場合、これらは別個のsuper-readsを生成するだろう。もちろん、super-readsはde Bruijn graphを使って簡単に計算できる。要点は、super-readsが作成されると、super-readsと接続するメイトペアと一緒に、de brujin graphを一括して置き換えることである。ペア情報は、以下に説明するOLCアセンブリ工程を用いて行われる。スーパーリードデータ計算の2つの最も重要な特性は、以下の通りである。

  • 各元のリードはsuper-readsに含まれているため(情報が失われていない)。
  • オリジナルリードの多くは同じsuper-readsになるため、super-readsを使用するとデータセットが大幅に削減される。

 

(複数段落省略)

言い換えれば、super-readsのセットは、リード内のすべての情報を含み、エラーを導入していない。この証拠はsuper-readsの構築に従っている。各リードはsuper-readsに含まれるため、データは失われない(一部略)。MaSuRCAアセンブラは、super-readsからコンティグとスキャホールドを作成するため、CABOGアセンブラアセンブリ技術を利用している。

 

 

ジョンズ ホプキンス大学 Center for Computational Biology(CCB)のツール一覧

http://ccb.jhu.edu/software.shtml

MaSuRCA (Maryland Super-Read Celera Assembler) ページ

http://masurca.blogspot.com

Slideshare

 

発表から時間が経ってますが、ハイパフォーマンスなアセンブリツールの1つなのでまとめておきます。

 

インストール

cent os6とubuntu18.04のdockerコンテナ内でテストした。

本体 Github

https://github.com/alekseyzimin/masurca

ビルドに失敗したので、dockerで導入する。複数イメージがアップロードされているが、sauloさんのイメージを使用した(https://hub.docker.com/r/sauloal/masurca/)。

docker pull sauloal/masurca

#イメージのhomeとカレントディレクトリを共有ディレクトリにして立ち上げる。メモリは30gとした(*2)
docker run --rm -m "30g" -itv ${PWD}/:/home/ sauloal/masurca

#移動
cd /MaSuRCA-2.2.1/bin/

./masurca 

# ./masurca 

USAGE: ./masurca <config_file> at ./masurca line 161.

root@bba6c6f0a8ff:/MaSuRCA-2.2.1/bin# ./masurca -h

Create the assembly script from a MaSuRCA configuration file. An

sample configuration file can be generated with the -g switch. The

assembly script will run the assembly proper.

 

Options:

 -o, --output              Assembly script (assemble.sh)

 -g, --generate            Generate example configuration file

 -p, --path                Prepend to PATH in assembly script

 -l, --ld-library-path     Prepend to LD_LIBRARY_PATH in assembly script

     --skip-path-checking  Skip checking availability of other executables

 -h, --help                This message

 

 

ラン

ランするにはconfigファイルが必要になる。configに記載することで、複数のライブラリを効率的に管理し、まとめてアセンブリできるようになっている。ここでは1組のペアエンドfastqだけランするとして、以下のように修正した。

 

configのマスターシートをコピーして編集する(*1)。

mkdir /MaSuRCA-2.2.1/test/
cd /MaSuRCA-2.2.1/test/
cp /MaSuRCA-2.2.1/sr_config_example.txt config.txt
vi /MaSuRCA-2.2.1/config.txt

cat config.txt

# cat config.txt 

# example configuration file 

 

# DATA is specified as type {PE,JUMP,OTHER} and 5 fields:

# 1)two_letter_prefix 2)mean 3)stdev 4)fastq(.gz)_fwd_reads

# 5)fastq(.gz)_rev_reads. The PE reads are always assumed to be

# innies, i.e. --->.<---, and JUMP are assumed to be outties

# <---.--->. If there are any jump libraries that are innies, such as

# longjump, specify them as JUMP and specify NEGATIVE mean. Reverse reads

# are optional for PE libraries and mandatory for JUMP libraries. Any

# OTHER sequence data (454, Sanger, Ion torrent, etc) must be first

# converted into Celera Assembler compatible .frg files (see

# http://wgs-assembler.sourceforge.com)

DATA

PE= pe 600 60  /MaSuRCA-2.2.1/test/pair1.fq /MaSuRCA-2.2.1/test/pair2.fq

END

 

PARAMETERS

#this is k-mer size for deBruijn graph values between 25 and 101 are supported, auto will compute the optimal size based on the read data and GC content

GRAPH_KMER_SIZE=auto

#set this to 1 for Illumina-only assemblies and to 0 if you have 1x or more long (Sanger, 454) reads, you can also set this to 0 for large data sets with high jumping clone coverage, e.g. >50x

USE_LINKING_MATES=1

#this parameter is useful if you have too many jumping library mates. Typically set it to 60 for bacteria and something large (300) for mammals

LIMIT_JUMP_COVERAGE = 60

#these are the additional parameters to Celera Assembler.  do not worry about performance, number or processors or batch sizes -- these are computed automatically. for mammals do not set cgwErrorRate above 0.15!!!

CA_PARAMETERS = ovlMerSize=30 cgwErrorRate=0.25 ovlMemory=4GB

#minimum count k-mers used in error correction 1 means all k-mers are used.  one can increase to 2 if coverage >100

KMER_COUNT_THRESHOLD = 1

#auto-detected number of cpus to use

NUM_THREADS= 20

#this is mandatory jellyfish hash size

JF_SIZE=100000000

#this specifies if we do (1) or do not (0) want to trim long runs of homopolymers (e.g. GGGGGGGG) from 3' read ends, use it for high GC genomes

DO_HOMOPOLYMER_TRIM=0

END

赤い部分修正したところ。イルミナのショートリードのみの場合。KMER_COUNT_THRESHOLD も1から2に変えてもいいかもしれない。サンプル情報部分は

PE= pe 600 60  /MaSuRCA-2.2.1/test/pair1.fq /MaSuRCA-2.2.1/test/pair2.fq

のように5カラムでデータの情報を記載する。1目のPEはペアエンド、2カラム目の600はインサートサイズ、3カラム目の60はインサートサイズのSD。4−5カラムにペアエンドファイルのパスを記載する。それ以外の行のパラメータは、条件に応じて変更する(コメントアウトされてない行がパラメータ)。

インサートサイズやそのSDをなるべく精度よく求めたければ、seqkit(紹介)などで10%ほどサンプリングしてminimap2(紹介)等でbamを作り、picard-toolsやgoleftなどで計算してください(リンク)。

 

MaSuRCAの実行。

../bin/masurca config.txt -o assemble.sh
#assemble.shができる

#de novoアセンブリ実行
./assemble.sh

エラーコレクションが終わると、そのまま自動でアセンブリが実行される。

[Sun Jun 10 01:33:04 UTC 2018] Gap closing

Gap close success. Output sequence is in CA/10-gapclose/genome.{ctg,scf}.fasta

ランが終わると、カレントに大量のtempoファイルとCAディレクトリができる。CA/10-gapclose/genome.{ctg,scf}.fastaが出力となる。

ホストとの共有ディレクトリにコピー。

cp CA/10-gapclose/genome.*fasta /home/

説明しやすいのでコンテナ中で作業しましたが、sauloさんが記載されているように、コンテナの外からもジョブは投げれます。

 

性能については、GAGE-Bやレビューを参考にしてください。GAGEは論文で用いられたハイカバレッジのシーケンスデータもダウンロードしてテストできるようになっています。

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3702249/

 

 

OLCとde brujin graphのハイブリッドde novoアセンブリツールも出てきています。

イカバレッジのバクテリアのde novoアセンブリツール。

Unicycler

引用

The MaSuRCA genome assembler.

Zimin AV, Marçais G, Puiu D, Roberts M, Salzberg SL, Yorke JA.

Bioinformatics. 2013 Nov 1;29(21):2669-77.

 

GIthubへの移行

http://masurca.blogspot.com/2018/01/masurca-is-now-on-github.html

 

 

参考

*1

emacsに慣れているなら、その場でインストールして使う。

apt-get update && apt-get install emacs

 

*2

何度もランするなら--rmは消してください。毎回configを作るは面倒です。