Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

R, join within a range vectorised

Tags:

merge

r

dplyr

sqldf

I'm trying to join two datasets where a variable (or position along a genome) in one dataset fits within a range in the second (gene start/stop position). However, positions are not unique, but nested within an additional column (chromosome). The same goes for the gene start/stop positions. My goal is to link each position with the corresponding annotation and effect.

For example:

library(sqldf)
set.seed(100)
a <- data.frame(
    annotation = sample(c("this", "that", "other"), 3, replace=TRUE),
    start = seq(1, 30, 10),
    chr = sample(1:3, 3, replace=TRUE)
  )
a$stop <- a$start + 10
b <- data.frame(
    chr = sample(1:3, 3, replace=TRUE),
    position = sample(1:15, 3, replace=TRUE),
    effect = sample(c("high", "low"), 3, replace=TRUE)
  )

An SQL inner join gets me part of the way there:

df<-sqldf("SELECT a.start, a.stop, a.annotation, b.effect, b.position
    FROM a, b
    inner JOIN a b on(b.position >= a.start and b.position <= a.stop);")

But this doesn't account for the repetition of position per chromosome. I'm having conceptual trouble wrapping this into a loop or apply function.

I'm not wedded to SQL, it's just how I tackled a simpler problem previously. I'm also not sure that making an additional index column is appropriate since I have thousands of chromosome values.

My desired output would look like the following:

    df$chr<-c("NA","2","2")
      start stop annotation effect position chr
1     1   11       this   high        3  NA
2     1   11       this   high       10  NA
3    11   21       this    low       14   2

Where each position has been placed between the start and stop points on the correct chr, or given NA where no points on a chr match.

like image 557
smm Avatar asked Aug 01 '26 22:08

smm


1 Answers

The development version of data.table introduces non-equi joins, allowing for:

library(data.table)
setDT(a) # converting to data.table in place
setDT(b)

b[a, on = .(position >= start, position <= stop), nomatch = 0,
  .(start, stop, annotation, effect, x.position, chr = ifelse(i.chr == x.chr, i.chr, NA))]
#   start stop annotation effect x.position chr
#1:     1   11       this   high          3  NA
#2:     1   11       this   high         10  NA
#3:    11   21       this    low         14   2
like image 83
eddi Avatar answered Aug 04 '26 16:08

eddi