搜索
查看: 1772|回复: 0

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

[复制链接]

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
发表于 2017-2-25 16:02:17 | 显示全部楼层 |阅读模式
本帖最后由 zckoo007 于 2017-2-25 16:03 编辑

【直播】我的基因组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的统计结果如下:

[AppleScript] 纯文本查看 复制代码
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的统计结果如下:

[AppleScript] 纯文本查看 复制代码
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,就可以进行基本的统计啦!



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

[AppleScript] 纯文本查看 复制代码
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即可:



  • [AppleScript] 纯文本查看 复制代码
    a=read.table('type.txt',stringsAsFactors = F) ##这个type.txt文件就是上面截图的数据
    library(ggplot2)
    ggplot(a,aes(x=V5,y=V1))+geom_bar(stat = 'identity',aes(fill=paste0(V4,V2)))+facet_wrap(~V3)



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


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

文:Jimmy

图文编辑:吃瓜群众




本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:【直播】我的基因组53:几个找变异的软件的效果比较
下一篇:【直播】我的基因组55:简单的PCA分析千人基因组的人群分布
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-16 15:38 , Processed in 0.031604 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.