improve distance binning for FR,RR,FR,RF pairs "scalings" in stats output

[sergpolly]

https://github.com/mirnylab/pairtools/blob/d1ddf9c39a336662f7fc725fa5a70ec68df9ba95/pairtools/pairtools_stats.py#L147

consider replacing it with something more readable and usable, e.g. @mimakaev 's robust bins:

# ~10 bins per order of magnitude
bins = np.logspace(0,9, num = 9*10+1,dtype=int)
bins = np.unique(bins)
bins = np.cumsum(np.sort(np.r_[1,np.diff(bins)]))

currently we have:

min_log10_dist=0
max_log10_dist=9
og10_dist_bin_step=0.25
bins = np.r_[0, np.round(10**np.arange(min_log10_dist, max_log10_dist+0.001, log10_dist_bin_step)).astype(np.int)]

which are also non-decreasing, but are too sparsely spaced ... - and code is hard to read

[sergpolly]

also should we expose number_of_bins_per_order_of_magnitude as a parameter in pairtools/distiller or , just stick with something

pairtools stats is more of a QC than the analysis tools at the end

[mimakaev]

So the problem is matching bins between cooler and pairtools. It's a little tricky. There are a few options.

• One is the current one: not to try to match them.
• Another is to match them at 1kb.
• A third is to match them at 100bp, 200bp, and 1kb, making 1kb a perfect resolution (and sacrificing 100,200 performance). (100,500,1000) is also an option. Not both (500/200 is 2.5).
• A fourth is to do a partial match (that pairs -> cooler bins are reduseable).

Current proposal for nice cooler bins: 1,2,3,4,5,6,8,10,13,16,20,25,32,40,50,63,79,100,

at 1kb those bins become: 1000,2000,3000,4000,5000,6000,8000,10000,13000,16000,20000,25000,32000,40000,50000,63000,79000...
at 100bp those bins become
100,200,300,400,500,600,800,1000,1300,1600,2000,2500,3200,4000,5000,6300,7900,10000,12600,15800,20000,25100,31600,39800,50100,63100,79400,100000,125900,158500,199500,251200,316200,398100,501200,631000,794300,1000000,1258900,1584900,1995300,2511900,3162300,3981100,5011900,6309600,7943300

The two are clearly different

pairtools bins not matched to cooler:
...126,158,200,251,316,398,501,631,794,1000,1259,1585,1995,2512,3162,3981,5012,6310,7943,10000,12589,15849,19953,25119,31623,

pairtools= bins matched to cooler at 1kb:
1,2,3,4,5,6,8,10,13,16,20,25,32,40,50,63,79,100,126,158,200,251,316,398,501,631,794,1000,2000,3000,4000,5000,6000,8000,10000,13000,16000,20000,25000,32000,40000,50000,63000,

pairtools bis matched to cooler at 1kb and 100bp/200bp if 100bp/200bp cooler uses modified bins:
1,2,3,4,5,6,8,10,13,16,20,25,32,40,50,63,79,100,200,400,600,1000,2000,3000,4000,5000,6000,8000,10000,13000,16000,20000,25000,32000,40000,50000,63000

And bins matched to cooler at 100, 200, 1000bp resolutions with extra bins for pairs.
1,2,3,4,5,6,8,10,13,16,20,25,32,40,50,63,79,100,126,159,200,252,318,400,490,600,693,800,1000,1260,_1590,_2000,3000,4000,5000,6000,6930,8000,10000,13000,16000,20000,25000,32000,40000,50000,63000 (bin in italics are the ones that were added)

I couldn't think of a more general solution. Powers of two obviously has one, but not here...

[mimakaev]

@golobor @sergpolly @agalitsyna - what do you guys think? We should probably decide on this before we merge in cooltools logbin_expected.

-- Should we aim at matching at one resolution, or at two, or matching at all?
-- Would there ever be a need in sub-kb scalings? 100bp resolution?

[golobor]

alternatively, we can let users decide between two options: (a) keep bins nice within all orders of magnitude, (b) do not use nice bins at all.

[sergpolly]

IMHO - ~100 bp is needed , at least for pair-level stuff because of DNase/MNase-based methods like microC, OmniC, and whateverC might happen "tomorrow"

100bp coolers for microC isn't a crazy thing to do, so perhaps it makes sense to match it like @mimakaev suggested:

And bins matched to cooler at 100, 200, 1000bp resolutions with extra bins for pairs.
1,2,3,4,5,6,8,10,13,16,20,25,32,40,50,63,79,100,126,159,200,252,318,400,490,600,693,800,1000,1260,_1590,_2000,3000,4000,5000,6000,6930,8000,10000,13000,16000,20000,25000,32000,40000,50000,63000 (bin in italics are the ones that were added)

but this would only work for high-resolution coolers and wouldn't be applicable to sparse data - <50-100M pairs of usable pairs in a cooler. So like @golobor is suggesting - this matching between bins for coolers and pairs could be optional

another IMHO - i don't think it is THAT crucial to match bins for pairtools stats with coolers, as it is just for QC, and it's not meant to be for used for "real" analysis (by-arm, individual chromosomes, etc etc)

[mimakaev]

Yeah, that would probably be ideal. I will a little better engineer that set and make sure it is actually matched.

[mimakaev]

These are ratios of neighboring pair bins in the current version of bins.

bins = [10,13,16,20,25,32,40,50,63,79,100,126,159,200,240,300,400,490,600,800,1000,1200,1600,2000,2400,3000,4000,5000,6000,8000,10000,13000,16000,20000,25000,32000,40000,50000,63000]

Bins for 100bp and 200bp (just without 100 and 300)
100,200,300,400,600,800,1000,1200,1600,2000,2400,3000,4000,5000,6000,8000,10000,13000,16000,20000,25000,32000,40000,50000,63000

[golobor]

i'll repeat this, but - why not making bins nice in all orders of magnitude? ( [1,2,3,4,5,6,8,10] + [1,2,3,4,5,6,8,10] * 10 + [1,2,3,4,5,6,8,10] * 100 ) What negative consequences would this have?

…

On Fri, 6 Mar 2020 at 01:55, Maksim Imakaev ***@***.***> wrote: [image: image] <https://user-images.githubusercontent.com/9454715/76039728-a7af5080-5f23-11ea-9501-8f288256dfa4.png> These are ratios of neighboring pair bins in the current version of bins. bins = [10,13,16,20,25,32,40,50,63,79,100,126,159,200,240,300,400,490,600,800,1000,1200,1600,2000,2400,3000,4000,5000,6000,8000,10000,13000,16000,20000,25000,32000,40000,50000,63000] Bins for 100bp and 200bp (just without 100 and 300) 100,200,300,400,600,800,1000,1200,1600,2000,2400,3000,4000,5000,6000,8000,10000,13000,16000,20000,25000,32000,40000,50000,63000 — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#81>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAG64CRFPJEFYQW2JFOCFDDRGBCZ5ANCNFSM4KZJNMVQ> .

[mimakaev]

ok, now I get it.

A large negative consequence is a two-fold jump from 1 to 2. Could have used 1 2 5 10 instead - that's at least even.

A partial remedy is to use these bins, and drop #2,3,5 in the first order of magnitude
1000,1300,1600,2000,2400,3000,4000,5000,6000,8000 which has slightly better bin size ratios.

[agalitsyna]

I will convert this to the discussion for now, but feel free to comment or open an issue if binning improvements are needed!

[agalitsyna]

Since 1.0.0 we'll have pairtools scaling utility, with fairly readable output. Kudos to @golobor for proposing transfer of pairlib to pairtools.
#117