ericsdavis / hictoolsr Goto Github PK

When using extractCounts() to pull from a bedpe file (read in as a data.frame and converted to a GInteractions object with as_ginteractions()), the resulting object contained no rows.

This turned out to be an issue with the chromosome names - in the loop bedpe file used they were listed as 1, 2, etc. However this portion of the extractCounts() function requires that they be in "chr1", "chr2" format:

## Split binned binnedEP by chromosome
  chrBedpe <- lapply(chroms, function(chr) {
    bedpe[seqnames(first(bedpe)) == paste0('chr',chr)]
  })

Simply pasting "chr" before the chromosomes in the original bedpe object fixed this issue!

binBedpe is used to assign anchors of a GInteractions object to a resolution-binned range. If this function is used to simultaneously expand anchor regions then keep in mind that the expansion will not be equal on both sides since it is then reassigning the range to a bin after the resize operation. Therefore, successive calls to binBedpe on the same dataset (either to shrink or expand the resolution range) will propagate error from the original ranges.

For example:

GInteractions(anchor1 = GRanges("chr1:5-6"),
              anchor2 = GRanges("chr1:15-16")) |>
    binBedpe(res = 5, a1Pos = 'center', a2Pos = 'center')

results in :

GInteractions object with 1 interaction and 0 metadata columns:
      seqnames1   ranges1     seqnames2   ranges2
          <Rle> <IRanges>         <Rle> <IRanges>
  [1]      chr1      5-10 ---      chr1     15-20
  -------
  regions: 2 ranges and 0 metadata columns
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

This is the correct 5-length bin from the center of each anchor range. However, if you try to work backwards to obtain the original range, you will end up with this:

>   GInteractions(anchor1 = GRanges("chr1:5-6"),
+                 anchor2 = GRanges("chr1:15-16")) |>
+     binBedpe(res = 5, a1Pos = 'center', a2Pos = 'center') |>
+     binBedpe(res = 1, a1Pos = 'center', a2Pos = 'center')
GInteractions object with 1 interaction and 0 metadata columns:
      seqnames1   ranges1     seqnames2   ranges2
          <Rle> <IRanges>         <Rle> <IRanges>
  [1]      chr1       7-8 ---      chr1     17-18
  -------
  regions: 2 ranges and 0 metadata columns
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

For this example you could set a1Pos and a2Pos to 'start', but this will not always result in the correct shifted range. Whenever you need to resize the range, be sure to work with the original GInteractions object.

Currently, the matrix is always set to "observed". I would like to extract observed/expected counts. Could you add a parameter to the extractCounts function so that it can be set to either "observed" or "oe"?

Add normalization argument to extractCounts()