变异位点的归一化(normalization of indel)

一、定义
二、indel的左对齐
三、多等位位点的左对齐

定义

在二代测序中，使用bcftools call或者 gatk HaplotypeCaller工具检出变异位点的时候，会得到原始的vcf文件，这些vcf文件中的indel一般是没有进行归一化的，这里的归一化主要针对两个方面:(1) indel的左对齐和(2) 多等位位点的左对齐。
如果不进行归一化，相同的变异可能会因为形式的不同而导致无法比较，例如下面的两个位点其实是相同的位点，如果不进行处理，对后续的注释和分析会带来很大的麻烦。

#CHROM  POS     ID      REF     ALT
1       900010  .       GC      GCC
1       900010  .       G      GC

indel的左对齐

如下面所示的vcf文件中检出的indel(这里仅展示前5列内容)

#CHROM  POS     ID      REF     ALT
1       900010  .       GC      GCC
1       1225506 .       TCCCCC  TCCC

通过左对齐后可以变成

#CHROM  POS     ID      REF     ALT
1       900010  .       G      GC
1       1225506 .       TCC    T

第一个indel可以进行左对齐简化，由G变成GC，表明在该位置发生了一个C的插入，而第一个indel可以去掉右边冗余的C，由 TCC变为T,发生了一个CC的缺失。
该过程我们可以使用bcftools norm来实现,关于该命令的详细用法请见bcftools文档,这里重点介绍几个参数:

bcftools norm的使用

基础使用方法
Usage: bcftools norm [options] file.vcf.gz
主要参数
1 -f/--fasta-ref 该参数接一个参考序列的fasta文件，例如人类的hg19，hg38的参考基因组。当出现该参数时，bcftools norm命令将会对vcf文件中的indel进行左对齐。
2 -m/multiallelics 该参数的值由两部分组成，第一部分的可选值为[+/-]分别对应于对多等位位点进行合并(+)或者拆分(-)。第二部分的可选值为[snps|indels|both|any]该部分指定了针对哪些类型的变异执行多等位位点的合并或者拆分，针对snps，或者indels或者两者都进行合并或者拆分both，或者对snp和indel合并为同一个记录any。

接下来使用bcftools norm来实现上述过程：

bcftools norm -f b37.fasta -m +both test.vcf -Ov -o test.normed.vcf

命令中的-Ov为 bcftools的通用命令,我们可以指定输出的文件类型，b为 compressed bcf(压缩的bcf文件)，u为uncompressed bcf(不压缩的bcf文件),v为 uncompressed vcf(未压缩的vcf文件),z为compressed vcf(压缩的vcf文件)。

多等位位点的左对齐

由于人类是二倍体生物，同源染色体都是以成对的方式出现的，因此会出现多等位变异。请观察如下的几个位点：

#CHROM  POS     ID      REF     ALT
1       10002   .       A       T,G
1       2332681 .       TCC     TCCCC,TCCC

上述的几个位点为多等位位点，两条同源染色体在相同的位置分别发生了不同的变异。通过左对齐可以变成如下的样子：

#CHROM  POS     ID      REF     ALT
1       10002   .       A       T,G
1       2332681 .       T     TCC,TC

snp的多等位位点不发生改变，而indel的多等位位点的REF和ALT可以同时去掉冗余的两个CC碱基。
上述的过程可以使用bcftools norm命令来实现

bcftools norm -f b37.fasta -m +both test.vcf -Ov -o test.normed.vcf

可以看到和普通的indel的归一化并没有什么区别，这里你可能有疑问既然两者的归一化的实现方式相同，为什么要将indel的归一化和多等位位点的归一化分开来分析呢？
让我们来看下面的例子：

#CHROM  POS     ID      REF     ALT
1       10003   .       AC      AG,ACT

上述的位点是一个多等位位点，两条同源染色体在10003这个位置上分别发生了AC到AG和AC到ACT的变异，这里可以注意到AC到AG可以去掉碱基A,而AC到ACT可以通过左对齐来变成C到CT,接下来让我们将该位点放到vcf文件中并使用命令bcftools norm -f b37.fasta -m +both test.vcf -Ov -o test.normed.vcf来进行归一化，结果如下：

#CHROM  POS     ID      REF     ALT
1       10004   .       C       G,CT

你会发现该位置的多等位位点变成了一个snp和一个indel！这显然不是我们想要的，为了解决这个问题，我们可以先对多等位位点进行拆分，在左对齐后，再进行合并，合并的时候，我们可以选择-m +both而不让snp和indel合并到一起（如果你使用-m +any参数，结果它们将再次合并到一起）。

bcftools norm -m -both test.vcf |
bcftools norm -m +both -f b37.fasta test.vcf -Ov -o test.normed.vcf

最后得到如下的结果

#CHROM  POS     ID      REF     ALT
1       10004   .       C       G
1       10004   .       C       CT

巧妙的借助管道符来解决这个问题，这条命令请及时做下笔记。