Gracefully handle variable number of input and output files

user_1114453 0
Viewed 250

You have a pipeline where some samples have one input file and produce one output file while other samples have two input files and produce two output files. The typical case for those in bioinformatics is NGS sequencing libraries where some samples are paired-end and other samples are single-end sequenced. If you need to trim reads and align them you have to account for the variable number of input/output files.

What is the most appropriate way to handle this? I feel using checkpoints is overkill since in my opinion checkpoints are a bit cryptic, but I may be wrong...

Here's how I would do it for the case where the number of input/output of files can be only 1 or 2: Have an if-else in the run directive based first on the number of input files. If the number of inputs is 1 touch the second files. For the following rules have again an if-else this time checking whether the second file has 0 bytes.

Here's an example (not tested but it should be about right):

import os

samples = {'S1': ['s1.R1.fastq'], 
           'S2': ['s2.R1.fastq', 's2.R2.fastq']}

rule all:
    input:
        expand('bam/{sample_id}.bam', sample_id= samples),

rule cutadapt:
    input:
        fastq= lambda wc: samples[wc.sample_id],
    output:
        r1= 'cutadapt/{sample_id}.R1.fastq',
        r2= touch('cutadapt/{sample_id}.R2.fastq'),
    run:
        if len(input.fastq) == 1:
            # {output.r1} created, {output.r2} touched
            shell('cutadapt ... -o {output.r1} {input.fastq}')
        else:
            shell('cutadapt ... -o {output.r1} -p {output.r2} {input.fastq}')

rule align:
    input:
        r1= 'cutadapt/{sample_id}.R1.fastq',
        r2= 'cutadapt/{sample_id}.R2.fastq',
    output:
        'bam/{sample_id}.bam',
    run:
        if os.path.getsize(input.r2) == 0:
        # or 
        # if len(samples[wildcards.sample_id]) == 1:
            shell('hisat2 ... {input.r1} > {output.bam}')
        else:
            shell('hisat2 ... -1 {input.r1} -2 {input.r2} > {output.bam}')

This works, but I find it awkward to artificially create the second file just to keep the workflow happy. Are there better solutions?

snakemake python
2 Answers

There shall be two separate rules for single ended and pair ended case:

rule single_end:
    input:
        fastq = 's{n}.R1.fastq'
    output:
        r1 = 'cutadapt/S{n}.R1.fastq'
    shell:
        'cutadapt ... -o {output.r1} {input.fastq}'

rule pair_ends:
    input:
        r1 = 's{n}.R1.fastq',
        r2 = 's{n}.R2.fastq'
    output:
        r1 = 'cutadapt/S{n}.R1.fastq',
        r2 = 'cutadapt/S{n}.R2.fastq'
    shell:
        'cutadapt ... -o {output.r1} -p {output.r2} {input}'

Now there is an ambiguity: whenever the pair_ends can be applied, the single_end may be applied too. This problem can be solved with can be ruleorder:

ruleorder pair_ends > single_end
user_112186870

I agree that checkpoints would be overkill here.

You can compose arbitrary inputs and command lines using input functions and parameter functions. The trick is optional r2 output from your cutadapt rule, which may or may not be present, which complicates the workflow's DAG. In this case it's probably appropriate to use Snakemake's directories as output functionality:

CUTADAPT_R1 = 'cutadapt/{sample_id}/R1.fastq'
CUTADAPT_R2 = 'cutadapt/{sample_id}/R2.fastq'

def get_cutadapt_outargs(wldc, input):
    r1 = CUTADAPT_R1.format(**wldc)
    r2 = CUTADAPT_R2.format(**wldc)
    ret = f"-o '{r1}'"
    if len(input.fastq) > 1:
        ret += f" -p '{r2}'"
    return ret

def get_align_inargs(wldc):
    r1 = CUTADAPT_R1.format(**wldc)
    r2 = CUTADAPT_R2.format(**wldc)
    if len(samples[wldc.sample_id]) > 1:
        return "-1 '{}' -2 '{}'".format(r1, r2)
    else:
        return "'{}'".format(r1)

rule all:
    input:
        expand('bam/{sample_id}.bam', sample_id= samples),

rule cutadapt:
    input:
        fastq= lambda wc: samples[wc.sample_id],
    output:
        fq = directory('cutadapt/{sample_id}'),
    params:
        out_args = get_cutadapt_outargs,
    shell:
        """
        cutadapt ... {params.out_args} {input.fastq}
        """

rule align:
    input:
        fq = rules.cutadapt.output.fq,
    output:
        bam = 'bam/{sample_id}.bam',
    params:
        input_args = get_align_inargs,
    shell:
        """
        hisat2 ... {params.input_args} > {output.bam}
        """
user_32037190
Related