Help us understand the problem. What is going on with this article?

変異のHGVS表記の相互変換(Ensembl variant_recoder)

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

今回はここまで:smile:

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away