R中的匹配和计数字符串（DNA的k-mer）答案

【问题标题】：matching and counting strings (k-mer of DNA) in RR中的匹配和计数字符串（DNA的k-mer）
【发布时间】：2014-10-28 04:07:21
【问题描述】：

我有一个字符串列表（DNA 序列），包括 A、T、C、G。我想找到所有匹配项并插入表中，其列都是这些 DNA 字母表的所有可能组合（4 ^ k；“k”是每个匹配项的长度 - K-mer - 并且必须由用户指定）并且行表示在列表中按顺序匹配。

假设我的列表包括 5 个成员：

DNAlst<-list("CAAACTGATTTT","GATGAAAGTAAAATACCG","ATTATGC","TGGA","CGCGCATCAA")

我想设置k=2 (2-mer) 所以4^2=16 组合可用，包括AA,AT,AC,AG,TA,TT,...

所以我的表将有5 rows 和16 columns。我想计算我的 k-mers 和列表成员之间的匹配数。

我想要的结果：df:

lstMemb AA AT AC AG TA TT TC ...
  1     2  1  1  0  0  3  0
  2       ...
  3
  4
  5

你能帮我在 R 中实现这个吗？

【问题讨论】：

我的数据库很大，所以效率在这里也很重要。谢谢。

标签： regex r string dna-sequence

【解决方案1】：

如果您正在寻找速度，显而易见的解决方案是 stringi 包。有stri_count_fixed 函数用于计数模式。现在，检查代码和基准测试！

DNAlst<-list("CAAACTGATTTT","GATGAAAGTAAAATACCG","ATTATGC","TGGA","CGCGCATCAA")
dna <- stri_paste(rep(c("A","C","G","T"),each=4),c("A","C","G","T"))
result <- t(sapply(DNAlst, stri_count_fixed,pattern=dna,overlap=TRUE))
colnames(result) <- dna
result
     AA AC AG AT CA CC CG CT GA GC GG GT TA TC TG TT
[1,]  2  1  0  1  1  0  0  1  1  0  0  0  0  0  1  3
[2,]  5  1  1  2  0  1  1  0  2  0  0  1  2  0  1  0
[3,]  0  0  0  2  0  0  0  0  0  1  0  0  1  0  1  1
[4,]  0  0  0  0  0  0  0  0  1  0  1  0  0  0  1  0
[5,]  1  0  0  1  2  0  2  0  0  2  0  0  0  1  0  0



fstri <- function(x){
    t(sapply(x, stri_count_fixed,dna,T))
}
fbio <- function(x){
    t(sapply(x, function(x){x1 <-  DNAString(x); oligonucleotideFrequency(x1,2)}))
}

all(fstri(DNAlst)==fbio(DNAlst)) #results are the same
[1] TRUE

longDNA <- sample(DNAlst,100,T)
microbenchmark(fstri(longDNA),fbio(longDNA))
Unit: microseconds
           expr        min         lq        mean     median         uq        max neval
 fstri(longDNA)    689.378    738.184    825.3014    766.862    793.134   6027.039   100
  fbio(longDNA) 118371.825 125552.401 129543.6585 127245.489 129165.711 359335.294   100
127245.489/766.862
## [1] 165.9301

Ca 速度快 165 倍 :)

【讨论】：

【解决方案2】：

这可能有帮助

 source("http://bioconductor.org/biocLite.R")
 biocLite("Biostrings")
 library(Biostrings)
 t(sapply(DNAlst, function(x){x1 <-  DNAString(x)
                   oligonucleotideFrequency(x1,2)}))
  #     AA AC AG AT CA CC CG CT GA GC GG GT TA TC TG TT
  #[1,]  2  1  0  1  1  0  0  1  1  0  0  0  0  0  1  3
  #[2,]  5  1  1  2  0  1  1  0  2  0  0  1  2  0  1  0
  #[3,]  0  0  0  2  0  0  0  0  0  1  0  0  1  0  1  1
  #[4,]  0  0  0  0  0  0  0  0  1  0  1  0  0  0  1  0
  #[5,]  1  0  0  1  2  0  2  0  0  2  0  0  0  1  0  0

或者按照@Arun 的建议，首先将list 转换为vector

   oligonucleotideFrequency(DNAStringSet(unlist(DNAlst)), 2L)
   #     AA AC AG AT CA CC CG CT GA GC GG GT TA TC TG TT
   #[1,]  2  1  0  1  1  0  0  1  1  0  0  0  0  0  1  3
   #[2,]  5  1  1  2  0  1  1  0  2  0  0  1  2  0  1  0
   #[3,]  0  0  0  2  0  0  0  0  0  1  0  0  1  0  1  1
   #[4,]  0  0  0  0  0  0  0  0  1  0  1  0  0  0  1  0
   #[5,]  1  0  0  1  2  0  2  0  0  2  0  0  0  1  0  0

【讨论】：

你应该这样做：oligonucleotideFrequency(DNAStringSet(x), 2L)，其中x 是unlist(DNAlist)！！
@Arun 谢谢，这样更好。

【解决方案3】：

我的回答没有@bartektartanus 快。但是，它也很快，我编写了代码......：D

与其他代码相比，我的代码的优点是：

不需要安装stri_count_fixed的未实现版本
可能stringi 包对于大型 k-mer 来说会变得非常慢，因为它必须为模式生成所有可能的组合，然后，检查它们在数据中的存在并计算它出现的次数。
它也适用于长 single 和 multiple 序列同样的输出速度非常快。
您可以为k 设置一个值，而不是创建模式字符串。
如果运行oligonucleotideFrequency 时k 大于12 在一个大的序列中，函数会因内存使用过多而冻结，并且 R 重新启动，而我的函数运行得非常快。

我的代码

sequence_kmers <- function(sequence, k){
    k_mers <- lapply(sequence,function(x){
        seq_loop_size <- length(DNAString(x))-k+1

        kmers <- sapply(1:seq_loop_size, function(z){
            y <- z + k -1
            kmer <- substr(x=x, start=z, stop=y)
            return(kmer)
        })
        return(kmers)
    })

    uniq <- unique(unlist(k_mers))
    ind <- t(sapply(k_mers, function(x){
        tabulate(match(x, uniq), length(uniq))
    }))
    colnames(ind) <- uniq

    return(ind)
}

我只使用Biostringspackage 来计算碱基...您可以使用其他选项，例如stringi 来计算... 如果您删除 k_mers lapply 和 return(k_mers) 下面的所有代码，它只返回列表...所有 k-mers 以及相应的重复向量

`sequence`这里是1000bp的序列

#same output for 1 or multiple sequences
> sequence_kmers(sequence,4)[,1:10]
GTCT TCTG CTGA TGAA GAAC AACG ACGC CGCG GCGA CGAG 
   4    4    3    4    4    8    6    4    5    5 
> sequence_kmers(c(sequence,sequence),4)[,1:10]
     GTCT TCTG CTGA TGAA GAAC AACG ACGC CGCG GCGA CGAG
[1,]    4    4    3    4    4    8    6    4    5    5
[2,]    4    4    3    4    4    8    6    4    5    5

使用我的函数完成的测试：

#super fast for 1 sequence
> system.time({sequence_kmers(sequence,13)})
  usuário   sistema decorrido 
     0.08      0.00      0.08 

#works fast for 1 sequence or 50 sequences of 1000bps
> system.time({sequence_kmers(rep(sequence,50),4)})
     user    system   elapsed
     3.61      0.00      3.61 

#same speed for 3-mers or 13-mers
> system.time({sequence_kmers(rep(sequence,50),13)})
     user    system   elapsed
     3.63      0.00      3.62

使用Biostrings 完成的测试：

#Slow 1 sequence 12-mers
> system.time({oligonucleotideFrequency(DNAString(sequence),12)})
     user    system   elapsed 
   150.11      1.14    151.37 

#Biostrings package freezes for a single sequence of 13-mers
> system.time({oligonucleotideFrequency(sequence,13)})  
freezes, used all my 8gb RAM

【讨论】：

【解决方案4】：

我们最近发布了作为 Bioconductor 的一部分的“烤肉串”包 3.0 发布。虽然这个包旨在提供序列内核用于分类、回归和其他任务，例如基于相似性的聚类，该软件包还包括有效计算 k-mer 频率的功能：

#installing kebabs:
#source("http://bioconductor.org/biocLite.R")
#biocLite(c("kebabs", "Biostrings"))
library(kebabs)

s1 <- DNAString("ATCGATCGATCGATCGATCGATCGACTGACTAGCTAGCTACGATCGACTG")
s1
s2 <- DNAString(paste0(rep(s1, 200), collate=""))
s2

sk13 <- spectrumKernel(k=13, normalized=FALSE)
system.time(kmerFreq <- drop(getExRep(s1, sk13)))
kmerFreq
system.time(kmerFreq <- drop(getExRep(s2, sk13)))
kmerFreq

所以你看到 k-mer 频率是作为显式获得的 k=13 的标准（非归一化）谱内核的特征向量。这个函数是用高效的 C++ 代码实现的前缀树，只考虑实际出现在序列（根据您的要求）。即使对于 k=13 和一个序列，您也会看到有数万个碱基，计算只需要几分之一一秒钟（在我们使用了 5 年的戴尔服务器上为 19 毫秒）。上述函数也适用于 DNAStringSets，但在这种情况下，您应该删除 drop() 得到一个 k-mer 频率矩阵。矩阵是默认的 sparse（类 'dgRMatrix'），但您也可以强制结果在标准密集矩阵格式（但是，仍然省略不发生在任何序列中）：

sv <- c(DNAStringSet(s1), DNAStringSet(s2))
system.time(kmerFreq <- getExRep(sv, sk13))
kmerFreq
system.time(kmerFreq <- getExRep(sv, sk13, sparse=FALSE))
kmerFreq

k-mers 的长度可能取决于您的系统。在我们的系统上， DNA 序列的限制似乎是 k=22。同样适用于 RNA 和氨基酸序列。然而，对于后者，在 k 方面的限制显着降低，因为特征空间明显很多相同的 k 更大。

#for the kebabs documentation please see:
browseVignettes("kebabs")

我希望这会有所帮助。如果您还有任何问题，请告诉我。

最好的问候，乌尔里希

【讨论】：

感谢您的明确回答。 spectrumKernel(k=13, normalized=TRUE) 中关于标准化的一个问题；它是按行还是按列标准化？
核归一化定义为 k'(x, y) = k(x, y) / sqrt(k(x, x) * k(y, y))。在显式表示的情况下，这对应于通过将每一行除以其欧几里得范数来进行行标准化。如果您需要按列归一化，请计算没有归一化的显式表示（如我上面的答案）并将 scale() 应用于矩阵。

【解决方案5】：

另一种方法：

DNAlst<-list("CAAACTGATTTT","GATGAAAGTAAAATACCG","ATTATGC","TGGA","CGCGCATCAA","ACACACACACCA")
len <- 4
stri_sub_fun <- function(x) table(stri_sub(x,1:(stri_length(x)-len+1),length = len))
sapply(DNAlst, stri_sub_fun)
[[1]]

AAAC AACT ACTG ATTT CAAA CTGA GATT TGAT TTTT 
   1    1    1    1    1    1    1    1    1 

[[2]]

AAAA AAAG AAAT AAGT AATA ACCG AGTA ATAC ATGA GAAA GATG GTAA TAAA TACC TGAA 
   1    1    1    1    1    1    1    1    1    1    1    1    1    1    1 

[[3]]

ATGC ATTA TATG TTAT 
   1    1    1    1 

[[4]]

TGGA 
   1 

[[5]]

ATCA CATC CGCA CGCG GCAT GCGC TCAA 
   1    1    1    1    1    1    1 

[[6]]

ACAC ACCA CACA CACC 
   4    1    3    1

【讨论】：

与其他代码相比，我的代码的优点是：

我的代码

sequence这里是1000bp的序列

`sequence`这里是1000bp的序列