Add cluster metrics viz

modules/nf-core/custom/clustermetrics/environment.yml

@@ -0,0 +1,12 @@
+---
+# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/environment-schema.json
+channels:
+  - conda-forge
+  - bioconda
+dependencies:
+  - conda-forge::matplotlib=3.10.9
+  - conda-forge::numpy=2.4.4
+  - conda-forge::pandas=3.0.3
+  - conda-forge::python=3.12.13
+  - conda-forge::pyyaml=6.0.3
+  - conda-forge::scikit-learn=1.8.0

modules/nf-core/custom/clustermetrics/main.nf

@@ -0,0 +1,46 @@
+process CUSTOM_CLUSTERMETRICS {
+    tag "$meta.id"
+    label 'process_medium'
+
+    conda "${moduleDir}/environment.yml"
+    container "${ workflow.containerEngine in ['singularity', 'apptainer'] && !task.ext.singularity_pull_docker_container ?
+        'https://community-cr-prod.seqera.io/docker/registry/v2/blobs/sha256/25/25129a5258522a434c386b800d3e2e3e6dc72d8a1171b7b10f21df3488526795/data' :
+        'community.wave.seqera.io/library/matplotlib_numpy_pandas_python_pruned:169e228afc7d3686' }"
+
+    input:
+    tuple val(meta), path(features), path(clusters)
+
+    output:
+    tuple val(meta), path("*.metrics.tsv")  , emit: metrics
+    tuple val(meta), path("*.k_sweep.csv")  , emit: k_sweep
+    tuple val(meta), path("*.selected.json"), emit: selected
+    tuple val(meta), path("*.png")          , emit: plots, optional: true
+    path "versions.yml"                     , emit: versions, topic: versions
+
+    when:
+    task.ext.when == null || task.ext.when
+
+    script:
+    template 'cluster_metrics.py'
+
+    stub:
+    def prefix = task.ext.prefix ?: "${meta.id}"
+    """
+    touch ${prefix}.metrics.tsv
+    touch ${prefix}.k_sweep.csv
+    touch ${prefix}.selected.json
+    touch ${prefix}.elbow.png
+    touch ${prefix}.silhouette.png
+    touch ${prefix}.davies_bouldin.png
+    touch ${prefix}.calinski_harabasz.png
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        python: \$(python3 --version | sed 's/Python //')
+        matplotlib: \$(python3 -c "from importlib.metadata import version; print(version('matplotlib'))")
+        numpy: \$(python3 -c "from importlib.metadata import version; print(version('numpy'))")
+        pandas: \$(python3 -c "from importlib.metadata import version; print(version('pandas'))")
+        scikit-learn: \$(python3 -c "from importlib.metadata import version; print(version('scikit-learn'))")
+    END_VERSIONS
+    """
+}

modules/nf-core/custom/clustermetrics/meta.yml

@@ -0,0 +1,96 @@
+name: "CUSTOM_CLUSTERMETRICS"
+description: "Computes clustering quality metrics (silhouette, Calinski-Harabasz,
+  Davies-Bouldin) and performs k-sweep analysis"
+keywords:
+  - clustering
+  - metrics
+  - silhouette
+  - calinski-harabasz
+  - davies-bouldin
+  - evaluation
+tools:
+  - "scikit-learn":
+      description: "Machine learning library for clustering metrics"
+      homepage: "https://scikit-learn.org/"
+      documentation: "https://scikit-learn.org/stable/modules/clustering.html"
+      licence:
+        - "BSD-3-Clause"
+      identifier: ""
+input:
+  - - meta:
+        type: map
+        description: |
+          Groovy Map containing sample information
+          e.g. `[ id:'sample1' ]`
+    - features:
+        type: file
+        description: |
+          Tab-separated feature matrix with a `sample_id` column and one
+          column per numeric feature (e.g. PCA scores).
+        pattern: "*.tsv"
+        ontologies:
+          - edam: http://edamontology.org/format_3475
+    - clusters:
+        type: file
+        description: |
+          Comma-separated cluster assignments with `sample_id` and integer
+          `cluster` columns. Label -1 is treated as DBSCAN noise.
+        pattern: "*.csv"
+        ontologies:
+          - edam: http://edamontology.org/format_3752
+output:
+  metrics:
+    - - meta:
+          type: map
+          description: Groovy Map containing sample information
+      - "*.metrics.tsv":
+          type: file
+          description: TSV with selected cluster quality metrics
+          pattern: "*.metrics.tsv"
+          ontologies:
+            - edam: http://edamontology.org/format_3475
+  k_sweep:
+    - - meta:
+          type: map
+          description: Groovy Map containing sample information
+      - "*.k_sweep.csv":
+          type: file
+          description: CSV with metrics for different values of k
+          pattern: "*.k_sweep.csv"
+          ontologies:
+            - edam: http://edamontology.org/format_3752
+  selected:
+    - - meta:
+          type: map
+          description: Groovy Map containing sample information
+      - "*.selected.json":
+          type: file
+          description: JSON with the selected/best metrics
+          pattern: "*.selected.json"
+          ontologies:
+            - edam: http://edamontology.org/format_3464
+  plots:
+    - - meta:
+          type: map
+          description: Groovy Map containing sample information
+      - "*.png":
+          type: file
+          description: Optional PNG plots (elbow, silhouette, etc.)
+          pattern: "*.png"
+          ontologies: []
+  versions:
+    - "versions.yml":
+        type: file
+        description: File containing software versions
+        pattern: "versions.yml"
+        ontologies:
+          - edam: http://edamontology.org/format_3750
+topics:
+  versions:
+    - versions.yml:
+        type: string
+        description: The name of the process
+authors:
+  - "@dbaku42"
+maintainers:
+  - "@dbaku42"

modules/nf-core/custom/clustermetrics/templates/cluster_metrics.py

@@ -0,0 +1,143 @@
+#!/usr/bin/env python3
+
+# Copyright (c) nf-core
+# This software is licensed under the MIT License.
+# SPDX-License-Identifier: MIT
+
+import argparse
+import json
+import platform
+import shlex
+
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import pandas as pd
+import sklearn
+import yaml
+from sklearn.cluster import KMeans
+from sklearn.metrics import (
+    calinski_harabasz_score,
+    davies_bouldin_score,
+    silhouette_score,
+)
+
+
+def load_features(path):
+    """Read a TSV of `sample_id` + numeric feature columns, indexed by sample_id."""
+    df = pd.read_csv(path, sep="\\t")
+    if "sample_id" not in df.columns:
+        raise ValueError(f"features file must have a 'sample_id' column. Found: {list(df.columns)}")
+    df["sample_id"] = df["sample_id"].astype(str)
+    return df.set_index("sample_id").apply(pd.to_numeric, errors="coerce").fillna(0.0)
+
+
+def load_clusters(path):
+    """Read a CSV of `sample_id` + `cluster`, returning a Series of int labels."""
+    df = pd.read_csv(path)
+    if "sample_id" not in df.columns or "cluster" not in df.columns:
+        raise ValueError(f"clusters file must have 'sample_id' and 'cluster' columns. Found: {list(df.columns)}")
+    df["sample_id"] = df["sample_id"].astype(str)
+    return df.set_index("sample_id")["cluster"].astype(int)
+
+
+def cluster_quality(x, labels):
+    """Silhouette / Calinski-Harabasz / Davies-Bouldin for given (x, labels).
+
+    Treats label -1 as DBSCAN noise and excludes those points. Returns None
+    for each score when fewer than 2 clusters of more than one point remain.
+    """
+    mask = labels != -1
+    x, labels = x[mask], labels[mask]
+    n = len(set(labels))
+    valid = 2 <= n < len(x)
+    return {
+        "silhouette": float(silhouette_score(x, labels)) if valid else None,
+        "calinski_harabasz": float(calinski_harabasz_score(x, labels)) if valid else None,
+        "davies_bouldin": float(davies_bouldin_score(x, labels)) if valid else None,
+    }
+
+
+def plot_curve(sweep_df, metric, title, ylabel, out_png):
+    plt.figure(figsize=(7, 4.5))
+    vals = sweep_df[metric].dropna()
+    ks = sweep_df.loc[vals.index, "k"]
+    plt.plot(ks, vals, marker="o")
+    plt.xticks(sweep_df["k"].tolist())
+    plt.title(title)
+    plt.xlabel("k")
+    plt.ylabel(ylabel)
+    plt.tight_layout()
+    plt.savefig(out_png, dpi=200)
+    plt.close()
+
+
+def main():
+    features = "$features"
+    clusters_path = "$clusters"
+    prefix = "${task.ext.prefix ?: meta.id}"
+
+    # Optional configuration via task.ext.args (nf-core convention).
+    raw_args = "$task.ext.args"
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--k-min", type=int, default=2)
+    parser.add_argument("--k-max", type=int, default=12)
+    opts = parser.parse_args(shlex.split(raw_args) if raw_args and raw_args != "null" else [])
+
+    joined = load_features(features).join(load_clusters(clusters_path), how="inner")
+    if len(joined) < 2:
+        raise ValueError(f"Need at least 2 samples with matching sample_id in both inputs. Got {len(joined)}.")
+
+    labels = joined["cluster"].values
+    x = joined.drop(columns=["cluster"]).to_numpy(dtype=float)
+
+    # Quality metrics on the supplied labels.
+    selected = {"n_clusters": len(set(labels) - {-1}), **cluster_quality(x, labels)}
+    pd.DataFrame([selected]).to_csv(f"{prefix}.metrics.tsv", sep="\\t", index=False)
+    with open(f"{prefix}.selected.json", "w") as fh:
+        json.dump(selected, fh, indent=2)
+
+    # KMeans k-sweep for downstream comparison.
+    rows = []
+    for k in range(opts.k_min, min(opts.k_max, len(x)) + 1):
+        model = KMeans(n_clusters=k, n_init=10, random_state=42).fit(x)
+        rows.append({"k": k, "inertia": float(model.inertia_), **cluster_quality(x, model.labels_)})
+
+    sweep_df = pd.DataFrame(rows)
+    sweep_df.to_csv(f"{prefix}.k_sweep.csv", index=False, float_format="%.10g")
+
+    if not sweep_df.empty:
+        plot_curve(sweep_df, "inertia", "Elbow method (KMeans inertia)", "inertia", f"{prefix}.elbow.png")
+        plot_curve(
+            sweep_df, "silhouette", "Silhouette score (higher is better)", "silhouette", f"{prefix}.silhouette.png"
+        )
+        plot_curve(
+            sweep_df,
+            "davies_bouldin",
+            "Davies-Bouldin index (lower is better)",
+            "davies_bouldin",
+            f"{prefix}.davies_bouldin.png",
+        )
+        plot_curve(
+            sweep_df,
+            "calinski_harabasz",
+            "Calinski-Harabasz index (higher is better)",
+            "calinski_harabasz",
+            f"{prefix}.calinski_harabasz.png",
+        )
+
+    versions = {
+        "${task.process}": {
+            "python": platform.python_version(),
+            "pandas": pd.__version__,
+            "scikit-learn": sklearn.__version__,
+            "matplotlib": matplotlib.__version__,
+        }
+    }
+    with open("versions.yml", "w") as fh:
+        yaml.dump(versions, fh, default_flow_style=False, sort_keys=False)
+
+
+if __name__ == "__main__":
+    main()

modules/nf-core/custom/clustermetrics/tests/main.nf.test

@@ -0,0 +1,70 @@
+nextflow_process {
+
+    name "Test Process CUSTOM_CLUSTERMETRICS"
+    script "../main.nf"
+    process "CUSTOM_CLUSTERMETRICS"
+
+    tag "modules"
+    tag "modules_nfcore"
+    tag "custom"
+    tag "custom/clustermetrics"
+
+    test("clustermetrics - features and clusters") {
+
+        when {
+            process {
+                """
+                input[0] = [
+                    [ id:'test' ],
+                    file(params.modules_testdata_base_path + "genomics/homo_sapiens/popgen/clustering/test_features.tsv", checkIfExists: true),
+                    file(params.modules_testdata_base_path + "genomics/homo_sapiens/popgen/clustering/test_clusters.csv", checkIfExists: true)
+                ]
+                """
+            }
+        }
+
+        then {
+            assertAll(
+                { assert process.success },
+                { assert snapshot(
+                    process.out.metrics,
+                    process.out.k_sweep,
+                    process.out.selected,
+                    process.out.versions,
+                    path(process.out.versions[0]).yaml
+                ).match() }
+            )
+        }
+    }
+
+    test("clustermetrics - features and clusters - stub") {
+
+        options "-stub"
+
+        when {
+            process {
+                """
+                input[0] = [
+                    [ id:'test' ],
+                    file(params.modules_testdata_base_path + "genomics/homo_sapiens/popgen/clustering/test_features.tsv", checkIfExists: true),
+                    file(params.modules_testdata_base_path + "genomics/homo_sapiens/popgen/clustering/test_clusters.csv", checkIfExists: true)
+                ]
+                """
+            }
+        }
+
+        then {
+            assertAll(
+                { assert process.success },
+                { assert snapshot(
+                    process.out.metrics,
+                    process.out.k_sweep,
+                    process.out.selected,
+                    process.out.plots,
+                    process.out.versions,
+                    path(process.out.versions[0]).yaml
+                ).match() }
+            )
+        }
+    }
+}