macでインフォマティクス

macでインフォマティクス

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

アライメントフリーでk-merデータベースから高速にバリアントを検出する FastGT

small indel k-mer 高速なツール ゲノム比較 (comparative genomics) human genome reference free human whole genome human exome download Nature Scientific Reports 2017

 

ゲノム変異の研究には、次世代シーケンシング（NGS）技術が広く使用されている。ヒトゲノムの変異は、通常、配列決定されたリードをマッピングし、次いでgenotypeのコールを行うことによって検出される（論文より ref.1-4)。標準的なパイプラインでは、rawシーケンスデータからマルチスレッドサーバー上のバリアントコールまで30倍のカバレッジでヒトゲノムを処理するには40〜50時間が必要となる。マッピングとコールは、多数の利用可能なソフトウェアオプションに精通した熟練者を必要とする最先端のプロセスである。異なるパイプラインがわずかに異なる遺伝子型のコールを生成することは驚くべきことではない（ref.5-9）。幸いにも、一貫性のないgenotypeコールは特定のゲノム領域にのみ関連しており、一方でゲノムの残りの80〜90％での遺伝子型決定は堅牢で信頼性が高い。

　コンピュータが大量のシーケンシングデータをより効率的に扱うことができるため、ゲノム解析におけるk-mers（長さkの部分文字列）の使用が増加した。例えば、すべての既知の細菌の系統樹は、それらのゲノムDNAからのk-merを用いて容易に構築することができる。株特異的なk-mers(ref.16-18)を検索することにより、メタゲノミックデータから細菌株を迅速に同定することができる。 K-mersは、rawリードにおけるシーケンスエラーを修正するためにも使用されている（ref.19-22）。 最近の論文の１つでは、k-mersの頻度をカウントすることに基づくアライメントフリーのSNV呼び出し方法を記載している（pubmed）。このメソッドは、未処理リードをBurrows-Wheelerトランスフォームに変換し、バリアントを囲む可変長の一意の部分文字列を使用してカウントすることによって遺伝子型を呼び出す。

　著者らは、ユニークなk-mersを数えることによって、NGSデータから既知の変異を直接決定する可能性を提供する新しい方法を開発した。この方法は、ゲノムの信頼できる領域のみを使用し、従来のマッピングベースの遺伝子型検出よりも約1〜2桁高速である。したがって、フルスケール分析が完了する前にマーカーのサブセットの迅速で予備的な分析に理想的である。

　このメソッドはC言語で実装され、FastGTソフトウェアパッケージとして利用できる。 FastGTは現在、既知のゲノム変異のコールに限定されている。なぜなら、特定のk-mersはすべての既知の対立遺伝子に対して予め選択されなければならないからである。したがって、従来のマッピングやバリアントコールの代わりにNGSベースのゲノム解析の特定の側面を容易にする補完的な方法ではない。実際、FastGTはNGSデータを入力として使用する大規模デジタルマイクロアレイに匹敵する。この方法は、（1）ユニークなk-merの選択のための手順、（2）FASTQファイルから直接的にk-mersを記憶し、カウントするためのカスタマイズされたデータ構造、および（3）最尤k-mer数から遺伝子型を推定するために特別に設計された方法である。

 

 

公式HP

http://bioinfo.ut.ee/FastGT/index.php?r=site/index

チュートリアル

http://bioinfo.ut.ee/FastGT/index.php?r=site/page&view=manual

k-merデータベース（ヒトゲノムのみ）

http://bioinfo.ut.ee/FastGT/index.php?r=site/page&view=kmers

 

 

公式サイト

http://bioinfo.ut.ee/FastGT/

 

インストール

依存

  • GenomeTester4

http://kazumaxneo.hatenablog.com/entry/2017/12/28/233029

本体 Github

git clone https://github.com/bioinfo-ut/GenomeTester4.git
cd GenomeTester4/
#まずGenomeTester4
cd src/
make clean 
make 

#次にFastGT
make gmer_counter 
make gmer_caller

 またはバイナリをダウンロードする(linux)。

wget http://bioinfo.ut.ee/FastGT/downloads/fastgt_binaries_1.0.tar.gz 
tar zxvf fastgt_binaries_1.0.tar.gz 
cd fastgt_1.0/

>  gmer_counter 

$ gmer_counter 

Nothing to do!

gmer_counter version 4.1 (unstable)

Usage:

  gmer_counter ARGUMENTS SEQUENCES...

Arguments:

    -v | --version   - Print version information and exit

    -db DATABASE     - SNP/KMER database file

    -dbb DBBINARY    - binary database file

    -w FILENAME      - write binary database to file

    -32              - use 32-bit integeres for counts (default 16-bit)

    --max_kmers NUM  - maximum number of kmers per node

    --silent         - do not output kmer counts (useful if only compiling db or index is needed

    --header         - print header row

    --total          - print the total number of kmers per node

    --unique         - print the number of nonzero kmers per node

    --kmers          - print individual kmer counts (default if no other output)

    --compile_index FILENAME - Add read index to database and write it to file

    --distribution NUM  - print kmer distribution (up to given number)

    --num_threads    - number of worker threads (default 24)

    --prefetch       - prefetch memory mapped files (faster on high-memory systems)

    --recover        - recover from FastA/FastQ errors (useful for corrupted streams)

    --stats          - print some statistics about sequence and kmers

    -D               - increase debug level

    -DDB             - increase database debug level

 パスを通しておく。

 

 

実行方法

チュートリアルに従ってバリアントの検出を行う。

 

 k-merデータベースのダウンロード。

#テキストバージョン
wget http://bioinfo.ut.ee/FastGT/downloads/kmer_list_WG30238282.db

#binaryバージョンはこちら
wget http://bioinfo.ut.ee/FastGT/downloads/kmer_db_WG30238282.dbb

 

NA12877の独立した３つのシーケンスデータをダウンロード。 

wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR194/ERR194146/ERR194146.fastq.gz 
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR194/ERR194146/ERR194146_1.fastq.gz 
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR194/ERR194146/ERR194146_2.fastq.gz
gunzip ERR194146*.fastq

 

genotypesをカウント。

gmer_counter -db kmer_list_WG30238282.db ERR194146*.fastq > ERR194146.counts

#ERR194146.fastqだけ使うなら
gmer_counter -db kmer_list_WG30238282.db ERR194146.fastq > ERR194146.counts
  • --num_threads      number of worker threads (default 24)
  • -db    SNP/KMER database file
  • -dbb     binary database file

 

genotypesをコール。テキストファイルが出力される。

gmer_caller ERR194146.counts > ERR194146.calls

先頭を開く。

head ERR194146.calls

$ head ERR194146.calls 

1:657698:rs565995692:C/T NC 0 0

1:658426:rs188842781:G/T NC 0 0

1:668346:rs115048193:G/A NC 2 0

1:668374:rs138476838:G/A NC 2 0

1:676127:rs150820983:C/T NC 2 0

1:693588:rs574459339:G/A NC 0 0

1:693747:rs532427839:A/G NC 0 0

1:697919:rs539728205:T/C NC 0 0

1:705475:rs564206378:G/A NC 0 0

1:706645:rs148085246:A/C NC 0 0

  •  カラム１: k-merデータベースのヒットしたマーカー（chr:position:rs_number:reference_allele/alternative_allele）
  • カラム２: genotype。AAはリファレンスと一致、Bは一致しない。
  • カラム３: genotypeの尤度
  • カラム４:リファレンスと一致したk-merの頻度
  • カラム５:リファレンスと一致しなかったk-merの頻度

 

リファレンスと一致しないならAB（ヘテロ）かBB（ホモ）がつくので、 これをgrepで抽出する (AA、BB、AB、上のNC (no call) がある )。

grep AB ERR194146.calls | head

AAやBBも同じ。 

 

A|Bだけでなく全ての遺伝子型をコールしたいなら以下のように打つ。

gmer_caller --non_canonical ERR194146.counts > ERR194146.all_calls

 

引用

FastGT: an alignment-free method for calling common SNVs directly from raw sequencing reads

Pajuste FD, Kaplinski L, Möls M, Puurand T, Lepamets M, Remm M.

Sci Rep. 2017 May 31;7(1):2537.