2
$\begingroup$

In Snakemake I want to use two input files for comparison. The one file has a sample name (e.g. 526846), with an experimental id (e.g. EF_S13' and 'MF_S5'). The file I want to compare it to has the same sample name only. I want to compare each sample with experiment ID (first files) to the file with corresponding sample name (second file). These files are in different directories. I thought I would be able to use two expand functions with two lists, as shown in the code below. Note: these are example lists and the lists in reality is muchlarger.

Relevant code:


SAMPLE_SET = ['12345', '56789', '54321']
ID_SET = ['EF_S13', 'MF_S5', 'EF_S11', 'MF_S3', 'EF_S10', 'MF_S2']

rule all:
    input:
        array = expand("results/{sample}_{samp_id}_array/{sample}_{samp_id}.summary.csv", zip, sample = SAMPLE_SET, samp_id = ID_SET)


rule happy_rtg_array:
    """ Normalizes and compares vcf files with array samples"""
    input:
        query_vcf = expand("WESP_run/results/variant_calling/{sample}_{samp_id}.vcf", zip, sample = SAMPLE_SET, samp_id = ID_SET),
        truth_vcf = expand(os.path.join(config["truth_dir_vcf_gsa"],"{sample}.vcf.gz"), sample = SAMPLE_SET),
    output: 
        "WESPValidation/results/{sample}_{samp_id}_array/{sample}_{samp_id}.summary.csv"
[more code]

However, when I look at the output log file with --dry-run, I get as input for the job a huge list with all the inputs in one go, rather then per file/comparison, and one output file (see relevant output), and that times the amount of the (right) combinations of the file name and id (counted from the output).

Also the wildcards are not what I thought they should be. I might be confused with the expand function, as I thought how I used them now should "lock in" the wildcards and become placeholders filled in by the lists.

Relevant output:

[Tue Oct  3 14:31:20 2023]
rule happy_rtg_array:
    input: 12345_EF_S13.vcf, 12345_MF_S5.vcf, 56789_EF_S11.vcf, 56789_MF_S3.vcf, 54321_EF_S10.vcf, 54321_MF_S2.vcf, 12345.vcf, 56789.vcf, 54321.vcf
    output: results/12345_EF_S13_array/12345_EF_S13.summary.csv
    log: happy_rtg_array/12345_EF_S13.log
    jobid: 6
    reason: Missing output files: WESPValidation/results/12345_EF_S13_array/12345_EF_S13.summary.csv
    wildcards: sample=12345_EF, samp_id=S13
    threads: 10
    resources: mem_mb=10000, disk_mb=7855, tmpdir=/tmp, threads=1, runtime_min=120

I doubt this is how it should be, and wonder what I am doing wrong. Is there a better/correct way to do this.

Thank you in advance

$\endgroup$

2 Answers 2

1
$\begingroup$

Using expand in your rule's input: causes snakemake to aggregate all vcfs into one list and leads to your undesired behaviour. By simply removing the expand calls from the input: you should get your desired result.

https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#the-expand-function

$\endgroup$
1
  • $\begingroup$ This is only partially the answer, I tried this before indeed, but still get the following wildcards "wildcards: sample=12345_EF, samp_id=S13" and they should be: "wildcards: sample=12345, samp_id=EF_S13", to get the right file name for the input "truth_vcf" $\endgroup$
    – Dandelion
    Oct 5, 2023 at 12:38
0
$\begingroup$

After some trial and error, I fixed it with the expand function only in the rule all and removing the expand from the input in the rule happy_rtg_array (thanks to @Pallie), and by changing the output file paths:

SAMPLE_SET = ['12345', '56789', '54321']
ID_SET = ['EF_S13', 'MF_S5', 'EF_S11', 'MF_S3', 'EF_S10', 'MF_S2']

rule all:
    input:
        array = expand("results/{sample}/{samp_id}/summary.csv", zip, sample = SAMPLE_SET, samp_id = ID_SET)


rule happy_rtg_array:
    """ Normalizes and compares vcf files with array samples"""
    input:
        query_vcf = "_".join(["WESP_run/results/variant_calling/{sample}","{samp_id}.vcf"]),
        truth_vcf = os.path.join(config["truth_dir_vcf_gsa"],"{sample}.vcf.gz"),
    output: 
        "WESPValidation/results/{sample}/{samp_id}/summary.csv"
[more code]

Basically, I replaced all the underscores with slashes. For some reason these underscores are messing up the wildcards when trying to connect them with it. Maybe someone more knowledgeable knows why this is.

[Tue Oct  5 12:24:10 2023]
rule happy_rtg_array:
    input: 12345_EF_S13.vcf
    output: results/12345/EF_S13/summary.csv
    log: happy_rtg_array/12345_EF_S13.log
    jobid: 6
    reason: Missing output files: WESPValidation/results/12345/EF_S13/summary.csv
    wildcards: sample=12345, samp_id=EF_S13
    threads: 10
    resources: mem_mb=10000, disk_mb=7855, tmpdir=/tmp, threads=1, runtime_min=120
$\endgroup$
2
  • 1
    $\begingroup$ For the reason why, check the part about ambiguous wildcards: snakemake.readthedocs.io/en/stable/snakefiles/…. I missed the amibguity of your wildcards when I wrote my answer :) $\endgroup$
    – Pallie
    Oct 5, 2023 at 13:19
  • $\begingroup$ Dank je @Pallie :-) $\endgroup$
    – Dandelion
    Oct 5, 2023 at 13:30

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.