Hi,

I have a set of bins in a GrangesList, and a set of regions I would like to mask. My goal is to get a list of vectors, where each vector indicates the fractions of the bins not overlapped with the masked regions, so I can use these vectors for further scaling other things.

I've implemented a small function to achieve my goal. Although it did work on small testing data, it won't work in real and much larger data. Below is the testing data and the function I wrote to illustrate my goal.

The set of bins in a GrangesList, bin size is set as 200bp here.

grp_bin1 <-
    GenomicRanges::GRanges(
        seqnames = S4Vectors::Rle(c("1"), c(4)),
        ranges = IRanges::IRanges(start = seq(from = 1, by = 200, length.out = 4),
                                  end = seq(from = 200, by = 200, length.out = 4)),
        strand = S4Vectors::Rle(BiocGenerics::strand(c("+")), c(4)))

grp_bin2 <-
    GenomicRanges::GRanges(
        seqnames = S4Vectors::Rle(c("2"), c(4)),
        ranges = IRanges::IRanges(start = seq(from = 1, by = 200, length.out = 4),
                                  end = seq(from = 200, by = 200, length.out = 4)),
        strand = S4Vectors::Rle(BiocGenerics::strand(c("-")), c(4)))

grp_bins <-
    GenomicRanges::GRangesList("bin1" = grp_bin1, "bin2" = grp_bin2)

> grp_bins
GRangesList object of length 2:
$bin1
GRanges object with 4 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]        1     1-200      +
  [2]        1   201-400      +
  [3]        1   401-600      +
  [4]        1   601-800      +
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

$bin2
GRanges object with 4 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]        2     1-200      -
  [2]        2   201-400      -
  [3]        2   401-600      -
  [4]        2   601-800      -
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

The set of regions I would like to mask

add_mask1 <-
    GenomicRanges::GRanges(
        seqnames = S4Vectors::Rle(c("1", "2"), c(2, 1)),
        ranges = IRanges::IRanges(start = c(51, 251, 151),
                                  end = c(200, 400, 250)),
        strand = S4Vectors::Rle(BiocGenerics::strand(c("*")), c(3)))

> add_mask1
GRanges object with 3 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]        1    51-200      *
  [2]        1   251-400      *
  [3]        2   151-250      *
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

Then function I wrote and the expected output

scale_additional_masks <- function(bins, add_mask, bin_size) {
    # Find overlap
    ovr <- IRanges::findOverlapPairs(add_mask, bins)
    # Compute the percentages of the bins overlap with the masks
    ovr_intersect <-
        IRanges::pintersect(ovr, drop.nohit.ranges = FALSE)
    ovr_val <- GenomicRanges::width(ovr_intersect) / bin_size
    # Convert the overlap fractions into vectors
    ovr_val_ls <- base::split(ovr_val, names(ovr_val))
    add_scale <- lapply(ovr_val_ls, function(x) {
        1 - colSums(as.matrix(x))
    })
    return(add_scale)
}

> scale_additional_masks(bins = grp_bins, add_mask = add_mask1, bin_size = 200)
$bin1
[1] 0.25 0.25 1.00 1.00

$bin2
[1] 0.75 0.75 1.00 1.00

Everthing works fine until I move on to my real data, which has ~18K groups of bins, and ~7M non-overlapped regions need to be masked.

Browse[2]> bins
GRangesList object of length 18101:
$`10000_-`
GRanges object with 250 ranges and 0 metadata columns:
        seqnames        ranges strand
           <Rle>     <IRanges>  <Rle>
    [1]       11 644233-644482      -
    [2]       11 644483-644732      -
    [3]       11 644733-644982      -
    [4]       11 644983-645232      -
    [5]       11 645233-645482      -
    ...      ...           ...    ...
  [246]       11 705483-705732      -
  [247]       11 705733-705982      -
  [248]       11 705983-706232      -
  [249]       11 706233-706482      -
  [250]       11 706483-706732      -
  -------
  seqinfo: 23 sequences from an unspecified genome; no seqlengths

...
<18100 more elements>

Browse[2]> add_mask
GRanges object with 7238684 ranges and 0 metadata columns:
            seqnames              ranges strand
               <Rle>           <IRanges>  <Rle>
        [1]        1            1-382433      *
        [2]        1       382512-382513      *
        [3]        1       382529-382530      *
        [4]        1       382631-611597      *
        [5]        1       611655-611656      *
        ...      ...                 ...    ...
  [7238680]        X 156010147-156010167      *
  [7238681]        X           156010554      *
  [7238682]        X 156011799-156011812      *
  [7238683]        X 156012646-156012649      *
  [7238684]        X 156013711-156040895      *
  -------
  seqinfo: 23 sequences from an unspecified genome; no seqlengths

When I run and debug the function, it pops an error at this line.

ovr <- IRanges::findOverlapPairs(add_mask, bins)

 Error in METHOD(x, i) : 
  Subsetting operation on CompressedGRangesList object 'x' produces a result that is too big to be
  represented as a CompressedList object. Please try to coerce 'x' to a SimpleList object first (with 'as(x,
  "SimpleList")').

I cannot convert my mask regions to a SimpleList as suggested by the error

> as(add_mask1,
+   "SimpleList")
Error in getListElement(x, i, ...) : 
  GRanges objects don't support [[, as.list(), lapply(), or unlist() at the moment

I also tried to truncated the mask regions, and it seems working ok when the number of rows is around 3e6, but no more than 4e6. When it runs at this scale, it takes more than 100 GB RAM so I have to run it on a server rather than my local computer.

I'm wondering if there is any solutions could help me when data come to this scale? We have high memory node on server (1.5TB) so it might be ok to use high memory for a short time, but I guess it's also my implementation too naive so it make things so costly.

Any suggestions are highly appreciated!

Best,

Yixin

You're generating a lot (billions) of overlaps. A more practical approach is to compute the positions of the genome that is not covered by any masks, and then compute the fraction of each bin consisting of those positions.

Something like:

all_bins <- stack(grp_bins, "grp")
cov <- coverage(add_mask1)
views <- Views(cov == 0L, all_bins)
viewSums(views) / width(views)

But please be aware that it's critical to specify the seqinfo representing the full genome (or whatever your universe is) on the masks before calling coverage().