I'm stuck on how to do this with Snakemake.
First, say my "all" rule is:
rule all:
input: "SA.txt", "SA_T1.txt", "SA_T2.txt",
"SB.txt", "SB_T1.txt", "SB_T2.txt", "SB_T3.txt"
Notice that SA
has two _T#
files while SB
has three such files, a crucial element of this.
Now I want to write a rule like this to generate these files:
rule X:
output: N="S{X}.txt", T="S{X}_T{Y}.txt"
(etc.)
But SnakeMake requires that both output templates have the same wildcards, which these don't. Further, even if SnakeMake could handle the multiple wildcards, it would presumably want to find a single filename match for the S{X}_T{Y}.txt
template, but I want that to match ALL files where {X}
matches the first template's {X}
, i.e. I want output.T
to be a list, not a single file. So it would seem the way to do it is:
def myExpand(T, wildcards):
T = T.replace("{X}", wildcards.X)
T = [T.replace("{Y}", S) for S in theXYs[wildcards.X]]
return T
rule X:
output: N="S{X}.txt", T=lambda wildcards: myExpand("S{X}_T{Y}.txt", wildcards)
(etc.)
But I can't do this, because a lambda function cannot be used in an output section.
How do I do it?
It seems to me this argues for supporting lambda functions on output statements, providing a wildcards dictionary filled with values from already-parsed sections of the output statement.
Additions responding to comments:
The value of wildcard Y is needed because other rules have inputs for those files that have the wildcard Y.
My rule knows the different values for Y (and X) that it needs to work with from data read from a database into python dictionaries.
There are many values for X, and from 2 to 6 values of Y for each value of X. I don't think it makes sense to use separate rules for each value of X. However, I might be wrong as I recently learned that one can put a rule inside a loop, and create multiple rules.
More info about the workflow: I am combining somatic variant VCF files for several tumor samples from one person together into a single VCF file, and doing it such that for each called variant in any one tumor, all tumors not calling that variant are analyzed to determine read depth at the variant, which is included in the merged VCF file.
The full process involves about 14 steps, which could perhaps be as many as 14 rules. I actually didn't want to use 14 rules, but preferred to just do it all in one rule.
However, I now think the solution is to indeed use lots of separate rules. I was avoiding this partly because of the large number of unnecessary intermediate files, but actually, these exist anyway, temporarily, within a single large rule. With multiple rules I can mark them temp() so Snakemake will delete them at the end.
For the sake of fleshing out this discussion, which I believe is a legitimate one, let's assume a simple situation that might arise. Say that for each of a number of persons, you have N (>=2) tumor VCF files, as I do, and that you want to write a rule that will produce N+1 output files, one output file per tumor plus one more file that is associated with the person. Use wildcard X for person ID and wildcard Y for tumor ID within person X. Say that the operation is to put all variants present in ALL tumor VCF files into the person output VCF file, and all OTHER variants into the corresponding tumor output files for the tumors in which they appear. Say a single program generates all N+1 files from the N input files. How do you write the rule?
You want this:
rule:
output: COMMON="{X}.common.vcf", INDIV="{X}.{Y}.indiv.vcf"
input: "{X}.{Y}.vcf"
shell: """
getCommonAndIndividualVariants --inputs {input} \
--common {output.COMMON} --indiv {output.INDIV}
"""
But that violates the rules for output wildcards.
The way I did it, which is less than satisfactory, but works, is to use two rules, the first of which has the output template with more wildcards, and the second the template with fewer wildcards, and having the second rule create temporary output files which are renamed to the final name by the first rule:
rule A:
output: "{X}.{Y}.indiv.vcf"
input: "{X}.common.vcf"
run: "for infile in {input}: os.system('mv '+infile+'.tmp'+' '+infile)"
rule B:
output: "{X}.common.vcf"
input: lambda wildcards: \
expand("{X}.{Y}.vcf", **wildcards, Y=getYfromDB(wildcards["X"]))
params: OUT=lambda wildcards: \
expand("{X}.{Y}.indiv.vcf.tmp", Y=getYfromDB(wildcards["X"]))
shell: """
getCommonAndIndividualVariants --inputs {input} \
--common {output} --indiv {params.OUT}
"""