milou User Manual¶

1. Prerequisites & Environment¶

milou requires a POSIX-compliant environment (Linux/macOS) with the following minimum software foundations:

Nextflow: Version >= 24.04.2 (tested up through 25.10.x).
Container Engine: Singularity / Apptainer (>= 3.8.0, recommended for HPC) or Docker (>= 20.10).
NVIDIA CUDA Driver & Toolkit (>= 12.0): Required only when executing the high-throughput GPU alignment track.

Nextflow JVM Memory Configuration

For large cohorts, ensure Nextflow has sufficient Java heap memory allocated in your shell environment (~/.bashrc or ~/.zshrc):

export NXF_OPTS="-Xms4g -Xmx16g"

Important Information

Almost all tools in milou run with Conda, but for GPU runs, it is recommended to use Singularity or Docker.

Conda implementation for GPU runs is not recommended.

2. Quick Start & Execution Profiles¶

2.1 Basic Execution Syntax¶

Run the pipeline using the base syntax:

nextflow run JD2112/milou \
    -r 1.2.0 \
    -profile singularity,clinical \
    --sample_sheet samplesheet.csv \
    --genome_fasta /data/genomes/GRCh38.fa \
    --diff_meth_method dss,edger \
    --outdir results/milou_clinical

Run Information

Key Arguments:

-r: Specifies the version or branch of the pipeline (e.g., -r main or -r 1.2.0).
-profile: Defines the execution environment and mode. You can combine multiple profiles using commas (e.g., singularity,gpu,clinical).
- singularity: Recommended for HPC environments.
- docker: Executes pipeline using Docker containers.
- conda: Executes pipeline using Conda environments (requires environment_main.yml).
- clinical: Enables clinical-style reporting mode.
- clinical_offline: Enables strict air-gapped clinical mode, disabling all external internet egress.
- research: Enables research mode (default).
- test_local: Executes pipeline with a minimal test dataset (6 samples).
- test_full: Executes pipeline with a test dataset mapping full differential methylation steps.
- gpu: Enables NVIDIA Parabricks acceleration.

Test the Pipeline

It is highly recommended to run the test profile on your system to ensure everything is configured properly before processing your own full datasets. The test runs require a GPU (from Parabricks) for maximum throughput and test real human genome data over varying sub-samples:

Local Test (runs quickly, uses small subset of data):

nextflow run main.nf -profile test_local,singularity,gpu,clinical

Full Pipeline Test (multi-group evaluation across 24 samples):

nextflow run main.nf -profile test_full,singularity,gpu

Running the test command will generate standard pipeline files in your working directory:

work/                               # Nextflow intermediate task executions
.nextflow.log                       # Detailed execution log
results/milou_gpu/                  # Published outputs, QC dashboards, and Quarto reports

2.2 Cohort Execution Archetypes¶

Whole-Genome Bisulfite (WGBS) - GPU

nextflow run JD2112/milou -r 1.2.0 \
    -profile gpu,singularity,clinical \
    --sample_sheet samplesheet_wgbs.csv \
    --genome_fasta GRCh38.primary_assembly.fa \
    --diff_meth_method all \
    --smoothing FALSE \
    --outdir results/wgbs_cohort

(Note: --smoothing FALSE is recommended for whole-genome cohorts to contain DSS memory under 50 GB RAM across all 28M human CpGs).

Enzymatic Methyl-Seq (EM-seq) - CPU

nextflow run JD2112/milou -r 1.2.0 \
    -profile singularity,clinical \
    --sample_sheet samplesheet_emseq.csv \
    --genome_fasta GRCh38.fa \
    --assay_type emseq \
    --diff_meth_method dss,edger \
    --outdir results/emseq_cohort

Targeted Hybridization Capture (Twist Panel)

nextflow run JD2112/milou -r 1.2.0 \
    -profile gpu,singularity,clinical \
    --sample_sheet samplesheet_twist.csv \
    --genome_fasta GRCh38.fa \
    --assay_type twist \
    --diff_meth_method all \
    --outdir results/twist_cohort

Running JD2112/milou with your data

Running the pipeline with your own data is as simple as running the base command with your:

Prepare a samplesheet,
Prepare reference files (Genome FASTA file, RefSeq BED file, and GTF file for annotation). If you are using Twist Target Region BED file, download it here.
Select workflows to run (select one or more methods).
Select a profile to run the pipeline (singularity, docker, conda, gpu, clinical, research).
Run the pipeline with the base command.

3. Sample Sheet & Reference Specifications¶

3.1 Sample Sheet Format (`samplesheet.csv`)¶

The pipeline requires a CSV samplesheet with a header row. This file maps your raw sequencing data to specific groups for comparison.

3.1.1 Required Column Schema¶

Column Header	Requirement	Description
`sample_id`	Required	Unique alphanumeric sample identifier (`A-Z`, `a-z`, `0-9`, `-`, `_`).
`group`	Required	Experimental group or clinical condition (e.g., `Control`, `Treated`). Minimum 3 replicates per group required for differential methylation.
`read1`	Required	Local filepath or remote URL to the forward paired-end FASTQ (`.fastq.gz`).
`read2`	Required	Local filepath or remote URL to the reverse paired-end FASTQ (`.fastq.gz`).
`assay_type`	Optional	Per-sample assay override: `wgbs`, `emseq`, or `twist`.

HIPAA Compliance & PII

The sample_id field is strictly validated via JSON schema to prevent the accidental inclusion of Patient Identifiable Information (PII). Spaces and special characters are forbidden. This ensures that patient names or medical record numbers are not leaked into filenames, logs, or clinical reports.

3.1.2 Minimal & Benchmark Cohort Examples¶

Minimal Example Samplesheet (DONOT RUN)

sample_id,group,read1,read2
SN09,Healthy,SN09_R1_001.fastq.gz,SN09_R2_001.fastq.gz
SN10,Disease,SN10_R1_001.fastq.gz,SN10_R2_001.fastq.gz

Full Example Samplesheet

Check out benchmarked Sample_sheet.csv

Sample sheet Information

The following data was downloaded from ENA webserver and the original data was published in Krumpolec et al (2024)

sample_id,group,read1,read2
ERR11284501,CS,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR112/001/ERR11284501/ERR11284501_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR112/001/ERR11284501/ERR11284501_2.fastq.gz
ERR11435646,VD,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/046/ERR11435646/ERR11435646_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/046/ERR11435646/ERR11435646_2.fastq.gz
ERR11435647,CS,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/047/ERR11435647/ERR11435647_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/047/ERR11435647/ERR11435647_2.fastq.gz
ERR11435651,VD,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/051/ERR11435651/ERR11435651_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/051/ERR11435651/ERR11435651_2.fastq.gz
ERR11435652,CS,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/052/ERR11435652/ERR11435652_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/052/ERR11435652/ERR11435652_2.fastq.gz
ERR11435654,VD,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/054/ERR11435654/ERR11435654_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/054/ERR11435654/ERR11435654_2.fastq.gz
ERR11435659,CS,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/059/ERR11435659/ERR11435659_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/059/ERR11435659/ERR11435659_2.fastq.gz
ERR11435661,VD,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/061/ERR11435661/ERR11435661_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/061/ERR11435661/ERR11435661_2.fastq.gz
ERR11435663,CS,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/063/ERR11435663/ERR11435663_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/063/ERR11435663/ERR11435663_2.fastq.gz
ERR11435640,VD,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/040/ERR11435640/ERR11435640_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/040/ERR11435640/ERR11435640_2.fastq.gz
ERR11435643,CS,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/043/ERR11435643/ERR11435643_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/043/ERR11435643/ERR11435643_2.fastq.gz
ERR11435644,VD,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/044/ERR11435644/ERR11435644_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/044/ERR11435644/ERR11435644_2.fastq.gz
ERR11435645,CS,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/045/ERR11435645/ERR11435645_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/045/ERR11435645/ERR11435645_2.fastq.gz
ERR11435649,VD,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/049/ERR11435649/ERR11435649_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/049/ERR11435649/ERR11435649_2.fastq.gz
ERR11435650,CS,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/050/ERR11435650/ERR11435650_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/050/ERR11435650/ERR11435650_2.fastq.gz
ERR11435653,VD,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/053/ERR11435653/ERR11435653_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/053/ERR11435653/ERR11435653_2.fastq.gz
ERR11435655,CS,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/055/ERR11435655/ERR11435655_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/055/ERR11435655/ERR11435655_2.fastq.gz
ERR11435657,VD,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/057/ERR11435657/ERR11435657_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/057/ERR11435657/ERR11435657_2.fastq.gz
ERR11435641,CS,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/041/ERR11435641/ERR11435641_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/041/ERR11435641/ERR11435641_2.fastq.gz
ERR11435642,VD,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/042/ERR11435642/ERR11435642_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/042/ERR11435642/ERR11435642_2.fastq.gz
ERR11435648,CS,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/048/ERR11435648/ERR11435648_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/048/ERR11435648/ERR11435648_2.fastq.gz
ERR11435656,VD,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/056/ERR11435656/ERR11435656_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/056/ERR11435656/ERR11435656_2.fastq.gz
ERR11435658,CS,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/058/ERR11435658/ERR11435658_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/058/ERR11435658/ERR11435658_2.fastq.gz
ERR11435660,VD,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/060/ERR11435660/ERR11435660_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/060/ERR11435660/ERR11435660_2.fastq.gz
ERR11435662,CS,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/062/ERR11435662/ERR11435662_1.fastq.gz,https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR114/062/ERR11435662/ERR11435662_2.fastq.gz

Differential Methylation Analysis Requirements

You must have minimum 3 samples per group to perform a statistically valid differential methylation analysis.

Replicating the analysis in the publication (Krumpolec et al 2024)

Download the fastq.gz files from https://www.ebi.ac.uk/ena/browser/view/ERP146869
Download TSV file with "sample_alias" (examples/filereport_read_run_ERP146869.tsv)
Run bin/map_era_samples.py to get the final csv file (examples/Sample_sheet_replicate.csv). Note- provide exact file path.
Run nextflow nextflow run main.nf -profile replicate_article,singularity,gpu

a. provide hg19.fasta file. https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000001405.13/

b. provide covered_targets_Twist_Methylome_hg19_annotated_collapsed_final.bed.zip file https://www.twistbioscience.com/resources/data-files/twist-human-methylome-panel-target-bed-file

3.2 Reference Files & Genome Builds¶

Genome FASTA (--genome_fasta): Uncompressed or .gz FASTA containing primary assembly chromosomes.
GTF Annotation (--gtf_file): Ensembl/GENCODE GTF for genomic feature mapping (promoters, exons, introns).
RefSeq BED (--refseq_file): Gene model coordinates for Gviz locus visualizations.
Target BED (--methylkit.bed_file): Required for targeted hybrid-capture panels (e.g. Twist Human Methylome).

Reference Files

The reference files should matched with the parameters configuration file in the pipeline. Check the parameters configuration file for more information.

Example: For human genome assembly GRCh38

Download the Genome Fasta file.
Download the RefSeq bed file.
Download the GTF file.
Download the Twist Target Region BED file.

4. Pipeline Architecture & Subworkflows¶

As mentioned above, milou has several subworkflows, each with its own set of parameters and outputs:

Running the pipeline with default CPU track (bismark) or GPU track (parabricks). To run with the GPU track, the user need to set the --profile gpu. By default, the pipeline will run with the CPU track.
User can choose to run the pipeline with directly with reference genome fasta file --genome_fasta [path/to/fasta] or bismark indexed file --bismark_index [path/to/index]. Default is --genome_fasta [path/to/fasta].
User also has the option to run the pipeline with aligned BAM files --aligned_bams [path/to/bams]. Default is false.
For the differential methylation analysis, the user can choose to run the pipeline with one or more methods: --diff_meth_method dss, --diff_meth_method edger or --diff_meth_method methylkit. Default is dss. If the user wants to run multiple methods, the user can provide a comma-separated list, e.g. --diff_meth_method dss,edger. They can also choose to skip the differential methylation analysis by setting the --skip_diff_meth true.

4.1 Quality Control & Adapter Trimming¶

Before alignment, the pipeline ensures your data is clean and high-quality. The pipeline uses fastqc to assess the raw sequencing quality and trim_galore to remove adapter sequences and clip low-quality base calls from the ends of reads. These are default parameters and user can adjust the --fastqc and --trim_galore arguments using --args to change the default parameters.

The pipelines also generates a post-alignments quality control report using qualimap and a multiqc report using multiqc. User can adjust the multiqc_config and multiqc_title arguments to change the default parameters. For qualimap extra arguments, check the Qualimap arguments. For multiqc extra arguments, check the MultiQC arguments.

FastQC arguments

Check this page for more information: FastQC. For fastqc extra arguments, check the FastQC arguments.

Trim Galore! arguments

Check this page for more information: Trim Galore!. For trim_galore extra arguments, check the Trim Galore! arguments.

Parameters	Description
`--multiqc_config`	multiqc config file. Default is `null`. `--multiqc_config [path/to/config]` to use multiqc config file.
`--multiqc_title`	multiqc title. Default is `null`. `--multiqc_title [title]` to use multiqc title.
`--qualimap_args`	qualimap arguments. Default is `null`. `--qualimap_args [args]` to use qualimap arguments.

4.2 Read Alignment & BAM Processing¶

Aligment is the default track used in the pipeline if not the --aligned_bams true is specified. The standard CPU track uses the Bismark suite for comprehensive bisulfite analysis.

4.2.1 Bismark Alignment (CPU Track)¶

Prepares the reference genome for bisulfite alignment by converting it in silico (C->T and G->A). The default parameter is --genome_fasta [path/to/fasta]. User can also use --bismark_index [path/to/index] to skip this step if they already have a pre-built index. Maps trimmed reads to the converted reference using Bowtie2. The default parameter is --aligner bismark. Removes PCR duplicates based on mapping position. The default process is BISMARK_DEDUPLICATE. It is essential for accurate methylation estimation, as duplicates can bias the counts. samtools is used to convert the Bismark output into sorted/indexed BAM files for downstream compatibility.

4.2.2 BWA-meth / Clara Parabricks (GPU Track)¶

Since bismark is quite slow and runs only on CPU, we recommend using the GPU track for larger datasets. User can modify the GPU runs by changing the conf/gpu.config file.

When the GPU profile is active, the pipeline automatically builds a BWA-meth index, which is required by NVIDIA Parabricks. Use --aligner bwameth to use the GPU track. A highly optimized, GPU-accelerated version of the alignment and deduplication process. The default parameter is --use_parabricks true and -profile gpu. User can also use --aligner bwameth to use the GPU track.

4.3 Methylation Calling & Extraction¶

This stage converts aligned BAM files into site-specific methylation counts.

4.3.1 Bismark Methylation Extractor (CPU Track)¶

The pipeline extracts CpG, CHG, and CHH methylation calls using the BISMARK_METHYLATION_EXTRACTOR process and generates .cov.gz files containing methylation percentages and coverage.

4.3.2 MethylDackel Extraction (GPU Track)¶

Used in the GPU track for rapid extraction of methylation metrics directly from Parabricks BAMs. The pipeline extracts CpG, CHG, and CHH methylation calls using the METHYLDACKEL_EXTRACT process and generates .cov.gz files containing methylation percentages and coverage.

Parameters	Description
`--profile`	use cpu or gpu track. Default is `cpu`. `--profile gpu` to use GPU track.
`--aligner`	aligner to use. Default is `bismark`. Change to `bwameth` to use BWA-meth.
`--use_parabricks`	use NVIDIA Parabricks. Default is `false`. `--use_parabricks true` to use NVIDIA Parabricks.
`--genome_fasta`	genome fasta file. Default is `false`. `--genome_fasta [path/to/fasta]` to use genome fasta file.
`--bismark_index`	bismark index file. Default is `false`. `--bismark_index [path/to/index]` to use bismark index file.
`--aligned_bams`	aligned bams file. Default is `false`. `--aligned_bams [path/to/bams]` to use aligned bams file.

5. Downstream Statistical Analysis¶

Once methylation calls are extracted, the pipeline performs statistical analysis to identify differentially methylated regions and genes, if not skipped by the user using --skip_diff_meth true. By default, the pipeline will run the differential methylation analysis using diff_meth_method parameter. The default parameter is --diff_meth_method dss. User can also use --diff_meth_method edger or --diff_meth_method methylkit to run the differential methylation analysis using EdgeR or MethylKit. If the user wants to run all methods, the user can set the --diff_meth_method all.

5.1 Coverage Filtering & Group Comparisons¶

Before the differential methylation analysis, CpG sites with low read depth are filtered out to ensure statistical power using --coverage_threshold (Default: 3). It also checks for the presence of samplesheet using --compare_str parameter. By default, the --compare_str all to ensure comparing different groups in the samplesheet.

For example:

If the user has a samplesheet with two groups, "control" and "treatment", the user can set the --compare_str "control,treatment" to compare the two groups. The output will be a CSV file with the differentially methylated regions and genes.
If the user has a design file with three groups, "control", "treatment1", and "treatment2", the user can set the --compare_str "control,treatment1,treatment2" to compare the three groups. The output will be a CSV file with the differentially methylated regions and genes.

Parameters	Description
`--coverage_threshold`	Minimum coverage threshold for filtering CpG sites. Default is `3`.
`--compare_str`	String of groups to compare (e.g., `control,disease`). Default is `all` to compare all groups in the samplesheet.
`--diff_meth_method`	Method to use for differential methylation analysis. Default is `dss`. Provide comma-separated list for multiple methods.
`--skip_diff_meth`	Skip differential methylation analysis. Default is `false`. `--skip_diff_meth true` to skip differential methylation analysis.

5.2 Multi-Method Consensus Scoring Framework (π-Value)¶

A central innovation of milou is its multi-method consensus scoring framework, designed to overcome the divergence between disparate differential methylation algorithms.

5.2.1 The Epigenetic Challenge¶

Distinct differential methylation callers rely on fundamentally different mathematical and statistical distributions:

DSS: Employs a Bayesian hierarchical β-binomial distribution with spatial smoothing across adjacent CpGs and empirical Bayes dispersion shrinkage.
edgeR: Models count data via generalized linear models (GLMs) with negative binomial quasi-likelihood F-tests.
methylKit: Fits logistic regression models with overdispersion correction.

Because each engine operates under distinct distributional assumptions, candidate gene lists from single tools frequently diverge, especially in cohorts with modest sequencing coverage or high biological variability.

5.2.2 Mathematical Formulation of the π-Value¶

To synthesize these orthogonal outputs into a deterministic, clinically actionable ranking, milou adapts the \(\pi\)-value framework (Xiao et al., 2014) for DNA methylation data:

\[\pi_{g} = \overline{\left| \log_{2}{(FC)}_{g} \right|} \times \left( - \log_{10}(P_{\min,g}) \right)\]

Where:

\(\overline{\left| \log_{2}{(FC)}_{g} \right|}\) is the mean absolute effect size (\(\log_2 \text{FC}\) or \(|\Delta\beta|\)) across all statistical engines that detected gene \(g\) as differentially methylated.
\(P_{\min,g}\) is the most significant (minimum) \(p\)-value observed across the callers for that gene.

5.2.3 Biological Prioritization vs. Classical Meta-Analysis¶

Standard meta-analytic methods, such as Fisher's Combined Probability Test:

\[\chi^{2} = - 2\sum_{i=1}^{k}{\ln(p_{i})}\]

rank candidate genes purely on cumulative statistical significance. In high-depth NGS cohorts, this often leads to a major clinical pitfall: biologically trivial methylation shifts (e.g., \(|\Delta\beta| < 0.05\)) can achieve extreme \(p\)-values simply due to large sample numbers or high read counts.

In contrast, the \(\pi\)-value:

Balances Statistical Rigor and Biological Magnitude: It directly weights biological effect size (\(|\Delta\beta|\)), ensuring that top-ranked biomarkers exhibit substantial epigenetic perturbation.
Suppresses Analytical Noise: Features detected by only a single engine with marginal significance are de-prioritized.
Strong Empirical Concordance: As demonstrated in our whole-genome Alzheimer's benchmark (Fetahu et al. cohort), \(\pi\)-value ranking correlates strongly with Fisher's test (\(r_s = 0.9269, p < 2.2 \times 10^{-16}\)) while correctly reprioritizing biologically critical targets (e.g., ABCA13, IMPG1) over candidates with negligible effect sizes.

5.2.4 Consensus Modes & Output Artifacts¶

--mode research (default): Retains all detected features while calculating cross-engine \(\pi\)-values and rank orders.
--mode clinical: Enforces strict majority consensus voting (requires a feature to be independently called by at least 2 out of 3 differential callers), filtering out single-caller artifacts before clinical reporting.

The resulting prioritized tables are output to: * results/unified_layer/consensus_pi_value_ranking.tsv * results/unified_layer/unified_consensus_ranking.csv * The interactive candidate tables and volcano plots in results/clinical_reporting/milou_clinical_report.html (and .pdf).

5.3 Post-Processing & Functional Enrichment Analysis¶

The pipeline performs post-processing and functional enrichment analysis on the differentially methylated regions and genes. For edgeR, it requires annotation with a GTF file (--gtf_file). For methylKit and DSS, it uses RefSeq data (--refseq_file). The post-processing steps produce a unified results layer that integrates statistical summaries of CpGs, DMRs, and genes.

5.3.1 Functional Enrichment (GO, KEGG, DisGeNET)¶

The pipeline performs advanced enrichment analysis using the clusterProfiler package:

GO Enrichment: Molecular Function, Cellular Component, and Biological Process.
KEGG Pathway Enrichment: Identification of significantly enriched biological pathways.
Disease Annotation: Automated disease association analysis using the DisGeNET database.

The pipeline also generates a Chord Diagram (using GOChord) and Dot Plots to visualize the relationship between genes and their associated pathways.

5.3.2 Key Parameters & Visualization Controls¶

Parameters	Description
`--top_n_genes`	Number of top genes to use for enrichment analysis. Default is `100`.
`--mode`	Pipeline mode: `research` (default) or `clinical` (unified reporting enabled)
`--run_clinical_report`	Force generation of Quarto clinical PDF report (default: false)
`--promoter_dist`	Distance from TSS to define a Promoter region (default: 2000)
`--enhancer_dist`	Distance from TSS to define an Enhancer/Distal region (default: 10000)

The post-processing also generates volcano plots for all comparisons. You can adjust the plotting aesthetics in conf/params.config.

Parameters	Description
`--logfc_cutoff`	Log fold change cutoff for significance. Default `0.5`.
`--pvalue_cutoff`	P-value cutoff for significance. Default `0.05`.
`--hyper_color`	Color for hypermethylated regions. Default `red`.
`--hypo_color`	Color for hypomethylated regions. Default `blue`.
`--nonsig_color`	Color for non-significant regions. Default `black`.

6. Performance, Reproducibility & Data Governance¶

6.1 Runtime Estimates & Hardware Benchmarks¶

Runtime varies based on dataset size and computational resources.

With the CPU track, processing can take several hours to days depending on the number of samples on a single node (dual Intel(R) Xeon(R) Platinum 8592+ CPU @3.90GHz). The most time-consuming step is the bismark alignment step, requiring approximately ~14–15 hours per sample. The bismark methylation calling step takes ~2–3 hours per sample.

On the GPU track (with 3 x NVIDIA L40S GPUs, 48GB VRAM), alignment takes only 15–20 minutes per sample, and a similar duration for methylation extraction. Memory ranges from <1 GB to 40 GB depending on the step.

Storage footprint is approximately 2 TB for a 24 paired-end sample cohort, including intermediate work files and final results.

Runtime Benchmarks Overview

Feature	CPU Track (Standard)	GPU Track (Parabricks)
Time/sample	16-18 hours	15-20 minutes
Hardware	12+ CPU Cores	NVIDIA GPU (16GB+ VRAM)
Memory	~10 GB	~40 GB
Storage	~2 TB	~2 TB

6.2 Bitwise Reproducibility & Statistical Determinism¶

milou is designed for Bitwise Reproducibility. Using the version tag, users can execute the pipeline with the identical codebase across any compute cluster (-r 1.2.0).

Furthermore, the pipeline enforces Statistical Determinism by injecting fixed random seeds (set.seed(42)) into all R-based modules (edgeR, methylKit, DSS, Enrichment). This guarantees that algorithms relying on stochastic sampling produce identical P-values and results across different compute environments.

6.3 Clinical Data Governance & Input Integrity¶

To meet rigorous IVDR Class C standards, the pipeline implements automated data governance:

Input Integrity (SHA256): The pipeline automatically calculates cryptographic SHA256 checksums for every input FASTQ file. Hashes are recorded in the results and printed in the Clinical Report, providing a verifiable audit trail against transfer corruption.
Output Sanity Checks: The reporting layer performs automated Pass/Fail validation, evaluating sample coverage thresholds and mapping quality. Failures are flagged prominently in the final report.
PHI Segregation: The pipeline intercepts patient identifiers right at the data ingress layer, dynamically stripping them and assigning secured pseudo-IDs (e.g., MILOU-SPEC-001).
Conversion Efficiency Gating: For bisulfite and enzymatic sequencing, the pipeline actively probes lambda phage spike-in controls. If non-conversion exceeds 1%, downstream statistical analysis halts and a diagnostic failure is flagged.
Offline Determinism: The clinical_offline profile locks the pipeline into zero-trust air-gapped execution. All disease ontologies, such as sovereign DisGeNET annotations, run entirely locally without external network calls.
Full Audit Trail: The Clinical Report embeds a complete snapshot of all resolved Nextflow parameters in JSON format alongside the execution cryptographic UUID.

6.4 Core Nextflow Arguments & JVM Memory Configuration¶

Core Nextflow Command-Line Arguments

These options are passed directly to Nextflow using a single hyphen.

-profile

Selects execution profiles (singularity, docker, conda, gpu, clinical, clinical_offline). Profiles can be combined with commas, e.g. -profile singularity,gpu,clinical.

-resume

Restarts the pipeline using cached results for any tasks whose inputs and code have not changed.

-c

Specifies a custom configuration file to override or extend pipeline settings.

-name

Assigns a custom execution name to the run for easy tracking in Nextflow Tower and local logs.

Nextflow Java virtual machines can request a large amount of memory for high-throughput cohorts. Set the following in your shell configuration (~/.bashrc or ~/.zshrc):

NXF_OPTS="-Xms4g -Xmx16g"

For a full list of flags, see the Parameters Reference.

7. Complete Parameter Reference¶

This section provides a complete reference for all command-line parameters available in milou. Default values are appropriately showcased.

To modify or adjust any parameters, please edit:

conf/params.config,
conf/resources.config for resource allocation, and
conf/gpu.config for GPU configuration.
conf/test_local.config for minimal test data run
conf/test_full.config for full test data run
conf/benchmark.config for performance tracking and resource usage benchmarking
conf/replicate_article.config for replicate article run

DAG rendering options can also be adjusted in conf/dag.config.

7.1 Input/Output Options¶

Parameter	Description	Type	Default	Required
`--sample_sheet`	Path to Samplesheet.csv with following headers: `sample_id, group, read1, read2, assay_type` (optional).	`string`	`null`	Yes
`--genome_fasta`	Path to the reference genome FASTA file.	`string`	`null`	Yes
`--assay_type`	Global chemistry override for clipping offsets (`'wgbs'`, `'emseq'`, `'twist'`).	`string`	`'wgbs'`	No
`--save_reference`	If `true`, the Bismark/BWA-meth index is saved to a persistent storeDir for reuse.	`boolean`	`false`	No
`--bismark_index`	Path to a pre-built Bismark index directory.	`string`	`false`	No
`--aligned_bams`	Start the pipeline from previously aligned BAM files instead of fastQ.	`boolean`	`false`	No
`--outdir`	Output directory where results will be stored.	`string`	`'results'`	Yes
`--design_file`	Path to a custom design matrix for specialized comparisons.	`string`	`null`	No

7.2 QC & Alignment Parameters¶

Parameter	Description	Type	Default	Required
`--aligner`	Aligner to use: `'bismark'` or `'bwameth'`.	`string`	`'bismark'`	No
`--use_parabricks`	Set to `true` to use NVIDIA Parabricks for alignment.	`boolean`	`false`	No
`--qualimap_args`	Additional arguments to pass to Qualimap.	`string`	`""`	No
`--multiqc_config`	multiqc config file. Default is `null`. `--multiqc_config [path/to/config]` to use multiqc config file.	`string`	`null`	No
`--multiqc_title`	multiqc title. Default is `null`. `--multiqc_title [title]` to use multiqc title.	`string`	`null`	No

7.3 Methylation Calling & Extraction Parameters¶

Parameter	Description	Type	Default	Required
`--coverage_threshold`	Minimum read coverage for analysis AND the threshold for automated Clinical Report Pass/Fail validation.	`integer`	`3`	Yes

7.4 Differential Methylation Analysis Parameters¶

Parameter	Description	Type	Default	Required
`--diff_meth_method`	Differential method(s) to use: `'dss'`, `'edger'`, `'methylkit'`, or a comma-separated list.	`string`	`'dss'`	Yes
`--compare_str`	The group comparison string limit (e.g., `'Healthy-Tumor'`).	`string`	`'all'`	Yes
`--smoothing`	Enable or disable moving-average spline smoothing in DSS. Set to `false` (`--smoothing FALSE`) for whole-genome WGBS cohorts (~28M CpGs) to reduce peak RAM from >250 GB to <50 GB.	`boolean`	`true`	No
`--skip_diff_meth`	Skip the differential methylation stage completely.	`boolean`	`false`	No
`--methylkit.assembly`	Assembly name for MethylKit context.	`string`	`'hg38'`	Yes
`--methylkit.diff`	Minimum methylation difference percentage for Methylkit.	`number`	`0.05`	Yes
`--methylkit.qvalue`	Maximum Q-value threshold for Methylkit.	`number`	`1`	Yes
`--methylkit.mc_cores`	Number of cores to use for MethylKit.	`integer`	`16`	Yes
`--methylkit.bed_file`	Path to the Twist Target region BED file for MethylKit context.	`string`	`null`	No

7.5 Functional Annotation & Enrichment Parameters¶

Parameter	Description	Type	Default	Required
`--gtf_file`	Path to a GTF annotation file (e.g., Gencode/Ensembl) for gene mapping.	`string`	`null`	Yes
`--refseq_file`	Path to a RefSeq BED file for annotation.	`string`	`null`	Yes
`--disgenet_db`	Path to a sovereign, local DisGeNET TSV release for clinical offline mapping.	`string`	`null`	No
`--post_processing`	Enable post-processing summaries and visualization.	`boolean`	`true`	Yes
`--logfc_cutoff`	Log2 Fold Change threshold for significance.	`number`	`0.5`	Yes
`--pvalue_cutoff`	P-value threshold for statistical significance.	`number`	`0.05`	Yes
`--top_n_genes`	Number of top genes to include in GO/KEGG enrichment analysis.	`integer`	`100`	Yes
`--hyper_color`	Color for hyper-methylated points in graphs.	`string`	`'red'`	Yes
`--hypo_color`	Color for hypo-methylated points in graphs.	`string`	`'blue'`	Yes
`--nonsig_color`	Color for non-significant points in graphs.	`string`	`'black'`	Yes

7.6 Clinical Reporting Layer Parameters¶

Parameter	Description	Type	Default	Required
`--mode`	Pipeline mode: `'research'` (standard) or `'clinical'` (enforces 2-out-of-3 consensus voting).	`string`	`'research'`	No
`--offline`	Disables internet access and API queries. Enforced automatically in `clinical_offline` profile.	`boolean`	`false`	No
`--run_clinical_report`	Force generation of automated Quarto clinical PDF/HTML report.	`boolean`	`false`	No
`--promoter_dist`	Distance from TSS (upstream) to define a Promoter region.	`integer`	`2000`	No
`--enhancer_dist`	Distance from TSS to define an Enhancer/Distal region.	`integer`	`10000`	No

7.7 Logging & Resource Management Parameters¶

Parameter	Description	Type	Default	Required
`--max_memory`	Maximum amount of RAM available for any single task.	`string`	`'128.GB'`	No
`--max_cpus`	Maximum number of CPUs available for any single task.	`integer`	`16`	No
`--max_time`	Maximum walltime for any single task.	`string`	`'240.h'`	No
`--dag.file`	Path to the output DAG file.	`string`	`'workflow_dag.dot'`	No
`--dag.overwrite`	Overwrite the output DAG file if it already exists.	`boolean`	`false`	No
`--dag.renderHTML`	Render the DAG as HTML.	`boolean`	`true`	No
`--dag.renderFormat`	Format to render the DAG.	`string`	`'png'`	No
`--dag.renderOptions`	Options to pass to the rendering tool.	`string`	`'-Tpng -Gdpi=300'`	No

7.8 Generic Pipeline Options¶

Parameter	Description	Type	Default	Required
`-profile`	Configuration profile: `docker`, `singularity`, `gpu`, `clinical`, `clinical_offline`, `benchmark`.	`string`	variable	Yes
`-resume`	Re-start the pipeline from where it left off.	`boolean`	`false`	No
`-w`	Custom working directory for intermediate files.	`string`	`'work/'`	No
`--help`	Display the pipeline help message.	`boolean`	`false`	No

8. Operational Troubleshooting¶

8.1 Memory Scaling in Whole-Genome Sequencing (WGBS)¶

Symptom: Process terminated with exit status 137 (Out of Memory) during DIFFERENTIAL_METHYLATION:DSS_ANALYSIS.
Root Cause: By default, DSS::DMLtest() computes 2D moving-average spline smoothing across adjacent genomic loci. On whole-genome human sequencing spanning all ~28–30 million CpGs, calculating continuous moving-average splines constructs large distance matrices that cause memory spikes exceeding 250 GB RAM.
Resolution: Pass --smoothing FALSE to disable the moving-average spline fitting while preserving empirical Bayes dispersion shrinkage and contiguous DMR detection via callDMR():
```
nextflow run main.nf -profile test_bisulfite_gpu,gpu,singularity \
    --smoothing FALSE
```
Bypassing spline smoothing reduces peak resident memory to ~30–47 GB RAM, allowing WGBS cohorts to run smoothly on standard 64–256 GB nodes.

8.2 Singularity / Apptainer User Namespace Collisions¶

Symptom: ERROR: Could not write info to setgroups: Permission denied or Error while waiting event for user namespace mappings.
Root Cause: Newer Apptainer versions (>=1.3+) colliding with host cluster user namespace permission restrictions.
Resolution:

Load the cluster's native Singularity-CE module:

module unload apptainer
module load singularity-ce

Set the global bind path explicitly:
```
export SINGULARITY_BINDPATH="/data"
```

8.3 Resuming Interrupted Workflows¶

Nextflow tracks all completed tasks via cryptographic input hashes. To resume an interrupted pipeline execution without re-computing upstream steps:

nextflow run JD2112/milou -profile singularity,gpu -resume