【直播】我的基因组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 条评论
登录 后参与评论

相关文章

来自专栏UAI人工智能

OpenAI/Universe-万物

2585
来自专栏SAP最佳业务实践

从SAP最佳业务实践看企业管理(128)-MM-139寄售库存的采购和消耗

MM139寄售库存的采购和消耗 寄售:在寄售处理中, 由供应商提供物料, 并将它们存储在你处。在贵公司将这些物料从寄存货物中提出之前, 该供应商一直是这些物料的...

3884
来自专栏刺客博客

SS详解(一):原理

3.9K4
来自专栏C语言及其他语言

想举办自己的编程比赛? 来这里

大家都知道C语言网(www.dotcpp.com)是一个集在线学习、训练、比赛、求职于一体的综合性编程学习网站。大家除了日常在网站上看视频学习、做题训练、写题...

1672
来自专栏程序员互动联盟

小菜学Chromium开发系列:OpenGL学习

风雨送春归,飞雪迎春到。已是悬崖百丈冰,犹有花枝俏。俏也不争春, 只把春来报。待到山花烂漫时,她在丛中笑。 这首《毛泽东诗词》·卜算子·咏...

2805
来自专栏嵌入式程序猿

怎么样用万用表检测IGBT

IGBT是变频器中非常重要的功率器件,这段时间给大家推送了多篇变频器拆解,设计的文章,而功率器件的好坏是检测中非常重要的一环,那么如何用万用表检测呢,我们以英飞...

3463
来自专栏逸鹏说道

2018年码文方向统计

昨天发图片,琢磨着发送成功了,就把图片素材删了(占空间),然后发现不能看了,今天还是发下吧:

1114
来自专栏牛客网

(19实习)51信用卡Java面经

1705
来自专栏BestSDK

一文看懂 iOS 11所有新功能:文件管理器、多任务处理、全新文件 App等

全新文件 App 通过这个App,可以把各类文件汇集在一起,浏览、搜索和整理文件。用户最近使用过的各种文件都会在一个专属的位置,以方便查看。除了 iPad 上的...

3228
来自专栏Y大宽

Cytoscape插件5:DisGeNET(2)应用实例

广泛的定义为“对于一个给定的基因集来说,哪些基因和某疾病相关?” 哮喘病是一个严重的慢性疾病,有时会威胁生命,这种病肺部发生病变。哮喘的病因很复杂包括遗传和环...

3003

扫码关注云+社区