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 条评论
登录 后参与评论

相关文章

来自专栏生信技能树

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

昨天我们学习了正常情况下,6种SNV(C>A, C>G, C>T, T>A, T>C, T>G)突变频谱的制作,但是如果考虑到突变的上下文,就可以变成96种(如...

3087
来自专栏JadePeng的技术博客

10个有用的jquery 图片插件

jquery的灵活性为我们所熟知并热爱。 给人映象最深的jquery应用通常与图片相关。 事实上,你可以借助jquery来处理图片达到给你的项目增添令人惊奇的功...

2484
来自专栏JackieZheng

AngularJS in Action读书笔记1——扫平一揽子专业术语

前(fei)言(hua):   数月前,以一个盲人摸象的姿态看了一些关于AngularJS的视频书籍,留下了我个人的一点或许是指点迷津或许是误人子弟的读后感。自...

1707
来自专栏HT

拓扑图弹力布局呈现Flickr图片搜索结果

十年前有值得分享的图片我都存在Flickr上,可惜yahoo收购了Flickr之后堕落​好多年,最近yahoo在梅姐带领下Flickr团队终于恢复了生机,个人免...

1779
来自专栏非著名程序员

Android应用程序优化注意事项

? 我们在开发过程中,如果不注意性能的优化,代码的优化等等,可能会导致应用程序的卡顿和效率极慢,所以开发过程中,注意细节,注意代码的编写和变量,常量的使用,可...

17810
来自专栏编程

实现微信朋友圈所有动态点赞的自动化用例

本人在是呀UiAutomator的过程中,突发奇想,写一个自动给朋友圈点赞的用例,经过尝试,终于成功,效果不错。这个方法用的是for循环,也可以用while循环...

26210
来自专栏Golang语言社区

使用Redis之前5个必须了解的事情

使用Redis开发应用程序是一个很愉快的过程,但是就像其他技术一样,基于Redis的应用程序设计你同样需要牢记几点。在之前,你可能已经对关系型数据库开发的那一整...

34410
来自专栏技术分享

聊下git pull --rebase

有一种场景是经常发生的。 大家都基于develop拉出分支进行并行开发,这里的分支可能是多到数十个。然后彼此在进行自己的逻辑编写,时间可能需要几天或者几周。在这...

1897
来自专栏猿人谷

C++ STL编程轻松入门基础

C++ STL编程轻松入门基础 1 初识STL:解答一些疑问 1.1 一个最关心的问题:什么是STL 1.2 追根溯源:STL的历史 1.3 千丝万缕的联系 ...

1869
来自专栏喔家ArchiSelf

一行Python代码

自从08年接触Python,就有爱不释手的感觉,逐渐地,有些不忍地疏远了Perl 和Shell编程,因为python 的优雅么? 不全是,主要是可以高效开发吧。

754

扫描关注云+社区