解析 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 记录下。例如，以下列表中每个 alternate_bases 下都有一个重复的记录字段 CSQ：

Name Type Mode 说明
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 的一项功能。

  • 您可以在 Variant Transforms 工具中使用 --minimal_vep_alt_matching 标志来模拟 --minimal 标志。不过，这可能导致不确定的匹配。这些不确定的匹配会在 Dataflow 控制台中进行计数和报告，并会添加到创建表格的新字段下。

    为获得最佳效果，请使用 --minimal，而不是 --minimal_vep_alt_matching。

使用 VEP 为 VCF 文件自动添加注释

通常，如果要解析注释字段并将其拆分为单独的字段，首先必须为 VCF 文件预先添加注释。不过，如果使用 Variant Transforms 工具，则无需为 VCF 文件预先添加注释。

Variant Transforms 可在将 VCF 文件加载到 BigQuery 中时自动为其添加注释。为此，您必须同时运行 Variant Transforms 工具和 VEP，这需要创建并运行包含预配置 Docker 映像的 Compute Engine 虚拟机实例。

同时运行 VEP 和 Variant Transforms 工具

如需同时运行 Variant Transforms 工具和 VEP，并在将 VCF 文件加载到 BigQuery 中时自动为其添加注释，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 False 启用自动注释。
--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 对单个文件运行 VEP 时使用的本地进程数。如需了解详情，请参阅 VEP 文档中的 --fork [num_forks] 标志。
--shard_variants 是 这用于确定在使用 VEP 向 VCF 文件添加注释之前是否将输入文件分片为较小的临时 VCF 文件。如果输入文件较小（例如，如果每个 VCF 文件包含的变体少于 50000 个），那么，将此标志设置为 True 可能会浪费计算资源。
--number_of_variants_per_shard 2000 这是 shard_variants 为 True 时写入每个分片的变体数量上限。如果您的数据集包含许多样本，不妨使用较小的值。如果提供较小的值，流水线可能需要更长的时间才能完成。在大多数情况下，建议使用该默认值。

问题排查

如果同时运行 VEP 和 Variant Transforms 工具，而该工具发生故障，系统会创建一个 logs 目录，提供来自运行进程的虚拟机的日志。logs 目录在 --annotation_output_dir 标志所用的路径中指定。除了标准的 Variant Transforms 工具日志之外，您还可以检查这些日志文件，以确定问题的原因。

当您同时运行 VEP 和 Variant Transforms 工具时，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

运行查询将返回以下内容：

行 reference_name start_position reference_bases alt allele Consequence IMPACT SYMBOL Gene Feature_type Feature BIOTYPE EXON INTRON
1 11 5246947 G GC C frameshift_variant HIGH HBB ENSG00000244734 Transcript ENST00000335295 protein_coding 3/3
2 11 5246957 T C C splice_acceptor_variant HIGH HBB ENSG00000244734 Transcript ENST00000335295 proten_coding 2/2
3 11 5246957 T C C splice_acceptor_variant&non_coding_transcript_variant HIGH HBB ENSG00000244734 Transcript ENST00000475226 retained_intron 1/1
4 11 5247805 C T T splice_donor_variant HIGH HBB ENSG00000244734 Transcript ENST00000335295 protein_coding 2/2
5 11 5247805 C T T splice_donor_variant&non_coding_transcript_variant HIGH HBB ENSG00000244734 Transcript ENST00000475226 retained_intron 1/1
6 11 5247991 CAAAG C - frameshift_variant HIGH HBB ENSG00000244734 Transcript ENST00000335295 protein_coding 2/3
7 11 5247991 CAAAG C - frameshift_variant HIGH HBB ENSG00000244734 Transcript ENST00000380315 protein_coding 4/4
8 11 5248003 G A A stop_gained HIGH HBB ENSG00000244734 Transcript ENST00000335295 protein_coding 2/3
9 11 5248003 G A A stop_gained HIGH HBB ENSG00000244734 Transcript ENST00000380315 protein_coding 4/4
10 11 5248199 T A A stop_gained HIGH HBB ENSG00000244734 Transcript ENST00000335295 protein_coding 1/3
11 11 5248199 T A A stop_gained HIGH HBB ENSG00000244734 Transcript ENST00000380315 protein_coding 3/4

gnomAD 数据集内没有任何检出或样本，但您可以扩展此查询，使选择的样本包含有较大影响的变体。

将数据与 BigQuery 公开数据库联接

以下示例展示了如何将数据与 BigQuery 中提供的其他数据库（如途径数据库）联接。可通过 ISB-CGC 获得的部分数据库包括：

以下示例展示了如何根据 GO（基因本体论）生物过程 (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

运行查询将返回以下内容：

行 reference_name start_position reference_bases alt Consequence IMPACT SYMBOL
1 5 917449 G T splice_acceptor_variant&non_coding_transcript_variant HIGH TRIP13
2 12 1022568 A C stop_gained HIGH RAD52
3 12 1022568 A C stop_gained HIGH RAD52
4 12 1022568 A C stop_gained HIGH RAD52
5 12 1023125 G A stop_gained HIGH RAD52
6 12 1023125 G A stop_gained HIGH RAD52
7 12 1023125 G A stop_gained HIGH RAD52
8 12 1023167 G A stop_gained HIGH RAD52
9 12 1023167 G A stop_gained HIGH RAD52
10 12 1023167 G A stop_gained HIGH RAD52
11 12 1023217 G T stop_gained HIGH RAD52
12 12 1023217 G T stop_gained HIGH RAD52
13 12 1023217 G T stop_gained HIGH RAD52
14 12 1025654 TGA T frameshift_variant HIGH RAD52
15 12 1025699 T G splice_acceptor_variant&non_coding_transcript_variant HIGH RAD52
16 12 1036004 GC G frameshift_variant HIGH RAD52
17 12 1038977 C CT frameshift_variant HIGH RAD52