从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 条评论
登录 后参与评论

相关文章

来自专栏李智的专栏

CK+表情数据库及使用

这个数据库是在 Cohn-Kanade Dataset 的基础上扩展来的,发布于2010年。这个数据库比起JAFFE 要大的多。而且也可以免费获取,包含表情的l...

1311
来自专栏Python中文社区

Github上影响力最大的十位Pythoner

10、Shipeng Feng [1] 来自:北京市 Fllowers:213 Stared:59 代表项目:plan [2] - 一个用Python...

1996
来自专栏python小白到大牛

Python终级教程!语音识别!大四学生实现语音识别技能!吊的不行

语音识别源于 20 世纪 50 年代早期在贝尔实验室所做的研究。早期语音识别系统仅能识别单个讲话者以及只有约十几个单词的词汇量。现代语音识别系统已经取得了很大进...

2162
来自专栏吉浦迅科技

为啥在Matlab上用NVIDIA Titan V训练的速度没有GTX1080快?

在Matlab官方论坛上看到这个帖子,希望给大家带来参考 有一天,有人在Matlab的论坛上发出了求救帖: ? 楼主说: 我想要加快我的神经网络训练,所以把G...

4658
来自专栏WindCoder

Python版方正教务系统查询

准确的说是Python+wxPython版河北师大方正教务系统查询。目前仅实现了信息查询、成绩查询、平均学分绩计算。

691
来自专栏FreeBuf

利用esp8266制作一个可随身携带的WiFi密码钓鱼器

孙子兵法:故上兵伐谋,其次伐交,其次伐兵,其下攻城;攻城之法为不得已。 正是如此靠机器硬破不如直接问对方。 下面这就介绍下WIFI密码社工的进阶玩法 玩过Wi...

4264
来自专栏章鱼的慢慢技术路

用ARM实现音乐电子相册

1312
来自专栏区块链源码分析

超级账本(Hyperledger Fabric)源码分析之一:总览

1)Go,注意设置好gopath(笔者安装的是go1.8.3,对应的源码是v1.0.0这个tag,版本不对可能会出现编译不过或者运行出现问题)

3265
来自专栏Python中文社区

Python中文社区开源项目计划:ImagePy

开源图像处理框架,插件式设计,可以快速集成任何基于或支持numpy的图像处理算法,经过简单配置,快速生成交互环境,可供非计算机专业人员使用,可以理解为是算法研发...

801
来自专栏Python中文社区

pyecharts(二):Python可视化利器

專 欄 ❈陈键冬,Python中文社区专栏作者 GitHub: https://github.com/chenjiandongx ❈ 恭喜本社区专栏作者陈键冬...

2717

扫码关注云+社区