前面給大家介紹了，自己不會寫R函數(shù)如何去“抄”高手寫好的函數(shù)，我們直接“拿來”用就可以了。有讀者反映為什么不直接用gdcVolcanoPlot這個(gè)函數(shù)，既然人家都已經(jīng)寫好了。這是一個(gè)很好的問題，這里我解答一下。原因有兩個(gè)

1. 要想直接用gdcVolcanoPlot這個(gè)函數(shù)，首先你必須在你的R環(huán)境里安裝GDCRNATools這個(gè)包，因?yàn)檫@個(gè)函數(shù)是這個(gè)包里面的。而GDCRNATools這個(gè)包有很多依賴的其他的包，安裝起來比較費(fèi)時(shí)費(fèi)力，安裝大概需要十到二十分鐘，并且網(wǎng)速要好，裝好大概有1G左右。如果你只想畫一個(gè)火山圖，實(shí)際上沒有必要把這個(gè)R包全部安裝了。有點(diǎn)高射炮打蚊子的感覺。
2. gdcVolcanoPlot這個(gè)函數(shù)，原作者在寫的時(shí)候考慮的不是很周全，有些參數(shù)設(shè)置的不是很靈活。小編在使用的時(shí)候，發(fā)現(xiàn)了一些小問題。今天小編就會給大家展示一下，如何站在巨人的肩膀上看的更遠(yuǎn)。即使是“抄”也要“抄”出水平來。

我們上次“抄”來的gdcVolcanoPlot函數(shù)如下

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

我們畫了所有差異表達(dá)基因的火山圖

load("DEGAll.rda")
#這里用到ggplot2這個(gè)包來畫圖
library(ggplot2)
gdcVolcanoPlot(DEGAll)

接下來我們用這個(gè)函數(shù)來畫差異表達(dá)miRNA的火山圖，在DEGAll.rda這套數(shù)據(jù)里面保存了兩個(gè)數(shù)據(jù)框?？梢酝ㄟ^ls()來查看。

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

DEGAll里面存放的是所有mRNA差異表達(dá)分析的結(jié)果，而DEGMIR里面存放的是miRNA差異表達(dá)分析的結(jié)果。我們還是用前面定義的gdcVolcanoPlot來畫miRNA的火山圖

gdcVolcanoPlot(DEGMIR)

我們會得到下面這張火山圖，初略看上去也沒啥問題。但是對于有強(qiáng)迫癥的小編來說，總覺的哪里不太對勁。原來是這張圖看上去點(diǎn)比較稀疏一點(diǎn)，讓人覺得點(diǎn)比較小，而mRNA的火山圖看上去比較密集一點(diǎn)。原因是人編碼蛋白的基因大概有2萬多，而人的miRNA大概就2600多個(gè)，點(diǎn)的數(shù)目本身就要少很多，所以會顯得比較稀疏。

仔細(xì)研究了原作者對函數(shù)的描述，一共只有三個(gè)參數(shù)，并沒有可以控制圓點(diǎn)大小的參數(shù)。瞬間感覺有點(diǎn)方

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

等冷靜下來，發(fā)現(xiàn)這個(gè)函數(shù)的源代碼都有了，“改”它。原作者沒有想到，我們可以改原作者的代碼?。∪绻f前面完全照“抄”是囫圇吞棗，那么今天我們就來細(xì)嚼慢咽。通讀函數(shù)所有代碼，找到了控制圓點(diǎn)大小的部分，就在size這里。原作者把這個(gè)參數(shù)寫死了，就是0.8。當(dāng)然我們可以簡單粗暴的把這個(gè)0.8改大一些，比如2，完全可以實(shí)現(xiàn)我們想要的效果。但是這樣改顯得沒有水平，因?yàn)橄麓稳绻阆氚腰c(diǎn)改的再大一些，比如4，你又得把函數(shù)全部重寫一遍。

volcano + geom_point(aes(color = threshold), alpha = 1, 
                       size = 0.8) + xlab("log2(Fold Change)") + ylab("-log10(FDR)")

前面我們講R函數(shù)的時(shí)候，提到過參數(shù)，以及默認(rèn)參數(shù)。

我們可以重新定義一個(gè)新的畫火山圖的函數(shù)gdcVolcanoPlot2，增加一個(gè)新的參數(shù)叫dotsize，來控制點(diǎn)的大小，我們把默認(rèn)值設(shè)置成0.8，把原來size=0.8的地方改成size=dotsize這個(gè)參數(shù)，這樣就一勞永逸了。

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

接下來我們來試試，“改”過的函數(shù)gdcVolcanoPlot2

#不指定dotsize，就用默認(rèn)值0.8來繪圖
gdcVolcanoPlot2(DEGMIR)
#指定了dotsize，就用指定值2來繪圖
gdcVolcanoPlot2(DEGMIR,dotsize=2)

這個(gè)點(diǎn)的大小看上去就好多了

參考文獻(xiàn)：

1. R函數(shù)
2. R的save，load函數(shù)和 .rda文件
3. R函數(shù)不會寫，"抄"總會吧！

參考下面這篇文章，獲取DEGAll.rda文件

R函數(shù)，如何“抄”出水平?mp.weixin.qq.com