比较并获得行之间的区间交点

Question

我有一个类似下面的数据库。

pos1<-c(5,15,25,40,80,5,18,22,38,84,5,16,50,92,31,50,20,30,50,70,27,50,60,50,90,20,40)
pos2<-c(10,17,30,42,90,10,20,24,42,87,10,19,52,100,40,70,25,32,60,90,30,60,71,60,100,25,50)
chr<-c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2)
n<-c(25,65,78,56,35,78,58,98,14,25,65,85,98,74,20,36,48,98,52,69,21,47,53,10,12,37,82)
pop<-c("A","A","A","A","A","B","B","B","B","C","C","C","C","C","D","D","A","A","A","A","B","B","B","C","C","D","D")
data<-data.frame(pos1,pos2,chr,pop,n)

位置 1 和位置 2 为每个 chr 和总体设计了一个区间的起点和终点。我的目的是获取流行音乐 A、B 和 C（不是 D）之间的哪些区间相交，以及每个群体的哪些区间是唯一的。

因此，对于唯一的间隔，我将有一个结果 data.frame，如下所示：

pos1.u<-c(25,50,92,20,30,27,90)
pos2.u<-c(30,52,100,25,32,30,100)
chr.u<-c(1,1,1,2,2,2,2)
pop.u<-c("A","B","C","A","A","B","C")
n.u<-c(78,98,74,48,98,21,12)
data.u<-data.frame(pos1.u,pos2.u,chr.u,pop.u,n.u)

对于这 3 个总体之间相交的区间，data.frame 如下所示：

pos1.c<-c(5,15,40,80,5,38,85,5,16,50,70,50,60,50)
pos2.c<-c(10,17,42,90,10,42,87,10,19,60,90,60,71,60)
chr.c<-c(1,1,1,1,1,1,1,1,1,2,2,2,2,2)
pop.c<-c("A","A","A","A","B","B","B","C","C","A","A","B","B","C")
n.c<-c(25,65,56,35,78,14,25,65,85,52,69,47,53,10)
data.c<-data.frame(pos1.c,pos2.c,chr.c,pop.c,n.c)

我不知道如何编写一个恰好可以做到这一点的脚本，你能帮我吗？

Answer 1

我认为以下代码可以满足您的要求，尽管它产生的结果与您的不同 - 所以请仔细检查！我认为差异在于开区间和闭区间的定义。以下假设不包括任何端点，而我怀疑这可能不是您的意思（否则 (15,18) 和 (17,19) 不会算作重叠，因为两者都没有整数值） .因此，您可能需要调整下面的打开/关闭定义。

pos1<-c(5,15,25,40,80,5,18,22,38,84,5,16,50,92,31,50,20,30,50,70,27,50,60,50,90,20,40)
pos2<-c(10,17,30,42,90,10,20,24,42,87,10,19,52,100,40,70,25,32,60,90,30,60,71,60,100,25,50)
chr<-c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2)
n<-c(25,65,78,56,35,78,58,98,14,25,65,85,98,74,20,36,48,98,52,69,21,47,53,10,12,37,82)
pop<-c("A","A","A","A","A","B","B","B","B","C","C","C","C","C","D","D","A","A","A","A","B","B","B","C","C","D","D")
data<-data.frame(pos1,pos2,chr,pop,n,stringsAsFactors = FALSE)

library(intervals)
data<-data[data$pop!="D",] #remove irrelevant D entries
rownames(data) <- seq_len(nrow(data)) #reset rownames to allow for removed Ds

#set ints as a list of intervals (as required by intervals package)
ints <- tapply(1:nrow(data),data$pop,function(v) 
         Intervals(as.matrix(data[v,c("pos1","pos2")]),
         closed=c(FALSE,FALSE), #this is where you adjust open/closed lower and upper ends of the intervals - TRUE means end value included
         type="Z")) #Z is integers
pops <- unique(data$pop) #unique values of pop
popidx <- lapply(pops,function(x) which(data$pop==x)) #list of indices of these values in data
names(popidx) <- pops

#sets is a df of all pairwise combinations to check
sets <- expand.grid(pops,pops,stringsAsFactors = FALSE) 
sets <- sets[sets$Var1!=sets$Var2,]

olap <- lapply(1:nrow(sets),function(i) 
        interval_overlap(ints[[sets$Var1[i]]],ints[[sets$Var2[i]]])) #list of overlaps
olap <- lapply(1:nrow(sets),function(i) {
  df<-as.data.frame(olap[[i]],stringsAsFactors=FALSE)
  df$pos1 <- as.numeric(rownames(df))
  df$pos2 <- sapply(1:nrow(df),function(j) popidx[[sets$Var2[i]]][df[j,1][[1]][1]])
  return(df)}) #tidy up as dfs, with correct indices in data (rather than in ints)
olap <- do.call(rbind,olap)[,-1] #join dataframes
olap$olaps <- !is.na(olap$pos2) #identify those with overlaps

#group by unique pos1 and identify max and min no of overlaps with other groups
olap <- data.frame(minoverlap=tapply(olap$olaps,olap$pos1,min),maxoverlap=tapply(olap$olaps,olap$pos1,max))
olap$rowno <- as.numeric(rownames(olap))

uniques <- data[olap$rowno[olap$maxoverlap==0],] #intervals appearing in just one pop
commons <- data[olap$rowno[olap$minoverlap>0],] #intervals with an overlap in all other pops