変異のHGVS表記の変換(Python hgvs)

# 変異のHGVS表記の相互変換

[Ensembl variant_recoder](https://qiita.com/percipere/items/d884aa2f72b1730f193d)を利用する方法は、変換時間がそれなりに必要なので、今回は[biocommons/hgvs](https://github.com/biocommons/hgvs) を利用して変異の記述の変換を行ってみよう

ここでもDockerを利用して作業を行う
hgvs の利用にはPython3.6以上が必要とされるため、dockerのイメージとしてpythonを利用する(記述時点ではPython 3.7.2)
python のイメージのダウンロード

```terminal
$ docker pull python
$ docker run -it python /bin/bash
```

##作業手順
この作業はpythonのコンテナ内部で実行する
主な手順は下記の通り

- docker内部でhgvs動作環境の構築
    - hgvsのインストール
    - seqrepoのインストール
    - Fasta配列及びDatabaseのダウンロード

```terminal
$ pip install hgvs
$ apt update
$ apt upgrade -y
$ apt install -y tabix rsync
$ pip install seqrepo
$ mkdir /usr/local/share/seqrepo
$ seqrepo pull
```

`seqrepo pull`でDB及びFastaデータのダウンロードが開始される
ダウンロード終了後のseqrepo 2018-11-26 のサイズはおよそ11GB

```terminal
$ du -sh /usr/local/share/seqrepo
11G	/usr/local/share/seqrepo
```

ダウンロードしたデータを確認する

```terminal
$ ls -d /usr/local/share/seqrepo/*
/usr/local/share/seqrepo/2018-11-26
```

**後でdocker コンテナを起動する際にこのディレクトリが必要なので書き留めておく**

## dockerイメージの作成
ここまでの手順を行ったdockerコンテナをイメージとして保存する

dockerコンテナをデタッチするためのキーストローク
Ctrl+p Ctrl+qを押す

ホストにカーソルが戻るので、現在のコンテナの状況を確認する

```terminal
$ docker ps -l
CONTAINER ID        IMAGE               COMMAND             CREATED             STATUS              PORTS               NAMES
c9342540d9b5        python              "/bin/bash"         23 hours ago        Up 23 hours                             modest_wilson
```

`CONTAINER ID`については動作毎に異なる事に注意

上記で確認したコンテナの停止

```terminal
$ docker stop c93
```

> `stop`の後に指定するIDは識別するのに必要な最低限の長さで良い

停止したdockerコンテナを利用してdockerイメージを作成する`c93`部分はこれまでと同じく`CONTAINER ID`なので、動作毎に異なっている
ここで取得した`hgvs`は1.2.5.post1だったので、名称を`hgvs:1.2.5.post1`とした
seqrepoのサイズが大きいのでイメージの作成には時間がかかる

```terminal
$ docker commit c93 hgvs:1.2.5.post1
```

```
$ docker images hgvs
REPOSITORY          TAG                 IMAGE ID            CREATED             SIZE
hgvs                1.2.5.post1         b82d7412ceeb        15 minutes ago      12.4GB
```

12.4GBという大きなサイズのコンテナが作成された

# hgvsを利用する

これまでの手順で準備ができたのでhgvsを利用してcdna配列の変異からゲノム上の変異に変換してみよう
次の様なプログラムをカレントディレクトリに準備しよう

```python3:convert_ctog.py
#!/usr/local/bin/python3
import sys
from hgvs import parser, variantmapper, assemblymapper
from hgvs.dataproviders import uta

target = 'NM_000603.4:c.894T>G'
if len(sys.argv) > 1:
    target = sys.argv[1]

#Set uta
hdp = uta.connect()

#prepare hgvs parser
parser = parser.Parser()
variant = parser.parse_hgvs_variant(target)
print(variant)

# prepare human genome mapper
mapper_grch37 = assemblymapper.AssemblyMapper(hdp, assembly_name='GRCh37', alt_aln_method='splign')

# map to genome
mapped = mapper_grch37.c_to_g(variant)
print(mapped)
```

コマンドの利用方法は下記の通り

```terminal
docker run -it \
    -u "$(id -u):$(id -g)" \
    -e HGVS_SEQREPO_DIR="seqrepoでダウンロードしたデータ及びシーケンスのディレクトリ" \
    -e UTA_DB_URL="UTAのpostgres SQLサーバーの接続先" \
    -v ${PWD}:/work \
    --workdir /work \
    hgvs:1.2.5.post1 \
    ./convert_ctog.py 
```


`NM_000603.4:c.894T>G`からGRCh37の変異の位置を求める例は次の通り

```terminal
$ docker run -it \
-u "$(id -u):$(id -g)" \
-e HGVS_SEQREPO_DIR="/usr/local/share/seqrepo/2018-11-26" \
-e UTA_DB_URL="postgresql://anonymous:anonymous@uta.biocommons.org/uta/uta_20180821" \
-v ${PWD}:/work \
--workdir /work \
hgvs:1.2.5.post1 ./convert_ctog.py `NM_000603.4:c.894T>G`
NM_000603.4:c.894T>G
NC_000007.13:g.150696111T>G
Closing connection; future mapping and validation will fail.
```

> `HGVS_SEQREPO_DIR`に設定するディレクトリは先程書き留めたものを利用する
> UTA_DB_URLについては、最新のUTA(Universal Transcript Archive)を設定
> [dl.biocommons.org](http://dl.biocommons.org/uta)を確認するとuta_20180821が最新
> UTA_DB_URL="postgresql://anonymous:anonymous@uta.biocommons.org/uta/uta_20180821"

`./convert_ctog.py` の後に引数を与える事もできる
'NM_182763.2:c.688+403C>T'を引数として与えてみる

```terminal
$ docker run -it \
-u "$(id -u):$(id -g)" \
-e HGVS_SEQREPO_DIR="/usr/local/share/seqrepo/2018-11-26" \
-e UTA_DB_URL="postgresql://anonymous:anonymous@uta.biocommons.org/uta/uta_20180821" \
-v ${PWD}:/work  \
--workdir /work \
hgvs:1.2.5.post1 ./convert_ctog.py 'NM_182763.2:c.688+403C>T'
NM_182763.2:c.688+403C>T
NC_000001.10:g.150550916G>A
Closing connection; future mapping and validation will fail.
```

> 与える文字列は必ずクオーテーションでくくること

[biocommons/hgvs](https://github.com/biocommons/hgvs)を利用して変異の変換が実行できた

なお、今回の実行時間は

```terminal
$ time docker run -it \
-u "$(id -u):$(id -g)" \
-e HGVS_SEQREPO_DIR="/usr/local/share/seqrepo/2018-11-26" \
-e UTA_DB_URL="postgresql://anonymous:anonymous@uta.biocommons.org/uta/uta_20180821" \
-v ${PWD}:/work  \
--workdir /work hgvs:1.2.5.post1 ./convert_ctog.py 'NM_182763.2:c.688+403C>T'
...省略...
real	0m8.473s
user	0m0.029s
sys	0m0.014s
```


今回はここまで:smile: