前面我們簡單的介紹了R函數(shù)清笨。有些人可能會說，我現(xiàn)在的R水平有限刃跛，還不足以寫出很高級的函數(shù)抠艾，該怎么辦？俗話說前人栽樹后人乘涼桨昙，他山之石可以攻玉检号，魯迅同志也提出過“拿來”主義。已經(jīng)有前人蛙酪，高手寫出了很多很實用齐苛，很強大的R函數(shù)，你直接拿來用就可以了桂塞。如果你很好學(xué)凹蜂，也可以把人家的函數(shù)源代碼拿來學(xué)習(xí)，其實這也是一種學(xué)習(xí)R的很好的方法阁危。你如果完全讀懂了原作者的函數(shù)玛痊，你還可以稍作修改用作他用，甚至可以讓這個函數(shù)功能更加強大狂打。
下面給大家舉個具體的例子擂煞，火山圖大家可能都不陌生，是一種展示差異表達分析結(jié)果的常用可視化方式趴乡。
在R的GDCRNATools包中就內(nèi)置了一個專門畫火山圖的函數(shù)对省，叫做gdcVolcanoPlot。我們有兩種方法可以獲取這個函數(shù)的源代碼晾捏。
1.通過下面的鏈接來獲取gdcVolcanoPlot的源代碼
https://rdrr.io/bioc/GDCRNATools/src/R/gdcDEGVisulization.R
2.從Bioconductor官網(wǎng)上去下載這個R包的所有源代碼蒿涎，
http://www.bioconductor.org/packages/release/bioc/html/GDCRNATools.html
注意一定要下載tar.gz格式的文件。.zip格式的文件是windows系統(tǒng)下的R安裝包惦辛，都是已經(jīng)編譯過的劳秋，你是無法看到源代碼的。
解壓之后你就看到所有函數(shù)的源代碼
我們要找的gdcVolcanoPlot的源代碼就在gdcDEGVisulization.R這個文件中。
我們照"抄"gdcVolcanoPlot這個函數(shù)俗批，接下來我們就可以用這個函數(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))
}
?
接下來我們來畫火山圖市怎，數(shù)據(jù)是從DEGAll.rda這個文件中來岁忘，具體如何生成這個文件和如何使用這個文件可以參考R的save，load函數(shù)和 .rda文件区匠。加載之后你就會有DEGall這個變量了干像，里面存放的是差異表達分析之后的結(jié)果。畫火山圖需要用到logFC驰弄，F(xiàn)DR麻汰。
load("DEGAll.rda")
#這里用到ggplot2這個包來畫圖
library(ggplot2)
gdcVolcanoPlot(DEGAll)
你就會得到下面這張火山圖，是不是很方便戚篙，不會寫函數(shù)一樣可以畫火山圖五鲫。
Reference：
2.R函數(shù)
DEGAll.rda文件的獲取方式請參考下面這篇文章
