【直播】我的基因组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 "$tmpt$_"}' 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
图文编辑:吃瓜群众
- Hadoop数据分析平台实战——010hadoop介绍安装
- Python为什么文件运行和在命令行运行同样语句但结果却不同?
- HDU 2034 人见人爱A-B
- Hadoop数据分析平台实战——030Hadoop Shell命令02(熟悉linux跳过)离线数据分析平台实战——030Hadoop Shell命令02
- Hadoop数据分析平台实战——070深入理解MapReduce 02(案例)离线数据分析平台实战——070深入理解MapReduce 02
- Hadoop数据分析平台实战——040HDFS介绍(熟悉基础概念跳过)离线数据分析平台实战——040HDFS&JAVA API(熟悉基础概念跳过)
- Hadoop数据分析平台实战——060深入理解MapReduce 01(案例)离线数据分析平台实战——060深入理解MapReduce 01(案例)
- 利用向量积(叉积)计算三角形的面积和多边形的面积
- HDU 1556 Color the ball
- Hadoop数据分析平台实战——080HBase介绍和安装离线数据分析平台实战——080HBase介绍和安装
- Hadoop数据分析平台实战——130Hive Shell命令介绍 02(熟悉Hive略过)离线数据分析平台实战——130Hive Shell命令介绍 02(熟悉Hive略过)
- ECJTUACM16 Winter vacation training #4 题解&源码
- Hadoop数据分析平台实战——090HBase shell客户端和Java Api介绍离线数据分析平台实战——090HBase shell客户端和Java Api介绍
- Hadoop数据分析平台实战——140Hive函数以及自定义函数讲解离线数据分析平台实战——140Hive函数以及自定义函数讲解
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- python 自动化测试(1):获取验证码图片,实现自动登录
- RocketMQ学习第一步之源码构建
- python 库学习之:openpyxl
- python 学习之:读取xml配置文件
- 我的C语言入门笔记~!
- 宜信OCR技术探索之版面分析业务实践|技术沙龙直播速记
- python 学习之:将字符串转换成变量,调用该变量实例对象的方法
- iOS 登录接口封装实践
- 机器学习模型部署—PMML
- 「文档数据库之争」MongoDB和CouchDB的比较
- 【Python】【爬虫】最近想买电脑,用Python爬取京东评论做个参考
- 新手如何快速入门Python
- 10行Python代码自动清理电脑内重复文件,解放双手!
- 模型评价指标—KS
- python 学习之:修饰器