I have a set of `GenomicRanges` and all I need is to expand it by certain distance on start and another distance on end (taking into account the strand). For example, from this set:

expanding by 20 nt at start and by 30 nt at the end, and considering chromosome length as 350, I want to get this:

Note that 1) for range #2 shift values are inversed because it's negative strand, 2) for range #3 resulting end position is 350 because it can't be longer than chromosome length.

Here is the code to generate this sample data:

mygr <- GRanges(seqnames  = Rle(rep('chr1', 3)),
               ranges    = IRanges(c(100, 200, 300), c(110, 220, 330)),
               strand    = c('+','-','+'),
               seqlengths= c(chr1=350))

So far I wrote the following function:

grexpand <- function(inputGR, befLen=0, aftLen=0){
  inputGR <- resize(inputGR,  width = width(inputGR)+befLen, fix = 'end')
  inputGR <- resize(inputGR,  width = width(inputGR)+aftLen, fix = 'start')
  return(trim(inputGR));
}

it does the job but I'm wondering if this functionality is already included in the package and I just don't know the right way to call it...

Hi,

AFAIK we don't have a function in GenomicRanges to perform this particular transformation. FWIW here is an alternative to your grexpand():

extend <- function(x, upstream=0, downstream=0)     
{
    if (any(strand(x) == "*"))
        warning("'*' ranges were treated as '+'")
    on_plus <- strand(x) == "+" | strand(x) == "*"
    new_start <- start(x) - ifelse(on_plus, upstream, downstream)
    new_end <- end(x) + ifelse(on_plus, downstream, upstream)
    ranges(x) <- IRanges(new_start, new_end)
    trim(x)
}

It avoids a collision you can get with grexpand() between the start and end of the intermediate GRanges object (i.e. the one produced by your first call to resize()) when befLen is negative:

grexpand(mygr, -12, 12)
# Error in .local(x, width, fix, use.names, ...) : 
#   'width' values must be non-negative

extend(mygr, -12, 12)
# GRanges object with 3 ranges and 0 metadata columns:
#       seqnames     ranges strand
#          <Rle>  <IRanges>  <Rle>
#   [1]     chr1 [112, 122]      +
#   [2]     chr1 [188, 208]      -
#   [3]     chr1 [312, 342]      +
#   -------
#   seqinfo: 1 sequence from an unspecified genome

More generally speaking, when downstream == -upstream, the transformation should be equivalent to a strand aware shift() and so should always work.

Cheers,

H.