I have the below bash script that I would like to convert into a snakefile:
mmseqs rbh flye_db megahit_db flye_megahit_rbh --min-seq-id 0.9 mmseq2_tmp --threads 12
mmseqs rbh flye_db metaspades_db flye_metaspades_rbh --min-seq-id 0.9 mmseq2_tmp --threads 12
mmseqs rbh megahit_db metaspades_db megahit_metaspades_rbh --min-seq-id 0.9 mmseq2_tmp --threads 12
I've managed to come up with the following, but was wondering if there's a way to use regex or expand to improve the code further:
rule mmseq2_compare:
input:
mm1=expand(os.path.join(RESULTS_DIR, "annotation/mmseq2/{assembler}_db"), assembler="flye"),
mm2=expand(os.path.join(RESULTS_DIR, "annotation/mmseq2/{assembler}_db"), assembler="megahit"),
mm3=expand(os.path.join(RESULTS_DIR, "annotation/mmseq2/{assembler}_db"), assembler="metaspades_hybrid")
output:
mo1=os.path.join(RESULTS_DIR, "annotation/mmseq2/flye_megahit_rbh"),
mo2=os.path.join(RESULTS_DIR, "annotation/mmseq2/flye_metaspades_hybrid_rbh"),
mo3=os.path.join(RESULTS_DIR, "annotation/mmseq2/megahit_metaspades_hybrid_rbh")
log: os.path.join(RESULTS_DIR, "annotation/mmseq2/compare.mmseq2.log")
conda: "cd-hit.yml"
shell:
"""
(date &&\
mmseqs rbh {input.mm1} {input.mm2} {output.mo1} --min-seq-id 0.9 mmseq2_tmp --threads 12 &&\
mmseqs rbh {input.mm1} {input.mm3} {output.mo2} --min-seq-id 0.9 mmseq2_tmp --threads 12 &&\
mmseqs rbh {input.mm2} {input.mm3} {output.mo3} --min-seq-id 0.9 mmseq2_tmp --threads 12 &&\
date) &> >(tee {log})
"""
With the 3 assemblers (flye, megahit and metaspades_hybrid) is there any way to remove the redundancy especially in the 'shell'?
Thank you!
Dry-run output
Building DAG of jobs...
Job counts:
count jobs
1 all
1 mmseq_compare
2
[Thu Mar 26 12:25:14 2020]
rule mmseq_compare:
input: results/annotation/mmseq2/flye_db, results/annotation/mmseq2/megahit_db
output: results/annotation/mmseq2/flye_megahit_rbh
jobid: 1
wildcards: assembler1=flye, assembler2=megahit
mmseqs rbh results/annotation/mmseq2/flye_db results/annotation/mmseq2/megahit_db results/annotation/mmseq2/flye_megahit_rbh --min-seq-id 0.9 mmseq2_tmp --threads 12
[Thu Mar 26 12:25:14 2020]
localrule all:
input: results/annotation/mmseq2/flye_megahit_rbh, results/annotation/mmseq2/flye_metaspades_hybrid_rbh, results/annotation/mmseq2/megahit_metaspades_hybrid_rbh
jobid: 0
Job counts:
count jobs
1 all
1 mmseq_compare
2
This was a dry-run (flag -n). The order of jobs does not reflect the order of execution.```
Firstly you don't need an expand
in your input. This is needed if you wish to create a list of filenames that have the same pattern.
Next, as long as you already use Unix-type slash in your paths, you may add the RESULTS_DIR
into f-strings for readability (but don't forget to double the braces for wildcards).
Finally, there is no need in having a pipeline of scripts separated with &&: that is what Snakemake is designed for.
My version of reworked script:
rule all:
input:
expand(f"{RESULTS_DIR}/annotation/mmseq2/{{assembler1}}__{{assembler2}}_rbh", zip,
assembler1=["flye", "flye", "megahit"],
assembler2=["megahit", "megahit_metaspades_hybrid", "megahit_metaspades_hybrid"])
rule mmseq_compare:
input:
f"{RESULTS_DIR}/annotation/mmseq2/{{assembler1}}_db",
f"{RESULTS_DIR}/annotation/mmseq2/{{assembler2}}_db"
output:
f"{RESULTS_DIR}/annotation/mmseq2/{{assembler1}}__{{assembler2}}_rbh"
conda:
"cd-hit.yml"
shell:
"mmseqs rbh {input[0]} {input[1]} {output} --min-seq-id 0.9 mmseq2_tmp --threads 12"
I've excluded the date
and logging. My solution has a limitaton that the order of execution of the comparisons is undefined: you need to reconsider the logging strategy in this case.