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

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

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

詳細は[Sequence Variant Nomenclature](http://varnomen.hgvs.org)を参照

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

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

Variant RecodeはEnsemblのVariant Effect Predictor [Ensembl/ensembl-vep](https://github.com/Ensembl/ensembl-vep)と共にgit上で公開されているツール

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

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

```terminal
$ docker pull ensemblorg/ensembl-vep 
```

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

```terminal
$ 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表記のデータが得られる

```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`とする

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

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

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

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

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

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

```text:refseq_nm.txt
NM_000603.4:c.894T>G
NM_000029.3:c.803T>C
```

```terminal
$ 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ディレクトリ内でツールを実行するよう指定している

結果は次の通り

```json
[
   {
      "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'を変換してみよう

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

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

```json
[
   {
      "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](https://www.ncbi.nlm.nih.gov/nuccore/NM_000603)の最新版はNM_000603.5なのだけれど、これはEnsemblのDBの現時点のそれとは異なっている

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

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

```text
-------------------- 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以外の情報も入っている

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

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

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

```terminal
$ 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: