使用 R 分割字符串和计算字符的更快方法?
我正在寻找一种更快的方法来计算从 FASTA 文件读取的 DNA 字符串的 GC 含量。这归结为获取一个字符串并计算字母“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
对于使此代码更快有什么建议吗?
如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱,如果其他用户或者作者回复了您的评论,将不能在第一时间通知您!
发布评论
评论(6)
最好根本不拆分,只计算匹配数:
这会快一个数量级。
一个仅迭代字符的小型 C 函数将会快一个数量级。
Better to not split at all, just count the matches:
That's an order of magnitude faster.
A small C function that just iterates over the characters would be yet another order of magnitude faster.
一班轮:
A one liner:
我不知道它是否更快,但您可能想查看 R 包 seqinR - http://pbil.univ-lyon1.fr/software/seqinr/home.php?lang=eng。它是一个优秀的通用生物信息学软件包,具有多种序列分析方法。它在 CRAN 中(在我写这篇文章时它似乎已经关闭了)。
GC 内容将是:
这是来自字符串,您还可以使用“read.fasta()”读取 fasta 文件。
I don't know that it's any faster, but you might want to look at the R package seqinR - http://pbil.univ-lyon1.fr/software/seqinr/home.php?lang=eng. It is an excellent, general bioinformatics package with many methods for sequence analysis. It's in CRAN (which seems to be down as I write this).
GC content would be:
That's from a string, you can also read in a fasta file using "read.fasta()".
这里不需要使用循环。
试试这个:
There's no need to use a loop here.
Try this:
尝试使用
stringi
包中的这个函数,或者您可以使用正则表达式版本来计算 g 和 G
,或者您可以先使用 tolower 函数,然后使用 stri_count
时间性能
另一个较长字符串的示例。 stri_dup 复制字符串 n 次
如您所见,对于较长的序列 stri_count 更快:)
Try this function from
stringi
packageor you can use regex version to count g and G
or you can use tolower function first and then stri_count
time performance
another example for longer string. stri_dup replicates string n-times
As you can see, for longer sequence stri_count is faster :)
感谢这篇文章的所有人,
为了优化我想要计算 200bp 的 100M 序列的 GC 含量的脚本,我最终测试了此处提出的不同方法。 Ken Williams 的方法表现最好(2.5 小时),优于 seqinr(3.6 小时)。使用 stringr str_count 减少到 1.5 小时。
最后我用 C++ 编写了它,并使用 Rcpp 调用它,这将计算时间减少到了 10 分钟!
这是 C++ 代码:
我从 R 键入中调用它:
Thanks to all for this post,
To optimize a script in which I want to calculate GC content of 100M sequences of 200bp, I ended up testing different methods proposed here. Ken Williams' method performed best (2.5 hours), better than seqinr (3.6 hours). Using stringr str_count reduced to 1.5 hour.
In the end I coded it in C++ and called it using Rcpp, which cuts the computation time down to 10 minutes!
here is the C++ code:
Which I call from R typing: