diff --git a/modules/nf-core/custom/clustermetrics/main.nf b/modules/nf-core/custom/clustermetrics/main.nf index 2eb096b08f40..1ec2bc9c0c97 100644 --- a/modules/nf-core/custom/clustermetrics/main.nf +++ b/modules/nf-core/custom/clustermetrics/main.nf @@ -33,6 +33,7 @@ process CUSTOM_CLUSTERMETRICS { touch ${prefix}.silhouette.png touch ${prefix}.davies_bouldin.png touch ${prefix}.calinski_harabasz.png + touch ${prefix}.k_distance.png cat <<-END_VERSIONS > versions.yml "${task.process}": diff --git a/modules/nf-core/custom/clustermetrics/meta.yml b/modules/nf-core/custom/clustermetrics/meta.yml index e364f55b6f75..9e0272f7d346 100644 --- a/modules/nf-core/custom/clustermetrics/meta.yml +++ b/modules/nf-core/custom/clustermetrics/meta.yml @@ -1,12 +1,13 @@ name: "CUSTOM_CLUSTERMETRICS" description: "Computes clustering quality metrics (silhouette, Calinski-Harabasz, - Davies-Bouldin) and performs k-sweep analysis" + Davies-Bouldin), performs k-sweep and k-distance analysis" keywords: - clustering - metrics - silhouette - calinski-harabasz - davies-bouldin + - k-distance - evaluation tools: - "scikit-learn": @@ -25,9 +26,11 @@ input: - 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" + Sample-by-feature matrix: TSV with a `sample_id` column, or a + PLINK2 `.eigenvec` file (`#FID IID PC1 ...` or `#IID PC1 ...`). + Family/source ID columns are dropped; remaining columns are + numeric features (e.g. PCA scores). + pattern: "*.{tsv,txt,eigenvec}" ontologies: - edam: http://edamontology.org/format_3475 - clusters: diff --git a/modules/nf-core/custom/clustermetrics/templates/cluster_metrics.py b/modules/nf-core/custom/clustermetrics/templates/cluster_metrics.py index fd8856fd41cf..6a19d786d0b5 100644 --- a/modules/nf-core/custom/clustermetrics/templates/cluster_metrics.py +++ b/modules/nf-core/custom/clustermetrics/templates/cluster_metrics.py @@ -13,6 +13,7 @@ matplotlib.use("Agg") import matplotlib.pyplot as plt +import numpy as np import pandas as pd import sklearn import yaml @@ -22,13 +23,28 @@ davies_bouldin_score, silhouette_score, ) +from sklearn.neighbors import NearestNeighbors def load_features(path): - """Read a TSV of `sample_id` + numeric feature columns, indexed by sample_id.""" - df = pd.read_csv(path, sep="\\t") + """Read sample-by-feature matrix, indexed by sample_id. + + Accepts a TSV with a sample_id column, or a PLINK2 .eigenvec file + (`#FID IID PC1 ...` or `#IID PC1 ...`). FID/SID are dropped. + """ + df = pd.read_csv(path, sep="\\t", dtype=str) + df.columns = [str(col).lstrip("#") for col in df.columns] + + if "IID" in df.columns: + df = df.rename(columns={"IID": "sample_id"}) + + drop = [col for col in ("FID", "SID") if col in df.columns] + if drop: + df = df.drop(columns=drop) + 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) @@ -73,6 +89,28 @@ def plot_curve(sweep_df, metric, title, ylabel, out_png): plt.close() +def plot_k_distance(x, k, out_png, eps=None): + """k-distance plot used to choose DBSCAN eps.""" + k = int(min(max(k, 1), len(x) - 1)) + nn = NearestNeighbors(n_neighbors=k + 1) + nn.fit(x) + distances, _ = nn.kneighbors(x) + k_distances = np.sort(distances[:, -1]) + + plt.figure(figsize=(8, 4)) + plt.plot(k_distances) + if eps is not None: + plt.axhline(y=eps, linestyle="--", label=f"eps = {eps}") + plt.legend() + plt.ylabel(f"Distance to {k}-th nearest neighbour") + plt.xlabel("Points sorted by distance") + plt.title("k-distance plot") + plt.grid(True) + plt.tight_layout() + plt.savefig(out_png, dpi=200) + plt.close() + + def main(): features = "$features" clusters_path = "$clusters" @@ -83,6 +121,8 @@ def main(): parser = argparse.ArgumentParser() parser.add_argument("--k-min", type=int, default=2) parser.add_argument("--k-max", type=int, default=12) + parser.add_argument("--k-neighbors", type=int, default=None) + parser.add_argument("--eps", type=float, default=None) 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") @@ -127,6 +167,9 @@ def main(): f"{prefix}.calinski_harabasz.png", ) + k_neighbors = opts.k_neighbors if opts.k_neighbors is not None else x.shape[1] + plot_k_distance(x, k_neighbors, f"{prefix}.k_distance.png", eps=opts.eps) + versions = { "${task.process}": { "python": platform.python_version(), diff --git a/modules/nf-core/custom/clustermetrics/tests/main.nf.test b/modules/nf-core/custom/clustermetrics/tests/main.nf.test index d5a42e682198..8b951c7e1cfe 100644 --- a/modules/nf-core/custom/clustermetrics/tests/main.nf.test +++ b/modules/nf-core/custom/clustermetrics/tests/main.nf.test @@ -37,6 +37,34 @@ nextflow_process { } } + test("clustermetrics - features from eigenvec") { + + when { + process { + """ + input[0] = [ + [ id:'test' ], + file(params.modules_testdata_base_path + "genomics/homo_sapiens/popgen/clustering/test.eigenvec", 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" diff --git a/modules/nf-core/custom/clustermetrics/tests/main.nf.test.snap b/modules/nf-core/custom/clustermetrics/tests/main.nf.test.snap index 0ec0df908bc5..6bdff3a47c01 100644 --- a/modules/nf-core/custom/clustermetrics/tests/main.nf.test.snap +++ b/modules/nf-core/custom/clustermetrics/tests/main.nf.test.snap @@ -78,6 +78,7 @@ "test.calinski_harabasz.png:md5,d41d8cd98f00b204e9800998ecf8427e", "test.davies_bouldin.png:md5,d41d8cd98f00b204e9800998ecf8427e", "test.elbow.png:md5,d41d8cd98f00b204e9800998ecf8427e", + "test.k_distance.png:md5,d41d8cd98f00b204e9800998ecf8427e", "test.silhouette.png:md5,d41d8cd98f00b204e9800998ecf8427e" ] ] @@ -95,10 +96,54 @@ } } ], - "timestamp": "2026-05-19T14:11:36.600888057", + "timestamp": "2026-09-08T11:43:04.630166242", "meta": { "nf-test": "0.9.5", - "nextflow": "25.09.0" + "nextflow": "26.04.6" + } + }, + "clustermetrics - features from eigenvec": { + "content": [ + [ + [ + { + "id": "test" + }, + "test.metrics.tsv:md5,15c36eab43e480e0311c4bcc3d511477" + ] + ], + [ + [ + { + "id": "test" + }, + "test.k_sweep.csv:md5,98635ea739c5e136ced833916ae5d931" + ] + ], + [ + [ + { + "id": "test" + }, + "test.selected.json:md5,d1d1b3788b7a111f38bf3b27d4bb1ab4" + ] + ], + [ + "versions.yml:md5,602aa5dfe6c0b807d758c4f5cf3fc5e4" + ], + { + "CUSTOM_CLUSTERMETRICS": { + "python": "3.12.13", + "pandas": "3.0.3", + "scikit-learn": "1.8.0", + "matplotlib": "3.10.9" + } + } + ], + "timestamp": "2026-09-08T15:41:30.682155747", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.6" } } } \ No newline at end of file diff --git a/modules/nf-core/custom/clustervisualization/meta.yml b/modules/nf-core/custom/clustervisualization/meta.yml index 62bed07814e5..7f3abf907165 100644 --- a/modules/nf-core/custom/clustervisualization/meta.yml +++ b/modules/nf-core/custom/clustervisualization/meta.yml @@ -29,9 +29,11 @@ input: - 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" + Sample-by-feature matrix: TSV with a `sample_id` column, or a + PLINK2 `.eigenvec` file (`#FID IID PC1 ...` or `#IID PC1 ...`). + Family/source ID columns are dropped; remaining columns are + numeric features (e.g. PCA scores). + pattern: "*.{tsv,txt,eigenvec}" ontologies: - edam: http://edamontology.org/format_3475 - clusters: diff --git a/modules/nf-core/custom/clustervisualization/templates/cluster_viz.py b/modules/nf-core/custom/clustervisualization/templates/cluster_viz.py index 6f6a60d5ac37..a46e4676ee11 100644 --- a/modules/nf-core/custom/clustervisualization/templates/cluster_viz.py +++ b/modules/nf-core/custom/clustervisualization/templates/cluster_viz.py @@ -29,10 +29,24 @@ def load_features(path): - """Read a TSV of `sample_id` + numeric feature columns, indexed by sample_id.""" - df = pd.read_csv(path, sep="\\t") + """Read sample-by-feature matrix, indexed by sample_id. + + Accepts a TSV with a sample_id column, or a PLINK2 .eigenvec file + (`#FID IID PC1 ...` or `#IID PC1 ...`). FID/SID are dropped. + """ + df = pd.read_csv(path, sep="\\t", dtype=str) + df.columns = [str(col).lstrip("#") for col in df.columns] + + if "IID" in df.columns: + df = df.rename(columns={"IID": "sample_id"}) + + drop = [col for col in ("FID", "SID") if col in df.columns] + if drop: + df = df.drop(columns=drop) + 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) @@ -54,9 +68,9 @@ def embed(x, method, umap_neighbors, tsne_perplexity): """ n = len(x) if method == "umap": - reducer = umap.UMAP(n_components=2, n_neighbors=min(umap_neighbors, max(2, n - 1)), random_state=42) + reducer = umap.UMAP(n_components=2, n_neighbors=min(umap_neighbors, max(2, n - 1)), n_jobs=1, random_state=42) elif method == "tsne": - reducer = TSNE(n_components=2, perplexity=min(tsne_perplexity, max(2, n - 1)), random_state=42) + reducer = TSNE(n_components=2, perplexity=min(tsne_perplexity, max(2, n - 1)), n_jobs=1, random_state=42) else: raise ValueError(f"Unknown method '{method}' (expected 'umap' or 'tsne')") return reducer.fit_transform(x) diff --git a/modules/nf-core/custom/clustervisualization/tests/main.nf.test b/modules/nf-core/custom/clustervisualization/tests/main.nf.test index 0334f0d09434..e33116c5ad7d 100644 --- a/modules/nf-core/custom/clustervisualization/tests/main.nf.test +++ b/modules/nf-core/custom/clustervisualization/tests/main.nf.test @@ -1,11 +1,14 @@ nextflow_process { + name "Test Process CUSTOM_CLUSTERVISUALIZATION" script "../main.nf" process "CUSTOM_CLUSTERVISUALIZATION" + tag "modules" tag "modules_nfcore" tag "custom" tag "custom/clustervisualization" + test("clustervisualization - features and clusters") { when { process { @@ -18,18 +21,60 @@ nextflow_process { """ } } + then { + def umap = path(process.out.umap_tsv[0][1]).csv(sep: "\t") + def tsne = path(process.out.tsne_tsv[0][1]).csv(sep: "\t") assertAll( { assert process.success }, + { assert umap.columnNames == ["sample_id", "Dim1", "Dim2", "cluster"] }, + { assert tsne.columnNames == ["sample_id", "Dim1", "Dim2", "cluster"] }, + { assert umap.rowCount == tsne.rowCount }, { assert snapshot( - process.out.umap_tsv, - process.out.tsne_tsv, + umap.columnNames, + umap.rowCount, + tsne.columnNames, + tsne.rowCount, + process.out.versions, + path(process.out.versions[0]).yaml + ).match() } + ) + } + } + + test("clustervisualization - features from eigenvec") { + when { + process { + """ + input[0] = [ + [ id:'test' ], + file(params.modules_testdata_base_path + "genomics/homo_sapiens/popgen/clustering/test.eigenvec", checkIfExists: true), + file(params.modules_testdata_base_path + "genomics/homo_sapiens/popgen/clustering/test_clusters.csv", checkIfExists: true) + ] + """ + } + } + + then { + def umap = path(process.out.umap_tsv[0][1]).csv(sep: "\t") + def tsne = path(process.out.tsne_tsv[0][1]).csv(sep: "\t") + assertAll( + { assert process.success }, + { assert umap.columnNames == ["sample_id", "Dim1", "Dim2", "cluster"] }, + { assert tsne.columnNames == ["sample_id", "Dim1", "Dim2", "cluster"] }, + { assert umap.rowCount == tsne.rowCount }, + { assert snapshot( + umap.columnNames, + umap.rowCount, + tsne.columnNames, + tsne.rowCount, process.out.versions, path(process.out.versions[0]).yaml ).match() } ) } } + test("clustervisualization - features and clusters - stub") { options "-stub" when { @@ -43,6 +88,7 @@ nextflow_process { """ } } + then { assertAll( { assert process.success }, diff --git a/modules/nf-core/custom/clustervisualization/tests/main.nf.test.snap b/modules/nf-core/custom/clustervisualization/tests/main.nf.test.snap index f180dfa8b993..3d113566219a 100644 --- a/modules/nf-core/custom/clustervisualization/tests/main.nf.test.snap +++ b/modules/nf-core/custom/clustervisualization/tests/main.nf.test.snap @@ -54,24 +54,58 @@ "nextflow": "25.09.0" } }, + "clustervisualization - features from eigenvec": { + "content": [ + [ + "sample_id", + "Dim1", + "Dim2", + "cluster" + ], + 200, + [ + "sample_id", + "Dim1", + "Dim2", + "cluster" + ], + 200, + [ + "versions.yml:md5,a00bfbb9b1b4145177ec0e8a7406caf9" + ], + { + "CUSTOM_CLUSTERVISUALIZATION": { + "python": "3.12.13", + "pandas": "3.0.3", + "matplotlib": "3.10.9", + "seaborn": "0.13.2", + "umap-learn": "0.5.12", + "scikit-learn": "1.8.0" + } + } + ], + "timestamp": "2026-09-08T15:44:25.637933809", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.6" + } + }, "clustervisualization - features and clusters": { "content": [ [ - [ - { - "id": "test" - }, - "test.umap.tsv:md5,d338c1ef0e979dbf653e21c3417e975e" - ] + "sample_id", + "Dim1", + "Dim2", + "cluster" ], + 200, [ - [ - { - "id": "test" - }, - "test.tsne.tsv:md5,eb00f27d82530f665552b158e8e3c8ff" - ] + "sample_id", + "Dim1", + "Dim2", + "cluster" ], + 200, [ "versions.yml:md5,a00bfbb9b1b4145177ec0e8a7406caf9" ], @@ -86,10 +120,10 @@ } } ], - "timestamp": "2026-05-19T17:08:32.73077357", + "timestamp": "2026-09-08T12:55:48.907052107", "meta": { "nf-test": "0.9.5", - "nextflow": "25.10.2" + "nextflow": "26.04.6" } } } \ No newline at end of file