r - 在 R 中有效地将多个字符串/关键字与多个文本匹配

Question

我正在尝试有效地将精确的肽(26 个字符字母表 A-Z¹ 中的氨基酸短序列)映射到蛋白质(同一字母表的较长序列)。据我所知，最有效的方法是使用 Aho-Corasick trie(其中肽是关键字)。不幸的是，我在 R 中找不到适用于非核苷酸字母表的 AC 版本(Biostrings 的 PDict 和 Starr 的 match_ac 都是 DNA 硬编码) .

作为拐杖，我一直在尝试并行化基本的 grep 方法。但我很难找到一种方法来做到这一点而不产生大量的 IO 开销。这是一个简短的示例:

peptides = c("FSSSGGGGGGGR","GAHLQGGAK","GGSGGSYGGGGSGGGYGGGSGSR","IISNASCTTNCLAPLAK")
if (!exists("proteins"))
{
  biocLite("biomaRt", ask=F, suppressUpdates=T, suppressAutoUpdate=T)
  library(biomaRt)
  ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
  proteins = getBM(attributes=c('peptide', 'refseq_peptide'), filters='refseq_peptide', values=c("NP_000217", "NP_001276675"), mart=ensembl)
  row.names(proteins) = proteins$refseq_peptide
}

library(snowfall)
library(Biostrings)
library(plyr)
sfInit(parallel=T, cpus=detectCores()-1)

allPeptideInstances = NULL
i=1
increment=100
count=nrow(proteins)
while(T)
{
  print(paste(i, min(count, i+increment), sep=":"))
  text_source = proteins[i:min(count, i+increment),]
  text = text_source$peptide

  #peptideInstances = sapply(peptides, regexpr, text, fixed=T, useBytes=T)
  peptideInstances = sfSapply(peptides, regexpr, text, fixed=T, useBytes=T)
  dimnames(peptideInstances) = list(text_source$refseq_peptide, colnames(peptideInstances))

  sparsePeptideInstances = alply(peptideInstances, 2, .fun = function(x) {x[x > 0]}, .dims = T)

  allPeptideInstances = c(allPeptideInstances, sparsePeptideInstances, recursive=T)
  if (i==count | nrow(text_source) < increment)
    break
  i = i+increment
}

sfStop()

这里有几个问题:

• peptideInstances 这里是一个稠密矩阵，所以从每个 worker 那里返回它是非常冗长的。我已经把它拆散了分成 block ，这样我就不会处理 40,000(蛋白质)x 60,000(肽)矩阵。
• 在肽上进行并行化，当它可以使并行化蛋白质更有意义，因为它们更大。但我对尝试通过蛋白质来做到这一点感到沮丧，因为:
• 如果 text_source 中只有一种蛋白质，此代码就会中断。

或者，如果有人知道 R 中有更好的解决方案，我很乐意使用它。我已经在这方面花费了足够的时间，实现 Aho-Corasick 可能会得到更好的服务。

¹ 其中一些是歧义代码，但为了简单起见，请忽略它们。

Answer 1

我学习了 Rcpp 并自己实现了 Aho-Corasick。现在CRAN有一个很好的通用多关键字搜索package .

以下是一些使用示例:

listEquals = function(a, b) { is.null(unlist(a)) && is.null(unlist(b)) || !is.null(a) && !is.null(b) && all(unlist(a) == unlist(b)) }

# simple search of multiple keywords in a single text
keywords = c("Abra", "cadabra", "is", "the", "Magic", "Word")
oneSearch = AhoCorasickSearch(keywords, "Is Abracadabra the Magic Word?")
stopifnot(listEquals(oneSearch[[1]][[1]], list(keyword="Abra", offset=4)))
stopifnot(listEquals(oneSearch[[1]][[2]], list(keyword="cadabra", offset=8)))
stopifnot(listEquals(oneSearch[[1]][[3]], list(keyword="the", offset=16)))
stopifnot(listEquals(oneSearch[[1]][[4]], list(keyword="Magic", offset=20)))
stopifnot(listEquals(oneSearch[[1]][[5]], list(keyword="Word", offset=26)))

# search a list of lists
# * sublists are accessed by index
# * texts are accessed by index
# * non-matched texts are kept (to preserve index order)
listSearch = AhoCorasickSearchList(keywords, list(c("What in", "the world"), c("is"), "secret about", "the Magic Word?"))
stopifnot(listEquals(listSearch[[1]][[1]], list()))
stopifnot(listEquals(listSearch[[1]][[2]][[1]], list(keyword="the", offset=1)))
stopifnot(listEquals(listSearch[[2]][[1]][[1]], list(keyword="is", offset=1)))
stopifnot(listEquals(listSearch[[3]], list()))
stopifnot(listEquals(listSearch[[4]][[1]][[1]], list(keyword="the", offset=1)))
stopifnot(listEquals(listSearch[[4]][[1]][[2]], list(keyword="Magic", offset=5)))
stopifnot(listEquals(listSearch[[4]][[1]][[3]], list(keyword="Word", offset=11)))

# named search of a list of lists
# * sublists are accessed by name
# * matched texts are accessed by name
# * non-matched texts are dropped
namedSearch = AhoCorasickSearchList(keywords, list(subject=c(phrase1="What in", phrase2="the world"),
                                                   verb=c(phrase1="is"),
                                                   predicate1=c(phrase1="secret about"),
                                                   predicate2=c(phrase1="the Magic Word?")))
stopifnot(listEquals(namedSearch$subject$phrase2[[1]], list(keyword="the", offset=1)))
stopifnot(listEquals(namedSearch$verb$phrase1[[1]], list(keyword="is", offset=1)))
stopifnot(listEquals(namedSearch$predicate1, list()))
stopifnot(listEquals(namedSearch$predicate2$phrase1[[1]], list(keyword="the", offset=1)))
stopifnot(listEquals(namedSearch$predicate2$phrase1[[2]], list(keyword="Magic", offset=5)))
stopifnot(listEquals(namedSearch$predicate2$phrase1[[3]], list(keyword="Word", offset=11)))

# named search of multiple texts in a single list with keyword grouping and aminoacid alphabet
# * all matches to a keyword are accessed by name
# * non-matched keywords are dropped
proteins = c(protein1="PEPTIDEPEPTIDEDADADARARARARAKEKEKEKEPEPTIDE",
             protein2="DERPADERPAPEWPEWPEEPEERAWRAWWARRAGTAGPEPTIDEKESEQUENCE")
peptides = c("PEPTIDE", "DERPA", "SEQUENCE", "KEKE", "PEPPIE")
peptideSearch = AhoCorasickSearch(peptides, proteins, alphabet="aminoacid", groupByKeyword=T)
stopifnot(listEquals(peptideSearch$PEPTIDE, list(list(keyword="protein1", offset=1),
                                                 list(keyword="protein1", offset=8),
                                                 list(keyword="protein1", offset=37),
                                                 list(keyword="protein2", offset=38))))
stopifnot(listEquals(peptideSearch$DERPA, list(list(keyword="protein2", offset=1),
                                               list(keyword="protein2", offset=6))))
stopifnot(listEquals(peptideSearch$SEQUENCE, list(list(keyword="protein2", offset=47))))
stopifnot(listEquals(peptideSearch$KEKE, list(list(keyword="protein1", offset=29),
                                              list(keyword="protein1", offset=31),
                                              list(keyword="protein1", offset=33))))
stopifnot(listEquals(peptideSearch$PEPPIE, NULL))

# grouping by keyword without text names: offsets are given without reference to the text
names(proteins) = NULL
peptideSearch = AhoCorasickSearch(peptides, proteins, groupByKeyword=T)
stopifnot(listEquals(peptideSearch$PEPTIDE, list(1, 8, 37, 38)))
stopifnot(listEquals(peptideSearch$DERPA, list(1, 6)))
stopifnot(listEquals(peptideSearch$SEQUENCE, list(47)))
stopifnot(listEquals(peptideSearch$KEKE, list(29, 31, 33)))