[discussion] snparcher CLI

[pdimens]

Hey @cademirch @tsackton, I saw on the v2 roadmap that you want to work towards a CLI for snparcher, which is a great idea. I'm trying to work towards a similar thing with our friends over at loco-pipe, which is conceptually similar to SNParcher with respect to one snakemake invocation doing a great many things/analyses. In the case of loco-pipe (a project I'm unaffiliated with, but am a friend and colleague of Nicolas'), I wanted the users at my uni to have some kind of setup automation (a frequent pain point), but without refactoring anything loco-pipe actually does. My WIP attempt is bundling loco-pipe (and its assets) as a minimal python package, then exposing loco-pipe init to setup a run directory, followed by loco-pipe start to do a basic directory validation and call snakemake on your behalf with all the configs present in the output dir from init, along with any user modifications since (hpc, etc.).

I think SNParcher's infrastructure would work well for that kind of design, i.e. snparcher init + snparcher start (or something similar), without compromising its strength of being a one-stop shop, in contrast to e.g. the discrete and narrow workflows in harpy, each with its own CLI. I'm still refining the loco-pipe approach, but you can see the barebones design here, where it uses a bit of simple python to copy loco-pipe configs/workflows, and fill in placeholder config text with format strings (stdlib and click, but argparse would suffice without adding extra deps). It would be ideal to run snakemake using the python API instead of as a subprocess, but I still can't wrap my head around it after the big API change between versions 7/8/9 and the inconsistent documentation associated with those changes.

So far, this CLI looks like (and is unlikely to become any more sophisticated):

➜ loco-pipe init   
Usage: loco-pipe init [OPTIONS] SAMPLES...

  Initialize a configured loco-pipe project

  Create a directory (--output) with all the components loco-pipe needs to
  run. Must include a `--reference` file. Provide the input BAM files and/or
  directories at the end of the command as individual files/folders, using
  shell wildcards (e.g. `data/*.bam`), or both.

Options:
  -o, --output TEXT     name of directory to create  [default: .]
  -r, --reference FILE  reference fasta file which BAMs were aligned to
                        [required]
  -s, --simple          Do not annotate the workflow YAML file

➜ loco-pipe init -o koyfish -r toyfish/reference/toy_refgen.fa -s toyfish/bams
Notices:
- Samples have all been assigned to a single group and need to have distinct classifications in koyfish/docs/samples.tsv.

➜ loco-pipe start  
Usage: loco-pipe start [OPTIONS] [DIRECTORY]

  Launch the loco-pipe snakemake workflow

  The target `directory` is expected to have been created with `loco-pipe init`.
  Make sure you have addressed any notices provided by the init command prior
  to launching loco-pipe. Will print snakemake output to the terminal.

Options:
  -@, --threads INTEGER RANGE  Number of threads to use  [default: 8;
                               1<=x<=999]
  -d, --dry                    Perform a snakemake dry run

[cademirch]

Hi @pdimens - thank you for sharing! Funnily enough I was just discussing an idea much like this with @tsackton at our last meeting. Happy to see you've already thought about this and implemented it! I really do like the idea of distributing a workflow this way - guiding the user setup their project directory and copying the workflow to that directory.

@tsackton what are your thoughts?

[pdimens]

@cademirch thanks, yeah it would be a bit of handholding where it matters most to at least get a massive pipeline like that started and running. Plus, the additional benefit of having a snparcher versioned conda (or container) install via bioconda, and a snparcher command associated with running the workflow.

[tsackton]

I think this is a good approach and it is nice to see your framework @pdimens for inspiration. We've also experimented with slightly less polished versions of this, e.g. in PhyloAcc we have a Python frontend that generates a Snakefile for users to run.

We are working to modularize snparcher a little, so that for example it is easier to start with just bams, or to run the QC step only, ultimately with the goal of exposing a variety of snparcher subcommands. But starting with an init command, maybe some kind of prep/download command that just fetches data, and then a run command makes sense to me.

[pdimens]

I see, I was unaware that you were considering this kind of modularization. If thats the case, I think the tooling in Harpy would be of value for that. Each workflow is a snakemake pipeline, and the commands essentially verify data and build out a properly configured directory and run snakemake from there. You're welcome to use as much of the existing harpy internals for that as you would like.

Overview

There are internal classes that facilitate common actions to reduce code duplication, these include:

• FASTQ that checks and validates FASTQ files
• SAM the same, but alignment files
• FASTA same ...
• VCF ...

Along with Workflow, which is probably the most useful of the bunch see the dev version here. The Workflow class is what takes in all the CLI arguments/options and turns them into an output folder, sets up the config.yaml, hpc.yaml, pulls scripts and reports, grabs the snakefile, and ultimately launches Snakemake. A realistic implementation of that in one of the subcommands looks like:

workflow = Workflow("qc", "qc.smk", output, container, clean, quiet)
workflow.setup_snakemake(threads, hpc, snakemake, no_temp)
workflow.notebook_files = ["qc_bx_stats.ipynb"]
workflow.conda = ["qc"]

## checks and validations....
fastq = FASTQ(inputs)

workflow.notebooks["skip"] = skip_reports
workflow.linkedreads["type"] = fastq.lr_type
workflow.input(fastq.files)
workflow.param(trim_adapters, "trim-adapters")
workflow.param(deduplicate, "deduplicate")
workflow.param(length[0], "min-len")
workflow.param(length[1], "max-len")
if extra_params:
    workflow.param(extra_params, "extra")
workflow.initialize(setup)

[tsackton]

Hmm, I'll have to look at this in more depth, thanks for sharing. We've been thinking about this as setting up Snakemake rules, e.g. if you want to just produce the callable sites output, you can use rule callable_sites that has as input only the bed and zarr files produced by that branch of the the DAG.

[pdimens]

I understand. In that case, you would be able to sneak in the name of the target non-wildcard rule as the final positional argument TARGET to snakemake (source)

targets

Targets to build. May be rules or files.

snakemake --cores ... callable_sites