今天的帖子，是对在R语言里面批量操作的总结：
事情的起源来自于临床师弟的需求
## 读入数据
test <- data.table::fread("test.csv",data.table = F,header = T)
因为excel的使用习惯，很多临床医生收集标本的时候，喜欢合并把同一个患者的信息，只写一次。但是，如果想要更好地分析，我们需要的是清洁数据：
行是观测，列是变量，也就是patient那一列不能有空，需要自动填充。
这个需求，我不知道在excel如何实现，也不知道在R如何方便地实现。但是我有底线：
一定能用R实现，实在不行就编程啊。
获取第一列的数据
获取非缺失部分的位置
计算出两个位置之间的重复次数
接下来只要把1重复1次，2重复4次，6重复3次，9重复5次，14重复6次就可以了。这个可以批量运行
得到了位置，填充也就十分简单
为了代码复用，我把让他们写成了一个函数，实际上，只要能清晰地定义每一次的操作，就可以方便地写出for循环，写出function。
fill_blank <- function(column){
index_T <- which(column != "")
index1 <- c(index_T[-1],length(column)+1)-index_T
index2 <- do.call(c,sapply(1:length(index_T), function(i){rep(index_T[i],index1[i])}))
column[index2]
}
给原数据增加新的一列
test$fill <- fill_blank(test$V1)
这是师弟在组会上问我的时候，我给出的方法。但是经过那一次技能大比拼之后，我知道了更多批量的方法。
R语言学习路上的忆苦思甜
其中一个就是mapply，他是apply家族的一员，批量的同时，接受多个参数。比如，我想把1重复1次，2重复2次，3重复3次，一次类推。mapply可以轻松实现
mapply(rep,1:10,1:10)
又比如，每次把不同字母重复不一样的次数
mapply(rep,LETTERS[1:4],c(2,6,5,3))
mapply返回的是list，可以通过unlist去掉
unlist(mapply(rep,LETTERS[1:4],c(2,6,5,3)))
我写的那个函数中有很长的一步就是处理这个事情，现在可以简化了
fill_blank2 <- function(column){
index_T <- which(column != "")
index1 <- c(index_T[-1],length(column)+1)-index_T
index2 <- unlist(mapply(rep,index_T,index1))
column[index2]
}
给原来的数据增加一列
test$fill2 <- fill_blank2(test$V1)
没有任何问题，此时我开始膨胀了，我隐隐感觉我能解决掉1年前熊给我提出的需求。
未能解决的来自大神熊的需求
有两个关于snp的数据，第一个里面有三列，第一列是染色体号，第列是具体位置，第三列是一个value值
数据名称叫big，是个数据框，15957行
第二个数据是，也是三列，表示的是染色体上按照100万分成的无数间隔（这种间隔成为bin）
这个数据的名称是window，很大
现在的需求是：
统计出染色体上每个间隔中，总共有多少个snp，以及所含有snp的平均value。
熊说了，在linux上实现起来相当容易，R语言里面处理起来速度太慢，当时我因为要面子，硬着头皮怒答一波：
然而并没有解决问题，速度还是太慢了。
近期，体内有种不知名的力量在涌动，所以，准备再来一次。
首先，我准备一个bin一个bin的统计，测试前1000个bin的时间
system.time(for(i in 1:1000){
## 从第一个bin开始计算
print(i)
## 限定snp的范围，在当前的染色体
big_f = big[big$chr==window$chr[i],]
rownames(big_f) =big_f$snp
## 看哪些snp在bin的范围中
index = intersect(big_f$snp,seq(window$start[i],window$end[i]))
## 统计个数
count = length(index)
## 计算平均值
mean = mean(big_f[as.character(index),"value"],na.rm = T)
})
大概27s，如果统计完所有的bin，大概是11个小时，显然是不符合要求的。
我喜欢用for循环的原因是，我可以通过限制某一行不运行来确定限速环节，最终发现是intersect那一行是限速行。
我把求交集的策略换成了比较大小
system.time(for(i in 1:1000){
print(i)
big_f = big[big$chr==window$chr[i],]
#rownames(big_f) =big_f$snp
index = big_f$snp < window$end[i] & big_f$snp > window$start[i]
count = sum(index)
mean = mean(big_f$value[index],na.rm = T)
})
这时候大概需要0.22s，计算一下完成所有的操作时间是5分钟
我问了一下当事人，大概多长时间能接受，他说的是1分钟以内。
于是，我就把for循环改成了函数
此时我知道了mapply的存在，所以就打算先把两个数据都按照染色体切割分类，然后写出一个函数，同时接受两个参数，最终再用futurue来提速。
函数是这样的
count_mean <- function(window_f,big_f){
do.call(rbind,lapply(1:nrow(window_f),function(i){
index = big_f$snp < window_f$end[i] & big_f$snp > window_f$start[i]
count = sum(index)
mean = mean(big_f$value[index],na.rm = T)
return(c(count=count,mean=mean))
}))
}
测试一下是否有效
test <- count_mean(split(window,window$chr)[[2]],split(big,big$chr)[[2]])
批量运行
library(future.apply)
plan(multiprocess)
system.time(test <- future_mapply(count_mean,split(window,window$chr),split(big,big$chr)))
最终消耗11秒，圆满完成需求。
转换数据的过程，时间可以忽略不计
results <- cbind(window,do.call(rbind,test))
思维总结
现在我的脑子里面对于批量运算有了小方案，不再惧怕批量运行了，我在R语言里面经常干的事情就是，批量把结果全部给我返回，然后我来挑最好的。
- 写出for循环，调整方案，测试速度
- apply+function 批量运行
- 如果速度慢，用future来并行化处理
关于批量，我还写过以下帖子，没事的时候可以看看。
8秒完成2万个基因的生存分析，人人都可以！
神奇的lapply
Reduce函数实现多个数据框批量合并
R语言学习路上的忆苦思甜
网友评论