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近くある。

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

横軸がリードの長さ(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のサイズ分布をヒストグラムで調べる。

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

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

アセンブル

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

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

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

E.coli=> 3 contig

S. cerevisiae=> 27 contig

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

グラフの左から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

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

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

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

spades

E.coli=> 763 contig

S. cerevisiae=> 1389 contig

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

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

spades

E.coli=>  602 contig

S. cerevisiae=>  1145 contig

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

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