VCF | 多等位基因位点如何拆分？InDel变异如何标准化？

统一标准

多等位基因位点 (Multiallelic sites) 的拆分，左对齐标准化InDel，对ANNOVAR注释、MAF文件转换、映射人群频率、映射致病性等在少数特定位点上有影响，除非不分析多等位基因 (可能是因为分析起来过于复杂，或者这些位点不是那么重要，很多研究和应用领域不考虑这些位点)。

工具：bcftools norm

bcftools norm

About:

Left-alignandnormalizeindels;

check if REF allelesmatchthe reference;

split multiallelic sitesinto multiple rows;

recover multiallelics frommultiple rows.

Usage: bcftools norm [options] <in.vcf.gz>

Options:

-c, --check-ref <e|w|x|s>check REFalleles and exit (e), warn (w), exclude (x), or set (s) bad sites [e]

-D, --remove-duplicatesremoveduplicate linesof the same type.

-d, --rm-dup <type> removeduplicate snps|indels|both|all|none

-f, --fasta-ref <file>reference sequence (MANDATORY)

-m, --multiallelics <-|+>[type]splitmultiallelics (-) orjoinbiallelics (+),type: snps|indels|both|any [both]

--no-version do not append version and command line to the header

-N, --do-not-normalizedo not normalize indels(with -m or -c s)

-o, --output <file> write output to a file [standard output]

-O, --output-type <type>'b' compressed BCF; 'u' uncompressed BCF; 'z' compressed VCF; 'v' uncompressed VCF [v]

-r, --regions <region>restrict to comma-separated list of regions

-R, --regions-file <file> restrict to regions listed in a file

-s, --strict-filter when merging (-m+), merged site is PASS only if all sites being merged PASS

-t, --targets <region>similar to -r but streams rather than index-jumps

-T, --targets-file <file> similar to -R but streams rather than index-jumps

--threads <int> number of extra (de)compression threads [0]

-w, --site-win <int> buffer for sorting lines which changed position during realignment [1000]

位点测试

位点1 - 含Del

#拆分前zcat ~/wgs/result/Genotype.cohort.dbSNP.g.vcf.gz | \grep"CHROM\|rs34364959" | sed 's/\t/ /g'# 拆分后bcftools norm -f /db/gatk/hg38/Homo_sapiens_assembly38.fasta \--multiallelics -both -Ov ~/wgs/result/Genotype.cohort.dbSNP.g.vcf.gz | \grep"CHROM\|rs34364959" | sed 's/\t/ /g'

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT CHB1 CHB2 GWD1 GWD2

chr112411436rs34364959CG TG,C652.11PASSAC=1,1;AF=0.125,0.125;AN=8;BaseQRankSum=2.92;DB;DP=227;ExcessHet=3.6798;FS=3.302;MLEAC=1,1;MLEAF=0.125,0.125;MQ=59.87;MQRankSum=0;QD=7.25;ReadPosRankSum=-1.212;SOR=1.121GT:AD:DP:GQ:PGT:PID:PL:PS0/1:28,22,0:50:99:.:.:473,0,498,554,564,1118:.0|2:33,0,7:40:99:0|1:2411398_C_A:195,294,1680,0,1386,1365:2411398 0/0:63,0,0:63:0:.:.:0,0,1220,0,1220,1220:. 0/0:71,0,0:71:59:.:.:0,59,1675,59,1675,1675:.

拆分后

Lines total/split/realigned/skipped:1439/51/9/0

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT CHB1 CHB2 GWD1 GWD2

chr112411436rs34364959C T652.11PASSAC=1;AF=0.125;AN=8;BaseQRankSum=2.92;DB;DP=227;ExcessHet=3.6798;FS=3.302;MLEAC=1;MLEAF=0.125;MQ=59.87;MQRankSum=0;QD=7.25;ReadPosRankSum=-1.212;SOR=1.121GT:AD:DP:GQ:PGT:PID:PL:PS0/1:28,22:50:99:.:.:473,0,498:.0|0:33,0:40:99:0|1:2411398_C_A:195,294,1680:2411398 0/0:63,0:63:0:.:.:0,0,1220:. 0/0:71,0:71:59:.:.:0,59,1675:.

chr112411436rs34364959CG C652.11PASSAC=1;AF=0.125;AN=8;BaseQRankSum=2.92;DB;DP=227;ExcessHet=3.6798;FS=3.302;MLEAC=1;MLEAF=0.125;MQ=59.87;MQRankSum=0;QD=7.25;ReadPosRankSum=-1.212;SOR=1.121GT:AD:DP:GQ:PGT:PID:PL:PS0/0:28,0:50:99:.:.:473,554,1118:.0|1:33,7:40:99:0|1:2411398_C_A:195,0,1365:2411398 0/0:63,0:63:0:.:.:0,0,1220:. 0/0:71,0:71:59:.:.:0,59,1675:.

还不错，基本该做到的、都做到了；难以做到的、也都没做到

位点2 - 纯SNP

# 拆分前zcat ~/wgs/result/Genotype.cohort.dbSNP.g.vcf.gz | \grep "CHROM\|rs752579667" | sed 's/\t/ /g'# 拆分后bcftools norm -f /db/gatk/hg38/Homo_sapiens_assembly38.fasta \--multiallelics -both -Ov ~/wgs/result/Genotype.cohort.dbSNP.g.vcf.gz | \grep "CHROM\|rs752579667" | sed 's/\t/ /g'

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT CHB1 CHB2 GWD1 GWD2

chr113361836rs752579667C T,A69.04PASS AC=1,1;AF=0.125,0.125;AN=8;BaseQRankSum=2.24;DB;DP=86;ExcessHet=3.6798;FS=12.725;MLEAC=1,1;MLEAF=0.125,0.125;MQ=52.03;MQRankSum=-1.799;QD=2.16;ReadPosRankSum=1.59;SOR=1.057 GT:AD:DP:GQ:PGT:PID:PL:PS 0/0:17,0,0:17:23:.:.:0,23,365,23,365,365:. 0/0:35,0,0:35:78:.:.:0,78,1170,78,1170,1170:.0/1:13,2,1:17:42:.:.:45,0,540,42,504,585:.0|2:14,0,2:16:42:0|1:3361828_A_G:42,84,657,0,573,567:3361828

拆分后

Lines total/split/realigned/skipped:1439/51/9/0

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT CHB1 CHB2 GWD1 GWD2

chr113361836rs752579667C T69.04 PASS AC=1;AF=0.125;AN=8;BaseQRankSum=2.24;DB;DP=86;ExcessHet=3.6798;FS=12.725;MLEAC=1;MLEAF=0.125;MQ=52.03;MQRankSum=-1.799;QD=2.16;ReadPosRankSum=1.59;SOR=1.057 GT:AD:DP:GQ:PGT:PID:PL:PS 0/0:17,0:17:23:.:.:0,23,365:. 0/0:35,0:35:78:.:.:0,78,1170:.0/1:13,2:17:42:.:.:45,0,540:. 0|0:14,0:16:42:0|1:3361828_A_G:42,84,657:3361828

chr113361836rs752579667C A69.04 PASS AC=1;AF=0.125;AN=8;BaseQRankSum=2.24;DB;DP=86;ExcessHet=3.6798;FS=12.725;MLEAC=1;MLEAF=0.125;MQ=52.03;MQRankSum=-1.799;QD=2.16;ReadPosRankSum=1.59;SOR=1.057 GT:AD:DP:GQ:PGT:PID:PL:PS 0/0:17,0:17:23:.:.:0,23,365:. 0/0:35,0:35:78:.:.:0,78,1170:. 0/0:13,1:17:42:.:.:45,42,585:.0|1:14,2:16:42:0|1:3361828_A_G:42,0,567:3361828

位点3 - 含Del

#拆分前zcat ~/wgs/result/Genotype.cohort.dbSNP.g.vcf.gz | \grep "CHROM\|rs1050228\|rs749512034\|rs558809956" | sed 's/\t/ /g'#拆分后bcftools norm -f /db/gatk/hg38/Homo_sapiens_assembly38.fasta \--multiallelics -both -Ov ~/wgs/result/Genotype.cohort.dbSNP.g.vcf.gz | \grep "CHROM\|rs1050228\|rs749512034\|rs558809956" | sed 's/\t/ /g'

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT CHB1 CHB2 GWD1 GWD2

chr116390705rs1050228;rs749512034;rs558809956TGCTGGCGCTGGCGCTGGCTGCTGGCGCTGGC,CGCTGGCGCTGGCGCTGGC,T2139.3 PASS AC=4,2,1;AF=0.5,0.25,0.125;AN=8;BaseQRankSum=0.105;DB;DP=126;ExcessHet=3.0103;FS=1.475;MLEAC=4,2,1;MLEAF=0.5,0.25,0.125;MQ=59.95;MQRankSum=0;QD=30.56;ReadPosRankSum=-1.026;SOR=0.399 GT:AD:DP:GQ:PL 1/2:0,6,7,0:13:99:372,166,199,220,0,254,387,198,250,433 1/2:0,4,3,0:7:86:245,86,114,168,0,159,252,109,168,270 0/1:10,14,0,0:24:99:546,0,331,576,375,951,576,375,951,951 1/3:0,13,0,13:26:99:1007,516,476,1048,548,1237,503,0,729,878

拆分后

Lines total/split/realigned/skipped: 1439/51/9/0

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT CHB1 CHB2 GWD1 GWD2

chr116390705rs1050228;rs749512034;rs558809956TGCTGGCT2139.3 PASS AC=4;AF=0.5;AN=8;BaseQRankSum=0.105;DB;DP=126;ExcessHet=3.0103;FS=1.475;MLEAC=4;MLEAF=0.5;MQ=59.95;MQRankSum=0;QD=30.56;ReadPosRankSum=-1.026;SOR=0.399 GT:AD:DP:GQ:PL 1/0:0,6:13:99:372,166,199 1/0:0,4:7:86:245,86,114 0/1:10,14:24:99:546,0,331 1/0:0,13:26:99:1007,516,476

chr116390705rs1050228;rs749512034;rs558809956TC2139.3 PASS AC=2;AF=0.25;AN=8;BaseQRankSum=0.105;DB;DP=126;ExcessHet=3.0103;FS=1.475;MLEAC=2;MLEAF=0.25;MQ=59.95;MQRankSum=0;QD=30.56;ReadPosRankSum=-1.026;SOR=0.399 GT:AD:DP:GQ:PL 0/1:0,7:13:99:372,220,254 0/1:0,3:7:86:245,168,159 0/0:10,0:24:99:546,576,951 0/0:0,0:26:99:1007,1048,1237

chr116390705rs1050228;rs749512034;rs558809956TGCTGGCGCTGGCGCTGGC T2139.3 PASS AC=1;AF=0.125;AN=8;BaseQRankSum=0.105;DB;DP=126;ExcessHet=3.0103;FS=1.475;MLEAC=1;MLEAF=0.125;MQ=59.95;MQRankSum=0;QD=30.56;ReadPosRankSum=-1.026;SOR=0.399 GT:AD:DP:GQ:PL 0/0:0,0:13:99:372,387,433 0/0:0,0:7:86:245,252,270 0/0:10,0:24:99:546,576,951 0/1:0,13:26:99:1007,503,878

位点4- 含Ins

# chr11:281137# 拆分前zcat ~/wgs/result/Genotype.cohort.dbSNP.g.vcf.gz | \grep -P 'CHROM|chr11\t281137' | sed 's/\t/ /g'# 拆分后bcftools norm -f /db/gatk/hg38/Homo_sapiens_assembly38.fasta \--multiallelics -both -Ov ~/wgs/result/Genotype.cohort.dbSNP.g.vcf.gz | \grep -P 'CHROM|chr11\t281137' | sed 's/\t/ /g'

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT CHB1 CHB2 GWD1 GWD2

chr11281137.TTTTTCTCGCAATACGGCCGTATCATCACCTCACGAATCCTGGTTGATCAAGTCACA,TTTTCTCGCAATAC1406.78 PASS AC=2,1;AF=0.25,0.125;AN=8;BaseQRankSum=-0.471;DP=417;ExcessHet=5.4407;FS=157.296;MLEAC=2,1;MLEAF=0.25,0.125;MQ=57.36;MQRankSum=-11.61;QD=4.08;ReadPosRankSum=-0.004;SOR=8.051 GT:AD:DP:GQ:PGT:PID:PL:PS0|2:25,0,26:51:99:0|1:281137_T_TTTTCTCGCAATAC:1017,1092,2108,0,1016,938:281137 0/0:56,0,0:56:90:.:.:0,90,1350,90,1350,1350:.0/1:143,15,0:158:99:.:.:195,0,5834,630,5879,6509:.0/1:122,14,0:136:99:.:.:215,0,5041,588,5083,5671:.

拆分后

Lines total/split/realigned/skipped:1439/51/9/0

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT CHB1 CHB2 GWD1 GWD2

chr11281137.TTTTTCTCGCAATACGGCCGTATCATCACCTCACGAATCCTGGTTGATCAAGTCACA 1406.78 PASS AC=2;AF=0.25;AN=8;BaseQRankSum=-0.471;DP=417;ExcessHet=5.4407;FS=157.296;MLEAC=2;MLEAF=0.25;MQ=57.36;MQRankSum=-11.61;QD=4.08;ReadPosRankSum=-0.004;SOR=8.051 GT:AD:DP:GQ:PGT:PID:PL:PS 0|0:25,0:51:99:0|1:281137_T_TTTTCTCGCAATAC:1017,1092,2108:281137 0/0:56,0:56:90:.:.:0,90,1350:.0/1:143,15:158:99:.:.:195,0,5834:.0/1:122,14:136:99:.:.:215,0,5041:.

chr11281137.TTTTTCTCGCAATAC1406.78 PASS AC=1;AF=0.125;AN=8;BaseQRankSum=-0.471;DP=417;ExcessHet=5.4407;FS=157.296;MLEAC=1;MLEAF=0.125;MQ=57.36;MQRankSum=-11.61;QD=4.08;ReadPosRankSum=-0.004;SOR=8.051 GT:AD:DP:GQ:PGT:PID:PL:PS0|1:25,26:51:99:0|1:281137_T_TTTTCTCGCAATAC:1017,0,938:281137 0/0:56,0:56:90:.:.:0,90,1350:. 0/0:143,0:158:99:.:.:195,630,6509:. 0/0:122,0:136:99:.:.:215,588,5671:.

使用方法

#变异标准化#多等位基因位点(Multiallelicsites)的拆分bcftools norm -f ${ref} --multiallelics -both ${result}/Genotype.cohort.dbSNP.g.vcf \-Oz-o${result}/Genotype.cohort.dbSNP.g.vcf.gz

变异标准化及左对齐等相关理论

下面的例子中，多核苷酸多态性（Multi Nucleotide Polymorphism,MNP）的前 3个核苷酸是多余的 (Representedsuperfluously)，而从第 4 个核苷酸开始表示时，是简约化的 (Parsimoniously)。

Variant normalization

当一个变异在左侧有多余的核苷酸时，它被定义为"非左简约化"，因为需要向左修剪 (Defined asnot beingleft parsimoniousas there is a need toleft trim)。相对于这个概念有：Right parsimony and trimming。

Parsimony 也适用于 Indel，及左对齐 (Left alignment)。

更多解释：

https://genome.sph.umich.edu/wiki/Variant_Normalization

Adrian Tan, Gonçalo R. Abecasis and Hyun Min Kang. ()Unified Representation of Genetic Variants. Bioinformatics.

如果觉得《VCF | 多等位基因位点如何拆分？InDel变异如何标准化？》对你有帮助，请点赞、收藏，并留下你的观点哦！