vcf文件中短序列插入缺失（INDEL）标准化

问题：INDEL在vcf文件中表示的非唯一特性

INDEL 在Wikipedia中的解释：

Indel is a molecular biology term for an insertion or deletion of bases in the genome of an organism. It is classified among small genetic variations, measuring from 1 to 10 000 base pairs in length...

由于牵涉基因组上多个位点的变化，同时变化情况多样，导致同一种序列的变化可以写成多种样子，而不同的软件也会造成这种情况，Berkeley的SMaSH页面详细说明了这个问题，我这里简单说明：

序列变化 第一种写法 第二种写法 第三种写法 第四种写法

由以上可以看出，针对一个序列变化，至少有4种写法，而在vcf文件中的不同写法可以理解成不同的变化，对于后期的数据解读简直是灾难。

这个问题并不是新发现。

• 密歇根大学的统计遗传学中心关于变异的标准化的wiki页面这样说：

  ...However, variant representation in VCF is non-unique for variants that have explicitly expressed reference and alternate sequences. A failure to recognize this will frequently result in inaccurate analyses...

  他们认为变异标准vcf表示应该包含两个部分：

  1. 简单性原则，即保证ref和alt的长度不为0的前提下用尽量少的碱基表示该变异；

  2. 以变异的长度和基因组位置的自然表达为前提，将其向左对齐。
     .
     该页面还提供了实现向左对齐的算法示例，如下图。

     
     向左对齐算法

• 非常有名的变异注释工具ANNOVAR，其作者王凯在该软件的文档中对该问题也作了简单描述和自己的看法。

  ...There is no community consensus yet on how to describe an indel in an unique way. Many users prefer to do a left-normalization...

  他建议的做法，也是ANNOVAR的做法是，将我们生成的vcf文件统一做左对齐的标准化，同时对注释使用的数据库也做左对齐的标准化，这样避免了很多潜在的问题。事实上，ANNOVAR下载的数据库包括1000g、ESP6500si、dbSNP等都作了预标准化工作。

我们怎么办？

我想通过以上简单的介绍，大家能做一个简单的事情，那就是向王凯那样对所有的vcf文件做向左对齐的标准化工作，同时对公开和自建的数据库也做向左对齐的标准化，再使用工具进行注释。

向大家介绍一个工具，就是密歇根大学的vt。vt工具是密歇根大学开发的用来处理二代测序产生的短序列变异的工具。向左对齐的标准化命令是

# normalize variants and write out to dbsnp.normalized.vcf
vt normalize dbsnp.vcf -r seq.fa -o dbsnp.normalized.vcf

P.S.

1. SpeedSeq产生的vcf文件是将SNP和INDEL柔和在一起的，可以通过不同的工具将两者分开成两个单独的vcf文件，只对INDEL的vcf文件做处理就可以了。

2. 很多时候我们vcf文件中单行表示了多种变异情况，同样可以使用vt工具将多种变异情况的单行记录拆分成多行，命令行如下。

# decomposes multiallelic variants into biallelic variants and write out to gatk.decomposed.vcf
vt decompose gatk.vcf -o gatk.decomposed.vcf 
# before decomposition
#CHROM  POS     ID   REF     ALT         QUAL   FILTER  INFO                  FORMAT    S1                                     S2                                                                          
1       3759889 .    TA      TAA,TAAA,T  .      PASS    AF=0.342,0.173,0.037  GT:DP:PL      1/2:81:281,5,9,58,0,115,338,46,116,809     0/0:86:0,30,323,31,365,483,38,291,325,567  

# after decomposition
#CHROM  POS     ID   REF     ALT         QUAL   FILTER  INFO                                        FORMAT   S1               S2             
1   3759889 .    TA      TAA     .      PASS    OLD_MULTIALLELIC=1:3759889:TA/TAA/TAAA/T    GT:PL    1/.:281,5,9      0/0:0,30,323    
1   3759889 .    TA      TAAA        .      .       OLD_MULTIALLELIC=1:3759889:TA/TAA/TAAA/T    GT:PL    ./1:281,58,115   0/0:0,31,483    
1   3759889 .    TA      T           .      .       OLD_MULTIALLELIC=1:3759889:TA/TAA/TAAA/T    GT:PL    ./.:281,338,809  0/0:0,38,567

参考：

1. INDEL Wikipedia：https://en.wikipedia.org/wiki/Indel；
2. Berkeley SMaSH: http://smash.cs.berkeley.edu/normalize.html；
3. Umich Center of Statistical Genetics Wiki: http://genome.sph.umich.edu/wiki/Variant_Normalization；
4. ANNOVAR documentation: http://annovar.openbioinformatics.org/en/latest/articles/VCF/；
5. Umich vt: http://genome.sph.umich.edu/wiki/Vt#Normalization。