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

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

要想直接用gdcVolcanoPlot这个函数，首先你必须在你的R环境里安装GDCRNATools这个包，因为这个函数是这个包里面的。而GDCRNATools这个包有很多依赖的其他的包，安装起来比较费时费力，安装大概需要十到二十分钟，并且网速要好，装好大概有1G左右。如果你只想画一个火山图，实际上没有必要把这个R包全部安装了。有点高射炮打蚊子的感觉。
gdcVolcanoPlot这个函数，原作者在写的时候考虑的不是很周全，有些参数设置的不是很灵活。小编在使用的时候，发现了一些小问题。今天小编就会给大家展示一下，如何站在巨人的肩膀上看（niao）的更远。即使是“抄”也要“抄”出水平来。

我们上次“抄”来的gdcVolcanoPlot函数如下

代码语言：javascript复制

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))
}

我们画了所有差异表达基因的火山图

代码语言：javascript复制

load("DEGAll.rda")
#这里用到ggplot2这个包来画图
library(ggplot2)
gdcVolcanoPlot(DEGAll)

接下来我们用这个函数来画差异表达miRNA的火山图，在DEGAll.rda这套数据里面保存了两个数据框。可以通过ls()来查看。如何生成和使用rda文件可以参考R的save，load函数和 .rda文件

代码语言：javascript复制

ls()
#[1] "DEGAll"         "DEGMIR"

DEGAll里面存放的是所有mRNA差异表达分析的结果，而DEGMIR里面存放的是miRNA差异表达分析的结果。我们还是用前面定义的gdcVolcanoPlot来画miRNA的火山图

代码语言：javascript复制

gdcVolcanoPlot(DEGMIR)

我们会得到下面这张火山图，初略看上去也没啥问题。但是对于有强迫症的小编来说，总觉的哪里不太对劲。原来是这张图看上去点比较稀疏一点，让人觉得点比较小，而mRNA的火山图看上去比较密集一点。原因是人编码蛋白的基因大概有2万多，而人的miRNA大概就2600多个，点的数目本身就要少很多，所以会显得比较稀疏。

仔细研究了原作者对函数的描述，一共只有三个参数，并没有可以控制圆点大小的参数。瞬间感觉有点方

代码语言：javascript复制

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，你又得把函数全部重写一遍。

代码语言：javascript复制

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这个参数，这样就一劳永逸了。

代码语言：javascript复制

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

代码语言：javascript复制

#不指定dotsize，就用默认值0.8来绘图
gdcVolcanoPlot2(DEGMIR)
#指定了dotsize，就用指定值2来绘图
gdcVolcanoPlot2(DEGMIR,dotsize=2)

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

参考文献：

R函数
R的save，load函数和 .rda文件
R函数不会写，"抄"总会吧！

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

r语言 nat NAT网关

0 人点赞