为data.table中的每个元素创建一个“索引”。

Question

我的数据按照V6中的ID进行分组，并按位置（V1:V3）排序：

dt
      V1      V2      V3 V4 V5                 V6
 1: chr1 3054233 3054733  .  + ENSMUSG00000090025
 2: chr1 3102016 3102125  .  + ENSMUSG00000064842
 3: chr1 3205901 3207317  .  - ENSMUSG00000051951
 4: chr1 3206523 3207317  .  - ENSMUSG00000051951
 5: chr1 3213439 3215632  .  - ENSMUSG00000051951
 6: chr1 3213609 3216344  .  - ENSMUSG00000051951
 7: chr1 3214482 3216968  .  - ENSMUSG00000051951
 8: chr1 3421702 3421901  .  - ENSMUSG00000051951
 9: chr1 3466587 3466687  .  + ENSMUSG00000089699
10: chr1 3513405 3513553  .  + ENSMUSG00000089699

我想做的是在每个V6组中，按位置添加一个额外的带索引的列。也就是说，每个V6组的第一个元素为“1”，第二个元素为“2”，以此类推。 我可以使用ddply和自定义函数实现这一点：

rankExons <- function(x){
   if(unique(x$V5) == "+"){ 
         x$index <- seq_len(nrow(x))}
   else{
         x$index <- rev(seq_len(nrow(x)))}
   x
}

indexed <- ddply(dt, .(V6), rankExons)
indexed
     V1      V2      V3 V4 V5                 V6 index
1  chr1 3205901 3207317  .  - ENSMUSG00000051951     6
2  chr1 3206523 3207317  .  - ENSMUSG00000051951     5
3  chr1 3213439 3215632  .  - ENSMUSG00000051951     4
4  chr1 3213609 3216344  .  - ENSMUSG00000051951     3
5  chr1 3214482 3216968  .  - ENSMUSG00000051951     2
6  chr1 3421702 3421901  .  - ENSMUSG00000051951     1
7  chr1 3102016 3102125  .  + ENSMUSG00000064842     1
8  chr1 3466587 3466687  .  + ENSMUSG00000089699     1
9  chr1 3513405 3513553  .  + ENSMUSG00000089699     2
10 chr1 3054233 3054733  .  + ENSMUSG00000090025     1

不幸的是，在完整数据集（约620k行）上运行非常缓慢，而且使用并行时会崩溃。

library(doMC)
registerDoMC(cores=6)
indexed <- ddply(dt, .(V6), rankExons, .parallel=TRUE)
Error: serialization is too large to store in a raw vector
Error: serialization is too large to store in a raw vector
Error: serialization is too large to store in a raw vector
Error: serialization is too large to store in a raw vector
Error: serialization is too large to store in a raw vector
Error: serialization is too large to store in a raw vector
Warning message:
In mclapply(argsList, FUN, mc.preschedule = preschedule, mc.set.seed = set.seed,  :
  all scheduled cores encountered errors in user code

所以，我尝试使用data.table，但无法使其正常工作。这是我尝试过的：

setkey(dt, "V6")

dt[,index:=rankExons(dt), by=V6]
dt[,rankExons(.sd), by=V6, .SDcols=c("V5, V6")]

两者都失败了。我该如何使用data.table重新创建我的ddply？

dt <- data.table::data.table(
  V1 = rep("chr1", 10L),
  V2 = c(
    3054233L, 3102016L, 3205901L, 3206523L, 3213439L, 3213609L,
    3214482L, 3421702L, 3466587L, 3513405L
  ),
  V3 = c(
    3054733L, 3102125L, 3207317L, 3207317L, 3215632L, 3216344L,
    3216968L, 3421901L, 3466687L, 3513553L
  ),
  V4 = rep(".", 10L),
  V5 = rep(c("+", "-", "+"), c(2L, 6L, 2L)),
  V6 = rep(
    c(
      "ENSMUSG00000090025", "ENSMUSG00000064842", "ENSMUSG00000051951",
      "ENSMUSG00000089699"
    ),
    c(1L, 1L, 6L, 2L)
  )
)

- fridaymeetssunday

---

“提出好问题，获得好答案”应该是 Stack Overflow 的座右铭 :) - fridaymeetssunday

Answer 1

作为一名生物信息学家同行，我经常遇到这种操作。这正是我欣赏 data.table 的 modify subset of rows by reference 功能的原因所在！

我会像这样做：

dt[V5 == "+", index := 1:.N, by=V6]
dt[V5 == "-", index := .N:1, by=V6]

不需要函数。这更有优势，因为它避免了在每个组内检查== "+"或"-"一次！相反，您可以首先使用+一次对所有组进行子集操作，然后按V6进行分组，并在原位修改只有这些行！

同样，您可以再次使用"-"进行操作。希望这能帮到您。

注意：.N是一个特殊变量，包含每个组的观察数。

- Arun

---

谢谢 - 我总是惊叹于 data.table 的速度。我也尝试过使用 .N，但从未接近解决方案。 - fridaymeetssunday

Answer 2

首先，我将加载您的样本数据到R中（目前无法使用dput()与data.table）：

df <- read.table(header = TRUE, stringsAsFactors = FALSE, text = "
V1      V2      V3 V4 V5                 V6
1  chr1 3205901 3207317  .  - ENSMUSG00000051951
2  chr1 3206523 3207317  .  - ENSMUSG00000051951
3  chr1 3213439 3215632  .  - ENSMUSG00000051951
4  chr1 3213609 3216344  .  - ENSMUSG00000051951
5  chr1 3214482 3216968  .  - ENSMUSG00000051951
6  chr1 3421702 3421901  .  - ENSMUSG00000051951
7  chr1 3102016 3102125  .  + ENSMUSG00000064842
8  chr1 3466587 3466687  .  + ENSMUSG00000089699
9  chr1 3513405 3513553  .  + ENSMUSG00000089699
10 chr1 3054233 3054733  .  + ENSMUSG00000090025")

您几乎可以用dplyr优雅地解决您的问题：

使用dplyr，您几乎可以优雅地解决问题：

library(dplyr)

df %>% 
  group_by(V6, V5) %>%
  mutate(index = row_number(V2))

我假设V2是您想要按索引的变量 - 我认为明确说明比依赖于行顺序更好

但是，您希望针对不同的子集使用不同的摘要，这在dplyr中目前并不容易。一种方法是拆分然后重新组合：

rbind_list(
  df %>% filter(V5 == "+") %>% mutate(index = row_number(V2)),
  df %>% filter(V5 == "-") %>% mutate(index = row_number(desc(V2)))
)

但是这种方法会相对较慢，因为您需要复制数据两次。

另一种方法是在摘要内使用 if：

df %>% 
  group_by(V6, V5) %>%
  mutate(index = row_number(if (V5[1] == "+") V2 else desc(V2)))

- hadley

---

感谢@hadley添加了一个dplyr解决方案。我还没有研究过这个包，但这可能是一个好的开始。 - fridaymeetssunday

如果你熟悉plyr，转换到dplyr应该很容易。@fridaymeetssunday - hadley