R函数,如何“抄”出水平
前面给大家介绍了,自己不会写R函数如何去“抄”高手写好的函数,我们直接“拿来”用就可以了。有读者反映为什么不直接用gdcVolcanoPlot这个函数,既然人家都已经写好了。这是一个很好的问题,这里我解答一下。原因有两个
- 要想直接用gdcVolcanoPlot这个函数,首先你必须在你的R环境里安装GDCRNATools这个包,因为这个函数是这个包里面的。而GDCRNATools这个包有很多依赖的其他的包,安装起来比较费时费力,安装大概需要十到二十分钟,并且网速要好,装好大概有1G左右。如果你只想画一个火山图,实际上没有必要把这个R包全部安装了。有点高射炮打蚊子的感觉。
- gdcVolcanoPlot这个函数,原作者在写的时候考虑的不是很周全,有些参数设置的不是很灵活。小编在使用的时候,发现了一些小问题。今天小编就会给大家展示一下,如何站在巨人的肩膀上看(niao)的更远。即使是“抄”也要“抄”出水平来。
我们上次“抄”来的gdcVolcanoPlot函数如下
gdcVolcanoPlot<-function (deg.all, fc = 2, pval = 0.01)
{
geneList <- deg.all
geneList$threshold <- c()
geneList$threshold[geneList$logFC > log(fc, 2) & geneList$FDR <
pval] <- 1
geneList$threshold[geneList$logFC >= -log(fc, 2) & geneList$logFC <=
log(fc, 2) | geneList$FDR >= pval] <- 2
geneList$threshold[geneList$logFC < -log(fc, 2) & geneList$FDR <
pval] <- 3
geneList$threshold <- as.factor(geneList$threshold)
lim <- max(max(geneList$logFC), abs(min(geneList$logFC))) +
0.5
volcano <- ggplot(data = geneList, aes(x = logFC,
y = -log10(FDR)))
volcano + geom_point(aes(color = threshold), alpha = 1,
size = 0.8) + xlab("log2(Fold Change)") + ylab("-log10(FDR)") +
scale_colour_manual(values = c("red", "black", "green3")) + xlim(c(-lim, lim)) +
geom_vline(xintercept = c(-log(fc, 2), log(fc, 2)), color = "darkgreen",
linetype = 3) + geom_hline(yintercept = -log(pval,
10), color = "darkgreen", linetype = 3) + theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.border = element_rect(colour = "black"),
panel.background = element_blank()) + theme(legend.position = "none") +
theme(axis.text = element_text(size = 14), axis.title = element_text(size = 16))
}
我们画了所有差异表达基因的火山图
load("DEGAll.rda")
#这里用到ggplot2这个包来画图
library(ggplot2)
gdcVolcanoPlot(DEGAll)
接下来我们用这个函数来画差异表达miRNA的火山图,在DEGAll.rda这套数据里面保存了两个数据框。可以通过ls()来查看。如何生成和使用rda文件可以参考R的save,load函数和 .rda文件
ls()
#[1] "DEGAll" "DEGMIR"
DEGAll里面存放的是所有mRNA差异表达分析的结果,而DEGMIR里面存放的是miRNA差异表达分析的结果。我们还是用前面定义的gdcVolcanoPlot来画miRNA的火山图
gdcVolcanoPlot(DEGMIR)
我们会得到下面这张火山图,初略看上去也没啥问题。但是对于有强迫症的小编来说,总觉的哪里不太对劲。原来是这张图看上去点比较稀疏一点,让人觉得点比较小,而mRNA的火山图看上去比较密集一点。原因是人编码蛋白的基因大概有2万多,而人的miRNA大概就2600多个,点的数目本身就要少很多,所以会显得比较稀疏。
仔细研究了原作者对函数的描述,一共只有三个参数,并没有可以控制圆点大小的参数。瞬间感觉有点方
Description
A volcano plot showing differentially expressed genes/miRNAs
Usage
gdcVolcanoPlot(deg.all, fc = 2, pval = 0.01)
Arguments
deg.all
a dataframe generated from gdcDEAnalysis containing all genes of analysis no matter they are differentially expressed or not
fc
a numeric value specifying the threshold of fold change
pval
a nuemric value specifying the threshold of p value
等冷静下来,发现这个函数的源代码都有了,“盘”它。原作者没有想到,我们可以改原作者的代码啊!如果说前面完全照“抄”是囫囵吞枣,那么今天我们就来细嚼慢咽。通读函数所有代码,找到了控制圆点大小的部分,就在size这里。原作者把这个参数写死了,就是0.8。当然我们可以简单粗暴的把这个0.8改大一些,比如2,完全可以实现我们想要的效果。但是这样改显得没有水平,因为下次如果你想把点改的再大一些,比如4,你又得把函数全部重写一遍。
volcano + geom_point(aes(color = threshold), alpha = 1,
size = 0.8) + xlab("log2(Fold Change)") + ylab("-log10(FDR)")
前面我们讲R函数的时候,提到过参数,以及默认参数。
我们可以重新定义一个新的画火山图的函数gdcVolcanoPlot2,增加一个新的参数叫dotsize,来控制点的大小,我们把默认值设置成0.8,把原来size=0.8的地方改成size=dotsize这个参数,这样就一劳永逸了。
gdcVolcanoPlot2<-function (deg.all, fc = 2, pval = 0.01, dotsize=0.8)
{
geneList <- deg.all
geneList$threshold <- c()
geneList$threshold[geneList$logFC > log(fc, 2) & geneList$FDR <
pval] <- 1
geneList$threshold[geneList$logFC >= -log(fc, 2) & geneList$logFC <=
log(fc, 2) | geneList$FDR >= pval] <- 2
geneList$threshold[geneList$logFC < -log(fc, 2) & geneList$FDR <
pval] <- 3
geneList$threshold <- as.factor(geneList$threshold)
lim <- max(max(geneList$logFC), abs(min(geneList$logFC))) +
0.5
volcano <- ggplot(data = geneList, aes(x = logFC,
y = -log10(FDR)))
volcano + geom_point(aes(color = threshold), alpha = 1,
size = dotsize) + xlab("log2(Fold Change)") + ylab("-log10(FDR)") +
scale_colour_manual(values = c("red", "black", "green3")) + xlim(c(-lim, lim)) +
geom_vline(xintercept = c(-log(fc, 2), log(fc, 2)), color = "darkgreen",
linetype = 3) + geom_hline(yintercept = -log(pval,
10), color = "darkgreen", linetype = 3) + theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.border = element_rect(colour = "black"),
panel.background = element_blank()) + theme(legend.position = "none") +
theme(axis.text = element_text(size = 14), axis.title = element_text(size = 16))
}
接下来我们来试试,“改”过的函数gdcVolcanoPlot2
#不指定dotsize,就用默认值0.8来绘图
gdcVolcanoPlot2(DEGMIR)
#指定了dotsize,就用指定值2来绘图
gdcVolcanoPlot2(DEGMIR,dotsize=2)
这个点的大小看上去就好多了
参考文献:
关注“生信交流平台”公众号,后台回复"火山图",获取DEGAll.rda文件。
- 百年老牌的创新之路:看可口可乐如何用AI、大数据颠覆传统营销
- 微信小游戏上线,小程序或将成为未来的营销工具
- WooCommerce 自定义商品价格显示HTML结构
- 借助Github 为第三方WordPress 主题/插件添加“自动更新”功能
- vue-cli#2.0 webpack 配置分析
- 短代码插件S-shortcodes 更新2.4版本:修复font icon的冲突问题
- 景驰落户广州 王劲称不知百度为何指控 四条回应两大疑点
- 全球各行业2020年将需要270万位数据科学家
- 解决iOS 版Safari 中浮动(float)导致页面右侧偏移的bug
- 学而思网校又玩大了:引入人工智能技术,办了一场“人机对话”英语赛事
- 全球首个机器人公民索菲亚亮相2017双12知商节 引爆全场知识产权新高潮
- WP Settings Generator:生成WordPress设置相关代码的工具
- Reactjs 入门基础(三)
- 小谈中文环境下中文排版的font-family 字体选择
- 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 数组属性和方法
- 聊聊dubbo-go的TokenFilter
- 面试 | 卡掉不少人的一道腾讯算法面试题,高手来试试?
- 面试 | 百度测试开发岗位面试题目回顾
- ESP8266简单介绍
- 基于MTCNN和MobileFaceNet实现的人脸识别
- 学习 | egg.js 从入门到精通
- 形式化分析工具AVISPA(三)学习User micro-manual of AVISPA
- 形式化分析工具AVISPA(三)2.学习User micro-manual of AVISPA
- s6中class的一些基础知识和es5语法的对比
- 在CentOS 8上使用Elastic Stack: Elasticsearch/Kibana 7.8的部署与认证配置
- 做一个简单的京东购物栏
- 解决Elasticsearch SQL命令行启动报错 ./x-pack-env: No such file or directory
- MySql 入门到精通-sql查询语句的执行过程,你真的知道吗?
- 用 float 存储金额,老板说损失从工资里扣!
- 分布式追踪实战