前言

之前学习 perl, 我最直观的感觉是,写到函数以后,涉及到函数的引用、面向对象编程时,% {} $ -> @_ $_ $$ $ \@ \% 这堆符号让人伤透了心,于是当我看见 python 没有这堆烦人的东西、并且可以支持画图、支持科学、统计计算以后,立即决定转投python。

相对python 而言,R 给我的最直观感觉是 for 循环什么的速度非常慢,比python 慢 5~10 倍吧,而且想自己写点能当library用的东西太难,一度有些想放弃。

然而确实是有些所谓的“毅种循环”吧,后来代码写多了,认识也提高了,很多之前觉得用R语言没法算的东西,发现其实是有相应的模块可以用的,只是这些模块跟perl 的那堆符号一样难以理解。又过了很长时间,自己写python 能写出经常复用的模块以后、懂了些map reduce以后,再回头看看这些R代码,才发现是这么回事(我也够笨的)。

今天的博客,主要介绍 R 语言中取代for循环的一系列东西。

1. 函数 apply, lapply, sapply的基本用法

参考自: http://stackoverflow.com/questions/3505701/r-grouping-functions-sapply-vs-lapply-vs-apply-vs-tapply-vs-by-vs-aggrega.

1.1 apply

1.1.1 基本

应用:对一个matrix 的每行、每列执行某种运算

M <- matrix(seq(1,160000), 400, 400)

# 行求和
M.rowSum <- apply(M, 1, sum)

# 列求和
M.colSum <- apply(M, 2, sum)

1.1.2 进阶

使用apply 可以使 R 的运算速度提高10倍,但仍然比python 慢几倍。若想进一步加速,可以考虑 RCpp。详见我之前的blog.http://blog.csdn.net/u012551177/article/details/25737303

1.2 lapply

1.2.1 基本

应用:对一个list 的每个元素,执行某种计算

x <- list(a = 1, b = 1:3, c = 10:100) 
lapply(x, FUN = length) 
## $a
## [1] 1
## 
## $b
## [1] 3
## 
## $c
## [1] 91
lapply(x, FUN = sum) 
## $a
## [1] 1
## 
## $b
## [1] 6
## 
## $c
## [1] 5005

1.2.2 进阶

个人认为,栈溢出给的例子远远不够。这里来个高难度的:

使用 iris 数据

head(iris)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa
unique(iris$Species)
## [1] setosa     versicolor virginica 
## Levels: setosa versicolor virginica

这里要求按照 Species: virginica setosa versicolor 顺序排列, Sepal.Length 从高到低排列

df_out <- lapply(c("virginica", "setosa", "versicolor"),function(i){
  tmp <- iris[iris$Species==i,]
  tmp[order(tmp$Sepal.Length, decreasing = T),]
})

分成了三个列表,每个列表叫做 “vireginica”, “setosa”, “versicolor”,然后多个列表逐行合并。

df_out_sort <- Reduce(rbind, df_out)

1.3 sapply

1.3.1 基本

x <- list(a = 1, b = 1:3, c = 10:100)

求 a b c 各自的个数

sapply(x, length)
##  a  b  c 
##  1  3 91

求 a b c 各自组的平均值

sapply(x, mean)
##  a  b  c 
##  1  2  55

1.3.2 进阶

strsplit 与 paste + sapply 相互结合,分割、合并字符串。 以分割路径为例:

join_str<-function(in_vec,join_str){
  out_str<-in_vec[1]
  if (length(in_vec) > 1){
    for (i in 2:length(in_vec))
    {
      out_str = paste(out_str,in_vec[i],sep=join_str  )  
    }
  }
  return(out_str)
}

get_args_infoV2<-function(args_in){
  split.arg<-strsplit(args_in,'/',perl=TRUE)
  prefix_strs<-split.arg[[1]][1:length(split.arg[[1]])]
  file_Info<-NULL
  file_Info$prefix<-join_str( prefix_strs[1:length(split.arg[[1]])-1],"/" )
  if ( is.na(file_Info$prefix) ){
    file_Info$prefix<-c("./")
  }
  file_Info$infile<-tail(prefix_strs,n=1)
  split.out<-strsplit(file_Info$infile,'\\.',perl=TRUE)
  prefix_out_strs<-split.out[[1]][1:length(split.out[[1]])]
  file_Info$outprefix<-join_str( prefix_out_strs[1:length(split.out[[1]])-1],"." )
  return( file_Info )
}

args_in <- c("/Users/hubq/Downloads/Project/github/Rmarkdown_setting")
get_args_infoV2(args_in)
## $prefix
## [1] "/Users/hubq/Downloads/Project/github"
## 
## $infile
## [1] "Rmarkdown_setting"
## 
## $outprefix
## [1] NA

lapply + sapply 相互结合算 pearson correlation 的 correlation 值 与 pvalue。

M_1 <- matrix(seq(1,16), 4, 4)
M_2 <- matrix(seq(16,1), 4, 4)
set.seed(1)
M_1 <- M_1 + matrix(rnorm(16), 4,4)
M_2 <- M_2 + matrix(rnorm(16), 4,4)

l_value <- lapply( c(1:dim(M_1)[1]), function(i){
  c <- cor.test(M_1[i,], M_2[i,])
  df<- data.frame("corr"=as.numeric(c$estimate), "pval"=as.numeric(c$p.value))
  df
})
sapply(l_value, "[[", "corr")
## [1] -0.9741776 -0.9971313 -0.9987768 -0.9697326
sapply(l_value, "[[", "pval")
## [1] 0.025822354 0.002868681 0.001223214 0.030267376

2. 多核运算

以 RNA-seq 的差异表达基因分析为例。如果是经典的 实验组-对照组 的实验设计,则对实验组、对照组分组完成后,找差异表达基因即可。

但如果不是呢?比如拿到了不同组织器官的RNA-seq数据,如果有N种,两两比较,就有 N(N-1)/2 次比较,而每次比较都要做一次 差异表达分析。

这种情况下,想加快分析,就需要使用更多的CPU。

一种方法是把一个任务分解成多个进程。这种方法如果样本少是不错,但如果N达到了20以上,则我们需要先执行200以上个程序,然后等结束后,执行一个程序整理结果。这样虽然快,但这个程序比较尴尬,差不多运行20分钟,想去干点别的也不好安排时间。如果在命令行中要求执行完200个自动执行下一个,又总觉得不太对劲。

能否一个程序,调动多个CPU并行分析,然后汇总结果?

这里提供一个例子:

args<-commandArgs(T)
if(length(args)<1)
{
  cat("Rscript DESeq.R merge.dexseq_clean.xls\n")
  q()
}

join_str<-function(in_vec,join_str){
  out_str<-in_vec[1]
  if (length(in_vec) > 1){
    for (i in 2:length(in_vec))
    {
      out_str = paste(out_str,in_vec[i],sep=join_str  )  
    }
  }
  return(out_str)
}

get_args_infoV1<-function(args_in){
  split.arg<-strsplit(args_in[1],'/',perl=TRUE)
  prefix_strs<-split.arg[[1]][1:length(split.arg[[1]])]
  file_Info<-NULL
  file_Info$prefix<-join_str( prefix_strs[1:length(split.arg[[1]])-1],"/" )
  if ( is.na(file_Info$prefix) ){
    file_Info$prefix<-c("./")
  }
  file_Info$infile<-tail(prefix_strs,n=1)
  split.out<-strsplit(file_Info$infile,'\\.',perl=TRUE)
  prefix_out_strs<-split.out[[1]][1:length(split.out[[1]])]
  file_Info$outprefix<-join_str( prefix_out_strs[1:length(split.out[[1]])-1],"." )
  return( file_Info )
}

DESeq_func<-function(comp,data_info,data,experimentDesign,file_Info){
  require(DESeq2)
  print(comp)
  stage1 = strsplit( comp ,"__vs__")[[1]][1]
  stage2 = strsplit( comp ,"__vs__")[[1]][2]
  print(stage1)
  print(stage2)
  all_idx<-rep(1:length(data_info$sample))
  idx_s1 = all_idx[data_info$stage==stage1]
  idx_s2 = all_idx[data_info$stage==stage2]
  
  countTable<-data[,c(idx_s1,idx_s2)]
  print(stage1)
  print(stage2)
  experimentDesign_less<-subset( experimentDesign, stage==stage1 | stage==stage2 )
  print( experimentDesign_less$stage  )
  experimentDesign_less$stage<-factor( experimentDesign_less$stage,level=unique(experimentDesign_less$stage) )
  
  print( experimentDesign_less  )
  dds = DESeqDataSetFromMatrix(countData=countTable,colData=experimentDesign_less,design=~stage)
  
  dds = estimateSizeFactors( dds )
  dds = estimateDispersions( dds )
  dds <- nbinomWaldTest(dds)
  
  res <- results( dds )
  resOrdered <- res[order(res$padj),]
  df_res<-as.data.frame(resOrdered)
  df_res$gene<-row.names(df_res)
  write.table( df_res,sep="\t",quote=FALSE,row.names=F, file=paste(file_Info$outprefix,"no_Mol.result",stage1,stage2,"xls",sep="."))  
}

library( "DESeq2" )
args<-c("CountRNA.table.xls")  # result of HTSeq
file_Info<-get_args_infoV1(args)
setwd(file_Info$prefix)
data<-read.table( file_Info$infile,header = TRUE )
rownames(data)<-data$Gene
data$Gene<-NULL
data_info=data.frame(
    sample = colnames(data),
    stage = paste(
                sapply( strsplit(colnames(data),"_"), "[[",1),
                sapply( strsplit(colnames(data),"_"), "[[",2),
                sapply( strsplit(colnames(data),"_"), "[[",3), sep="_" 
            )
    )
row.names(data_info)<-data_info$sample
data_info$stage <-factor( data_info$stage, level=unique(data_info$stage)  )

experimentDesign = data.frame(
    stage = as.character(paste(
            sapply( strsplit(colnames(data),"_"), "[[",1),
            sapply( strsplit(colnames(data),"_"), "[[",2),
            sapply( strsplit(colnames(data),"_"), "[[",3),
        sep="_")),
    row.names=colnames(data)
)
experimentDesign$stage<-factor( experimentDesign[,],level=unique(experimentDesign$stage) )

stages<-as.vector(unique(experimentDesign$stage))

数据准备完毕。下面生成两两比较的结果:

cll<-combn(stages,2)

使用 parallel 库并行运算。需要注意的一点是,前面 lapply 里面可以写全局变量,即lapply内部的function里面的变量可以外部定义,但是在多线程状态下就必须在内部定义。因此这里写了一长串参数。

此外上个月跟一大神学习了MPI编程,发现这个其实也可以用R MPI分析。还没想好怎么做,以后会写个专题专门研究。

library(parallel)

comp<-paste(t(cll)[,1],t(cll)[,2],sep="__vs__")
cl <- makeCluster(mc <- getOption("cl.cores", 8))
dfd<-parLapply(cl,comp,function(i,DESeq_func,data_info,data,experimentDesign,file_Info){
  print( i )
  x<-DESeq_func(i,data_info,data,experimentDesign,file_Info)
},DESeq_func,data_info,data,experimentDesign,file_Info
)
stopCluster(cl)



本文作者Boqiang Hu, 欢迎评论、交流。
转载请务必标注出处: R语言进阶之 apply lapply