【直播】我的基因组54:把我的variation跟dbSNP数据库相比较

问:有许多朋友后台留言问我,为什么找变异的步骤不用GATK的best practice呢?而是选取bcftools,freebayes这种小众软件。

答:我这里统一回复一下,不同软件找到的variation的差别我前面已经说了,它们小众,并不代表不堪大用,对一个真正的variation来说,不管是什么软件,都是可以找到的,对一个模棱两可的variation,我一定会去去IGV用肉眼查看的,毕竟这可事关我的健康呀,不能马虎的!不同软件就是在sensitivity和specificity之间找平衡,而我早期并不需要有多么精的sensitivity,那些模棱两可的位点,干脆就不要报告给我,本来位点几百万就够我头疼的了,先把这些有把握的位点给探索一遍吧,等将来有空了我再回过头来看看是不是我的基因组还一些待挖掘的细节。

前面考虑到X,Y染色体的variation不是很可靠,本次分析我就统统排除掉了,直接从常染色体的高质量(bcftools quality>20)的variation开始分析起。

我们很简单统计一下杂合变异位点跟纯和变异位点的分布,用的是perl和shell脚本。

对INDEL的统计结果如下:

grep INDEL autochr.highQuali.dbsnp.vcf |perl -alne '{@tmp=split/:/,$F[9];print $tmp[0]}' |sort |uniq -c

结果如下:

288977 0/1

272779 1/1

39671 1/2

杂合纯合比例差不多,说明本次call variation的步骤还算合理。

对SNV的统计结果如下:

grep -v INDEL autochr.highQuali.dbsnp.vcf |perl -alne '{@tmp=split/:/,$F[9];print $tmp[0]}' |sort |uniq -c

结果如下:

2260576 0/1

1540114 1/1

1739 1/2

同时也统计了千人基因组计划(20130502版本)的2504个人的杂合纯合比例情况!证明了我的杂合纯合比例是符合大众规律的。

也下载了dbSNP(b147_GRCh37p13版本),并且把我的VCF文件注释到了dbSNP,就可以进行基本的统计啦!

有了这些信息,就可以进行下面的统计了!

代码如下:

cat autochr.highQuali.dbsnp.vcf |perl -alne 'next if /^#/;$type="SNV",$KGPhase3="NO";$type="INDEL" if /INDEL/;$KGPhase3="KGPhase3" if /KGPhase3/;$h=(split/:/,$F[9])[0];print substr($F[2],0,2),"\t$type\t$KGPhase3\t$h"'|sort |uniq -c

统计结果如下:

带rs标记的说明这个位点在dbSNP里面有记录,带有KGPhase3的说明在千人基因组计划里面有记录!在千人基因组计划里面发现了的snp一定在dbSNP里面有记录!

这个表格看起来不够直观,需要可视化如下:

每次画图我都很头疼,我简单说明一下上面的图吧!

3种颜色,NO代表着dbSNP(b147_GRCh37p13版本)和千人基因组计划(20130502版本)都没有记载,是我本人的全新突变!!而NOrs代表着在dbSNP有,在千人里面没有。而KGPhase3rs代表着在dbSNP和千人都有啦!

左边是INDEL,右边是SNV咯,可以看出来INDEL的全新太多了,有两个可能

第一,每个人之间的INDEL差异真的很大

第二,我这个找INDEL真的很不准确!

至于0/1,1/1,1/2就不用我多说了,分别是杂合,纯合的突变,其中1/2这种情况不是很好理解,暂时不需要深究!

代码很简单,就是把上面的数据导入R里面,用ggplot即可:

  1. a=read.table('type.txt',stringsAsFactors = F) ##这个type.txt文件就是上面截图的数据
  2. library(ggplot2)
  3. ggplot(a,aes(x=V5,y=V1))+geom_bar(stat = 'identity',aes(fill=paste0(V4,V2)))+facet_wrap(~V3)

其实正常的科研用图应该是下面这个,(⊙﹏⊙)b,高下立判!

图来源于:https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3737162/figure/fig1/

文:Jimmy

图文编辑:吃瓜群众

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-01-28

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏java一日一条

健康知识图谱,阿里工程师如何实现?

阿里妹导读:“The world is not made of strings , but is made of things.” 大千世界,万物相联。保险领...

812
来自专栏生信宝典

ChIP-seq基本分析流程

前年在中科院做培训时,整理了一套ChIP-seq分析流程,截选实战部分,略作修改,分享出来,希望大家指正。 那次培训内容比较杂,还有关于TCGA、ICGC、Pr...

35210
来自专栏ThoughtWorks

其实你早就学过响应式编程 | TW洞见

今日洞见 文章作者/图片 来自ThoughtWorks:佟达。 本文所有内容,包括文字、图片和音视频资料,版权均属ThoughtWorks公司所有,任何媒体、网...

33312
来自专栏FreeBuf

网络蜜罐技术探讨

*原创作者bt0sea,本文属FreeBuf原创奖励计划,未经许可禁止转载 之前在FreeBuf发表的第一篇有关蜜罐文章,引起了业界不小的轰动,但是...

3169
来自专栏ThoughtWorks

DDD实战篇:分层架构的代码结构

不同于其它的架构方法,领域驱动设计DDD(Domain Driven Design)提出了从业务设计到代码实现一致性的要求,不再对分析模型和实现模型进行区分。也...

2704
来自专栏生信宝典

小学生都学Python了,你还不知道怎么开始

最近Python又火了一把,一是我大山东省小学六年级的教材中加入了Python的内容;二是从2018年起,Python也将成为浙江高考的内容之一;三是计算机二级...

1809
来自专栏程序人生

谈谈分布式系统

先讲个黑色笑话: 半年前,一个谁也没见过的日本浪人推出的理财产品突然在七侠镇火爆起来,据说买上点屯着,不出几月就能把同福客栈,甚至龙门镖局都盘下。我们家小六的七...

31510
来自专栏生信宝典

测序数据可视化 (四)- Epigenomebrowser

Epigenomebrowser是华盛顿大学王艇教授团队开发的强大的基因组浏览器,可以显示Hi-C和ChIA-PET等三维基因组结构的数据。同时还带有比较多的小...

1836
来自专栏安恒信息

解析阻止机器学习的十种网络攻击

即使是瑟曦.兰尼斯特的阴谋诡计或者乔拉.莫尔蒙爵士父亲般的保护(译注:两者都是HBO剧集《权力的游戏》中的人物)也无法阻止攻击者攻破HBO的网络并窃取了1.5T...

2667
来自专栏玉树芝兰

如何高效使用外脑?

我们总有一种错觉,认为人类很聪明、很擅长记忆东西。认知科学家没有你这么乐观。他们提出人脑工作记忆容量很小,只有7±2个组块而已。

683

扫码关注云+社区