【直播】我的基因组46: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-22

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏大数据

数据挖掘敲门砖-Python爬虫入门

WHAT ? 数据挖掘是一门综合的技术,随着Ai的兴起,在国内的需求日渐增大。 数据挖掘的职业方向通常有三个,顺便概要地提一下所需的技能(不仅于此) 数据分析方...

1869
来自专栏大数据文摘

解密千万密码:透过密码看人性

1232
来自专栏玉树芝兰

Python编程遇问题,文科生怎么办?

敲黑板了啊,答疑时间到。如果你没有良好的Python编程基础,在尝试应用数据科学方法时遇到了问题和困难,又不知道该如何有效解决,那么这篇文章就是为你写的。请务必...

542
来自专栏吉浦迅科技

Gromacs也开始支持OpenCL啦!

GROMACS 是目前最常用的分子动力学开源软件。主要用于蛋白、高分子化学和碳纳米管模拟。 荷兰一家OpenCL技术服务公司StreamComputin...

2816
来自专栏DT数据侠

4500个热门景点数据,告诉你国庆长假的正确打开姿势

国庆出游,确实是个让人头痛的问题。今天这位数据侠,不仅用数据告诉你国庆如何成功避开“people mountain people sea”,还手把手带你用Pyt...

670
来自专栏CDA数据分析师

Python告诉你:单词软件火了,但真的有那么多人在背单词吗?

0x00 前言 你想知道背单词软件有大概多少人注册第一天都没有背完嘛? 你想知道背单词软件这么火,这么多人在使用,真的有多少人真的在背诵嘛? 别急,Python...

1787
来自专栏程序员的碎碎念

数据挖掘敲门砖--Python爬虫入门

? WHAT 数据挖掘是一门综合的技术,随着Ai的兴起,在国内的需求日渐增大。 数据挖掘的职业方向通常有三个,顺便概要地提一下所需的技能(不仅于此) 数据分...

4398
来自专栏CDA数据分析师

大数据告诉你:土豪们都用哪些密码?

【摘要】你的密码为什么老被盗?土豪们都喜欢用哪些密码? 对于密码,我们已经知道了不少。比如,多数密码短小、简单、且容易破解。但我们对一个人选择某个密码的心理原因...

1795
来自专栏Crossin的编程教室

【我问Crossin】Python 入门之后难以提高,该如何解决?

报错 EOL 大多都是因为代码中的引号没有成对。或者其中有引号被转义,导致没起到引号的作用。

36913
来自专栏我和PYTHON有个约会

大牧夜话——爬虫篇-预告片PYTHON爬虫-江湖夜话

应大家的要求,最近打算整理一下PYTHON爬虫的东东,希望能对入门的童鞋们有所助益!本人技术一般水平有限,如有不妥请联系或者私信本人,互相进步。 内容会同步在...

632

扫描关注云+社区