macでインフォマティクス

macでインフォマティクス

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

ナノポアのアセンブルデータのキュレーション及び変異の検出 nanopolish

 

ナノポアリードでアセンブルしたcontigのエラーを修復し精度を上げるためのプログラム。変異のコールや推定メチル化サイトの検出の行うことができる。

インストールから動作まで見ていく。

 

誤りが見つかったため、初投稿からいくつか内容を修正しています。

R7のデータがnanopolishで解析できなかったため、 テストデータについてもR9で読んだデータを使うように修正しました。

 

ダウンロード

Github HP

GitHub - jts/nanopolish: Signal-level algorithms for MinION data

 

インストー

 依存するもの

解凍したディレクトリのルートでmakeするだけでbiopython以外の依存するパッケージは自動導入してくれる。

make 

mac ではソースコードからのビルドでエラーを起こすことが多いらしく、brewで導入する方が良いと思われる(エラーについてはbiostarの関連スレッドを参照)。他に並列実行するためparallelを使う。brewでインストールしておく。2017年6月にbrewで導入されたバージョンは0.6.0だった。

 

 

ラン 

ミスマッチが修復されるだけでなく、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が入っている(一定の間隔でデータが区切られて出力されている)。

f:id:kazumaxneo:20170625002046j:plain

前置きが長くなったが、ここからが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ゲノムにアライメントしてみる。

f:id:kazumaxneo:20170704123616j:plain

 上半分がPoretoolsで抽出したリードで、下半分がnanopolishで抽出したリードである。nanopolishの方がリードは少ないが、綺麗である。抽出されたリードの数になぜここまで差がつくのか、今のところ分からない。

 

 

次にナノポアのロングリードをアセンブルしてcontigを作る。

  • 2、de novo assembly

ナノポアリードをアセンブルしてcontigを出力する。ナノポアリードのアセンブルについては、以下のエントリをー参考にしてください。

 

  • エラーコレクションはないが、高速なMiniasmでアセンブルする。

他にもspadesを使ってショートリードとナノポアのロングリードをhybrid assemblyする方法も論文でたびたび見かける。

今回はcanuでアセンブルした。contigが1つだけ出力されたのでblastをかけると、chloroplastゲノムに高い相同性を示した。このことから、ナズナのクロロプラストゲノムがアセンブルされたことがわかった。説明がないので不明だが、何らかの方法でクロロプラストが単離されてDNA抽出されているのか、クロロプラストゲノムのコピー数が多いためクロロプラストゲノムが大半を占めてしまっているのかもしれない(120Mのゲノムをアセンブルするにはカバレッジが足りない)。

 

 

ナノポアのロングリードのアライメントには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

 

 

 

 

2, Oxford Nanopore sequencing, hybrid error correction, and de novo assembly of a eukaryotic genome

Sara Goodwin,1 James Gurtowski,1 Scott Ethe-Sayers, Panchajanya Deshpande, Michael C. Schatz, and W. Richard McCombie

Genome Res. 2015 Nov; 25(11): 1750–1756. doi: 10.1101/gr.191395.115

 

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