cli wrapper for expected

[sergpolly]

Simple CLI wrapper for expected.py:

cooltools compute_expected --cis in.cool
cooltools compute_expected --trans in.cool

Using bunch of cores:

cooltools compute_expected -n 8 --cis in.cool
cooltools compute_expected -n 8 --trans in.cool (would raise NotImplemented)

NOTES:

Wrapper works fine with the latest dask(0.17.0) and distributed(1.21.0).
Only stdout output is implemented, just like for diamond_insulation.

TODO:

• Should probably add ignore_diags as a click option - will do later.
• Improve my usage of click options - i provided links to click features I tried to use, but help looks ugly a bit.
• do we have to explicitly shutdown distributed cluster ?

[sergpolly]

Forgot to mention: the wrapper computes expected only for ALL chromosomes (as provided by cooler.chromnames) both for --cis and for --trans, so we don't exploit the flexibility of expected.py to specify smaller sub-chromosome-regions (--cis) or a subset of chromosomes (--trans).

I feel like it's good enough for CLI wrapper - it's going to be used for the entire genome in 99% of cases.

[Phlya]

What do you mean by all chromosomes - there will be an expected value for each distance for each chromosome, or one value for the whole genome?
If it's the latter, I think excluding chrM and chrY, for example, would be very useful! I saw some artefacts on scaling curves computed from expected values coming from them - e.g. a peak at the end of chrY.

[sergpolly]

In the case of --cis there will be single output that would contain expected for all chromosomes, one by one, one after another:

chrom	diag	n_valid	balanced.sum	balanced.avg
chr1	0	223	nan	nan
chr1	1	218	nan	nan
chr1	2	217	6.334600193162759	0.029191705959275387
chr1	3	215	3.179269773687505	0.01478730127296514
...
chr1    247     1       0.0003275742729092722   0.0003275742729092722
chr1    248     0       0.0     nan
chr1    249     0       0.0     nan
chr2    0       237     nan     nan
chr2    1       233     nan     nan
chr2    2       232     7.965345378195444       0.03433338525084243
chr2    3       231     3.906494527809641       0.016911231722119658
...

So there will be no option (for now at least) to compute expected only for chr10 for example, or for an even more limited region, like chr10:0-10000000. expected.py allows you to do all of that, but do we need it in a CLI wrapper ?!

In a --trans case, output would look like that:

chrom1  chrom2  balanced        n_valid
chr1    chr10   7.356826602356368       33811
chr1    chr11   7.344436681470667       33811
chr1    chr12   7.482985420556106       33419
chr1    chr13   5.300837695437331       28460
...
chr9    chrY    0.6505189901154699      7233
chrX    chrM    0.04497483626805778     156
chrX    chrY    2.821101613527302       8970
chrY    chrM    0.00488750315282725     60

So, again , you'd get an expected for ALL pairs of chromosomes, without an option to compute expected for a subset in inter-chromosomal regions, say ['chr5', 'chr6', 'chr7'] - which is again possible by using expected.py in your custom scripts.

Does that sound reasonable, or do you think we need to add some of that flexibility to the initial version of the wrapper ?

[sergpolly]

Another thing we should do right away, is to figure out naming:
for --cis: balanced.sum is the total number of pairs (normalized, aka balanced) and balanced.avg is the actual expected balanced.avg = balanced.sum / n_valid. Nomenclature seems ok for me .

Should we do the same for --trans case as well ?! rename balanced to balanced.sum, because that's what it is, and add another column balanced.avg = balanced.sum / n_valid, so that everything seems consistent ?

Other naming preferences ?

[Phlya]

Awesome, thanks!

For what would one need the .sum values?

[sergpolly]

they we there in the beginning ... - so history as the reason N1 ))
it's nice sometimes to see n_valid and total signal , i.e. balanced.sum, and not just the ratio.

For --trans case , there is another thing: in order to get a single trans-expected for the entire genome, you'd need to sum(balanced.sum)/sum(n_valid), which is not the same as sum( balanced.sum/n_valid ).

[sergpolly]

(1) Changed interface from

cooltools compute_expected --cis in.cool
cooltools compute_expected --trans in.cool

to

cooltools compute_expected --contact-type cis in.cool
cooltools compute_expected --contact-type trans in.cool

correspondingly. Seems more intuitive after all, and help prints out better.

(2) added --drop-diags flag to control # of diags ignored during cis-calculations

(3) changed the expected.py itself a bit, to make trans_expected return a dataframe with column-names more in line with those of cis_expected:
chrom1 chrom2 balanced.sum n_valid balanced.avg
Meaning is preserved.

(4) tested - works fine, no stupid mistakes

[nvictus]

probably good practice, but we can leave it commented for now

[nvictus]

Looks good so far!

[nvictus]

The next thing to do is to provide an API to do finer partitions in cis than full chromosomes, i.e. chromosome arms. Will open a new issue.

[nvictus]

Replying to some of the discussion above, note that having the .sum and n_valid values in the table allows you to generate a genome-wide expected via groupby.

[nvictus]

hmm...this suggests that if someone has launched their own dask scheduler as a daemon process, we could provide an option to pass in the <host>:<port> to the scheduler instead of starting a new cluster ourselves.