以前の記事でリファレンス配列についてだらだらと書かせて頂きました。
www.kimoton.com
手持ちのデータは最新のhg19でマッピングされているけれど論文の情報はhg38でのマッピングが行われている場合など、異なるアセンブリのデータを取得したい場合、本来はそのアセンブリのリファレンスを用いて再解析を行う必要があります。
一方で手持ちのデータのポジションを必要なアセンブリのポジションに変換することで対応できる場合も少なくありません。
このような処理を「リフトオーバー」と呼びます。
今回はbedファイル及びVCFファイルについてリフトオーバーの行い方を見ていきます。
公開されているスクリプトを使用すれば、ポジションのみならず、rs番号を変換することもできたりします。
ファイルフォーマット
まずは今回使用するファイルフォーマットについてざっくり説明していきます。
BEDファイル
ゲノムのポジションに関する情報は、大抵BEDというフォーマットによって表されます。 これは、Browser Extensible Dataの略で、恐らくUCSC Genome Browserでゲノム情報を可視化する際に使用したことが始まりだと思われます。
BEDファイルでは、以下のように1行ごとに位置情報 (chrom, chromStart, chromEnd) 及びそれに対応する情報が記載されています。
フォーマットの詳細はUCSCのFAQで書かれているので参照してください。
chr1 743267 743268 rs3115860 chr1 766408 766409 rs12124819 chr1 773885 773886 rs17160939
なお、3列のみ（すなわち位置情報のみ）のBEDファイルを3BED、
4列のみ（すなわち位置情報 + 情報1つ）のBEDファイルを4BEDと呼んだりします。
だからどうしたって感じですが。。
VCFファイル
Variant Call Format、即ち変異コールに関する情報を含めたフォーマットです。 samtoolsのページでVCFフォーマットに関する詳細なドキュメントが公開されているので時間がある方は見てみてください。
以下のようなフォーマットになっています。
「#」から始まるヘッダー行の下に、
CHROM、POS、ID、REF、ALT、QUAL、FILTER、INFO...と、変異に関する情報がつらつらと書かれています。
例えば1行目のデータに関して言えば、20番染色体の10001019番目のT（チミン）がG（グアニン）に変異したSNPを示しています。1塩基変異のみならず、SV（構造多型）なんかを表現することもできます。
[HEADER LINES] #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA12878 20 10001019 . T G 364.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=0.699;ClippingRankSum=0.00;DP=34;ExcessHet=3.0103;FS=3.064;MLEAC=1;MLEAF=0.500;MQ=42.48;MQRankSum=-3.219e+00;QD=11.05;ReadPosRankSum=-6.450e-01;SOR=0.537 GT:AD:DP:GQ:PL 0/1:18,15:33:99:393,0,480 20 10001298 . T A 884.77 . AC=2;AF=1.00;AN=2;DP=30;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=29.49;SOR=1.765 GT:AD:DP:GQ:PL 1/1:0,30:30:89:913,89,0 20 10001436 . A AAGGCT 1222.73 . AC=2;AF=1.00;AN=2;DP=29;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=25.36;SOR=0.836 GT:AD:DP:GQ:PL 1/1:0,28:28:84:1260,84,0 20 10001474 . C T 843.77 . AC=2;AF=1.00;AN=2;DP=27;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=31.25;SOR=1.302 GT:AD:DP:GQ:PL 1/1:0,27:27:81:872,81,0 20 10001617 . C A 493.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=1.63;ClippingRankSum=0.00;DP=38;ExcessHet=3.0103;FS=1.323;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.00;QD=12.99;ReadPosRankSum=0.170;SOR=1.179 GT:AD:DP:GQ:PL 0/1:19,19:38:99:522,0,480
GATK | Tool Documentation Indexより引用
chainファイル
リフトオーバーを行うには、対応するchainファイルと呼ばれるファイルが必要になります。
これはリファレンスポジションの対応表のようなもので、以下からダウンロードすることができます。
Index of /goldenPath/hg19/liftOver
hg19からhg38への変換に必要なchain fileは以下のコマンドからダウンロード可能です。
$ wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/liftOver/hg19ToHg38.over.chain.gz $ gunzip hg19ToHg38.over.chain.gz
中身はただのテキストファイルとなっています。
$ head hg19ToHg38.over.chain chain 20851231461 chr1 249250621 + 10000 249240621 chr1 248956422 + 10000 248946422 2 167376 50041 80290 40302 253649 288020 1044699 1 2 3716 0 3 1134 4 18 3377 0 1 7258 1 1 27 1 1 1275 1680 5595
リフトオーバー
BEDファイルの変換
BEDファイルの変換にはliftOverを使用します。
このツールはブラウザ上でも使うことができます。
ブラウザ版の場合、リフトオーバーを行いたいファイルをアップロード、若しくはその中身をペーストした後にSubmitボタンを押すだけです。
バイナリ版を使用する際は以下のコマンドによりダウンロードしてください。
$ wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/liftOver $ chmod +x ./liftOver
例えば以下のように使用します。
$ liftOver input.bed hg19ToHg38.over.chain output.bed unlifted.bed
unlifted.bedには変換先のリファレンスに対応するポジションが存在しなかったゲノムポジションが格納されます。
古いゲノムビルドで存在していたcontigが新しいビルドでは存在しなくなった場合、このファイルに格納されることになります。
LiftoverVcf
VCFのリフトオーバーには、LiftoverVcfを使用します。
Picardsのコマンドの1つであり、Javaで動作します。
Javaのインストールはここを参照してください。
Picardはここから最新版をダウンロード。
テストデータにはMITライセンスで公開されている以下のデータを使いましょう。dbSNP build138のVCFファイルみたいですね。
LiftoverVcfを使用する際には、まずdictファイルを作成する必要があります。
TARGET genomeのリファレンスファイルが置かれたディレクトリにCreateSequenceDictionary関数により.dictファイルを作成しましょう。
java -jar picard.jar CreateSequenceDictionary \ R=/home/kimoton/genome/hg38.fa \ O=/home/kimoton/genome/hg38.dict
お目当てのLiftoverVcf関数を実行します。
hg19リファレンス基準のVCFファイルを変換する際は以下のようにコマンドを指定します。
java -jar picard.jar LiftoverVcf \ I=dbSNP_138.hg19.vcf \ O=lifted_over.vcf \ CHAIN=hg19ToHg38.over.chain \ REJECT=rejected_variants.vcf \ R=/home/kimoton/genome/hg38.fa
bedファイルの変換時同様に、変換先のリファレンスに対応するポジションが存在しなかった場合は別ファイルとしてrejected_variants.vcfに格納されることになります。
