diff --git a/aviary/modules/assembly/assembly.smk b/aviary/modules/assembly/assembly.smk index fa06fea4..07bb5caa 100644 --- a/aviary/modules/assembly/assembly.smk +++ b/aviary/modules/assembly/assembly.smk @@ -69,6 +69,28 @@ rule flye_assembly: script: "scripts/run_flye.py" +rule metamdbg_assembly: + input: + fastq = "data/long_reads.fastq.gz" + output: + fasta = "data/metamdbg/contigs.fasta.gz", + log = "data/metamdbg/metaMDBG.log", # Some info here is surplus to snakemake log file. + junk1 = temp(directory("data/metamdbg/tmp/")), + params: + long_read_type = config["long_read_type"] + threads: + config["max_threads"] + resources: # TODO: Currently these are copied from flye_assembly, but they should be adjusted + mem_mb = lambda wildcards, attempt: min(int(config["max_memory"])*1024, 512*1024*attempt), + runtime = lambda wildcards, attempt: 24*60 + 24*60*attempt, + log: + "logs/metamdbg_assembly.log" + conda: + "envs/metamdbg.yaml" + benchmark: + "benchmarks/metamdbg_assembly.benchmark.txt" + script: + "scripts/run_metamdbg.py" # Polish the long reads assembly with Racon or Medaka rule polish_metagenome_flye: diff --git a/aviary/modules/assembly/envs/metamdbg.yaml b/aviary/modules/assembly/envs/metamdbg.yaml new file mode 100644 index 00000000..fe951800 --- /dev/null +++ b/aviary/modules/assembly/envs/metamdbg.yaml @@ -0,0 +1,5 @@ +channels: + - conda-forge + - bioconda +dependencies: + - metamdbg = 1.0 \ No newline at end of file diff --git a/aviary/modules/assembly/scripts/run_metamdbg.py b/aviary/modules/assembly/scripts/run_metamdbg.py new file mode 100755 index 00000000..1417e609 --- /dev/null +++ b/aviary/modules/assembly/scripts/run_metamdbg.py @@ -0,0 +1,70 @@ +from subprocess import run, STDOUT +import os +from pathlib import Path + + +def dummy_fasta(output_fasta: str): + with open(output_fasta, "w") as out: + out.write(">dummy\n") + out.write("ACGT\n") + + +def create_dummy_output(output_dir: str): + # need to create dummy output files + os.makedirs(output_dir, exist_ok=True) + output_fasta = f"{output_dir}/assembly.fasta" + dummy_fasta(output_fasta) + # Path.joinpath(Path(output_dir), "assembly_info.txt").touch() + # Path.joinpath(Path(output_dir), "assembly_graph.gfa").touch() + # os.makedirs(f"{output_dir}/00-assembly", exist_ok=True) + # os.makedirs(f"{output_dir}/10-consensus", exist_ok=True) + # os.makedirs(f"{output_dir}/20-repeat", exist_ok=True) + # os.makedirs(f"{output_dir}/30-contigger", exist_ok=True) + # os.makedirs(f"{output_dir}/40-polishing", exist_ok=True) + + +def run_metamdbg( + long_read_type: str, + input_fastq: str, + output_dir: str, + threads: int, + log: str, +): + # check if input_fastq has any reads + if os.path.getsize(input_fastq) == 0: + with open(log, "a") as logf: + logf.write(f"Input fastq file {input_fastq} is empty\n") + logf.write(f"Skipping metamdbg assembly\n") + create_dummy_output(output_dir) + return + + if long_read_type == 'ont': + read_type = "--in-ont" + elif long_read_type == 'ont_hq': + read_type = "--in-ont" + elif long_read_type == 'ccs' or long_read_type == 'hifi': + read_type = "--in-hifi" + else: + raise ValueError(f"Invalid long read type: {long_read_type}") + + cmd = f"metaMDBG asm {read_type} {input_fastq} --out-dir {output_dir} --threads {threads}".split() + with open(log, "a") as logf: + run(cmd, stdout=logf, stderr=STDOUT) + + +if __name__ == '__main__': + long_read_type = snakemake.params.long_read_type + input_fastq = snakemake.input.fastq + output_dir = "data/metamdbg" + threads = snakemake.threads + log = snakemake.log[0] + + with open(log, "w") as logf: pass + + run_metamdbg( + long_read_type=long_read_type, + input_fastq=input_fastq, + output_dir=output_dir, + threads=threads, + log=log, + ) diff --git a/test/data/ont_badread_10.4.1_GCA_000503915.1_ASM50391v1_genomic.fastq.gz b/test/data/ont_badread_10.4.1_GCA_000503915.1_ASM50391v1_genomic.fastq.gz new file mode 100644 index 00000000..bf084e45 Binary files /dev/null and b/test/data/ont_badread_10.4.1_GCA_000503915.1_ASM50391v1_genomic.fastq.gz differ