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

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

More than 1 year has passed since last update.

変異の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:

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