专栏首页用户7627119的专栏R函数,如何“抄”出水平

R函数,如何“抄”出水平

前面给大家介绍了,自己不会写R函数如何去“抄”高手写好的函数,我们直接“拿来”用就可以了。有读者反映为什么不直接用gdcVolcanoPlot这个函数,既然人家都已经写好了。这是一个很好的问题,这里我解答一下。原因有两个

  1. 要想直接用gdcVolcanoPlot这个函数,首先你必须在你的R环境里安装GDCRNATools这个包,因为这个函数是这个包里面的。而GDCRNATools这个包有很多依赖的其他的包,安装起来比较费时费力,安装大概需要十到二十分钟,并且网速要好,装好大概有1G左右。如果你只想画一个火山图,实际上没有必要把这个R包全部安装了。有点高射炮打蚊子的感觉。
  2. 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)

这个点的大小看上去就好多了

参考文献:

  1. R函数
  2. R的save,load函数和 .rda文件
  3. R函数不会写,"抄"总会吧!

关注“生信交流平台”公众号,后台回复"火山图",获取DEGAll.rda文件。

本文分享自微信公众号 - 生信交流平台(gh_d04ce007f7b8),作者:生信交流平台

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-10-07

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • R函数不会写,"抄"总会吧!

    前面我们简单的介绍了R函数。有些人可能会说,我现在的R水平有限,还不足以写出很高级的函数,该怎么办?俗话说前人栽树后人乘凉,他山之石可以攻玉,鲁迅同志也...

    生信交流平台
  • 小白学数据 | 28张小抄表大放送:Python,R,大数据,机器学习

    大数据文摘
  • 「R」如何计算几何平均数

    刚遇到一个有意思的问题,如何用R计算几何平均数。如果数字少,简单,计算很容易,直观上,先用prod函数连乘,然后开方即可。

    王诗翔呀
  • 如何搞定数据库水平切分?

    本文将介绍数据库架构设计中的一些基本概念,常见问题以及对应解决方案,为了便于读者理解,将以“用户中心”为例,讲解数据库架构设计的常见玩法。

    用户1737318
  • 增长乏力、违法成本低,互联网抄袭为何屡禁不止?

    今天,尽管国人对于知识产权越发重视,但是互联网抄袭之风却从未真正停止过,甚至还有愈演愈烈之势。日前,老牌K歌APP唱吧被媒体曝出抄袭唱鸭独家研发的弹唱功能,再次...

    刘旷
  • 如何看出一个程序员的技术水平?

    如何看出一个程序员的技术能力和水平?这个题目是比较复杂的,它包含的东西比较多,认真讨论能写几万字。如果是专业研究,能写一本书了。

    Java技术栈
  • 【译文】怎样学习R(下)

    何品言翻译,广东科技学院大学生,喜欢R语言和数据科学。 王陆勤审核,从事数据挖掘工作,专注机器学习研究与应用。 英文链接:http://www.r-blogge...

    小莹莹
  • 学编程一开始就值得坚持的习惯

    学习任何一门技术,在一开始就养成优秀的习惯,这是非常重要的。 1 看官方文档 遇到不清楚或不懂的知识点,先去看官方文档! 很多官方文档是英文的,硬着头皮也要看...

    老九君
  • 【R语言经典实例8】如何定义一个R函数。

    使用关键字function,并在其后跟随函数参数列表和函数主体。其基本形式如下: function(param1, ...., paramN) expr

    统计学家
  • R问题|如何查看函数的源代码

    最近有读者问我,如何查看R语言某包中某函数的源代码呢?我第一时间给出了自己比较常用的方法(见方法一),今天打算做个这方面的推文,于是又查了些资料,才发现原来水好...

    庄闪闪
  • 两个神奇的R包介绍,外加实用小抄

    认识Tidy Data1.Reshape Data2.Handle Missing Values3.Expand Tables4.split cells一、测...

    生信技能树
  • 零基础学编程039:生成群文章目录(2)

    每个月的月底,“分享与成长群”要汇总所有成员的原创文章,这次我改用了水滴微信平台把数据采集到一个电子表格文件中。在《零基础学编程019:生成群文章目录》这一节里...

    申龙斌
  • 我的5年Python7年R,述说她们的差异在哪里?

    首次接触R语言是在2012年读研的时候,有一门课程是统计分析与R语言,清晰地记得期末考试时,由于把答案给同学抄,最终落了个重考的后果(重考92分)。那个时候真的...

    1480
  • 如何通过外表看出应聘的程序员的水平?

    在软件行业混的时间比较长,前前后后面试过多个程序员,程序员的技术面试很难以貌取人,技术属于内在的东西不是靠长相,真要看外貌是不是能够是个做编程的料,倒是有个比较...

    程序员互动联盟
  • 新手绘图一站式R包ggstatsplot

    但绝大部分小伙伴仍然是选择躺平,不愿意动手实战,提高自己。对这样的小白来说,各种拥有操作界面的软件可能是更适合,比如orgin和prism等等,其实R里面也有类...

    生信技能树
  • CV学习笔记(十五):直线检测

    在这一篇文章中我们将学习使用OpenCV中的 HoughLines 函数和 HoughLinesP 函数来检测图像中的直线.

    云时之间
  • typescript高级用法之infer的理解与使用

    以前一直不会用infer,要么直接就是returnType,压根不需要用infer,网上那些教程只给示例不给具体场景就无法让人很好理解这玩意。

    徐小夕
  • java有参构造函数如何输出

    不知道你说的什么意思,你建几个public变量给它存起来不就能用了么。如果是序列化,那么把所有属性序列化就可以了,传递的什么参数该做的改变它都做掉了,序列化不需...

    用户7886150
  • 新手绘图一站式R包之ggpubr

    但绝大部分小伙伴仍然是选择躺平,不愿意动手实战,提高自己。对这样的小白来说,各种拥有操作界面的软件可能是更适合,比如orgin和prism等等,其实R里面也有类...

    生信技能树

扫码关注云+社区

领取腾讯云代金券