从字符串生成所有的有间隔k-mer序列

Question

我想制作一个从序列中生成所有间隔k-mer的程序。其中，间隔k-mer被定义为长度为k的序列与长度为k的另一条序列之间最多相隔m个位置的序列。例如，“序列CAGAT”, 当核对k=1和m=2的间隔k-mer时，可以找到具有零到两个不相关位置的单体对。也就是说，它可以找到特征CA，C.G，C..A，AG，A.A，A..T，GA，G.T和AT。

replacefxn <- function(x, k, m) {
  substr(x, k + 1, k + m) <- paste(rep("X", m), collapse = "")
  return(x)
}
gappedkmersfxn <- function(x, k, m) {
  n <- (2 * k + m)
  subseq <-
    substring(x, seq(from = 1, to = (nchar(x) - n + 1)), seq(from = n, to = nchar(x)))
  return(sapply(subseq, replacefxn, k, m))
}



allgappedkmersfxn <- function(x, k, m) {
  kmers <- list()
  for (i in 0:m) {
    kmers[[i]] <- gappedkmersfxn(x, k, i)
  }
  kmers <- unlist(kmers)
  return(kmers)
}

allgappedkmersfxn是我目前实现的方式，但它不会添加没有间隙的特征（m是最大间隙，但从0到m），因此无法给出所有所需的特征（例如“CAGAT”的示例）。此外，在一次处理数百万个序列时，它非常缓慢和低效。它的代码也很糟糕，但由于我在R方面经验有限，我不知道如何改进它。

在确保输出包含所有期望子序列（例如：对于k=1，m=2，CAGAT -> CA、C.G、C..A、AG、A.A、A..T、GA、G.T和AT）的情况下，最有效的方法是什么？

谢谢！

- Jack Arnestad

Answer 1

你可以查看Bioconductor kebabs package中gappy pair kernel的实现。

安装方法：

## try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
biocLite("kebabs")

生成一个k = 1，m = 2内核的方法如下：

library(kebabs)
gappyK1M2 <- gappyPairKernel(k = 1, m = 2)

从DNA序列生成明确的表示：

dnaseqs <- DNAStringSet("CAGAT")
dnaseqsrep <- getExRep(dnaseqs, gappyK1M2)

dnaseqsrep@Dimnames[[2]]
[1] "A.A"  "AG"   "AT"   "A..T" "CA"   "C..A" "C.G"  "GA"   "G.T" 

- neilfws