macでインフォマティクス

macでインフォマティクス

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

ONTリードのレイアウトを計算しコンセンサス配列を出力する spectrassembler

 

 イルミナのようなハイスループットシーケンシング技術は、リード長は犠牲になるがシーケンシングコストは大幅に減らすことができる。しかし、ゲノムにはリピート配列が含まれているため、デノボアセンブリで問題になる。PacBioのSMRTやOxford Nanopore Technology(ONT)などの最近のロングリードシーケンシング技術は、10 kb以上の長さのリード(論文より Koren and Phillippy、2015)を生み出すため、デノボアセンブリルネサンスを促した。しかし、高いエラー率(〜15%)は、アセンブリ作業を困難にし、複雑かつ計算集約的なパイプラインを必要とする。ロングリードアセンブリのほとんどのアプローチ、例えばCanu (Koren et al., 2016)(以前のCelera Assembler (Myers et al., 2000)) は、アセンブリを実行する前にリードを修正することでこの問題に対処する。

 ハイブリッド技術は、ショートリードおよびロングリードを組み合わせる。この方法は、PacBio(Koren et al、2012)およびONT(Goodwin et al、2015)データを用いて原核生物および真核生物のゲノムを首尾よく組み立てることが示された。階層的アセンブリマッピングとコンセンサスの原則に従う。十分なカバレッジが存在し、シーケンシングエラーがランダムであれば、例えばゲノム上のある位置の塩基が、50リードのうち8リードが間違っている場合、多数票で正しい塩基がもたらされる。階層的方法は、ロングリードをお互いにマップし、各リードについて、それと重複するすべてリードに基づくコンセンサス配列を導く。このようなアプローチは、PacBio SMRTデータを組み立てるためにHGAP(Chin et al。、2013)に実装され、さらに最近ではLomanら(2015)のONTデータを用いた大腸菌の完全なデノボアセンブリの例がある。

 最近、Li(2016)は、 all-versus-all のrawロングリードの minimapによるマッピングとminiasmによるアセンブルの2ステップだけでノイズの多いリードのデノボアセンブリが実行可能であることを示している 。このアセンブラはCelera Assemblerからインスパイアされ、アセンブリー・グラフの構築を通じてユニットを生成する。制限は、エラー率が未処理のリードと同じオーダーのドラフトであることにある。spectrassemblerは、raw ONTリードのレイアウトを計算するための新しい方法である。GooglePageRank(Page et al、1999)に似た単純なスペクトルアルゴリズムに基づいている。スペクトラムアルゴリズムを使用して、ペアワイズオーバラップ情報からリードの順序を見つける。これはゲノムの物理的マッピング(Atkins and Middendorf、1996)、先祖ゲノムの再構成(Jones et al、2012)、またはlocus ordering problem (Cheema et al、2010)に首尾よく適用され、デノボのアセンブル上の問題には適用されていない。この方法は、カバレッジがあまり高くない、いくつかのバクテリアノムにおいて、ONTリードから質の高いゲノムサイズのコンティグ(〜99%の同一性)を作ることが示されている。

Githubでは、このツールは実験的なものであり、メインストリームのアセンブラを置き換えるものではない、と記載されている)。

 

 

インストール 

cent OSに導入した。

依存

pythonのライブラリはpipで導入する。"pip install biopython"

pip install numpy
pip install scipy
pip install biopython

minimapはcloneしてビルドする。

git clone https://github.com/lh3/minimap && (cd minimap && make)

本体 Github

https://github.com/antrec/spectrassembler

git clone https://github.com/antrec/spectrassembler
cd spectrassembler && srcd=`pwd` && git submodule init && git submodule update &&
cd tools/spoa && git submodule init && git submodule update && make
cd ../../../

Bioconda環境ならcondaでも導入できる(リンク)。(linux only)

 > python spectrassembler.py -h

$ python spectrassembler.py -h

usage: spectrassembler.py [-h] [-r ROOT] -f READS_FN -m MINIMAPFN

                          [--min_cc_len MIN_CC_LEN] [--w_len W_LEN]

                          [--w_ovl_len W_OVL_LEN] [--len_thr LEN_THR]

                          [--sim_qtile SIM_QTILE] [-v]

                          [--ref_pos_csvf REF_POS_CSVF] [--spoapath SPOAPATH]

                          [--nproc NPROC] [--margin MARGIN]

                          [--trim_margin TRIM_MARGIN] [--julia JULIA]

 

De novo experimental assemblerbased on a spectral algorithm to reorder the

reads

 

optional arguments:

  -h, --help            show this help message and exit

  -r ROOT, --root ROOT  directory where to store layout and consensus files.

  -f READS_FN, --READS_FN READS_FN

                        path to reads file (fasta or fastq)

  -m MINIMAPFN, --minimapfn MINIMAPFN

                        overlap file path (from minimap in PAF format).

  --min_cc_len MIN_CC_LEN

                        minimum number of reads for a contig to be considered

  --w_len W_LEN         length of consensus windows for POA.

  --w_ovl_len W_OVL_LEN

                        overlap length between two successive consensus

                        windows.

  --len_thr LEN_THR     threshold on length of overlaps (similarity matrix

                        preprocessing).

  --sim_qtile SIM_QTILE

                        quantile threshold on overlap score (similarity matrix

                        preprocessing.)0.5 means you keep only overlaps with

                        num_match > quantile(num_matches, 0.5)

  -v, --verbosity       verbosity level (-v, -vv or none)

  --ref_pos_csvf REF_POS_CSVF

                        csv file with position of reads (in same order as in

                        READS_FN)obtained from BWA, in order to plot reads

                        position found vs reference.

  --spoapath SPOAPATH   path to spoa executable

  --nproc NPROC         number of parallel processes

  --margin MARGIN       number of bases to add to current consensus to make

                        sure it overlaps next window

  --trim_margin TRIM_MARGIN

                        length to cut in beginning and end of consensus

                        sequences from spoa (where the consensus isless good)

  --julia JULIA         path to Julia (optional,though eigenvector

                        computations are clearly faster in Julia than in

                        Python)

またはcondaでも導入できる(リンク)。

 

ラン

テストデータをダウンロードしてランしてみる。

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

mkdir oxford-test && cd oxford-test
curl -L -o oxford.fasta http://nanopore.s3.climb.ac.uk/MAP006-PCR-1_2D_pass.fasta

 

2、minimapでアセンブリ

minimap -Sw5 -L100 -m0 oxford.fasta oxford.fasta > oxford.mini.paf

 

3、spectrassemblerをラン。

python spectrassembler.py -f oxford.fasta -m oxford.mini.paf -r temp -v --nproc 12 > contigs.fasta
  • -f path to reads file (fasta or fastq)
  • -m overlap file path (from minimap in PAF format).
  • -r  directory where to store layout and consensus files.
  • -v  verbosity level (-v, -vv, -vvv or none), default -v]
  • --nproc    number of parallel processes
  • --sim_qtile quantile threshold on overlap score (similarity matrix preprocessing.) 0.5 means you keep only overlaps with um_match > quantile(num_matches, 0.5)

デフォルトは40%の最小値(--sim_qtile 0.4)を削除する)を使用して、重複からの類似性マトリックスを構築する(カバレッジがあまり高くないバクテリアゲノムのONTデータ)。コンティグはstdoutに書き込まれる。高カバレッジ(> 80x)の真核生物のゲノムまたはPacbioデータの場合は、--sim_qtile 0.9オプションを使用して値の90%をゼロに設定する方がよいと書かれている。

 

 

引用

A spectral algorithm for fast de novo layout of uncorrected long nanopore reads.

Recanati A, Brüls T, d'Aspremont A

Bioinformatics. 2017 Oct 15;33(20):3188-3194.