【R语言】从Function到lapply,再到future.a

写在前面
上一节介绍了For循环在科研上的小例子，本节来讲述一下自己编写函数的乐趣。

1.引入

结构

function.name <- function(arguments){
  computations_on_the_arguments
}

举例一：
传入两个参数，计算x的y次方加上1。

> yourSecondfun <- function(a,b){
  a^b+1
  }
>  yourSecondfun(5,3)
>  126

举例二：
计算乘方和乘积，同时返回。默认只返回最后一个参数，如果想返回多个参数，要使用return()函数。

doubleget <- function(x,y){
  a <- x^y
  b <- x*y
  return(c(a,b))
}
doubleget(2,3)

2.For循环和函数的相互转换

For循环根据TCGA_id第14和15位的数字来判断是tumor还是sample。

rm(list = ls())
load(file = "data/TCGA_ggplot.Rdata")
metadata <- data.frame(TCGA_id=rownames(exprSet)) #转换成数据框
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")

改写成函数，使其输入一个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.写一个作图函数

> 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

改写成函数的形式，使其输入一个基因名，就可以拿到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

基因和免疫浸润的相关性分析

rm(list = ls())
load(file = "data/expr_immu_BRCA.Rdata")
gene <- "KLF5"
y <- as.numeric(expr_data[,gene])
### 批量操作的具体实现过程：
### 1.设定容器,最终生成的数据放在什么地方？
correlation <- data.frame()
### 2.批量把数据导出到容器
for(i in 1:length(colnames(immu_data))){
  ## 1.指示
  print(i)
  ## 2.计算
  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")

将上述代码转换成lapply+function的模式来计算。
lapply 以及do.call 的介绍：https://mp.weixin.qq.com/s/VvbZOXlefGGEcekYWl0ZbQ
三步走战略。
第一步：写出单次处理的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)
}
### 测试函数功能
mycor("Activated B cell")

第二步：lapply批量作用于函数，返回list

lapplylist = lapply(colnames(immu_data),mycor)

第三步：do.call 转换list。lapply填充的是列表+函数；do.call填充的是函数+列表。

cor_data <- do.call(rbind,lapplylist)

通过画图来展示数据结果。

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. 单基因的批量相关性分析

加载数据

rm(list = ls())
load(file = "data/TCGA_BRCA_exprSet_plot.Rdata")
test <- exprSet[1:10,1:10]

A.使用for循环来完成。

##1.设定容器,最终生成的数据放在什么地方？
correlation <- data.frame()

##2.准备数据
data <- exprSet[,-c(1,2,3)]
test <- data[1:10,1:10]

##3.获取批量操作的范围,应该是个向量
genelist <- colnames(data)

##4.开始for循环,批量执行单次操作
gene <- "FOXA1"
genedata <- as.numeric(data[,gene])
for(i in 1:length(genelist)){
  ## 1.指示
  print(i)
  ## 2.计算
  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)
}
### 测试函数功能
mycor("GATA3")

### 第2步lapply批量作用于函数，返回list
system.time(lapplylist <- lapply(genelist,mycor))

### 第3步do.call 转换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的不同之处在于，前者可以调用电脑的多线程来工作，而后者只能用单线程。如果计算量不大的话，可能多线程工作还不如单线程省时。
加载R包。

install.packages("future.apply")
library(future.apply)
plan(multisession) #查看电脑总共有多少线程

三步走战略

### 第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)
}
### 测试函数功能
mycor("GATA3")

### 第2步lapply批量作用于函数，返回list
### system.time(lapplylist <- lapply(genelist,mycor))
system.time(lapplylist <- future_lapply(genelist,mycor))

### 第3步do.call 转换list
cor_data <- do.call(rbind,lapplylist)