macでインフォマティクス

macでインフォマティクス

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

Oxford Nanoporeリードのアセンブリ パフォーマンス比較

 

S. cerevisiaeとE.coli K12をilluminaとnanoporeで読んだ論文のデータ が公開されている。

http://schatzlab.cshl.edu/data/nanocorr/

 

このデータを使い、2017年6月現在のナノポアリードのパフォーマンスを調べてみる。

上記URLからMinionとilluminaのシーケンスデータ (fastq) をダウンロードした。 

 

 

E.coliのリード数などは以下の通り。

・Oxford nanopore Minion: 59008 reads、total_length 240Mb、average 4068 bp、max 43000bp、min 120bp

・illumina Miseq (301x2): 22Mx2 reads、total_length 13.4Gb、max 301 bp

ナノポアのリードは平均4000bpを超えているが、リードのトータルサイズが2.4億bpで、カバレッジは48xにしかならない。一方でMiseqのトータルサイズは301bp計算で67億bpを超えている!。カバレッジを計算すると2680xで、かなりディープに読まれたデータだ。gzから解凍したファイルサイズはそれぞれ15GB近くある。

ロングリードのサイズ分布をヒストグラム解析して調べてみた。

f:id:kazumaxneo:20170622150029j:plain

横軸がリードの長さ(bp)で、縦軸がその長さのリードが出てくる頻度である。(右端は単位を5000bpに変えていることに注意してください)。一番多いのは1000bp程度のリードであるが、綺麗な二項分布にはなっていない。6000-7000付近にも弱いピークがあるように見える。

やはりとてつもなく長い。1リードで5000bp以上の長さがあれば、rrnaのクラスターもいきなりアセンブルできる可能性が高い。

 

S. cerevisiaeのリード数は以下の通り。

・Oxford nanopore Minion: 88218 reads、total_length 526Mb、average 5969 bp、max 72879bp、min 103bp

・illumina Miseq (301x2): 25Mx2 reads、total_length 15Gb、max 301 bp

S. cerevisiaeのゲノムサイズを大雑把に12.5Mbとした時、ナノポアのカバレッジは44x、Miseqのカバレッジは1300xとなる。

ナノポアのlong readのサイズ分布をヒストグラムで調べる。

f:id:kazumaxneo:20170622150005j:plain

横軸がリードの長さ(bp)で、縦軸がその長さのリードが出てくる頻度である。(同じく右端は単位を5000bpに変えていることに注意してください)。一番多いのは1000bp程度のリードであるが、7000付近にも強いピークがある。

推測だが、7000bp付近のリードが本来のナノポアMinionのパフォーマンスで、左端のピークはDNA抽出時に物理的に切れてしまったDNAではないだろうか?(海外研究者のナノポアのワークショップの投稿で、λ HInd III cutをナノポアで読むと、等カバレッジにならずサイズが長くなるほどカバレッジが下がという話がある)。

 

 

アセンブル

以下はE.coliの例であるが、 S. cerevisiaeもゲノムサイズ以外は同じパラメータで実行した。 

1、canuを使いONTリードだけでアセンブル。

canu -p test -d long_read genomeSize=4.8m -nanopore-raw oxt.fq

E.coli=> 3 contig

S. cerevisiae=> 27 contig

contigの長さをグラフにすると以下のようになる。

f:id:kazumaxneo:20170623112456j:plain

グラフの左からcontigの長さ順に示す。縦軸がcontigの長さ(bp)である。

 

S.cerevisiaeは17クロモソーム持つので、クロモソームのサイズとcontigの長さを比較した。下の表は左側がS. cerevisiaeのゲノムをサイズ順に並べたもので、右側がナノポアリードをアセンブルして作ったcontigをサイズ順に並べたものである。

chr Size (bp)   contig Size (bp)
IV 1530000   tig00000001  1546985
VII 1090000   tig00000003  1091026
XV 1090000   tig00000007  1086779
XII 1080000   tig00000006  939389
XVI 948066   tig00000013  805656
XIII 924431   tig00000022  737318
II 813184   tig00013668  678294
XIV 784333   tig00000024  644519
X 745751   tig00000053  573971
XI 666816   tig00013665  519780
V 576874   tig00000181  487116
VIII 562643   tig00000055  462121
IX 439888   tig00000058  429563
III 316620   tig00000072  271056
VI 270161   tig00000185  259055
I 230218   tig00000066  239610
MT 85779   tig00000183  237080
      tig00000050  190760
      tig00000075  188868
      tig00000035  187180
      tig00000184  162871
      tig00000086  155346
      tig00000042  140091
      tig00000036  110210
      tig00013670  86062
      tig00000009  52963
      tig00000046  25807

おそらく多くのクロモソームは1つにアセンブルされている。quastでエラーチェックも行なったが、リファレンスゲノムと比べて明確な構造変化は起こっていない。つまり、あとはPCRを総当たりすれば、ゲノムの構造は決まってしまう。かなり衝撃的な結果である。(ただしエラーチェックは入念に行う必要がある!)

E.coli K12の方は、3つのcontigにアセンブルされていた(E.coli K12 MG1655のゲノムは4.64Mbのクロモソームのみを持つ。プラスミドはない)。

 

  

2、spadesを使いilluminaのペアリードをアセンブル。

spades

E.coli=> 763 contig

S. cerevisiae=> 1389 contig

f:id:kazumaxneo:20170623111700j:plain

 

lengthのtop100に絞って表示する。縦軸がcontigの長さ(bp)である。

 

 

3、spadesを使いONTリードとilluminaリード両方でhybirdアセンブル。

spades

E.coli=>  602 contig

S. cerevisiae=>  1145 contig

E.coliのcontigは一番長いもので4.64Mbで、残りは2kb以下のcontigだった。グラフにすると下のようになる。

f:id:kazumaxneo:20170623112142j:plain

lengthのtop30に絞って表示する。縦軸がcontigの長さ(bp)である。