HGVSとは
ゲノム等のリファレンス配列に対する違いを表現する表記としてHGVS(Human Genome Variation Society)が推奨している表記で、変異を表す際の標準となっている表記法
例えばcDNA上での変異であれば下記のように表現される
NM_000603.4:c.894T>G
詳細はSequence Variant Nomenclatureを参照
変異の表現方法の変換
cDNA上で表現された変異とその位置をゲノム上や、発現しているアミノ酸の位置を求めるにはどうすれば良いだろうか
ここではensemblで提供されているVariant Recoderを利用して変換を行ってみよう
Variant RecodeはEnsemblのVariant Effect Predictor Ensembl/ensembl-vepと共にgit上で公開されているツール
本体はPerlで記述されているので、インストールは比較的容易だけれど、Docker hub上でDockerイメージも提供されているので、ここはDockerを利用する事でインストールの手間を省こう
まずは、Docker イメージの取得
$ docker pull ensemblorg/ensembl-vep
Dockerのイメージが取得できたら後は実行してみよう
$ 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表記のデータが得られる
[
{
"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とする
$ docker run ensemblorg/ensembl-vep variant_recoder --input_data 'NM_000603.4:c.894T>G' --grch37 --fields hgvsg --pretty
これで、次の結果が得られる
[
{
"hgvsg" : [
"NC_000007.13:g.150696111T>G"
],
"input" : "NM_000603.4:c.894T>G"
}
]
複数の項目を指定する場合には、--fieldsオプションの後に必要なキーをカンマ区切りで指定すれば良い
$ docker run ensemblorg/ensembl-vep variant_recoder --input_data 'NM_000603.4:c.894T>G' --grch37 --fields hgvsg,hgvsc --pretty
複数の対象の変換
ここまでは一つの変異を対象に処理をおこなっていたが、対象が複数ある場合には、変異の一覧を作成して処理する方が効率が良い
変換する対象を一行ごとに区切って格納したファイルを作成し、--input-fileオプションで指定しよう
NM_000603.4:c.894T>G
NM_000029.3:c.803T>C
$ 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ディレクトリ内でツールを実行するよう指定している
結果は次の通り
[
{
"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'を変換してみよう
$ docker run ensemblorg/ensembl-vep variant_recoder --input_data 'NM_000603.5:c.894T>G' --grch37 --fields hgvsg --pretty
下記の様に、キーとしてwarningsが出力される
[
{
"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の最新版はNM_000603.5なのだけれど、これはEnsemblのDBの現時点のそれとは異なっている
バージョン番号を含めてがEnsemblのDBと同じでなければエラーになってしまう
*エラーが発生するとエラーログに Ensembl API version 番号が表示される
-------------------- 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以外の情報も入っている
$ cut -f 3-4 xref.txt | grep "^NP_" > np_map.tsv
$ cut -f 3-4 xref.txt | grep "^NM_" > nm_map.tsv
上記のような方法で作成したデータを利用してRefseqのIDの対応をとってから変換をおこなうと良いだろう
残念なところ
変換機能に関しては期待どおりに動作するけれど、変換に時間がかかる事は否めない
$ 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
今回はここまで