<p>Ok, so I feel like a right PawPaw. This is <strong>NOT a snakemake problem</strong>. The issue was multi-part but it was me, python and good ol' PBS/Torque Cluster environments.</p>
<h3>Empty input files:</h3>
<p>I was using a wrapper python script by the name of <a href="https://github.com/sritchie73/liftOverPlink" rel="nofollow noreferrer">liftOverPlink.py</a> to convert my files from PLINK-1.9's binary file format into the format needed by the software liftOver by UCSC. In my infinite wisdom, I did make some custom alterations. Turns out in python, there is a big difference between <code>&quot;w+&quot;</code> and <code>&quot;a+&quot;</code> for the <code>open()</code> function <em>(<code>&quot;w+&quot;</code> = write mode and create if needed, <code>&quot;a+&quot;</code> = append mode and create if needed)</em>. <code>&quot;w+&quot;</code> opens or creates a file but places the cursor at the beginning of the file, effectively wiping the file if it already had content inside, wheras <code>&quot;a+&quot;</code> puts the cursor at the end of the files content (which is beginning of file if it had to create file). I used <code>&quot;w+&quot;</code> which was effectively erasing my input files content...</p>
<h3>Why alter the script in the first place?</h3>
<p>As it turns out, our PBS/Torque cluster was not properly serving up the software at runtime (software in question being liftOver). From what I have been able to gather, normally liftOverPlink would create intermediate files, liftOver would run and create its files where necessary, and liftOverPlink would then step in and perform format conversion back into PLINK-1.9 format. Since liftOver was not firing, the files did not exist and liftOverPlink kicked up errors since it was never designed to handle that (functionality to detect missing files would have been redundant). I attempted to ugly-fix it by using</p>
<pre class="lang-py prettyprint-override"><code>filename = open(filename, &quot;w+&quot;)
filename.close()
</code></pre>
<p>to open and close the file, effectively forcefully creating it if needed. Unfortunately that let to the above problem. Providing an absolute path to the liftOver executable fixed runtime issues and my custom code was no-longer necessary.</p>

CC BY-SA 4.0
2020-04-08T12:18:11.087Z
2020-04-08T12:18:11.087Z
2020-06-18T08:30:41.743Z
Graeme Ford
<h1>Problem:</h1>
<p>When I run my Snakefile on our cluster, all the processes complete successfully. I can see from the logs that my shell scripts are being run correctly and are claiming to complete successfully. About half-way through when analysis steps begins, the workflow crashes due to empty input files (which is normal for <a href="https://www.cog-genomics.org/plink/1.9/" rel="nofollow noreferrer">PLINK-1.9</a>, what's not normal is that the input files are empty). I went back and checked, most of the data-containing files (<code>*.map</code> and <code>*.ped</code>, in some cases/processes Plinks binary file set of <code>*.bed</code>, <code>*.bim</code> and <code>*.fam</code>) are completely empty since the beginning of the workflow? (Logs still claim processes all completed and created output files successfully?)</p>
<h1>Nuances:</h1>
<ul>
<li>When I was playing around, I accidentally messed up the dynamic naming variables (escaped the variable interpret syntax, causing literal usage) and caused creation of single file with actual variable names, not dataset names. When I did this however, the file was not empty?</li>
<li>Some 'side' files produced by the software (E.g. PLINK-1.9 produces a <code>*.nosex</code> file indicating samples with missing sex information) are not empty but the actuall data-containing output is empty?</li>
<li>Had a problem earlier using <code>cut -f 4 {params.prefix}.unlifted | sed '/^#/d' &gt;&gt; {params.prefix3}</code> <em>(Intention was to collate list of entries ID's to exclude in next workflow step)</em>, where the use of &gt;&gt; with Snakemake on cluster environment simply would not be allowed (Permission denied error despite it being created line before by the same script?). Had to fix by adding in <code>chmod a+rwx {params.prefix3}</code> just before. Might this be a symptom of the same problem?</li>
<li>This is being executed on a PBS/Torque cluster environment where the scripts are run as user who submitted them so permissions are not <em>supposed</em> to be an issue. Cluster profile was built for Snakemake using <a href="https://github.com/snakemake-profiles/pbs-torque" rel="nofollow noreferrer">this PBS/Torque Profile</a>. <em>(I just stripped out dependency flags because they were breaking things)</em></li>
</ul>
<h1>Code:</h1>
<p>Here is the first three processes. So far, it crashes at the third <code>CLEAN</code> rule as the software PLINK-1.9 does not play nice with empty inputs.</p>
<pre><code>rule TRIM:
    input:
        expand(&quot;rawData/{{sample}}.{extension}&quot;, extension=['bed', 'bim', 'fam'])

    output:
        expand(&quot;Results/{{sample}}_{{location}}_TRIMMED.{fileExtensions}&quot;, fileExtensions=[&quot;map&quot;, &quot;ped&quot;])

    params:
        fromBP = lambda wildcards: config[&quot;locations&quot;][wildcards.location][&quot;GRCh37&quot;][&quot;from&quot;],
        toBP = lambda wildcards: config[&quot;locations&quot;][wildcards.location][&quot;GRCh37&quot;][&quot;to&quot;],
        chr = lambda wildcards: config[&quot;locations&quot;][wildcards.location][&quot;GRCh37&quot;][&quot;chromosome&quot;]

    shell:
        &quot;&quot;&quot;
        module load plink-1.9
        plink --bfile rawData/{wildcards.sample} --chr {params.chr} --set-missing-var-ids @_# --recode --keep-allele-order --from-bp {params.fromBP} --to-bp {params.toBP} --out Results/{wildcards.sample}_{wildcards.location}_TRIMMED
        &quot;&quot;&quot;


#/********* LiftOverPlink (LIFTOVER) *********/

rule LIFTOVER:
    input:
        expand(&quot;Results/{{sample}}_{{location}}_TRIMMED.{fileExtensions}&quot;, fileExtensions=[&quot;map&quot;, &quot;ped&quot;])

    output:
        outMap=&quot;Results/{sample}_{location}_LIFTED.map&quot;,
        outPed=&quot;Results/{sample}_{location}_LIFTED.ped&quot;,
        exclusion=&quot;Results/{sample}_{location}_EXCLUDE.dat&quot;

    params:
        prefix=lambda wildcards: &quot;Results/{}_{}_PRE_FILTER&quot;.format(wildcards.sample, wildcards.location),
        prefix2=lambda wildcards: &quot;Results/{}_{}_POST_FILTER&quot;.format(wildcards.sample, wildcards.location),
        prefix3=lambda wildcards: &quot;Results/{}_{}_EXCLUDE.dat&quot;.format(wildcards.sample, wildcards.location),
        chainFile=&quot;Binaries/hg19ToHg38.over.chain&quot;,
        LiftOver=&quot;Binaries/liftOverPlink.py&quot;,
        rmBadLifts=&quot;Binaries/rmBadLifts.py&quot;,
        sexes=&quot;rawData/1000g.sexes&quot;

    run:
        print(&quot;Determining Liftover requirements now...&quot;)
        if config['samples'][wildcards.sample]['refGenome'] != &quot;GRCh38&quot;:
            print(&quot;Liftover required. Dataset {} is mapped to {}&quot;.format(wildcards.sample, config['samples'][wildcards.sample]['refGenome'])),
            shell(&quot;module load liftover&quot;),
            if config['samples'][wildcards.sample]['refGenome'] == &quot;GRCh37&quot; or config['samples'][wildcards.sample]['refGenome'] == &quot;Hg19&quot;:
                print(&quot;Lifting from GRCh37 to GRCh38.&quot;),
                shell(&quot;python {params.LiftOver} -m Results/{wildcards.sample}_{wildcards.location}_TRIMMED.map -p Results/{wildcards.sample}_{wildcards.location}_TRIMMED.ped -o {params.prefix} -c {params.chainFile}&quot;),
                shell(&quot;ls Results&quot;),
                shell(&quot;python {params.rmBadLifts} --map {params.prefix}.map --out {output.outMap} --log {params.prefix2}.log&quot;),
                shell(&quot;cut -f 2 {params.prefix2}.log &gt; {params.prefix3}&quot;),
                shell(&quot;chmod a+rwx {params.prefix3}&quot;),
                shell(&quot;cut -f 4 {params.prefix}.unlifted | sed '/^#/d' &gt;&gt; {params.prefix3}&quot;),
                shell(&quot;less {params.prefix3} | {output.exclusion}&quot;),
                shell(&quot;awk 'NR==FNR{{a[<span class="math-container">$1]=$</span>4;next}}{{ if (a[<span class="math-container">$1]!= 'NULL') {{ if (a[$1] == 'female') {{ $5=2; print }} else if (a[$1] == 'male') {{ $5=1; print }} else {{$</span>5=0; print}}}}}}' {params.sexes} {params.prefix}.ped &gt; Results/{wildcards.sample}_{wildcards.location}_LIFTED.ped&quot;)
                # ToDo: Add conditionals for other human reference genome builds
        else:
            print(&quot;No liftover required. Dataset {} is already mapped to GRCh38.&quot;.format(wildcards.sample)),
            shell(&quot;cp Results/{wildcards.sample}_{wildcards.location}_TRIMMED.map Results/{wildcards.sample}_{wildcards.location}_LIFTED.map&quot;),
            shell(&quot;cp Results/{wildcards.sample}_{wildcards.location}_TRIMMED.ped Results/{wildcards.sample}_{wildcards.location}_LIFTED.ped&quot;),
            shell(&quot;touch Resuts/{wildcards.sample}_{wildcards.location}_EXCLUDE.dat&quot;)


rule CLEAN:
    input:
        mapFile=&quot;Results/{sample}_{location}_LIFTED.map&quot;,
        pedFile=&quot;Results/{sample}_{location}_LIFTED.ped&quot;,
        exclude=&quot;Results/{sample}_{location}_EXCLUDE.dat&quot;

    output:
        expand(&quot;Results/{{sample}}_{{location}}_CLEANED.{extension}&quot;, extension=['bed', 'bim', 'fam'])

    shell:
        &quot;&quot;&quot;
        module load plink-1.9
        plink --map {input.mapFile} --ped {input.pedFile} --out Results/{wildcards.sample}_{wildcards.location}_CLEANED --make-bed --keep-allele-order --exclude {input.exclude}
        &quot;&quot;&quot;
<span class="math-container">```</span>
</code></pre>

CC BY-SA 4.0
2020-04-07T15:01:43.663Z
#Update:
I just re-enabled the dependency flags and now all outputs are looking normal again? (When I don't remove it, my cluster management software (PBS/Torque) kicks up an error, saying `Invalid Job Dependancy`) Why would dependancy flags on `qsub` matter or impact output file write operations?
Graeme Ford
CC BY-SA 4.0
2020-04-07T15:53:39.967Z
https://stackoverflow.com/questions/54964785/the-touch-of-death
This is extremely enlightening. Turns out it may be a Snakemake problem?
Graeme Ford
CC BY-SA 4.0
2020-04-07T16:15:52.640Z
Ok, so I have done some more debugging. It seems if I interrupt my Workflow between the rule `TRIM` and `LIFTOVER`, I can see that it is in fact creating valid outputs. As soon as I re-run the pipeline and it carries on, it wipes the input files and then as a result, generates empty output files? Why does Snakemake do this?
Graeme Ford
CC BY-SA 4.0
2020-04-07T16:57:19.260Z
What command are you using to submit your jobs to cluster?
Phoenix Mu
CC BY-SA 4.0
2020-04-07T17:01:14.210Z
It’s a PBS/Torque cluster so uses the qsub command
Graeme Ford
CC BY-SA 4.0
2020-04-07T17:26:41.973Z
snakeMake isn't generating empty output files unless you have it `touch()` them. It will delete them if it thinks an input file may have changes (e.g., due to being `temp()` and therefore having been previously removed).
Devon Ryan
CC BY-SA 4.0
2020-04-07T17:31:14.863Z
Ok but my problem is the inputs. It’s wiping the inputs? Also I didn’t use touch or temp (mostly bc I didn’t know touch existed and bc I want the intermediate files as well XD)
Graeme Ford
<snakemake><workflow-management>
Snakemake on cluster produces empty output files
# Snakemake on cluster produces empty output files

(Asked by: username_0 on 2020-04-07)

Problem:
When I run my Snakefile on our cluster, all the processes complete successfully. I can see from the logs that my shell scripts are being run correctly and are claiming to complete successfully. About half-way through when analysis steps begins, the workflow crashes due to empty input files (which is normal for PLINK-1.9, what's not normal is that the input files are empty). I went back and checked, most of the data-containing files (*.map and *.ped, in some cases/processes Plinks binary file set of *.bed, *.bim and *.fam) are completely empty since the beginning of the workflow? (Logs still claim processes all completed and created output files successfully?)
Nuances:

When I was playing around, I accidentally messed up the dynamic naming variables (escaped the variable interpret syntax, causing literal usage) and caused creation of single file with actual variable names, not dataset names. When I did this however, the file was not empty?
Some 'side' files produced by the software (E.g. PLINK-1.9 produces a *.nosex file indicating samples with missing sex information) are not empty but the actuall data-containing output is empty?
Had a problem earlier using cut -f 4 {params.prefix}.unlifted | sed '/^#/d' >> {params.prefix3} (Intention was to collate list of entries ID's to exclude in next workflow step), where the use of >> with Snakemake on cluster environment simply would not be allowed (Permission denied error despite it being created line before by the same script?). Had to fix by adding in chmod a+rwx {params.prefix3} just before. Might this be a symptom of the same problem?
This is being executed on a PBS/Torque cluster environment where the scripts are run as user who submitted them so permissions are not supposed to be an issue. Cluster profile was built for Snakemake using this PBS/Torque Profile. (I just stripped out dependency flags because they were breaking things)

Code:
Here is the first three processes. So far, it crashes at the third CLEAN rule as the software PLINK-1.9 does not play nice with empty inputs.
rule TRIM:
    input:
        expand("rawData/{{sample}}.{extension}", extension=['bed', 'bim', 'fam'])

    output:
        expand("Results/{{sample}}_{{location}}_TRIMMED.{fileExtensions}", fileExtensions=["map", "ped"])

    params:
        fromBP = lambda wildcards: config["locations"][wildcards.location]["GRCh37"]["from"],
        toBP = lambda wildcards: config["locations"][wildcards.location]["GRCh37"]["to"],
        chr = lambda wildcards: config["locations"][wildcards.location]["GRCh37"]["chromosome"]

    shell:
        """
        module load plink-1.9
        plink --bfile rawData/{wildcards.sample} --chr {params.chr} --set-missing-var-ids @_# --recode --keep-allele-order --from-bp {params.fromBP} --to-bp {params.toBP} --out Results/{wildcards.sample}_{wildcards.location}_TRIMMED
        """


#/********* LiftOverPlink (LIFTOVER) *********/

rule LIFTOVER:
    input:
        expand("Results/{{sample}}_{{location}}_TRIMMED.{fileExtensions}", fileExtensions=["map", "ped"])

    output:
        outMap="Results/{sample}_{location}_LIFTED.map",
        outPed="Results/{sample}_{location}_LIFTED.ped",
        exclusion="Results/{sample}_{location}_EXCLUDE.dat"

    params:
        prefix=lambda wildcards: "Results/{}_{}_PRE_FILTER".format(wildcards.sample, wildcards.location),
        prefix2=lambda wildcards: "Results/{}_{}_POST_FILTER".format(wildcards.sample, wildcards.location),
        prefix3=lambda wildcards: "Results/{}_{}_EXCLUDE.dat".format(wildcards.sample, wildcards.location),
        chainFile="Binaries/hg19ToHg38.over.chain",
        LiftOver="Binaries/liftOverPlink.py",
        rmBadLifts="Binaries/rmBadLifts.py",
        sexes="rawData/1000g.sexes"

    run:
        print("Determining Liftover requirements now...")
        if config['samples'][wildcards.sample]['refGenome'] != "GRCh38":
            print("Liftover required. Dataset {} is mapped to {}".format(wildcards.sample, config['samples'][wildcards.sample]['refGenome'])),
            shell("module load liftover"),
            if config['samples'][wildcards.sample]['refGenome'] == "GRCh37" or config['samples'][wildcards.sample]['refGenome'] == "Hg19":
                print("Lifting from GRCh37 to GRCh38."),
                shell("python {params.LiftOver} -m Results/{wildcards.sample}_{wildcards.location}_TRIMMED.map -p Results/{wildcards.sample}_{wildcards.location}_TRIMMED.ped -o {params.prefix} -c {params.chainFile}"),
                shell("ls Results"),
                shell("python {params.rmBadLifts} --map {params.prefix}.map --out {output.outMap} --log {params.prefix2}.log"),
                shell("cut -f 2 {params.prefix2}.log > {params.prefix3}"),
                shell("chmod a+rwx {params.prefix3}"),
                shell("cut -f 4 {params.prefix}.unlifted | sed '/^#/d' >> {params.prefix3}"),
                shell("less {params.prefix3} | {output.exclusion}"),
                shell("awk 'NR==FNR{{a[$1]=$4;next}}{{ if (a[$1]!= 'NULL') {{ if (a[$1] == 'female') {{ $5=2; print }} else if (a[$1] == 'male') {{ $5=1; print }} else {{$5=0; print}}}}}}' {params.sexes} {params.prefix}.ped > Results/{wildcards.sample}_{wildcards.location}_LIFTED.ped")
                # ToDo: Add conditionals for other human reference genome builds
        else:
            print("No liftover required. Dataset {} is already mapped to GRCh38.".format(wildcards.sample)),
            shell("cp Results/{wildcards.sample}_{wildcards.location}_TRIMMED.map Results/{wildcards.sample}_{wildcards.location}_LIFTED.map"),
            shell("cp Results/{wildcards.sample}_{wildcards.location}_TRIMMED.ped Results/{wildcards.sample}_{wildcards.location}_LIFTED.ped"),
            shell("touch Resuts/{wildcards.sample}_{wildcards.location}_EXCLUDE.dat")


rule CLEAN:
    input:
        mapFile="Results/{sample}_{location}_LIFTED.map",
        pedFile="Results/{sample}_{location}_LIFTED.ped",
        exclude="Results/{sample}_{location}_EXCLUDE.dat"

    output:
        expand("Results/{{sample}}_{{location}}_CLEANED.{extension}", extension=['bed', 'bim', 'fam'])

    shell:
        """
        module load plink-1.9
        plink --map {input.mapFile} --ped {input.pedFile} --out Results/{wildcards.sample}_{wildcards.location}_CLEANED --make-bed --keep-allele-order --exclude {input.exclude}
        """
```

**Tags:** snakemake, workflow-management
> **Comments:**
> **username_0**: #Update: I just re-enabled the dependency flags and now all outputs are looking normal again? (When I don't remove it, my cluster management software (PBS/Torque) kicks up an error, saying `Invalid Job Dependancy`) Why would dependancy flags on `qsub` matter or impact output file write operations?  
> **username_0**: https://stackoverflow.com/questions/54964785/the-touch-of-death This is extremely enlightening. Turns out it may be a Snakemake problem?  
> **username_0**: Ok, so I have done some more debugging. It seems if I interrupt my Workflow between the rule `TRIM` and `LIFTOVER`, I can see that it is in fact creating valid outputs. As soon as I re-run the pipeline and it carries on, it wipes the input files and then as a result, generates empty output files? Why does Snakemake do this?  
> **username_1**: What command are you using to submit your jobs to cluster?  
> **username_0**: It’s a PBS/Torque cluster so uses the qsub command  
> **username_2**: snakeMake isn't generating empty output files unless you have it `touch()` them. It will delete them if it thinks an input file may have changes (e.g., due to being `temp()` and therefore having been previously removed).  
> **username_0**: Ok but my problem is the inputs. It’s wiping the inputs? Also I didn’t use touch or temp (mostly bc I didn’t know touch existed and bc I want the intermediate files as well XD)  

---

## 1 answers

### Answer by username_0 on 2020-04-08. Score: 1

Ok, so I feel like a right PawPaw. This is NOT a snakemake problem. The issue was multi-part but it was me, python and good ol' PBS/Torque Cluster environments.
Empty input files:
I was using a wrapper python script by the name of liftOverPlink.py to convert my files from PLINK-1.9's binary file format into the format needed by the software liftOver by UCSC. In my infinite wisdom, I did make some custom alterations. Turns out in python, there is a big difference between "w+" and "a+" for the open() function ("w+" = write mode and create if needed, "a+" = append mode and create if needed). "w+" opens or creates a file but places the cursor at the beginning of the file, effectively wiping the file if it already had content inside, wheras "a+" puts the cursor at the end of the files content (which is beginning of file if it had to create file). I used "w+" which was effectively erasing my input files content...
Why alter the script in the first place?
As it turns out, our PBS/Torque cluster was not properly serving up the software at runtime (software in question being liftOver). From what I have been able to gather, normally liftOverPlink would create intermediate files, liftOver would run and create its files where necessary, and liftOverPlink would then step in and perform format conversion back into PLINK-1.9 format. Since liftOver was not firing, the files did not exist and liftOverPlink kicked up errors since it was never designed to handle that (functionality to detect missing files would have been redundant). I attempted to ugly-fix it by using
filename = open(filename, "w+")
filename.close()

to open and close the file, effectively forcefully creating it if needed. Unfortunately that let to the above problem. Providing an absolute path to the liftOver executable fixed runtime issues and my custom code was no-longer necessary.


---

