diff --git a/docker/minfi/1.48.0-1/Dockerfile b/docker/minfi/1.48.0-1/Dockerfile
new file mode 100644
index 00000000..cc376784
--- /dev/null
+++ b/docker/minfi/1.48.0-1/Dockerfile
@@ -0,0 +1,8 @@
+FROM quay.io/biocontainers/bioconductor-minfi:1.48.0--r43hdfd78af_0
+
+RUN R --no-save <<SCRIPT
+    if (!require("BiocManager", quietly = TRUE))
+        install.packages("BiocManager", repos='http://cran.us.r-project.org')
+    BiocManager::install("IlluminaHumanMethylationEPICmanifest")
+    BiocManager::install("IlluminaHumanMethylationEPICanno.ilm10b4.hg19")
+SCRIPT
diff --git a/docker/python-plotting/1.0.0/Dockerfile b/docker/python-plotting/1.0.0/Dockerfile
new file mode 100644
index 00000000..9ab805c1
--- /dev/null
+++ b/docker/python-plotting/1.0.0/Dockerfile
@@ -0,0 +1,3 @@
+FROM python:3.13.0
+
+RUN pip install pandas numpy matplotlib seaborn plotly
\ No newline at end of file
diff --git a/docker/umap/0.5.7-1/Dockerfile b/docker/umap/0.5.7-1/Dockerfile
new file mode 100644
index 00000000..058d0b3f
--- /dev/null
+++ b/docker/umap/0.5.7-1/Dockerfile
@@ -0,0 +1,3 @@
+FROM python:3.12.7
+
+RUN pip install "umap-learn==0.5.7"
\ No newline at end of file
diff --git a/workflows/methylation/methylation-cohort.wdl b/workflows/methylation/methylation-cohort.wdl
new file mode 100644
index 00000000..b34b099e
--- /dev/null
+++ b/workflows/methylation/methylation-cohort.wdl
@@ -0,0 +1,257 @@
+version 1.1
+
+workflow methylation_cohort {
+    meta {
+        description: "Process methylation data for a cohort of samples"
+        outputs: {
+            combined_beta: "Combined beta values for all samples"
+            filtered_beta: "Filtered beta values for all samples"
+            umap_embedding: "UMAP embedding for all samples"
+            umap_plot: "UMAP plot for all samples"
+        }
+    }
+
+    parameter_meta {
+        unfiltered_normalized_beta: "Array of unfiltered normalized beta values for each sample"
+    }
+
+    input {
+        Array[File] unfiltered_normalized_beta
+    }
+
+    call combine_data { input:
+        unfiltered_normalized_beta,
+    }
+
+    call filter_probes { input:
+        combined_beta_values = combine_data.combined_beta,
+    }
+
+    call generate_umap { input:
+        filtered_beta_values = filter_probes.filtered_beta_values,
+    }
+
+    call plot_umap { input:
+        umap = generate_umap.umap,
+    }
+
+    output {
+        File combined_beta = combine_data.combined_beta
+        File filtered_beta = filter_probes.filtered_beta_values
+        File umap_embedding = generate_umap.umap
+        File umap_plot = plot_umap.umap_plot
+    }
+}
+
+task combine_data {
+    meta {
+        description: "Combine data from multiple samples"
+        outputs: {
+            combined_beta: "Combined beta values for all samples"
+        }
+    }
+
+    parameter_meta {
+        unfiltered_normalized_beta: "Array of unfiltered normalized beta values for each sample"
+    }
+
+    input {
+        Array[File] unfiltered_normalized_beta
+    }
+
+    Int disk_size_gb = ceil(size(unfiltered_normalized_beta, "GiB") * 2)
+
+    command <<<
+        echo "Combining data"
+        ln -s ~{sep(" ", unfiltered_normalized_beta)} .
+
+        cat <<SCRIPT > run.py
+        import sys
+        import argparse
+        import pandas as pd
+
+        def get_args():
+            parser = argparse.ArgumentParser(
+                description="Combine CSV files.")
+            parser.add_argument(
+                "csvs", type=str, nargs="+", help="List of CSV files.")
+
+            args = parser.parse_args()
+
+            return args
+        
+        def read(filename):
+            df = pd.read_csv(filename)
+            df.rename(columns={df.columns[0]: "probe"}, inplace=True)
+            df.set_index("probe", inplace=True)
+            return df
+
+        if __name__ == "__main__":
+            args = get_args()
+
+            # Combine data
+            df = pd.concat([read(f) for f in args.csvs], axis=1, join="inner")
+            df.to_csv("combined_beta.csv", index=True)
+        SCRIPT
+
+        python run.py *.beta_swan_norm_unfiltered.csv
+
+    >>>
+
+    output {
+        File combined_beta = "combined_beta.csv"
+    }
+
+    runtime {
+        container: "quay.io/biocontainers/pandas:2.2.1"
+        memory: "8 GB"
+        cpu: 1
+        disks: "~{disk_size_gb} GB"
+        maxRetries: 1
+    }
+}
+
+task filter_probes {
+    meta {
+        description: "Filter probes based on standard deviation"
+        outputs: {
+            filtered_beta_values: "Filtered beta values for all samples"
+        }
+    }
+
+    parameter_meta {
+        combined_beta_values: "Combined beta values for all samples"
+    }
+
+    input {
+        File combined_beta_values
+    }
+
+    Int disk_size_gb = ceil(size(combined_beta_values, "GiB") * 2)
+
+    command <<<
+        ln -s ~{combined_beta_values} beta.csv
+
+        python <<SCRIPT
+        import pandas as pd
+
+        # Read beta values
+        beta = pd.read_csv("beta.csv", index_col=0)
+
+        # Calculate standard deviation for each probe
+        beta["sd"] = beta.std(axis='columns')
+
+        # Filter to only the top 10,000 probes with the highest standard deviation
+        beta.sort_values(by="sd", ascending=False).head(10000).drop(columns={"sd"}).to_csv("filtered_beta.csv")
+        SCRIPT
+    >>>
+
+    output {
+        File filtered_beta_values = "filtered_beta.csv"
+    }
+
+    runtime {
+        container: "quay.io/biocontainers/pandas:2.2.1"
+        memory: "8 GB"
+        cpu: 1
+        disks: "~{disk_size_gb} GB"
+        maxRetries: 1
+    }
+}
+
+task generate_umap {
+    meta {
+        description: "Generate UMAP embedding"
+        outputs: {
+            umap: "UMAP embedding for all samples"
+        }
+    }
+
+    parameter_meta {
+        filtered_beta_values: "Filtered beta values for all samples"
+    }
+
+    input {
+        File filtered_beta_values
+    }
+
+    Int disk_size_gb = ceil(size(filtered_beta_values, "GiB") * 2)
+
+    command <<<
+        ln -s ~{filtered_beta_values} beta.csv
+
+        python <<SCRIPT
+        import pandas as pd
+        import umap
+
+        # Read beta values
+        beta = pd.read_csv("beta.csv", index_col=0)
+
+        # Perform UMAP
+        embedding = umap.UMAP().fit_transform(beta.T)
+
+        # Save UMAP embedding
+        umap = pd.DataFrame(data=embedding, index=beta.T.index, columns=["UMAP1", "UMAP2"])
+        umap.index_name = "sample"
+        umap.to_csv("umap.csv")
+        SCRIPT
+    >>>
+
+    output {
+        File umap = "umap.csv"
+    }
+
+    runtime {
+        container: "adthrasher/umap:0.5.7-1"
+        memory: "8 GB"
+        cpu: 1
+        disks: "~{disk_size_gb} GB"
+        maxRetries: 1
+    }
+}
+
+task plot_umap {
+    meta {
+        description: "Plot UMAP embedding"
+        outputs: {
+            umap_plot: "UMAP plot for all samples"
+        }
+    }
+
+    parameter_meta {
+        umap: "UMAP embedding for all samples"
+    }
+
+    input {
+        File umap
+    }
+
+    command <<<
+        python <<SCRIPT
+        import pandas as pd
+        import matplotlib.pyplot as plt
+
+        # Read UMAP embedding
+        umap = pd.read_csv("umap.csv", index_col=0)
+
+        # Plot UMAP
+        plt.scatter(umap["UMAP1"], umap["UMAP2"])
+        plt.xlabel("UMAP1")
+        plt.ylabel("UMAP2")
+        plt.title("UMAP Embedding")
+        plt.savefig("umap.png")    
+        SCRIPT
+    >>>
+
+    output {
+        File umap_plot = "umap.png"
+    }
+
+    runtime {
+        container: "adthrasher/python-plotting:1.0.0"
+        memory: "4 GB"
+        cpu: 1
+        disks: "4 GB"
+        maxRetries: 1
+    }
+}
\ No newline at end of file
diff --git a/workflows/methylation/methylation-preprocess.wdl b/workflows/methylation/methylation-preprocess.wdl
new file mode 100644
index 00000000..5dab66e8
--- /dev/null
+++ b/workflows/methylation/methylation-preprocess.wdl
@@ -0,0 +1,164 @@
+version 1.1
+
+workflow methylation_preprocess {
+    meta {
+
+    }
+
+    parameter_meta {
+        idats: "Array of raw IDAT files from the Illumina methlyation array"
+    }
+
+    input {
+        Pair[File, File] idats
+    }
+
+    call process_raw_idats { input:
+        idats,
+    }
+
+    output {
+        File beta_swan_norm_unfiltered = process_raw_idats.beta_swan_norm_unfiltered
+        File mset = process_raw_idats.mset
+        File rgset = process_raw_idats.rgset
+        File rset = process_raw_idats.rset
+        File annotation = process_raw_idats.annotation
+        File beta = process_raw_idats.beta
+        File cn_values = process_raw_idats.cn_values
+        File m_values = process_raw_idats.m_values
+        File pheno = process_raw_idats.pheno
+        File phenoData = process_raw_idats.phenoData
+        File probe_names = process_raw_idats.probe_names
+        File sample_names = process_raw_idats.sample_names
+    }
+}
+
+task process_raw_idats {
+    meta {
+        description: "Process raw IDAT files from the Illumina methylation array for a single sample"
+        outputs: {
+
+        }
+    }
+
+    parameter_meta {
+        idats: "Array of raw IDAT files from the Illumina methlyation array"
+    }
+
+    input {
+        Pair[File, File] idats
+        Int disk_size_gb = 10
+    }
+
+    String idat_base = basename(idats.left)
+    String out_base = sub(idat_base, "(_Grn|_Red).idat", "")
+
+    command <<<
+        echo "Processing IDAT files"
+        ln -s ~{idats.left} ~{idats.right} .
+
+        R --no-save <<SCRIPT
+            options(error = function() traceback(3))
+
+            library(minfi)
+            library(dplyr)
+            library(data.table)
+            library(IlluminaHumanMethylationEPICmanifest)
+            library(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
+
+            set.seed(1)
+
+            dir=getwd()
+
+            RGSet <- read.metharray("~{idat_base}", verbose = TRUE, force = TRUE)
+            saveRDS(RGSet, "~{out_base}.RGSet.rds")
+            phenoData <- pData(RGSet)
+            write.csv(phenoData, "~{out_base}.phenoData.csv")
+
+            # The manifest is needed by preprocessRAW
+            manifest <- getManifest(RGSet)
+            manifest
+
+            # Load raw data into a MethylSet object be converting red/green channels into a matrix of methlyated and unmethylated signals.
+            MSet <- preprocessRaw(RGSet)
+            saveRDS(MSet, "~{out_base}.MSet.rds")
+
+            # Convert to a RatioSet
+            RSet <- ratioConvert(MSet, what = "both", keepCN = TRUE)
+            saveRDS(RSet, "~{out_base}.RSet.rds")
+
+            # Add genomic coordinates to each probe (plus additional annotation)
+            GRset <- mapToGenome(RSet)
+            saveRDS(GRSet, "~{out_base}.GRSet.rds")
+
+            # Take the genomic mapped RatioSet and fill Beta values (non-normalized).
+            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            #Get the NON-normalized beta values:
+            beta <- getBeta(GRset)
+            write.csv(beta, "~{out_base}.beta.csv")
+            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            
+            # Get M-value matrix and copy-number matrix
+            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            #Get M and CN vals if interested:
+            M <- getM(GRset)
+            write.csv(M, "~{out_base}.m_values.csv")
+            CN <- getCN(GRset)
+            write.csv(CN, "~{out_base}.cn_values.csv")
+            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            # Get sample names, probe names, and phenotype data
+            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            sampleNames <- sampleNames(GRset)
+            write.csv(sampleNames, "~{out_base}.sampleNames.csv")
+            probeNames <- featureNames(GRset)
+            write.csv(probeNames, "~{out_base}.probeNames.csv")
+            pheno <- pData(GRset)
+            write.csv(pheno, "~{out_base}.pheno.csv")
+
+            gr <- granges(GRset)
+            write.csv(gr, "~{out_base}.gr.csv")
+
+            annotation <- getAnnotation(GRset)
+            write.csv(annotation, "~{out_base}.annotation.csv")
+            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            #Perform SWAN normalization on beta values
+            GRset.swan_norm <- preprocessSWAN(RGSet)
+            write.csv(GRset.swan_norm, "~{out_base}.GRset.swan_norm.csv")
+            beta_swan_norm <- getBeta(GRset.swan_norm)
+            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            #Write the normalized beta-values that have NOT yet had low-variance probes filtered out
+            write.csv(beta_swan_norm, "~{out_base}.beta_swan_norm_unfiltered.csv")
+            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        SCRIPT
+    >>>
+
+    output {
+        File beta_swan_norm_unfiltered = out_base + ".beta_swan_norm_unfiltered.csv"
+        File mset = out_base + ".MSet.rds"
+        File rgset = out_base + ".RGSet.rds"
+        File rset = out_base + ".RSet.rds"
+        File annotation = out_base + ".annotation.csv"
+        File beta = out_base + ".beta.csv"
+        File cn_values = out_base + ".cn_values.csv"
+        File m_values = out_base + ".m_values.csv"
+        File pheno = out_base + ".pheno.csv"
+        File phenoData = out_base + ".phenoData.csv"
+        File probe_names = out_base + ".probeNames.csv"
+        File sample_names = out_base + ".sampleNames.csv"
+    }
+
+    runtime {
+        #container: "quay.io/biocontainers/bioconductor-minfi:1.48.0--r43hdfd78af_0"
+        container: "adthrasher/minfi:dev"
+        memory: "8 GB"
+        cpu: 1
+        disks: "~{disk_size_gb} GB"
+        maxRetries: 1
+    }
+}
\ No newline at end of file