Repository: AroneyS/binchicken (language: python)

Issue:
Fix QC
- [x] Update minimum read length to 80 (some SRA are 96x96 reads)
- [x] QC even without unmapping
- [ ] Check QC'd file size: if zero then warning, if a coassemble sample, then error ~ Not easy since QC is now in main Snakemake pipeline. Change to 80 min len should help.
- [ ] Kingfisher SRA download with download files as output so they are removed with failure ~ Can't do since output can be *1.fastq.gz + *2.fastq.gz or just *.fastq.gz


Discussion / hints:

Fix patch:
diff --git a/ibis/config/template_coassemble.yaml b/ibis/config/template_coassemble.yaml
index e64c2a26..323ac0e3 100644
--- a/ibis/config/template_coassemble.yaml
+++ b/ibis/config/template_coassemble.yaml
@@ -5,7 +5,7 @@ genomes:
 taxa_of_interest: ""
 max_coassembly_size: 
 assemble_unmapped: false
-read_qc: false
+run_qc: false
 sra: false
 prodigal_meta: false
 single_assembly: false
diff --git a/ibis/ibis.py b/ibis/ibis.py
index 3aa17b4e..569da4b1 100755
--- a/ibis/ibis.py
+++ b/ibis/ibis.py
@@ -157,8 +157,13 @@ def download_sra(args):
     expected_reverse = [sra_dir + f + "_2" + SRA_SUFFIX for f in args.forward]
 
     # Fix single file outputs if interleaved or error if unpaired
-    single_ended = {}
-    if not args.dryrun and args.sra != "build":
+    if os.path.isfile(sra_dir + "single_ended.tsv"):
+        with open(sra_dir + "single_ended.tsv") as f:
+            single_ended = {line.split("\t")[0]: line.split("\t")[1].strip() for line in f.readlines()[1:]}
+    else:
+        single_ended = {}
+
+    if not single_ended and not args.dryrun and args.sra != "build":
         for f, r in zip(expected_forward, expected_reverse):
             if os.path.isfile(f) & os.path.isfile(r):
                 continue
@@ -1108,8 +1113,6 @@ def coassemble_argument_verification(args):
             raise Exception("Input SingleM query (--sample-query) requires SingleM otu tables (--sample-singlem) for coverage")
         if args.assemble_unmapped and args.single_assembly:
             raise Exception("Assemble unmapped is incompatible with single-sample assembly")
-        if args.run_qc and not args.assemble_unmapped:
-            raise Exception("Run QC requires unmapping (--assemble-unmapped)")
         if args.assemble_unmapped and not args.genomes and not args.genomes_list:
             raise Exception("Reference genomes must be provided to assemble unmapped reads")
         if not args.singlem_metapackage and not os.environ['SINGLEM_METAPACKAGE_PATH'] and not args.sample_query and not args.sample_query_list:
diff --git a/ibis/workflow/coassemble.smk b/ibis/workflow/coassemble.smk
index a538fbc0..efa055de 100644
--- a/ibis/workflow/coassemble.smk
+++ b/ibis/workflow/coassemble.smk
@@ -16,6 +16,9 @@ logs_dir = output_dir + "/logs"
 mapped_reads_1 = {read: output_dir + f"/mapping/{read}_unmapped.1.fq.gz" for read in config["reads_1"]}
 mapped_reads_2 = {read: output_dir + f"/mapping/{read}_unmapped.2.fq.gz" for read in config["reads_2"]}
 
+qc_reads_1 = {read: output_dir + f"/qc/{read}_1.fastq.gz" for read in config["reads_1"]}
+qc_reads_2 = {read: output_dir + f"/qc/{read}_2.fastq.gz" for read in config["reads_2"]}
+
 def get_genomes(wildcards, version=None):
     version = version if version else wildcards.version
     if version == "":
@@ -37,6 +40,11 @@ def get_reads(wildcards, forward=True, version=None):
             return mapped_reads_1
         else:
             return mapped_reads_2
+    elif version == "qc_":
+        if forward:
+            return qc_reads_1
+        else:
+            return qc_reads_2
     else:
         raise ValueError("Version should be empty, 'whole' or 'unmapped'")
 
@@ -56,6 +64,8 @@ def get_reads_coassembly(wildcards, forward=True, recover=False):
 
     if config["assemble_unmapped"]:
         version = "unmapped_"
+    elif config["run_qc"]:
+        version = "qc_"
     else:
         version = "whole"
 
@@ -430,7 +440,7 @@ rule qc_reads:
     params:
         quality_cutoff = 15,
         unqualified_percent_limit = 40,
-        min_length = 100,
+        min_length = 80,
     threads: 16
     resources:
         mem_mb=125*1000,
@@ -555,20 +565,29 @@ rule finish_mapping:
     shell:
         "touch {output}"
 
+rule finish_qc:
+    input:
+        qc_reads_1.values(),
+        qc_reads_2.values(),
+    output:
+        output_dir + "/qc/done"
+    shell:
+        "touch {output}"
+
 ##############################
 ### Create Aviary commands ###
 ##############################
 rule aviary_commands:
     input:
-        output_dir + "/mapping/done" if config["assemble_unmapped"] else [],
+        output_dir + "/mapping/done" if config["assemble_unmapped"] else output_dir + "/qc/done" if config["run_qc"] else [],
         elusive_clusters = output_dir + "/target/elusive_clusters.tsv",
     output:
         coassemble_commands = output_dir + "/commands/coassemble_commands.sh",
         recover_commands = output_dir + "/commands/recover_commands.sh"
     threads: 8
     params:
-        reads_1 = mapped_reads_1 if config["assemble_unmapped"] else config["reads_1"],
-        reads_2 = mapped_reads_2 if config["assemble_unmapped"] else config["reads_2"],
+        reads_1 = mapped_reads_1 if config["assemble_unmapped"] else qc_reads_1 if config["run_qc"] else config["reads_1"],
+        reads_2 = mapped_reads_2 if config["assemble_unmapped"] else qc_reads_2 if config["run_qc"] else config["reads_2"],
         dir = output_dir,
         memory = config["aviary_memory"],
         threads = config["aviary_threads"],
@@ -583,7 +602,7 @@ rule aviary_commands:
 #########################################
 rule aviary_assemble:
     input:
-        output_dir + "/mapping/done" if config["assemble_unmapped"] else [],
+        output_dir + "/mapping/done" if config["assemble_unmapped"] else output_dir + "/qc/done" if config["run_qc"] else [],
         elusive_clusters = output_dir + "/target/elusive_clusters.tsv",
     output:
         dir = directory(output_dir + "/coassemble/{coassembly}/assemble"),


Test patch:
diff --git a/test/test_coassemble.py b/test/test_coassemble.py
index 7fb81df6..b126fe04 100644
--- a/test/test_coassemble.py
+++ b/test/test_coassemble.py
@@ -235,6 +235,7 @@ def test_coassemble_taxa_of_interest(self):
                 f"--output test "
                 f"--taxa-of-interest \"p__Actinobacteriota\" "
                 f"--aviary-speed comprehensive "
+                f"--run-qc "
                 f"--conda-prefix {path_to_conda} "
             )
             extern.run(cmd)
@@ -275,6 +276,21 @@ def test_coassemble_taxa_of_interest(self):
             with open(cluster_path) as f:
                 self.assertEqual(expected, f.read())
 
+            qc_sample_1F_path = os.path.join("test", "coassemble", "qc", "sample_1_1.fastq.gz")
+            self.assertTrue(os.path.exists(qc_sample_1F_path))
+            qc_sample_1R_path = os.path.join("test", "coassemble", "qc", "sample_1_2.fastq.gz")
+            self.assertTrue(os.path.exists(qc_sample_1R_path))
+
+            qc_sample_2F_path = os.path.join("test", "coassemble", "qc", "sample_2_1.fastq.gz")
+            self.assertTrue(os.path.exists(qc_sample_2F_path))
+            qc_sample_2R_path = os.path.join("test", "coassemble", "qc", "sample_2_2.fastq.gz")
+            self.assertTrue(os.path.exists(qc_sample_2R_path))
+
+            qc_sample_3F_path = os.path.join("test", "coassemble", "qc", "sample_3_1.fastq.gz")
+            self.assertTrue(os.path.exists(qc_sample_3F_path))
+            qc_sample_3R_path = os.path.join("test", "coassemble", "qc", "sample_3_2.fastq.gz")
+            self.assertTrue(os.path.exists(qc_sample_3R_path))
+
             coassemble_path = os.path.join("test", "coassemble", "commands", "coassemble_commands.sh")
             self.assertTrue(os.path.exists(coassemble_path))
             test_dir = os.path.abspath("test")
@@ -282,11 +298,11 @@ def test_coassemble_taxa_of_interest(self):
                 [
                     " ".join([
                         "aviary assemble --coassemble -1",
-                        os.path.join(os.path.abspath(path_to_data), "sample_1.1.fq"),
-                        os.path.join(os.path.abspath(path_to_data), "sample_3.1.fq"),
+                        os.path.join(test_dir, "coassemble", "qc", "sample_1_1.fastq.gz"),
+                        os.path.join(test_dir, "coassemble", "qc", "sample_3_1.fastq.gz"),
                         "-2",
-                        os.path.join(os.path.abspath(path_to_data), "sample_1.2.fq"),
-                        os.path.join(os.path.abspath(path_to_data), "sample_3.2.fq"),
+                        os.path.join(test_dir, "coassemble", "qc", "sample_1_2.fastq.gz"),
+                        os.path.join(test_dir, "coassemble", "qc", "sample_3_2.fastq.gz"),
                         "--output", os.path.join(test_dir, "coassemble", "coassemble", "coassembly_0", "assemble"),
                         "-n 64 -t 64 -m 500 --skip-qc &>",
                         os.path.join(test_dir, "coassemble", "coassemble", "logs", "coassembly_0_assemble.log"),
@@ -305,11 +321,11 @@ def test_coassemble_taxa_of_interest(self):
                     " ".join([
                         "aviary recover --assembly", os.path.join(test_dir, "coassemble", "coassemble", "coassembly_0", "assemble", "assembly", "final_contigs.fasta"),
                         "-1",
-                        os.path.join(os.path.abspath(path_to_data), "sample_1.1.fq"),
-                        os.path.join(os.path.abspath(path_to_data), "sample_3.1.fq"),
+                        os.path.join(test_dir, "coassemble", "qc", "sample_1_1.fastq.gz"),
+                        os.path.join(test_dir, "coassemble", "qc", "sample_3_1.fastq.gz"),
                         "-2",
-                        os.path.join(os.path.abspath(path_to_data), "sample_1.2.fq"),
-                        os.path.join(os.path.abspath(path_to_data), "sample_3.2.fq"),
+                        os.path.join(test_dir, "coassemble", "qc", "sample_1_2.fastq.gz"),
+                        os.path.join(test_dir, "coassemble", "qc", "sample_3_2.fastq.gz"),
                         "--output", os.path.join(test_dir, "coassemble", "coassemble", "coassembly_0", "recover"),
                         "-n 32 -t 32 -m 250 --skip-qc &>",
                         os.path.join(test_dir, "coassemble", "coassemble", "logs", "coassembly_0_recover.log"),
diff --git a/test/test_update.py b/test/test_update.py
index 0a2d29c0..42019dbb 100644
--- a/test/test_update.py
+++ b/test/test_update.py
@@ -341,6 +341,54 @@ def test_update_sra_download_mock_filter_single(self):
             with open(elusive_clusters_path) as f:
                 self.assertEqual(expected, f.read())
 
+    def test_update_sra_download_mock_filter_single_rerun(self):
+        with in_tempdir():
+            os.makedirs("test/coassemble/sra")
+            with open("test/coassemble/sra/single_ended.tsv", "w") as f:
+                f.write("\t".join(["sra", "reason"]) + "\n")
+                f.write("\t".join(["SRR3309137", "Consecutive reads do not match (1/10)"]) + "\n")
+
+            cmd = (
+                f"ibis update "
+                f"--forward SRR3309137 SRR8334323 SRR8334324 SRR3309137_mismatched "
+                f"--sra "
+                f"--genomes {GENOMES} "
+                f"--coassemble-unbinned {MOCK_UNBINNED_SRA} "
+                f"--coassemble-binned {MOCK_BINNED_SRA} "
+                f"--coassemble-targets {MOCK_TARGETS} "
+                f"--coassemble-elusive-edges {MOCK_ELUSIVE_EDGES} "
+                f"--coassemble-elusive-clusters {MOCK_ELUSIVE_CLUSTERS_SRA_MOCK} "
+                f"--coassemble-summary {MOCK_SUMMARY} "
+                f"--output test "
+                f"--conda-prefix {path_to_conda} "
+                f"--snakemake-args \" --config mock_sra=True run_qc=False\" "
+            )
+            extern.run(cmd)
+
+            single_ended_path = os.path.join("test", "coassemble", "sra", "single_ended.tsv")
+            self.assertTrue(os.path.exists(single_ended_path))
+            expected = "\n".join(
+                [
+                    "\t".join(["sra", "reason"]),
+                    "\t".join(["SRR3309137", "Consecutive reads do not match (1/10)"]),
+                    ""
+                ]
+            )
+            with open(single_ended_path) as f:
+                self.assertEqual(expected, f.read())
+
+            elusive_clusters_path = os.path.join("test", "coassemble", "target", "elusive_clusters.tsv")
+            self.assertTrue(os.path.exists(elusive_clusters_path))
+            expected = "\n".join(
+                [
+                    "\t".join(["samples", "length", "total_targets", "total_size", "recover_samples", "coassembly"]),
+                    "\t".join(["SRR8334324,SRR8334323", "2", "2", "0", "SRR8334323,SRR8334324", "coassembly_0"]),
+                    ""
+                ]
+            )
+            with open(elusive_clusters_path) as f:
+                self.assertEqual(expected, f.read())
+
     def test_update_sra_download_mock_fail(self):
         with in_tempdir():
             cmd = (