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

相关文章

来自专栏生信技能树

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

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

2785
来自专栏新智元

TensorFlow正式发布1.5.0,支持CUDA 9和cuDNN 7,双倍提速

来源:Github 编译:费欣欣 【新智元导读】TensorFlow今天正式发布了1.5.0版本,支持CUDA 9和cuDNN 7,进一步提速。并且,从1.6版...

2966
来自专栏炉边夜话

利用Oprofile对多核多线程进行性能分析

在对应用程序不断调优的过程中,除了制定完备的测试基准(Benchmark)外,还需要一把直中要害的利器——性能分析工具。

632
来自专栏Y大宽

ChIP-seq详细分析流程

参考生信技能树1以及生信技能树2 只记录从数据下载,到最终结果展示,具体生物学知识请自行查阅 稍后关于ChIP-seq的背景知识我会再发布一篇文章。 数据...

831
来自专栏向治洪

ARKit 简介

ARKit 简介 苹果在AR一直布局VR,最近的苹果开发者大会上,果家终于放出大招:iOS移动端ARKit平台以及VR兼容新桌面操作系统macOS High S...

2206
来自专栏Java成长之路

Solr理论基础

传统数据库是为了解决结构化存储而产生的,如关系型数据库、键值存储、操作磁盘文件的map-reduce(映射-规约)引擎,图引擎等。 传统型数据库的缺点:

773
来自专栏机器之心

腾讯Angel 1.0正式版发布:基于Java与Scala的机器学习高性能计算平台

机器之心报道 Tencent 深度学习是近些年来人工智能技术发展的核心,伴随而来的机器学习框架平台也层出不穷。到现在,一家科技巨头没有一个主导的机器学习平台都不...

2875
来自专栏大数据挖掘DT机器学习

小白一天之内玩转机器学习!

很多朋友都对机器学习心存各种敬畏之心。实际上,机器学习更多的也不过是我们“统计学习”的扩展延伸和行业实现的具体化。无非是通过样本数据发现规律性的东西...

3469
来自专栏生信技能树

找个motif嘛,简单

2387
来自专栏数据科学与人工智能

【陆勤践行】小白一天之内玩转机器学习!

很多朋友都对机器学习心存各种敬畏之心。实际上,机器学习更多的也不过是我们“统计学习”的扩展延伸和行业实现的具体化。无非是通过样本数据发现规律性的东西而已。何况“...

1775

扫码关注云+社区