R语言入门13：apply系列函数（apply、lapply、sapply、tapply）

本期目录：

• apply

• lapply

• sapply

• tapply

• 其他apply函数

• Reduce和do.call

• Reduce

• do.call

for循环是一个元素一个元素的操作，在R语言中这种做法是比较低效的，更好的做法是向量化操作，也就是同时对一整行/列进行操作，不用逐元素操作，这样可以大大加快运行速度。

apply函数家族就是这样的一组函数，专门实现向量化操作，可替代for循环。

先举个简单的例子，说明下什么是向量化。假如你有如下一个向量a，你想让其中的每个元素都加1，你不用把每个元素单独拎出来加1：

a <- c(1,2,3,NA)
a + 1 # 直接加1即可，是不是很方便？
## [1]  2  3  4 NA

再举个例子，下面这个数据框中有一些NA除此之外还有一些空白，或者空格。如何批量替换这些值？

tmp <- data.frame(a = c(1,1,3,4),
                  b = c("one","two","three","four"),
                  d = c(""," ",NA,90),
                  e = c(" ",NA, "",20)
                  )
tmp
##   a     b    d    e
## 1 1   one          
## 2 1   two      <NA>
## 3 3 three <NA>     
## 4 4  four   90   20

比如，让NA都变成999。

常规的做法是：检查每一个值，确认它是不是NA，如果是，就改成999，如果不是，就不改。

向量化的做法是：

tmp[is.na(tmp)] <- 999
# tmp[tmp == NA] <- 999 # 错误的做法
tmp
##   a     b   d   e
## 1 1   one        
## 2 1   two     999
## 3 3 three 999    
## 4 4  four  90  20

再比如，让空白的地方变成NA：

tmp[tmp == ""] <- NA
tmp
##   a     b    d    e
## 1 1   one <NA>     
## 2 1   two       999
## 3 3 three  999 <NA>
## 4 4  four   90   20

为什么还有一些空白？因为有的空白是真空白，有的则是空格！

tmp[tmp == " "] <- NA
tmp
##   a     b    d    e
## 1 1   one <NA> <NA>
## 2 1   two <NA>  999
## 3 3 three  999 <NA>
## 4 4  four   90   20

以上示例旨在告诉大家，有很多时候并不需要逐元素循环，向量化是更好的方式。

apply

对数据框（或矩阵）按行或者按列执行某个操作。

下面使用一个例子演示。示例数据是从TCGA官网下载的COAD的mrna的表达矩阵，一共有1000行，100列，每一行表示一个基因，每一列表示一个样本。

load(file = "datasets/coad_mran_df.rdata")

dim(coad_mrna_df)
## [1] 1000  100
class(coad_mrna_df)
## [1] "data.frame"
coad_mrna_df[1:4,1:3]
##        TCGA-5M-AAT6-01A-11R-A41B-07 TCGA-AA-3552-01A-01R-0821-07
## MT-CO2                     28026.23                     32915.04
## MT-CO3                     29725.85                     30837.60
## MT-ND4                     19509.82                     22026.42
## MT-CO1                     23193.16                     20924.84
##        TCGA-AA-3867-01A-01R-1022-07
## MT-CO2                     21030.00
## MT-CO3                     21997.99
## MT-ND4                     17171.58
## MT-CO1                     15485.43

如果要对表达矩阵进行log2转换，无需单独对每个元素进行log2，直接对整个数据框进行log2即可：

coad_mrna_df <- log2(coad_mrna_df + 1)

如果要计算每一个基因在所有样本中的平均表达量，也就是计算每一行的平均值，使用apply就非常简单：

# apply主要是3个参数
# 第1个是你的数据框
# 第2个是选择行或者列，1表示行，2表示列
# 第3个是要执行的操作，可以是R自带函数，也可以是自编函数
# 自带函数不用加括号，直接写名字即可
tmp <- apply(coad_mrna_df, 1, mean)
head(tmp)
##   MT-CO2   MT-CO3   MT-ND4   MT-CO1  MT-ATP6   MT-ND3 
## 14.59276 14.43845 14.01330 14.04316 13.57397 13.40406

如果使用for循环，就会显得很麻烦，运行时间也会长一点：

tmp <- vector("numeric", nrow(coad_mrna_df))
for(i in 1:nrow(coad_mrna_df)){
  tmp[i] <- mean(as.numeric(coad_mrna_df[i,]))
}
head(tmp)
## [1] 14.59276 14.43845 14.01330 14.04316 13.57397 13.40406

除了3个主要的参数，apply还有一个...参数，它表示：如果你要执行的操作中还有其他参数，可以直接往后写。比如mean()这个函数有一个na.rm参数，表示要不要在计算时去除缺失值，你可以直接把这个参数写在后面：

tmp <- apply(coad_mrna_df, 1, mean, na.rm = TRUE) # na.rm是mean的参数
head(tmp)
##   MT-CO2   MT-CO3   MT-ND4   MT-CO1  MT-ATP6   MT-ND3 
## 14.59276 14.43845 14.01330 14.04316 13.57397 13.40406

如果要计算每一列的平均值，第2个参数就写2即可：

# 1是行，2是列
tmp <- apply(coad_mrna_df, 2, mean, na.rm = TRUE)
head(tmp)
## TCGA-5M-AAT6-01A-11R-A41B-07 TCGA-AA-3552-01A-01R-0821-07 
##                     7.754459                     7.921157 
## TCGA-AA-3867-01A-01R-1022-07 TCGA-AD-6895-01A-11R-1928-07 
##                     8.131564                     8.198273 
## TCGA-AA-3560-01A-01R-0821-07 TCGA-CM-6676-01A-11R-1839-07 
##                     7.917137                     8.056527

上面的示例只是为了演示apply的用法，实际上在计算某一行/列的均值/加和时，R自带了几个函数，比如计算每一行的均值：

tmp <- rowMeans(coad_mrna_df)
head(tmp)
##   MT-CO2   MT-CO3   MT-ND4   MT-CO1  MT-ATP6   MT-ND3 
## 14.59276 14.43845 14.01330 14.04316 13.57397 13.40406

其他几个类似函数：

rowMeans()
rowSums()
colMeans()
colSums()

下面比较一下3种方法的运行时间：

system.time({  # 最慢
  tmp <- vector("numeric", nrow(coad_mrna_df))
  for(i in 1:nrow(coad_mrna_df)){
    tmp[i] <- mean(as.numeric(coad_mrna_df[i,]))
    }
})
##    user  system elapsed 
##    0.39    0.00    0.40
system.time(tmp <- apply(coad_mrna_df, 1, mean))
##    user  system elapsed 
##    0.01    0.00    0.00
system.time(tmp <- rowMeans(coad_mrna_df)) # 最快
##    user  system elapsed 
##       0       0       0

要执行的操作除了可以是R自带的函数外，还可以是自编函数。比如：筛选在所有样本中的表达量的加和大于800的基因：

# 对每一行执行1个操作
# 计算每一行的加和，并和800进行比较 
tmp <- apply(coad_mrna_df, 1, function(x){sum(x)>800})
head(tmp)
##  MT-CO2  MT-CO3  MT-ND4  MT-CO1 MT-ATP6  MT-ND3 
##    TRUE    TRUE    TRUE    TRUE    TRUE    TRUE
table(tmp)
## tmp
## FALSE  TRUE 
##   650   350
#coad_mrna_df[tmp,]

当然上面只是为了演示如何在apply中使用自编函数，实际使用时还是用rowSums更快更简单：

tmp <- rowSums(coad_mrna_df) > 800
head(tmp)
##  MT-CO2  MT-CO3  MT-ND4  MT-CO1 MT-ATP6  MT-ND3 
##    TRUE    TRUE    TRUE    TRUE    TRUE    TRUE
table(tmp)
## tmp
## FALSE  TRUE 
##   650   350

再举个例子，选择方差大于1的行（方差小说明这个基因在所有样本中表达量都很接近，这种基因没有意义）

tmp <- coad_mrna_df[apply(coad_mrna_df,1,function(x){var(x)>1}),]
dim(tmp)
## [1] 178 100

lapply

对list的每一个对象执行某个操作，或者对data.frame的每一列执行某个操作，输出结果是list。lapply的首字母就是list的首字母。

使用方法：

lapply(X, FUN, ...)

# x是你的数据框或者列表
# FUN是你要执行的操作
# ...和apply中的...一样

比如，选择方差大于1的列：

# ?lapply
# 和apply非常像，但是不用选择行或列，默认就是列
tmp <- lapply(coad_mrna_df, function(x){var(x)>1})

class(tmp)
## [1] "list"
length(tmp)
## [1] 100

# coad_mrna_df[tmp,]

计算每一列的中位数：

tmp <- lapply(coad_mrna_df, median)
class(tmp)
## [1] "list"
length(tmp)
## [1] 100

展开列表：

class(unlist(tmp))
## [1] "numeric"

查看列表中每个对象的长度：

# 创建一个列表
g <- "My First List" # 字符串
h <- c(25, 26, 18, 39) # 数值型向量
j <- matrix(1:10, nrow=5) # 矩阵
k <- c("one", "two", "three") # 字符型向量
l <- list("apple",1,TRUE) # 列表

mylist <- list(title=g, ages=h, j, k, l)

查看每个对象的长度：

lapply(mylist, length)
## $title
## [1] 1
## 
## $ages
## [1] 4
## 
## [[3]]
## [1] 10
## 
## [[4]]
## [1] 3
## 
## [[5]]
## [1] 3
unlist(lapply(mylist, length))
## title  ages                   
##     1     4    10     3     3

多个数据框的批量保存，lapply版本：

df1 <- data.frame(
  patientID = c("甲","乙","丙","丁"),
  age = c(23,43,45,34),
  gender = c("男","女","女","男")
)

df2 <- data.frame(
  patientID = c("甲","乙","戊","几","庚","丁"),
  hb = c(110,124,138,142,108,120),
  wbc = c(3.7,4.6,6.4,4.2,5.6,5.2)
)

df3 <- data.frame(
  patientID = c("丙","乙","几","庚","丁"),
  rbc = c(4.5,4.3,4.5,3.4,4.2),
  plt = c(180,250,360,120,220))

df4 <- data.frame(
  patientID = c("丙","乙","几","庚","丁","甲","戊"),
  a = rnorm(7, 20),
  b = rnorm(7,10)
)

df5 <- data.frame(
  patientID = c("丙","乙","甲","戊"),
  d = rnorm(4, 2),
  e = rnorm(4,1)
)

df6 <- data.frame(
  patientID = c("乙","几","庚","丁"),
  f = rnorm(4, 2),
  g = rnorm(4,1)
)

使用lapply的方式和for循环非常像。

先把这些数据框放到一个列表中：

dataframes <- list(df1,df2,df3,df4,df5,df6)

然后批量保存，和前面的for循环比较一下，是不是基本一样？

lapply(1:length(dataframes), function(x){ 
  write.csv(dataframes[[x]],
            file = paste0("datasets/csvs/","df",x,".csv"),
            quote = F,row.names = F)
})
## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL
## 
## [[6]]
## NULL

如果列表中的对象有名字，也可以像下面这样实现，还是和for循环基本一样：

dataframes <- list(df1,df2,df3,df4,df5,df6)  # 放到1个列表中
names(dataframes) <- c("df1","df2","df3","df4","df5","df6") # 添加名字
names(dataframes) # 查看名字
## [1] "df1" "df2" "df3" "df4" "df5" "df6"

lapply(names(dataframes), function(x){
    write.csv(dataframes[[x]],
            file = paste0("datasets/csvs/",x,".csv"),
            quote = F,row.names = F)
})
## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL
## 
## [[6]]
## NULL

多个数据框的批量读取：

allfiles <- list.files("datasets/csvs",full.names = T)
allfiles
## [1] "datasets/csvs/df1.csv" "datasets/csvs/df2.csv" "datasets/csvs/df3.csv"
## [4] "datasets/csvs/df4.csv" "datasets/csvs/df5.csv" "datasets/csvs/df6.csv"

# 1行代码解决，可以和前面的for循环对比下
dfs <- lapply(allfiles, read.csv)

dfs[[1]]
##   patientID age gender
## 1        甲  23     男
## 2        乙  43     女
## 3        丙  45     女
## 4        丁  34     男

如果你没有使用全名，需要自己构建文件路径+文件名，借助paste0即可：

allfiles <- list.files("datasets/csvs")
allfiles
## [1] "df1.csv" "df2.csv" "df3.csv" "df4.csv" "df5.csv" "df6.csv"

# 自己写个函数即可
dfs <- lapply(allfiles,
              function(x){read.csv(paste0("datasets/csvs/",x))})

dfs[[1]]
##   patientID age gender
## 1        甲  23     男
## 2        乙  43     女
## 3        丙  45     女
## 4        丁  34     男

此时的x就代指df1.csv、df2.csv这些名字。

sapply

lapply的简化版本，输出结果不是list。如果simplify=FALSE和USE.NAMES=FALSE，那么sapply函数就等于lapply函数了。不如lapply使用广泛。

选择方差大于1的列：

tmp <- sapply(coad_mrna_df, function(x){var(x)>1})

# coad_mrna_df[tmp,]

计算每一列的中位数：

tmp <- sapply(coad_mrna_df, median)
class(tmp)
## [1] "numeric"
length(tmp)
## [1] 100

head(tmp)
## TCGA-5M-AAT6-01A-11R-A41B-07 TCGA-AA-3552-01A-01R-0821-07 
##                     7.632902                     7.631332 
## TCGA-AA-3867-01A-01R-1022-07 TCGA-AD-6895-01A-11R-1928-07 
##                     7.882883                     8.042666 
## TCGA-AA-3560-01A-01R-0821-07 TCGA-CM-6676-01A-11R-1839-07 
##                     7.730625                     7.873826

tapply

分组操作。根据某一个条件进行分组，然后对每一个组进行某种操作，最后进行汇总。这种数据处理思想是非常出名的：split-apply-combine。

brca_clin <- read.csv("datasets/brca_clin.csv",header = T)
dim(brca_clin)
## [1] 20  9
brca_clin[,4:5]
##            sample_type initial_weight
## 1  Solid Tissue Normal            260
## 2  Solid Tissue Normal            220
## 3  Solid Tissue Normal            130
## 4  Solid Tissue Normal            260
## 5  Solid Tissue Normal            200
## 6  Solid Tissue Normal             60
## 7  Solid Tissue Normal            320
## 8  Solid Tissue Normal            310
## 9  Solid Tissue Normal            100
## 10 Solid Tissue Normal            250
## 11       Primary Tumor            130
## 12       Primary Tumor            110
## 13       Primary Tumor            470
## 14       Primary Tumor             90
## 15       Primary Tumor            200
## 16       Primary Tumor             70
## 17       Primary Tumor            130
## 18       Primary Tumor            770
## 19       Primary Tumor            200
## 20       Primary Tumor            250

分别计算normal组和tumor组的weight的平均值：

# 主要是3个参数
tapply(X = brca_clin$initial_weight, 
       INDEX = brca_clin$sample_type, #组别是分类变量，不能数值型
       FUN = mean)
##       Primary Tumor Solid Tissue Normal 
##                 242                 211

分别计算normal组和tumor组的age的中位数：

tapply(brca_clin$age_at_index, 
       brca_clin$sample_type, 
       median)
##       Primary Tumor Solid Tissue Normal 
##                55.0                59.5

还有几个类似的函数，比如：aggregate和by。

# 和tapply基本一样，但是第2个参数必须是list
# 并支持根据多个变量进行分组
aggregate(brca_clin$age_at_index, 
          list(brca_clin$sample_type), 
          median)
##               Group.1    x
## 1       Primary Tumor 55.0
## 2 Solid Tissue Normal 59.5

aggregate(brca_clin$age_at_index, 
          list(brca_clin$sample_type
               ,brca_clin$ajcc_pathologic_stage), 
          median)
##                Group.1    Group.2    x
## 1        Primary Tumor    Stage I 56.0
## 2  Solid Tissue Normal    Stage I 68.5
## 3        Primary Tumor   Stage IA 49.0
## 4  Solid Tissue Normal   Stage IA 63.0
## 5        Primary Tumor  Stage IIA 67.5
## 6  Solid Tissue Normal  Stage IIA 78.0
## 7        Primary Tumor  Stage IIB 63.0
## 8  Solid Tissue Normal  Stage IIB 54.0
## 9        Primary Tumor Stage IIIA 47.0
## 10 Solid Tissue Normal Stage IIIA 39.0
## 11       Primary Tumor Stage IIIC 36.0

by也是一样的用法：组别需要是因子型或者列表：

by(brca_clin$age_at_index, 
   list(brca_clin$sample_type), 
   median)
## : Primary Tumor
## [1] 55
## ------------------------------------------------------------ 
## : Solid Tissue Normal
## [1] 59.5

by(brca_clin$age_at_index, 
   list(brca_clin$sample_type
        ,brca_clin$ajcc_pathologic_stage), 
   median)
## : Primary Tumor
## : Stage I
## [1] 56
## ------------------------------------------------------------ 
## : Solid Tissue Normal
## : Stage I
## [1] 68.5
## ------------------------------------------------------------ 
## : Primary Tumor
## : Stage IA
## [1] 49
## ------------------------------------------------------------ 
## : Solid Tissue Normal
## : Stage IA
## [1] 63
## ------------------------------------------------------------ 
## : Primary Tumor
## : Stage IIA
## [1] 67.5
## ------------------------------------------------------------ 
## : Solid Tissue Normal
## : Stage IIA
## [1] 78
## ------------------------------------------------------------ 
## : Primary Tumor
## : Stage IIB
## [1] 63
## ------------------------------------------------------------ 
## : Solid Tissue Normal
## : Stage IIB
## [1] 54
## ------------------------------------------------------------ 
## : Primary Tumor
## : Stage IIIA
## [1] 47
## ------------------------------------------------------------ 
## : Solid Tissue Normal
## : Stage IIIA
## [1] 39
## ------------------------------------------------------------ 
## : Primary Tumor
## : Stage IIIC
## [1] 36
## ------------------------------------------------------------ 
## : Solid Tissue Normal
## : Stage IIIC
## [1] NA

组别是因子型也可以（实测字符型也可以），比如：

# 可以看到sample_type是字符型
str(brca_clin)
## 'data.frame':    20 obs. of  9 variables:
##  $ barcode               : chr  "TCGA-BH-A1FC-11A-32R-A13Q-07" "TCGA-AC-A2FM-11B-32R-A19W-07" "TCGA-BH-A0DO-11A-22R-A12D-07" "TCGA-E2-A1BC-11A-32R-A12P-07" ...
##  $ patient               : chr  "TCGA-BH-A1FC" "TCGA-AC-A2FM" "TCGA-BH-A0DO" "TCGA-E2-A1BC" ...
##  $ sample                : chr  "TCGA-BH-A1FC-11A" "TCGA-AC-A2FM-11B" "TCGA-BH-A0DO-11A" "TCGA-E2-A1BC-11A" ...
##  $ sample_type           : chr  "Solid Tissue Normal" "Solid Tissue Normal" "Solid Tissue Normal" "Solid Tissue Normal" ...
##  $ initial_weight        : int  260 220 130 260 200 60 320 310 100 250 ...
##  $ ajcc_pathologic_stage : chr  "Stage IIA" "Stage IIB" "Stage I" "Stage IA" ...
##  $ days_to_last_follow_up: int  NA NA 1644 501 660 3247 NA NA 1876 707 ...
##  $ gender                : chr  "female" "female" "female" "female" ...
##  $ age_at_index          : int  78 87 78 63 41 59 60 39 54 51 ...
class(brca_clin$sample_type)
## [1] "character"

by(brca_clin$age_at_index, 
   brca_clin$sample_type, # 字符型也可以
   median)
## brca_clin$sample_type: Primary Tumor
## [1] 55
## ------------------------------------------------------------ 
## brca_clin$sample_type: Solid Tissue Normal
## [1] 59.5

先把sample_type变成因子型也可以：

brca_clin$sample_type <- factor(brca_clin$sample_type)
class(brca_clin$sample_type) # 变成因子型了
## [1] "factor"

# 也OK
by(brca_clin$age_at_index, 
   brca_clin$sample_type, # 字符型也可以
   median)
## brca_clin$sample_type: Primary Tumor
## [1] 55
## ------------------------------------------------------------ 
## brca_clin$sample_type: Solid Tissue Normal
## [1] 59.5

其他apply函数

还有vapply、mapply、rapply、eapply，用的很少，不再介绍。

vapply类似于sapply，提供了FUN.VALUE参数，用来控制返回值的行名，这样可以让程序更清晰易懂。

Reduce和do.call

Reduce

对多个对象进行累积操作。

比如，累加：

Reduce("+", 1:100)
## [1] 5050

再比如，多个数据框的merge，merge函数只能对两个数据框进行合并，但是如果有多个数据框需要合并怎么办？有100个怎么办？

批量读取多个数据框：

# 6个数据框
allfiles <- list.files("datasets/csvs",full.names = T)
allfiles
## [1] "datasets/csvs/df1.csv" "datasets/csvs/df2.csv" "datasets/csvs/df3.csv"
## [4] "datasets/csvs/df4.csv" "datasets/csvs/df5.csv" "datasets/csvs/df6.csv"

# 1行代码解决
dfs <- lapply(allfiles, read.csv)

# 查看其中1个
dfs[[2]]
##   patientID  hb wbc
## 1        甲 110 3.7
## 2        乙 124 4.6
## 3        戊 138 6.4
## 4        几 142 4.2
## 5        庚 108 5.6
## 6        丁 120 5.2

6个数据框的merge：

Reduce(merge, dfs)
##   patientID age gender  hb wbc rbc plt      a        b        d        e
## 1        乙  43     女 124 4.6 4.3 250 19.664 10.51165 2.474508 1.372298
##          f         g
## 1 2.862749 -0.384265

如果想要使用merge里面的参数怎么办？自己写函数即可：

# 这个函数只能有两个参数
Reduce(function(x,y){merge(x,y, by = "patientID")},
       dfs)
##   patientID age gender  hb wbc rbc plt      a        b        d        e
## 1        乙  43     女 124 4.6 4.3 250 19.664 10.51165 2.474508 1.372298
##          f         g
## 1 2.862749 -0.384265

do.call

使用场景：你有很多个数据框，而且每个数据框的内容都一样，你想把这些数据框拼接到一起。

df1 <- data.frame(
  patientID = 1:4,
  aa = rnorm(4,10),
  bb = rnorm(4,16)
)
df2 <- data.frame(
  patientID = 5:8,
  aa = rnorm(4,10),
  bb = rnorm(4,16)
)
df3 <- data.frame(
  patientID = 9:12,
  aa = rnorm(4,10),
  bb = rnorm(4,16)
)
df4 <- data.frame(
  patientID = 13:16,
  aa = rnorm(4,10),
  bb = rnorm(4,16)
)

不断地重复写rbind？没有必要。

ll <- list(df1,df2,df3,df4)

do.call(rbind, ll)
##    patientID        aa       bb
## 1          1  9.574481 15.24356
## 2          2  9.933919 15.83192
## 3          3 10.675271 15.60532
## 4          4 11.130001 16.94735
## 5          5 10.068181 15.07117
## 6          6  9.832190 18.76410
## 7          7  8.944788 15.92174
## 8          8 10.282279 17.53555
## 9          9  9.580775 16.12769
## 10        10  9.956511 16.31920
## 11        11 11.776207 15.37159
## 12        12 11.313994 14.55692
## 13        13 10.306852 16.04596
## 14        14  9.194999 13.03253
## 15        15  8.295845 17.77535
## 16        16  9.482168 16.35076

其实这种场景下使用Reduce也可以，但是数据量比较大的话还是do.call更快。

Reduce(rbind, ll)
##    patientID        aa       bb
## 1          1  9.574481 15.24356
## 2          2  9.933919 15.83192
## 3          3 10.675271 15.60532
## 4          4 11.130001 16.94735
## 5          5 10.068181 15.07117
## 6          6  9.832190 18.76410
## 7          7  8.944788 15.92174
## 8          8 10.282279 17.53555
## 9          9  9.580775 16.12769
## 10        10  9.956511 16.31920
## 11        11 11.776207 15.37159
## 12        12 11.313994 14.55692
## 13        13 10.306852 16.04596
## 14        14  9.194999 13.03253
## 15        15  8.295845 17.77535
## 16        16  9.482168 16.35076