bschilder
commented
2 years ago @NathanSkene Got it! The key was going back to the papers and understanding how each peak caller filters their peaks.
I've written a function (plot_precision_recall, and its subfunction precision_recall) which can generate precision-recall plots (and F1 plots as a bonus) for any combination of peakfiles and reference files.
It works with peaks called with any of the following: MACS2, MACS3, HOMER, and SEACR
In short, it works by iterating through increasingly stringent thresholds (the exact thresholds can be adjusted by users with several function args) and computing percentiles for the relevant metrics in each GRanges object, then computes precision-recall between the peakfiles vs. reference files (all combinations). I also parallelized some steps to speed things up further for many comparisons.
NOTE: It cannot compute precision-recall for thresholds that are less stringent that the one used to call the peaks originally. This is because those peaks that didn't reach the peak caller's threshold criterion were dropped from the files entirely (thus we don't have access to them in EpiCompare). Of course, anyone can always go back and rerun their peak calling with the most lenient thresholds and then rerun plot_precision_recall across the full range of thresholds.
By default, SEACR uses a threshold of 0.5 when you don't supply an IgG control file. Basically, it just takes the top 50% strongest peaks:
...a numeric threshold n between 0 and 1 returns the top n fraction of peaks based on total signal within peaks.
In light of this, it makes a lot of sense that our SEACR peaks have a max recall of almost exactly 50%. I suspect if we used a less stringent threshold when calling peaks (or perhaps used an IgG control) we would see greater recall (perhaps at the cost of some precision). Much like the plots in the initial post in this Issue indicate.
Example
data("CnR_H3K27ac")
data("CnT_H3K27ac")
data("encode_H3K27ac")
peakfiles <- list(CnR_H3K27ac=CnR_H3K27ac, CnT_H3K27ac=CnT_H3K27ac)
reference <- list("encode_H3K27ac" = encode_H3K27ac)
pr_out <- plot_precision_recall(peakfiles = peakfiles,
reference = reference)
In addition to returning both plots in a list, it also returns the processed data used to make the plots:
NathanSkene
For peak files where the necessary info is included, we can make precision/recall plots like this, using the reference as the "ground truth":
@leylaabbasova could you point us to where you wrote the code for this?