Cách nhanh hơn để chia chuỗi và đếm ký tự bằng R?

Question

Tôi đang tìm cách tính toán nội dung GC cho chuỗi DNA được đọc từ tệp FASTA nhanh hơn. Điều này nhún xuống để lấy một chuỗi và đếm số lần chữ 'G' hoặc 'C' xuất hiện. Tôi cũng muốn chỉ định phạm vi các ký tự cần xem xét.

Tôi có chức năng hoạt động khá chậm và gây tắc nghẽn trong mã của tôi. Nó trông giống như thế này:

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

Chạy Rprof mang lại cho tôi kết quả như sau:

> 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 

Bất cứ lời khuyên để làm mã này nhanh hơn?

Answer 1

Tốt hơn để không chia ở tất cả, chỉ cần đếm các trận đấu:

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

Đó là một thứ tự cường độ nhanh .

Chức năng C nhỏ chỉ lặp qua các ký tự sẽ là một thứ tự cường độ nhanh hơn.

Answer 2

Không cần sử dụng vòng lặp tại đây.

Hãy thử điều này:

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

Answer 3

Một lót:

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

Answer 4

Tôi không biết rằng nó nhanh hơn, nhưng bạn có thể muốn xem gói R seqinR - http://pbil.univ-lyon1.fr/software/seqinr/home.php?lang=eng. Đây là một gói tin sinh học tổng quát, tuyệt vời với nhiều phương pháp phân tích trình tự. Đó là trong CRAN (mà dường như là xuống như tôi viết này).

nội dung GC sẽ là:

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

Đó là từ một chuỗi, bạn cũng có thể đọc trong một tập tin FASTA sử dụng "read.fasta()".

Answer 5

Hãy thử chức năng này từ stringi gói

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

hoặc bạn có thể sử dụng phiên bản regex để đếm g và G

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

hoặc bạn có thể sử dụng ToLower hoạt động đầu tiên và sau đó stri_count

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

hiệu suất thời gian

> 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 

một ví dụ khác cho chuỗi dài hơn.stri_dup tái tạo chuỗi n-lần

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

Như bạn có thể thấy, trong thời gian dài chuỗi stri_count nhanh :)

> 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

Nhờ tất cả cho bài này,

Để tối ưu hóa một kịch bản trong đó Tôi muốn tính hàm lượng GC của 100 triệu chuỗi 200bp, tôi đã kết thúc thử nghiệm các phương pháp khác nhau được đề xuất ở đây. Phương pháp của Ken Williams được thực hiện tốt nhất (2,5 giờ), tốt hơn so với seqinr (3,6 giờ). Sử dụng stringr str_count giảm xuống còn 1,5 giờ.

Cuối cùng, tôi đã mã hóa nó bằng C++ và gọi nó bằng Rcpp, giúp giảm thời gian tính toán xuống còn 10 phút!

ở đây là 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; 
} 

Mà tôi gọi từ R gõ:

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