从WGS测序得到的VCF文件里面提取位于外显子区域的【直播】我的基因组84

首先要下载并且得到人类基因组的外显子坐标记录文件

这里我用的参考基因组版本仍然是hg19,所以去CCDS数据库里面下载对应版本,并且格式化成BED文件。

wget  ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/archive/Hs37.3/CCDS.20110907.txtcat CCDS.20110907.txt |perl -alne '{/[(.*?)]/;next unless $1;$exons=$1;$exons=~s/\s//g;$exons=~s/-/\t/g;print "$F[0]\t$" foreach split/,/,$exons;}' >hg19exon.bed

制作好的bed格式的人类全部的exon区域坐标文件如下:

1    801942  8024331    861321  8613921    865534  8657151    866418  8664681    871151  8712751    874419  8745081    874654  8748391    876523  8766851    877515  8776301    877789  877867

从VCF文件里面根据BED文件进行抽提

这里就不自己造轮子了,用现成的工具,而且是我们用过很多次的SnpEff套件,代码如下

cat snp.vcf | java -jar  ~/biosoft/SnpEff/snpEff/SnpSift.jar intervals hg19exon.bed >hg19exon.snp.vcfcat indel.vcf | java -jar  ~/biosoft/SnpEff/snpEff/SnpSift.jar intervals hg19exon.bed >hg19exon.indel.vcf

可以把我经由GATK best practice流程得到的SNP/INDEL记录的VCF文件都进行提取,用代码 wc -l *vcf简单统计一下提取的效果,如下:

      1042 hg19_exon.indel.vcf     25067 hg19_exon.snp.vcf    754755 indel.vcf   3784343 snp.vcf

很明显可以看到,位于外显子区域的mutation毕竟是少数,这时候还可以继续看看那些在外显子上面却没有被dbSNP数据库记录的mutation还有多少:

cat hg19_exon.snp.vcf |grep -v "^#" |cut -f 3 |grep '\.' |wc

仍然有2315个SNV在外显子区域,却没有被dbSNP数据库记录,可能是我的家族特异性的位点,属于正常的基因型多样性,也有极小的可能性这些位点是后发突变,也就是通常癌症研究领域的somatic mutation 。

用下面的代码可以简单浏览一下这些在外显子上面的未知突变的情况。

cat hg19_exon.snp.vcf |perl -alne '{print if $F[2] eq "."}' |less -Scat hg19_exon.indel.vcf |perl -alne '{print if $F[2] eq "."}' |less -S

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

原文发表时间:2017-08-23

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏逸鹏说道

Uncaught RangeError: Maximum call stack size exceeded 调试日记

异常处理汇总-前端系列 http://www.cnblogs.com/dunitian/p/4523015.html 开发道路上不是解决问题最重要,而是解决问题...

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

SAP最佳业务实践:ETO–项目装配(240)-6审批 WBS 要素

image.png CJ20N审批 WBS 要素 为确保预先采购(下个步骤),WBS 要素的状态必须为 已释放。 角色项目经理 后勤®项目系统®项目®项目构造...

3106
来自专栏HBStream流媒体与音视频技术

借用PortAudio采集和播放音频,实现双路混音器

2975
来自专栏我是业余自学C/C++的

Mac部分按键失灵问题解决

1914
来自专栏HBStream流媒体与音视频技术

借助开源项目,又好又快的实现视频文件”剧情连拍(剧情截图)”功能

2927
来自专栏腾讯Bugly的专栏

WebVR如此近 - three.js的WebVR示例程序解析

关于WebVR 最近VR的发展十分吸引人们的眼球,很多同学应该也心痒痒的想体验VR设备,然而现在的专业硬件价格还比较高,入手一个估计就要吃土了。但是,对于我们...

3318
来自专栏Y大宽

RNA-seq(4):下载参考基因组及基因注释

那下载哪个基因组呢?先了解一下: https://bitesizebio.com/38335/get-to-know-your-reference-genom...

1025
来自专栏生信技能树

生信菜鸟团博客2周年精选文章集(6)三个最基础生信软件教程

其实我现在已经不写软件教程了! fastqc对原始测序reads质控 NCBI的blast++软件使用说明书 SRA工具sratoolkit把原始测序数据转为...

40511
来自专栏cnblogs

vue原来可以这样上手

       今儿与一群友讨论vue相关问题让我思量极深,1.我们是否在争对性解决问题或者说是帮助别人;2.我们是否在炫耀自己的技能。以下是被戏剧化的对白: ...

1919
来自专栏FreeBuf

黑了记者:写个恶意软件玩玩(二)

该篇是3篇系列中的第2篇(可在此读第1篇 http://www.freebuf.com/articles/others-articles/34277.html)...

1809

扫描关注云+社区