アメリエフの技術ブログ

Amelieff Staff Blog

LiftoverVcfでVCFのゲノムビルドを変換する

ReseqRNA-seq他、多くのNGS解析において欠かせないのが、reference genome、参照ゲノム塩基配列です。

NGS解析のワークフローはおおよそこんな感じです。

f:id:kubo-m:20211006114008p:plain
解析フロー
https://biosciencedbc.jp/gadget/human/20160725_amelieff_20160803.pdf

まずはシークエンスデータをゲノム参照配列にマッピングしてから、変異や発現量定量などを行うという流れは、どんな解析でも共通しています。
このマッピング以降の解析では、どのゲノムビルドバージョンで作成されたデータを扱っているのかは非常に重要な情報になります。

生物のゲノム配列は、一度公開された後様々な修正や更新を経て、バージョンアップします。各バージョンのことをゲノムビルドバージョンと呼びます。

同じ生物種のデータの同じ変異や遺伝子、ゲノム領域でも、異なるゲノムビルド間ではgenome coordinateが全く異なり、そのままではgenome coordinateを用いた比較はできません。たとえばヒトの最新のゲノムビルドバージョンはGRCh38/hg38であり、そのひとつ前はGRCh37/hg19でしたが、同じ遺伝子についてもそれぞれでまったく異なるgenome coordinateで表現され、例としてAPOB遺伝子は

  • hg19:chr2:21,224,301-21,266,945
  • hg38:hr2:21,001,429-21,044,073

となります。
そのため、異なるゲノムビルドで作成されたデータを統合したり比較したい場合は、どちらかのgenome coordinateをもう片方に合わせて変換する必要があります。

liftOverでgenome coordinateの変換

BEDファイルのゲノムビルドは liftOver というソフトウェアで変換することが可能です。

このツールはBEDファイルにのみ対応しているので、他のファイル形式のデータを変換したいときは、適用することができません。
VCFのゲノムビルドを変換する際はGATKのLiftoverVcfを使用することができます。

https://gatk.broadinstitute.org/hc/en-us/articles/360037060932-LiftoverVcf-Picard-

用意するもの

  1. GATK
    • 本記事ではバージョン4.1.0.0を使用(現在の最新バージョンは4.2.2.0ですが、おそらく使い方は同じです)
  2. 変換したいVCF
    • 本記事では、hg19を参照ゲノムに作成されたVCFをhg38に変換することにします
  3. chainファイル
    変換するゲノムビルド同士のgenome coordinate対応に関する情報を含むファイルです。オリジナルのliftOverツールで使用するのと同じものを使用します
    • UCSC Genome BrowserのHPから、
      Downloads > Genome Data > [変換元のゲノムビルド(ここではhg19)] > LiftOver files
      f:id:kubo-m:20211007113634p:plain
      変換元のゲノムビルドから、その他のゲノムに変換するためのchainファイルが用意されています。
      本記事では、hg19からhg38に変換したいので、hg19ToHg38.over.chain.gzをダウンロードします。
      ところで、hg19ToMm10.over.chain.gz のように他生物種に変換するためのchainファイルもたくさん用意されているのですが、そんな簡単に異なる生物種間のゲノムを変換できるんですか……すごいですね。
  4. 変換する先の参照ゲノムの塩基配列fastaファイル ここではhg19→hg38に変換するのでhg38のfastaを用意します。

実行

gatk LiftoverVcf -I input.vcf -O output.vcf --CHAIN hg19ToHg38.over.chain.gz --REJECT rejected_variants.vcf -R /path/to/hg38/genome.fa

-Oで指定したファイル、ここではoutput.vcfがhg38に変換された結果です。

  • liftOverでは、変異のgenome coordinateを100%変換できるわけではありません。何らかの理由で変換できなかった変異は --REJECT で指定したファイルに書き出されます。変換したら目的の変異がなくなった!!などと慌てる前にrejected_variants.vcfを確認しましょう(自戒)。
  • 変換できない理由 → https://genome.sph.umich.edu/wiki/LiftOver#Various_reasons_that_lift_over_could_fail