神奇的多参数批量函数mapply

今天的帖子，是对在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语言学习路上的忆苦思甜