HGVSとは

ゲノム等のリファレンス配列に対する違いを表現する表記としてHGVS(Human Genome Variation Society)が推奨している表記で、変異を表す際の標準となっている表記法

例えばcDNA上での変異であれば下記のように表現される
NM_000603.4:c.894T>G

詳細はSequence Variant Nomenclatureを参照

変異の表現方法の変換

cDNA上で表現された変異とその位置をゲノム上や、発現しているアミノ酸の位置を求めるにはどうすれば良いだろうか

ここではensemblで提供されているVariant Recoderを利用して変換を行ってみよう

Variant RecodeはEnsemblのVariant Effect Predictor Ensembl/ensembl-vepと共にgit上で公開されているツール

本体はPerlで記述されているので、インストールは比較的容易だけれど、Docker hub上でDockerイメージも提供されているので、ここはDockerを利用する事でインストールの手間を省こう

まずは、Docker イメージの取得

$ docker pull ensemblorg/ensembl-vep

Dockerのイメージが取得できたら後は実行してみよう

$ docker run ensemblorg/ensembl-vep variant_recoder --input_data "NM_000603.4:c.894T>G" --grch37  --pretty

オプション	説明	コメント
--input_data	変換元のhgvsで表記された文字列	ダブルクオーテーション又はシングルクオーテーションで区切る
--grch37	ヒトゲノムのバージョン指定	無指定であればhg38が利用される
--pretty	JSON表記を見やすく整形表示する

###変換結果
変換には若干の時間が必要だけれど、この例では次のようなJSON表記のデータが得られる

[
   {
      "hgvsc" : [
         "ENST00000297494.3:c.894T>G",
         "ENST00000461406.1:c.276T>G",
         "ENST00000467517.1:c.894T>G",
         "ENST00000484524.1:c.894T>G",
         "NM_000603.4:c.894T>G",
         "NM_001160109.1:c.894T>G",
         "NM_001160110.1:c.894T>G",
         "NM_001160111.1:c.894T>G"
      ],
      "input" : "NM_000603.4:c.894T>G",
      "id" : [
         "rs1799983",
         "CM981388"
      ],
      "hgvsp" : [
         "ENSP00000297494.3:p.Asp298Glu",
         "ENSP00000417143.1:p.Asp92Glu",
         "ENSP00000420551.1:p.Asp298Glu",
         "ENSP00000420215.1:p.Asp298Glu",
         "NP_000594.2:p.Asp298Glu",
         "NP_001153581.1:p.Asp298Glu",
         "NP_001153582.1:p.Asp298Glu",
         "NP_001153583.1:p.Asp298Glu"
      ],
      "hgvsg" : [
         "NC_000007.13:g.150696111T>G"
      ]
   }
]

戻り値のJSONのキーと値は次の通り

キー	意味
hgvsc	cDNA上の変異
hgvsp	アミノ酸配列上の変異
hgvsg	ゲノム配列上の変異
input	問い合わせした変異文字列
id	変異のID:dbSNPのrsIDやHGMDのAccession等

問い合わせの変異情報を表す文字列としては、ここで利用しているHGVS表記の他に、VCFの表記、VEPの表記、dbSNPのrsID、HGMDのAccession番号が利用できる

###項目の絞り込み
取得したい項目のみ指定して変換する事もできる
例えば、ゲノム上の位置のみが必要であれば、--fields hgvsgとする

$ docker run ensemblorg/ensembl-vep variant_recoder --input_data 'NM_000603.4:c.894T>G' --grch37 --fields hgvsg --pretty

これで、次の結果が得られる

[
   {
      "hgvsg" : [
         "NC_000007.13:g.150696111T>G"
      ],
      "input" : "NM_000603.4:c.894T>G"
   }
]

複数の項目を指定する場合には、--fieldsオプションの後に必要なキーをカンマ区切りで指定すれば良い

$ docker run ensemblorg/ensembl-vep variant_recoder --input_data 'NM_000603.4:c.894T>G' --grch37 --fields hgvsg,hgvsc --pretty

複数の対象の変換

ここまでは一つの変異を対象に処理をおこなっていたが、対象が複数ある場合には、変異の一覧を作成して処理する方が効率が良い
変換する対象を一行ごとに区切って格納したファイルを作成し、--input-file オプションで指定しよう

refseq_nm.txt

NM_000603.4:c.894T>G
NM_000029.3:c.803T>C

$ docker run -v${PWD}:/work -w=/work ensemblorg/ensembl-vep variant_recoder --input_file refseq_nm.txt --grch37 --pretty

オプション	説明
--input_file	変異の一覧を記入したファイルを指定

カレントディレクトリのrefseq_nm.txtファイルをdocker内で読み込むために、dockerに-v${PWD}:/work -w=/work というオプションを指定している
このオプションで、ホストのカレントディレクトリをdocker内部の/workディレクトリとしてマウントして、/workディレクトリ内でツールを実行するよう指定している

結果は次の通り

[
   {
      "hgvsp" : [
         "ENSP00000297494.3:p.Asp298Glu",
         "ENSP00000417143.1:p.Asp92Glu",
         "ENSP00000420551.1:p.Asp298Glu",
         "ENSP00000420215.1:p.Asp298Glu",
         "NP_000594.2:p.Asp298Glu",
         "NP_001153581.1:p.Asp298Glu",
         "NP_001153582.1:p.Asp298Glu",
         "NP_001153583.1:p.Asp298Glu"
      ],
      "hgvsg" : [
         "NC_000007.13:g.150696111T>G"
      ],
      "hgvsc" : [
         "ENST00000297494.3:c.894T>G",
         "ENST00000461406.1:c.276T>G",
         "ENST00000467517.1:c.894T>G",
         "ENST00000484524.1:c.894T>G",
         "NM_000603.4:c.894T>G",
         "NM_001160109.1:c.894T>G",
         "NM_001160110.1:c.894T>G",
         "NM_001160111.1:c.894T>G"
      ],
      "input" : "NM_000603.4:c.894T>G",
      "id" : [
         "rs1799983",
         "CM981388"
      ]
   },
   {
      "input" : "NM_000029.3:c.803T>C",
      "hgvsc" : [
         "ENST00000366667.4:c.803T>C",
         "NM_000029.3:c.803T>C"
      ],
      "id" : [
         "rs699",
         "CM920010"
      ],
      "hgvsg" : [
         "NC_000001.10:g.230845794A>G"
      ],
      "hgvsp" : [
         "ENSP00000355627.4:p.Met268Thr",
         "NP_000020.1:p.Met268Thr"
      ]
   }
]

##変換時に何かの問題が発生した場合

例えば'NM_000603.5:c.894T>G'を変換してみよう

$ docker run ensemblorg/ensembl-vep variant_recoder --input_data 'NM_000603.5:c.894T>G' --grch37 --fields hgvsg --pretty

下記の様に、キーとしてwarningsが出力される

[
   {
      "warnings" : [
         "Unable to parse HGVS notation 'NM_000603.5:c.894T>G': Could not get a Transcript object for 'NM_000603.5': Could not get a Transcript object for 'NM_000603.5'"
      ]
   }
]

Refseq のバージョン番号の解決

上記の通り、RefseqのIDのバージョン番号が利用しているEnsebmlのDBのそれと一致していないと変換が実行できない
例えば、これを記述している時点でのNM_000603の最新版はNM_000603.5なのだけれど、これはEnsemblのDBの現時点のそれとは異なっている

バージョン番号を含めてがEnsemblのDBと同じでなければエラーになってしまう

*エラーが発生するとエラーログに　Ensembl API version 番号が表示される

-------------------- EXCEPTION --------------------
(省略)
Ensembl API version = 95
---------------------------------------------------

Ensembl内のRefseqのバージョン番号を取得する

この問題を解決するためにEnsemblのref.txtファイルをftpで取得する(それぞれのサイズは28MB程)

Ensemblのサイトで最新のバージョンを確認してほしい
これを記述している時点での最新版はrelease-95

release-95のヒトゲノムのリファレンスGRCh37の場合のEnsemblのファイルは下記
ftp://ftp.ensembl.org/pub/grch37/release-95/mysql/homo_sapiens_core_95_37/xref.txt.gz.bz2

release-95のヒトゲノムのリファレンスGRCh38の場合のEnsemblのファイルは下記
ftp://ftp.ensembl.org/pub/release-95/mysql/homo_sapiens_core_95_38/ref.txt.gz

これらのファイルを展開して下記のカラムを抽出する

カラム	コメント
カラム3	Version番号を含まないRefseqのID
カラム4	Version番号を含むRefseqのID (Version番号無しも有る)

実際のところここにはRefseq以外の情報も入っている

$ cut -f 3-4 xref.txt | grep "^NP_" > np_map.tsv
$ cut -f 3-4 xref.txt | grep "^NM_" > nm_map.tsv

上記のような方法で作成したデータを利用してRefseqのIDの対応をとってから変換をおこなうと良いだろう

残念なところ

変換機能に関しては期待どおりに動作するけれど、変換に時間がかかる事は否めない

$ time docker run ensemblorg/ensembl-vep variant_recoder --input_data 'NP_000594.2:p.Asp298Glu' --grch37 --pretty 
...(省略)
real	2m14.798s
user	0m0.032s
sys	0m0.015s

今回はここまで