3
$\begingroup$

I'm running a custom perl script using Snakemake. I use this rule:

rule complexity_20mer_counter:
    input:
       os.path.join(fastq_trimmed_dir, '{sample}_{read}' + fastq_extension.split('.')[0] + '_val_' + '{read}'.strip('R') + '.fq.gz')
    output:
        os.path.join(custom_qc_dir, '{sample}_{read}_20mer_counter.txt')
    shell:
        '''
        #!/bin/bash
        perl scripts/20mer_counter.pl {input} > {output}
        '''

It fails (I'll get to why later) with this error:

Building DAG of jobs...
MissingInputException in line 1 of /rst1/2017-0205_illuminaseq/scratch/swo-406/2019-01-24/rules/20mer_counter.smk:
Missing input files for rule complexity_20mer_counter:
fastq_trimmed/xs0067_P2018SEQE27R01_S1_R1_001_val_R1.fq.gz

... this is because I use the output of trim_galore, this rule looks like this:

rule trim_galore:
    input:
       [os.path.join(fastq_raw_dir, '{sample}_'+ read + fastq_extension)
       for read in reads]
    output:
        timmed_fastq_files = [os.path.join(fastq_trimmed_dir, '{sample}_' + read + fastq_extension.split('.')[0] + '_val_' + read.strip('R') + '.fq.gz')
        for read in reads],
        trimming_reports = [os.path.join(trim_galore_qc_dir, '{sample}_'+ read + fastq_extension + '_trimming_report.txt')
        for read in reads]
    conda:
        "../envs/fastqc.yaml"
    shell:
        '''
        #!/bin/bash
        trim_galore --fastqc --gzip -o {fastq_trimmed_dir} --paired --retain_unpaired {input}
        # Move all qc reports to the fastq_trimmed_dir
        mv {fastq_trimmed_dir}/*.zip \
           {fastq_trimmed_dir}/*.html \
           {fastq_trimmed_dir}/*.txt \
           {trim_galore_qc_dir}
        '''

As you can see, trim_galore not only uses my {read} parameter in the output filename but again adds it at the end, only now without the "R".

I tried solving it in my first rule by applying a python .strip('R') on my {read} wildcard, but this does not seem to work. How can I correctly get '1' and '2' from 'R1' and 'R2' (which are my {read} wildcards)?

(Extra bonus points for a nice output.trim_galore.something as input for my complexity_20mer_counter rule, which does not work because the output is a list (it needs 2 fastq files at the same time), but my complexity_20mer_counter rule can be run in parallel again.)

Improve this question
$\endgroup$
3
  • $\begingroup$ For me it's quite hard to give you a final solution, without having a look on the whole Snakefile or at least a runnable minimal example. I would suggest that the values for read are already without the R. Than you can easily decide when a R is needed in the filename and when not. $\endgroup$
    – finswimmer
    Commented Jan 31, 2019 at 13:55
  • 1
    $\begingroup$ This Snakefile is quite difficult to follow, and its complexity means that subtle mistakes might easily creep in. For example, in your trim_galore rule you have timmed_fastq_files instead of trimmed_fastq_files as the label for your output files. I don't know if this is a problem, but little things like that can definitely mess up the whole pipeline. Also, the shebang at the beginning of each "shell" directive is unnecessary. :-) $\endgroup$
    – Daniel Standage
    Commented Jan 31, 2019 at 14:13
  • $\begingroup$ Yes, I already caught the timmed error, it was not it. I did try putting reads to ['1', '2'] and making _R part of the string but may people use my pipelines and I want to make the default what every is used most. Doing what @finswimmer suggests will lead to a lot of code to make the "R" optional somehow. Although one could probably figure something out. I would just like to do string operations in the wildcard, like {read}.strip('R'). I'll try to look into making something minimal working but that is a lot of effort in this large pipeline. $\endgroup$
    – Freek
    Commented Jan 31, 2019 at 14:59

1 Answer 1

Reset to default
2
$\begingroup$

It seems that it is not possible to apply strip() directly on the wildcard placeholders. Luckily it is possible to executed functions in the input directive.

So try changing your input of the complexity_20mer_counter rule to this:

lambda wc: os.path.join(fastq_trimmed_dir, wc.sample + '_' + wc.read + fastq_extension.split('.')[0] + '_val_' + wc.read.strip('R') + '.fq.gz')

The wildcard object is passed to our anonymous function. So have access to the wildcards value via wc.sample and wc.read, which are now handled as normal strings and we can apply strip on them.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Yes, this is what I was looking for, it works, thank you! (I did call "wildcards", not "wc", not sure if that is supported, I call "wildcards" everywhere. $\endgroup$
    – Freek
    Commented Feb 4, 2019 at 13:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.