Schnellerer Weg, um eine Zeichenfolge zu teilen und Zeichen mit R zu zählen?

Question

Ich suche einen schnelleren Weg, um den GC -Inhalt für DNA -Strings aus einer Fasta -Datei zu berechnen. Dies läuft darauf hinaus, eine Zeichenfolge zu nehmen und die Häufigkeit zu zählen, in der der Buchstabe 'G' oder 'C' erscheint. Ich möchte auch den zu berücksichtigenden Zeichenbereich angeben.

Ich habe eine Arbeitsfunktion, die ziemlich langsam ist, und sie verursacht einen Engpass in meinem Code. Es sieht aus wie das:

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

Das Ausführen von RPROF gibt mir die folgende Ausgabe:

> 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

Irgendwelche Ratschläge, um diesen Code schneller zu machen?

Answer 1

Besser nicht aufgeteilt, nur die Übereinstimmungen zählen:

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

Das ist eine Größenordnung schneller.

Eine kleine C -Funktion, die nur über die Zeichen iteriert, wäre eine weitere Größenordnung schneller.

Answer 2

Ein One Liner:

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

Answer 3

Ich weiß nicht, dass es schneller ist, aber vielleicht möchten Sie sich das R -Paket seqinr - ansehen - http://pbil.univ-lyon1.fr/software/seqinr/home.php?lang=eng. Es ist ein ausgezeichnetes allgemeines Bioinformatikpaket mit vielen Methoden zur Sequenzanalyse. Es ist in Cran (was scheinbar unten zu sein scheint, wenn ich das schreibe).

GC -Inhalt wäre:

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

Das stammt aus einer Zeichenfolge, Sie können auch in einer Fasta -Datei mit "read.fasta ()" lesen.

Answer 4

Hier müssen keine Schleife verwendet werden.

Versuche dies:

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

Answer 5

Versuchen Sie diese Funktion von stringi Paket

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

Oder Sie können die Regex -Version verwenden, um g und g zu zählen

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

Oder Sie können zuerst die Tolower -Funktion und dann STRI_COUNT verwenden

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

Zeitleistung

    > 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

Ein weiteres Beispiel für eine längere Zeichenfolge. stri_dup repliziert die String n-z-Zeit

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

Wie Sie sehen können, ist für eine längere Sequenz STRI_COUNT schneller :)

> 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

Vielen Dank an alle für diesen Beitrag,

Um ein Skript zu optimieren, in dem ich den GC -Inhalt von 100 m -Sequenzen von 200 bp berechnen möchte, habe ich hier verschiedene Methoden getestet, die hier vorgeschlagen wurden. Die Methode von Ken Williams zeigte sich am besten (2,5 Stunden), besser als Seqinr (3,6 Stunden). Verwenden von Stringr str_count auf 1,5 Stunden reduziert.

Am Ende habe ich es in C ++ codiert und es mit RCPP aufgerufen, was die Rechenzeit auf 10 Minuten verkürzt!

Hier ist der C ++ - Code:

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

Was ich vom R -Tippen rufe:

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