SnpSift学习笔记(三)
本篇主要介绍caseControl, rmRefGen, tstv, rmInfo, gt, vcfcheck这6个命令的用法。
1. caseControl
caseControl命令用来计算不同样本分组中的纯合突变和杂合突变的样本数，以及突变位点的总数，在表示样本分组时，可以采用命令行参数的形式，也可以采用TFAM格式的文件来指定。
用命令行参数表示样本分组时，+代表case组，-代表control组，0代表忽略该样本，用法如下
java -jar SnpSift.jar caseControl "++++0-----" cc.vcfcc.vcf中包含了10个样本的基因型信息，其中1-4号样本为case组，5号样本忽略不计，6-10号样本为control组。输出结果会在INFO字段中添加如下信息
INFO FORMAT Sample_01 Sample_02 Sample_03 Sample_04 Sample_05 Sample_06 Sample_07 Sample_08 Sample_09 Sample_10
AF=0.01;Cases=1,2,4;Controls=2,2,6 GT 0/1 1/1 1/0 0/0 0/0 0/1 1/1 1/1 1/0 0/0case组的4个样本对应的基因型分别为0/1,1/1,1/0,0/0,其中纯合突变的样本数为1个，杂合突变的样本数为2个，突变位点总数为1个纯合的2个突变位点加上2个杂合的2个突变位点，一共是4个突变位点；control组5个样本对应的基因型分别为0/1, 1/1, 1/1, 1/0,0/0,其中纯合突变的样本数为2个，杂合突变的样本数为2个，突变位点总数为2 * 2 + 2 = 6个。
通过命令行参数表示样本分组时，必须按照VCF文件中样本的顺序依次进行指定。如果采用TFAM文件指定样本分组，则不需要考虑样本顺序，用法如下
java -jar SnpSift.jar caseControl -tfam samples.tfam cc.vcf2. rmRefGen
rmRefGen命令用于删除VCF文件中没有发生突变的样本，用法如下
cat file.vcf | java -jar SnpSift.jar rmRefGen > file_noref.vcf原始的file.vcf文件内容如下
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT M1 M2 X1 X2
2L 426906 . C G 53.30 . DP=169 GT:PL:GQ 0/1:7,0,255:4 0/1:7,0,255:4 0/0:0,0,0:6 0/0:0,30,255:35
2L 601171 . C A 999.00 . DP=154 GT:PL:GQ 0/1:81,0,141:78 0/1:42,0,251:39 0/0:0,0,0:4 0/0:0,33,255:36
2L 648611 . A T 999.00 . DP=225 GT:PL:GQ 0/1:52,0,42:47 1/1:75,21,0:14 0/0:0,0,0:3 0/0:0,60,255:61
2L 807373 . A G 106.00 . DP=349 GT:PL:GQ 0/1:14,0,65:12 0/1:60,0,42:50 0/0:0,0,0:4 0/0:0,69,255:72
2L 816766 . G T 999.00 . DP=411 GT:PL:GQ 0/1:108,0,45:53 0/1:7,0,255:6 0/0:0,0,0:4 0/0:0,57,255:59可以看到X1和X2两个样本的基因型都为0/0,没有发生突变。rmRefGen处理之后，输出结果如下
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT M1 M2 X1 X2
2L 426906 . C G 53.30 . DP=169 GT:PL:GQ 0/1:7,0,255:4 0/1:7,0,255:4 . .
2L 601171 . C A 999.00 . DP=154 GT:PL:GQ 0/1:81,0,141:78 0/1:42,0,251:39 . .
2L 648611 . A T 999.00 . DP=225 GT:PL:GQ 0/1:52,0,42:47 1/1:75,21,0:14 . .
2L 807373 . A G 106.00 . DP=349 GT:PL:GQ 0/1:14,0,65:12 0/1:60,0,42:50 . .
2L 816766 . G T 999.00 . DP=411 GT:PL:GQ 0/1:108,0,45:53 0/1:7,0,255:6 . .可以看到，X1和X2样本对应的记录都变成了.。
3. tstv
tstv命令用于计算突变位点中，转换和颠换的比例，用法如下
java -jar SnpSift.jar tstv hom s.vcf输出结果如下
Sample : 1 2 3 4 5 6 7 8 9 10 11 12 Total
Transitions : 150488 150464 158752 156674 152936 160356 152276 155314 156484 149276 151182 153468 1847670
Transversions : 70878 70358 73688 72434 70828 76150 72030 71958 72960 69348 70180 71688 862500
Ts/Tv : 2.123 2.139 2.154 2.163 2.159 2.106 2.114 2.158 2.145 2.153 2.154 2.141 2.1424. rmInfo
rmInfo命令用于删除VCF文件中INFO中的sub fields的信息，用法如下
java -jar SnpSift.jar rmInfo test.snpeff.vcf EFF这里的EFF代表你想要删除的字段的名字。
原始VCF文件内容如下
cat test.snpeff.vcf
#CHROM POS ID REF ALT QUAL FILTER INFO
1 734462 1032 G A . s50 AC=348;EFF=DOWNSTREAM(MODIFIER|||||RP11-206L10.8|processed_transcript|NON_CODING|ENST00000447500||1),INTRON(MODIFIER|||||RP11-206L10.6|processed_transcript|NON_CODING|ENST00000429505|1|1)删除之后的内容如下
#CHROM POS ID REF ALT QUAL FILTER INFO
1 734462 1032 G A . s50 AC=3485. gt
gt命令用于压缩VCF文件，主要是压缩其中各个样本的基因型信息，以减少VCF文件的大小。这个文件压缩的方法效率是非常高的，对于1000多个样本的VCF文件，可以由原来的1TB压缩到1G。
其压缩的方式也很简单，类似稀疏矩阵的存储方式，由于大部分样本都是没有突变的，基因型都为0/0, 所以只记录其中发生了突变位点的样本数，在实际处理时，只记录下面3种类型的样本数
- HO
- HE
- NA
HO代表纯合突变，基因型为1/1;HE代表杂合突变，基因型为0/1; NA代表没有基因型信息。
用法如下
java -jar SnpSift.jar gt test.vcf | tee test.gt.vcf输出结果如下：
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample_1 Sample_2 Sample_3 Sample_4 Sample_5 Sample_6 Sample_7 Sample_8 Sample_9 Sample_10 Sample_11 Sample_12 Sample_13 Sample_14 Sample_15
1 861276 . A G . PASS AC=1;HO=1所有样本的基因型信息都没有了，取而代之的是INFO字段中，记录的HO纯合突变的样本数。
除了压缩之外，通过添加-u参数，也可以解压缩，还原出之前的基因型信息，用法如下
java -jar SnpSift.jar gt -u test.gt.vcf6. vcfCheck
vcfCheck命令用于检测VCF文件的格式是否有问题，用法如下
java -jar SnpSift.jar vcfCheck bad.vcf当VCF存在一些常见错误时，会打印报错信息，可以当做一个debug的工具来使用。