这是 Arun 使用 data.table 的相当优雅的答案的另一种方法。我决定发布它是因为它包含两个额外的方面,这些方面是您问题中的重要考虑因素:
浮点比较:比较以查看浮点值是否在区间内,需要在计算区间时考虑舍入误差。这是比较实数的浮点表示的一般问题。参见R上下文中的this和this。下面在函数in.interval中实现了这种比较。
多个匹配:如果间隔重叠,您的间隔匹配标准可能会导致多个匹配。以下假设您只需要第一个匹配项(相对于每个txt.import.matrix 矩阵的增加行)。这在函数match.interval 中实现,并在后面的注释中解释。如果您想获得符合标准的区域的平均值,则需要其他逻辑。
为了从txt.import.matrix 中为矩阵reduced.list.pre.filtering 中的每一行找到矩阵中的匹配行,以下代码将比较函数在@987654329 之间的所有枚举行对的空间上的应用向量化@ 和来自txt.import.matrix 的矩阵。对于此应用程序的功能,这与 Arun 使用 data.table 的 non-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
注意事项:
-
这部分构造数据,以便在reduced.list.pre.filtering 和来自txt.import.matrix 的矩阵之间的所有枚举行对的空间上实现比较函数应用的矢量化。要构造的数据是四个数组,它们是在比较标准中使用的两列的复制(或扩展),reduced.list.pre.filtering 在来自txt.import.matrix 的每个矩阵的行维度中以及两列的复制,用于比较标准,来自txt.import.matrix 的每个矩阵在reduced.list.pre.filtering 的行维度中。这里,术语数组是指二维矩阵或一维向量。得到的四个数组是:
-
r.rt 是reduced.list.pre.filtering 的RT. 列(即r[,1])在t 的行维度中的复制
-
t.rt 是矩阵的RT. 列从txt.import.matrix(即t[,2])在r 的行维度中的复制
-
r.mz 是reduced.list.pre.filtering 的m.z. 列(即r[,2])在t 的行维度中的复制
-
t.mz 是矩阵的m.z. 列从txt.import.matrix(即t[,4])在r 的行维度中的复制
重要的是每个数组的索引以相同的方式枚举r 和t 中的所有行对。具体来说,将这些数组视为大小为M by N 的二维矩阵,其中M=nrow(t) 和N=nrow(r),行索引对应于t 的行,列索引对应于@987654369 的行@。因此,i-th 行和j-th 列(四个数组中的每一个)的数组值(所有四个数组)是j-th 之间的比较标准中使用的值r 的行和t 的i-th 行。此复制过程的实现使用 R 函数rep。例如,在计算r.rt 时,使用rep 和each=M,其作用是将其数组输入r[,1] 视为行向量,并将该行复制M 次以形成M 行。结果是,与r 中的一行相对应的每一列都具有来自r 相应行的RT. 值,并且该值对于r.rt 的所有行(该列的)都相同,每一个都对应t中的一行。这意味着在将r 中的该行与t 中的任何行进行比较时,将使用r 中该行的RT. 的值。相反,在计算t.rt 时,使用rep 和times=N,其作用是将其数组输入视为列向量并将该列复制N 次以形成N 列。结果是,t.rt 中的每一行(对应于t 中的一行)具有来自t 相应行的RT. 值,并且该值对于(该行的)所有列都是相同的t.rt 中的每一行对应r 中的一行。这意味着在将t 中的该行与r 中的任何行进行比较时,将使用t 中该行的RT. 的值。同样,r.mz 和 t.mz 的计算分别使用来自 r 和 t 的 m.z. 列。
这将执行矢量化比较,生成 M 与 N 的逻辑矩阵,其中 i-th 行和 j-th 列是 TRUE,如果 j-th r 的行与i 的第t 行匹配条件,否则FALSE。 which() 的输出是此逻辑比较结果矩阵的数组索引数组,其中元素为TRUE。我们希望将这些数组索引转换为比较结果矩阵的行索引和列索引,以引用r 和t 的行。下一行从数组索引中提取列索引。请注意,变量名称是r.ind,表示这些对应于r 的行。我们首先提取它,因为它对于检测r 中的一行的多个匹配项很重要。
这部分为r 中的给定行处理t 中可能的多个匹配项。多个匹配项将在r.ind 中显示为重复值。如上所述,这里的逻辑只保留t 中增加行的第一个匹配项。函数duplicated 返回数组中所有重复值的索引。因此删除这些元素将做我们想要的。代码首先从r.ind 中删除它们,然后从ind 中删除它们,最后使用修剪后的ind 和@987654438 计算与t 的行相对应的比较结果矩阵的列索引@。 match.interval 返回的是一个矩阵,其行是匹配的行索引对,其第一列是r 的行索引,第二列是t 的行索引。
get.area.matched 函数仅使用来自match.ind 的结果从t 中提取Area 以获取所有匹配项。请注意,返回的结果是一个(列)向量,其长度等于r 中的行数,并初始化为NA。这样,r 中与 t 不匹配的行将返回 Area 的 NA。
这使用lapply 在列表txt.import.matrix 上应用函数get.area.matched,并将返回的匹配Area 结果附加到reduced.list.pre.filtering 作为列向量。同样,相应的列名也被附加并设置在结果res中。
编辑:使用foreach 包的替代实现
事后看来,更好的实现是使用foreach 包来对比较进行矢量化处理。在此实现中,需要 foreach 和 magrittr 包
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))
}
这更容易解析并反映正在执行的内容。请注意,嵌套 foreach 与 cbind 组合返回的内容再次是一个逻辑矩阵。最后,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,"]"))}))
希望这会有所帮助。