Build a Nextflow cfDNA Variant Pipeline

Build a Nextflow DSL2 pipeline that classifies cancer variants in cfDNA data.

Introduction

30 Second Summary

Liquid biopsies detect cancer signals by analyzing tiny fragments of tumor DNA floating in a patient's blood. The challenge is distinguishing real cancer mutations from sequencing noise, and machine learning is how modern labs solve it.

In this project, you will build a Nextflow bioinformatics pipeline that classifies cell-free DNA (cfDNA) variants as true somatic mutations or sequencing artifacts using a Random Forest classifier.

What You'll Build

A containerized Nextflow DSL2 pipeline that orchestrates synthetic cfDNA data generation, model training with class imbalance handling, and precision-recall evaluation.

By the end of this project, you'll have:

  • A Nextflow DSL2 pipeline with three chained processes that take variant data from generation through training to evaluation.
  • A Docker container packaging your Python ML environment for reproducible execution on any machine.
  • A nextflow.config with switchable profiles for local Docker execution and HPC Slurm submission.
  • Secret Mission: Add multi-sample parallelization so the pipeline processes multiple patient samples concurrently and aggregates results.

Are there any prerequisites?

You should be comfortable with basic Python (functions, pandas DataFrames, pip) and familiar with using the terminal. You also need Docker installed and running on your machine before starting.

Set Up the Project

Before writing any pipeline code, you need a clean project structure that follows Nextflow conventions. Nextflow pipelines use specific directories to organize executable scripts, configurations, and workflow definitions.

In this step, you will create the directory structure for your pipeline and open the project in VS Code so you are ready to start coding.

In this step, get ready to:
  • Create the bin/ directory for your pipeline scripts.
  • Open VS Code at the project root.
Create the bin/ directory

Nextflow has a special convention for the bin/ directory. Any executable scripts placed inside it are automatically added to the $PATH for all pipeline tasks.

  • Switch back to your terminal. Make sure you are inside your cfdna-pipeline folder (the project directory from the previous step).
  • Create the bin/ directory by running this command:
mkdir bin

What is the bin/ directory in Nextflow?

In a Nextflow pipeline repository, the bin/ directory has special meaning. Nextflow automatically adds it to the $PATH for every task, so scripts inside it can be called like regular commands from any process.

If your tasks run on a remote machine or in a container, Nextflow copies the bin/ directory there automatically. This keeps your scripts versioned alongside your workflow.

  • Confirm the directory exists by running:
ls

You should see bin listed in the output alongside your venv directory.

Don't see the bin directory?

Make sure you ran mkdir bin from inside your cfdna-pipeline folder. If you are in the wrong directory, the folder was created elsewhere.

Run pwd to check your current location. You should see a path ending in cfdna-pipeline.

Help me find where my bin directory was created.

Open VS Code at the project root

With the directory structure in place, open the project in VS Code so you have a code editor ready for the next steps.

  • Press Cmd+Space (macOS) or the Windows key (Windows) to open your search bar.
  • Type Visual Studio Code and press Enter to open it.
  • In VS Code, click File in the top menu bar.
  • Click Open Folder.
  • Navigate to your Desktop and select the cfdna-pipeline folder, then click Open.

You should see the bin folder and venv folder listed in the Explorer sidebar on the left.

Trust the workspace

VS Code may show a dialog asking if you trust the authors of this folder. Click Yes, I trust the authors to enable all editor features.

Don't see the Explorer sidebar?

Press Cmd+Shift+E (macOS) or Ctrl+Shift+E (Windows) to toggle the Explorer panel.

Help me set up my VS Code workspace for this pipeline project.

Your project structure is ready. Next up, you will write the Python scripts that generate synthetic cfDNA data, train the classifier, and evaluate performance.

Write the Python Scripts

Your environment is ready with Nextflow, Python, and Docker all verified. Now it is time to write the core logic that your pipeline will orchestrate.

A real cfDNA variant-calling workflow needs to generate training data, fit a classifier, and evaluate the model's performance. You will create three standalone Python scripts that handle each of these stages. Later, Nextflow will chain them together as pipeline processes.

In this step, get ready to:
  • Create a data generation script that outputs synthetic cfDNA variant features as CSV.
  • Create a model training script that fits a Random Forest classifier with cross-validation.
  • Create an evaluation script that produces precision-recall plots and feature importance charts.
Create the data generation script

Real cfDNA sequencing datasets are multi-gigabyte BAM files that require specialized bioinformatics tools to process. Synthetic data with realistic feature distributions lets you build and validate the full pipeline without any of that overhead.

Why synthetic data instead of real cfDNA samples?

Downloading real cfDNA data from repositories like TCGA requires 1.7GB+ files and BAM processing tools. The synthetic approach generates features with distributions modeled on published cfDNA research, teaching the same ML and pipeline engineering concepts in a fraction of the time.

The data generation script has two parts. First, a function that builds the feature matrix. Then, a CLI interface that writes it to CSV.

  • In VS Code, create a new file called generate_data.py inside the bin/ folder.
  • Paste the following code, which defines the generate_variants() function:
#!/usr/bin/env python3
"""Generate synthetic cfDNA variant feature data for classifier training."""

import argparse
import numpy as np
import pandas as pd


def generate_variants(n_variants=2000, true_somatic_frac=0.10, seed=42):
    """Generate synthetic cfDNA variant features with class imbalance."""
    rng = np.random.default_rng(seed)

    n_true = int(n_variants * true_somatic_frac)
    n_artifact = n_variants - n_true

    # True somatic variants: higher depth, lower AF, moderate strand bias
    true_features = {
        'read_depth': rng.normal(150, 40, n_true).clip(10, 500),
        'allele_frequency': rng.beta(2, 20, n_true).clip(0.005, 0.5),
        'strand_bias': rng.beta(5, 5, n_true),
        'fragment_length': rng.normal(167, 20, n_true).clip(100, 250),
        'mapping_quality': rng.normal(55, 5, n_true).clip(20, 60),
        'gc_content': rng.normal(0.45, 0.10, n_true).clip(0.2, 0.8),
        'in_cosmic': rng.binomial(1, 0.7, n_true),
        'in_dbsnp': rng.binomial(1, 0.1, n_true),
    }

    # Artifacts: lower depth, variable AF, extreme strand bias
    artifact_features = {
        'read_depth': rng.normal(80, 50, n_artifact).clip(5, 500),
        'allele_frequency': rng.beta(1, 5, n_artifact).clip(0.005, 0.5),
        'strand_bias': rng.beta(1, 3, n_artifact),
        'fragment_length': rng.normal(140, 35, n_artifact).clip(80, 300),
        'mapping_quality': rng.normal(40, 12, n_artifact).clip(0, 60),
        'gc_content': rng.normal(0.50, 0.15, n_artifact).clip(0.1, 0.9),
        'in_cosmic': rng.binomial(1, 0.05, n_artifact),
        'in_dbsnp': rng.binomial(1, 0.6, n_artifact),
    }

    true_df = pd.DataFrame(true_features)
    true_df['label'] = 1

    artifact_df = pd.DataFrame(artifact_features)
    artifact_df['label'] = 0

    df = pd.concat([true_df, artifact_df], ignore_index=True)
    df = df.sample(frac=1, random_state=seed).reset_index(drop=True)

    return df

What does this code do?

  • np.random.default_rng(seed) creates a reproducible random number generator. Using a fixed seed means the same data is generated every time.
  • True somatic variants (10% of the data) get higher read depth, lower allele frequency, and a 70% chance of appearing in the COSMIC database. These mimic real cancer mutations.
  • Artifacts (90% of the data) get lower depth, more extreme strand bias, and a 60% chance of appearing in dbSNP. These mimic sequencing noise that should be filtered out.
  • The 10/90 class imbalance reflects the real-world ratio where true somatic variants are rare compared to artifacts in cfDNA data.

You will not see a change yet. The next chunk adds the CLI interface that actually runs the function and writes output.

  • Below the generate_variants function in the same file, add the main() function and entry point:

def main():
    parser = argparse.ArgumentParser(description='Generate synthetic cfDNA variant data')
    parser.add_argument('--output', type=str, default='variants.csv',
                        help='Output CSV file path')
    parser.add_argument('--n-variants', type=int, default=2000,
                        help='Number of variants to generate')
    parser.add_argument('--seed', type=int, default=42,
                        help='Random seed for reproducibility')
    args = parser.parse_args()

    df = generate_variants(n_variants=args.n_variants, seed=args.seed)
    df.to_csv(args.output, index=False)
    print(f"Generated {len(df)} variants ({df['label'].sum()} true somatic) -> {args.output}")


if __name__ == '__main__':
    main()

What does this code do?

  • argparse defines a CLI with three flags: --output (CSV destination), --n-variants (how many to generate), and --seed (for reproducibility).
  • The if __name__ == '__main__' guard ensures the script runs as a standalone command but can also be imported as a module by other scripts.
  • Save bin/generate_data.py.
  • Run the script to verify it produces output:
python3 bin/generate_data.py --output variants.csv
  • You should see output like: Generated 2000 variants (200 true somatic) -> variants.csv.

Seeing ModuleNotFoundError?

  • Make sure your virtual environment is active. You should see (venv) at the start of your terminal prompt.
  • If it is not active, run source venv/bin/activate from your project root, then try again.

help me fix a ModuleNotFoundError when running my generate_data.py script

✔️ Awesome, I've got everything!

Great. Make sure you have saved bin/generate_data.py before moving on.

ⓧ I'd like to double check the full code

#!/usr/bin/env python3
"""Generate synthetic cfDNA variant feature data for classifier training."""

import argparse
import numpy as np
import pandas as pd


def generate_variants(n_variants=2000, true_somatic_frac=0.10, seed=42):
    """Generate synthetic cfDNA variant features with class imbalance."""
    rng = np.random.default_rng(seed)

    n_true = int(n_variants * true_somatic_frac)
    n_artifact = n_variants - n_true

    # True somatic variants: higher depth, lower AF, moderate strand bias
    true_features = {
        'read_depth': rng.normal(150, 40, n_true).clip(10, 500),
        'allele_frequency': rng.beta(2, 20, n_true).clip(0.005, 0.5),
        'strand_bias': rng.beta(5, 5, n_true),
        'fragment_length': rng.normal(167, 20, n_true).clip(100, 250),
        'mapping_quality': rng.normal(55, 5, n_true).clip(20, 60),
        'gc_content': rng.normal(0.45, 0.10, n_true).clip(0.2, 0.8),
        'in_cosmic': rng.binomial(1, 0.7, n_true),
        'in_dbsnp': rng.binomial(1, 0.1, n_true),
    }

    # Artifacts: lower depth, variable AF, extreme strand bias
    artifact_features = {
        'read_depth': rng.normal(80, 50, n_artifact).clip(5, 500),
        'allele_frequency': rng.beta(1, 5, n_artifact).clip(0.005, 0.5),
        'strand_bias': rng.beta(1, 3, n_artifact),
        'fragment_length': rng.normal(140, 35, n_artifact).clip(80, 300),
        'mapping_quality': rng.normal(40, 12, n_artifact).clip(0, 60),
        'gc_content': rng.normal(0.50, 0.15, n_artifact).clip(0.1, 0.9),
        'in_cosmic': rng.binomial(1, 0.05, n_artifact),
        'in_dbsnp': rng.binomial(1, 0.6, n_artifact),
    }

    true_df = pd.DataFrame(true_features)
    true_df['label'] = 1

    artifact_df = pd.DataFrame(artifact_features)
    artifact_df['label'] = 0

    df = pd.concat([true_df, artifact_df], ignore_index=True)
    df = df.sample(frac=1, random_state=seed).reset_index(drop=True)

    return df


def main():
    parser = argparse.ArgumentParser(description='Generate synthetic cfDNA variant data')
    parser.add_argument('--output', type=str, default='variants.csv',
                        help='Output CSV file path')
    parser.add_argument('--n-variants', type=int, default=2000,
                        help='Number of variants to generate')
    parser.add_argument('--seed', type=int, default=42,
                        help='Random seed for reproducibility')
    args = parser.parse_args()

    df = generate_variants(n_variants=args.n_variants, seed=args.seed)
    df.to_csv(args.output, index=False)
    print(f"Generated {len(df)} variants ({df['label'].sum()} true somatic) -> {args.output}")


if __name__ == '__main__':
    main()
Create the model training script

The training script takes the CSV from the previous stage and fits a Random Forest classifier. It uses stratified cross-validation to estimate performance on held-out data, then saves the trained model and metrics.

Why Random Forest with class_weight='balanced'?

In cfDNA data, true somatic variants are rare (around 10% of calls). A naive classifier would get 90% accuracy by labeling everything as an artifact. Setting class_weight='balanced' tells scikit-learn to automatically upweight the minority class during training, producing a model that is sensitive to true mutations.

This approach is modeled on Maruzani et al. 2025, which used Random Forest with resampling to predict high-confidence somatic variants in cfDNA whole-exome sequencing data.

  • Create a new file called train_model.py inside the bin/ folder.
  • Paste the following code:
#!/usr/bin/env python3
"""Train a Random Forest classifier on cfDNA variant features."""

import argparse
import json
import pickle

import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import StratifiedKFold, cross_validate
from sklearn.metrics import make_scorer, precision_score, recall_score, f1_score


def main():
    parser = argparse.ArgumentParser(description='Train cfDNA variant classifier')
    parser.add_argument('--input', type=str, required=True,
                        help='Input CSV with variant features')
    parser.add_argument('--model-output', type=str, default='model.pkl',
                        help='Output path for trained model')
    parser.add_argument('--metrics-output', type=str, default='metrics.json',
                        help='Output path for metrics JSON')
    args = parser.parse_args()

    df = pd.read_csv(args.input)
    feature_cols = ['read_depth', 'allele_frequency', 'strand_bias',
                    'fragment_length', 'mapping_quality', 'gc_content',
                    'in_cosmic', 'in_dbsnp']
    X = df[feature_cols]
    y = df['label']

    clf = RandomForestClassifier(
        n_estimators=100,
        class_weight='balanced',
        random_state=42,
        n_jobs=-1
    )

    scoring = {
        'precision': make_scorer(precision_score, zero_division=0),
        'recall': make_scorer(recall_score, zero_division=0),
        'f1': make_scorer(f1_score, zero_division=0),
    }

    cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    cv_results = cross_validate(clf, X, y, cv=cv, scoring=scoring)

    clf.fit(X, y)

    metrics = {
        'cv_precision_mean': float(np.mean(cv_results['test_precision'])),
        'cv_recall_mean': float(np.mean(cv_results['test_recall'])),
        'cv_f1_mean': float(np.mean(cv_results['test_f1'])),
        'cv_precision_std': float(np.std(cv_results['test_precision'])),
        'cv_recall_std': float(np.std(cv_results['test_recall'])),
        'cv_f1_std': float(np.std(cv_results['test_f1'])),
        'n_samples': len(df),
        'n_positive': int(y.sum()),
        'n_negative': int(len(y) - y.sum()),
        'feature_importances': dict(zip(feature_cols,
                                        [float(x) for x in clf.feature_importances_])),
    }

    with open(args.model_output, 'wb') as f:
        pickle.dump(clf, f)

    with open(args.metrics_output, 'w') as f:
        json.dump(metrics, f, indent=2)

    print(f"Model trained. CV F1: {metrics['cv_f1_mean']:.3f} +/- {metrics['cv_f1_std']:.3f}")
    print(f"Saved model -> {args.model_output}")
    print(f"Saved metrics -> {args.metrics_output}")


if __name__ == '__main__':
    main()

What does this code do?

  • Loads the variant CSV and extracts the 8 feature columns plus the label column as the target.
  • Creates a RandomForestClassifier with 100 trees, balanced class weights, and a fixed random state for reproducibility.
  • Runs 5-fold StratifiedKFold cross-validation, scoring precision, recall, and F1. Stratified folds preserve the 10/90 class ratio in each split.
  • Fits the final model on all data after CV, then saves both the model (as a pickle file) and the metrics (as JSON).
  • Save bin/train_model.py.
  • Run the training script, pointing it at the CSV you generated earlier:
python3 bin/train_model.py --input variants.csv
  • You should see output showing the cross-validation F1 score and confirmation that model.pkl and metrics.json were saved.

Seeing a FileNotFoundError?

  • Make sure you ran generate_data.py first to produce variants.csv in your project root.
  • Check that your terminal is in the project root directory (the folder containing bin/), not inside bin/ itself.

help me fix a FileNotFoundError in my train_model.py script

✔️ Awesome, I've got everything!

Great. Make sure you have saved bin/train_model.py before moving on.

ⓧ I'd like to double check the full code

#!/usr/bin/env python3
"""Train a Random Forest classifier on cfDNA variant features."""

import argparse
import json
import pickle

import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import StratifiedKFold, cross_validate
from sklearn.metrics import make_scorer, precision_score, recall_score, f1_score


def main():
    parser = argparse.ArgumentParser(description='Train cfDNA variant classifier')
    parser.add_argument('--input', type=str, required=True,
                        help='Input CSV with variant features')
    parser.add_argument('--model-output', type=str, default='model.pkl',
                        help='Output path for trained model')
    parser.add_argument('--metrics-output', type=str, default='metrics.json',
                        help='Output path for metrics JSON')
    args = parser.parse_args()

    df = pd.read_csv(args.input)
    feature_cols = ['read_depth', 'allele_frequency', 'strand_bias',
                    'fragment_length', 'mapping_quality', 'gc_content',
                    'in_cosmic', 'in_dbsnp']
    X = df[feature_cols]
    y = df['label']

    clf = RandomForestClassifier(
        n_estimators=100,
        class_weight='balanced',
        random_state=42,
        n_jobs=-1
    )

    scoring = {
        'precision': make_scorer(precision_score, zero_division=0),
        'recall': make_scorer(recall_score, zero_division=0),
        'f1': make_scorer(f1_score, zero_division=0),
    }

    cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    cv_results = cross_validate(clf, X, y, cv=cv, scoring=scoring)

    clf.fit(X, y)

    metrics = {
        'cv_precision_mean': float(np.mean(cv_results['test_precision'])),
        'cv_recall_mean': float(np.mean(cv_results['test_recall'])),
        'cv_f1_mean': float(np.mean(cv_results['test_f1'])),
        'cv_precision_std': float(np.std(cv_results['test_precision'])),
        'cv_recall_std': float(np.std(cv_results['test_recall'])),
        'cv_f1_std': float(np.std(cv_results['test_f1'])),
        'n_samples': len(df),
        'n_positive': int(y.sum()),
        'n_negative': int(len(y) - y.sum()),
        'feature_importances': dict(zip(feature_cols,
                                        [float(x) for x in clf.feature_importances_])),
    }

    with open(args.model_output, 'wb') as f:
        pickle.dump(clf, f)

    with open(args.metrics_output, 'w') as f:
        json.dump(metrics, f, indent=2)

    print(f"Model trained. CV F1: {metrics['cv_f1_mean']:.3f} +/- {metrics['cv_f1_std']:.3f}")
    print(f"Saved model -> {args.model_output}")
    print(f"Saved metrics -> {args.metrics_output}")


if __name__ == '__main__':
    main()
Create the evaluation script

The evaluation script takes the trained model and the original data, then produces visual outputs. A precision-recall curve shows how well the classifier distinguishes true somatic variants from artifacts at different confidence thresholds. Feature importance plots reveal which biological signals the model relies on most.

  • Create a new file called evaluate_model.py inside the bin/ folder.
  • Paste the following code:
#!/usr/bin/env python3
"""Evaluate the trained cfDNA variant classifier and produce plots."""

import argparse
import json
import pickle

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from sklearn.metrics import precision_recall_curve, average_precision_score
from sklearn.model_selection import cross_val_predict, StratifiedKFold


def main():
    parser = argparse.ArgumentParser(description='Evaluate cfDNA classifier')
    parser.add_argument('--model', type=str, required=True,
                        help='Path to trained model pickle')
    parser.add_argument('--data', type=str, required=True,
                        help='Path to variant CSV data')
    parser.add_argument('--output-dir', type=str, default='.',
                        help='Directory for output plots and report')
    args = parser.parse_args()

    with open(args.model, 'rb') as f:
        clf = pickle.load(f)

    df = pd.read_csv(args.data)
    feature_cols = ['read_depth', 'allele_frequency', 'strand_bias',
                    'fragment_length', 'mapping_quality', 'gc_content',
                    'in_cosmic', 'in_dbsnp']
    X = df[feature_cols]
    y = df['label']

    cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    y_scores = cross_val_predict(clf, X, y, cv=cv, method='predict_proba')[:, 1]

    precision, recall, thresholds = precision_recall_curve(y, y_scores)
    ap_score = average_precision_score(y, y_scores)

    fig, ax = plt.subplots(figsize=(8, 6))
    ax.plot(recall, precision, linewidth=2, color='#2563eb')
    ax.set_xlabel('Recall (Sensitivity)')
    ax.set_ylabel('Precision')
    ax.set_title(f'Precision-Recall Curve (AP = {ap_score:.3f})')
    ax.set_xlim([0, 1])
    ax.set_ylim([0, 1])
    ax.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig(f'{args.output_dir}/precision_recall_curve.png', dpi=150)
    plt.close()

    importances = clf.feature_importances_
    indices = np.argsort(importances)[::-1]

    fig, ax = plt.subplots(figsize=(8, 5))
    ax.barh(range(len(feature_cols)),
            importances[indices],
            color='#2563eb', alpha=0.8)
    ax.set_yticks(range(len(feature_cols)))
    ax.set_yticklabels([feature_cols[i] for i in indices])
    ax.set_xlabel('Feature Importance')
    ax.set_title('Random Forest Feature Importances')
    ax.invert_yaxis()
    plt.tight_layout()
    plt.savefig(f'{args.output_dir}/feature_importance.png', dpi=150)
    plt.close()

    report = {
        'average_precision': float(ap_score),
        'n_variants': len(df),
        'class_balance': f"{y.sum()}/{len(y)} positive ({100*y.mean():.1f}%)",
        'top_features': [feature_cols[i] for i in indices[:3]],
    }

    with open(f'{args.output_dir}/evaluation_report.json', 'w') as f:
        json.dump(report, f, indent=2)

    print(f"Average Precision: {ap_score:.3f}")
    print(f"Top features: {report['top_features']}")
    print(f"Plots saved to {args.output_dir}/")


if __name__ == '__main__':
    main()

What does this code do?

  • matplotlib.use('Agg') forces matplotlib to render plots to files without needing a display. This is essential for running inside Docker containers later.
  • cross_val_predict generates out-of-fold probability scores, giving an honest estimate of how well the model generalizes to unseen data.
  • The precision-recall curve plots the tradeoff between catching true variants (recall) and avoiding false positives (precision). Average precision summarizes the curve as a single number.
  • Feature importance reveals which variant features the model relies on most, helping bioinformaticians interpret what drives classification decisions.
  • Save bin/evaluate_model.py.
  • Run the evaluation script, pointing it at the model and data you generated:
python3 bin/evaluate_model.py --model model.pkl --data variants.csv
  • You should see the average precision score and a list of top features printed to the terminal.
  • Confirm the output files exist by running:
ls *.png *.json
  • You should see precision_recall_curve.png, feature_importance.png, evaluation_report.json, and metrics.json listed.

Script runs but no PNG files appear?

  • Check that matplotlib is installed in your virtual environment. Run pip list | grep matplotlib to confirm.
  • If you see an error about Agg backend, make sure matplotlib.use('Agg') appears before import matplotlib.pyplot as plt in your script.

help me debug why my evaluate_model.py script is not producing PNG plot files

✔️ Awesome, I've got everything!

Great. Make sure all three scripts are saved in your bin/ folder.

ⓧ I'd like to double check the full code

#!/usr/bin/env python3
"""Evaluate the trained cfDNA variant classifier and produce plots."""

import argparse
import json
import pickle

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from sklearn.metrics import precision_recall_curve, average_precision_score
from sklearn.model_selection import cross_val_predict, StratifiedKFold


def main():
    parser = argparse.ArgumentParser(description='Evaluate cfDNA classifier')
    parser.add_argument('--model', type=str, required=True,
                        help='Path to trained model pickle')
    parser.add_argument('--data', type=str, required=True,
                        help='Path to variant CSV data')
    parser.add_argument('--output-dir', type=str, default='.',
                        help='Directory for output plots and report')
    args = parser.parse_args()

    with open(args.model, 'rb') as f:
        clf = pickle.load(f)

    df = pd.read_csv(args.data)
    feature_cols = ['read_depth', 'allele_frequency', 'strand_bias',
                    'fragment_length', 'mapping_quality', 'gc_content',
                    'in_cosmic', 'in_dbsnp']
    X = df[feature_cols]
    y = df['label']

    cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    y_scores = cross_val_predict(clf, X, y, cv=cv, method='predict_proba')[:, 1]

    precision, recall, thresholds = precision_recall_curve(y, y_scores)
    ap_score = average_precision_score(y, y_scores)

    fig, ax = plt.subplots(figsize=(8, 6))
    ax.plot(recall, precision, linewidth=2, color='#2563eb')
    ax.set_xlabel('Recall (Sensitivity)')
    ax.set_ylabel('Precision')
    ax.set_title(f'Precision-Recall Curve (AP = {ap_score:.3f})')
    ax.set_xlim([0, 1])
    ax.set_ylim([0, 1])
    ax.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig(f'{args.output_dir}/precision_recall_curve.png', dpi=150)
    plt.close()

    importances = clf.feature_importances_
    indices = np.argsort(importances)[::-1]

    fig, ax = plt.subplots(figsize=(8, 5))
    ax.barh(range(len(feature_cols)),
            importances[indices],
            color='#2563eb', alpha=0.8)
    ax.set_yticks(range(len(feature_cols)))
    ax.set_yticklabels([feature_cols[i] for i in indices])
    ax.set_xlabel('Feature Importance')
    ax.set_title('Random Forest Feature Importances')
    ax.invert_yaxis()
    plt.tight_layout()
    plt.savefig(f'{args.output_dir}/feature_importance.png', dpi=150)
    plt.close()

    report = {
        'average_precision': float(ap_score),
        'n_variants': len(df),
        'class_balance': f"{y.sum()}/{len(y)} positive ({100*y.mean():.1f}%)",
        'top_features': [feature_cols[i] for i in indices[:3]],
    }

    with open(f'{args.output_dir}/evaluation_report.json', 'w') as f:
        json.dump(report, f, indent=2)

    print(f"Average Precision: {ap_score:.3f}")
    print(f"Top features: {report['top_features']}")
    print(f"Plots saved to {args.output_dir}/")


if __name__ == '__main__':
    main()

All three scripts run independently and produce the expected outputs. Next up, you will package this Python environment into a Docker container so the pipeline runs identically on any machine.

Containerize with Docker

Your Python scripts run perfectly on your local machine. But when this pipeline moves to a cloud cluster or a colleague's workstation, missing packages or version mismatches can silently break results.

Docker solves this by packaging your scripts, their exact dependencies, and the runtime into a single portable image. In this step, you will write a Dockerfile, build the image, and verify the container works.

In this step, get ready to:
  • Write a Dockerfile that defines the container environment.
  • Build and tag the Docker image.
  • Test that the container produces expected output.
Write the Dockerfile

A Dockerfile is a recipe that tells Docker exactly how to build your image. Each line adds a layer: the base OS, installed packages, copied files, and environment configuration.

  • In VS Code, create a new file called Dockerfile in the root of your project folder (the same level as the bin/ directory).
  • Paste the following into Dockerfile:
FROM python:3.11-slim

WORKDIR /app

# Install the exact package versions used in the project
RUN pip install --no-cache-dir \
    scikit-learn==1.6.1 \
    pandas==2.2.3 \
    numpy==2.2.3 \
    matplotlib==3.10.1

# Copy the pipeline scripts into the container
COPY bin/ /app/bin/

# Make the scripts executable
RUN chmod +x /app/bin/*.py

# Add the scripts directory to PATH
ENV PATH="/app/bin:${PATH}"

What does this Dockerfile do?

  • FROM python:3.11-slim starts from a minimal Python 3.11 image. The slim variant keeps the image small by excluding build tools you don't need.
  • RUN pip install --no-cache-dir installs the exact pinned versions of your dependencies. The --no-cache-dir flag keeps the image size down.
  • COPY bin/ /app/bin/ brings your three Python scripts into the container's filesystem.
  • ENV PATH adds /app/bin to the container's PATH so Nextflow can call your scripts directly.
  • Save Dockerfile.
  • Confirm the file sits next to the bin/ folder in your project root by running:
ls Dockerfile bin/

You should see both Dockerfile and the bin/ directory listed.

Don't see Dockerfile?

Make sure you created the file at the project root, not inside the bin/ folder. The file must be named exactly Dockerfile with a capital D and no file extension.

Help me find where my Dockerfile ended up.

Why pin exact package versions?

Without pinned versions, pip install scikit-learn grabs whatever version is latest when the image is built. If a breaking change ships between runs, your pipeline could silently produce different results.

Pinning versions like scikit-learn==1.6.1 guarantees identical behaviour every time you rebuild the image. This is essential for reproducible science.

Build the Docker image

With the Dockerfile written, you can now build the image. This installs all the packages and bakes your scripts into a portable artifact.

  • Build and tag your image by running this command:
docker build -t cfdna-classifier:latest .

What does this command do?

  • docker build reads the Dockerfile and executes each instruction as a layer.
  • -t cfdna-classifier:latest tags the resulting image with a name and version. Nextflow will reference this tag later.
  • . tells Docker to use the current directory as the build context (so COPY bin/ can find your scripts).

The build takes 1-2 minutes as Docker downloads the base image and installs packages. You should see output ending with a line containing your image tag.

  • Verify the image was created by running:
docker images

You should see cfdna-classifier in the REPOSITORY column with the latest tag.

Build failed with a pip error?

If you see a network timeout or package resolution failure, check that Docker Desktop is running and has internet access. Try running the build command again.

If the error mentions a missing Dockerfile, make sure your terminal is in the project root (the folder containing Dockerfile).

Help me debug my Docker build failure.

Test the container

A successful build does not guarantee the scripts work inside the container. A quick test run confirms everything is wired up correctly.

  • Run the data generation script inside the container by running this command:
docker run --rm cfdna-classifier:latest python /app/bin/generate_data.py --output /dev/stdout | head

What does this command do?

  • docker run --rm starts a new container and removes it automatically when the command finishes.
  • --output /dev/stdout writes the CSV to standard output instead of a file, so you can see it immediately.
  • | head shows only the first 10 lines so you can confirm the CSV header and a few data rows without flooding your terminal.

You should see a CSV header row (read_depth,allele_frequency,strand_bias,...) followed by numeric data rows. This confirms your scripts and all their dependencies are working correctly inside the container.

Seeing a ModuleNotFoundError?

This means a package did not install properly during the build. Rebuild the image to ensure the RUN pip install layer completes without errors:

docker build --no-cache -t cfdna-classifier:latest .

The --no-cache flag forces Docker to reinstall all packages from scratch.

Help me fix my container import error.

✔️ Awesome, I've got everything!

Great. Double check you have saved your Dockerfile and that docker images shows cfdna-classifier:latest.

ⓧ I'd like to double check the full code

Here is the complete Dockerfile:

FROM python:3.11-slim

WORKDIR /app

RUN pip install --no-cache-dir \
    scikit-learn==1.6.1 \
    pandas==2.2.3 \
    numpy==2.2.3 \
    matplotlib==3.10.1

COPY bin/ /app/bin/

RUN chmod +x /app/bin/*.py

ENV PATH="/app/bin:${PATH}"

Your ML environment is now locked inside a portable container. Next up, you will wire everything together into a Nextflow pipeline that orchestrates data generation, training, and evaluation as chained processes.

Build the Nextflow Pipeline

You have a Docker container with all three Python scripts packaged inside. Now you need something to orchestrate them in the correct order, passing outputs from one script as inputs to the next.

That orchestrator is Nextflow. In this step, you will write a DSL2 pipeline definition that chains your data generation, model training, and evaluation into a single reproducible workflow. You will also configure execution profiles so the same pipeline can run locally with Docker or on a Slurm HPC cluster.

In this step, get ready to:
  • Create main.nf with three DSL2 processes wired together via channels.
  • Create nextflow.config with local Docker and HPC Slurm profiles.
  • Run the full pipeline and verify all outputs are published.
Create the Nextflow pipeline file

Nextflow DSL2 uses processes as self-contained units of work and channels to pass data between them. Each process declares its container, inputs, outputs, and a shell script to execute.

Why Nextflow over Snakemake?

Nextflow's dataflow model maps naturally to cfDNA processing where samples flow through stages independently. The nf-core ecosystem provides 300+ curated bioinformatics pipelines built on Nextflow, making it the dominant workflow manager in clinical genomics.

Switching between local Docker execution and HPC Slurm is a single command-line flag change. No rewriting of the pipeline logic is needed.

  • Create a new file called main.nf in your project root (the same folder that contains bin/ and your Dockerfile).
  • Add the pipeline parameters and the first process by pasting this code:
params.n_variants = 2000
params.seed = 42
params.outdir = 'results'

process GENERATE_DATA {
    container 'cfdna-classifier:latest'

    output:
    path 'variants.csv', emit: variants

    script:
    """
    python3 /app/bin/generate_data.py \\
        --output variants.csv \\
        --n-variants ${params.n_variants} \\
        --seed ${params.seed}
    """
}

What does this code do?

The three params lines define pipeline-level parameters that can be overridden at runtime with --n-variants 5000 on the command line.

The GENERATE_DATA process declares which Docker container to use, what file it produces (variants.csv), and the shell command to run inside the container. The emit: variants label names this output channel so downstream processes can reference it.

  • Below the GENERATE_DATA process, add the TRAIN_MODEL process:
process TRAIN_MODEL {
    container 'cfdna-classifier:latest'

    input:
    path variants

    output:
    path 'model.pkl', emit: model
    path 'metrics.json', emit: metrics

    script:
    """
    python3 /app/bin/train_model.py \\
        --input ${variants} \\
        --model-output model.pkl \\
        --metrics-output metrics.json
    """
}

What does this code do?

The input: block declares that this process expects a file path on its input channel. Nextflow automatically stages the file from the upstream process into this process's working directory.

The process produces two named outputs: model (the trained pickle) and metrics (the JSON report). Both are available as separate channels downstream.

  • Below TRAIN_MODEL, add the EVALUATE_MODEL process:
process EVALUATE_MODEL {
    container 'cfdna-classifier:latest'

    input:
    path model
    path variants

    output:
    path 'precision_recall_curve.png', emit: pr_plot
    path 'feature_importance.png', emit: fi_plot
    path 'evaluation_report.json', emit: report

    script:
    """
    python3 /app/bin/evaluate_model.py \\
        --model ${model} \\
        --data ${variants} \\
        --output-dir .
    """
}

What does this code do?

EVALUATE_MODEL takes two inputs: the trained model from TRAIN_MODEL and the original variant data from GENERATE_DATA. It produces three output files: a precision-recall plot, a feature importance plot, and a JSON evaluation report.

Each output gets its own emit label so you can publish them to different directories in the output block.

  • Below the three processes, add the workflow and output blocks that wire everything together:
workflow {
    main:
    GENERATE_DATA()
    TRAIN_MODEL(GENERATE_DATA.out.variants)
    EVALUATE_MODEL(TRAIN_MODEL.out.model, GENERATE_DATA.out.variants)

    publish:
    data = GENERATE_DATA.out.variants
    model = TRAIN_MODEL.out.model
    metrics = TRAIN_MODEL.out.metrics
    pr_plot = EVALUATE_MODEL.out.pr_plot
    fi_plot = EVALUATE_MODEL.out.fi_plot
    report = EVALUATE_MODEL.out.report
}

output {
    data {
        path 'data'
    }
    model {
        path 'model'
    }
    metrics {
        path 'metrics'
    }
    pr_plot {
        path 'plots'
    }
    fi_plot {
        path 'plots'
    }
    report {
        path 'reports'
    }
}

What does this code do?

The workflow block's main: section calls each process and passes channel outputs using the .out.label syntax. This defines the execution graph: GENERATE_DATA runs first, then TRAIN_MODEL consumes its output, then EVALUATE_MODEL consumes outputs from both.

The publish: section maps named channels to published outputs. The output block declares where each published channel writes its files (into subdirectories like data/, plots/, reports/).

  • Save main.nf.

Seeing a syntax error in your editor?

Nextflow uses Groovy-based syntax. Your editor may not recognize .nf files by default. The syntax is valid as long as you matched the structure above. You can install a Nextflow syntax extension for VS Code if you want highlighting.

Help me debug my main.nf syntax error

✔️ Awesome, I've got everything!

Great. Double check you have saved main.nf in your project root.

ⓧ I'd like to double check the full code

params.n_variants = 2000
params.seed = 42
params.outdir = 'results'

process GENERATE_DATA {
    container 'cfdna-classifier:latest'

    output:
    path 'variants.csv', emit: variants

    script:
    """
    python3 /app/bin/generate_data.py \\
        --output variants.csv \\
        --n-variants ${params.n_variants} \\
        --seed ${params.seed}
    """
}

process TRAIN_MODEL {
    container 'cfdna-classifier:latest'

    input:
    path variants

    output:
    path 'model.pkl', emit: model
    path 'metrics.json', emit: metrics

    script:
    """
    python3 /app/bin/train_model.py \\
        --input ${variants} \\
        --model-output model.pkl \\
        --metrics-output metrics.json
    """
}

process EVALUATE_MODEL {
    container 'cfdna-classifier:latest'

    input:
    path model
    path variants

    output:
    path 'precision_recall_curve.png', emit: pr_plot
    path 'feature_importance.png', emit: fi_plot
    path 'evaluation_report.json', emit: report

    script:
    """
    python3 /app/bin/evaluate_model.py \\
        --model ${model} \\
        --data ${variants} \\
        --output-dir .
    """
}

workflow {
    main:
    GENERATE_DATA()
    TRAIN_MODEL(GENERATE_DATA.out.variants)
    EVALUATE_MODEL(TRAIN_MODEL.out.model, GENERATE_DATA.out.variants)

    publish:
    data = GENERATE_DATA.out.variants
    model = TRAIN_MODEL.out.model
    metrics = TRAIN_MODEL.out.metrics
    pr_plot = EVALUATE_MODEL.out.pr_plot
    fi_plot = EVALUATE_MODEL.out.fi_plot
    report = EVALUATE_MODEL.out.report
}

output {
    data {
        path 'data'
    }
    model {
        path 'model'
    }
    metrics {
        path 'metrics'
    }
    pr_plot {
        path 'plots'
    }
    fi_plot {
        path 'plots'
    }
    report {
        path 'reports'
    }
}
Create the Nextflow configuration file

Nextflow separates pipeline logic from execution settings. The configuration file defines profiles that control where and how processes run, without changing a single line of main.nf.

  • Create a new file called nextflow.config in your project root (same folder as main.nf).
  • Paste the following profiles configuration:
profiles {
    docker_local {
        process.executor = 'local'
        docker.enabled = true
    }
    slurm_hpc {
        process.executor = 'slurm'
        process.queue = 'short'
        process.memory = '8 GB'
        process.cpus = 4
        process.time = '1h'
        process.resourceLimits = [
            memory: 750.GB,
            cpus: 200,
            time: 30.d
        ]
        singularity.enabled = true
    }
}

What does this code do?

The docker_local profile runs processes on your machine using Docker containers. This is what you will use locally.

The slurm_hpc profile configures Nextflow to submit each process as a Slurm job on an HPC cluster. It requests 4 CPUs, 8 GB of memory, and 1 hour of wall time per job. The resourceLimits cap prevents accidental over-requests. Singularity replaces Docker on most HPC systems since Docker requires root access.

Switching from local to HPC execution is a single flag: -profile slurm_hpc instead of -profile docker_local.

  • Save nextflow.config.

Unsure about the Slurm profile?

You do not need a Slurm cluster to continue. The slurm_hpc profile is configured but only activates when explicitly selected with -profile slurm_hpc. You will only use docker_local in this project.

Explain how Nextflow profiles work

✔️ Awesome, I've got everything!

Great. Double check you have saved nextflow.config in your project root.

ⓧ I'd like to double check the full code

profiles {
    docker_local {
        process.executor = 'local'
        docker.enabled = true
    }
    slurm_hpc {
        process.executor = 'slurm'
        process.queue = 'short'
        process.memory = '8 GB'
        process.cpus = 4
        process.time = '1h'
        process.resourceLimits = [
            memory: 750.GB,
            cpus: 200,
            time: 30.d
        ]
        singularity.enabled = true
    }
}
Run the pipeline locally with Docker

Everything is in place. Your main.nf defines the data flow, your nextflow.config tells Nextflow to use Docker, and your container image has the Python scripts baked in.

  • Run the full pipeline by executing this command in your project root:
nextflow run main.nf -profile docker_local

What should I see?

Nextflow prints a progress bar for each process. You should see all three processes show status COMPLETED:

  • GENERATE_DATA (1) - COMPLETED
  • TRAIN_MODEL (1) - COMPLETED
  • EVALUATE_MODEL (1) - COMPLETED

The entire run takes about 30-60 seconds depending on your machine.

  • Confirm the output files were published by running:
ls results/

You should see subdirectories data/, metrics/, model/, plots/, and reports/ containing the pipeline outputs.

Pipeline failed or processes not completing?

  • Make sure Docker Desktop is running. Nextflow needs the Docker daemon to pull and run containers.
  • Check that your Docker image is named exactly cfdna-classifier:latest by running docker images. The container name in main.nf must match.
  • Look at the error log Nextflow prints. It usually shows which process failed and the exact command that errored inside the container.

Help me debug my Nextflow pipeline failure

Your pipeline is running end-to-end. You have a fully containerized, reproducible workflow that generates data, trains a model, and evaluates results in a single command. Next up, you will take this further with multi-sample parallelization in the Secret Mission.

Secret mission

Multi-Sample Parallelization

Real bioinformatics pipelines don't process one patient at a time. In this secret mission, you'll modify your pipeline to process five patient samples concurrently using Nextflow's native channel parallelization and aggregate the results into a single summary report.

Clean Up Your Resources

Clean Up Your Resources

This project runs entirely locally with Docker, so there are no ongoing cloud costs. Decide whether to keep your resources for further experimentation or delete them to reclaim disk space.

Resources you used:

  • Docker image cfdna-classifier:latest (local, no cost).
  • Nextflow work directory work/ containing cached task outputs.
  • Nextflow cache files .nextflow/ and .nextflow.log*.
  • Project folder on your Desktop with all pipeline source files.

Keep everything running

No action needed. Choose this if you plan to keep iterating on your pipeline, experimenting with different features, or extending the Secret Mission further.

  • Your Docker image stays available for instant pipeline re-runs.
  • The work/ directory caches completed tasks, so Nextflow skips already-finished processes on the next run.
  • No cloud resources are consuming credits in the background.

Pause - I'll come back to this later

Since everything in this project is local files and a Docker image, there are no running processes to pause. Your resources sit idle on disk without consuming CPU or memory.

  • If Docker Desktop is using memory, you can quit it from the menu bar. The image will still be there next time you start Docker.
  • Your project files remain on your Desktop ready for you to pick up later.

Delete - I don't want to use this again

Remove all project resources to reclaim disk space. This deletes everything you built in this project.

  • Remove the Docker image by running:
  • Remove the Nextflow cache and work directory by running:
  • Delete the entire project folder from your Desktop.

Remove the Docker image:

docker rmi cfdna-classifier:latest

Remove the Nextflow work directory and cache files:

rm -rf work/ .nextflow/ .nextflow.log*

Delete the project folder (adjust the folder name if yours is different):

cd ~/Desktop
rm -rf cfdna-pipeline

What about Nextflow and Java?

This cleanup removes the project files and Docker image only. Nextflow and Java remain installed on your system. If you want to keep using them for other bioinformatics projects, leave them in place.

Nice Work!

Nice Work!

You just built a production-style Nextflow bioinformatics pipeline that classifies somatic variants in cell-free DNA data using machine learning. This is the exact workflow management, containerization, and applied ML intersection that bioinformatics engineering roles demand.

What you learned:

  • Built a Nextflow DSL2 pipeline with three chained processes orchestrating data generation, Random Forest model training, and precision-recall evaluation for cfDNA variant classification.
  • Containerized the entire ML environment with Docker for reproducibility, and configured switchable execution profiles (local Docker vs Slurm HPC) in nextflow.config.
  • Trained a scikit-learn classifier with class imbalance handling, 5-fold stratified cross-validation, and generated feature importance analysis and precision-recall plots as pipeline outputs.
  • Secret Mission: Parallelized the pipeline across multiple patient samples using Nextflow channels and channel.of(), aggregating per-sample metrics into a summary report.

Ready to quiz yourself?