2
$\begingroup$

All, I have a workflow that in some rules requires the sample's read-lengths. For this I always ran a script that would parse all fastq files and would add the read lengths to the sample table. Rules could then take the read lengths for each sample from there during execution.

However, I recently made a workflow/pipeline that downloads files from SRA, so I need to compute the read lengths at a later stage during workflow execution (after prefetch/fasterq-dump).

I added a step to update the sample table, however, I think the sample table is only read at the very beginning of the workflow.

What would be the nicest way, the Snakemake preferred way, to re-read the samples.tsv file before a rule is executed?

In general I like storing meta data in the sample table so, I'll probably use this mechanism more often in the future, if I find a nice methods to re-parse the samples table.

Right now I do this in rules:

rule deletion_profile:
    input:
        rseqc_input_bam,
        rules.add_max_read_length.output.done_file
    output:
        os.path.join(rseqc_dir, '{sample}.deletion_profile.txt')
    conda:
        "../envs/rseqc.yaml"
    params:
        read_length = lambda wildcards: samples.loc[wildcards.sample, 'read_length']
    shell:
        '''
        deletion_profile.py -i {input} -l {params.read_length} -o {rseqc_dir}/{wildcards.sample}

Thanx.

$\endgroup$

1 Answer 1

3
$\begingroup$

You're right, it doesn't re-run things again for each job, so your original definition of samples is still the one it's working with.

It seems like it goes against the grain of Snakemake to overwrite a file in-place that acts as input to a rule, so I think the most Snakemake-y way you could do it would be to treat your samples.tsv as part of the input/output graph like everything else. You could have a rule with samples.tsv as input and samples.somethingelse.tsv as output where that latter one has read lengths and whatever else gets automatically added. (Alternatively, you could just have your updated samples.tsv read within the rule itself rather than relying on a shared copy loaded earlier, which is kind of halfway to the above idea.)

Either way, maybe something like:

rule deletion_profile:
    input:
        bam=rseqc_input_bam,
        samples="samples.extended.tsv"
    output:
        os.path.join(rseqc_dir, '{sample}.deletion_profile.txt')
    conda:
        "../envs/rseqc.yaml"
    run:
        samples = whatever_loads_samples(input.samples)
        read_length = samples.loc[wildcards.sample, 'read_length']
        shell("deletion_profile.py -i {input.bam} -l {read_length} -o {output}")

(On a side note, you can use {output} rather than repeating yourself for that part.)

$\endgroup$

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.