2019 3/5 インストール追記、コマンドのわかりにくい部分を修正
2019 5/14 リンク追加
2019 5/27 docker追加、オプションヘルプ追加
2019 8/27 twitter追記
ハイスループットシーケンシングにより、新規ゲノムのシーケンシングが日常的に可能になっている。しかしながら、これらのゲノムの最も基本的な特徴、例えばサイズまたはヘテロ接合率などは、最初は未知であり、例えばリードマッパー、デノボアセンブラ、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プロファイルをプロットし、ユーザが原因を調べることができる。
インストール
ubuntu16.04でテストした。
依存
#ubuntuならapt-getで1.1を導入できる
apt update && apt install -y jellyfish
git clone https://github.com/schatzlab/genomescope.git
cd genomescope/
> Rscript genomescope.R
$ Rscript genomescope.R -h
USAGE: genomescope.R histogram_file k-mer_length read_length output_dir [kmer_max] [verbose]
> jellyfish count -h
$ jellyfish count -h
Usage: jellyfish count [options] file:path+
Count k-mers or qmers in fasta or fastq files
Options (default value in (), *required):
-m, --mer-len=uint32 *Length of mer
-s, --size=uint64 *Hash size
-t, --threads=uint32 Number of threads (1)
-o, --output=string Output prefix (mer_counts)
-c, --counter-len=Length in bits Length of counting field (7)
--out-counter-len=Length in bytes Length of counter field in output (4)
-C, --both-strands Count both strand, canonical representation (false)
-p, --reprobes=uint32 Maximum number of reprobes (62)
-r, --raw Write raw database (false)
-q, --quake Quake compatibility mode (false)
--quality-start=uint32 Starting ASCII for quality values (64)
--min-quality=uint32 Minimum quality. A base with lesser quality becomes an N (0)
-L, --lower-count=uint64 Don't output k-mer with count < lower-count
-U, --upper-count=uint64 Don't output k-mer with count > upper-count
--invalid-char=warn|ignore|error How to treat invalid characters. The char is changed to a N. (warn)
--matrix=Matrix file Hash function binary matrix
--timing=Timing file Print timing information
--stats=Stats file Print stats
--usage Usage
-h, --help This message
--full-help Detailed help
-V, --version Version
dockerイメージ
https://hub.docker.com/r/greatfireball/ime_genomescope
ラン
まずfastqファイルを分析する。
jellyfish count -C -m 21 -s 1000000000 -t 12 *fq -o reads.jf
- -m Length of mer
- -s Initial hash size
- -t Number of threads (1)
ヒストグラムファイルを出力する。
jellyfish histo -t 12 reads.jf* > reads.histogram
k-merスペクトラムのグラフを出力する。k-mer_maxの推奨値は1000になっている。
cd genomescope/
Rscript genomescope.R reads.histogram 21 301 output_dir
出力
cat outdir/summary.txt
GenomeScope version 1.0
k = 31
property min max
Heterozygosity 0.000695161% 0.0158001%
Genome Haploid Length 5,157,352 bp 5,166,793 bp
Genome Repeat Length 2,071,904 bp 2,075,697 bp
Genome Unique Length 3,085,448 bp 3,091,096 bp
Model Fit 71.0875% 73.9962%
Read Error Rate 0.32858% 0.32858%
下に載せたのはサンプルデータ（シロイヌナズナのF1のシーケンスデータ）の分析結果。縦軸の頻度は分かりやすいようlogになっている。カバレッジの異なる２つのピークが確認できる。duplicationを無視すると、カバレッジが低い方のピークはヘテロ接合のゲノム領域由来のピーク、カバレッジが高い方のピークはホモ接合のゲノム領域由来のピークと考えられる。カバレッジは極端に少ないが、頻度は非常に多いグラフの左端のk-merの分布は、シーケンスエラーとしてフィッティングされている。
上のグラフは片対数。下は両対数。印象が違うのは、両グラフで縦軸の目盛り間隔が異なるため。
dockerイメージも利用できる。
docker run --rm -itv $PWD:/data/ -w /data \ greatfireball/ime_genomescope reads.histogram 21 300 output
histファイルを分析するwebサーバーも用意されている。
GenomeScope online
ヒストグラムテキストファイルを指定すれば分析できる。
引用
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.
追記
GenomeScope 2.0 and Smudgeplots: Reference-free profiling of polyploid genomes https://t.co/XVBuc0F0R5 #biorxiv_bioinfo
— bioRxiv Bioinfo (@biorxiv_bioinfo) August 26, 2019
