10

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

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

我想我可以像这样使用 IRange 包：

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

但是后来我无法获得折叠的列：我尝试过使用 findOvarlaps 但这

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

但我认为这是不对的。

任何帮助将不胜感激。

4

2 回答 2

12

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：这里折叠的变量名不是字符串而是因子列表。例如，这比使用粘贴的折叠字符更有效且更易于访问。

OP 澄清后编辑，我将通过 chrom 创建组变量。我的意思是现在为每个 chrom 组调用 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
于 2013-06-06T09:19:06.923 回答
5

对数据进行排序后，您可以轻松测试一个区间是否与前一个区间重叠，并为每组重叠区间分配一个标签。一旦你有了这些标签，你就可以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
于 2013-06-06T09:11:50.160 回答