【直播】我的基因组 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 条评论
登录 后参与评论

相关文章

来自专栏生信技能树

bedtools 用法大全(一文就够吧)

bedtools等工具号称是可以代替普通的生物信息学数据处理工程师的!我这里用一个专题来讲解它的用法,其实它能实现的需求,我们写脚本都是可以做的,而且我强烈建议...

1.6K6
来自专栏Y大宽

Enrichment Map User guide用户指南

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

3373
来自专栏Android小菜鸡

Andorid pcm转码wav

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

2762
来自专栏生信技能树

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

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

4458
来自专栏简书专栏

Python数据分析及可视化-小测验

本文中测验需要的文件夹下载链接: https://pan.baidu.com/s/1OqFM2TNY75iOST6fBlm6jw 密码: rmbt 下载压缩包...

3132
来自专栏Python

Python任务调度模块 – APScheduler,Flask-APScheduler实现定时任务

  看代码,定义一个函数,然后定义一个scheduler类型,添加一个job,然后执行,就可以了,代码是不是超级简单,而且非常清晰。看看结果吧。

1K0
来自专栏利炳根的专栏

学习笔记TF065: TensorFlowOnSpark

Hadoop生态大数据系统分为Yam、 HDFS、MapReduce计算框架。TensorFlow分布式相当于MapReduce计算框架,Kubernetes相...

1K0
来自专栏生信技能树

转录组数据拼接之应用篇

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

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

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

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

1295
来自专栏高爽的专栏

登录之验证码

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

2790

扫码关注云+社区

领取腾讯云代金券