diff --git a/.gitignore b/.gitignore index 524f096..9f585e9 100644 --- a/.gitignore +++ b/.gitignore @@ -1,24 +1,17 @@ -# Compiled class file -*.class +.nextflow* +work/ +data/ +results/ +.DS_Store +testing/ +testing* +*.pyc +null/ +\#* -# Log file -*.log +# NF-test +.nf-test/ +.nf-test/* +.nf-test* -# BlueJ files -*.ctxt - -# Mobile Tools for Java (J2ME) -.mtj.tmp/ - -# Package Files # -*.jar -*.war -*.nar -*.ear -*.zip -*.tar.gz -*.rar - -# virtual machine crash logs, see http://www.java.com/en/download/help/error_hotspot.xml -hs_err_pid* -replay_pid* +modules/nf-core diff --git a/modules/UMCUGenetics/ancestry_knn/calc/main.nf b/modules/UMCUGenetics/ancestry_knn/calc/main.nf new file mode 100644 index 0000000..4701eb3 --- /dev/null +++ b/modules/UMCUGenetics/ancestry_knn/calc/main.nf @@ -0,0 +1,34 @@ +process ANCESTRY_KNN { + tag "${meta.id}" + label "process_medium" + + container "ghcr.io/astral-sh/uv:python3.13-bookworm" + + input: + tuple val(meta), path(eigenvec) + tuple val(meta2), path(ref_metadata) + + output: + tuple val(meta), path("*_knn.tsv"), emit: knn_tsv + tuple val(meta), path("*_knn_pca.png"), emit: knn_pca_plot, optional: true + tuple val("${task.process}"), val('ancestry_knn'), eval('echo 1.0.0'), emit: versions_ancestry_knn, topic: versions + + script: + def prefix = task.ext.prefix ?: meta.id + def args = task.ext.args ?: "" + """ + knn.py \ + --eig ${eigenvec} \\ + --labels ${ref_metadata} \\ + ${args} \\ + --plot-output ${prefix}_knn_pca.png \\ + --output ${prefix}_knn.tsv + """ + + stub: + def prefix = task.ext.prefix ?: meta.id + """ + touch ${prefix}_knn.tsv + touch ${prefix}_knn_pca.png + """ +} diff --git a/modules/UMCUGenetics/ancestry_knn/calc/resources/usr/bin/knn.py b/modules/UMCUGenetics/ancestry_knn/calc/resources/usr/bin/knn.py new file mode 100755 index 0000000..b33400c --- /dev/null +++ b/modules/UMCUGenetics/ancestry_knn/calc/resources/usr/bin/knn.py @@ -0,0 +1,424 @@ +#!/usr/bin/env -S uv run --script --no-cache +# /// script +# requires-python = ">=3.12" +# dependencies = [ +# "matplotlib~=3.10.0", +# "numpy==2.4.6", +# "pandas==3.0.3", +# "typer==0.26.7", +# ] +# /// +"""Run K-nearest-neighbor ancestry prediction from PCA coordinates.""" + +from __future__ import annotations + +from pathlib import Path +from typing import Annotated + +import matplotlib +import numpy as np +import pandas as pd +import typer + +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +app = typer.Typer(add_completion=False, help="KNN ancestry prediction from PCA coordinates.") + + +def read_input_tables( + eig_path: Path, + labels_path: Path, + sep: str, + id_col: str, + label_col: str, +) -> pd.DataFrame: + """Read PCA and label tables, validate required columns, and merge by ID. + + Args: + eig_path: Path to the PCA eigenvector file. + labels_path: Path to the sample label file. + sep: Delimiter used in both input files. + id_col: Column name for sample identifiers. + label_col: Column name for group labels. + + Returns: + Merged DataFrame with PCA features and labels joined on ``id_col``. + + Raises: + typer.BadParameter: If required columns are missing from either file. + """ + eig = pd.read_csv(eig_path, sep=sep) + labels = pd.read_csv(labels_path, sep=sep) + + required_eig_cols = {id_col} + required_label_cols = {id_col, label_col} + missing_eig = required_eig_cols - set(eig.columns) + missing_lab = required_label_cols - set(labels.columns) + + if missing_eig: + raise typer.BadParameter( + f"Missing required column(s) in eig file: {', '.join(sorted(missing_eig))}" + ) + if missing_lab: + raise typer.BadParameter( + f"Missing required column(s) in labels file: {', '.join(sorted(missing_lab))}" + ) + + return eig.merge(labels[[id_col, label_col]], on=id_col, how="left") + + +def get_pc_columns(df: pd.DataFrame) -> list[str]: + """Return PCA feature columns named like PC1, PC2, etc. + + Args: + df: DataFrame expected to contain one or more PC columns. + + Returns: + List of column names whose names start with ``PC`` (case-insensitive). + + Raises: + typer.BadParameter: If no PC columns are found in ``df``. + """ + pc_cols = [column for column in df.columns if column.upper().startswith("PC")] + if not pc_cols: + raise typer.BadParameter( + "No PC columns found. Expected columns like PC1, PC2, ..." + ) + return pc_cols + + +def sort_pc_columns(pc_cols: list[str]) -> list[str]: + """Sort PC columns by their numeric suffix when present (PC1, PC2, PC10, ...). + + Columns with a non-numeric suffix are sorted after numeric ones. + + Args: + pc_cols: List of PC column names to sort. + + Returns: + Sorted list of PC column names. + """ + def _pc_sort_key(column: str) -> tuple[int, str]: + suffix = column[2:] + return (int(suffix), column) if suffix.isdigit() else (10**9, column) + + return sorted(pc_cols, key=_pc_sort_key) + + +def split_train_predict( + merged_df: pd.DataFrame, + label_col: str, +) -> tuple[pd.DataFrame, pd.DataFrame]: + """Split merged data into labeled training rows and unlabeled prediction rows. + + Args: + merged_df: DataFrame containing both labeled reference and unlabeled samples. + label_col: Column name that distinguishes labeled (non-NaN) from unlabeled rows. + + Returns: + A ``(train, predict)`` tuple where ``train`` contains rows with a known label + and ``predict`` contains rows whose label is NaN. + + Raises: + typer.BadParameter: If there are no labeled or no unlabeled samples after splitting. + """ + train = merged_df[merged_df[label_col].notna()].copy() + predict = merged_df[merged_df[label_col].isna()].copy() + + if train.empty: + raise typer.BadParameter( + "No labeled samples found after merge. Check if IDs match between files." + ) + if predict.empty: + raise typer.BadParameter( + "No unlabeled samples found. Nothing to predict." + ) + + return train, predict + + +def standardize_train_query( + train_features: np.ndarray, + query_features: np.ndarray, +) -> tuple[np.ndarray, np.ndarray]: + """Standardize features using training-set mean and standard deviation. + + Columns with zero variance in the training set are left unchanged (std set to 1). + + Args: + train_features: Feature matrix for labeled training samples ``(n_train, n_pcs)``. + query_features: Feature matrix for samples to predict ``(n_query, n_pcs)``. + + Returns: + A ``(train_std, query_std)`` tuple of standardized feature matrices. + """ + mean = train_features.mean(axis=0) + std = train_features.std(axis=0) + std[std == 0] = 1.0 + return (train_features - mean) / std, (query_features - mean) / std + + +def predict_knn( + train_features: np.ndarray, + train_labels: np.ndarray, + query_features: np.ndarray, + k: int, +) -> tuple[np.ndarray, np.ndarray]: + """Predict class labels and confidence scores using a simple KNN vote. + + Confidence is the fraction of the *k* nearest neighbors that agree on the + winning label. When ``k`` exceeds the number of training samples, the + effective neighbourhood is capped at ``len(train_labels)``. + + Args: + train_features: Feature matrix for labeled training samples ``(n_train, n_pcs)``. + train_labels: Label array aligned with ``train_features`` ``(n_train,)``. + query_features: Feature matrix for samples to predict ``(n_query, n_pcs)``. + k: Number of nearest neighbors to consider. + + Returns: + A ``(predicted_labels, confidences)`` tuple, both of shape ``(n_query,)``. + + Raises: + typer.BadParameter: If ``k`` is less than 1. + """ + if k < 1: + raise typer.BadParameter("k must be >= 1.") + + effective_k = min(k, len(train_labels)) + predicted_labels: list[str] = [] + confidences: list[float] = [] + + for query_row in query_features: + distances = np.sum((train_features - query_row) ** 2, axis=1) + neighbor_idx = np.argpartition(distances, effective_k - 1)[:effective_k] + neighbor_labels = train_labels[neighbor_idx] + + labels, counts = np.unique(neighbor_labels, return_counts=True) + winner_index = int(np.argmax(counts)) + predicted_labels.append(str(labels[winner_index])) + confidences.append(float(counts[winner_index] / effective_k)) + + return np.array(predicted_labels), np.array(confidences) + + +def write_predictions( + sample_ids: np.ndarray, + predicted_labels: np.ndarray, + confidences: np.ndarray, + output_path: Path, +) -> None: + """Write KNN predictions to a tab-separated output file. + + Args: + sample_ids: Array of sample identifiers ``(n_query,)``. + predicted_labels: Predicted group label for each sample ``(n_query,)``. + confidences: KNN confidence score for each prediction ``(n_query,)``. + output_path: Destination path for the output TSV. + """ + out_df = pd.DataFrame( + { + "#IID": sample_ids, + "pred_group": predicted_labels, + "knn_conf": confidences, + } + ) + out_df.to_csv(output_path, sep="\t", index=False) + + +def write_pca_plot( + train_df: pd.DataFrame, + predict_df: pd.DataFrame, + predicted_labels: np.ndarray, + pc_x: str, + pc_y: str, + label_col: str, + output_path: Path, +) -> None: + """Generate a PCA scatter plot colored by cluster/group labels. + + Reference samples are drawn as small filled circles; predicted samples are + drawn as larger ``X`` markers with a black edge, using the same color as + their assigned cluster. + + Args: + train_df: DataFrame of labeled reference samples including ``pc_x``, ``pc_y``, + and ``label_col`` columns. + predict_df: DataFrame of samples to predict including ``pc_x`` and ``pc_y`` columns. + predicted_labels: Predicted group label for each row in ``predict_df`` ``(n_query,)``. + pc_x: Column name to use as the x-axis. + pc_y: Column name to use as the y-axis. + label_col: Column name holding the reference group label in ``train_df``. + output_path: Destination path for the output PNG. + """ + plot_train = train_df[[pc_x, pc_y, label_col]].copy() + plot_train["cluster"] = plot_train[label_col].astype(str) + plot_train["source"] = "reference" + + plot_predict = predict_df[[pc_x, pc_y]].copy() + plot_predict["cluster"] = predicted_labels + plot_predict["source"] = "predicted" + + plot_df = pd.concat( + [ + plot_train[[pc_x, pc_y, "cluster", "source"]], + plot_predict[[pc_x, pc_y, "cluster", "source"]], + ], + ignore_index=True, + ) + + clusters = sorted(plot_df["cluster"].unique()) + cmap = plt.cm.get_cmap("tab20", max(len(clusters), 1)) + color_map = {cluster: cmap(i) for i, cluster in enumerate(clusters)} + + fig, ax = plt.subplots(figsize=(8, 6)) + + for cluster in clusters: + cluster_train = plot_df[(plot_df["cluster"] == cluster) & (plot_df["source"] == "reference")] + cluster_pred = plot_df[(plot_df["cluster"] == cluster) & (plot_df["source"] == "predicted")] + color = color_map[cluster] + + if not cluster_train.empty: + ax.scatter( + cluster_train[pc_x], + cluster_train[pc_y], + s=28, + alpha=0.75, + color=color, + edgecolors="none", + label=str(cluster), + ) + + if not cluster_pred.empty: + ax.scatter( + cluster_pred[pc_x], + cluster_pred[pc_y], + s=80, + alpha=1.0, + marker="X", + color=color, + edgecolors="black", + linewidths=0.6, + ) + + ax.set_xlabel(pc_x) + ax.set_ylabel(pc_y) + ax.set_title("PCA Clusters with KNN Predictions") + + handles, labels = ax.get_legend_handles_labels() + if handles: + by_label = dict(zip(labels, handles)) + ax.legend(by_label.values(), by_label.keys(), title="Cluster", bbox_to_anchor=(1.02, 1), loc="upper left") + + fig.tight_layout() + fig.savefig(output_path, dpi=200) + plt.close(fig) + + +@app.command() +def main( + eig: Annotated[ + Path, + typer.Option( + "--eig", + "-e", + exists=True, + dir_okay=False, + readable=True, + help="PCA table with ID and PC columns.", + ), + ], + labels: Annotated[ + Path, + typer.Option( + "--labels", + "-l", + exists=True, + dir_okay=False, + readable=True, + help="Label table with ID and group label columns.", + ), + ], + output: Annotated[ + Path, + typer.Option("--output", "-o", help="Output TSV path."), + ] = Path("knn_pred.tsv"), + k: Annotated[int, typer.Option("--k", "-k", min=1, help="Number of neighbors.")] = 5, + id_col: Annotated[str, typer.Option("--id-col", help="Sample ID column name.")] = "#IID", + label_col: Annotated[str, typer.Option("--label-col", help="Label column name.")] = "SuperPop", + sep: Annotated[str, typer.Option("--sep", help="Input delimiter used by both files.")] = "\t", + normalize: Annotated[ + bool, typer.Option("--normalize/--no-normalize", help="Standardize PC columns.") + ] = True, + conf_threshold: Annotated[float, typer.Option('--conf_threshold', help="Confidence threshold for plot outputting")] = 0.6, + req_label: Annotated[str, typer.Option("--required_superpop", help="Expected super population for ancestry flagging")] = "EUR", + plot_output: Annotated[ + Path, typer.Option("--plot-output", help="Output PNG path for PCA scatter plot.") + ] = Path("knn_pca.png"), +) -> None: + """Run the KNN workflow and write prediction output. + + Reads PCA and label files, trains a KNN classifier on labeled reference + samples, predicts ancestry for unlabeled samples, and writes a TSV of + predictions. A PCA scatter plot is written when the first sample's + confidence falls at or below ``conf_threshold`` or its predicted label + differs from ``req_label``. + """ + merged_df = read_input_tables( + eig_path=eig, + labels_path=labels, + sep=sep, + id_col=id_col, + label_col=label_col, + ) + pc_cols = sort_pc_columns(get_pc_columns(merged_df)) + if len(pc_cols) < 2: + raise typer.BadParameter( + "At least two PC columns are required." + ) + train_df, predict_df = split_train_predict(merged_df, label_col=label_col) + + train_features = train_df[pc_cols].to_numpy(dtype=float) + query_features = predict_df[pc_cols].to_numpy(dtype=float) + train_labels = train_df[label_col].to_numpy(dtype=str) + + if normalize: + train_features, query_features = standardize_train_query( + train_features=train_features, + query_features=query_features, + ) + + predicted_labels, confidences = predict_knn( + train_features=train_features, + train_labels=train_labels, + query_features=query_features, + k=k, + ) + + write_predictions( + sample_ids=predict_df[id_col].to_numpy(), + predicted_labels=predicted_labels, + confidences=confidences, + output_path=output, + ) + + first_confidence = confidences[0] + first_predicted_label = predicted_labels[0] + + if first_confidence <= conf_threshold or first_predicted_label != req_label: + write_pca_plot( + train_df=train_df, + predict_df=predict_df, + predicted_labels=predicted_labels, + pc_x=pc_cols[0], + pc_y=pc_cols[1], + label_col=label_col, + output_path=plot_output, + ) + typer.echo(f"Wrote {plot_output}") + typer.echo(f"Wrote {output}") + + +if __name__ == "__main__": + app() diff --git a/modules/UMCUGenetics/ancestry_knn/calc/resources/usr/bin/test_knn.py b/modules/UMCUGenetics/ancestry_knn/calc/resources/usr/bin/test_knn.py new file mode 100644 index 0000000..a1448cd --- /dev/null +++ b/modules/UMCUGenetics/ancestry_knn/calc/resources/usr/bin/test_knn.py @@ -0,0 +1,152 @@ +"""Unit tests for knn.py helper functions.""" + +import numpy as np +import pandas as pd +import pytest +import typer + +from knn import ( + get_pc_columns, + predict_knn, + sort_pc_columns, + split_train_predict, + standardize_train_query, +) + + +# --------------------------------------------------------------------------- +# sort_pc_columns +# --------------------------------------------------------------------------- + +def test_sort_pc_columns_numeric_order(): + assert sort_pc_columns(["PC10", "PC2", "PC1"]) == ["PC1", "PC2", "PC10"] + + +def test_sort_pc_columns_non_numeric_suffix_last(): + result = sort_pc_columns(["PCx", "PC2", "PC1"]) + assert result[:2] == ["PC1", "PC2"] + assert result[-1] == "PCx" + + +def test_sort_pc_columns_single_element(): + assert sort_pc_columns(["PC1"]) == ["PC1"] + + +# --------------------------------------------------------------------------- +# get_pc_columns +# --------------------------------------------------------------------------- + +def test_get_pc_columns_returns_matching_columns(): + df = pd.DataFrame(columns=["#IID", "PC1", "PC2", "SuperPop"]) + assert set(get_pc_columns(df)) == {"PC1", "PC2"} + + +def test_get_pc_columns_case_insensitive(): + df = pd.DataFrame(columns=["#IID", "pc1", "PC2"]) + assert set(get_pc_columns(df)) == {"pc1", "PC2"} + + +def test_get_pc_columns_raises_when_none_found(): + df = pd.DataFrame(columns=["#IID", "SuperPop"]) + with pytest.raises(typer.BadParameter): + get_pc_columns(df) + + +# --------------------------------------------------------------------------- +# split_train_predict +# --------------------------------------------------------------------------- + +def _make_merged_df(labeled_n: int, unlabeled_n: int) -> pd.DataFrame: + rows = [] + for i in range(labeled_n): + rows.append({"#IID": f"ref_{i}", "PC1": float(i), "SuperPop": "EUR"}) + for i in range(unlabeled_n): + rows.append({"#IID": f"sample_{i}", "PC1": float(i + 100), "SuperPop": np.nan}) + return pd.DataFrame(rows) + + +def test_split_train_predict_sizes(): + df = _make_merged_df(labeled_n=10, unlabeled_n=3) + train, predict = split_train_predict(df, label_col="SuperPop") + assert len(train) == 10 + assert len(predict) == 3 + + +def test_split_train_predict_raises_no_train(): + df = _make_merged_df(labeled_n=0, unlabeled_n=3) + with pytest.raises(typer.BadParameter, match="No labeled"): + split_train_predict(df, label_col="SuperPop") + + +def test_split_train_predict_raises_no_predict(): + df = _make_merged_df(labeled_n=5, unlabeled_n=0) + with pytest.raises(typer.BadParameter, match="Nothing to predict"): + split_train_predict(df, label_col="SuperPop") + + +# --------------------------------------------------------------------------- +# standardize_train_query +# --------------------------------------------------------------------------- + +def test_standardize_zero_mean_unit_std(): + train = np.array([[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]]) + query = np.array([[3.0, 4.0]]) + train_std, query_std = standardize_train_query(train, query) + np.testing.assert_allclose(train_std.mean(axis=0), 0.0, atol=1e-10) + np.testing.assert_allclose(train_std.std(axis=0), 1.0, atol=1e-10) + # query center-of-training → all zeros after standardization + np.testing.assert_allclose(query_std, [[0.0, 0.0]], atol=1e-10) + + +def test_standardize_zero_variance_column_unchanged(): + """Columns with zero variance must not cause division-by-zero.""" + train = np.array([[1.0, 5.0], [1.0, 7.0]]) # column 0 has std=0 + query = np.array([[1.0, 6.0]]) + train_std, _ = standardize_train_query(train, query) + # zero-variance column is left unchanged (divided by 1 after clamp) + np.testing.assert_allclose(train_std[:, 0], 0.0, atol=1e-10) + + +# --------------------------------------------------------------------------- +# predict_knn +# --------------------------------------------------------------------------- + +def _simple_train(): + """Two clusters far apart: A around 0, B around 10.""" + train_features = np.array( + [[0.0], [0.1], [0.2], [10.0], [10.1], [10.2]] + ) + train_labels = np.array(["A", "A", "A", "B", "B", "B"]) + return train_features, train_labels + + +def test_predict_knn_clear_majority(): + train_features, train_labels = _simple_train() + query = np.array([[0.05], [10.05]]) + labels, confs = predict_knn(train_features, train_labels, query, k=3) + assert list(labels) == ["A", "B"] + np.testing.assert_allclose(confs, [1.0, 1.0]) + + +def test_predict_knn_confidence_range(): + train_features, train_labels = _simple_train() + query = np.array([[5.0]]) # equidistant — confidence will be < 1 + labels, confs = predict_knn(train_features, train_labels, query, k=6) + assert 0.0 <= confs[0] <= 1.0 + + +def test_predict_knn_k_exceeds_training_size(): + """k larger than training set should not raise, just use all neighbors.""" + train_features = np.array([[0.0], [1.0]]) + train_labels = np.array(["A", "A"]) + query = np.array([[0.5]]) + labels, confs = predict_knn(train_features, train_labels, query, k=100) + assert labels[0] == "A" + assert confs[0] == 1.0 + + +def test_predict_knn_k_less_than_1_raises(): + train_features = np.array([[0.0]]) + train_labels = np.array(["A"]) + with pytest.raises(typer.BadParameter, match="k must be"): + predict_knn(train_features, train_labels, np.array([[0.0]]), k=0) diff --git a/modules/UMCUGenetics/ancestry_knn/calc/tests/main.nf.test b/modules/UMCUGenetics/ancestry_knn/calc/tests/main.nf.test new file mode 100644 index 0000000..786e682 --- /dev/null +++ b/modules/UMCUGenetics/ancestry_knn/calc/tests/main.nf.test @@ -0,0 +1,50 @@ +nextflow_process { + name "Test ancestry_knn/calc" + script "../main.nf" + process "ANCESTRY_KNN" + + tag "modules/UMCUGenetics" + tag "modules" + tag "ancestry" + tag "prs" + + test("Test Ancestry KNN calculations") { + when { + process { + """ + input[0] = [[id: "eigenvec"], + file(params.modules_umcu_testdata_base_path + "ancestry_knn/pca_eig.tsv", + checkIfExists: true) + ] + input[1] = [[id: 'labels'], + file(params.modules_umcu_testdata_base_path + "ancestry_knn/labels.tsv", + checkIfExists: true) + ] + """ + } + + } + then { + assertAll( + { assert process.success }, + + { assert path(process.out.knn_tsv[0][1]).readLines().size() == 5 }, + { assert path(process.out.knn_tsv[0][1]).readLines()[0] + == "#IID\tpred_group\tknn_conf" }, + + { assert path(process.out.knn_tsv[0][1]).text.contains("QRY_001\tAFR\t0.6") }, + { assert path(process.out.knn_tsv[0][1]).text.contains("QRY_003\tEAS\t1.0") }, + + { assert path(process.out.knn_pca_plot[0][1]).exists() }, + { assert path(process.out.knn_pca_plot[0][1]).size() > 0 }, + { assert process.out.knn_pca_plot[0][1].endsWith(".png") }, + + { assert snapshot( + process.out.knn_tsv, + process.out.findAll { key, val -> key.startsWith('versions') } + ).match() + } + ) + } + } +} diff --git a/modules/UMCUGenetics/ancestry_knn/calc/tests/main.nf.test.snap b/modules/UMCUGenetics/ancestry_knn/calc/tests/main.nf.test.snap new file mode 100644 index 0000000..792a96d --- /dev/null +++ b/modules/UMCUGenetics/ancestry_knn/calc/tests/main.nf.test.snap @@ -0,0 +1,28 @@ +{ + "Test Ancestry KNN calculations": { + "content": [ + [ + [ + { + "id": "eigenvec" + }, + "eigenvec_knn.tsv:md5,bab411a7128384b8a5db6729cf144888" + ] + ], + { + "versions_ancestry_knn": [ + [ + "ANCESTRY_KNN", + "ancestry_knn", + "1.0.0" + ] + ] + } + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "timestamp": "2026-07-13T13:25:42.562506" + } +} \ No newline at end of file diff --git a/modules/UMCUGenetics/ancestry_knn/merge/main.nf b/modules/UMCUGenetics/ancestry_knn/merge/main.nf new file mode 100644 index 0000000..29c64fb --- /dev/null +++ b/modules/UMCUGenetics/ancestry_knn/merge/main.nf @@ -0,0 +1,20 @@ +process ANCESTRY_KNN_MERGE { + tag "ANCESTRY_KNN_MERGE" + label "process_low" + + input: + path(knn_tsvs) + + output: + path("ancestry_knn_mqc.tsv"), emit: knn_mqc_tsv + tuple val("${task.process}"), val('ancestry_knn'), eval('echo 1.0.0'), emit: versions_ancestry_knn, topic: versions + + script: + """ + echo "Sample ID\tPrediction Group\tConfidence" > ancestry_knn_mqc.tsv + + for f in ${knn_tsvs}; do + tail -n 1 \$f >> ancestry_knn_mqc.tsv + done + """ +} diff --git a/modules/UMCUGenetics/ancestry_knn/merge/tests/main.nf.test b/modules/UMCUGenetics/ancestry_knn/merge/tests/main.nf.test new file mode 100644 index 0000000..916479c --- /dev/null +++ b/modules/UMCUGenetics/ancestry_knn/merge/tests/main.nf.test @@ -0,0 +1,34 @@ +nextflow_process { + name "Test ancestry_knn/merge" + script "../main.nf" + process "ANCESTRY_KNN_MERGE" + + tag "modules/UMCUGenetics" + tag "modules" + tag "ancestry" + tag "prs" + + test("Test ancestry KNN merging") { + when { + process { + """ + input[0] = [ + file(params.modules_umcu_testdata_base_path + "ancestry_knn/ancestry_result.tsv", + checkIfExists: true), + file(params.modules_umcu_testdata_base_path + "ancestry_knn/ancestry_result2.tsv", + checkIfExists: true) + ] + """ + } + } + then { + assertAll( + { assert process.success }, + { assert snapshot( + process.out.knn_mqc_tsv, + process.out.findAll { key, val -> key.startsWith('versions') } + ).match() } + ) + } + } +} diff --git a/modules/UMCUGenetics/ancestry_knn/merge/tests/main.nf.test.snap b/modules/UMCUGenetics/ancestry_knn/merge/tests/main.nf.test.snap new file mode 100644 index 0000000..6dd8573 --- /dev/null +++ b/modules/UMCUGenetics/ancestry_knn/merge/tests/main.nf.test.snap @@ -0,0 +1,23 @@ +{ + "Test ancestry KNN merging": { + "content": [ + [ + "ancestry_knn_mqc.tsv:md5,d79057c19aada0ce7f002428174f16d9" + ], + { + "versions_ancestry_knn": [ + [ + "ANCESTRY_KNN_MERGE", + "ancestry_knn", + "1.0.0" + ] + ] + } + ], + "meta": { + "nf-test": "0.9.3", + "nextflow": "25.10.4" + }, + "timestamp": "2026-07-13T13:25:46.370844" + } +} \ No newline at end of file diff --git a/modules/UMCUGenetics/gatk4/haplotypecaller_alleles/environment.yml b/modules/UMCUGenetics/gatk4/haplotypecaller_alleles/environment.yml new file mode 100644 index 0000000..67e0eb8 --- /dev/null +++ b/modules/UMCUGenetics/gatk4/haplotypecaller_alleles/environment.yml @@ -0,0 +1,10 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/environment-schema.json +channels: + - conda-forge + - bioconda +dependencies: + # renovate: datasource=conda depName=bioconda/gatk4 + - bioconda::gatk4=4.6.2.0 + # renovate: datasource=conda depName=bioconda/gcnvkernel + - bioconda::gcnvkernel=0.9 diff --git a/modules/UMCUGenetics/gatk4/haplotypecaller_alleles/main.nf b/modules/UMCUGenetics/gatk4/haplotypecaller_alleles/main.nf new file mode 100644 index 0000000..22723ee --- /dev/null +++ b/modules/UMCUGenetics/gatk4/haplotypecaller_alleles/main.nf @@ -0,0 +1,82 @@ +process GATK4_HAPLOTYPECALLER { + tag "${meta.id}" + label 'process_low' + + conda "${moduleDir}/environment.yml" + container "${workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container + ? 'https://community-cr-prod.seqera.io/docker/registry/v2/blobs/sha256/b2/b28daf5d9bb2f0d129dcad1b7410e0dd8a9b087aaf3ec7ced929b1f57624ad98/data' + : 'community.wave.seqera.io/library/gatk4_gcnvkernel:e48d414933d188cd'}" + + input: + tuple val(meta), path(input), path(input_index), path(intervals), path(dragstr_model) + tuple val(meta2), path(fasta) + tuple val(meta3), path(fai) + tuple val(meta4), path(dict) + tuple val(meta5), path(dbsnp) + tuple val(meta6), path(dbsnp_tbi) + tuple val(meta7), path(alleles), path(alleles_tbi) + + + output: + tuple val(meta), path("*.vcf.gz"), emit: vcf + tuple val(meta), path("*.tbi"), optional: true, emit: tbi + tuple val(meta), path("*.realigned.bam"), optional: true, emit: bam + path "versions.yml", emit: versions + + when: + task.ext.when == null || task.ext.when + + script: + def args = task.ext.args ?: '' + def prefix = task.ext.prefix ?: "${meta.id}" + def dbsnp_command = dbsnp ? "--dbsnp ${dbsnp}" : "" + def interval_command = intervals ? "--intervals ${intervals}" : "" + def dragstr_command = dragstr_model ? "--dragstr-params-path ${dragstr_model}" : "" + def bamout_command = args.contains("--bam-writer-type") ? "--bam-output ${prefix.replaceAll('.g\\s*$', '')}.realigned.bam" : "" + def alleles_command = alleles ? "--alleles ${alleles}" : "" + + def avail_mem = 3072 + if (!task.memory) { + log.info('[GATK HaplotypeCaller] Available memory not known - defaulting to 3GB. Specify process memory requirements to change this.') + } + else { + avail_mem = (task.memory.mega * 0.8).intValue() + } + """ + gatk --java-options "-Xmx${avail_mem}M -XX:-UsePerfData" \\ + HaplotypeCaller \\ + --input ${input} \\ + --output ${prefix}.vcf.gz \\ + --reference ${fasta} \\ + --native-pair-hmm-threads ${task.cpus} \\ + ${dbsnp_command} \\ + ${interval_command} \\ + ${dragstr_command} \\ + ${bamout_command} \\ + ${alleles_command} \\ + --tmp-dir . \\ + ${args} + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + gatk4: \$(echo \$(gatk --version 2>&1) | sed 's/^.*(GATK) v//; s/ .*\$//') + END_VERSIONS + """ + + stub: + def args = task.ext.args ?: '' + def prefix = task.ext.prefix ?: "${meta.id}" + def bamout_command = args.contains("--bam-writer-type") ? "--bam-output ${prefix.replaceAll('.g\\s*$', '')}.realigned.bam" : "" + + def stub_realigned_bam = bamout_command ? "touch ${prefix.replaceAll('.g\\s*$', '')}.realigned.bam" : "" + """ + touch ${prefix}.vcf.gz + touch ${prefix}.vcf.gz.tbi + ${stub_realigned_bam} + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + gatk4: \$(echo \$(gatk --version 2>&1) | sed 's/^.*(GATK) v//; s/ .*\$//') + END_VERSIONS + """ +} diff --git a/modules/UMCUGenetics/gatk4/haplotypecaller_alleles/meta.yml b/modules/UMCUGenetics/gatk4/haplotypecaller_alleles/meta.yml new file mode 100644 index 0000000..5f9dbd5 --- /dev/null +++ b/modules/UMCUGenetics/gatk4/haplotypecaller_alleles/meta.yml @@ -0,0 +1,139 @@ +name: gatk4_haplotypecaller +description: Call germline SNPs and indels via local re-assembly of haplotypes +keywords: + - gatk4 + - haplotype + - haplotypecaller +tools: + - gatk4: + description: | + Developed in the Data Sciences Platform at the Broad Institute, the toolkit offers a wide variety of tools + with a primary focus on variant discovery and genotyping. Its powerful processing engine + and high-performance computing features make it capable of taking on projects of any size. + homepage: https://gatk.broadinstitute.org/hc/en-us + documentation: https://gatk.broadinstitute.org/hc/en-us/categories/360002369672s + doi: 10.1158/1538-7445.AM2017-3590 + licence: ["Apache-2.0"] + identifier: "" +input: + - - meta: + type: map + description: | + Groovy Map containing sample information + e.g. [ id:'test', single_end:false ] + - input: + type: file + description: BAM/CRAM file from alignment + pattern: "*.{bam,cram}" + ontologies: [] + - input_index: + type: file + description: BAI/CRAI file from alignment + pattern: "*.{bai,crai}" + ontologies: [] + - intervals: + type: file + description: Bed file with the genomic regions included in the library (optional) + ontologies: [] + - dragstr_model: + type: file + description: Text file containing the DragSTR model of the used BAM/CRAM file + (optional) + pattern: "*.txt" + ontologies: [] + - - meta2: + type: map + description: | + Groovy Map containing reference information + e.g. [ id:'test_reference' ] + - fasta: + type: file + description: The reference fasta file + pattern: "*.fasta" + ontologies: [] + - - meta3: + type: map + description: | + Groovy Map containing reference information + e.g. [ id:'test_reference' ] + - fai: + type: file + description: Index of reference fasta file + pattern: "fasta.fai" + ontologies: [] + - - meta4: + type: map + description: | + Groovy Map containing reference information + e.g. [ id:'test_reference' ] + - dict: + type: file + description: GATK sequence dictionary + pattern: "*.dict" + ontologies: [] + - - meta5: + type: map + description: | + Groovy Map containing dbsnp information + e.g. [ id:'test_dbsnp' ] + - dbsnp: + type: file + description: VCF file containing known sites (optional) + ontologies: [] + - - meta6: + type: map + description: | + Groovy Map containing dbsnp information + e.g. [ id:'test_dbsnp' ] + - dbsnp_tbi: + type: file + description: VCF index of dbsnp (optional) + ontologies: [] +output: + vcf: + - - meta: + type: map + description: | + Groovy Map containing sample information + e.g. [ id:'test', single_end:false ] + - "*.vcf.gz": + type: file + description: Compressed VCF file + pattern: "*.vcf.gz" + ontologies: + - edam: http://edamontology.org/format_3989 # GZIP format + tbi: + - - meta: + type: map + description: | + Groovy Map containing sample information + e.g. [ id:'test', single_end:false ] + - "*.tbi": + type: file + description: Index of VCF file + pattern: "*.vcf.gz.tbi" + ontologies: [] + bam: + - - meta: + type: map + description: | + Groovy Map containing sample information + e.g. [ id:'test', single_end:false ] + - "*.realigned.bam": + type: file + description: Assembled haplotypes and locally realigned reads + pattern: "*.realigned.bam" + ontologies: [] + versions: + - versions.yml: + type: file + description: File containing software versions + pattern: "versions.yml" + ontologies: + - edam: http://edamontology.org/format_3750 # YAML +authors: + - "@suzannejin" + - "@FriederikeHanssen" +maintainers: + - "@suzannejin" + - "@FriederikeHanssen" diff --git a/modules/UMCUGenetics/gatk4/haplotypecaller_alleles/tests/main.nf.test b/modules/UMCUGenetics/gatk4/haplotypecaller_alleles/tests/main.nf.test new file mode 100644 index 0000000..18d35f4 --- /dev/null +++ b/modules/UMCUGenetics/gatk4/haplotypecaller_alleles/tests/main.nf.test @@ -0,0 +1,154 @@ +// nf-core modules test gatk4/haplotypecaller +nextflow_process { + + name "Test Process GATK4_HAPLOTYPECALLER" + script "../main.nf" + process "GATK4_HAPLOTYPECALLER" + + tag "modules" + tag "modules_nfcore" + tag "gatk4" + tag "gatk4/haplotypecaller" + + test("homo_sapiens - [bam, bai] - fasta - fai - dict") { + + when { + process { + """ + input[0] = [ + [ id:'test_bam' ], // meta map + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/bam/test.paired_end.sorted.bam', checkIfExists: true), + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/bam/test.paired_end.sorted.bam.bai', checkIfExists: true), + [], + [] + ] + input[1] = [ [ id:'test_fa' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta', checkIfExists: true) ] + input[2] = [ [ id:'test_fai' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta.fai', checkIfExists: true) ] + input[3] = [ [ id:'test_dict' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.dict', checkIfExists: true) ] + input[4] = [ [], [] ] + input[5] = [ [], [] ] + + """ + } + } + + then { + assertAll( + { assert process.success }, + { assert snapshot( + file(process.out.vcf[0][1]).name, + file(process.out.tbi[0][1]).name, + process.out.versions + ).match() + } + ) + } + + } + + test("homo_sapiens - [cram, crai] - fasta - fai - dict") { + + when { + process { + """ + input[0] = [ + [ id:'test_cram' ], // meta map + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/cram/test.paired_end.sorted.cram', checkIfExists: true), + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/cram/test.paired_end.sorted.cram.crai', checkIfExists: true), + [], + [] + ] + input[1] = [ [ id:'test_fa' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta', checkIfExists: true) ] + input[2] = [ [ id:'test_fai' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta.fai', checkIfExists: true) ] + input[3] = [ [ id:'test_dict' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.dict', checkIfExists: true) ] + input[4] = [ [], [] ] + input[5] = [ [], [] ] + """ + } + } + + then { + assertAll( + { assert process.success }, + { assert snapshot( + file(process.out.vcf[0][1]).name, + file(process.out.tbi[0][1]).name, + process.out.versions + ).match() + } + ) + } + + } + + test("homo_sapiens - [cram, crai] - fasta - fai - dict - sites - sites_tbi") { + + when { + process { + """ + input[0] = [ + [ id:'test_cram_sites' ], // meta map + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/cram/test.paired_end.sorted.cram', checkIfExists: true), + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/cram/test.paired_end.sorted.cram.crai', checkIfExists: true), + [], + [] + ] + input[1] = [ [ id:'test_fa' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta', checkIfExists: true) ] + input[2] = [ [ id:'test_fai' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta.fai', checkIfExists: true) ] + input[3] = [ [ id:'test_dict' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.dict', checkIfExists: true) ] + input[4] = [ [ id:'test_sites' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/vcf/dbsnp_146.hg38.vcf.gz', checkIfExists: true) ] + input[5] = [ [ id:'test_sites_tbi' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/vcf/dbsnp_146.hg38.vcf.gz.tbi', checkIfExists: true) ] + """ + } + } + + then { + assertAll( + { assert process.success }, + { assert snapshot( + file(process.out.vcf[0][1]).name, + file(process.out.tbi[0][1]).name, + process.out.versions + ).match() + } + ) + } + + } + + test("homo_sapiens - [cram, crai, dragstr_model] - fasta - fai - dict - sites - sites_tbi") { + + when { + process { + """ + input[0] = [ + [ id:'test_cram_sites_dragstr' ], // meta map + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/cram/test.paired_end.sorted.cram', checkIfExists: true), + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/cram/test.paired_end.sorted.cram.crai', checkIfExists: true), + [], + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/gatk/test_paired_end_sorted_dragstrmodel.txt', checkIfExists: true) + ] + input[1] = [ [ id:'test_fa' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta', checkIfExists: true) ] + input[2] = [ [ id:'test_fai' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta.fai', checkIfExists: true) ] + input[3] = [ [ id:'test_dict' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.dict', checkIfExists: true) ] + input[4] = [ [ id:'test_sites' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/vcf/dbsnp_146.hg38.vcf.gz', checkIfExists: true) ] + input[5] = [ [ id:'test_sites_tbi' ], file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/vcf/dbsnp_146.hg38.vcf.gz.tbi', checkIfExists: true) ] + """ + } + } + + then { + assertAll( + { assert process.success }, + { assert snapshot( + file(process.out.vcf[0][1]).name, + file(process.out.tbi[0][1]).name, + process.out.versions + ).match() + } + ) + } + + } + +} diff --git a/modules/UMCUGenetics/gatk4/haplotypecaller_alleles/tests/main.nf.test.snap b/modules/UMCUGenetics/gatk4/haplotypecaller_alleles/tests/main.nf.test.snap new file mode 100644 index 0000000..afcebb0 --- /dev/null +++ b/modules/UMCUGenetics/gatk4/haplotypecaller_alleles/tests/main.nf.test.snap @@ -0,0 +1,58 @@ +{ + "homo_sapiens - [cram, crai] - fasta - fai - dict": { + "content": [ + "test_cram.vcf.gz", + "test_cram.vcf.gz.tbi", + [ + "versions.yml:md5,1481e5d37e1e19e76def11cf7ca38683" + ] + ], + "meta": { + "nf-test": "0.9.2", + "nextflow": "25.04.7" + }, + "timestamp": "2025-09-15T15:26:35.293382853" + }, + "homo_sapiens - [cram, crai] - fasta - fai - dict - sites - sites_tbi": { + "content": [ + "test_cram_sites.vcf.gz", + "test_cram_sites.vcf.gz.tbi", + [ + "versions.yml:md5,1481e5d37e1e19e76def11cf7ca38683" + ] + ], + "meta": { + "nf-test": "0.9.2", + "nextflow": "25.04.7" + }, + "timestamp": "2025-09-15T15:26:51.176703555" + }, + "homo_sapiens - [bam, bai] - fasta - fai - dict": { + "content": [ + "test_bam.vcf.gz", + "test_bam.vcf.gz.tbi", + [ + "versions.yml:md5,1481e5d37e1e19e76def11cf7ca38683" + ] + ], + "meta": { + "nf-test": "0.9.2", + "nextflow": "25.04.7" + }, + "timestamp": "2025-09-15T15:26:19.52126577" + }, + "homo_sapiens - [cram, crai, dragstr_model] - fasta - fai - dict - sites - sites_tbi": { + "content": [ + "test_cram_sites_dragstr.vcf.gz", + "test_cram_sites_dragstr.vcf.gz.tbi", + [ + "versions.yml:md5,1481e5d37e1e19e76def11cf7ca38683" + ] + ], + "meta": { + "nf-test": "0.9.2", + "nextflow": "25.04.7" + }, + "timestamp": "2025-09-15T15:27:07.223882838" + } +} \ No newline at end of file diff --git a/modules/UMCUGenetics/pgscatalog/combine/main.nf b/modules/UMCUGenetics/pgscatalog/combine/main.nf new file mode 100644 index 0000000..a651c19 --- /dev/null +++ b/modules/UMCUGenetics/pgscatalog/combine/main.nf @@ -0,0 +1,27 @@ +process PGSCATALOG_COMBINE { + tag "${meta.id}" + + container "${workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container + ? 'https://depot.galaxyproject.org/singularity/pgscatalog-utils:1.4.4--pyhdfd78af_0' + : 'biocontainers/pgscatalog-utils:1.4.4--pyhdfd78af_0'}" + + input: + tuple val(meta), path(scoring_file) + val assembly_version + + output: + tuple val(meta), path("*_normalised.txt.gz"), emit: normalised_model + path "versions.yml", emit: versions + + script: + def prefix = task.ext.prefix ?: meta.id + """ + pgscatalog-combine \\ + -s ${scoring_file} \\ + -t ${assembly_version} \\ + -o ${prefix}_normalised.txt.gz + + + echo "pgscatalog-combine: 1.4.4" > versions.yml + """ +} diff --git a/modules/UMCUGenetics/pgscatalog/combine/tests/main.nf.test b/modules/UMCUGenetics/pgscatalog/combine/tests/main.nf.test new file mode 100644 index 0000000..748609f --- /dev/null +++ b/modules/UMCUGenetics/pgscatalog/combine/tests/main.nf.test @@ -0,0 +1,27 @@ +nextflow_process { + name "Test Process PGScatalog-combine" + script "../main.nf" + process "PGSCATALOG_COMBINE" + + tag "modules/local" + + + test("Test Correct model formatting "){ + when{ + process{ + """ + input[0] = [[id: 'model'], + file("${projectDir}/assets/models/BCAC_313_PRS_GRCh38.txt", + checkifExists: true)] + input[1] = channel.value("GRCh38") + """ + } + + then { + assertAll( + {assert process.success}, + ) + } + } + } +} diff --git a/modules/UMCUGenetics/pgscatalog/match/main.nf b/modules/UMCUGenetics/pgscatalog/match/main.nf new file mode 100644 index 0000000..8d2db1b --- /dev/null +++ b/modules/UMCUGenetics/pgscatalog/match/main.nf @@ -0,0 +1,42 @@ +process PGSCATALOG_MATCH { + tag "${meta.id}" + + container "${workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container + ? 'https://depot.galaxyproject.org/singularity/pgscatalog-utils:1.4.4--pyhdfd78af_0' + : 'biocontainers/pgscatalog-utils:1.4.4--pyhdfd78af_0'}" + + input: + tuple val(meta), path(pvar) + tuple val(meta2), path(scoring_file) + + + output: + tuple val(meta), path("*_summary.csv"), emit: summary + tuple val(meta), path("*.scorefile.gz"), emit: scorefile + tuple val(meta), path("*_log.csv.gz"), emit: log + path "versions.yml", emit: versions + + script: + def prefix = task.ext.prefix ?: meta.id + def args = task.ext.args ?: "" + """ + pgscatalog-match \\ + ${args} \\ + --dataset ${prefix} \\ + --scorefiles ${scoring_file} \\ + --target ${pvar} \\ + --outdir ./ + + echo "pgscatalog-match: 1.4.4" > versions.yml + """ + + stub: + def prefix = task.ext.prefix ?: "Cohort" + """ + touch ${prefix}_summary.csv + touch ${prefix}.scorefile.gz + touch ${prefix}_log.csv.gz + + echo "pgscatalog-match: 1.4.4" > versions.yml + """ +} diff --git a/modules/UMCUGenetics/pgscatalog/match/tests/main.nf.test b/modules/UMCUGenetics/pgscatalog/match/tests/main.nf.test new file mode 100644 index 0000000..8987caf --- /dev/null +++ b/modules/UMCUGenetics/pgscatalog/match/tests/main.nf.test @@ -0,0 +1,60 @@ +nextflow_process { + name "Test Process pgscatalog/match" + script "../main.nf" + process "PGSCATALOG_MATCH" + + + + tag "modules" + tag "pgscatalog" + tag "pgscatalog/match" + + test("Test SNV ambiguity detection"){ + setup{ + run("PLINK2_VCF") { + script "../../../../nf-core/plink2/vcf/main.nf" + + process{ + """ + input[0] = [[id: "vcf"], file("${projectDir}/assets/test-data/test.vcf",checkIfExists: true)] + """ + } + } + } + + when { + process { + """ + input[0] = [[id: "model"], + file( + "${projectDir}/assets/test-data/norm_test_model_subworkflow.txt.gz", + checkIfExists: true)] + input[1] = PLINK2_VCF.out.pvar_zst + + """ + } + } + then { + // Analysing the csv is necessary as the snapshots do not always match + // (row order differs) + def csv = path(process.out.summary[0][1]).csv().table + def matched = csv.stringColumn("match_status").isEqualTo("matched") + + def n_ambiguous_SNV = csv.where( + csv.booleanColumn("ambiguous").isTrue().and(matched) + ).row(0).getInt("count") + + def n_unambiguous_SNV = csv.where( + csv.booleanColumn("ambiguous").isFalse().and(matched) + ).row(0).getInt("count") + + assertAll ( + {assert process.success}, + {assert n_ambiguous_SNV == 1}, + {assert n_unambiguous_SNV == 3} + ) + } + + + } +} diff --git a/nf-test.config b/nf-test.config index 609c4c9..239ef6e 100644 --- a/nf-test.config +++ b/nf-test.config @@ -11,4 +11,7 @@ config { // run all test with the defined docker profile from the main nextflow.config profile = "" + plugins { + load "nft-utils@1.0.0" + } } diff --git a/subworkflows/UMCUGenetics/bam_haplotypecaller_norm/main.nf b/subworkflows/UMCUGenetics/bam_haplotypecaller_norm/main.nf new file mode 100644 index 0000000..1685dc6 --- /dev/null +++ b/subworkflows/UMCUGenetics/bam_haplotypecaller_norm/main.nf @@ -0,0 +1,59 @@ +include { GATK4_HAPLOTYPECALLER } from '../../../modules/UMCUGenetics/gatk4/haplotypecaller_alleles/main' +include { BCFTOOLS_NORM } from '../../../modules/nf-core/bcftools/norm/main' + +workflow BAM_HAPLOTYPECALLER_NORM { + take: + ch_samplesheet + ch_genome_fasta + ch_genome_index + ch_genome_dict + ch_dbsnp + ch_dbsnp_index + ch_snp_list + ch_snp_vcf + + + main: + + // [model_meta, snplist, vcf, tbi] + ch_per_model = ch_snp_list.join(ch_snp_vcf) + + ch_hc = ch_samplesheet.combine(ch_per_model) + .multiMap { sm, bam, bai, mm, snplist, vcf, tbi -> + def meta = [ + id: "${sm.id}__${mm.id}", + sample_id: sm.id, + model_id: mm.id, + mu: mm.mu, + sd: mm.sd, + alpha: mm.alpha, + ] + input: [meta, bam, bai, snplist, []] + alleles: [[id: meta.id], vcf, tbi] + } + + GATK4_HAPLOTYPECALLER( + ch_hc.input, + ch_genome_fasta, + ch_genome_index, + ch_genome_dict, + ch_dbsnp, + ch_dbsnp_index, + ch_hc.alleles + ) + + BCFTOOLS_NORM( + GATK4_HAPLOTYPECALLER.out.vcf + .join(GATK4_HAPLOTYPECALLER.out.tbi), + ch_genome_fasta + ) + + // Collate software versions + ch_versions = Channel.empty() + ch_versions = ch_versions.mix(GATK4_HAPLOTYPECALLER.out.versions) + + emit: + vcf = BCFTOOLS_NORM.out.vcf + tbi = BCFTOOLS_NORM.out.index + ch_versions = ch_versions +} diff --git a/subworkflows/UMCUGenetics/bam_haplotypecaller_norm/meta.yml b/subworkflows/UMCUGenetics/bam_haplotypecaller_norm/meta.yml new file mode 100644 index 0000000..d869065 --- /dev/null +++ b/subworkflows/UMCUGenetics/bam_haplotypecaller_norm/meta.yml @@ -0,0 +1,12 @@ +name: "bam_haplotypecaller_norm" +description: "PRS tailored haplotypecaller + bcftools norm implementation" +keywords: + - prs + - gatk + - variant_calling + - bcftools + - bam +components: + - gatk4/haplotypecaller + - bcftools/norm: + git_remote: https://github.com/nf-core/modules.git diff --git a/subworkflows/UMCUGenetics/bam_haplotypecaller_norm/tests/main.nf.test b/subworkflows/UMCUGenetics/bam_haplotypecaller_norm/tests/main.nf.test new file mode 100644 index 0000000..fc9d2fe --- /dev/null +++ b/subworkflows/UMCUGenetics/bam_haplotypecaller_norm/tests/main.nf.test @@ -0,0 +1,78 @@ +nextflow_workflow { + name "Test Subworkflow Variant calling" + script "../main.nf" + config "./nextflow.config" + + workflow "BAM_HAPLOTYPECALLER_NORM" + + tag "subworkflows" + tag "subworkflows/local" + tag "subworkflows_UMCUGenetics" + tag "subworkflows/bam_haplotypecaller_norm" + tag "local" + tag "gatk4/haplotypecaller_alleles" + tag "subworkflows/../../modules/nf-core/bcftools/norm" + + + setup{ + nfcoreInitialise("${launchDir}/library/") + nfcoreInstall( + "${launchDir}/library/", + [ + "bcftools/norm" + ] + ) + nfcoreLink("${launchDir}/library/", "${baseDir}/modules/") + } + test("Validate Variant calling subworkflow") { + options "-stub" + + when { + workflow { + """ + input[0] = Channel.of([ + [id: "sample1"], + file("https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/modules/data/genomics/homo_sapiens/illumina/bam/test.paired_end.sorted.bam", + checkIfExists: true), + file("https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/modules/data/genomics/homo_sapiens/illumina/bam/test.paired_end.sorted.bam.bai", + checkIfExists: true) + ]) + input[1] = Channel.of([[id: "genome_fa"], + file("https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/modules/data/genomics/homo_sapiens/genome/genome.fasta", checkIfExists: true)]).first() + input[2] = Channel.of([[id: "genome_fai"], + file("https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/modules/data/genomics/homo_sapiens/genome/genome.fasta.fai", checkIfExists: true)]).first() + input[3] = Channel.of([[id: "genome_dict"], + file("https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/modules/data/genomics/homo_sapiens/genome/genome.dict", checkIfExists: true)]).first() + input[4] = Channel.of([[id: "dbsnp"], + file("https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/modules/data/genomics/homo_sapiens/genome/vcf/dbsnp_146.hg38.vcf.gz",checkIfExists: true)]).first() + input[5] = Channel.of([[id: "dbsnp"], + file("https://raw.githubusercontent.com/nf-core/test-datasets/refs/heads/modules/data/genomics/homo_sapiens/genome/vcf/dbsnp_146.hg38.vcf.gz.tbi",checkIfExists: true)]).first() + input[6] = Channel.of( + [ [id: "modelA"], + file(params.subworkflows_umcu_testdata_base_path + "bam_haplotypecaller_norm/BCAC_313_PRS_GRCh38_snplist.list", checkIfExists: true)], + + [[id: "modelB"], + file(params.subworkflows_umcu_testdata_base_path + "bam_haplotypecaller_norm/BCAC_313_PRS_GRCh38_snplist.list", checkIfExists: true)] + ) + input[7] = Channel.of( + [[id: "modelA"], + file(params.subworkflows_umcu_testdata_base_path + "bam_haplotypecaller_norm/BCAC_313_PRS_GRCh38_genotypes_sorted.vcf.gz", checkIfExists: true), + file(params.subworkflows_umcu_testdata_base_path + "bam_haplotypecaller_norm/BCAC_313_PRS_GRCh38_genotypes_sorted.vcf.gz.tbi", checkIfExists: true) + ], + [[id: "modelB"], + file(params.subworkflows_umcu_testdata_base_path + "bam_haplotypecaller_norm/BCAC_313_PRS_GRCh38_genotypes_sorted.vcf.gz", checkIfExists: true), + file(params.subworkflows_umcu_testdata_base_path + "bam_haplotypecaller_norm/BCAC_313_PRS_GRCh38_genotypes_sorted.vcf.gz.tbi", checkIfExists: true) + ] + + ) + """ + } + } + + then { + assertAll( + { assert workflow.success } + ) + } + } +} diff --git a/subworkflows/UMCUGenetics/bam_haplotypecaller_norm/tests/nextflow.config b/subworkflows/UMCUGenetics/bam_haplotypecaller_norm/tests/nextflow.config new file mode 100644 index 0000000..4526ee8 --- /dev/null +++ b/subworkflows/UMCUGenetics/bam_haplotypecaller_norm/tests/nextflow.config @@ -0,0 +1,11 @@ +process { + + withName: BCFTOOLS_NORM { + ext.args = [ + "-m -", // Split multi-allelic sites + "--output-type z", + "--write-index=tbi" + ].join(' ').trim() + ext.prefix = { "${meta.id}_norm"} + } +} diff --git a/subworkflows/yaml-schema.json b/subworkflows/yaml-schema.json new file mode 100644 index 0000000..619135e --- /dev/null +++ b/subworkflows/yaml-schema.json @@ -0,0 +1,121 @@ +{ + "$schema": "http://json-schema.org/draft-07/schema", + "title": "Meta yaml", + "description": "Validate the meta yaml file for an nf-core subworkflow", + "type": "object", + "properties": { + "name": { + "type": "string", + "description": "Name of the subworkflow" + }, + "description": { + "type": "string", + "description": "Description of the subworkflow" + }, + "authors": { + "type": "array", + "description": "Authors of the subworkflow", + "items": { + "type": "string" + } + }, + "maintainers": { + "type": "array", + "description": "Maintainers of the subworkflow", + "items": { + "type": "string" + } + }, + "components": { + "type": "array", + "description": "Modules and subworkflows used in the subworkflow", + "items": { + "type": ["object", "string"] + }, + "minItems": 0 + }, + "keywords": { + "type": "array", + "description": "Keywords for the module", + "items": { + "type": "string" + }, + "minItems": 3 + }, + "input": { + "type": "array", + "description": "Input channels for the subworkflow", + "items": { + "type": "object", + "patternProperties": { + ".*": { + "type": "object", + "properties": { + "type": { + "type": "string", + "description": "Type of the input channel" + }, + "description": { + "type": "string", + "description": "Description of the input channel" + }, + "pattern": { + "type": "string", + "description": "Pattern of the input channel, given in Java glob syntax" + }, + "default": { + "type": ["string", "number", "boolean", "array", "object"], + "description": "Default value for the input channel" + }, + "enum": { + "type": "array", + "description": "List of allowed values for the input channel", + "items": { + "type": ["string", "number", "boolean", "array", "object"] + }, + "uniqueItems": true + } + }, + "required": ["description"] + } + } + } + }, + "output": { + "type": "array", + "description": "Output channels for the subworkflow", + "items": { + "type": "object", + "patternProperties": { + ".*": { + "type": "object", + "properties": { + "type": { + "type": "string", + "description": "Type of the output channel" + }, + "description": { + "type": "string", + "description": "Description of the output channel" + }, + "pattern": { + "type": "string", + "description": "Pattern of the input channel, given in Java glob syntax" + }, + "enum": { + "type": "array", + "description": "List of allowed values for the output channel", + "items": { + "type": ["string", "number", "boolean", "array", "object"] + }, + "uniqueItems": true + } + }, + "required": ["description"] + } + } + } + } + }, + "required": ["name", "description", "keywords", "authors", "output", "components"] +} diff --git a/tests/config/nf-test.config b/tests/config/nf-test.config index 3a1cca8..f946df0 100644 --- a/tests/config/nf-test.config +++ b/tests/config/nf-test.config @@ -1,8 +1,12 @@ +nextflow.enable.moduleBinaries = true + params { publish_dir_mode = "copy" singularity_pull_docker_container = false test_data_base = 'https://raw.githubusercontent.com/nf-core/test-datasets/modules' modules_testdata_base_path = 'https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/' + modules_umcu_testdata_base_path = 'https://raw.githubusercontent.com/UMCUGenetics/DxNextflowTestData/main/Modules/' + subworkflows_umcu_testdata_base_path = 'https://raw.githubusercontent.com/UMCUGenetics/DxNextflowTestData/main/Subworkflows/' } process {