【问题标题】:match with an interval and extract values between two matrix R匹配区间并提取两个矩阵 R 之间的值
【发布时间】:2016-07-17 23:16:19
【问题描述】:

我在一个列表中有 n 个矩阵和一个附加矩阵,其中包含我想在矩阵列表中找到的值。

要获取矩阵列表,我使用以下代码:

setwd("C:\\~\\Documents\\R") 


import.multiple.txt.files<-function(pattern=".txt",header=T)
{
list.1<-list.files(pattern=".txt")
list.2<-list()
for (i in 1:length(list.1))
{
list.2[[i]]<-read.delim(list.1[i])
}
names(list.2)<-list.1
list.2

}


txt.import.matrix<-cbind(txt.import)

我的列表如下所示:(我只展示了一个 n=2 的示例)。每个数组的行数不同(这里我只取 5 行和 6 行来简化,但我的真实数据超过 500)。

txt.import.matrix[1]

    [[1]]
     X.     RT.     Area.  m.z.      
1     1     1.01   2820.1  358.9777  
2     2     1.03   9571.8  368.4238  
3     3     2.03   6674.0  284.3294  
4     4     2.03   5856.3  922.0094  
5     5     3.03   27814.6 261.1299  


txt.import.matrix[2]

    [[2]]
     X.     RT.     Area.  m.z.      
1     1     1.01    7820.1 358.9777  
2     2     1.06    8271.8 368.4238  
3     3     2.03   12674.0 284.3294  
4     4     2.03    5856.6 922.0096  
5     5     2.03   17814.6 261.1299
6     6     3.65    5546.5 528.6475  

我想在矩阵列表中找到另一个值数组。该数组是通过将列表中的所有数组组合成一个数组并删除重复项获得的。

reduced.list.pre.filtering

     RT.   m.z.
1    1.01  358.9777
2    1.07  368.4238
3    2.05  284.3295
4    2.03  922.0092
5    3.03  261.1299
6    3.56  869.4558

我想获得一个新矩阵,其中写入了列表中所有矩阵匹配的RT. ± 0.02m.z. ± 0.0002Area. 结果。输出可能是这样的。

     RT.   m.z.        Area.[1]      Area.[2]
1    1.01  358.9777    2820.1        7820.1
2    1.07  368.4238                  8271.8      
3    2.05  284.3295    6674.0        12674.0
4    2.03  922.0092    5856.3             
5    3.03  261.1299    27814.6            
6    3.65  528.6475    

我只知道如何只匹配一个数组中的一个精确值。这里的困难是在数组列表中找到值,并且需要找到值±一个区间。如果您有任何建议,我将不胜感激。

【问题讨论】:

    标签: r find match intervals


    【解决方案1】:

    使用来自 data.table 的当前开发版本 v1.9.7 的 non-equi 联接(请参阅 installation instructions),它允许为联接提供非相等条件:

    require(data.table) # v1.9.7
    names(ll) = c("Area1", "Area2")
    A = rbindlist(lapply(ll, as.data.table), idcol = "id")           ## (1)
    
    B = as.data.table(mat)
    B[, c("RT.minus", "RT.plus") := .(RT.-0.02, RT.+0.02)]
    B[, c("m.z.minus", "m.z.plus") := .(m.z.-0.0002, m.z.+0.0002)]   ## (2)
    
    ans = A[B, .(id, X., RT. = i.RT., m.z. = i.m.z., Area.), 
               on = .(RT. >= RT.minus, RT. <= RT.plus, 
                      m.z. >= m.z.minus, m.z. <= m.z.plus)]          ## (3)
    
    dcast(ans, RT. + m.z. ~ id)                                      ## (4)
    # or dcast(ans, RT. + m.z. ~ id, fill = 0)
    #     RT.     m.z.   Area1   Area2
    # 1: 1.01 358.9777  2820.1  7820.1
    # 2: 1.07 368.4238      NA  8271.8
    # 3: 2.03 922.0092  5856.3      NA
    # 4: 2.05 284.3295  6674.0 12674.0
    # 5: 3.03 261.1299 27814.6      NA
    

    [1] 命名矩阵列表(此处称为ll)并使用lapply() 将它们中的每一个转换为data.table,并使用rbindlist 将它们逐行绑定,并将名称添加为额外的列 (idcol)。叫它A

    [2] 将第二个矩阵(此处称为mat)也转换为data.table。添加与您要搜索的范围/间隔相对应的其他列(因为on= 参数,正如我们将在下一步中看到的那样,无法处理表达式)。叫它B

    [3] 执行条件连接/子集。对于B 中的每一行,在A 中找到与提供给on= 参数的条件相对应的匹配行,并为那些匹配的行索引提取列id, X., R.T. and m.z.

    [4] 最好保留在 [3]。但是,如果您希望它如您的答案所示,我们可以将其重塑为宽格式。 fill = 0 会将结果中的 NAs 替换为 0

    【讨论】:

      【解决方案2】:

      这是 Arun 使用 data.table 的相当优雅的答案的另一种方法。我决定发布它是因为它包含两个额外的方面,这些方面是您问题中的重要考虑因素:

      1. 浮点比较:比较以查看浮点值是否在区间内,需要在计算区间时考虑舍入误差。这是比较实数的浮点表示的一般问题。参见R上下文中的thisthis。下面在函数in.interval中实现了这种比较。

      2. 多个匹配:如果间隔重叠,您的间隔匹配标准可能会导致多个匹配。以下假设您只需要第一个匹配项(相对于每个txt.import.matrix 矩阵的增加行)。这在函数match.interval 中实现,并在后面的注释中解释。如果您想获得符合标准的区域的平均值,则需要其他逻辑。

      为了从txt.import.matrix 中为矩阵reduced.list.pre.filtering 中的每一行找到矩阵中的匹配行,以下代码将比较函数在@987654329 之间的所有枚举行对的空间上的应用向量化@ 和来自txt.import.matrix 的矩阵。对于此应用程序的功能,这与 Arun 使用 data.tablenon-equi 连接的解决方案相同;但是,non-equi 连接功能更通用,data.table 实现很可能针对此应用程序的内存使用和速度进行了更好的优化。

      in.interval <- function(x, center, deviation, tol = .Machine$double.eps^0.5) {
        return (abs(x-center) <= (deviation + tol))
      }
      
      match.interval <- function(r, t) {
        r.rt <- rep(r[,1], each=nrow(t))
        t.rt <- rep(t[,2], times=nrow(r))
        r.mz <- rep(r[,2], each=nrow(t))
        t.mz <- rep(t[,4], times=nrow(r))                                       ## 1.
      
        ind <- which(in.interval(r.rt, t.rt, 0.02) & 
                     in.interval(r.mz, t.mz, 0.0002))
        r.ind <- floor((ind - 1)/nrow(t)) + 1                                   ## 2.
      
        dup <- duplicated(r.ind)
        r.ind <- r.ind[!dup]
        t.ind <- ind[!dup] - (r.ind - 1)*nrow(t)                                ## 3.
        return(cbind(r.ind,t.ind))                       
      }
      
      get.area.matched <- function(r, t) {
        match.ind <- match.interval(r, t)
        area <- rep(NA,nrow(r))
        area[match.ind[,1]] <- t[match.ind[,2], 3]                              ## 4.
        return(area)
      }
      
      res <- cbind(reduced.list.pre.filtering,
                   do.call(cbind,lapply(txt.import.matrix, 
                                        get.area.matched, 
                                        r=reduced.list.pre.filtering)))         ## 5.
      colnames(res) <- c(colnames(reduced.list.pre.filtering), 
                         sapply(seq_len(length(txt.import.matrix)), 
                                function(i) {return(paste0("Area.[",i,"]"))}))  ## 6.
      print(res)
      ##      RT.     m.z. Area.[1] Area.[2]
      ##[1,] 1.01 358.9777   2820.1   7820.1
      ##[2,] 1.07 368.4238       NA   8271.8
      ##[3,] 2.05 284.3295   6674.0  12674.0
      ##[4,] 2.03 922.0092   5856.3       NA
      ##[5,] 3.03 261.1299  27814.6       NA
      ##[6,] 3.56 869.4558       NA       NA
      

      注意事项:

      1. 这部分构造数据,以便在reduced.list.pre.filtering 和来自txt.import.matrix 的矩阵之间的所有枚举行对的空间上实现比较函数应用的矢量化。要构造的数据是四个数组,它们是在比较标准中使用的两列的复制(或扩展),reduced.list.pre.filtering 在来自txt.import.matrix 的每个矩阵的行维度中以及两列的复制,用于比较标准,来自txt.import.matrix 的每个矩阵在reduced.list.pre.filtering 的行维度中。这里,术语数组是指二维矩阵或一维向量。得到的四个数组是:

        • r.rtreduced.list.pre.filteringRT. 列(即r[,1])在t 的行维度中的复制
        • t.rt 是矩阵的RT. 列从txt.import.matrix(即t[,2])在r 的行维度中的复制
        • r.mzreduced.list.pre.filteringm.z. 列(即r[,2])在t 的行维度中的复制
        • t.mz 是矩阵的m.z. 列从txt.import.matrix(即t[,4])在r 的行维度中的复制

        重要的是每个数组的索引以相同的方式枚举rt 中的所有行对。具体来说,将这些数组视为大小为M by N 的二维矩阵,其中M=nrow(t)N=nrow(r),行索引对应于t 的行,列索引对应于@987654369 的行@。因此,i-th 行和j-th 列(四个数组中的每一个)的数组值(所有四个数组)是j-th 之间的比较标准中使用的值r 的行和ti-th 行。此复制过程的实现使用 R 函数rep。例如,在计算r.rt 时,使用repeach=M,其作用是将其数组输入r[,1] 视为行向量,并将该行复制M 次以形成M 行。结果是,与r 中的一行相对应的每一列都具有来自r 相应行的RT. 值,并且该值对于r.rt 的所有行(该列的)都相同,每一个都对应t中的一行。这意味着在将r 中的该行与t 中的任何行进行比较时,将使用r 中该行的RT. 的值。相反,在计算t.rt 时,使用reptimes=N,其作用是将其数组输入视为列向量并将该列复制N 次以形成N 列。结果是,t.rt 中的每一行(对应于t 中的一行)具有来自t 相应行的RT. 值,并且该值对于(该行的)所有列都是相同的t.rt 中的每一行对应r 中的一行。这意味着在将t 中的该行与r 中的任何行进行比较时,将使用t 中该行的RT. 的值。同样,r.mzt.mz 的计算分别使用来自 rtm.z. 列。

      2. 这将执行矢量化比较,生成 MN 的逻辑矩阵,其中 i-th 行和 j-th 列是 TRUE,如果 j-th r 的行与i 的第t 行匹配条件,否则FALSEwhich() 的输出是此逻辑比较结果矩阵的数组索引数组,其中元素为TRUE。我们希望将这些数组索引转换为比较结果矩阵的行索引和列索引,以引用rt 的行。下一行从数组索引中提取列索引。请注意,变量名称是r.ind,表示这些对应于r 的行。我们首先提取它,因为它对于检测r 中的一行的多个匹配项很重要。

      3. 这部分为r 中的给定行处理t 中可能的多个匹配项。多个匹配项将在r.ind 中显示为重复值。如上所述,这里的逻辑只保留t 中增加行的第一个匹配项。函数duplicated 返回数组中所有重复值的索引。因此删除这些元素将做我们想要的。代码首先从r.ind 中删除它们,然后从ind 中删除它们,最后使用修剪后的ind 和@987654438 计算与t 的行相对应的比较结果矩阵的列索引@。 match.interval 返回的是一个矩阵,其行是匹配的行索引对,其第一列是r 的行索引,第二列是t 的行索引。

      4. get.area.matched 函数仅使用来自match.ind 的结果从t 中提取Area 以获取所有匹配项。请注意,返回的结果是一个(列)向量,其长度等于r 中的行数,并初始化为NA。这样,r 中与 t 不匹配的行将返回 AreaNA

      5. 这使用lapply 在列表txt.import.matrix 上应用函数get.area.matched,并将返回的匹配Area 结果附加到reduced.list.pre.filtering 作为列向量。同样,相应的列名也被附加并设置在结果res中。

      编辑:使用foreach 包的替代实现

      事后看来,更好的实现是使用foreach 包来对比较进行矢量化处理。在此实现中,需要 foreachmagrittr

      require("magrittr")  ## for %>%
      require("foreach")
      

      然后是match.interval中的代码,用于向量化比较

      r.rt <- rep(r[,1], each=nrow(t))
      t.rt <- rep(t[,2], times=nrow(r))
      r.mz <- rep(r[,2], each=nrow(t))
      t.mz <- rep(t[,4], times=nrow(r))                       # 1.
      
      ind <- which(in.interval(r.rt, t.rt, 0.02) & 
                   in.interval(r.mz, t.mz, 0.0002))
      

      可以替换为

      ind <- foreach(r.row = 1:nrow(r), .combine=cbind) %:% 
               foreach(t.row = 1:nrow(t)) %do% 
                 match.criterion(r.row, t.row, r, t) %>% 
                   as.logical(.) %>% which(.)
      

      match.criterion 定义为

      match.criterion <- function(r.row, t.row, r, t) {
        return(in.interval(r[r.row,1], t[t.row,2], 0.02) & 
               in.interval(r[r.row,2], t[t.row,4], 0.0002))
      }
      

      这更容易解析并反映正在执行的内容。请注意,嵌套 foreachcbind 组合返回的内容再次是一个逻辑矩阵。最后,get.area.matched 函数对列表txt.import.matrix 的应用也可以使用foreach 执行:

      res <- foreach(i = 1:length(txt.import.matrix), .combine=cbind) %do% 
               get.area.matched(reduced.list.pre.filtering, txt.import.matrix[[i]]) %>%
                 cbind(reduced.list.pre.filtering,.)
      

      使用foreach的完整代码如下:

      require("magrittr")
      require("foreach")
      
      in.interval <- function(x, center, deviation, tol = .Machine$double.eps^0.5) {
        return (abs(x-center) <= (deviation + tol))
      }
      
      match.criterion <- function(r.row, t.row, r, t) {
        return(in.interval(r[r.row,1], t[t.row,2], 0.02) & 
           in.interval(r[r.row,2], t[t.row,4], 0.0002))
      }
      
      match.interval <- function(r, t) {
        ind <- foreach(r.row = 1:nrow(r), .combine=cbind) %:% 
             foreach(t.row = 1:nrow(t)) %do% 
           match.criterion(r.row, t.row, r, t) %>% 
             as.logical(.) %>% which(.)
        # which returns 1-D indices (row-major),
        # convert these to 2-D indices in (row,col)
        r.ind <- floor((ind - 1)/nrow(t)) + 1                   ## 2.
        # detect duplicates in r.ind and remove them from ind
        dup <- duplicated(r.ind)
        r.ind <- r.ind[!dup]
        t.ind <- ind[!dup] - (r.ind - 1)*nrow(t)                ## 3.
      
        return(cbind(r.ind,t.ind))                       
      }
      
      get.area.matched <- function(r, t) {
        match.ind <- match.interval(r, t)
        area <- rep(NA,nrow(r))
        area[match.ind[,1]] <- t[match.ind[,2], 3]
        return(area)
      }
      
      res <- foreach(i = 1:length(txt.import.matrix), .combine=cbind) %do% 
           get.area.matched(reduced.list.pre.filtering, txt.import.matrix[[i]]) %>%
             cbind(reduced.list.pre.filtering,.)
      
      colnames(res) <- c(colnames(reduced.list.pre.filtering), 
                 sapply(seq_len(length(txt.import.matrix)), 
                    function(i) {return(paste0("Area.[",i,"]"))}))
      

      希望这会有所帮助。

      【讨论】:

      • 谢谢爱超!确实,您的代码对我来说有点复杂,但尝试理解它是一个很好的挑战。它运行得很好,我得到了我想要的输出。是的,这很有帮助!
      • @Vanbell:谢谢。我还在编辑帖子。明天再回来看看,希望它会更清楚。当然,如果您有任何问题,请发表评论。同样,我认为我应该发帖指出您的问题中我认为很重要的两个可能方面。
      • 好的,我在等你的帖子。我的事情,我会问你一些问题,因为它太容易复制和过去;)
      • 我仔细检查了结果,但没有得到预期的结果。我的数据有点复杂。事实上,对于构成txt.import.matrix 的所有矩阵,我没有相同的行数。我想我必须更改这部分代码sapply(seq_len(length(txt.import.matrix)), function(i) {return(paste0("Area[",i,"]"))})) 我试过function(i) { sapply(seq_len(length(txt.import.matrix[i])), return(paste0("Area[",i,"]"))})) 但它不起作用。我不明白,因为它应该考虑到列表中每个矩阵的长度。
      • @Vanbell:你能用更多说明问题的数据更新你的问题吗?请使用dput(),以便我可以快速重构它。与此同时,我已经更新了这篇文章,这样你可能会更清楚。我还打算再次编辑以使用 foreach 包添加另一个实现。
      【解决方案3】:

      这是一种快速粗略的方法,如果我得到您想要做的事情,可能会有所帮助。

      从两个矩阵的每个变量中取消列出值

      areas <- unlist(lapply(txt.import.matrix, function(x) x$Area.))
      rts <- unlist(lapply(txt.import.matrix, function(x) x$RT.))
      mzs <- unlist(lapply(txt.import.matrix, function(x) x$m.z.))
      

      查找那些 RT 和 m.z 值的索引。最接近第三个矩阵/df 中的值:

       rtmins <- lapply(reduced.list.pre.filtering$RT., function(x) which(abs(rts-x)==min(abs(rts-x))))
      mzmins <- lapply(reduced.list.pre.filtering$m.z., function(x) which(abs(mzs-x)==min(abs(mzs-x))))
      

      使用purrr 快速计算两者中的索引(即每个索引的最小差异):

      inboth <- purrr::map2(rtmins,mzmins,intersect)
      

      获取对应的面积值:

      vals<-lapply(inboth, function(x) areas[x])
      

      使用reshape2转成宽格式:

      vals2 <- reshape2::melt(vals)
      vals2$number <- ave(vals2$L1, vals2$L1, FUN = seq_along)
      vals.wide <-reshape2::dcast(vals2, L1 ~ number, value.var="value")
      
      cbind(reduced.list.pre.filtering, vals.wide)
      
      #   RT.     m.z. L1       1       2
      #1 1.01 358.9777  1  2820.1  7820.1
      #2 1.07 368.4238  2  8271.8      NA
      #3 2.05 284.3295  3  6674.0 12674.0
      #4 2.03 922.0092  4  5856.3      NA
      #5 3.03 261.1299  5 27814.6      NA
      

      这可能会给你一些想法。可以很容易地适应检查共享最小值是否超过 +/- 一个值。

      【讨论】:

      • 谢谢,申请成功后第一时间通知您。
      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2021-10-09
      • 1970-01-01
      相关资源
      最近更新 更多