减少脚本运行的时间（基于循环）

Question

我编写了一个基于两个 for 循环的脚本，我想对其进行优化以加快其运行时间。

以下是我用我在自己的数据上使用的代码简化的可重现数据。

nuc 是一个具有 101 个“位置”的向量，并且 tel 是坐标“aa”和“bb”不同的数据框

目的是为每个位置计算每个位置包含在每个 aa 和 bb 坐标之间的次数。例如位置 111 包含在 3 对坐标之间：G、I 和 J

#data

tel=data.frame(aa=c(153,113,163,117,193,162,110,109,186,103),
               bb=c(189,176,185,130,200,189,156,123,198,189),
               ID=c("A", "B", "C", "D", "E", "F", "G", "H", "I", "J"))

> tel
    aa  bb ID
1  153 189  A
2  113 176  B
3  163 185  C
4  117 130  D
5  193 200  E
6  162 189  F
7  110 156  G
8  109 123  H
9  186 198  I
10 103 189  J

nuc2=100:200

# Loop

count_occ=0
count_occ_int=NULL
count_occ_fin=NULL

for (j in 1:length(nuc2)){
  for (i in 1:nrow(tel)) {
    
    if (nuc2[j]< tel$bb[i] & nuc2[j]>tel$aa[i]) 
    {count_occ=count_occ+1}
  }
  count_occ_int=count_occ
  count_occ_fin=c(count_occ_fin,count_occ_int)
  count_occ=0
}

nuc_occ=data.frame(nuc=nuc2, occ=count_occ_fin)

> head(nuc_occ,20)
   nuc occ
1  100   0
2  101   0
3  102   0
4  103   0
5  104   1
6  105   1
7  106   1
8  107   1
9  108   1
10 109   1
11 110   2
12 111   3
13 112   3
14 113   3
15 114   4
16 115   4
17 116   4
18 117   4
19 118   5
20 119   5

在我的数据中，我的 nuc 向量的长度是 9304567，坐标对的数量是 53（我很快就会有几百个），运行代码花了 60 多个小时！！

有什么想法可以帮助我加快这段代码的速度吗？

我虽然使用了 apply 函数，但我不确定如何将这两个 for 循环操作结合起来。

Answer 1

这是一个 tidyverse 解决方案：

lapply(
  100:200,
  \(x) tel %>% 
         filter(aa <= x & x <= bb) %>% 
         summarise(occ=n(), .groups="drop") %>% 
         add_column(nuc=x, .before=1)
) %>% 
bind_rows() %>% 
as_tibble()
# A tibble: 101 × 2
     nuc   occ
   <int> <int>
 1   100     0
 2   101     0
 3   102     0
 4   103     1
 5   104     1
 6   105     1
 7   106     1
 8   107     1
 9   108     1
10   109     2
# … with 91 more rows

使用microbenchmark 来评估性能，这给出了

Unit: nanoseconds
     expr min lq mean median uq max neval
   lapply   7  9  8.8      9  9   9    10
 original   8  9 23.8      9  9 158    10

换句话说，速度降低了大约三分之二。 tidyverse 并不以速度着称。基础 R 解决方案可能会更快。

Answer 2

您可以像这样使用 data.table 非等连接：

library(data.table)
setDT(tel)[SJ(v=nuc2), on=.(aa<=v, bb>=v)][,.(occ = sum(!is.na(ID))), by=.(nuc=aa)]

解释：

• setDT(tel) 将 tel data.frame 设置为 data.table 类
• SJ(v=nuc2) 是将向量转换为 data.table 的便捷函数；在这种情况下，将nuc2 转换为具有一列v 的data.table。我这样做是因为我想加入两个 data.tables，一个是 tel（包含列 aa、bb 和 v），另一个是包含单个列 v 的值nuc2
• 连接条件在setDT(tel)[...]子句的on=..参数中；这里的连接条件是 v 值必须 >= aa 值并且必须是 bb 值
• 最后一步（即下一个链接的 data.table 操作）仅按nuc 值 (by=.(nuc=aa)) 计算 ID 不为 NA 的行数

输出：

       nuc   occ
     <int> <int>
  1:   100     0
  2:   101     0
  3:   102     0
  4:   103     1
  5:   104     1
 ---            
 97:   196     2
 98:   197     2
 99:   198     2
100:   199     1
101:   200     1