a <- c(1,2,3,NA)
a + 1 # 直接加1即可，是不是很方便？
## [1] 2 3 4 NA11 apply系列
for循环是一个元素一个元素的操作，在R语言中这种做法是比较低效的，更好的做法是向量化操作，也就是同时对一整行/列进行操作，不用逐元素操作，这样可以大大加快运行速度。
apply函数家族就是这样的一组函数，专门实现向量化操作，可替代for循环。
先举个简单的例子，说明下什么是向量化。假如你有如下一个向量a，你想让其中的每个元素都加1，你不用把每个元素单独拎出来加1：
再举个例子，下面这个数据框中有一些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以上示例旨在告诉大家，有很多时候并不需要逐元素循环，向量化是更好的方式。
11.1 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.38 0.03 0.40
system.time(tmp <- apply(coad_mrna_df, 1, mean))
## user system elapsed
## 0.00 0.00 0.02
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 10011.2 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这些名字。
11.3 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.87382611.4 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.0by也是一样的用法：组别需要是因子型或者列表：
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.511.5 其他apply函数
还有vapply、mapply、rapply、eapply，用的很少，不再介绍。
vapply类似于sapply，提供了FUN.VALUE参数，用来控制返回值的行名，这样可以让程序更清晰易懂。
11.6 Reduce和do.call
11.6.1 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.26个数据框的merge：
Reduce(merge, dfs)
## patientID age gender hb wbc rbc plt a b d e
## 1 乙 43 女 124 4.6 4.3 250 19.54459 11.54804 0.6582348 0.4118437
## f g
## 1 1.465577 0.5200834如果想要使用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.54459 11.54804 0.6582348 0.4118437
## f g
## 1 1.465577 0.520083411.6.2 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()：
rbind(df1,df2,df3,df4)
## patientID aa bb
## 1 1 9.279656 16.08581
## 2 2 10.413684 14.81278
## 3 3 8.740267 16.61363
## 4 4 12.043236 15.67159
## 5 5 9.530871 16.08338
## 6 6 9.559042 15.21700
## 7 7 10.589750 15.29838
## 8 8 10.405087 18.42527
## 9 9 9.551934 15.77828
## 10 10 9.499350 16.40950
## 11 11 9.321850 16.47779
## 12 12 9.810859 17.11033
## 13 13 9.815785 15.54795
## 14 14 10.639532 12.96574
## 15 15 9.112032 16.42428
## 16 16 10.102449 17.18567但是通常我们在进行数据分析时，经常会使用循环或者向量化运行，这样得到的结果是列表，数据框可能是列表中的元素，此时直接使用rbind()是不行的：
# 数据框在列表中
ll <- list(df1,df2,df3,df4)
# 再用rbind是不能有同样效果的
#rbind(ll)
#lapply(ll, rbind)不断地重复写rbind？没有必要。
ll <- list(df1,df2,df3,df4)
do.call(rbind, ll)
## patientID aa bb
## 1 1 9.279656 16.08581
## 2 2 10.413684 14.81278
## 3 3 8.740267 16.61363
## 4 4 12.043236 15.67159
## 5 5 9.530871 16.08338
## 6 6 9.559042 15.21700
## 7 7 10.589750 15.29838
## 8 8 10.405087 18.42527
## 9 9 9.551934 15.77828
## 10 10 9.499350 16.40950
## 11 11 9.321850 16.47779
## 12 12 9.810859 17.11033
## 13 13 9.815785 15.54795
## 14 14 10.639532 12.96574
## 15 15 9.112032 16.42428
## 16 16 10.102449 17.18567其实这种场景下使用Reduce也可以，但是数据量比较大的话还是do.call更快。
Reduce(rbind, ll)
## patientID aa bb
## 1 1 9.279656 16.08581
## 2 2 10.413684 14.81278
## 3 3 8.740267 16.61363
## 4 4 12.043236 15.67159
## 5 5 9.530871 16.08338
## 6 6 9.559042 15.21700
## 7 7 10.589750 15.29838
## 8 8 10.405087 18.42527
## 9 9 9.551934 15.77828
## 10 10 9.499350 16.40950
## 11 11 9.321850 16.47779
## 12 12 9.810859 17.11033
## 13 13 9.815785 15.54795
## 14 14 10.639532 12.96574
## 15 15 9.112032 16.42428
## 16 16 10.102449 17.18567