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

Question

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

6

我有一个字符串列表（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行和16列。我想要计算我的k-mers与列表成员之间的匹配数。

我的期望结果：df：

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

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

- Cina

我的数据库非常庞大，因此效率在这里也非常重要。谢谢。 - Cina

5个回答

7

如果您想要更快的速度，显然的解决方案是使用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

比之前快了165倍 :)

- bartektartanus

6

我们最近在Bioconductor 3.0发布中发布了我们的“kebabs”软件包。尽管该软件包旨在为分类、回归和基于相似性的聚类等任务提供序列核函数，但该软件包还包括高效计算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-mer（正如您所要求的）。即使对于k = 13和具有成千上万个碱基的序列，计算也只需花费几分之一秒钟（我们5年前的Dell服务器上为19毫秒）。上述函数也适用于DNAStringSets，但在这种情况下，您应该删除drop（）以获取k-mer频率矩阵。默认情况下，矩阵是稀疏的（类'dgRMatrix'），但您也可以强制将结果格式化为标准密集矩阵格式（但仍然省略在任何序列中都不出现的k-mer）。

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

k-mer的长度可能取决于您的系统。在我们的系统上，DNA序列的限制似乎为k=22。RNA和氨基酸序列同样适用。然而，对于后者，由于特征空间显然更大，因此在k方面的限制要显著较低。

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

我希望这能帮到您。如果您还有任何问题，请告诉我。

最好的问候， Ulrich

- UBod

谢谢你明确的回答。关于spectrumKernel(k=13, normalized=TRUE)中的归一化，我有一个问题：它是按行归一化还是按列归一化？ - Cina

核规范化被定义为 k'(x, y) = k(x, y) / sqrt(k(x, x) * k(y, y))。在显式表示的情况下，这相当于通过将每行除以其欧几里德范数进行逐行归一化。如果需要列向归一化，则计算不带规范化的显式表示（如上面的答案）并对矩阵应用scale()。 - UBod

5

我的回答不像@bartektartanus那么快。但是，它也非常快，而且我写了代码... :D

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

不需要安装未实现版本的stri_count_fixed
对于大的k-mer，可能stringi包会变得非常慢，因为它必须为模式生成所有可能的组合，然后检查它们在数据中的存在并计算出它们出现的次数。
它还可以对长的单个和多个序列进行快速处理，并具有相同的输出。
您可以为k设置一个值，而不是创建一个模式字符串。
如果您在大序列中使用k大于12运行oligonucleotideFrequency函数，则由于超额内存使用而导致函数冻结并重新启动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)
}

我只使用Biostrings包来计算碱基数量...你可以使用其他选项，如stringi来进行计数... 如果你删除k_mers lapply和return(k_mers)下面的所有代码，它将返回一个列表...包含所有k-mer及其对应的重复向量

`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

- Rodrigo de Alexandre

2

另一种方法是这样的：

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

- user4415867

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- akrun · Accepted Answer

也许这会有所帮助。

 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

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

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

我的代码

sequence在这里是一个长度为1000bp的序列

`sequence`在这里是一个长度为1000bp的序列