Getting a stretch of genomic ranges from a dataframe/granges object based on metadata column

Question

I have a "test" data.frame object in R, which is basically a small subset of a 66000 row dataframe, which looks as follows:

chr   start    end      PC1     comp
chr1  3360001  3400000  -0.009  B
chr1  3400001  3440000  -0.004  B
chr1  3440001  3480000  -0.001  B
chr1  3480001  3520000  0.003  A
chr1  3520001  3560000  0.002  A
chr1  3560001  3600000  -0.001  B
chr1  3600001  3640000  0.002  A
chr1  3640001  3680000  0.004  A
chr1  3720001  3760000  0.003  A
chr1  3760001  3800000  0.003  A
chr1  3800001  3840000  0.007  A
chr1  3840001  3880000  0.006  A
chr1  3880001  3920000  -0.004  B
chr1  3920001  3960000  -0.017  B

As you can see, the interval between each start and end position is 40kb. The column comp can have only two values, A or B. Going from row 1 to the last row of my dataframe and based on the values of the column comp, I want to get the whole stretch/block of the genomic region, once the compartment(character in column comp) is changed, and perform a mean of the PC1 values in that stretch.

So in this example, I basically want to shrink my test dataframe as follows:

chr1  3360001  3480000  mean(c(-0.009,-0.004,-0.001))  B
chr1  3480001  3560000  mean(c(0.003,0.002))  A
chr1  3560001  3600000  mean(c(-0.001))  B  <---- #notice this is a 'singleton'
chr1  3600001  3880000  mean(c(0.002,0.004,0.003,0.003,0.007,0.006))  A
chr1  3880001  3960000  mean(c(-0.004,-0.017))  B

I tried with a simple code, but it only runs until the first compartment value change(so row 4 in this test), and I am unable to figure out how to make it run for the whole dataset of 66000 rows. Also the use of multiple for loops, which is what I tried, may not be the prettiest way to do this.

I guess there exists a Bioconductor package which performs this task very easily. If anyone can point me to that direction, it'd be great.

Answer 1

The GenomicRanges package in bioconductor will do this. It lets you do 'set' operations on genomic ranges, such as concatenating them, finding overlaps, etc. Hopefully something in the following code will be inspiring, even if I haven't solved the exact problem for you here:

library(GenomicRanges)

# make a GenomicRanges set corresponding to the original table
test.gr = makeGRangesFromDataFrame(df = test, start.field = "start", end.field = "end", seqnames.field = "chr")

# make separate GenomicRanges sets for each of the two partitions of the genome
test_A.gr = makeGRangesFromDataFrame(df = test[test$comp=="A",], start.field = "start", end.field = "end", seqnames.field = "chr")
test_B.gr = makeGRangesFromDataFrame(df = test[test$comp=="B",], start.field = "start", end.field = "end", seqnames.field = "chr")

# apply the reduce() function to concatenate segments which are adjoining each other into single segments
reduced_test_A.gr = reduce(test_A.gr)
reduced_test_B.gr = reduce(test_B.gr)

# group PC1 values by segment classification
test[queryHits(findOverlaps(test.gr, reduced_test_A.gr)), "PC1"]
test[queryHits(findOverlaps(test.gr, reduced_test_B.gr)), "PC1"]

Answer 2

Take a look at the rle function:

Compute the lengths and values of runs of equal values in a vector – or the reverse operation.

In this way you can mark the moment in which A changes to B and the length of the streches of A and B.

With that I think you can reduce the number of for loops you have to use (I can not think a way of doing this without at least one loop, but for sure it may exist)

Answer 3

What is proposed by @Jonathan works, just adding another possible solution here.

library(GenomicRanges)

test.gr = makeGRangesFromDataFrame(test,keep.extra.columns=TRUE)
# as long as the previous comp does not agree with current, we make it a new segment
test.gr$seg = c(0,cumsum(diff(as.numeric(test.gr$comp)) != 0))

# simple function to collapse the range
fn = function(u){
  out = reduce(u)
  out$comp = unique(u$comp)
  out$mean = mean(u$PC1)
  out
}

reduced_range = unlist(GRangesList(lapply(split(test.gr,test.gr$seg),fn))))


```