dplyr 中的快速字符串计数

Question

我需要对小标题中的列中保存的字符串中的模式（固定）实例进行计数。mutate(x = str_count(col, pattern))正是我想要的，但它的速度不够快，无法处理我必须评估的字符串量。

简化的小标题

test = tibble(
  id=c(1,2,3),
  seq=c("ATCG","ATTT","CGCG")
)

有效，但在微基准测试中速度太慢

test %>% mutate(CpG = str_count(seq, "CG"))

理论上在单个序列上更快，但在我的 tibble 列上不起作用

这只是给出第一行计数的单个值

test %>% mutate(CpG = sum(gregexpr("CG", seq, fixed=TRUE)[[1]] > 0))

我试图让 purrr::map 工作，但我失败了......这是我无法运行的尝试：

test %>% mutate(CpG = map(~sum(gregexpr("CG", seq, fixed=TRUE)[[1]] > 0)))
test %>% mutate(CpG = map(seq, ~sum(gregexpr("CG", .x, fixed=TRUE)[[1]] > 0)))

编辑：时间基准

看起来stringi::stri_count_fixed是测试集上的最佳选择。

microbenchmark(
  mutate(test, CpG = stringr::str_count(seq, "CG")),
  mutate(test, CpG = purrr::map_int(seq, ~sum(gregexpr("CG", .x, fixed=TRUE)[[1]] > 0))),
  mutate(test, CpG = Biostrings::vcountPattern("CG", DNAStringSet(seq))),
  mutate(test, CpG = lengths(regmatches(seq, gregexpr("CG", seq)))),
  mutate(test, CpG = stringi::stri_count_fixed(seq,"CG"))
)  
Unit: microseconds
                                                                                          expr     min       lq     mean   median       uq      max neval
                                             mutate(test, CpG = stringr::str_count(seq, "CG")) 125.502 133.9995 159.0844 147.9045 164.2925  441.242   100
 mutate(test, CpG = purrr::map_int(seq, ~sum(gregexpr("CG", .x,      fixed = TRUE)[[1]] > 0))) 167.210 183.4695 215.7746 202.8890 230.5275  450.177   100
                        mutate(test, CpG = Biostrings::vcountPattern("CG", DNAStringSet(seq))) 789.889 827.3810 955.8919 887.3400 989.7415 1845.212   100
                             mutate(test, CpG = lengths(regmatches(seq, gregexpr("CG", seq)))) 150.315 159.4750 178.9039 169.1840 185.1345  316.479   100
                                      mutate(test, CpG = stringi::stri_count_fixed(seq, "CG")) 112.015 120.6615 131.4281 125.0470 135.4330  213.184   100

而且，以下是我的一部分真实数据（约 20,000 个 1500 nt 序列的实例）在远程服务器上运行的结果，该服务器的可用 RAM 和内核约为初始运行的 3 倍。我仍然对 的stringi实现印象深刻，它以某种方式在 20 毫秒左右的时间里读完了所有内容。

Unit: milliseconds
                                                                                          expr        min         lq       mean     median        uq        max neval
                                             mutate(test, CpG = stringr::str_count(seq, "CG"))  402.55513  408.04101  419.82213  414.31085  425.4113  481.08128   100
 mutate(test, CpG = purrr::map_int(seq, ~sum(gregexpr("CG", .x,      fixed = TRUE)[[1]] > 0)))  429.84148  436.02402  474.08864  438.84198  447.3058 1054.48253   100
                        mutate(test, CpG = Biostrings::vcountPattern("CG", DNAStringSet(seq)))  301.26062  309.76674  310.97071  310.08229  310.8837  336.12834   100
                             mutate(test, CpG = lengths(regmatches(seq, gregexpr("CG", seq)))) 1981.78309 1990.84206 2050.41928 1999.07389 2020.3675 2980.62566   100
                                      mutate(test, CpG = stringi::stri_count_fixed(seq, "CG"))   23.33313   23.55215   23.90881   23.66268   23.8916   27.40499   100

因为我可以在本机 DNAStringSet 上运行它，因为我的数据以 GRanges 对象开始，所以我也对序列向量进行了微基准测试，并得到了 240-250 毫秒范围内的范围。

我有另一个用例，其中包含 1,000-100,000 nt 长的序列和更长的搜索模式，并计划在我回到该项目时更新这些结果。

Answer 1

您可以使用Biostrings::vcountPattern

library(Biostrings)
test %>% mutate(CpG = vcountPattern("CG", DNAStringSet(seq)))
#  # A tibble: 3 x 3
#     id seq     CpG
#  <dbl> <chr> <int>
#1     1 ATCG      1
#2     2 ATTT      0
#3     3 CGCG      2

Biostrings::vcountPatternstr_count比较长 (DNA) 字符串向量的解决方案快得多；这是一些microbenchmark结果

# Sample data
df <- tibble(
    id = 1:100,
    seq = replicate(
        100,
        paste0(sample(c("A", "C", "G", "T"), 100000, replace = T), collapse = "")))


library(microbenchmark)
res <- microbenchmark(
    str_count = df %>% mutate(CpG = str_count(seq, "CG")),
    vmatchPattern = df %>% mutate(CpG = vcountPattern("CG", DNAStringSet(seq))))
#Unit: milliseconds
#          expr       min       lq     mean   median       uq      max neval
#     str_count 156.37642 162.1629 167.2005 164.4958 169.2456 251.3426   100
# vmatchPattern  95.68998 102.4176 108.2428 105.2501 108.5511 322.3143   100

library(ggplot2)
autoplot(res)