macでインフォマティクス

macでインフォマティクス

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

変異はどんな遺伝子に起きているのか? SnpEffを使ってindel検出結果のアノテーションを行う

 多くのindel検出ツールは変異のポジションしか出力しないため、その変異がどのようなアミノ酸変化を引き起こすか、サイレントなのかは別途調べる必要がある。ただし手動でやるのはしんどいし、間違いの元になる。snpEffはこうした作業をサポートするプログラムである。バリアント検出結果のvcfを入力として、データベースを元にアミノ酸変化などをコールしてくれる。SnpEffは2500以上のモデル生物のアノテーションデータを保持している。たいていの生物ですぐにアノテーションをできるのも広く使われている理由の1つと思われる。

 

8/24追記

 

ダウンロード

SnpEff

公式マニュアル

http://snpeff.sourceforge.net/SnpEff_manual.html#databases

最近brewでも導入可能になった。 

 

ダウンロードしたら解凍し、snpEff_latest_coresnpEff/内に入る。まずは、アノテーションをつけたい生物のデータベースが存在するかコマンドを打って確認する。

cd snpEff_latest_core/snpEff/
java -jar snpEff.jar databases | grep -i 7002|cut -f 1 # 下の*1を参照

ここではSynechococcus sp. PCC7002というシアノバクテリアを検索している(まずは7002だけで検索)。上記のコマンドを打つと

Campylobacter_jejuni_subsp_jejuni_icdccj07002 
Neisseria_meningitidis_70021
Streptococcus_pneumoniae_gca_001170025
Synechococcus_sp_pcc_7002
Treponema_medium_atcc_700293

 がヒットした。4つ目が今回の生物に該当する。

 

生物名が不明な場合、上のコマンドのパイプでつないだcutを"cut -f 1,2" に変えて実行する。データベースで説明されている"GRCh38.76"を検索すると

 >java -jar snpEff.jar databases | grep -i GRCh38|cut -f 1,2
GRCh38.86 Homo_sapiens
GRCh38.p7.RefSeq Human genome GRCh38 using RefSeq transcripts

人のデータとわかる。rnaも登録されていてアノテーションできるようだ。

 

またはsnpEff.configを開き、どんな生物が登録されているのか確認する。

snpEff databases html > list.html
open list.html

f:id:kazumaxneo:20170821205125j:plain

 

 

 

 

データベース名がわかったら、それをダウンロードする。

java -jar snpEff.jar download -v Synechococcus_sp_pcc_7002 -c /Users/user/local/snippy/etc/snpeff.config

-v: 進捗状況の表示。

-c: snpeff.configのパスを指定 (-cで指定していないと途中で止まってしまう)

snpeff.configの場所は各自検索してください。上記はbrewでインストールした場合です。

 

これでアノテーションする準備ができた。ランは以下のコマンドで行う。

java -Xmx4g -jar snpEff.jar Synechococcus_sp_pcc_7002 indel.vcf > annotation.vcf

-Xmx4g: javaのコマンド。メモリを4Gまでに制限。

indel.vcf: GATKなどの変異検出ソフトで出力されたファイル。

databasesコマンドで表示された名前をデータベース名として入力する。

  

入力のindelファイルはVariant Call Format(VCF)に対応している。GATK以外に多くのindel検出ソフトの出力はVCFフォーマットの規則に則っているので、そのままsnpEffで変異部位のアノテーションはできると思われる。VCFの詳細はsnpFff (http://snpeff.sourceforge.net/SnpEff_manual.html)

のVCF filesの項目を参照する。

 

 

 

データベースとリファレンスゲノムをセットで拾う例も紹介する。

例えばクラミドモナスなら

>snpEff databases | grep -i Chlamydomonas|cut -f 1,2
v3.1.29                                                   Chlamydomonas_reinhardtii 

ゲノムのversion3.1.29をアノテーションに使っていることがわかった。植物ゲノムが登録されているRSATデータベースのgenomeフォルダ内のここから.faファイルをダウンロードできる。

並行してsnpEffのデータベースもダウンロードする。

java -jar snpEff.jar download -v v3.1.29 -c /Users/user/local/snippy/etc/snpeff.config

ランが正常に終わると、まとめのhtmlファイルも作成される。 

f:id:kazumaxneo:20170824141212j:plain

f:id:kazumaxneo:20170824141056j:plain

 

 

f:id:kazumaxneo:20170824141147j:plain

f:id:kazumaxneo:20170824141155j:plain

f:id:kazumaxneo:20170824141207j:plain

*1 vcfファイルのコメント行を捨ててしまうと、htmlのサマリーはエラーになる。

 

*2 brewでインストールした人は最初を端折って"snpEff"だけで実行する。

 

*3 データベースとリファレンスゲノムの名前がマッチせずエラーになる場合は、data/にあるsnpeffのデータベースを確認し、それと一致する名前に変えればいい。例えば、data/Chlamydomonas_reinhardtii/snpEffectPredictor.binを解凍してsnpEffectPredictorを開くと、先頭に

user$ head /Users/user/local/snpEff_latest_core/snpEff/data/Chlamydomonas_reinhardtii/snpEffectPredictor 

SnpEff 4.3

CHROMOSOME 2 1 0 8222 DS496638 false

CHROMOSOME 3 1 0 8212 DS496637 false

CHROMOSOME 4 1 0 8189 DS496639 false

CHROMOSOME 5 1 0 8224 DS496634 false

CHROMOSOME 6 1 0 8224 DS496633 false

CHROMOSOME 7 1 0 8381 DS496636 false

CHROMOSOME 8 1 0 8219 DS496635 false

CHROMOSOME 9 1 0 8290 DS496630 false

CHROMOSOME 10 1 0 8229 DS496632 false

ヘッダーの部分に当たる名前が載っている。これと同じ名前をfastaに使う。例えば

CHROMOSOME 2なら、DS496638をfastaのヘッダーにする。

>DS496638

fastaを修正して変異解析をやり直しても良いが、vcfファイルの1カラム目を修正した方が早い。

 

 

 

参考

https://www.biostars.org/p/124885/