Methurator: New module (MultiQC#3447)

Author: ewels

multiqc/modules/methurator/__init__.py

@@ -0,0 +1,3 @@
+from .methurator import MultiqcModule
+
+__all__ = ["MultiqcModule"]

multiqc/modules/methurator/methurator.py

@@ -0,0 +1,368 @@
+import logging
+from typing import Dict, Any
+
+import yaml
+
+from multiqc.base_module import BaseMultiqcModule, ModuleNoSamplesFound
+from multiqc.plots import linegraph
+from multiqc.plots.linegraph import LinePlotConfig
+
+log = logging.getLogger(__name__)
+
+
+def _make_alert(message: str, level: str = "info") -> str:
+    """Create an HTML alert div for empty state messages."""
+    return f'<div class="alert alert-{level}"><strong>{message}</strong></div>'
+
+
+class MultiqcModule(BaseMultiqcModule):
+    """
+    Methurator is a Python package designed to estimate sequencing saturation for
+    reduced-representation bisulfite sequencing (RRBS) data.
+
+    The module parses the `methurator_summary.yml` file generated by the `methurator downsample`
+    command. This file contains:
+
+    - **Reads summary**: Read counts at each downsampling percentage
+    - **CpGs summary**: CpG site counts at different coverage thresholds
+    - **Saturation analysis**: Model fit results including asymptote and saturation estimates
+
+    The module displays key metrics in the General Statistics table and creates saturation
+    curve plots showing how CpG detection changes with sequencing depth.
+    """
+
+    def __init__(self):
+        super(MultiqcModule, self).__init__(
+            name="Methurator",
+            anchor="methurator",
+            href="https://github.com/VIBTOBIlab/methurator",
+            info="Estimates sequencing saturation for reduced-representation bisulfite sequencing (RRBS) data.",
+            # doi="",  # No DOI available
+        )
+
+        # Store parsed data
+        self.methurator_data: Dict[str, Dict[str, Any]] = {}
+
+        # Find and parse methurator summary files
+        for f in self.find_log_files("methurator"):
+            try:
+                parsed = yaml.load(f["f"], Loader=yaml.SafeLoader)
+            except Exception as e:
+                log.warning(f"Could not parse methurator YAML file '{f['fn']}': {e}")
+                continue
+
+            # Extract sample data from the parsed YAML
+            samples_data = self._extract_sample_data(parsed)
+            if not samples_data:
+                log.warning(f"No sample data found in {f['fn']}")
+                continue
+
+            for sample_name, sample_data in samples_data.items():
+                s_name = self.clean_s_name(sample_name, f)
+                if s_name in self.methurator_data:
+                    log.debug(f"Duplicate sample name found! Overwriting: {s_name}")
+                self.methurator_data[s_name] = sample_data
+                self.add_data_source(f, s_name=s_name)
+
+                # Add software version if available
+                version = parsed.get("methurator_summary", {}).get("metadata", {}).get("methurator_version")
+                self.add_software_version(version, s_name)
+
+        # Filter ignored samples
+        self.methurator_data = self.ignore_samples(self.methurator_data)
+
+        # Raise exception if no data found
+        if not self.methurator_data:
+            raise ModuleNoSamplesFound
+
+        log.info(f"Found {len(self.methurator_data)} samples")
+
+        # Add general stats
+        self._add_general_stats()
+
+        # Add saturation curve plot (CpGs vs reads)
+        self._add_saturation_plot()
+
+        # Add saturation percentage plot (% saturation vs reads)
+        self._add_saturation_pct_plot()
+
+        # Write data file (must be last)
+        self.write_data_file(self.methurator_data, "methurator")
+
+    def _extract_sample_data(self, parsed: Dict) -> Dict[str, Dict[str, Any]]:
+        """Extract per-sample data from the parsed YAML structure."""
+        samples_data: Dict[str, Dict[str, Any]] = {}
+
+        summary = parsed.get("methurator_summary", {})
+        reads_summary = summary.get("reads_summary", [])
+        cpgs_summary = summary.get("cpgs_summary", [])
+        saturation_analysis = summary.get("saturation_analysis", [])
+
+        # Build a mapping of sample names to their data
+        # The structure is: [{sample_name: [[pct, value], ...]}, ...]
+
+        # Parse reads summary
+        reads_by_sample: Dict[str, list] = {}
+        for item in reads_summary:
+            if isinstance(item, dict):
+                for sample_name, data in item.items():
+                    reads_by_sample[sample_name] = data
+
+        # Parse CpGs summary
+        cpgs_by_sample: Dict[str, Dict[int, list]] = {}
+        for item in cpgs_summary:
+            if isinstance(item, dict):
+                for sample_name, coverage_data in item.items():
+                    cpgs_by_sample[sample_name] = {}
+                    if isinstance(coverage_data, list):
+                        for cov_entry in coverage_data:
+                            if isinstance(cov_entry, dict):
+                                min_cov = cov_entry.get("minimum_coverage")
+                                data = cov_entry.get("data", [])
+                                if min_cov is not None:
+                                    cpgs_by_sample[sample_name][min_cov] = data
+
+        # Parse saturation analysis
+        saturation_by_sample: Dict[str, Dict[int, Dict[str, Any]]] = {}
+        for item in saturation_analysis:
+            if isinstance(item, dict):
+                for sample_name, sat_data in item.items():
+                    saturation_by_sample[sample_name] = {}
+                    if isinstance(sat_data, list):
+                        for sat_entry in sat_data:
+                            if isinstance(sat_entry, dict):
+                                min_cov = sat_entry.get("minimum_coverage")
+                                if min_cov is not None:
+                                    saturation_by_sample[sample_name][min_cov] = sat_entry
+
+        # Combine all data for each sample
+        all_samples = set(reads_by_sample.keys()) | set(cpgs_by_sample.keys()) | set(saturation_by_sample.keys())
+
+        for sample_name in all_samples:
+            sample_data: Dict[str, Any] = {
+                "reads": reads_by_sample.get(sample_name, []),
+                "cpgs": cpgs_by_sample.get(sample_name, {}),
+                "saturation_analysis": saturation_by_sample.get(sample_name, {}),
+            }
+
+            # Extract key metrics for general stats (use minimum_coverage=1 as default)
+            # Get total reads at 100% (pct=1.0)
+            for pct, reads in sample_data["reads"]:
+                if pct == 1.0:
+                    sample_data["total_reads"] = reads
+                    break
+
+            # Get CpGs and saturation metrics for the lowest minimum coverage
+            if sample_data["cpgs"]:
+                min_coverage = min(sample_data["cpgs"].keys())
+                sample_data["min_coverage_used"] = min_coverage
+
+                # Get total CpGs at 100%
+                for pct, cpgs in sample_data["cpgs"].get(min_coverage, []):
+                    if pct == 1.0:
+                        sample_data["total_cpgs"] = cpgs
+                        break
+
+                # Get saturation metrics
+                sat_data = sample_data["saturation_analysis"].get(min_coverage, {})
+                if sat_data.get("fit_success"):
+                    sample_data["asymptote"] = sat_data.get("asymptote")
+                    sample_data["saturation_pct"] = sat_data.get("saturation")
+                    sample_data["fit_success"] = True
+                else:
+                    sample_data["fit_success"] = False
+                    sample_data["fit_error"] = sat_data.get("fit_error")
+
+            samples_data[sample_name] = sample_data
+
+        return samples_data
+
+    def _add_general_stats(self):
+        """Add columns to the general statistics table."""
+        general_stats_data: Dict[str, Dict[str, Any]] = {}
+
+        for s_name, data in self.methurator_data.items():
+            general_stats_data[s_name] = {
+                "total_reads": data.get("total_reads"),
+                "total_cpgs": data.get("total_cpgs"),
+                "asymptote": data.get("asymptote"),
+                "saturation_pct": data.get("saturation_pct"),
+            }
+
+        headers = {
+            "saturation_pct": {
+                "title": "Saturation",
+                "description": "CpG sequencing saturation (percentage of theoretical maximum CpGs detected)",
+                "suffix": "%",
+                "max": 100,
+                "min": 0,
+                "scale": "RdYlGn",
+            },
+            "asymptote": {
+                "title": "Asymptote",
+                "description": "Theoretical maximum number of CpG sites (model estimate)",
+                "format": "{:,.0f}",
+                "scale": "BuPu",
+            },
+            "total_cpgs": {
+                "title": "CpGs",
+                "description": "Total number of CpG sites detected at full sequencing depth",
+                "format": "{:,.0f}",
+                "scale": "Greens",
+            },
+            "total_reads": {
+                "title": "Reads",
+                "description": "Total number of reads in the BAM file",
+                "format": "{:,.0f}",
+                "scale": "Blues",
+                "hidden": True,
+            },
+        }
+
+        self.general_stats_addcols(general_stats_data, headers)
+
+    def _add_saturation_plot(self):
+        """Add the CpG saturation curve plot showing CpGs vs number of reads."""
+        # Collect plot data for each sample, one dataset per minimum coverage level
+        plot_data_by_coverage: Dict[int, Dict[str, Dict[float, float]]] = {}
+
+        for s_name, data in self.methurator_data.items():
+            saturation_analysis = data.get("saturation_analysis", {})
+            # Build a lookup from percentage to reads count
+            reads_data = data.get("reads", [])
+            pct_to_reads = {pct: reads for pct, reads in reads_data}
+
+            for min_cov, sat_data in saturation_analysis.items():
+                if min_cov not in plot_data_by_coverage:
+                    plot_data_by_coverage[min_cov] = {}
+
+                # Extract data points: [downsampling_pct, cpgs, saturation_pct, is_extrapolated]
+                sat_points = sat_data.get("data", [])
+                sample_curve: Dict[float, float] = {}
+                for point in sat_points:
+                    if len(point) >= 2:
+                        pct = point[0]
+                        cpgs = point[1]
+                        # Map percentage to reads count if available
+                        if pct in pct_to_reads:
+                            sample_curve[float(pct_to_reads[pct])] = float(cpgs)
+
+                if sample_curve:
+                    plot_data_by_coverage[min_cov][s_name] = sample_curve
+
+        # Filter out coverage levels with no valid data and sort
+        coverage_levels = sorted([cov for cov, data in plot_data_by_coverage.items() if data])
+
+        if not coverage_levels:
+            self.add_section(
+                name="CpG Saturation Curve",
+                anchor="methurator_saturation",
+                description="Saturation curves showing the number of CpG sites detected at each sequencing depth.",
+                content=_make_alert(
+                    "No CpG saturation data available. "
+                    "No valid saturation analysis data was found for any coverage level."
+                ),
+            )
+            return
+
+        # Build datasets and labels (works for single or multiple coverage levels)
+        datasets = [plot_data_by_coverage[cov] for cov in coverage_levels]
+        data_labels = [
+            {
+                "name": f"Min Coverage {cov}x",
+                "ylab": "Number of CpGs",
+                "title": f"Methurator: CpG Saturation Curve (min. coverage {cov}x)",
+            }
+            for cov in coverage_levels
+        ]
+
+        pconfig = LinePlotConfig(
+            id="methurator_saturation_plot",
+            title=f"Methurator: CpG Saturation Curve (min. coverage {coverage_levels[0]}x)",
+            xlab="Number of Reads",
+            ylab="Number of CpGs",
+            xmin=0,
+            ymin=0,
+            tt_label="<b>{point.x:,.0f}</b> reads: {point.y:,.0f} CpGs",
+            data_labels=data_labels,
+        )
+
+        self.add_section(
+            name="CpG Saturation Curve",
+            anchor="methurator_saturation",
+            description="Saturation curves showing the number of CpG sites detected at each sequencing depth.",
+            helptext="A flattening curve indicates saturation, where additional sequencing yields diminishing returns.",
+            plot=linegraph.plot(datasets, pconfig),
+        )
+
+    def _add_saturation_pct_plot(self):
+        """Add a plot showing saturation percentage vs percentage of reads."""
+        # Collect plot data for each sample, one dataset per minimum coverage level
+        plot_data_by_coverage: Dict[int, Dict[str, Dict[float, float]]] = {}
+
+        for s_name, data in self.methurator_data.items():
+            saturation_analysis = data.get("saturation_analysis", {})
+            for min_cov, sat_data in saturation_analysis.items():
+                if min_cov not in plot_data_by_coverage:
+                    plot_data_by_coverage[min_cov] = {}
+
+                # Extract data points: [downsampling_pct, cpgs, saturation_pct, is_extrapolated]
+                sat_points = sat_data.get("data", [])
+                sample_curve: Dict[float, float] = {}
+                for point in sat_points:
+                    if len(point) >= 3 and point[2] is not None:
+                        sample_curve[point[0] * 100] = float(point[2])
+
+                if sample_curve:
+                    plot_data_by_coverage[min_cov][s_name] = sample_curve
+
+        # Filter out coverage levels with no valid data and sort
+        coverage_levels = sorted([cov for cov, data in plot_data_by_coverage.items() if data])
+
+        if not coverage_levels:
+            self.add_section(
+                name="Saturation Percentage",
+                anchor="methurator_saturation_pct",
+                description="Saturation percentage curves showing the fraction of theoretical maximum "
+                "CpG sites detected at each sequencing depth.",
+                content=_make_alert(
+                    "No saturation percentage data available. "
+                    "The saturation model failed to fit for all coverage levels, "
+                    "so saturation percentages could not be calculated. "
+                    "This typically occurs when there are too few CpG sites detected."
+                ),
+            )
+            return
+
+        # Build datasets and labels (works for single or multiple coverage levels)
+        datasets = [plot_data_by_coverage[cov] for cov in coverage_levels]
+        data_labels = [
+            {
+                "name": f"Min Coverage {cov}x",
+                "ylab": "CpG Saturation (%)",
+                "title": f"Methurator: Saturation Percentage (min. coverage {cov}x)",
+            }
+            for cov in coverage_levels
+        ]
+
+        pconfig = LinePlotConfig(
+            id="methurator_saturation_pct_plot",
+            title=f"Methurator: Saturation Percentage (min. coverage {coverage_levels[0]}x)",
+            xlab="% Sequenced Reads",
+            ylab="CpG Saturation (%)",
+            xmin=0,
+            ymin=0,
+            ymax=100,
+            tt_label="<b>{point.x:.0f}%</b> depth: {point.y:.1f}% saturation",
+            data_labels=data_labels,
+        )
+
+        self.add_section(
+            name="Saturation Percentage",
+            anchor="methurator_saturation_pct",
+            description="Saturation percentage curves showing the fraction of theoretical maximum "
+            "CpG sites detected at each sequencing depth. Only coverage levels with successful model fits are shown.",
+            helptext="100% on the x-axis represents the actual sequencing depth, with extrapolation beyond. "
+            "100% saturation on the y-axis would mean all detectable CpG sites have been found.",
+            plot=linegraph.plot(datasets, pconfig),
+        )

multiqc/search_patterns.yaml

@@ -572,6 +572,8 @@ megahit:
 metaphlan:
   fn: "*.txt"
   contents: "#clade_name\tNCBI_tax_id\trelative_abundance\t"
+methurator:
+  fn: "*methurator_summary.yml"
 methylqa:
   fn: "*.report"
   shared: true

pyproject.toml

@@ -189,6 +189,7 @@ malt = "multiqc.modules.malt:MultiqcModule"
 mapdamage = "multiqc.modules.mapdamage:MultiqcModule"
 megahit = "multiqc.modules.megahit:MultiqcModule"
 metaphlan = "multiqc.modules.metaphlan:MultiqcModule"
+methurator = "multiqc.modules.methurator:MultiqcModule"
 methylqa = "multiqc.modules.methylqa:MultiqcModule"
 mgikit = "multiqc.modules.mgikit:MultiqcModule"
 minionqc = "multiqc.modules.minionqc:MultiqcModule"