変異の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配列の変異からゲノム上の変異に変換してみよう
次の様なプログラムをカレントディレクトリに準備しよう
# !/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
今回はここまで