2019 7/26 condaインストール追記、流れ修正
ゲノムのデノボアセンブリは、いわゆるロングリードシークエンシング技術の導入により劇的に恩恵を受けている。 PacbioによるSMRTシーケンシングやOxford Nanopore Technologiesによるナノポアシークエンシングプラットフォームなどのこれらの技術は、典型的には10キロベース読める。これらのリードは、以前は解明できなかったゲノムの反復性または低複雑性領域をカバーすることができる。しかし残念なことに、これらの技術の比較的高いエラー率は、新しい情報および分析上の課題をもたらす。エラーを修正するための効果的かつ効率的な方法が必要である[論文より ref.1-4]。
ロングリードのデノボアセンブリは、ペアワイズオーバーラップステップによって駆動される。これらのロングリードを互いに重複させるのに必要な時間は、リードの数に対して2次的に増加する。 MHAP [ref.5]やMinimap [ref.6]などの新しい方法ではこれを改善することを目指しているが、実際には必要な計算時間とメモリは非常に多くなる。プリアセンブリの修正により、オーバーラップステップで考慮する必要のある分散を減らすことによって、後続のオーバーラップレイアウトが大幅に簡素化される。特に、エラー訂正を受けたロングリードなら、長い同一の伸長を共有する可能性が高くなる。修正がより長く、より正確であれば、より迅速かつ正確に長いリードをアセンブルすることができる。
ロングリードのエラー訂正アルゴリズムは、自己訂正アルゴリズムまたはハイブリッド訂正アルゴリズムのいずれかに大別できる。自己訂正アルゴリズムは、他のロングリードシーケンスのみを使用してロングリードを修正する。 Sprai [ref.7]、LoRMA [ref.8]、HGAP [ref.1]、およびPBcR [ref.3]は、ロングリードを互いに整列させ、コンセンサスシーケンスを生成する。正確なコンセンサスを生成するために、これらの方法は比較的高いカバレッジを必要とする。対照的に、ハイブリッド補正アルゴリズムは、同じサンプル由来のショートリードシーケンシングを使用してロングリードをエラー訂正する。ショートリードシークエンシングではシークエンシングエラーが少なく、シーケンコストも低いため、ヌクレオチドあたりのコストはずっと低くなる。多くのハイブリッドエラー訂正手法は、ショートリードデータを最初にプリアセンブルし、次いでショートリードとロングリードとの間のアラインメントを行う。これらのアプローチは合理的だが、2つの問題が発生する。第1に、ショートリードのサイズを超えるリピートが確実に解決できないという点である。第2に、ショートリードアセンブリとそれに続くペアワイズアラインメント/ショートリードコンティグによるロングリードのオーバラップは、直接的なロングリードのエラー訂正手法よりも著しく遅い。ハイブリッド補正アルゴリズムは、LoRDEC [ref.4]、ECTools [ref.9]、Jabba [ref.10]、CoLoRMap [ref.11]、Nanocorr [ref.12]などがある。 Cerulean [ref.13]、DBG2OLC [ref.14]、hybridSPAdes [ref.15]などの方法では、ロングリードとショートリードのハイブリッドアセンブリを実行するが、ロングリードのシーケンスエラーは明示的に訂正しない。これらのハイブリッド法は実質的により低コストで、より正確で連続したアセンブリをしばしば構築できる。 ECTools [ref.9]およびNanocorr [ref.12]は、同じ基本的な方法論に基づいているが、Pacbioおよびナノポアシーケンスのために設計されている。これらの手法はショートリードとロングリードの間の完全な位置合わせを行うが、微生物ゲノムよりも大きなものを実行するにはかなり長い時間がかかり、それ以上考慮されていない（現在は推奨されていない）。
生物医学的または経済的に重要な、大規模で複雑なゲノムのエラー訂正やアセンブリの場合、可能な限り迅速に、可能な限り少ない計算資源を使用して、できるだけ正確にアセンブリを実行することが重要な課題となる。現在の方法は、しばしば、非常に大きな計算資源を必要とする。アセンブリに適切なパラメータを見つけることはしばしば反復的なプロセスであることを考えると、これらの高い計算コストは障壁になる。
FMLRC（FM index Long Read Corrector）はロングリードのシーケンスエラーを修正するためのハイブリッドメソッドである。 本手法の主な貢献は、ショートリードシーケンシングデータのマルチストリングBurrows-Wheeler Transform（BWT）[ref.16]から構築されたFM indexの適用である。 FM indexは、k-mer頻度をFMLRCがO（k）ステップで検索することを可能にする。 他のデータ構造とは異なり、kの長さはFM indexの構築時には固定されず、実行時に選択される。 結果として、FMLRCはFM indexを使用して、ショートリードシーケンシングデータセットのすべてのde Bruijnグラフ[ref.17]を暗示的に表現する。 次に、これらのde Bruijnグラフを使用して、ショートリードシーケンシングデータでサポートされていないロングリードを修正する。論文中で、FMLRCが既存のハイブリッドエラー訂正方法よりもかなり高品質のエラー訂正シーケンスを効率的に生成することを示している。また、多くのケースで、既存のハイブリッドおよびロングリード単独のデノボアセンブリ方法よりも、より長いアセンブリを生成することも示している。
Quick start
https://github.com/holtjma/fmlrc/wiki/Quick-start-test
インストール
依存
- msbwt
- ropebwt2
#ラン前にシーケンシングデータのBWTを構築する必要がある。msbwtとropebwt2を導入しておく。
easy_install msbwt
git clone https://github.com/lh3/ropebwt2.git
cd ropebwt2/
make
./ropebwt2 -h
#パスの通ったディレクトリに移動しておく。
本体 Github
https://github.com/holtjma/fmlrc
git clone https://github.com/holtjma/fmlrc.git
cd fmlrc/
make
./fmlrc -h
#bioconda (link)
conda install -c bioconda -y fmlrc
$ ./fmlrc -h
Usage: fmlrc [options] <comp_msbwt.npy> <long_reads.fa> <corrected_reads.fa>
Options: -h print help menu
-v print version number and exit
-k INT small k-mer size (default: 21)
-K INT large K-mer size (default: 59), set K=k for single pass
-p INT number of correction threads
-b INT index of read to start with (default: 0)
-e INT index of read to end with (default: end of file)
-m INT absolute minimum count to consider a path (default: 5)
-f FLOAT dynamic minimum fraction of median to consider a path (default: .10)
-B INT set branch limit to <INT>*<k or K> (default: 4)
-i build a sampled FM-index instead of bit arrays
-F INT FM-index is sampled every 2**<INT> values (default: 8); requires -i
-V verbose output
> ropebwt2 -h
$ ./ropebwt2 -h
./ropebwt2: illegal option -- h
Usage: ropebwt2-r187 [options] <in.fq.gz>
Options: -l INT leaf block length [512]
-n INT max number children per internal node [64]
-s build BWT in the reverse lexicographical order (RLO)
-r build BWT in RCLO, overriding -s
-m INT batch size for multi-string indexing; 0 for single-string [10g]
-P always use a single thread
-M INT switch to single thread when < INT strings remain in a batch [1000]
-i FILE read existing index in the FMR format from FILE, overriding -s/-r [null]
-L input in the one-sequence-per-line format
-F skip forward strand
-R skip reverse strand
-N skip sequences containing ambiguous bases
-x INT cut at ambiguous bases and discard segment with length <INT [0]
-C cut one base if forward==reverse
-q INT hard mask bases with QUAL<INT [0]
-o FILE write output to FILE [stdout]
-b dump the index in the binary FMR format
-d dump the index in fermi's FMD format
-T output the index in the Newick format (for debugging)
ラン
１、BWTを構築。
gunzip -c short_read.fq.gz | awk 'NR % 4 == 2' | sort | tr NT TN | ropebwt2 -LR | tr NT TN | msbwt convert output
outout/にcomp_msbwt.npyができる。
２、エラー訂正を実行する。
fmlrc comp_msbwt.npy long_reads.fa corrected_reads.fa
参考
メモリ効率がよくないとされるがもmsbwt（リンク）も使える。
git clone https://github.com/holtjma/msbwt
easy_install msbwt
msbwt -h
$ msbwt -h
usage: msbwt [-h] [-V]
{cffq,pp,cfpp,merge,query,massquery,compress,decompress,convert}
...
A multi-string BWT package for DNA and RNA.
positional arguments:
{cffq,pp,cfpp,merge,query,massquery,compress,decompress,convert}
cffq create a MSBWT from FASTQ files (pp + cfpp)
pp pre-process FASTQ files before BWT creation
cfpp create a MSBWT from pre-processed sequences and offsets
merge merge many MSBWTs into a single MSBWT
query search for a sequence in an MSBWT, prints sequence and seqID
massquery search for many sequences in an MSBWT
compress compress a MSBWT from byte/base to RLE
decompress decompress a MSBWT from RLE to byte/base
convert convert from a raw text input to RLE
optional arguments:
-h, --help show this help message and exit
-V, --version show program's version number and exit
BWTを構築。
msbwt cffq outout/ pair1.fq pair2.fq
出力されるoffsets.npyを使う。
引用
FMLRC: Hybrid long read error correction using an FM-index.
Wang JR, Holt J, McMillan L, Jones CD
BMC Bioinformatics. 2018 Feb 9;19(1):50.
比較のプレプリントが出ています。
https://www.biorxiv.org/content/biorxiv/early/2019/05/29/519330.full.pdf
