Подсчет уникальных вхождений в текстовой строке в R

Я использую R, и у меня есть кадр данных, содержащий строки из 4 уникальных букв (ДНК). Меня интересует подсчет количества раз, когда в этих строках встречаются определенные уникальные комбинации букв. Один из возможных сценариев — определить, сколько раз я вижу одну и ту же букву подряд.

Я столкнулся с несколькими возможными способами добиться этого с помощью регулярных выражений и таких пакетов, как stringr, но у меня все еще есть одна проблема.

Эти методы не перебирают подстроку (буква за буквой) и рассматривают следующую букву в строке как соблюдение. Это проблема, когда одна и та же буква повторяется более 2 раз.

Пример (где я хочу подсчитать количество раз, когда появляется «CC», а столбец true_count является моим желаемым результатом):

sequence  stringr_count  true_count
ACCTACGT      1             1
CCCCCCCC      4             7
ACCCGCCT      2             3

person user1607359    schedule 19.06.2015    source источник
comment
это должно помочь stackoverflow.com/questions/7878992/ lengths(gregexpr("(?=CC)", x, perl=T) )   -  person user20650    schedule 19.06.2015
comment
использование lengths(gregexpr((?=CC), x, perl=T) в строках, где нет вхождений определенной комбинации, приводит к возвращению 1. Почему это так и каков хороший обходной путь? ищите подстроку AA в моем примере выше.   -  person user1607359    schedule 19.06.2015
comment
Да, извините, предыдущий комментарий неверен: попробуйте sapply(gregexpr("(?=CC)", x, perl=T) , function(i) length(i[i>0]))   -  person user20650    schedule 19.06.2015
comment
Это помогло. Спасибо за помощь!   -  person user1607359    schedule 19.06.2015
comment
добро пожаловать ... обратите внимание, что ответ с нулями, кажется, также вытаскивает правильные подсчеты   -  person user20650    schedule 19.06.2015
comment
может так лучше sapply(gregexpr("(?=AA)", x, perl=T) , function(i) sum(i>0))   -  person user20650    schedule 19.06.2015
comment
Хотя использование gregexpr работает, это просто смущающе медленно. Пожалуйста, проверьте обновленную версию моего ответа.   -  person zero323    schedule 20.06.2015


Ответы (1)


Я бы рекомендовал использовать stringi::stri_count_fixed следующим образом:

> library(stringi)
> seqs <- data.frame(sequence=c('ACCTACGT', 'CCCCCCCC', 'ACCCGCCT'))
> opts <- stri_opts_fixed(overlap=TRUE)
> seqs$true_count <- stri_count_fixed(str=seqs$sequence, pattern='CC', opts_fixed=opts)
> seqs
  sequence true_count
1 ACCTACGT          1
2 CCCCCCCC          7
3 ACCCGCCT          3

С фиксированным шаблоном stringi выполняется на порядок быстрее, чем при использовании gregexpr:

library(microbenchmark)

# Answer provided by @user20650 in the comments
f1 <- function(x) sapply(gregexpr('(?=CC)', x, perl=T) , function(i) sum(i>0))

f2 <- function(x) stri_count_fixed(
    str=x, pattern='CC',
    opts_fixed=stri_opts_fixed(overlap=TRUE))

# Generate random sequences
sequence <- stri_rand_strings(n=10000, length=1000, pattern='[ATGC]')

Результаты микробенчмарка:

> microbenchmark(f1(sequence), f2(sequence))
Unit: milliseconds
         expr       min        lq      mean    median       uq       max neval
 f1(sequence) 290.90393 304.87107 329.11392 313.39819 327.9860 437.10229   100
 f2(sequence)  20.99733  21.12559  21.39206  21.26017  21.4377  27.68867   100

Вы также можете ознакомиться с библиотекой Biostrings. По моему опыту, обычно это медленнее, чем работа с stringi, и требует некоторых дополнительных шагов, но предоставляет множество полезных функций, предназначенных для работы с биологическими последовательностями, включая countPattern:

library(Biostrings)

bsequence <- DNAStringSet(sequence)
f3 <- function(x) vcountPattern('CC', x)

Результаты микробенчмарка:

> microbenchmark(f2(sequence), f3(bsequence))
Unit: milliseconds
          expr      min       lq     mean   median       uq      max neval
  f2(sequence) 20.83336 21.11473 21.36759 21.25088 21.45000 23.80708   100
 f3(bsequence) 86.95430 89.10023 89.51665 89.37103 89.87699 91.88203   100

И просто для уверенности:

> identical(f1(seqs$sequence), f2(seqs$sequence),  f3(BStringSet(seqs$sequence))
[1] TRUE
person zero323    schedule 19.06.2015