SnpSift学习笔记(三)-CSDN博客

欢迎关注"生信修炼手册"！

本篇主要介绍caseControl, rmRefGen, tstv, rmInfo, gt, vcfcheck这6个命令的用法。

1. caseControl

caseControl命令用来计算不同样本分组中的纯合突变和杂合突变的样本数，以及突变位点的总数，在表示样本分组时，可以采用命令行参数的形式，也可以采用TFAM格式的文件来指定。

用命令行参数表示样本分组时，+代表case组，-代表control组，0代表忽略该样本，用法如下

java -jar SnpSift.jar caseControl "++++0-----" cc.vcf

cc.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/0

case组的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.vcf

2. 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.142

4. 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=348

5. gt

gt命令用于压缩VCF文件，主要是压缩其中各个样本的基因型信息，以减少VCF文件的大小。这个文件压缩的方法效率是非常高的，对于1000多个样本的VCF文件，可以由原来的1TB压缩到1G。

其压缩的方式也很简单，类似稀疏矩阵的存储方式，由于大部分样本都是没有突变的，基因型都为0/0, 所以只记录其中发生了突变位点的样本数，在实际处理时，只记录下面3种类型的样本数

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参数，也可以解压缩，还原出之前的基因型信息，用法如下