誤りが見つかったため、初投稿からいくつか内容を修正しています。
R7のデータがnanopolishで解析できなかったため、 テストデータについてもR9で読んだデータを使うように修正しました。
オックスフォードNanoporeシーケンシングデータのシグナルレベル解析のためのツール。 ナノポリッシュは、ドラフトゲノムアセンブリのコンセンサス配列を計算したり、塩基修飾を検出し、参照ゲノムおよびそれ以上についてSNPおよびindelsを呼び出すことができる。ショートリードを使わずに、ロングリード自身でナノポアリードでアセンブルしたcontigのエラーを修復し精度を上げるコンセンサス配列出力ツールとしてよく用いられる。
追記
- インストール
- fast5からindexを作成するstepなどが追加されているようなので、時間があるときテストして追記します。
インストール
依存
ビルド
git clone --recursive https://github.com/jts/nanopolish.git
cd nanopolish
make
またはbioconda環境ならcondaで簡単インストールできる(リンク) 。
conda install -c bioconda -y nanopolish
bioconda導入についてはこちらの方の説明を参考にしてください。分かりやすく説明されています(リンク)。
> nanopolish -h
# nanopolish -h
usage: nanopolish [command] [options]
valid commands:
--help
--version
call-methylation
eventalign
extract
getmodel
help
index
methyltrain
phase-reads
scorereads
variants
for help on given command, type nanopolish command --help
実行方法
ミスマッチが修復されるだけでなく、contigが伸びたりスキャホールドが長くなることもある。ドラフトゲノムをpolishしてクオリティを上げるワークフロ-は以下のようになる。
1、FAST5からfasta(fastq)を抽出
↓
2、de novo assembly
↓
3、nanoporeのリードをcontigにマッピングする
↓
4、nanopolishでゲノムを50kbのセグメントに分け、コンセンサス配列を分析。
↓
5、nanopolishで各セグメントをマージし、polishしたゲノム (.fa) を出力。
論文化されているナノポアのデータを使ってnanopolishの機能を確認していく。ArabidopsisをシーケンスしたデータがENAに上がっているので、これを使う。
http://www.ebi.ac.uk/ena/data/view/PRJEB8681
Minionは電圧をかけている間だけポアをDNAが通ってデータが読まれていくので、任意の時間だけランできる(つまりポアが生きている限り何度でもランできる)。Minionで何時間シーケンスされたか情報はなかったが、ダウンロードしたraw_readディレクトリには、passに9560のfast5ファイル(29.1 GB)、failに36000以上のfast5(60.5 GB)があった。
ダウンロードと解凍には丸一日かかった。
解凍するとフォルダの中には生データのfast5が入っている(一定の間隔でデータが区切られて出力されている)。
前置きが長くなったが、ここからがnanopolsihの出番である。
- 1、FAST5からfasta(fastq)を抽出
まずは nanopolishのextractを使い、FAST5をfastq、またはfastaに変換する作業から始める。
nanopolish extract -t any flowcell_17/downloads/ -o reads.fa
- -v display verbose output
- -r recurse into subdirectories
- -q extract fastq (default: fasta)
- -t read type: template, complement, 2d, 2d-or-template, any (default: 2d-or-template)
- -oFILE write output to FILE (default: stdout)
生データは非常にサイズが大きく捨てたくなるが、後で参照するので捨ててはいけない。
nanopolishのanyオプション付きと、nanoporeの抽出でファイルサイズに6倍くらい差があった。
nanopolish extract -t any => 19118 line、117 MB、117074353 charactor (wc -m)
poretools fasta => 57354 line、 326 MB、326362638 charactor (wc -m)
nanoplosihとは関係ないが、nanopolishとporetoolsでfast5から抽出されたリードに違いが出た理由を知りたい。bwa memを使いchloroplstゲノムにアライメントしてみる。
上半分がPoretoolsで抽出したリードで、下半分がnanopolishで抽出したリードである。nanopolishの方がリードは少ないが、綺麗である。抽出されたリードの数になぜここまで差がつくのか、今のところ分からない。
次にナノポアのロングリードをアセンブルしてcontigを作る。
- 2、de novo assembly
ナノポアリードをアセンブルしてcontigを出力する。ナノポアリードのアセンブルについては、以下のエントリをー参考にしてください。
- エラーコレクションのあるcanuでアセンブルする。
- エラーコレクションはないが、高速なMiniasmでアセンブルする。
他にもspadesを使ってショートリードとナノポアのロングリードをhybrid assemblyする方法も論文でたびたび見かける。
今回はcanuでアセンブルした。contigが1つだけ出力されたのでblastをかけると、chloroplastゲノムに高い相同性を示した。このことから、ナズナのクロロプラストゲノムがアセンブルされたことがわかった。説明がないので不明だが、何らかの方法でクロロプラストが単離されてDNA抽出されているのか、クロロプラストゲノムのコピー数が多いためクロロプラストゲノムが大半を占めてしまっているのかもしれない(120Mのゲノムをアセンブルするにはカバレッジが足りない)。
- 3、nanoporeのリードをcontigにマッピングする
ナノポアのロングリードのアライメントにはbwa mem の"-x on2d"を使う。
bwa index ecoli_draft.fa #1 contig.faのindexを作る
bwa mem -x ont2d -t 12 ecoli_draft.fa reads.fa | samtools sort -o reads.sorted.bam -T reads.tmp - #2
samtools view -H reads.sorted.bam|head -1 #3 sortの確認
samtools index reads.sorted.bam #4 bamのindexを作る
- #1 e_coli_draft.fa => アセンブルして作ったcontig
- #2の前半 -tでbwa memのスレッド数を指定
- #2のパイプ 今回はbwa memの結果をパイプでsamtools sortにつないでいる (samtools viewは必要ない。samtools sortはsamを扱える)。samを一回出力してsortしても同じである。
- #2の後半 マルチfastaにアライメントしたsam/bamをsamtools sortすると、フォルダに一時的に中間ファイルが出力される。上では-Tで中間ファイルが画面に出力されないように工夫している。"-T file_name"がなくても結果は同じ。
- #2の後半 ハイフン(-)は"|"でパイプした時は標準入力として扱われる。わかりやすいよう書かれているだけで、無くても結果は同じである(参考)。また、右端に書かれているのは、"-o"の手前だと"- -o"となり可読性が悪いからである。
- #3 samtools view -Hでsam/bamのヘッダー部分を確認できる。"SO:coordinate"という記載が1行目にあればソートできている。"-n"をつけてソートすると"SO:queryname"になる。
- #4 #2で出力されたsort.bamのidne作り
以上でアライメント作業は終了。
- 4、nanopolishでゲノムを50kbのセグメントに分け、コンセンサス配列を出力。
ここがnanopolishのコアになろうドラフトゲノムのキュレーションのステップである。
python nanopolish_makerange.py ecoli_draft.fa | parallel --results nanopolish.results -P 8 nanopolish variants --consensus polished.{1}.fa -w {1} -r reads.fa -b reads.sorted.bam -g draft.fa -t 4 --min-candidate-frequency 0.1
パイプの前半
- このコマンドnanopolish_makerange.pyはbrewでインストールすると見つからないかもしれない。ない人はgitからソースコードをダウンロードして、src/に入っているこのコマンドをフルパス指定して走らせる。
- nanopolish_makerange.pyはおそらくリファレンスのサイズを50kbp(正確には500200 bp?)ごとに区切る数値を出すのコマンドである。
パイプの後半
- parallelという作業を並列化するコマンドを使っている。parallelの後に並列化させるコマンドを打ち込む。
長くて見にくいコマンドだが、パイプ後半のparallelは作業を並列化させるために使っている。よって後半は実質、下のようなコマンドになる。
nanopolish variants --consensus polished.{1}.fa -w {1} -r reads.fa -b reads.sorted.bam -g draft.fa -t 4 --min-candidate-frequency 0.1
nanopolish variantsはドラフトゲノムからpolishされたコンセンサス配列を出力するコマンド。ただし並列化が難しく1ファイルだけで長い配列を計算させると時間がかかるので、50kbpごとにわけ(200bp重複させ)、それをパラレルに計算させることで高速化を図っている。必要な引数はドラフトゲノム(-g draft.fa)、ロングリード(-r reads.fa)、3で作ったbamファイル(-b reads.sorted.bam)、である。version 0.6.0以前は--modelファイルも必要だったらしい(biostarログ)。
上記では--min-candidate-frequency 0.1(バリアントの割合が0.1=10%以上の部位を出力。defaultでは20%)、-t 4(use NUM threads)もつけている。
4のコマンドは全体として不安定で、biostarにはエラーの質問が多く上がっている。原因はparallelをnanopolishがうまく扱えないために起こることが多いらしい。私もエラーを吐いたので色々調べたところ、使っていたHDF5の最新バージョンに微妙なバグがあるらしいことがわかった(見つけたリンク)。またこれとは別に、mInionのR7でシーケンスされたデータがnanopolishでpolishできないというスレッドも上がっていた(今回のナズナのデータはMinionの R9でシーケンスされたため、その症状は回避できた)。Minionの古いデータを使う場合は、githubのバージョンで残っているnanopolishの0.5以下の古いバージョンをテストするとうまく行くかもしれない。
上記のコマンドが終わると、カレントディレクトリに50kb単位でpolishされたfastaファイルが出力される。nanopolish.resultsディレクトリにはログが残る。最後にマージする。リファレレンスゲノムがはっきりしないため分からないが、今回のpolishのテストでは、300箇所以上の塩基が置換され、3000箇所以上でgap修正が行われた。ロングリードのエラーがランダムなことを利用しており、精度の良い解析にはある程度のカバレッジが必要。
他にも以下のような機能がある。
- メチル化サイトをコールする。
nanopolish call-methylation: predict genomic bases that may be methylated
- SNVとindelをコールする。
nanopolish variants: detect SNPs and indels with respect to a reference genome
追記
Sparcも参考にしてください。
引用(修正 )
A complete bacterial genome assembled de novo using only nanopore sequencing data
Loman NJ, Quick J, Simpson JT
Nat Methods. 2015 Aug;12(8):733-5