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

相关文章

来自专栏编程

无数据库权限下载文献攻略大全

阅读完本篇文章你就会学会了在家里,在路上,在可以连接到网络的任何地点都轻松下载您所需要的文献!!! 在之前的内容中,我们为大家推送过两篇关于如何在没有权限的情况...

25380
来自专栏Coding01

推荐一个 PHP 图像处理操作插件 Intervention Image

无论 Web 前端,还是 APP 开发,都避免不了和图像处理打交道,对于前端来说,图像处理都还好说,也比较简单。

33120
来自专栏嵌入式程序猿

mscan VS flexcan

在嵌入式程序猿的公众号里,曾多次介绍过NXP的flexcan以及基于flexcan的一些其他协议和开发,最近在用NXP的另外一款片子,使用的是mscan,这两种...

50390
来自专栏后端云

compute node ha 主流开源实现

nova evacuate和热迁移很像。都是想实例从一个节点转移到另外一个节点。区别主要是热迁移在正常状态下进行的,疏散时在异常状态下进行的。用一个形象的比如就...

22630
来自专栏生信技能树

【直播】我的基因组 43:简单粗糙的WGS数据分析流程

前面我们扯到bam文件的各种操作,vcf文件的各种操作,基础知识不牢固的同学可能已经云里雾里了。这次我们来讲一个简单的。就是拿到了fastq的测序数据,如何把全...

50790
来自专栏坚毅的PHP

zookeeper学习系列:四、Paxos算法和zookeeper的关系

一、问题起源 淘宝搜索的博客 http://www.searchtb.com/2011/01/zookeeper-research.html  提到Paxos是...

44340
来自专栏SDNLAB

如何像Facebook一样构建数据中心 – BGP在大规模数据中心中的应用(3)

作者简介:史梦晨,曾就职于国内金牌集成商, 现就职于EANTC( 欧洲高级网络测试中心),研究方向:网络架构,测试,运维(大规模数据中心,SD-WAN,EVPN...

24010
来自专栏数据魔术师

词云图:论一个精致猪猪男孩的数据修养

“词云”就是对网络文本中出现频率较高的“关键词”予以视觉上的突出,形成“关键词云层”或“关键词渲染”,从而过滤掉大量的文本信息,使浏览网页者只要一眼扫过...

17240
来自专栏生信技能树

【直播】我的基因组 31:vcf文件标记dbSNP的rsID号

vcf文件标记dbSNP的rsID号的这个问题非常多的人问过,大部分的variation calling软件给出的vcf文件里面第3列都是一个纯粹的do...

38380
来自专栏铭毅天下

吃透 | Elasticsearch filter和query的不同

除了确定文档是否匹配外,查询子句还计算了表示文档与其他文档相比匹配程度的_score。

17720

扫码关注云+社区

领取腾讯云代金券