ReseqやRNA-seq他、多くのNGS解析において欠かせないのが、reference genome、参照ゲノム塩基配列です。
NGS解析のワークフローはおおよそこんな感じです。 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 というソフトウェアで変換することが可能です。
- webツール
Lift Genome Annotations - コマンドラインツール(非営利のアカデミアや個人目的での使用のみ無料で利用可能)
https://genome-store.ucsc.edu/
このツールはBEDファイルにのみ対応しているので、他のファイル形式のデータを変換したいときは、適用することができません。
VCFのゲノムビルドを変換する際はGATKのLiftoverVcfを使用することができます。
https://gatk.broadinstitute.org/hc/en-us/articles/360037060932-LiftoverVcf-Picard-
用意するもの
- GATK
- 本記事ではバージョン4.1.0.0を使用(現在の最新バージョンは4.2.2.0ですが、おそらく使い方は同じです)
- 変換したいVCF
- 本記事では、hg19を参照ゲノムに作成されたVCFをhg38に変換することにします
- chainファイル
変換するゲノムビルド同士のgenome coordinate対応に関する情報を含むファイルです。オリジナルのliftOverツールで使用するのと同じものを使用します- UCSC Genome BrowserのHPから、
Downloads > Genome Data > [変換元のゲノムビルド(ここではhg19)] > LiftOver files
変換元のゲノムビルドから、その他のゲノムに変換するためのchainファイルが用意されています。
本記事では、hg19からhg38に変換したいので、hg19ToHg38.over.chain.gz
をダウンロードします。
ところで、hg19ToMm10.over.chain.gz
のように他生物種に変換するためのchainファイルもたくさん用意されているのですが、そんな簡単に異なる生物種間のゲノムを変換できるんですか……すごいですね。
- UCSC Genome BrowserのHPから、
- 変換する先の参照ゲノムの塩基配列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