Macs peak calling

[immanuelazn]

Description

Add macs2/macs3 peak calling Utilize a multi-threaded approach separated by cluster within prep_macs_inputs(). Allow users to choose a specific step to run (ie "prep-inputs", "run-macs", "read-outputs") to allow for executing macs calls in external slurm cluster.
Also change behaviour in InsertionIterator to not keep insertions in memory after moving to next fragment.

Tests

• Added test cases for prep_macs() parameterized by cluster existence, threading, and insertion_mode

[bnprks]

Here's some initial comments I came up with from a read-through. Overall looks like a pretty good structure and approach. I'm focusing on some higher-level comments here.

There are also several small things that I think could be tweaked around choices of naming, docs, default arguments, organization etc. Rather than dictate solutions to those while you're still drafting, I'll let you finish up and do your own polishing pass before we discuss those.

C++ comments:

• I'm pretty sure that writeInsertionBedByPseudobulk should write out duplicate coordinates rather than explicitly skipping printing repeats
• Once we get a chance to do a performance profile, we should check if gzprintf is taking up a bunch of time (and comparing with/without a .gz file ending to isolate if the compression or the printf part is slow)
  • We might also want to take a quick look to see if macs3 can support any other kinds of compressed input that would be faster to write
• We will definitely need to debug the memory issue you identified in InsertionIterator as well

R comments

• With the option to split up execution using step, I'm a bit concerned about errors that might happen if users edit the input folder between steps, or try calling a later step with arguments that are mismatched/incompatible with their earlier step
  • This might mean making cell_names non-nullable, and maybe even fragments, though that might not be worth it given the usability trade-off.
• prep-inputs should probably write the macs command to a .sh file (one-per-cluster). Also returning the commands is as a vector is fine, but later steps could then read from the .sh files to simplify
  • Also we usually have a few more arguments present, including e.g. --shift. You might want to check the ArchR source a bit more for that as I think arguments are listed in at least two separate places
• For mclapply we probably also want to always use the argument mc.preschedule=FALSE

[immanuelazn]

There are a few more tests to be added that directly compare writing insertions from fragments in ArchR vs BPCells. I might just copy the output files to eliminate the dependence of building an entire ArchR project.
Otherwise, I think I've addressed most of your comments, and the feature is in a relatively polished state.

[bnprks]

Looking pretty good! Most remaining comments are style-related, though I did end up with a large number of those in call_macs_peaks(). Trying out a bulleted list here in rough order of the code rather than using Github's inline comment feature since it's a bit easier for me to prep.

Good job setting up the priority queue InsertionIterator -- I'm fine leaving it as a stopgap for this PR but I'd really like to get the radix sort algorithm fixed up and running soon after we merge this.

InsertionIterator

• Very cool that you were able to compare performance between the priority queue and the sorting version!
• Minor style note but *(starts + i) is the same as starts[i]

writeInsertionBedByPseudobulk()

• Given that this writes out just one pseudobulk, I'm not sure why both the cells and group_of_interest arguments are needed. Arguably, you could even just pass in fragments to already be subset to just the cells of interest from the R side with select_cells()
  • At least filtering prior to InsertionIterator will save some sorting work
• Unclear if keep_dups argument is needed (more discussion below). This also might eliminate the need for the cleanup at end of chromosome section
• Error messages shouldn't be prefixed with writeBedgraph: anymore
• I'd suggest just calling this function just writeInsertionBed for simplicity. As a cleanup task I can edit writeInsertionBedgraph to just write one pseudobulk at a time anyhow as I think that will make the most sense
• Either this file name will need an update or this function should go into its own file because it doesn't handle Bedgraph files. (I'd be happy to rename to maybe BedWriter since bedgraph is kinda just another type of bed file)

prep_macs_inputs()

• Does this need to be a separate exported function? Seems like it could mostly slot in to call_macs_peaks or at least be private

call_macs_peaks()

1. Make docs and parameter names consistent with call_peaks_tile (cell_groups and effective_genome_size arguments). Can use @inheritParams to save on re-typing docs if the text will be identical.
   • You might also consider copying the default initialization of cell_groups, though I'm not sure how that would work if we allow fragments to be NULL
2. path docs. Put directory structure in @details section. Specify inputs are <path>/input/<cluster>.bed(.gz) and shell commands are <path>/input/<cluster>.sh
3. insertion_mode docs don't quite read naturally. If you want consistency with write_insertion_bedgraph I'd update both to say maybe "for coverage calculation" rather than "for insertion counts calculation"
4. Do we need the keep_dups option? I was just looking through the ArchR code and I don't think it would ever de-duplicate (and as mentioned previously I think it's likely undesirable to remove duplicate insertion sites in a peak calling context). Though if your testing indicates otherwise we can discuss
5. step argument should default to all
6. macs_version could probably be macs_executable instead, and allow the user to pass a full path if desired. (And I'd say default to macs3 at this point, or make the default NULL and try to auto-detect macs2 vs macs3 with test system() calls)
   • Either way, we probably want to check if macs can be run before prepping any inputs in the case of step = all
7. additionalParams should be snake case like you have it in the docs. It also might be nice if the name made it clear that these are parameters getting passed to MACS
   • Still don't have a --shift argument being passed, which might make sense to add here
8. quiet: If we multi-thread run-macs, then we might want to write the macs stdout/stderr to a log file in outputs/ and have this just print when we start writing each input file and when we start running each copy of macs (maybe also printing when each finishes?). Also, for consistency this should flip and be named verbose
9. use_gz docs would be helpful if they could give some absolute file sizes + time to write with /without compression for a small-ish example dataset (e.g. one of the public 10x pbmc datasets)
10. threads argument in call_macs_peaks should probably also be used in the run-macs section
11. Return value should be a single data frame, with an extra column specifying the group that the cluster was called in. Try to place the group column early in the ordering (probably either first or right before name?). A list of data frames is a little awkward to deal with as output
12. Will the system2 call work properly when there are spaces in cluster names? Might want to add some quoting to macs_call_template. Also maybe get rid of any / character in a cluster names?
13. In docs, bedfile should be split up as two words for style

tests

• Looks reasonable to me. I like your tests that the insertion bed file has the right contents. If prep_macs_inputs stops being a public/separate function, you could adjust to just call call_macs_peaks(step="prep-inputs") in that test

[immanuelazn]

I think I'm happy with the current state, but please let me know if you see any more glaring changes I should be making. It will be a little bit of a journey attempting to extract insertion beds from ArchR, while removing the hardwired pseudobulk replicates and tn5 bias steps. To make sure I can continue progressing, I suggest that I table creating a testcase from comparing outputs from ArchR. Getting ArchR to skip those steps will require a decent amount of software engineering time by itself. Would that be reasonable?

[immanuelazn]

Ok it should be ready! I added in the test comparing Archr and BPCells results. The only change required is setting devtools::load_all() to target your bpcells/archr install.

[bnprks]

Overall, I like how this is going and fundamentally it's quite close to being ready. A lot of the details I asked for from the last round are fixed. I like the ordering you chose for the call_macs_peaks parameters, and making write_insertion_bed have a matching interface to write_insertion_bedgraph. The performance testing and comparison against ArchR outputs I imagine were tricky and time-consuming, but good job on them.

There are several more detail-oriented issues remaining though, several from my last round of suggestions which were either forgotten or only partially addressed. I've also included a few new ones that came up as I was able to try out the function more.

Sorry for having taken so much of the week to get back to you on this, I hope it hasn't been too much of a blocker. I'd be happy to do a call early next week to discuss any questions/clarification about this round of comments.

Follow-ups still to be addressed from last comment round

writeInsertionBed

• After removing the keep_dups argument, there's no reason to have the variable last_base and hence the end-of-chromosome cleanup step.

call_macs_peaks()

• path docs in details section.
  • For input files, I'd describe contents/purpose of the files. E.g just "shell commands" doesn't say what they are for. "Bed files" doesn't explain what the contents are.
  • For outputs, for the files BPCells doesn't use, we can just link to the macs3 documentation on output files.
  • Very useful to say we only read the narrowPeaks file, but you need an extra line break so it doesn't get swallowed into the last bullet when rendered as a webpage
  • In the docs for the path argument, you should add a (see Details for directory structure) or equivalent message
• macs_version was not renamed in the function declaration, only in the docs. Also, using match.arg means that it's not possible for the user to pass a custom path for a macs executable
• use_gz docs would be helpful to give numbers for the time trade-off as well. Most straightforward would be to give both time + input folder size for with and without compression.
• threads argument in call_macs_peaks should probably also be used in the run-macs section
• Return value "cluster" column would ideally not be the last column
• Will the system2 call work properly when there are spaces in cluster names? Might want to add some quoting to macs_call_template. Also maybe get rid of any / character in a cluster names?

New comments:

Verbose logging setup:

• log_message variable will not be defined in "run-macs" section if run standalone (and I also don't think it's necessary to include that information in the log file)

• In write_insertion_bed(), the logging should be inside the mclapply call so it prints at the correct time. (Despite the theoretical race condition, I think this basically always works out fine even in very tight loops)

  • The style of log you use in run-macs is better also, rather than printing out cluster names on individual lines

  • Example of the current (sub-optimal) output

    [1] "2024-08-23 16:10:15.508477 Writing bed file for clusters:"
    [1] "consonant"
    [1] "vowel"
    [1] "2024-08-23 16:10:15.532154 Finished writing bed files"
    

• I'd suggest format(Sys.time()) to avoid having the fractional seconds get printed

• I think using message() rather than print() will be better to avoid the [1] prefixes getting printed

run-macs and error checking:

• macs_executable now has a slightly weird combination of both auto-detect and a default value.

  I'd say the best option is to default NULL, and in the case of NULL we can do your auto-detect procedure. If not null, we just check that macs can be found else we error out

  • Also, the existing stop_run and log_message variables can probably be eliminated since I don't think there's a need / easy way to hang on to log_message for use in later steps.

• The existing check on macs_executable should check the output of macs3 --version to make sure it contains the text macs -- otherwise the check will pass for almost any executable even if it's not macs

• In run-macs, it's probably better for system2 to write stdout/stderr directly to the log file so that logs are written as the command is running rather than all at the end. It's OK if this prevents writing the command as it appears MACS already prints the input command at the top of the log file

• Would be good to check in run-macs that the command ran without error (had a 0 exit code)

Other:

• The test for macs_e2e_works now introduces a system dependency on having macs2 installed. Can we make it skip instead if macs is not installed?
• --call-summits probably can go in --additional-params, as it seems to just affect the peak calling method used by macs while leaving the output file format consistent
• in write_insertion_bed, cell_groups should probably come after path for consistency with write_insertion_bedgraph, and cell_groups should maybe default to NULL for consistency.

Questions/discussion areas (can do via call/slack)

• Is there a better name possible for use_gz? I'm just realizing the name doesn't really give any hints about what is getting compressed, so it would be nice to have a name that would be more specific to the macs bed insertion inputs. Let me know if you have a suggestion here
• What do you think of the adjusting column names/values in the output? Reading the [UCSC narrowPeak format docs](https://genome.ucsc.edu/FAQ/FAQformat.html#format12), I think it's kinda weird that pValue is given as -log10(p), and peak/pointSource is given relative to start. On the one hand it's nice to maintain a direct correspondence with the macs output, but some parts of it are kind of non-intuitive