snakemake: How to use glob_wildcards for newly created files?

user_10299385 0
Viewed 3258

The problem:

I have a large workflow which creates at some point an arbitrary number of files per {sample}, named e.g. test1.txt, test2.txt, etc.

I then need to use these files for further processing. The input files for the next rule are then {sample}/test1.txt, {sample}/test2.txt, etc. Thus test1, test2, etc become wildcards.

The data structure is:

---data
 ---sample1
   ---test1.txt
   ---test2.txt
   ---test3.txt
 ---sample2
   ---test1.txt
   ---test2.txt
Snakefile

I am struggling how snakemake can be used for such problems. I have looked into the function glob_wildcards, but couldn't figure out how to use it.

Intuitively, I would have done something like this:

samples = ['sample1', 'sample2']

rule append_hello:
  input:
    glob_wildcards('data/{sample}/{id}.txt')
  output:
    'data/{sample}/{id}_2.txt'
  shell:
    " echo {input} 'hello' >> {output} "

I have two questions:

  1. How can this problem be handled in Snakemkae?
  2. How would you construct a rule all in order to run this.

Any inputs or any hints into further reading would be appreciated.

Edit

I think it has to do with wildcard constraints. When I run:

assemblies = []
for filename in glob_wildcards(os.path.join("data/{sample}", "{i}.txt")):
    assemblies.append(filename)
print(assemblies)

I get two lists where the corresponding index matches:

[['sample1', 'sample1', 'sample1', 'sample2', 'sample2'], ['test1', 'test2', 'test3', 'test5', 'test4']]

Now I basically only need to tell snakemake to use corresponding wildcards values.

snakemake python
1 Answers

You problem is that glob_wildcards is evaluated only once before any rule is executed, so the system doesn't know how many files would that rule generate.

What you need is a checkpoint. This feature allows you to stop Snakemake on some point and reevaluate the DAG.

samples = ["sample1", "sample2"]

rule all:
    input:
        expand("data/{sample}/processed.txt", sample=samples)


checkpoint generate_arbitrary:
    output:
        directory("data/{sample}/arbitrary")
    run:
        if wildcards.sample == "sample1":
            n = 3
        else:
            n = 2

        shell("mkdir {output}")
        for id in range(1, n + 1):
            shell(f"echo '{{id}}' > data/{wildcards.sample}/arbitrary/{id}.txt")



def aggregate_input(wildcards):
    checkpoints.generate_arbitrary.get(sample=wildcards.sample)
    ids = glob_wildcards(f"data/{wildcards.sample}/arbitrary/{{id}}.txt").id
    return expand(f"data/{wildcards.sample}/arbitrary/{{id}}.txt", id=ids)


rule append_hello:
    input:
        aggregate_input
    output:
        "data/{sample}/processed.txt"
    shell:
        "echo {input} 'hello' > {output}"
user_112186870
Related