macでインフォマティクス

macでインフォマティクス

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

シーケンスデータからk-merスペクトラム分析を行う GenomeScope

 

 ハイスループットシーケンシングにより、新規ゲノムのシーケンシングが日常的に可能になっている。しかしながら、これらのゲノムの最も基本的な特徴、例えばサイズまたはヘテロ接合率などは、最初は未知であり、例えばリードマッパー、デノボアセンブラ、SNPコーラーなどについて適切な分析方法を選択することが困難である(論文より Smolka et al、2015)。これらの特性をあらかじめ決定すれば、変異の数を過小報告したり、ゲノムのかなりの部分をアセンブルできないなど、解析方法がゲノムの複雑性を捕捉をできていないかどうかを明らかにすることができる。これらの特性のいくつかを測定する実験が利用可能であるが、かなりの費用と労力を必要とする可能性がある。

 アセンブルされていないシークエンシングリードからゲノムサイズを推定するための計算手法は現在ほとんどない(Chikhi and Medvedev、2014; Melsted and Halldorsson、2014)。いくつかのゲノムアセンブラアルゴリズムを導くために関連統計を内部的に計算する(Bankevich et al、2012; Gnerre et al、2011)。これらの方法は、BACの長さをライブラリのショットガンサンガーシークエンシングから推測する以前の研究に従っている(Li and Waterman、2003)。しかしながら、ヘテロ接合性の割合などのより複雑な特性を測定するための方法はごくわずかしかなく、これらの方法は使用または解釈が難しい場合がある。 Simpson(2014)は、デノボアセンブリ技術を用いてシーケンシングリードからこれらの特性のいくつかを推定するための計算方法を提案した。しかしながら、この方法は計算集約的であり、より直接的なヘテロ接合率ではなく、変異誘発分枝率のような結果がアセンブリグラフに関連して報告されるので、解釈が難しい場合がある。 GCE法(Liu et al、2013)もまた、ゲノムサイズおよびヘテロ接合率を決定しようと試みるが、完全に自動化されておらず、ユーザがいくつかのカットオフを手動で指定する必要がある。また、実際のシークエンシングデータで推定値が低くなる可能性のあるポアソンカバレッジモデルを使用し、ヘテロ接合性モデルは反復配列のないゲノムに限定されている。

 ここでは、ゲノム特性(全ゲノム長と一倍体ゲノム長、リピート含量とヘテロ接合率のパーセンテージ)、raw シーケンシングデータから全体的なリード特性(カバレッジ、重複とエラー率)を参照ゲノムを必要とせず、k-merプロフィールの統計的分析によって自動的に推測するGenomeScopeを紹介する。 k-merプロファイル(k-merスペクトラムとも呼ばれる)は、k-mers; 長さkの部分文字列がシーケンスリードでどのくらいの頻度で発生するかを測定するもので(一部略)、プロファイルはゲノムの複雑さを反映している:ホモ接合ゲノムは単純なポアソンプロファイルを有するが、ヘテロ接合型ゲノムは特徴的な二峰性のプロファイルを有する(Kajitani et al、2014)。反復はさらに高いk-merスペクトラムでピークを追加するが、シークエンシングエラーとリードのduplicationは低いスペクトラムのfalse k-merと分散の増加を伴いプロファイルを歪ませる(Kelley et al、2010; Miller et al、2011)。

 これらの潜在的複雑さを認識して、本ツールは、ヘテロ接合およびホモ接合、ユニークおよび2コピーの配列の相対存在量を測定するために、k-merプロファイルに4つの等間隔の負の二項分布の混合モデルを適合させる(論文 supplementary S2)。 ポアソン分布に比べて実際のシークエンシングデータはより分散していることが多いため、GenomeScopeはポアソン分布より負の二項分布の混合モデルを使用する(Miller et al、2011)。フィッティングは、Rのnls関数によって実装された非線形最小二乗推定値を使用して計算される(Bates and Watts、1988)。モデルフィッティングをより強固にするために、ゲノムスコープは、シークエンシングエラーによって引き起こされる可能性のある低頻度k-merの異なる部分を除いて、数回のモデルフィッティングを試み、正しいヘテロ接合性およびホモ接合性ピークを決定する際のあいまいさを調整する。パラメータの最終セットは、観察されたk-merプロファイルに対するモデルの残差平方和誤差(RSSE)を最小にするパラメータとして選択される。その後、シーケンスエラーおよびより高いコピー反復は、モデル範囲外のk-mersによって同定され、ホモ接合配列の平均カバレッジ値に対して観察されたk-mer頻度を正規化することによって、全ゲノムサイズを推定する。

GenomeScopeは、コマンドラインのRアプリケーションとしても、また使いやすいWebアプリケーションとしても利用できる。GenomeScopeのコマンドラインまたはオンラインバージョンは、一般的に1分未満で完了し、必要なメモリ量は少なく、推定ゲノム特性を持つテキストファイルも出力する。モデルが収束しない場合、通常、カバレッジが低いか、低品質なため、モデルパラメータを表示せずにk-merプロファイルをプロットし、ユーザが原因を調べることができる。

 

インストール

依存

 

本体  Github

https://github.com/schatzlab/genomescope

git clone https://github.com/schatzlab/genomescope.git
cd genomescope/

Rscript genomescope.R 

$ Rscript genomescope.R 

USAGE: genomescope.R histogram_file k-mer_length read_length output_dir [kmer_max] [verbose]

 

ラン

まずfastqファイルを分析する。

jellyfish count -C -m 21 -s 1000000000 -t 10 *fq -o reads.jf

ヒストグラムファイルを出力する。

jellyfish histo -t 10 reads.jf > reads.histo

k-merスペクトラムのグラフを出力する。k-mer_maxの推奨値は1000になっている。

Rscript genomescope.R histogram_file 21 301 output_dir [kmer_max]

 

 

サンプルデータ(シロイヌナズナのF1のシーケンスデータ)分析結果。縦軸の頻度は分かりやすいようlogになっている。カバレッジの異なる2つのピークが確認できる。duplicationを無視すると、カバレッジが低い方のピークはヘテロ接合のゲノム領域由来のピーク、カバレッジが高い方のピークはホモ接合のゲノム領域由来のピークと考えられる。カバレッジは極端に少ないが、頻度は非常に多いグラフの左端のk-merの分布は、シーケンスエラーとしてフィッティングされている。

f:id:kazumaxneo:20180428224309j:plain

上のグラフは片対数。下は両対数。印象が違うのは、両グラフで縦軸の目盛り間隔が異なるため。

f:id:kazumaxneo:20180428224319j:plain

 

 

histファイルを分析するwebサーバーも用意されている。

GenomeScope online

f:id:kazumaxneo:20180428210201j:plain

ヒストグラムテキストファイルを指定すれば分析できる。

 

 

引用

GenomeScope: fast reference-free genome profiling from short reads.

Vurture GW, Sedlazeck FJ, Nattestad M, Underwood CJ, Fang H, Gurtowski J, Schatz MC

Bioinformatics. 2017 Jul 15;33(14):2202-2204.