Более быстрый способ разбить строку и посчитать символы, используя R?

Question

Я ищу более быстрый способ расчета содержания GC для строк ДНК, считанных из файла FASTA. Это сводится к тому, чтобы взять строку и посчитать, сколько раз появляется буква «G» или «C». Я также хочу указать диапазон символов для рассмотрения.

У меня есть работающая функция, которая работает довольно медленно, и это вызывает узкое место в моем коде. Это выглядит так:

##
## count the number of GCs in the characters between start and stop
##
gcCount <-  function(line, st, sp){
  chars = strsplit(as.character(line),"")[[1]]
  numGC = 0
  for(j in st:sp){
    ##nested ifs faster than an OR (|) construction
    if(chars[[j]] == "g"){
      numGC <- numGC + 1
    }else if(chars[[j]] == "G"){
      numGC <- numGC + 1
    }else if(chars[[j]] == "c"){
      numGC <- numGC + 1
    }else if(chars[[j]] == "C"){
      numGC <- numGC + 1
    }
  }
  return(numGC)
}

Запуск Rprof дает мне следующий вывод:

> a = "GCCCAAAATTTTCCGGatttaagcagacataaattcgagg"
> Rprof(filename="Rprof.out")
> for(i in 1:500000){gcCount(a,1,40)};
> Rprof(NULL)
> summaryRprof(filename="Rprof.out")

                   self.time self.pct total.time total.pct
"gcCount"          77.36     76.8     100.74     100.0
"=="               18.30     18.2      18.30      18.2
"strsplit"          3.58      3.6       3.64       3.6
"+"                 1.14      1.1       1.14       1.1
":"                 0.30      0.3       0.30       0.3
"as.logical"        0.04      0.0       0.04       0.0
"as.character"      0.02      0.0       0.02       0.0

$by.total
               total.time total.pct self.time self.pct
"gcCount"          100.74     100.0     77.36     76.8
"=="                18.30      18.2     18.30     18.2
"strsplit"           3.64       3.6      3.58      3.6
"+"                  1.14       1.1      1.14      1.1
":"                  0.30       0.3      0.30      0.3
"as.logical"         0.04       0.0      0.04      0.0
"as.character"       0.02       0.0      0.02      0.0

$sampling.time
[1] 100.74

Какой совет, как сделать этот код быстрее?

Answer 1

Лучше вообще не разбивать, просто посчитай совпадения:

gcCount2 <-  function(line, st, sp){
  sum(gregexpr('[GCgc]', substr(line, st, sp))[[1]] > 0)
}

Это на порядок быстрее.

Маленькая функция C, которая просто перебирает символы, будет еще на порядок быстрее.

Answer 2

Один вкладыш:

table(strsplit(toupper(a), '')[[1]])

Answer 3

Я не знаю, что это быстрее, но вы можете посмотреть на пакет R seqinR - http://pbil.univ -lyon1.fr / software / seqinr / home.php? Lang = eng, Это отличный общий пакет биоинформатики с множеством методов для анализа последовательности. Это в CRAN (который, кажется, не работает, когда я пишу это).

Содержание GC будет:

mysequence <- s2c("agtctggggggccccttttaagtagatagatagctagtcgta")
    GC(mysequence)  # 0.4761905

Это из строки, вы также можете прочитать в файле fasta, используя read.fasta ().

Answer 4

Попробуйте эту функцию из stringi пакета

> stri_count_fixed("GCCCAAAATTTTCCGG",c("G","C"))
[1] 3 5

или вы можете использовать регулярное выражение для подсчета g и G

> stri_count_regex("GCCCAAAATTTTCCGGggcc",c("G|g|C|c"))
[1] 12

или вы можете сначала использовать функцию tolower, а затем stri_count

> stri_trans_tolower("GCCCAAAATTTTCCGGggcc")
[1] "gcccaaaattttccggggcc"

время исполнения

    > microbenchmark(gcCount(x,1,40),gcCount2(x,1,40), stri_count_regex(x,c("[GgCc]")))
Unit: microseconds
                             expr     min     lq  median      uq     max neval
                gcCount(x, 1, 40) 109.568 112.42 113.771 116.473 146.492   100
               gcCount2(x, 1, 40)  15.010  16.51  18.312  19.213  40.826   100
 stri_count_regex(x, c("[GgCc]"))  15.610  16.51  18.912  20.112  61.239   100

еще один пример для более длинной строки. stri_dup копирует строку n раз

> stri_dup("abc",3)
[1] "abcabcabc"

Как видите, для более длинной последовательности stri_count работает быстрее:)

> y <- stri_dup("GCCCAAAATTTTCCGGatttaagcagacataaattcgagg",100)
    > microbenchmark(gcCount(y,1,40*100),gcCount2(y,1,40*100), stri_count_regex(y,c("[GgCc]")))
    Unit: microseconds
                                 expr       min         lq     median        uq       max neval
              gcCount(y, 1, 40 * 100) 10367.880 10597.5235 10744.4655 11655.685 12523.828   100
             gcCount2(y, 1, 40 * 100)   360.225   369.5315   383.6400   399.100   438.274   100
     stri_count_regex(y, c("[GgCc]"))   131.483   137.9370   151.8955   176.511   221.839   100

Answer 5

Здесь нет необходимости использовать цикл.

Попробуйте это:

gcCount <-  function(line, st, sp){
  chars = strsplit(as.character(line),"")[[1]][st:sp]
  length(which(tolower(chars) == "g" | tolower(chars) == "c"))
}

Answer 6

Спасибо всем за этот пост,

Чтобы оптимизировать сценарий, в котором я хочу рассчитать содержание GC для 100M последовательностей по 200 бп, я закончил тестировать различные методы, предложенные здесь. Метод Кена Уильямса показал себя лучше (2,5 часа), лучше, чем Seqinr (3,6 часа). Использование stringr str_count сокращено до 1,5 часа.

В конце я написал код на C ++ и вызвал его с помощью Rcpp, что сокращает время вычислений до 10 минут!

вот код C ++:

#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
float pGC_cpp(std::string s) {
  int count = 0;

  for (int i = 0; i < s.size(); i++) 
    if (s[i] == 'G') count++;
    else if (s[i] == 'C') count++;

  float pGC = (float)count / s.size();
  pGC = pGC * 100;
  return pGC;
}

Который я звоню из R набрав:

sourceCpp("pGC_cpp.cpp")
pGC_cpp("ATGCCC")