Repository: AroneyS/binchicken (language: python)

Issue:
Remove bbduk qc filtering
fastp filtering added to Aviary


Discussion / hints:

Fix patch:
diff --git a/ibis/ibis.py b/ibis/ibis.py
index 4a10036f..f10b1610 100755
--- a/ibis/ibis.py
+++ b/ibis/ibis.py
@@ -105,7 +105,7 @@ def run_workflow(config, workflow, output_dir, cores=16, dryrun=False,
         snakefile=os.path.join(os.path.dirname(__file__), "workflow", workflow),
         config=config,
         output_dir=output_dir,
-        jobs=f"--jobs {cores}" if cores is not None else "",
+        jobs=f"--cores {cores}" if cores is not None else "--cores 1",
         profile="" if not profile else f"--profile {profile}",
         retries="" if (cluster_retries is None) else f"--retries {cluster_retries}",
         conda_frontend=f"--conda-frontend {conda_frontend}" if conda_frontend is not None else "",
@@ -263,42 +263,16 @@ def download_sra(args):
         except FileNotFoundError:
             raise Exception("Single-ended reads detected. Remove from input SRA and try again.")
 
-    config_items = {
-        "sra": [s for s in args.forward if s not in single_ended],
-        "reads_1": {},
-        "reads_2": {},
-        "snakemake_profile": args.snakemake_profile,
-        "cluster_retries": args.cluster_retries,
-    }
-
-    config_path = make_config(
-        importlib.resources.files("ibis.config").joinpath("template_coassemble.yaml"),
-        args.output,
-        config_items
-        )
-
-    run_workflow(
-        config = config_path,
-        workflow = "coassemble.smk",
-        output_dir = args.output,
-        cores = args.cores,
-        dryrun = args.dryrun,
-        profile = args.snakemake_profile,
-        cluster_retries = args.cluster_retries,
-        conda_prefix = args.conda_prefix,
-        snakemake_args = args.snakemake_args + " -- compile_sra_qc" if args.snakemake_args else "-- compile_sra_qc",
-    )
-
     # Need to convince Snakemake that the SRA data predates the completed outputs
     # Otherwise it will rerun all the rules with reason "updated input files"
     # Using Jan 1 2000 as pre-date
-    sra_qc_dir = args.output + "/coassemble/sra_qc/"
-    os.makedirs(sra_qc_dir, exist_ok=True)
+    sra_dir = args.output + "/coassemble/sra/"
+    os.makedirs(sra_dir, exist_ok=True)
     for sra in [f + s for f in args.forward for s in ["_1", "_2"]]:
-        subprocess.check_call(f"touch -t 200001011200 {sra_qc_dir + sra + SRA_SUFFIX}", shell=True)
+        subprocess.check_call(f"touch -t 200001011200 {sra_dir + sra + SRA_SUFFIX}", shell=True)
 
-    forward = sorted([sra_qc_dir + f for f in os.listdir(sra_qc_dir) if f.endswith("_1" + SRA_SUFFIX)])
-    reverse = sorted([sra_qc_dir + f for f in os.listdir(sra_qc_dir) if f.endswith("_2" + SRA_SUFFIX)])
+    forward = sorted([sra_dir + f for f in os.listdir(sra_dir) if f.endswith("_1" + SRA_SUFFIX)])
+    reverse = sorted([sra_dir + f for f in os.listdir(sra_dir) if f.endswith("_2" + SRA_SUFFIX)])
 
     return forward, reverse
 
diff --git a/ibis/workflow/coassemble.smk b/ibis/workflow/coassemble.smk
index 19375d42..debbff9b 100644
--- a/ibis/workflow/coassemble.smk
+++ b/ibis/workflow/coassemble.smk
@@ -3,7 +3,7 @@
 #############
 ruleorder: no_genomes > query_processing > update_appraise > singlem_appraise
 ruleorder: mock_download_sra > download_sra
-localrules: all, summary, singlem_summarise_genomes, singlem_appraise_filtered, no_genomes, mock_download_sra, compile_sra_qc, collect_genomes, finish_mapping, aviary_commands, aviary_recover, aviary_combine
+localrules: all, summary, singlem_summarise_genomes, singlem_appraise_filtered, no_genomes, mock_download_sra, collect_genomes, finish_mapping, aviary_commands, aviary_recover, aviary_combine
 
 import os
 import pandas as pd
@@ -403,43 +403,6 @@ rule mock_download_sra:
         "mkdir -p {output} && "
         "cp {params.sra_u} {params.sra_f} {params.sra_r} {output}"
 
-rule sra_qc:
-    output:
-        out1 = output_dir + "/sra_qc/{sra}_1.fastq.gz",
-        out2 = output_dir + "/sra_qc/{sra}_2.fastq.gz",
-        stats = output_dir + "/sra_qc/{sra}_qc.stats",
-    params:
-        in1 = output_dir + "/sra/{sra}_1.fastq.gz",
-        in2 = output_dir + "/sra/{sra}_2.fastq.gz",
-    threads: 16
-    resources:
-        mem_mb=125*1000,
-        runtime = "24h",
-    conda:
-        "env/bbtools.yml"
-    log:
-        logs_dir + "/sra_qc/{sra}_qc.log"
-    shell:
-        "bbduk.sh "
-        "in1={params.in1} in2={params.in2} "
-        "out1={output.out1} out2={output.out2} "
-        "ref=adapters,phix "
-        "ktrim=r k=23 mink=11 hdist=1 "
-        "qtrim=r trimq=10 "
-        "minlen=30 "
-        "stats={output.stats} statscolumns=3 "
-        "t={threads} "
-        "&> {log} "
-
-rule compile_sra_qc:
-    input:
-        expand(output_dir + "/sra_qc/{sra}_1.fastq.gz", sra=config["sra"]),
-        expand(output_dir + "/sra_qc/{sra}_2.fastq.gz", sra=config["sra"]),
-    output:
-        done = output_dir + "/sra_qc/done"
-    shell:
-        "touch {output.done} "
-
 #####################################
 ### Map reads to matching genomes ###
 #####################################
@@ -463,10 +426,11 @@ rule map_reads:
         genomes = output_dir + "/mapping/{read}_reference.fna",
     output:
         dir = temp(directory(output_dir + "/mapping/{read}_coverm")),
+    group: "unmapping"
     threads: 16
     resources:
         mem_mb=125*1000,
-        runtime = "24h",
+        runtime = "12h",
     log:
         logs_dir + "/mapping/{read}_coverm.log",
     conda:
@@ -485,6 +449,7 @@ rule filter_bam_files:
         output_dir + "/mapping/{read}_coverm",
     output:
         temp(output_dir + "/mapping/{read}_unmapped.bam"),
+    group: "unmapping"
     params:
         genomes = "{read}_reference.fna",
         reads_1 = lambda wildcards: os.path.basename(config["reads_1"][wildcards.read]),
@@ -493,7 +458,7 @@ rule filter_bam_files:
     threads: 16
     resources:
         mem_mb=125*1000,
-        runtime = "24h",
+        runtime = "4h",
     log:
         logs_dir + "/mapping/{read}_filter.log",
     conda:
@@ -514,10 +479,11 @@ rule bam_to_fastq:
     output:
         reads_1 = output_dir + "/mapping/{read}_unmapped.1.fq.gz",
         reads_2 = output_dir + "/mapping/{read}_unmapped.2.fq.gz",
+    group: "unmapping"
     threads: 16
     resources:
         mem_mb=125*1000,
-        runtime = "24h",
+        runtime = "4h",
     log:
         logs_dir + "/mapping/{read}_fastq.log",
     conda:
diff --git a/ibis/workflow/env/bbtools.yml b/ibis/workflow/env/bbtools.yml
deleted file mode 100644
index 8e761172..00000000
--- a/ibis/workflow/env/bbtools.yml
+++ /dev/null
@@ -1,7 +0,0 @@
-channels:
-  - conda-forge
-  - bioconda
-  - AgBiome
-dependencies:
-  - python>=3.7
-  - bbtools


Test patch:
diff --git a/test/test_manual.py b/test/test_manual.py
index 9f317e29..ac610622 100644
--- a/test/test_manual.py
+++ b/test/test_manual.py
@@ -76,18 +76,6 @@ def test_update_sra_download_real(self):
         self.assertTrue(os.path.exists(os.path.join(output_dir, "coassemble", "sra", "SRR8334324_1.fastq.gz")))
         self.assertTrue(os.path.exists(os.path.join(output_dir, "coassemble", "sra", "SRR8334324_2.fastq.gz")))
 
-        sra_1_path = os.path.join(output_dir, "coassemble", "sra_qc", "SRR8334323_1.fastq.gz")
-        self.assertTrue(os.path.exists(sra_1_path))
-        with gzip.open(sra_1_path) as f:
-            file = f.readline().decode()
-            self.assertTrue("@SRR8334323.1 HS2:487:H80UEADXX:1:1101:1148:1986/1" in file)
-            self.assertTrue("@SRR8334323.11 HS2:487:H80UEADXX:1:1101:2569:1958/1" not in file)
-            self.assertTrue("@SRR8334323.15 HS2:487:H80UEADXX:1:1101:3302:1937/1" not in file)
-
-        self.assertTrue(os.path.exists(os.path.join(output_dir, "coassemble", "sra_qc", "SRR8334323_2.fastq.gz")))
-        self.assertTrue(os.path.exists(os.path.join(output_dir, "coassemble", "sra_qc", "SRR8334324_1.fastq.gz")))
-        self.assertTrue(os.path.exists(os.path.join(output_dir, "coassemble", "sra_qc", "SRR8334324_2.fastq.gz")))
-
     def test_update_aviary_run_real(self):
         output_dir = os.path.join("example", "test_update_aviary_run_real")
         self.setup_output_dir(output_dir)
diff --git a/test/test_update.py b/test/test_update.py
index 6a4aec0f..5b9406b6 100644
--- a/test/test_update.py
+++ b/test/test_update.py
@@ -220,17 +220,9 @@ def test_update_sra_download(self):
 
             output_sra = output_comb.split("Building DAG of jobs...")[1]
             self.assertTrue("download_sra" in output_sra)
-            self.assertTrue("sra_qc" not in output_sra)
-            self.assertTrue("compile_sra_qc" not in output_sra)
             self.assertTrue("aviary_commands" not in output_sra)
 
-            output_sra = output_comb.split("Building DAG of jobs...")[2]
-            self.assertTrue("download_sra" not in output_sra)
-            self.assertTrue("sra_qc" in output_sra)
-            self.assertTrue("compile_sra_qc" in output_sra)
-            self.assertTrue("aviary_commands" not in output_sra)
-
-            output = output_comb.split("Building DAG of jobs...")[3]
+            output = output_comb.split("Building DAG of jobs...")[2]
             self.assertTrue("singlem_pipe_reads" not in output)
             self.assertTrue("genome_transcripts" not in output)
             self.assertTrue("singlem_pipe_genomes" not in output)
@@ -265,27 +257,27 @@ def test_update_sra_download_mock(self):
             )
             extern.run(cmd)
 
-            sra_f0_path = os.path.join("test", "coassemble", "sra_qc", "SRR3309137_1.fastq.gz")
+            sra_f0_path = os.path.join("test", "coassemble", "sra", "SRR3309137_1.fastq.gz")
             self.assertTrue(os.path.exists(sra_f0_path))
             with gzip.open(sra_f0_path) as f:
                 file = f.readline().decode()
                 self.assertTrue("@SRR3309137.1 HISEQ06:195:D1DRHACXX:5:1101:1597:2236/1" in file)
                 self.assertTrue("@SRR3309137.2 HISEQ06:195:D1DRHACXX:5:1101:1597:2236/1" not in file)
 
-            sra_r0_path = os.path.join("test", "coassemble", "sra_qc", "SRR3309137_2.fastq.gz")
+            sra_r0_path = os.path.join("test", "coassemble", "sra", "SRR3309137_2.fastq.gz")
             self.assertTrue(os.path.exists(sra_r0_path))
             with gzip.open(sra_r0_path) as f:
                 file = f.readline().decode()
                 self.assertTrue("@SRR3309137.2 HISEQ06:195:D1DRHACXX:5:1101:1597:2236/1" in file)
                 self.assertTrue("@SRR3309137.1 HISEQ06:195:D1DRHACXX:5:1101:1597:2236/1" not in file)
 
-            sra_f1_path = os.path.join("test", "coassemble", "sra_qc", "SRR8334323_1.fastq.gz")
+            sra_f1_path = os.path.join("test", "coassemble", "sra", "SRR8334323_1.fastq.gz")
             self.assertTrue(os.path.exists(sra_f1_path))
             with gzip.open(sra_f1_path) as f:
                 file = f.readline().decode()
-                self.assertTrue("@SEQ_ID.1" not in file)
+                self.assertTrue("@SEQ_ID.1" in file)
 
-            sra_f2_path = os.path.join("test", "coassemble", "sra_qc", "SRR8334324_1.fastq.gz")
+            sra_f2_path = os.path.join("test", "coassemble", "sra", "SRR8334324_1.fastq.gz")
             self.assertTrue(os.path.exists(sra_f2_path))
             with gzip.open(sra_f2_path) as f:
                 file = f.readline().decode()
@@ -295,12 +287,12 @@ def test_update_sra_download_mock(self):
             self.assertTrue(os.path.exists(recover_path))
             with open(recover_path) as f:
                 file = f.readline()
-                self.assertTrue("sra_qc/SRR3309137_1.fastq.gz" in file)
-                self.assertTrue("sra_qc/SRR3309137_2.fastq.gz" in file)
-                self.assertTrue("sra_qc/SRR8334323_1.fastq.gz" in file)
-                self.assertTrue("sra_qc/SRR8334323_2.fastq.gz" in file)
-                self.assertTrue("sra_qc/SRR8334324_1.fastq.gz" in file)
-                self.assertTrue("sra_qc/SRR8334324_2.fastq.gz" in file)
+                self.assertTrue("sra/SRR3309137_1.fastq.gz" in file)
+                self.assertTrue("sra/SRR3309137_2.fastq.gz" in file)
+                self.assertTrue("sra/SRR8334323_1.fastq.gz" in file)
+                self.assertTrue("sra/SRR8334323_2.fastq.gz" in file)
+                self.assertTrue("sra/SRR8334324_1.fastq.gz" in file)
+                self.assertTrue("sra/SRR8334324_2.fastq.gz" in file)
 
     def test_update_sra_download_mock_filter_single(self):
         with in_tempdir():
@@ -391,17 +383,9 @@ def test_update_aviary_run(self):
 
             output_sra = output_comb.split("Building DAG of jobs...")[1]
             self.assertTrue("download_sra" in output_sra)
-            self.assertTrue("sra_qc" not in output_sra)
-            self.assertTrue("compile_sra_qc" not in output_sra)
-            self.assertTrue("aviary_commands" not in output_sra)
-
-            output_sra = output_comb.split("Building DAG of jobs...")[2]
-            self.assertTrue("download_sra" not in output_sra)
-            self.assertTrue("sra_qc" in output_sra)
-            self.assertTrue("compile_sra_qc" in output_sra)
             self.assertTrue("aviary_commands" not in output_sra)
 
-            output = output_comb.split("Building DAG of jobs...")[3]
+            output = output_comb.split("Building DAG of jobs...")[2]
             self.assertTrue("singlem_pipe_reads" not in output)
             self.assertTrue("genome_transcripts" not in output)
             self.assertTrue("singlem_pipe_genomes" not in output)
@@ -515,12 +499,6 @@ def test_update_duplicate_contig_names(self):
                 f"--snakemake-args \" --notemp\" "
             )
             extern.run(cmd)
-            # try:
-            #     import subprocess
-            #     subprocess.run(cmd, shell=True, check=True)
-            # except:
-            #     pass
-            # import pdb; pdb.set_trace()
 
             bins_reference_path = os.path.join("test", "coassemble", "mapping", "sample_1_reference.fna")
             self.assertTrue(os.path.exists(bins_reference_path))