Snakemake on cluster produces empty output files

Question

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}
        """
```

Answer 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.