使用 ANNOVAR 之前，你应该知道

2020 年 3 月 30 日
筆記

翻译自 ANNOVAR 文档：https://doc-openbio.readthedocs.io/projects/annovar/en/latest/articles/VCF/

本文介绍了有关 VCF 预处理的一些想法，以确保VCF文件的功能更加准确/可靠。

VCF 文件简介及其复杂性

最初开发 ANNOVAR 时，几乎所有 call 突变的软件都有自己的一套输出格式（SamTools，SOAPSNP，SOLiD BioScope，Illumina CASAVA，CG ASM-var，CG ASM-masterVAR 等），因此 ANNOVAR 就决定采用一种最简单的格式（仅包含 chr, start, end, ref, alt 以及 optional fields）作为输入。现将其称为 avinput 文件。我们也在 ANNOVAR 软件包中提供了 convert2annovar.pl 程序，以方便进行格式转换。

后来，VCF 逐渐成为描述突变的主流格式。有关其格式规范的详细信息，请参见此处[1]。

如今，几乎每个进行突变分析的研究人员都会使用 VCF 或 MAF 格式，这极大地促进了研究人员之间的交流和沟通。ANNOVAR 也提供了从 VCF/MAF 转换 avinput 格式的功能，以便用户可以注释 VCF 文件。

但许多用户可能不完全了解什么是 VCF 文件，因此，我收到了许多用户 email 来询问有关这方面的问题。与一些 ANNOVAR 用户进行交流之后，我决定写一个简单的文章来描述 VCF 文件的一般准则，以帮助 ANNOVAR 用户了解如何更好地注释 VCF 文件。

首先你需要知道的一些基本事实：

1.VCF 是一种用于描述基因座的格式；从技术上讲，尽管其名称为“ Variant Call Format”，但它并不用于描述突变或基因型。也可以仅包含基因型（或突变）突变信息，对于许多非二倍体物种或许多其他情况（例如线粒体或人类癌症），有时甚至没有意义进行基因型注释。虽然大多数软件的目的是生成基因型信息，它们使用 VCF 作为输出文件的格式，但这并不意味着 VCF 就是被设计用于存储突变信息。2.由于 VCF 是 locus descriptor，这样就会产生许多后果。首先，突变和位置并不是一一对应的。多个突变可以位于同一基因座中，因此在有基因型信息时，VCF文件中的一行，原则上可以描述多个突变（包括野生型非突变等位基因）和多种类型的基因型。例如，看下面的示例 VCF 记录。它包含八个制表符分隔的列。在 ALT 列中，有几个逗号分隔的替代等位基因。因此，在一行中，会同时存在数个插入和缺失以及一个单核苷酸突变（SNV）。

1 112240038 . CTTT CTTTT,CTTTTT,CTTA,CTT,CT,C . PASS AC=986,3,1249,3,127,3;AF=0.196885,0.000599042,0.249401,0.000599042,0.0253594,0.000599042

许多用户喜欢在 VCF 文件中包含突变的注释信息（ INFO 字段）。因此，在上述情况下，我们需要在同一行的 INFO 字段中为这六个等位基因都添加注释，并确保用户知道哪个注释对应于哪个等位基因。（在 ANNOVAR 中可通过 table_annovar.pl 正确处理此类问题）

1.VCF 可能会劫持你的突变，将 SNV 变成多核苷酸突变，并将简单的 indel 变成复杂的描述符。这将对注释产生挑战，因为理想情况下，一个突变在特定参考基因组中有且仅有一种注释信息。以上面的 CTTT->CTTA 的突变为例，它应该只是 T->A 的单碱基突变，但 indel 劫持了基因座，因此将其写为 CTTT->CTTA，而不是 T->A。考虑到一些数据库（例如千人基因组计划数据库）仅包含 T->A 突变，而没有CTTT->CTTA，所以注释软件就会忽略这种突变，作为 private variant，而事实上在千人基因组中也存在这种突变。

同样，CTTT->CTTTT 突变可能不会记录在数据库中，因为 C->CT 是记录这个突变更合适的方式。

1.关于如何描述 indel，目前尚无共识。许多人喜欢进行左归一化（left-normalization）。左归一化是指将突变的开始位置移至最小。但是，HGVS 明确指定将在 cDNA（mRNA）坐标上执行左归一化，这意味着人类基因组中一半基因需要进行右归一化。2.阅读这些事实之后，现在的问题是，我们应该如何对 VCF 文件进行更准确的注释？

由于左归一化越来越流行，因此我的建议是只使用左归一化，并且数据库和用户都使用这种做法以便我们可以将两种数据进行比较。我的第二个建议是，每条 VCF 记录仅描述一个单一的突变，这样 indel 就不会劫持 SNP，以确保与数据库进行1对1匹配。

因此，作为 ANNOVAR 开发人员，我决定重新处理所有千人基因组数据库以及 ESP6500si 数据库和 dbSNP 数据库，实现每一行仅包含一个突变，并且使每个突变都进行左归一化。更新的数据库已于 2014 年 12 月公布。

因此，作为用户，你应该：

•第一步：拆分 VCF，使之每行只有一个突变•第二步：对所有 VCF 行进行左归一化•第三步：用 ANNOVAR 进行注释。

例如，假设输入是 ex1.vcf.gz，你应该执行以下命令：

bcftools norm -m-both -o ex1.step1.vcf ex1.vcf.gz    bcftools norm -f human_g1k_v37.fasta -o ex1.step2.vcf ex1.step1.vcf

第一个命令表示把包含多个等位基因突变的行拆分成单独的行，第二个命令执行左归一化。（有时候进行第一步操作时会遇到 bug，即存在这样的突变但命令却无法拆分，遇到这样的情况可以尝试 vt[2]。在完成预处理步骤之后，你可以开始用 ANNOVAR 注释 ex1.step2.vcf 了。

1.你应牢记上述方法存在几个问题。

首先，正链基因很可能在不同参考基因组中被转化为负链基因（例如，同一参考基因组的不同版本，例如 hg17/hg18/hg19/hg38，或者来自不同种族群体的参考基因组，例如高加索人/非洲人/中国人/朝鲜人），因此左归一化会导致一些不一致的蛋白质注释；但如果我们采用 HGVS 标准，这就不是问题了。无论如何，由于当今人类基因组已经相对成熟，我认为这至少对于人类基因组来说这是一个相对较小的问题。

其次，由于 VCF 的设计方式，进行左归一化和拆分的软件工具并不像你想象的那样聪明，并且 INFO 字段可能无法正确拆分。我们以一个 ESP6500 文件中的记录为例： EA_AC = 76,129,1560 可能出现在 VCF 文件的 INFO 字段中，但它表示替代等位基因 1 的计数，替代等位基因 2 的计数，参考等位基因的计数（但是，诸如 bcftools 之类的软件并不知道这些隐藏含义，也不知道等位基因的确切顺序）。现在，如果你拆分这个 VCF 并对它进行左归一化，任何软件工具都不会足够聪明地重新生成正确的记录，导致你无法再正确解释 INFO 字段。为了解决这个问题，你可能会在同一位置重新加入多个突变并生成一个新的 VCF 文件。

当前，ANNOVAR 中的以下数据库是左归一化的，你可以直接使用它们与你的左归一化 VCF 文件进行比较：

•avsnp138•avsnp142•clinvar_20150330•1000g2014oct•exac03•esp6500siv2

引用链接

[1] http://www.1000genomes.org/wiki/Analysis/Variant%2520Call%2520Format/vcf-variant-call-format-version-41 [2] https://genome.sph.umich.edu/wiki/Vt

生信技能树目前已经公开了三个生信知识库，记得来关注哦~

每周文献分享

https://www.yuque.com/biotrainee/weeklypaper

肿瘤外显子分析指南

https://www.yuque.com/biotrainee/wes

生物统计从理论到实践

https://www.yuque.com/biotrainee/biostat