¿Forma más rápida de dividir una cadena y contar caracteres usando R?

Question

Estoy buscando una forma más rápida de calcular el contenido de GC para cadenas de ADN leídas desde un archivo FASTA. Esto se reduce a tomar una cadena y contar el número de veces que aparece la letra 'G' o 'C'. También quiero especificar el rango de caracteres a considerar.

Tengo una función de trabajo que es bastante lenta y está causando un cuello de botella en mi código. Se ve así:

## 
## 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) 
} 

Correr Rprof me da el siguiente resultado:

> 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 

Algún consejo para hacer el código más rápido?

Answer 1

Mejor divididos en absoluto, simplemente contar los partidos

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

Eso es un orden de magnitud más rápido .

Una pequeña función C que simplemente itera sobre los caracteres sería otro orden de magnitud más rápido.

Answer 2

No es necesario utilizar un bucle aquí.

Prueba esto:

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

Answer 3

Un chiste:

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

Answer 4

No sé si es más rápido, pero es posible que desee ver el paquete R seqinR - http://pbil.univ-lyon1.fr/software/seqinr/home.php?lang=eng. Es un excelente paquete de bioinformática general con muchos métodos para el análisis de secuencias. Está en CRAN (que parece estar caído mientras escribo esto).

contenido de GC sería:

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

Eso es de una cadena, también se puede leer en un archivo FASTA usando "read.fasta()".

Answer 5

Prueba esta función desde stringi paquete

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

o puede utilizar la versión de expresiones regulares para contar gy G

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

o puede utilizar tolower funcionar primero y luego stri_count

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

tiempo de rendimiento

> 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 

otro ejemplo de cadena más larga.stri_dup replica cadena de n-veces

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

Como se puede ver, por más stri_count secuencia es más rápido :)

> 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 6

Gracias a todos por este post,

Para optimizar un guión en el que Quiero calcular el contenido de GC de secuencias 100M de 200 pb, terminé probando diferentes métodos propuestos aquí. El método de Ken Williams se comportó mejor (2.5 horas), mejor que el seqinr (3.6 horas). Usar stringr str_count reducido a 1.5 hora.

¡Al final lo codifiqué en C++ y lo llamé usando Rcpp, que reduce el tiempo de cálculo a 10 minutos!

aquí es el código 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; 
} 

que llamo de R escribiendo:

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