Error trying to use multiple file as input in Snakemake

Viewed 34

I'm trying to use some files as input for a bioinformatics tool. My basic code follows:

configfile: "config.yaml"


WORK_DATA = config["WORK_DATA"]

rule all:
  input:
    expand(WORK_DATA + "{sample}_{read}_quality-pass.fastq", sample=config["samples"], read=[1, 2])
     

rule filter:
  input:
    R1 = expand(WORK_DATA + "{sample}_1.fastq", sample=config["samples"]),
    R2 = expand(WORK_DATA + "{sample}_2.fastq", sample=config["samples"])
  
  output:
    expand(WORK_DATA + "{sample}_{read}_quality-pass.fastq", sample=config["samples"], read=[1, 2])
 
  params:
    outname1 = expand("{sample}_1", sample=config["samples"]),
    outname2 = expand("{sample}_2", sample=config["samples"])

  shell:
    "FilterSeq.py quality -s {input.R1} -q 20 --outname {params.outname1} --log FS1.log\n"
    "FilterSeq.py quality -s {input.R2} -q 20 --outname {params.outname2} --log FS2.log"

However, I get an error looking something like this:

Error in rule filter:
    jobid: 1
    output: /home/path/SRR1383456_1_quality-pass.fastq, /home/path/SRR1383456_2_quality-pass.fastq, /home/path/SRR1383457_1_quality-pass.fastq, /home/path/SRR1383457_2_quality-pass.fastq
    shell:
FilterSeq.py quality -s /home/path/SRR1383456_1.fastq /home/path/SRR1383457_1.fastq -q 20 --outname SRR1383456_1 SRR1383457_1 --log FS1.log
FilterSeq.py quality -s /home/path/SRR1383456_2.fastq /home/path/SRR1383457_2.fastq -q 20 --outname SRR1383456_2 SRR1383457_2 --log FS2.log
        (one of the commands exited with non-zero exit code; note that snakemake uses bash strict mode!)

Shutting down, this might take some time.
Exiting because a job execution failed. Look above for error message
Complete log: .snakemake/log/2022-04-28T113611.757898.snakemake.log

I don't know exactly the cause of this error, but it seems that those files are being placed as input (input.R1/R2) at the same time. Consequently, the FilterSeq.py program get's confused since it accepts only one input per run. I wonder how could I solve this issue, so that I could run all my intended files smoothly. Any help on this will be appreciated!

1 Answers

I think the issue lies in the multiple uses of expand. In the all rule, it looks good. It will evaluate to a list that looks something like

[sampleA_1_quality-pass.fastq, sampleA_2_quality-pass.fastq,
sampleB_1_quality-pass.fastq, sampleB_2_quality-pass.fastq, ....]

These should be a list of the files you expect as output. The general idea behind snakemake is that you can define abstract rules that rely on wildcard to "figure out" when they need to run to produce a desired output. In your filter rule however, the use of expand is explicitly defining the required input and output based on any file that matches in the filesystem.

This could create a problem because when the filter rule is trying to produce output for sampleA_1_quality-pass.fastq, it could be also attempting to include sampleB files as input or output during the expand step. You can see this in the command that ultimately run that includes two different samples here --outname SRR1383456_1 SRR1383457_1. I'm not sure what your FilterSeq.py program is doing but I believe that could be the problem.

Consider the following suggestion:

configfile: "config.yaml"

WORK_DATA = config["WORK_DATA"]

rule all:
  input:
    expand(WORK_DATA + "{sample}_{read}_quality-pass.fastq", sample=config["samples"], read=[1, 2])
     
rule filter:
  input:
    WORK_DATA+"{sample}_{read}.fastq"      
  output:
    WORK_DATA+"{sample}_{read}_quality-pass.fastq"
  params:
    outname="{sample}_{read}"
  shell:
    "FilterSeq.py quality -s {input} -q 20 --outname {params.outname} --log FS{wildcards.read}.log"

This way, the filter rule will run once for every sample/read combination specified in the all rule. As a final note, you did not include WORK_DATA as a prefix to your params so I didn't here as well. But I suspect you might need to in order to have the output appear where it is expected to.

Related