【直播】我的基因组 45:SNV突变(6种)频谱的制作

突变频谱呢,就是对含有SNV的VCF格式的文件进行一个统计。

全基因组SNP突变可以分成6类(C>A, C>G, C>T, A>C, A>G, A>T)。肯定会有人问为什么是六类?

以A:T>C:G为例,此种类型SNP突变包括A>C和T>G。由于测序数据既可比对到参考基因组的正链,也可比对到参考基因组的负链,当T>C类型突变出现在参考基因组正链上,A>G类型突变即在参考基因组负链的相同位置,所以将T>C和A>G划分成一类,换句话说我们只考虑正链的突变形式,参考碱基只允许有C或者T,因为它们等价于G或者A。所以全基因组SNP突变可以分成这6类。

很明显,我们只需要考虑VCF文件的第4,5行即可!

cat realign.vcf |grep -v INDEL |grep -v "^#" |cut -f 1-5 |head

cat realign.vcf |grep -v INDEL |grep -v "^#" |cut -f 4,5|sort |uniq -c |grep -v ","

我们过滤掉了多种变异形式的SNV,比如T,突变成G或者C!最后的结果如下:

一般来说,是要可视化一下的,我用R语言的ggplot来画一个呗~

  1. dat <- data.frame(type=c('C>A(G>T)','C>T(G>A)','C>G(G>C)','T>A(A>T)','T>G(A>C)','T>C(A>G)'),
  2. counts=c(180515+181567,698322+697568,184176+185144,148387+148580,177215+177415,676816+675821)
  3. )
  4. library(ggplot2)
  5. p=ggplot(dat,aes( x=type,y=counts))+geom_bar(stat="identity")
  6. print(p)

当然,mutation spectrum这个画图代码只能出一个最简单的条形图,如果你想达到下面的效果,需要学习ggplot啦!

画条形图请参考:http://docs.ggplot2.org/0.9.3.1/geom_bar.html

如果要区分染色体,可以重新考虑第1行,拿去可视化!

http://www.bio-info-trainee.com/1619.html

文:Jimmy

图文编辑:吃瓜群众

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

原文发表时间:2017-01-16

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏生信宝典

分子对接简明教程 (4)

文件格式解释 PDB文件 (详细格式描述) 基本信息部分 HEADER记录: 包括分子的分类、提交日期、PDB ID TITLE记录: 为该结构的描述,如果有多...

2837
来自专栏Android小菜鸡

Andorid pcm转码wav

参考文章:https://blog.csdn.net/hesong1120/article/details/79043482

1762
来自专栏生信宝典

如何获取目标基因的转录因子(上)——Biomart下载基因和motif位置信息

科研过程中我们经常会使用Ensembl(http://asia.ensembl.org/index.html) 网站来获取物种的参考基因组,其中BioMart工...

6943
来自专栏Python小屋

Win10系统配置Python3.6+OpenGL环境详细步骤

1、首先登录https://www.opengl.org/resources/libraries/glut/,下载下图箭头所指的文件 ? 2、解压缩,如下图所示...

3787
来自专栏吉浦迅科技

DAY46:阅读Surface Reference API

reads the CUDA array bound to the one-dimensional surface reference surfRef usin...

1015
来自专栏哈雷彗星撞地球

Objective-C 中如何测量代码的效率背景

因此,我们不可避免的要用到一些方法来计算代码的执行效率。计算代码的执行效率可以使用的API有:

945
来自专栏高爽的专栏

登录之验证码

       产生验证码,MakeCertPic.java: public class MakeCertPic { // 验证码图片中可以出现的字符集,可...

2750
来自专栏生信技能树

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

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

4268
来自专栏生信技能树

转录组数据拼接之应用篇

前前后后接触了一些基因组和转录组拼接的工作,而且后期还会持续进行。期间遇到了各种各样莫名其妙的坑,也尝试了一些不同的方法和软件,简单做一个阶段性小结。上周的今天...

5136
来自专栏Y大宽

Enrichment Map User guide用户指南

http://www.baderlab.org/Software/EnrichmentMap/UserManual#rnk

2583

扫码关注云+社区