Find the intersection of overlapping ranges in two tables using data.table function foverlaps
@Seth provided the fastest way to solve the problem of intersection overlaps using the data.table foverlaps function. However, this solution did not take into account the fact that the input bed files may have overlapping ranges that needed to be reduced into single regions. @Martin Morgan solved that with his solution using the GenomicRanges package, that did both the intersecting and range reducing. However, Martin's solution didn't use the foverlaps function. @Arun pointed out that the overlapping ranges in different rows within a table was not currently possible using foverlaps. Thanks to the answers provided, and some additional research on stackoverflow, I came up with this hybrid solution.
Create example BED files without overlapping regions within each file.
chr <- c(1:22,"X","Y","MT")
#bedA contains 5 million rows
bedA <- data.table(
CHROM = as.vector(sapply(chr, function(x) rep(x,200000))),
START = rep(as.integer(seq(1,200000000,1000)),25),
STOP = rep(as.integer(seq(500,200000000,1000)),25),
key = c("CHROM","START","STOP")
)
#bedB contains 500 thousand rows
bedB <- data.table(
CHROM = as.vector(sapply(chr, function(x) rep(x,20000))),
START = rep(as.integer(seq(200,200000000,10000)),25),
STOP = rep(as.integer(seq(600,200000000,10000)),25),
key = c("CHROM","START","STOP")
)
Now create a new bed file containing the intersecting regions in bedA and bedB.
#This solution uses foverlaps
system.time(tmpA <- intersectBedFiles.foverlaps(bedA,bedB))
user system elapsed
1.25 0.02 1.37
#This solution uses GenomicRanges
system.time(tmpB <- intersectBedFiles.GR(bedA,bedB))
user system elapsed
12.95 0.06 13.04
identical(tmpA,tmpB)
[1] TRUE
Now, modify bedA and bedB such that they contain overlapping regions:
#Create overlapping ranges
makeOverlaps <- as.integer(c(0,0,600,0,0,0,600,0,0,0))
bedC <- bedA[, STOP := STOP + makeOverlaps, by=CHROM]
bedD <- bedB[, STOP := STOP + makeOverlaps, by=CHROM]
Test time to intersect bed files with overlapping ranges using either the foverlaps or GenomicRanges fucntions.
#This solution uses foverlaps to find the intersection and then run GenomicRanges on the result
system.time(tmpC <- intersectBedFiles.foverlaps(bedC,bedD))
user system elapsed
1.83 0.05 1.89
#This solution uses GenomicRanges
system.time(tmpD <- intersectBedFiles.GR(bedC,bedD))
user system elapsed
12.95 0.04 12.99
identical(tmpC,tmpD)
[1] TRUE
The winner: foverlaps!
FUNCTIONS USED
This is the function based upon foverlaps, and will only call the GenomicRanges function (reduceBed.GenomicRanges) if there are overlapping ranges (which are checked for using the rowShift function).
intersectBedFiles.foverlaps <- function(bed1,bed2) {
require(data.table)
bedKey <- c("CHROM","START","STOP")
if(nrow(bed1)>nrow(bed2)) {
bed <- foverlaps(bed1, bed2, nomatch = 0)
} else {
bed <- foverlaps(bed2, bed1, nomatch = 0)
}
bed[, START := pmax(START, i.START)]
bed[, STOP := pmin(STOP, i.STOP)]
bed[, `:=`(i.START = NULL, i.STOP = NULL)]
if(!identical(key(bed),bedKey)) setkeyv(bed,bedKey)
if(any(bed[, STOP+1 >= rowShift(START), by=CHROM][,V1], na.rm = T)) {
bed <- reduceBed.GenomicRanges(bed)
}
return(bed)
}
rowShift <- function(x, shiftLen = 1L) {
#Note this function was described in this thread:
#http://stackoverflow.com/questions/14689424/use-a-value-from-the-previous-row-in-an-r-data-table-calculation
r <- (1L + shiftLen):(length(x) + shiftLen)
r[r<1] <- NA
return(x[r])
}
reduceBed.GenomicRanges <- function(bed) {
setnames(bed,colnames(bed),bedKey)
if(!identical(key(bed),bedKey)) setkeyv(bed,bedKey)
grBed <- makeGRangesFromDataFrame(bed,
seqnames.field = "CHROM",start.field="START",end.field="STOP")
grBed <- reduce(grBed)
grBed <- data.table(
CHROM=as.character(seqnames(grBed)),
START=start(grBed),
STOP=end(grBed),
key = c("CHROM","START","STOP"))
return(grBed)
}
This function strictly used the GenomicRanges package, produces the same result, but is about 10 fold slower that the foverlaps funciton.
intersectBedFiles.GR <- function(bed1,bed2) {
require(data.table)
require(GenomicRanges)
bed1 <- makeGRangesFromDataFrame(bed1,
seqnames.field = "CHROM",start.field="START",end.field="STOP")
bed2 <- makeGRangesFromDataFrame(bed2,
seqnames.field = "CHROM",start.field="START",end.field="STOP")
grMerge <- suppressWarnings(intersect(bed1,bed2))
resultTable <- data.table(
CHROM=as.character(seqnames(grMerge)),
START=start(grMerge),
STOP=end(grMerge),
key = c("CHROM","START","STOP"))
return(resultTable)
}
An additional comparison using IRanges
I found a solution to collapse overlapping regions using IRanges but it is more than 10 fold slower than GenomicRanges.
reduceBed.IRanges <- function(bed) {
bed.tmp <- bed
bed.tmp[,group := {
ir <- IRanges(START, STOP);
subjectHits(findOverlaps(ir, reduce(ir)))
}, by=CHROM]
bed.tmp <- bed.tmp[, list(CHROM=unique(CHROM),
START=min(START),
STOP=max(STOP)),
by=list(group,CHROM)]
setkeyv(bed.tmp,bedKey)
bed[,group := NULL]
return(bed.tmp[, -(1:2)])
}
system.time(bedC.reduced <- reduceBed.GenomicRanges(bedC))
user system elapsed
10.86 0.01 10.89
system.time(bedD.reduced <- reduceBed.IRanges(bedC))
user system elapsed
137.12 0.14 137.58
identical(bedC.reduced,bedD.reduced)
[1] TRUE
Find overlap of multiple ranges in data.table
Here's a one-line example:
library(data.table)
dt = data.table(a = c(3,4,5), b = c(13,12,19))
dt[, c("o.l", "o.r") := as.list(range(Reduce(intersect, mapply(seq, a, b, 1))))]
dt
# a b o.l o.r
# 1: 3 13 5 12
# 2: 4 12 5 12
# 3: 5 19 5 12
Where the core of the problem is
dt = data.table(a = c(3,4,5), b = c(13,12,19))
dt[, Reduce(intersect, mapply(seq, a, b, 1))]
# [1] 5 6 7 8 9 10 11 12
Finding all overlaps in one iteration of foverlap in R's data.table
The IRanges package on Bioconductor from which data.table's foverlaps() was inspired has some handy functions for questions like this.
Perhaps, reduce() might be the function you are looking for to merge all overlapping periods:
library(data.table)
dt = data.table(start = c(1,2,4,6,8,10),
end = c(2,3,6,8,10,12))
library(IRanges)
ir <- IRanges(dt$start, dt$end)
ir
IRanges object with 6 ranges and 0 metadata columns:
start end width
<integer> <integer> <integer>
[1] 1 2 2
[2] 2 3 2
[3] 4 6 3
[4] 6 8 3
[5] 8 10 3
[6] 10 12 3
reduce(ir, min.gapwidth = 0L)
IRanges object with 2 ranges and 0 metadata columns:
start end width
<integer> <integer> <integer>
[1] 1 3 3
[2] 4 12 9
as.data.table(reduce(ir, min.gapwidth = 0L))
start end width
1: 1 3 3
2: 4 12 9
On Bioconductor, there is a comprehensive Introduction to IRanges available.
Edit: The OP has supplied a second sample data set which includes an id column and is asking if IRanges has support for joining intervals by id.
Adding data to IRanges seems to specialize quickly into the area of genome research which is terra incognita to me. However, I've found the following approach using IRanges:
Grouping with IRanges
library(data.table)
# 2nd sample data set provided by the OP
dt = data.table(id = c('A','A','A','A','A','B','B','B'),
eventStart = c(1,2,4,6,8,10,11,15),
eventEnd = c(2,3,6,8,10,12,14,16))
library(IRanges)
# set names when constructing IRanges object
ir <- IRanges(dt$eventStart, dt$eventEnd, names = dt$id)
lapply(split(ir, names(ir)), reduce, min.gapwidth = 0L)
$A
IRanges object with 2 ranges and 0 metadata columns:
start end width
<integer> <integer> <integer>
[1] 1 3 3
[2] 4 10 7
$B
IRanges object with 2 ranges and 0 metadata columns:
start end width
<integer> <integer> <integer>
[1] 10 14 5
[2] 15 16 2
Converting this back to data.table leads to a rather clumsy piece of code:
ir <- IRanges(dt$eventStart, dt$eventEnd, names = dt$id)
rbindlist(lapply(split(ir, names(ir)),
function(x) as.data.table(reduce(x, min.gapwidth = 0L))),
idcol = "id")
id start end width
1: A 1 3 3
2: A 4 10 7
3: B 10 14 5
4: B 15 16 2
Grouping within data.table
We can get the same result with a less convoluted code if we group within data.table and apply reduce() on the single chunks:
dt[, as.data.table(reduce(IRanges(eventStart, eventEnd), min.gapwidth = 0L)), id]
R Overlapping period for event with two identifiers
I had a hard time wrapping my head around the solution for the other scenario,
so I'm not sure how to adjust that version,
but I believe you can achieve the result you want in this case with an inner non-equi join:
dt <- data.table(df, key = c("FirmID", "ID", "start", "end"))
dt[, firm_total := .N, by = "FirmID"][
dt, .(FirmID, firm_total), on = .(FirmID, ID < ID, start <= end, end >= start), nomatch = NULL, mult = "all"][
, .(n = .N, pct = .N / (firm_total[1L] * (firm_total[1L] - 1L) / 2L)), by = "FirmID"]
We avoid redundant rows in the join with ID < ID
(and note that here that should be interpreted as left-hand side column's ID < right-hand side column's ID).
Redundancy can happen because,
if ID x overlaps with ID y,y overlaps with x.
If you think of all pairs of IDs and put them in a matrix,
the maximum number of overlaps would be the number of elements in the lower triangular,
which can be calculated with n * (n - 1) / 2,
that's why we initially add firm_total.
I didn't do extensive testing but this version may be better for this scenario.
The documentation of foverlaps states that it's mainly targeted at joins where one table is much smaller than the other one,
and explains why it may be an expensive operation.
You're doing a self-join,
so both tables are the same size.
And a table.express version of the solution because why not:
library(table.express)
dt %>%
group_by(FirmID) %>%
mutate(firm_total = .N) %>%
inner_join(dt, FirmID, ID < ID, start <= end, end >= start, .expr = TRUE) %>%
select(FirmID, firm_total) %>%
group_by(FirmID) %>%
summarize(n = .N, pct = .N / (firm_total[1L] * (firm_total[1L] - 1L) / 2))
EDIT: and if you want your overlap column you could compute that with size of lower triangular - n.
Find overlaps in different datatables
I think for this you want to use set operations.
Set Operations
This should work:
dat1 <- data.frame(gen = c("aaaaa", "1494_f_at", "1111", "!!!!"), stringsAsFactors = FALSE)
dat2 <- data.frame(gen = c("1494_f_at", "cccccc", "!!!!","999"), stringsAsFactors = FALSE)
dat3 <- data.frame(gen = c( "!!!!","1494_f_at", "999", "dddddd"), stringsAsFactors = FALSE)
intersect(intersect(dat1[,1], dat2[,1]), dat3[,1])
How to join 2 data tables by time interval and summarize overlapping and non-overlapping time periods by factor variable
Not the entire answer (see last paragraph).. but I think this will get you what you want.
library( data.table )
library( lubridate )
set.seed(13)
EffortType = sample(c("A","B","C"), 100, replace = TRUE)
On = sample(seq(as.POSIXct('2016/01/01 01:00:00'), as.POSIXct('2016/01/03 01:00:00'), by = "15 mins"), 100, replace=T)
Off = On + minutes(sample(1:60, 100, replace=T))
Effort1 = data.table(EffortType, On, Off)
EffortType2 = sample(c("A","B","C"), 100, replace = TRUE)
On = sample(seq(as.POSIXct('2016/01/01 12:00:00'), as.POSIXct('2016/01/03 12:00:00'), by = "15 mins"), 100, replace=T)
Off = On + minutes(sample(1:60, 100, replace=T))
Effort2 = data.table(EffortType2, On, Off)
#create DT of minutes, spanning your entire period.
dt.minutes <- data.table( On = seq(as.POSIXct('2016/01/01 01:00:00'), as.POSIXct('2016/01/03 12:00:00'), by = "1 mins"),
Off = seq(as.POSIXct('2016/01/01 01:00:00'), as.POSIXct('2016/01/03 12:00:00'), by = "1 mins") + 60 )
#prep for using foverlaps
setkey(Effort1, On, Off)
setkey(Effort2, On, Off)
#overlap join both efforts on the dt.minutes. note the use of "within" an "nomatch" to throw away minutes without events.
m1 <- foverlaps(dt.minutes, Effort1 ,type="within",nomatch=0L)
m2 <- foverlaps(dt.minutes, Effort2 ,type="within",nomatch=0L)
#bind together
result <- rbindlist(list(m1,m2))[, `:=`(On=i.On, Off = i.Off)][, `:=`(i.On = NULL, i.Off = NULL)]
#cast the result
result.cast <- dcast( result, On + Off ~ EffortType, value.var = "EffortType")
results in
head( result.cast, 10)
# On Off A B C
# 1: 2016-01-01 01:00:00 2016-01-01 01:01:00 1 0 1
# 2: 2016-01-01 01:01:00 2016-01-01 01:02:00 1 0 1
# 3: 2016-01-01 01:02:00 2016-01-01 01:03:00 1 0 1
# 4: 2016-01-01 01:03:00 2016-01-01 01:04:00 1 0 1
# 5: 2016-01-01 01:04:00 2016-01-01 01:05:00 1 0 1
# 6: 2016-01-01 01:05:00 2016-01-01 01:06:00 1 0 1
# 7: 2016-01-01 01:06:00 2016-01-01 01:07:00 1 0 1
# 8: 2016-01-01 01:07:00 2016-01-01 01:08:00 1 0 1
# 9: 2016-01-01 01:08:00 2016-01-01 01:09:00 1 0 1
# 10: 2016-01-01 01:09:00 2016-01-01 01:10:00 1 0 1
Sometimes a event occurs 2-3 times within the same minute, like
# On Off A B C
#53: 2016-01-02 14:36:00 2016-01-02 14:37:00 2 2 3
Not sure on how you want to sum that...
If you can treat them as a single minute, then:
> sum( result.cast[A>0 & B==0, C==0, ] )
[1] 476
> sum( result.cast[A==0 & B>0, C==0, ] )
[1] 386
> sum( result.cast[A==0 & B==0, C>0, ] )
[1] 504
> sum( result.cast[A>0 & B>0, C==0, ] )
[1] 371
> sum( result.cast[A==0 & B>0, C>0, ] )
[1] 341
> sum( result.cast[A>0 & B==0, C>0, ] )
[1] 472
> sum( result.cast[A>0 & B>0, C>0, ] )
[1] 265
will do the trick to get duration in minutes, I think (although this can probably be done in a much smarter way)
Find non-overlapping values from two tables in R
We can use fsetdiff from data.table by including only the common columns in both datasets
fsetdiff(Input, FDate[ , names(Input), with = FALSE])
# Date Cycle
#1: 10 320
Or a join as @Frank mentioned
Input[!FDate, on=.(Date)]
# Date Cycle
#1: 10 320
In the OP's code,
setdiff(FDate$Date,Input$Date)
the first argument is from the 'Date' column from 'FDate' All of the elements in that column is also in the master data 'Input$Date'. So, it returns integer(0)). If we do the reverse, it would return 10
Finding overlapping ranges between two interval data
In general, it's very appropriate to use the bioconductor package IRanges to deal with problems related to intervals. It does so efficiently by implementing interval tree. GenomicRanges is another package that builds on top of IRanges, specifically for handling, well, "Genomic Ranges".
require(GenomicRanges)
gr1 = with(dtFrags, GRanges(Rle(factor(chr,
levels=c("1", "2", "X", "Y"))), IRanges(start, end)))
gr2 = with(dtCoords, GRanges(Rle(factor(chr,
levels=c("1", "2", "X", "Y"))), IRanges(coord, coord)))
olaps = findOverlaps(gr2, gr1)
dtCoords[, grp := seq_len(nrow(dtCoords))]
dtFrags[subjectHits(olaps), grp := queryHits(olaps)]
setkey(dtCoords, grp)
setkey(dtFrags, grp)
dtFrags[, list(grp, id, type)][dtCoords]
grp id type id.1 chr coord
1: 1 1 exon 10 1 150
2: 2 2 intron 20 2 300
3: 2 4 exon 20 2 300
4: 3 NA NA 30 Y 500
Related Topics
Saving Grid.Arrange() Plot to File
Rstudio Rmarkdown: Both Portrait and Landscape Layout in a Single PDF
Add Multiple Columns to R Data.Table in One Function Call
Run a for Loop in Parallel in R
Remove Rows from Data Frame Where a Row Matches a String
How to Flatten a List of Lists
Align Ggplot2 Plots Vertically
Reverse Order of Discrete Y Axis in Ggplot2
How to Read Only Lines That Fulfil a Condition from a CSV into R
Proper Idiom for Adding Zero Count Rows in Tidyr/Dplyr
Re-Ordering Factor Levels in Data Frame
How to Create a Loop That Includes Both a Code Chunk and Text with Knitr in R
R Convert Zipcode or Lat/Long to County
If/Else Constructs Inside and Outside Functions
R Split Numeric Vector at Position
Subsetting Data.Table by 2Nd Column Only of a 2 Column Key, Using Binary Search Not Vector Scan