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

Merge pull request 'dev/PreprocessingInUI' (#1) from dev/PreprocessingInUI into main

Browse Source 
Reviewed-on: https://codeberg.org/pixelguy14/JuliaMSI/pulls/1
This commit is contained in:
José Julián Sierra Álvarez 2026-01-20 19:26:54 +01:00
commit 4414671a30
36 changed files with 18914 additions and 2021 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
.gitignore vendored
View File

@ -1,6 +1,10 @@
public/*
public/css/imgOver.png
!public/css/
!public/css/autogenerated.css
!public/css/LABI_logo.png
!public/css/masks/static.txt
public/css/masks/*
log/*
R_original_scripts/
test/results/
test/binarization_results/
Manifest.toml

1364
Manifest.toml
View File

File diff suppressed because it is too large Load Diff

82
Project.toml
View File

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

3
README.md
View File

@ -6,6 +6,7 @@ https://codeberg.org/LabABI/JuliaMSI
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).
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
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,7 +23,7 @@ https://codeberg.org/LabABI/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):
```
julia --project=. start_MSI_GUI.jl
julia --threads auto --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.

3786
app.jl
View File

File diff suppressed because it is too large Load Diff

838
app.jl.html
View File

File diff suppressed because it is too large Load Diff

3
config/env/global.jl vendored
View File

@ -1 +1,2 @@
ENV["GENIE_ENV"] = "prod"
ENV["GENIE_ENV"] = "dev"
#ENV["GENIE_ENV"] = "prod"

6
config/env/prod.jl vendored Normal file
View File

@ -0,0 +1,6 @@
# 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 = []

1192
julia_imzML_visual.jl
View File

File diff suppressed because it is too large Load Diff

1021
mask.jl Normal file
View File

File diff suppressed because it is too large Load Diff

171
mask.jl.html Normal file
View File

@ -0,0 +1,171 @@
<header id="header">
<img src="/css/LABI_logo.png" alt="Labi Logo Icon" id="imgLogo">
<div>
<h4>JuliaMSI - Mask Editor&nbsp;</h4>
</div>
</header>
<div id="extDivStyle" class="row col-12 q-pa-xl">
<div class="row col-4">
<!-- Left Panel: Controls -->
<div class="st-col col-12 st-module q-pa-md">
<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" v-on:focus="refetch_folders = true"></q-select>
</div>
<p class="text-center" v-html="msgimg"></p>
<!-- Step 2: Provide Mask Input -->
<div class="text-subtitle1 q-mt-md">Step 2: Create/Upload Mask</div>
<div class="row justify-around">
<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>
<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>
<q-btn class="q-ma-sm btn-style" v-on:click="btn_upload_premade_mask = true" label="Upload Pre-made Mask"
:disable="!is_browsing_slices"></q-btn>
<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>
</div>
<!-- Step 3: Automated Processing & Overlay -->
<div class="text-subtitle1 q-mt-md">Step 3: Adjust & Verify</div>
<div class="q-mt-sm">
<p class="text-caption text-center">Adjust automated processing for the initial mask:</p>
<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 class="row col-6">
<!-- Right Panel: Image Displays -->
<div id="intDivStyle-right" class="st-col col-12 st-module">
<div class="text-h6">Mask Preview</div>
<div class="row items-center">
<q-btn icon="arrow_back" class="q-my-sm btn-style" v-on:click="btn_img_minus=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 v-if="!show_verification_plot" class="text-center text-grey q-pa-xl">
<q-icon name="image" size="5em" class="q-mb-md" />
<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>
</div>
<div v-if="show_verification_plot">
<plotly :data="plotdata_verify" :layout="plotlayout_verify" class="q-pa-none q-ma-none pixelated-plot">
</plotly>
<div class="q-mt-md text-center">
<q-btn color="primary" class="btn-style" label="Edit Manually"
v-on:click="btn_edit_manually=true" />
</div>
</div>
<p class="text-center" v-html="msgimg"></p>
</div>
</div>
</div>
<q-dialog v-model="show_editor">
<q-card style="min-width: 80vw; min-height: 80vh;">
<q-card-section class="row items-center q-pb-none">
<div class="text-h6">Manual Mask Editor</div>
<q-space />
<q-btn icon="close" flat round dense v-on:click="syncMaskToServer()" v-close-popup />
</q-card-section>
<q-card-section class="scroll" style="max-height: 70vh;">
<!-- Editor Tools -->
<div class="q-pa-md q-gutter-sm row justify-center items-center">
<!-- Tool Selection -->
<q-btn-group>
<q-btn label="Brush" v-on:click="current_tool = 'brush'"
: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="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-separator vertical inset />
<!-- Tool Settings -->
<q-input dense filled label="Brush Size" type="number" v-model.number="brush_size"
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-separator vertical inset />
<!-- Image Manipulation -->
<q-input dense filled label="Rotate (deg)" type="number" v-model.number="rotate_degrees"
style="max-width: 150px" :step="90" />
<q-select dense filled label="Flip" v-model="flip_direction" :options="['horizontal', 'vertical']"
style="max-width: 150px" />
<q-btn label="Apply Flip" v-on:click="btn_flip_mask = !btn_flip_mask" />
</div>
<!-- Drag Tool Controls -->
<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;">
<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_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_forward" v-on:click="move_mask_payload = { direction: 'right', t: Date.now() }" />
</div>
<!-- Interactive Canvas -->
<div class="row justify-center q-pa-md">
<div style="position: relative; display: inline-block; overflow: auto; max-width: 100%;">
<!-- Added for scrolling very large images -->
<canvas id="maskCanvas" :width="canvas_width" :height="canvas_height" :style="{
width: canvas_width * editor_scale + 'px',
height: canvas_height * editor_scale + 'px',
border: '2px solid #ccc',
cursor: current_tool === 'drag' ? 'move' : 'crosshair',
imageRendering: 'pixelated'
}" v-on:mousedown.prevent="startDrawing" v-on:mousemove.prevent="draw"
v-on:mouseup.prevent="stopDrawing()" v-on:mouseleave.prevent="stopDrawing()">
</canvas>
</div>
</div>
</q-card-section>
<!--<q-card-actions align="right" class="q-pa-md">
<q-btn label="Cancel" color="negative" v-on:click="show_editor = false"/>
</q-card-actions>-->
</q-card>
</q-dialog>

27
ongoing_tasks.txt
View File

@ -1,27 +0,0 @@
PENDING
Rmsi & julia coherence with image creation and colorbar
colorbar revamp CURRENTLY ONGOING
Create function that makes the imzML from mzML files ?
Possibility to add multiple imzML to process at once ?
Add comparison image (rotate, transform, translate, transparency) CURRENTLY ONGOING
DONE
Quicker start for julia GUI
Better UI
refined style of page
improved responsiveness
UI functionality added
Code cleanse
better code readability
Better blind search for images
3d topology plots for images
2d contour plots for images
Image flip upside down
Measure for the time it takes creating image or plot
Rotate vertical images
Easier input for files
TrIQ default values
Comparative for two views
Even faster initial boot and subsectuential boot
Multiple spectra plot types
Plot creation for spectra Per pixel of image

129
public/css/genieapp.css
View File

@ -1,11 +1,11 @@
#btnStyle {
background: #009f90;
color: rgb(229, 236, 246);
.btn-style {
background: #009f90 !important;
color: rgb(229, 236, 246) !important;
transition: background-color 0.3s ease;
}
#btnStyle:hover {
background-color: #007f72;
.btn-style:hover {
background-color: #007f72 !important;
}
#imgInt {
@ -41,10 +41,129 @@
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 {
color: #009f90;
}
#intDivStyle-left, #intDivStyle-right {
border-radius: 10px;
padding: 10px;
background-color: #fdfdfd; /* Add a background color to distinguish from the outer div */
}
/* Tab styling to match your theme */
.q-tabs {
background-color: transparent !important;
}
.q-tab {
color: #7f8389 !important;
font-weight: 500;
}
.q-tab--active {
color: #009f90 !important;
background-color: rgba(0, 159, 144, 0.1) !important;
border-radius: 5px;
}
.q-tabs[active-color="primary"] .q-tab--active {
color: #009f90 !important;
}
/* Indicator color to match your primary color */
.q-tab__indicator {
background-color: #009f90 !important;
}
/* Dense tabs adjustment */
.q-tabs--dense .q-tab {
padding: 8px 12px;
font-size: 0.9rem;
}
#intDivStyle-left .q-tab-panel, #intDivStyle-right .q-tab-panel {
height: auto; /* Let content define height */
overflow-y: hidden; /* Remove scrollbar */
}
.pixelated-plot svg image {
image-rendering: pixelated;
image-rendering: -moz-crisp-edges; /* Firefox */
image-rendering: crisp-edges;
}
*{
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;
}

1
public/css/masks/static.txt Normal file
View File

@ -0,0 +1 @@
masks

1
public/css/themes/static.txt Normal file
View File

@ -0,0 +1 @@
themes

3
routes.jl Normal file
View File

@ -0,0 +1,3 @@
# Include all your app modules
include("app.jl")
include("mask.jl")

255
src/BloomFilters.jl Normal file
View File

@ -0,0 +1,255 @@
# 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
"""
StreamingBloomFilter
A memory-efficient Bloom filter that doesn't store all values at once.
"""
mutable struct StreamingBloomFilter
bits::BitVector
size::Int
hash_count::Int
seed::UInt64
count::Int
# Streaming state
current_chunk::Vector{Int}
chunk_size::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`)
"""
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
"""
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
- `p::Float64`: Desired false positive rate
"""
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.
"""
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.
"""
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.
"""
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.
"""
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.
"""
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

217
src/Common.jl Normal file
View File

@ -0,0 +1,217 @@
# src/Common.jl - Updated with BloomFilter
using Base.Threads
# --- Buffer Pooling ---
"""
BufferPool
A thread-safe pool of `Vector{UInt8}` buffers, categorized by their size.
This helps reduce memory allocations and garbage collection overhead by reusing buffers.
# Fields
- `lock`: A `ReentrantLock` to ensure thread-safe access to the pool.
- `buffers`: A dictionary mapping buffer size (in bytes) to a list of available buffers of that size.
"""
mutable struct BufferPool
lock::ReentrantLock
buffers::Dict{Int, Vector{Vector{UInt8}}}
function BufferPool()
new(ReentrantLock(), Dict{Int, Vector{Vector{UInt8}}}())
end
end
"""
get_buffer!(pool::BufferPool, 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 `BufferPool` to retrieve the buffer from.
- `size`: The minimum desired size of the buffer in bytes.
# Returns
- A `Vector{UInt8}` buffer.
"""
function get_buffer!(pool::BufferPool, size::Int)::Vector{UInt8}
lock(pool.lock) do
# Find an existing buffer that is large enough
for (buffer_size, buffer_list) in pool.buffers
if buffer_size >= size && !isempty(buffer_list)
buffer = pop!(buffer_list)
# Resize if necessary (and if it's significantly larger than needed)
if length(buffer) > size * 2 # Heuristic: if buffer is more than twice the requested size
return Vector{UInt8}(undef, size) # Allocate new smaller buffer
else
return buffer
end
end
end
# No suitable buffer found, allocate a new one
return Vector{UInt8}(undef, size)
end
end
"""
release_buffer!(pool::BufferPool, buffer::Vector{UInt8})
Returns a `Vector{UInt8}` buffer to the pool for future reuse.
# Arguments
- `pool`: The `BufferPool` to return the buffer to.
- `buffer`: The `Vector{UInt8}` buffer to release.
"""
function release_buffer!(pool::BufferPool, buffer::Vector{UInt8})
lock(pool.lock) do
buffer_size = length(buffer)
if !haskey(pool.buffers, buffer_size)
pool.buffers[buffer_size] = Vector{Vector{UInt8}}()
end
push!(pool.buffers[buffer_size], buffer)
end
return nothing
end
"""
SimpleBufferPool
A dramatically simplified buffer pool that avoids complex locking.
"""
mutable struct SimpleBufferPool
buffers::Dict{Int, Vector{Vector{UInt8}}}
max_pool_size::Int
end
SimpleBufferPool() = SimpleBufferPool(Dict{Int, Vector{Vector{UInt8}}}(), 50)
function get_buffer!(pool::SimpleBufferPool, size::Int)::Vector{UInt8}
# Check for existing buffers of exact size first
if haskey(pool.buffers, size) && !isempty(pool.buffers[size])
return pop!(pool.buffers[size])
end
# No suitable buffer found, allocate new one
return Vector{UInt8}(undef, size)
end
function release_buffer!(pool::SimpleBufferPool, buffer::Vector{UInt8})
size = length(buffer)
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
# 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
const DEFAULT_NUM_BINS = 2000
const DEFAULT_BLOOM_FILTER_SIZE = 10000 # Default size for empty spectra
const DEFAULT_FALSE_POSITIVE_RATE = 0.01
"""
validate_spectrum_data(mz, intensity, id)
Checks a spectrum for common data integrity issues.
"""
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.
"""
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.
"""
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.
"""
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

109
src/ImageProcessing.jl Normal file
View File

@ -0,0 +1,109 @@
# src/ImageProcessing
using Images
using ImageBinarization
using ImageMorphology
using ImageComponentAnalysis
using Colors # For converting to grayscale
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
# ===================================================================
"""
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;
otsu_scale=1.0,
noise_size_percent=0.1,
hole_size_percent=0.05,
smoothing=2)
# --- Step 1: Otsu Binarization ---
otsu_threshold = find_threshold(gray_img, Otsu())
adjusted_threshold = otsu_threshold * otsu_scale
binary_img = gray_img .>= adjusted_threshold
# --- Smart Parameter Scaling ---
image_area = length(gray_img)
noise_size_pixels = round(Int, image_area * noise_size_percent)
hole_size_pixels = round(Int, image_area * hole_size_percent)
# --- Step 2: Remove Small White Regions (Noise) ---
# area_opening is the correct morphological operation for this.
noise_removed_img = area_opening(binary_img, min_area=noise_size_pixels)
# --- Step 3: Fill Small Black Holes ---
# area_closing is the dual of area_opening and fills holes.
holes_filled_img = area_closing(noise_removed_img, min_area=hole_size_pixels)
# --- Step 4: Smooth Edges ---
# A morphological closing with a small disk smooths outlines.
smoothing_kernel = ones(Bool, (smoothing, smoothing))
smoothed_img = closing(holes_filled_img, smoothing_kernel)
# --- Return all intermediate steps for visualization ---
return binary_img, noise_removed_img, holes_filled_img, smoothed_img
end

1550
src/MSIData.jl Normal file
View File

File diff suppressed because it is too large Load Diff

108
src/MSI_src.jl Normal file
View File

@ -0,0 +1,108 @@
# MSI_src.jl
module MSI_src
# Export the public MSI API
export OpenMSIData,
GetSpectrum,
IterateSpectra,
ImportMzmlFile,
plot_slice,
get_total_spectrum,
get_average_spectrum,
precompute_analytics,
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
# Export the public preprocessing & precalculations API
export run_preprocessing_analysis,
main_precalculation,
FeatureMatrix,
Calibration,
Smoothing,
BaselineCorrection,
Normalization,
PeakPicking,
PeakBinning,
get_masked_spectrum_indices,
detect_peaks_profile_core,
detect_peaks_centroid_core,
smooth_spectrum_core,
apply_baseline_correction_core,
apply_normalization_core,
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
# Include all source files directly into the main module
include("BloomFilters.jl")
include("Common.jl")
include("MSIData.jl")
include("ParserHelpers.jl")
include("mzML.jl")
include("imzML.jl")
include("MzmlConverter.jl")
include("Preprocessing.jl")
include("ImageProcessing.jl")
include("Precalculations.jl")
include("PreprocessingPipeline.jl")
using Setfield # For immutable struct updates
# --- Main Entry Point --- #
"""
OpenMSIData(filepath::String; cache_size=100)
Opens a .mzML or .imzML file and prepares it for data access.
This is the main entry point for the new data access API.
"""
function OpenMSIData(filepath::String; cache_size=300, spectrum_type_map::Union{Dict{Int, Symbol}, Nothing}=nothing)
local msi_data
if endswith(lowercase(filepath), ".mzml")
msi_data = load_mzml_lazy(filepath, cache_size=cache_size)
elseif endswith(lowercase(filepath), ".imzml")
msi_data = load_imzml_lazy(filepath, cache_size=cache_size)
else
error("Unsupported file type: $filepath. Please provide a .mzML or .imzML file.")
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 # module MSI_src

873
src/MzmlConverter.jl Normal file
View File

@ -0,0 +1,873 @@
# src/MzmlConverter.jl
"""
This file contains the workflow for converting .mzML files (with one spectrum per pixel)
into a proper .imzML/.ibd file pair, using a separate synchronization file.
It replicates the functionality of the original R scripts that use MALDIquant.
"""
using DataFrames, Printf, CSV, UUIDs, ProgressMeter
# 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
# A struct to hold the processed pixel data before exporting
"""
ProcessedPixel
A temporary struct to hold the data for a single, fully rendered pixel before
it is written to the final `.ibd` file.
# Fields
- `coords`: A tuple `(x, y)` of the pixel's spatial coordinates.
- `mz`: The m/z array for the pixel.
- `intensity`: The calculated intensity array for the pixel.
"""
struct ProcessedPixel
coords::Tuple{Int, Int}
mz::Vector{Float64} # Assuming m/z is consistent, can be optimized later
intensity::Vector{Float32}
end
"""
BinaryMetadata
A struct to hold the file offset and length for a spectrum's binary data arrays
(m/z and intensity) after they have been written to the `.ibd` file.
# Fields
- `mz_offset`, `mz_length`: Byte offset and length for the m/z array.
- `int_offset`, `int_length`: Byte offset and length for the intensity array.
"""
struct BinaryMetadata
mz_offset::UInt64
mz_length::UInt64
int_offset::UInt64
int_length::UInt64
end
"""
GetMzmlScanTime_linebyline(fileName::String)
Parses a `.mzML` file line-by-line to extract the scan start time for each
spectrum. This is a slow and memory-intensive fallback method used only when
the faster, index-based `GetMzmlScanTime` fails.
# Arguments
- `fileName`: Path to the `.mzML` file.
# Returns
- A `Matrix{Int64}` where each row is `[spectrum_index, time_in_milliseconds]`.
"""
function GetMzmlScanTime_linebyline(fileName::String)
times = Tuple{Int64, Int64}[]
try
file_size = filesize(fileName)
estimated_spectra = max(1000, file_size ÷ 10000)
sizehint!(times, estimated_spectra)
catch
end
open(fileName, "r") do stream
state = :outside_spectrum
current_index = 0
current_time = nothing
for line in eachline(stream)
if state == :outside_spectrum
if occursin("<spectrum ", line)
state = :in_spectrum
current_index = 0
current_time = nothing
idx_match = match(r"index=\"(\d+)\"", line)
if idx_match !== nothing
current_index = parse(Int, idx_match.captures[1]) + 1
end
end
elseif state == :in_spectrum
if occursin("<cvParam", line)
if occursin("MS:1000016", line) || occursin("MS:1000015", line)
time_match = match(r"value=\"([\d\.]+)\"", line)
if time_match !== nothing
time_val = parse(Float64, time_match.captures[1])
unit_scale = occursin("MS:1000016", line) ? 60000 : 1000
current_time = round(Int64, time_val * unit_scale)
end
end
end
if occursin("</spectrum>", line)
state = :outside_spectrum
if current_index > 0 && current_time !== nothing
push!(times, (current_index, current_time))
end
end
end
end
end
result = Matrix{Int64}(undef, length(times), 2)
for (i, t) in enumerate(times)
result[i, 1] = t[1]
result[i, 2] = t[2]
end
return result
end
"""
GetMzmlScanTime(fileName::String)
Parses a .mzML file to extract the scan start time for each spectrum.
Uses the indexed part of the .mzML file for fast access, falling back to
a slower line-by-line parse if the index is not present.
# Arguments
* `fileName`: Path to the .mzML file.
# Returns
- A `Matrix{Int64}` where each row is `[spectrum_index, time_in_milliseconds]`.
"""
function GetMzmlScanTime(fileName::String)
times = Tuple{Int64, Int64}[]
try
open(fileName, "r") do stream
# 1. Find and parse the spectrum index offsets
seekend(stream)
end_chunk_size = min(filesize(stream), 8192)
seek(stream, filesize(stream) - end_chunk_size)
footer = read(stream, String)
index_offset_match = match(r"<indexListOffset>(\d+)</indexListOffset>", footer)
if index_offset_match === nothing
@warn "No <indexListOffset> found. Falling back to slow line-by-line parsing for scan times. This may be memory intensive."
return GetMzmlScanTime_linebyline(fileName)
end
index_offset = parse(Int64, index_offset_match.captures[1])
seek(stream, index_offset)
# The find_tag function is defined in ParserHelpers.jl
if find_tag(stream, r"<index\s+name=\"spectrum\"") === nothing
@warn "Could not find spectrum index. Falling back to slow line-by-line parsing."
return GetMzmlScanTime_linebyline(fileName)
end
# The parse_offset_list function is defined in mzML.jl
spectrum_offsets = parse_offset_list(stream)
# 2. Iterate through offsets and parse time for each spectrum
for (idx, offset) in enumerate(spectrum_offsets)
seek(stream, offset)
current_time = nothing
# Read a limited number of lines to find the time
for _ in 1:100
if eof(stream) break end
line = readline(stream)
if occursin("MS:1000016", line) || occursin("MS:1000015", line) # scan start time
time_match = match(r"value=\"([\d\.]+)\"", line)
if time_match !== nothing
time_val = parse(Float64, time_match.captures[1])
unit_scale = occursin("MS:1000016", line) ? 60000 : 1000 # minutes vs seconds
current_time = round(Int64, time_val * unit_scale)
break
end
end
if occursin("<binary>", line) || occursin("</spectrum>", line)
break
end
end
if current_time !== nothing
push!(times, (idx, current_time))
end
end
end
catch e
@error "Failed to parse scan times with indexed method. Falling back to line-by-line." exception=(e, catch_backtrace())
return GetMzmlScanTime_linebyline(fileName)
end
# The index is already sorted by spectrum index, so no need to sort `times`.
# Normalize times relative to the first scan
if !isempty(times)
first_time = times[1][2]
for i in eachindex(times)
times[i] = (times[i][1], times[i][2] - first_time)
end
end
result = Matrix{Int64}(undef, length(times), 2)
for (i, t) in enumerate(times)
result[i, 1] = t[1]
result[i, 2] = t[2]
end
return result
end
"""
MatchAcquireTime(sync_file_path::String, scans::Matrix{Int64})
Correlates pixel acquisition times from a synchronization file with scan
acquisition times from an mzML file.
This is a Julia implementation of the `MatchAcquireTime` function from the R scripts.
# Arguments
* `sync_file_path`: Path to the synchronization file (.txt).
* `scans`: A matrix of scan times, as returned by `GetMzmlScanTime`.
# Returns
- A `Matrix{Int64}` where each row is `[pixel_index, pixel_time_ms, first_scan_index, last_scan_index]`.
"""
function MatchAcquireTime(sync_file_path::String, scans::Matrix{Int64}; img_width::Int=0, img_height::Int=0)
if !isfile(sync_file_path)
error("Synchronization file not found: $sync_file_path")
end
pixel_df = CSV.read(sync_file_path, DataFrame, header=false, skipto=3)
pixel_matrix = Matrix(pixel_df)
num_pixels_original = size(pixel_matrix, 1)
if num_pixels_original == 0
return zeros(Int64, 0, 5)
end
# Determine image dimensions and generate coordinates
current_img_width = img_width
current_img_height = img_height
if current_img_width == 0 || current_img_height == 0
if size(pixel_matrix, 2) >= 3
@info "Sync file contains X, Y coordinates"
coordinates = convert(Matrix{Int}, pixel_matrix[:, 1:2])
pixel_times = convert(Vector{Int64}, pixel_matrix[:, 3])
current_img_width = maximum(coordinates[:, 1])
current_img_height = maximum(coordinates[:, 2])
# Apply R's pixel truncation logic
num_pixels = (floor(Int, (num_pixels_original - 1) / current_img_width)) * current_img_width
if num_pixels <= 0
error("Calculated num_pixels to process is zero or negative: $num_pixels")
end
# Truncate arrays
final_coordinates = coordinates[1:num_pixels, :]
final_pixel_times = pixel_times[1:num_pixels]
else
@info "Sync file contains Index and Time; generating coordinates"
pixel_times = convert(Vector{Int64}, pixel_matrix[:, 2])
# R's width detection logic
diffs = pixel_matrix[2:end, 1] .- pixel_matrix[1:end-1, 1]
width_indices = findall(x -> x != 1, diffs)
if !isempty(width_indices)
current_img_width = width_indices[1]
else
current_img_width = num_pixels_original
end
current_img_height = num_pixels_original ÷ current_img_width
# R's pixel truncation logic
num_pixels = (floor(Int, (num_pixels_original - 1) / current_img_width)) * current_img_width
if num_pixels <= 0
error("Calculated num_pixels to process is zero or negative: $num_pixels")
end
# Generate coordinates for truncated pixels
final_coordinates = zeros(Int, num_pixels, 2)
final_pixel_times = zeros(Int64, num_pixels)
for i in 1:num_pixels
idx = pixel_matrix[i, 1]
final_coordinates[i, 1] = ((idx - 1) % current_img_width) + 1
final_coordinates[i, 2] = fld(idx - 1, current_img_width) + 1
final_pixel_times[i] = pixel_times[i]
end
end
else
@info "Using provided dimensions: $(current_img_width)x$(current_img_height)"
pixel_times = convert(Vector{Int64}, pixel_matrix[:, 2])
# R's pixel truncation logic (FIXED: consistent formula)
num_pixels = (floor(Int, (num_pixels_original - 1) / current_img_width)) * current_img_width
if num_pixels <= 0
error("Calculated num_pixels to process is zero or negative: $num_pixels")
end
# Generate coordinates
final_coordinates = zeros(Int, num_pixels, 2)
final_pixel_times = zeros(Int64, num_pixels)
for i in 1:num_pixels
idx = pixel_matrix[i, 1]
final_coordinates[i, 1] = ((idx - 1) % current_img_width) + 1
final_coordinates[i, 2] = fld(idx - 1, current_img_width) + 1
final_pixel_times[i] = pixel_times[i]
end
end
# Normalize pixel times
if !isempty(final_pixel_times)
min_pixel_time = minimum(final_pixel_times)
final_pixel_times .-= min_pixel_time
end
num_scans = size(scans, 1)
if num_scans == 0
return hcat(final_coordinates, final_pixel_times, zeros(Int64, num_pixels, 2))
end
# FIXED R's time matching algorithm with bounds checking
first_idx = 1
last_idx = 1
index_matrix = zeros(Int, num_pixels, 2)
for i_pixel in 1:num_pixels
pixel_time = final_pixel_times[i_pixel]
# Find the first scan that reaches or exceeds pixel time
while last_idx <= num_scans && scans[last_idx, 2] < pixel_time
last_idx += 1
end
# Ensure valid indices
if last_idx > num_scans
# No more scans available for remaining pixels
index_matrix[i_pixel:end, 1] .= num_scans + 1 # Invalid index
index_matrix[i_pixel:end, 2] .= num_scans # Invalid index
break
end
# Assign scan indices with bounds checking
start_scan = max(1, first_idx)
end_scan = max(1, last_idx - 1)
# Ensure start_scan <= end_scan
if start_scan > end_scan
start_scan = end_scan
end
index_matrix[i_pixel, 1] = start_scan
index_matrix[i_pixel, 2] = end_scan
# Update for next iteration (R's algorithm)
first_idx = last_idx - 1
last_idx = first_idx
# Ensure first_idx doesn't go below 1
if first_idx < 1
first_idx = 1
last_idx = 1
end
end
return hcat(final_coordinates, final_pixel_times, index_matrix)
end
"""
RenderPixel(pixel_info, scans, msi_data, scan_time_deltas, pixel_time_deltas)
Reconstructs the spectrum for a single pixel by combining intensities from the
raw MS scans that occurred during the pixel's acquisition time. It uses a
weighted interpolation scheme based on the relative timing of scans and pixels.
This function handles three cases:
1. A single scan falls entirely within the pixel's time window.
2. The pixel's time window is covered by two partial scans.
3. The pixel's time window covers one or more full scans plus two partial scans.
# Arguments
- `pixel_info`: A row from the timing matrix containing the pixel's time and scan indices.
- `scans`: The matrix of scan times.
- `msi_data`: The `MSIData` object for the source `.mzML` file.
- `scan_time_deltas`: Pre-calculated time durations for each scan.
- `pixel_time_deltas`: Pre-calculated time durations for each pixel.
# Returns
- A tuple `(mz_array, intensity_array)` for the rendered pixel spectrum.
"""
function RenderPixel(
intensity_buffer::Vector{Float32},
pixel_info::AbstractVector{Int64},
scans::AbstractMatrix{Int64},
msi_data::MSIData,
scan_time_deltas::AbstractVector{Int64},
pixel_time_deltas::AbstractVector{Int64}
)
pixel_time = pixel_info[3]
first_scan = pixel_info[4]
last_scan = pixel_info[5]
num_actions = last_scan - first_scan
# Get reference m/z array from the first scan involved.
mz_array, _ = GetSpectrum(msi_data, first_scan)
# If the first spectrum was empty, we can't do anything else.
if isempty(mz_array)
return (mz_array, view(intensity_buffer, 0:-1), UNKNOWN)
end
# Determine the mode for the output spectrum. If any contributing scan is profile, the result is profile.
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
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
end
if num_actions == 0 # Single scan contributes to the pixel
process_spectrum(msi_data, first_scan) do _, intensity
# SAFETY: Ensure positive scaling
scale = max(pixel_time_deltas[first_scan] / scan_time_deltas[first_scan], 0.0f0)
new_intensity .= intensity .* scale
end
elseif num_actions == 1 # Two partial scans contribute
# First partial scan
process_spectrum(msi_data, first_scan) do _, intensity1
scale1 = max((scans[last_scan, 2] - pixel_time) / scan_time_deltas[first_scan], 0.0f0)
new_intensity .+= intensity1 .* scale1
end
# Second partial scan
process_spectrum(msi_data, last_scan) do _, intensity2
next_pixel_time = pixel_time + pixel_time_deltas[first_scan]
scale2 = max((next_pixel_time - scans[last_scan, 2]) / scan_time_deltas[last_scan], 0.0f0)
new_intensity .+= intensity2 .* scale2
end
elseif num_actions > 1 # Multiple scans contribute
# First partial scan
process_spectrum(msi_data, first_scan) do _, intensity1
scale1 = max((scans[first_scan + 1, 2] - pixel_time) / scan_time_deltas[first_scan], 0.0f0)
new_intensity .+= intensity1 .* scale1
end
# Full scans in the middle
for i in (first_scan + 1):(last_scan - 1)
process_spectrum(msi_data, i) do _, intensity_middle
new_intensity .+= intensity_middle
end
end
# Last partial scan
process_spectrum(msi_data, last_scan) do _, intensity2
next_pixel_time = pixel_time + pixel_time_deltas[first_scan]
scale2 = max((next_pixel_time - scans[last_scan, 2]) / scan_time_deltas[last_scan], 0.0f0)
new_intensity .+= intensity2 .* scale2
end
end
# FINAL SAFETY: Clamp any negative values to zero
new_intensity = max.(new_intensity, 0.0f0)
return (mz_array, new_intensity, final_mode)
end
"""
ConvertMzmlToImzml(source_file, target_ibd_file, timing_matrix, scans)
Orchestrates the conversion of spectra from a `.mzML` file into a binary `.ibd`
file. It iterates through each pixel defined in the `timing_matrix`, calls
`RenderPixel` to reconstruct the pixel's spectrum, and writes the resulting
m/z and intensity arrays to the `.ibd` file in little-endian byte order.
# Arguments
- `source_file`: Path to the source `.mzML` file.
- `target_ibd_file`: Path for the output `.ibd` binary file.
- `timing_matrix`: The output from `MatchAcquireTime`, mapping pixels to scans.
- `scans`: The matrix of scan times from `GetMzmlScanTime`.
# Returns
- A tuple `(binary_meta_vec, coords_vec, (width, height))` containing:
- A vector of `BinaryMetadata` for each spectrum.
- A vector of `(x, y)` coordinate tuples.
- A tuple of the final image dimensions.
"""
function ConvertMzmlToImzml(source_file::String, target_ibd_file::String, timing_matrix::Matrix{Int64}, scans::Matrix{Int64})
if size(timing_matrix, 1) == 0
# Create an empty .ibd file if there's nothing to process
open(target_ibd_file, "w") do ibd_stream
write(ibd_stream, zeros(UInt8, 16)) # UUID placeholder
end
return BinaryMetadata[], Tuple{Int, Int}[], (0, 0), SpectrumMode[], uuid4()
end
width = maximum(timing_matrix[:, 1])
height = maximum(timing_matrix[:, 2])
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))
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))
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
return binary_meta_vec, coords_vec, (width, height), pixel_modes, ibd_uuid, msi_data.instrument_metadata
end
"""
ExportImzml(target_file, binary_meta, coords, dims)
Generates the `.imzML` metadata file. This XML file contains all the necessary
metadata to interpret the corresponding `.ibd` binary file, including references
to external data offsets, image dimensions, and CV parameters describing the
experiment and data format.
# Arguments
- `target_file`: The path for the output `.imzML` file.
- `binary_meta`: A vector of `BinaryMetadata` structs with offset and length info.
- `coords`: A vector of `(x, y)` coordinates for each spectrum.
- `dims`: A tuple `(width, height)` of the final image dimensions.
# Returns
- `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)
ibd_file = replace(target_file, r"\.imzML$"i => ".ibd")
if isempty(binary_meta)
@warn "No binary metadata to export; creating empty imzML file."
# Still create a valid, empty imzML file
end
try
# The .ibd file is now written by ConvertMzmlToImzml.
# This function is only responsible for the .imzML XML metadata file.
open(target_file, "w") do imzml_stream
# XML Header
write(imzml_stream, """<?xml version="1.0" encoding="ISO-8859-1"?>
<mzML xmlns="http://psi.hupo.org/ms/mzml" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://psi.hupo.org/ms/mzml http://psidev.info/files/ms/mzML/xsd/mzML1.1.0.xsd" version="1.1" id="$(splitext(basename(target_file))[1])">
""")
# CV List, File Description, etc. (static parts)
write(imzml_stream, """ <cvList count="3">
<cv id="MS" fullName="Proteomics Standards Initiative Mass Spectrometry Ontology" version="1.3.1" URI="http://psidev.info/ms/mzML/psi-ms.obo"/>
<cv id="UO" fullName="Unit Ontology" version="1.15" URI="http://obo.cvs.sourceforge.net/obo/obo/ontology/phenotype/unit.obo"/>
<cv id="IMS" fullName="Imaging MS Ontology" version="0.9.1" URI="http://www.maldi-msi.org/download/imzml/imagingMS.obo"/>
</cvList>
""")
write(imzml_stream, """ <fileDescription>
<fileContent>
<cvParam cvRef="MS" accession="MS:1000579" name="MS1 spectrum"/>
<cvParam cvRef="IMS" accession="IMS:1000080" name="mass spectrum"/>
<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>
</fileDescription>
""")
write(imzml_stream, """ <referenceableParamGroupList count="2">
<referenceableParamGroup id="mzArray">
<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:1000523" name="64-bit float"/>
</referenceableParamGroup>
<referenceableParamGroup id="intensityArray">
<cvParam cvRef="MS" accession="MS:1000576" name="no compression"/>
<cvParam cvRef="MS" accession="MS:1000515" name="intensity array" unitCvRef="MS" unitAccession="MS:1000131" unitName="number of counts"/>
<cvParam cvRef="MS" accession="MS:1000521" name="32-bit float"/>
</referenceableParamGroup>
</referenceableParamGroupList>
""")
write(imzml_stream, """ <sampleList count="1">
<sample id="sample1" name="ImagingSample">
<cvParam cvRef="MS" accession="MS:1000001" name="sample number" value="1"/>
</sample>
</sampleList>
""")
write(imzml_stream, """ <softwareList count="1">
<software id="MSIConverter" version="1.0">
<cvParam cvRef="MS" accession="MS:1000799" name="custom unreleased software tool" value="mzML to imzML converter"/>
</software>
</softwareList>
""")
write(imzml_stream, """ <scanSettingsList count="1">
<scanSettings id="scanSettings1">
<cvParam cvRef="IMS" accession="IMS:1000042" name="max count of pixel x" value="$(dims[1])"/>
<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: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>
</scanSettingsList>
""")
write(imzml_stream, """ <instrumentConfigurationList count="1">
<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">
<source order="1">
<cvParam cvRef="MS" accession="MS:1000075" name="MALDI source"/>
</source>
<analyzer order="2">
<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>
<detector order="3">
<cvParam cvRef="MS" accession="MS:1000253" name="electron multiplier"/>
</detector>
</componentList>
<softwareRef ref="MSIConverter"/>
</instrumentConfiguration>
</instrumentConfigurationList>
""")
write(imzml_stream, """ <dataProcessingList count="$(1 + (instrument_meta.vendor_preprocessing_steps !== nothing ? length(instrument_meta.vendor_preprocessing_steps) : 0))">
<dataProcessing id="conversionProcessing">
<processingMethod order="1" softwareRef="MSIConverter">
<cvParam cvRef="MS" accession="MS:1000544" name="Conversion to imzML"/>
</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>
</dataProcessingList>
""")
# Run and Spectrum List
spectrum_offsets = UInt64[]
write(imzml_stream, """ <run defaultInstrumentConfigurationRef="instrument1" id="run1" sampleRef="sample1">
<spectrumList count="$(length(binary_meta))" defaultDataProcessingRef="conversionProcessing">
""")
# Write each spectrum's metadata
for (i, meta) in enumerate(binary_meta)
x, y = coords[i]
spectrum_start = position(imzml_stream)
push!(spectrum_offsets, spectrum_start)
# Calculate number of points from byte length
mz_points = meta.mz_length ÷ sizeof(Float64)
int_points = meta.int_length ÷ sizeof(Float32)
write(imzml_stream, """ <spectrum id="Scan=$(i)" defaultArrayLength="$(mz_points)" index="$(i-1)">
<cvParam cvRef="MS" accession="MS:1000511" name="ms level" value="1"/>
""")
if modes[i] == CENTROID
write(imzml_stream, """ <cvParam cvRef="MS" accession="MS:1000127" name="centroid spectrum"/>
""")
elseif modes[i] == PROFILE
write(imzml_stream, """ <cvParam cvRef="MS" accession="MS:1000128" name="profile spectrum"/>
""")
end
write(imzml_stream, """ <scanList count="1">
<scan instrumentConfigurationRef="instrument1">
<cvParam cvRef="IMS" accession="IMS:1000050" name="position x" value="$(x)"/>
<cvParam cvRef="IMS" accession="IMS:1000051" name="position y" value="$(y)"/>
</scan>
</scanList>
<binaryDataArrayList count="2">
<binaryDataArray encodedLength="0">
<referenceableParamGroupRef ref="mzArray"/>
<cvParam cvRef="IMS" accession="IMS:1000103" name="external array length" value="$(mz_points)"/>
<cvParam cvRef="IMS" accession="IMS:1000102" name="external offset" value="$(meta.mz_offset)"/>
<binary/>
</binaryDataArray>
<binaryDataArray encodedLength="0">
<referenceableParamGroupRef ref="intensityArray"/>
<cvParam cvRef="IMS" accession="IMS:1000103" name="external array length" value="$(int_points)"/>
<cvParam cvRef="IMS" accession="IMS:1000102" name="external offset" value="$(meta.int_offset)"/>
<binary/>
</binaryDataArray>
</binaryDataArrayList>
</spectrum>
""")
end
write(imzml_stream, """ </spectrumList>
</run>
</mzML>
""")
end
println("Successfully created: $target_file")
println("Successfully created: $ibd_file")
return true
catch e
@error "Failed to export imzML metadata file" exception=(e, catch_backtrace())
# Clean up partial .imzML file
isfile(target_file) && rm(target_file, force=true)
# Do not delete the .ibd file as it might be useful for debugging
return false
end
end
"""
ImportMzmlFile(source_file::String, sync_file::String, target_file::String)
Main workflow function to convert a .mzML file to an .imzML file.
# Arguments
* `source_file`: Path to the input .mzML file.
* `sync_file`: Path to the synchronization text file.
* `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_height`: height dimention for the creation of the y axis
"""
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...")
scans = GetMzmlScanTime(source_file)
println("Step 2: Matching acquisition times...")
timing_matrix = MatchAcquireTime(sync_file, scans; img_width=img_width, img_height=img_height)
println("Step 3: Converting spectra and writing .ibd file...")
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)
# Flip image vertically to match R script output
flipped_coords = [(x, height - y + 1) for (x, y) in coords]
println("Step 4: Exporting .imzML metadata file...")
success = ExportImzml(target_file, binary_meta, flipped_coords, (width, height), pixel_modes, ibd_uuid, instrument_meta)
if success
println("Conversion successful: $target_file")
else
println("Conversion failed.")
end
return success
end

135
src/ParserHelpers.jl Normal file
View File

@ -0,0 +1,135 @@
# src/ParserHelpers.jl
"""
This file provides common helper functions for parsing mzML and imzML files.
"""
# ============================================================================
#
# Data Structures and Low-Level Parser Helpers
#
# ============================================================================
"""
SpecDim
A struct to hold configuration for a spectral data axis (e.g., m/z or intensity).
"""
mutable struct SpecDim
Format::Type
Packed::Bool
Axis::Int
Skip::Int
Mode::SpectrumMode
end
"""
find_tag(stream, regex::Regex)
Reads a stream line-by-line until a line matches the provided regex.
# Returns
- A `RegexMatch` object if a match is found, otherwise throws an error.
"""
function find_tag(stream, regex::Regex)
while !eof(stream)
line = readline(stream)
isTag = match(regex, line)
if isTag !== nothing
return isTag
end
end
error("Tag not found for regex: $regex")
end
"""
get_attribute(source::String, tag::String = "([^=]+)")
Retrieves an attribute's value from an XML tag string.
# Returns
- A `RegexMatch` object containing the attribute and its value.
"""
function get_attribute(source::AbstractString, tag::String = "([^=]+)")
# Construct the regex pattern string
pattern_str = "\\s*" * tag * "=\"([^\"]*)\""
regStr = Regex(pattern_str)
return match(regStr, source)
end
"""
configure_spec_dim(stream)
Parses a block of `<cvParam>` tags from an XML stream to configure a `SpecDim`
struct. It reads accessions to determine the data format (e.g., `Float32`),
compression status (`zlib`), and axis type (m/z vs. intensity).
# Arguments
- `stream`: An IO stream positioned at the start of the `cvParam` block.
# Returns
- A `SpecDim` struct populated with the parsed configuration.
"""
function configure_spec_dim(stream)
axis = SpecDim(Float64, false, 1, 0, UNKNOWN) # Add UNKNOWN as default mode
while !eof(stream)
currLine = readline(stream)
matchInfo = match(r"^\s*<(cvParam)", currLine)
if matchInfo === nothing
# Check for end of cvParam block
if match(r"^\s*</", currLine) !== nothing
return axis
end
continue
end
index = length(matchInfo.captures[1])
attr_match = get_attribute(currLine[index:end], "accession")
if attr_match !== nothing
accession = attr_match.captures[1]
if accession == "MS:1000515" # intensity array
axis.Axis = 2
elseif accession == "MS:1000519" # 32-bit integer
axis.Format = Int32
elseif accession == "MS:1000521" # 32-bit float
axis.Format = Float32
elseif accession == "MS:1000522" # 64-bit integer
axis.Format = Int64
elseif accession == "MS:1000574" # zlib compression
axis.Packed = true
elseif accession == "MS:1000127" # centroid spectrum
axis.Mode = CENTROID
elseif accession == "MS:1000128" # profile spectrum
axis.Mode = PROFILE
end
end
end
return axis
end
# ============================================================================
#
# Data Structures and Helpers for mzML lazy parsing
#
# ============================================================================
"""
CVParams
A mutable struct used during the parsing of mzML files to temporarily hold
CV parameter information for a data array before it is converted into a
`SpectrumAsset`.
# Fields
- `format`: The data type (e.g., `Float64`).
- `is_compressed`: A boolean indicating if the data is compressed.
- `axis_type`: A symbol indicating the array type (`:mz` or `:intensity`).
"""
mutable struct CVParams
format::Type
is_compressed::Bool
axis_type::Symbol
end

1759
src/Precalculations.jl Normal file
View File

File diff suppressed because it is too large Load Diff

1246
src/Preprocessing.jl Normal file
View File

File diff suppressed because it is too large Load Diff

456
src/PreprocessingPipeline.jl Normal file
View File

@ -0,0 +1,456 @@
# 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)
baseline = apply_baseline_correction_core(s.intensity; method=method, iterations=iterations, window=window)
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)
s.intensity = transform_intensity_core(s.intensity; method=method)
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)
smoothed_intensity = max.(0.0, smooth_spectrum_core(s.intensity; method=method, window=window, order=order))
s.intensity = 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)
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
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)
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())
s.mz = itp(s.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)
s.mz = alignment_func.(s.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)
s.intensity = apply_normalization_core(s.intensity; method=method)
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

1536
src/imzML.jl Normal file
View File

File diff suppressed because it is too large Load Diff

462
src/mzML.jl Normal file
View File

@ -0,0 +1,462 @@
# src/mzML.jl
# This file is responsible for parsing metadata from .mzML files.
# It has been refactored to produce a unified MSIData object.
# Constants for CV parameter accessions - defined once for performance
const MZ_AXIS_ACCESSION = "MS:1000514"
const INTENSITY_AXIS_ACCESSION = "MS:1000515"
const COMPRESSION_ACCESSION = "MS:1000574"
const NO_COMPRESSION_ACCESSION = "MS:1000576"
# Data format accessions as constants
const DATA_FORMAT_ACCESSIONS = Dict{String, DataType}(
"MS:1000518" => Int16,
"MS:1000519" => Int32,
"MS:1000520" => Float64,
"MS:1000521" => Float32,
"MS:1000522" => Int64,
"MS:1000523" => Float64
)
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)
Parses a `<binaryDataArray>` block within an mzML file to extract metadata
for a single data array (e.g., m/z or intensity). It reads CV parameters to
determine the data type, compression, and axis type.
# Arguments
- `stream`: An IO stream positioned at the beginning of a `<binaryDataArray>` block.
# Returns
- A `SpectrumAsset` struct containing the parsed metadata, including the binary
data offset, encoded length, format, and compression status.
"""
function get_spectrum_asset_metadata(stream::IO)
start_pos = position(stream)
#println("DEBUG: Entering get_spectrum_asset_metadata to parse binaryDataArray...")
bda_tag = find_tag(stream, r"<binaryDataArray\s+encodedLength=\"(\d+)\"")
if bda_tag === nothing
throw(FileFormatError("Cannot find binaryDataArray"))
end
encoded_length = parse(Int32, bda_tag.captures[1])
#println("DEBUG: Encoded length: $encoded_length")
# Initialize parameters as separate variables with concrete types
data_format::DataType = Float64
compression_flag::Bool = false
axis::Symbol = :mz
while !eof(stream)
line = readline(stream)
if occursin("</binaryDataArray>", line)
break
end
if occursin("<cvParam", line)
accession = get_attribute(line, "accession")
if accession === nothing
continue
end
acc_str = accession.captures[1]
# Use constant comparisons and dictionary lookup for better performance
if acc_str == MZ_AXIS_ACCESSION
axis = :mz
#println("DEBUG: Axis type identified as: m/z")
elseif acc_str == INTENSITY_AXIS_ACCESSION
axis = :intensity
#println("DEBUG: Axis type identified as: intensity")
elseif haskey(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
compression_flag = true
#println("DEBUG: Compression: true")
elseif acc_str == NO_COMPRESSION_ACCESSION
compression_flag = false
#println("DEBUG: Compression: false")
end
end
end
seek(stream, start_pos)
readuntil(stream, "<binary>")
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
readuntil(stream, "</binaryDataArray>")
#println("DEBUG: Exiting get_spectrum_asset_metadata.")
# Create SpectrumAsset directly from the variables
return SpectrumAsset(data_format, compression_flag, binary_offset, encoded_length, axis)
end
# This function is updated to return the generic SpectrumMetadata struct
"""
parse_spectrum_metadata(stream::IO, offset::Int64)
Parses an entire `<spectrum>` block from an mzML file, given a starting offset.
It extracts the spectrum ID and calls `get_spectrum_asset_metadata` to parse
the m/z and intensity array metadata.
# Arguments
- `stream`: An IO stream for the mzML file.
- `offset`: The byte offset where the `<spectrum>` block begins.
# Returns
- A `SpectrumMetadata` struct containing the parsed metadata for one spectrum.
"""
function parse_spectrum_metadata(stream::IO, offset::Int64)
seek(stream, offset)
# Read the whole spectrum block to parse mode
spectrum_start_pos = position(stream)
line = ""
spectrum_buffer = IOBuffer()
while !eof(stream)
line = readline(stream)
write(spectrum_buffer, line)
if occursin("</spectrum>", line)
break
end
end
spectrum_xml = String(take!(spectrum_buffer))
seek(stream, spectrum_start_pos) # Reset for other parsing
id_match = match(r"<spectrum\s+index=\"\d+\"\s+id=\"([^\"]+)", spectrum_xml)
id = id_match === nothing ? "" : id_match.captures[1]
#println("DEBUG: Parsing spectrum ID: $id")
# Determine mode from the XML block
mode = UNKNOWN
if occursin("MS:1000127", spectrum_xml)
mode = CENTROID
#println("DEBUG: Spectrum mode: CENTROID")
elseif occursin("MS:1000128", spectrum_xml)
mode = PROFILE
#println("DEBUG: Spectrum mode: PROFILE")
end
# Find where the binary data list starts to parse assets
binary_list_match = findfirst("<binaryDataArrayList", spectrum_xml)
if binary_list_match !== nothing
seek(stream, spectrum_start_pos + binary_list_match.start - 1)
end
asset1 = get_spectrum_asset_metadata(stream)
asset2 = get_spectrum_asset_metadata(stream)
# Determine which asset is m/z and which is intensity
mz_asset, int_asset = if asset1.axis_type == :mz
(asset1, asset2)
else
(asset2, asset1)
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
# 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)
end
"""
parse_offset_list(stream::IO)
Parses the `<index name="spectrum">` block in an indexed mzML file to extract
the byte offsets for each spectrum.
# Arguments
- `stream`: An IO stream positioned at the start of the `<index>` block.
# Returns
- A `Vector{Int64}` containing the byte offsets for all spectra.
"""
function parse_offset_list(stream::IO)
offsets = Int64[]
offset_regex = r"<offset[^>]*>(\d+)</offset>"
# Actively search for offset tags, ignoring other lines until the end of the index is found.
while !eof(stream)
line = readline(stream)
# First, check for the end condition.
if occursin("</index>", line) || occursin("</indexedmzML>", line)
break
end
# If it's not the end, see if it's an offset tag.
m = match(offset_regex, line)
if m !== nothing
push!(offsets, parse(Int64, m.captures[1]))
end
# If it's neither, ignore the line and continue the loop.
end
return offsets
end
"""
find_index_offset(stream::IO)::Int64
Finds the index offset in an mzML file by reading from the end.
Optimized version with better buffer management.
"""
function find_index_offset(stream::IO)::Int64
file_size = filesize(stream)
seekend(stream)
# Read larger chunk for better chance of finding the offset
chunk_size = min(8192, file_size)
seek(stream, file_size - chunk_size)
footer = read(stream, String)
index_offset_match = match(r"<indexListOffset>(\d+)</indexListOffset>", footer)
if index_offset_match === nothing
throw(FileFormatError("Could not find <indexListOffset>. File may not be an indexed mzML."))
end
return parse(Int64, index_offset_match.captures[1])
end
# This is the main lazy-loading function for mzML, now returning an MSIData object.
"""
load_mzml_lazy(file_path::String; cache_size::Int=100)
Lazily loads an indexed `.mzML` file by parsing only the metadata. It reads the
spectrum index from the end of the file to get the offsets of each spectrum,
then parses the metadata for each spectrum without loading the binary data.
# Arguments
- `file_path`: The path to the `.mzML` file.
- `cache_size`: The number of spectra to hold in an LRU cache for faster access.
# Returns
- An `MSIData` object ready for lazy data access.
"""
function load_mzml_lazy(file_path::String; cache_size::Int=100)
println("DEBUG: Opening file stream for $file_path")
ts_stream = ThreadSafeFileHandle(file_path, "r")
try
# --- NEW: Parse instrument metadata from header ---
println("DEBUG: Parsing instrument metadata from header...")
instrument_meta = parse_instrument_metadata_mzml(ts_stream.handle)
println("--- Extracted Instrument Metadata ---")
println("Resolution: ", instrument_meta.resolution)
println("Acquisition Mode (pre-check): ", instrument_meta.acquisition_mode)
println("Calibration Status: ", instrument_meta.calibration_status)
println("Instrument Model: ", instrument_meta.instrument_model)
println("Mass Accuracy (ppm): ", instrument_meta.mass_accuracy_ppm)
println("Laser Settings: ", instrument_meta.laser_settings)
println("Polarity: ", instrument_meta.polarity)
println("------------------------------------")
seekstart(ts_stream.handle) # Reset stream after header parsing
println("DEBUG: Finding index offset...")
index_offset = find_index_offset(ts_stream.handle)
println("DEBUG: Seeking to index list at offset $index_offset.")
seek(ts_stream.handle, index_offset)
println("DEBUG: Searching for '<index name=\"spectrum\">'.")
if find_tag(ts_stream.handle, r"<index\s+name=\"spectrum\"") === nothing
throw(FileFormatError("Could not find spectrum index."))
end
println("DEBUG: Found spectrum index tag.")
println("DEBUG: Parsing spectrum offsets...")
spectrum_offsets = parse_offset_list(ts_stream.handle)
if isempty(spectrum_offsets)
throw(FileFormatError("No spectrum offsets found."))
end
num_spectra = length(spectrum_offsets)
println("DEBUG: Found $num_spectra spectrum offsets.")
println("DEBUG: Parsing metadata for each spectrum...")
# Pre-allocate the metadata vector for better performance
spectra_metadata = Vector{SpectrumMetadata}(undef, num_spectra)
# Use @inbounds for faster indexing in the loop
@inbounds for i in 1:num_spectra
spectra_metadata[i] = parse_spectrum_metadata(ts_stream.handle, spectrum_offsets[i])
# Progress reporting for large files
if i % 1000 == 0
println("DEBUG: Processed $i/$num_spectra spectra")
end
end
println("DEBUG: Metadata parsing complete for all $num_spectra spectra.")
# Assuming uniform data formats, take from the first spectrum
first_meta = spectra_metadata[1]
mz_format = first_meta.mz_asset.format
intensity_format = first_meta.int_asset.format
println("DEBUG: Inferred global m/z format: $mz_format")
println("DEBUG: Inferred global intensity format: $intensity_format")
# --- NEW: Determine overall acquisition mode ---
modes = [meta.mode for meta in spectra_metadata]
num_centroid = count(m -> m == CENTROID, modes)
num_profile = count(m -> m == PROFILE, modes)
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
println("DEBUG: Inferred overall acquisition mode: $acq_mode_symbol (Centroid: $num_centroid, Profile: $num_profile)")
final_instrument_meta = InstrumentMetadata(
instrument_meta.resolution,
acq_mode_symbol, # Update with parsed mode
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 # Add this new field
)
source = MzMLSource(ts_stream, mz_format, intensity_format)
println("DEBUG: Creating MSIData object.")
return MSIData(source, spectra_metadata, final_instrument_meta, (0, 0), nothing, cache_size)
catch e
close(ts_stream) # Ensure stream is closed on error
rethrow(e)
end
end

18
start_MSI_GUI.jl
View File

@ -1,14 +1,32 @@
using Pkg
sTime=time()
Pkg.activate(".")
ENV["GENIE_ENV"] = "dev"
#ENV["GENIE_ENV"] = "prod"
manifest_path = joinpath(@__DIR__, "Manifest.toml")
# Only instantiate in development mode
if get(ENV, "GENIE_ENV", "dev") != "prod" || !isfile(manifest_path)
@info "Development environment detected. Instantiating packages..."
if !isfile(manifest_path)
@info "Manifest.toml not found. Generating it based on Project.toml..."
end
Pkg.resolve()
Pkg.instantiate()
Pkg.gc()
end
using Genie
# Load and configure Genie
Genie.loadapp()
# Remove html parser error discrepancy
redirect_stderr(devnull)
# Start the Genie server
@async begin
up(host="127.0.0.1", port=1481)

932
test/benchmark.jl Normal file
View File

@ -0,0 +1,932 @@
# 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/imzML_LTP/ltpmsi-chilli.imzML",420,0.1, "Chilli Pepper"), # Chilli
BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_DESI/ColAd_Individual/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/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"), # 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

81
test/readme.md Normal file
View File

@ -0,0 +1,81 @@
# JuliaMSI Test Suite
This document provides instructions on how to set up and run the test suite for the `JuliaMSI` package. The tests validate the core functionality of the data processing workflows, including loading, converting, and visualizing mass spectrometry data.
## Local Installation and Setup
### 1. Prerequisites
Make sure you have Julia installed (at least version 1.6). If not, we recommend using `juliaup` for installation.
- **juliaup**: [https://github.com/JuliaLang/juliaup](https://github.com/JuliaLang/juliaup)
- **Official binaries**: [https://julialang.org/downloads/](https://julialang.org/downloads/)
### 2. Download the Project
Download and decompress the repository. You can get it from the official Codeberg page or clone it using Git.
- **Download ZIP**: [https://codeberg.org/LabABI/JuliaMSI/archive/main.zip](https://codeberg.org/LabABI/JuliaMSI/archive/main.zip)
- **Git URL**: `https://codeberg.org/LabABI/JuliaMSI.git`
Ensure you know the directory where the `JuliaMSI` files are located.
### 3. Download Test Data
The test script requires example mass spectrometry data files. You can download the required test data from Zenodo:
- **Example Data**: [https://doi.org/10.5281/zenodo.10084132](https://doi.org/10.5281/zenodo.10084132)
Download the data and place it in a known location on your computer.
### 4. Configure the Test Script
Before running the tests, you must edit the `test/run_tests.jl` file to point to the test data you downloaded.
1. Open `test/run_tests.jl` in a text editor.
2. Locate the `CONFIG: PLEASE FILL IN YOUR FILE PATHS HERE` section.
3. Update the constant variables (e.g., `TEST_MZML_FILE`, `CONVERSION_SOURCE_MZML`, `CONVERSION_SYNC_FILE`) with the absolute paths to the corresponding files on your system.
## Running the Tests
1. **Navigate to the Project Directory**:
Open your terminal and use the `cd` command to navigate to the root folder of the `JuliaMSI` project.
```bash
cd /path/to/your/JuliaMSI
```
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.
```bash
julia --threads auto --project=. test/run_tests.jl
```
```bash
julia --threads auto --project=. test/run_preprocessing.jl
```
3. **Check the Results**:
The script will print its progress to the console. Any generated images (plots and image slices) will be saved in the `test/results/` directory.
## Test Case Configuration
You can customize the test run by editing the variables in `test/run_tests.jl` or `test/run_preprocessing.jl` respectively.
### Enabling and Disabling Test Cases
You can run or skip specific test cases of run_tests by setting the corresponding boolean variables to `true` or `false`.
```julia
test1 = true # Runs Test Case 1
test2 = true # Runs Test Case 2
test3 = true # Runs Test Case 3
```
- **Test Case 1**: Validates a standard `.mzML` file and generates plots for a single spectrum, the total spectrum, and the average spectrum.
- **Test Case 2**: Tests the conversion of a `.mzML` file (and its corresponding `.txt` sync file) into the `.imzML` format. It then validates the newly created file.
- **Test Case 3**: Validates an existing `.imzML` file. It generates plots for a spectrum at specific coordinates, the total spectrum, the average spectrum, and an ion image slice.
### Input Configuration Variables
All configuration variables are located in the `CONFIG` section of `test/run_tests.jl`.
- `TEST_MZML_FILE`: Path to the standard `.mzML` file for Test Case 1.
- `SPECTRUM_TO_PLOT`: The index of the spectrum to plot from the `.mzML` file in Test Case 1.
- `CONVERSION_SOURCE_MZML`: Path to the imaging `.mzML` file to be converted in Test Case 2.
- `CONVERSION_SYNC_FILE`: Path to the `.txt` synchronization file corresponding to `CONVERSION_SOURCE_MZML`.
- `CONVERSION_TARGET_IMZML`: The output path for the `.imzML` file generated in Test Case 2. This is also used as the default input for Test Case 3.
- `TEST_IMZML_FILE`: Path to the `.imzML` file to be tested in Test Case 3.
- `MZ_VALUE_FOR_SLICE`: The m/z value to use for creating an ion image slice in Test Case 3.
- `MZ_TOLERANCE`: The tolerance (+/-) for the `MZ_VALUE_FOR_SLICE` when generating the image.
- `COORDS_TO_PLOT`: An `(X, Y)` tuple specifying the coordinates of the spectrum to plot from the `.imzML` file in Test Case 3.
- `RESULTS_DIR`: The directory where output images will be saved.

323
test/run_descriptive_preprocessing.jl Normal file
View File

@ -0,0 +1,323 @@
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 Normal file
View File

@ -0,0 +1,157 @@
# 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()

281
test/run_precalculation_example.jl Normal file
View File

@ -0,0 +1,281 @@
# 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 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
=#

711
test/run_preprocessing.jl Normal file
View File

@ -0,0 +1,711 @@
# run_preprocessing.jl
using Printf
import Pkg
using CairoMakie
using DataFrames # For creating dataframes
using CSV
using Statistics
using Interpolations
# --- Load the MSI_src Module ---
Pkg.activate(joinpath(@__DIR__, ".."))
using MSI_src
# ===================================================================
# CONFIG: PLEASE FILL IN YOUR FILE PATHS HERE
# ===================================================================
# const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Thricoderma_etc/Imaging_interaccion_trichoderma_vs_streptomyces.imzML"
# const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/set de datos MS/Atropina_tuneo_fraq_20ev.mzML"
const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML"
# const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML"
const MASK_ROUTE = "/home/pixel/Documents/Cinvestav_2025/JuliaMSI/public/css/masks/Stomach_DHB_uncompressed.png"
# const MASK_ROUTE = ""
const OUTPUT_DIR = "./test/results/preprocessing_results"
# Reference peaks for calibration and alignment
const reference_peaks = Dict(
137.0244 => "DHB_fragment",
155.0349 => "DHB_M+H",
177.0168 => "DHB_M+Na",
496.3398 => "PC_16:0_16:0",
520.3398 => "PC_16:0_18:1",
760.5851 => "PC_16:0_18:1_Na",
391.2843 => "PDMS",
413.2662 => "PDMS_Na",
842.5092 => "Protein_standard",
1045.532 => "Protein_standard",
290.1747 => "Atropine [M+H]+",
304.1903 => "Scopolamine [M+H]+",
124.0393 => "Tropine [M+H]+",
)
const PIPELINE_STP = [
"stabilization",
#"baseline_correction",
"smoothing",
"peak_picking",
"peak_selection",
#"calibration",
"peak_alignment",
"normalization",
"peak_binning"
]
# ===================================================================
# USER OVERRIDES: Manually specify parameters here
# ===================================================================
# This dictionary allows you to override any auto-detected parameters.
# The structure should match the output of `main_precalculation`.
# Example: Force a less aggressive peak prominence threshold.
const USER_OVERRIDES = Dict(
:Stabilization => Dict(
:method => :sqrt # Default stabilization method
),
:PeakPicking => Dict(
#:snr_threshold => 1.5,
#:half_window => 5,
:snr_threshold => 8.0,
#:merge_peaks_tolerance => 2.5,
#:half_window => 2
#:half_window => 3
),
# :Smoothing => Dict(
# :window => 11
# )
#:PeakAlignment => Dict(
# :tolerance => 0.01
#)
#:PeakSelection => Dict(
#:frequency_threshold => 0,
#:min_shape_r2 => 0.5,
#:max_fwhm_ppm => 30.0,
#:min_fwhm_ppm => 2.0
#),
#:PeakBinning => Dict(
#:tolerance => 30.0,
#:frequency_threshold => 0
#)
)
# ===================================================================
# HELPER FUNCTIONS
# ===================================================================
function ensure_output_dir()
if !isdir(OUTPUT_DIR)
mkdir(OUTPUT_DIR)
println("Created output directory: $OUTPUT_DIR")
end
end
function plot_spectrum_step(mz, intensity, step_name; peaks=nothing, spectrum_index=1)
fig = Figure(size=(1200, 600))
ax = Axis(fig[1, 1], xlabel="m/z", ylabel="Intensity", title="Step: $step_name (Spectrum $spectrum_index)")
lines!(ax, mz, intensity, color=:blue, linewidth=1)
if peaks !== nothing && !isempty(peaks)
peak_mz = [p.mz for p in peaks]
peak_intensity = [p.intensity for p in peaks]
scatter!(ax, peak_mz, peak_intensity, color=:red, markersize=5)
end
filename = joinpath(OUTPUT_DIR, "spectrum_$(step_name)_index_$spectrum_index.png")
save(filename, fig)
println(" - Saved plot: $(basename(filename))")
end
function print_step_header(step_name::String)
println("\n" * ">"^50)
println("APPLYING STEP: $step_name")
println(">"^50)
end
# ===================================================================
# PREPROCESSING PIPELINE FUNCTIONS (IN-PLACE)
# ===================================================================
function apply_baseline_correction_core(spectra::Vector{MutableSpectrum}, params::Dict, msi_data::MSIData)
print_step_header("Baseline Correction")
method = get(params, :method, :snip)
iterations = get(params, :iterations, 100)
window = get(params, :window, 20)
println(" Method: $method, Iterations: $iterations, Window: $window")
# Safely plot the first spectrum
if !isempty(spectra)
s = spectra[1]
if validate_spectrum(s.mz, s.intensity)
baseline = MSI_src.apply_baseline_correction_core(s.intensity; method=method, iterations=iterations, window=window)
corrected_intensity = max.(0.0, s.intensity .- baseline)
plot_spectrum_step(s.mz, corrected_intensity, "baseline_correction")
end
end
Threads.@threads for s in spectra
if validate_spectrum(s.mz, s.intensity)
baseline = MSI_src.apply_baseline_correction_core(s.intensity; method=method, iterations=iterations, window=window)
s.intensity = max.(0.0, s.intensity .- baseline)
end
end
end
function apply_smoothing(spectra::Vector{MutableSpectrum}, params::Dict)
print_step_header("Smoothing")
method = get(params, :method, :savitzky_golay)
window = get(params, :window, 9)
order = get(params, :order, 2)
println(" Method: $method, Window: $window, Order: $order")
# Safely plot the first spectrum
if !isempty(spectra)
s = spectra[1]
if validate_spectrum(s.mz, s.intensity)
smoothed_intensity = max.(0.0, smooth_spectrum_core(s.intensity; method=method, window=window, order=order))
plot_spectrum_step(s.mz, smoothed_intensity, "smoothing", spectrum_index=s.id)
end
end
Threads.@threads for s in spectra
if validate_spectrum(s.mz, s.intensity)
smoothed_intensity = max.(0.0, smooth_spectrum_core(s.intensity; method=method, window=window, order=order))
s.intensity = smoothed_intensity
end
end
end
function apply_peak_picking(spectra::Vector{MutableSpectrum}, params::Dict)
print_step_header("Peak Picking")
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)
println(" Method: $method, SNR: $snr_threshold, Half Window: $half_window")
found_invalid_spectrum = Threads.Atomic{Bool}(false)
Threads.@threads for s in spectra
is_valid = validate_spectrum(s.mz, s.intensity)
if !is_valid && !found_invalid_spectrum[]
if Threads.atomic_xchg!(found_invalid_spectrum, true) == false
# Debug printing for the first invalid spectrum found
end
end
if is_valid
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(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
else
s.peaks = []
end
end
# Safely plot the first spectrum
if !isempty(spectra)
s1 = findfirst(s -> s.id == 1, spectra)
if s1 !== nothing
plot_spectrum_step(spectra[s1].mz, spectra[s1].intensity, "peak_picking", peaks=spectra[s1].peaks)
println(" - Detected $(length(spectra[s1].peaks)) peaks in spectrum 1")
end
end
end
function apply_peak_selection(spectra::Vector{MutableSpectrum}, params::Dict)
print_step_header("Peak Selection")
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` values from params, default to non-filtering values
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
println(" - Min SNR: $min_snr")
println(" - FWHM Range (ppm): [$min_fwhm, $max_fwhm]")
println(" - Min Shape R²: $min_r2")
total_peaks_before = sum(s -> length(s.peaks), spectra)
Threads.@threads for s in spectra
if !isempty(s.peaks)
s.peaks = filter(p ->
p.snr >= min_snr &&
(min_fwhm <= p.fwhm <= max_fwhm) &&
p.shape_r2 >= min_r2,
s.peaks
)
end
end
total_peaks_after = sum(s -> length(s.peaks), spectra)
println(" - Peaks before: $total_peaks_before, Peaks after: $total_peaks_after")
# Safely plot the first spectrum to show effect of filtering
if !isempty(spectra)
s1_idx = findfirst(s -> s.id == 1, spectra)
if s1_idx !== nothing
s1 = spectra[s1_idx]
plot_spectrum_step(s1.mz, s1.intensity, "peak_selection", peaks=s1.peaks, spectrum_index=s1.id)
println(" - Spectrum 1 now has $(length(s1.peaks)) peaks after selection.")
end
end
end
function apply_calibration(spectra::Vector{MutableSpectrum}, params::Dict, reference_peaks::Dict)
print_step_header("Calibration")
method = get(params, :method, :none)
ppm_tolerance = get(params, :ppm_tolerance, 20.0)
println(" Method: $method, PPM Tolerance: $ppm_tolerance")
if method == :none || isempty(reference_peaks)
println(" - Skipping calibration (no method or reference peaks)")
return
end
reference_masses = collect(keys(reference_peaks))
calibration_info = Vector{String}(undef, length(spectra))
Threads.@threads for i in 1:length(spectra)
s = spectra[i]
info_message = ""
if validate_spectrum(s.mz, s.intensity)
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())
s.mz = itp(s.mz) # Modify mz-axis in-place
info_message = " - Spectrum $(s.id): calibrated using $(length(matched_peaks)) reference peaks"
else
info_message = " - Spectrum $(s.id): insufficient reference peaks ($(length(matched_peaks)) found), skipping"
end
end
calibration_info[i] = info_message
end
# Print summary of results
printed_count = 0
for info in calibration_info
if !isempty(info) && printed_count < 5
println(info)
printed_count += 1
end
end
end
function apply_peak_alignment(spectra::Vector{MutableSpectrum}, params::Dict)
print_step_header("Peak Alignment")
method = get(params, :method, :none)
tolerance = get(params, :tolerance, 0.002)
tolerance_unit = get(params, :tolerance_unit, :mz)
println(" Method: $method, Tolerance: $tolerance $tolerance_unit")
if method == :none
println(" - Skipping peak alignment")
return
end
ref_find_idx = findfirst(s -> !isempty(s.peaks), spectra)
if ref_find_idx === nothing
println(" - 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]
println(" - Using spectrum $(ref_spectrum.id) as reference with $(length(ref_peaks_mz)) 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)
s.mz = alignment_func.(s.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
function apply_normalization_core(spectra::Vector{MutableSpectrum}, params::Dict)
print_step_header("Normalization")
method = get(params, :method, :tic)
println(" Method: $method")
# Safely plot the first spectrum
if !isempty(spectra)
s1 = findfirst(s -> s.id == 1, spectra)
if s1 !== nothing
s = spectra[s1]
if validate_spectrum(s.mz, s.intensity)
normalized_intensity = MSI_src.apply_normalization_core(s.intensity; method=method)
plot_spectrum_step(s.mz, normalized_intensity, "normalization")
end
end
end
Threads.@threads for s in spectra
if validate_spectrum(s.mz, s.intensity)
s.intensity = MSI_src.apply_normalization_core(s.intensity; method=method)
end
end
end
function apply_peak_binning(spectra::Vector{MutableSpectrum}, params::Dict)
print_step_header("Peak Binning")
method = get(params, :method, :adaptive)
tolerance = get(params, :tolerance, 20.0)
tolerance_unit = get(params, :tolerance_unit, :ppm)
min_peak_per_bin = get(params, :min_peak_per_bin, 3)
println(" Method: $method, Tolerance: $tolerance $tolerance_unit")
if isempty(spectra) || all(s -> isempty(s.peaks), spectra)
@warn "No peaks found for binning. Returning empty feature matrix."
return nothing, nothing
end
println(" - Binning peaks from $(length(spectra)) spectra")
# Collect all peaks with their intensities
all_peaks = Vector{Tuple{Float64, Float64}}() # (mz, intensity)
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])
# Create bins - just store mz_center and intensity
bin_centers = Float64[]
bin_intensities = Float64[]
i = 1
while i <= length(all_peaks)
current_bin_start = i
current_peak = all_peaks[i]
# Find all peaks in this bin
j = i + 1
while j <= length(all_peaks)
next_peak = all_peaks[j]
# Calculate tolerance
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
# Check if we have enough peaks in this bin
if bin_size >= min_peak_per_bin
bin_peaks_core = all_peaks[current_bin_start:current_bin_end]
# Calculate m/z center and average intensity
mz_sum = 0.0
intensity_sum = 0.0
for peak in bin_peaks_core
mz_sum += peak[1]
intensity_sum += peak[2]
end
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 # Move to next potential bin
end
# Create the 2-row matrix
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
println(" - Created feature matrix: 2 × $n_bins")
println(" - Number of bins created: $n_bins")
println(" - m/z range: $(round(bin_centers[1], digits=4)) - $(round(bin_centers[end], digits=4))")
# Return bin info as a vector of tuples for compatibility
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
function save_feature_matrix(feature_matrix::Matrix{Float64}, bin_info)
print_step_header("Saving Results")
# 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 header
write(io, "mz,intensity\n")
# Write data: m/z values in first column, intensities in second
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
println(" - 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 header with m/z values as column names
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 the aggregated intensity values
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
println(" - Saved standard format matrix: $csv_path_standard")
return csv_path, csv_path_standard
end
# ===================================================================
# MAIN PREPROCESSING PIPELINE
# ===================================================================
function run_preprocessing_pipeline()
println("="^80)
println("MSI PREPROCESSING PIPELINE")
println("="^80)
# Ensure output directory exists
ensure_output_dir()
# Load data
println("\nLoading data: $(basename(TEST_FILE))")
msi_data = OpenMSIData(TEST_FILE)
# Precompute analytics
println("\nPrecomputing analytics...")
precompute_analytics(msi_data)
# Run precalculation to get automatic parameters
println("\nRunning precalculation for automatic parameter determination...")
auto_params = main_precalculation(msi_data, reference_peaks=reference_peaks, mask_path=MASK_ROUTE)
# --- Apply User Overrides ---
if !isempty(USER_OVERRIDES)
println("\nApplying user overrides...")
for (step, params) in USER_OVERRIDES
if haskey(auto_params, step)
for (param, value) in params
@info "OVERRIDE: For step `:$step`, setting `:$param` to `$value` (was `$(get(auto_params[step], param, "not set"))`)"
end
merge!(auto_params[step], params)
else
@warn "User override for non-existent step `:$step` ignored."
end
end
end
# --- Final Parameter Summary ---
println("\n" * "-"^60)
println("FINAL PARAMETERS FOR PIPELINE RUN")
println("-"^60)
for (step, params) in sort(collect(pairs(auto_params)), by=p -> p.first)
println(" - Step: $step")
if isempty(params)
println(" (No parameters)")
continue
end
for (param, value) in sort(collect(pairs(params)), by=p -> p.first)
println(" - $(rpad(param, 25)): $value")
end
end
# Define pipeline steps
pipeline_steps = PIPELINE_STP
println("\nPipeline steps: $(join(pipeline_steps, " -> "))")
# Initialize `current_spectra` as a Vector of mutable structs for in-place modification
println("\nInitializing spectra data structure...")
local spectrum_indices_to_process::AbstractVector{Int}
if !isempty(MASK_ROUTE)
println("Applying mask from: $(MASK_ROUTE)")
try
mask_matrix = MSI_src.load_and_prepare_mask(MASK_ROUTE, msi_data.image_dims)
masked_indices_set = MSI_src.get_masked_spectrum_indices(msi_data, mask_matrix)
spectrum_indices_to_process = collect(masked_indices_set)
println("Mask applied. $(length(spectrum_indices_to_process)) spectra are within the masked region.")
catch e
@error "Failed to load or apply mask: $e. Proceeding without mask."
spectrum_indices_to_process = 1:length(msi_data.spectra_metadata)
end
else
spectrum_indices_to_process = 1:length(msi_data.spectra_metadata)
end
if isempty(spectrum_indices_to_process)
@warn "No spectra available for processing after applying mask/filter. Exiting pipeline."
close(msi_data)
return
end
num_spectra_to_process = length(spectrum_indices_to_process)
current_spectra = Vector{MutableSpectrum}(undef, num_spectra_to_process)
Threads.@threads for i in 1:num_spectra_to_process
original_idx = spectrum_indices_to_process[i]
mz, intensity = MSI_src.GetSpectrum(msi_data, original_idx)
current_spectra[i] = MSI_src.MutableSpectrum(original_idx, mz, intensity, [])
# Plot raw spectrum without masks (only the first processed one)
if i == 1
plot_spectrum_step(mz, intensity, "raw_unmasked_spectrum", spectrum_index=original_idx)
end
end
# These variables will be populated by the pipeline steps
feature_matrix = nothing
bin_definitions = nothing
# Apply pipeline steps, modifying `current_spectra` in-place
for step in pipeline_steps
println("\n" * "-"^60)
println("PROCESSING STEP: $step")
println("-"^60)
if step == "stabilization"
print_step_header("Intensity Transformation (Stabilization)")
method = get(auto_params[:Stabilization], :method, :sqrt)
println(" Method: $method")
# Safely plot the first spectrum before transformation
if !isempty(current_spectra)
s = current_spectra[1]
if validate_spectrum(s.mz, s.intensity)
# Use a temporary spectrum to plot before and after
initial_intensity = deepcopy(s.intensity)
plot_spectrum_step(s.mz, initial_intensity, "stabilization_before", spectrum_index=s.id)
end
end
@time apply_intensity_transformation(current_spectra, auto_params[:Stabilization])
# Safely plot the first spectrum after transformation
if !isempty(current_spectra)
s = current_spectra[1]
if validate_spectrum(s.mz, s.intensity)
plot_spectrum_step(s.mz, s.intensity, "stabilization_after", spectrum_index=s.id)
end
end
elseif step == "baseline_correction"
@time apply_baseline_correction_core(current_spectra, auto_params[:BaselineCorrection], msi_data)
elseif step == "smoothing"
apply_smoothing(current_spectra, auto_params[:Smoothing])
elseif step == "peak_picking"
@time apply_peak_picking(current_spectra, auto_params[:PeakPicking])
elseif step == "peak_selection"
@time apply_peak_selection(current_spectra, auto_params[:PeakSelection])
elseif step == "calibration"
@time apply_calibration(current_spectra, auto_params[:Calibration], reference_peaks)
elseif step == "peak_alignment"
@time apply_peak_alignment(current_spectra, auto_params[:PeakAlignment])
elseif step == "normalization"
@time apply_normalization_core(current_spectra, auto_params[:Normalization])
elseif step == "peak_binning"
# This step is different as it generates the final matrix, not modifying spectra in-place
feature_matrix, bin_info = @time apply_peak_binning(current_spectra, auto_params[:PeakBinning])
if feature_matrix !== nothing
@time save_feature_matrix(feature_matrix, bin_info)
end
else
@warn "Unknown step: $step, skipping"
end
# Print progress
println("✓ Completed step: $step")
end
# Clean up
close(msi_data)
GC.gc()
println("\n" * "="^80)
println("PREPROCESSING PIPELINE COMPLETED SUCCESSFULLY!")
println("Results saved to: $OUTPUT_DIR")
println("="^80)
end
# ===================================================================
# EXECUTE PIPELINE
# ===================================================================
if abspath(PROGRAM_FILE) == @__FILE__
@time run_preprocessing_pipeline()
end

347
test/run_tests.jl Normal file
View File

@ -0,0 +1,347 @@
# test/run_tests.jl
# ===================================================================
# Test Environment for JuliaMSI Package
# ===================================================================
# This script validates the core functionality of the data processing
# workflows, including loading, converting, and visualizing mass
# spectrometry data.
#
# 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/run_tests.jl
# 3. Check the `test/results/` folder for the output images.
# ===================================================================
using Printf
using CairoMakie
import Pkg
# --- Load the MSI_src Module ---
# Activate the project environment at the parent directory of this test script
Pkg.activate(joinpath(@__DIR__, ".."))
using MSI_src
# ===================================================================
# CONFIG: PLEASE FILL IN YOUR FILE PATHS HERE
# ===================================================================
# --- Test Case 1: Standard .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/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 prueba Roya 1/Roya.mzML"
const SPECTRUM_TO_PLOT = 1 # Which spectrum to plot from the file
# --- Test Case 2: .mzML + Sync File for Conversion ---
# The special .mzML file with one spectrum per pixel.
# const CONVERSION_SOURCE_MZML = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML"
const CONVERSION_SOURCE_MZML = TEST_MZML_FILE
# The corresponding synchronization text file.
const CONVERSION_SYNC_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.txt"
# const CONVERSION_SYNC_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging prueba Roya 1/Synchro.txt"
# const CONVERSION_SOURCE_MZML = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging_paper_spray/Imaging_paper_spray.mzML"
# const CONVERSION_SYNC_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging_paper_spray/Imaging_paper_spray.txt"
# The desired output path for the new .imzML file.
const CONVERSION_TARGET_IMZML = "test/results/converted_mzml.imzML"
# --- Test Case 3: Standard .imzML file ---
# An existing imzML file (can be the one generated from Case 2).
# const TEST_IMZML_FILE = CONVERSION_TARGET_IMZML # The output from case 2
# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML"
# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging_paper_spray/Imaging_paper_spray.imzML" #profile
# 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/Stomach_DHB_compressed.imzML" # centroid compressed
const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML" # centroid
# The m/z value to use for creating an image slice.
# const MZ_VALUE_FOR_SLICE = 309.06 # BF
# const MZ_VALUE_FOR_SLICE = 896.0 # HR2MSI
# const MZ_VALUE_FOR_SLICE = 76.03 # I PS
# const MZ_VALUE_FOR_SLICE = 313 # ROYA
const MZ_VALUE_FOR_SLICE = 100 # advanced processing
# const MZ_TOLERANCE = 0.1
# const MZ_TOLERANCE = 1
const MZ_TOLERANCE = 0.1
# Coordinates to plot a specific spectrum from imzML
const COORDS_TO_PLOT = (50, 50) # Example coordinates (X, Y)
# --- Output Directory ---
const RESULTS_DIR = "test/results"
test1 = true
test2 = false
test3 = true
# ===================================================================
# DATA VALIDATION UTILITY
# ===================================================================
"""
validate_msi_data(filepath::String)
Performs a series of checks on a .mzML or .imzML file using the MSIData API.
"""
function validate_msi_data(filepath::String)
println("" * "-"^10 * " Running validation for $filepath " * "-"^10)
if !isfile(filepath)
println("SKIPPED VALIDATION: File not found: $filepath")
return false
end
try
# 1. Basic structure validation
println("Opening file with OpenMSIData...")
msi_data = @time OpenMSIData(filepath)
# 2. Compare spectrum counts
num_spectra = length(msi_data.spectra_metadata)
println("Found $num_spectra spectra")
@assert num_spectra > 0 "No spectra found in file."
# 3. Test random access
println("Testing random access to spectra...")
test_indices = unique([1, max(1, num_spectra ÷ 2), num_spectra])
println("Testing indices: $test_indices")
for idx in test_indices
print("Fetching spectrum #$idx... ")
@time process_spectrum(msi_data, idx) do mz, intensity
@assert length(mz) == length(intensity) "Spectrum $idx: mz/intensity length mismatch. Got $(length(mz)) mz values and $(length(intensity)) intensity values."
println("OK, $(length(mz)) points.")
end
end
# 4. Test iteration
println("Testing iteration over all spectra...")
count = 0
iter_time = @elapsed for (idx, (mz, intensity)) in IterateSpectra(msi_data)
count += 1
# Basic data validation
@assert all(isfinite, mz) "Non-finite mz values in spectrum $idx"
@assert all(>=(0), intensity) "Negative intensities in spectrum $idx"
end
println("Iterated over $count spectra in $iter_time seconds.")
@assert count == num_spectra "Iteration count mismatch: expected $num_spectra, got $count."
println("VALIDATION SUCCESSFUL for $filepath")
return true
catch e
println("VALIDATION FAILED for $filepath.")
showerror(stdout, e, catch_backtrace())
println()
return false
end
end
function debug_xml_parsing(file_path::String)
println("=== DEBUG XML PARSING ===")
stream = open(file_path, "r")
# Find and print the first spectrum
while !eof(stream)
line = readline(stream)
if occursin("<spectrum", line)
println("FOUND FIRST SPECTRUM:")
spectrum_xml = line
# Read until end of spectrum
while !eof(stream) && !occursin("</spectrum>", line)
line = readline(stream)
spectrum_xml *= line
end
println("SPECTRUM XML:")
println(spectrum_xml)
break
end
end
close(stream)
println("=== END DEBUG ===")
end
# ===================================================================
# TEST RUNNER
# ===================================================================
function run_test()
println("Starting MSI_src Test Suite...")
# --- Test Case 1: Process a standard .mzML file ---
println("" * "="^20 * " Test Case 1: Processing .mzML " * "="^20)
if test1 == true
# Run new, stronger validation
validate_msi_data(TEST_MZML_FILE)
# Also run original plotting test to ensure visualization still works
if isfile(TEST_MZML_FILE)
try
println("Plotting a sample spectrum from $TEST_MZML_FILE...")
msi_data = @time OpenMSIData(TEST_MZML_FILE)
process_spectrum(msi_data, SPECTRUM_TO_PLOT) do mz, intensity
fig = Figure(size = (800, 600))
ax = Axis(fig[1, 1], xlabel="m/z", ylabel="Intensity", title="Spectrum #$SPECTRUM_TO_PLOT from $(basename(TEST_MZML_FILE))")
lines!(ax, mz, intensity)
output_path = joinpath(RESULTS_DIR, "test_mzml_spectrum.png")
save(output_path, fig)
println("SUCCESS: Spectrum plot saved to $output_path")
end
# Get the summed spectrum data
mz, intensity = get_total_spectrum(msi_data)
# Plot the data
println("Plotting total spectrum...")
fig = Figure(size = (800, 600))
ax = Axis(fig[1, 1], xlabel="m/z", ylabel="Total Intensity", title="Total Spectrum from $(basename(TEST_MZML_FILE))")
lines!(ax, mz, intensity)
# Saving the output
output_path = joinpath(RESULTS_DIR, "test_mzml_total_spectrum.png")
save(output_path, fig)
println("SUCCESS: Total spectrum plot saved to $output_path")
# Get the averaged spectrum data
mz, intensity = get_average_spectrum(msi_data)
# Plot the data
println("Plotting averaged spectrum...")
fig = Figure(size = (800, 600))
ax = Axis(fig[1, 1], xlabel="m/z", ylabel="Average Intensity", title="Average Spectrum from $(basename(TEST_MZML_FILE))")
lines!(ax, mz, intensity)
# Saving the output
output_path = joinpath(RESULTS_DIR, "test_mzml_average_spectrum.png")
save(output_path, fig)
println("SUCCESS: Total spectrum plot saved to $output_path")
catch e
println("ERROR during plotting in Test Case 1: $e")
end
end
else
println("SKIPPED Test Case 1.")
end
# --- Test Case 2: Convert .mzML + .txt to .imzML ---
println("" * "="^20 * " Test Case 2: Converting to .imzML " * "="^20)
if isfile(CONVERSION_SOURCE_MZML) && isfile(CONVERSION_SYNC_FILE) && test2 == true
try
println("Running conversion process (with profiling)...")
success = @time ImportMzmlFile(CONVERSION_SOURCE_MZML, CONVERSION_SYNC_FILE, CONVERSION_TARGET_IMZML)
if success
println("SUCCESS: Conversion process completed.")
# Validate the newly created imzML file
validate_msi_data(CONVERSION_TARGET_IMZML)
else
println("FAILURE: Conversion process failed.")
end
catch e
println("ERROR in Test Case 2: $e")
end
else
println("SKIPPED: Files not found for Test Case 2.")
println(" - mzML: $CONVERSION_SOURCE_MZML")
println(" - Sync: $CONVERSION_SYNC_FILE")
end
# --- Test Case 3: Process an existing .imzML file ---
println("" * "="^20 * " Test Case 3: Processing .imzML " * "="^20)
if test3 == true
# Run new, stronger validation for imzML
validate_msi_data(TEST_IMZML_FILE)
# Also run tests for plotting spectrum and image slice
if isfile(TEST_IMZML_FILE)
debug_xml_parsing(TEST_IMZML_FILE)
msi_data = @time OpenMSIData(TEST_IMZML_FILE)
precompute_analytics(msi_data)
# Add spectrum plotting for imzML to match Test Case 1
try
# Get the msi data from the imzml
println("Plotting a sample spectrum from $TEST_IMZML_FILE...")
# Use coordinates from the first spectrum in the metadata
first_spectrum_meta = msi_data.spectra_metadata[1]
x_coord = first_spectrum_meta.x
y_coord = first_spectrum_meta.y
println("DIAGNOSTIC_READ: Reading from coordinate_map[$x_coord, $y_coord]. Value is $(msi_data.coordinate_map[x_coord, y_coord]).")
# Get the x y coordinate spectrum data and plot it using the function barrier
process_spectrum(msi_data, Int(x_coord), Int(y_coord)) do mz, intensity
# Plot the data
fig = Figure(size = (800, 600))
ax = Axis(fig[1, 1], xlabel="m/z", ylabel="Intensity", title="Spectrum at ($x_coord, $y_coord) from $(basename(TEST_IMZML_FILE))")
lines!(ax, mz, intensity)
# Saving the output
output_path = joinpath(RESULTS_DIR, "test_imzml_spectrum.png")
save(output_path, fig)
println("SUCCESS: Spectrum plot saved to $output_path")
end
# Get the summed spectrum data
mz, intensity = get_total_spectrum(msi_data)
# Plot the data
println("Plotting total spectrum...")
fig = Figure(size = (800, 600))
ax = Axis(fig[1, 1], xlabel="m/z", ylabel="Total Intensity", title="Total Spectrum from $(basename(TEST_IMZML_FILE))")
lines!(ax, mz, intensity)
# Saving the output
output_path = joinpath(RESULTS_DIR, "test_imzml_total_spectrum.png")
save(output_path, fig)
println("SUCCESS: Total spectrum plot saved to $output_path")
# Get the averaged spectrum data
mz, intensity = get_average_spectrum(msi_data)
# Plot the data
println("Plotting averaged spectrum...")
fig = Figure(size = (800, 600))
ax = Axis(fig[1, 1], xlabel="m/z", ylabel="Average Intensity", title="Average Spectrum from $(basename(TEST_IMZML_FILE))")
lines!(ax, mz, intensity)
# Saving the output
output_path = joinpath(RESULTS_DIR, "test_imzml_average_spectrum.png")
save(output_path, fig)
println("SUCCESS: Total spectrum plot saved to $output_path")
# println("No spectrums tested on this try.")
catch e
println("ERROR during spectrum plotting in Test Case 3: $e")
end
# Test the plot_slice function
try
println("Testing plot_slice function on $TEST_IMZML_FILE...")
# Use the base peak m/z from the first spectrum for the slice
# Or a more general approach: use the global max intensity m/z
# For now, let's use the base peak m/z of the first spectrum
slice_mz_value = msi_data.spectrum_stats_df.BasePeakMZ[1]
@time plot_slice(msi_data, slice_mz_value, MZ_TOLERANCE, RESULTS_DIR, stage_name="test_imzml_single_slice")
# The success message is now inside plot_slice
catch e
println("ERROR during plot_slice test in Test Case 3: $e")
end
end
else
println("SKIPPED Test Case 3.")
end
println("Tests for all 3 cases is finished.")
end
# --- Execute ---
# Ensure the results directory exists
mkpath(RESULTS_DIR)
@time run_test()