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

相关文章

来自专栏生信技能树

lncRNA实战项目-第三步-了解参考基因组及注释文件

下载原始测序数据: 在GEO数据库搜索GSE87182, 这里没有直接给出ftp地址,需要先从BioProject找到SRA号,可以得到RNA-Seq的SRA的...

5125
来自专栏肖洒的博客

【爬虫】(七)Python数据存储之MySQL(下)

上一篇关于Python和MySQL的简单联调做了学习。 这次主要是将这个过程再优化扩大点。 对教务处需要的数据都进行了处理存进数据库了。 也是对bug问题的总结...

841
来自专栏木宛城主

SharePoint 2013 Designer工作流——Parallel Block的应用

参考目录 安装和配置SharePoint 2013 Workflow SharePoint 2013 实现多级审批工作流 在自定义Workflow...

24010
来自专栏张善友的专栏

Entity Framework Code First 支持存储过程

存储过程(Stored Procedure)不仅仅是将多得简直荒唐的业务逻辑塞入数据库的一种方式;它还是避免将多得简直荒唐的存储逻辑塞入应用程序层(applic...

1988
来自专栏生信技能树

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

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

3838
来自专栏北京马哥教育

余生只够写50行代码,这么写绝对赚翻了

学Python最简单的方法是什么?推荐阅读:Python开发工程师成长魔法 假如有一天死神来找你,警告你最多只能再写50行代码,然后就得随他而去,应该写点什么...

3158
来自专栏张善友的专栏

在MongoDB中实现聚合函数

随着组织产生的数据爆炸性增长,从GB到TB,从TB到PB,传统的数据库已经无法通过垂直扩展来管理如此之大数据。传统方法存储和处理数据的成本将会随着数据量增长而显...

2727
来自专栏积累沉淀

linux学习之硬盘的存储原理和内部架构

首先,让我们看一下硬盘的发展史: 1956年9月13日,IBM的IBM 350 RAMAC(Random Access Method of Accounting...

3476
来自专栏施炯的IoT开发专栏

PhoneFinder--寻找丢失的手机

    手机丢了怎么办?那就打电话给手机,如果运气好的话,捡到的好心人能够把手机还给你。如果手机是被偷的,那就没有办法了,即使手机开着,估计小偷也不会接电话。当...

2834
来自专栏芋道源码1024

中国程序员容易发音错误的单词

1133

扫码关注云+社区