SNV突变(96种)频谱的制作

昨天我们学习了正常情况下,6种SNV(C>A, C>G, C>T, T>A, T>C, T>G)突变频谱的制作,但是如果考虑到突变的上下文,就可以变成96种(如下图所示)!(如果你还不了解mutation siganures,请自行复制http://www.bio-info-trainee.com/1619.html或查看原文)

The mutational spectrum of a set of SNVs was determined by classifying all SNVs contained in the set by their type of mutation (C . A, C . G, C . T, T . A, T . C, T . G) and the sequence context (i.e., the preceding and the following base). The resulting count matrix with dimensions 4 · 4 · 6 (with each cell representing a mutation of one base triplet into another) was then normalized for the observed frequency of each source base triplet in the genome that the calls were made against. An additional conversion into percentage was performed to allow for comparison of SNV sets with different sizes.

这里我们可以自己小脚本来做,也可以直接使用程序,我这里还是用号称可以替代生物学工程师的强大的bedtools软件。

http://bedtools.readthedocs.io/en/latest/content/tools/getfasta.html

简单的阅读该软件的说明书就知道了,需要把vcf文件转为3列的bed格式,就染色体号,起始终止坐标即可!

cat realign.vcf |grep -v INDEL |grep -v "^#" |cut -f 1-5|grep -v "," |perl -alne '{print join("\t",$F[0],$F[1]-2,$F[1]+1,"$F[3]:$F[4]")}' >vcf.bed

注意VCF文件的坐标不仅仅是上下文3个碱基,起始坐标应该左移,这是bed文件的特性,从0开始的!

还有第4列是突变形式,在下面的bedtools里面会用得到!

然后调用bedtools即可,代码如下:

~/biosoft/bedtools/bedtools2/bin/bedtools getfasta -fi ~/reference/genome/human_g1k_v37/human_g1k_v37.fasta -bed vcf.bed -tab -name -fo vcf.fasta

结果显示共4131526行数据!

这个结果可以用一个网页工具来检查一下:http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chr1:10248,10251

接下来我们就是要对上面的四百多万行数据进行统计咯,左边一列是突变形式,右边是上下文,我们还是跟上一讲一样,突变形式只考虑6种!【直播】我的基因组 45:SNV突变(6种)频谱的制作

代码如下:

perl -alne '{$tmp=$_;s/A:C/T:G/; s/A:T/T:A/; s/A:G/T:C/; s/G:A/C:T/; s/G:C/C:G/; s/G:T/C:A/; print "$tmp\t$_"}' vcf.fasta |cut -f 3,4 |sort |uniq -c >96context.counts

结果如下:

可视化如下,其实应该有更好的展现方式的,而且我的代码稍微有点复杂了:

dat=read.table('96context.counts')colnames(dat)=c('counts','mut','context')dat$percent = 100*dat$counts/sum(dat$counts)library(ggplot2)p=ggplot(dat,aes( x=1:nrow(dat),y=percent,fill=mut))+geom_bar(stat="identity")p=p+ylab('Mutation type probabilty')+ xlab('context')+ggtitle("Mutation signature")p=p+scale_x_continuous(breaks=1:192,labels = dat$context, expand = c(0, 0))+scale_y_continuous(expand = c(0, 0))p=p+theme_set(theme_set(theme_bw(base_size=20)))p=p+theme(text=element_text(face='bold'),axis.text.x=element_text(angle=30,hjust=1,size =6),plot.title = element_text(hjust = 0.5) ,panel.grid = element_blank(),#panel.border = element_blank())print(p)

http://www.cookbook-r.com/Graphs/Axes_(ggplot2)/

http://stackoverflow.com/questions/40675778/center-plot-title-in-ggplot2

http://stackoverflow.com/questions/22945651/how-to-remove-space-between-axis-area-plot-in-ggplot2

如果要自己写脚本,请参考生信技能树论坛,我发的帖子:

http://www.biotrainee.com/thread-666-1-1.html

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

http://cancer.sanger.ac.uk/cosmic/signatures

文:Jimmy

图文编辑:吃瓜群众

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

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

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏林德熙的博客

win10 支持默认把触摸提升鼠标事件 打开 Pointer 消息

在 WPF 经常需要重写一套触摸事件,没有UWP的Pointer那么好用。 如果一直都觉得 WPF 的触摸做的不好,或想解决 WPF 的触摸问题,但是没有方法,...

871
来自专栏ASP.NET MVC5 后台权限管理系统

ASP.NET MVC5+EF6+EasyUI 后台管理系统(45)-工作流设计-设计步骤

步骤设计很重要,特别是规则的选择。 我这里分为几个规则 1.按自行选择(在起草时候自行选审批人,比较灵活) 2.按上级(无需指定,当时需要知道用户的上司是谁,可...

2097
来自专栏鹅厂优文

开发效率太低?您可能没看这篇文章

还记得刚参加工作的时候, 有位开发的同事软件使用效率奇高. 我曾亲眼目睹他在几秒之内打开开发软件, 优雅地调出隐藏的功能, 输入数据输出结果的过程行云流水, 一...

5302
来自专栏恰同学骚年

Unity3D游戏开发初探—1.跨平台的游戏引擎让.NET程序员新生

  Unity是由Unity Technologies开发的一个让轻松创建诸如三维视频游戏、建筑可视化、实时三维动画等类型互动内容的多平台的综合型游戏开发工具,...

593
来自专栏张戈的专栏

教你如何去掉友荐和无觅的隐藏外链和版权链接

最近,为了做无觅的 APP,装上了无觅的相关推荐,结果果断不给力,打包了 2 个星期还在打包,还能再坑点么? ? 蜗牛般的效率,暂且就不吐槽了。偶尔用站长工具检...

4148
来自专栏Android-JessYan

MVPArms官方首发一键生成组件化,体验纯傻瓜式组件化开发

原文地址: https://www.jianshu.com/p/2452ea776a45

923
来自专栏张戈的专栏

替换WordPress默认搜索为百度站内搜索(知更鸟主题可照搬)

今天,中国博客联盟 QQ 群里的【58 说】博友提到百度站长平台推出绿色收录通道了。连忙登陆站长平台看了下,意外的发现张戈博客已开通了站内搜索功能。之前确实给管...

3784
来自专栏姬小光

姬小光前端兴趣班【第008期】- 真正的切图大法

上一期我们使用了 word 转存大法来生成网页,主要是为了帮助大家理解表格布局的原理。那么今天我们就来学习一下真正的切图大法。

762
来自专栏用户2442861的专栏

跟我一起学习VIM - The Life Changing Editor

前两天同事让我在小组内部分享一下VIM,于是我花了一点时间写了个简短的教程。虽然准备有限,但分享过程中大家大多带着一种惊叹的表情,原来编辑器可以这样强大,这算是...

862
来自专栏Alice

ios 浅谈一下UITextFiled UITextView 在tableview的cell上边展示

最近在项目中。要做到在tableview的cell上边加一个输入框。允许用户输入。 1.我首先选的是在uitextView  然后在通知键盘出现的时候,将tab...

1875

扫码关注云+社区