macでインフォマティクス

macでインフォマティクス

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

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

2017 8/24追記

2019 4/24 インストール追記

 

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

 

ダウンロード

SnpEff

公式マニュアル

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

最近brewでも導入可能になった。 またcondaでも導入できる。

#anaconda環境を使っているなら
conda install -c bioconda -y snpeff

 > snpEff

$ snpEff

SnpEff version SnpEff 4.3t (build 2017-11-24 10:18), by Pablo Cingolani

Usage: snpEff [command] [options] [files]

 

Run 'java -jar snpEff.jar command' for help on each specific command

 

Available commands: 

[eff|ann]                    : Annotate variants / calculate effects (you can use either 'ann' or 'eff', they mean the same). Default: ann (no command or 'ann').

build                        : Build a SnpEff database.

buildNextProt                : Build a SnpEff for NextProt (using NextProt's XML files).

cds                          : Compare CDS sequences calculated form a SnpEff database to the one in a FASTA file. Used for checking databases correctness.

closest                      : Annotate the closest genomic region.

count                        : Count how many intervals (from a BAM, BED or VCF file) overlap with each genomic interval.

databases                    : Show currently available databases (from local config file).

download                     : Download a SnpEff database.

dump                         : Dump to STDOUT a SnpEff database (mostly used for debugging).

genes2bed                    : Create a bed file from a genes list.

len                          : Calculate total genomic length for each marker type.

pdb                          : Build interaction database (based on PDB data).

protein                      : Compare protein sequences calculated form a SnpEff database to the one in a FASTA file. Used for checking databases correctness.

seq                          : Show sequence (from command line) translation.

show                         : Show a text representation of genes or transcripts coordiantes, DNA sequence and protein sequence.

translocReport               : Create a translocations report (from VCF file).

 

Generic options:

-c , -config                 : Specify config file

-configOption name=value     : Override a config file option

-d , -debug                  : Debug mode (very verbose).

-dataDir <path>              : Override data_dir parameter from config file.

-download                    : Download a SnpEff database, if not available locally. Default: true

-nodownload                  : Do not download a SnpEff database, if not available locally.

-h , -help                   : Show this help and exit

-noLog                       : Do not report usage statistics to server

-t                           : Use multiple threads (implies '-noStats'). Default 'off'

-q , -quiet                  : Quiet mode (do not show any messages or errors)

-v , -verbose                : Verbose mode

-version                     : Show version number and exit

 

Database options:

-canon                       : Only use canonical transcripts.

-canonList <file>            : Only use canonical transcripts, replace some transcripts using the 'gene_id  transcript_id' entries in <file>.

-interaction                 : Annotate using inteactions (requires interaciton database). Default: true

-interval <file>             : Use a custom intervals in TXT/BED/BigBed/VCF/GFF file (you may use this option many times)

-maxTSL <TSL_number>         : Only use transcripts having Transcript Support Level lower than <TSL_number>.

-motif                       : Annotate using motifs (requires Motif database). Default: true

-nextProt                    : Annotate using NextProt (requires NextProt database).

-noGenome                    : Do not load any genomic database (e.g. annotate using custom files).

-noExpandIUB                 : Disable IUB code expansion in input variants

-noInteraction               : Disable inteaction annotations

-noMotif                     : Disable motif annotations.

-noNextProt                  : Disable NextProt annotations.

-onlyReg                     : Only use regulation tracks.

-onlyProtein                 : Only use protein coding transcripts. Default: false

-onlyTr <file.txt>           : Only use the transcripts in this file. Format: One transcript ID per line.

-reg <name>                  : Regulation track to use (this option can be used add several times).

-ss , -spliceSiteSize <int>  : Set size for splice sites (donor and acceptor) in bases. Default: 2

-spliceRegionExonSize <int>  : Set size for splice site region within exons. Default: 3 bases

-spliceRegionIntronMin <int> : Set minimum number of bases for splice site region within intron. Default: 3 bases

-spliceRegionIntronMax <int> : Set maximum number of bases for splice site region within intron. Default: 8 bases

-strict                      : Only use 'validated' transcripts (i.e. sequence has been checked). Default: false

-ud , -upDownStreamLen <int> : Set upstream downstream interval length (in bases)

 

 

データベース

ダウンロードしたら解凍し、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カラム目を修正した方が早い。

 

 

引用

A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3
Cingolani P, Platts A, Wang le L, Coon M, Nguyen T, Wang L, Land SJ, Lu X, Ruden DM.

Fly (Austin). 2012 Apr-Jun;6(2):80-92. doi: 10.4161/fly.19695.

 

 参考

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