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

相关文章

来自专栏IT技术精选文摘

本地IDC机房数据库容灾解决方案

风险无处不在,包括自然灾害以及突发事件等,有时候我们无法预测到一些风险,比如天津港爆炸事件。IT领域也一样,总是有意想不到的事情,风险具有不可预测性,万全之策就...

942
来自专栏IT技术精选文摘

Qzone视频下载如何做到多快好省?

Qzone的日均视频播放量已经突破了10亿,其中Android端的播放量在总播放量中的占比超过70%,相比年初,播放量的增长了超过10倍。视频下载是整个视频播放...

19810
来自专栏携程技术中心

干货 | 携程高可用架构的演变和迭代——应用开发者视角

作者简介 周源,携程技术中心基础业务研发部高级研发经理,从事软件开发10余年。2012年加入携程,先后参与支付、营销、客服、用户中心的设计和研发。 前言 携程的...

2505
来自专栏Rindew的iOS技术分享

寻找成套的 App SDK 服务

1575
来自专栏企鹅号快讯

C经典类库 需要的收藏

现实中,C++的库门类繁多,解决的问题也是极其广泛,库从轻量级到重量级的都有。本文为你介绍了十一种类库,有我们常见的,也有不常见的,一起来看。 C++类库介绍 ...

1917
来自专栏腾讯大数据的专栏

全民拥抱Docker云--Lhotse系统经验分享

前言 “只要站在风口,猪也能飞起来”,这碗心灵鸡汤不知道激励了多少英雄豪杰踏上寻风口之路。而现如今,Docker这阵龙卷风呼啸来袭,更让众人生起迎风而上、直冲云...

1969
来自专栏日志易的专栏

日志易:IT 运维分析及海量日志搜索的实践之路(上)

IT运维分析(IT Operation Analytics, ITOA)是近年兴起的把大数据技术应用于分析IT运维产生的大量数据,数据来源主要有日志、网络流量、...

3211
来自专栏QQ空间开发团队的专栏

播放量突破 10 亿,Qzone 视频下载如何做到多快好省?

基于 Qzone 的日均视频播放量已增长超过 10 倍,突破了 10 亿,我们将视频下载总结为“多快好省”四个方面,以下载成功率、首次缓冲时长和缓冲概率为主要的...

5101
来自专栏即时通讯技术

移动端IM系统的协议选型:UDP还是TCP?

对于有过网络编程经验的开发者来说,使用何种数据传输层协议来实现数据的通信,是个非常基础的问题,它涉及到你的第一行代码�该如何编写。

571
来自专栏数据库新发现

Oracle11g的新特性-11g New Features

随着这几天Oracle OpenWorld大会的召开,Oracle11g的新特性越来越多的被展现出来。

674

扫描关注云+社区