r - 使用 apply、stringr、stringi 和 rbind 使函数运行得更快

Question

背景:我将提供此代码的应用背景和编程背景。希望两者都有帮助。我从事基因组学计算工作。是的——只是另一位冒充计算机科学家的生物学家。我正在编写一个脚本，该脚本将允许我按照人类基因组中的每个位置整合一堆数据集。这意味着数据帧超过 30 亿行 x 12 列。作为测试数据集，我正在使用酵母基因组构建分析管道，这将生成一个包含约 2500 万行和 12 列的数据帧。

问题:我当前的代码工作正常，但速度非常慢。例如，我在 45 分钟前启动了我的管道，其大约 1/3 穿过了酵母基因组。这意味着完成一份酵母样本可能需要 135 分钟，或者完成一份人体样本需要 270 小时……现在将其乘以我准备分析的 90 个人体样本，您有望看到我的问题。我需要加快速度。我将对此进行并行化，但即便如此，我认为它本身的代码也太笨重了。我需要帮助使我现有的功能变得更快。请不要告诉我我需要并行化它(这会得到否决票)。

示例数据:

chrom <- c("chr1", "chr1", "chr1", "chr1")
start <- c("0","1","2","6")
stop <- c("1","2","6","7")
sequence <- c("a", "t", "tcag", "a")
seqData <- data.frame(chrom, start, stop, sequence)

示例输出:

chrom_out <- c("chr1", "chr1", "chr1", "chr1", "chr1", "chr1", "chr1")
start_out <- c("0", "1", "2", "3", "4", "5", "6")
stop_out <- c("1", "2", "3", "4", "5", "6", "7")
sequence_out <- c("a", "t", "t", "c", "a", "g", "a")
out_seqdata <- data.frame(chrom_out, start_out, stop_out, sequence_out)

当前代码:

library(dplyr)
library(stringi)
library(stringr) 


wl = function(x){

  length<- stri_length(x["sequence"])
  if(length ==1){
    tmpseq<- x["sequence"]
    tmpstart <- as.numeric(x["start"])
    tmpstop <- as.numeric(x["stop"])
    tmpchrom <- x["chrom"]
    tmpdf <- data.frame(tmpseq, tmpstart, tmpstop, tmpchrom)
    colnames(tmpdf)<- c("tmpseq", "tmpstart", "tmpstop", "tmpchrom")
    print(tmpdf)
  }else{
    tmpseq<- strsplit(x["sequence"], "(?<=.{1})", perl = TRUE)
    tmpstart <- as.numeric(x["start"])+(1:length-1)
    tmpstop<- as.numeric(x["start"])+(1:length)
    tmpdf <- data.frame(tmpseq, tmpstart, tmpstop)
    tmpdf$tmpchrom <- x["chrom"]
    colnames(tmpdf)<- c("tmpseq", "tmpstart", "tmpstop", "tmpchrom")
    print(tmpdf)
  }
}

代码说明:我使用 apply 来迭代数据帧的每一行。数据框是坐标列表以及这些坐标的基因组序列。 Chrom = 染色体，start = 染色体上的起始位置，stop = 停止位置，sequence 为实际序列。数据当前采用压缩格式，以第三行数据为例。我想扩展这些数据，使每个基因组字母成为自己的行，然后适当调整坐标范围。函数 wl(代表从宽到长)执行此操作。它首先确定序列的字符串长度。如果长度等于 1，则将该行作为数据帧返回，无需进一步操作；否则，它将字符串分解为单独的字符，确定每个字符的坐标，并返回此数据帧。结果是一个数据帧列表，然后将其绑定(bind)在一起，生成示例输出数据。

我需要什么:我将对基因组进行分块，创建一个列表，从而允许我并行化该列表。这些 block 将产生一系列长度约为 2500 万行的数据帧。我也将并行化多个样本。并行化中的并行化……听起来像是使集群崩溃的好方法。我知道如何做到这一点(既编写此代码又使集群崩溃)。我需要帮助的是使实际功能更快。使用我当前的函数处理 2500 万行仍然需要很长时间。任何想法将不胜感激。请编辑我的函数或推荐一种新方法 - 欢迎所有想法。除了增加更多马力之外，我不知道还有更快的方法。

Answer 1

您可以矢量化所有操作:

# Generate vector of start positions
# Goes from 0 (minimal position in given data) to maximum base position in chromosome
foo <- 0:max(as.numeric(as.character(seqData$start)))
# Split sequence into a character vector
bar <- unlist(strsplit(as.character(seqData$sequence), ""))
# Generate final data frame
data.frame(start = foo, end = foo + 1, seq = bar)
#   start end seq
# 1     0   1   a
# 2     1   2   t
# 3     2   3   t
# 4     3   4   c
# 5     4   5   a
# 6     5   6   g
# 7     6   7   a

您可以使用此代码一次迭代一个染色体。

自定义函数和易于并行的 foreach 循环可能如下所示:

wl <- function(data, chr) {
    startPos <- 0:max(as.numeric(as.character(data$start)))
    nucs     <- unlist(strsplit(as.character(data$sequence), ""))
    data.frame(chr, start = startPos, end = startPos + 1, seq = nucs)
}
library(foreach)
# use dopar for parallel computations 
foreach(i = unique(seqData$chr), .combine = rbind) %do% {
    wl(subset(seqData, chrom == i), i)
}

PS:我永远不会使用基因组坐标作为字符向量。此外，创建 end 列只是浪费空间，因为您知道它的位置距离 start 为 1。