Быстрый подсчет строк в 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"))

Теоретически быстрее на отдельных последовательностях, но не работает на моем столбце тиблов

Это просто дает единственное значение числа первой строки

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

Я пытался заставить работать purrr :: map, но у меня ничего не получается ... вот мои попытки, которые победили ' t run:

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 нт), запущенных на удаленном сервере с ~ 3X доступной оперативной памятью и ядрами моего начального запуска. Я по-прежнему впечатлен реализацией 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 мс.

У меня есть еще один вариант использования с 1000-100000 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::vcountPattern значительно быстрее, чем решение str_count для вектора более длинных (ДНК) строк; Вот некоторые 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)

Answer 2

Попробуйте использовать stri_count_fixed из библиотеки stringi, которая, как известно, работает быстрее, чем ее stringr аналоги.

test %>% mutate(CpG = stringi::stri_count_fixed(seq, 'CG'))

#     id seq     CpG
#  <dbl> <chr> <int>
#1     1 ATCG      1
#2     2 ATTT      0
#3     3 CGCG      2

---

Кроме того, вы можете использовать map_int, который будет возвращать целочисленный вектор по сравнению с map, который возвращает список для использования с gregexpr.

library(dplyr)

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

Или другой вариант - regmatches с gregexpr.

test %>% mutate(CpG = lengths(regmatches(seq, gregexpr("CG", seq))))

Answer 3

Используя набор тестовых данных @MauritsEvers df, вы можете использовать oligoNucleotideFrequency в Biostrings.

library(Biostrings)
head(oligonucleotideFrequency(DNAStringSet(df$seq),2))
       AA   AC   AG   AT   CA   CC   CG   CT   GA   GC   GG   GT   TA   TC   TG
[1,] 6402 6147 6351 6190 6129 6081 6180 6232 6311 6222 6323 6341 6249 6172 6342
[2,] 6305 6320 6295 6217 6379 6253 6139 6291 6227 6279 6239 6163 6225 6210 6235
[3,] 6284 6352 6305 6133 6211 6293 6356 6231 6374 6286 6280 6202 6205 6160 6201
[4,] 6217 6253 6265 6300 6189 6296 6302 6175 6379 6275 6207 6195 6251 6137 6282
[5,] 6346 6315 6217 6291 6289 6060 6254 6223 6269 6203 6313 6259 6264 6248 6260
[6,] 6066 6216 6282 6286 6376 6339 6223 6156 6216 6235 6360 6241 6193 6304 6186
       TT
[1,] 6327
[2,] 6222
[3,] 6126
[4,] 6276
[5,] 6188
[6,] 6320

Мы просто получаем столбец "CG":

library(microbenchmark)
res <- microbenchmark(
    str_count = df %>% mutate(CpG = str_count(seq, "CG")),
    vmatchPattern = df %>% mutate(CpG = vcountPattern("CG", DNAStringSet(seq))),
    oligo=df %>% mutate(CpG = oligonucleotideFrequency(DNAStringSet(seq),2)[,"CG"])
)

> res
Unit: milliseconds
          expr       min        lq      mean    median        uq      max neval
     str_count 129.73744 135.05308 140.01877 139.11458 144.47085 158.9793   100
 vmatchPattern  82.30784  88.48554  92.13986  93.01389  95.53329 107.3482   100
         oligo  57.25595  62.11515  67.63134  64.50559  66.99784 349.3443   100
 cld
   c
  b 
 a  

Вы уже знаете, что оснований или динуклеотидов, чтобы ожидать, так что это вопрос подсчета, а не сопоставления строк, что делает это быстрее.