2019 7/26 追記
生物のDNA配列の正確な決定、すなわち、DNA分子中のヌクレオチドA、C、GおよびTの正確な順序を確立することは、生物学における基本的かつ挑戦的な問題である。本質的にこのプロセスは2つのステップから成っている：（1）ケミカルプロセスによってDNAをシークエンシングし、多数のリードを生じさせ、（2）リードを処理して完全なDNA配列を再構築するゲノムアセンブリを行う。すべてのシーケンシング技術ではシーケンスにエラーが含まれ、エラープロファイルはプラットフォーム間で大きく異なる。第2世代のシーケンサーと第3世代のシーケンサーとの間には明確な区別があり、後者の方がはるかに高いエラー率ではあるが大幅に改善されたリード長を特徴とする。
第2世代シーケンシングでは、主にイルミナのプラットフォームを検討する。イルミナの技術は、高い精度（< 2％以下のエラー率、主に置換）で短い（100-300ヌクレオチド）シーケンスを高スループットかつ低コストに生成する。 de Bruijnグラフに基づく新しいアルゴリズムは、膨大な量の第2世代シークエンシングデータのアセンブリに効率的に対処するために特別に開発された。 k-merを共有するリード間、すなわち長さkの部分文字列の間で、短いリード間の重複が線形時間内に確立される。しかしながら、de Bruijnグラフの繰り返し分解能は、第2世代データの非常に短いリード長によって著しく妨げられる。
最近、第3世代配列決定技術（Pacific Biosciences、2013; Oxford Nano Technologies、2014）が出現し始めた。 Pacific Biosciences社のSMRTシークエンシングは、高いエラー率（論文執筆次点 最大15％、主に挿入および欠失およびより少ない程度の置換）であるにもかかわらず、はるかに長いリード（平均> 5000ヌクレオチド）をもたらす。 この高いエラー率にもかかわらず、非常に高いコンセンサス精度が達成される。なぜなら、エラーは読み取りにわたって均一に分散されるからである。 カバレッジが十分に高く、読み取り間のオーバーラップが正しく確立されている場合、この一様なエラー分布は非常に正確なコンセンサス呼び出しを可能にする。 これらの重なりを計算することは、エラー率が高いと誤ったk-mersが過剰になるため、de Bruijn graphによって効率的に達成することはできない。 したがって、第3世代のリード間のペアワイズアライメントを計算するための他の効率的な方法が開発されている[論文より ref.1,2]。
シーケンシングリードの処理では、通常、相互に潜在的なオーバーラップを確立するためにリードをアライメントさせるか、リファレンスゲノムにマッピングする。リードのエラーは、これらのアライメントにノイズを導入し、対応するエラーのないリードよりも弱いアライメントにつながる。より低い評価のアラインメントは、その後の分析のために廃棄され、重大な情報を廃棄する可能性がある。これは、カバレッジの低い領域で低品質のリードを処理する場合に特に問題になる。このシーケンシングノイズを処理するために、エラー訂正を適用する。リードにおけるエラーを訂正することにより、最適なアライメントをより正確に同定し、より適切に定格化することができる。第2世代のリードを訂正するアルゴリズムは3つのタイプに分類されている[ref.4]。 k-merスペクトラムベースの方法[ref.5,6]は、k-merが実際のDNA配列の一部を表すかどうかを決定するためにカバレッジ閾値に依存している[ref.5,6]。サフィックスツリーベースの方法[ref.7,8]は、複数のk値を一度に処理することによってkスペクトラム法を一般化する。最後に、複数の配列アライメントベースの方法[ref.9]は、いくつかの同様のリードを整列させた後にリードを修正する。
第3世代のリードを訂正するために、互いにリードをアライメントさせ、重複するリード間のコンセンサス配列を計算する方法を用いることができる。しかし、第3世代のリードの高精度コンセンサスに基づく補正に必要なカバレッジは、多くのシーケンシングプロジェクトにとって非常に高いファイナンシャルコストにつながる。ハイブリッドエラー訂正方法は代替手段を提供する。目標は、第2世代のリードに含まれるより正確なシーケンス情報を使用して、長い第3世代のリードを修正することである。第3世代データのカバレッジにかかわらず、（比較的安価な）第2世代のデータセットでは、長いリードを修正するのに十分であるという考えがある。これは、低カバレッジの第3世代のデータで十分かもしれないので、シーケンシングのためのファイナンシャルコストの削減をもたらす可能性がある。ハイブリッドエラー訂正方法はまた、長いリード間のペアワイズ比較を回避し、したがって二次計算の複雑さを回避するので、計算上の観点から魅力的であるように見える。第1のタイプのハイブリッドエラー訂正方法LSC [ref.10]、PacBioToCA [ref.11]およびproovread [ref.12]は、ショートリードからロングリードへのマッピングに依存し、次いで、この複数のアライメントからコンセンサス配列を呼び出す。しかしながら、そのような方法は、ショートリードを個別にマッピングし、ショートリードが生じるコンテキストを利用しない。より最近のハイブリッドエラー訂正方法であるLoRDECは、最初にショートリードからde Bruijn graphを構築し、このgraph上の長いリードをマッピングする。ロングリードが揃うgraph内の経路によって暗示されるシーケンスは、訂正されたリードを表す。 de Bruijn graphの使用には、ショートリード間のオーバーラップがロングリードにマッピングする前に確立されるという利点がある。 [ref.13]では、LoRDECが他のエラー訂正方法と同様の精度を達成したが、ランタイムが大幅に改善されていることが示された。 LoRDECは、すべてのシードがグラフのノードに対応するk-merインデックスを使用する。
ここでは第3世代のリードのためのハイブリッドエラー訂正方法であるJabbaを紹介する。 Jabbaでは、第3世代のリードは、第2世代のリードから構築されたde Bruijn graph[ref.14]にマップされ、シード・アンド・エクステンション手法に基づく擬似的なアライメント手法を使用する。graphの結果として得られる経路は、リード訂正を指示する。シードは、個々のリードとgraphのノードとの間の最大完全一致（MEM）である。シードとしてのMEMの使用は、それらがLoRDECで使用されるので、k-mersよりもいくつかの利点を有する。第一に、シードはより長くなり得る。たとえ長いシードがまれにしか発生しないとしても、いくつかのより長いシードで、リードをどのようにしてgraphにアライメントさせるかを概算することができる。これをさらに洗練するために、より短いシードを使用することができる。第2に、拡張添字配列[ref.15]が与えられると、このインデックスを再構築する必要なしに、任意の長さのシードを探すことができる。これは、k-merインデックス（ハッシュテーブルなど）では当てはまらない。アラインメント処理中に異なる値のkを使用する必要がある場合、グラフのさまざまなk-merインデックスを構築する必要がある。最後に、MEMの使用は、de Bruijn graphを構築するためにkの任意の値を使用することを可能にする。第3世代のリードのエラー率が高いことは、シードサイズの最小値を限定する要因であるため、ハイブリッドエラー訂正技術の現状をはるかに上回る。シードサイズとk-merサイズのこのような切り離しは、より大きな値のkを使用してde Bruijn graphを構築することを可能にし、その結果、Bruijnグラフはより複雑になる。 de Bruijn graphのk-merサイズは、第2世代データのエラー率によって制限される。このようにして、graphを構築する前にショートリードを修正し、MEMをシードとして使用することで、de Bruijn graphの大きなk値が可能になり、多くの小さな繰り返しが効果的に解決される。
マニュアル
https://github.com/biointec/jabba/wiki
インストール
cent OS6でテストした。
依存
- CMake 2.6
- GCC 4.7
git clone https://github.com/biointec/jabba.git
cd jabba/
mkdir build
cd build/
cmake ../
make -j
cd src/
> ./jabba
$ ./jabba
Jabba: Hybrid Error Correction.
Jabba v.1.0.0
Copyright 2014, 2015 Giles Miclotte (giles.miclotte@intec.ugent.be)
This is free software; see the source for copying conditions. There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
Usage: Jabba [options] [file_options] file1 [[file_options] file2]...
Corrects sequence reads in file(s)
[options]
-h --help display help page
-i --info display information page
[options arg]
-l --length minimal seed size [default = 20]
-k --dbgk de Bruijn graph k-mer size
-e --essak sparseness factor of the enhance suffix array [default = 1]
-t --threads number of threads [default = available cores]
-p --passes maximal number of passes per read [default = 2]
-m --outputmode short (do not extend the reads) or long (maximally extend reads) [default = short]
[file_options file_name]
-o --output output directory [default = Jabba_output]
-fastq fastq input files
-g --graph graph input file [default = DBGraph.fasta]
examples:
./Jabba --dbgk 31 --graph DBGraph.txt -fastq reads.fastq
./Jabba -o Jabba -l 20 -k 31 -p 2 -e 12 -g DBGraph.txt -fastq reads.fastq
ラン
de brujin graphは3-rdパーティのツールで作る必要がある。ここではBrownieを使う（マニュアル）。
https://github.com/biointec/brownie
sparsehashが必要になる。なければ"sudo yum install sparsehash-devel"で導入。またはここからrpmをダウンロードして"sudo rpm -ivh sparsehash-devel-1.12-3.sdl7.x86_64.rpm"を実行（2018 05実行時のバージョン）。
git clone https://github.com/biointec/brownie.git
cd brownie/
mkdir build
cd build
cmake ..
make install
> brownie -h
$ brownie -h
Usage: brownie command [options] [file_options] file1 [[file_options] file2]...
Corrects sequence reads in file(s)
command
readCorrection Correct input reads [default]
graphConstruction Build de bruijn graph only
graphCorrection Correct constructed de Bruijn graph only
[options]
-h --help display help page
-i --info display information page
-s --singlestranded enable single stranded DNA [default = false]
[options arg]
-k --kmersize kmer size [default = 31]
-t --threads number of threads [default = available cores]
-v --visits maximum number of visited nodes during bubble detection [default = 1000]
-d --depth maximum number of visited nodes during read correction [default = 1000]
-e --essa sparseness factor of the enhanced sparse suffix array [default = 1]
-c --cutoff value to separate true and false nodes based on their coverage [default = calculated based on poisson mixture model]
-p --pathtotmp path to directory to store temporary files [default = current directory]
[file_options]
-o --output output file name [default = inputfile.corr]
examples:
./brownie readCorrection inputA.fastq
./brownie readCorrection -k 29 -t 4 -o outputA.fasta inputA.fasta -o outputB.fasta inputB.fastq
(anaconda2-4.2.0) [uesaka@cyano src]$
de brujin graph作成前にはエラー訂正もしておく必要がある。karectを使う（インストール紹介）。
karect -correct -threads=20 -matchtype=hamming -celltype=haploid -inputfile=./pair1.fastq -inputfile=./pair2.fastq -kmer=9
brownieを使いde brujin graphを構築する。
#ペアエンドはインターレースにしておく。（seqtk）
seqtk mergepe karect_pair1.fastq karect_pair2.fastq > karect_interlace.fastq
mkdir brownie_data
brownie graphCorrection -p brownie_data -k 31 karect_interlace.fastq
brownie_data/にDBGraph.fasta等が出力される。
$ head brownie_data/DBGraph.fasta
>NODE 1 85 1 -730 1 -436
GTGTATACAGAGTAAAATTATCAGACGTCAAGCTTTTTTGTCAATTTGAATTTTGATGATAATTGCGTTGAGATTTGATACGATA
>NODE 2 61 1 -9475 1 -7634
GAGGAATGCAAATTTCAAAGAAAAGAAATGATCGTTGTACAAAGTCTGGTTTATTGCTTGA
>NODE 3 133 1 2887 1 2887
TTCAAATTCTAGTTATATATTTATACATACGATACAACGCAGCTTGGTCAAATATATCATAAGTATGAAAATATATGCAAATCCTGAACTCATGCGCGCGCAAGCGCGCGTAAATTCAAATTCTAGTTATATA
>NODE 4 61 1 8641 1 3729
AAAAAATCTTTTTGGAAACATTCCATGTTGTGATGATTCCCGCGTCGTTTTACACCAGACC
>NODE 5 32 2 4982 3452 1 -3507
AGTACAATTTATTACTTATTTAAGTAACGAGA
最後にJabbaでエラーコレクションを行う。
jabba -o Jabba -l 20 -k 31 -p 2 -e 1 -g DBGraph.fasta -fastq karect_pair1.fastq karect_pair2.fastq -fasta long_reads.fasta
ショートリード、ロングリードそれぞれについて、エラー訂正されたリードのファイルと、エラー訂正されなかったリードのファイルが出力される。
ショートリードは１度目のエラーコレクションツール（karect）とJabbaで２回エラー訂正が行われる。2018年２月のNature Communicationsにpublishされた、シロイヌナズナのゲノムのONTシーケンスの論文（ PCR-free paired-end readsでpolish）（https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5803254/）のデータでテストしてみた。
illumina paired-end fastq: 10GBx2（ERR2173372）
karectの解析でユニークなk-mer は2億から5千万まで減っている。それからjabbaでハイブリッドエラーコレクションを行うことでユニークk-merは120万まで低下している。シーケンスデータの40%ほどしかjabbaでは出力されなかったので、正規化すると差は縮まるものの、シーケンスエラー由来と考えられるユニークk-merは大きく減った。
ONT fastq : 8GB（ERR2173373）
ONTリードは30%ほどしかエラー訂正できなかったが、エラー訂正後は、ユニークk-merが大きく減り、raw fastqでか確認できなかった真のピークらしきk-mer頻度が出てきている。
引用
Jabba: hybrid error correction for long sequencing reads
Giles Miclotte, Mahdi Heydari, Piet Demeester, Stephane Rombauts, Yves Van de Peer, Pieter Audenaert, Jan Fostier
Algorithms Mol Biol. 2016; 11: 10
比較のプレプリントが出ています。
https://www.biorxiv.org/content/biorxiv/early/2019/05/29/519330.full.pdf
