在R中折叠相交区域

Question

我试图找到一种方法来折叠具有相交范围的行，用“开始”和“停止”列表示，并将折叠值记录到新列中。例如，我有这个数据帧：

my.df<- data.frame(chrom=c(1,1,1,1,14,16,16), name=c("a","b","c","d","e","f","g"), start=as.numeric(c(0,70001,70203,70060, 40004, 50000872, 50000872)), stop=as.numeric(c(71200,71200,80001,71051, 42004, 50000890, 51000952))) 


chrom name start stop 
1 a  0 71200 
1 b 70001 71200 
1 c 70203 80001 
1 d 70060 71051 
14 e 40004 42004 
16 f 50000872 50000890 
16 g 50000872 51000952 

，我试图找到重叠的范围，并记录在“开始”和“停止”的坍塌重叠行和折叠的行的名称所覆盖的最大范围，所以我会得到这样的：

chrom start stop  name 
1 70001 80001 a,b,c,d 
14 40004 42004 e 
16 50000872 51000952 f,g 

我想我可以用包IRanges这样的：

library(IRanges) 
ranges <- split(IRanges(my.df$start, my.df$stop), my.df$chrom) 

但后来我有麻烦塌陷列：我有线索d与findOvarlaps但这

ov <- findOverlaps(ranges, ranges, type="any") 

但我不认为这是正确的。

任何帮助将不胜感激。

谢谢！ -fra

Answer 1

排序完数据后，您可以轻松测试间隔是否与前一个重叠，并为每组重叠间隔分配一个标签。 一旦你有这些标签，你可以使用ddply来聚合数据。

d <- data.frame(
    chrom = c(1,1,1,14,16,16), 
    name = c("a","b","c","d","e","f"), 
    start = as.numeric(c(70001,70203,70060, 40004, 50000872, 50000872)), 
    stop = as.numeric(c(71200,80001,71051, 42004, 50000890, 51000952)) 
) 

# Make sure the data is sorted 
d <- d[ order(d$start), ] 

# Check if a record should be linked with the previous 
d$previous_stop <- c(NA, d$stop[-nrow(d)]) 
d$previous_stop <- cummax(ifelse(is.na(d$previous_stop),0,d$previous_stop)) 
d$new_group <- is.na(d$previous_stop) | d$start >= d$previous_stop 

# The number of the current group of records is the number of times we have switched to a new group 
d$group <- cumsum(d$new_group) 

# We can now aggregate the data 
library(plyr) 
ddply( 
    d, "group", summarize, 
    start=min(start), stop=max(stop), name=paste(name,collapse=",") 
) 
# group start  stop name 
# 1  1  0 80001 a,d,c,b 
# 2  2 50000872 51000952  e,f 

但是这忽略了chrom柱：考虑到它，你可以做同样的事，每个染色体，分别。

d <- d[ order(d$chrom, d$start), ] 
d <- ddply(d, "chrom", function(u) { 
    x <- c(NA, u$stop[-nrow(u)]) 
    y <- ifelse(is.na(x), 0, x) 
    y <- cummax(y) 
    y[ is.na(x) ] <- NA 
    u$previous_stop <- y 
    u 
}) 
d$new_group <- is.na(d$previous_stop) | d$start >= d$previous_stop 
d$group <- cumsum(d$new_group) 
ddply( 
    d, .(chrom,group), summarize, 
    start=min(start), stop=max(stop), name=paste(name,collapse=",") 
) 
# chrom group start  stop name 
# 1  1  1  0 80001 a,c,b 
# 2 14  2 40004 42004  d 
# 3 16  3 50000872 51000952 e,f 

Answer 2

IRanges是一个很好的候选人这样的工作。不需要使用chrom变量。

ir <- IRanges(my.df$start, my.df$stop) 
## I create a new grouping variable Note the use of reduce here(performance issue) 
my.df$group2 <- subjectHits(findOverlaps(ir, reduce(ir))) 
# chrom name start  stop group2 
# 1  1 a 70001 71200  2 
# 2  1 b 70203 80001  2 
# 3  1 c 70060 71051  2 
# 4 14 d 40004 42004  1 
# 5 16 e 50000872 50000890  3 
# 6 16 f 50000872 51000952  3 

新的group2变量是范围指示符。现在，使用data.table我无法将数据转换为所需的输出：

library(data.table) 
DT <- as.data.table(my.df) 
DT[, list(start=min(start),stop=max(stop), 
     name=list(name),chrom=unique(chrom)), 
       by=group2] 

# group2 start  stop name chrom 
# 1:  2 70001 80001 a,b,c  1 
# 2:  1 40004 42004  d 14 
# 3:  3 50000872 51000952 e,f 16 

PS：这里倒塌的变量名不是字符串，但一个名单因素的。这比使用粘贴的collapased角色更高效，更易于访问。

编辑因为OP的说明，我会通过chrom创建组varibale。我的意思是现在为每个染色体组调用Iranges代码。我稍微修改你的数据，创建一组同一染色体的区间。

my.df<- data.frame(chrom=c(1,1,1,1,14,16,16), 
        name=c("a","b","c","d","e","f","g"), 
        start=as.numeric(c(0,3000,70203,70060, 40004, 50000872, 50000872)), 
        stop=as.numeric(c(1,5000,80001,71051, 42004, 50000890, 51000952))) 

library(data.table) 
DT <- as.data.table(my.df) 

## find interval for each chromsom 
DT[,group := { 
     ir <- IRanges(start, stop); 
     subjectHits(findOverlaps(ir, reduce(ir))) 
     },by=chrom] 

## Now I group by group and chrom 
DT[, list(start=min(start),stop=max(stop),name=list(name),chrom=unique(chrom)), 
    by=list(group,chrom)] 

    group chrom start  stop name chrom 
1:  1  1  0  1 a  1 
2:  2  1  3000  5000 b  1 
3:  3  1 70060 80001 c,d  1 
4:  1 14 40004 42004 e 14 
5:  1 16 50000872 51000952 f,g 16