更快的方式将一个字符串分解和,使用R计算字符数?
-
20-09-2019 - |
题
我在找一个更快的方法来计算GC含量的DNA串从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
用于制备本代码更快任何意见?
解决方案
不如不是在所有分割,就计数匹配:
gcCount2 <- function(line, st, sp){
sum(gregexpr('[GCgc]', substr(line, st, sp))[[1]] > 0)
}
这是大小更快的顺序。
一个小的C函数,仅仅在字符迭代将大小的又一顺序更快。
其他提示
一个一个衬里:
table(strsplit(toupper(a), '')[[1]])
我不知道这是任何更快,但你可能想看看将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()”。
有没有必要在这里使用一个循环。
尝试这种情况:
gcCount <- function(line, st, sp){
chars = strsplit(as.character(line),"")[[1]][st:sp]
length(which(tolower(chars) == "g" | tolower(chars) == "c"))
}
这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
由于所有此信息,
要优化中,我要计算的200bp的序列100M的GC含量的脚本,我在这里结束了提出的测试不同的方法。根·威廉斯的方法进行最佳(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")
不隶属于 StackOverflow