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

相关文章

来自专栏生信宝典

一文教会你查找基因的启动子、UTR、TSS等区域以及预测转录因子结合位点

本文授权转载自科研小助手(ID:SciRes)斜体小一号字体为生信宝典的备注或校正。

1231
来自专栏Y大宽

RNA-seq分析简洁版

Tumor:SRR316214,SRR316215 Adjacent Normal Liver:SRR316212,SRR316213

562
来自专栏生信技能树

【直播】我的基因组62:用Delly检测SV

人类单体型(Haplotype)及单核苷酸多态性位点(Single Nucleotide Polymorphism, SNP),能够揭示对药物和环境因子的个体反...

4168
来自专栏生信宝典

分子对接简明教程 (一)

分子对接(Molecular Docking)理论 所谓分子对接就是两个或多个分子之间通过几何匹配和能量匹配相互识别找到最佳匹配模式的过程。分子对接对酶学研究和...

2379
来自专栏生信宝典

不是原配也可以-对接非原生配体

Docking非原生配体 在前面的例子中,AutoDock Vina能把配体构象调整到几乎原生的构象,验证了这一预测方法的准确度。下面,我们尝试docking另...

1868
来自专栏生信宝典

生信宝典之傻瓜式(六)查找转录因子的靶基因

1917
来自专栏Y大宽

RNA-seq(10):KEGG通路可视化:gage和pathview

开始用gage包进行富集分析,gage()函数需要fold change 和Entrez gene IDs

662
来自专栏FreeBuf

利用开源工具TempestSDR实现屏显内容远程窃取

你可能还不了解“TEMPEST”,它是用来窃取远程视频信息的一种基于软件定义的无线电平台技术,可被当做间谍工具包使用,用来针对某些目标电子设备的射频信号(声音和...

1905
来自专栏生信技能树

得到一个物种所有基因的TSS(转录起始位点)区域的bed文件。

首先在UCSC的table browser 里面下载下面这个文件: ? 可以看到我这里选择的mm10的refseq系统的所有基因,共有29037个不同的tss,...

2908
来自专栏生信技能树

找个motif嘛,简单

2187

扫描关注云+社区