在R中存储长字符串（DNA序列）

Question

我已经编写了一个函数，用于查找长DNA序列中子序列的索引。当我的较长DNA序列长度小于约4000个字符时，它可以正常工作。然而，当我尝试将相同的函数应用于更长的序列时，控制台会给出一个+而不是一个>... 这使我相信问题在于字符串的长度。

例如：当较长的序列为："GATATATGCATATACTT"，而子序列为："ATAT"时，我得到的索引是"1, 3, 9"（基于0）

dnaMatch <- function(dna, sequence) {
  ret <- list()
  k <- str_length(sequence)
  c <- str_length(dna) - k 
  for(i in 1:(c+1)) {
    ret[i] = str_sub(dna, i, i+k-1)
  }
  ret <- unlist(ret)
  TFret <- lapply (ret, identical, sequence)
  TFret <- which(unlist(TFret), arr.ind = TRUE) -1
  print(TFret)
}

基本上，我的问题是...在字符串类中有没有绕过字符限制的方法？

- Zoe

---

控制台给我一个 + 而不是 > - 这通常是由于缺少 ), ", ' 或 } 导致的。你有没有仔细检查过代码中的拼写错误？ - nrussell

我明白了，我通过运行 cat(paste0(rep("abcdef",1000),collapse="")) 复制输出并尝试将其分配给对象 x <- "[paste copied text]"，成功复现了这个问题。很可能是由于某种缓冲区限制导致的。 - nrussell

我可以重现你的例子，nrussell，但是以下代码可以正确赋值：x<-paste0(rep("abcdef",1000),collapse="")。 - bjoseph

她在这种情况下也可以有效地使用正则表达式匹配，例如 ?gregxpr 或 ?regexpr。详见 http://rfunction.com/archives/1719 和 http://stat.ethz.ch/R-manual/R-devel/library/base/html/regex.html。 - bjoseph

使用dnaMatch(x,'ATATAT')和gregexpr('ATATAT',x)在一个长度为100k的字符串上匹配，分别需要49秒和0.003秒。抱歉，我没有耐心再运行更长时间的循环。 - rawr

Answer 1

我可以复制nrussell的示例，但这会正确地分配 x<-paste0(rep("abcdef",1000),collapse="") -- 一种潜在的解决方法是将字符字符串写入.txt文件并直接读取.txt文件到R中：

test.txt是一个6,000个字符长的字符串。

`test<-read.table('test.txt',stringsAsFactors = FALSE)
 length(class(test[1,1]))
[1] 1
class(test[1,1])
[1] "character"
 nchar(test[1,1])
[1] 6000`

- bjoseph

---

谢谢 - 这就是我最终做的！ - Zoe

太棒了！你应该研究一下正则表达式匹配。在这种情况下，它可以拯救生命。 - bjoseph

Answer 2

不要自己编写函数，为什么不使用seqinr包中的words.pos函数呢？它似乎可以处理长度达到一百万个碱基对的字符串。

例如，

library(seqinr)
data(ec999)
myseq <- paste(ec999[[1]], collapse="")
myseq <- paste(rep(myseq,100), collapse="")
words.pos("atat", myseq)

- J. Win.

---

https://github.com/cran/seqinr/blob/master/R/words.pos.R 该代码使用了 while 循环、对象增长和 substr 函数。可能不是最快的方法，而且只是 regexpr 函数的包装器。 - rawr

不是最快的，但她的例子表明她想找到非不相交的匹配。 - J. Win.

基本上，我的问题是...在字符串类中有没有绕过字符限制的方法？我认为她的问题实际上是关于R解释器缓冲区的问题，即它一次只能处理有限量的文本输入到控制台，而不是关于她的函数性能的问题。 - nrussell