【直播】我的基因组68:看看哪些基因的突变较多,哪些较少

全基因组分析后的vcf突变文件记录了四百多万个位点,前面我们讲到了如何把它们注释到dbSNP数据库ID,一般来说有注释的位点也就顺便注释到了基因,所以可以简单写一个程序来看看哪些基因的突变位点最多:

cat autochr.highQuali.dbsnp.vcf |perl -alne '{/GENEINFO=(.*?)[;,\|]/;$h{$1}++ if $1 }END{print "$_\t$h{$_}" foreach keys %h}' >~/tmp/vcf.gene.stat

简单统计结果如下;

当然,其实并不需要注释到dbSNP数据后再进行统计,可以直接对vcf文件进行基因注释,因为vcf文件本身就记录着坐标,把vcf文件按照bed格式稍微转换一下,就可以用bedtools来进行注释啦。

首先制备好基因的坐标文件,染色体号,基因定位的起始终止坐标即可,比如下面这个SPIN1基因:

可以看到, 有10个突变位点注释到了这个基因,可以其中只有4个是dbSNP数据库记录的,所以最开始统计的基因的突变个数排行不是很准确。

那么我们就针对这个bedtools closest 结果进行统计吧:

可以看到几乎每个基因的突变个数都增加了,因为不需要被dbSNP数据库收录啦。

再看看基因突变个数的个数的变化:

之前突变个数为1的那些基因有1324个,但是现在只剩下了712个!

如果可视化在图上,就能看到条形图很明显的右移,不过这不是重点。

FAQ

为什么有712个基因仅仅发现一个突变呢?是这个基因太保守了吗?还是这个基因长度太短了?同理,那些突变异常多的基因又有什么特征呢?

我选取了那712个只有一个变异位点的基因,还有超过400个变异位点的909个基因。

很明显,从长度来解释是一个很好的角度~~

以上的变异位点,都应该改名叫做多态性位点。

文:Jimmy

图文编辑:吃瓜群众

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

原文发表时间:2017-04-25

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏生信技能树

比对到hg19和hg38对somatic变异的寻找影响很大

其中B是正常组织的WES数据,使用varscan找somatic mutation的时候作为normal,然后对另外两个样本(D和T)计算。 从这个bam文件可...

1213
来自专栏美团技术团队

ReactiveCocoa核心元素与信号流

概述 ReactiveCocoa(以下简称“RAC”)是一个函数响应式编程框架,它能让我们脱离Cocoa API的束缚,给我们提供另外一套编码的思路与可能性,它...

3654
来自专栏令仔很忙

设计模式六大原则---依赖倒置原则(DIP)

    依赖倒置原则告诉我们:细节是多变的,而抽象是相对稳定的。所以我们编程的时候要注重抽象的编程,而非细节编程。

682
来自专栏生信技能树

使用R语言的cgdsr包获取TCGA数据

前些天被TCGA的终结新闻刷屏,但是一直比较忙,还没来得及仔细研读,但是笔记本躺着的一些TCGA教程快发霉了,借此契机好好整理一下吧,预计二十篇左右的笔记

1263
来自专栏林冠宏的技术文章

通俗地介绍下---数据结构之堆

(出处:https://cloud.tencent.com/developer/user/1148436/activities) 前序:   堆是基础数据结构中...

1837
来自专栏落影的专栏

使用AudioToolbox编码AAC

前言 使用VideoToolbox硬编码H.264 使用VideoToolbox硬解码H.264 这次在编码H.264视频流的同时,录制并编码AAC音频流。...

3647
来自专栏点滴积累

使用Python以优雅的方式实现根据shp数据对栅格影像进行切割

一、前言        前面一篇文章(使用Python实现子区域数据分类统计)讲述了通过geopandas库实现对子区域数据的分类统计,说白了也就是如何根据一个...

56610
来自专栏DannyHoo的专栏

iOS开发中使用算法之二分搜索算法

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/u010105969/article/details/...

592
来自专栏生信技能树

【直播】我的基因组77:批量计算每个蛋白编码基因的测序深度及覆盖度

目前我使用的仍然是hg19系统的参考基因组,所以就在gencode数据库里面下载了基于hg19的gtf注释文件,并格式化如下: head ~/reference...

3728
来自专栏生信技能树

使用R包genefu来根据基因集进行表达谱分类

教程略微有点复杂:https://rdrr.io/bioc/genefu/f/inst/doc/genefu.pdf

1024

扫码关注云+社区