Repository: AroneyS/binchicken (language: python)

Issue:
Update SingleM
Switch to conda release v0.15.0


Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index 105ba1e3..c0dfcf6d 100644
--- a/README.md
+++ b/README.md
@@ -46,6 +46,7 @@ The taxa of the considered sequences can be filtered to target a specific taxon
 The graph is clustered using the Girvan-Newman algorithm to provide sample groupings.
 Aviary assemble/recover commands are generated based on proposed coassemblies.
 Optionally, reads can be mapped to the matched bins with only unmapped reads being assembled.
+Paired end reads of form reads_1.1.fq, reads_1_1.fq and reads_1_R1.fq are automatically detected and matched to their basename.
 
 ```bash
 # Example: cluster reads into proposed coassemblies
diff --git a/ibis.yml b/ibis.yml
index 60975b46..2d83a0eb 100644
--- a/ibis.yml
+++ b/ibis.yml
@@ -12,7 +12,7 @@ dependencies:
   - bird_tool_utils_python=0.4.*
   - extern=0.4.*
   - ruamel.yaml=0.17.*
-  - polars=0.18.*
+  - polars=0.19.7
   - pigz=2.3.*
   - pyarrow=12.0.*
   - parallel=20230522
diff --git a/ibis/ibis.py b/ibis/ibis.py
index 7e5bbf5b..3aa17b4e 100755
--- a/ibis/ibis.py
+++ b/ibis/ibis.py
@@ -744,7 +744,7 @@ def iterate(args):
 
             if comb_cluster.height > 0:
                 _ = comb_cluster.select(
-                    pl.col("coassembly").apply(lambda x: logging.warn(f"{x} has been previously suggested"))
+                    pl.col("coassembly").map_elements(lambda x: logging.warn(f"{x} has been previously suggested"))
                     )
     elif not args.exclude_coassemblies:
         logging.warn("Suggested coassemblies may match those from previous iterations. To check, use `--elusive-clusters`. To exclude manually, use `--exclude-coassembles`")
@@ -1079,7 +1079,7 @@ def add_coassemble_arguments(parser):
     logging.info(f"Ibis v{__version__}")
     logging.info(f"Command: {' '.join(['ibis'] + sys.argv[1:])}")
 
-    os.environ["POLARS_MAX_THREADS"] = str(args.cores)
+    os.environ["POLARS_MAX_THREADS"] = str(args.local_cores)
     import polars as pl
 
     args.output = os.path.abspath(args.output)
diff --git a/ibis/workflow/coassemble.smk b/ibis/workflow/coassemble.smk
index 6ff88dde..a538fbc0 100644
--- a/ibis/workflow/coassemble.smk
+++ b/ibis/workflow/coassemble.smk
@@ -6,6 +6,7 @@ ruleorder: mock_download_sra > download_sra
 localrules: all, summary, singlem_summarise_genomes, singlem_appraise_filtered, no_genomes, download_sra, mock_download_sra, collect_genomes, finish_mapping, aviary_commands, aviary_recover, aviary_combine
 
 import os
+os.environ['OPENBLAS_NUM_THREADS'] = '1'
 import pandas as pd
 from ibis.ibis import FAST_AVIARY_MODE
 
diff --git a/ibis/workflow/env/singlem.yml b/ibis/workflow/env/singlem.yml
index d25b3b7e..4dff6012 100644
--- a/ibis/workflow/env/singlem.yml
+++ b/ibis/workflow/env/singlem.yml
@@ -3,34 +3,4 @@ channels:
   - bioconda
 dependencies:
   - python>=3.7
-  - diamond>2.0.11,<2.1.3
-  - biopython
-  - hmmer
-  - orfm
-  - mfqe
-  - extern
-  - graftm>=0.14.0
-  - krona
-  - pplacer
-  - sra-tools
-  - ncbi-ngs-sdk
-  - sqlite
-  - squarify
-  - nose
-  - mafft
-  - seqmagick
-  - expressbetadiversity
-  - cd-hit
-  - fasttree
-  - fastalite
-  - jinja2
-  - pip
-  - sqlalchemy
-  - pandas
-  - bird_tool_utils_python>=0.4.1
-  - zenodo_backpack
-  - smafa>=0.7.0
-  - git
-  - git-lfs
-  - pip:
-    - git+https://github.com/wwood/singlem.git@8074a40
+  - singlem>=0.15.*
diff --git a/ibis/workflow/evaluate.smk b/ibis/workflow/evaluate.smk
index 8456cc7a..8824ffb1 100644
--- a/ibis/workflow/evaluate.smk
+++ b/ibis/workflow/evaluate.smk
@@ -1,8 +1,9 @@
 #############
 ### Setup ###
 #############
-import pandas as pd
 import os
+os.environ['OPENBLAS_NUM_THREADS'] = '1'
+import pandas as pd
 import re
 
 output_dir = os.path.abspath("evaluate")
diff --git a/ibis/workflow/scripts/cluster_graph.py b/ibis/workflow/scripts/cluster_graph.py
index 65feaa35..9fa0f1a3 100644
--- a/ibis/workflow/scripts/cluster_graph.py
+++ b/ibis/workflow/scripts/cluster_graph.py
@@ -42,7 +42,7 @@ def join_list_subsets(df1, df2):
             how="left",
         )
         .filter(pl.col("samples").list.set_intersection("right_samples").list.lengths() >= pl.col("length"))
-        .groupby("samples_hash")
+        .group_by("samples_hash")
         .agg(pl.col("extra_targets").flatten())
         .select(
             pl.col("samples_hash"),
@@ -79,12 +79,12 @@ def find_recover_candidates(df, samples_df, MAX_RECOVERY_SAMPLES=20):
         .select("samples_hash", "target_ids")
         .explode("target_ids")
         .join(samples_df, on="target_ids", how="left")
-        .groupby("samples_hash", "recover_candidates")
+        .group_by("samples_hash", "recover_candidates")
         .agg(pl.count())
         .sort("count", descending=True)
-        .groupby("samples_hash")
+        .group_by("samples_hash")
         .head(MAX_RECOVERY_SAMPLES)
-        .groupby("samples_hash")
+        .group_by("samples_hash")
         .agg("recover_candidates")
     )
 
@@ -154,7 +154,7 @@ def pipeline(
             clusters = [
                 elusive_edges
                 .explode("samples")
-                .groupby("samples")
+                .group_by("samples")
                 .agg(pl.col("target_ids").flatten())
                 .with_columns(
                     pl.concat_list(pl.col("samples")),
@@ -180,7 +180,7 @@ def pipeline(
                     # Prevent combinatorial explosion (also, large clusters are less useful for distinguishing between clusters)
                     .filter(pl.col("samples").list.lengths() < MAX_SAMPLES_COMBINATIONS)
                     .with_columns(
-                        samples_combinations = pl.struct(["length", "cluster_size"]).apply(
+                        samples_combinations = pl.struct(["length", "cluster_size"]).map_elements(
                             lambda x: [i for i in itertools.combinations(range(x["length"]), x["cluster_size"])],
                             return_dtype=pl.List(pl.List(pl.Int64)),
                             )
@@ -189,7 +189,7 @@ def pipeline(
                     .with_columns(pl.col("samples").list.take(pl.col("samples_combinations")))
                     .with_columns(samples_hash = pl.col("samples").list.sort().hash())
                     .select("samples_hash", "samples", "target_ids", "cluster_size")
-                    .groupby("samples_hash")
+                    .group_by("samples_hash")
                     .agg(
                         pl.first("samples"),
                         pl.col("target_ids").flatten(),
@@ -205,7 +205,7 @@ def pipeline(
             .explode("recover_candidates")
             .explode("target_ids")
             .unique()
-            .groupby("recover_candidates")
+            .group_by("recover_candidates")
             .agg("target_ids")
         )
 
@@ -216,7 +216,7 @@ def pipeline(
             .with_columns(length = pl.col("samples").list.lengths())
             .explode("samples")
             .join(read_size, on="samples", how="left")
-            .groupby("samples_hash")
+            .group_by("samples_hash")
             .agg(
                 pl.col("samples"),
                 pl.first("length"),
@@ -235,7 +235,7 @@ def pipeline(
             .with_columns(
                 unique_samples = 
                     pl.col("samples")
-                    .map(accumulate_clusters, return_dtype=pl.Boolean),
+                    .map_batches(accumulate_clusters, return_dtype=pl.Boolean),
                 )
             .filter(pl.col("unique_samples"))
             .drop("unique_samples")
diff --git a/ibis/workflow/scripts/collect_reference_bins.py b/ibis/workflow/scripts/collect_reference_bins.py
index d5ef55dc..c8062671 100755
--- a/ibis/workflow/scripts/collect_reference_bins.py
+++ b/ibis/workflow/scripts/collect_reference_bins.py
@@ -51,7 +51,7 @@ def pipeline(appraise_binned, appraise_unbinned, sample, MIN_APPRAISED=0.1, TRIM
         "found_in"
     ).with_columns(
         pl.col("found_in").str.replace("_protein$", "")
-    ).groupby(
+    ).group_by(
         ["gene", "found_in"]
     ).agg(
         pl.col("coverage").sum()
@@ -60,7 +60,7 @@ def pipeline(appraise_binned, appraise_unbinned, sample, MIN_APPRAISED=0.1, TRIM
     ).melt(
         id_vars="gene", variable_name="found_in", value_name="coverage"
     ).fill_null(0
-    ).groupby(
+    ).group_by(
         "found_in"
     ).agg(
         (pl.col("coverage").len() * TRIM_FRACTION).floor().cast(int).alias("cut"),
diff --git a/ibis/workflow/scripts/evaluate.py b/ibis/workflow/scripts/evaluate.py
index 62691af9..6ad8c375 100755
--- a/ibis/workflow/scripts/evaluate.py
+++ b/ibis/workflow/scripts/evaluate.py
@@ -65,7 +65,7 @@ def evaluate(target_otu_table, binned_otu_table, elusive_clusters, elusive_edges
             )
         .explode("samples")
         .join(sample_coassemblies, on="samples")
-        .groupby("samples_hash", "coassembly")
+        .group_by("samples_hash", "coassembly")
         .agg(
             pl.first("target_ids"),
             pl.count(),
@@ -86,7 +86,7 @@ def evaluate(target_otu_table, binned_otu_table, elusive_clusters, elusive_edges
     # Create otu table with original sequence, samples present, cluster id, target id and associated coassemblies
     sample_edges = (
         elusive_edges
-        .groupby("coassembly", "target")
+        .group_by("coassembly", "target")
         .agg(
             source_samples = pl.col("samples")
                 .flatten()
@@ -122,7 +122,7 @@ def evaluate(target_otu_table, binned_otu_table, elusive_clusters, elusive_edges
             ])
         .join(sample_coassemblies, left_on="sample", right_on="samples", how="left")
         .drop_nulls("coassembly")
-        .groupby(["gene", "sequence", "coassembly"])
+        .group_by(["gene", "sequence", "coassembly"])
         .agg([
             pl.first("taxonomy"),
             pl.first("found_in"),
@@ -200,7 +200,7 @@ def summarise_stats(matches, combined_otu_table, recovered_bins):
                 .otherwise(pl.lit("nonmatch"))
                 .alias("status")
             )
-        .groupby([
+        .group_by([
             "coassembly", "status"
             ])
         .agg(
@@ -216,7 +216,7 @@ def summarise_stats(matches, combined_otu_table, recovered_bins):
                         pl.col("found_in").is_not_null()
                     ).then(pl.lit("match")).otherwise(pl.lit("nonmatch")).alias("status")
                     )
-                .groupby(["coassembly", "status"])
+                .group_by(["coassembly", "status"])
                 .agg(pl.col("sequence").len().alias("nontarget_bin_sequences")),
             on=["coassembly", "status"], how="outer"
             )
@@ -228,7 +228,7 @@ def summarise_stats(matches, combined_otu_table, recovered_bins):
                         pl.all_horizontal(pl.col(["target", "found_in"]).is_null()) & (pl.col("source_samples").is_not_null())
                     ).then(pl.lit("match")).otherwise(pl.lit("nonmatch")).alias("status")
                     )
-                .groupby(["coassembly", "status"])
+                .group_by(["coassembly", "status"])
                 .agg(pl.col("sequence").len().alias("nontarget_unbin_sequences")),
             on=["coassembly", "status"], how="outer"
             )
@@ -240,7 +240,7 @@ def summarise_stats(matches, combined_otu_table, recovered_bins):
                         pl.all_horizontal(pl.col(["target", "found_in", "source_samples"]).is_null())
                     ).then(pl.lit("match")).otherwise(pl.lit("nonmatch")).alias("status")
                     )
-                .groupby(["coassembly", "status"])
+                .group_by(["coassembly", "status"])
                 .agg(pl.col("sequence").len().alias("novel_sequences")),
             on=["coassembly", "status"], how="outer"
             )
@@ -255,7 +255,7 @@ def summarise_stats(matches, combined_otu_table, recovered_bins):
             pl.col("name").str.split("-").list.to_struct(upper_bound=2, fields=["coassembly", "genome"])
             )
         .unnest("name")
-        .groupby("coassembly")
+        .group_by("coassembly")
         .agg(pl.count().alias("n"))
         .with_columns(
             pl.lit("match").alias("status"),
diff --git a/ibis/workflow/scripts/query_processing.py b/ibis/workflow/scripts/query_processing.py
index 19c2f32b..005930ef 100755
--- a/ibis/workflow/scripts/query_processing.py
+++ b/ibis/workflow/scripts/query_processing.py
@@ -43,7 +43,7 @@ def processing(
         "taxonomy", "num_hits", "coverage"
     ]).join(
         pipe_read, on=["gene", "sample", "sequence"], how="inner"
-    ).groupby(
+    ).group_by(
         ["gene", "sample", "sequence", "num_hits", "coverage", "taxonomy", "divergence"]
     ).agg(
         pl.col("found_in").sort().str.concat(",")
diff --git a/ibis/workflow/scripts/separate_SingleM_seq.py b/ibis/workflow/scripts/separate_SingleM_seq.py
index 99db5b96..9176cda7 100755
--- a/ibis/workflow/scripts/separate_SingleM_seq.py
+++ b/ibis/workflow/scripts/separate_SingleM_seq.py
@@ -3,6 +3,8 @@
 ###############################
 # Author: Samuel Aroney
 
+import os
+os.environ['OPENBLAS_NUM_THREADS'] = '1'
 import pandas as pd
 import re
 from Bio.Seq import Seq
diff --git a/ibis/workflow/scripts/summarise_coassemblies.py b/ibis/workflow/scripts/summarise_coassemblies.py
index 60c03b2f..4c462eb9 100644
--- a/ibis/workflow/scripts/summarise_coassemblies.py
+++ b/ibis/workflow/scripts/summarise_coassemblies.py
@@ -3,6 +3,8 @@
 #################################
 # Author: Samuel Aroney
 
+import os
+os.environ['OPENBLAS_NUM_THREADS'] = '1'
 import pandas as pd
 
 elusive_clusters = pd.read_csv(snakemake.input.elusive_clusters, sep="\t")
diff --git a/ibis/workflow/scripts/target_elusive.py b/ibis/workflow/scripts/target_elusive.py
index 05ea3629..0ab87f94 100755
--- a/ibis/workflow/scripts/target_elusive.py
+++ b/ibis/workflow/scripts/target_elusive.py
@@ -44,7 +44,7 @@ def pipeline(
         .with_columns(
             pl.when(pl.col("sample").is_in(samples))
             .then(pl.col("sample"))
-            .otherwise(pl.col("sample").str.replace(r"(_|\.)1$", ""))
+            .otherwise(pl.col("sample").str.replace(r"(_|\.)R?1$", ""))
             )
         .filter(pl.col("sample").is_in(samples))
         .drop("found_in")
@@ -58,6 +58,10 @@ def pipeline(
             )
     )
 
+    if unbinned.height == 0:
+        logging.warning("No SingleM sequences found for the given samples")
+        return unbinned, pl.DataFrame(schema=EDGES_COLUMNS)
+
     def process_groups(df):
         if df.height == 1:
             return pl.DataFrame(schema={"style": str, "cluster_size": pl.Int64, "samples": str, "target": str})
@@ -90,7 +94,7 @@ def process_groups(df):
             df
             .join(cluster_sizes, how="cross")
             .filter(pl.col("coverage") > float(MIN_COASSEMBLY_COVERAGE) / pl.col("cluster_size").cast(float))
-            .groupby("target", "cluster_size")
+            .group_by("target", "cluster_size")
             .agg(pl.concat_list("samples").flatten())
             .filter(pl.col("samples").list.lengths() >= pl.col("cluster_size"))
             .select(
@@ -109,11 +113,11 @@ def process_groups(df):
         .select(
             "target",
             "coverage",
-            samples = pl.col("sample").apply(lambda x: [x], return_dtype=pl.List(pl.Utf8)),
+            samples = pl.col("sample").map_elements(lambda x: [x], return_dtype=pl.List(pl.Utf8)),
             )
-        .groupby("target")
-        .apply(process_groups)
-        .groupby(["style", "cluster_size", "samples"])
+        .group_by("target")
+        .map_groups(process_groups)
+        .group_by(["style", "cluster_size", "samples"])
         .agg(target_ids = pl.col("target").sort().str.concat(","))
     )
 


Test patch:
diff --git a/test/test_target_elusive.py b/test/test_target_elusive.py
index 1c6aa403..7a9767f0 100644
--- a/test/test_target_elusive.py
+++ b/test/test_target_elusive.py
@@ -75,7 +75,7 @@ def test_target_elusive_suffix(self):
         self.assertDataFrameEqual(expected_targets, observed_targets)
         self.assertDataFrameEqual(expected_edges, observed_edges)
 
-    def test_target_elusive_suffix(self):
+    def test_target_elusive_suffix_underscore(self):
         unbinned = pl.DataFrame([
             ["S3.1", "sample_1_1", "AAA", 5, 10, "Root", ""],
             ["S3.1", "sample_2_1", "AAA", 5, 10, "Root", ""],
@@ -94,6 +94,25 @@ def test_target_elusive_suffix(self):
         self.assertDataFrameEqual(expected_targets, observed_targets)
         self.assertDataFrameEqual(expected_edges, observed_edges)
 
+    def test_target_elusive_suffix_R(self):
+        unbinned = pl.DataFrame([
+            ["S3.1", "sample_1_R1", "AAA", 5, 10, "Root", ""],
+            ["S3.1", "sample_2_R1", "AAA", 5, 10, "Root", ""],
+        ], schema=APPRAISE_COLUMNS)
+        samples = set(["sample_1", "sample_2"])
+
+        expected_targets = pl.DataFrame([
+            ["S3.1", "sample_1", "AAA", 5, 10, "Root", "0"],
+            ["S3.1", "sample_2", "AAA", 5, 10, "Root", "0"],
+        ], schema=TARGETS_COLUMNS)
+        expected_edges = pl.DataFrame([
+            ["match", 2, "sample_1,sample_2", "0"],
+        ], schema=EDGES_COLUMNS)
+
+        observed_targets, observed_edges = pipeline(unbinned, samples)
+        self.assertDataFrameEqual(expected_targets, observed_targets)
+        self.assertDataFrameEqual(expected_edges, observed_edges)
+
     def test_target_elusive_empty_input(self):
         unbinned = pl.DataFrame([
         ], schema=APPRAISE_COLUMNS)
@@ -108,6 +127,22 @@ def test_target_elusive_empty_input(self):
         self.assertDataFrameEqual(expected_targets, observed_targets)
         self.assertDataFrameEqual(expected_edges, observed_edges)
 
+    def test_target_elusive_no_matching_samples(self):
+        unbinned = pl.DataFrame([
+            ["S3.1", "sample_1.1", "AAA", 5, 10, "Root", ""],
+            ["S3.1", "sample_2.1", "AAA", 5, 10, "Root", ""],
+        ], schema=APPRAISE_COLUMNS)
+        samples = set(["sample_3", "sample_4"])
+
+        expected_targets = pl.DataFrame([
+        ], schema=TARGETS_COLUMNS)
+        expected_edges = pl.DataFrame([
+        ], schema=EDGES_COLUMNS)
+
+        observed_targets, observed_edges = pipeline(unbinned, samples)
+        self.assertDataFrameEqual(expected_targets, observed_targets)
+        self.assertDataFrameEqual(expected_edges, observed_edges)
+
     def test_target_elusive_low_coverage(self):
         unbinned = pl.DataFrame([
             ["S3.1", "sample_1", "AAA", 5, 5, "Root", ""],