搜索
查看: 2619|回复: 1

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

[复制链接]

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
发表于 2017-2-7 20:02:36 | 显示全部楼层 |阅读模式
【直播】我的基因组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软件。

[url=]http://bedtools.readthedocs.io/en/latest/content/tools/getfasta.html[/url]

简单的阅读该软件的说明书就知道了,需要把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行数据!

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



接下来我们就是要对上面的四百多万行数据进行统计咯,左边一列是突变形式,右边是上下文,我们还是跟上一讲一样,突变形式只考虑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)




[url=]http://www.cookbook-r.com/Graphs/Axes_(ggplot2)/[/url]
[url=]http://stackoverflow.com/questions/40675778/center-plot-title-in-ggplot2[/url]
[url=]http://stackoverflow.com/questions/22945651/how-to-remove-space-between-axis-area-plot-in-ggplot2[/url]



如果要自己写脚本,请参考生信技能树论坛,我发的帖子:
[url=]http://www.biotrainee.com/thread-666-1-1.html[/url]

[url=]http://www.bio-info-trainee.com/1623.html[/url]
[url=]http://cancer.sanger.ac.uk/cosmic/signatures[/url]

文:Jimmy

图文编辑:吃瓜群众



本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:【直播】我的基因组 45:SNV突变(6种)频谱的制作
下一篇:生信人的windows电脑必须安装的软件
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复

使用道具 举报

0

主题

3

帖子

38

积分

新手上路

Rank: 1

积分
38
发表于 2017-2-12 17:03:23 | 显示全部楼层
bedtools get!
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-16 15:38 , Processed in 0.034723 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.