saddles in cis by distance - flexible way

[sergpolly]

intro: saddles plots discrete/continuous by distance are useful and informative ...

"By distance" means - taking into account interactions within a certain distance range, e.g. 1MB<D<5MB when calculating "the saddle" ...

Currently, the "by distance" part is done by filling anything that does not belong to a selected distance band - i.e. outside of min_diag and max_diag with NaNs:

# snippet from inside the _accumulate subfunction of saddle module

# matrix is an array with HiC data in it
# cis, corresponding to 1 chromosome or an arm
for d in np.arange(-min_diag + 1, min_diag):
    numutils.set_diag(matrix, np.nan, d)
if max_diag >= 0:
    for d in np.append(
        np.arange(-matrix.shape[0], -max_diag),
        np.arange(max_diag + 1, matrix.shape[0]),
    ):
    numutils.set_diag(matrix, np.nan, d)

Instead one can generate saddles for every diagonal without significant performance penalty:

So, instead of existing _accumulate one could use something like that:

def _accumulate_dist( S, C, getmatrix, digitized, reg1, reg2, max_dist ):
    """
    Helper function to aggregate across region pairs. Does it by distance !
    max_dist is the bin_id of 
    """
   # S, C now have to be a 3-dimensional arrays of n_bins X n_bins X max_dist
   # S(D,n,n) - first dim correspond to a certain distance (in bins) and the latter 2 are the typical saddles n*n
   # C(D,n,n) - same thing here as well ...

    n_bins = S.shape[-1]
    # should work only in cis - i.e. reg1 or reg2, or inter-arm
    matrix = getmatrix(reg1, reg2)

    # Toeplitz matrix with distance for every pixel ...
    _dist_vec = np.arange(matrix.shape[0])
    dist_mat = np.abs(_dist_vec[None,:] - _dist_vec[:,None])

    for i in range(n_bins):
        row_mask = digitized[reg1] == i
        for j in range(n_bins):
            col_mask = digitized[reg2] == j
            data = matrix[row_mask, :][:, col_mask]
            dist = dist_mat[row_mask, :][:, col_mask]
            is_finite_mask = np.isfinite(data)
            data = data[is_finite_mask]
            dist = dist[is_finite_mask]
            # unleash the power of bincounter here :
            # S unrolled by distances - inefficient memory access - isn't it ?
            S[:, i, j] += np.bincount(dist, weights=data, minlength=max_dist)
            # C unrolled by distances
            C[:, i, j] += np.bincount(dist, minlength=max_dist).astype(float)

This way one would accumulate saddles (sum and counts first) into 3D-stacked arrays, ala pileups/snipping: S[diagonals, saddle_bins, saddle_bins] - such that it would be easy to do by distance saddles like so:

saddle_mat_1to5MB = S[bin_1MB:bin_5Mb+1, :, :].sum(axis=0)/C[bin_1MB:bin_5Mb+1, :, :].sum(axis=0)
saddle_mat_5to25MB = S[bin_5MB:bin_25Mb, :, :].sum(axis=0)/C[bin_5MB:bin_25Mb, :, :].sum(axis=0)

Full prototype solution below:
https://gist.github.com/sergpolly/ee39a452c1e30f12d5100b28f35f4ee0

[sergpolly]

There is a step missing in the provided solution - symmetrising the final interaction matrix (saddle) sums + sums.T ...
So for the 3D sums and count arrays this can be achieved like so:

# transposing 2nd and 3rd dimensions, leaving 1st dimension alone !
_sX = _sumX + _sumX.transpose(0,2,1)
_cX = _countX + _countX.transpose(0,2,1)

it is unclear though why that symmetrising is there in the first place and when does it come into play ...

Potentially useful for future fully sparse implementation of saddle

[golobor]

related PR: #484
@nvictus @gfudenberg @Phlya @agalitsyna

[golobor]

to me, this seems like a very useful extension of the saddle functionality!
Though, i do have a question - how heavy could the final objects be, potentially?

A simple back-of-the-envelope calculation:

• in human genomes, max separation is s_max = ~250Mb.
• at 1Kb, resolution, max_diag = 2.5e8 / 1e3 = 2.5e5 diagonals
• if digitized has ~50 values, each matrix has 50x50=2.5e3 elements.
• thus, if aggregated across all chromosomes, the final object will contain ~6.3e8 elements, which is manageable, although borderline.
• splitting into chromosomes or chromosomal arms would push this towards 1e10 elements, which would be impossible to keep in memory at once and difficult to analyse.
  Using 10kb resolution would only reduce the memory requirements to 1e9, which is still borderline.

Are these estimates correct?

I can think of two solutions to mitigate these issues:

• use distance binning. This makes the API and algorithms a bit uglier, but would greatly help with performance.
• use out-of-core computations, combined with some kind of binary format for saddle storage (hdf5, zarr, etc). Could be done too, but may potentially be computationally costly due to I/O.

[sergpolly]

Are these estimates correct?

100%, yes ! this would explode in memory very very quickly ... definitely .

Yes, this could become another way of shooting yourself in the foot - but not the first one:

• pairwise pileups in trans are pretty heavy - would be enabled by Refactor snipping to enable trans, inter-arm, etc pileups #475
• eigs_trans is still dense and explodes in memory as well
• something else ? - dense matrix extraction from cooler itself can go >>1e10 easily

But enough of that ...

Practical solutions:

1. "use distance binning." is of course practical, but ...
   • would require those 3D stacks to turn into labeled arrays/dataframe/xarrays for interpretation of distance/diagonal intervals
   • OR can also supply user-defined distance bins dist_bins=[0, 1000, 5000, 20000], this is logical extension for existing min_dist/max_dist
2. storing intermediate stacks into npz/hdf5 - also practical - that's what I do when I'm actually encountering "heavy" intermediate arrays that are a pain to recalculate
   • the only problem is - not sure how to implement in a clean way inside API
   • maybe, like, each _saddle_per_region- function (the _accumulate replacement) could store each 3D stack into npz/hdf5 ?
3. MY favorite solution ! cause - lazy and need by_distance
   • Enable all of the plumbing internally for by region dictionary of giant 3D stacks, BUT:
   • keep existing user saddle-API in saddle - still return it only in the aggregated form - i.e. 2D array for all regions combined
   • add unrestricted saddle_by_region_by_dist API into sandbox - see how often is this used, refine API and ultimately see if (1) or (2) is the way

[sergpolly]

such ultimately flexible saddle_by_dist is also a flashback to existing (in sandbox now) contrast_by_diagonal function by Johannes https://github.com/open2c/cooltools/blob/master/cooltools/sandbox/contrast.py -> there is a lot going on there, but I believe they are exactly the same conceptually

[sergpolly]

some examples just in case - here is how a saddle by distance looks like for WT-like and deltaRad21-like samples - biology aside - short range does look quite different (saddles are done using their own eigen vectors @25kb binsize, intra-arm only)

so, yes, - we don't need all of those individual diagonals per se - but they are nice to have - to aggregate them into distance bins after the fact

I personally haven't even check different chroms separately - that was purely theoretical excercise - I think still it would be useful