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

2020-10-23 11:08:18 浏览数 (1)

前面我们简单的介绍了R函数。有些人可能会说,我现在的R水平有限,还不足以写出很高级的函数,该怎么办?俗话说前人栽树后人乘凉,他山之石可以攻玉,鲁迅同志也提出过“拿来”主义。已经有前人,高手写出了很多很实用,很强大的R函数,你直接拿来用就可以了。如果你很好学,也可以把人家的函数源代码拿来学习,其实这也是一种学习R的很好的方法。你如果完全读懂了原作者的函数,你还可以稍作修改用作他用,甚至可以让这个函数功能更加强大。

下面给大家举个具体的例子,火山图大家可能都不陌生,是一种展示差异表达分析结果的常用可视化方式。

在R的GDCRNATools包中就内置了一个专门画火山图的函数,叫做gdcVolcanoPlot。我们有两种方法可以获取这个函数的源代码。

1.通过下面的链接来获取gdcVolcanoPlot的源代码

https://rdrr.io/bioc/GDCRNATools/src/R/gdcDEGVisulization.R

2.从Bioconductor官网上去下载这个R包的所有源代码,

http://www.bioconductor.org/packages/release/bioc/html/GDCRNATools.html

注意一定要下载tar.gz格式的文件。.zip格式的文件是windows系统下的R安装包,都是已经编译过的,你是无法看到源代码的。

解压之后你就看到所有函数的源代码

我们要找的gdcVolcanoPlot的源代码就在gdcDEGVisulization.R这个文件中。

我们照""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))
}

接下来我们来画火山图,数据是从DEGAll.rda这个文件中来,具体如何生成这个文件和如何使用这个文件可以参考R的save,load函数和 .rda文件。加载之后你就会有DEGall这个变量了,里面存放的是差异表达分析之后的结果。画火山图需要用到logFC,FDR。

代码语言:javascript复制
load("DEGAll.rda")
#这里用到ggplot2这个包来画图
library(ggplot2)
gdcVolcanoPlot(DEGAll)

你就会得到下面这张火山图,是不是很方便,不会写函数一样可以画火山图。

关注公众号,后台回复"火山图",获取DEGAll.rda文件。

Reference:

1.R的save,load函数和 .rda文件

2.R函数

0 人点赞