In 4.4 Peak-calling in the interaction space in the diffHiC users guide, it is suggested to use a "modest" minimum enrichment threshold of 0.5. It is also suggested that most parent bins should have low numbers (<5) of nested bin pairs, and that larger proportions require more aggressive filtering. A more aggressive minimum enrichment of 0 is suggested in the 2nd example.

While the user guide has suggested that a lower threshold is more stringent, instinctually, I would expect a higher minimum enrichment to be more stringent than a lower one, so for my own data set I looked at several thresholds above and below 0.5. At a resolution of 20kb, a lower minimum enrichment actually does result in an increase in nested bin pairs, contrary to what is shown in the user guide example. For the most part, increasing the min enrichment seems to lower the number of nested bins pairs, though some do increase depending on the threshold. However at lower resolutions, the number of bin pairs seems to go up or down depending on how far I go in either direction, but even decreasing the threshold to an extreme negative number cannot produce bin numbers of <5 for most parent bins. It's a bit unclear whether I should be increasing or decreasing the min enrichment threshold or what might be a good value to use.

How exactly is this minimum enrichment threshold used in peak calling? Does the sign (+ vs -) affect the results? What is the best way to determine what threshold to use?

Below are some examples of different thresholds, using the code suggested in the users guide. (I've copied it only once below -the only thing I changed for the output was the min.enrich value or the data set since I used two different resolutions).
Thank you, to anyone who might have some insight.

> peak.keep_1e6_stringent <- filterPeaks(original_1e6, enrichments_1e6, min.enrich=0.5, min.count=5, min.diag=2L)
> sum(peak.keep_1e6_stringent)
[1] 27464
> neighborhood_1e6_stringent <- (2*flank.width + 1) * metadata(original_1e6)$width
> boxed_1e6_stringent <- boxPairs(data[peak.keep_1e6_stringent], reference=neighborhood_1e6_stringent)
> out_1e6_stringent <- tabulate(tabulate(boxed_1e6_stringent$indices[[1]]))
> setNames(round(c(out_1e6_stringent[1:5], sum(out_1e6_stringent[5:length(out_1e6_stringent)]))/sum(out_1e6_stringent)*100, 1), c(1:5, ">5"))

1e6 min.enrich = 0.5:

1e6 min.enrich= 0.8 (this seems to be the most reasonable threshold, with most parent values <5):

1e6 min.enrich = 0 (numbers change marginally when the enrichment score is increasingly negative, all the way down to -9):

20kb min.enrich = 0.5 (only 3 parents are <5, but perhaps that's enough?):

20kb min.enrich = 0:

20kb min.enrich = 1 (increasing this:

20kb min.enrich = -1:

For good measure, the summary of enrichment scores of the data at both resolutions:

> summary(enrichments_1e6)

> summary(enrichments_20kb)

There might be some misunderstanding with the wording in the user's guide:

Larger proportions suggest that more aggressive filtering on the enrichment score is required, e.g., in the second example.

What I meant was that there are larger "> 5" proportions in the second example, so more aggressive filtering is required in that situation by increasing min.enrich. Your results seem to be fairly consistent with what is expected, i.e., increasing min.enrich increases the isolation of bin pairs corresponding to putative peaks, resulting in a decrease in the proportion of parents containing many internal bin pairs.

I guess I could change the wording to be more explicit. In any case, we're working on a more statistically rigorous method of defining the "peakness" cut-off, which will hopefully supercede the need to fiddle with min.enrich. (Though having said that, it would require some other threshold, e.g., based on the FDR, so you'll have to fiddle with something.)