如何使用 R 在大型数据帧中加速 forloop grep

Question

我需要帮助。

我的脚本适用于许多数据帧，即使这需要几个小时（在集群上：> 100 GB 内存）。对于一些大型数据帧（> 300 万行），即使运行两天后，for 循环也不起作用。 因此，如果有办法加快 for 循环或用 R 中更多的加速器函数替换脚本，我需要帮助。

这是对我的脚本/数据的简短描述：

    snp1 <- c("R0100004", "R0100009", "R0100044", "R0100061", "R0100066","R0100067") # 3 million SNPs
    snp2 <- c("R0100039", "R0100152", "R0100066", "R0100067", "R0100068", "R0100082") # 3 million SNPs
    
    blocks <- c("R0100004|R0100009|R0100190|R0100015|R0100016|R0100017|R0100018|R0100021|R0100022|R0100024|R0100025",
                "R0100039|R0100038|R0100037|R0100036|R0100043|R0100044",                                   
                "R0100220|R0100052|R0100053|R0100054|R0100055|R0100057|R0100058|R0100059",                         
                "R0100061|R0100066|R0100067",                                                                      
                "R0100068|R0100069|R0100071|R0100072|R0100073|R0100074|R0133440|R0100076|R0100077|R0100078",         
                "R0100079|R0100081|R0100082") # 50000 haplotype block: each block contain > 2 SNPs. The SNP could be in 2 or more blocks.

# THE OBJECTIVE: For each SNP (snp1 and snp2) find his haplotype block    
    # This is my forloop
    I <- length(snp1) # 3000000
    res1 <- list()
    res2 <- list()
    for(j in 1:I){ 
        myres1 <- list(grep(snp1[j], blocks, value=T))
        myres2 <- list(grep(snp2[j], blocks, value=T))
        res1[j] <- myres1
        res2[j] <- myres2
    }

for 循环适用于中等数据帧，但对于具有大行的数据帧则需要几天时间。 如何替换或加快这个 for 循环？

提前致谢。

Answer 1

这会很快，但可能不会完全符合您的要求。

我根据示例数据假设了一些事情 - 例如每个snp1 或snp2 值最多只出现在一个块中，并且blocks 比snp1、snp2 短得多。

library(data.table)
y <- data.table(block = blocks)[,
  .(Rnum = unlist(strsplit(block, "\\|"))),
  .(block)]
res1 <- as.list(y[snp1, on = 'Rnum', block])
res2 <- as.list(y[snp2, on = 'Rnum', block])

如果我的假设有误，请告诉我。

但是，我强烈怀疑，如果您描述您实际尝试做的事情，将会有比这更好的解决方案。

Answer 2

我尝试了完全不同的方法，并故意将您的数据大小增加了 500。我提出的解决方案是对 500 倍大的数据集快 5 倍，如果数据大小增加，应该会越来越快。

library(tidyverse)

scale_factor <- 500 # to basically make your objects larger

# recreating your objects, but scaling them using the scaling_factor
snp1 <- c("R0100004", "R0100009", "R0100044", "R0100061", "R0100066","R0100067") %>% rep(scale_factor)
snp2 <- c("R0100039", "R0100152", "R0100066", "R0100067", "R0100068", "R0100082") %>% rep(scale_factor)

blocks <- c("R0100004|R0100009|R0100190|R0100015|R0100016|R0100017|R0100018|R0100021|R0100022|R0100024|R0100025",
            "R0100039|R0100038|R0100037|R0100036|R0100043|R0100044",                                   
            "R0100220|R0100052|R0100053|R0100054|R0100055|R0100057|R0100058|R0100059",                         
            "R0100061|R0100066|R0100067",                                                                      
            "R0100068|R0100069|R0100071|R0100072|R0100073|R0100074|R0133440|R0100076|R0100077|R0100078",         
            "R0100079|R0100081|R0100082")

# this is your original code
original <- function(snp1, snp2, blocks){
  # This is my forloop
  I <- length(snp1) # 3000000
  res1 <- list()
  res2 <- list()
  for(j in 1:I){ 
    myres1 <- list(grep(snp1[j], blocks, value=T))
    myres2 <- list(grep(snp2[j], blocks, value=T))
    res1[j] <- myres1
    res2[j] <- myres2
  }
  
}

tuned <- function(snp1, snp2, blocks){
  
  df1 <- data.frame(snp1)
  df2 <- data.frame(snp2)
  
  df_blocks <- data.frame(blocks) %>%
    mutate(blocks2 = blocks %>% str_split("\\|")) %>%
    unnest(cols = c(blocks2))
  
  res1 <- left_join(df1, df_blocks, by = c("snp1" = "blocks2")) %>%
    pull(blocks) %>%
    map(~.x)
  
  res2 <- left_join(df2, df_blocks, by = c("snp2" = "blocks2")) %>%
    pull(blocks) %>%
    map(~.x)

  
}

# this benchmarks both solutions
library(microbenchmark)
microbenchmark(original(snp1, snp2, blocks), tuned(snp1, snp2, blocks))
# Unit: milliseconds
# expr     min       lq     mean  median       uq      max neval
# original(snp1, snp2, blocks) 63.0723 65.65775 71.18956 69.3941 73.14755 102.5951   100
# tuned(snp1, snp2, blocks) 13.3522 14.39015 16.24824 15.5053 16.17280  33.6679   100