寫在前面
上一節(jié)介紹了For循環(huán)在科研上的小例子轨香，本節(jié)來講述一下自己編寫函數(shù)的樂趣。
1.引入
結(jié)構(gòu)
function.name <- function(arguments){
computations_on_the_arguments
}
舉例一：
傳入兩個(gè)參數(shù)幼东，計(jì)算x的y次方加上1臂容。
> yourSecondfun <- function(a,b){
a^b+1
}
> yourSecondfun(5,3)
> 126
舉例二：
計(jì)算乘方和乘積，同時(shí)返回根蟹。默認(rèn)只返回最后一個(gè)參數(shù)脓杉，如果想返回多個(gè)參數(shù)，要使用return()函數(shù)简逮。
doubleget <- function(x,y){
a <- x^y
b <- x*y
return(c(a,b))
}
doubleget(2,3)
2.For循環(huán)和函數(shù)的相互轉(zhuǎn)換
For循環(huán)根據(jù)TCGA_id第14和15位的數(shù)字來判斷是tumor還是sample球散。
rm(list = ls())
load(file = "data/TCGA_ggplot.Rdata")
metadata <- data.frame(TCGA_id=rownames(exprSet)) #轉(zhuǎn)換成數(shù)據(jù)框
for (i in 1:length(metadata[,1])) {
num <- as.numeric(substring(metadata[i,1],14,15))
if (num %in% seq(1,9)) {metadata[i,2] = "Tumor"}
if (num %in% seq(10,29)) {metadata[i,2] = "Normal"}
}
colnames(metadata) <- c("TCGA_id","sample")
改寫成函數(shù)，使其輸入一個(gè)TCGA_id散庶，就能判斷ID屬于是tumor還是sample蕉堰。
tell_me_sample <- function(TCGAID){
num <- as.numeric(substring(TCGAID,14,15))
if (num %in% seq(1,9)) {sample = "Tumor"}
if (num %in% seq(10,29)) {sample = "Normal"}
return(sample)
}
tell_me_sample("TCGA-AC-A3OD-01B-06R-A22O-07")
tell_me_sample("TCGA-BH-A0B5-11A-23R-A12P-07")
tell_me_sample(metadata[1,1])
tell_me_sample(metadata[3,1])
3.寫一個(gè)作圖函數(shù)
> rm(list = ls())
> load(file = "data/TCGA_ggplot.Rdata")
> my_comparisons <- list(
c("LumA", "Normal"),
c("LumB", "Normal"),
c("Her2", "Normal"),
c("Basal", "Normal")
)
> library(ggpubr)
> ggboxplot(
exprSet, x = "subgroup", y = "ESR1",
color = "subgroup", palette = c("#00AFBB", "#E7B800", "#FC4E07","#A687D8", "#89E4A4"),
add = "jitter"
)+
stat_compare_means(comparisons = my_comparisons, method = "t.test")
畫圖代碼來源:https://mp.weixin.qq.com/s/AE3Rt3IoWAkhvPYtEiwWQg
改寫成函數(shù)的形式，使其輸入一個(gè)基因名悲龟，就可以拿到boxplot屋讶。
> subgroup_plot <- function(gene){
ggboxplot(
exprSet,x='subgroup',y=gene,
color = 'subgroup',palette = c("#00AFBB","#E7B800","#FC4E07","#A687D8","#89E4A4"),
add = "jitter"
)+
stat_compare_means(comparisons = my_comparisons,method = "t.test")
}
> subgroup_plot("SAMD11")
> subgroup_plot("BRCA1")
4.lapply
基因和免疫浸潤的相關(guān)性分析
rm(list = ls())
load(file = "data/expr_immu_BRCA.Rdata")
gene <- "KLF5"
y <- as.numeric(expr_data[,gene])
### 批量操作的具體實(shí)現(xiàn)過程：
### 1.設(shè)定容器,最終生成的數(shù)據(jù)放在什么地方？
correlation <- data.frame()
### 2.批量把數(shù)據(jù)導(dǎo)出到容器
for(i in 1:length(colnames(immu_data))){
## 1.指示
print(i)
## 2.計(jì)算
dd = cor.test(as.numeric(immu_data[,i]),y,method="spearman")
## 3.填充
correlation[i,1] = colnames(immu_data)[i]
correlation[i,2] = dd$estimate
correlation[i,3] = dd$p.value
}
### 修改名稱
colnames(correlation) <- c("cell","cor","p.value")
將上述代碼轉(zhuǎn)換成lapply+function的模式來計(jì)算躲舌。
lapply 以及do.call 的介紹：https://mp.weixin.qq.com/s/VvbZOXlefGGEcekYWl0ZbQ
三步走戰(zhàn)略丑婿。
第一步：寫出單次處理的function
mycor = function(x){
dd = cor.test(as.numeric(immu_data[, x]),y,method ="spearman",exact=FALSE)
data.frame(cell=x,cor=dd$estimate,p.value=dd$p.value)
}
### 測試函數(shù)功能
mycor("Activated B cell")
第二步：lapply批量作用于函數(shù)，返回list
lapplylist = lapply(colnames(immu_data),mycor)
第三步：do.call 轉(zhuǎn)換list没卸。lapply填充的是列表+函數(shù)羹奉；do.call填充的是函數(shù)+列表。
cor_data <- do.call(rbind,lapplylist)
通過畫圖來展示數(shù)據(jù)結(jié)果约计。
library(dplyr)
library(ggplot2)
cor_data %>%
filter(p.value <0.05) %>%
ggplot(aes(cor,forcats::fct_reorder(cell,cor)))+
geom_segment(aes(xend=0,yend=cell))+
geom_point(aes(col=p.value,size=abs(cor)))+
scale_colour_gradientn(colours=c("#7fc97f","#984ea3"))+
#scale_color_viridis_c(begin = 0.5, end = 1)+
scale_size_continuous(range =c(2,8))+
theme_bw()+
ylab(NULL)
如果熟練了诀拭，三步走可以變成一步。
cor_data <- do.call(rbind,lapply(colnames(immu_data),function(x){
dd = cor.test(as.numeric(immu_data[,x]),y,method ="spearman",exact=FALSE)
data.frame(cell=x,cor=dd$estimate,p.value=dd$p.value)
}))
5. 單基因的批量相關(guān)性分析
加載數(shù)據(jù)
rm(list = ls())
load(file = "data/TCGA_BRCA_exprSet_plot.Rdata")
test <- exprSet[1:10,1:10]
A.使用for循環(huán)來完成煤蚌。
##1.設(shè)定容器,最終生成的數(shù)據(jù)放在什么地方耕挨？
correlation <- data.frame()
##2.準(zhǔn)備數(shù)據(jù)
data <- exprSet[,-c(1,2,3)]
test <- data[1:10,1:10]
##3.獲取批量操作的范圍,應(yīng)該是個(gè)向量
genelist <- colnames(data)
##4.開始for循環(huán),批量執(zhí)行單次操作
gene <- "FOXA1"
genedata <- as.numeric(data[,gene])
for(i in 1:length(genelist)){
## 1.指示
print(i)
## 2.計(jì)算
dd = cor.test(genedata,as.numeric(data[,i]),method="spearman")
## 3.填充
correlation[i,1] = gene
correlation[i,2] = genelist[i]
correlation[i,3] = dd$estimate
correlation[i,4] = dd$p.value
}
colnames(correlation) <- c("gene1","gene2","cor","p.value")
B.使用function和lapply來完成细卧。
### 第1,寫出單次處理的function
mycor = function(i){
dd = cor.test(genedata,as.numeric(data[,i]),method="spearman")
data.frame(gene1=gene,gene2=i,cor=dd$estimate,p.value=dd$p.value)
}
### 測試函數(shù)功能
mycor("GATA3")
### 第2步lapply批量作用于函數(shù)，返回list
system.time(lapplylist <- lapply(genelist,mycor))
### 第3步do.call 轉(zhuǎn)換list
cor_data <- do.call(rbind,lapplylist)
### 熟悉了可以寫成1步
cor_data <- do.call(rbind,lapply(genelist,function(i){
dd = cor.test(genedata,as.numeric(data[,i]),method="spearman")
data.frame(gene1=gene,gene2=i,cor=dd$estimate,p.value=dd$p.value)
}))
C.使用future_lapply + funcxtion
future_lapply和lapply的不同之處在于筒占，前者可以調(diào)用電腦的多線程來工作贪庙，而后者只能用單線程。如果計(jì)算量不大的話翰苫，可能多線程工作還不如單線程省時(shí)止邮。
加載R包。
install.packages("future.apply")
library(future.apply)
plan(multisession) #查看電腦總共有多少線程
三步走戰(zhàn)略
### 第1,寫出單次處理的function
mycor = function(i){
dd = cor.test(genedata,as.numeric(data[,i]),method="spearman")
data.frame(gene1=gene,gene2=i,cor=dd$estimate,p.value=dd$p.value)
}
### 測試函數(shù)功能
mycor("GATA3")
### 第2步lapply批量作用于函數(shù)奏窑，返回list
### system.time(lapplylist <- lapply(genelist,mycor))
system.time(lapplylist <- future_lapply(genelist,mycor))
### 第3步do.call 轉(zhuǎn)換list
cor_data <- do.call(rbind,lapplylist)
