Best practices for handling single-end and paired-end data in a Snakemake pipeline

Question

I am fairly new to Snakemake and am building an NGS data processing pipeline to align ChIP-seq data and call peaks on it, plus some other analyses. The data can be single-end or paired-end, depending on the experiment.

Originally, I was intending to handle single-end and paired-end data using a sample_sp wildcard for the tools that require different paramters, and impose a regex constraint on it. That way, all my paired-end data files could bear a _pe suffix, and would only be processed by the paired-end rules in the pipeline like so:

rule bowtie2_map_se:
    input:
        "DATA/FASTQ/{sample_sp}.filtered.fq.gz"
    output:
        temp("DATA/BAM/UNSORTED/{sample_sp}.sam")
    ...


rule bowtie2_map_pe:
    input:
        fq1="DATA/FASTQ/{sample_sp}_fw.filtered.fq.gz",
        fq2="DATA/FASTQ/{sample_sp}_rv.filtered.fq.gz"
    output:
        temp("DATA/BAM/UNSORTED/{sample_sp}_pe.sam")
    ...

By constraining sample_sp to not end with {"_fw", "_rv", "_pe"}, then snakemake will not attempt to use those data files for single-end specific rules.

However, I encountered issues when trying to implement these regex constraints.

What are the best practices for handling single-end and paired-end data in your snakemake pipelines? Are there ways to implement these regex constraints, or am I approaching the issue from the wrong point of view?

Answer 1

The way I solve that is with a metaData table where I define what is what. This is not exacly what you ask for but, I find it more flexible since the pipeline does not crash if a file is named differently.

metaData.csv

path;sample_sp;pairedStatus
input/A.fq.gz;A;SE
input/B-fw.fq.gz;B;PE
input/B-rv.fq.gz;B;PE

snakefile

import pandas as pd
samples = pd.read_table("metaData.csv",sep=';').set_index("path", drop=False)


rule all:
    input:
        expand("out/{sample_sp}.sam",sample_sp=list(set([f"{a}_{b}" for a, b in zip(samples.sample_sp.iloc, samples.pairedStatus)])) )


rule test_PE:
    input:
        lambda wc: samples[(samples.sample_sp == wc.id) & (samples.pairedStatus == "PE")].path
    output:
        temp("out/{id}_PE.sam")
    shell:
        "test"

rule test_SE:
    input:
        lambda wc: samples[(samples.sample_sp == wc.id) & (samples.pairedStatus == "SE")].path
    output:
        temp("out/{id}_SE.sam")
    shell:
        "test"

Answer 2

Not sure about best practice, but this is basically also how i would do it. Having two independent rules handling PE and SE separately. So my guess would be that the error stems from a later rule which combines both PE and SE outputs again, or associating the right input to the correct rule.

Could you provide a minimal example, this might help to spot the error.

Answer 3

The negative lookbehind does seem to work for me. With just the below in a test Snakefile:

wildcard_constraints:
    sample_sp = r".+(?<!_pe)"

rule bowtie2_map_se:
    input:
        "DATA/FASTQ/{sample_sp}.filtered.fq.gz"
    output:
        temp("DATA/BAM/UNSORTED/{sample_sp}.sam")

rule bowtie2_map_pe:
    input:
        fq1="DATA/FASTQ/{sample_sp}_fw.filtered.fq.gz",
        fq2="DATA/FASTQ/{sample_sp}_rv.filtered.fq.gz"
    output:
        temp("DATA/BAM/UNSORTED/{sample_sp}_pe.sam")

Then...

$ snakemake -n DATA/BAM/UNSORTED/foo_pe.sam
Building DAG of jobs...
MissingInputException in rule bowtie2_map_pe in file /home/tbooth2/workspace/snakemake_scratch/stack_overflow/Snakefile, line 10:
Missing input files for rule bowtie2_map_pe:
    output: DATA/BAM/UNSORTED/foo_pe.sam
    wildcards: sample_sp=foo
    affected files:
        DATA/FASTQ/foo_fw.filtered.fq.gz
        DATA/FASTQ/foo_rv.filtered.fq.gz

If you comment out the wildcard constraints, the error from Snakemake will indicates that it considers both rules to be candidates to generate this output. I could touch the input files to get rid of this error, or use --debug-dag, but there is no point as the error message tells me everything I need.

Remember that wildcard_constraints have no effect on rule inputs - they are only used by Snakemake when trying to determine if a given rule can create some requested output.

While this can work, I'd not call it "best practice". My inclination would be to have the data path reflect the type of BAM being created, eg:

rule bowtie2_map_se:
   input: ...
   output:
        temp("DATA/BAM/UNSORTED_SE/{sample_sp}.sam")

rule bowtie2_map_pe:
    input: ...
    output:
        temp("DATA/BAM/UNSORTED_PE/{sample_sp}_pe.sam")

Now whatever rule deals with the unsorted BAM files (presumably using samtools or picard to sort the BAM) is going to have to do an if/else to decide if it wants to get input from UNSORTED_PE or UNSORTED_SE, so you'll need an input function there, but this can implement whatever logic you like, based on the sample name, or on a metadata table as suggested by @kEks, or just a global config item if you are not mixing se/pe in one run.