SnpSift学习笔记(三)

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

本篇主要介绍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

可以看到X1X2两个样本的基因型都为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 . .

可以看到,X1X2样本对应的记录都变成了.

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种类型的样本数

  1. HO

  2. HE

  3. 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.vcf

6. vcfCheck

vcfCheck命令用于检测VCF文件的格式是否有问题,用法如下

java -jar SnpSift.jar vcfCheck bad.vcf

当VCF存在一些常见错误时,会打印报错信息,可以当做一个debug的工具来使用。

扫描关注微信号,更多精彩内容等着你!


©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 216,372评论 6 498
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,368评论 3 392
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 162,415评论 0 353
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,157评论 1 292
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,171评论 6 388
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,125评论 1 297
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,028评论 3 417
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,887评论 0 274
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,310评论 1 310
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,533评论 2 332
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,690评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,411评论 5 343
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,004评论 3 325
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,659评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,812评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,693评论 2 368
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,577评论 2 353

推荐阅读更多精彩内容