Julian/JuliaMSI
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Compare commits

base: Julian:main
Branches Tags
Julian:main
Julian:dev/Memory_pooling_adjustments
Julian:dev/Adjustments
Julian:hotfix/bugfixes_12_03
Julian:dev/PreprocessingInUI
Julian:dev/ProcessingWithMasks
Julian:dev/AdvancedProcessingAPI
Julian:dev/FasterXMLDataAccess
..
compare: Julian:dev/AdvancedProcessingAPI
Branches Tags
Julian:main
Julian:dev/Memory_pooling_adjustments
Julian:dev/Adjustments
Julian:hotfix/bugfixes_12_03
Julian:dev/PreprocessingInUI
Julian:dev/ProcessingWithMasks
Julian:dev/AdvancedProcessingAPI
Julian:dev/FasterXMLDataAccess

No commits in common. "main" and "dev/AdvancedProcessingAPI" have entirely different histories.
main ... dev/Advanc

46 changed files with 6869 additions and 16509 deletions
Show Stats Expand all files Collapse all files

32
.github/workflows/ci.yml vendored
View File

@ -1,32 +0,0 @@
name: CI
on:
push:
branches:
- main
pull_request:
jobs:
test:
name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
runs-on: ${{ matrix.os }}
strategy:
fail-fast: false
matrix:
version:
- '1.10'
os:
- ubuntu-latest
- windows-latest
- macOS-latest
arch:
- x64
steps:
- uses: actions/checkout@v4
- uses: julia-actions/setup-julia@v2
with:
version: ${{ matrix.version }}
arch: ${{ matrix.arch }}
- uses: julia-actions/cache@v2
- uses: julia-actions/julia-buildpkg@v1
- uses: julia-actions/julia-runtest@v1

3
.gitignore vendored
View File

@ -7,6 +7,3 @@ log/*
R_original_scripts/ R_original_scripts/
test/results/ test/results/
test/binarization_results/ test/binarization_results/
Manifest.toml
MSI_sysimage.dll
MSI_sysimage.so

31
Dockerfile.headless
View File

@ -1,31 +0,0 @@
# Dockerfile.headless
# This Dockerfile creates a minimal, high-speed container solely intended
# for running JuliaMSI scripts via the pre-compiled headless sysimage.
# It explicitly excludes the Genie UI and graphical routing to optimize size and start speed.
FROM julia:1.11
# System dependencies for core MSI math libraries
RUN apt-get update && apt-get install -y \
build-essential \
zlib1g-dev \
&& rm -rf /var/lib/apt/lists/*
WORKDIR /app
# Only copy essential library files (no UI or app.jl required)
COPY Project.toml Manifest.toml ./
COPY src/ ./src/
COPY build_sysimage.jl precompile_script.jl ./
COPY test/ ./test/
# Instantiate the environment
RUN julia --project=. -e 'import Pkg; Pkg.instantiate(); Pkg.precompile()'
# Build the headless sysimage
# This command strips out Genie and plot libraries, creating a pure data-engine .so
RUN julia --project=. build_sysimage.jl --headless
# The resulting image is heavily optimized for zero-JIT execution.
# Default command simply enters a fast REPL.
CMD ["julia", "--project=.", "--threads", "auto", "-J", "sys_msi_headless.so"]

217
JuliaMSI_MainGuide.md
View File

@ -1,217 +0,0 @@
# JuliaMSI Platform User Guide
## Core Workflow and Interactive Analysis
---
### Table of Contents
1. [Introduction](#introduction)
2. [Loading a Dataset](#loading-a-dataset)
3. [Converting mzML to imzML](#converting-mzml-to-imzml)
4. [Batch Processing](#batch-processing)
5. [Spectral Visualization](#spectral-visualization)
6. [Imaging Slice Generation](#imaging-slice-generation)
7. [Topography and Surface Plots](#topography-and-surface-plots)
8. [Small Features and Tools](#small-features-and-tools)
* [Metadata Viewing](#metadata-viewing)
* [Optical Image Overlay](#optical-image-overlay)
* [Comparison of Datasets](#comparison-of-datasets)
* [Registry and Safe Processing](#registry-and-safe-processing)
9. [Mask Generation and Usage (mask.jl)](#mask-generation-and-usage-maskjl)
10. [Preprocessing](#preprocessing)
11. [Practical Tips for Analysts](#practical-tips-for-analysts)
* [Julia Version Compatibility](#julia-version-compatibility)
* [First-Time Startup](#first-time-startup)
* [Accessing the Web Interface](#accessing-the-web-interface)
* [Memory Management](#memory-management)
* [Registry and File Locations](#registry-and-file-locations)
* [Mask Handling](#mask-handling)
12. [Troubleshooting Common Issues](#troubleshooting-common-issues)
13. [Citation and Contact](#citation-and-contact)
---
## Introduction
This document serves as an entry point for analysts using the **JuliaMSI** platform. It focuses on the user-facing modules: `app.jl` (the main application entry), `julia_imzml_visual.jl` (visualization and heavy procedures), and `mask.jl` (mask generation and filtering). The goal is to explain why each step requires specific inputs and what those inputs mean scientifically.
**JuliaMSI** is a user-friendly Graphical User Interface (GUI) developed in Julia for MSI data analyses. Our framework allows for rapid **imzML** MSI data analysis and **mzML** spectra plot loading, significantly improving analytical workflows. High-resolution instruments produce massive datasets, requiring efficient memory management and parallel processing—reasons why JuliaMSI is written in Julia, which offers high performance.
### Core Objectives:
* **Simplification**: Reducing command-line dependency through intuitive visualizations and easy parameter input.
* **Accessibility**: Breaking the proprietary software barrier across fields like plant sciences, biomedical research, and environmental studies.
* **Interactive Exploration**: Using the `Plotly.jl` package for dynamic visualizations. Users can explore multiple interpretations of images, select specific coordinates for spectra analysis, or pick an *m/z* value directly from the displayed spectrum.
### Key Capabilities:
* **Smart data handling**: Processes large files in manageable sections without overwhelming memory.
* **Batch processing**: Analyzes multiple samples simultaneously for high-throughput studies.
* **Interactive exploration**: Web interface enables real-time data visualization and analysis.
---
## Loading a Dataset
**Action**: Click the **Search** button (`@onbutton btnSearch`)
Load an **imzML** or **mzML** file into the workspace, extract metadata, and prepare the dataset for visualization and processing.
The platform checks whether a pre-processed version of the same file exists in the registry (**fast load**). If not, it performs a **full load** using `OpenMSIData`. The GUI handles multiple file uploads for MSI slice generation; however, for spectra and preprocessing, it will only use the last selected file.
---
## Converting mzML to imzML
**Action**: Click **Convert/Process** (`@onbutton convert_process`)
This tool converts an **mzML** file (from a non-imaging mass spectrometry run) into an **imzML** file suitable for imaging analysis. This is necessary when spectra were acquired in a raster pattern but saved as a continuous **mzML** file without spatial coordinates.
### Inputs Required:
* **mzML file**: The raw mass spectrometry data file.
* **Sync file**: A text file that maps each spectrum in the **mzML** to a specific pixel coordinate (*x, y*). The sync file typically contains one line per spectrum with the format `x y` or similar.
> **Why these inputs?**: **mzML** files lack spatial metadata; they contain a list of spectra in acquisition order. Without coordinate information, the data cannot be reconstructed as an image. **Sync files** provide the missing spatial mapping, often generated by instrument software from a motor stage log. The conversion process embeds these coordinates into the **imzML** structure.
---
## Batch Processing
**Action**: Click **Main Process** (`@onbutton mainProcess`)
This is the core image generation routine. It takes a list of comma-separated *m/z* values and produces ion images (and optionally **TrIQ Quantized** images) for each file loaded.
### Required Inputs:
* **m/z values**: A comma-separated list of *m/z* values to visualize spatial distribution. These are the core of MSI: each image maps the abundance of a particular ion across the sample. Without a target, there is no image.
* **Tolerance**: The mass window around each target ion. Necessary because no mass spectrometer measures masses with infinite precision. The signal from a given ion is spread over a small *m/z* range; integrating over that range (by summing intensities of all data points within ±tolerance) captures the true signal.
* **Color level**: Number of discrete color levels in the output BMP image. A low number (e.g., 16) gives a coarse view; a high number (256) shows subtle intensity variations but can be noisy.
* **TrIQ**: A toggle to enable **Threshold Intensity Quantization**. This produces images with enhanced contrast and reduced influence of extreme outliers, making spatial patterns more visible. It is a contrast normalization and quantization technique.
* **TrIQ Threshold**: Probability threshold for outlier removal in TrIQ.
* **Median filter**: Replaces each pixel’s intensity with the median of its neighbors, effectively removing isolated "hot pixels" while preserving sharp boundaries.
* **Mask**: Toggle to focus on biologically relevant regions (e.g., tissue vs. background), avoiding wasted computational resources on irrelevant pixels.
> **Note**: The mask path is determined inside the JSON file, which can be manually altered, or using the specific section in `mask.jl`.
---
## Imaging Slice Generation (Part of Main Process)
In MSI, each pixel contains a full mass spectrum. To create an ion image (or "slice"), we must decide which data points in that spectrum belong to the ion of interest.
* **Tolerance**: Defines the *m/z* window. Too narrow may miss part of the peak; too wide includes noise or interfering peaks. Optimal tolerance depends on instrument resolving power and calibration stability.
* **Colorbar Generation**: Each image is accompanied by a colorbar PNG showing the mapping between intensity and colors. The bounds are determined during quantization (TrIQ or min-max scaling) and passed to `generate_colorbar_image`. This ensures consistent color interpretation with labels indicating the intensity range and scaling factors (e.g., "Intensity ×10³").
---
## Spectral Visualization
**Functions**: `createMeanPlot`, `createSumPlot`, `createXYPlot`, `createNSpectrumPlot`
Displays spectra from the loaded dataset, optionally filtered by a mask. This generates a Plotly layout plot for quality control, identification, and *m/z* selection.
### Inputs & Features:
* **Dataset**: The last **imzML** / **mzML** file loaded in the GUI.
* **Mask**: Toggle for mask filtering. The spectrum is computed only from pixels inside the mask.
* **User click**: If a user clicks on an already generated **slice** (2D image), two vertical lines appear and values are filled (one for *x*-axis, one for *y*-axis) to generate spectrum plot indicating the position of the click in the image. This allows the creation of a spectrum plot correlating to specific spatial selection.
* **NSpectrumPlot**: Allows the user to select specific spectra based on the **ID** (the index of the spectra in the dataset).
### Spectral Types:
* **Mean spectrum**: The average of all pixel spectra. Gives an overview of ions present in the entire sample and their relative abundances. Used to identify peaks for further imaging.
* **Sum spectrum**: Total ion current summed across all pixels. Emphasizes ions that are abundant across many pixels.
* **Per-pixel spectra (XY and N)**: Reveal the distribution of ions in a specific location, useful for comparing tissue regions (e.g., tumor vs. normal).
> **Why these inputs?**: Spectra are the raw data of MSI. Visualizing them is the first step: checking for unexpected peaks, assessing signal-to-noise ratio, and verifying mass calibration. Applying a mask (e.g., only tumor region) allows researchers to extract region-specific molecular signatures.
---
## Topography and Surface Plots
**Functions**: `image3dPlot`, `imageCPlot`, `triq3dPlot`, `triqCPlot`
Convert a 2D ion image into a 3D surface plot where the height (*z*-axis) represents ion intensity.
* **Inputs**: Current ion image and optional Mask.
* **Why these inputs?**: 3D surfaces and 2D Contour plots make intensity variations more tangible, especially for features with large dynamic range. They are useful for understanding the topography of ion distributions (e.g., a lipid enriched in a specific tissue layer). The 3D plot is interactive, allowing for rotation and zoom.
---
## Small Features and Tools
### Metadata Viewing (`showMetadataBtn`)
Display a summary of the dataset’s metadata: instrument settings, acquisition parameters, image dimensions, number of pixels, *m/z* range, etc.
### Optical Image Overlay (`btnOptical`)
Overlay a histological or photographic image onto the MSI image to correlate molecular distributions with tissue anatomy.
* **Inputs**: Current ion image, Transparency slider (`imgTrans`), and Optical image file (PNG, BMP, JPG, JPEG).
* **Output**: The main image tab updates to show the MSI image (or TrIQ image) blended with the optical background.
### Comparison of Datasets
Functionality to compare two datasets (e.g., control vs. treated). This involves side-by-side visualization or ratio image generation to highlight differences in molecular abundance. It allows visualization of slices with or without TrIQ filtering.
### Registry and Safe Processing
Batch processing is orchestrated by `process_file_safely`. The registry (`registry.json`) stores critical info: source path, image dimensions, summary statistics, and mask path. This file can be edited manually following its generated structure.
---
## Mask Generation and Usage (`mask.jl`)
Create binary masks that define **Regions of Interest (ROI)** in the sample. Masks filter spectra, preprocessing, and images to include only pixels inside the ROI.
### User-Defined Inputs:
* **Drawing tools**: Draw freehand shapes on the ion image or optical overlay.
* **Otsu Scale**: Automatically generates a mask based on the Otsu thresholding algorithm with the slider percentage.
* **Noise Size**: Removes noise from the mask based on the slider percentage.
* **Hole Size**: Removes holes in the mask based on the slider percentage.
* **Smoothing**: Smooths the mask based on the slider percentage.
* **Thresholding**: Set intensity thresholds to automatically segment regions.
* **Import**: Load a pre-existing mask image (e.g., from a segmentation tool).
### Outputs:
* A binary mask matrix (same dimensions as the MSI image) saved in the public/css/masks folder as a `.png` file with the name of the dataset.
* The mask is registered in the dataset’s metadata entry for automatic application in subsequent analyses.
> **Why these inputs?**: Often only a portion of the imaged area contains tissue; the rest could be background. Masking removes background pixels, improving statistical power and avoiding artifacts. Manual delineation based on optical images is standard practice in histology-directed MSI. It can also work for a selection of an area of interest (e.g. only the tumor region) in samples with multiple areas of interest.
---
## Preprocessing
To understand the deep structure of the preprocessing pipeline, please refer to the dedicated guide: **[JuliaMSI Preprocessing Guide](JuliaMSI_Preprocessing.md)**
This document describes the rationale, methods (TIC/PQN, SNIP, Wavelets), and implementation choices in detail.
---
## Practical Tips for Analysts
### Julia Version Compatibility
* The platform is tested and compatible with **Julia 1.11**.
* **Recommendation**: Use Julia 1.11 to ensure all dependencies work correctly. Use **Juliaup** to manage multiple versions easily. Downgrading your system Julia is unnecessary; you can keep multiple versions and launch the platform with the specific executable path.
### First-Time Startup
* The script will detect a missing `Manifest.toml` and automatically instantiate packages. This may take a few minutes. Subsequent startups will be much faster.
* **Caution**: Do not close the terminal/console while the server is running; it hosts the web interface.
### Accessing the Web Interface
* The platform opens your browser to `http://127.0.0.1:1481`.
* If port 1481 is in use, change this segment in the startup script: `up(host="127.0.0.1", port=1481)`.
### Memory Management
* MSI datasets can be large (gigabytes). JuliaMSI includes automatic **garbage collection** and, on Linux, calls `malloc_trim` to return freed memory to the OS.
* **OOM Troubleshooting**: Process fewer *m/z* values at once, use a mask to exclude background, or increase system swap space.
### Registry and File Locations
* `registry.json` stores paths to your datasets. If you move original **imzML** files, update the paths in the JSON or delete the registry file to let the platform rebuild it.
* Processed images and colorbars are saved in the `public/` folder, organized by dataset name.
### Mask Handling
* Ensure manually created masks have the **exact same dimensions** (width × height) as the MSI image.
* Masks are automatically applied during batch processing if the "Mask" toggle is enabled. If you move the mask file, you must update its path in the registry.
---
## Troubleshooting Common Issues
* **"Error loading active file"**: Check that the file is a valid **imzML** or **mzML** and that the corresponding `.ibd` file is in the same directory.
* **Images appear all black**: This can happen if the tolerance is too narrow (no signal captured) or if the mask is incorrectly applied. Verify *m/z* and tolerance.
* **Slow performance**: Reduce the number of *m/z* values or disable the median filter. Ensure you have enough RAM.
* **Browser interface not responding**: Refresh the page. Check the console for errors; you may need to restart the Julia process.
---
## Citation and Contact
Any revision to the code, bug, or error should be reported to: **julian.sierrag@icloud.com**
*Note: This guide is based on collective knowledge of these algorithms; please perform your own research and cite properly. To cite the JuliaMSI platform, please use:*
> José Julián Sierra-Álvarez, Martín Orlando Camargo-Escalante, Carlos Daniel Sierra-Álvarez, Carmelo Hernández-Caricio, Juan Francisco Moreno-Luna, Isabel Buendía-Corona, Robert Winkler,
> **JuliaMSI: A high-performance graphical platform for mass spectrometry imaging data analysis**,
> *Analytica Chimica Acta*, Volume 1377, 2025, 344613, ISSN 0003-2670,
> [https://doi.org/10.1016/j.aca.2025.344613](https://doi.org/10.1016/j.aca.2025.344613)

354
JuliaMSI_Preprocessing.md
View File

@ -1,354 +0,0 @@
# A Comprehensive Guide to Mass Spectrometry Imaging Preprocessing
## Rationale, Methods, and Implementation
---
### Table of Contents
1. [Introduction](#introduction)
2. [Profile vs. Centroid Mode Data](#profile-vs-centroid-mode-data)
3. [General Principles of MSI Preprocessing](#general-principles-of-msi-preprocessing)
4. [Step-by-Step Breakdown of Preprocessing Steps](#step-by-step-breakdown-of-preprocessing-steps)
* [Validation and Quality Control](#validation-and-quality-control)
* [Intensity Transformation (Stabilization)](#intensity-transformation-stabilization)
* [Smoothing](#smoothing)
* [Baseline Correction](#baseline-correction)
* [Normalization](#normalization)
* [Peak Detection](#peak-detection)
* [Peak Selection (Filtering)](#peak-selection-filtering)
* [Calibration](#calibration)
* [Peak Alignment](#peak-alignment)
* [Peak Binning and Feature Matrix Generation](#peak-binning-and-feature-matrix-generation)
5. [Example Walkthrough: From Raw Spectra to Feature Matrix](#example-walkthrough-from-raw-spectra-to-feature-matrix)
* [Phase 1 – Instrument & Data Characteristics](#phase-1--instrument--data-characteristics)
* [Phase 2 – Noise & Signal Quality](#phase-2--noise--signal-quality)
* [Phase 3 – Mass Accuracy Analysis](#phase-3--mass-accuracy-analysis)
* [Phase 4 – Spatial Region Analysis](#phase-4--spatial-region-analysis)
* [Phase 5 – Peak Characteristics](#phase-5--peak-characteristics)
* [Phase 6 – Preprocessing Recommendations](#phase-6--preprocessing-recommendations)
6. [Extra Tips for the Analyst](#extra-tips-for-the-analyst)
* [Isotopic and Charge State Matters](#isotopic-and-charge-state-matters)
* [Always Visualize the Average Spectrum First](#always-visualize-the-average-spectrum-first)
* [Preprocessing Is Not a One-Shot Process](#preprocessing-is-not-a-one-shot-process)
---
## Introduction
Mass Spectrometry Imaging (MSI) generates thousands of mass spectra, each acquired from a distinct pixel location. In the literature, before any interpretation can be drawn from these data, a series of preprocessing steps must be applied. These steps correct for instrumental artifacts, reduce noise, and align data across pixels to produce a clean, comparable feature matrix suitable for statistical analysis and image reconstruction.
Inside the **JuliaMSI** code (`Preprocessing.jl` and `PreprocessingPipeline.jl`), we have a preprocessing workflow inspired by established R packages like **MALDIquant** and **Cardinal**.
This document explains the rationale behind each step, the required inputs, the scientific meaning of parameters, and the context behind key algorithms such as **TIC normalization**, **PQN**, **wavelet-based peak detection**, and **LOWESS alignment**.
---
## Profile vs. Centroid Mode Data
Mass spectrometry data can be acquired in two primary modes: **profile** and **centroid**. Understanding the difference is crucial because preprocessing steps must be adapted accordingly.
* **Profile mode**:
The instrument records intensity at every digitized point along the *m/z* axis, producing a continuous waveform. Peaks appear as Gaussian-like shapes. This mode retains full peak shape information, which is valuable for accurate mass measurement, resolution assessment, and advanced peak detection. However, file sizes are large and noise is present between peaks.
* **Centroid mode**:
The instrument performs real-time peak detection and reduces each peak to a single centroid (*m/z*, intensity) pair, discarding the underlying profile. This dramatically reduces file size but loses peak shape details. Centroid data are essentially a list of detected peaks with no information about peak width or baseline.
In **profile mode**, algorithms must operate on the full *m/z* and intensity vectors, performing tasks like smoothing, baseline correction, and peak picking on continuous data. In **centroid mode**, these steps are often unnecessary or simplified: smoothing and baseline correction are not applicable because the data are already peak lists; peak picking is replaced by filtering based on intensity thresholds.
> [!IMPORTANT]
> The JuliaMSI code internally distinguishes this difference in some steps by acquiring the mode from the metadata on the files, but it's up to the analyst to decide skipping smoothing/baseline steps for centroid data, for example.
---
## General Principles of MSI Preprocessing
All preprocessing steps share a common goal: **to end up with an aligned feature matrix where rows correspond to pixels (spectra) and columns correspond to *m/z* bins (features).**
The steps are applied sequentially, often in a specific order to avoid compounding errors. For example, **baseline correction should precede normalization** because baseline offsets can distort total ion current calculations. Similarly, **peak alignment should follow peak detection** because alignment requires a set of common peaks across spectra.
---
## Step-by-Step Breakdown of Preprocessing Steps
### Validation and Quality Control
*(This step doesn't require analyst input)*
* **Functions in code**: `validate_spectrum`, `qc_is_empty`, `qc_is_regular`
* **Purpose**: Ensure that each input spectrum is well-formed and contains usable data before any processing.
* **Inputs required**:
* `mz::AbstractVector{<:Real}`: The *m/z* values.
* `intensity::AbstractVector{<:Real}`: The corresponding intensities.
> **Why these inputs?**
> Every subsequent algorithm depends on the integrity of these vectors. Checking for non-emptiness, equal length, finite values, non-negativity, and monotonic *m/z* order prevents cryptic errors later.
---
### Intensity Transformation (Stabilization)
* **Function**: `transform_intensity_core`
* **Purpose**: Apply a variance-stabilizing transformation to make the data more homoscedastic (constant variance across intensity range).
* **Inputs required**:
* `intensity::AbstractVector{<:Real}`
* `method::Symbol` (e.g., `:sqrt`, `:log1p`)
> **Why these inputs?**
> The transformation operates on the intensity vector only; *m/z* values remain unchanged. The choice of method depends on the data characteristics (e.g., presence of zeros).
>
> Stabilization uses functions like **square root**, **log**, **log2**, **log10** and **log1p** helps the analyst stabilize variance for count data, or handle quadratical variance. The analyst must know which one of these methods to use according to the data they are visualizing (e.g., by the mean spectrum plot, or a plot from a specific coordinate).
>
> In MSI in general, **square root stabilization** is common because ion counting statistics approximate Poisson.
>
> **Example**: A spectrum with intense peaks and low baseline noise: after square root, the dynamic range is compressed, making low-intensity features more comparable to high-intensity ones.
---
### Smoothing
* **Function**: `smooth_spectrum_core`, `moving_average_smooth`
* **Purpose**: Reduce high-frequency noise while preserving peak shapes.
* **Inputs required**:
* `y::AbstractVector{<:Real}`: Intensity vector.
* `method::Symbol` (`:savitzky_golay`, `:moving_average`)
* `window::Int`: Size of the smoothing window.
* `order::Int`: Polynomial order for Savitzky–Golay.
> **Why these inputs?**
> Smoothing operates on the intensity vector only. The window size must be an odd integer for symmetry; the polynomial order must be less than the window. These constraints ensure the filter behaves properly (e.g., Savitzky–Golay requires the window bigger than the order parameter).
>
> **Savitzky–Golay** fits a low-degree polynomial to the local window, preserving moments of the peak (height, width) better than a simple moving average, which can distort peak shapes.
>
> **Example**: A noisy peak with 5% random noise: after Savitzky–Golay, the peak becomes smoother while its centroid and height remain accurate.
---
### Baseline Correction
* **Function**: `apply_baseline_correction_core`, `_snip_baseline_impl`, `convex_hull_baseline`, `median_baseline`
* **Purpose**: Estimate and subtract the background signal (baseline) from each spectrum.
* **Inputs required**:
* `y::AbstractVector{<:Real}`: Intensity vector.
* `method::Symbol` (`:snip`, `:convex_hull`, `:median`)
* `iterations::Int` (for SNIP)
* `window::Int` (for median)
> **Why these inputs?**
> Baseline correction requires only the intensity vector. The method-specific parameters control the aggressiveness and shape of the baseline estimate.
>
> * **SNIP (Sensitive Nonlinear Iterative Peak clipping)**: Iteratively replaces each point with the minimum of itself and the average of its neighbors. After many iterations, only the baseline remains. SNIP is the default in **MALDIquant** due to its effectiveness.
> * **Convex hull**: Constructs the upper convex hull of the spectrum (or lower hull for baseline) and interpolates.
> * **Moving median**: Estimates baseline as the local median. Simple but can be biased by dense peaks.
>
> **Example**: In a spectrum with a rising baseline after SNIP correction, the baseline is removed, revealing peaks sitting on a flat background.
---
### Normalization
* **Function**: `apply_normalization_core`, `tic_normalize`, `pqn_normalize`, `median_normalize`, `rms_normalize`
* **Purpose**: Correct for differences in total ion abundance between spectra, making intensities comparable across pixels.
* **Inputs required**:
* `y::AbstractVector{<:Real}` (per-spectrum) or a matrix `M` (for PQN).
* `method::Symbol`
> **Why these inputs?**
> Normalization operates on intensities. **PQN** requires a matrix of spectra (rows = *m/z* bins, columns = spectra) because it computes a reference spectrum across all samples.
>
> * **TIC (Total Ion Current)**: Divides each intensity by the sum of all intensities in the spectrum. A simple algorithm but sensitive to a few very intense peaks (e.g., matrix peaks).
> * **Median normalization**: Divides by the median intensity. Robust to outliers but can be unstable if many intensities are zero.
> * **RMS normalization**: Divides by root-mean-square, preserving the shape of the spectrum while scaling to unit energy.
> * **PQN (Probabilistic Quotient Normalization)**: Calculates a reference spectrum and scales each spectrum by the median of the quotients (intensity/reference). This method is preferred when some peaks vary drastically because it downweights their influence.
>
> **Example**: Two spectra from different tissue regions: one has a huge lipid peak, the other has lower overall signal. TIC would suppress the lipid peak in the first spectrum, possibly losing biological information. PQN would scale both spectra more equitably.
---
### Peak Detection
* **Functions**: `detect_peaks_profile_core`, `detect_peaks_wavelet_core`, `detect_peaks_centroid_core`
* **Purpose**: Identify peaks (signals of interest) in each spectrum.
* **Inputs required**:
* *m/z*, intensity vectors.
* Method-specific parameters: `snr_threshold`, `half_window`, `min_peak_prominence`, `merge_peaks_tolerance`, etc.
> **Why these inputs?**
> Peak detection uses both *m/z* and intensity to locate local maxima and evaluate peak quality. Parameters control sensitivity and filtering. A peak is a local maximum that rises significantly above the local noise. Detection algorithms vary by data mode:
>
> * **Profile mode (local maxima)**: Scans for points that are highest in a sliding window (`half_window`). Then filters by **SNR** (signal-to-noise ratio) and prominence (height relative to surrounding minima). Close peaks are merged if within `merge_peaks_tolerance`.
> * **Wavelet-based**: Uses **Continuous Wavelet Transform (CWT)** to identify peaks at multiple scales, making it robust to varying peak widths. Local maxima in the CWT magnitude indicate peak positions. Wavelet is often used for complex spectra with overlapping peaks. CWT-based peak detection is implemented in **MALDIquant** as the default method for profile data.
> * **Centroid mode**: Simply filters the existing peak list by SNR, as peaks are already identified.
>
> **Example**: A spectrum with a small peak on the shoulder of a large peak: wavelet detection may resolve both, while simple local maxima might miss the shoulder.
---
### Peak Selection (Filtering)
* **Function**: `apply_peak_selection`
* **Purpose**: Remove peaks that do not meet quality criteria, such as minimum SNR, acceptable width, or good Gaussian shape.
* **Inputs required**:
* peaks list (from peak detection).
* Thresholds: `min_snr`, `min_fwhm_ppm`, `max_fwhm_ppm`, `min_shape_r2`.
> **Why these inputs?**
> These thresholds reflect the analyst's confidence in peak quality. They are typically derived from instrument specifications or empirical data.
>
> Not all local maxima are true peaks. Noise spikes, baseline artifacts, or electronic interference can produce spurious peaks. Filtering by **SNR** removes noise. **FWHM** (full width at half maximum) in **ppm** ensures that peaks have physically reasonable widths (e.g., isotopic peaks have predictable spacing). **Shape R²** from a Gaussian fit assesses how well the peak resembles a theoretical peak; poor shape may indicate saturation or interference.
>
> **Example**: A peak with SNR < 3 is likely noise and is discarded. A peak with FWHM > 100 *ppm* might be background rather than a specific ion.
---
### Calibration
*(This step requires that you have a list of Internal Standards inside JuliaMSI preprocessing tab filled in)*
* **Function**: `calibrate_spectra_core`, `find_calibration_peaks_core`
* **Purpose**: Correct systematic mass errors by aligning the *m/z* axis to known reference masses (internal standards).
* **Inputs required**:
* spectra (list of `MutableSpectrum`).
* `internal_standards`: dictionary of theoretical *m/z* values.
* `ppm_tolerance`: matching tolerance.
* `fit_order`: polynomial order for calibration curve.
> **Why these inputs?**
> Calibration requires a set of known masses that are present in the spectra. The algorithm matches detected peaks to these references, then fits a transformation (e.g., linear) to map measured to theoretical *m/z*.
>
> Mass spectrometers drift; this step is made to correct this drift, ensuring that the same *m/z* value corresponds to the same ion across the entire imaging dataset. **Linear or quadratic calibration** is standard. The code currently uses linear interpolation between matched peaks, which is sufficient for most MSI applications where drift is small.
>
> **Example**: A reference peak expected at **760.585 Da** (a lipid) is found at **760.590 Da** in a spectrum. Calibration shifts the entire axis so that the peak aligns with the theoretical value.
---
### Peak Alignment
* **Function**: `align_peaks_lowess_core`
* **Purpose**: Correct residual *m/z* shifts between spectra after calibration, ensuring that the same analyte peak appears at the same *m/z* across all pixels.
* **Inputs required**:
* Reference spectrum's peak *m/z* list.
* Target spectrum's peak *m/z* list.
* Parameters: method (`:lowess`, `:linear`, `:ransac`), `tolerance`, `tolerance_unit`, `max_shift_ppm`.
> **Why these inputs?**
> Alignment operates on peak lists, not full spectra, to be efficient. It finds common peaks between a reference spectrum (e.g., the one with most peaks) and each target, then fits a warping function.
>
> Even after calibration, small nonlinear shifts may persist due to sample topography or matrix effects. Alignment brings all spectra onto a common *m/z* scale, which is essential for constructing a feature matrix.
>
> **LOWESS (Locally Weighted Scatterplot Smoothing)**: Fits a smooth curve through the matched peak pairs, allowing for nonlinear corrections. It is robust to outliers and does not assume a global polynomial shape. It has been used in **MALDIquant**. **RANSAC (Random Sample Consensus)** is a robust alternative. The code defaults to LOWESS because it is well-suited for smooth, continuous drift.
>
> **Example**: Peaks in a target spectrum are systematically shifted by **+0.01 Da** at low *m/z* and **+0.005 Da** at high *m/z*. LOWESS fits a curve that corrects each peak individually, producing aligned *m/z* values.
---
### Peak Binning and Feature Matrix Generation
* **Function**: `bin_peaks_core`
* **Purpose**: Group peaks from all spectra into common *m/z* bins, producing a final feature matrix.
* **Inputs required**:
* spectra (with detected peaks).
* `params::PeakBinning`: method (`:adaptive` or `:uniform`), `tolerance`, `frequency_threshold`, etc.
> **Why these inputs?**
> Binning requires the list of peaks from each spectrum. The method determines how bins are defined.
>
> * **Adaptive binning**: Clusters peaks based on their *m/z* proximity (using tolerance). Bins are created where peaks are dense, and bins without enough peaks are discarded.
> * **Uniform binning**: Divides the total *m/z* range into a fixed number of equally spaced bins.
>
> After alignment, peaks from different spectra still have slightly different *m/z* values due to residual noise. Binning assigns them to discrete bins, creating a matrix suitable for multivariate analysis. The **frequency threshold** removes bins that appear in too few spectra (likely noise).
>
> **Example**: Peaks around *m/z* **760.58** from 100 spectra are grouped into one bin. The intensity for each spectrum in that bin is the maximum intensity of peaks assigned to it.
---
## Example Walkthrough: From Raw Spectra to Feature Matrix
### with Intelligent Parameter Recommendations
To illustrate the entire pipeline, consider a mass spectrometry imaging dataset of a mouse brain section acquired in **profile mode**. The analyst has a list of internal standards for lipids (e.g., known phospholipids at *m/z* **760.585**, **798.540**, **834.528**) to aid in mass calibration.
Before running the full preprocessing, the **JuliaMSI** code executes the intelligent pre-analysis pipeline (`main_precalculation` which is a wrapper for `run_preprocessing_analysis`). This function samples a representative subset of spectra (e.g., 100 pixels) and performs a battery of diagnostic tests to suggest optimal parameters.
#### The pre-analysis yields a report:
* **Phase 1 – Instrument & Data Characteristics**:
* Acquisition mode: **profile** (continuous waveforms)
* *m/z* axis type: **regular** (constant step of ~0.01 Da)
* Dynamic range: ~4 orders of magnitude
* **Phase 2 – Noise & Signal Quality**:
* Estimated noise level (median absolute deviation): 45 counts
* Suggested SNR threshold: **3.0**
* TIC coefficient of variation: **35%** → indicates moderate heterogeneity, so robust normalization (e.g., **PQN**) may be beneficial.
* **Phase 3 – Mass Accuracy Analysis** (using the provided lipid standards):
* Mean PPM error: **12.3 ppm**
* Standard deviation: **8.1 ppm**
* Suggested bin tolerance (mean + 3σ): **36.6 ppm** (capped between 10 and 100 ppm).
* **Phase 4 – Spatial Region Analysis**:
* *(if region masks were supplied) – skipped in this example.*
* **Phase 5 – Peak Characteristics**:
* Mean FWHM: **22.5 ppm** (converted from Δm at *m/z* 500)
* Mean Gaussian R²: **0.85** (peaks are well-shaped)
* Mean peaks per spectrum: **187**
* **Phase 6 – Preprocessing Recommendations**:
* The system synthesizes all information into a dictionary of suggested parameters for each step. Parameters that cannot be determined automatically (e.g., frequency threshold for binning, minimum matched peaks) are left as `nothing`; the analyst must supply them based on experimental goals.
> [!NOTE]
> If you have internal standards, remember to re-run the preprocessing recommendations after you have loaded them. This is only for precalculation and is not a required step for analysts that know about their data.
#### Execution of the Pipeline:
The analyst can now execute the pipeline, possibly overriding a few parameters (e.g., increasing SNR threshold to 4 for more stringent peak picking).
1. **Validation**: All spectra pass QC; one empty spectrum is flagged and removed.
2. **Intensity transformation (sqrt)**: Stabilizes variance, making low-intensity features more comparable.
3. **Baseline correction (SNIP, 100 iterations)**: Removes background using the recommended window of 200 points.
4. **Smoothing (Savitzky–Golay, window 9, order 2)**: Reduces noise while preserving peak shape.
5. **Peak detection (profile mode)**: Detects peaks with SNR ≥ 3, prominence ≥ 0.002, and uses the recommended half-window of 5 points. Peaks closer than 18 *ppm* are merged.
6. **Peak selection**: Filters peaks: keeps only those with SNR ≥ 3, FWHM between 11 and 55 *ppm*, and shape R² ≥ 0.68.
7. **Calibration**: Using the internal standards, each spectrum’s *m/z* axis is linearly corrected based on matched peaks within 36.6 *ppm* tolerance.
8. **Peak alignment (LOWESS)**: Aligns all spectra to a reference spectrum (the one with most peaks) with a span of 0.6, using the same 36.6 *ppm* tolerance. This corrects residual nonlinear drift.
9. **Normalization (RMS)**: Divides each spectrum by its root-mean-square intensity, chosen because TIC variation was high.
10. **Peak binning (adaptive)**: Groups peaks within 36.6 *ppm* into bins. Bins containing fewer than 3 peaks or appearing in <20% of spectra are discarded. The final bin centers are intensity-weighted averages.
11. **Feature matrix**: A **999 × 425 matrix** is produced, where rows are pixels and columns are *m/z* bins. This matrix is ready for PCA, segmentation, or statistical testing.
#### How the Analyst Uses the Recommendations:
* The suggested values are data-driven, but the analyst can override them based on domain knowledge. Always verify; experimenting is encouraged as data integrity is kept.
* Parameters like `frequency_threshold` in peak selection and binning require the analyst’s judgment. A threshold of **0.2** means that a bin must contain a peak in at least 20% of all pixels to be retained. This reduces noise but may discard rare but real signals.
* This is an iterative refinement; the analyst may rerun the analysis with adjusted parameters (e.g., a lower bin tolerance).
* JuliaMSI aims for a tailored preprocessing workflow that respects both the instrument’s characteristics and the biological question.
---
## Extra Tips for the Analyst
The preprocessing pipeline is powerful, but its success ultimately depends on the analyst’s informed choices. Here are practical tips to guide you through the process, especially when working with internal standards and interpreting intermediate results.
> [!WARNING]
> This is a really resource-intensive pipeline dependent on the sample, so run your precautions.
### Isotopic and Charge State Matters
When providing a dictionary of reference masses for calibration, ensure that the theoretical *m/z* values correspond to the same charge state and isotopic composition as the peaks you expect to see in your data.
* **Example**: If your instrument detects lipids mainly as **[M+H]⁺** ions, the reference mass must be the monoisotopic mass plus a proton. Using the neutral mass will introduce a systematic offset of **1.0078 Da**, leading to poor calibration. (Make the same adjustment for negative ion mode).
* If isotopic peaks are prominent, you may want to include multiple isotopes (e.g., **M**, **M+1**, **M+2**) as separate references, but be aware that their relative intensities vary. The calibration algorithm will match the closest detected peak. This is common for MSI preprocessing.
* The pre-analysis suggests a bin tolerance based on observed mass error. A tolerance that is too large may match the wrong peak (e.g., an isobaric interference); a tolerance that is too small may miss valid matches.
### Always Visualize the Average Spectrum First
Before running any preprocessing, generate a sum or average spectrum. This single plot tells you:
1. **Baseline shape**: If it's flat, rising, or undulating, it can guide the baseline correction method (**SNIP** works well for most shapes; **convex hull** is faster but may fail on curvy baselines).
2. **Noise level**: By zooming in.
3. ***m/z* Scatter**: How scattered are the spectra *m/z*.
4. **Peak density & Intensity range**: Influence the peak detection thresholds (SNR, prominence) and binning strategy.
* If the average spectrum shows a **rising baseline** from low to high *m/z*, SNIP with ~100 iterations will handle it well. If extremely wavy, increase iterations or try a different method.
* If running from the console (`run_preprocessing.jl`), the output of `run_preprocessing_analysis` provides:
* **TIC coefficient of variation (CV)**: A CV > 30% suggests substantial pixel-to-pixel variation, favoring **RMS** or **PQM** over TIC.
* **Mean PPM error**: If >20 *ppm*, calibration is essential. If <5 *ppm* and drift is linear, you might skip alignment.
* **FWHM and Gaussian R²**: Indicate peak quality. If R² < 0.6, peaks may be poorly shaped; consider increasing the `min_peak_shape_r2` threshold.
* **Acquisition mode**: For **centroid data**, never apply smoothing or baseline correction because those steps assume continuous waveforms.
* The pipeline will set `:method = :centroid` in peak picking and use permissive filters.
* Use a low SNR threshold (e.g., 2.0).
* Disable shape-based filtering (`min_peak_shape_r2 = 0.0`).
* Centroid data do not retain peak width information; set `min_fwhm_ppm = 0.0` and `max_fwhm_ppm = 500.0` to avoid discarding real peaks.
* Adaptive binning still works; tolerance should be based on expected mass accuracy.
### Preprocessing Is Not a One-Shot Process
The recommended parameters are a starting point, not a final answer.
* **Reproducibility**: Keep a record of all parameters used. The `params` dictionary can be saved as a **JSON**. This is essential for publications.
* **Biological Context**: If looking for rare metabolites appearing in few pixels, lower the `frequency_threshold` in binning. If interested only in high-abundance lipids, increase the SNR threshold.
* **Modularity**: The pipeline is modular; you can always rerun parts of it with adjusted parameters without starting from scratch.
---
*Any revision to the code, bug, error, report to julian.sierrag@icloud.com, if you we're to cite this file, remember i provide no citation and its based on internet knowledge of the web about these algorithms, so do your own research and cite properly. to cite JuliaMSI platform, please cite:*
José Julián Sierra-Álvarez, Martín Orlando Camargo-Escalante, Carlos Daniel Sierra-Álvarez, Carmelo Hernández-Caricio, Juan Francisco Moreno-Luna, Isabel Buendía-Corona, Robert Winkler,
JuliaMSI: A high-performance graphical platform for mass spectrometry imaging data analysis,
Analytica Chimica Acta,
Volume 1377,
2025,
344613,
ISSN 0003-2670,
https://doi.org/10.1016/j.aca.2025.344613

2653
Manifest.toml Normal file
View File

File diff suppressed because it is too large Load Diff

66
Project.toml
View File

@ -1,5 +1,4 @@
name = "MSI_src" name = "MSI_src"
uuid = "8f395f85-3b9c-4f7f-bf7e-12345678abcd"
authors = ["JJSA"] authors = ["JJSA"]
version = "0.1.0" version = "0.1.0"
@ -8,10 +7,8 @@ Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
CodecBase = "6c391c72-fb7b-5838-ba82-7cfb1bcfecbf"
ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
ContinuousWavelets = "96eb917e-2868-4417-9cb6-27e7ff17528f"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
@ -22,13 +19,10 @@ HistogramThresholding = "2c695a8d-9458-5d45-9878-1b8a99cf7853"
ImageBinarization = "cbc4b850-ae4b-5111-9e64-df94c024a13d" ImageBinarization = "cbc4b850-ae4b-5111-9e64-df94c024a13d"
ImageComponentAnalysis = "d9b9e9a0-1569-11e9-2cb5-bbca914b0e89" ImageComponentAnalysis = "d9b9e9a0-1569-11e9-2cb5-bbca914b0e89"
ImageContrastAdjustment = "f332f351-ec65-5f6a-b3d1-319c6670881a" ImageContrastAdjustment = "f332f351-ec65-5f6a-b3d1-319c6670881a"
ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5" ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
ImageMorphology = "787d08f9-d448-5407-9aad-5290dd7ab264" ImageMorphology = "787d08f9-d448-5407-9aad-5290dd7ab264"
ImageSegmentation = "80713f31-8817-5129-9cf8-209ff8fb23e1" ImageSegmentation = "80713f31-8817-5129-9cf8-209ff8fb23e1"
Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0" Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
Libz = "2ec943e9-cfe8-584d-b93d-64dcb6d567b7" Libz = "2ec943e9-cfe8-584d-b93d-64dcb6d567b7"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@ -36,70 +30,10 @@ Loess = "4345ca2d-374a-55d4-8d30-97f9976e7612"
Mmap = "a63ad114-7e13-5084-954f-fe012c677804" Mmap = "a63ad114-7e13-5084-954f-fe012c677804"
NativeFileDialog = "e1fe445b-aa65-4df4-81c1-2041507f0fd4" NativeFileDialog = "e1fe445b-aa65-4df4-81c1-2041507f0fd4"
NaturalSort = "c020b1a1-e9b0-503a-9c33-f039bfc54a85" NaturalSort = "c020b1a1-e9b0-503a-9c33-f039bfc54a85"
PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5" PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
SavitzkyGolay = "c4bf5708-b6a6-4fbe-bcd0-6850ed671584" SavitzkyGolay = "c4bf5708-b6a6-4fbe-bcd0-6850ed671584"
Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
StipplePlotly = "ec984513-233d-481d-95b0-a3b58b97af2b" StipplePlotly = "ec984513-233d-481d-95b0-a3b58b97af2b"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
[compat]
Accessors = "0.1"
Base64 = "1.11"
CSV = "0.10"
CairoMakie = "0.13"
CodecBase = "0.3"
ColorSchemes = "3.30"
Colors = "0.12"
ContinuousWavelets = "1.1"
DataFrames = "1.7"
Dates = "1.11"
FileIO = "1.17"
GLMakie = "0.11"
Genie = "5.31"
GenieFramework = "2.8"
HistogramThresholding = "0.3"
ImageBinarization = "0.3"
ImageComponentAnalysis = "0.2"
ImageContrastAdjustment = "0.3"
ImageCore = "0.10"
ImageFiltering = "0.7"
ImageMorphology = "0.4"
ImageSegmentation = "1.9"
Images = "0.26"
Interpolations = "0.15"
JLD2 = "0.6"
JSON = "0.21"
Libz = "1.0"
LinearAlgebra = "1.11"
Loess = "0.6"
Mmap = "1.11"
NativeFileDialog = "0.2"
NaturalSort = "1.0"
Parameters = "0.12"
PlotlyBase = "0.8"
Printf = "1.11"
ProgressMeter = "1.11"
SavitzkyGolay = "0.9"
Serialization = "1.11"
Setfield = "1.1"
Statistics = "1.11"
StatsBase = "0.34"
StipplePlotly = "0.13"
Test = "1.11.0"
UUIDs = "1.11"
julia = "1.11, 1.12"
[extras]
Test = "8dfed614-e60c-5e00-aba5-a33c2d43236e"
[targets]
test = ["Test"]

56
README.md
View File

@ -3,10 +3,9 @@
https://codeberg.org/LabABI/JuliaMSI https://codeberg.org/LabABI/JuliaMSI
## Local Installation ## Local Installation
1. Make sure you have Julia installed and at least Julia version 1.6 but not higher than 1.11, If not, follow the next installation guide which contains the latest files and instructions for [juliaup](https://github.com/JuliaLang/juliaup) (recommended) alternatively, you can follow the official Julia installation instructions at [julia](https://julialang.org/downloads/platform/). 1. Make sure you have Julia installed and at least Julia version 1.6, If not, follow the next installation guide which contains the latest files and instructions for [juliaup](https://github.com/JuliaLang/juliaup) (recommended) alternatively, you can follow the official Julia installation instructions at [julia](https://julialang.org/downloads/platform/).
2. Download and decompress the following repository (You can download the repository from [download](https://codeberg.org/LabABI/JuliaMSI/archive/main.zip) or locate the option manually at the top right of the page, in the box with [•••] (More operations). For a more advanced method using Git, use this link: https://codeberg.org/LabABI/JuliaMSI.git). 2. Download and decompress the following repository (You can download the repository from [download](https://codeberg.org/LabABI/JuliaMSI/archive/main.zip) or locate the option manually at the top right of the page, in the box with [•••] (More operations). For a more advanced method using Git, use this link: https://codeberg.org/LabABI/JuliaMSI.git).
3. Ensure that you have decompressed the ZIP file and know the working directory where it is located (the folder on your computer where you saved the JuliaMSI files). 3. Ensure that you have decompressed the ZIP file and know the working directory where it is located (the folder on your computer where you saved the JuliaMSI files).
4. Julia Version 1.11 is required. You might need to downgrade your whole Julia installation to support it.
## Load User Interface ## Load User Interface
1. Set the working directory to where JuliaMSI is located (the specific location of the folder on your computer where you decompressed the ZIP) in your terminal.<br> 1. Set the working directory to where JuliaMSI is located (the specific location of the folder on your computer where you decompressed the ZIP) in your terminal.<br>
@ -22,8 +21,8 @@ https://codeberg.org/LabABI/JuliaMSI
Example of a correct route:<br> Example of a correct route:<br>
~/Downloads/JuliaMSI-main/juliamsi ~/Downloads/JuliaMSI-main/juliamsi
2. Without entering the Julia environment, launch the project in your terminal with the following command (which works for all operating systems): 2. Without entering the Julia environment, launch the project in your terminal with the following command (which works for all operating systems):
```bash ```
julia --threads auto --project=. start_MSI_GUI.jl julia --project=. start_MSI_GUI.jl
``` ```
3. After the script has finished loading, you can open a [page](http://127.0.0.1:1481/) in your browser with the web app running. 3. After the script has finished loading, you can open a [page](http://127.0.0.1:1481/) in your browser with the web app running.
@ -34,47 +33,6 @@ Minimum system requirements: 4 core processor, 8 GB RAM<br>
JuliaMSI is a Graphical User Interface for a library of MSI tools in Julia: https://github.com/CINVESTAV-LABI/julia_mzML_imzML JuliaMSI is a Graphical User Interface for a library of MSI tools in Julia: https://github.com/CINVESTAV-LABI/julia_mzML_imzML
## Build the System Image (One-time)
Alternatively, you can generate custom pre-compiled `.so` / `.dll` system images by running the build script in your directory. This bakes the heavy mass-spectrometry kernels and UI dependencies into a single machine-code binary for ultra-fast startup times.
### 1. Build the GUI Sysimage (for full App usage)
```bash
julia --project=. build_sysimage.jl
```
This may take 5–15 minutes. The resulting `sys_msi_gui.so` (or `.dll` on Windows) file will be around **300MB–600MB** because it contains the pre-compiled machine code for your entire graphical environment.
*Note: A sysimage built on Linux (.so) will not work on Windows. You must run the build script once on each target operating system.*
Once `sys_msi_gui.so` is created in your directory, adapt your launch command:
```bash
julia --project=. -e 'using Pkg; Pkg.precompile()'
# Adjust the .so extension to .dll if on Windows or .dylib if on macOS
julia --threads auto --project=. --sysimage sys_msi_gui.so start_MSI_GUI.jl
```
### 2. Build the Headless Sysimage (for Scripts & Data Scientists)
If you only want to run scripts or the core `MSI_src` engine without the overhead of the Genie web server and Plotly, you can build a stripped-down, high-speed system image:
```bash
julia --project=. build_sysimage.jl --headless
```
Launch your automated tests or custom processing scripts using the headless image to bypass virtually all compilation overhead:
```bash
julia --threads auto --project=. --sysimage sys_msi_headless.so test/test_streaming_pipeline.jl
```
### 3. Docker Option for High-Speed Processing
For maximum portability and execution speed without clutter, we provide a `Dockerfile.headless`. This option creates a minimal, ultra-fast container that completely omits `app.jl` and the web interface.
It automatically builds and bundles the headless sysimage, giving researchers an instant-start engine ready for heavy-duty `.imzML` batch pipelines with zero JIT latency.
To build and run the Docker image:
```bash
docker build -t juliamsi-headless -f Dockerfile.headless .
# Run a script mapped from your local data folder:
docker run -it -v /my/local/data:/data juliamsi-headless julia -J sys_msi_headless.so /data/my_processing_script.jl
```
## License ## License
JuliaMSI is published under the terms of the MIT License. JuliaMSI is published under the terms of the MIT License.
@ -82,10 +40,4 @@ JuliaMSI is published under the terms of the MIT License.
## Example data ## Example data
Robert Winkler. (2023). mzML mass spectrometry and imzML mass spectrometry imaging test data [Data set]. Robert Winkler. (2023). mzML mass spectrometry and imzML mass spectrometry imaging test data [Data set].
Zenodo. <https://doi.org/10.5281/zenodo.10084132> Zenodo. <https://doi.org/10.5281/zenodo.10084132>
## User Guides
For a general guide on how to use the JuliaMSI platform, please refer to the [JuliaMSI Main Guide](JuliaMSI_MainGuide.md).
For a more detailed guide on the preprocessing pipeline, please refer to the [JuliaMSI Preprocessing Guide](JuliaMSI_Preprocessing.md).

3573
app.jl
View File

File diff suppressed because it is too large Load Diff

1565
app.jl.html
View File

File diff suppressed because it is too large Load Diff

49
build_sysimage.jl
View File

@ -1,49 +0,0 @@
# build_sysimage.jl
using PackageCompiler
headless = "--headless" in ARGS
println("--- Starting Custom System Image Build ---")
println("Mode: ", headless ? "Headless (Scripting/Docker)" : "Full GUI (Genie App)")
println("This process can take 5-15 minutes.")
ext = Sys.iswindows() ? ".dll" : Sys.isapple() ? ".dylib" : ".so"
sysimage_name = headless ? "sys_msi_headless" * ext : "sys_msi_gui" * ext
sysimage_path = joinpath(@__DIR__, sysimage_name)
packages_to_include = Symbol[
:MSI_src, # <--- Bake our core library
:DataFrames,
:SparseArrays,
:Mmap,
:Libz,
:Serialization,
:Printf,
:JSON
]
if !headless
append!(packages_to_include, [
:Genie,
:GenieFramework,
:PlotlyBase
])
end
create_sysimage(
packages_to_include,
sysimage_path = sysimage_path,
precompile_execution_file = "precompile_script.jl",
incremental = true,
filter_stdlibs = false
)
println("--- System Image Build Successful! ---")
println("Created: $(sysimage_path)")
if headless
println("To use in a headles script:")
println("julia --threads auto --project=. --sysimage $(sysimage_name) my_script.jl")
else
println("To start Julia with the GUI image, use:")
println("julia --threads auto --project=. --sysimage $(sysimage_name) start_MSI_GUI.jl")
end

6
config/env/prod.jl vendored
View File

@ -1,6 +0,0 @@
# config/env/prod.jl
# This file is loaded when GENIE_ENV is set to "prod"
# Define PLUGINS_WITH_ASSETS to prevent UndefVarError in production mode
# It's typically an empty vector if no specific plugins with assets are tracked.
const PLUGINS_WITH_ASSETS = []

1108
julia_imzML_visual.jl
View File

File diff suppressed because it is too large Load Diff

244
mask.jl
View File

@ -9,14 +9,14 @@ using Statistics, NaturalSort, LinearAlgebra, StipplePlotly
using Base.Filesystem: mv using Base.Filesystem: mv
using MSI_src using MSI_src
using .MSI_src: MSIData, process_image_pipeline, REGISTRY_LOCK using .MSI_src: MSIData
# Plot Handling # Plot Handling
if !@isdefined(increment_image) include("./julia_imzML_visual.jl")
include("./julia_imzML_visual.jl")
end
# Image Processing Pipeline # Image Processing Pipeline
include("src/ImageProcessing.jl")
using .ImageProcessing
using ImageBinarization using ImageBinarization
function load_and_binarize_mask(path) function load_and_binarize_mask(path)
@ -46,7 +46,7 @@ function alter_image(img_path, otsu_scale, noise_size_percent, hole_size_percent
gray_img = Float32.(ensure_grayscale(original_img)) gray_img = Float32.(ensure_grayscale(original_img))
binary, noise_removed, holes_filled, smoothed = binary, noise_removed, holes_filled, smoothed =
process_image_pipeline(gray_img; ImageProcessing.process_image_pipeline(gray_img;
otsu_scale=otsu_scale, noise_size_percent=noise_size_percent, otsu_scale=otsu_scale, noise_size_percent=noise_size_percent,
hole_size_percent=hole_size_percent, smoothing=smoothing_level) hole_size_percent=hole_size_percent, smoothing=smoothing_level)
@ -152,9 +152,7 @@ end
@in selected_folder_main = "" @in selected_folder_main = ""
@out available_folders = String[] @out available_folders = String[]
@out image_available_folders = String[] @out image_available_folders = String[]
@out registry_path = abspath(joinpath(@__DIR__, "public", "registry.json"))
@private registry_init_done = false @private registry_init_done = false
@in refetch_folders = false
@out imgInt = "" # Path to current slice @out imgInt = "" # Path to current slice
@out current_msi = "" @out current_msi = ""
@ -193,7 +191,6 @@ end
@in btn_flip_mask = false @in btn_flip_mask = false
@in btn_save_final_mask = false @in btn_save_final_mask = false
@in btn_move_mask = false @in btn_move_mask = false
@in btn_upload_premade_mask = false
# --- Manual Editor State --- # --- Manual Editor State ---
@in brush_size = 10 @in brush_size = 10
@ -253,6 +250,7 @@ end
is_browsing_slices = true is_browsing_slices = true
is_editing_mask = false is_editing_mask = false
folder_path = joinpath("public", selected_folder_main) folder_path = joinpath("public", selected_folder_main)
println("Selected folder: $selected_folder_main")
if !isdir(folder_path) if !isdir(folder_path)
imgInt = "" imgInt = ""
msgimg = "Folder not found." msgimg = "Folder not found."
@ -275,33 +273,6 @@ end
msgimg = "No MSI images found in this dataset." msgimg = "No MSI images found in this dataset."
show_verification_plot = false show_verification_plot = false
end end
# --- Check for mask in registry and update message ---
reg_path = abspath(joinpath(@__DIR__, "public", "registry.json"))
if isfile(reg_path)
registry = JSON.parsefile(reg_path) # Assuming JSON.jl is loaded
if haskey(registry, selected_folder_main)
dataset_entry = registry[selected_folder_main]
if get(dataset_entry, "has_mask", false) && haskey(dataset_entry, "mask_path")
mask_filepath = dataset_entry["mask_path"]
mask_filename = basename(mask_filepath)
mask_editor_message = "Dataset has an existing mask: $(mask_filename)"
# Optionally, load this mask if the user is currently editing
if is_editing_mask && isempty(smoothed_mask_path)
# Only auto-load if not already editing a mask and smoothed_mask_path is empty
smoothed_mask_path = get_timestamped_path("/css/masks/$(mask_filename)")
plotdata_verify, plotlayout_verify, show_verification_plot = update_main_plot(imgInt, smoothed_mask_path, imgTrans)
mask_editor_message *= " (loaded for editing)"
end
else
mask_editor_message = "No existing mask found for this dataset."
end
else
mask_editor_message = "Dataset entry not found in registry."
end
else
mask_editor_message = "Registry file (registry.json) not found."
end
end end
end end
@ -376,86 +347,6 @@ end
end end
end end
@onbutton btn_upload_premade_mask begin
picked_path = pick_file(; filterlist="png,bmp,jpg,jpeg")
if !isempty(picked_path)
target_dir = joinpath("public", "css", "masks")
mkpath(target_dir)
uploaded_filename = "premade_mask.png"
destination_path = joinpath(target_dir, uploaded_filename)
local_uploaded_path = "" # Define here for finally block cleanup
try
cp(picked_path, destination_path; force=true)
local_uploaded_path = destination_path # Full path to the copied image
# --- Binary Check ---
img = load(local_uploaded_path)
gray_img = ensure_grayscale(img)
# Check if image is effectively binary (only two distinct values, or close enough)
# Convert to raw pixel values for robust checking
pixel_values = Float64.(collect(Iterators.flatten(channelview(gray_img))))
unique_pixel_values = unique(pixel_values)
# Heuristic: if more than 2 distinct values exist, it's not truly binary
# Allow for floating point inaccuracies
if length(unique_pixel_values) > 2
# Find min and max for range check
min_val, max_val = extrema(unique_pixel_values)
# If values are not clustered at 0 and 1, it's likely not binary
if !(min_val < 0.1 && max_val > 0.9 && all(v -> v < 0.1 || v > 0.9, unique_pixel_values))
mask_editor_message = "Warning: Uploaded image is not purely binary (found $(length(unique_pixel_values)) distinct pixel values). Please use 'Upload to Create Mask' for processing or ensure your image is strictly black and white."
mask_editor_warning = true
rm(local_uploaded_path; force=true) # Clean up non-binary file
return
end
end
# If it has 2 values, or only one value, ensure it's a true BitMatrix
# This handles cases like 0-255 images being converted to 0.0-1.0 and then to true/false
binary_img = gray_img .>= 0.5 # Ensure pure binary (true/false)
save(local_uploaded_path, binary_img)
# --- Resizing to Slice (if imgInt is present) ---
slice_path_cleaned = replace(imgInt, r"\?.*" => "")
slice_full_path = joinpath("public", lstrip(slice_path_cleaned, '/'))
if !isempty(imgInt) && isfile(slice_full_path) # Only resize if a slice is present
slice_img = load(slice_full_path)
mask_img = load(local_uploaded_path) # Reload potentially re-binarized image to get correct type/size
if size(slice_img) != size(mask_img)
@info "Resizing pre-made mask to match slice dimensions."
resized_mask = imresize(mask_img, size(slice_img))
save(local_uploaded_path, resized_mask)
end
end
is_browsing_slices = false
is_editing_mask = true
# Update smoothed_mask_path to the new file, with timestamp for cache busting
relative_path_for_web = "/css/masks/$(uploaded_filename)"
smoothed_mask_path = get_timestamped_path(relative_path_for_web)
mask_editor_message = "Pre-made mask uploaded and loaded successfully!"
mask_editor_warning = false
# Trigger plot update
plotdata_verify, plotlayout_verify, show_verification_plot = update_main_plot(imgInt, smoothed_mask_path, imgTrans)
catch e
@error "Failed to upload pre-made mask" exception=(e, catch_backtrace())
mask_editor_message = "Error uploading pre-made mask: $(sprint(showerror, e))"
mask_editor_warning = true
if !isempty(local_uploaded_path) && isfile(local_uploaded_path) rm(local_uploaded_path; force=true) end # Clean up
end
else
mask_editor_message = "No file selected for upload."
mask_editor_warning = true
end
end
@onbutton btn_change_slice begin @onbutton btn_change_slice begin
is_browsing_slices = true is_browsing_slices = true
is_editing_mask = false is_editing_mask = false
@ -783,17 +674,15 @@ end
# Update registry directly in the handler # Update registry directly in the handler
reg_path = abspath(joinpath(@__DIR__, "public", "registry.json")) reg_path = abspath(joinpath(@__DIR__, "public", "registry.json"))
registry = isfile(reg_path) ? JSON.parsefile(reg_path) : Dict{String, Any}() registry = isfile(reg_path) ? JSON.parsefile(reg_path) : Dict{String, Any}()
full_mask_path = abspath(final_path)
# Update the registry entry # Update the registry entry
if haskey(registry, selected_folder_main) if haskey(registry, selected_folder_main)
registry[selected_folder_main]["mask_path"] = full_mask_path registry[selected_folder_main]["mask_path"] = "/css/masks/$(final_mask_name)"
registry[selected_folder_main]["has_mask"] = true registry[selected_folder_main]["has_mask"] = true
else else
# Create a new entry if folder doesn't exist in registry # Create a new entry if folder doesn't exist in registry
registry[selected_folder_main] = Dict( registry[selected_folder_main] = Dict(
"mask_path" => full_mask_path, "mask_path" => "/css/masks/$(final_mask_name)",
"has_mask" => true, "has_mask" => true,
"is_imzML" => true, "is_imzML" => true,
"processed_date" => string(Dates.now()) "processed_date" => string(Dates.now())
@ -801,7 +690,9 @@ end
end end
# Save the updated registry # Save the updated registry
save_registry(reg_path, registry) open(reg_path, "w") do f
JSON.print(f, registry, 4)
end
@info "Registry updated with mask: $(final_mask_name)" @info "Registry updated with mask: $(final_mask_name)"
@ -821,74 +712,61 @@ end
end end
end end
@onbutton refetch_folders begin
# First, re-run the logic to populate the folder list
registry = load_registry(registry_path)
all_folders = sort(collect(keys(registry)), lt=natural)
img_folders = filter(folder -> get(get(registry, folder, Dict()), "is_imzML", false), all_folders)
available_folders = deepcopy(all_folders)
image_available_folders = deepcopy(img_folders)
# Apply the logic to select the first item if none is selected
if !isempty(image_available_folders) && isempty(selected_folder_main)
# This ensures the list is up-to-date
image_available_folders = deepcopy(img_folders)
selected_folder_main = first(image_available_folders)
end
end
@onchange isready begin @onchange isready begin
if isready && !registry_init_done if isready && !registry_init_done
sleep(1.0) # Give frontend time to initialize @async begin # Run asynchronously to not block startup
try sleep(1.0) # Give frontend time to initialize
println("Synchronizing registry for mask editor...") try
reg_path = abspath(joinpath(@__DIR__, "public", "registry.json")) println("Synchronizing registry for mask editor...")
# Assuming load_registry is available from MSI_src or a similar utility file reg_path = abspath(joinpath(@__DIR__, "public", "registry.json"))
# For now, handle its absence gracefully if it's not explicitly defined here. # Assuming load_registry is available from MSI_src or a similar utility file
registry = isfile(reg_path) ? load_registry(reg_path) : Dict{String, Any}() # For now, handle its absence gracefully if it's not explicitly defined here.
registry = isfile(reg_path) ? load_registry(reg_path) : Dict{String, Any}()
public_dirs = isdir("public") ? readdir("public") : [] public_dirs = isdir("public") ? readdir("public") : []
ignored_dirs = ["css", "masks"] ignored_dirs = ["css", "masks"]
dataset_dirs = filter(d -> isdir(joinpath("public", d)) && !(d in ignored_dirs), public_dirs) dataset_dirs = filter(d -> isdir(joinpath("public", d)) && !(d in ignored_dirs), public_dirs)
registry_keys = Set(keys(registry)) registry_keys = Set(keys(registry))
folder_set = Set(dataset_dirs) folder_set = Set(dataset_dirs)
new_folders = setdiff(folder_set, registry_keys) new_folders = setdiff(folder_set, registry_keys)
for folder in new_folders for folder in new_folders
println("Found new folder: $folder") println("Found new folder: $folder")
registry[folder] = Dict( registry[folder] = Dict(
"source_path" => "unknown (manually added)", "source_path" => "unknown (manually added)",
"processed_date" => "unknown", "processed_date" => "unknown",
"metadata" => Dict(), "metadata" => Dict(),
"is_imzML" => true # Assume folder contains images if found this way "is_imzML" => true # Assume folder contains images if found this way
) )
end
removed_folders = setdiff(registry_keys, folder_set)
for folder in removed_folders
delete!(registry, folder)
end
if !isempty(new_folders) || !isempty(removed_folders)
println("Registry changed, saving...")
open(reg_path, "w") do f
JSON.print(f, registry, 4)
end
end
all_folders = sort(collect(keys(registry)), lt=natural)
img_folders = filter(folder -> get(get(registry, folder, Dict()), "is_imzML", false), all_folders)
available_folders = deepcopy(all_folders)
image_available_folders = deepcopy(img_folders)
println("Mask editor UI lists updated. All: $(length(available_folders)), Images: $(length(image_available_folders))")
catch e
@warn "Mask editor registry synchronization failed: $e"
available_folders = []
image_available_folders = []
finally
registry_init_done = true
end end
removed_folders = setdiff(registry_keys, folder_set)
for folder in removed_folders
delete!(registry, folder)
end
if !isempty(new_folders) || !isempty(removed_folders)
println("Registry changed, saving...")
save_registry(reg_path, registry)
end
all_folders = sort(collect(keys(registry)), lt=natural)
img_folders = filter(folder -> get(get(registry, folder, Dict()), "is_imzML", false), all_folders)
available_folders = deepcopy(all_folders)
image_available_folders = deepcopy(img_folders)
println("Mask editor UI lists updated. All: $(length(available_folders)), Images: $(length(image_available_folders))")
catch e
@warn "Mask editor registry synchronization failed: $e"
available_folders = []
image_available_folders = []
finally
registry_init_done = true
end end
end end
end end

166
mask.jl.html
View File

@ -1,93 +1,72 @@
<header id="header"> <header id="header">
<img src="/css/LABI_logo.png" alt="Labi Logo Icon" id="imgLogo"> <img src="/css/LABI_logo.png" alt="Labi Logo Icon" id="imgLogo">
<div> <div>
<h4>JuliaMSI - Mask Editor&nbsp;</h4> <h4>JuliaMSI - Mask Editor&nbsp;</h4>
</div> </div>
</header> </header>
<div id="extDivStyle" class="row col-12 q-pa-xl"> <div id="extDivStyle" class="row col-12 q-pa-xl">
<div class="row col-4"> <div class="row col-4">
<!-- Left Panel: Controls --> <!-- Left Panel: Controls -->
<div class="st-col col-12 st-module q-pa-md"> <div class="st-col col-12 st-module q-pa-md">
<div class="text-h6">Mask Generation Workflow</div> <div class="text-h6">Mask Generation Workflow</div>
<!-- Step 1: Select Slice -->
<div class="text-subtitle1 q-mt-md">Step 1: Select Slice</div>
<div class="row items-center">
<q-select v-model="selected_folder_main" :options="image_available_folders" label="Select Dataset" class="q-ma-sm col" standout="custom-standout" :disable="is_editing_mask"></q-select>
</div>
<p class="text-center" v-html="msgimg"></p>
<!-- Step 1: Select Slice --> <!-- Step 2: Provide Mask Input -->
<div class="text-subtitle1 q-mt-md">Step 1: Select Slice</div> <div class="text-subtitle1 q-mt-md">Step 2: Create/Upload Mask</div>
<div class="row items-center"> <div class="row justify-around">
<q-select v-model="selected_folder_main" :options="image_available_folders" label="Select Dataset" <q-btn class="q-ma-sm btn-style" v-on:click="btn_use_slice_as_mask = true" label="Use Current Slice" :disable="!is_browsing_slices || !imgInt"></q-btn>
class="q-ma-sm col" standout="custom-standout" :disable="is_editing_mask" v-on:focus="refetch_folders = true"></q-select> <q-btn class="q-ma-sm btn-style" v-on:click="btn_upload_mask = true" label="Upload to Create Mask" :disable="!is_browsing_slices || !imgInt"></q-btn>
</div> <q-btn class="q-ma-sm btn-style" v-on:click="btn_change_slice = true" label="Change Mask Blueprint" :disable="!is_editing_mask"></q-btn>
<p class="text-center" v-html="msgimg"></p> </div>
<!-- Step 2: Provide Mask Input --> <!-- Step 3: Automated Processing & Overlay -->
<div class="text-subtitle1 q-mt-md">Step 2: Create/Upload Mask</div> <div class="text-subtitle1 q-mt-md">Step 3: Adjust & Verify</div>
<div class="row justify-around"> <div class="q-mt-sm">
<q-btn class="q-ma-sm btn-style" v-on:click="btn_use_slice_as_mask = true" label="Use Current Slice" <p class="text-caption text-center">Adjust automated processing for the initial mask:</p>
:disable="!is_browsing_slices || !imgInt"></q-btn> <q-slider color="primary" label-always v-model="otsu_scale" :label-value="'Otsu Scale: ' + otsu_scale.toFixed(2)" :step="0.01" :min="0.1" :max="2.0" :disable="!is_editing_mask"></q-slider>
<q-btn class="q-ma-sm btn-style" v-on:click="btn_upload_mask = true" label="Upload to Create Mask" <q-slider color="primary" label-always v-model="noise_size_percent" :label-value="'Noise Size: ' + (noise_size_percent * 100).toFixed(3) + '%'" :step="0.001" :min="0.001" :max="0.5" :disable="!is_editing_mask"></q-slider>
:disable="!is_browsing_slices || !imgInt"></q-btn> <q-slider color="primary" label-always v-model="hole_size_percent" :label-value="'Hole Size: ' + (hole_size_percent * 100).toFixed(4) + '%'" :step="0.0005" :min="0.0005" :max="0.5" :disable="!is_editing_mask"></q-slider>
<q-btn class="q-ma-sm btn-style" v-on:click="btn_upload_premade_mask = true" label="Upload Pre-made Mask" <q-slider color="primary" label-always v-model="smoothing_level" :label-value="'Smoothing: ' + smoothing_level + 'px'" :step="2" :min="1" :max="9" :disable="!is_editing_mask"></q-slider>
:disable="!is_browsing_slices"></q-btn> </div>
<q-btn class="q-ma-sm btn-style" v-on:click="btn_change_slice = true" label="Change Mask Blueprint" <div class="q-mt-md">
:disable="!is_editing_mask"></q-btn> <p class="text-caption text-center">Adjust slice overlay transparency:</p>
</div> <q-slider color="black" v-model="imgTrans" :min="0.0" :max="1.0" :step="0.05" label-always :label-value="'Transparency: ' + imgTrans.toFixed(2)" :disable="!is_editing_mask" />
</div>
<!-- Step 3: Automated Processing & Overlay --> <p class="q-mt-md" :class="{'text-negative': mask_editor_warning}">{{ mask_editor_message }}</p>
<div class="text-subtitle1 q-mt-md">Step 3: Adjust & Verify</div> <div class="row">
<div class="q-mt-sm"> <q-btn :loading="progress" class="q-ma-sm btn-style" :disable="!is_editing_mask" v-on:click="btn_save_final_mask = !btn_save_final_mask" padding="lg" icon="save" label="Save Final Mask"/>
<p class="text-caption text-center">Adjust automated processing for the initial mask:</p> <q-btn class="q-ma-sm btn-style" icon="arrow_back" label="Return" href="/" ></q-btn>
<q-slider color="primary" label-always v-model="otsu_scale"
:label-value="'Otsu Scale: ' + otsu_scale.toFixed(2)" :step="0.01" :min="0.1" :max="2.0"
:disable="!is_editing_mask"></q-slider>
<q-slider color="primary" label-always v-model="noise_size_percent"
:label-value="'Noise Size: ' + (noise_size_percent * 100).toFixed(3) + '%'" :step="0.001"
:min="0.001" :max="0.5" :disable="!is_editing_mask"></q-slider>
<q-slider color="primary" label-always v-model="hole_size_percent"
:label-value="'Hole Size: ' + (hole_size_percent * 100).toFixed(4) + '%'" :step="0.0005"
:min="0.0005" :max="0.5" :disable="!is_editing_mask"></q-slider>
<q-slider color="primary" label-always v-model="smoothing_level"
:label-value="'Smoothing: ' + smoothing_level + 'px'" :step="2" :min="1" :max="9"
:disable="!is_editing_mask"></q-slider>
</div>
<div class="q-mt-md">
<p class="text-caption text-center">Adjust slice overlay transparency:</p>
<q-slider color="black" v-model="imgTrans" :min="0.0" :max="1.0" :step="0.05" label-always
:label-value="'Transparency: ' + imgTrans.toFixed(2)" :disable="!is_editing_mask" />
</div>
<p class="q-mt-md" :class="{'text-negative': mask_editor_warning}">{{ mask_editor_message }}</p>
<div class="row">
<q-btn :loading="progress" class="q-ma-sm btn-style" :disable="!is_editing_mask"
v-on:click="btn_save_final_mask = !btn_save_final_mask" padding="lg" icon="save"
label="Save Final Mask" />
<q-btn class="q-ma-sm btn-style" icon="arrow_back" label="Return" href="/"></q-btn>
</div>
</div> </div>
</div> </div>
</div>
<div class="row col-6"> <div class="row col-6">
<!-- Right Panel: Image Displays --> <!-- Right Panel: Image Displays -->
<div id="intDivStyle-right" class="st-col col-12 st-module"> <div id="intDivStyle-right" class="st-col col-12 st-module">
<div class="text-h6">Mask Preview</div> <div class="text-h6">Mask Preview</div>
<div class="row items-center"> <div class="row items-center">
<q-btn icon="arrow_back" class="q-my-sm btn-style" v-on:click="btn_img_minus=true" <q-btn icon="arrow_back" class="q-my-sm btn-style" v-on:click="btn_img_minus=true" :disable="is_editing_mask"></q-btn>
:disable="is_editing_mask"></q-btn> <q-btn icon="arrow_forward" class="q-my-sm on-right btn-style" v-on:click="btn_img_plus=true" :disable="is_editing_mask"></q-btn>
<q-btn icon="arrow_forward" class="q-my-sm on-right btn-style" v-on:click="btn_img_plus=true"
:disable="is_editing_mask"></q-btn>
</div> </div>
<div v-if="!show_verification_plot" class="text-center text-grey q-pa-xl"> <div v-if="!show_verification_plot" class="text-center text-grey q-pa-xl">
<q-icon name="image" size="5em" class="q-mb-md" /> <q-icon name="image" size="5em" class="q-mb-md" />
<p>Please select a dataset and provide a mask input to begin.</p> <p>Please select a dataset and provide a mask input to begin.</p>
<p class="text-caption">Use the controls on the left to load a slice or upload a mask image.</p> <p class="text-caption">Use the controls on the left to load a slice or upload a mask image.</p>
</div> </div>
<div v-if="show_verification_plot"> <div v-if="show_verification_plot">
<plotly :data="plotdata_verify" :layout="plotlayout_verify" class="q-pa-none q-ma-none pixelated-plot"> <plotly :data="plotdata_verify" :layout="plotlayout_verify" class="q-pa-none q-ma-none pixelated-plot"></plotly>
</plotly>
<div class="q-mt-md text-center"> <div class="q-mt-md text-center">
<q-btn color="primary" class="btn-style" label="Edit Manually" <q-btn color="primary" class="btn-style" label="Edit Manually" v-on:click="btn_edit_manually=true" />
v-on:click="btn_edit_manually=true" />
</div> </div>
</div> </div>
<p class="text-center" v-html="msgimg"></p> <p class="text-center" v-html="msgimg"></p>
@ -108,59 +87,54 @@
<div class="q-pa-md q-gutter-sm row justify-center items-center"> <div class="q-pa-md q-gutter-sm row justify-center items-center">
<!-- Tool Selection --> <!-- Tool Selection -->
<q-btn-group> <q-btn-group>
<q-btn label="Brush" v-on:click="current_tool = 'brush'" <q-btn label="Brush" v-on:click="current_tool = 'brush'" :color="current_tool === 'brush' ? 'primary' : 'white'" text-color="black"/>
:color="current_tool === 'brush' ? 'primary' : 'white'" text-color="black" /> <q-btn label="Eraser" v-on:click="current_tool = 'eraser'" :color="current_tool === 'eraser' ? 'primary' : 'white'" text-color="black"/>
<q-btn label="Eraser" v-on:click="current_tool = 'eraser'" <q-btn label="Bucket" v-on:click="current_tool = 'bucket'" :color="current_tool === 'bucket' ? 'primary' : 'white'" text-color="black"/>
:color="current_tool === 'eraser' ? 'primary' : 'white'" text-color="black" /> <q-btn label="Drag" v-on:click="current_tool = 'drag'" :color="current_tool === 'drag' ? 'primary' : 'white'" text-color="black"/>
<q-btn label="Bucket" v-on:click="current_tool = 'bucket'"
:color="current_tool === 'bucket' ? 'primary' : 'white'" text-color="black" />
<q-btn label="Drag" v-on:click="current_tool = 'drag'"
:color="current_tool === 'drag' ? 'primary' : 'white'" text-color="black" />
</q-btn-group> </q-btn-group>
<q-separator vertical inset /> <q-separator vertical inset />
<!-- Tool Settings --> <!-- Tool Settings -->
<q-input dense filled label="Brush Size" type="number" v-model.number="brush_size" <q-input dense filled label="Brush Size" type="number" v-model.number="brush_size" style="max-width: 120px" :min="1"/>
style="max-width: 120px" :min="1" /> <q-input dense filled label="Zoom" type="number" step="0.1" v-model.number="editor_scale" style="max-width: 120px" :min="0.1"/>
<q-input dense filled label="Zoom" type="number" step="0.1" v-model.number="editor_scale"
style="max-width: 120px" :min="0.1" />
<q-separator vertical inset /> <q-separator vertical inset />
<!-- Image Manipulation --> <!-- Image Manipulation -->
<q-input dense filled label="Rotate (deg)" type="number" v-model.number="rotate_degrees" <q-input dense filled label="Rotate (deg)" type="number" v-model.number="rotate_degrees" style="max-width: 150px" :step="90"/>
style="max-width: 150px" :step="90" /> <q-select dense filled label="Flip" v-model="flip_direction" :options="['horizontal', 'vertical']" style="max-width: 150px"/>
<q-select dense filled label="Flip" v-model="flip_direction" :options="['horizontal', 'vertical']" <q-btn label="Apply Flip" v-on:click="btn_flip_mask = !btn_flip_mask"/>
style="max-width: 150px" />
<q-btn label="Apply Flip" v-on:click="btn_flip_mask = !btn_flip_mask" />
</div> </div>
<!-- Drag Tool Controls --> <!-- Drag Tool Controls -->
<div v-if="current_tool === 'drag'" class="q-pa-sm q-gutter-sm row justify-center items-center" <div v-if="current_tool === 'drag'" class="q-pa-sm q-gutter-sm row justify-center items-center" style="border: 1px solid #ccc; border-radius: 4px;">
style="border: 1px solid #ccc; border-radius: 4px;"> <q-input dense filled label="Move Pixels" type="number" v-model.number="move_pixels" style="max-width: 120px" :min="1"/>
<q-input dense filled label="Move Pixels" type="number" v-model.number="move_pixels"
style="max-width: 120px" :min="1" />
<q-btn icon="arrow_upward" v-on:click="move_mask_payload = { direction: 'up', t: Date.now() }" /> <q-btn icon="arrow_upward" v-on:click="move_mask_payload = { direction: 'up', t: Date.now() }" />
<q-btn icon="arrow_downward" v-on:click="move_mask_payload = { direction: 'down', t: Date.now() }" /> <q-btn icon="arrow_downward" v-on:click="move_mask_payload = { direction: 'down', t: Date.now() }" />
<q-btn icon="arrow_back" v-on:click="move_mask_payload = { direction: 'left', t: Date.now() }" /> <q-btn icon="arrow_back" v-on:click="move_mask_payload = { direction: 'left', t: Date.now() }" />
<q-btn icon="arrow_forward" v-on:click="move_mask_payload = { direction: 'right', t: Date.now() }" /> <q-btn icon="arrow_forward" v-on:click="move_mask_payload = { direction: 'right', t: Date.now() }" />
</div> </div>
<!-- Interactive Canvas --> <!-- Interactive Canvas -->
<div class="row justify-center q-pa-md"> <div class="row justify-center q-pa-md">
<div style="position: relative; display: inline-block; overflow: auto; max-width: 100%;"> <div style="position: relative; display: inline-block; overflow: auto; max-width: 100%;"> <!-- Added for scrolling very large images -->
<!-- Added for scrolling very large images --> <canvas id="maskCanvas"
<canvas id="maskCanvas" :width="canvas_width" :height="canvas_height" :style="{ :width="canvas_width"
:height="canvas_height"
:style="{
width: canvas_width * editor_scale + 'px', width: canvas_width * editor_scale + 'px',
height: canvas_height * editor_scale + 'px', height: canvas_height * editor_scale + 'px',
border: '2px solid #ccc', border: '2px solid #ccc',
cursor: current_tool === 'drag' ? 'move' : 'crosshair', cursor: current_tool === 'drag' ? 'move' : 'crosshair',
imageRendering: 'pixelated' imageRendering: 'pixelated'
}" v-on:mousedown.prevent="startDrawing" v-on:mousemove.prevent="draw" }"
v-on:mouseup.prevent="stopDrawing()" v-on:mouseleave.prevent="stopDrawing()"> v-on:mousedown.prevent="startDrawing"
</canvas> v-on:mousemove.prevent="draw"
</div> v-on:mouseup.prevent="stopDrawing()"
v-on:mouseleave.prevent="stopDrawing()">
</canvas>
</div>
</div> </div>
</q-card-section> </q-card-section>
@ -168,4 +142,4 @@
<q-btn label="Cancel" color="negative" v-on:click="show_editor = false"/> <q-btn label="Cancel" color="negative" v-on:click="show_editor = false"/>
</q-card-actions>--> </q-card-actions>-->
</q-card> </q-card>
</q-dialog> </q-dialog>

173
precompile_script.jl
View File

@ -1,173 +0,0 @@
# precompile_script.jl
# This script exercises the core MSI processing kernels (imzML and mzML)
# to ensure they are precompiled into the custom system image.
using MSI_src
using Mmap
using Base64
using PlotlyBase
# --- DYNAMIC WARM-UP WORKLOAD --- #
function create_warmup_datasets(dir, T_mz::Type, T_int::Type)
# 1. Create minimal imzML/ibd
imzml_path = joinpath(dir, "warmup_$(T_mz)_$(T_int).imzML")
ibd_path = joinpath(dir, "warmup_$(T_mz)_$(T_int).ibd")
n_points = 10
n_pixels = 1
mz_bytes = n_points * sizeof(T_mz)
int_bytes = n_points * sizeof(T_int)
# Write some dummy binary data (m/z must be sorted for validation)
open(ibd_path, "w") do f
write(f, sort(rand(T_mz, n_points)))
write(f, rand(T_int, n_points))
end
# Minimal but VALID imzML structure that passes axes_config_img
imzml_content = """<?xml version="1.0" encoding="utf-8"?>
<mzML xmlns="http://psi.hupo.org/ms/mzML" version="1.1.0">
<cvList count="2">
<cv id="MS" fullName="Proteomics Standards Initiative Mass Spectrometry Ontology" version="3.30.0" URI="http://psidev.cvs.sourceforge.net/*checkout*/psidev/ms/mzML/ontology/psi-ms.obo"/>
<cv id="IMS" fullName="Imaging MS Ontology" version="0.9.1" URI="http://www.imzml.org/ontology/imzML1.1.0.obo"/>
</cvList>
<fileDescription>
<fileContent>
<cvParam cvRef="IMS" accession="IMS:1000031" name="processed" value=""/>
<cvParam cvRef="IMS" accession="IMS:1000080" name="universally unique identifier" value="00000000-0000-0000-0000-000000000000"/>
<cvParam cvRef="IMS" accession="IMS:1000091" name="ibd MD5 HTTP" value="00000000000000000000000000000000"/>
</fileContent>
</fileDescription>
<referenceableParamGroupList count="2">
<referenceableParamGroup id="mz_array_settings">
<cvParam cvRef="MS" accession="MS:1000514" name="m/z array" value=""/>
<cvParam cvRef="MS" accession="MS:$(T_mz == Float64 ? "1000523" : "1000521")" name="$(T_mz == Float64 ? "64-bit float" : "32-bit float")" value=""/>
<cvParam cvRef="MS" accession="MS:1000574" name="zlib compression" value=""/>
</referenceableParamGroup>
<referenceableParamGroup id="intensity_array_settings">
<cvParam cvRef="MS" accession="MS:1000515" name="intensity array" value=""/>
<cvParam cvRef="MS" accession="MS:$(T_int == Float64 ? "1000523" : "1000521")" name="$(T_int == Float64 ? "64-bit float" : "32-bit float")" value=""/>
<cvParam cvRef="MS" accession="MS:1000574" name="zlib compression" value=""/>
</referenceableParamGroup>
</referenceableParamGroupList>
<run id="run_1">
<spectrumList count="1">
<spectrum index="0" id="scan=1" defaultArrayLength="$n_points">
<cvParam cvRef="IMS" accession="IMS:1000050" name="coordinate x" value="1"/>
<cvParam cvRef="IMS" accession="IMS:1000051" name="coordinate y" value="1"/>
<binaryDataArrayList count="2">
<binaryDataArray encodedLength="0">
<referenceableParamGroupRef ref="mz_array_settings"/>
<cvParam cvRef="IMS" accession="IMS:1000102" name="external data" value=""/>
<cvParam cvRef="IMS" accession="IMS:1000103" name="external offset" value="0"/>
<cvParam cvRef="IMS" accession="IMS:1000104" name="external array length" value="$n_points"/>
<cvParam cvRef="IMS" accession="IMS:1000101" name="external encoded length" value="$mz_bytes"/>
<binary/>
</binaryDataArray>
<binaryDataArray encodedLength="0">
<referenceableParamGroupRef ref="intensity_array_settings"/>
<cvParam cvRef="IMS" accession="IMS:1000102" name="external data" value=""/>
<cvParam cvRef="IMS" accession="IMS:1000103" name="external offset" value="$mz_bytes"/>
<cvParam cvRef="IMS" accession="IMS:1000104" name="external array length" value="$n_points"/>
<cvParam cvRef="IMS" accession="IMS:1000101" name="external encoded length" value="$int_bytes"/>
<binary/>
</binaryDataArray>
</binaryDataArrayList>
</spectrum>
</spectrumList>
</run>
</mzML>
"""
write(imzml_path, imzml_content)
# 2. Create minimal mzML (Base64)
mzml_path = joinpath(dir, "warmup_$(T_mz)_$(T_int).mzML")
b64_mz = base64encode(sort(rand(T_mz, n_points)))
b64_int = base64encode(rand(T_int, n_points))
mzml_content = """<?xml version="1.0" encoding="utf-8"?>
<mzML xmlns="http://psi.hupo.org/ms/mzML" version="1.1.0">
<run id="run_1">
<spectrumList count="1">
<spectrum index="0" id="scan=1" defaultArrayLength="$n_points">
<binaryDataArrayList count="2">
<binaryDataArray encodedLength="$(length(b64_mz))">
<cvParam cvRef="MS" accession="MS:1000514" name="m/z array" value=""/>
<cvParam cvRef="MS" accession="MS:$(T_mz == Float64 ? "1000523" : "1000521")" name="" value=""/>
<binary>$b64_mz</binary>
</binaryDataArray>
<binaryDataArray encodedLength="$(length(b64_int))">
<cvParam cvRef="MS" accession="MS:1000515" name="intensity array" value=""/>
<cvParam cvRef="MS" accession="MS:$(T_int == Float64 ? "1000523" : "1000521")" name="" value=""/>
<binary>$b64_int</binary>
</binaryDataArray>
</binaryDataArrayList>
</spectrum>
</spectrumList>
</run>
<indexList count="1"><index name="spectrum"><offset idRef="scan=1">PLACEHOLDER</offset></index></indexList>
<indexListOffset>ID_OFFSET</indexListOffset>
</mzML>
"""
spec_start = findfirst("<spectrum", mzml_content).start - 1
index_start = findfirst("<indexList", mzml_content).start - 1
mzml_content = replace(mzml_content, "PLACEHOLDER" => string(spec_start))
mzml_content = replace(mzml_content, "ID_OFFSET" => string(index_start))
write(mzml_path, mzml_content)
return imzml_path, mzml_path
end
println("Starting robust precompilation warmup (imzML + mzML)...")
try
mktempdir() do tmp_dir
for (T_mz, T_int) in [(Float32, Float32), (Float64, Float32)]
println("Exercising pathways for $T_mz and $T_int...")
imzml_path, mzml_path = create_warmup_datasets(tmp_dir, T_mz, T_int)
try
# Exercise imzML pathway
imzml_data = MSI_src.load_imzml_lazy(imzml_path; use_mmap=true)
MSI_src.get_multiple_mz_slices(imzml_data, [10.0, 20.0], 0.1)
# Exercise Sprint 2 Streaming Pipeline
config = MSI_src.PipelineConfig(
steps=[
MSI_src.StreamingStep(:baseline_correction, Dict(:method => :snip, :iterations => 5)),
MSI_src.StreamingStep(:normalization, Dict(:method => :tic)),
MSI_src.StreamingStep(:peak_picking, Dict(:method => :profile, :snr_threshold => 3.0))
],
num_bins=2000
)
try
MSI_src.process_dataset!(imzml_data, config)
catch
# Ignore matrix mismatch errors from tiny random data
end
# Exercise mzML pathway
mzml_data = MSI_src.load_mzml_lazy(mzml_path)
MSI_src.GetSpectrum(mzml_data, 1)
catch e
@warn "Warmup failed for $T_mz/$T_int: $e"
end
end
# --- COMMON KERNELS ---
println("Precompiling image processing kernels...")
dummy_slice = rand(Float32, 4, 4)
MSI_src.TrIQ(dummy_slice, 256, 1.0)
MSI_src.save_bitmap(joinpath(tmp_dir, "d.bmp"), zeros(UInt8, 4, 4), MSI_src.ViridisPalette)
# Plotly Warmup
p = PlotlyBase.Plot(PlotlyBase.scatter(x=1:2, y=[1,2]))
end
catch e
if !(e isa InterruptException)
@warn "Global Warmup failed: $e"
end
end
println("Precompilation workload completed.")

84
public/css/genieapp.css
View File

@ -41,46 +41,13 @@
color: rgb(0, 0, 0) !important; color: rgb(0, 0, 0) !important;
} }
/* Placeholder color for custom-standout q-inputs */
.custom-standout .q-field__native::placeholder {
color: #CFD8DC !important; /* Lighter grey for placeholder text */
opacity: 1 !important; /* Ensure full visibility */
}
.custom-standout .q-field__control::before {
border-color: #7f8389 !important; /* Keep the original border color */
}
/* Placeholder color for various browsers */
.custom-standout input::placeholder {
color: #CFD8DC !important;
opacity: 1 !important;
}
.custom-standout input::-webkit-input-placeholder { /* WebKit, Blink, Edge */
color: #CFD8DC !important;
opacity: 1 !important;
}
.custom-standout input::-moz-placeholder { /* Mozilla Firefox 19+ */
color: #CFD8DC !important;
opacity: 1 !important;
}
.custom-standout input:-ms-input-placeholder { /* Internet Explorer 10-11 */
color: #CFD8DC !important;
opacity: 1 !important;
}
.custom-standout input::-ms-input-placeholder { /* Microsoft Edge */
color: #CFD8DC !important;
opacity: 1 !important;
}
#tabHeader { #tabHeader {
color: #009f90; color: #009f90;
} }
#intDivStyle-left, #intDivStyle-right { #intDivStyle-left {
border-radius: 10px; border-radius: 10px;
padding: 10px; padding: 10px;
background-color: #fdfdfd; /* Add a background color to distinguish from the outer div */
} }
/* Tab styling to match your theme */ /* Tab styling to match your theme */
@ -114,9 +81,13 @@
font-size: 0.9rem; font-size: 0.9rem;
} }
#intDivStyle-left .q-tab-panel, #intDivStyle-right .q-tab-panel { #intDivStyle-left .q-tab-panels {
height: auto; /* Let content define height */ height: 700px; /* Set this to accommodate your tallest content */
overflow-y: hidden; /* Remove scrollbar */ }
#intDivStyle-left .q-tab-panel {
height: 100%;
overflow-y: auto; /* Add scroll if content overflows */
} }
.pixelated-plot svg image { .pixelated-plot svg image {
@ -127,43 +98,4 @@
*{ *{
font-family: 'Roboto', 'Lato', sans-serif; font-family: 'Roboto', 'Lato', sans-serif;
}
.q-card {
border-radius: 8px !important;
box-shadow: 0 2px 4px rgba(0,0,0,0.1) !important;
margin-bottom: 16px !important;
}
.q-card__section--dark {
background: #f7f7f7 !important;
}
.q-radio__label {
font-size: 1rem !important;
}
.q-option-group > div {
margin-bottom: 8px;
}
/* Loading Overlay Styles */
.loading-overlay {
position: fixed;
top: 0;
left: 0;
width: 100%;
height: 100%;
background-color: rgba(0, 0, 0, 0.7);
display: flex;
justify-content: center;
align-items: center;
z-index: 9999;
text-align: center;
}
.loading-content .q-spinner {
color: white;
}
.loading-content .text-h6 {
color: white;
} }

279
src/BloomFilters.jl
View File

@ -1,279 +0,0 @@
# src/BloomFilter.jl
"""
BloomFilter{T}
A simple, efficient Bloom filter implementation for probabilistic set membership testing.
# Fields
- `bits::BitVector`: The underlying bit array
- `size::Int`: Number of bits in the filter
- `hash_count::Int`: Number of hash functions to use
- `seed::UInt64`: Random seed for hash functions
- `count::Int`: Number of elements added (for monitoring)
"""
mutable struct BloomFilter{T}
bits::BitVector
size::Int
hash_count::Int
seed::UInt64
count::Int
end
"""
BloomFilter{T}(expected_elements::Int, false_positive_rate::Float64=0.01; kwargs...) -> Return type
Creates a Bloom filter optimized for the expected number of elements and desired false positive rate.
# Arguments
- `expected_elements::Int`: Argument description
- `false_positive_rate::Float64`: Argument description
(**Default**: `0.01`)
# Keywords
- `seed::Union{UInt32,UInt64}`: Keyword description
(**Default**: `0x12345678`)
# Returns
- A `BloomFilter{T}`.
"""
function BloomFilter{T}(expected_elements::Int, false_positive_rate::Float64=0.01; seed::Union{UInt32,UInt64}=0x12345678) where T
# Convert seed to UInt64 for consistency
seed_uint64 = UInt64(seed)
# Calculate optimal parameters using standard formulas
size = optimal_bit_size(expected_elements, false_positive_rate)
hash_count = optimal_hash_count(expected_elements, size)
bits = falses(size)
return BloomFilter{T}(bits, size, hash_count, seed_uint64, 0)
end
"""
optimal_bit_size(n::Int, p::Float64) -> Int
Calculates the optimal number of bits for a Bloom filter
# Arguments
- `n::Int`: Expected number of elements
- `p::Float64`: Desired false positive rate
# Returns
- An `Int` to determine the optimal number of bits for a Bloom filter.
"""
function optimal_bit_size(n::Int, p::Float64)::Int
if p <= 0.0 || p >= 1.0
throw(ArgumentError("False positive rate must be between 0 and 1"))
end
# m = - (n * ln(p)) / (ln(2)^2)
m = ceil(Int, - (n * log(p)) / (log(2)^2))
return max(m, 1)
end
"""
optimal_hash_count(n::Int, m::Int) -> Int
Calculates the optimal number of hash functions for a Bloom filter.
# Arguments
- `n::Int`: Expected number of elements
- `m::Int`: Number of bits in the filter
# Returns
- An `Int` to determine the optimal number of hash functions for a Bloom filter.
"""
function optimal_hash_count(n::Int, m::Int)::Int
if n <= 0 || m <= 0
return 1
end
# k = (m / n) * ln(2)
k = max(1, round(Int, (m / n) * log(2)))
return min(k, 8) # Practical limit to avoid too many hashes
end
"""
hash_functions(item::T, count::Int, size::Int, seed::UInt64) -> Vector{Int}
Generates multiple hash values for an item using double hashing technique.
# Arguments
- `item::T`: Argument description
- `count::Int`: Argument description
- `size::Int`: Argument description
- `seed::UInt64`: Argument description
# Returns
- A `Vector{Int}` of hash values.
"""
function hash_functions(item::T, count::Int, size::Int, seed::UInt64) where T
# Use Julia's built-in hash with different seeds
hashes = Vector{Int}(undef, count)
# First hash with the main seed
h1 = hash(item, seed)
h2 = hash(item, seed + 1)
for i in 1:count
# Double hashing: h_i = h1 + i * h2
combined_hash = UInt64(h1) + UInt64(i) * UInt64(h2)
hashes[i] = (combined_hash % UInt64(size)) + 1 # 1-based indexing
end
return hashes
end
"""
Base.push!(bf::BloomFilter{T}, item::T)
Adds an element to the Bloom filter.
# Arguments:
- `bf::BloomFilter{T}`: Argument description
- `item::T`: Argument description
"""
function Base.push!(bf::BloomFilter{T}, item::T) where T
hashes = hash_functions(item, bf.hash_count, bf.size, bf.seed)
for h in hashes
bf.bits[h] = true
end
bf.count += 1
return bf
end
"""
Base.in(item::T, bf::BloomFilter{T}) -> Bool
Checks whether an element is possibly in the Bloom filter.
Returns `true` if the element might be in the set, `false` if it's definitely not.
# Arguments:
- `item::T`: Argument description
- `bf::BloomFilter{T}`: Argument description
"""
function Base.in(item::T, bf::BloomFilter{T})::Bool where T
hashes = hash_functions(item, bf.hash_count, bf.size, bf.seed)
for h in hashes
if !bf.bits[h]
return false # Definitely not in the set
end
end
return true # Possibly in the set
end
"""
contains(bf::BloomFilter{T}, item::T) -> Bool
Alias for `in()` for compatibility with your existing code.
# Arguments:
- `bf::BloomFilter{T}`: Argument description
- `item::T`: Argument description
"""
contains(bf::BloomFilter{T}, item::T) where T = item in bf
"""
add!(bf::BloomFilter{T}, item::T)
Alias for `push!()` for compatibility with your existing code.
# Arguments:
- `bf::BloomFilter{T}`: Argument description
- `item::T`: Argument description
"""
add!(bf::BloomFilter{T}, item::T) where T = push!(bf, item)
"""
false_positive_rate(bf::BloomFilter) -> Float64
Estimates the current false positive rate of the Bloom filter.
# Arguments:
- `bf::BloomFilter`: Argument description
"""
function false_positive_rate(bf::BloomFilter)::Float64
if bf.count == 0
return 0.0
end
# Theoretical false positive rate: (1 - e^(-k * n / m)) ^ k
k = bf.hash_count
n = bf.count
m = bf.size
return (1 - exp(-k * n / m)) ^ k
end
"""
fill_ratio(bf::BloomFilter)::Float64 return count(bf.bits) / length(bf.bits) end -> Return type
Returns the fraction of bits that are set to 1.
# Arguments
- `bf::BloomFilter`: Argument description
"""
function fill_ratio(bf::BloomFilter)::Float64
return count(bf.bits) / length(bf.bits)
end
"""
is_empty(bf::BloomFilter)::Bool return bf.count == 0 end -> Return type
Checks if the Bloom filter is empty (no elements added).
# Arguments
- `bf::BloomFilter`: Argument description
"""
function is_empty(bf::BloomFilter)::Bool
return bf.count == 0
end
"""
reset!(bf::BloomFilter) fill!(bf.bits, false) bf.count = 0 return bf end -> Return type
Clears the Bloom filter, removing all elements.
# Arguments
- `bf::BloomFilter`: Argument description
"""
function reset!(bf::BloomFilter)
fill!(bf.bits, false)
bf.count = 0
return bf
end
# Specialized constructor for empty Bloom filters
"""
BloomFilter{T}(; kwargs...) -> Return type
Description of the function
# Keywords
- `size::Int`: Keyword description
(**Default**: `100`)
- `hash_count::Int`: Keyword description
(**Default**: `3`)
- `seed::Union{UInt32,UInt64}`: Keyword description
(**Default**: `0x12345678`)
"""
function BloomFilter{T}(;size::Int=100, hash_count::Int=3, seed::Union{UInt32,UInt64}=0x12345678) where T
seed_uint64 = UInt64(seed)
bits = falses(size)
return BloomFilter{T}(bits, size, hash_count, seed_uint64, 0)
end

245
src/Common.jl
View File

@ -1,245 +0,0 @@
using Base.Threads
using Mmap
# POSIX madvise constants
const MADV_NORMAL = 0
const MADV_RANDOM = 1
const MADV_SEQUENTIAL = 2
const MADV_WILLNEED = 3
const MADV_DONTNEED = 4
"""
posix_madvise(buffer::AbstractArray, advice::Integer)
A safe wrapper for the OS `madvise` system call. Signals the kernel about the
access pattern for a memory-mapped region. Currently supports Linux/Unix systems.
"""
function posix_madvise(buffer::AbstractArray, advice::Integer)
@static if Sys.isunix()
try
ptr = pointer(buffer)
len = sizeof(buffer)
# ccall(:madvise, return_type, (arg_types...), args...)
ret = ccall(:madvise, Int32, (Ptr{Cvoid}, Csize_t, Int32), ptr, len, Int32(advice))
return ret == 0
catch
return false
end
else
return false # Not supported on this OS
end
end
# --- Buffer Pooling ---
"""
SimpleBufferPool
A dramatically simplified buffer pool that avoids complex locking.
# Arguments:
- `buffers::Dict{Int, Vector{Vector{UInt8}}}`: A dictionary mapping buffer size (in bytes) to a list of available buffers of that size.
- `max_pool_size::Int`: The maximum number of buffers to keep in the pool.
- `lock::ReentrantLock`: A lock to ensure thread-safe access to the pool.
"""
mutable struct SimpleBufferPool
buffers::Dict{Int, Vector{Vector{UInt8}}}
max_pool_size::Int
lock::ReentrantLock # Add lock for thread safety
end
"""
SimpleBufferPool()
Creates a new `SimpleBufferPool` with a default maximum pool size of 50.
"""
SimpleBufferPool() = SimpleBufferPool(Dict{Int, Vector{Vector{UInt8}}}(), 50, ReentrantLock())
"""
get_buffer!(pool::SimpleBufferPool, size::Int) -> Vector{UInt8}
Retrieves a `Vector{UInt8}` buffer of at least `size` bytes from the pool.
If no suitable buffer is available, a new one is allocated.
# Arguments:
- `pool`: The `SimpleBufferPool` to retrieve the buffer from.
- `size`: The minimum desired size of the buffer in bytes.
# Returns:
- A `Vector{UInt8}` buffer.
"""
function get_buffer!(pool::SimpleBufferPool, size::Int)::Vector{UInt8}
buf = lock(pool.lock) do
# Check for existing buffers of exact size first
if haskey(pool.buffers, size) && !isempty(pool.buffers[size])
return pop!(pool.buffers[size])
end
return nothing
end
if buf !== nothing
return buf
end
# No suitable buffer found, allocate new one (outside lock to reduce contention)
return Vector{UInt8}(undef, size)
end
"""
release_buffer!(pool::SimpleBufferPool, buffer::Vector{UInt8})
Returns a `Vector{UInt8}` buffer to the pool for future reuse.
# Arguments:
- `pool`: The `SimpleBufferPool` to return the buffer to.
- `buffer`: The `Vector{UInt8}` buffer to release.
"""
function release_buffer!(pool::SimpleBufferPool, buffer::Vector{UInt8})
size = length(buffer)
lock(pool.lock) do
if !haskey(pool.buffers, size)
pool.buffers[size] = Vector{Vector{UInt8}}()
end
# Limit pool size to prevent memory bloat
if length(pool.buffers[size]) < pool.max_pool_size
push!(pool.buffers[size], buffer)
end
end
# If pool is full, let buffer get GC'd
end
# --- Unified Error Types ---
abstract type MSIError <: Exception end
struct FileFormatError <: MSIError
msg::String
end
struct SpectrumNotFoundError <: MSIError
id
end
struct InvalidSpectrumError <: MSIError
id
reason::String
end
# --- Shared Constants ---
const DEFAULT_CACHE_SIZE = 100 # Default cache size for spectra
const DEFAULT_NUM_BINS = 2000 # Default number of bins for spectra
const DEFAULT_BLOOM_FILTER_SIZE = 10000 # Default size for empty spectra
const DEFAULT_FALSE_POSITIVE_RATE = 0.01 # Default false positive rate for bloom filters
"""
validate_spectrum_data(mz, intensity, id)
Checks a spectrum for common data integrity issues.
# Arguments:
- `mz`: The m/z values of the spectrum.
- `intensity`: The intensity values of the spectrum.
- `id`: The ID of the spectrum.
# Returns:
- `true` if the spectrum is valid.
# Throws:
- `InvalidSpectrumError` if the spectrum is invalid.
"""
function validate_spectrum_data(mz::AbstractVector, intensity::AbstractVector, id)
if length(mz) != length(intensity)
throw(InvalidSpectrumError(id, "m/z and intensity arrays have different lengths"))
end
if !all(isfinite, mz) || !all(isfinite, intensity)
throw(InvalidSpectrumError(id, "Spectrum contains NaN or Inf values"))
end
if !issorted(mz)
throw(InvalidSpectrumError(id, "m/z array is not sorted"))
end
return true
end
"""
create_bloom_filter_for_spectrum(mz::AbstractVector{<:Real}) -> BloomFilter{Int}
Memory-optimized version that processes data in chunks to reduce temporary allocations
and avoid holding large intermediate arrays.
# Arguments:
- `mz`: The m/z values of the spectrum.
# Returns:
- A `BloomFilter{Int}` for the spectrum.
"""
function create_bloom_filter_for_spectrum(mz::AbstractVector{<:Real})::BloomFilter{Int}
if isempty(mz)
return empty_bloom_filter()
end
# Size Bloom filter appropriately
n_points = length(mz)
bf = BloomFilter{Int}(n_points, DEFAULT_FALSE_POSITIVE_RATE)
discretization_factor = 100.0
chunk_size = 10000 # Process in chunks to reduce memory pressure
# Process spectrum data in chunks to avoid large temporary arrays
for chunk_start in 1:chunk_size:n_points
chunk_end = min(chunk_start + chunk_size - 1, n_points)
@inbounds for i in chunk_start:chunk_end
discretized_val = round(Int, mz[i] * discretization_factor)
push!(bf, discretized_val)
end
end
return bf
end
"""
empty_bloom_filter() -> BloomFilter{Float64}
Creates an empty Bloom filter for spectra with no data.
# Returns:
- An empty `BloomFilter{Float64}`.
"""
function empty_bloom_filter()::BloomFilter{Int}
return BloomFilter{Int}(size=DEFAULT_BLOOM_FILTER_SIZE, hash_count=3)
end
"""
AtomicFlag
A thread-safe boolean flag using atomic operations.
# Fields:
- `value::Base.Threads.Atomic{Int}`: The atomic value of the flag.
# Arguments:
- `value::Base.Threads.Atomic{Int}`: The atomic value of the flag.
"""
struct AtomicFlag
value::Base.Threads.Atomic{Int}
end
AtomicFlag() = AtomicFlag(Base.Threads.Atomic{Int}(0))
function set!(flag::AtomicFlag)
Base.Threads.atomic_add!(flag.value, 1)
end
function is_set(flag::AtomicFlag)::Bool
return Base.Threads.atomic_add!(flag.value, 0) > 0
end
function reset!(flag::AtomicFlag)
Base.Threads.atomic_xchg!(flag.value, 0)
end

68
src/ImageProcessing.jl
View File

@ -1,80 +1,16 @@
# src/ImageProcessing module ImageProcessing
using Images using Images
using ImageBinarization using ImageBinarization
using ImageMorphology using ImageMorphology
using ImageComponentAnalysis using ImageComponentAnalysis
using Colors # For converting to grayscale
export process_image_pipeline export process_image_pipeline
export load_and_prepare_mask # Export the new function
"""
load_and_prepare_mask(mask_path::String, target_dims::Tuple{Int, Int})
Loads a PNG image mask, converts it to a binary (Boolean) matrix,
and resizes it to the specified `target_dims`. White pixels in the mask
are considered `true` (part of the ROI), and black pixels are `false`.
# Arguments
- `mask_path`: Absolute path to the PNG mask file.
- `target_dims`: A tuple `(width, height)` representing the desired output dimensions.
# Returns
- A `BitMatrix` of `target_dims` where `true` indicates the ROI.
"""
function load_and_prepare_mask(mask_path::String, target_dims::Tuple{Int, Int})
if !isfile(mask_path)
error("Mask file not found: $(mask_path)")
end
# Load the image
img = Images.load(mask_path)
# Convert to grayscale if it's a color image
gray_img = Gray.(img)
# Binarize using Otsu's method - white regions become 'true'
binary_img = binarize(gray_img, Otsu())
# Resize to target dimensions
resized_img = imresize(binary_img, (target_dims[2], target_dims[1]))
# Ensure the output is a BitMatrix, as imresize can change the type
return resized_img .> 0.5
end
# =================================================================== # ===================================================================
# CORE PROCESSING PIPELINE # CORE PROCESSING PIPELINE
# =================================================================== # ===================================================================
"""
process_image_pipeline(gray_img; otsu_scale=1.0, noise_size_percent=0.1, hole_size_percent=0.05, smoothing=2)
Applies a multi-step image processing pipeline to a grayscale image to segment regions of interest.
The pipeline consists of:
1. **Binarization**: An adjusted Otsu's threshold is used to create a binary image.
2. **Noise Removal**: Small white regions (noise) are removed using an area opening operation.
3. **Hole Filling**: Small black regions (holes) within larger objects are filled.
4. **Edge Smoothing**: The edges of the final regions are smoothed using a morphological closing operation.
# Arguments
- `gray_img`: The input grayscale image (`Matrix{<:Gray}`).
# Keyword Arguments
- `otsu_scale`: A factor to scale the automatically determined Otsu threshold. Values > 1.0 make the threshold stricter (less white), < 1.0 make it more lenient (more white). Default: `1.0`.
- `noise_size_percent`: The percentage of the total image area used as a threshold to remove small white noise components. Default: `0.1`.
- `hole_size_percent`: The percentage of the total image area used as a threshold to fill black holes in white components. Default: `0.05`.
- `smoothing`: The size of the kernel for the final edge smoothing (closing) operation. Default: `2`.
# Returns
- A tuple containing four images representing the intermediate steps of the pipeline:
1. `binary_img`: The result of the initial binarization.
2. `noise_removed_img`: The image after noise removal.
3. `holes_filled_img`: The image after filling holes.
4. `smoothed_img`: The final smoothed image.
"""
function process_image_pipeline(gray_img; function process_image_pipeline(gray_img;
otsu_scale=1.0, otsu_scale=1.0,
noise_size_percent=0.1, noise_size_percent=0.1,
@ -107,3 +43,5 @@ function process_image_pipeline(gray_img;
# --- Return all intermediate steps for visualization --- # --- Return all intermediate steps for visualization ---
return binary_img, noise_removed_img, holes_filled_img, smoothed_img return binary_img, noise_removed_img, holes_filled_img, smoothed_img
end end
end

1564
src/MSIData.jl
View File

File diff suppressed because it is too large Load Diff

145
src/MSI_src.jl
View File

@ -6,99 +6,43 @@ module MSI_src
export OpenMSIData, export OpenMSIData,
GetSpectrum, GetSpectrum,
IterateSpectra, IterateSpectra,
ImportMzmlFile, ImportMzmlFile,
load_slices,
plot_slices,
plot_slice, plot_slice,
get_total_spectrum, get_total_spectrum,
get_average_spectrum, get_average_spectrum,
LoadMzml,
precompute_analytics, precompute_analytics,
process_spectrum, process_spectrum
generate_colorbar_image,
process_image_pipeline,
load_and_prepare_mask,
set_global_mz_range!,
get_global_mz_range,
MSIData,
_iterate_spectra_fast,
validate_spectrum,
get_mz_slice,
REGISTRY_LOCK,
SpectrumMetadata,
SpectrumAsset,
SpectrumMode,
CENTROID,
PROFILE,
MzMLSource,
ImzMLSource,
SimpleBufferPool,
get_buffer!,
release_buffer!
# Define shared registry lock # Export the public Preprocessing API
const REGISTRY_LOCK = ReentrantLock() export FeatureMatrix,
run_preprocessing_pipeline,
# Export the public preprocessing & precalculations API qc_is_empty,
export run_preprocessing_analysis, qc_is_regular,
main_precalculation, transform_intensity,
FeatureMatrix, smooth_spectrum,
Calibration, snip_baseline,
Smoothing, tic_normalize,
BaselineCorrection, pqn_normalize,
Normalization, detect_peaks_profile,
PeakPicking, align_peaks_lowess,
PeakBinning, bin_peaks,
get_masked_spectrum_indices, plot_stage_spectrum,
detect_peaks_profile_core, calculate_ppm_error,
detect_peaks_centroid_core, calculate_resolution_fwhm,
smooth_spectrum_core, analyze_mass_accuracy,
apply_baseline_correction_core, generate_qc_report,
apply_normalization_core, get_common_calibration_standards
bin_peaks_core,
PeakSelection,
PeakAlignment,
find_calibration_peaks_core,
align_peaks_lowess_core,
MutableSpectrum,
transform_intensity_core,
calibrate_spectra_core
export apply_baseline_correction,
apply_smoothing,
apply_peak_picking,
apply_calibration,
apply_peak_alignment,
apply_normalization,
apply_peak_binning,
apply_intensity_transformation,
save_feature_matrix
# Sprint 2: Streaming Pipeline API
export process_dataset!,
PipelineConfig,
StreamingStep,
normalize_inplace!,
transform_inplace!,
smooth_inplace!,
baseline_subtract_inplace!,
detect_peaks_streaming,
calibrate_inplace!
# Include all source files directly into the main module # Include all source files directly into the main module
include("BloomFilters.jl")
include("Common.jl")
include("ResourcePool.jl")
include("MSIData.jl") include("MSIData.jl")
include("ParserHelpers.jl") include("ParserHelpers.jl")
include("mzML.jl") include("mzML.jl")
include("imzML.jl") include("imzML.jl")
include("MzmlConverter.jl") include("MzmlConverter.jl")
include("Preprocessing.jl") include("Preprocessing.jl")
include("ImageProcessing.jl")
include("Precalculations.jl")
include("PreprocessingPipeline.jl")
include("StreamingKernels.jl")
include("StreamingPipeline.jl")
using Setfield # For immutable struct updates
# --- Main Entry Point --- # # --- Main Entry Point --- #
@ -109,42 +53,15 @@ using Setfield # For immutable struct updates
Opens a .mzML or .imzML file and prepares it for data access. Opens a .mzML or .imzML file and prepares it for data access.
This is the main entry point for the new data access API. This is the main entry point for the new data access API.
# Arguments
- `filepath::String`: Path to the .mzML or .imzML file to open.
- `cache_size::Int`: Number of spectra to cache in memory.
- `spectrum_type_map::Union{Dict{Int, Symbol}, Nothing}`: Optional mapping of spectrum indices to types.
# Returns
- An `MSIData` object.
""" """
function OpenMSIData(filepath::String; cache_size=300, spectrum_type_map::Union{Dict{Int, Symbol}, Nothing}=nothing) function OpenMSIData(filepath::String; cache_size=300)
# Apply standard path normalization for cross-platform compatibility if endswith(lowercase(filepath), ".mzml")
# Ensure Windows backslashes are converted to OS-native separators return load_mzml_lazy(filepath, cache_size=cache_size)
norm_filepath = normpath(replace(filepath, "\\" => "/")) elseif endswith(lowercase(filepath), ".imzml")
return load_imzml_lazy(filepath, cache_size=cache_size)
local msi_data
if endswith(lowercase(norm_filepath), ".mzml")
msi_data = load_mzml_lazy(norm_filepath, cache_size=cache_size)
elseif endswith(lowercase(norm_filepath), ".imzml")
msi_data = load_imzml_lazy(norm_filepath, cache_size=cache_size)
else else
error("Unsupported file type: $norm_filepath. Please provide a .mzML or .imzML file.") error("Unsupported file type: $filepath. Please provide a .mzML or .imzML file.")
end end
# Apply spectrum type map if provided
if spectrum_type_map !== nothing
for (idx, type_symbol) in spectrum_type_map
if 1 <= idx <= length(msi_data.spectra_metadata)
msi_data.spectra_metadata[idx] = @set msi_data.spectra_metadata[idx].type = type_symbol
else
@warn "Spectrum index $idx out of bounds for spectrum_type_map. Skipping."
end
end
end
return msi_data
end end
end # module MSI_src end # module MSI_src

257
src/MzmlConverter.jl
View File

@ -5,7 +5,7 @@ into a proper .imzML/.ibd file pair, using a separate synchronization file.
It replicates the functionality of the original R scripts that use MALDIquant. It replicates the functionality of the original R scripts that use MALDIquant.
""" """
using DataFrames, Printf, CSV, UUIDs, ProgressMeter using DataFrames, Printf, CSV
# This file assumes that the main application file (e.g., app.jl) has already included # This file assumes that the main application file (e.g., app.jl) has already included
# the necessary source files: MSIData.jl, mzML.jl, imzML.jl # the necessary source files: MSIData.jl, mzML.jl, imzML.jl
@ -396,7 +396,6 @@ This function handles three cases:
- A tuple `(mz_array, intensity_array)` for the rendered pixel spectrum. - A tuple `(mz_array, intensity_array)` for the rendered pixel spectrum.
""" """
function RenderPixel( function RenderPixel(
intensity_buffer::Vector{Float32},
pixel_info::AbstractVector{Int64}, pixel_info::AbstractVector{Int64},
scans::AbstractMatrix{Int64}, scans::AbstractMatrix{Int64},
msi_data::MSIData, msi_data::MSIData,
@ -409,32 +408,19 @@ function RenderPixel(
num_actions = last_scan - first_scan num_actions = last_scan - first_scan
# Get reference m/z array from the first scan involved. # Get reference m/z array from the first scan involved.
# This call remains type-unstable, but its impact is now isolated and only paid once.
mz_array, _ = GetSpectrum(msi_data, first_scan) mz_array, _ = GetSpectrum(msi_data, first_scan)
# If the first spectrum was empty, we can't do anything else. # If the first spectrum was empty, we can't do anything else.
if isempty(mz_array) if isempty(mz_array)
return (mz_array, view(intensity_buffer, 0:-1), UNKNOWN) return (mz_array, Float32[])
end end
# Determine the mode for the output spectrum. If any contributing scan is profile, the result is profile. new_intensity = zeros(Float32, length(mz_array))
final_mode = CENTROID
for i in first_scan:last_scan
if msi_data.spectra_metadata[i].mode == PROFILE
final_mode = PROFILE
break
end
end
# Use the provided buffer instead of allocating a new one.
if length(intensity_buffer) < length(mz_array)
error("Provided intensity buffer is too small for spectrum of length $(length(mz_array))")
end
new_intensity = view(intensity_buffer, 1:length(mz_array))
fill!(new_intensity, 0.0f0)
# SAFETY: Ensure we have valid scan indices # SAFETY: Ensure we have valid scan indices
if first_scan < 1 || last_scan > size(scans, 1) || first_scan > last_scan if first_scan < 1 || last_scan > size(scans, 1) || first_scan > last_scan
return (mz_array, new_intensity, final_mode) # Return zero intensity for invalid ranges return (mz_array, new_intensity) # Return zero intensity for invalid ranges
end end
if num_actions == 0 # Single scan contributes to the pixel if num_actions == 0 # Single scan contributes to the pixel
@ -483,7 +469,7 @@ function RenderPixel(
# FINAL SAFETY: Clamp any negative values to zero # FINAL SAFETY: Clamp any negative values to zero
new_intensity = max.(new_intensity, 0.0f0) new_intensity = max.(new_intensity, 0.0f0)
return (mz_array, new_intensity, final_mode) return (mz_array, new_intensity)
end end
""" """
@ -512,96 +498,84 @@ function ConvertMzmlToImzml(source_file::String, target_ibd_file::String, timing
open(target_ibd_file, "w") do ibd_stream open(target_ibd_file, "w") do ibd_stream
write(ibd_stream, zeros(UInt8, 16)) # UUID placeholder write(ibd_stream, zeros(UInt8, 16)) # UUID placeholder
end end
return BinaryMetadata[], Tuple{Int, Int}[], (0, 0), SpectrumMode[], uuid4() return BinaryMetadata[], Tuple{Int, Int}[], (0, 0), UNKNOWN
end end
width = maximum(timing_matrix[:, 1]) width = maximum(timing_matrix[:, 1])
height = maximum(timing_matrix[:, 2]) height = maximum(timing_matrix[:, 2])
msi_data = OpenMSIData(source_file) msi_data = OpenMSIData(source_file)
local binary_meta_vec, coords_vec, pixel_modes, ibd_uuid
try
precompute_analytics(msi_data)
max_points = isempty(msi_data.spectrum_stats_df.NumPoints) ? 0 : maximum(msi_data.spectrum_stats_df.NumPoints)
intensity_buffer = zeros(Float32, max_points)
scan_time_deltas = zeros(Int64, size(scans, 1)) source_mode = UNKNOWN
if size(scans, 1) > 1 if !isempty(msi_data.spectra_metadata)
for i in 1:(size(scans, 1) - 1) source_mode = msi_data.spectra_metadata[1].mode
delta = scans[i+1, 2] - scans[i, 2]
scan_time_deltas[i] = max(1, delta)
end
scan_time_deltas[end] = max(1, scan_time_deltas[end-1])
end
pixel_time_deltas = zeros(Int64, size(timing_matrix, 1))
if size(timing_matrix, 1) > 1
for i in 1:(size(timing_matrix, 1) - 1)
delta = timing_matrix[i+1, 3] - timing_matrix[i, 3]
pixel_time_deltas[i] = max(1, delta)
end
pixel_time_deltas[end] = max(1, pixel_time_deltas[end-1])
end
binary_meta_vec = BinaryMetadata[]
sizehint!(binary_meta_vec, size(timing_matrix, 1))
coords_vec = Tuple{Int, Int}[]
sizehint!(coords_vec, size(timing_matrix, 1))
pixel_modes = SpectrumMode[]
sizehint!(pixel_modes, size(timing_matrix, 1))
empty_pixel_count = 0
open(target_ibd_file, "w") do ibd_stream
# Generate and write a valid UUID
ibd_uuid = uuid4()
write(ibd_stream, htol(ibd_uuid.value))
p = Progress(size(timing_matrix, 1), 1, "Converting pixels... ")
for i in 1:size(timing_matrix, 1)
pixel_info = timing_matrix[i, :]
x, y = pixel_info[1], pixel_info[2]
push!(coords_vec, (x, y))
first_scan = pixel_info[4]
last_scan = pixel_info[5]
if first_scan > last_scan || first_scan < 1 || last_scan > size(scans, 1)
empty_pixel_count += 1
current_pos = position(ibd_stream)
push!(binary_meta_vec, BinaryMetadata(current_pos, 0, current_pos, 0))
push!(pixel_modes, UNKNOWN) # Mode for empty pixel
next!(p)
continue
end
mz, intensity, pixel_mode = RenderPixel(intensity_buffer, pixel_info, scans, msi_data, scan_time_deltas, pixel_time_deltas)
push!(pixel_modes, pixel_mode)
# Write m/z array (as Float64)
mz_offset = position(ibd_stream)
write(ibd_stream, htol.(mz)) # mz is Vector{Float64}
mz_length = position(ibd_stream) - mz_offset
# Write intensity array (as Float32)
int_offset = position(ibd_stream)
write(ibd_stream, htol.(intensity)) # intensity is Vector{Float32}
int_length = position(ibd_stream) - int_offset
push!(binary_meta_vec, BinaryMetadata(mz_offset, mz_length, int_offset, int_length))
next!(p)
end
end
@info "Found and processed $empty_pixel_count empty pixels out of $(size(timing_matrix, 1)) total."
finally
close(msi_data)
end end
return binary_meta_vec, coords_vec, (width, height), pixel_modes, ibd_uuid, msi_data.instrument_metadata scan_time_deltas = zeros(Int64, size(scans, 1))
if size(scans, 1) > 1
for i in 1:(size(scans, 1) - 1)
delta = scans[i+1, 2] - scans[i, 2]
scan_time_deltas[i] = max(1, delta)
end
scan_time_deltas[end] = max(1, scan_time_deltas[end-1])
end
pixel_time_deltas = zeros(Int64, size(timing_matrix, 1))
if size(timing_matrix, 1) > 1
for i in 1:(size(timing_matrix, 1) - 1)
delta = timing_matrix[i+1, 3] - timing_matrix[i, 3]
pixel_time_deltas[i] = max(1, delta)
end
pixel_time_deltas[end] = max(1, pixel_time_deltas[end-1])
end
binary_meta_vec = BinaryMetadata[]
sizehint!(binary_meta_vec, size(timing_matrix, 1))
coords_vec = Tuple{Int, Int}[]
sizehint!(coords_vec, size(timing_matrix, 1))
empty_pixel_count = 0
open(target_ibd_file, "w") do ibd_stream
write(ibd_stream, zeros(UInt8, 16)) # UUID placeholder
for i in 1:size(timing_matrix, 1)
pixel_info = timing_matrix[i, :]
x, y = pixel_info[1], pixel_info[2]
push!(coords_vec, (x, y))
first_scan = pixel_info[4]
last_scan = pixel_info[5]
if first_scan > last_scan || first_scan < 1 || last_scan > size(scans, 1)
empty_pixel_count += 1
# For empty pixels, offsets point to the current end of file, with zero length
current_pos = position(ibd_stream)
push!(binary_meta_vec, BinaryMetadata(current_pos, 0, current_pos, 0))
continue
end
mz, intensity = RenderPixel(pixel_info, scans, msi_data, scan_time_deltas, pixel_time_deltas)
# Write m/z array
mz_offset = position(ibd_stream)
for val in mz
write(ibd_stream, htol(Float32(val)))
end
mz_length = position(ibd_stream) - mz_offset
# Write intensity array
int_offset = position(ibd_stream)
for val in intensity
write(ibd_stream, htol(Float32(val)))
end
int_length = position(ibd_stream) - int_offset
push!(binary_meta_vec, BinaryMetadata(mz_offset, mz_length, int_offset, int_length))
end
end
@info "Found and processed $empty_pixel_count empty pixels out of $(size(timing_matrix, 1)) total."
return binary_meta_vec, coords_vec, (width, height), source_mode
end end
""" """
@ -621,7 +595,7 @@ experiment and data format.
# Returns # Returns
- `true` on success, `false` on failure. - `true` on success, `false` on failure.
""" """
function ExportImzml(target_file::String, binary_meta::Vector{BinaryMetadata}, coords::Vector{Tuple{Int, Int}}, dims::Tuple{Int, Int}, modes::Vector{SpectrumMode}, ibd_uuid::UUID, instrument_meta::InstrumentMetadata) function ExportImzml(target_file::String, binary_meta::Vector{BinaryMetadata}, coords::Vector{Tuple{Int, Int}}, dims::Tuple{Int, Int}, mode::SpectrumMode)
ibd_file = replace(target_file, r"\.imzML$"i => ".ibd") ibd_file = replace(target_file, r"\.imzML$"i => ".ibd")
if isempty(binary_meta) if isempty(binary_meta)
@ -651,16 +625,6 @@ function ExportImzml(target_file::String, binary_meta::Vector{BinaryMetadata}, c
<cvParam cvRef="MS" accession="MS:1000579" name="MS1 spectrum"/> <cvParam cvRef="MS" accession="MS:1000579" name="MS1 spectrum"/>
<cvParam cvRef="IMS" accession="IMS:1000080" name="mass spectrum"/> <cvParam cvRef="IMS" accession="IMS:1000080" name="mass spectrum"/>
<cvParam cvRef="IMS" accession="IMS:1000031" name="processed"/> <cvParam cvRef="IMS" accession="IMS:1000031" name="processed"/>
<cvParam cvRef="IMS" accession="IMS:1000081" name="ibd uuid" value="$(ibd_uuid)"/>
$(
if instrument_meta.polarity == :positive
""" <cvParam cvRef="MS" accession="MS:1000130" name="positive scan"/>"""
elseif instrument_meta.polarity == :negative
""" <cvParam cvRef="MS" accession="MS:1000129" name="negative scan"/>"""
else
""
end
)
</fileContent> </fileContent>
</fileDescription> </fileDescription>
""") """)
@ -668,7 +632,7 @@ $(
<referenceableParamGroup id="mzArray"> <referenceableParamGroup id="mzArray">
<cvParam cvRef="MS" accession="MS:1000576" name="no compression"/> <cvParam cvRef="MS" accession="MS:1000576" name="no compression"/>
<cvParam cvRef="MS" accession="MS:1000514" name="m/z array" unitCvRef="MS" unitAccession="MS:1000040" unitName="m/z"/> <cvParam cvRef="MS" accession="MS:1000514" name="m/z array" unitCvRef="MS" unitAccession="MS:1000040" unitName="m/z"/>
<cvParam cvRef="MS" accession="MS:1000523" name="64-bit float"/> <cvParam cvRef="MS" accession="MS:1000521" name="32-bit float"/>
</referenceableParamGroup> </referenceableParamGroup>
<referenceableParamGroup id="intensityArray"> <referenceableParamGroup id="intensityArray">
<cvParam cvRef="MS" accession="MS:1000576" name="no compression"/> <cvParam cvRef="MS" accession="MS:1000576" name="no compression"/>
@ -695,42 +659,17 @@ $(
<cvParam cvRef="IMS" accession="IMS:1000043" name="max count of pixel y" value="$(dims[2])"/> <cvParam cvRef="IMS" accession="IMS:1000043" name="max count of pixel y" value="$(dims[2])"/>
<cvParam cvRef="IMS" accession="IMS:1000046" name="pixel size x" value="1" unitCvRef="UO" unitAccession="UO:0000017" unitName="micrometer"/> <cvParam cvRef="IMS" accession="IMS:1000046" name="pixel size x" value="1" unitCvRef="UO" unitAccession="UO:0000017" unitName="micrometer"/>
<cvParam cvRef="IMS" accession="IMS:1000047" name="pixel size y" value="1" unitCvRef="UO" unitAccession="UO:0000017" unitName="micrometer"/> <cvParam cvRef="IMS" accession="IMS:1000047" name="pixel size y" value="1" unitCvRef="UO" unitAccession="UO:0000017" unitName="micrometer"/>
<cvParam cvRef="IMS" accession="IMS:1000092" name="line scan sequence" value="line scan"/>
<cvParam cvRef="IMS" accession="IMS:1000093" name="spotsize" value="1"/>
</scanSettings> </scanSettings>
</scanSettingsList> </scanSettingsList>
""") """)
write(imzml_stream, """ <instrumentConfigurationList count="1"> write(imzml_stream, """ <instrumentConfigurationList count="1">
<instrumentConfiguration id="instrument1"> <instrumentConfiguration id="instrument1">
$(
if instrument_meta.instrument_model != "Not Available"
""" <cvParam cvRef="MS" accession="MS:1000031" name="instrument model" value="$(instrument_meta.instrument_model)"/>"""
else
""
end
)
$(
if instrument_meta.resolution !== nothing
""" <cvParam cvRef="MS" accession="MS:1001496" name="mass resolving power" value="$(instrument_meta.resolution)"/>"""
else
""
end
)
<componentList count="3"> <componentList count="3">
<source order="1"> <source order="1">
<cvParam cvRef="MS" accession="MS:1000075" name="MALDI source"/> <cvParam cvRef="MS" accession="MS:1000075" name="MALDI source"/>
</source> </source>
<analyzer order="2"> <analyzer order="2">
<cvParam cvRef="MS" accession="MS:1000084" name="time-of-flight"/> <cvParam cvRef="MS" accession="MS:1000084" name="time-of-flight"/>
$(
if instrument_meta.acquisition_mode == :centroid
""" <cvParam cvRef="MS" accession="MS:1000127" name="centroid spectrum"/>"""
elseif instrument_meta.acquisition_mode == :profile
""" <cvParam cvRef="MS" accession="MS:1000128" name="profile spectrum"/>"""
else
""
end
)
</analyzer> </analyzer>
<detector order="3"> <detector order="3">
<cvParam cvRef="MS" accession="MS:1000253" name="electron multiplier"/> <cvParam cvRef="MS" accession="MS:1000253" name="electron multiplier"/>
@ -740,23 +679,11 @@ $(
</instrumentConfiguration> </instrumentConfiguration>
</instrumentConfigurationList> </instrumentConfigurationList>
""") """)
write(imzml_stream, """ <dataProcessingList count="$(1 + (instrument_meta.vendor_preprocessing_steps !== nothing ? length(instrument_meta.vendor_preprocessing_steps) : 0))"> write(imzml_stream, """ <dataProcessingList count="1">
<dataProcessing id="conversionProcessing"> <dataProcessing id="conversionProcessing">
<processingMethod order="1" softwareRef="MSIConverter"> <processingMethod order="1" softwareRef="MSIConverter">
<cvParam cvRef="MS" accession="MS:1000544" name="Conversion to imzML"/> <cvParam cvRef="MS" accession="MS:1000544" name="Conversion to imzML"/>
</processingMethod> </processingMethod>
$(
if instrument_meta.vendor_preprocessing_steps !== nothing
join([
""" <processingMethod order="$(i + 1)" softwareRef="MSIConverter">
<cvParam cvRef="MS" accession="MS:1000589" name="data processing" value="$(step)"/>
</processingMethod>"""
for (i, step) in enumerate(instrument_meta.vendor_preprocessing_steps)
], "\n")
else
""
end
)
</dataProcessing> </dataProcessing>
</dataProcessingList> </dataProcessingList>
""") """)
@ -775,16 +702,16 @@ $(
push!(spectrum_offsets, spectrum_start) push!(spectrum_offsets, spectrum_start)
# Calculate number of points from byte length # Calculate number of points from byte length
mz_points = meta.mz_length ÷ sizeof(Float64) mz_points = meta.mz_length ÷ sizeof(Float32)
int_points = meta.int_length ÷ sizeof(Float32) int_points = meta.int_length ÷ sizeof(Float32)
write(imzml_stream, """ <spectrum id="Scan=$(i)" defaultArrayLength="$(mz_points)" index="$(i-1)"> write(imzml_stream, """ <spectrum id="Scan=$(i)" defaultArrayLength="$(mz_points)" index="$(i-1)">
<cvParam cvRef="MS" accession="MS:1000511" name="ms level" value="1"/> <cvParam cvRef="MS" accession="MS:1000511" name="ms level" value="1"/>
""") """)
if modes[i] == CENTROID if mode == CENTROID
write(imzml_stream, """ <cvParam cvRef="MS" accession="MS:1000127" name="centroid spectrum"/> write(imzml_stream, """ <cvParam cvRef="MS" accession="MS:1000127" name="centroid spectrum"/>
""") """)
elseif modes[i] == PROFILE else
write(imzml_stream, """ <cvParam cvRef="MS" accession="MS:1000128" name="profile spectrum"/> write(imzml_stream, """ <cvParam cvRef="MS" accession="MS:1000128" name="profile spectrum"/>
""") """)
end end
@ -842,22 +769,8 @@ Main workflow function to convert a .mzML file to an .imzML file.
* `target_file`: Path for the output .imzML file (the .ibd will be named accordingly). * `target_file`: Path for the output .imzML file (the .ibd will be named accordingly).
* `img_width`: width dimention for the creation of the x axis * `img_width`: width dimention for the creation of the x axis
* `img_height`: height dimention for the creation of the y axis * `img_height`: height dimention for the creation of the y axis
# Returns
- `true` if the conversion was successful.
- `false` if the conversion failed.
# Example
```julia
ImportMzmlFile("test.mzML", "test.txt", "test.imzML")
```
""" """
function ImportMzmlFile(source_file::String, sync_file::String, target_file::String; img_width::Int=0, img_height::Int=0) function ImportMzmlFile(source_file::String, sync_file::String, target_file::String; img_width::Int=0, img_height::Int=0)
if !isfile(source_file)
throw(ArgumentError("Source mzML file not found: $source_file"))
end
println("Step 1: Getting scan times from .mzML file...") println("Step 1: Getting scan times from .mzML file...")
scans = GetMzmlScanTime(source_file) scans = GetMzmlScanTime(source_file)
@ -866,13 +779,13 @@ function ImportMzmlFile(source_file::String, sync_file::String, target_file::Str
println("Step 3: Converting spectra and writing .ibd file...") println("Step 3: Converting spectra and writing .ibd file...")
ibd_file = replace(target_file, r"\.imzML$"i => ".ibd") ibd_file = replace(target_file, r"\.imzML$"i => ".ibd")
binary_meta, coords, (width, height), pixel_modes, ibd_uuid, instrument_meta = ConvertMzmlToImzml(source_file, ibd_file, timing_matrix, scans) binary_meta, coords, (width, height), source_mode = ConvertMzmlToImzml(source_file, ibd_file, timing_matrix, scans)
# Flip image vertically to match R script output # Flip image vertically to match R script output
flipped_coords = [(x, height - y + 1) for (x, y) in coords] flipped_coords = [(x, height - y + 1) for (x, y) in coords]
println("Step 4: Exporting .imzML metadata file...") println("Step 4: Exporting .imzML metadata file...")
success = ExportImzml(target_file, binary_meta, flipped_coords, (width, height), pixel_modes, ibd_uuid, instrument_meta) success = ExportImzml(target_file, binary_meta, flipped_coords, (width, height), source_mode)
if success if success
println("Conversion successful: $target_file") println("Conversion successful: $target_file")

28
src/ParserHelpers.jl
View File

@ -28,18 +28,8 @@ end
Reads a stream line-by-line until a line matches the provided regex. Reads a stream line-by-line until a line matches the provided regex.
# Arguments
- `stream::IO`: The stream to parse the metadata from.
- `regex::Regex`: The regex to match.
# Returns # Returns
- A `RegexMatch` object if a match is found, otherwise throws an error. - A `RegexMatch` object if a match is found, otherwise throws an error.
# Example
```julia
find_tag(open("test.mzML"), r"<spectrum")
```
""" """
function find_tag(stream, regex::Regex) function find_tag(stream, regex::Regex)
while !eof(stream) while !eof(stream)
@ -57,18 +47,8 @@ end
Retrieves an attribute's value from an XML tag string. Retrieves an attribute's value from an XML tag string.
# Arguments
- `source::String`: The string to parse the attribute from.
- `tag::String`: The tag to match.
# Returns # Returns
- A `RegexMatch` object containing the attribute and its value. - A `RegexMatch` object containing the attribute and its value.
# Example
```julia
get_attribute("<spectrum index=\"1\">", "index")
```
""" """
function get_attribute(source::AbstractString, tag::String = "([^=]+)") function get_attribute(source::AbstractString, tag::String = "([^=]+)")
# Construct the regex pattern string # Construct the regex pattern string
@ -85,16 +65,10 @@ struct. It reads accessions to determine the data format (e.g., `Float32`),
compression status (`zlib`), and axis type (m/z vs. intensity). compression status (`zlib`), and axis type (m/z vs. intensity).
# Arguments # Arguments
- `stream::IO`: An IO stream positioned at the start of the `cvParam` block. - `stream`: An IO stream positioned at the start of the `cvParam` block.
# Returns # Returns
- A `SpecDim` struct populated with the parsed configuration. - A `SpecDim` struct populated with the parsed configuration.
# Example
```julia
configure_spec_dim(open("test.mzML"))
```
""" """
function configure_spec_dim(stream) function configure_spec_dim(stream)
axis = SpecDim(Float64, false, 1, 0, UNKNOWN) # Add UNKNOWN as default mode axis = SpecDim(Float64, false, 1, 0, UNKNOWN) # Add UNKNOWN as default mode

1794
src/Precalculations.jl
View File

File diff suppressed because it is too large Load Diff

1818
src/Preprocessing.jl
View File

File diff suppressed because it is too large Load Diff

510
src/PreprocessingPipeline.jl
View File

@ -1,510 +0,0 @@
# src/PreprocessingPipeline.jl
using Base.Threads # For multithreading
using Printf # For @sprintf
using Interpolations # For linear_interpolation
using DataFrames # For saving feature matrix
using CSV # For saving feature matrix
# This file provides a set of functions that apply preprocessing steps to a vector of
# `MutableSpectrum` objects. Each function takes the vector of spectra and a dictionary
# of parameters, modifying the spectra in-place where appropriate. This mirrors the
# logic from `test/run_preprocessing.jl` but is intended for use in the main application.
# ===================================================================
# PREPROCESSING PIPELINE FUNCTIONS (IN-PLACE)
# ===================================================================
"""
apply_baseline_correction(spectra::Vector{MutableSpectrum}, params::Dict)
Applies baseline correction to the intensity data of each spectrum. This function
modifies the `.intensity` field of each `MutableSpectrum` object in-place.
# Parameters from `params` Dict:
- `:method` (Symbol): The algorithm to use. Supports `:snip`, `:convex_hull`, `:median`. Defaults to `:snip`.
- `:iterations` (Int): The number of iterations for the SNIP algorithm. Defaults to 100.
- `:window` (Int): The window size for the Median algorithm. Defaults to 20.
"""
function apply_baseline_correction(spectra::Vector{MutableSpectrum}, params::Dict)
method = get(params, :method, :snip)
iterations = get(params, :iterations, 100)
window = get(params, :window, 20)
Threads.@threads for s in spectra
if validate_spectrum(s.mz, s.intensity)
original_length = length(s.intensity)
baseline = apply_baseline_correction_core(s.intensity; method=method, iterations=iterations, window=window)
# CRITICAL: Ensure baseline has the same length as intensity
if length(baseline) != original_length
@warn "Baseline correction: length mismatch for spectrum $(s.id). baseline=$(length(baseline)), intensity=$original_length. Skipping this spectrum."
continue
end
s.intensity = max.(0.0, s.intensity .- baseline)
end
end
end
"""
apply_intensity_transformation(spectra::Vector{MutableSpectrum}, params::Dict)
Applies an intensity transformation to the intensity data of each spectrum. This function
modifies the `.intensity` field of each `MutableSpectrum` object in-place.
# Parameters from `params` Dict:
- `:method` (Symbol): The transformation to apply. Supports `:sqrt`, `:log`, `:log2`, `:log10`, `:log1p`. Defaults to `:sqrt`.
"""
function apply_intensity_transformation(spectra::Vector{MutableSpectrum}, params::Dict)
method = get(params, :method, :sqrt)
Threads.@threads for s in spectra
if validate_spectrum(s.mz, s.intensity)
original_length = length(s.intensity)
transformed = transform_intensity_core(s.intensity; method=method)
# CRITICAL: Ensure transformation preserves array length
if length(transformed) != original_length
@warn "Intensity transformation: length mismatch for spectrum $(s.id). transformed=$(length(transformed)), original=$original_length. Skipping this spectrum."
continue
end
s.intensity = transformed
end
end
end
"""
apply_smoothing(spectra::Vector{MutableSpectrum}, params::Dict)
Applies a smoothing filter to the intensity data of each spectrum. This function
modifies the `.intensity` field of each `MutableSpectrum` object in-place.
# Parameters from `params` Dict:
- `:method` (Symbol): The smoothing algorithm. Supports `:savitzky_golay`, `:moving_average`. Defaults to `:savitzky_golay`.
- `:window` (Int): The size of the smoothing window. Defaults to 9.
- `:order` (Int): The polynomial order for the Savitzky-Golay filter. Defaults to 2.
"""
function apply_smoothing(spectra::Vector{MutableSpectrum}, params::Dict)
method = get(params, :method, :savitzky_golay)
window = get(params, :window, 9)
order = get(params, :order, 2)
Threads.@threads for s in spectra
if validate_spectrum(s.mz, s.intensity)
original_length = length(s.intensity)
smoothed_intensity = smooth_spectrum_core(s.intensity; method=method, window=window, order=order)
# CRITICAL: Ensure smoothing preserves array length
if length(smoothed_intensity) != original_length
@warn "Smoothing: length mismatch for spectrum $(s.id). smoothed=$(length(smoothed_intensity)), original=$original_length, method=$method, window=$window. Skipping this spectrum."
continue
end
s.intensity = max.(0.0, smoothed_intensity)
end
end
end
"""
apply_peak_picking(spectra::Vector{MutableSpectrum}, params::Dict)
Detects peaks in each spectrum and stores them in the `.peaks` field of each
`MutableSpectrum` object, modifying it in-place.
# Parameters from `params` Dict:
- `:method` (Symbol): The peak detection algorithm. Supports `:profile`, `:wavelet`, `:centroid`. Defaults to `:profile`.
- `:snr_threshold` (Float64): Signal-to-Noise Ratio threshold.
- `:half_window` (Int): Half-window size for local maxima detection.
- `:min_peak_prominence` (Float64): Minimum required prominence for a peak.
- `:merge_peaks_tolerance` (Float64): m/z tolerance to merge adjacent peaks.
"""
function apply_peak_picking(spectra::Vector{MutableSpectrum}, params::Dict)
method = get(params, :method, :profile)
snr_threshold = get(params, :snr_threshold, 3.0)
half_window = get(params, :half_window, 10)
min_peak_prominence = get(params, :min_peak_prominence, 0.1)
merge_peaks_tolerance = get(params, :merge_peaks_tolerance, 0.002)
Threads.@threads for s in spectra
if validate_spectrum(s.mz, s.intensity)
# old_len = length(s.peaks)
if method == :profile
s.peaks = detect_peaks_profile_core(s.mz, s.intensity; snr_threshold=snr_threshold, half_window=half_window, min_peak_prominence=min_peak_prominence, merge_peaks_tolerance=merge_peaks_tolerance)
elseif method == :wavelet
s.peaks = detect_peaks_wavelet_core(s.mz, s.intensity; snr_threshold=snr_threshold, half_window=half_window)
elseif method == :centroid
s.peaks = detect_peaks_centroid_core(s.mz, s.intensity; snr_threshold=snr_threshold)
else
s.peaks = detect_peaks_profile_core(s.mz, s.intensity; snr_threshold=snr_threshold, half_window=half_window)
end
# @info "Spectrum $(s.id): detected $(length(s.peaks)) peaks"
else
s.peaks = []
end
end
end
"""
apply_peak_selection(spectra::Vector{MutableSpectrum}, params::Dict)
Filters peaks within each spectrum based on quality criteria. This step removes peaks
that do not meet the specified thresholds for signal-to-noise ratio (SNR),
full width at half maximum (FWHM), and peak shape.
This function modifies the `.peaks` field of each `MutableSpectrum` object in the `spectra` vector in-place.
# Parameters from `params` Dict:
- `:min_snr` (Float64): The minimum Signal-to-Noise Ratio required for a peak to be kept.
- `:min_fwhm_ppm` (Float64): The minimum FWHM (in ppm) for a peak.
- `:max_fwhm_ppm` (Float64): The maximum FWHM (in ppm) for a peak.
- `:min_shape_r2` (Float64): The minimum R² value from a Gaussian fit, measuring peak shape quality.
"""
function apply_peak_selection(spectra::Vector{MutableSpectrum}, params::Dict)
min_snr = get(params, :min_snr, 0.0)
min_fwhm = get(params, :min_fwhm_ppm, 0.0)
max_fwhm = get(params, :max_fwhm_ppm, Inf)
min_r2 = get(params, :min_shape_r2, 0.0)
# Handle `nothing` from params, which can happen if pre-calculation fails.
min_snr = isnothing(min_snr) ? 0.0 : min_snr
min_fwhm = isnothing(min_fwhm) ? 0.0 : min_fwhm
max_fwhm = isnothing(max_fwhm) ? Inf : max_fwhm
min_r2 = isnothing(min_r2) ? 0.0 : min_r2
Threads.@threads for s in spectra
if !isempty(s.peaks)
filter!(p ->
p.snr >= min_snr &&
(min_fwhm <= p.fwhm <= max_fwhm) &&
p.shape_r2 >= min_r2,
s.peaks
)
end
end
end
"""
apply_calibration(spectra::Vector{MutableSpectrum}, params::Dict, reference_peaks::Dict)
Performs mass calibration on each spectrum using a list of internal standards.
This function modifies the `.mz` axis of each `MutableSpectrum` object in-place.
# Parameters from `params` Dict:
- `:method` (Symbol): The calibration method. Only `:internal_standards` is currently meaningful.
- `:ppm_tolerance` (Float64): The tolerance in PPM for matching detected peaks to reference masses.
- `:fit_order` (Int): The polynomial order for the calibration fit (not yet used in this implementation, defaults to linear).
"""
function apply_calibration(spectra::Vector{MutableSpectrum}, params::Dict, reference_peaks::Dict)
method = get(params, :method, :none)
ppm_tolerance = get(params, :ppm_tolerance, 20.0)
if method == :none || isempty(reference_peaks)
return
end
reference_masses = collect(keys(reference_peaks))
Threads.@threads for i in 1:length(spectra)
s = spectra[i]
if validate_spectrum(s.mz, s.intensity)
original_length = length(s.mz)
matched_peaks = find_calibration_peaks_core(s.mz, s.intensity, reference_masses; ppm_tolerance=ppm_tolerance)
if length(matched_peaks) >= 2
measured = sort(collect(values(matched_peaks)))
theoretical = sort(collect(keys(matched_peaks)))
itp = linear_interpolation(measured, theoretical, extrapolation_bc=Line())
new_mz = itp(s.mz)
# CRITICAL: Ensure m/z axis preserves array length
if length(new_mz) != original_length
@warn "Calibration: length mismatch for spectrum $(s.id). new_mz=$(length(new_mz)), original=$original_length. Skipping this spectrum."
continue
end
s.mz = new_mz # Modify mz-axis in-place
else
@warn "Spectrum $(s.id): insufficient reference peaks ($(length(matched_peaks)) found), skipping calibration."
end
end
end
end
"""
apply_peak_alignment(spectra::Vector{MutableSpectrum}, params::Dict)
Aligns the m/z axis of all spectra to a chosen reference spectrum. This function
modifies both the `.mz` axis and the m/z values within the `.peaks` field of each
`MutableSpectrum` object in-place.
# Parameters from `params` Dict:
- `:method` (Symbol): The alignment algorithm. Supports `:lowess`, `:linear`, `:ransac`.
- `:tolerance` (Float64): The tolerance for matching peaks between spectra.
- `:tolerance_unit` (Symbol): The unit for tolerance, `:mz` or `:ppm`.
"""
function apply_peak_alignment(spectra::Vector{MutableSpectrum}, params::Dict)
method = get(params, :method, :none)
tolerance = get(params, :tolerance, 0.002)
tolerance_unit = get(params, :tolerance_unit, :mz)
if method == :none
return
end
ref_find_idx = findfirst(s -> !isempty(s.peaks), spectra)
if ref_find_idx === nothing
@warn "Insufficient spectra with peaks for alignment. Skipping."
return
end
ref_spectrum = spectra[ref_find_idx]
ref_peaks_mz = [p.mz for p in ref_spectrum.peaks]
Threads.@threads for s in spectra
if s.id == ref_spectrum.id || isempty(s.peaks)
continue
end
current_peaks_mz = [p.mz for p in s.peaks]
alignment_func = align_peaks_lowess_core(ref_peaks_mz, current_peaks_mz; method=method, tolerance=tolerance, tolerance_unit=tolerance_unit)
original_length = length(s.mz)
new_mz = alignment_func.(s.mz)
# CRITICAL: Ensure alignment preserves array length
if length(new_mz) != original_length
@warn "Peak alignment: m/z length mismatch for spectrum $(s.id). new_mz=$(length(new_mz)), original=$original_length. Skipping this spectrum."
continue
end
s.mz = new_mz # Update m/z axis
# Update peak m/z values
for i in 1:length(s.peaks)
old_peak = s.peaks[i]
aligned_peak_mz = alignment_func(old_peak.mz)
s.peaks[i] = (mz=aligned_peak_mz, intensity=old_peak.intensity, fwhm=old_peak.fwhm, shape_r2=old_peak.shape_r2, snr=old_peak.snr, prominence=old_peak.prominence)
end
end
end
"""
apply_normalization(spectra::Vector{MutableSpectrum}, params::Dict)
Applies intensity normalization to each spectrum. This function modifies the
`.intensity` field of each `MutableSpectrum` object in-place.
# Parameters from `params` Dict:
- `:method` (Symbol): The normalization method. Supports `:tic`, `:median`, `:rms`, `:none`.
"""
function apply_normalization(spectra::Vector{MutableSpectrum}, params::Dict)
method = get(params, :method, :tic)
Threads.@threads for s in spectra
if validate_spectrum(s.mz, s.intensity)
original_length = length(s.intensity)
normalized = apply_normalization_core(s.intensity; method=method)
# CRITICAL: Ensure normalization preserves array length
if length(normalized) != original_length
@warn "Normalization: length mismatch for spectrum $(s.id). normalized=$(length(normalized)), original=$original_length. Skipping this spectrum."
continue
end
s.intensity = normalized
end
end
end
function apply_peak_binning(spectra::Vector{MutableSpectrum}, params::Dict)
tolerance = get(params, :tolerance, 20.0)
tolerance_unit = get(params, :tolerance_unit, :ppm)
min_peak_per_bin = get(params, :min_peak_per_bin, 3)
if isempty(spectra) || all(s -> isempty(s.peaks), spectra)
@warn "No peaks found for binning. Returning empty feature matrix."
return nothing, nothing
end
all_peaks = Vector{Tuple{Float64, Float64}}()
for s in spectra
for p in s.peaks
push!(all_peaks, (p.mz, p.intensity))
end
end
if isempty(all_peaks)
@warn "No peaks collected for binning."
return nothing, nothing
end
sort!(all_peaks, by=x->x[1])
bin_centers = Float64[]
bin_intensities = Float64[]
i = 1
while i <= length(all_peaks)
current_bin_start = i
current_peak = all_peaks[i]
j = i + 1
while j <= length(all_peaks)
next_peak = all_peaks[j]
tol = (tolerance_unit == :ppm) ? (current_peak[1] * tolerance / 1e6) : tolerance
if (next_peak[1] - current_peak[1]) <= tol
j += 1
else
break
end
end
current_bin_end = j - 1
bin_size = current_bin_end - current_bin_start + 1
if bin_size >= min_peak_per_bin
bin_peaks = all_peaks[current_bin_start:current_bin_end]
mz_sum = sum(p[1] for p in bin_peaks)
intensity_sum = sum(p[2] for p in bin_peaks)
mz_center = mz_sum / bin_size
avg_intensity = intensity_sum / bin_size
push!(bin_centers, mz_center)
push!(bin_intensities, avg_intensity)
end
i = j
end
if !isempty(bin_centers)
n_bins = length(bin_centers)
feature_matrix = Matrix{Float64}(undef, 2, n_bins)
for i in 1:n_bins
feature_matrix[1, i] = bin_centers[i]
feature_matrix[2, i] = bin_intensities[i]
end
bin_info = [(bin_centers[i], bin_intensities[i]) for i in 1:n_bins]
return feature_matrix, bin_info
else
@warn "No bins created after filtering"
return nothing, nothing
end
end
"""
save_feature_matrix(feature_matrix::Matrix{Float64}, bin_info, output_dir::String) -> Tuple{String, String}
Saves the aggregated `2 x n_bins` feature matrix into two different CSV formats.
1. **Simple Format (`feature_matrix_simple.csv`):** A two-column CSV with "mz" and "intensity".
2. **Standard Format (`feature_matrix_standard.csv`):** A row-based format where m/z values are headers and there is a single data row for the aggregated spectrum.
# Arguments
- `feature_matrix::Matrix{Float64}`: The `2 x n_bins` matrix from `apply_peak_binning`.
- `bin_info`: The associated bin information (currently unused but kept for compatibility).
- `output_dir::String`: The directory where the output CSV files will be saved.
# Returns
- A tuple containing the paths to the two saved files.
"""
function save_feature_matrix(feature_matrix::Matrix{Float64}, bin_info, output_dir::String)
# Save as simple CSV with m/z and intensity rows
csv_path = joinpath(output_dir, "feature_matrix_simple.csv")
open(csv_path, "w") do io
write(io, "mz,intensity\n")
for i in 1:size(feature_matrix, 2)
mz = feature_matrix[1, i]
intensity = feature_matrix[2, i]
write(io, "$mz,$intensity\n")
end
end
@info "Saved simple feature matrix: $csv_path"
# Also save in a more standard format for MSI
csv_path_standard = joinpath(output_dir, "feature_matrix_standard.csv")
open(csv_path_standard, "w") do io
write(io, "sample_type,")
mz_headers = [@sprintf("mz_%.4f", feature_matrix[1, i]) for i in 1:size(feature_matrix, 2)]
write(io, join(mz_headers, ",") * "\n")
write(io, "aggregated_spectrum,")
intensity_values = [feature_matrix[2, i] for i in 1:size(feature_matrix, 2)]
write(io, join(string.(intensity_values), ",") * "\n")
end
@info "Saved standard format matrix: $csv_path_standard"
return csv_path, csv_path_standard
end
function execute_full_preprocessing(progress_callback::Function, # Now first positional
spectra::Vector{MutableSpectrum}, params::Dict,
pipeline_steps::Vector{String}, reference_peaks::Dict,
mask_path::Union{String, Nothing} # Positional argument
)
println("Starting preprocessing pipeline with $(length(spectra)) spectra")
println("Steps: $(join(pipeline_steps, " -> "))")
# These variables will be populated by the pipeline steps
feature_matrix = nothing
bin_definitions = nothing
# Apply pipeline steps, modifying `spectra` in-place
for step in pipeline_steps
progress_callback(step)
println("\n" * "-"^60)
println("PROCESSING STEP: $step")
println("-"^60)
if step == "stabilization"
println(" Applying intensity transformation (stabilization)")
apply_intensity_transformation(spectra, get(params, :Stabilization, Dict()))
elseif step == "baseline_correction"
println(" Applying baseline correction")
apply_baseline_correction(spectra, get(params, :BaselineCorrection, Dict()))
elseif step == "smoothing"
println(" Applying smoothing")
apply_smoothing(spectra, get(params, :Smoothing, Dict()))
elseif step == "peak_picking"
println(" Applying peak picking")
apply_peak_picking(spectra, get(params, :PeakPicking, Dict()))
elseif step == "peak_selection"
println(" Applying peak selection")
apply_peak_selection(spectra, get(params, :PeakSelection, Dict()))
elseif step == "calibration"
println(" Applying calibration")
apply_calibration(spectra, get(params, :Calibration, Dict()), reference_peaks)
elseif step == "peak_alignment"
println(" Applying peak alignment")
apply_peak_alignment(spectra, get(params, :PeakAlignment, Dict()))
elseif step == "normalization"
println(" Applying normalization")
apply_normalization(spectra, get(params, :Normalization, Dict()))
elseif step == "peak_binning"
println(" Applying peak binning")
feature_matrix, bin_definitions = apply_peak_binning(spectra, get(params, :PeakBinning, Dict()))
else
@warn "Unknown step: $step, skipping"
end
println("✓ Completed step: $step")
end
return feature_matrix, bin_definitions
end

79
src/ResourcePool.jl
View File

@ -1,79 +0,0 @@
# src/ResourcePool.jl
using Base.Threads
"""
ResourcePool{T}
A thread-safe pool for reusing objects of type `T` to minimize allocations and GC pressure.
Specifically designed for high-performance computing tasks where large buffers are needed
repeatedly across multiple threads.
# Fields:
- `pool::Vector{T}`: The underlying storage for idle resources.
- `lock::ReentrantLock`: Ensures thread-safe access to the pool.
- `max_size::Int`: Maximum number of resources to hold in the pool.
- `constructor::Function`: A function to create a new resource if the pool is empty.
"""
mutable struct ResourcePool{T}
pool::Vector{T}
lock::ReentrantLock
max_size::Int
constructor::Function
end
"""
ResourcePool{T}(constructor::Function; max_size::Int=2 * nthreads())
Creates a new `ResourcePool` for resources of type `T`.
"""
function ResourcePool{T}(constructor::Function; max_size::Int=2 * nthreads()) where T
return ResourcePool{T}(T[], ReentrantLock(), max_size, constructor)
end
"""
acquire(pool::ResourcePool{T}) -> T
Retrieves a resource from the pool. If the pool is empty, a new resource is created
using the constructor.
"""
function acquire(pool::ResourcePool{T}) where T
lock(pool.lock) do
if !isempty(pool.pool)
return pop!(pool.pool)
end
end
# Create new resource outside of lock to minimize contention
return pool.constructor()
end
"""
release!(pool::ResourcePool{T}, resource::T)
Returns a resource to the pool for later reuse. If the pool is already at `max_size`,
the resource is allowed to be garbage collected.
"""
function release!(pool::ResourcePool{T}, resource::T) where T
lock(pool.lock) do
if length(pool.pool) < pool.max_size
push!(pool.pool, resource)
end
end
return nothing
end
"""
with_resource(f::Function, pool::ResourcePool{T})
Acquires a resource from the pool, executes the function `f(resource)`, and
automatically releases the resource back to the pool when finished.
"""
function with_resource(f::Function, pool::ResourcePool{T}) where T
resource = acquire(pool)
try
return f(resource)
finally
release!(pool, resource)
end
end

317
src/StreamingKernels.jl
View File

@ -1,317 +0,0 @@
# src/StreamingKernels.jl
# ============================================================================
# In-Place Spectral Kernels for the Streaming Pipeline
#
# These functions operate on raw (mz, intensity) views from the Sprint 1
# Mmap engine. They write results back to the input buffers using .= to
# achieve zero-allocation processing per spectrum.
#
# Design contract:
# - All !-suffixed functions modify their arguments in-place
# - If a kernel needs temporary storage, it borrows from data.resource_pool
# - No function creates MutableSpectrum objects
# ============================================================================
using Statistics: mean, median
# =============================================================================
# Category A: Purely Streamable Kernels
# =============================================================================
"""
normalize_inplace!(intensity::AbstractVector{<:Real}, method::Symbol)
Normalizes intensity values in-place. Supports :tic, :median, :rms.
Zero-allocation for the normalization itself.
"""
@inline function normalize_inplace!(intensity::AbstractVector{<:Real}, method::Symbol)
if method === :tic
s = sum(intensity)
if s > 0
intensity ./= s
end
elseif method === :median
m = median(intensity)
if m > 0
intensity ./= m
end
elseif method === :rms
s = sqrt(sum(abs2, intensity) / length(intensity))
if s > 0
intensity ./= s
end
end
return intensity
end
"""
transform_inplace!(intensity::AbstractVector{Float64}, method::Symbol)
Applies intensity transformation in-place. Supports :sqrt, :log1p, :log, :log2, :log10.
"""
@inline function transform_inplace!(intensity::AbstractVector{Float64}, method::Symbol)
if method === :sqrt
@inbounds @simd for i in eachindex(intensity)
intensity[i] = sqrt(max(0.0, intensity[i]))
end
elseif method === :log1p
@inbounds @simd for i in eachindex(intensity)
intensity[i] = log1p(max(0.0, intensity[i]))
end
elseif method === :log
@inbounds @simd for i in eachindex(intensity)
intensity[i] = log(max(eps(Float64), intensity[i]))
end
elseif method === :log2
@inbounds @simd for i in eachindex(intensity)
intensity[i] = log2(max(eps(Float64), intensity[i]))
end
elseif method === :log10
@inbounds @simd for i in eachindex(intensity)
intensity[i] = log10(max(eps(Float64), intensity[i]))
end
end
return intensity
end
"""
smooth_inplace!(intensity::AbstractVector{Float64}, data::MSIData;
method::Symbol=:savitzky_golay, window::Int=9, order::Int=2)
Smooths intensity in-place using a temporary buffer from the resource pool.
The SavitzkyGolay library allocates internally, but we copy the result back
to the original buffer and return the pool buffer.
"""
function smooth_inplace!(intensity::AbstractVector{Float64}, scratch::AbstractVector{Float64}, data::MSIData;
method::Symbol=:savitzky_golay, window::Int=9, order::Int=2)
n = length(intensity)
if n < 3
return intensity
end
if method === :savitzky_golay
win = isodd(window) ? window : window + 1
if n < win
return intensity
end
# SavitzkyGolay handles its own math but causes mild allocation.
res = SavitzkyGolay.savitzky_golay(collect(intensity), win, order)
@inbounds for i in eachindex(intensity)
intensity[i] = max(0.0, res.y[i])
end
elseif method === :moving_average
copyto!(scratch, intensity)
half_w = div(window, 2)
@inbounds for i in 1:n
s_idx = max(1, i - half_w)
e_idx = min(n, i + half_w)
s = 0.0
@simd for j in s_idx:e_idx
s += scratch[j]
end
intensity[i] = max(0.0, s / (e_idx - s_idx + 1))
end
end
return intensity
end
"""
baseline_subtract_inplace!(intensity::AbstractVector{Float64}, data::MSIData;
method::Symbol=:snip, iterations::Int=100, window::Int=20)
Subtracts baseline from intensity in-place. Uses two pool buffers for the
SNIP ping-pong iteration to avoid any heap allocation in the hot loop.
"""
function baseline_subtract_inplace!(intensity::AbstractVector{Float64}, scratch::AbstractVector{Float64}, data::MSIData;
method::Symbol=:snip, iterations::Int=100, window::Int=20)
n = length(intensity)
if n < 3
return intensity
end
if method === :snip
copyto!(scratch, intensity)
for k in 1:iterations
prev_val = scratch[1]
scratch[1] = min(scratch[1], scratch[2])
@inbounds for i in 2:n-1
curr_val = scratch[i]
scratch[i] = min(curr_val, 0.5 * (prev_val + scratch[i+1]))
prev_val = curr_val
end
scratch[n] = min(scratch[n], prev_val)
end
@inbounds @simd for i in 1:n
intensity[i] = max(0.0, intensity[i] - scratch[i])
end
elseif method === :convex_hull
baseline = convex_hull_baseline(intensity)
@inbounds @simd for i in eachindex(intensity)
intensity[i] = max(0.0, intensity[i] - baseline[i])
end
elseif method === :median
baseline = median_baseline(intensity; window=window)
@inbounds @simd for i in eachindex(intensity)
intensity[i] = max(0.0, intensity[i] - baseline[i])
end
end
return intensity
end
"""
detect_peaks_streaming(mz::AbstractVector, intensity::AbstractVector;
method::Symbol=:profile, snr_threshold::Float64=3.0,
half_window::Int=10, min_peak_prominence::Float64=0.1,
merge_peaks_tolerance::Float64=0.002)
Detects peaks and returns a vector of (mz, intensity) tuples.
This delegates to existing _core functions but returns a lightweight format
suitable for sparse accumulation (no NamedTuple overhead in the hot path).
"""
function detect_peaks_streaming(callback::Function, mz::AbstractVector{Float64}, intensity::AbstractVector{Float64}, scratch::AbstractVector{Float64};
method::Symbol=:profile, snr_threshold::Float64=3.0,
half_window::Int=10, min_peak_prominence::Float64=0.1,
merge_peaks_tolerance::Float64=0.002)
n = length(intensity)
if n < 3
return
end
if method === :profile || method === :wavelet
# Zero-allocation noisy estimation (using mean of bottom half)
sum_i = 0.0
@simd for i in 1:n
sum_i += intensity[i]
end
mean_i = sum_i / n
sum_noise = 0.0
count_noise = 0
@inbounds for i in 1:n
if intensity[i] < mean_i
sum_noise += intensity[i]
count_noise += 1
end
end
# Use * 1.5 as an approximation to MAD
noise_level = count_noise > 0 ? (sum_noise / count_noise) * 1.5 + eps(Float64) : mean_i + eps(Float64)
# We will use the scratch buffer to store candidate indices to avoid allocating `Int[]`
# Because scratch is Float64, we can safely store integer indices up to 2^53 exactly.
num_candidates = 0
@inbounds for i in 2:n-1
if intensity[i] > snr_threshold * noise_level
left = max(1, i - half_window)
right = min(n, i + half_window)
is_max = true
for j in left:right
if intensity[j] > intensity[i]
is_max = false
break
end
end
if is_max
# Compute prominence
min_left = intensity[i]
for j in left:i
if intensity[j] < min_left
min_left = intensity[j]
end
end
min_right = intensity[i]
for j in i:right
if intensity[j] < min_right
min_right = intensity[j]
end
end
prominence = intensity[i] - max(min_left, min_right)
if prominence > min_peak_prominence * intensity[i]
num_candidates += 1
scratch[num_candidates] = i
end
end
end
end
# Merge close peaks
if num_candidates > 0
if merge_peaks_tolerance > 0
last_idx = trunc(Int, scratch[1])
# We emit the first peak lazily down below, so let's compact them in place
num_merged = 1
for i in 2:num_candidates
idx = trunc(Int, scratch[i])
if (mz[idx] - mz[last_idx]) > merge_peaks_tolerance
num_merged += 1
scratch[num_merged] = idx
last_idx = idx
elseif intensity[idx] > intensity[last_idx]
scratch[num_merged] = idx
last_idx = idx
end
end
num_candidates = num_merged
end
# Emit merged peaks
for i in 1:num_candidates
idx = trunc(Int, scratch[i])
callback(mz[idx], intensity[idx])
end
end
elseif method === :centroid
# Just use any value over snr_threshold * mean_noise
sum_i = sum(intensity)
mean_i = sum_i / n
noise_level = mean_i + eps(Float64)
@inbounds for i in 1:n
if intensity[i] > snr_threshold * noise_level
callback(mz[i], intensity[i])
end
end
end
end
# =============================================================================
# Category B: Conditionally Streamable Kernels (Fixed-Reference)
# =============================================================================
"""
calibrate_inplace!(mz::Vector{Float64}, intensity::AbstractVector,
reference_masses::Vector{Float64}; ppm_tolerance::Float64=20.0)
Calibrates the m/z axis in-place using a fixed dictionary of internal standard
reference masses. This is streamable because the reference is constant.
Returns `true` if calibration was applied, `false` if insufficient peaks were found.
"""
function calibrate_inplace!(mz::Vector{Float64}, intensity::AbstractVector,
reference_masses::Vector{Float64}; ppm_tolerance::Float64=20.0)
matched_peaks = find_calibration_peaks_core(mz, intensity, reference_masses;
ppm_tolerance=ppm_tolerance)
if length(matched_peaks) < 2
return false # Insufficient reference peaks
end
measured = sort(collect(values(matched_peaks)))
theoretical = sort(collect(keys(matched_peaks)))
itp = linear_interpolation(measured, theoretical, extrapolation_bc=Line())
# Apply calibration in-place
@inbounds for i in eachindex(mz)
mz[i] = itp(mz[i])
end
return true
end

402
src/StreamingPipeline.jl
View File

@ -1,402 +0,0 @@
# src/StreamingPipeline.jl
# ============================================================================
# The Streaming Pipeline Executor
#
# This module provides `process_dataset!`, the Sprint 2 master function that
# streams spectral data through an in-place kernel chain and accumulates
# results into a SparseMatrixCSC without ever holding more than 1 spectrum
# per thread in RAM.
#
# Architecture:
# 1. _iterate_spectra_fast → Mmap zero-copy views
# 2. copyto!(writable_buf, view) → makes mutable copy for kernels
# 3. Kernel chain: smooth! → baseline! → peaks → bin
# 4. Thread-local (I, J, V) sparse accumulators
# 5. Final sparse(I, J, V, num_bins, num_spectra) assembly
#
# This works alongside the existing execute_full_preprocessing in
# PreprocessingPipeline.jl — it does NOT replace the app.jl integration.
# ============================================================================
using SparseArrays
using Printf
# =============================================================================
# Configuration Structs
# =============================================================================
"""
StreamingStep
Represents a single step in the streaming pipeline.
"""
struct StreamingStep
name::Symbol
params::Dict{Symbol, Any}
end
"""
PipelineConfig
Holds the complete configuration for a streaming pipeline execution.
# Fields
- `steps::Vector{StreamingStep}` — ordered sequence of processing steps
- `reference_peaks::Vector{Float64}` — fixed m/z values for calibration (Category B)
- `num_bins::Int` — number of bins for the output feature matrix
- `min_peaks_per_bin::Int` — minimum peak count to keep a bin
- `frequency_threshold::Float64` — minimum fraction of spectra a bin must appear in (0.0-1.0)
# Example
```julia
config = PipelineConfig(
steps = [
StreamingStep(:smoothing, Dict(:method => :savitzky_golay, :window => 9, :order => 2)),
StreamingStep(:baseline_correction, Dict(:method => :snip, :iterations => 100)),
StreamingStep(:normalization, Dict(:method => :tic)),
StreamingStep(:peak_picking, Dict(:method => :profile, :snr_threshold => 3.0)),
],
num_bins = 2000
)
```
"""
struct PipelineConfig
steps::Vector{StreamingStep}
reference_peaks::Vector{Float64}
num_bins::Int
min_peaks_per_bin::Int
frequency_threshold::Float64
end
# Convenience constructor with defaults
function PipelineConfig(; steps::Vector{StreamingStep}=StreamingStep[],
reference_peaks::Vector{Float64}=Float64[],
num_bins::Int=2000,
min_peaks_per_bin::Int=3,
frequency_threshold::Float64=0.0)
return PipelineConfig(steps, reference_peaks, num_bins, min_peaks_per_bin, frequency_threshold)
end
# =============================================================================
# Sparse Accumulator (Thread-Local)
# =============================================================================
"""
SparseAccumulator
Thread-local accumulator for sparse matrix construction.
Collects (row, col, val) triplets that will be assembled into
a SparseMatrixCSC at the end of the pipeline.
"""
mutable struct SparseAccumulator
I::Vector{Int} # Row indices (bin indices)
J::Vector{Int} # Column indices (spectrum indices)
V::Vector{Float64} # Values (intensities)
lck::Base.Threads.SpinLock
function SparseAccumulator(capacity_hint::Int=10000)
acc = new(
Vector{Int}(undef, 0),
Vector{Int}(undef, 0),
Vector{Float64}(undef, 0),
Base.Threads.SpinLock()
)
sizehint!(acc.I, capacity_hint)
sizehint!(acc.J, capacity_hint)
sizehint!(acc.V, capacity_hint)
return acc
end
end
"""
accumulate!(acc::SparseAccumulator, spectrum_idx::Int, bin_indices::AbstractVector{Int},
intensities::AbstractVector{Float64})
Appends peak data for one spectrum into the sparse accumulator.
"""
@inline function accumulate!(acc::SparseAccumulator, spectrum_idx::Int,
bin_indices::AbstractVector{Int},
intensities::AbstractVector{Float64})
n = length(bin_indices)
for k in 1:n
@inbounds begin
push!(acc.I, bin_indices[k])
push!(acc.J, spectrum_idx)
push!(acc.V, intensities[k])
end
end
end
# =============================================================================
# The Pipeline Executor
# =============================================================================
"""
process_dataset!(data::MSIData, config::PipelineConfig;
progress_callback::Union{Function, Nothing}=nothing,
masked_indices::Union{AbstractVector{Int}, Nothing}=nothing)
The Sprint 2 master streaming function. Processes an entire MSI dataset through
a kernel chain without holding more than 1 spectrum per thread in RAM.
# Returns
- `SparseMatrixCSC{Float64, Int}`: The feature matrix (bins × spectra)
- `Vector{Float64}`: The m/z bin centers
# Architecture
1. Ensures analytics are computed (for global m/z range)
2. Creates thread-local SparseAccumulators
3. Streams spectra via `_iterate_spectra_fast`
4. Per spectrum: copy view → kernel chain → peak detect → bin → accumulate
5. Merges accumulators → `sparse(I, J, V)`
"""
function process_dataset!(data::MSIData, config::PipelineConfig;
progress_callback::Union{Function, Nothing}=nothing,
masked_indices::Union{AbstractVector{Int}, Nothing}=nothing)
# --- Step 1: Ensure analytics are computed (provides global m/z range) ---
if !is_set(data.analytics_ready)
println("Pre-computing analytics for streaming pipeline...")
precompute_analytics(data)
end
# Determine global m/z range for binning
global_min_mz = Base.Threads.atomic_add!(data.global_min_mz, 0.0)
global_max_mz = Base.Threads.atomic_add!(data.global_max_mz, 0.0)
if !isfinite(global_min_mz) || !isfinite(global_max_mz) || global_min_mz >= global_max_mz
@warn "Invalid global m/z range: [$global_min_mz, $global_max_mz]. Cannot bin peaks."
return spzeros(0, 0), Float64[]
end
num_bins = config.num_bins
bin_edges = range(global_min_mz, stop=global_max_mz, length=num_bins + 1)
bin_centers = [(bin_edges[i] + bin_edges[i+1]) / 2 for i in 1:num_bins]
inv_bin_width = 1.0 / step(bin_edges)
num_spectra = length(data.spectra_metadata)
indices_to_process = masked_indices === nothing ? nothing : masked_indices
# --- Step 2: Create thread-local accumulators ---
n_threads = Base.Threads.nthreads()
accumulators = [SparseAccumulator(num_spectra * 10) for _ in 1:n_threads]
spectra_processed = Base.Threads.Atomic{Int}(0)
# NEW: Create dedicated workspace buffers for each thread.
# This completely eliminates the need for acquire/release and prevents deadlocks.
workspaces_mz = [Vector{Float64}(undef, 0) for _ in 1:n_threads]
workspaces_int = [Vector{Float64}(undef, 0) for _ in 1:n_threads]
workspaces_scratch = [Vector{Float64}(undef, 0) for _ in 1:n_threads]
# Pre-parse step configuration for fast dispatch in the hot loop
has_smoothing = false
has_baseline = false
has_normalization = false
has_transform = false
has_peak_picking = false
has_calibration = false
smooth_params = Dict{Symbol, Any}()
baseline_params = Dict{Symbol, Any}()
norm_params = Dict{Symbol, Any}()
transform_params = Dict{Symbol, Any}()
peak_params = Dict{Symbol, Any}()
for s in config.steps
if s.name === :smoothing
has_smoothing = true
smooth_params = s.params
elseif s.name === :baseline_correction
has_baseline = true
baseline_params = s.params
elseif s.name === :normalization
has_normalization = true
norm_params = s.params
elseif s.name === :stabilization || s.name === :intensity_transformation
has_transform = true
transform_params = s.params
elseif s.name === :peak_picking
has_peak_picking = true
peak_params = s.params
elseif s.name === :calibration
has_calibration = true
end
end
reference_masses = config.reference_peaks
# --- Step 3: Stream and process ---
start_time = time_ns()
# Use let block to capture all variables cleanly for the closure
let data=data, accumulators=accumulators, spectra_processed=spectra_processed,
bin_edges=bin_edges, num_bins=num_bins, inv_bin_width=inv_bin_width,
global_min_mz=global_min_mz,
workspaces_mz=workspaces_mz, workspaces_int=workspaces_int, workspaces_scratch=workspaces_scratch,
has_smoothing=has_smoothing, has_baseline=has_baseline,
has_normalization=has_normalization, has_transform=has_transform,
has_peak_picking=has_peak_picking, has_calibration=has_calibration,
smooth_params=smooth_params, baseline_params=baseline_params,
norm_params=norm_params, transform_params=transform_params,
peak_params=peak_params, reference_masses=reference_masses
_iterate_spectra_fast(data, indices_to_process) do idx, mz_view, int_view
thread_id = Base.Threads.threadid()
acc = accumulators[thread_id]
# --- Grab Thread-Local Workspaces ---
# No locking, no blocking, guaranteed to be available
mz_buf = workspaces_mz[thread_id]
int_buf = workspaces_int[thread_id]
scratch_buf = workspaces_scratch[thread_id]
resize!(mz_buf, length(mz_view))
resize!(int_buf, length(int_view))
resize!(scratch_buf, length(int_view))
copyto!(mz_buf, mz_view)
copyto!(int_buf, int_view)
# --- Kernel Chain (in pipeline order) ---
# Category B: Fixed-reference calibration
if has_calibration && !isempty(reference_masses)
calibrate_inplace!(mz_buf, int_buf, reference_masses)
end
# Category A: Intensity transformation
if has_transform
transform_inplace!(int_buf, get(transform_params, :method, :sqrt))
end
# Category A: Smoothing
if has_smoothing
smooth_inplace!(int_buf, scratch_buf, data;
method=get(smooth_params, :method, :savitzky_golay),
window=get(smooth_params, :window, 9),
order=get(smooth_params, :order, 2))
end
# Category A: Baseline correction
if has_baseline
baseline_subtract_inplace!(int_buf, scratch_buf, data;
method=get(baseline_params, :method, :snip),
iterations=get(baseline_params, :iterations, 100),
window=get(baseline_params, :window, 20))
end
# Category A: Normalization
if has_normalization
normalize_inplace!(int_buf, get(norm_params, :method, :tic))
end
# --- Peak Detection & Binning ---
if has_peak_picking
detect_peaks_streaming(mz_buf, int_buf, scratch_buf;
method=get(peak_params, :method, :profile),
snr_threshold=Float64(get(peak_params, :snr_threshold, 3.0)),
half_window=Int(get(peak_params, :half_window, 10)),
min_peak_prominence=Float64(get(peak_params, :min_peak_prominence, 0.1)),
merge_peaks_tolerance=Float64(get(peak_params, :merge_peaks_tolerance, 0.002))) do peak_mz, peak_int
# Bin each discovered peak directly
bin_idx = trunc(Int, (peak_mz - global_min_mz) * inv_bin_width) + 1
bin_idx = clamp(bin_idx, 1, num_bins)
push!(acc.I, bin_idx)
push!(acc.J, idx)
push!(acc.V, peak_int)
end
else
# No peak picking: bin raw intensity directly
@inbounds for i in eachindex(mz_buf)
bin_idx = trunc(Int, (mz_buf[i] - global_min_mz) * inv_bin_width) + 1
bin_idx = clamp(bin_idx, 1, num_bins)
push!(acc.I, bin_idx)
push!(acc.J, idx)
push!(acc.V, int_buf[i])
end
end
Base.Threads.atomic_add!(spectra_processed, 1)
end
end
# --- Step 4: Merge thread-local accumulators ---
total_entries = sum(length(acc.I) for acc in accumulators)
merged_I = Vector{Int}(undef, total_entries)
merged_J = Vector{Int}(undef, total_entries)
merged_V = Vector{Float64}(undef, total_entries)
offset = 0
for acc in accumulators
n = length(acc.I)
if n > 0
copyto!(merged_I, offset + 1, acc.I, 1, n)
copyto!(merged_J, offset + 1, acc.J, 1, n)
copyto!(merged_V, offset + 1, acc.V, 1, n)
offset += n
end
end
# --- Step 5: Assemble sparse matrix ---
# Use max combiner: when multiple peaks map to the same bin for same spectrum,
# keep the maximum intensity
feature_matrix = sparse(merged_I, merged_J, merged_V, num_bins, num_spectra, max)
# --- Step 6: Apply frequency threshold if configured ---
if config.frequency_threshold > 0.0
# Count how many spectra have a non-zero value in each bin
bin_presence = vec(sum(feature_matrix .> 0, dims=2))
min_count = ceil(Int, config.frequency_threshold * num_spectra)
keep_bins = findall(bin_presence .>= min_count)
feature_matrix = feature_matrix[keep_bins, :]
bin_centers = bin_centers[keep_bins]
end
duration = (time_ns() - start_time) / 1e9
n_processed = spectra_processed[]
n_nonzeros = nnz(feature_matrix)
sparsity = 1.0 - n_nonzeros / (size(feature_matrix, 1) * size(feature_matrix, 2) + 1)
@printf "Streaming pipeline complete: %d spectra processed in %.2f seconds.\n" n_processed duration
@printf "Feature matrix: %d bins × %d spectra, %d non-zeros (%.1f%% sparse)\n" size(feature_matrix, 1) size(feature_matrix, 2) n_nonzeros sparsity * 100
@printf "RAM: %.1f MB (vs %.1f MB dense)\n" (n_nonzeros * 16) / 1e6 (size(feature_matrix, 1) * size(feature_matrix, 2) * 8) / 1e6
if progress_callback !== nothing
progress_callback(1.0)
end
return feature_matrix, collect(Float64, bin_centers)
end
"""
save_sparse_matrix(matrix::SparseMatrixCSC, output_path::String)
Exports a highly optimized SparseMatrixCSC array to disk using the standard
Matrix Market Coordinate format (`.mtx`), guaranteeing no bottleneck or OOM crashes
for extremely large MS dataset persistence.
"""
function save_sparse_matrix(matrix::SparseMatrixCSC{Float64, Int}, output_path::String)
m, n = size(matrix)
nnz_val = nnz(matrix)
# Use streaming I/O with a large buffer for ultra-fast persistence
open(output_path, "w") do io
# Write Matrix Market Header
write(io, "%%MatrixMarket matrix coordinate real general\n")
write(io, "$m $n $nnz_val\n")
# Directly extract CSC properties (O(1) memory, zero allocation)
row_indices = rowvals(matrix)
values_array = nonzeros(matrix)
@inbounds for filter_j in 1:n
# nzrange returns the index bounds for non-zero elements in column 'j'
for idx in nzrange(matrix, filter_j)
i = row_indices[idx]
v = values_array[idx]
write(io, "$i $filter_j $v\n")
end
end
end
end

1671
src/imzML.jl
View File

File diff suppressed because it is too large Load Diff

381
src/mzML.jl
View File

@ -19,144 +19,6 @@ const DATA_FORMAT_ACCESSIONS = Dict{String, DataType}(
"MS:1000523" => Float64 "MS:1000523" => Float64
) )
"""
parse_instrument_metadata_mzml(stream::IO)
Parses instrument metadata from an mzML file.
# Arguments
- `stream::IO`: The stream to parse the metadata from.
# Returns
- `InstrumentMetadata`: The instrument metadata.
# Example
```julia
parse_instrument_metadata_mzml(open("test.mzML"))
```
"""
function parse_instrument_metadata_mzml(stream::IO)
println("DEBUG: Starting mzML instrument metadata parsing...")
# Initialize with default values from the InstrumentMetadata constructor
instrument_meta = InstrumentMetadata()
# Create temporary variables to hold parsed values
resolution = instrument_meta.resolution
instrument_model = instrument_meta.instrument_model
mass_accuracy_ppm = instrument_meta.mass_accuracy_ppm
polarity = instrument_meta.polarity
calibration_status = instrument_meta.calibration_status
laser_settings = Dict{String, Any}() # Use Any for heterogeneous values
vendor_preprocessing_steps = String[] # Initialize vendor_preprocessing_steps
try
# Reset stream and read a sufficiently large header block to find metadata.
seekstart(stream)
header_block = ""
header_read_limit = 50000 # Read up to 50KB of header
bytes_read = 0
while !eof(stream) && bytes_read < header_read_limit
line_pos = position(stream)
line = readline(stream)
bytes_read += (position(stream) - line_pos)
# Stop at the start of the main data section
if occursin("<run>", line) || occursin("<spectrumList>", line)
break
end
header_block *= line * "\n"
end
header_io = IOBuffer(header_block)
while !eof(header_io)
line = readline(header_io)
if occursin("<cvParam", line)
accession_match = get_attribute(line, "accession")
value_match = get_attribute(line, "value")
name_match = get_attribute(line, "name") # For laser attributes and vendor_preprocessing names
if accession_match !== nothing
acc = accession_match.captures[1]
val = (value_match !== nothing) ? value_match.captures[1] : ""
name = (name_match !== nothing) ? name_match.captures[1] : "" # NEW: Get name for logging
if acc == "MS:1000031" # instrument model
instrument_model = val
#println("DEBUG: Instrument model: $instrument_model")
elseif acc == "MS:1001496" # mass resolving power (more specific)
resolution = tryparse(Float64, val)
#println("DEBUG: Resolution (specific): $resolution")
elseif acc == "MS:1000011" && resolution === nothing # resolution (less specific)
resolution = tryparse(Float64, val)
#println("DEBUG: Resolution (less specific): $resolution")
elseif acc == "MS:1000016" # mass accuracy (ppm)
mass_accuracy_ppm = tryparse(Float64, val)
#println("DEBUG: Mass accuracy (ppm): $mass_accuracy_ppm")
elseif acc == "MS:1000130" # positive scan
polarity = :positive
#println("DEBUG: Polarity: positive")
elseif acc == "MS:1000129" # negative scan
polarity = :negative
#println("DEBUG: Polarity: negative")
elseif acc == "MS:1000592" # external calibration
calibration_status = :external
#println("DEBUG: Calibration status: external")
elseif acc == "MS:1000593" # internal calibration
calibration_status = :internal
#println("DEBUG: Calibration status: internal")
elseif acc == "MS:1000747" && calibration_status == :uncalibrated # instrument specific calibration
calibration_status = :internal # Assume as a form of internal calibration
#println("DEBUG: Calibration status: instrument specific (internal)")
elseif acc == "MS:1000867" # laser wavelength
laser_settings["wavelength_nm"] = tryparse(Float64, val)
#println("DEBUG: Laser wavelength: $(laser_settings["wavelength_nm"]) nm")
elseif acc == "MS:1000868" # laser fluence
laser_settings["fluence"] = tryparse(Float64, val)
#println("DEBUG: Laser fluence: $(laser_settings["fluence"])")
elseif acc == "MS:1000869" # laser repetition rate
laser_settings["repetition_rate_hz"] = tryparse(Float64, val)
#println("DEBUG: Laser repetition rate: $(laser_settings["repetition_rate_hz"]) Hz")
# NEW: Vendor Preprocessing terms
elseif acc == "MS:1000579" # baseline correction
push!(vendor_preprocessing_steps, "Baseline Correction")
#println("DEBUG: Vendor preprocessing step: Baseline Correction")
elseif acc == "MS:1000580" # smoothing
push!(vendor_preprocessing_steps, "Smoothing")
#println("DEBUG: Vendor preprocessing step: Smoothing")
elseif acc == "MS:1000578" # data transformation (e.g., centroiding)
push!(vendor_preprocessing_steps, "Data Transformation: $(name)")
#println("DEBUG: Vendor preprocessing step: Data Transformation: $(name)")
elseif acc == "MS:1000800" # deisotoping
push!(vendor_preprocessing_steps, "Deisotoping")
#println("DEBUG: Vendor preprocessing step: Deisotoping")
end
end
end
end
catch e
@warn "Could not fully parse instrument metadata from mzML header. Using defaults. Error: $e"
end
println("DEBUG: Finished mzML instrument metadata parsing.")
# Always return a valid object
return InstrumentMetadata(
resolution,
instrument_meta.acquisition_mode, # To be determined by load_mzml_lazy
instrument_meta.mz_axis_type,
calibration_status,
instrument_model,
mass_accuracy_ppm,
isempty(laser_settings) ? nothing : laser_settings,
polarity,
isempty(vendor_preprocessing_steps) ? nothing : vendor_preprocessing_steps # NEW field
)
end
""" """
get_spectrum_asset_metadata(stream::IO) get_spectrum_asset_metadata(stream::IO)
@ -174,15 +36,11 @@ determine the data type, compression, and axis type.
function get_spectrum_asset_metadata(stream::IO) function get_spectrum_asset_metadata(stream::IO)
start_pos = position(stream) start_pos = position(stream)
#println("DEBUG: Entering get_spectrum_asset_metadata to parse binaryDataArray...") bda_tag = find_tag(stream, r"<binaryDataArray\s+encodedLength=\"(\d+)\"" )
bda_tag = find_tag(stream, r"<binaryDataArray\s+encodedLength=\"(\d+)\"")
if bda_tag === nothing if bda_tag === nothing
throw(FileFormatError("Cannot find binaryDataArray")) error("Cannot find binaryDataArray")
end end
encoded_length = parse(Int32, bda_tag.captures[1]) encoded_length = parse(Int32, bda_tag.captures[1])
#println("DEBUG: Encoded length: $encoded_length")
# Initialize parameters as separate variables with concrete types # Initialize parameters as separate variables with concrete types
data_format::DataType = Float64 data_format::DataType = Float64
@ -204,19 +62,14 @@ function get_spectrum_asset_metadata(stream::IO)
# Use constant comparisons and dictionary lookup for better performance # Use constant comparisons and dictionary lookup for better performance
if acc_str == MZ_AXIS_ACCESSION if acc_str == MZ_AXIS_ACCESSION
axis = :mz axis = :mz
#println("DEBUG: Axis type identified as: m/z")
elseif acc_str == INTENSITY_AXIS_ACCESSION elseif acc_str == INTENSITY_AXIS_ACCESSION
axis = :intensity axis = :intensity
#println("DEBUG: Axis type identified as: intensity")
elseif haskey(DATA_FORMAT_ACCESSIONS, acc_str) elseif haskey(DATA_FORMAT_ACCESSIONS, acc_str)
data_format = DATA_FORMAT_ACCESSIONS[acc_str] data_format = DATA_FORMAT_ACCESSIONS[acc_str]
#println("DEBUG: Data format identified as: $data_format")
elseif acc_str == COMPRESSION_ACCESSION elseif acc_str == COMPRESSION_ACCESSION
compression_flag = true compression_flag = true
#println("DEBUG: Compression: true")
elseif acc_str == NO_COMPRESSION_ACCESSION elseif acc_str == NO_COMPRESSION_ACCESSION
compression_flag = false compression_flag = false
#println("DEBUG: Compression: false")
end end
end end
end end
@ -224,14 +77,12 @@ function get_spectrum_asset_metadata(stream::IO)
seek(stream, start_pos) seek(stream, start_pos)
readuntil(stream, "<binary>") readuntil(stream, "<binary>")
binary_offset = position(stream) binary_offset = position(stream)
#println("DEBUG: Binary data offset: $binary_offset")
# Move stream to the end of the binary data array for the next iteration # Move stream to the end of the binary data array for the next iteration
readuntil(stream, "</binaryDataArray>") readuntil(stream, "</binaryDataArray>")
#println("DEBUG: Exiting get_spectrum_asset_metadata.")
# Create SpectrumAsset directly from the variables # Create SpectrumAsset directly from the variables
return SpectrumAsset(data_format, compression_flag, binary_offset, encoded_length, axis, 0.0, 0.0) return SpectrumAsset(data_format, compression_flag, binary_offset, encoded_length, axis)
end end
# This function is updated to return the generic SpectrumMetadata struct # This function is updated to return the generic SpectrumMetadata struct
@ -268,16 +119,13 @@ function parse_spectrum_metadata(stream::IO, offset::Int64)
id_match = match(r"<spectrum\s+index=\"\d+\"\s+id=\"([^\"]+)", spectrum_xml) id_match = match(r"<spectrum\s+index=\"\d+\"\s+id=\"([^\"]+)", spectrum_xml)
id = id_match === nothing ? "" : id_match.captures[1] id = id_match === nothing ? "" : id_match.captures[1]
#println("DEBUG: Parsing spectrum ID: $id")
# Determine mode from the XML block # Determine mode from the XML block
mode = UNKNOWN mode = UNKNOWN
if occursin("MS:1000127", spectrum_xml) if occursin("MS:1000127", spectrum_xml)
mode = CENTROID mode = CENTROID
#println("DEBUG: Spectrum mode: CENTROID")
elseif occursin("MS:1000128", spectrum_xml) elseif occursin("MS:1000128", spectrum_xml)
mode = PROFILE mode = PROFILE
#println("DEBUG: Spectrum mode: PROFILE")
end end
# Find where the binary data list starts to parse assets # Find where the binary data list starts to parse assets
@ -296,13 +144,9 @@ function parse_spectrum_metadata(stream::IO, offset::Int64)
(asset2, asset1) (asset2, asset1)
end end
#println("DEBUG: m/z Asset - Format: $(mz_asset.format), Compressed: $(mz_asset.is_compressed), Offset: $(mz_asset.offset), Encoded Length: $(mz_asset.encoded_length)")
#println("DEBUG: Intensity Asset - Format: $(int_asset.format), Compressed: $(int_asset.is_compressed), Offset: $(int_asset.offset), Encoded Length: $(int_asset.encoded_length)")
#println("DEBUG: Finished parsing spectrum ID: $id metadata.")
# Create the new unified metadata object # Create the new unified metadata object
# For mzML, x and y coordinates are not applicable, so we use 0. # For mzML, x and y coordinates are not applicable, so we use 0.
return SpectrumMetadata(Int32(0), Int32(0), id, :sample, mode, mz_asset, int_asset) return SpectrumMetadata(Int32(0), Int32(0), id, mode, mz_asset, int_asset)
end end
""" """
@ -346,19 +190,6 @@ end
Finds the index offset in an mzML file by reading from the end. Finds the index offset in an mzML file by reading from the end.
Optimized version with better buffer management. Optimized version with better buffer management.
# Arguments
- `stream::IO`: The stream to parse the metadata from.
# Returns
- `Int64`: The index offset.
# Example
```julia
find_index_offset(open("test.mzML"))
```
""" """
function find_index_offset(stream::IO)::Int64 function find_index_offset(stream::IO)::Int64
file_size = filesize(stream) file_size = filesize(stream)
@ -371,7 +202,7 @@ function find_index_offset(stream::IO)::Int64
index_offset_match = match(r"<indexListOffset>(\d+)</indexListOffset>", footer) index_offset_match = match(r"<indexListOffset>(\d+)</indexListOffset>", footer)
if index_offset_match === nothing if index_offset_match === nothing
throw(FileFormatError("Could not find <indexListOffset>. File may not be an indexed mzML.")) error("Could not find <indexListOffset>. File may not be an indexed mzML.")
end end
return parse(Int64, index_offset_match.captures[1]) return parse(Int64, index_offset_match.captures[1])
@ -394,106 +225,170 @@ then parses the metadata for each spectrum without loading the binary data.
""" """
function load_mzml_lazy(file_path::String; cache_size::Int=100) function load_mzml_lazy(file_path::String; cache_size::Int=100)
println("DEBUG: Opening file stream for $file_path") println("DEBUG: Opening file stream for $file_path")
stream = open(file_path, "r")
# --- Handle Pool Optimization ---
# Open multiple handles to the .mzML file to avoid lock contention
num_handles = Threads.nthreads()
mzml_handles = [open(file_path, "r") for _ in 1:num_handles]
# Use the first handle for initial parsing
primary_handle = mzml_handles[1]
try try
# --- NEW: Parse instrument metadata from header ---
println("DEBUG: Parsing instrument metadata from header...")
instrument_meta = parse_instrument_metadata_mzml(primary_handle)
seekstart(primary_handle) # Reset stream after header parsing
println("DEBUG: Finding index offset...") println("DEBUG: Finding index offset...")
index_offset = find_index_offset(primary_handle) index_offset = find_index_offset(stream)
# --- NEW: Mmap Optimization with RAM Safety ---
mmap_data = nothing
try
file_size = filesize(file_path)
println("DEBUG: Memory mapping .mzML file...")
seekstart(primary_handle) # Anchor Mmap to the beginning of the file to prevent overflow
mmap_data = Mmap.mmap(primary_handle, Vector{UInt8}, (file_size,))
println("DEBUG: .mzML file mmapped successfully.")
catch e
@warn "Memory mapping failed for mzML, falling back to standard I/O: $e"
end
println("DEBUG: Seeking to index list at offset $index_offset.") println("DEBUG: Seeking to index list at offset $index_offset.")
seek(primary_handle, index_offset) seek(stream, index_offset)
println("DEBUG: Searching for '<index name=\"spectrum\">'.") println("DEBUG: Searching for '<index name=\"spectrum\">'.")
if find_tag(primary_handle, r"<index\s+name=\"spectrum\"") === nothing if find_tag(stream, r"<index\s+name=\"spectrum\"") === nothing
throw(FileFormatError("Could not find spectrum index.")) error("Could not find spectrum index.")
end end
println("DEBUG: Found spectrum index tag.")
println("DEBUG: Parsing spectrum offsets...") println("DEBUG: Parsing spectrum offsets...")
spectrum_offsets = parse_offset_list(primary_handle) spectrum_offsets = parse_offset_list(stream)
if isempty(spectrum_offsets) if isempty(spectrum_offsets)
throw(FileFormatError("No spectrum offsets found.")) error("No spectrum offsets found.")
end end
num_spectra = length(spectrum_offsets) num_spectra = length(spectrum_offsets)
println("DEBUG: Found $num_spectra spectrum offsets.") println("DEBUG: Found $num_spectra spectrum offsets.")
println("DEBUG: Parsing metadata for each spectrum...") println("DEBUG: Parsing metadata for each spectrum...")
# Pre-allocate the metadata vector for better performance
spectra_metadata = Vector{SpectrumMetadata}(undef, num_spectra) spectra_metadata = Vector{SpectrumMetadata}(undef, num_spectra)
# Use @inbounds for faster indexing in the loop
@inbounds for i in 1:num_spectra @inbounds for i in 1:num_spectra
spectra_metadata[i] = parse_spectrum_metadata(primary_handle, spectrum_offsets[i]) spectra_metadata[i] = parse_spectrum_metadata(stream, spectrum_offsets[i])
# Progress reporting for large files
if i % 1000 == 0 if i % 1000 == 0
println("DEBUG: Processed $i/$num_spectra spectra") println("DEBUG: Processed $i/$num_spectra spectra")
end end
end end
println("DEBUG: Metadata parsing complete for all $num_spectra spectra.") println("DEBUG: Metadata parsing complete.")
# Inferred global formats from first spectrum # Assuming uniform data formats, take from the first spectrum
first_meta = spectra_metadata[1] first_meta = spectra_metadata[1]
mz_format = first_meta.mz_asset.format mz_format = first_meta.mz_asset.format
intensity_format = first_meta.int_asset.format intensity_format = first_meta.int_asset.format
# Determine overall acquisition mode ... source = MzMLSource(stream, mz_format, intensity_format)
num_centroid = count(m -> m.mode == CENTROID, spectra_metadata)
num_profile = count(m -> m.mode == PROFILE, spectra_metadata)
acq_mode_symbol = if num_centroid > 0 && num_profile == 0
:centroid
elseif num_profile > 0 && num_centroid == 0
:profile
elseif num_centroid > 0 && num_profile > 0
:mixed
else
:unknown
end
final_instrument_meta = InstrumentMetadata(
instrument_meta.resolution,
acq_mode_symbol,
instrument_meta.mz_axis_type,
instrument_meta.calibration_status,
instrument_meta.instrument_model,
instrument_meta.mass_accuracy_ppm,
instrument_meta.laser_settings,
instrument_meta.polarity,
instrument_meta.vendor_preprocessing_steps
)
source = MzMLSource(mzml_handles, mz_format, intensity_format, mmap_data)
println("DEBUG: Creating MSIData object.") println("DEBUG: Creating MSIData object.")
return MSIData(source, spectra_metadata, final_instrument_meta, (0, 0), nothing, cache_size) return MSIData(source, spectra_metadata, (0, 0), nothing, cache_size)
catch e catch e
# Close all handles in the pool if initialization fails close(stream) # Ensure stream is closed on error
for h in mzml_handles
isopen(h) && close(h)
end
rethrow(e) rethrow(e)
end end
end end
"""
LoadMzml(fileName::String)
Eagerly loads all spectra from a .mzML file into memory.
This function now uses the new lazy-loading
architecture internally but presents the data in the old format.
# Arguments
- `fileName`: The path to the `.mzML` file.
# Returns
- A `2xN` matrix where `N` is the number of spectra. The first row contains
m/z arrays and the second row contains intensity arrays.
"""
function LoadMzml(fileName::String)
# Use the lazy loader to get the MSIData object
msi_data = load_mzml_lazy(fileName, cache_size=0) # No need to cache if we load all
num_spectra = length(msi_data.spectra_metadata)
# FIXED: Use concrete typed arrays instead of Array{Any}
spectra_matrix = Vector{Tuple{Vector{Float64}, Vector{Float64}}}(undef, num_spectra)
# Pre-allocate and use bounds checking optimization
@inbounds for i in 1:num_spectra
# Use the new GetSpectrum API
mz, intensity = GetSpectrum(msi_data, i)
spectra_matrix[i] = (mz, intensity)
end
# Convert to the expected 2xN format if needed by downstream code
# Note: This maintains the original interface but with better typing
result_matrix = Array{Any}(undef, (2, num_spectra))
@inbounds for i in 1:num_spectra
result_matrix[1, i] = spectra_matrix[i][1]
result_matrix[2, i] = spectra_matrix[i][2]
end
return result_matrix
end
"""
load_mzml_batch(file_path::String, spectrum_indices::AbstractVector{Int})
Loads a specific batch of spectra from an mzML file efficiently.
Useful for parallel processing or when only specific spectra are needed.
# Arguments
- `file_path`: Path to the mzML file
- `spectrum_indices`: Indices of spectra to load
# Returns
- Vector of (mz_array, intensity_array) tuples
"""
function load_mzml_batch(file_path::String, spectrum_indices::AbstractVector{Int})
msi_data = load_mzml_lazy(file_path)
num_to_load = length(spectrum_indices)
# Pre-allocate result with concrete types
results = Vector{Tuple{Vector{Float64}, Vector{Float64}}}(undef, num_to_load)
@inbounds for (i, idx) in enumerate(spectrum_indices)
mz, intensity = GetSpectrum(msi_data, idx)
results[i] = (mz, intensity)
end
return results
end
"""
get_mzml_summary(file_path::String)::NamedTuple
Quickly extracts summary information from an mzML file without loading all metadata.
# Returns
- Named tuple with: num_spectra, mz_range, intensity_range, data_formats
"""
function get_mzml_summary(file_path::String)::NamedTuple
stream = open(file_path, "r")
try
index_offset = find_index_offset(stream)
seek(stream, index_offset)
if find_tag(stream, r"<index\s+name=\"spectrum\"") === nothing
error("Could not find spectrum index.")
end
spectrum_offsets = parse_offset_list(stream)
num_spectra = length(spectrum_offsets)
# Sample first few spectra to determine formats and ranges
sample_size = min(10, num_spectra)
sample_metadata = [parse_spectrum_metadata(stream, spectrum_offsets[i]) for i in 1:sample_size]
# Extract formats from sample
mz_format = sample_metadata[1].mz_asset.format
intensity_format = sample_metadata[1].int_asset.format
return (num_spectra=num_spectra,
mz_format=mz_format,
intensity_format=intensity_format,
sample_spectra=sample_size)
finally
close(stream)
end
end
# Performance optimization: Cache frequently used regex patterns
const PRE_COMPILED_REGEX = (
encoded_length = r"<binaryDataArray\s+encodedLength=\"(\d+)\"",
spectrum_id = r"<spectrum\s+index=\"\d+\"\s+id=\"([^\"]+)",
offset = r"<offset[^>]*>(\d+)</offset>",
index_list = r"<index\s+name=\"spectrum\"",
index_offset = r"<indexListOffset>(\d+)</indexListOffset>"
)

66
start_MSI_GUI.jl
View File

@ -1,76 +1,18 @@
using Pkg using Pkg
sTime=time() sTime=time()
Pkg.activate(".") Pkg.activate(".")
ENV["GENIE_ENV"] = "dev" # Only instantiate in development mode
#ENV["GENIE_ENV"] = "prod" if get(ENV, "GENIE_ENV", "dev") != "prod"
@info "Development environment detected. Instantiating packages..."
manifest_path = joinpath(@__DIR__, "Manifest.toml")
# Selective instantiation for faster startup
if get(ENV, "GENIE_ENV", "dev") != "prod" && !isfile(manifest_path)
@info "Development environment detected and Manifest.toml missing. Instantiating packages..."
Pkg.resolve()
Pkg.instantiate() Pkg.instantiate()
Pkg.gc()
elseif get(ENV, "GENIE_ENV", "dev") != "prod"
@info "Manifest.toml found. Skipping Pkg.instantiate() for faster boot. Delete Manifest.toml if you need to re-instantiate."
ENV["GENIE_ENV"] = "prod"
end end
Pkg.gc()
using Genie using Genie
# --- Cross-Platform Startup Cleanup ---
# Remove orphaned GenieSessionFileSession directories from previous runs.
# These accumulate in the OS temp directory as jl_XXXXXX folders containing
# serialized session files (64-char hex filenames). Over long sessions or
# after crashes, they can consume gigabytes of disk space.
function cleanup_orphaned_sessions()
tmp = Base.tempdir()
cleaned_count = 0
cleaned_bytes = 0
for entry in readdir(tmp; join=false)
# Only target directories matching Julia's temp naming pattern
startswith(entry, "jl_") || continue
full_path = joinpath(tmp, entry)
isdir(full_path) || continue
# Validate: a Genie session dir contains files with 64-char hex names
try
contents = readdir(full_path)
isempty(contents) && continue
# Check if at least one file matches the 64-char hex session ID pattern
is_session_dir = any(contents) do f
length(f) == 64 && all(c -> c in "0123456789abcdef", f)
end
is_session_dir || continue
# Safe to remove — this is an orphaned Genie session directory
dir_size = sum(filesize(joinpath(full_path, f)) for f in contents; init=0)
rm(full_path; recursive=true, force=true)
cleaned_count += 1
cleaned_bytes += dir_size
catch e
@debug "Skipping $entry during cleanup: $e"
end
end
if cleaned_count > 0
size_mb = round(cleaned_bytes / (1024^2), digits=1)
@info "Startup cleanup: removed $cleaned_count orphaned session dir(s), freed $(size_mb) MB"
end
end
cleanup_orphaned_sessions()
# Load and configure Genie # Load and configure Genie
Genie.loadapp() Genie.loadapp()
# Remove html parser error discrepancy
redirect_stderr(devnull)
# Start the Genie server # Start the Genie server
@async begin @async begin
up(host="127.0.0.1", port=1481) up(host="127.0.0.1", port=1481)

936
test/benchmark.jl
View File

@ -1,936 +0,0 @@
# test/benchmark.jl
# ===================================================================
# Performance Benchmark Suite for JuliaMSI
# ===================================================================
# This script compares the performance of the new `JuliaMSI` pipeline
# against the old `julia_mzML_imzML` library.
#
# It measures the time and memory allocation for loading an `.imzML`
# file and generating an image slice multiple times.
#
# Instructions:
# 1. Fill in the placeholder paths in the "CONFIG" section below.
# 2. Run the script from the project's root directory:
# julia --project=. test/benchmark.jl
# 3. Check `test/results/` for `benchmark_results.csv` and plots.
# ===================================================================
using Pkg
Pkg.activate(joinpath(@__DIR__, "..")) # Activate main project environment
using DataFrames
using CSV
using CairoMakie
using ColorSchemes
using Statistics
using Printf
using JSON
using MSI_src # Import the new library
using .MSI_src: get_mz_slice, quantize_intensity, save_bitmap, get_multiple_mz_slices
using julia_mzML_imzML
# Remember, to run this, you must enter pkg mode in julia and add this library:
# add github.com/CINVESTAV-LABI/julia_mzML_imzML
# ===================================================================
# CONFIG: PLEASE FILL IN YOUR FILE PATHS AND PARAMETERS
# ===================================================================
struct BenchmarkCase
filepath::String
mz_value::Float64
mz_tolerance::Float64
name::String # New field for a descriptive name
end
const BENCHMARK_CASES = [
# Add your benchmark cases here. Example:
# BenchmarkCase("/path/to/your/small_file.imzML", 309.06, 0.1),
# BenchmarkCase("/path/to/your/medium_file.imzML", 896.0, 1.0),
# BenchmarkCase("/path/to/your/large_file.imzML", 100.0, 0.1),
# BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_LA-ESI/180817_NEG_Thaliana_Leaf_bottom_1_0841.imzML",116.07,0.1, "Thaliana Leaf"),
BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/Chilli/ltpmsi-chilli.imzML",420,0.1, "Chilli Pepper"), # Chilli
BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_DESI/40TopL,10TopR,30BottomL,20BottomR/40TopL,10TopR,30BottomL,20BottomR-centroid.imzML",885.5,0.1, "Colon Cancer Human"), # Human Cancer
BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML", 716.053,0.1, "Mouse Urinary Bladder"), # Mouse bladder
BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/Leafs/CE1_Leaf_R3.imzML",306.1,0.1, "Leaf"),
BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/Liv2_imzML_TIMSConvert-selected/Liv2.imzML",796.18,0.1, "Liver Cut"), #Lib2
BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML",804.3,0.1, "Mouse Stomach 2GB"), # Mouse Stomach
BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/salida_Tims/Stomach_DHB.imzML",804.3,0.1, "Mouse Stomach 4GB"), # Mouse Stomach
]
const NUM_REPETITIONS = 50 # Number of times to generate the image for averaging
const BULK_MZ_VARIATIONS = 50 # Number of m/z variations to generate for the bulk benchmark (e.g., mz, mz+0.1, ..., mz+(N-1)*0.1)
const RESULTS_DIR = joinpath(@__DIR__, "results")
# Helper function to generate m/z variations for bulk processing
function get_mz_variation_vector(base_mz::Float64, num_variations::Int)
return [base_mz + i * 0.1 for i in 0:(num_variations - 1)]
end
# Helper function to get the combined size of .imzML and .ibd files
function get_total_file_size_mb(filepath::String)
imzml_size = isfile(filepath) ? filesize(filepath) : 0
ibd_path = replace(filepath, r"\.(imzML|imzml)$" => ".ibd")
ibd_size = isfile(ibd_path) ? filesize(ibd_path) : 0
return round((imzml_size + ibd_size) / 1024^2, digits=2)
end
# ===================================================================
# BENCHMARK: JuliaMSI (New Library)
# ===================================================================
"""
benchmark_julianew(case::BenchmarkCase)
Runs the benchmark for the modern JuliaMSI library.
"""
function benchmark_julianew(case::BenchmarkCase)
println("--- Benchmarking JuliaMSI (New) on $(basename(case.filepath)) ---")
# 1. Load the data
load_stats = @timed OpenMSIData(case.filepath)
msi_data = load_stats.value
load_time = load_stats.time
load_mem = load_stats.bytes / 1024^2 # Convert to MB
image_times = []
image_mems = []
local final_image # Declare final_image outside the loop
# 2. Generate image multiple times
for i in 1:NUM_REPETITIONS
print("\rRepetition $i/$NUM_REPETITIONS...")
image_stats = @timed get_mz_slice(msi_data, case.mz_value, case.mz_tolerance)
final_image = image_stats.value
push!(image_times, image_stats.time)
push!(image_mems, image_stats.bytes / 1024^2)
end
println("\nDone.")
# 3. Save one resulting image for validation
output_path_bitmap = joinpath(RESULTS_DIR, "benchmark_JuliaMSI_$(basename(case.filepath)).bmp")
# Quantize and save the image
quantized_image, _ = quantize_intensity(final_image, 20) # 256 levels
save_bitmap(output_path_bitmap, quantized_image, MSI_src.ViridisPalette)
println("Saved validation image to $output_path_bitmap")
# 4. Calculate statistics
println("DEBUG (New): NUM_REPETITIONS = $NUM_REPETITIONS")
println("DEBUG (New): load_time = $load_time")
println("DEBUG (New): image_times = $(image_times)")
println("DEBUG (New): length(image_times) = $(length(image_times))")
avg_image_time = mean(image_times)
avg_image_mem = mean(image_mems) # Added missing calculation
println("DEBUG (New): avg_image_time = $avg_image_time")
println("DEBUG (New): sum(image_times) = $(sum(image_times))")
total_time = load_time + sum(image_times)
println("DEBUG (New): calculated total_time = $total_time")
return (
library="JuliaMSI",
filepath=basename(case.filepath),
file_size_mb=get_total_file_size_mb(case.filepath), # Use helper function
load_time_s=load_time,
avg_image_time_s=avg_image_time,
total_time_s=total_time,
load_memory_mb=load_mem,
avg_image_memory_mb=avg_image_mem,
name=case.name # Added case name to results
)
end
"""
benchmark_julianew_bulk(case::BenchmarkCase)
Runs the benchmark for the modern JuliaMSI library's bulk slice generation.
"""
function benchmark_julianew_bulk(case::BenchmarkCase)
println("--- Benchmarking JuliaMSI (Bulk) on $(basename(case.filepath)) ---")
# Generate multiple m/z values for bulk processing
mz_values_for_bulk = get_mz_variation_vector(case.mz_value, BULK_MZ_VARIATIONS)
println(" Generating $(length(mz_values_for_bulk)) m/z variations: $(mz_values_for_bulk)")
# 1. Load the data
load_stats = @timed OpenMSIData(case.filepath)
msi_data = load_stats.value
load_time = load_stats.time
load_mem = load_stats.bytes / 1024^2 # Convert to MB
image_times = []
image_mems = []
local final_slices_dict # Declare outside the loop
# 2. Generate multiple images in bulk (repeated for averaging)
for i in 1:1
print("\rRepetition $i/$NUM_REPETITIONS for bulk...")
# Assume MSI_src.get_multiple_mz_slices exists and returns Dict{Real, Matrix{Float64}}
image_stats = @timed MSI_src.get_multiple_mz_slices(msi_data, mz_values_for_bulk, case.mz_tolerance)
final_slices_dict = image_stats.value
push!(image_times, image_stats.time)
# Calculate memory for all slices combined
total_slices_mem = sum(Base.summarysize(slice) for slice in values(final_slices_dict)) / 1024^2
push!(image_mems, total_slices_mem)
end
println("\nDone.")
# No validation image saving for bulk due to multiple slices
# 3. Calculate statistics
avg_image_time = mean(image_times)
avg_image_mem = mean(image_mems)
total_time = load_time + sum(image_times) # Total time includes all repetitions
return (
library="JuliaMSIBulk",
filepath=basename(case.filepath),
file_size_mb=get_total_file_size_mb(case.filepath), # Use helper function
load_time_s=load_time,
avg_image_time_s=avg_image_time,
total_time_s=total_time,
load_memory_mb=load_mem,
avg_image_memory_mb=avg_image_mem,
name=case.name # Added case name to results
)
end
"""
benchmark_juliaold(case::BenchmarkCase)
Runs the benchmark for the old library using dynamic environment switching.
"""
function benchmark_juliaold(case::BenchmarkCase)
println("--- Benchmarking julia_mzML_imzML (Old) on $(basename(case.filepath)) ---")
try
# --- 1. Load Data ---
load_stats = @timed LoadImzml(case.filepath)
msi_data = load_stats.value
load_time = load_stats.time
load_mem = load_stats.bytes / 1024^2 # Convert to MB
image_times = []
image_mems = []
local last_image_stats_value # Declare outside the loop
# --- 2. Generate Image Repeatedly ---
for i in 1:NUM_REPETITIONS
print("\rRepetition $i/$NUM_REPETITIONS...")
image_stats = @timed GetSlice(msi_data, case.mz_value, case.mz_tolerance)
last_image_stats_value = image_stats.value # Store the value
push!(image_times, image_stats.time)
push!(image_mems, image_stats.bytes / 1024^2)
end
println("\nDone.")
# --- 3. Save Validation Image ---
final_image = last_image_stats_value
output_path_bitmap = joinpath(RESULTS_DIR, "benchmark_old_lib_$(basename(case.filepath)).bmp")
SaveBitmap(output_path_bitmap, IntQuant(final_image), ViridisPalette)
println("Saved validation image to $output_path_bitmap")
# --- 4. Calculate Averages and Totals ---
println("DEBUG (Old): NUM_REPETITIONS = $NUM_REPETITIONS")
println("DEBUG (Old): load_time = $load_time")
println("DEBUG (Old): image_times = $(image_times)")
println("DEBUG (Old): length(image_times) = $(length(image_times))")
avg_image_time = mean(image_times)
println("DEBUG (Old): avg_image_time = $avg_image_time")
println("DEBUG (Old): sum(image_times) = $(sum(image_times))")
total_time = load_time + sum(image_times)
println("DEBUG (Old): calculated total_time = $total_time")
avg_image_mem = mean(image_mems) # Added missing calculation
return (
library="julia_mzML_imzML",
filepath=basename(case.filepath),
file_size_mb=get_total_file_size_mb(case.filepath), # Use helper function
load_time_s=load_time,
avg_image_time_s=avg_image_time,
total_time_s=total_time,
load_memory_mb=load_mem,
avg_image_memory_mb=avg_image_mem,
name=case.name # Added case name to results
)
catch e
println("ERROR benchmarking old library: $e")
println("Make sure julia_mzML_imzML is installed and OLD_LIB_ENV_PATH is set correctly.")
return nothing
end
end
function benchmark_multi_file_walkthrough(cases::Vector{BenchmarkCase})
println("--- Multi-File Walkthrough Benchmark ---")
println("Simulating workflow: load file -> generate 1 image -> close -> next file")
results = []
for case in cases
if !isfile(case.filepath)
@warn "File not found: $(case.filepath). Skipping."
continue
end
println("\nProcessing: $(case.name)")
# --- NEW LIBRARY (JuliaMSI) ---
println(" [JuliaMSI]")
# Time includes metadata loading + single image generation
new_total_stats = @timed begin
data = OpenMSIData(case.filepath)
img = get_mz_slice(data, case.mz_value, case.mz_tolerance)
close(data) # Explicitly release resources
end
# --- OLD LIBRARY (julia_mzML_imzML) ---
println(" [julia_mzML_imzML]")
old_total_stats = try
@timed begin
data = LoadImzml(case.filepath) # Loads ENTIRE dataset into RAM
img = GetSlice(data, case.mz_value, case.mz_tolerance)
# Note: No explicit close function in old library
end
catch e
@warn "Old library failed on $(basename(case.filepath)): $e"
(time=Inf, bytes=Inf, value=nothing)
end
push!(results, (
name=case.name,
file_size_mb=get_total_file_size_mb(case.filepath),
julianew_time_s=new_total_stats.time,
julianew_memory_mb=new_total_stats.bytes/1024^2,
juliaold_time_s=old_total_stats.time,
juliaold_memory_mb=old_total_stats.bytes/1024^2,
speedup=old_total_stats.time / new_total_stats.time
))
# Force garbage collection between files to level the playing field
GC.gc()
if Sys.islinux()
ccall(:malloc_trim, Int32, (Int32,), 0)
end
end
return results
end
function benchmark_real_world_pipeline(cases::Vector{BenchmarkCase})
println("--- REAL-WORLD PIPELINE BENCHMARK ---")
println("Simulating: Load → Multiple Slices → Process → Save → Next File")
# Realistic m/z values (similar to your UI)
mz_values = [116.07, 420.0, 885.5, 716.053, 796.18, 804.3]
tolerance = 0.1
results = []
for case in cases
if !isfile(case.filepath)
@warn "File not found: $(case.filepath). Skipping."
continue
end
println("\nProcessing: $(case.name)")
# --- NEW LIBRARY (JuliaMSI) with bulk processing ---
println(" [JuliaMSI - Bulk Processing]")
new_stats = try
@timed begin
# 1. Load (lazy, metadata only)
data = OpenMSIData(case.filepath)
# 2. Generate ALL slices in one bulk operation
slice_dict = get_multiple_mz_slices(data, mz_values, tolerance)
# 3. Process each slice (simulate quantization, saving, etc.)
for (mass, slice) in slice_dict
# Simulate minimal processing for fair comparison
# In real pipeline: TrIQ, median filter, save bitmap
processed = quantize_intensity(slice, 256)[1]
end
# 4. Close and cleanup
close(data)
end
catch e
@warn "New library failed: $e"
(time=Inf, bytes=Inf, value=nothing)
end
# --- OLD LIBRARY (sequential processing) ---
println(" [Old Library - Sequential]")
old_stats = try
@timed begin
# 1. Load (everything into RAM)
data = LoadImzml(case.filepath)
# 2. Generate slices one by one
for mass in mz_values
slice = GetSlice(data, mass, tolerance)
# Same processing simulation
processed = IntQuant(slice)
end
# 3. Old library doesn't have explicit close
# Data stays in RAM until GC
end
catch e
@warn "Old library failed: $e"
(time=Inf, bytes=Inf, value=nothing)
end
# Memory cleanup between files (mirroring your pipeline)
GC.gc(true)
if Sys.islinux()
ccall(:malloc_trim, Int32, (Int32,), 0)
end
push!(results, (
name=case.name,
file_size_mb=get_total_file_size_mb(case.filepath),
n_slices=length(mz_values),
julianew_time_s=new_stats.time,
julianew_memory_mb=new_stats.bytes/1024^2,
juliaold_time_s=old_stats.time,
juliaold_memory_mb=old_stats.bytes/1024^2,
speedup=old_stats.time / new_stats.time,
# NEW METRICS for your pipeline:
memory_reduction_pct=100 * (old_stats.bytes - new_stats.bytes) / old_stats.bytes,
slices_per_second_new=length(mz_values) / new_stats.time,
slices_per_second_old=length(mz_values) / old_stats.time
))
end
return results
end
# ===================================================================
# Main Runner
function run_all_benchmarks(skip_benchmark::Bool=false)
results_csv_path = joinpath(RESULTS_DIR, "benchmark_results.csv")
multi_file_walkthrough_csv_path = joinpath(RESULTS_DIR, "multi_file_walkthrough.csv")
real_world_pipeline_csv_path = joinpath(RESULTS_DIR, "real_world_pipeline.csv")
local all_results::DataFrame
local walkthrough_df::DataFrame = DataFrame() # Initialize as empty
local pipeline_df::DataFrame = DataFrame() # Initialize as empty
if skip_benchmark && isfile(results_csv_path) && isfile(multi_file_walkthrough_csv_path) && isfile(real_world_pipeline_csv_path)
@info "Skipping benchmark run. Loading results from $results_csv_path"
all_results = CSV.read(results_csv_path, DataFrame)
walkthrough_df = CSV.read(multi_file_walkthrough_csv_path, DataFrame)
pipeline_df = CSV.read(real_world_pipeline_csv_path, DataFrame)
else
if isempty(BENCHMARK_CASES)
@warn "No benchmark cases found. Please fill in the `BENCHMARK_CASES` constant in `test/benchmark.jl`."
return
end
all_results = DataFrame()
for case in BENCHMARK_CASES
if !isfile(case.filepath)
@warn "File not found: $(case.filepath). Skipping."
continue
end
# Run for new library
new_results = benchmark_julianew(case)
push!(all_results, new_results, cols=:union)
GC.gc() # Trigger garbage collection
if Sys.islinux()
ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure julia returns the freed memory to OS
end
# Run for new library (Bulk)
new_bulk_results = benchmark_julianew_bulk(case)
push!(all_results, new_bulk_results, cols=:union)
GC.gc() # Trigger garbage collection
if Sys.islinux()
ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure julia returns the freed memory to OS
end
# Run for old library
old_results = benchmark_juliaold(case)
if old_results !== nothing
push!(all_results, old_results, cols=:union)
end
GC.gc() # Trigger garbage collection
if Sys.islinux()
ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure julia returns the freed memory to OS
end
end
# Save results to CSV (only if benchmark was actually run)
CSV.write(results_csv_path, all_results)
println("\nBenchmark results saved to $results_csv_path")
# --- Multi-File Walkthrough Benchmark ---
println("\n" * "="^60)
println("MULTI-FILE WALKTHROUGH BENCHMARK")
println("="^60)
walkthrough_results = benchmark_multi_file_walkthrough(BENCHMARK_CASES)
# Create a summary DataFrame
walkthrough_df = DataFrame(
name=[r.name for r in walkthrough_results],
file_size_mb=[r.file_size_mb for r in walkthrough_results],
JuliaMSI_time_s=[r.julianew_time_s for r in walkthrough_results],
Old_time_s=[r.juliaold_time_s for r in walkthrough_results],
speedup=[r.speedup for r in walkthrough_results],
JuliaMSI_mem_mb=[r.julianew_memory_mb for r in walkthrough_results],
Old_mem_mb=[r.juliaold_memory_mb for r in walkthrough_results]
)
println("\nMulti-file walkthrough results:")
show(walkthrough_df, allrows=true)
# Save to CSV
CSV.write(multi_file_walkthrough_csv_path, walkthrough_df)
# --- REAL-WORLD PIPELINE BENCHMARK ---
println("\n" * "="^60)
println("REAL-WORLD PIPELINE BENCHMARK")
println("="^60)
pipeline_results = benchmark_real_world_pipeline(BENCHMARK_CASES)
# Create a summary DataFrame
pipeline_df = DataFrame(
name=[r.name for r in pipeline_results],
file_size_mb=[r.file_size_mb for r in pipeline_results],
n_slices=[r.n_slices for r in pipeline_results],
JuliaMSI_time_s=[r.julianew_time_s for r in pipeline_results],
Old_time_s=[r.juliaold_time_s for r in pipeline_results],
speedup=[r.speedup for r in pipeline_results],
JuliaMSI_mem_mb=[r.julianew_memory_mb for r in pipeline_results],
Old_mem_mb=[r.juliaold_memory_mb for r in pipeline_results],
memory_reduction_pct=[r.memory_reduction_pct for r in pipeline_results],
slices_per_second_new=[r.slices_per_second_new for r in pipeline_results],
slices_per_second_old=[r.slices_per_second_old for r in pipeline_results]
)
println("\nReal-world pipeline results:")
show(pipeline_df, allrows=true)
# Save to CSV
CSV.write(real_world_pipeline_csv_path, pipeline_df)
end
# Generate and save plots (always happens)
generate_plots(all_results)
plot_walkthrough_results(walkthrough_df)
plot_real_world_pipeline_results(pipeline_df)
end
# ===================================================================
# PLOTTING - Professional comparative visualizations
# ===================================================================
function generate_plots(df::DataFrame)
if isempty(df)
println("Cannot generate plots from empty dataframe.")
return
end
unique_names = unique(df.name) # Use the new name field
num_files = length(unique_names)
# Calculate positions for dodged bars
positions = 1:num_files # Define positions here
# Generate A-Z labels for X-axis
alphabet_labels = [string(Char('A' + i - 1)) for i in 1:num_files]
# Use a scientific color palette
colors = Dict(
"JuliaMSI" => "#2E86AB", # Professional blue
"julia_mzML_imzML" => "#A23B72", # Professional magenta
"JuliaMSIBulk" => "#2CA02C" # Professional green for bulk
)
# Shorter labels for display
library_labels = Dict(
"JuliaMSI" => "New",
"julia_mzML_imzML" => "Old",
"JuliaMSIBulk" => "New Bulk"
)
# Define specific colors for the speedup metrics for each library
speedup_colors_new = ["#1E88E5", "#43A047", "#FFB300", "#E53935"] # Blue, Green, Yellow, Red for JuliaMSI
speedup_colors_bulk = ["#0056B3", "#008000", "#FFA500", "#CC0000"] # Slightly darker/different hues for JuliaMSIBulk
# Set modern theme
set_theme!(Theme(
fontsize = 72, # General font size (1.5x of 48)
font = "Helvetica", # Changed to Helvetica
Axis = (
backgroundcolor = :white,
spinewidth = 1.5,
xgridvisible = false,
ygridvisible = true,
ygridcolor = (:gray, 0.2),
ygridwidth = 0.5,
xticklabelrotation = π/8,
xticklabalign = (:right, :center),
yticklabalign = (:right, :center),
xlabelsize = 72, # Axis label size (1.5x of 48)
ylabelsize = 72, # Axis label size (1.5x of 48)
titlesize = 72, # Axis title size (1.5x of 48)
xticklabelsize = 60, # Tick label size (1.5x of 40)
yticklabelsize = 60, # Tick label size (1.5x of 40)
titlefont = :bold, # Ensured bold
),
Legend = (
backgroundcolor = (:white, 0.95),
framecolor = :gray,
framewidth = 1,
padding = (10, 10, 5, 5),
labelsize = 72, # Legend label size (1.5x of 48)
titlesize = 72, # Legend title size (1.5x of 48)
)
))
# ===================================================================
# Plot 1: Main Performance Dashboard
# ===================================================================
fig1 = Figure(size = (3600, 2400)) # Explicitly set figure size
# Define grid layout with specific areas
g1 = fig1[2, 1:2] = GridLayout()
# Total Time
ax1 = Axis(g1[1, 1],
title = "Total Time",
ylabel = "Seconds", # Changed to linear scale, removed "(log scale)"
xticks = (positions, alphabet_labels) # Use alphabetical labels for x-axis
)
# Average Image Time
ax2 = Axis(g1[1, 2],
title = "Image Generation Time",
ylabel = "Seconds",
yticks = LinearTicks(5),
xticks = (positions, alphabet_labels) # Use alphabetical labels for x-axis
)
# Loading Memory
ax3 = Axis(g1[2, 1],
title = "Loading Memory",
ylabel = "MB",
yticks = LinearTicks(5),
xticks = (positions, alphabet_labels) # Use alphabetical labels for x-axis
)
# Image Generation Memory
ax4 = Axis(g1[2, 2],
title = "Image Generation Memory",
ylabel = "MB",
yticks = LinearTicks(5),
xticks = (positions, alphabet_labels) # Use alphabetical labels for x-axis
)
bar_width = 0.25 # Reduced bar width
x_offsets = [-0.3, 0, 0.3] # Offsets for 3 bars to avoid overlap
# Arrays to store speedup factors
total_time_speedup_new = zeros(num_files)
image_time_speedup_new = zeros(num_files)
load_mem_speedup_new = zeros(num_files)
img_mem_speedup_new = zeros(num_files)
total_time_speedup_bulk = zeros(num_files)
image_time_speedup_bulk = zeros(num_files)
load_mem_speedup_bulk = zeros(num_files)
img_mem_speedup_bulk = zeros(num_files)
# Initialize a list to hold rows for the summary DataFrame
summary_df_rows = []
for (i, name_for_df_filter) in enumerate(unique_names) # Iterate through unique names
# Get data for this file
new_data = df[(df.name .== name_for_df_filter) .& (df.library .== "JuliaMSI"), :]
bulk_data = df[(df.name .== name_for_df_filter) .& (df.library .== "JuliaMSIBulk"), :]
old_data = df[(df.name .== name_for_df_filter) .& (df.library .== "julia_mzML_imzML"), :]
# Plot total time
barplot!(ax1, [i + x_offsets[1], i + x_offsets[2], i + x_offsets[3]],
[new_data.total_time_s[1], bulk_data.total_time_s[1], old_data.total_time_s[1]],
color = [colors["JuliaMSI"], colors["JuliaMSIBulk"], colors["julia_mzML_imzML"]],
width = bar_width
)
# Plot average image time
barplot!(ax2, [i + x_offsets[1], i + x_offsets[2], i + x_offsets[3]],
[new_data.avg_image_time_s[1], bulk_data.avg_image_time_s[1], old_data.avg_image_time_s[1]],
color = [colors["JuliaMSI"], colors["JuliaMSIBulk"], colors["julia_mzML_imzML"]],
width = bar_width
)
# Plot loading memory
barplot!(ax3, [i + x_offsets[1], i + x_offsets[2], i + x_offsets[3]],
[new_data.load_memory_mb[1], bulk_data.load_memory_mb[1], old_data.load_memory_mb[1]],
color = [colors["JuliaMSI"], colors["JuliaMSIBulk"], colors["julia_mzML_imzML"]],
width = bar_width
)
# Plot image generation memory
barplot!(ax4, [i + x_offsets[1], i + x_offsets[2], i + x_offsets[3]],
[new_data.avg_image_memory_mb[1], bulk_data.avg_image_memory_mb[1], old_data.avg_image_memory_mb[1]],
color = [colors["JuliaMSI"], colors["JuliaMSIBulk"], colors["julia_mzML_imzML"]],
width = bar_width
)
# Calculate speedup factors (higher is better for new library)
if nrow(new_data) == 1 && nrow(old_data) == 1
total_time_speedup_new[i] = old_data.total_time_s[1] / new_data.total_time_s[1]
image_time_speedup_new[i] = old_data.avg_image_time_s[1] / new_data.avg_image_time_s[1]
load_mem_speedup_new[i] = old_data.load_memory_mb[1] / new_data.load_memory_mb[1]
img_mem_speedup_new[i] = old_data.avg_image_memory_mb[1] / new_data.avg_image_memory_mb[1]
end
if nrow(bulk_data) == 1 && nrow(old_data) == 1
total_time_speedup_bulk[i] = old_data.total_time_s[1] / bulk_data.total_time_s[1]
image_time_speedup_bulk[i] = old_data.avg_image_time_s[1] / bulk_data.avg_image_time_s[1]
load_mem_speedup_bulk[i] = old_data.load_memory_mb[1] / bulk_data.load_memory_mb[1]
img_mem_speedup_bulk[i] = old_data.avg_image_memory_mb[1] / bulk_data.avg_image_memory_mb[1]
end
# Collect data for the summary table CSV
current_case_name = unique_names[i]
case_row = df[(df.name .== current_case_name) .& (df.library .== "JuliaMSI"), :] # Any library row will have file_size_mb
file_size = nrow(case_row) == 1 ? case_row.file_size_mb[1] : NaN
push!(summary_df_rows, Dict(
"DS" => unique_names[i],
"File Size (MB)" => file_size,
"JMSI-TT Speedup" => round(total_time_speedup_new[i], digits=2),
"JMSI-IT Speedup" => round(image_time_speedup_new[i], digits=2),
"JMSI-LM Speedup" => round(load_mem_speedup_new[i], digits=2),
"JMSI-IM Speedup" => round(img_mem_speedup_new[i], digits=2),
"JBulk-TT Speedup" => round(total_time_speedup_bulk[i], digits=2),
"JBulk-IT Speedup" => round(image_time_speedup_bulk[i], digits=2),
"JBulk-LM Speedup" => round(load_mem_speedup_bulk[i], digits=2),
"JBulk-IM Speedup" => round(img_mem_speedup_bulk[i], digits=2)
))
end
# Create DataFrame from collected rows
summary_df = DataFrame(summary_df_rows)
# Save the summary table to CSV
CSV.write(joinpath(RESULTS_DIR, "benchmark_summary_table.csv"), summary_df)
println("\nBenchmark summary table saved to $(joinpath(RESULTS_DIR, "benchmark_summary_table.csv"))")
# Add overall title
Label(fig1[1, :], "MSI Library Performance Benchmark", fontsize = 72, font = :bold)
# Add a legend
Legend(fig1[3, :], # Moved to row 3, spanning all columns
[PolyElement(color = colors["JuliaMSI"]), PolyElement(color = colors["JuliaMSIBulk"]), PolyElement(color = colors["julia_mzML_imzML"])],
[library_labels["JuliaMSI"], library_labels["JuliaMSIBulk"], library_labels["julia_mzML_imzML"]],
"Library",
orientation = :horizontal,
tellwidth = false, tellheight = true,
halign = :center, valign = :top, # Centered horizontally, aligned to top
margin = (10, 10, 10, 10)
)
# Adjust layout spacing
colgap!(g1, 1, 30)
rowgap!(g1, 1, 20)
# Ensure all rows and columns in g1 expand to fill available space
rowsize!(g1, 1, Auto())
rowsize!(g1, 2, Auto())
colsize!(g1, 1, Auto())
colsize!(g1, 2, Auto())
# Ensure the row containing the legend expands in the top-level figure layout
rowsize!(fig1.layout, 3, Auto())
# Save plots
save(joinpath(RESULTS_DIR, "performance_dashboard.png"), fig1, px_per_unit = 1)
println("\nEnhanced visualizations saved to $RESULTS_DIR")
end
function plot_walkthrough_results(walkthrough_df)
fig = Figure(size=(3600, 2400))
# Plot 1: Total time per file
ax1 = Axis(fig[1, 1],
title="Time per File (Load + Generate 1 Image)",
ylabel="Time (seconds)",
xticks=(collect(1:nrow(walkthrough_df)), walkthrough_df.name),
xticklabelrotation=π/8)
# Grouped bars
barplot!(ax1, collect(1:nrow(walkthrough_df)) .- 0.2, walkthrough_df.JuliaMSI_time_s,
width=0.4, color="#2E86AB", label="JuliaMSI (New)")
barplot!(ax1, collect(1:nrow(walkthrough_df)) .+ 0.2, walkthrough_df.Old_time_s,
width=0.4, color="#A23B72", label="julia_mzML_imzML (Old)")
# Plot 2: Speedup factor
ax2 = Axis(fig[1, 2],
title="Speedup Factor (Old/New Time)",
ylabel="Speedup (Higher = Better)",
xticks=(collect(1:nrow(walkthrough_df)), walkthrough_df.name),
xticklabelrotation=π/8)
# Bars above 1 = new library is faster
colors = [x >= 1 ? "#2CA02C" : "#E53935" for x in walkthrough_df.speedup]
barplot!(ax2, collect(1:nrow(walkthrough_df)), walkthrough_df.speedup, color=colors)
hlines!(ax2, [1.0], color=:black, linestyle=:dash, linewidth=2)
# Plot 3: Memory comparison
ax3 = Axis(fig[2, 1],
title="Memory Usage per File",
ylabel="Memory (MB)",
xticks=(collect(1:nrow(walkthrough_df)), walkthrough_df.name),
xticklabelrotation=π/8)
barplot!(ax3, collect(1:nrow(walkthrough_df)) .- 0.2, walkthrough_df.JuliaMSI_mem_mb,
width=0.4, color="#2E86AB", label="JuliaMSI")
barplot!(ax3, collect(1:nrow(walkthrough_df)) .+ 0.2, walkthrough_df.Old_mem_mb,
width=0.4, color="#A23B72", label="Old Library")
# Plot 4: Time vs File Size (scatter)
ax4 = Axis(fig[2, 2],
title="Processing Time vs File Size",
xlabel="File Size (MB)",
ylabel="Time (seconds)")
scatter!(ax4, walkthrough_df.file_size_mb, walkthrough_df.JuliaMSI_time_s,
color="#2E86AB", markersize=20, label="JuliaMSI")
scatter!(ax4, walkthrough_df.file_size_mb, walkthrough_df.Old_time_s,
color="#A23B72", markersize=20, label="Old Library")
# Add linear trend lines
if nrow(walkthrough_df) > 1
linreg_new = hcat(fill(1.0, nrow(walkthrough_df)), walkthrough_df.file_size_mb) \ walkthrough_df.JuliaMSI_time_s
linreg_old = hcat(fill(1.0, nrow(walkthrough_df)), walkthrough_df.file_size_mb) \ walkthrough_df.Old_time_s
x_vals = [minimum(walkthrough_df.file_size_mb), maximum(walkthrough_df.file_size_mb)]
lines!(ax4, x_vals, linreg_new[1] .+ linreg_new[2] .* x_vals, color="#2E86AB", linewidth=2)
lines!(ax4, x_vals, linreg_old[1] .+ linreg_old[2] .* x_vals, color="#A23B72", linewidth=2)
end
# Add legends
Legend(fig[0, :], [PolyElement(color="#2E86AB"), PolyElement(color="#A23B72")],
["JuliaMSI (New)", "julia_mzML_imzML (Old)"],
orientation=:horizontal, tellwidth=false, framevisible=true)
# Adjust layout
rowgap!(fig.layout, 1)
colgap!(fig.layout, 1)
save(joinpath(RESULTS_DIR, "multi_file_walkthrough.png"), fig, px_per_unit=1)
println("\nMulti-file walkthrough plot saved.")
end
function plot_real_world_pipeline_results(pipeline_df::DataFrame)
fig = Figure(size=(3600, 2400)) # Larger figure for more detailed plots
# Filter out rows with Inf or NaN values in critical plotting columns
# Create a copy to avoid modifying the original DataFrame
filtered_df = deepcopy(pipeline_df)
# Filter for finite values in key metrics for plotting
filter!(row -> isfinite(row.JuliaMSI_time_s) && isfinite(row.Old_time_s) &&
isfinite(row.speedup) && isfinite(row.JuliaMSI_mem_mb) &&
isfinite(row.Old_mem_mb) && isfinite(row.memory_reduction_pct) &&
isfinite(row.slices_per_second_new) && isfinite(row.slices_per_second_old),
filtered_df)
if isempty(filtered_df)
println("No finite data points for real-world pipeline plotting after filtering. Skipping plot generation.")
save(joinpath(RESULTS_DIR, "real_world_pipeline.png"), fig, px_per_unit=1)
return
end
num_rows_filtered = nrow(filtered_df)
x_positions = collect(1:num_rows_filtered)
# Plot 1: Total Time Comparison
ax1 = Axis(fig[1, 1],
title="Real-World Pipeline: Total Time",
ylabel="Time (seconds)",
xticks=(x_positions, filtered_df.name),
xticklabelrotation=π/8)
barplot!(ax1, x_positions .- 0.2, filtered_df.JuliaMSI_time_s,
width=0.4, color="#2E86AB", label="JuliaMSI (New)")
barplot!(ax1, x_positions .+ 0.2, filtered_df.Old_time_s,
width=0.4, color="#A23B72", label="Old Library")
# Plot 2: Speedup Factor
ax2 = Axis(fig[1, 2],
title="Real-World Pipeline: Speedup Factor (Old/New Time)",
ylabel="Speedup (Higher = Better)",
xticks=(x_positions, filtered_df.name),
xticklabelrotation=π/8)
colors_speedup = [x >= 1 ? "#2CA02C" : "#E53935" for x in filtered_df.speedup]
barplot!(ax2, x_positions, filtered_df.speedup, color=colors_speedup)
hlines!(ax2, [1.0], color=:black, linestyle=:dash, linewidth=2)
# Plot 3: Memory Usage Comparison
ax3 = Axis(fig[2, 1],
title="Real-World Pipeline: Memory Usage",
ylabel="Memory (MB)",
xticks=(x_positions, filtered_df.name),
xticklabelrotation=π/8)
barplot!(ax3, x_positions .- 0.2, filtered_df.JuliaMSI_mem_mb,
width=0.4, color="#2E86AB", label="JuliaMSI")
barplot!(ax3, x_positions .+ 0.2, filtered_df.Old_mem_mb,
width=0.4, color="#A23B72", label="Old Library")
# Plot 4: Memory Reduction Percentage
ax4 = Axis(fig[2, 2],
title="Real-World Pipeline: Memory Reduction vs Old (%)",
ylabel="Memory Reduction (%)",
xticks=(x_positions, filtered_df.name),
xticklabelrotation=π/8)
colors_mem_red = [x > 0 ? "#2CA02C" : "#E53935" for x in filtered_df.memory_reduction_pct]
barplot!(ax4, x_positions, filtered_df.memory_reduction_pct, color=colors_mem_red)
hlines!(ax4, [0.0], color=:black, linestyle=:dash, linewidth=2)
# Plot 5: Slices Per Second
ax5 = Axis(fig[3, 1],
title="Real-World Pipeline: Slices per Second",
ylabel="Slices/sec (Higher = Better)",
xticks=(x_positions, filtered_df.name),
xticklabelrotation=π/8)
barplot!(ax5, x_positions .- 0.2, filtered_df.slices_per_second_new,
width=0.4, color="#2E86AB", label="JuliaMSI")
barplot!(ax5, x_positions .+ 0.2, filtered_df.slices_per_second_old,
width=0.4, color="#A23B72", label="Old Library")
# Overall Legend
Legend(fig[0, :], [PolyElement(color="#2E86AB"), PolyElement(color="#A23B72")],
["JuliaMSI (New)", "julia_mzML_imzML (Old)"],
orientation=:horizontal, tellwidth=false, framevisible=true)
# Adjust layout
rowgap!(fig.layout, 1)
colgap!(fig.layout, 1)
save(joinpath(RESULTS_DIR, "real_world_pipeline.png"), fig, px_per_unit=1)
println("\nReal-world pipeline plot saved.")
end
# --- Execute ---
mkpath(RESULTS_DIR)
run_all_benchmarks(true) # true: skip processing files # false: create new matrix

215
test/benchmark_v3.jl
View File

@ -1,215 +0,0 @@
# test/benchmark_v3.jl
# ===================================================================
# High-Precision Performance Benchmark Suite (v3)
# ===================================================================
# This script integrates the robust statistical sampling of `BenchmarkTools`
# with the visual comparative mechanics against the legacy library.
# It captures the paradigm shift from "Time per slice" to "Pipeline Throughput".
# ===================================================================
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
using BenchmarkTools
using DataFrames
using CSV
using CairoMakie
using Statistics
using MSI_src
using julia_mzML_imzML
struct BenchmarkCase
filepath::String
mz_value::Float64
mz_tolerance::Float64
name::String
end
const RESULTS_DIR = joinpath(@__DIR__, "results")
# List to benchmark
const BENCHMARK_CASES = [
BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/Chilli/ltpmsi-chilli.imzML", 420.0, 0.1, "Chilli Pepper"),
BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_DESI/40TopL,10TopR,30BottomL,20BottomR/40TopL,10TopR,30BottomL,20BottomR-centroid.imzML", 885.5, 0.1, "Colon Cancer Human"),
BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML", 716.053, 0.1, "Mouse Urinary Bladder"),
BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/Leafs/CE1_Leaf_R3.imzML", 306.1, 0.1, "Leaf"),
BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/Liv2_imzML_TIMSConvert-selected/Liv2.imzML", 796.18, 0.1, "Liver Cut"),
BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/salida_Tims/Stomach_DHB.imzML", 804.3, 0.1, "Mouse Stomach 4GB")
]
function get_total_file_size_mb(filepath::String)
imzml_size = isfile(filepath) ? filesize(filepath) : 0
ibd_path = replace(filepath, r"\.(imzML|imzml)$" => ".ibd")
ibd_size = isfile(ibd_path) ? filesize(ibd_path) : 0
return round((imzml_size + ibd_size) / 1024^2, digits=2)
end
function flush_memory()
GC.gc(true)
if Sys.islinux()
# Force the OS to reclaim memory from the glibc allocator
ccall(:malloc_trim, Int32, (Int32,), 0)
end
end
function run_v3_benchmarks()
mkpath(RESULTS_DIR)
results = DataFrame()
println("="^60)
println("STARTING ENTERPRISE BENCHMARK SUITE (v3)")
println("="^60)
for case in BENCHMARK_CASES
if !isfile(case.filepath)
@warn "File not found: $(case.filepath). Skipping."
continue
end
println("\n--- Target: $(case.name) ---")
file_size = get_total_file_size_mb(case.filepath)
# ------------------------------------------------------------
# 1. JuliaMSI (New Architecture)
# ------------------------------------------------------------
println("[JuliaMSI - New Engine]")
flush_memory()
# Load Phase (Metadata Only / Memory Mapping)
load_stats_new = @timed OpenMSIData(case.filepath)
msi_data = load_stats_new.value
load_time_new_s = load_stats_new.time
mem_load_new_mb = load_stats_new.bytes / 1024^2
# Slicing Phase (High Precision)
b_slice_new = @benchmark get_mz_slice($msi_data, $(case.mz_value), $(case.mz_tolerance)) samples=10 seconds=5
mean_time_new_s = mean(b_slice_new.times) / 1e9
# Calculate Amortized Throughput (10 slices)
# Includes the 'Loading Wall' penalty
total_time_10_new = load_time_new_s + (10 * mean_time_new_s)
amortized_sps_new = 10.0 / total_time_10_new
close(msi_data)
flush_memory()
# ------------------------------------------------------------
# 2. julia_mzML_imzML (Legacy Architecture)
# ------------------------------------------------------------
println("[julia_mzML_imzML - Legacy]")
load_stats_old = try
@timed LoadImzml(case.filepath)
catch e
@warn "Legacy load failed: $e"
(time=Inf, bytes=Inf, value=nothing)
end
old_data = load_stats_old.value
mem_load_old_mb = load_stats_old.bytes / 1024^2
load_time_old_s = Inf
mean_time_old_s = Inf
amortized_sps_old = 0.0
if old_data !== nothing
b_slice_old = try
# Legacy can be extremely slow, constrain it heavily
@benchmark GetSlice($old_data, $(case.mz_value), $(case.mz_tolerance)) samples=3 seconds=10
catch e
@warn "Legacy slice failed: $e"
nothing
end
if b_slice_old !== nothing
mean_time_old_s = mean(b_slice_old.times) / 1e9
# Accurately reflect the massive load time in the throughput
load_time_old_s = load_stats_old.time
total_time_10_old = load_time_old_s + (10 * mean_time_old_s)
amortized_sps_old = 10.0 / total_time_10_old
end
end
flush_memory()
# ------------------------------------------------------------
# Record
# ------------------------------------------------------------
push!(results, (
Dataset = case.name,
FileSize_MB = file_size,
# Legacy Metrics
Legacy_LoadMem_MB = mem_load_old_mb,
Legacy_LoadTime_s = load_stats_old.time,
Legacy_Amortized_SPS_10 = amortized_sps_old,
# New Metrics
New_LoadMem_MB = mem_load_new_mb,
New_LoadTime_s = load_time_new_s,
New_Amortized_SPS_10 = amortized_sps_new,
# Competitive Deltas
RAM_Reduction_Pct = isfinite(mem_load_old_mb) ? ((mem_load_old_mb - mem_load_new_mb) / mem_load_old_mb)*100 : NaN,
UX_Speedup_Factor = isfinite(amortized_sps_old) ? (amortized_sps_new / amortized_sps_old) : NaN
))
println(" > RAM Reduction: $(round(results[end, :RAM_Reduction_Pct], digits=2))%")
println(" > Amortized Throughput (10 Slices): $(round(amortized_sps_old, digits=2)) -> $(round(amortized_sps_new, digits=2)) slices/sec")
end
csv_path = joinpath(RESULTS_DIR, "v3_enterprise_benchmarks.csv")
CSV.write(csv_path, results)
println("\nData saved to $csv_path")
plot_v3_results(results)
return results
end
function plot_v3_results(df::DataFrame)
if isempty(df) return end
fig = Figure(size=(1800, 1200), fontsize=24)
x_pos = 1:nrow(df)
labels = df.Dataset
# ---------- Plot 1: The RAM Revolution ----------
ax1 = Axis(fig[1, 1],
title="Initial RAM Cost (Loading & Mmap)",
ylabel="Memory (MB)",
xticks=(x_pos, labels), xticklabelrotation=π/8)
barplot!(ax1, x_pos .- 0.2, df.New_LoadMem_MB, color="#10b981", width=0.4, label="JuliaMSI (Mmap Lazy)")
barplot!(ax1, x_pos .+ 0.2, df.Legacy_LoadMem_MB, color="#ef4444", width=0.4, label="Legacy (Dense Matrix)")
axislegend(ax1, position=:lt)
# ---------- Plot 2: Throughput Leap ----------
ax2 = Axis(fig[1, 2],
title="Amortized Throughput (10 Slices - Higher is Better)",
ylabel="Slices / Second (Inc. Load Time)",
xticks=(x_pos, labels), xticklabelrotation=π/8)
barplot!(ax2, x_pos .- 0.2, df.New_Amortized_SPS_10, color="#10b981", width=0.4, label="JuliaMSI")
barplot!(ax2, x_pos .+ 0.2, df.Legacy_Amortized_SPS_10, color="#ef4444", width=0.4, label="Legacy")
axislegend(ax2, position=:lt)
# ---------- Plot 3: The Scaling Wall (File Size vs Throughput) ----------
ax3 = Axis(fig[2, 1:2],
title="Performance Scaling: UX Throughput vs Dataset Size",
xlabel="Dataset File Size (MB)",
ylabel="Amortized Throughput (Slices/sec)")
scatterlines!(ax3, df.FileSize_MB, df.New_Amortized_SPS_10, color="#10b981", markersize=15, linewidth=4, label="JuliaMSI Engine")
scatterlines!(ax3, df.FileSize_MB, df.Legacy_Amortized_SPS_10, color="#ef4444", markersize=15, linewidth=4, label="Legacy Engine")
axislegend(ax3, position=:rt)
save(joinpath(RESULTS_DIR, "v3_enterprise_dashboard.png"), fig)
println("\nDashboards saved successfully.")
end
# Ensure we process if run as the main script
if abspath(PROGRAM_FILE) == @__FILE__
run_v3_benchmarks()
end

48
test/new_benchmark_mmap.jl
View File

@ -1,48 +0,0 @@
using BenchmarkTools
using MSI_src
using Statistics
using DataFrames
# ===================================================================
# HIGH-PRECISION COMPARATIVE SUITE
# ===================================================================
function run_advanced_benchmark(path, mz, tol)
println("\n" * "="^40)
println("TARGET: $(basename(path))")
println("="^40)
# 1. NEW LIBRARY: Metadata Load (The "Control Tower" startup)
# This measures how fast the Mmap and Cache system works
t_load_new = @belapsed OpenMSIData($path)
# 2. NEW LIBRARY: Slice Generation (The "Streaming" speed)
msi_new = OpenMSIData(path)
# We use @benchmark to get a distribution (min, mean, max)
b_slice_new = @benchmark get_mz_slice($msi_new, $mz, $tol)
# --- Metrics Table ---
results = DataFrame(
Metric = ["Metadata Load", "Slice Gen (Min)", "Slice Gen (Mean)", "Allocations"],
JuliaMSI = [
"$(round(t_load_new * 1000, digits=2)) ms",
"$(round(minimum(b_slice_new.times)/1e6, digits=2)) ms",
"$(round(mean(b_slice_new.times)/1e6, digits=2)) ms",
"$(b_slice_new.allocs) allocs"
]
)
println(results)
# --- The "Throughput" Test ---
# How many slices per second can we handle?
throughput_new = 1.0 / mean(b_slice_new.times/1e9)
println("\nThroughput: $(round(throughput_new, digits=1)) slices/sec")
return results
end
# Example Run
@time run_advanced_benchmark("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML", 716.053, 0.1)
# For multithread: julia --threads auto --project=. test/new_benchmark_mmap.jl

4
test/readme.md
View File

@ -40,10 +40,10 @@ Before running the tests, you must edit the `test/run_tests.jl` file to point to
2. **Execute the Different Test Scripts**: 2. **Execute the Different Test Scripts**:
Run the following command from the project's root directory. This will install the necessary dependencies and run the tests. Run the following command from the project's root directory. This will install the necessary dependencies and run the tests.
```bash ```bash
julia --threads auto --project=. test/run_tests.jl julia --project=. test/run_tests.jl
``` ```
```bash ```bash
julia --threads auto --project=. test/run_preprocessing.jl julia --project=. test/run_preprocessing.jl
``` ```
3. **Check the Results**: 3. **Check the Results**:

323
test/run_descriptive_preprocessing.jl
View File

@ -1,323 +0,0 @@
using CairoMakie
using Statistics
using Colors
"""
Plots the effect of baseline correction, showing how different numbers of SNIP
iterations affect the estimated baseline.
"""
function plot_baseline_correction_details()
fig = Figure(size=(1200, 700), fontsize=14)
ax = Axis(fig[1, 1], title="Detailed View: Baseline Correction (SNIP)", xlabel="m/z", ylabel="Intensity")
mz_range = range(200, 1000, length=1500)
true_signal = 1.2 .* exp.(-0.001 .* (mz_range .- 450).^2) .+ 0.8 .* exp.(-0.002 .* (mz_range .- 750).^2)
baseline_comp = 0.15 .+ 0.1 .* sin.(mz_range ./ 60) .+ 0.0001 .* mz_range
noise = 0.04 .* randn(length(mz_range))
raw_signal = true_signal .+ baseline_comp .+ noise
# Simplified SNIP simulation for visualization
function simple_snip(y, iterations)
b = copy(y)
for _ in 1:iterations
for i in 2:length(b)-1
b[i] = min(b[i], 0.5 * (b[i-1] + b[i+1]))
end
end
return b
end
baseline_iter_20 = simple_snip(raw_signal, 20)
baseline_iter_200 = simple_snip(raw_signal, 200)
corrected_signal = raw_signal .- baseline_iter_200
lines!(ax, mz_range, raw_signal, color=(:grey, 0.6), label="Raw Signal")
lines!(ax, mz_range, corrected_signal, color=:green, linewidth=2.5, label="Corrected Signal")
l1 = lines!(ax, mz_range, baseline_iter_20, color=(:red, 0.7), linestyle=:dash, linewidth=2, label="Baseline (iterations: 20)")
l2 = lines!(ax, mz_range, baseline_iter_200, color=:red, linewidth=2.5, label="Baseline (iterations: 200)")
text!(ax, "Parameter: `iterations`\nMore iterations result in a more aggressive baseline that follows the signal floor more closely.",
position=Point2f(0.05, 0.95), space=:relative, align=(:left, :top), fontsize=12)
axislegend(ax, position=:rt)
save("descriptive_plot_baseline_correction.png", fig)
println("Saved: descriptive_plot_baseline_correction.png")
end
"""
Plots the effect of smoothing, comparing different window sizes.
"""
function plot_smoothing_details()
fig = Figure(size=(1200, 700), fontsize=14)
ax = Axis(fig[1, 1], title="Detailed View: Smoothing (Savitzky-Golay)", xlabel="m/z", ylabel="Intensity")
mz_range = range(400, 700, length=1000)
true_signal = 0.8 .* exp.(-0.002 .* (mz_range .- 500).^2) .+ 0.6 .* exp.(-0.001 .* (mz_range .- 600).^2)
noisy_signal = true_signal .+ 0.1 .* randn(length(mz_range))
function simple_moving_average(y, window)
smoothed = similar(y)
for i in 1:length(y)
start_idx = max(1, i - window ÷ 2)
end_idx = min(length(y), i + window ÷ 2)
smoothed[i] = mean(y[start_idx:end_idx])
end
return smoothed
end
smoothed_small_window = simple_moving_average(noisy_signal, 5)
smoothed_large_window = simple_moving_average(noisy_signal, 21)
lines!(ax, mz_range, noisy_signal, color=(:red, 0.4), label="Noisy Signal")
lines!(ax, mz_range, true_signal, color=:black, linestyle=:dash, linewidth=2, label="True Signal")
lines!(ax, mz_range, smoothed_small_window, color=:blue, linewidth=2, label="Smoothed (window: 5)")
lines!(ax, mz_range, smoothed_large_window, color=:purple, linewidth=2, label="Smoothed (window: 21)")
text!(ax, "Parameter: `window`\nA larger window increases smoothing but may broaden peaks.",
position=Point2f(0.05, 0.95), space=:relative, align=(:left, :top), fontsize=12)
axislegend(ax, position=:rt)
save("descriptive_plot_smoothing.png", fig)
println("Saved: descriptive_plot_smoothing.png")
end
"""
Visualizes the peak picking process for profile-mode data, illustrating the
effects of SNR threshold and peak prominence.
"""
function plot_peak_picking_profile_details()
fig = Figure(size=(1200, 700), fontsize=14)
ax = Axis(fig[1, 1], title="Detailed View: Peak Picking (Profile Mode)", xlabel="m/z", ylabel="Intensity")
mz = 1:200
base_signal = 10 .* exp.(-((mz .- 50).^2) ./ (2*3^2)) .+ 7 .* exp.(-((mz .- 120).^2) ./ (2*5^2)) .+ 2
small_peak_signal = zeros(200)
for i in 80:90
small_peak_signal[i] = 3 * exp(-((i - 85)^2) / 2.0)
end
noise = 0.5 .* randn(200)
intensity = base_signal .+ small_peak_signal .+ noise
noise_level = median(abs.(intensity .- median(intensity))) * 1.4826 # MAD
snr_threshold_val = 3.0
intensity_threshold = noise_level * snr_threshold_val
picked_peaks_mz = [50, 85, 120]
picked_peaks_intensity = intensity[picked_peaks_mz]
hlines!(ax, [noise_level], color=:gray, linestyle=:dot, label="Est. Noise Level")
hlines!(ax, [intensity_threshold], color=:orange, linestyle=:dash, label="SNR Threshold (snr_threshold = 3.0)")
lines!(ax, mz, intensity, color=:blue, label="Profile Spectrum")
scatter!(ax, picked_peaks_mz, picked_peaks_intensity, color=:green, markersize=15, strokewidth=2, label="Peaks passing SNR")
scatter!(ax, [25], [intensity[25]], color=:red, marker=:x, markersize=15, label="Local max below SNR")
# Illustrate prominence
arrows!(ax, [120, 120], [intensity[135], intensity[120]], [0, 0], [intensity[120]-intensity[135], 0], color=:purple)
text!(ax, 125, (intensity[120]+intensity[135])/2, text="Prominence", color=:purple)
axislegend(ax, position=:rt)
save("descriptive_plot_peakpicking_profile.png", fig)
println("Saved: descriptive_plot_peakpicking_profile.png")
end
"""
Visualizes peak picking (filtering) for centroid-mode data based on an SNR
(intensity) threshold.
"""
function plot_peak_picking_centroid_details()
fig = Figure(size=(1200, 700), fontsize=14)
ax = Axis(fig[1, 1], title="Detailed View: Peak Picking (Centroid Mode)", xlabel="m/z", ylabel="Intensity")
mz = [100, 150, 200, 250, 300, 350, 400]
intensity = [10, 5, 25, 8, 3, 18, 12]
snr_threshold_val = 10.0
stem!(ax, mz, intensity, color=:gray, label="Input Centroids")
selected_mask = intensity .>= snr_threshold_val
stem!(ax, mz[selected_mask], intensity[selected_mask], color=:green, trunkwidth=3, label="Selected Peaks (intensity >= 10)")
stem!(ax, mz[.!selected_mask], intensity[.!selected_mask], color=:red, trunkwidth=3, label="Rejected Peaks (intensity < 10)")
hlines!(ax, [snr_threshold_val], color=:orange, linestyle=:dash, label="Intensity Threshold (snr_threshold)")
text!(ax, "Parameter: `snr_threshold`\nIn centroid mode, this acts as a direct intensity filter.",
position=Point2f(0.05, 0.95), space=:relative, align=(:left, :top), fontsize=12)
axislegend(ax, position=:rt)
save("descriptive_plot_peakpicking_centroid.png", fig)
println("Saved: descriptive_plot_peakpicking_centroid.png")
end
"""
Plots the effect of different normalization methods on a set of spectra.
"""
function plot_normalization_details()
fig = Figure(size=(1200, 700), fontsize=14)
mz = range(300, 400, length=500)
spec1 = 1.5 .* exp.(-((mz .- 350).^2) ./ 50)
spec2 = 0.8 .* exp.(-((mz .- 350).^2) ./ 50)
ax1 = Axis(fig[1, 1], title="Before Normalization", ylabel="Absolute Intensity")
lines!(ax1, mz, spec1, label="Spectrum A (High TIC)")
lines!(ax1, mz, spec2, label="Spectrum B (Low TIC)")
axislegend(ax1)
ax2 = Axis(fig[1, 2], title="After Normalization (TIC)", ylabel="Relative Intensity")
lines!(ax2, mz, spec1 ./ sum(spec1), label="Spectrum A (Normalized)")
lines!(ax2, mz, spec2 ./ sum(spec2), label="Spectrum B (Normalized)")
text!(ax2, "Effect: Spectra are scaled to have the same total area, making their intensities comparable.",
position=Point2f(0.05, 0.95), space=:relative, align=(:left, :top), fontsize=12, justification=:left)
save("descriptive_plot_normalization.png", fig)
println("Saved: descriptive_plot_normalization.png")
end
"""
Illustrates mass calibration, showing how a calibration curve corrects measured
m/z values based on reference peaks.
"""
function plot_calibration_details()
fig = Figure(size=(1200, 700), fontsize=14)
ax = Axis(fig[1, 1], title="Detailed View: Calibration", xlabel="Measured m/z", ylabel="m/z Error (Measured - Reference)")
ref_mz = [200, 400, 600, 800, 1000]
measured_mz = ref_mz .+ [0.1, 0.15, 0.2, 0.25, 0.3] .+ 0.02 .* randn(5)
errors = measured_mz .- ref_mz
# Fit a linear model to the error
A = [ones(5) measured_mz]
coeffs = A \ errors
correction_func(m) = m - (coeffs[1] .+ coeffs[2] .* m)
fit_line = coeffs[1] .+ coeffs[2] .* measured_mz
scatter!(ax, measured_mz, errors, color=:red, markersize=15, label="Measured Error")
lines!(ax, measured_mz, fit_line, color=:blue, label="Calibration Curve (fit_order=1)")
text!(ax, "Parameter: `fit_order`\nA curve is fit to the error of known reference peaks.\nThis curve is then used to correct all m/z values.",
position=Point2f(0.05, 0.95), space=:relative, align=(:left, :top), fontsize=12)
axislegend(ax, position=:rb)
save("descriptive_plot_calibration.png", fig)
println("Saved: descriptive_plot_calibration.png")
end
"""
Visualizes peak alignment by showing multiple spectra with misaligned peaks
before and after the alignment process.
"""
function plot_alignment_details()
fig = Figure(size=(1600, 600), fontsize=14)
ax1 = Axis(fig[1, 1], title="Before Alignment", xlabel="m/z", yticklabelsvisible=false, ygridvisible=false)
ax2 = Axis(fig[1, 2], title="After Alignment", xlabel="m/z", yticklabelsvisible=false, ygridvisible=false)
mz = range(490, 510, length=1000)
shifts = [-0.5, 0.0, 0.8]
colors = [:blue, :green, :purple]
for (i, shift) in enumerate(shifts)
peak_center = 500 + shift
spectrum = exp.(-((mz .- peak_center).^2) ./ 0.1)
lines!(ax1, mz, spectrum .+ i, color=colors[i])
vlines!(ax1, [peak_center], color=(colors[i], 0.5), linestyle=:dash)
# After alignment, all peaks are at 500
aligned_spectrum = exp.(-((mz .- 500).^2) ./ 0.1)
lines!(ax2, mz, aligned_spectrum .+ i, color=colors[i])
end
vlines!(ax2, [500], color=:red, linestyle=:dash, label="Reference m/z")
axislegend(ax2)
save("descriptive_plot_alignment.png", fig)
println("Saved: descriptive_plot_alignment.png")
end
"""
Illustrates peak selection by filtering a population of peaks based on
FWHM and SNR criteria.
"""
function plot_peak_selection_details()
fig = Figure(size=(1200, 700), fontsize=14)
ax = Axis(fig[1, 1], title="Detailed View: Peak Selection", xlabel="FWHM (ppm)", ylabel="Signal-to-Noise Ratio (SNR)")
n_peaks = 100
fwhm = rand(n_peaks) .* 150
snr = rand(n_peaks) .* 20
min_fwhm_ppm = 20.0
max_fwhm_ppm = 100.0
min_snr = 5.0
selected_mask = (fwhm .>= min_fwhm_ppm) .& (fwhm .<= max_fwhm_ppm) .& (snr .>= min_snr)
scatter!(ax, fwhm[.!selected_mask], snr[.!selected_mask], color=(:red, 0.5), label="Rejected Peaks")
scatter!(ax, fwhm[selected_mask], snr[selected_mask], color=:green, label="Selected Peaks")
vlines!(ax, [min_fwhm_ppm, max_fwhm_ppm], color=:blue, linestyle=:dash, label="FWHM bounds")
hlines!(ax, [min_snr], color=:orange, linestyle=:dash, label="SNR bound")
poly!(ax, BBox(min_fwhm_ppm, max_fwhm_ppm, min_snr, 22), color=(:green, 0.1))
text!(ax, "Selection Region", position=(60, 12), color=:green, fontsize=14)
axislegend(ax)
save("descriptive_plot_peak_selection.png", fig)
println("Saved: descriptive_plot_peak_selection.png")
end
"""
Visualizes the adaptive peak binning process, showing how peaks from different
spectra are grouped into a common bin based on a PPM tolerance.
"""
function plot_peak_binning_details()
fig = Figure(size=(1200, 700), fontsize=14)
ax = Axis(fig[1, 1], title="Detailed View: Adaptive Peak Binning", xlabel="m/z", yticklabelsvisible=false)
ref_mz = 500.0
tolerance_ppm = 50.0
tol_mz = ref_mz * tolerance_ppm / 1e6
bin_start = ref_mz - tol_mz/2
bin_end = ref_mz + tol_mz/2
peaks_mz = [ref_mz - 0.01, ref_mz + 0.005, ref_mz + 0.02, ref_mz - 0.015]
peak_intensities = [0.8, 1.0, 0.9, 0.7]
peak_colors = [:blue, :green, :purple, :orange]
vspan!(ax, bin_start, bin_end, color=(:gray, 0.2), label="Bin (tolerance: 50 ppm)")
stem!(ax, peaks_mz, peak_intensities, color=peak_colors, markersize=15)
# Show bin center
bin_center = mean(peaks_mz)
vlines!(ax, [bin_center], color=:red, linestyle=:dash, label="Calculated Bin Center")
text!(ax, "Parameter: `tolerance`\nPeaks from different spectra within the tolerance window are grouped into a single feature.",
position=Point2f(0.05, 0.95), space=:relative, align=(:left, :top), fontsize=12)
axislegend(ax, position=:rt)
xlims!(ax, ref_mz - tol_mz*2, ref_mz + tol_mz*2)
save("descriptive_plot_peak_binning.png", fig)
println("Saved: descriptive_plot_peak_binning.png")
end
"""
Main function to generate and save all descriptive plots.
"""
function create_and_save_all_plots()
println("Generating detailed descriptive plots for preprocessing steps...")
plot_baseline_correction_details()
plot_smoothing_details()
plot_peak_picking_profile_details()
plot_peak_picking_centroid_details()
plot_normalization_details()
plot_calibration_details()
plot_alignment_details()
plot_peak_selection_details()
plot_peak_binning_details()
println("\nAll descriptive plots have been saved in the current directory.")
end
# Execute the plot generation
if abspath(PROGRAM_FILE) == @__FILE__
create_and_save_all_plots()
end

157
test/run_precalculation_dict.jl
View File

@ -1,157 +0,0 @@
# test/run_precalculation_example.jl
using Printf
import Pkg
# --- Load the MSI_src Module ---
Pkg.activate(joinpath(@__DIR__, ".."))
using MSI_src
# ===================================================================
# CONFIG: PLEASE FILL IN YOUR FILE PATHS HERE
# ===================================================================
# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML"
#const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/mzML/Col_1.mzML"
const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Thricoderma_etc/Imaging_interaccion_trichoderma_vs_streptomyces.mzML"
# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.imzML"
#const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML"
const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Thricoderma_etc/Imaging_interaccion_trichoderma_vs_streptomyces.imzML"
#const MASK_ROUTE = "/home/pixel/Documents/Cinvestav_2025/JuliaMSI/public/css/masks/Stomach_DHB_uncompressed.png"
const MASK_ROUTE = ""
const reference_peaks = Dict(
# DHB Matrix peaks (should be present)
137.0244 => "DHB_fragment",
155.0349 => "DHB_M+H",
177.0168 => "DHB_M+Na",
# Common lipids in your mass range
496.3398 => "PC_16:0_16:0",
520.3398 => "PC_16:0_18:1",
760.5851 => "PC_16:0_18:1_Na",
# Common contaminants
391.2843 => "PDMS",
413.2662 => "PDMS_Na",
# Add some high mass peaks
842.5092 => "Protein_standard",
1045.532 => "Protein_standard",
290.1747 => "Atropine [M+H]+",
304.1903 => "Scopolamine [M+H]+",
124.0393 => "Tropine [M+H]+",
)
# ===================================================================
# HELPER FUNCTIONS FOR PRINTING
# ===================================================================
function print_header(title::String)
println("\n" * "="^80)
println("$(title)")
println("="^80)
end
function check_data_range(msi_data::MSIData)
println("\n--- Data Range Analysis ---")
min_mz, max_mz = get_global_mz_range(msi_data)
if isfinite(min_mz) && isfinite(max_mz) && min_mz < max_mz
println("Global m/z range: [$(min_mz), $(max_mz)]")
else
println("Global m/z range: Not yet determined or invalid (initial: [$(min_mz), $(max_mz)])")
end
# Check a few spectra to see actual m/z values
println("\nChecking first few spectra for actual m/z values:")
for i in 1:min(3, length(msi_data.spectra_metadata))
try
mz, intensity = GetSpectrum(msi_data, i)
if !isempty(mz)
println("Spectrum $i: m/z range [$(minimum(mz)), $(maximum(mz))], length=$(length(mz))")
# Print first and last few m/z values
if length(mz) > 10
println(" First 5 m/z: $(mz[1:5])")
println(" Last 5 m/z: $(mz[end-4:end])")
end
end
catch e
println("Spectrum $i: Error - $e")
end
end
end
# ===================================================================
# MAIN EXAMPLE RUNNER
# ===================================================================
function run_precalculation_example()
# --- Process mzML file ---
print_header("Processing mzML File: $(basename(TEST_MZML_FILE))")
if !isfile(TEST_MZML_FILE)
println("SKIPPING: mzML file not found at $(TEST_MZML_FILE)")
else
try
msi_data_mzml = @time OpenMSIData(TEST_MZML_FILE)
check_data_range(msi_data_mzml)
analysis_results_mzml = main_precalculation(msi_data_mzml, reference_peaks=reference_peaks)
println("\n" * "*"^80)
println("MZML PREPROCESSING ANALYSIS RESULTS")
println("*"^80)
for (step_name, params) in analysis_results_mzml
println("\n--- $(uppercase(string(step_name))) Parameters ---")
if isempty(params)
println(" No recommended parameters.")
else
for (param_name, param_value) in params
println(" $param_name: $param_value")
end
end
end
close(msi_data_mzml) # Close file handles
catch e
println("ERROR processing mzML file: $e")
showerror(stdout, e, catch_backtrace())
end
end
# --- Process imzML file ---
print_header("Processing imzML File: $(basename(TEST_IMZML_FILE))")
if !isfile(TEST_IMZML_FILE)
println("SKIPPING: imzML file not found at $(TEST_IMZML_FILE)")
else
try
msi_data_imzml = @time OpenMSIData(TEST_IMZML_FILE)
check_data_range(msi_data_imzml)
analysis_results_imzml = main_precalculation(msi_data_imzml, reference_peaks=reference_peaks, mask_path=MASK_ROUTE)
println("\n" * "*"^80)
println("IMZML PREPROCESSING ANALYSIS RESULTS")
println("*"^80)
for (step_name, params) in analysis_results_imzml
println("\n--- $(uppercase(string(step_name))) Parameters ---")
if isempty(params)
println(" No recommended parameters.")
else
for (param_name, param_value) in params
println(" $param_name: $param_value")
end
end
end
close(msi_data_imzml) # Close file handles
catch e
println("ERROR processing imzML file: $e")
showerror(stdout, e, catch_backtrace())
end
end
end
# --- Execute ---
@time run_precalculation_example()

282
test/run_precalculation_example.jl
View File

@ -1,282 +0,0 @@
# test/run_precalculation_example.jl
using Printf
import Pkg
# --- Load the MSI_src Module ---
Pkg.activate(joinpath(@__DIR__, ".."))
using MSI_src
# ===================================================================
# CONFIG: PLEASE FILL IN YOUR FILE PATHS HERE
# ===================================================================
# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML"
# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/mzML/Col_1.mzML"
const TEST_MZML_FILE = ""
# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.imzML"
#const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML"
# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Thricoderma_etc/Imaging_interaccion_trichoderma_vs_streptomyces.imzML"
#const MASK_ROUTE = "/home/pixel/Documents/Cinvestav_2025/JuliaMSI/public/css/masks/Stomach_DHB_uncompressed.png"
const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML"
const MASK_ROUTE = ""
#=
const reference_peaks = Dict(
# Common ESI positive mode reference compounds
121.0509 => "Purine",
149.0233 => "HP-921",
322.0481 => "Hexakis(1H,1H,3H-tetrafluoropropoxy)phosphazine",
622.0290 => "Hexakis(2,2-difluoroethoxy)phosphazine",
# Atropine and related compounds in positive mode
290.1747 => "Atropine [M+H]+",
304.1903 => "Scopolamine [M+H]+",
124.0393 => "Tropine [M+H]+",
# Common contaminants and lock masses
391.2843 => "Polydimethylcyclosiloxane [M+H]+",
413.2662 => "Polydimethylcyclosiloxane [M+Na]+",
429.2402 => "Polydimethylcyclosiloxane [M+K]+"
)
=#
const reference_peaks = Dict(
# DHB Matrix peaks (should be present)
137.0244 => "DHB_fragment",
155.0349 => "DHB_M+H",
177.0168 => "DHB_M+Na",
# Common lipids in your mass range
496.3398 => "PC_16:0_16:0",
520.3398 => "PC_16:0_18:1",
760.5851 => "PC_16:0_18:1_Na",
# Common contaminants
391.2843 => "PDMS",
413.2662 => "PDMS_Na",
# Add some high mass peaks
842.5092 => "Protein_standard",
1045.532 => "Protein_standard",
)
# ===================================================================
# HELPER FUNCTIONS FOR PRINTING
# ===================================================================
function print_header(title::String)
println("\n" * "="^80)
println("$(title)")
println("="^80)
end
function print_subheader(title::String)
println("\n" * "-"^80)
println("$(title)")
println("-"^80)
end
function print_param(key, value)
if value === nothing || (isa(value, Number) && isnan(value))
println(" " * "" * rpad(key, 30) * ": unable to get, user needs to input manually")
elseif isa(value, Symbol) && startswith(string(key), "method") # Heuristic for method selection
println(" " * "" * rpad(key, 30) * ": automatically detected this as the most optimal method: $(value)")
else
println(" " * "" * rpad(key, 30) * ": $(value)")
end
end
function print_recommendations(recommendations::Dict)
for (step, params) in recommendations
print_subheader("Recommendations for $(String(step))")
for (key, value) in params
print_param(key, value)
end
end
end
function check_data_range(msi_data::MSIData)
println("\n--- Data Range Analysis ---")
min_mz, max_mz = get_global_mz_range(msi_data)
if isfinite(min_mz) && isfinite(max_mz) && min_mz < max_mz
println("Global m/z range: [$(min_mz), $(max_mz)]")
else
println("Global m/z range: Not yet determined or invalid (initial: [$(min_mz), $(max_mz)])")
end
# Check a few spectra to see actual m/z values
println("\nChecking first few spectra for actual m/z values:")
for i in 1:min(3, length(msi_data.spectra_metadata))
try
mz, intensity = GetSpectrum(msi_data, i)
if !isempty(mz)
println("Spectrum $i: m/z range [$(minimum(mz)), $(maximum(mz))], length=$(length(mz))")
# Print first and last few m/z values
if length(mz) > 10
println(" First 5 m/z: $(mz[1:5])")
println(" Last 5 m/z: $(mz[end-4:end])")
end
end
catch e
println("Spectrum $i: Error - $e")
end
end
end
# ===================================================================
# MAIN EXAMPLE RUNNER
# ===================================================================
function run_precalculation_example()
# --- Process mzML file ---
print_header("Processing mzML File: $(basename(TEST_MZML_FILE))")
if !isfile(TEST_MZML_FILE)
println("SKIPPING: mzML file not found at $(TEST_MZML_FILE)")
else
try
msi_data_mzml = @time OpenMSIData(TEST_MZML_FILE)
check_data_range(msi_data_mzml)
analysis_results_mzml = run_preprocessing_analysis(msi_data_mzml, reference_peaks=reference_peaks)
println("\n" * "*"^80)
println("MZML PREPROCESSING ANALYSIS RESULTS")
println("*"^80)
for (phase, results) in analysis_results_mzml
if phase == :recommendations
print_subheader("Generated Preprocessing Recommendations")
print_recommendations(results)
else
print_subheader("Phase: $(String(phase))")
if isa(results, Dict)
for (key, value) in results
print_param(key, value)
end
elseif isa(results, NamedTuple)
for field in fieldnames(typeof(results))
print_param(field, getfield(results, field))
end
else
println(" $(String(phase)) results: $(results)")
end
end
end
close(msi_data_mzml) # Close file handles
catch e
println("ERROR processing mzML file: $e")
showerror(stdout, e, catch_backtrace())
end
end
# --- Process imzML file ---
print_header("Processing imzML File: $(basename(TEST_IMZML_FILE))")
if !isfile(TEST_IMZML_FILE)
println("SKIPPING: imzML file not found at $(TEST_IMZML_FILE)")
else
try
msi_data_imzml = @time OpenMSIData(TEST_IMZML_FILE)
check_data_range(msi_data_imzml)
analysis_results_imzml = run_preprocessing_analysis(msi_data_imzml, reference_peaks=reference_peaks, mask_path=MASK_ROUTE)
println("\n" * "*"^80)
println("IMZML PREPROCESSING ANALYSIS RESULTS")
println("*"^80)
for (phase, results) in analysis_results_imzml
if phase == :recommendations
print_subheader("Generated Preprocessing Recommendations")
print_recommendations(results)
else
print_subheader("Phase: $(String(phase))")
if isa(results, Dict)
for (key, value) in results
print_param(key, value)
end
elseif isa(results, NamedTuple)
for field in fieldnames(typeof(results))
print_param(field, getfield(results, field))
end
else
println(" $(String(phase)) results: $(results)")
end
end
end
close(msi_data_imzml) # Close file handles
catch e
println("ERROR processing imzML file: $e")
showerror(stdout, e, catch_backtrace())
end
end
end
# --- Execute ---
@time run_precalculation_example()
# =============================================================================
# Example Usage
# =============================================================================
#=
"""
example_preanalysis_workflow(msi_data::MSIData)
Demonstrates how to run the pre-analysis pipeline.
"""
function example_preanalysis_workflow(msi_data::MSIData)
# Load your MSIData object (this would come from your actual data loading)
# msi_data = load_imzml_dataset("path/to/your/data.imzML") # Assuming msi_data is passed
# Define reference peaks for mass accuracy analysis
reference_peaks = Dict(
89.04767 => "Alanin",
147.07642 => "Lysin",
189.12392 => "Unknown",
524.26496 => "PC(34:1) [M+H]+"
)
# Define region masks if you have spatial annotations
region_masks = Dict{Symbol, BitMatrix}()
# region_masks[:tumor] = load_mask("tumor_mask.png") # Not defined, comment out
# region_masks[:stroma] = load_mask("stroma_mask.png") # Not defined, comment out
println("Starting comprehensive pre-analysis...")
# Run the complete pre-analysis pipeline
analysis_results = run_preprocessing_analysis(
msi_data, # Your MSIData object
reference_peaks=reference_peaks,
region_masks=region_masks,
sample_size=200 # Adjust based on dataset size
)
# Access the recommendations
recommendations = analysis_results[:recommendations]
println("\n" * "="^60)
println("PREPROCESSING RECOMMENDATIONS")
println("="^60)
for (step, params) in recommendations
println("\n$step:")
for (key, value) in params
println(" - $key: $value")
end
end
return analysis_results
end
# You can also run individual analysis steps:
function run_targeted_analysis(msi_data::MSIData)
# Just analyze signal quality and peak characteristics
signal_analysis = analyze_signal_quality(msi_data, sample_size=100)
peak_analysis = analyze_peak_characteristics(msi_data, sample_size=50)
return (signal_analysis, peak_analysis)
end
=#

1035
test/run_preprocessing.jl
View File

File diff suppressed because it is too large Load Diff

14
test/run_tests.jl
View File

@ -31,10 +31,9 @@ using MSI_src
# --- Test Case 1: Standard .mzML file --- # --- Test Case 1: Standard .mzML file ---
# A regular, non-imaging mzML file. # A regular, non-imaging mzML file.
# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/mzML/T9_A1.mzML" # const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/mzML/T9_A1.mzML"
# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML" const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML"
# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging_paper_spray/Imaging_paper_spray.mzML" # const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging_paper_spray/Imaging_paper_spray.mzML"
# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging prueba Roya 1/Roya.mzML" # const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging prueba Roya 1/Roya.mzML"
const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/mzML"
const SPECTRUM_TO_PLOT = 1 # Which spectrum to plot from the file const SPECTRUM_TO_PLOT = 1 # Which spectrum to plot from the file
# --- Test Case 2: .mzML + Sync File for Conversion --- # --- Test Case 2: .mzML + Sync File for Conversion ---
@ -59,17 +58,16 @@ const CONVERSION_TARGET_IMZML = "test/results/converted_mzml.imzML"
# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging prueba Roya 1/royaimg.imzML" # const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging prueba Roya 1/royaimg.imzML"
# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/ltpmsi-chilli.imzML" # centroid aparently? # const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/ltpmsi-chilli.imzML" # centroid aparently?
# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_compressed.imzML" # centroid compressed # const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_compressed.imzML" # centroid compressed
# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML" # centroid const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML" # centroid
const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML"
# The m/z value to use for creating an image slice. # The m/z value to use for creating an image slice.
# const MZ_VALUE_FOR_SLICE = 309.06 # BF # const MZ_VALUE_FOR_SLICE = 309.06 # BF
const MZ_VALUE_FOR_SLICE = 896.0 # HR2MSI # const MZ_VALUE_FOR_SLICE = 896.0 # HR2MSI
# const MZ_VALUE_FOR_SLICE = 76.03 # I PS # const MZ_VALUE_FOR_SLICE = 76.03 # I PS
# const MZ_VALUE_FOR_SLICE = 313 # ROYA # const MZ_VALUE_FOR_SLICE = 313 # ROYA
# const MZ_VALUE_FOR_SLICE = 100 # advanced processing const MZ_VALUE_FOR_SLICE = 100 # advanced processing
# const MZ_TOLERANCE = 0.1 # const MZ_TOLERANCE = 0.1
# const MZ_TOLERANCE = 1 # const MZ_TOLERANCE = 1
const MZ_TOLERANCE = 0.2 const MZ_TOLERANCE = 0.1
# Coordinates to plot a specific spectrum from imzML # Coordinates to plot a specific spectrum from imzML
const COORDS_TO_PLOT = (50, 50) # Example coordinates (X, Y) const COORDS_TO_PLOT = (50, 50) # Example coordinates (X, Y)
@ -77,7 +75,7 @@ const COORDS_TO_PLOT = (50, 50) # Example coordinates (X, Y)
# --- Output Directory --- # --- Output Directory ---
const RESULTS_DIR = "test/results" const RESULTS_DIR = "test/results"
test1 = true test1 = false
test2 = false test2 = false
test3 = true test3 = true

58
test/runtests.jl
View File

@ -1,58 +0,0 @@
using Test
# We need to simulate loading the MSI_src module directly
include("../src/MSI_src.jl")
using .MSI_src
@testset "JuliaMSI Core Systems" begin
@testset "Basic Instantiation" begin
# 1. Test basic metadata structure initialization
meta = SpectrumMetadata(
Int32(1), Int32(2),
"test_id",
:sample,
CENTROID,
SpectrumAsset(Float64, false, Int64(0), 100, :mz, 0.0, 0.0),
SpectrumAsset(Float32, true, Int64(100), 50, :intensity, 0.0, 0.0)
)
@test meta.x == 1
@test meta.y == 2
@test meta.mode == CENTROID
@test meta.mz_asset.format == Float64
@test meta.int_asset.format == Float32
@test meta.int_asset.is_compressed == true
# 2. Test cache pool and MSIData structural stability
source = MzMLSource([], Float64, Float32, nothing) # Empty source for structural testing
# Test constructor doesn't throw
# Constructor signature: (source, metadata, instrument_meta, dims, coordinate_map, cache_size)
msi_data = MSIData(
source,
[meta],
nothing, # instrument_meta
(100, 100), # dims
nothing, # coord map
10 # cache size
)
@test msi_data.image_dims == (100, 100)
@test length(msi_data.spectra_metadata) == 1
@test msi_data.cache_size == 10
end
@testset "Buffer & Cache Subsystems" begin
pool = SimpleBufferPool()
# Test basic retrieval
buf = get_buffer!(pool, 1024)
@test length(buf) == 1024
# Test release
release_buffer!(pool, buf)
@test length(pool.buffers[1024]) == 1
# Test reuse
buf2 = get_buffer!(pool, 1024)
@test buf === buf2 # Should return the EXACT same buffer object
end
end

32
test/test_buffer_pool_concurrency.jl
View File

@ -1,32 +0,0 @@
using MSI_src
using Test
using Base.Threads
@testset "SimpleBufferPool Concurrency Stress Test" begin
pool = MSI_src.SimpleBufferPool()
n_threads = Threads.nthreads()
n_iterations = 1000
buffer_size = 1024
println("Running buffer pool stress test with $n_threads threads...")
# Parallel stress test
Threads.@threads for i in 1:(n_threads * n_iterations)
# get_buffer! and release_buffer! are now thread-safe
buf = MSI_src.get_buffer!(pool, buffer_size)
# Simulate some work
fill!(buf, UInt8(i % 256))
MSI_src.release_buffer!(pool, buf)
end
# After stress test, the dictionary should be coherent
sizes = collect(keys(pool.buffers))
if !isempty(sizes)
@test buffer_size ∈ sizes
@test length(pool.buffers[buffer_size]) <= pool.max_pool_size
end
println("Buffer pool stress test PASSED.")
end

146
test/test_pipeline_concurrency.jl
View File

@ -1,146 +0,0 @@
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
using MSI_src
using Test
using Base.Threads
# --- Test Configuration ---
const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Stomach/Stomach_DHB_uncompressed.imzML"
println("Starting concurrency test with file: $TEST_FILE")
if !isfile(TEST_FILE)
error("Test file not found: $TEST_FILE")
end
# --- Simulation of the Concurrency Issue ---
# Global "UI State" variable simulates the app's msi_data
global_msi_data = nothing
const GLOBAL_LOCK = ReentrantLock()
function ui_load_file(path)
global global_msi_data
lock(GLOBAL_LOCK) do
if global_msi_data !== nothing
close(global_msi_data)
end
println("UI: Loading file...")
global_msi_data = OpenMSIData(path)
println("UI: File loaded.")
end
end
function ui_close_file()
global global_msi_data
lock(GLOBAL_LOCK) do
if global_msi_data !== nothing
println("UI: Closing file...")
close(global_msi_data)
global_msi_data = nothing
println("UI: File closed.")
end
end
end
# Mock pipeline function that mimics `run_full_pipeline` in app.jl
# Crucially, it uses its own local `pipeline_msi_data` as per the fix
function run_mock_pipeline_isolated(path)
println("Pipeline: Starting isolated pipeline...")
# 1. Open ISOLATED instance
pipeline_msi_data = OpenMSIData(path)
println("Pipeline: Opened isolated MSIData instance.")
try
# 2. Simulate reading spectra
indices = 1:min(100, length(pipeline_msi_data.spectra_metadata)) # Read first 100 spectra
# Artificial delay to allow "user" interaction
sleep(0.5)
println("Pipeline: Reading spectra...")
Threads.@threads for i in indices
# Use local instance
mz, int = GetSpectrum(pipeline_msi_data, i)
# Simulate processing work
sum(int)
end
println("Pipeline: Finished reading spectra successfully.")
return true
catch e
println("Pipeline: CRASHED with error: $e")
return false
finally
close(pipeline_msi_data)
println("Pipeline: Closed isolated MSIData instance.")
end
end
# Mock pipeline that uses GLOBAL instance (The BUGGY version)
function run_mock_pipeline_buggy()
println("Buggy Pipeline: Starting...")
# Uses global_msi_data directly
try
global global_msi_data
if global_msi_data === nothing
println("Buggy Pipeline: No data loaded!")
return false
end
local_ref = global_msi_data # Still points to same object
indices = 1:min(100, length(local_ref.spectra_metadata))
sleep(0.5)
println("Buggy Pipeline: Reading spectra from shared object...")
Threads.@threads for i in indices
# This will fail if ui_close_file() happens concurrently
mz, int = GetSpectrum(local_ref, i)
sum(int)
end
println("Buggy Pipeline: Success (Unexpected if concurrency worked)")
return true
catch e
println("Buggy Pipeline: CRASHED as expected: $e")
return false
end
end
# --- execute checks ---
@testset "Concurrency Crash Fix Verification" begin
# 1. Setup: Load file initially
ui_load_file(TEST_FILE)
# 2. Test the FIX: Isolated Pipeline
println("\n--- Testing Fixed (Isolated) Pipeline ---")
t_pipeline = Threads.@spawn run_mock_pipeline_isolated(TEST_FILE)
# Simulate user closing/reloading file while pipeline runs
sleep(0.2)
ui_close_file()
# Wait for pipeline
success = fetch(t_pipeline)
@test success == true
println("Fixed pipeline result: ", success ? "PASSED" : "FAILED")
# 3. Test the BUG: Global Pipeline (Optional, to prove it crashes without fix)
# Uncomment to verify the bug exists if needed, but we assume it does based on user report.
# println("\n--- Testing Buggy (Shared) Pipeline ---")
# ui_load_file(TEST_FILE)
# t_buggy = Threads.@spawn run_mock_pipeline_buggy()
# sleep(0.2)
# ui_close_file()
# buggy_success = fetch(t_buggy)
# println("Buggy pipeline result: ", buggy_success ? "PASSED (No crash?)" : "FAILED (Crashed as expected)")
end

124
test/test_streaming_pipeline.jl
View File

@ -1,124 +0,0 @@
#!/usr/bin/env julia
# test/test_streaming_pipeline.jl
# ============================================================================
# Validation test for Sprint 2: The Streaming Pipeline
#
# This test exercises process_dataset! against the HR2MSI mouse bladder
# dataset and verifies:
# 1. Correct sparse matrix creation
# 2. Non-zero peak population
# 3. RAM savings vs dense equivalent
# 4. Allocation count and throughput
# ============================================================================
using Pkg
Pkg.activate(".")
using MSI_src
using SparseArrays
# =============================================================================
# Configuration
# =============================================================================
const IMZML_PATH = "/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML"
function main()
println("=" ^ 60)
println("SPRINT 2: Streaming Pipeline Validation")
println("=" ^ 60)
if !isfile(IMZML_PATH)
println("SKIPPED: Dataset not found at $IMZML_PATH")
return
end
# --- 1. Load dataset ---
println("\n--- Step 1: Loading dataset ---")
data = OpenMSIData(IMZML_PATH)
println("Loaded: $(length(data.spectra_metadata)) spectra")
# --- 2. Configure the streaming pipeline ---
println("\n--- Step 2: Configuring pipeline ---")
config = PipelineConfig(
steps = [
StreamingStep(:baseline_correction, Dict{Symbol,Any}(:method => :snip, :iterations => 50)),
StreamingStep(:normalization, Dict{Symbol,Any}(:method => :tic)),
StreamingStep(:peak_picking, Dict{Symbol,Any}(
:method => :profile,
:snr_threshold => 3.0,
:half_window => 10,
:min_peak_prominence => 0.1,
:merge_peaks_tolerance => 0.002
)),
],
num_bins = 2000,
frequency_threshold = 0.01 # Bins must appear in at least 1% of spectra
)
println("Steps: $(join([s.name for s in config.steps], " → "))")
println("Bins: $(config.num_bins), Frequency threshold: $(config.frequency_threshold)")
# --- 3. Run the streaming pipeline ---
println("\n--- Step 3: Running streaming pipeline ---")
stats = @timed begin
feature_matrix, bin_centers = process_dataset!(data, config)
end
feature_matrix = stats.value[1]
bin_centers = stats.value[2]
println("\n--- Results ---")
println(" Feature matrix size: $(size(feature_matrix))")
println(" Non-zeros: $(nnz(feature_matrix))")
println(" Bin centers: $(length(bin_centers))")
println(" Time: $(round(stats.time, digits=2))s")
println(" Allocations: $(stats.bytes ÷ 1_000_000) MB")
println(" GC time: $(round(stats.gctime, digits=2))s")
# --- 4. Validate ---
println("\n--- Step 4: Validation ---")
passed = true
# Check matrix dimensions
if size(feature_matrix, 1) > 0 && size(feature_matrix, 2) > 0
println(" ✓ Matrix has valid dimensions")
else
println(" ✗ Matrix has invalid dimensions: $(size(feature_matrix))")
passed = false
end
# Check non-zeros
if nnz(feature_matrix) > 0
println(" ✓ Matrix has $(nnz(feature_matrix)) non-zero entries")
else
println(" ✗ Matrix is completely empty")
passed = false
end
# Check sparsity savings
dense_mb = size(feature_matrix, 1) * size(feature_matrix, 2) * 8 / 1e6
sparse_mb = nnz(feature_matrix) * 16 / 1e6 # index + value per entry
if dense_mb > 0
savings = (1.0 - sparse_mb / dense_mb) * 100
println(" ✓ RAM savings: $(round(savings, digits=1))% ($(round(sparse_mb, digits=1)) MB vs $(round(dense_mb, digits=1)) MB dense)")
end
# Check bin centers alignment
if length(bin_centers) == size(feature_matrix, 1)
println(" ✓ Bin centers match matrix rows")
else
println(" ✗ Bin center count ($(length(bin_centers))) != matrix rows ($(size(feature_matrix, 1)))")
passed = false
end
println("\n" * "=" ^ 60)
if passed
println("ALL VALIDATIONS PASSED ✓")
else
println("SOME VALIDATIONS FAILED ✗")
end
println("=" ^ 60)
end
@time main()