【直播】我的基因组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，高下立判！

图来源于：

文：Jimmy

图文编辑：吃瓜群众