VCF ファイルの解析とアノテーション

このページでは、次の方法について説明します。

入力 VCF ファイルのアノテーションフィールドを解析します。
オープンソースのアノテーションツールと Variant Effect Predictor（VEP）などのデータベースを使用して、インポート中の VCF ファイルにアノテーションを付けます。

アノテーションフィールドの解析

Variant Transforms ツールは、VCF 形式のバリアントアノテーションで指定されたアノテーションフィールド形式をサポートしています。VEP はこれと同じ形式を使用します。次の例は、VEP によってアノテーションが付けられたファイルの VCF ヘッダーを示しています。

##INFO=<ID=CSQ,Number=.,Type=String,Description="Consequence annotations from Ensembl VEP. Format: Allele|Consequence|IMPACT|SYMBOL|Gene|Feature_type|Feature|BIOTYPE|EXON|INTRON|ALLELE_NUM|...>

REF=AAC と ALT=TAC,A のサンプルバリアントレコードの対応するアノテーションは、次のように示すことができます。

CSQ=-|downstream_gene_variant|MODIFIER|WASH7P|ENSG00000227232|Transcript|ENST00000423562|unprocessed_pseudogene|||2,...,T|upstream_gene_variant|MODIFIER|DDX11L1|ENSG00000223972|Transcript|ENST00000518655|transcribed_unprocessed_pseudogene|||1

カンマで区切られた複数のアノテーションリストが、バリアントレコードごとに指定されます。各セットの最初のフィールドは、代替対立遺伝子（ALT）に対応しています。上記の例では、ALT ごとにサブセットのみを使用しています。実際のファイルには、より多くのフィールドがあり、それぞれの代替対立遺伝子は異なるアノテーションで複数回繰り返すことができます。

Variant Transforms ツールを実行する際に、--annotation_fields を渡すことができます。このツールは、アノテーションリストを分割し、対応する alternate_bases レコードの配下に移動します。たとえば、以下の表ではレコードフィールド CSQ が各 alternate_bases 下に繰り返し見られることを確認できます。

名前	型	モード	説明
alternate_bases.CSQ	RECORD	REPEATED	この代替の CSQ アノテーションのリスト。
alternate_bases.CSQ.allele	STRING	NULLABLE	アノテーションフィールドの ALT 部分。
alternate_bases.CSQ.Consequence	STRING	Nullable	なし
alternate_bases.CSQ.IMPACT	STRING	Nullable	なし
alternate_bases.CSQ.SYMBOL	STRING	Nullable	なし
alternate_bases.CSQ.Gene	STRING	Nullable	なし
alternate_bases.CSQ.Feature_type	STRING	Nullable	なし
alternate_bases.CSQ.Feature	STRING	Nullable	なし
alternate_bases.CSQ.BIOTYPE	STRING	Nullable	なし
alternate_bases.CSQ.EXON	STRING	Nullable	なし
alternate_bases.CSQ.INTRON	STRING	Nullable	なし

代替対立遺伝子の一致

上記の例では、CSQ 列に表示されている代替対立遺伝子はバリアントレコードの ALT 文字列とは異なります。たとえば、AAC->C バリアントは「-」で示され、削除を表します。これは、VEP が --minimal フラグを持つデータセットで実行されるためです。Variant Transforms は、この minimal モードでの ALT フィールドの照合をサポートしています。

ALT 照合の次のモードがサポートされています。

デフォルトでは、ALT 一致フラグが設定されていない場合は、VCF 形式のバリアントアノテーションが使用されます。

注: がんサンプルはサポートされていませんが、サポートが開発中です。
上記の例で示すように、--use_allele_num フラグを使用すると、Variant Transforms ツールは ALLELE_NUM アノテーションフィールドの値を検索します。これは VEP の機能です。
--minimal_vep_alt_matching フラグを使用して、--minimal フラグを Variant Transforms ツールでシミュレートできます。ただし、一致があいまいな場合もあります。あいまいな一致は、Dataflow Console でカウントされ、報告されます。作成されたテーブルの新しいフィールドの配下にも追加されます。

--minimal を使用すると（--minimal_vep_alt_matching ではなく）最適な結果が得られます。

VEP を使用して VCF ファイルに自動的にアノテーションを付ける

通常、アノテーションフィールドを解析して個別のフィールドに分割するには、まず VCF ファイルにアノテーションを追加する必要があります。ただし、Variant Transforms ツールを使用する際、VCF ファイルに事前にアノテーションを追加する必要はありません。

Variant Transforms は、VCF ファイルが BigQuery に読み込まれる際に、自動的にアノテーションを付けることができます。そのためには、Variant Transforms ツールを使用して VEP を実行する必要があります。その場合は、事前に構成された Docker イメージを含む Compute Engine 仮想マシンインスタンスを作成し実行する必要があります。

Variant Transforms ツールを使用した VEP の実行

Variant Transforms ツールを VEP で実行し、BigQuery に読み込まれる際に VCF ファイルに自動的にアノテーションを付けるには、少なくとも次のフラグを使用して Variant Transforms ツールを実行する必要があります。

--run_annotation_pipeline
--max_num_workers MAX_NUM_WORKERS または --num_workers NUM_WORKERS を指定します。
--annotation_output_dir gs://CLOUD_STORAGE_BUCKET

CLOUD_STORAGE_BUCKET 変数に対して選択したバケットは、プロジェクトで所有している必要があります。

VEP の実行時に、Variant Transforms ツールに渡すことができるフラグの全一覧は、次のとおりです。

フラグ	デフォルト	説明
`--run_annotation_pipeline`	偽	自動アノテーションを有効にします。
`--max_num_workers` または `--num_workers`	自動的に決定	それぞれ、Cloud Dataflow が自動スケーリングするジョブのワーカーの最大数、または最初にジョブに割り当てられたワーカーの数を示します。
`--annotation_output_dir`	なし	VEP 出力ファイル用の Cloud Storage パス。入力ファイルの階層はこの場所に複製されますが、新しいファイルに `_vep_output.vcf` が追加されます。ディレクトリがすでに存在する場合、Variant Transforms ツールは失敗します。
`--vep_image_uri`	`gcr.io/gcp-variant-annotation/vep_91`	Cloud Genomics VEP Docker イメージ。Google では、デフォルト画像を保持しています（VEP バージョン 91）。
`--vep_cache_path`	`gs://gcp-variant-annotation-vep-cache/vep_cache_homo_sapiens_GRCh38_91.tar.gz`（GRCh38 参照配列に類似したヒトゲノムに推奨）	VEP キャッシュの圧縮バージョンの Cloud Storage パス。このファイルは、`--build_vep_cache.sh` スクリプトを使用して作成できます。
`--vep_info_field`	`CSQ_VT`	新規のアノテーションを保持する新しい INFO フィールドの名前。
`--vep_num_fork`	2	1 つのファイルに対して VEP を実行する際に使用される、ローカルのプロセス数。詳細については、VEP ドキュメントの `--fork [num_forks]` フラグをご覧ください。
`--shard_variants`	True	VEP を使用して VCF ファイルにアノテーションを付ける前に、入力ファイルをより小さい一時 VCF ファイルに分割するかどうかを指定します。入力ファイルが小さい場合（たとえば各 VCF ファイルのバリアントが 50,000 個未満の場合など）、このフラグを `True` に設定すると、コンピューティングリソースが浪費される可能性があります。
`--number_of_variants_per_shard`	2,000	`shard_variants` が `True` の場合に、各シャードに書き込まれるバリアントの最大数。サンプル数が多いデータセットの場合は、小さい値を使用することもできます。小さい値を指定すると、パイプラインの終了に時間がかかることがあります。ほとんどの場合に、デフォルト値が推奨されます。

トラブルシューティング

VEP を Variant Transforms ツールで実行し、ツールが失敗した場合は、プロセスが実行された仮想マシンのログを示す logs ディレクトリが作成されます。logs ディレクトリは、--annotation_output_dir フラグに使用するパスで指定します。これらのログファイルと、標準の Variant Transforms ツールのログを調べて、問題の原因を特定できます。

Variant Transforms ツールを使用して VEP を実行すると、VEP フラグ --check_ref が自動的に有効になります。そのため、logs ディレクトリのファイルには、参照と一致しないバリアントのリストが含まれます。正しい VEP キャッシュを使用していれば、不一致は自動的に破棄され、BigQuery の出力表には表示されないため、不一致は表示されません。

サンプルクエリ

影響の大きいバリアントの抽出

次のサンプルは、ヘモグロビンの β（ベータ）鎖をコードする HBB 遺伝子の影響の大きいバリアントを抽出する方法を示しています。

#standardSQL
SELECT reference_name, start_position, reference_bases, ALT.alt, CSQ.*
FROM vcf_imports_external.gnomad_genomes_chr_hg19 AS T, T.alternate_bases AS ALT, ALT.CSQ AS CSQ
WHERE CSQ.SYMBOL = "HBB" AND CSQ.IMPACT = "HIGH"
ORDER BY start_position

このクエリを実行すると、次の結果が返されます。

Row	reference_name	start_position	reference_bases	alt	allele	結果	影響	記号	遺伝子	Feature_type	機能	BIOTYPE	Exon	INTRON
1	11	5246947	G	GC	C	frameshift_variant	高	HBB	ENSG00000244734	文字起こし	ENST00000335295	protein_coding	3/3
2	11	5246957	T	C	C	splice_acceptor_variant	高	HBB	ENSG00000244734	文字起こし	ENST00000335295	proten_coding		2/2
3	11	5246957	T	C	C	splice_acceptor_variant&non_coding_transcript_variant	高	HBB	ENSG00000244734	文字起こし	ENST00000475226	retained_intron		1/1
4	11	5247805	C	T	T	splice_donor_variant	高	HBB	ENSG00000244734	文字起こし	ENST00000335295	protein_coding		2/2
5	11	5247805	C	T	T	splice_donor_variant&non_coding_transcript_variant	高	HBB	ENSG00000244734	文字起こし	ENST00000475226	retained_intron		1/1
6	11	5247991	CAAAG	C	-	frameshift_variant	高	HBB	ENSG00000244734	文字起こし	ENST00000335295	protein_coding	2/3
7	11	5247991	CAAAG	C	-	frameshift_variant	高	HBB	ENSG00000244734	文字起こし	ENST00000380315	protein_coding	4/4
8	11	5248003	G	A	A	stop_gained	高	HBB	ENSG00000244734	文字起こし	ENST00000335295	protein_coding	2/3
9	11	5248003	G	A	A	stop_gained	高	HBB	ENSG00000244734	文字起こし	ENST00000380315	protein_coding	4/4
10	11	5248199	T	A	A	stop_gained	高	HBB	ENSG00000244734	文字起こし	ENST00000335295	protein_coding	1/3
11	11	5248199	T	A	A	stop_gained	高	HBB	ENSG00000244734	文字起こし	ENST00000380315	protein_coding	3/4

gnomAD データセットには呼び出しまたはサンプルはありませんが、このクエリを拡張して影響の大きいバリアントがあるサンプルを選択できます。

公開 BigQuery データベースとのデータ結合

次のサンプルは、BigQuery で使用可能な pathway データベースなどの他のデータベースとデータを結合する方法を示しています。ISB-CGC では、たとえば次のようなデータベースを使用できます。

Reactome（BigQuery テーブル）
WikiPathways（BigQuery テーブル）

次のサンプルは、GO（Gene Ontology）生物学的プロセス（GO：0006302）で定義される二本鎖修復経路における高影響バリアントの抽出方法を示しています。

#standardSQL
SELECT
  reference_name, start_position, reference_bases, ALT.alt,
  CSQ.Consequence, CSQ.Impact, CSQ.SYMBOL
FROM
  `vcf_imports_external.gnomad_genomes_chr_hg19` AS T,
  T.alternate_bases AS ALT, ALT.CSQ AS CSQ
WHERE
  # Note: Matching based on symbol is "best effort" as the names may not be
  # standardized across sources.
  CSQ.SYMBOL IN (SELECT DB_Object_Symbol
                 FROM `isb-cgc.genome_reference.GO_Annotations`
                 WHERE GO_ID = 'GO:0006302')
  AND CSQ.IMPACT = "HIGH"
ORDER BY
  start_position