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

Question

我有兴趣从一个序列中生成所有有间隙的 kmers，其中有间隙的 kmer 定义为一个长度为 k 的序列，与另一个长度为 k 的序列最多分隔 m 个位置。因此，例如，“序列 CAGAT 间隙对 k = 1 和 m = 2 的内核找到其间具有零到两个不相关位置的单体对。即它找到特征 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 方面的经验有限，我不确定如何改进它。

执行此操作的最有效方法是什么，同时确保所有预期的子序列(例如:CAGAT -> CA、C.G、C..A、AG、A.A、A..T、GA、G.T 和 AT k=1, m=2) 包含在输出中？

谢谢!

Answer 1

您可能想在 Bioconductor kebabs package 中查看 gappy 对内核的实现。 .

安装:

## 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)

k-mers 存储在 Dimnames 槽中:

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