Conditional execution of multiplexed analysis with snakemake

Question

I've some troubles with Snakemake, up to now I didn’t found pertinent informations in the documentation (or somewhere else). In fact, I've a big file with different samples (multiplexed analyses) and I would like to stop the execution of the pipeline for some sample according to result found after rules.

I've already tried to change this value out of a rule definition (using a checkpoint or a def), to make conditional input for folowing rules and to considere wildcards as a simple list to delete one item. Below is an example of what I want to do (the conditional if is only indicative here) :

# Import the config file(s)
configfile: "../PATH/configfile.yaml"

# Wildcards
sample = config["SAMPLE"]
lauch = config["LAUCH"]

# Rules

rule all:
    input:
        expand("PATH_TO_OUTPUT/{lauch}.{sample}.output", lauch=lauch, sample=sample)


rule one:
    input:
        "PATH_TO_INPUT/{lauch}.{sample}.input"
    output:
        temp("PATH_TO_OUTPUT/{lauch}.{sample}.output.tmp")
    shell:
        """
        somescript.sh {input} {output}
        """

rule two:
    input:
        "PATH_TO_OUTPUT/{lauch}.{sample}.output.tmp"
    output:
        "PATH_TO_OUTPUT/{lauch}.{sample}.output"
    shell:
        """
        somecheckpoint.sh {input}       # Print a message and write in the log file for now

        if [ file_dont_pass_checkpoint ]; then
            # Delete the correspondant sample to the wildcard {sample}
            # to continu the analysis only with samples who are pass the validation
        fi


        somescript2.sh {input} {output}
        """

If someone has an idea I'm interested. Thank you in advance for your answers.

Answer 1

I think this is an interesting situation if I understand it correctly. If a sample passes some checks, then keep analysing it. Otherwise, stop early.

At the end of the pipeline, every sample must have a PATH_TO_OUTPUT/{lauch}.{sample}.output since this what the rule all asks for regardless of the check results.

You could have the rule(s) performing the checks writing a file containing a flag indicating whether for that sample the checks passed or not (say flag PASS or FAIL). Then according to that flag, the rule(s) doing the analysis either go for the full analysis (if PASS) or write an empty file (or whathever) if the flag is FAIL. Here's the gist:

rule all:
    input:
        expand('{sample}.output', sample= samples),

rule checker:
    input:
        '{sample}.input',
    output:
        '{sample}.check',
    shell:
        r"""
        if [ some_check_is_ok ]
        then
            echo "PASS" > {output}
        else
            echo "FAIL" > {output}
        fi
        """

rule do_analysis:
    input:
        chk= '{sample}.check',
        smp= '{sample}.input',
    output:
        '{sample}.output',
    shell:
        r"""
        if [ {input.chk} contains "PASS"]:
            do_long_analysis.sh {input.smp} > {output}
        else:
            > {output} # Do nothing: empty file
        """

If you don't want to see the failed, empty output files at all, you could use the onsuccess directive to get rid of them at the end of the pipeline:

onsuccess:
    for x in expand('{sample}.output', sample= samples):
        if os.path.getsize(x) == 0:
            print('Removing failed sample %s' % x)
            os.remove(x)

Answer 2

The canonical solution to problems like this is to use checkpoints. Consider the following example:

import pandas as pd

def get_results(wildcards):
    qc = pd.read_csv(checkpoints.qc.get().output[0].open(), sep="\t")
    return expand(
        "results/processed/{sample}.txt", 
        sample=qc[qc["some-qc-criterion"] > config["qc-threshold"]]["sample"]
    )


rule all:
    input:
        get_results


checkpoint qc:
    input:
        expand("results/preprocessed/{sample}.txt", sample=config["samples"])
    output:
        "results/qc.tsv"
    shell:
        "perfom-qc {input} > {output}"


rule process:
    input:
        "results/preprocessed/{sample}.txt"
    output:
        "results/processed/{sample.txt}"
    shell:
        "process {input} > {output}"

The idea is the following: at some point in your pipeline, after some (let's say) preprocessing, you add a checkpoint rule, which aggregates over all samples and generates some kind of QC table. Then, downstream of that, there is a rule that aggregates over samples (e.g. the rule all, or some other aggregation inside of the workflow). Let's say in that aggregation you only want to consider samples that pass the QC. For that, you let the required files ("results/processed/{sample}.txt") be determined via an input function, which reads the QC table generated by the checkpoint rule. Snakemake's checkpoint mechanism ensures that this input function is evaluated after the checkpoint has been executed, so that you can actually read the table results and base your decision about the samples on the qc criteria contained in that table. Any intermediate rules (like here the process rule) will then be automatically applied by Snakemake when re-evaluating the DAG.