使用R分割字符串和計數字符的更快方法？

我正在尋找更快速的方式來計算從FASTA文件中讀取的DNA字符串的GC內容。這可以歸結爲一個字符串並計算字母「G」或「C」出現的次數。我也想指定要考慮的字符範圍。使用R分割字符串和計數字符的更快方法？

我有一個相當慢的工作函數，它在我的代碼中造成瓶頸。它看起來像這樣：

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

使這一代碼快有什麼建議？

來源

2010-03-15 chrisamiller

對於它的價值，我最終決定的是，R是要處理約3十億乾脆太慢來自人類基因組的鹼基對，並使用一個小的Perl腳本來代替。 – chrisamiller 2011-11-04 16:27:09

更好

不拆所有，只是算上比賽：

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

這是一個數量級的速度更快。

只是迭代字符的一個小C函數會更快一個數量級。

來源

2010-03-15 18:20:40

更好（約快7倍）。謝謝！ – chrisamiller 2010-03-15 18:54:06

此優惠的重要補充 - 請注意，如果子字符串不包含G \ C，則長度函數可能會返回（-1）而不是0，因此需要對其進行檢查。 – dan12345 2011-08-16 13:05:44

感謝dan12345 - @ user2265478剛剛建議編輯修復此問題，我將其納入（即使該編輯被拒絕[不是由我]）。 – 2016-02-19 20:57:52

這裏沒有必要使用循環。

試試這個：

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

來源

2010-03-15 17:34:54 Shane

已投票。比我的回答更好（它會是這樣的：在BioPerl中做:-)） – 2010-03-15 17:39:24

非常感謝。這個速度大約快4倍，這僅僅是我圍繞Rajarshi的代碼構建的函數。你可以說我仍然在學習R - 很難擺脫這種多年來我一直使用的以循環爲中心的思維方式。 – chrisamiller 2010-03-15 17:57:54

您可能會嘗試的另一件事是：％c（「g」，「c」）中的tolower（chars）％。不知道哪個更快，儘管我懷疑OR'|'運算符比'％in％'快。 – Shane 2010-03-15 18:11:05

一個一個班輪：

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

來源

2010-03-15 17:37:32 rguha

我不知道它有什麼更快，但你可能想看看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（）」。

來源

2010-03-16 08:13:35 neilfws

從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

來源

2014-03-14 10:13:13 bartektartanus

感謝所有爲這個職位，

要優化腳本中我想計算200bp的100M序列的GC含量，我最終測試了這裏提出的不同方法。肯威廉姆斯的方法表現最好（2.5小時），優於賽勤爾（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")

來源

2017-10-16 20:01:47

使用R分割字符串和計數字符的更快方法？

回答

相關問題