diff --git a/MONSDA/Versions.py b/MONSDA/Versions.py index 47d5b3dc..599d5ddd 100644 --- a/MONSDA/Versions.py +++ b/MONSDA/Versions.py @@ -359,6 +359,25 @@ def write_versions(config, steps, envdir, monsda_version, configfile=None, path= rows = collect_versions(config, steps, envdir) + body = "" + if configfile: + body += "# config: {c}\n".format(c=str(configfile)) + body += ( + "# versions as pinned in the conda env yamls of the configured tools, not queried from the running binaries\n" + ) + body += "#tool\tversion\tmultiqc_module\tsteps\n" + for row in rows: + body += "\t".join(row) + "\n" + + if os.path.isfile(path): + with open(path, "r") as existing: + old_body = "".join( + line for line in existing if not line.startswith("# MONSDA ") + ) + if old_body == body: + log.info(logid + "Tool versions unchanged, keeping " + str(path)) + return path + with open(path, "w") as out: out.write( "# MONSDA {v}\t{ts}\n".format( @@ -366,14 +385,7 @@ def write_versions(config, steps, envdir, monsda_version, configfile=None, path= ts=datetime.datetime.now().strftime("%Y%m%d_%H_%M_%S"), ) ) - if configfile: - out.write("# config: {c}\n".format(c=str(configfile))) - out.write( - "# versions as pinned in the conda env yamls of the configured tools, not queried from the running binaries\n" - ) - out.write("#tool\tversion\tmultiqc_module\tsteps\n") - for row in rows: - out.write("\t".join(row) + "\n") + out.write(body) log.info(logid + "Tool versions written to " + str(path)) return path diff --git a/MONSDA/Workflows.py b/MONSDA/Workflows.py index 7d686347..1eef410b 100755 --- a/MONSDA/Workflows.py +++ b/MONSDA/Workflows.py @@ -2356,6 +2356,13 @@ def nf_fetch_params( retconf["REFDIR"] = REFDIR retconf["IP"] = IP + # Resolve all annotation paths to absolute paths once, so every nf + # process can read them reliably regardless of its own work directory + # (mirrors the equivalent hardening in workflows/header.smk). + for key in list(retconf.keys()): + if key.endswith("ANNO") and isinstance(retconf[key], str) and retconf[key]: + retconf[key] = os.path.abspath(retconf[key]) + return retconf diff --git a/tests/golden/config_DTU/nf/allconditions_fastqc-trimgalore-picard-star_DTU_dexseq_subflow.nf b/tests/golden/config_DTU/nf/allconditions_fastqc-trimgalore-picard-star_DTU_dexseq_subflow.nf index 240b5324..74b21753 100644 --- a/tests/golden/config_DTU/nf/allconditions_fastqc-trimgalore-picard-star_DTU_dexseq_subflow.nf +++ b/tests/golden/config_DTU/nf/allconditions_fastqc-trimgalore-picard-star_DTU_dexseq_subflow.nf @@ -152,6 +152,25 @@ process salmon_quant{ """ } } +process decompress_annotation_dtu{ + conda "/envs/${params.gCOUNTENV}"+".yaml" + container "oras://ghcr.io/jfallmann/monsda:"+"${params.gCOUNTENV}"+"-VERSION" + cpus 1 + cache 'lenient' + publishDir "${workflow.workDir}/../" , mode: 'link', + saveAs: {filename -> + if (filename.indexOf("log") > 0) "LOGS/${params.gSCOMBO}/DTU/"+"${params.gCOUNTENV}"+"/decompress_annotation.log" + } + input: + path anno + output: + path "annotation.gtf", emit: anno_plain + path "log", emit: log + script: + """ + zcat -f $anno > annotation.gtf 2> log; if [ ! -s annotation.gtf ]; then echo "ERROR: decompressed annotation from $anno is empty, check source file and filesystem" >> log; exit 1; fi + """ +} process oarfish_dtu_quant{ conda "/envs/${params.gCOUNTENV}"+".yaml" container "oras://ghcr.io/jfallmann/monsda:"+"${params.gCOUNTENV}"+"-VERSION" @@ -163,18 +182,20 @@ process oarfish_dtu_quant{ else "DTU/${params.gSCOMBO}/"+"${params.gCOUNTENV}"+"/${params.gCONDITION}/${file(filename).getName()}" } input: - path bam + path fls output: path "*.gz", includeInputs:false, emit: counts path "*.log", emit: logs script: + anno = fls[0] + bam = fls[1] fn = file(bam).getSimpleName().replaceAll(/\Q_mapped_sorted\E/,"") lf = "${params.gCOUNTENV}_"+fn+".log" prefix = fn+"/"+params.gCOUNTENV oz = fn+"/quant.sf.gz" ol = fn+"_counts.gz" """ - mkdir -p $fn; samtools collate -@ ${task.cpus} -o namecollated.bam $bam &>> $lf; zcat -f ${params.gDTUANNO} > $fn/annotation.gtf 2>> $lf; ${params.gCOUNTBIN} --threads ${task.cpus} --genome-alignments namecollated.bam --annotation $fn/annotation.gtf --genome-fasta ${params.gDTUREF} ${params.gCOUNTPARAMS} --output $prefix &>> $lf && rm -f namecollated.bam $fn/annotation.gtf && python3 ${params.gBINS}/Analysis/oarfish_to_salmon.py -i ${prefix}.quant -o $oz &>> $lf && ln -fs $oz $ol + mkdir -p $fn; : > $lf; if [ ! -s $anno ]; then echo "ERROR: decompressed annotation $anno is empty, check source file and filesystem" >> $lf; exit 1; fi; samtools collate -@ ${task.cpus} -o namecollated.bam $bam &>> $lf; ${params.gCOUNTBIN} --threads ${task.cpus} --genome-alignments namecollated.bam --annotation $anno --genome-fasta ${params.gDTUREF} ${params.gCOUNTPARAMS} --output $prefix &>> $lf && rm -f namecollated.bam && python3 ${params.gBINS}/Analysis/oarfish_to_salmon.py -i ${prefix}.quant -o $oz &>> $lf && ln -fs $oz $ol """ } process create_summary_snippet{ @@ -345,7 +366,9 @@ workflow DTU{ element -> return "${workflow.workDir}/../MAPPED/${params.gCOMBO}/"+element+"_mapped_sorted.bam" } mapsamples_ch = Channel.fromPath(MAPPEDSAMPLES.sort()) - oarfish_dtu_quant(mapsamples_ch.collate(1)) + decompress_annotation_dtu(annofile) + annogtf_dtu = decompress_annotation_dtu.out.anno_plain + oarfish_dtu_quant(annogtf_dtu.combine(mapsamples_ch.collate(1))) quantcounts = oarfish_dtu_quant.out.counts } prepare_dtu_annotation() diff --git a/tests/golden/config_DTU/nf/allconditions_fastqc-trimgalore-picard-star_DTU_drimseq_subflow.nf b/tests/golden/config_DTU/nf/allconditions_fastqc-trimgalore-picard-star_DTU_drimseq_subflow.nf index afbb4b21..129f93e8 100644 --- a/tests/golden/config_DTU/nf/allconditions_fastqc-trimgalore-picard-star_DTU_drimseq_subflow.nf +++ b/tests/golden/config_DTU/nf/allconditions_fastqc-trimgalore-picard-star_DTU_drimseq_subflow.nf @@ -152,6 +152,25 @@ process salmon_quant{ """ } } +process decompress_annotation_dtu{ + conda "/envs/${params.gCOUNTENV}"+".yaml" + container "oras://ghcr.io/jfallmann/monsda:"+"${params.gCOUNTENV}"+"-VERSION" + cpus 1 + cache 'lenient' + publishDir "${workflow.workDir}/../" , mode: 'link', + saveAs: {filename -> + if (filename.indexOf("log") > 0) "LOGS/${params.gSCOMBO}/DTU/"+"${params.gCOUNTENV}"+"/decompress_annotation.log" + } + input: + path anno + output: + path "annotation.gtf", emit: anno_plain + path "log", emit: log + script: + """ + zcat -f $anno > annotation.gtf 2> log; if [ ! -s annotation.gtf ]; then echo "ERROR: decompressed annotation from $anno is empty, check source file and filesystem" >> log; exit 1; fi + """ +} process oarfish_dtu_quant{ conda "/envs/${params.gCOUNTENV}"+".yaml" container "oras://ghcr.io/jfallmann/monsda:"+"${params.gCOUNTENV}"+"-VERSION" @@ -163,18 +182,20 @@ process oarfish_dtu_quant{ else "DTU/${params.gSCOMBO}/"+"${params.gCOUNTENV}"+"/${params.gCONDITION}/${file(filename).getName()}" } input: - path bam + path fls output: path "*.gz", includeInputs:false, emit: counts path "*.log", emit: logs script: + anno = fls[0] + bam = fls[1] fn = file(bam).getSimpleName().replaceAll(/\Q_mapped_sorted\E/,"") lf = "${params.gCOUNTENV}_"+fn+".log" prefix = fn+"/"+params.gCOUNTENV oz = fn+"/quant.sf.gz" ol = fn+"_counts.gz" """ - mkdir -p $fn; samtools collate -@ ${task.cpus} -o namecollated.bam $bam &>> $lf; zcat -f ${params.gDTUANNO} > $fn/annotation.gtf 2>> $lf; ${params.gCOUNTBIN} --threads ${task.cpus} --genome-alignments namecollated.bam --annotation $fn/annotation.gtf --genome-fasta ${params.gDTUREF} ${params.gCOUNTPARAMS} --output $prefix &>> $lf && rm -f namecollated.bam $fn/annotation.gtf && python3 ${params.gBINS}/Analysis/oarfish_to_salmon.py -i ${prefix}.quant -o $oz &>> $lf && ln -fs $oz $ol + mkdir -p $fn; : > $lf; if [ ! -s $anno ]; then echo "ERROR: decompressed annotation $anno is empty, check source file and filesystem" >> $lf; exit 1; fi; samtools collate -@ ${task.cpus} -o namecollated.bam $bam &>> $lf; ${params.gCOUNTBIN} --threads ${task.cpus} --genome-alignments namecollated.bam --annotation $anno --genome-fasta ${params.gDTUREF} ${params.gCOUNTPARAMS} --output $prefix &>> $lf && rm -f namecollated.bam && python3 ${params.gBINS}/Analysis/oarfish_to_salmon.py -i ${prefix}.quant -o $oz &>> $lf && ln -fs $oz $ol """ } process create_summary_snippet{ @@ -346,7 +367,9 @@ workflow DTU{ element -> return "${workflow.workDir}/../MAPPED/${params.gCOMBO}/"+element+"_mapped_sorted.bam" } mapsamples_ch = Channel.fromPath(MAPPEDSAMPLES.sort()) - oarfish_dtu_quant(mapsamples_ch.collate(1)) + decompress_annotation_dtu(annofile) + annogtf_dtu = decompress_annotation_dtu.out.anno_plain + oarfish_dtu_quant(annogtf_dtu.combine(mapsamples_ch.collate(1))) quantcounts = oarfish_dtu_quant.out.counts } prepare_dtu_annotation() diff --git a/tests/golden/config_DTU/nf/allconditions_fastqc-trimgalore-picard-star_DTU_spit_subflow.nf b/tests/golden/config_DTU/nf/allconditions_fastqc-trimgalore-picard-star_DTU_spit_subflow.nf index 2e0ac723..ccf4b98a 100644 --- a/tests/golden/config_DTU/nf/allconditions_fastqc-trimgalore-picard-star_DTU_spit_subflow.nf +++ b/tests/golden/config_DTU/nf/allconditions_fastqc-trimgalore-picard-star_DTU_spit_subflow.nf @@ -152,6 +152,25 @@ process salmon_quant{ """ } } +process decompress_annotation_dtu{ + conda "/envs/${params.gCOUNTENV}"+".yaml" + container "oras://ghcr.io/jfallmann/monsda:"+"${params.gCOUNTENV}"+"-VERSION" + cpus 1 + cache 'lenient' + publishDir "${workflow.workDir}/../" , mode: 'link', + saveAs: {filename -> + if (filename.indexOf("log") > 0) "LOGS/${params.gSCOMBO}/DTU/"+"${params.gCOUNTENV}"+"/decompress_annotation.log" + } + input: + path anno + output: + path "annotation.gtf", emit: anno_plain + path "log", emit: log + script: + """ + zcat -f $anno > annotation.gtf 2> log; if [ ! -s annotation.gtf ]; then echo "ERROR: decompressed annotation from $anno is empty, check source file and filesystem" >> log; exit 1; fi + """ +} process oarfish_dtu_quant{ conda "/envs/${params.gCOUNTENV}"+".yaml" container "oras://ghcr.io/jfallmann/monsda:"+"${params.gCOUNTENV}"+"-VERSION" @@ -163,18 +182,20 @@ process oarfish_dtu_quant{ else "DTU/${params.gSCOMBO}/"+"${params.gCOUNTENV}"+"/${params.gCONDITION}/${file(filename).getName()}" } input: - path bam + path fls output: path "*.gz", includeInputs:false, emit: counts path "*.log", emit: logs script: + anno = fls[0] + bam = fls[1] fn = file(bam).getSimpleName().replaceAll(/\Q_mapped_sorted\E/,"") lf = "${params.gCOUNTENV}_"+fn+".log" prefix = fn+"/"+params.gCOUNTENV oz = fn+"/quant.sf.gz" ol = fn+"_counts.gz" """ - mkdir -p $fn; samtools collate -@ ${task.cpus} -o namecollated.bam $bam &>> $lf; zcat -f ${params.gDTUANNO} > $fn/annotation.gtf 2>> $lf; ${params.gCOUNTBIN} --threads ${task.cpus} --genome-alignments namecollated.bam --annotation $fn/annotation.gtf --genome-fasta ${params.gDTUREF} ${params.gCOUNTPARAMS} --output $prefix &>> $lf && rm -f namecollated.bam $fn/annotation.gtf && python3 ${params.gBINS}/Analysis/oarfish_to_salmon.py -i ${prefix}.quant -o $oz &>> $lf && ln -fs $oz $ol + mkdir -p $fn; : > $lf; if [ ! -s $anno ]; then echo "ERROR: decompressed annotation $anno is empty, check source file and filesystem" >> $lf; exit 1; fi; samtools collate -@ ${task.cpus} -o namecollated.bam $bam &>> $lf; ${params.gCOUNTBIN} --threads ${task.cpus} --genome-alignments namecollated.bam --annotation $anno --genome-fasta ${params.gDTUREF} ${params.gCOUNTPARAMS} --output $prefix &>> $lf && rm -f namecollated.bam && python3 ${params.gBINS}/Analysis/oarfish_to_salmon.py -i ${prefix}.quant -o $oz &>> $lf && ln -fs $oz $ol """ } process create_summary_snippet{ @@ -344,7 +365,9 @@ workflow DTU{ element -> return "${workflow.workDir}/../MAPPED/${params.gCOMBO}/"+element+"_mapped_sorted.bam" } mapsamples_ch = Channel.fromPath(MAPPEDSAMPLES.sort()) - oarfish_dtu_quant(mapsamples_ch.collate(1)) + decompress_annotation_dtu(annofile) + annogtf_dtu = decompress_annotation_dtu.out.anno_plain + oarfish_dtu_quant(annogtf_dtu.combine(mapsamples_ch.collate(1))) quantcounts = oarfish_dtu_quant.out.counts } prepare_dtu_annotation() diff --git a/tests/golden/config_DTU/smk/allconditions_fastqc-trimgalore-picard-star_DTU_dexseq_subsnake.smk b/tests/golden/config_DTU/smk/allconditions_fastqc-trimgalore-picard-star_DTU_dexseq_subsnake.smk index e1aab12d..868fbe7a 100644 --- a/tests/golden/config_DTU/smk/allconditions_fastqc-trimgalore-picard-star_DTU_dexseq_subsnake.smk +++ b/tests/golden/config_DTU/smk/allconditions_fastqc-trimgalore-picard-star_DTU_dexseq_subsnake.smk @@ -293,6 +293,10 @@ if 'FUSIONS' in config: ANNOTATION = FANNO else: ANNOTATION = ANNO.get('GTF') if 'GTF' in ANNO and ANNO.get('GTF') != '' else ANNO.get('GFF') # by default GTF format will be used +if ANNOTATION and isinstance(ANNOTATION, str): + # Resolve once to an absolute path so every rule's shell command can read + # it reliably regardless of its own working directory (e.g. containers). + ANNOTATION = os.path.abspath(ANNOTATION) combo = '' if ANNOTATION: log.info(logid+f'Using COMBO: {combo}, REF: {REFERENCE}, ANNO: {ANNOTATION}, DECOY: {DECOY}, INDEX: {INDEX}, INDEX2: {INDEX2}, PREFIX: {PREFIX}, IP: {IP}, SETUP: {SETUP}, SEQUENCING: {SEQUENCING}, STRANDED: {stranded}, PAIRED: {paired}, RUNDEDUP: {rundedup}, PREDEDUP: {prededup}, SAMPLES: {SAMPLES}, BINS: {BINS}, MAXTHREAD: {MAXTHREAD}') @@ -404,9 +408,17 @@ else: # is needed, output is normalized to the salmon quant.sf schema so # downstream DTU tooling (tximport type="salmon", build_DTU_table.py, # SPIT.py) does not need to change. + rule decompress_annotation: + input: anno = ANNOTATION + output: anno = temp("DTU/{combo}/"+COUNTENV+"/annotation.gtf") + log: "LOGS/{combo}/DTU/"+COUNTENV+"/decompress_annotation.log" + conda: "/envs/"+COUNTENV+".yaml" + container: "oras://ghcr.io/jfallmann/monsda:"+COUNTENV+""+"-VERSION" + threads: 1 + shell: "zcat -f {input.anno} > {output.anno} 2> {log}; if [ ! -s {output.anno} ]; then echo 'ERROR: decompressed annotation from {input.anno} is empty, check source file and filesystem' >> {log}; exit 1; fi" rule mapping: input: bam = expand("MAPPED/{scombo}/{{file}}_mapped_sorted.bam", scombo=scombo), - anno = ANNOTATION, + anno = rules.decompress_annotation.output.anno, fa = REFERENCE output: cnts = report("DTU/{combo}/"+COUNTENV+"/{file}_counts.sf.gz", category="COUNTING"), ctsdir = report(directory("DTU/{combo}/"+COUNTENV+"/{file}"), category="COUNTING") @@ -417,11 +429,10 @@ else: params: cpara = lambda wildcards: tool_params(wildcards.file, None, config, 'DTU', DTUENV)['OPTIONS'].get('QUANT', ""), mapp = COUNTBIN, nbam = lambda wildcards, output: os.path.join(str(output.ctsdir), "namecollated.bam"), - anno = lambda wildcards, output: os.path.join(str(output.ctsdir), "annotation.gtf"), prefix = lambda wildcards, output: os.path.join(str(output.ctsdir), COUNTENV), linksf = lambda wildcards, output: str(os.path.abspath(output.ctsdir)), bins = BINS - shell: "set +euo pipefail; mkdir -p {output.ctsdir}; samtools collate -@ {threads} -o {params.nbam} {input.bam} &>> {log}; zcat -f {input.anno} > {params.anno} 2>> {log}; {params.mapp} --threads {threads} --genome-alignments {params.nbam} --annotation {params.anno} --genome-fasta {input.fa} {params.cpara} --output {params.prefix} &>> {log} && rm -f {params.nbam} {params.anno} && python3 {params.bins}/Analysis/oarfish_to_salmon.py -i {params.prefix}.quant -o {output.ctsdir}/quant.sf.gz &>> {log} ; ln -fs {params.linksf}/quant.sf.gz {output.cnts} &>> {log}" + shell: "set +euo pipefail; mkdir -p {output.ctsdir}; samtools collate -@ {threads} -o {params.nbam} {input.bam} &>> {log}; {params.mapp} --threads {threads} --genome-alignments {params.nbam} --annotation {input.anno} --genome-fasta {input.fa} {params.cpara} --output {params.prefix} &>> {log} && rm -f {params.nbam} && python3 {params.bins}/Analysis/oarfish_to_salmon.py -i {params.prefix}.quant -o {output.ctsdir}/quant.sf.gz &>> {log} ; ln -fs {params.linksf}/quant.sf.gz {output.cnts} &>> {log}" if runterminus: rule terminus_group: input: ctsdir = rules.mapping.output.ctsdir diff --git a/tests/golden/config_DTU/smk/allconditions_fastqc-trimgalore-picard-star_DTU_drimseq_subsnake.smk b/tests/golden/config_DTU/smk/allconditions_fastqc-trimgalore-picard-star_DTU_drimseq_subsnake.smk index c116c2f3..1cc1c0b5 100644 --- a/tests/golden/config_DTU/smk/allconditions_fastqc-trimgalore-picard-star_DTU_drimseq_subsnake.smk +++ b/tests/golden/config_DTU/smk/allconditions_fastqc-trimgalore-picard-star_DTU_drimseq_subsnake.smk @@ -293,6 +293,10 @@ if 'FUSIONS' in config: ANNOTATION = FANNO else: ANNOTATION = ANNO.get('GTF') if 'GTF' in ANNO and ANNO.get('GTF') != '' else ANNO.get('GFF') # by default GTF format will be used +if ANNOTATION and isinstance(ANNOTATION, str): + # Resolve once to an absolute path so every rule's shell command can read + # it reliably regardless of its own working directory (e.g. containers). + ANNOTATION = os.path.abspath(ANNOTATION) combo = '' if ANNOTATION: log.info(logid+f'Using COMBO: {combo}, REF: {REFERENCE}, ANNO: {ANNOTATION}, DECOY: {DECOY}, INDEX: {INDEX}, INDEX2: {INDEX2}, PREFIX: {PREFIX}, IP: {IP}, SETUP: {SETUP}, SEQUENCING: {SEQUENCING}, STRANDED: {stranded}, PAIRED: {paired}, RUNDEDUP: {rundedup}, PREDEDUP: {prededup}, SAMPLES: {SAMPLES}, BINS: {BINS}, MAXTHREAD: {MAXTHREAD}') @@ -414,9 +418,17 @@ else: # is needed, output is normalized to the salmon quant.sf schema so # downstream DTU tooling (tximport type="salmon", build_DTU_table.py, # SPIT.py) does not need to change. + rule decompress_annotation: + input: anno = ANNOTATION + output: anno = temp("DTU/{combo}/"+COUNTENV+"/annotation.gtf") + log: "LOGS/{combo}/DTU/"+COUNTENV+"/decompress_annotation.log" + conda: "/envs/"+COUNTENV+".yaml" + container: "oras://ghcr.io/jfallmann/monsda:"+COUNTENV+""+"-VERSION" + threads: 1 + shell: "zcat -f {input.anno} > {output.anno} 2> {log}; if [ ! -s {output.anno} ]; then echo 'ERROR: decompressed annotation from {input.anno} is empty, check source file and filesystem' >> {log}; exit 1; fi" rule mapping: input: bam = expand("MAPPED/{scombo}/{{file}}_mapped_sorted.bam", scombo=scombo), - anno = ANNOTATION, + anno = rules.decompress_annotation.output.anno, fa = REFERENCE output: cnts = report("DTU/{combo}/"+COUNTENV+"/{file}_counts.sf.gz", category="COUNTING"), ctsdir = report(directory("DTU/{combo}/"+COUNTENV+"/{file}"), category="COUNTING") @@ -427,11 +439,10 @@ else: params: cpara = lambda wildcards: tool_params(wildcards.file, None, config, 'DTU', DTUENV)['OPTIONS'].get('QUANT', ""), mapp = COUNTBIN, nbam = lambda wildcards, output: os.path.join(str(output.ctsdir), "namecollated.bam"), - anno = lambda wildcards, output: os.path.join(str(output.ctsdir), "annotation.gtf"), prefix = lambda wildcards, output: os.path.join(str(output.ctsdir), COUNTENV), linksf = lambda wildcards, output: str(os.path.abspath(output.ctsdir)), bins = BINS - shell: "set +euo pipefail; mkdir -p {output.ctsdir}; samtools collate -@ {threads} -o {params.nbam} {input.bam} &>> {log}; zcat -f {input.anno} > {params.anno} 2>> {log}; {params.mapp} --threads {threads} --genome-alignments {params.nbam} --annotation {params.anno} --genome-fasta {input.fa} {params.cpara} --output {params.prefix} &>> {log} && rm -f {params.nbam} {params.anno} && python3 {params.bins}/Analysis/oarfish_to_salmon.py -i {params.prefix}.quant -o {output.ctsdir}/quant.sf.gz &>> {log} ; ln -fs {params.linksf}/quant.sf.gz {output.cnts} &>> {log}" + shell: "set +euo pipefail; mkdir -p {output.ctsdir}; samtools collate -@ {threads} -o {params.nbam} {input.bam} &>> {log}; {params.mapp} --threads {threads} --genome-alignments {params.nbam} --annotation {input.anno} --genome-fasta {input.fa} {params.cpara} --output {params.prefix} &>> {log} && rm -f {params.nbam} && python3 {params.bins}/Analysis/oarfish_to_salmon.py -i {params.prefix}.quant -o {output.ctsdir}/quant.sf.gz &>> {log} ; ln -fs {params.linksf}/quant.sf.gz {output.cnts} &>> {log}" if runterminus: rule terminus_group: input: ctsdir = rules.mapping.output.ctsdir diff --git a/tests/golden/config_DTU/smk/allconditions_fastqc-trimgalore-picard-star_DTU_spit_subsnake.smk b/tests/golden/config_DTU/smk/allconditions_fastqc-trimgalore-picard-star_DTU_spit_subsnake.smk index 9e592ea8..4c4f8440 100644 --- a/tests/golden/config_DTU/smk/allconditions_fastqc-trimgalore-picard-star_DTU_spit_subsnake.smk +++ b/tests/golden/config_DTU/smk/allconditions_fastqc-trimgalore-picard-star_DTU_spit_subsnake.smk @@ -293,6 +293,10 @@ if 'FUSIONS' in config: ANNOTATION = FANNO else: ANNOTATION = ANNO.get('GTF') if 'GTF' in ANNO and ANNO.get('GTF') != '' else ANNO.get('GFF') # by default GTF format will be used +if ANNOTATION and isinstance(ANNOTATION, str): + # Resolve once to an absolute path so every rule's shell command can read + # it reliably regardless of its own working directory (e.g. containers). + ANNOTATION = os.path.abspath(ANNOTATION) combo = '' if ANNOTATION: log.info(logid+f'Using COMBO: {combo}, REF: {REFERENCE}, ANNO: {ANNOTATION}, DECOY: {DECOY}, INDEX: {INDEX}, INDEX2: {INDEX2}, PREFIX: {PREFIX}, IP: {IP}, SETUP: {SETUP}, SEQUENCING: {SEQUENCING}, STRANDED: {stranded}, PAIRED: {paired}, RUNDEDUP: {rundedup}, PREDEDUP: {prededup}, SAMPLES: {SAMPLES}, BINS: {BINS}, MAXTHREAD: {MAXTHREAD}') @@ -411,9 +415,17 @@ else: # is needed, output is normalized to the salmon quant.sf schema so # downstream DTU tooling (tximport type="salmon", build_DTU_table.py, # SPIT.py) does not need to change. + rule decompress_annotation: + input: anno = ANNOTATION + output: anno = temp("DTU/{combo}/"+COUNTENV+"/annotation.gtf") + log: "LOGS/{combo}/DTU/"+COUNTENV+"/decompress_annotation.log" + conda: "/envs/"+COUNTENV+".yaml" + container: "oras://ghcr.io/jfallmann/monsda:"+COUNTENV+""+"-VERSION" + threads: 1 + shell: "zcat -f {input.anno} > {output.anno} 2> {log}; if [ ! -s {output.anno} ]; then echo 'ERROR: decompressed annotation from {input.anno} is empty, check source file and filesystem' >> {log}; exit 1; fi" rule mapping: input: bam = expand("MAPPED/{scombo}/{{file}}_mapped_sorted.bam", scombo=scombo), - anno = ANNOTATION, + anno = rules.decompress_annotation.output.anno, fa = REFERENCE output: cnts = report("DTU/{combo}/"+COUNTENV+"/{file}_counts.sf.gz", category="COUNTING"), ctsdir = report(directory("DTU/{combo}/"+COUNTENV+"/{file}"), category="COUNTING") @@ -424,11 +436,10 @@ else: params: cpara = lambda wildcards: tool_params(wildcards.file, None, config, 'DTU', DTUENV)['OPTIONS'].get('QUANT', ""), mapp = COUNTBIN, nbam = lambda wildcards, output: os.path.join(str(output.ctsdir), "namecollated.bam"), - anno = lambda wildcards, output: os.path.join(str(output.ctsdir), "annotation.gtf"), prefix = lambda wildcards, output: os.path.join(str(output.ctsdir), COUNTENV), linksf = lambda wildcards, output: str(os.path.abspath(output.ctsdir)), bins = BINS - shell: "set +euo pipefail; mkdir -p {output.ctsdir}; samtools collate -@ {threads} -o {params.nbam} {input.bam} &>> {log}; zcat -f {input.anno} > {params.anno} 2>> {log}; {params.mapp} --threads {threads} --genome-alignments {params.nbam} --annotation {params.anno} --genome-fasta {input.fa} {params.cpara} --output {params.prefix} &>> {log} && rm -f {params.nbam} {params.anno} && python3 {params.bins}/Analysis/oarfish_to_salmon.py -i {params.prefix}.quant -o {output.ctsdir}/quant.sf.gz &>> {log} ; ln -fs {params.linksf}/quant.sf.gz {output.cnts} &>> {log}" + shell: "set +euo pipefail; mkdir -p {output.ctsdir}; samtools collate -@ {threads} -o {params.nbam} {input.bam} &>> {log}; {params.mapp} --threads {threads} --genome-alignments {params.nbam} --annotation {input.anno} --genome-fasta {input.fa} {params.cpara} --output {params.prefix} &>> {log} && rm -f {params.nbam} && python3 {params.bins}/Analysis/oarfish_to_salmon.py -i {params.prefix}.quant -o {output.ctsdir}/quant.sf.gz &>> {log} ; ln -fs {params.linksf}/quant.sf.gz {output.cnts} &>> {log}" if runterminus: rule terminus_group: input: ctsdir = rules.mapping.output.ctsdir diff --git a/tests/golden/config_FUSIONS/smk/allconditions_fastqc-trimgalore-picard-star_FUSIONS_starfusion_subsnake.smk b/tests/golden/config_FUSIONS/smk/allconditions_fastqc-trimgalore-picard-star_FUSIONS_starfusion_subsnake.smk index cde1eadf..6aba6a5b 100644 --- a/tests/golden/config_FUSIONS/smk/allconditions_fastqc-trimgalore-picard-star_FUSIONS_starfusion_subsnake.smk +++ b/tests/golden/config_FUSIONS/smk/allconditions_fastqc-trimgalore-picard-star_FUSIONS_starfusion_subsnake.smk @@ -293,6 +293,10 @@ if 'FUSIONS' in config: ANNOTATION = FANNO else: ANNOTATION = ANNO.get('GTF') if 'GTF' in ANNO and ANNO.get('GTF') != '' else ANNO.get('GFF') # by default GTF format will be used +if ANNOTATION and isinstance(ANNOTATION, str): + # Resolve once to an absolute path so every rule's shell command can read + # it reliably regardless of its own working directory (e.g. containers). + ANNOTATION = os.path.abspath(ANNOTATION) combo = '' if ANNOTATION: log.info(logid+f'Using COMBO: {combo}, REF: {REFERENCE}, ANNO: {ANNOTATION}, DECOY: {DECOY}, INDEX: {INDEX}, INDEX2: {INDEX2}, PREFIX: {PREFIX}, IP: {IP}, SETUP: {SETUP}, SEQUENCING: {SEQUENCING}, STRANDED: {stranded}, PAIRED: {paired}, RUNDEDUP: {rundedup}, PREDEDUP: {prededup}, SAMPLES: {SAMPLES}, BINS: {BINS}, MAXTHREAD: {MAXTHREAD}') diff --git a/tests/golden/config_Test/smk/Test_fastqc-bbduk-picard-bwa_subsnake.smk b/tests/golden/config_Test/smk/Test_fastqc-bbduk-picard-bwa_subsnake.smk index 7d78489f..f30e41e7 100644 --- a/tests/golden/config_Test/smk/Test_fastqc-bbduk-picard-bwa_subsnake.smk +++ b/tests/golden/config_Test/smk/Test_fastqc-bbduk-picard-bwa_subsnake.smk @@ -293,6 +293,10 @@ if 'FUSIONS' in config: ANNOTATION = FANNO else: ANNOTATION = ANNO.get('GTF') if 'GTF' in ANNO and ANNO.get('GTF') != '' else ANNO.get('GFF') # by default GTF format will be used +if ANNOTATION and isinstance(ANNOTATION, str): + # Resolve once to an absolute path so every rule's shell command can read + # it reliably regardless of its own working directory (e.g. containers). + ANNOTATION = os.path.abspath(ANNOTATION) combo = '' if ANNOTATION: log.info(logid+f'Using COMBO: {combo}, REF: {REFERENCE}, ANNO: {ANNOTATION}, DECOY: {DECOY}, INDEX: {INDEX}, INDEX2: {INDEX2}, PREFIX: {PREFIX}, IP: {IP}, SETUP: {SETUP}, SEQUENCING: {SEQUENCING}, STRANDED: {stranded}, PAIRED: {paired}, RUNDEDUP: {rundedup}, PREDEDUP: {prededup}, SAMPLES: {SAMPLES}, BINS: {BINS}, MAXTHREAD: {MAXTHREAD}') diff --git a/tests/golden/config_Test/smk/Test_fastqc-bbduk-picard-rustar_subsnake.smk b/tests/golden/config_Test/smk/Test_fastqc-bbduk-picard-rustar_subsnake.smk index e6c2b076..de09320b 100644 --- a/tests/golden/config_Test/smk/Test_fastqc-bbduk-picard-rustar_subsnake.smk +++ b/tests/golden/config_Test/smk/Test_fastqc-bbduk-picard-rustar_subsnake.smk @@ -293,6 +293,10 @@ if 'FUSIONS' in config: ANNOTATION = FANNO else: ANNOTATION = ANNO.get('GTF') if 'GTF' in ANNO and ANNO.get('GTF') != '' else ANNO.get('GFF') # by default GTF format will be used +if ANNOTATION and isinstance(ANNOTATION, str): + # Resolve once to an absolute path so every rule's shell command can read + # it reliably regardless of its own working directory (e.g. containers). + ANNOTATION = os.path.abspath(ANNOTATION) combo = '' if ANNOTATION: log.info(logid+f'Using COMBO: {combo}, REF: {REFERENCE}, ANNO: {ANNOTATION}, DECOY: {DECOY}, INDEX: {INDEX}, INDEX2: {INDEX2}, PREFIX: {PREFIX}, IP: {IP}, SETUP: {SETUP}, SEQUENCING: {SEQUENCING}, STRANDED: {stranded}, PAIRED: {paired}, RUNDEDUP: {rundedup}, PREDEDUP: {prededup}, SAMPLES: {SAMPLES}, BINS: {BINS}, MAXTHREAD: {MAXTHREAD}') diff --git a/tests/golden/config_Test/smk/Test_fastqc-bbduk-picard-star_subsnake.smk b/tests/golden/config_Test/smk/Test_fastqc-bbduk-picard-star_subsnake.smk index a97ad49a..b80e8001 100644 --- a/tests/golden/config_Test/smk/Test_fastqc-bbduk-picard-star_subsnake.smk +++ b/tests/golden/config_Test/smk/Test_fastqc-bbduk-picard-star_subsnake.smk @@ -293,6 +293,10 @@ if 'FUSIONS' in config: ANNOTATION = FANNO else: ANNOTATION = ANNO.get('GTF') if 'GTF' in ANNO and ANNO.get('GTF') != '' else ANNO.get('GFF') # by default GTF format will be used +if ANNOTATION and isinstance(ANNOTATION, str): + # Resolve once to an absolute path so every rule's shell command can read + # it reliably regardless of its own working directory (e.g. containers). + ANNOTATION = os.path.abspath(ANNOTATION) combo = '' if ANNOTATION: log.info(logid+f'Using COMBO: {combo}, REF: {REFERENCE}, ANNO: {ANNOTATION}, DECOY: {DECOY}, INDEX: {INDEX}, INDEX2: {INDEX2}, PREFIX: {PREFIX}, IP: {IP}, SETUP: {SETUP}, SEQUENCING: {SEQUENCING}, STRANDED: {stranded}, PAIRED: {paired}, RUNDEDUP: {rundedup}, PREDEDUP: {prededup}, SAMPLES: {SAMPLES}, BINS: {BINS}, MAXTHREAD: {MAXTHREAD}') diff --git a/tests/golden/config_Test/smk/Test_fastqc-cutadapt-picard-bwa_subsnake.smk b/tests/golden/config_Test/smk/Test_fastqc-cutadapt-picard-bwa_subsnake.smk index c75b659a..c77aa534 100644 --- a/tests/golden/config_Test/smk/Test_fastqc-cutadapt-picard-bwa_subsnake.smk +++ b/tests/golden/config_Test/smk/Test_fastqc-cutadapt-picard-bwa_subsnake.smk @@ -293,6 +293,10 @@ if 'FUSIONS' in config: ANNOTATION = FANNO else: ANNOTATION = ANNO.get('GTF') if 'GTF' in ANNO and ANNO.get('GTF') != '' else ANNO.get('GFF') # by default GTF format will be used +if ANNOTATION and isinstance(ANNOTATION, str): + # Resolve once to an absolute path so every rule's shell command can read + # it reliably regardless of its own working directory (e.g. containers). + ANNOTATION = os.path.abspath(ANNOTATION) combo = '' if ANNOTATION: log.info(logid+f'Using COMBO: {combo}, REF: {REFERENCE}, ANNO: {ANNOTATION}, DECOY: {DECOY}, INDEX: {INDEX}, INDEX2: {INDEX2}, PREFIX: {PREFIX}, IP: {IP}, SETUP: {SETUP}, SEQUENCING: {SEQUENCING}, STRANDED: {stranded}, PAIRED: {paired}, RUNDEDUP: {rundedup}, PREDEDUP: {prededup}, SAMPLES: {SAMPLES}, BINS: {BINS}, MAXTHREAD: {MAXTHREAD}') diff --git a/tests/golden/config_Test/smk/Test_fastqc-cutadapt-picard-rustar_subsnake.smk b/tests/golden/config_Test/smk/Test_fastqc-cutadapt-picard-rustar_subsnake.smk index f630ccd6..15ec53d0 100644 --- a/tests/golden/config_Test/smk/Test_fastqc-cutadapt-picard-rustar_subsnake.smk +++ b/tests/golden/config_Test/smk/Test_fastqc-cutadapt-picard-rustar_subsnake.smk @@ -293,6 +293,10 @@ if 'FUSIONS' in config: ANNOTATION = FANNO else: ANNOTATION = ANNO.get('GTF') if 'GTF' in ANNO and ANNO.get('GTF') != '' else ANNO.get('GFF') # by default GTF format will be used +if ANNOTATION and isinstance(ANNOTATION, str): + # Resolve once to an absolute path so every rule's shell command can read + # it reliably regardless of its own working directory (e.g. containers). + ANNOTATION = os.path.abspath(ANNOTATION) combo = '' if ANNOTATION: log.info(logid+f'Using COMBO: {combo}, REF: {REFERENCE}, ANNO: {ANNOTATION}, DECOY: {DECOY}, INDEX: {INDEX}, INDEX2: {INDEX2}, PREFIX: {PREFIX}, IP: {IP}, SETUP: {SETUP}, SEQUENCING: {SEQUENCING}, STRANDED: {stranded}, PAIRED: {paired}, RUNDEDUP: {rundedup}, PREDEDUP: {prededup}, SAMPLES: {SAMPLES}, BINS: {BINS}, MAXTHREAD: {MAXTHREAD}') diff --git a/tests/golden/config_Test/smk/Test_fastqc-cutadapt-picard-star_subsnake.smk b/tests/golden/config_Test/smk/Test_fastqc-cutadapt-picard-star_subsnake.smk index 8048518c..af644e01 100644 --- a/tests/golden/config_Test/smk/Test_fastqc-cutadapt-picard-star_subsnake.smk +++ b/tests/golden/config_Test/smk/Test_fastqc-cutadapt-picard-star_subsnake.smk @@ -293,6 +293,10 @@ if 'FUSIONS' in config: ANNOTATION = FANNO else: ANNOTATION = ANNO.get('GTF') if 'GTF' in ANNO and ANNO.get('GTF') != '' else ANNO.get('GFF') # by default GTF format will be used +if ANNOTATION and isinstance(ANNOTATION, str): + # Resolve once to an absolute path so every rule's shell command can read + # it reliably regardless of its own working directory (e.g. containers). + ANNOTATION = os.path.abspath(ANNOTATION) combo = '' if ANNOTATION: log.info(logid+f'Using COMBO: {combo}, REF: {REFERENCE}, ANNO: {ANNOTATION}, DECOY: {DECOY}, INDEX: {INDEX}, INDEX2: {INDEX2}, PREFIX: {PREFIX}, IP: {IP}, SETUP: {SETUP}, SEQUENCING: {SEQUENCING}, STRANDED: {stranded}, PAIRED: {paired}, RUNDEDUP: {rundedup}, PREDEDUP: {prededup}, SAMPLES: {SAMPLES}, BINS: {BINS}, MAXTHREAD: {MAXTHREAD}') diff --git a/tests/golden/config_Test/smk/Test_fastqc-trimgalore-picard-bwa_subsnake.smk b/tests/golden/config_Test/smk/Test_fastqc-trimgalore-picard-bwa_subsnake.smk index c31378a3..514ac158 100644 --- a/tests/golden/config_Test/smk/Test_fastqc-trimgalore-picard-bwa_subsnake.smk +++ b/tests/golden/config_Test/smk/Test_fastqc-trimgalore-picard-bwa_subsnake.smk @@ -293,6 +293,10 @@ if 'FUSIONS' in config: ANNOTATION = FANNO else: ANNOTATION = ANNO.get('GTF') if 'GTF' in ANNO and ANNO.get('GTF') != '' else ANNO.get('GFF') # by default GTF format will be used +if ANNOTATION and isinstance(ANNOTATION, str): + # Resolve once to an absolute path so every rule's shell command can read + # it reliably regardless of its own working directory (e.g. containers). + ANNOTATION = os.path.abspath(ANNOTATION) combo = '' if ANNOTATION: log.info(logid+f'Using COMBO: {combo}, REF: {REFERENCE}, ANNO: {ANNOTATION}, DECOY: {DECOY}, INDEX: {INDEX}, INDEX2: {INDEX2}, PREFIX: {PREFIX}, IP: {IP}, SETUP: {SETUP}, SEQUENCING: {SEQUENCING}, STRANDED: {stranded}, PAIRED: {paired}, RUNDEDUP: {rundedup}, PREDEDUP: {prededup}, SAMPLES: {SAMPLES}, BINS: {BINS}, MAXTHREAD: {MAXTHREAD}') diff --git a/tests/golden/config_Test/smk/Test_fastqc-trimgalore-picard-rustar_subsnake.smk b/tests/golden/config_Test/smk/Test_fastqc-trimgalore-picard-rustar_subsnake.smk index 3c0358ab..7ba6f00f 100644 --- a/tests/golden/config_Test/smk/Test_fastqc-trimgalore-picard-rustar_subsnake.smk +++ b/tests/golden/config_Test/smk/Test_fastqc-trimgalore-picard-rustar_subsnake.smk @@ -293,6 +293,10 @@ if 'FUSIONS' in config: ANNOTATION = FANNO else: ANNOTATION = ANNO.get('GTF') if 'GTF' in ANNO and ANNO.get('GTF') != '' else ANNO.get('GFF') # by default GTF format will be used +if ANNOTATION and isinstance(ANNOTATION, str): + # Resolve once to an absolute path so every rule's shell command can read + # it reliably regardless of its own working directory (e.g. containers). + ANNOTATION = os.path.abspath(ANNOTATION) combo = '' if ANNOTATION: log.info(logid+f'Using COMBO: {combo}, REF: {REFERENCE}, ANNO: {ANNOTATION}, DECOY: {DECOY}, INDEX: {INDEX}, INDEX2: {INDEX2}, PREFIX: {PREFIX}, IP: {IP}, SETUP: {SETUP}, SEQUENCING: {SEQUENCING}, STRANDED: {stranded}, PAIRED: {paired}, RUNDEDUP: {rundedup}, PREDEDUP: {prededup}, SAMPLES: {SAMPLES}, BINS: {BINS}, MAXTHREAD: {MAXTHREAD}') diff --git a/tests/golden/config_Test/smk/Test_fastqc-trimgalore-picard-star_subsnake.smk b/tests/golden/config_Test/smk/Test_fastqc-trimgalore-picard-star_subsnake.smk index e336b7fe..800e9ebd 100644 --- a/tests/golden/config_Test/smk/Test_fastqc-trimgalore-picard-star_subsnake.smk +++ b/tests/golden/config_Test/smk/Test_fastqc-trimgalore-picard-star_subsnake.smk @@ -293,6 +293,10 @@ if 'FUSIONS' in config: ANNOTATION = FANNO else: ANNOTATION = ANNO.get('GTF') if 'GTF' in ANNO and ANNO.get('GTF') != '' else ANNO.get('GFF') # by default GTF format will be used +if ANNOTATION and isinstance(ANNOTATION, str): + # Resolve once to an absolute path so every rule's shell command can read + # it reliably regardless of its own working directory (e.g. containers). + ANNOTATION = os.path.abspath(ANNOTATION) combo = '' if ANNOTATION: log.info(logid+f'Using COMBO: {combo}, REF: {REFERENCE}, ANNO: {ANNOTATION}, DECOY: {DECOY}, INDEX: {INDEX}, INDEX2: {INDEX2}, PREFIX: {PREFIX}, IP: {IP}, SETUP: {SETUP}, SEQUENCING: {SEQUENCING}, STRANDED: {stranded}, PAIRED: {paired}, RUNDEDUP: {rundedup}, PREDEDUP: {prededup}, SAMPLES: {SAMPLES}, BINS: {BINS}, MAXTHREAD: {MAXTHREAD}') diff --git a/workflows/deseq2_DE.nf b/workflows/deseq2_DE.nf index f15d570c..3f1e9c13 100644 --- a/workflows/deseq2_DE.nf +++ b/workflows/deseq2_DE.nf @@ -24,6 +24,30 @@ QUANTPARAMS = get_always('deseq2_DE_params_QUANT') ?: '' //DE PROCESSES +process decompress_annotation_deseq{ + conda "$COUNTENV"+".yaml" + container "oras://jfallmann/monsda:"+"$COUNTENV" + cpus 1 + cache 'lenient' + + publishDir "${workflow.workDir}/../" , mode: 'link', + saveAs: {filename -> + if (filename.indexOf("log") > 0) "LOGS/${SCOMBO}/DE/deseq2/decompress_annotation.log" + } + + input: + path anno + + output: + path "annotation.gtf", emit: anno_plain + path "log", emit: log + + script: + """ + zcat -f $anno > annotation.gtf 2> log; if [ ! -s annotation.gtf ]; then echo "ERROR: decompressed annotation from $anno is empty, check source file and filesystem" >> log; exit 1; fi + """ +} + process featurecount_deseq{ conda "$COUNTENV"+".yaml" container "oras://jfallmann/monsda:"+"$COUNTENV" @@ -69,7 +93,7 @@ process featurecount_deseq{ stranded = '' } """ - mkdir -p TMP; zcat -f $anno > annotation.gtf 2> $ol; $COUNTBIN -T ${task.cpus} $COUNTPARAMS $pair $stranded -a annotation.gtf -o tmpcts $reads 2>> $ol && rm -f annotation.gtf && head -n2 tmpcts |gzip > $oc && export LC_ALL=C; tail -n+3 tmpcts|sort --parallel=${task.cpus} -S $sortmem -T TMP -k1,1 -k2,2n -k3,3n -u |gzip >> $oc 2>> $ol && mv tmpcts.summary $os + mkdir -p TMP; : > $ol; if [ ! -s $anno ]; then echo "ERROR: decompressed annotation $anno is empty, check source file and filesystem" >> $ol; exit 1; fi; $COUNTBIN -T ${task.cpus} $COUNTPARAMS $pair $stranded -a $anno -o tmpcts $reads 2>> $ol && head -n2 tmpcts |gzip > $oc && export LC_ALL=C; tail -n+3 tmpcts|sort --parallel=${task.cpus} -S $sortmem -T TMP -k1,1 -k2,2n -k3,3n -u |gzip >> $oc 2>> $ol && mv tmpcts.summary $os """ } @@ -103,7 +127,7 @@ process oarfish_gene_count{ os = fn+".counts.summary" ol = fn+".log" """ - samtools collate -@ ${task.cpus} -o namecollated.bam $bam &>> $ol; zcat -f $anno > annotation.gtf 2>> $ol; $COUNTBIN --threads ${task.cpus} --genome-alignments namecollated.bam --annotation annotation.gtf --genome-fasta $DEREF $QUANTPARAMS --output $fn &>> $ol && rm -f namecollated.bam annotation.gtf && python3 $BINS/Analysis/oarfish_to_featurecounts.py -i ${fn}.quant -g $anno -s ${fn}_mapped_sorted.bam -o $oc 2>> $ol && touch $os + : > $ol; if [ ! -s $anno ]; then echo "ERROR: decompressed annotation $anno is empty, check source file and filesystem" >> $ol; exit 1; fi; samtools collate -@ ${task.cpus} -o namecollated.bam $bam &>> $ol; $COUNTBIN --threads ${task.cpus} --genome-alignments namecollated.bam --annotation $anno --genome-fasta $DEREF $QUANTPARAMS --output $fn &>> $ol && rm -f namecollated.bam && python3 $BINS/Analysis/oarfish_to_featurecounts.py -i ${fn}.quant -g $anno -s ${fn}_mapped_sorted.bam -o $oc 2>> $ol && touch $os """ } @@ -282,6 +306,9 @@ workflow DE{ annofile = Channel.fromPath(DEANNO) + decompress_annotation_deseq(annofile) + annogtf = decompress_annotation_deseq.out.anno_plain + if (COUNTTOOL == 'featureCounts'){ if (RUNDEDUP){ MAPPEDSAMPLES = LONGSAMPLES.collect{ @@ -296,7 +323,7 @@ workflow DE{ mapsamples_ch = Channel.fromPath(MAPPEDSAMPLES.sort()) //mapsamples_ch.subscribe { println "MAP: $it \t COMBO: ${COMBO} SCOMBO: ${SCOMBO} LONG: ${LONGSAMPLES}" } - featurecount_deseq(annofile.combine(mapsamples_ch.collate(1))) + featurecount_deseq(annogtf.combine(mapsamples_ch.collate(1))) countcts = featurecount_deseq.out.fc_cts }else{ MAPPEDSAMPLES = LONGSAMPLES.collect{ @@ -304,7 +331,7 @@ workflow DE{ } mapsamples_ch = Channel.fromPath(MAPPEDSAMPLES.sort()) - oarfish_gene_count(annofile.combine(mapsamples_ch.collate(1))) + oarfish_gene_count(annogtf.combine(mapsamples_ch.collate(1))) countcts = oarfish_gene_count.out.fc_cts } diff --git a/workflows/deseq2_DE.smk b/workflows/deseq2_DE.smk index 296d55db..2d48f6da 100644 --- a/workflows/deseq2_DE.smk +++ b/workflows/deseq2_DE.smk @@ -30,9 +30,19 @@ rule themall: Rmd = expand("REPORTS/SUMMARY/RmdSnippets/{combo}.Rmd", combo=combo), postde = expand("POSTDE/{combo}", combo=combo) if postde_enabled else [] +rule decompress_annotation: + input: anno = ANNOTATION + output: anno = temp("DE/{combo}/Featurecounts/annotation.gtf") + log: "LOGS/{combo}/DE/deseq2/decompress_annotation.log" + conda: ""+COUNTENV+".yaml" + container: "oras://jfallmann/monsda:"+COUNTENV+"" + threads: 1 + shell: "zcat -f {input.anno} > {output.anno} 2> {log}; if [ ! -s {output.anno} ]; then echo 'ERROR: decompressed annotation from {input.anno} is empty, check source file and filesystem' >> {log}; exit 1; fi" + if COUNTTOOL == 'featureCounts': rule featurecount_unique: - input: reads = expand("MAPPED/{scombo}/{{file}}_mapped_sorted_unique.bam", scombo=scombo) if not usededup else expand("MAPPED/{scombo}/{{file}}_mapped_sorted_unique_dedup.bam", scombo=scombo) + input: reads = expand("MAPPED/{scombo}/{{file}}_mapped_sorted_unique.bam", scombo=scombo) if not usededup else expand("MAPPED/{scombo}/{{file}}_mapped_sorted_unique_dedup.bam", scombo=scombo), + anno = rules.decompress_annotation.output.anno output: tmp = temp("DE/{combo}/Featurecounts/{file}_tmp.counts"), tmph = temp("DE/{combo}/Featurecounts/{file}_tmp.head.gz"), tmpc = temp("DE/{combo}/Featurecounts/{file}_tmp.count.gz"), @@ -42,30 +52,30 @@ if COUNTTOOL == 'featureCounts': container: "oras://jfallmann/monsda:"+COUNTENV+"" threads: MAXTHREAD params: countb = COUNTBIN, - anno = ANNOTATION, cpara = lambda wildcards: tool_params(wildcards.file, None, config, "DE", DEENV.split('_')[0])['OPTIONS'].get('COUNT', ""), paired = lambda x: '-p' if paired == 'paired' else '', stranded = lambda x: '-s 1' if stranded == 'fr' else '-s 2' if stranded == 'rf' else '', sortmem = get_sortmem, prefix = lambda wildcards: "DE/"+combo+"/Featurecounts/"+wildcards.file - shell: "mkdir -p $(dirname {params.prefix}); zcat -f {params.anno} > {params.prefix}_annotation.gtf 2> {log}; {params.countb} -T {threads} {params.cpara} {params.paired} {params.stranded} -a {params.prefix}_annotation.gtf -o {output.tmp} {input.reads} 2>> {log} && rm -f {params.prefix}_annotation.gtf && head -n2 {output.tmp} |gzip > {output.tmph} && export LC_ALL=C; tail -n+3 {output.tmp}|sort --parallel={threads} -S {params.sortmem}G -T TMP -k1,1 -k2,2n -k3,3n -u |gzip >> {output.tmpc} && zcat {output.tmph} {output.tmpc} |gzip > {output.cts} && mv {output.tmp}.summary {output.cts}.summary" + shell: "mkdir -p $(dirname {params.prefix}); {params.countb} -T {threads} {params.cpara} {params.paired} {params.stranded} -a {input.anno} -o {output.tmp} {input.reads} 2> {log} && head -n2 {output.tmp} |gzip > {output.tmph} && export LC_ALL=C; tail -n+3 {output.tmp}|sort --parallel={threads} -S {params.sortmem}G -T TMP -k1,1 -k2,2n -k3,3n -u |gzip >> {output.tmpc} && zcat {output.tmph} {output.tmpc} |gzip > {output.cts} && mv {output.tmp}.summary {output.cts}.summary" else: # Long-read quantifiers (e.g. oarfish) quantify directly from the # genome-mapped BAM and are aggregated to gene level afterwards. rule featurecount_unique: - input: reads = expand("MAPPED/{scombo}/{{file}}_mapped_sorted.bam", scombo=scombo) + input: reads = expand("MAPPED/{scombo}/{{file}}_mapped_sorted.bam", scombo=scombo), + anno = rules.decompress_annotation.output.anno output: cts = "DE/{combo}/Featurecounts/{file}_mapped_sorted_unique.counts.gz" if not usededup else "DE/{combo}/Featurecounts/{file}_mapped_sorted_unique_dedup.counts.gz" log: "LOGS/{combo}/{file}/DE/deseq2/"+COUNTENV+"_deseq2_unique.log" conda: ""+COUNTENV+".yaml" container: "oras://jfallmann/monsda:"+COUNTENV+"" threads: MAXTHREAD params: countb = COUNTBIN, - anno = ANNOTATION, + origanno = ANNOTATION, ref = REFERENCE, cpara = lambda wildcards: tool_params(wildcards.file, None, config, "DE", DEENV.split('_')[0])['OPTIONS'].get('QUANT', ""), bins = BINS, prefix = lambda wildcards: "DE/"+combo+"/Featurecounts/"+wildcards.file - shell: "mkdir -p $(dirname {params.prefix}); samtools collate -@ {threads} -o {params.prefix}_namecollated.bam {input.reads} 2> {log}; zcat -f {params.anno} > {params.prefix}_annotation.gtf 2>> {log}; {params.countb} --threads {threads} --genome-alignments {params.prefix}_namecollated.bam --annotation {params.prefix}_annotation.gtf --genome-fasta {params.ref} {params.cpara} --output {params.prefix} 2>> {log} && rm -f {params.prefix}_namecollated.bam {params.prefix}_annotation.gtf && python3 {params.bins}/Analysis/oarfish_to_featurecounts.py -i {params.prefix}.quant -g {params.anno} -s {wildcards.file}_mapped_sorted.bam -o {output.cts} 2>> {log}" + shell: "mkdir -p $(dirname {params.prefix}); samtools collate -@ {threads} -o {params.prefix}_namecollated.bam {input.reads} 2> {log}; {params.countb} --threads {threads} --genome-alignments {params.prefix}_namecollated.bam --annotation {input.anno} --genome-fasta {params.ref} {params.cpara} --output {params.prefix} 2>> {log} && rm -f {params.prefix}_namecollated.bam && python3 {params.bins}/Analysis/oarfish_to_featurecounts.py -i {params.prefix}.quant -g {params.origanno} -s {wildcards.file}_mapped_sorted.bam -o {output.cts} 2>> {log}" rule prepare_count_table: input: cnd = expand(rules.featurecount_unique.output.cts, combo=combo, file=samplecond(SAMPLES, config)) diff --git a/workflows/dexseq_DTU.nf b/workflows/dexseq_DTU.nf index 33f0e7d0..7e268524 100644 --- a/workflows/dexseq_DTU.nf +++ b/workflows/dexseq_DTU.nf @@ -27,7 +27,7 @@ PCOMBO = get_always('COMBO') ?: 'none' //DTU PROCESSES -include { salmon_quant; oarfish_dtu_quant; create_summary_snippet; terminus_collapse } from "dtu_base.nf" +include { salmon_quant; decompress_annotation_dtu; oarfish_dtu_quant; create_summary_snippet; terminus_collapse } from "dtu_base.nf" process salmon_idx{ conda "$COUNTENV"+".yaml" @@ -182,7 +182,9 @@ workflow DTU{ element -> return "${workflow.workDir}/../MAPPED/${COMBO}/"+element+"_mapped_sorted.bam" } mapsamples_ch = Channel.fromPath(MAPPEDSAMPLES.sort()) - oarfish_dtu_quant(mapsamples_ch.collate(1)) + decompress_annotation_dtu(annofile) + annogtf_dtu = decompress_annotation_dtu.out.anno_plain + oarfish_dtu_quant(annogtf_dtu.combine(mapsamples_ch.collate(1))) quantcounts = oarfish_dtu_quant.out.counts } diff --git a/workflows/drimseq_DTU.nf b/workflows/drimseq_DTU.nf index 89cabc59..bad83c5e 100644 --- a/workflows/drimseq_DTU.nf +++ b/workflows/drimseq_DTU.nf @@ -27,7 +27,7 @@ PCOMBO = get_always('COMBO') ?: 'none' //DTU PROCESSES -include { salmon_quant; oarfish_dtu_quant; create_summary_snippet; terminus_collapse } from "dtu_base.nf" +include { salmon_quant; decompress_annotation_dtu; oarfish_dtu_quant; create_summary_snippet; terminus_collapse } from "dtu_base.nf" process salmon_idx{ conda "$COUNTENV"+".yaml" @@ -183,7 +183,9 @@ workflow DTU{ element -> return "${workflow.workDir}/../MAPPED/${COMBO}/"+element+"_mapped_sorted.bam" } mapsamples_ch = Channel.fromPath(MAPPEDSAMPLES.sort()) - oarfish_dtu_quant(mapsamples_ch.collate(1)) + decompress_annotation_dtu(annofile) + annogtf_dtu = decompress_annotation_dtu.out.anno_plain + oarfish_dtu_quant(annogtf_dtu.combine(mapsamples_ch.collate(1))) quantcounts = oarfish_dtu_quant.out.counts } diff --git a/workflows/dtu_base.nf b/workflows/dtu_base.nf index 0948f95e..8ce085b7 100644 --- a/workflows/dtu_base.nf +++ b/workflows/dtu_base.nf @@ -62,6 +62,30 @@ process salmon_quant{ } } +process decompress_annotation_dtu{ + conda "$COUNTENV"+".yaml" + container "oras://jfallmann/monsda:"+"$COUNTENV" + cpus 1 + cache 'lenient' + + publishDir "${workflow.workDir}/../" , mode: 'link', + saveAs: {filename -> + if (filename.indexOf("log") > 0) "LOGS/${SCOMBO}/DTU/"+"${COUNTENV}"+"/decompress_annotation.log" + } + + input: + path anno + + output: + path "annotation.gtf", emit: anno_plain + path "log", emit: log + + script: + """ + zcat -f $anno > annotation.gtf 2> log; if [ ! -s annotation.gtf ]; then echo "ERROR: decompressed annotation from $anno is empty, check source file and filesystem" >> log; exit 1; fi + """ +} + process oarfish_dtu_quant{ conda "$COUNTENV"+".yaml" container "oras://jfallmann/monsda:"+"$COUNTENV" @@ -75,20 +99,22 @@ process oarfish_dtu_quant{ } input: - path bam + path fls output: path "*.gz", includeInputs:false, emit: counts path "*.log", emit: logs script: + anno = fls[0] + bam = fls[1] fn = file(bam).getSimpleName().replaceAll(/\Q_mapped_sorted\E/,"") lf = "${COUNTENV}_"+fn+".log" prefix = fn+"/"+COUNTENV oz = fn+"/quant.sf.gz" ol = fn+"_counts.gz" """ - mkdir -p $fn; samtools collate -@ ${task.cpus} -o namecollated.bam $bam &>> $lf; zcat -f $DTUANNO > $fn/annotation.gtf 2>> $lf; $COUNTBIN --threads ${task.cpus} --genome-alignments namecollated.bam --annotation $fn/annotation.gtf --genome-fasta $DTUREF $COUNTPARAMS --output $prefix &>> $lf && rm -f namecollated.bam $fn/annotation.gtf && python3 $BINS/Analysis/oarfish_to_salmon.py -i ${prefix}.quant -o $oz &>> $lf && ln -fs $oz $ol + mkdir -p $fn; : > $lf; if [ ! -s $anno ]; then echo "ERROR: decompressed annotation $anno is empty, check source file and filesystem" >> $lf; exit 1; fi; samtools collate -@ ${task.cpus} -o namecollated.bam $bam &>> $lf; $COUNTBIN --threads ${task.cpus} --genome-alignments namecollated.bam --annotation $anno --genome-fasta $DTUREF $COUNTPARAMS --output $prefix &>> $lf && rm -f namecollated.bam && python3 $BINS/Analysis/oarfish_to_salmon.py -i ${prefix}.quant -o $oz &>> $lf && ln -fs $oz $ol """ } diff --git a/workflows/dtu_base.smk b/workflows/dtu_base.smk index f8daf4af..61142ff2 100644 --- a/workflows/dtu_base.smk +++ b/workflows/dtu_base.smk @@ -73,9 +73,18 @@ else: # is needed, output is normalized to the salmon quant.sf schema so # downstream DTU tooling (tximport type="salmon", build_DTU_table.py, # SPIT.py) does not need to change. + rule decompress_annotation: + input: anno = ANNOTATION + output: anno = temp("DTU/{combo}/"+COUNTENV+"/annotation.gtf") + log: "LOGS/{combo}/DTU/"+COUNTENV+"/decompress_annotation.log" + conda: ""+COUNTENV+".yaml" + container: "oras://jfallmann/monsda:"+COUNTENV+"" + threads: 1 + shell: "zcat -f {input.anno} > {output.anno} 2> {log}; if [ ! -s {output.anno} ]; then echo 'ERROR: decompressed annotation from {input.anno} is empty, check source file and filesystem' >> {log}; exit 1; fi" + rule mapping: input: bam = expand("MAPPED/{scombo}/{{file}}_mapped_sorted.bam", scombo=scombo), - anno = ANNOTATION, + anno = rules.decompress_annotation.output.anno, fa = REFERENCE output: cnts = report("DTU/{combo}/"+COUNTENV+"/{file}_counts.sf.gz", category="COUNTING"), ctsdir = report(directory("DTU/{combo}/"+COUNTENV+"/{file}"), category="COUNTING") @@ -86,11 +95,10 @@ else: params: cpara = lambda wildcards: tool_params(wildcards.file, None, config, 'DTU', DTUENV)['OPTIONS'].get('QUANT', ""), mapp = COUNTBIN, nbam = lambda wildcards, output: os.path.join(str(output.ctsdir), "namecollated.bam"), - anno = lambda wildcards, output: os.path.join(str(output.ctsdir), "annotation.gtf"), prefix = lambda wildcards, output: os.path.join(str(output.ctsdir), COUNTENV), linksf = lambda wildcards, output: str(os.path.abspath(output.ctsdir)), bins = BINS - shell: "set +euo pipefail; mkdir -p {output.ctsdir}; samtools collate -@ {threads} -o {params.nbam} {input.bam} &>> {log}; zcat -f {input.anno} > {params.anno} 2>> {log}; {params.mapp} --threads {threads} --genome-alignments {params.nbam} --annotation {params.anno} --genome-fasta {input.fa} {params.cpara} --output {params.prefix} &>> {log} && rm -f {params.nbam} {params.anno} && python3 {params.bins}/Analysis/oarfish_to_salmon.py -i {params.prefix}.quant -o {output.ctsdir}/quant.sf.gz &>> {log} ; ln -fs {params.linksf}/quant.sf.gz {output.cnts} &>> {log}" + shell: "set +euo pipefail; mkdir -p {output.ctsdir}; samtools collate -@ {threads} -o {params.nbam} {input.bam} &>> {log}; {params.mapp} --threads {threads} --genome-alignments {params.nbam} --annotation {input.anno} --genome-fasta {input.fa} {params.cpara} --output {params.prefix} &>> {log} && rm -f {params.nbam} && python3 {params.bins}/Analysis/oarfish_to_salmon.py -i {params.prefix}.quant -o {output.ctsdir}/quant.sf.gz &>> {log} ; ln -fs {params.linksf}/quant.sf.gz {output.cnts} &>> {log}" if runterminus: diff --git a/workflows/edger_DE.nf b/workflows/edger_DE.nf index 9b8f4b37..5966e26b 100644 --- a/workflows/edger_DE.nf +++ b/workflows/edger_DE.nf @@ -24,6 +24,30 @@ COUNTENV = (COUNTTOOL == 'featureCounts') ? 'countreads_de' : COUNTTOOL //DE PROCESSES +process decompress_annotation_edger{ + conda "$COUNTENV"+".yaml" + container "oras://jfallmann/monsda:"+"$COUNTENV" + cpus 1 + cache 'lenient' + + publishDir "${workflow.workDir}/../" , mode: 'link', + saveAs: {filename -> + if (filename.indexOf("log") > 0) "LOGS/${SCOMBO}/DE/edger/decompress_annotation.log" + } + + input: + path anno + + output: + path "annotation.gtf", emit: anno_plain + path "log", emit: log + + script: + """ + zcat -f $anno > annotation.gtf 2> log; if [ ! -s annotation.gtf ]; then echo "ERROR: decompressed annotation from $anno is empty, check source file and filesystem" >> log; exit 1; fi + """ +} + process featurecount_edger{ conda "$COUNTENV"+".yaml" container "oras://jfallmann/monsda:"+"$COUNTENV" @@ -69,7 +93,7 @@ process featurecount_edger{ stranded = '' } """ - mkdir -p TMP; zcat -f $anno > annotation.gtf 2> $ol; $COUNTBIN -T ${task.cpus} $COUNTPARAMS $pair $stranded -a annotation.gtf -o tmpcts $reads 2>> $ol && rm -f annotation.gtf && head -n2 tmpcts |gzip > $oc && export LC_ALL=C; tail -n+3 tmpcts|sort --parallel=${task.cpus} -S $sortmem -T TMP -k1,1 -k2,2n -k3,3n -u |gzip >> $oc 2>> $ol && mv tmpcts.summary $os + mkdir -p TMP; : > $ol; if [ ! -s $anno ]; then echo "ERROR: decompressed annotation $anno is empty, check source file and filesystem" >> $ol; exit 1; fi; $COUNTBIN -T ${task.cpus} $COUNTPARAMS $pair $stranded -a $anno -o tmpcts $reads 2>> $ol && head -n2 tmpcts |gzip > $oc && export LC_ALL=C; tail -n+3 tmpcts|sort --parallel=${task.cpus} -S $sortmem -T TMP -k1,1 -k2,2n -k3,3n -u |gzip >> $oc 2>> $ol && mv tmpcts.summary $os """ } @@ -103,7 +127,7 @@ process oarfish_gene_count{ os = fn+".counts.summary" ol = fn+".log" """ - samtools collate -@ ${task.cpus} -o namecollated.bam $bam &>> $ol; zcat -f $anno > annotation.gtf 2>> $ol; $COUNTBIN --threads ${task.cpus} --genome-alignments namecollated.bam --annotation annotation.gtf --genome-fasta $DEREF $QUANTPARAMS --output $fn &>> $ol && rm -f namecollated.bam annotation.gtf && python3 $BINS/Analysis/oarfish_to_featurecounts.py -i ${fn}.quant -g $anno -s ${fn}_mapped_sorted.bam -o $oc 2>> $ol && touch $os + : > $ol; if [ ! -s $anno ]; then echo "ERROR: decompressed annotation $anno is empty, check source file and filesystem" >> $ol; exit 1; fi; samtools collate -@ ${task.cpus} -o namecollated.bam $bam &>> $ol; $COUNTBIN --threads ${task.cpus} --genome-alignments namecollated.bam --annotation $anno --genome-fasta $DEREF $QUANTPARAMS --output $fn &>> $ol && rm -f namecollated.bam && python3 $BINS/Analysis/oarfish_to_featurecounts.py -i ${fn}.quant -g $anno -s ${fn}_mapped_sorted.bam -o $oc 2>> $ol && touch $os """ } @@ -283,6 +307,9 @@ workflow DE{ annofile = Channel.fromPath(DEANNO) //annofile.subscribe { println "ANNO: $it \t COMBO: ${COMBO} SCOMBO: ${SCOMBO} LONG: ${LONGSAMPLES}" } + decompress_annotation_edger(annofile) + annogtf = decompress_annotation_edger.out.anno_plain + if (COUNTTOOL == 'featureCounts'){ if (RUNDEDUP){ MAPPEDSAMPLES = LONGSAMPLES.collect{ @@ -297,7 +324,7 @@ workflow DE{ mapsamples_ch = Channel.fromPath(MAPPEDSAMPLES.sort()) //mapsamples_ch.subscribe { println "MAP: $it \t COMBO: ${COMBO} SCOMBO: ${SCOMBO} LONG: ${LONGSAMPLES}" } - featurecount_edger(annofile.combine(mapsamples_ch.collate(1))) + featurecount_edger(annogtf.combine(mapsamples_ch.collate(1))) countcts = featurecount_edger.out.fc_cts }else{ MAPPEDSAMPLES = LONGSAMPLES.collect{ @@ -305,7 +332,7 @@ workflow DE{ } mapsamples_ch = Channel.fromPath(MAPPEDSAMPLES.sort()) - oarfish_gene_count(annofile.combine(mapsamples_ch.collate(1))) + oarfish_gene_count(annogtf.combine(mapsamples_ch.collate(1))) countcts = oarfish_gene_count.out.fc_cts } diff --git a/workflows/edger_DE.smk b/workflows/edger_DE.smk index 94c02d26..5fdcda3f 100644 --- a/workflows/edger_DE.smk +++ b/workflows/edger_DE.smk @@ -29,9 +29,19 @@ rule themall: Rmd = expand("REPORTS/SUMMARY/RmdSnippets/{combo}.Rmd", combo=combo), postde = expand("POSTDE/{combo}", combo=combo) if postde_enabled else [] +rule decompress_annotation: + input: anno = ANNOTATION + output: anno = temp("DE/{combo}/Featurecounts/annotation.gtf") + log: "LOGS/{combo}/DE/edger/decompress_annotation.log" + conda: ""+COUNTENV+".yaml" + container: "oras://jfallmann/monsda:"+COUNTENV+"" + threads: 1 + shell: "zcat -f {input.anno} > {output.anno} 2> {log}; if [ ! -s {output.anno} ]; then echo 'ERROR: decompressed annotation from {input.anno} is empty, check source file and filesystem' >> {log}; exit 1; fi" + if COUNTTOOL == 'featureCounts': rule featurecount_unique: - input: reads = expand("MAPPED/{scombo}/{{file}}_mapped_sorted_unique.bam", scombo=scombo) if not usededup else expand("MAPPED/{scombo}/{{file}}_mapped_sorted_unique_dedup.bam", scombo=scombo) + input: reads = expand("MAPPED/{scombo}/{{file}}_mapped_sorted_unique.bam", scombo=scombo) if not usededup else expand("MAPPED/{scombo}/{{file}}_mapped_sorted_unique_dedup.bam", scombo=scombo), + anno = rules.decompress_annotation.output.anno output: tmp = temp("DE/{combo}/Featurecounts/{file}_tmp.counts"), tmph = temp("DE/{combo}/Featurecounts/{file}_tmp.head.gz"), tmpc = temp("DE/{combo}/Featurecounts/{file}_tmp.count.gz"), @@ -41,30 +51,30 @@ if COUNTTOOL == 'featureCounts': container: "oras://jfallmann/monsda:"+COUNTENV+"" threads: MAXTHREAD params: countb = COUNTBIN, - anno = ANNOTATION, cpara = lambda wildcards: tool_params(wildcards.file, None, config, "DE", DEENV.split('_')[0])['OPTIONS'].get('COUNT', ""), paired = lambda x: '-p' if paired == 'paired' else '', stranded = lambda x: '-s 1' if stranded == 'fr' else '-s 2' if stranded == 'rf' else '', sortmem = get_sortmem, prefix = lambda wildcards: "DE/"+combo+"/Featurecounts/"+wildcards.file - shell: "mkdir -p $(dirname {params.prefix}); zcat -f {params.anno} > {params.prefix}_annotation.gtf 2> {log}; {params.countb} -T {threads} {params.cpara} {params.paired} {params.stranded} -a {params.prefix}_annotation.gtf -o {output.tmp} {input.reads} 2>> {log} && rm -f {params.prefix}_annotation.gtf && head -n2 {output.tmp} |gzip > {output.tmph} && export LC_ALL=C; tail -n+3 {output.tmp}|sort --parallel={threads} -S {params.sortmem}G -T TMP -k1,1 -k2,2n -k3,3n -u |gzip >> {output.tmpc} && zcat {output.tmph} {output.tmpc} |gzip > {output.cts} && mv {output.tmp}.summary {output.cts}.summary" + shell: "mkdir -p $(dirname {params.prefix}); {params.countb} -T {threads} {params.cpara} {params.paired} {params.stranded} -a {input.anno} -o {output.tmp} {input.reads} 2> {log} && head -n2 {output.tmp} |gzip > {output.tmph} && export LC_ALL=C; tail -n+3 {output.tmp}|sort --parallel={threads} -S {params.sortmem}G -T TMP -k1,1 -k2,2n -k3,3n -u |gzip >> {output.tmpc} && zcat {output.tmph} {output.tmpc} |gzip > {output.cts} && mv {output.tmp}.summary {output.cts}.summary" else: # Long-read quantifiers (e.g. oarfish) quantify directly from the # genome-mapped BAM and are aggregated to gene level afterwards. rule featurecount_unique: - input: reads = expand("MAPPED/{scombo}/{{file}}_mapped_sorted.bam", scombo=scombo) + input: reads = expand("MAPPED/{scombo}/{{file}}_mapped_sorted.bam", scombo=scombo), + anno = rules.decompress_annotation.output.anno output: cts = "DE/{combo}/Featurecounts/{file}_mapped_sorted_unique.counts.gz" if not usededup else "DE/{combo}/Featurecounts/{file}_mapped_sorted_unique_dedup.counts.gz" log: "LOGS/{combo}/{file}/DE/edger/"+COUNTENV+"_edger_unique.log" conda: ""+COUNTENV+".yaml" container: "oras://jfallmann/monsda:"+COUNTENV+"" threads: MAXTHREAD params: countb = COUNTBIN, - anno = ANNOTATION, + origanno = ANNOTATION, ref = REFERENCE, cpara = lambda wildcards: tool_params(wildcards.file, None, config, "DE", DEENV.split('_')[0])['OPTIONS'].get('QUANT', ""), bins = BINS, prefix = lambda wildcards: "DE/"+combo+"/Featurecounts/"+wildcards.file - shell: "mkdir -p $(dirname {params.prefix}); samtools collate -@ {threads} -o {params.prefix}_namecollated.bam {input.reads} 2> {log}; zcat -f {params.anno} > {params.prefix}_annotation.gtf 2>> {log}; {params.countb} --threads {threads} --genome-alignments {params.prefix}_namecollated.bam --annotation {params.prefix}_annotation.gtf --genome-fasta {params.ref} {params.cpara} --output {params.prefix} 2>> {log} && rm -f {params.prefix}_namecollated.bam {params.prefix}_annotation.gtf && python3 {params.bins}/Analysis/oarfish_to_featurecounts.py -i {params.prefix}.quant -g {params.anno} -s {wildcards.file}_mapped_sorted.bam -o {output.cts} 2>> {log}" + shell: "mkdir -p $(dirname {params.prefix}); samtools collate -@ {threads} -o {params.prefix}_namecollated.bam {input.reads} 2> {log}; {params.countb} --threads {threads} --genome-alignments {params.prefix}_namecollated.bam --annotation {input.anno} --genome-fasta {params.ref} {params.cpara} --output {params.prefix} 2>> {log} && rm -f {params.prefix}_namecollated.bam && python3 {params.bins}/Analysis/oarfish_to_featurecounts.py -i {params.prefix}.quant -g {params.origanno} -s {wildcards.file}_mapped_sorted.bam -o {output.cts} 2>> {log}" rule prepare_count_table: input: cnd = expand(rules.featurecount_unique.output.cts, combo=combo, file=samplecond(SAMPLES, config)) diff --git a/workflows/header.smk b/workflows/header.smk index 29496a2b..3fec00b6 100644 --- a/workflows/header.smk +++ b/workflows/header.smk @@ -326,6 +326,11 @@ if 'FUSIONS' in config: ANNOTATION = ANNO.get('GTF') if 'GTF' in ANNO and ANNO.get('GTF') != '' else ANNO.get('GFF') # by default GTF format will be used +if ANNOTATION and isinstance(ANNOTATION, str): + # Resolve once to an absolute path so every rule's shell command can read + # it reliably regardless of its own working directory (e.g. containers). + ANNOTATION = os.path.abspath(ANNOTATION) + combo = '' if ANNOTATION: diff --git a/workflows/oarfish.nf b/workflows/oarfish.nf index 05857a35..0fb889c9 100644 --- a/workflows/oarfish.nf +++ b/workflows/oarfish.nf @@ -10,6 +10,30 @@ COUNTPARAMS = get_always('oarfish_params_COUNT') ?: '' //COUNTING PROCESSES +process decompress_annotation_oarfish{ + conda "$COUNTENV"+".yaml" + container "oras://jfallmann/monsda:"+"$COUNTENV" + cpus 1 + cache 'lenient' + + publishDir "${workflow.workDir}/../" , mode: 'link', + saveAs: {filename -> + if (filename.indexOf("log") > 0) "LOGS/${COMBO}/COUNTING/oarfish/decompress_annotation.log" + } + + input: + path anno + + output: + path "annotation.gtf", emit: anno_plain + path "log", emit: log + + script: + """ + zcat -f $anno > annotation.gtf 2> log; if [ ! -s annotation.gtf ]; then echo "ERROR: decompressed annotation from $anno is empty, check source file and filesystem" >> log; exit 1; fi + """ +} + process oarfish_quant{ conda "$COUNTENV"+".yaml" container "oras://jfallmann/monsda:"+"$COUNTENV" @@ -23,7 +47,7 @@ process oarfish_quant{ } input: - path bam + path fls output: path "*.gz", includeInputs:false, emit: counts @@ -31,12 +55,14 @@ process oarfish_quant{ path "*", includeInputs:false, emit: rest script: + anno = fls[0] + bam = fls[1] fn = file(bam).getSimpleName().replaceAll(/\Q_mapped_sorted\E/,"") lf = "oarfish_"+fn+".log" oz = fn+"/quant.sf.gz" ol = fn+"_counts.gz" """ - mkdir -p $fn; samtools collate -@ ${task.cpus} -o namecollated.bam $bam &>> $lf; zcat -f $COUNTANNO > $fn/annotation.gtf 2>> $lf; $COUNTBIN --threads ${task.cpus} --genome-alignments namecollated.bam --annotation $fn/annotation.gtf --genome-fasta $COUNTREF $COUNTPARAMS --output $fn/oarfish &>> $lf && rm -f namecollated.bam $fn/annotation.gtf && python3 $BINS/Analysis/oarfish_to_salmon.py -i $fn/oarfish.quant -o $oz &>> $lf && ln -fs $oz $ol + mkdir -p $fn; : > $lf; if [ ! -s $anno ]; then echo "ERROR: decompressed annotation $anno is empty, check source file and filesystem" >> $lf; exit 1; fi; samtools collate -@ ${task.cpus} -o namecollated.bam $bam &>> $lf; $COUNTBIN --threads ${task.cpus} --genome-alignments namecollated.bam --annotation $anno --genome-fasta $COUNTREF $COUNTPARAMS --output $fn/oarfish &>> $lf && rm -f namecollated.bam && python3 $BINS/Analysis/oarfish_to_salmon.py -i $fn/oarfish.quant -o $oz &>> $lf && ln -fs $oz $ol """ } @@ -50,8 +76,11 @@ workflow COUNTING{ } mapsamples_ch = Channel.fromPath(MAPPEDSAMPLES.sort()) + annofile = Channel.fromPath(COUNTANNO) + decompress_annotation_oarfish(annofile) + annogtf = decompress_annotation_oarfish.out.anno_plain - oarfish_quant(mapsamples_ch.collate(1)) + oarfish_quant(annogtf.combine(mapsamples_ch.collate(1))) emit: counts = oarfish_quant.out.counts diff --git a/workflows/oarfish.smk b/workflows/oarfish.smk index 52396b44..d7818bb6 100644 --- a/workflows/oarfish.smk +++ b/workflows/oarfish.smk @@ -8,9 +8,18 @@ unik = get_dict_hash(keydict) rule themall: input: expand("COUNTS/{combo}/{file}_counts.sf.gz", combo=combo, file=samplecond(SAMPLES, config)) +rule decompress_annotation: + input: anno = ANNOTATION + output: anno = temp("COUNTS/{combo}/annotation.gtf") + log: "LOGS/{combo}/COUNTING/oarfish/decompress_annotation.log" + conda: ""+COUNTENV+".yaml" + container: "oras://jfallmann/monsda:"+COUNTENV+"" + threads: 1 + shell: "zcat -f {input.anno} > {output.anno} 2> {log}; if [ ! -s {output.anno} ]; then echo 'ERROR: decompressed annotation from {input.anno} is empty, check source file and filesystem' >> {log}; exit 1; fi" + rule mapping: input: bam = expand("MAPPED/{scombo}/{{file}}_mapped_sorted.bam", scombo=scombo), - anno = ANNOTATION, + anno = rules.decompress_annotation.output.anno, fa = REFERENCE output: cnts = report("COUNTS/{combo}/{file}_counts.sf.gz", category="COUNTING"), ctsdir = report(directory("COUNTS/{combo}/{file}"), category="COUNTING") @@ -21,8 +30,7 @@ rule mapping: params: cpara = lambda wildcards: tool_params(wildcards.file, None, config, 'COUNTING', COUNTENV)['OPTIONS'].get('COUNT', ""), mapp = COUNTBIN, nbam = lambda wildcards, output: os.path.join(str(output.ctsdir), "namecollated.bam"), - anno = lambda wildcards, output: os.path.join(str(output.ctsdir), "annotation.gtf"), prefix = lambda wildcards, output: os.path.join(str(output.ctsdir), "oarfish"), linksf = lambda wildcards, output: str(os.path.abspath(output.ctsdir)), bins = BINS - shell: "set +euo pipefail; mkdir -p {output.ctsdir}; samtools collate -@ {threads} -o {params.nbam} {input.bam} &>> {log}; zcat -f {input.anno} > {params.anno} 2>> {log}; {params.mapp} --threads {threads} --genome-alignments {params.nbam} --annotation {params.anno} --genome-fasta {input.fa} {params.cpara} --output {params.prefix} &>> {log} && rm -f {params.nbam} {params.anno} && python3 {params.bins}/Analysis/oarfish_to_salmon.py -i {params.prefix}.quant -o {output.ctsdir}/quant.sf.gz &>> {log} ; ln -fs {params.linksf}/quant.sf.gz {output.cnts} &>> {log}" + shell: "set +euo pipefail; mkdir -p {output.ctsdir}; samtools collate -@ {threads} -o {params.nbam} {input.bam} &>> {log}; {params.mapp} --threads {threads} --genome-alignments {params.nbam} --annotation {input.anno} --genome-fasta {input.fa} {params.cpara} --output {params.prefix} &>> {log} && rm -f {params.nbam} && python3 {params.bins}/Analysis/oarfish_to_salmon.py -i {params.prefix}.quant -o {output.ctsdir}/quant.sf.gz &>> {log} ; ln -fs {params.linksf}/quant.sf.gz {output.cnts} &>> {log}" diff --git a/workflows/spit_DTU.nf b/workflows/spit_DTU.nf index 4800f025..2ae920e5 100644 --- a/workflows/spit_DTU.nf +++ b/workflows/spit_DTU.nf @@ -27,7 +27,7 @@ PCOMBO = get_always('COMBO') ?: 'none' //DTU PROCESSES -include { salmon_quant; oarfish_dtu_quant; create_summary_snippet; terminus_collapse } from "dtu_base.nf" +include { salmon_quant; decompress_annotation_dtu; oarfish_dtu_quant; create_summary_snippet; terminus_collapse } from "dtu_base.nf" process salmon_idx{ conda "$COUNTENV"+".yaml" @@ -181,7 +181,9 @@ workflow DTU{ element -> return "${workflow.workDir}/../MAPPED/${COMBO}/"+element+"_mapped_sorted.bam" } mapsamples_ch = Channel.fromPath(MAPPEDSAMPLES.sort()) - oarfish_dtu_quant(mapsamples_ch.collate(1)) + decompress_annotation_dtu(annofile) + annogtf_dtu = decompress_annotation_dtu.out.anno_plain + oarfish_dtu_quant(annogtf_dtu.combine(mapsamples_ch.collate(1))) quantcounts = oarfish_dtu_quant.out.counts }