0
$\begingroup$

I'm trying to run a command using snakemake and a TSV config file. Here is my sample config file:

sample_id       FQ1_id  FQ2_id
MyCleanName-1       MyFqName-1_R1.fastq.gz        MyFqName-1_R2.fastq.gz
MyCleanName-2       MyFqName-2_R1.fastq.gz        MyFqName-2_R2.fastq.gz

Here is my snakefile:

import pandas as pd
samples = pd.read_table("samples.tsv").set_index("sample_id", drop = False)

rule all:
    input:
        expand("{sample}_{num}.fastq.gz", sample = samples.sample_id, num = [1,2])

rule my_job:
    input:
        fq1 = expand("{fq1_id}", fq1_id = samples.FQ1_id),
        fq2 = expand("{fq2_id}", fq2_id = samples.FQ2_id),
        logdir = expand("{sample}.logs", sample = samples.sample_id)
    output:
        "{sample}_1.fastq.gz",
        "{sample}_2.fastq.gz",
    message:
        "Running MyCmd on {wildcards.sample}"
    shell:
        """
        MyCmd --sample-name {wildcards.sample} -i1 {input.fq1} -i2 {input.fq2}
        """


rule create_logdir:
    output:
        directory("{sample}.logs")
    message:
        "Creating log directories"
    shell:
        "mkdir -pv {wildcards.sample}.logs/{{MyCmd,MyOtherCmd}}"

When I run snakemake -np MyCleanName-1.fastq.gz, I get the following:

Building DAG of jobs...
Job counts:
        count   jobs
        2       create_logdir
        1       my_job
        3

[Fri Dec 24 16:38:54 2021]
Job 1: Creating log directories

mkdir -pv MyCleanName-1.logs/{MyCmd,MyOtherCmd}

[Fri Dec 24 16:38:54 2021]
Job 2: Creating log directories

mkdir -pv MyCleanName-2.logs/{MyCmd,MyOtherCmd}

[Fri Dec 24 16:38:54 2021]
Job 0: Running MyCmd on MDA-IBC-3


        MyCmd         --sample-name MyCleanName-1 -i1 MyFqName-1_R1.fastq.gz MyFqName-2_R1.fastq.gz -i2 MyFqName-1_R2.fastq.gz MyFqName-2_R2.fastq.gz

Job counts:
        count   jobs
        2       create_logdir
        1       my_job
        3
This was a dry-run (flag -n). The order of jobs does not reflect the order of execution.

How can I restrict this so the run only creates log files for the Clean Name being used and uses only the corresponding FASTQs? I feel like I'm missing something trivial here. Should I not be using the config file variable samples directly in the rule my_job, but somehow wiring up the corresponding sample_id and FQ_ids on the all rule?

$\endgroup$
3
  • $\begingroup$ expand is a Snakemake help function that tends to create confusion for beginners. The goal is to generate lists of file names from patterns, usually to generate the input of a "gathering" / "summarizing" rule. If you want a single file name, I would recommend to use standard Python (>=3.6) f-strings. $\endgroup$
    – bli
    Dec 26, 2021 at 9:41
  • $\begingroup$ However, in your case, I think that the use of a function as input will likely be part of the solution: snakemake.readthedocs.io/en/stable/snakefiles/… Maybe you can have a look at this answer for some ideas of how to proceed: stackoverflow.com/a/64660447/1878788 $\endgroup$
    – bli
    Dec 26, 2021 at 9:45
  • $\begingroup$ Thank you, @bli - I'll look into this - functions as input might make the most sense for my case. $\endgroup$
    – Ram RS
    Dec 27, 2021 at 15:38

1 Answer 1

1
$\begingroup$

Using pointers provided by @bli, I made a few changes and the Snakefile started working as expected.

Instead of using expand for the rule's input, I switched to a function that accepts the wildcard sample as input and returns a dict with the actual inputs.

 # Could be a one-liner but I wrote it in a verbose
 # manner so it's easier to read
 def get_input_files(in_sample_id):
     fq1 = samples.loc[str(in_sample_id), 'FQ1_id']
     fq2 = samples.loc[str(in_sample_id), 'FQ2_id']
     logdir = "{in_sample}.logs".format(in_sample=str(in_sample_id))
     return_dict = { 'fq1' : fq1, 'fq2' : fq2, 'logdir' : logdir }
     return(return_dict)

And the rule my_job is now changed (just the input part):

rule my_job:
    input:
        unpack(get_input_files)
    output:
        "{sample}_1.fastq.gz",
        "{sample}_2.fastq.gz",
...
...

Hope this helps someone looking for a similar solution in the future!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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