2014-02-09 42 views
8

我的數據是由V6的ID分組而成的羣中的每個元素的「指數」,並下令位置(V1:V3):創建用於data.table

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中的每個組的第一個元素是「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?

dput(dt) 
structure(list(V1 = c("chr1", "chr1", "chr1", "chr1", "chr1", 
"chr1", "chr1", "chr1", "chr1", "chr1"), 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 = c(".", 
".", ".", ".", ".", ".", ".", ".", ".", "."), V5 = c("+", "+", 
"-", "-", "-", "-", "-", "-", "+", "+"), V6 = c("ENSMUSG00000090025", 
"ENSMUSG00000064842", "ENSMUSG00000051951", "ENSMUSG00000051951", 
"ENSMUSG00000051951", "ENSMUSG00000051951", "ENSMUSG00000051951", 
"ENSMUSG00000051951", "ENSMUSG00000089699", "ENSMUSG00000089699" 
)), .Names = c("V1", "V2", "V3", "V4", "V5", "V6"), class = c("data.table", 
"data.frame"), row.names = c(NA, -10L), .internal.selfref = <pointer: 0x1de6a88>) 
+0

+1非常很好的框架問題。 – Arun

+0

「提出好問題,獲得好答案」應該是stackoverflow的座右銘:) – fridaymeetssunday

回答

14

作爲一位生物信息學家,我經常遇到這種手術。這是我喜歡data.table修改行的子集引用功能!

我會做這樣的:

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

不需要的功能。這更有利一些,因爲它避免了必須爲每個組檢查=="+""-"一次!相反,你可以第一子集中所有組,+一次然後按V6和修改的地方只是那些行

同樣,你再次爲"-"做。希望有所幫助。

注意:.N是一個特殊的變量,它包含每個組的觀察次數。

+1

謝謝 - 我總是驚訝data.table速度有多快。我也和.N一起玩,但從來沒有接近解決方案。 – fridaymeetssunday

3

首先,我會加載示例數據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您的問題:

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))) 
) 

但是這樣做會相對較慢,因爲您必須製作兩份數據。

另一種方法是使用一個,如果彙總內:

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

謝謝@hadley添加dplyr解決方案。我還沒有看看包裝,但這可能是一個開始。 – fridaymeetssunday

+0

@fridaymeetssunday如果你熟悉plyr,轉換到dplyr應該很容易。 – hadley