snpSift的吊诡之处

1. 用作测试的数据库vcf文件： database_test.vcf

##fileformat=VCFv4.1
##fileDate=2021-11-06
##source=Inhouse
##reference=GRCh37
##ID=<Description="Variation ID">
##INFO=<ID=TEST,Number=.,Type=String,Description="the code representing an annotation">
#CHROM    POS    ID    REF    ALT    QUAL    FILTER    INFO
chr17    41243479    .    CTTGA    C    .    .    TEST=anno1
chr17    41243481    .    TGATT    T    .    .    TEST=anno2

然后将其压缩并建立index

bgzip -c database_test.vcf >database_test.vcf.gz
tabix -p vcf database_test.vcf.gz

 用作测试的待注释vcf文件： test.vcf

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  21S03396220_HuangSuWen_DX0622
chr17   41243479        .       CTTGA   C       .       germline        CONTQ=93;DP=2289;ECNT=1;GERMQ=1;MBQ=20,20;MFRL=56,126;MMQ=60,60;MPOS=22;POPAF=0.225;ROQ=93;SEQQ=93;STRANDQ=93;TLOD=6581.37     GT:AD:AF:DP:F1R2:F2R1:SB        0/1:8,2225:0.996:2233:2,1385:6,829:3,5,1142,1083
chr17   41243481        .       TGATT   T       .       germline;haplotype      CONTQ=93;DP=1789;ECNT=2;GERMQ=1;MBQ=9,20;MFRL=137,128;MMQ=60,60;MPOS=23;POPAF=0.035;ROQ=93;SEQQ=93;STRANDQ=93;TLOD=6252.32     GT:AD:AF:DP:F1R2:F2R1:PGT:PID:PS:SB     0|1:2,1731:0.999:1733:0,995:0,709:0|1:11854457_G_A:11854457:0,2,867,864
 

开始注释：

java -Xmx4g -jar ./snpEff/snpEff/SnpSift.jar annotate ./05.tests/snpSift/database_test.vcf.gz -noId -name testSift_ -exists InTestSift ./05.tests/snpSift/test.vcf >./05.tests/snpSift/test.snpSift.vcf

 结果文件（仅展示有效数据行）：

#CHROM    POS    ID    REF    ALT    QUAL    FILTER    INFO    FORMAT    21S03396220_HuangSuWen_DX0622
chr17    41243479    .    CTTGA    C    .    germline    CONTQ=93;DP=2289;ECNT=1;GERMQ=1;MBQ=20,20;MFRL=56,126;MMQ=60,60;MPOS=22;POPAF=0.225;ROQ=93;SEQQ=93;STRANDQ=93;TLOD=6581.37;testSift_TEST=anno1;InTestSift    GT:AD:AF:DP:F1R2:F2R1:SB    0/1:8,2225:0.996:2233:2,1385:6,829:3,5,1142,1083
chr17    41243481    .    TGATT    T    .    germline;haplotype    CONTQ=93;DP=1789;ECNT=2;GERMQ=1;MBQ=9,20;MFRL=137,128;MMQ=60,60;MPOS=23;POPAF=0.035;ROQ=93;SEQQ=93;STRANDQ=93;TLOD=6252.32    GT:AD:AF:DP:F1R2:F2R1:PGT:PID:PS:SB    0|1:2,1731:0.999:1733:0,995:0,709:0|1:11854457_G_A:11854457:0,2,867,864

可以看出，第二个突变：“chr17   41243481        .       TGATT   T ” 没有注释出来。为什么呢？

最后发现，原来是数据库文件中对应这个突变的那一行（已经是最后一行）没有换行符，也就是第二个突变这一行没有换行符。将数据库文件的最后一行加上换行符以后，再次压缩，建立index，再注释，就注释成功了。

#CHROM    POS    ID    REF    ALT    QUAL    FILTER    INFO    FORMAT    21S03396220_HuangSuWen_DX0622
chr17    41243479    .    CTTGA    C    .    germline    CONTQ=93;DP=2289;ECNT=1;GERMQ=1;MBQ=20,20;MFRL=56,126;MMQ=60,60;MPOS=22;POPAF=0.225;ROQ=93;SEQQ=93;STRANDQ=93;TLOD=6581.37;testSift_TEST=anno1;InTestSift    GT:AD:AF:DP:F1R2:F2R1:SB    0/1:8,2225:0.996:2233:2,1385:6,829:3,5,1142,1083
chr17    41243481    .    TGATT    T    .    germline;haplotype    CONTQ=93;DP=1789;ECNT=2;GERMQ=1;MBQ=9,20;MFRL=137,128;MMQ=60,60;MPOS=23;POPAF=0.035;ROQ=93;SEQQ=93;STRANDQ=93;TLOD=6252.32;testSift_TEST=anno2;InTestSift    GT:AD:AF:DP:F1R2:F2R1:PGT:PID:PS:SB    0|1:2,1731:0.999:1733:0,995:0,709:0|1:11854457_G_A:11854457:0,2,867,864
 

原则上，要求数据库文件和输入文件都要排序才行。

对输入文件是否排序，要求不是很严格。虽然会提出warning，仍然可以注释成功。

但对数据库文件要求严格，如果不排序，index都无法建立