Edited at

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


変異のHGVS表記の相互変換

Ensembl variant_recoderを利用する方法は、変換時間がそれなりに必要なので、今回はbiocommons/hgvs を利用して変異の記述の変換を行ってみよう

ここでもDockerを利用して作業を行う

hgvs の利用にはPython3.6以上が必要とされるため、dockerのイメージとしてpythonを利用する(記述時点ではPython 3.7.2)

python のイメージのダウンロード

$ docker pull python

$ docker run -it python /bin/bash


作業手順

この作業はpythonのコンテナ内部で実行する

主な手順は下記の通り


  • docker内部でhgvs動作環境の構築


    • hgvsのインストール

    • seqrepoのインストール

    • Fasta配列及びDatabaseのダウンロード



$ 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

$ du -sh /usr/local/share/seqrepo

11G /usr/local/share/seqrepo

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

$ ls -d /usr/local/share/seqrepo/*

/usr/local/share/seqrepo/2018-11-26

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


dockerイメージの作成

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

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

Ctrl+p Ctrl+qを押す

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

$ 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については動作毎に異なる事に注意

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

$ docker stop c93


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


停止したdockerコンテナを利用してdockerイメージを作成するc93部分はこれまでと同じくCONTAINER IDなので、動作毎に異なっている

ここで取得したhgvsは1.2.5.post1だったので、名称をhgvs:1.2.5.post1とした

seqrepoのサイズが大きいのでイメージの作成には時間がかかる

$ 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配列の変異からゲノム上の変異に変換してみよう

次の様なプログラムをカレントディレクトリに準備しよう


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)


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

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の変異の位置を求める例は次の通り

$ 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を確認するとuta_20180821が最新

UTA_DB_URL="postgresql://anonymous:anonymous@uta.biocommons.org/uta/uta_20180821"


./convert_ctog.py の後に引数を与える事もできる

'NM_182763.2:c.688+403C>T'を引数として与えてみる

$ 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を利用して変異の変換が実行できた

なお、今回の実行時間は

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