The peakCounts function in PEPATACr.R should probably need fixing

[zhongzheng1999]

When I used looper run and looper runp to obtain a consensus peak set, I specified this consensus peak set as frip_ref_peaks and ran looper run andlooper runp again. Confusingly, looper runp ignored the looper run results generated based on the consensus peak set, failing to produce a Count table with a PEPATAC-produced consensus peak set. Examining the log files, I identified the issue with PEPATACr. I discovered that in the peakcounts function:

    # check if coverage files are compressed
    if (any(file.exists(file.path(results_subdir,
                        sample_names, paste0("peak_calling_", genomes),
                        paste0(sample_names, "_ref_peaks_coverage.bed.gz"))))) {
        ext <- ".bed.gz"

These lines verify the existence of *_ref_peaks_coverage.bed.gz files, but the output files from looper run based on the consensus peak set are named *_ref_peaks_coverage.bed.

Therefore, when I modified the code:

if (any(file.exists(file.path(results_subdir,
                        sample_names, paste0("peak_calling_", genomes),
                        paste0(sample_names, "_ref_peaks_coverage.bed"))))) {
        ext <- ".bed"

It successfully produced the Count table with a consensus peak set generated by PEPATAC.

I believe I need to confirm with you whether the result obtained after the second run of looper run should be _ref_peaks_coverage.bed or _ref_peaks_coverage.bed.gz? The reason for my concern is that after the second looper run, I found temporary "tmp" files in each sample folder. As a result, I cannot determine if the issue lies with PEPATACr or if there is an incompleteness in my pipeline that is preventing the final Count Table generation.

[zhongzheng1999]

This issue has been reported before and I suspect it may be for similar reasons! #219