写一个自己的R函数

在实际的工作中，经常需要做大量相同的数据处理，重复而又冗长的操作很容易让人抓狂。因此，懒惰促使我学习了写一个自己的R函数，以及使用循环语句。

懒惰使人进步

（1）画一个火山图

以画差异表达基因的火山图为例，我们通过DESeq2获得了每个基因在每组比较中的log2 fold change和padj，并整理好如下：

> head(res.test)
      cMvscF.log2FC cMvscF.padj mFvscF.log2FC  mFvscF.padj mMvscM.log2FC  mMvscM.padj
gene1    -0.0776042  1.00000000     0.8235784 6.335025e-03   0.857999861 3.282493e-03
gene2    -0.3658063  0.00533996     0.7614857 1.594686e-11   0.832097420 7.320881e-14
gene3    -0.3389348  0.32209618     1.1255354 3.141754e-07   1.122211863 2.378666e-07

先画一个火山图试一试：

> #绘图数据
> volcanodata <- data.frame(LFC = res.test$cMvscF.log2FC,
+                           padj = res.test$cMvscF.padj,
+                           mean = res.test$baseMean)
> #上调下调
> volcanodata$Condition=ifelse(volcanodata$LFC>1 & volcanodata$padj<=0.1,"up",
+                              ifelse(volcanodata$LFC<(-1) & volcanodata$padj<=0.1,"down",
+                                     "normal"))
> volcanodata$Condition <- factor(volcanodata$Condition,levels=c("up","down","normal"))
> volcanodata <- volcanodata[order(volcanodata$padj,decreasing=T),]
> #数量统计
> print(table(volcanodata$Condition))
    up   down normal 
    96     19    284

> #绘图
> #library(ggplot2)
> p <- ggplot(data=volcanodata, aes(x=LFC, y=-log10(padj), colour=Condition)) + 
+   geom_point(alpha=0.8,shape=16)  +
+   xlab("log2 fold change") + 
+   ylab("-log10 padj")+
+   geom_vline(xintercept=c(1,-1),color="black",linetype=2) +
+   geom_vline(xintercept=0,color="orangered",linetype=1) +
+   geom_hline(yintercept=-log10(0.1),color="black",linetype=2) +
+   scale_color_manual(values=c('up'='red','down'='green','normal'='gray50')) +
+   theme_bw()
> #保存
> ggsave(p,filename="test.plot1.png",width=6, height=3)

得到了一张图片：

test.plot1.png

（2）写一个自己的函数

每一个R函数都包括三个部分：函数名，程序主体以及参数集合。在编写自定义R函数时，使用function函数将这三个部分各自储存在一个R对象中，形如： my_function<-function(){}。函数将返回最后一行的运行输出结果，如果最后一行不输出结果，整个函数也将不会有返回值。function()的括号内可以填入参数名称。如果只填参数名称，调用函数时不写参数（如果需要用到这个参数）会报错，因此可以预先设置一个初始默认值给括号内的参数，即缺省值。...参数是一种特殊的参数，表明一些可以传递给另一个函数的参数。常用于需要扩展另一个函数，而又不想复制原函数的整个参数列表时。

我们把第一部分的操作写成一个函数。把代码复制粘贴到大括号里，然后把需要用到的变量替换成参数名称。

plot.volcano.v1 <- function(LFC,padj,pcut=0.1,LFCcut=0,filename){
  #绘图数据
  volcanodata <- data.frame(LFC = LFC, padj = padj)
  #上调下调
  volcanodata$Condition=ifelse(volcanodata$LFC>LFCcut & volcanodata$padj<=pcut,"up",
                               ifelse(volcanodata$LFC<(-LFCcut) & volcanodata$padj<=pcut,"down",
                                      "normal"))
  volcanodata$Condition <- factor(volcanodata$Condition,levels=c("up","down","normal"))
  volcanodata <- volcanodata[order(volcanodata$padj,decreasing=T),]
  #数量统计
  print(table(volcanodata$Condition))
  #绘图
  p <- ggplot(data=volcanodata, aes(x=LFC, y=-log10(padj), colour=Condition)) + 
    geom_point(alpha=0.8,shape=16)  +
    xlab("log2 fold change") + 
    ylab("-log10 padj")+
    geom_vline(xintercept=c(LFCcut,-LFCcut),color="black",linetype=2) +
    geom_vline(xintercept=0,color="orangered",linetype=1) +
    geom_hline(yintercept=-log10(pcut),color="black",linetype=2) +
    scale_color_manual(values=c('up'='red','down'='green','normal'='gray50')) +
    theme_bw()
  #保存图片
  ggsave(p,filename=filename,width=6, height=3)
}

测试一下自定义函数：

> plot.volcano.v1(res.test$mFvscF.log2FC,res.test$mFvscF.padj,
+                 pcut=0.1,LFCcut=1,
+                 filename="test.plot2.png")
    up   down normal 
    80      8    311

test.plot2.png

把这个函数的代码复制到一个新的R Script，保存到当前工作路径下，命名为plot.volcano.v1.R。然后一键清空，source重新调用这个函数：

> rm(list=ls())
> source("plot.volcano.v1.R")
> load("res.test.RData")
> plot.volcano.v1(res.test$mMvscM.log2FC,res.test$mMvscM.padj,
+                 pcut=0.1,LFCcut=1,
+                 filename="test.plot3.png")
    up   down normal 
    93     12    294

test.plot3.png

（3）写一个循环

编程中减少代码重复的两个工具，一是循环，一是函数。有了上面那个高度个性化的自定义函数，R代码就变的简洁多了。如果想进一步解决重复操作的问题，就需要写一个循环。

for循环重复地执行一个语句，直到某个变量的值不再包含在序列seq中为止。

> for (i in 1:5) print("hello")
[1] "hello"
[1] "hello"
[1] "hello"
[1] "hello"
[1] "hello"

while循环重复地执行一个语句，直到条件不为真为止。

> i <- 5
> while (i>0) {print("hello"); i <- i-1}
[1] "hello"
[1] "hello"
[1] "hello"
[1] "hello"
[1] "hello"

if语句用来进行条件控制，以执行不同的语句。if (conditon) {expr1} else {expr2}。若condition条件为真，则执行expr1，否则执行expr2。else要紧跟着第一个语句的大括号}，不能另起一行，否则会认为if (conditon) {expr1}已经运行完了。ifelse会更加高效。

> if (i==5) {
+   print("is 5")
+ } else {
+   print("isn't 5")}
[1] "isn't 5"
> ifelse(i==1,"is 5","isn't 5")
[1] "isn't 5"

现在我们写一个for循环，把上面三个图一次性都画出来：

> for (i in 1:3) {
+   plot.volcano.v1(res.test[,(2*i-1)],res.test[,(2*i)],
+                   pcut=0.1,LFCcut=1,
+                   filename=paste0("test.plot.",i,".png"))
+ }
    up   down normal 
    96     19    284 
    up   down normal 
    80      8    311 
    up   down normal 
    93     12    294

变的简洁的代码看上去真是让一个强迫症患者感到开心呢。