R的字符串处理很慢?

yihui

放眼望去，遍地循环，能快得了么：）

你用substr()一个一个字符看当然慢啊，那么多次重复操作同一个字符串，何不一下拆开，向量化一口气替换完毕……

我不懂什么ATCG，高中学的早忘光了，但意思应该是AT互换，CG互换吧。我拿一个长度为三百万的字符串试验了一下，完成互换消耗的时间为0：

seqStr = paste(sample(c("A", "T", "C", "G"), 3e+06, <br />
    TRUE), collapse = "")<br />
system.time({<br />
    seqChar = strsplit(seqStr, "")[[1]]<br />
    seqName = c(A = 1, T = 2, C = 3, G = 4)<br />
    seqNum = seqName[seqChar]<br />
    seqSwitch = ifelse(seqNum%%2 == 0, seqNum - 1, seqNum + 1)<br />
    seqResChar = names(seqName)[seqSwitch]<br />
    seqResStr = paste(seqResChar, collapse = "")<br />
}) <br />
# 消耗时间<br />
#   user  system elapsed <br />
#      0       0       0 <br />
# 比如<br />
#> seqStr<br />
#[1] "AACAGTTCCGACACAGAGGC"<br />
#> seqResStr<br />
#[1] "TTGTCAAGGCTGTGTCTCCG"

如果你看不懂这段代码，可以把第一句的3e+06改小一些，比如20，也就是长度为20的字符串看看结果seqResStr是否正确。这个问题的关键在于，向量可以通过整数下标的方式取子集，也可以用字符名字的方式取子集。

yanlinlin82

谢老大的电脑果然比我的快，能跑到0！

<br />
seqStr <- paste(sample(c("A", "T", "C", "G"), 3e6, TRUE), collapse = "")<br />
system.time({<br />
    seqChar = strsplit(seqStr, "")[[1]]<br />
    seqName = c(A = 1, T = 2, C = 3, G = 4)<br />
    seqNum = seqName[seqChar]<br />
    seqSwitch = ifelse(seqNum%%2 == 0, seqNum - 1, seqNum + 1)<br />
    seqResChar = names(seqName)[seqSwitch]<br />
    seqResStr = paste(seqResChar, collapse = "")<br />
}) <br />
#   user  system elapsed <br />
#  5.860   0.210   6.459 <br />

我也学谢老大来变变“魔术”（经老大提醒，此处代码有错误，不能正确处理反转，正确代码见9楼后文）：

<br />
reverseAndComplement <- function(strseq, doReverse = T, doComplement = T) {<br />
    if (!doReverse) rev <- c<br />
    ifelse(doComplement, paste(c(A = 'T', T = 'A', C = 'G', G = 'C')[rev(strsplit(toupper(strseq),"")[[1]])], collapse = ""), rev(s))<br />
} <br />
system.time(reverseAndComplement(seqStr))<br />
#   user  system elapsed <br />
#  0.900   0.000   0.902 <br />

似乎能稍微快点。

iiiiiiiiiii

谢谢楼上两位,

我马上试一下. 我原来的写法也让我很不安,只是当时确实没啥其它想法.

yihui

你这个if语句实在是写得惊世骇俗！

下面的魔术看得我更加目瞪口呆，跟看刘谦的魔术似的……高，实在是高

另：ifelse()中的rev(s)疑为笔误……而且如果只反转字符串的话，结果就不对了，rev()只能反转向量，对单个的字符串中的字符无法反转。

iiiiiiiiiii

 <br />
 <br />
reverseAndComplement <- function(strseq) <br />
{ <br />
  paste(c(A='T',T='A',C='G',G='C')[rev(strsplit(toupper(strseq),'')[[1]])],collapse='')<br />
}<br />
seqStr <- paste(sample(c("A", "T", "C", "G"), 3e6, TRUE), collapse = "")<br />
system.time(reverseAndComplement(seqStr))<br />
<br />
这个机子的系统是最烂的vista.4G<br />
用户 系统 流逝 <br />
1.57 0.12 1.47 <br />
<br />
cluster上的比我自己的机子还慢..<br />
   user  system elapsed <br />
  1.983   0.237   2.222 <br />
<br />
<br />
再试了下1楼我写的慢程序:<br />
system.time(reverseAndComplement(seqStr))<br />
真的很慢~~ <br />
(好吧,我回完帖了,这结果还没出来,我就不帖了 -.- )<br />
<br />

处理小数据的时候,没啥感觉代码效率. 真到大数据的时候,真的很重要啊.

谢谢各位了. 学到很多~~

iiiiiiiiiii

那个if我完全没看懂,照搬来用着解决实际问题了先~~

ryusukekenji

谢老师和yanlinlin82 前辈的程序完全看不懂，太强了！！！

<br />
>seqChar = strsplit(seqStr, "")[[1]]<br />

看了好几遍，运行了一下还是摸不着头脑，seqStr经过strsplit后是list 吗？为什么要加上[[1]] 呢？

<br />
reverseAndComplement <- function(strseq, doReverse = T,<br />
    doComplement = T) {<br />
    if (!doReverse)<br />
        rev <- c<br />
    ifelse(doComplement, paste(c(A = "T", T = "A", C = "G", G = "C")[rev(strsplit(toupper(strseq),<br />
        "")[[1]])], collapse = ""), rev(s))<br />
}<br />
system.time(reverseAndComplement(seqStr))<br />
#   user  system elapsed<br />
#  0.900   0.000   0.902<br />

rev <- c 和 rev(s) 也看不懂？！ @_@

yanlinlin82

幸好被你指出，果然是笔误，导致无法正确处理反转。

改进的代码如下（这次还加了个names的技巧）：

<br />
reverseAndComplement <- function(strseq, doReverse = T, doComplement = T) {<br />
    map <- c('A','C','G','T')<br />
    names(map) <- ifelse(rep(doComplement, length(map)), rev(map), map)<br />
    if (!doReverse) rev <- c<br />
    paste(map[rev(strsplit(toupper(strseq), "")[[1]])], collapse = "")<br />
}<br />
strseq <- "CAGAACATTTAGCGGTAAAA"<br />
reverseAndComplement(strseq)       # "TTTTACCGCTAAATGTTCTG"<br />
reverseAndComplement(strseq, F)    # "GTCTTGTAAATCGCCATTTT"<br />
reverseAndComplement(strseq, T, F) # "AAAATGGCGATTTACAAGAC"<br />
reverseAndComplement(strseq, F, F) # "CAGAACATTTAGCGGTAAAA"<br />

结果看起来似乎很好。

那个“惊世骇俗”的“if (!doReverse) rev <- c”，我在开始写时就特别想笑，这纯粹是钻了R语言的空子，用一个局部的rev变量遮盖掉系统的rev函数，从而使反转失效。此法纯属搞笑，如非万不得已，尤其是在头脑不清醒时，请勿模仿使用，否则后果自负，哈哈 [s:11]

yanlinlin82

是的，seqStr经过strsplit后是list，所以用[[1]]来提取第一项。strsplit其实能够同时切分多个字符串，因此要用列表来存放每个字符串的切分结果。

ryusukekenji

[quote]引用第9楼yanlinlin82于2010-02-20 01:19发表的回 7楼(ryusukekenji) 的帖子 :

是的，seqStr经过strsplit后是list，所以用[[1]]来提取第一项。strsplit其实能够同时切分多个字符串，因此要用列表来存放每个字符串的切分结果。

[/quote]

> test <- c('A','B','C')<br />
> strsplit(test,'')<br />
[[1]]<br />
[1] "A"<br />
<br />
[[2]]<br />
[1] "B"<br />
<br />
[[3]]<br />
[1] "C"<br />

if (!doReverse) rev <- c 哈哈哈哈，了解了。strsplit 也了解了。谢谢yanlinlin82版主。

yihui

R的随意性跟C比起来就像个村夫，C一向都是西装笔挺，R灰头土脸云游四方

yanlinlin82

其实C也有夸张的要命的写法：

http://www.cnblogs.com/exclm/archive/2008/12/18/1357586.html

http://en.wikipedia.org/wiki/International_Obfuscated_C_Code_Contest

随意性与语言无关，而与人有关。所以才有那句话：好的代码是先给人读，然后才给机器读的。

记得老大曾经不时地给出一些关于写R程序的规则，诸如避免循环、避免attach之类的，确实为金玉良言。现在，为避免误导公众，我想“rev<-c”也该被收录列为禁忌之一了。 [s:13]

yihui

发现一个简单的调换字符的方法：

> chartr("ATCG", "TAGC", "CAGAACATTTAGCGGTAAAA")<br />
[1] "GTCTTGTAAATCGCCATTTT"<br />

yanlinlin82

回复第15楼的谢益辉：又领教了！一比较，前面的所有代码都显得很傻很天真了……

ryusukekenji

回复第14楼的 yanlinlin82：

#记得老大曾经不时地给出一些关于写R程序的规则，诸如避免循环、避免attach之类的，确实为金玉良言....

哦，attach原来也会拖慢速度吗？在R-help中还常看到一些编码attach了数据才运算。

回复第15楼的谢益辉：看了一下chartr的代码，才几个字竟然就那么好用，厉害。

yihui

attach()的问题不在于速度，而是它会把命名空间搞乱，可能会给一些人带来混淆。例如：

> x = sleep<br />
> attach(x)<br />
# 对x中的变量操作<br />
> extra = extra + 1<br />
# 但x本身的变量并没有受影响<br />
> identical(x$extra, extra)<br />
[1] FALSE<br />
> detach()<br />

为了不要把环境弄乱，一般尽量用with(data, {...})，比如

with(sleep, plot(extra ~ group))

这样能拥有attach()的优点（不必写dataname$varname，只需写varname），而且不会对命名空间产生影响。

当然，萝卜白菜各有所爱，也有很多人愿意用attach()，怎么方便怎么来。

ryusukekenji

回复第18楼的谢益辉：哦，原来是这样子，看到一些洋人研究人员都爱使用attach。这样的话不就出现另一个问题了？举个下列的例子，那么我们是否就该避免使用同样的变量名了？

>#例子<br />
>> x <- c(10:15)<br />
> smp <- data.frame(x=c(1:5),y=c(6:10))<br />
> x<br />
[1] 10 11 12 13 14 15<br />
> attach(smp)</p>
<p>        The following object(s) are masked _by_ .GlobalEnv :</p>
<p>         x<br />

guo.bailing

回复第18楼的谢益辉：正被这个问题困扰，感谢谢老大

月珥

完了我觉得我的C思想已经根深蒂固了。看到这贴首先想到的是：

直接把每个ATCG直接储存成两位二进制数，比如说令A=00，T=10，C=01，G=11，然后合并起来，得到一个00100111...这样的二进制链，最后每隔一位取反就行……

不知R里头做类似的工作方便不……

suckbunny

回复第21楼的月珥：用1减偶数项就好了。