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

finished test script, included instructions in test/readme.md, the test performs mzml, imzml and conversion tests, along with throwing reference plots and images which result from the files

Browse Source
This commit is contained in:
Pixelguy14 2025-09-29 10:57:11 -06:00
parent b1d7838a60
commit cd1115368f
17 changed files with 2145 additions and 2084 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -5,3 +5,4 @@ public/css/imgOver.png
!public/css/LABI_logo.png
log/*
R_original_scripts/
test/results/

1109
Manifest.toml
View File

File diff suppressed because it is too large Load Diff

10
Project.toml
View File

@ -1,6 +1,15 @@
name = "MSI_src"
authors = ["JJSA"]
version = "0.1.0"
[deps]
Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
Genie = "c43c736e-a2d1-11e8-161f-af95117fbd1e"
GenieFramework = "a59fdf5c-6bf0-4f5d-949c-a137c9e2f353"
Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
@ -9,6 +18,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
NativeFileDialog = "e1fe445b-aa65-4df4-81c1-2041507f0fd4"
NaturalSort = "c020b1a1-e9b0-503a-9c33-f039bfc54a85"
PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StipplePlotly = "ec984513-233d-481d-95b0-a3b58b97af2b"
julia_mzML_imzML = "38eb50d3-2fb6-4afa-992a-964ed8562ed9"

6
app.jl
View File

@ -6,7 +6,7 @@ using Libz
using PlotlyBase
using CairoMakie
using Colors
# using julia_mzML_imzML
using julia_mzML_imzML
using Statistics
using NaturalSort
using Images
@ -15,7 +15,7 @@ using NativeFileDialog # Opens the file explorer depending on the OS
using StipplePlotly
using Base.Filesystem: mv # To rename files in the system
include("./julia_imzML_visual.jl")
include("./src/imzML.jl")
@genietools
# == Reactive code ==
@ -363,7 +363,7 @@ include("./src/imzML.jl")
if isfile(full_route) && Nmass > 0 && Tol > 0 && Tol <=1 && colorLevel > 1 && colorLevel < 257
msg="File exists, Nmass=$(Nmass) Tol=$(Tol). Loading file will begin, please be patient."
try
spectra=load_imzml(full_route)
spectra=LoadImzml(full_route)
msg="File loaded. Creating spectra with the specific mass and tolerance, please be patient."
slice=GetMzSliceJl(spectra,Nmass,Tol)
fig=CairoMakie.Figure(size=(150, 250)) # Container

399
src/Bitmap.jl
View File

@ -1,399 +0,0 @@
# ********************************************************************
# SaveBitmap
# ********************************************************************
"""
SaveBitmap( name, pixMap, colorTable )
Save a discretized imzML slice as a bitmap file
# Arguments
* ` name`: Full path name of target bitmap
* ` pixMap`: [UInt8] Bidimensional matrix with image info
* `colorTable`: [UInt32] Vector with RGB colors for each gray level
# Examples
```julia
# Saves a bitmap with black and white alternate image squares
img = zeros(UInt8, 32, 32)
img[17:32,1:16] .= 1
img[1:16,17:32] .= 1
SaveBitmap( "test.bmp", img, UInt32[0, 0x00FFFFFF] )
```
The following image will be created on your hard disk
![](./Test.bmp)
"""
function SaveBitmap( name,
pixMap::Array{UInt8,2},
colorTable::Array{UInt32,1} )
# Get image dimensions
dim = size( pixMap )
if length( dim ) != 2
return 0
end
# Compute row padding
padding = ( 4 - dim[1] & 0x3 ) & 0x3
# Compute file dimensions. Header = 14 + 40 + ( 256 * 4 ) = 1078
offset = 1078
imgBytes = dim[2] * ( dim[1] + padding )
# Create file
stream = open( name, "w" )
# Save file header
write( stream, UInt16( 0x4D42 ) )
write( stream, UInt32[ offset + imgBytes, 0 , offset ] )
# Save info header
write( stream, UInt32[ 40, dim[1], dim[2], 0x80001, 0 ] )
write( stream, UInt32[ imgBytes, 0, 0, 256, 0 ] )
# Save color table
write( stream, colorTable )
if length( colorTable ) < 256
fixTable = zeros( UInt32, 256 - length( colorTable ) )
write( stream, fixTable )
end
# Save image pixels
if padding == 0
for i = 1:dim[2]
write( stream, pixMap[:,i] )
end
else
zeroPad = zeros( UInt8, padding )
for i in 1:dim[2]
write( stream, pixMap[:,i] )
write( stream, zeroPad )
end
end
# Close file
close( stream )
end
# ********************************************************************
# FindMass
# massVector: mz Vector sorted in ascending order
# mass: target mz value
# tolerance: bi-axial tolerance to find mass
# ********************************************************************
function FindMass( massVector, mass, tolerance )
index = Int( 0 )
lower = Int( 1 )
higher = Int( length( massVector ) )
while lower <= higher
# Compute mid element
index = ( lower + higher ) ÷ 2
# Go to lower portion ?
if massVector[ index ] > ( mass + tolerance )
higher = index - 1
continue
end
# Go to Higher portion?
if massVector[ index ] < ( mass - tolerance )
lower = index + 1
continue
end
# mass was found
return index
end
# mass was not found
return 0
end
# ********************************************************************
# GetSlice
# imzML: imzML data
# mass: mz target value
# tolerance: bi-axial tolerance to find mass
# ********************************************************************
"""
GetSlice( imzML, mass, tolerance )
Extract a mz-image from an imzML image data array loaded with
`LoadImzml` function. The resulting image array preserves the same
data type of the y-axis vectors in the file.
# Arguments
* ` imzML`: Image array loaded with `LoadImzml` function
* ` mass`: mz value for image extraction
* `tolerance`: Maximum lateral tolerance for a valid mz search
# Examples
```julia
# Load DESI MSI Carcinoma 885.55 mz slice
spectra = LoadImzml( "80TopL, 50TopR, 70BottomL, 60BottomR-centroid.imzML" )
slice = GetSlice( spectra, 885.55, 0.005 )
```
"""
function GetSlice( imzML, mass, tolerance )
# Alloc space for slice
width = maximum( imzML[1,:] )
height = maximum( imzML[2,:] )
image = zeros( Float64, width, height )
for i in 1:size( imzML )[2]
index = FindMass( imzML[3,i], mass, tolerance )
if index != 0
image[ imzML[1,i], imzML[2,i] ] = imzML[4,i][index]
end
end
return image
end
# ********************************************************************
# IntQuant, Discretize image amplitude in 0:255 range
# slice: Image matriz in measured arbitrary units
# ********************************************************************
"""
IntQuant( slice )
Zero Memory intensity quantizer. Discretize the image?s continuos
intensity range in discrete bins from 0 to 255 gray levels.
Returns a UInt8 bidimensional matrix
# Arguments
* `slice`: slice returned by the `GetSlice` function
# Examples
```julia
# Save DESI MSI Carcinoma 885.55 mz slice as a bitmap file
spectra = LoadImzml( "80TopL, 50TopR, 70BottomL, 60BottomR-centroid.imzML" )
slice = GetSlice( spectra, 885.55, 0.005
SaveBitmap( "Slice.bmp", IntQuant( slice ), ViridisPalette )
```
The following image will be created on your hard disk
![](./Slice.bmp)
"""
function IntQuant( slice )
# Compute scale factor for amplitude discretization
lower = minimum( slice )
scale = 255 / maximum( slice )
dim = size( slice )
image = zeros( UInt8, dim[1], dim[2] )
for i in 1:length( slice )
image[i] = convert( UInt8, floor( slice[i] * scale + 0.5 ) )
end
return image
end
# ********************************************************************
# TrIQ, Discretize image amplitude in 0:255 range, grouping outliers
# in the highest bin
# slice: Image matrix in measured arbitrary units
# ********************************************************************
function GetOutlierThres( slice, prob )
# Get bin width
low = minimum( slice )
upp = maximum( slice )
# Compute histogram's bin count & separation
if ( upp - low + 1 ) >= 100
bins = convert( Int32, 100 )
step = ( upp + (upp-low)/(bins-1) - low ) / bins
else
bins = convert( Int32, ceil( upp ) - floor( low ) + 1 )
step = 1
end
# Compute histogram
histCount = zeros( bins, 1 )
for k in slice
index = convert( Int64, floor( ( k - low ) / step ) + 1 )
histCount[ index ] += 1
end
# Get bin index, that accounts for desired probability
prob = 0.95
delta = cumsum( histCount, dims=1 ) / sum( histCount ) .- prob
index = findmin( broadcast( abs, delta ) )[2][1]
# Find max intensity image value
key = low + step * index
upp = -1
for k in slice[:]
if key > k && k > upp
upp = k
end
end
return [ low, upp ]
end
# ********************************************************************
# Discretize image's intensity using a given intensity range
# slice: Image matrix in measured arbitrary units
# bounds: max and min intensity discretizing range
# depth: number of intensity discrete steps in final image
# ********************************************************************
function SetPixelDepth( slice, bounds, depth )
# Compute intensity bins
bins = depth - 1
limits = collect(
range( bounds[1], stop = bounds[2], length = depth ) )
# Reserve memory for output image
imgBytes = zeros( UInt8, size( slice ) )
nPixels = length( slice )
step = limits[2] - limits[1]
# Set intensity depth
for k in 1:nPixels
# Find bin
i::Int64 = floor( ( slice[k] - limits[1] ) / step ) + 1
if i < bins
imgBytes[k] = i + ( slice[k] > limits[i] ) - 1
else
imgBytes[k] = bins
end
end
return imgBytes
end
# ********************************************************************
# Discretize image's intensity removing intensity outliers
# slice: Image matrix in measured arbitrary units
# depth: Number of intensity discrete steps in final image
# prob: Proportion of pixesl to take into account in discretization
# ********************************************************************
"""
TrIQ( slice, depth, prob = 0.98 )
TrIQ intensity quantizer as described in _DOI 10.7717/peerj-cs.585_.
Discretize the image?s continuos intensity range in discrete bins
from 0 to 255 gray levels, groupping intensity outiers in the
highehst bin. Returns a UInt8 bidimensional matrix
# Arguments
* `slice`: slice returned by the `GetSlice` function
* `depth`: Number ob discrete bins for intensity quantization
* ` prob`: Cummulative distribute function cutting value
# Examples
```julia
# Save DESI MSI Carcinoma 885.55 mz slice as a bitmap file
spectra = LoadImzml( "80TopL, 50TopR, 70BottomL, 60BottomR-centroid.imzML" )
slice = GetSlice( spectra, 885.55, 0.005 )
SaveBitmap( "TrIQ.bmp", TrIQ( slice, 256, 0.95 ), ViridisPalette )
```
The following image will be created on your hard disk
![](./TrIQ.bmp)
"""
function TrIQ( slice, depth, prob )
return SetPixelDepth(
slice,
GetOutlierThres( slice, prob ),
depth
)
end
# ********************************************************************
# Viridis color palette 256 color levels
# ********************************************************************
ViridisPalette = UInt32[
0x440154, 0x440256, 0x450457, 0x450559,
0x46075A, 0x46085C, 0x460A5D, 0x460B5E,
0x470D60, 0x470E61, 0x471063, 0x471164,
0x471365, 0x481467, 0x481668, 0x481769,
0x48186A, 0x481A6C, 0x481B6D, 0x481C6E,
0x481D6F, 0x481F70, 0x482071, 0x482173,
0x482374, 0x482475, 0x482576, 0x482677,
0x482878, 0x482979, 0x472A7A, 0x472C7A,
0x472D7B, 0x472E7C, 0x472F7D, 0x46307E,
0x46327E, 0x46337F, 0x463480, 0x453581,
0x453781, 0x453882, 0x443983, 0x443A83,
0x443B84, 0x433D84, 0x433E85, 0x423F85,
0x424086, 0x424186, 0x414287, 0x414487,
0x404588, 0x404688, 0x3F4788, 0x3F4889,
0x3E4989, 0x3E4A89, 0x3E4C8A, 0x3D4D8A,
0x3D4E8A, 0x3C4F8A, 0x3C508B, 0x3B518B,
0x3B528B, 0x3A538B, 0x3A548C, 0x39558C,
0x39568C, 0x38588C, 0x38598C, 0x375A8C,
0x375B8D, 0x365C8D, 0x365D8D, 0x355E8D,
0x355F8D, 0x34608D, 0x34618D, 0x33628D,
0x33638D, 0x32648E, 0x32658E, 0x31668E,
0x31678E, 0x31688E, 0x30698E, 0x306A8E,
0x2F6B8E, 0x2F6C8E, 0x2E6D8E, 0x2E6E8E,
0x2E6F8E, 0x2D708E, 0x2D718E, 0x2C718E,
0x2C728E, 0x2C738E, 0x2B748E, 0x2B758E,
0x2A768E, 0x2A778E, 0x2A788E, 0x29798E,
0x297A8E, 0x297B8E, 0x287C8E, 0x287D8E,
0x277E8E, 0x277F8E, 0x27808E, 0x26818E,
0x26828E, 0x26828E, 0x25838E, 0x25848E,
0x25858E, 0x24868E, 0x24878E, 0x23888E,
0x23898E, 0x238A8D, 0x228B8D, 0x228C8D,
0x228D8D, 0x218E8D, 0x218F8D, 0x21908D,
0x21918C, 0x20928C, 0x20928C, 0x20938C,
0x1F948C, 0x1F958B, 0x1F968B, 0x1F978B,
0x1F988B, 0x1F998A, 0x1F9A8A, 0x1E9B8A,
0x1E9C89, 0x1E9D89, 0x1F9E89, 0x1F9F88,
0x1FA088, 0x1FA188, 0x1FA187, 0x1FA287,
0x20A386, 0x20A486, 0x21A585, 0x21A685,
0x22A785, 0x22A884, 0x23A983, 0x24AA83,
0x25AB82, 0x25AC82, 0x26AD81, 0x27AD81,
0x28AE80, 0x29AF7F, 0x2AB07F, 0x2CB17E,
0x2DB27D, 0x2EB37C, 0x2FB47C, 0x31B57B,
0x32B67A, 0x34B679, 0x35B779, 0x37B878,
0x38B977, 0x3ABA76, 0x3BBB75, 0x3DBC74,
0x3FBC73, 0x40BD72, 0x42BE71, 0x44BF70,
0x46C06F, 0x48C16E, 0x4AC16D, 0x4CC26C,
0x4EC36B, 0x50C46A, 0x52C569, 0x54C568,
0x56C667, 0x58C765, 0x5AC864, 0x5CC863,
0x5EC962, 0x60CA60, 0x63CB5F, 0x65CB5E,
0x67CC5C, 0x69CD5B, 0x6CCD5A, 0x6ECE58,
0x70CF57, 0x73D056, 0x75D054, 0x77D153,
0x7AD151, 0x7CD250, 0x7FD34E, 0x81D34D,
0x84D44B, 0x86D549, 0x89D548, 0x8BD646,
0x8ED645, 0x90D743, 0x93D741, 0x95D840,
0x98D83E, 0x9BD93C, 0x9DD93B, 0xA0DA39,
0xA2DA37, 0xA5DB36, 0xA8DB34, 0xAADC32,
0xADDC30, 0xB0DD2F, 0xB2DD2D, 0xB5DE2B,
0xB8DE29, 0xBADE28, 0xBDDF26, 0xC0DF25,
0xC2DF23, 0xC5E021, 0xC8E020, 0xCAE11F,
0xCDE11D, 0xD0E11C, 0xD2E21B, 0xD5E21A,
0xD8E219, 0xDAE319, 0xDDE318, 0xDFE318,
0xE2E418, 0xE5E419, 0xE7E419, 0xEAE51A,
0xECE51B, 0xEFE51C, 0xF1E51D, 0xF4E61E,
0xF6E620, 0xF8E621, 0xFBE723, 0xFDE725 ]

94
src/Common.jl
View File

@ -1,94 +0,0 @@
# *******************************************************************
# VectorConfig, structure with axis decoding instructions
# *******************************************************************
mutable struct SpecDim
Format::DataType # Float32, Float64
Packed::Bool # 0: none 1: zlib compression
Axis::Int32 # 1: m/z 2: Amplitude
Skip::Int64 # Bytes to skip
end
# *******************************************************************
# Fill VectorConfig fields from "cvParam" tags
# *******************************************************************
function ConfigureSpecDim( stream )
# Vector configuration
axis = SpecDim( Float64, false, 1, 0 )
# Initial values
offset = position( stream )
currLine = ""
matchInfo = RegexMatch
while true
# Next field
currLine = readline( stream )
matchInfo = match( r"^\s*<(cvParam)", currLine )
if matchInfo === nothing
matchInfo = match( r"^\s*", currLine )
axis.Skip = position( stream ) - offset -
length( currLine ) + length( matchInfo.match )
return axis
end
index = length( matchInfo.captures[1] )
matchInfo = GetAttribute( currLine[index:end], "accession" )
if matchInfo.captures[1] == "MS:1000515" # intensity array
axis.Axis = 2
continue
elseif matchInfo.captures[1] == "MS:1000519" # 32-bit integer
axis.Format = Int32
continue
elseif matchInfo.captures[1] == "MS:1000521" # 32-bit float
axis.Format = Float32
continue
elseif matchInfo.captures[1] == "MS:1000522" # 64-bit integer
axis.Format = Int64
continue
elseif matchInfo.captures[1] == "MS:1000574" # zlib compresion
axis.Packed = true
continue
end
end
end
# *******************************************************************
# Read lines in file until regex match
# stream: Source file stream
# regex: target tag e.g. "binaryDataArray"
# *******************************************************************
function FindTag( stream, regex )
while true
isTag = match( regex, readline( stream ) )
if isTag !== nothing
return isTag
end
end
end
# *******************************************************************
# Retrieve attribute/value pair
# source: String with atr="value" pairs
# tag: Optional attribute, e.g. "\\spectrum"
# *******************************************************************
function GetAttribute( source, tag = "([^=]+)" )
# Build regex and retrieve attribute/value pair
regStr = Regex( "\\s" * tag * "=\"([^\"]*)\"" )
return( match( regStr, source ) )
end

82
src/DataAccess.jl
View File

@ -1,82 +0,0 @@
# src/DataAccess.jl
"""
load_spectra(fileName::String)
Eagerly loads all spectra from a .mzML or .imzML file into a 2xN matrix.
This function follows the eager-loading pattern of the original `LoadMzml`
for backward compatibility and for analyses that require all data in memory.
Each column in the returned matrix represents a single spectrum:
- Row 1: m/z array
- Row 2: Intensity array
For .imzML files, the (x,y) spatial coordinates are discarded.
# Arguments
* `fileName`: Full path to the .mzML or .imzML file.
# Returns
- An `Array{Any, 2}` containing all spectra.
"""
function load_spectra(fileName::String)
if endswith(lowercase(fileName), ".mzml")
# Detected .mzML file. Using the existing eager loader.
return LoadMzml(fileName)
elseif endswith(lowercase(fileName), ".imzml")
# Detected .imzML file. Eagerly loading all spectra from .ibd.
# Lazily load metadata to get offsets and counts
imzml_data = load_imzml(fileName)
try
num_spectra = length(imzml_data.spectra_metadata)
if num_spectra == 0
return Array{Any}(undef, (2, 0))
end
spectra_matrix = Array{Any}(undef, (2, num_spectra))
hIbd = imzml_data.ibd_handle
mz_format = imzml_data.mz_format
int_format = imzml_data.intensity_format
# For performance, pre-allocate one large buffer for reading.
# This avoids re-allocating memory for each spectrum inside the loop.
max_points = maximum(meta.mz_count for meta in imzml_data.spectra_metadata)
mz_array_buffer = Array{mz_format}(undef, max_points)
intensity_array_buffer = Array{int_format}(undef, max_points)
for i in 1:num_spectra
meta = imzml_data.spectra_metadata[i]
# Create a view into the buffer with the correct size for this spectrum
current_mz_array = view(mz_array_buffer, 1:meta.mz_count)
current_int_array = view(intensity_array_buffer, 1:meta.intensity_count)
# Read binary data directly into the sized views
seek(hIbd, meta.mz_offset)
read!(hIbd, current_mz_array)
seek(hIbd, meta.intensity_offset)
read!(hIbd, current_int_array)
# Store a copy in the final matrix. A copy is necessary because
# the buffer is overwritten in the next iteration.
spectra_matrix[1, i] = copy(current_mz_array)
spectra_matrix[2, i] = copy(current_int_array)
end
return spectra_matrix
finally
# Ensure the .ibd file handle is closed, as this is an eager load.
if isopen(imzml_data.ibd_handle)
close(imzml_data.ibd_handle)
end
end
else
error("Unsupported file type. Please provide a .mzML or .imzML file.")
end
end

369
src/MSIData.jl Normal file
View File

@ -0,0 +1,369 @@
# src/MSIData.jl
"""
This file defines the unified MSIData object and the associated data access layer,
including caching and iteration logic, for handling large mzML and imzML datasets
efficiently.
"""
using Base64, Libz # For reading binary data
# Abstract type for different data sources (e.g., mzML, imzML)
# This allows dispatching to the correct binary reading logic.
abstract type MSDataSource end
# Concrete type for imzML data sources
struct ImzMLSource <: MSDataSource
ibd_handle::IO
mz_format::Type
intensity_format::Type
end
# Concrete type for mzML data sources
struct MzMLSource <: MSDataSource
file_handle::IO
mz_format::Type
intensity_format::Type
end
# Struct to hold info about a binary data array (mz or intensity)
struct SpectrumAsset
format::Type
is_compressed::Bool
offset::Int64
encoded_length::Int32
# For mzML, axis_type is needed to distinguish mz from intensity.
# For imzML, this can be ignored as the order is fixed.
axis_type::Symbol
end
# Metadata for a single spectrum, common to both formats
struct SpectrumMetadata
# For imzML
x::Int32
y::Int32
# For mzML
id::String
# Common binary data info
mz_asset::SpectrumAsset
int_asset::SpectrumAsset
end
# The main unified data object
mutable struct MSIData
source::MSDataSource
spectra_metadata::Vector{SpectrumMetadata}
image_dims::Tuple{Int, Int} # (width, height) for imaging data
coordinate_map::Union{Matrix{Int}, Nothing} # Maps (x,y) to linear index for imzML
# LRU Cache implementation
cache::Dict{Int, Tuple{Vector, Vector}}
cache_order::Vector{Int} # Stores indices, with most recently used at the end
cache_size::Int # Max number of spectra in cache
function MSIData(source, metadata, dims, coordinate_map, cache_size)
obj = new(source, metadata, dims, coordinate_map, Dict(), [], cache_size)
# Ensure file handles are closed when the object is garbage collected
finalizer(obj) do o
if o.source isa ImzMLSource && isopen(o.source.ibd_handle)
close(o.source.ibd_handle)
elseif o.source isa MzMLSource && isopen(o.source.file_handle)
close(o.source.file_handle)
end
end
return obj
end
end
# --- Internal function for reading binary data --- #
function read_binary_vector(io::IO, asset::SpectrumAsset)
seek(io, asset.offset)
raw_b64 = read(io, asset.encoded_length)
decoded_bytes = base64decode(strip(String(raw_b64)))
bytes_io = IOBuffer(asset.is_compressed ? Libz.inflate(decoded_bytes) : decoded_bytes)
n_elements = Int(bytes_io.size / sizeof(asset.format))
out_array = Array{asset.format}(undef, n_elements)
read!(bytes_io, out_array)
return out_array
end
# Overload for different source types
function read_spectrum_from_disk(source::ImzMLSource, meta::SpectrumMetadata)
# For imzML, the binary data is raw, not base64 encoded.
# The `encoded_length` field in this case holds the number of points.
mz = Array{source.mz_format}(undef, meta.mz_asset.encoded_length)
intensity = Array{source.intensity_format}(undef, meta.int_asset.encoded_length)
seek(source.ibd_handle, meta.mz_asset.offset)
read!(source.ibd_handle, mz)
seek(source.ibd_handle, meta.int_asset.offset)
read!(source.ibd_handle, intensity)
return mz, intensity
end
function read_spectrum_from_disk(source::MzMLSource, meta::SpectrumMetadata)
# For mzML, data is Base64 encoded within the XML
mz = read_binary_vector(source.file_handle, meta.mz_asset)
intensity = read_binary_vector(source.file_handle, meta.int_asset)
return mz, intensity
end
# --- Public API --- #
"""
GetSpectrum(data::MSIData, index::Int)
Retrieves a single spectrum by its index, utilizing a cache for performance.
This function is the core of the "Indexed" and "Cache" access patterns.
"""
function GetSpectrum(data::MSIData, index::Int)
if index < 1 || index > length(data.spectra_metadata)
error("Spectrum index $index out of bounds.")
end
# Phase 1: Check the cache
if haskey(data.cache, index)
# Cache Hit: Move item to the end of the LRU list and return from cache
filter!(x -> x != index, data.cache_order)
push!(data.cache_order, index)
return data.cache[index]
end
# Phase 2: Cache Miss - Read from disk
meta = data.spectra_metadata[index]
spectrum = read_spectrum_from_disk(data.source, meta)
# Phase 3: Update cache
if data.cache_size > 0
if length(data.cache) >= data.cache_size
# Evict the least recently used item (at the front of the list)
lru_index = popfirst!(data.cache_order)
delete!(data.cache, lru_index)
end
data.cache[index] = spectrum
push!(data.cache_order, index)
end
return spectrum
end
"""
GetSpectrum(data::MSIData, x::Int, y::Int)
Retrieves a single spectrum by its (x, y) coordinates for imaging data.
Utilizes a coordinate map for efficient lookup and then the cache.
"""
function GetSpectrum(data::MSIData, x::Int, y::Int)
if data.coordinate_map === nothing
error("Coordinate map not available. This method is only for imaging data loaded from .imzML files.")
end
width, height = data.image_dims
if x < 1 || x > width || y < 1 || y > height
error("Coordinates ($x, $y) out of bounds for image dimensions ($width, $height).")
end
index = data.coordinate_map[x, y]
if index == 0
error("No spectrum found at coordinates ($x, $y).")
end
return GetSpectrum(data, index) # Call the existing method
end
"""
get_total_spectrum(msi_data::MSIData; num_bins::Int=20000) -> Tuple{Vector{Float64}, Vector{Float64}}
Calculates the sum of all spectra in the dataset by binning.
The m/z range is determined dynamically by finding the global min/max m/z values
across the entire dataset.
Returns a tuple containing two vectors: the binned m/z axis and the summed intensities.
"""
function get_total_spectrum(msi_data::MSIData; num_bins::Int=20000)
println("Calculating total spectrum...")
# 1. First Pass: Find the global m/z range
println(" Pass 1: Finding global m/z range...")
global_min_mz = Inf
global_max_mz = -Inf
_iterate_spectra_fast(msi_data) do idx, mz, _
if !isempty(mz)
local_min, local_max = extrema(mz)
global_min_mz = min(global_min_mz, local_min)
global_max_mz = max(global_max_mz, local_max)
end
end
if !isfinite(global_min_mz)
@warn "Could not determine a valid m/z range. All spectra might be empty."
return (Float64[], Float64[])
end
println(" Global m/z range found: [$(global_min_mz), $(global_max_mz)]")
# 2. Define Bins
mz_bins = range(global_min_mz, stop=global_max_mz, length=num_bins)
intensity_sum = zeros(Float64, num_bins)
bin_step = step(mz_bins)
# 3. Second Pass: Sum intensities into bins
println(" Pass 2: Summing intensities into $num_bins bins...")
_iterate_spectra_fast(msi_data) do idx, mz, intensity
if isempty(mz)
return # continue equivalent
end
for i in eachindex(mz)
# Find the correct bin for the current m/z value
bin_index = clamp(round(Int, (mz[i] - global_min_mz) / bin_step) + 1, 1, num_bins)
intensity_sum[bin_index] += intensity[i]
end
end
println("Total spectrum calculation complete.")
return (collect(mz_bins), intensity_sum)
end
"""
get_average_spectrum(msi_data::MSIData; num_bins::Int=20000) -> Tuple{Vector{Float64}, Vector{Float64}}
Calculates the average of all spectra in the dataset by binning.
This is effectively the total ion chromatogram (TIC) divided by the number of spectra.
Returns a tuple containing two vectors: the binned m/z axis and the averaged intensities.
"""
function get_average_spectrum(msi_data::MSIData; num_bins::Int=20000)
# This function uses the exact same logic as get_total_spectrum...
mz_bins, intensity_sum = get_total_spectrum(msi_data, num_bins=num_bins)
if isempty(intensity_sum)
return (mz_bins, intensity_sum)
end
# ...with one final step: dividing by the number of spectra to get the average.
println("Averaging spectrum...")
num_spectra = length(msi_data.spectra_metadata)
average_intensity = intensity_sum ./ num_spectra
return (mz_bins, average_intensity)
end
# --- Iterator Implementation --- #
struct MSIDataIterator
data::MSIData
end
"""
IterateSpectra(data::MSIData)
Returns an iterator that yields each spectrum, processing the file sequentially
with minimal memory overhead.
This function is the core of the "Event-driven" access pattern.
"""
function IterateSpectra(data::MSIData)
return MSIDataIterator(data)
end
# Define the iteration interface for the custom iterator
Base.length(it::MSIDataIterator) = length(it.data.spectra_metadata)
Base.eltype(it::MSIDataIterator) = Tuple{Int, Tuple{Vector, Vector}} # Index, (mz, intensity)
function Base.iterate(it::MSIDataIterator, state=1)
if state > length(it.data.spectra_metadata)
return nothing # End of iteration
end
# GetSpectrum uses the cache, so iteration is also cached
spectrum = GetSpectrum(it.data, state)
# Yield the index and the spectrum
return ((state, spectrum), state + 1)
end
# --- High-performance Internal Iterator --- #
"""
_iterate_spectra_fast(f::Function, data::MSIData)
An internal, high-performance iterator for bulk processing.
It avoids the overhead of `GetSpectrum` by reading directly from the file stream
and reusing pre-allocated buffers.
- `f`: A function to call for each spectrum, with signature `f(index, mz, intensity)`.
- `data`: The `MSIData` object.
"""
function _iterate_spectra_fast_impl(f::Function, data::MSIData, source::ImzMLSource)
# This implementation is for imzML and is optimized for performance by
# using pre-allocated buffers.
max_points = 0
for meta in data.spectra_metadata
# For imzML, encoded_length is the number of points
max_points = max(max_points, meta.mz_asset.encoded_length)
end
# If all spectra are empty, just call f with empty arrays
if max_points == 0 && !isempty(data.spectra_metadata)
for i in 1:length(data.spectra_metadata)
f(i, source.mz_format[], source.intensity_format[])
end
return
end
mz_buffer = Vector{source.mz_format}(undef, max_points)
int_buffer = Vector{source.intensity_format}(undef, max_points)
for i in 1:length(data.spectra_metadata)
meta = data.spectra_metadata[i]
nPoints = meta.mz_asset.encoded_length
if nPoints == 0
f(i, source.mz_format[], source.intensity_format[])
continue
end
mz_view = view(mz_buffer, 1:nPoints)
int_view = view(int_buffer, 1:nPoints)
seek(source.ibd_handle, meta.mz_asset.offset)
read!(source.ibd_handle, mz_view)
seek(source.ibd_handle, meta.int_asset.offset)
read!(source.ibd_handle, int_view)
f(i, mz_view, int_view)
end
end
function _iterate_spectra_fast_impl(f::Function, data::MSIData, source::MzMLSource)
# This implementation is for mzML. It iterates through each spectrum,
# decodes it, and then calls the function `f`. It's less performant than
# the imzML version because we cannot pre-allocate buffers of a known size,
# but it is correct and still faster than using GetSpectrum due to bypassing the cache.
for i in 1:length(data.spectra_metadata)
meta = data.spectra_metadata[i]
mz = read_binary_vector(source.file_handle, meta.mz_asset)
intensity = read_binary_vector(source.file_handle, meta.int_asset)
f(i, mz, intensity)
end
end
function _iterate_spectra_fast(f::Function, data::MSIData)
# Dispatch to the correct implementation based on the source type
_iterate_spectra_fast_impl(f, data, data.source)
end

33
src/MSI_src.jl Normal file
View File

@ -0,0 +1,33 @@
module MSI_src
# Export the public API
export OpenMSIData, GetSpectrum, IterateSpectra, ImportMzmlFile, load_slices, plot_slices, plot_slice, get_total_spectrum, get_average_spectrum, LoadMzml, LoadSpectra
# Include all source files directly into the main module
include("MSIData.jl")
include("ParserHelpers.jl")
include("mzML.jl")
include("imzML.jl")
include("MzmlConverter.jl")
# --- 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=100)
if endswith(lowercase(filepath), ".mzml")
return load_mzml_lazy(filepath, cache_size=cache_size)
elseif endswith(lowercase(filepath), ".imzml")
return load_imzml_lazy(filepath, cache_size=cache_size)
else
error("Unsupported file type: $filepath. Please provide a .mzML or .imzML file.")
end
end
end # module MSI_src

605
src/MzmlConverter.jl
View File

@ -1,16 +1,14 @@
# 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 project's original R scripts that use MALDIquant.
It replicates the functionality of the original R scripts that use MALDIquant.
"""
using CSV, DataFrames
using DataFrames, Printf, CSV
# 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
# Note: This file assumes that `LoadSpectra` from DataAccess.jl and the necessary
# structs/functions from ParserHelpers.jl, mzML.jl, and imzML.jl are available
# in the execution context (e.g., included in the main app.jl).
# A struct to hold the processed pixel data before exporting
struct ProcessedPixel
@ -19,11 +17,18 @@ struct ProcessedPixel
intensity::Vector{Float32}
end
# This struct will hold the metadata needed for the XML file
struct BinaryMetadata
mz_offset::UInt64
mz_length::UInt64
int_offset::UInt64
int_length::UInt64
end
"""
GetMzmlScanTime(fileName::String)
Parses a .mzML file to extract the scan start time for each spectrum.
This is a Julia implementation of the `GetMzmlScanTime` function from the R scripts.
# Arguments
* `fileName`: Path to the .mzML file.
@ -32,48 +37,75 @@ This is a Julia implementation of the `GetMzmlScanTime` function from the R scri
- A `Matrix{Int64}` where each row is `[spectrum_index, time_in_milliseconds]`.
"""
function GetMzmlScanTime(fileName::String)
spec_count = 0
times = Tuple{Int64, Int64}[]
# Pre-allocate based on an estimate to reduce re-allocations
try
file_size = filesize(fileName)
estimated_spectra = max(1000, file_size ÷ 10000) # Heuristic
sizehint!(times, estimated_spectra)
catch
# Ignore if filesize fails, sizehint! is just an optimization
end
open(fileName, "r") do stream
file_content = read(stream, String)
state = :outside_spectrum
current_index = 0
current_time = nothing
# First, find the total number of spectra to pre-allocate
spectrum_list_match = match(r"<spectrumList count=\"(\d+)\"", file_content)
if spectrum_list_match === nothing
error("Could not find <spectrumList> tag or count attribute.")
for line in eachline(stream)
if state == :outside_spectrum
if occursin("<spectrum ", line)
state = :in_spectrum
current_index = 0
current_time = nothing
# Process index from the start tag itself
idx_match = match(r"index=\"(\d+)\"", line)
if idx_match !== nothing
current_index = parse(Int, idx_match.captures[1]) + 1
end
total_spectra = parse(Int, spectrum_list_match.captures[1])
sizehint!(times, total_spectra)
# Iterate over each spectrum block
for spectrum_match in eachmatch(r"<spectrum index=\"(\d+)\"[^>]*>.*?<\/spectrum>", file_content, overlay=true)
spec_block = spectrum_match.match
# The R script appears to treat the 0-based mzML index as 1-based, so we add 1.
spec_index = parse(Int, spectrum_match.captures[1]) + 1
# Find the scan start time within the spectrum block
time_match = match(r"<cvParam accession=\"MS:1000016\".*?value=\"([\d\.]+)\"", spec_block)
end
elseif state == :in_spectrum
# Look for time cvParam
if occursin("<cvParam", line)
if occursin("MS:1000016", line) || occursin("MS:1000015", line)
time_match = match(r"value=\"([\d\.]+)\"", line)
if time_match !== nothing
# The value is in minutes. Convert to milliseconds and round.
time_in_minutes = parse(Float64, time_match.captures[1])
time_in_ms = round(Int64, time_in_minutes * 60000)
push!(times, (spec_index, time_in_ms))
time_val = parse(Float64, time_match.captures[1])
# MS:1000016 is minutes, MS:1000015 is seconds
unit_scale = occursin("MS:1000016", line) ? 60000 : 1000
current_time = round(Int64, time_val * unit_scale)
end
end
end
# Normalize times so the first scan is at t=0
# Check for end of spectrum
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
# Sort by spectrum index to ensure correct order before normalization
sort!(times, by = x -> x[1])
# Normalize times relative to the first scan
if !isempty(times)
sort!(times, by = x -> x[1]) # Ensure scans are sorted by index
first_time = times[1][2]
for i in eachindex(times)
times[i] = (times[i][1], times[i][2] - first_time)
end
end
end # open
# Convert vector of tuples to a matrix
if isempty(times)
return Matrix{Int64}(undef, 0, 2)
end
return permutedims(hcat(collect.(times)...))
end
@ -92,146 +124,493 @@ This is a Julia implementation of the `MatchAcquireTime` function from the R scr
# 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})
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
# Load xyz stage pixel delays from the sync file, skipping the first 2 header lines.
pixel_df = CSV.read(sync_file_path, DataFrame, header=false, skipto=3)
pixel = Matrix(pixel_df)
pixel_matrix = Matrix(pixel_df)
# Normalize pixel times so the first is at t=0
pixel[:, 2] .-= pixel[1, 2]
num_pixels = size(pixel, 1)
num_scans = size(scans, 1)
index_matrix = zeros(Int64, num_pixels, 2)
scan_row = 1
for pixel_row in 1:(num_pixels - 1)
start_scan = scan_row
# Find the scan that ends after the *next* pixel starts
while scan_row <= num_scans && scans[scan_row, 2] < pixel[pixel_row + 1, 2]
scan_row += 1
num_pixels_original = size(pixel_matrix, 1)
if num_pixels_original == 0
return zeros(Int64, 0, 5)
end
index_matrix[pixel_row, 1] = start_scan
index_matrix[pixel_row, 2] = scan_row - 1
# Determine image dimensions and generate coordinates
current_img_width = img_width
current_img_height = img_height
# The next pixel's scans start from the one that crossed the boundary
scan_row = max(1, scan_row - 1)
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
# Assign scans for the last pixel
index_matrix[num_pixels, 1] = scan_row
index_matrix[num_pixels, 2] = num_scans
# 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)
# Handle cases where some pixels have no scans by creating an empty range
for i in 1:num_pixels
if index_matrix[i, 2] < index_matrix[i, 1]
index_matrix[i, 2] = index_matrix[i, 1] - 1
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
return hcat(pixel, index_matrix)
# 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
function RenderPixel(pixel_info, scans, spectra, scan_time_deltas, pixel_time_deltas)
first_scan = pixel_info[3]
last_scan = pixel_info[4]
function RenderPixel(pixel_info, scans, msi_data::MSIData, scan_time_deltas, pixel_time_deltas)
pixel_time = pixel_info[3]
first_scan = pixel_info[4]
last_scan = pixel_info[5]
num_actions = last_scan - first_scan
# Use the m/z array of the first scan as the reference
mz_array = spectra[1, first_scan]
# Get reference m/z array from the first scan involved
mz_array, _ = GetSpectrum(msi_data, first_scan)
new_intensity = zeros(Float32, length(mz_array))
# 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) # Return zero intensity for invalid ranges
end
if num_actions == 0 # Single scan contributes to the pixel
scale = pixel_time_deltas[pixel_info[1]] / scan_time_deltas[first_scan]
new_intensity .+= spectra[2, first_scan] .* scale
_, intensity = GetSpectrum(msi_data, first_scan)
# SAFETY: Ensure positive scaling
scale = max(pixel_time_deltas[first_scan] / scan_time_deltas[first_scan], 0.0f0)
new_intensity .= intensity .* scale
elseif num_actions == 1 # Two partial scans contribute
# First partial scan
scale1 = (scans[last_scan, 2] - pixel_info[2]) / scan_time_deltas[first_scan]
new_intensity .+= spectra[2, first_scan] .* scale1
_, intensity1 = GetSpectrum(msi_data, first_scan)
scale1 = max((scans[last_scan, 2] - pixel_time) / scan_time_deltas[first_scan], 0.0f0)
new_intensity .+= intensity1 .* scale1
# Second partial scan
next_pixel_time = pixel_info[2] + pixel_time_deltas[pixel_info[1]]
scale2 = (next_pixel_time - scans[last_scan, 2]) / scan_time_deltas[last_scan]
new_intensity .+= spectra[2, last_scan] .* scale2
_, intensity2 = GetSpectrum(msi_data, last_scan)
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
elseif num_actions > 1 # Multiple scans contribute
# First partial scan
scale1 = (scans[first_scan + 1, 2] - pixel_info[2]) / scan_time_deltas[first_scan]
new_intensity .+= spectra[2, first_scan] .* scale1
_, intensity1 = GetSpectrum(msi_data, first_scan)
scale1 = max((scans[first_scan + 1, 2] - pixel_time) / scan_time_deltas[first_scan], 0.0f0)
new_intensity .+= intensity1 .* scale1
# Full scans in the middle
for i in (first_scan + 1):(last_scan - 1)
new_intensity .+= spectra[2, i]
_, intensity_middle = GetSpectrum(msi_data, i)
new_intensity .+= intensity_middle
end
# Last partial scan
next_pixel_time = pixel_info[2] + pixel_time_deltas[pixel_info[1]]
scale2 = (next_pixel_time - scans[last_scan, 2]) / scan_time_deltas[last_scan]
new_intensity .+= spectra[2, last_scan] .* scale2
_, intensity2 = GetSpectrum(msi_data, last_scan)
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
# FINAL SAFETY: Clamp any negative values to zero
new_intensity = max.(new_intensity, 0.0f0)
return (mz_array, new_intensity)
end
function ConvertMzmlToImzml(source_file::String, timing_matrix::Matrix{Int64}, scans::Matrix{Int64})
# Get image width from timing info
width = findfirst(i -> timing_matrix[i, 1] - timing_matrix[i-1, 1] != 1, 2:size(timing_matrix, 1))
height = size(timing_matrix, 1) ÷ width
if size(timing_matrix, 1) == 0
return ProcessedPixel[], (0, 0)
end
# Load the full mzML data
spectra = LoadSpectra(source_file)
# Get image dimensions from timing_matrix (assuming columns 1 and 2 are X and Y)
width = maximum(timing_matrix[:, 1])
height = maximum(timing_matrix[:, 2])
# Pre-calculate time deltas
scan_time_deltas = diff(scans[:, 2])
pixel_time_deltas = diff(timing_matrix[:, 2])
# Append a final delta for the last element
push!(scan_time_deltas, scan_time_deltas[end])
push!(pixel_time_deltas, pixel_time_deltas[end])
# Open the mzML file using the new memory-efficient API
msi_data = OpenMSIData(source_file)
# Pre-calculate time deltas robustly
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) # Ensure delta is at least 1
end
scan_time_deltas[end] = max(1, scan_time_deltas[end-1]) # Use previous delta for last scan, ensure at least 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)
# Assumes time is in the 3rd column of the timing_matrix
delta = timing_matrix[i+1, 3] - timing_matrix[i, 3]
pixel_time_deltas[i] = max(1, delta) # Ensure delta is at least 1
end
pixel_time_deltas[end] = max(1, pixel_time_deltas[end-1]) # Use previous delta for last pixel, ensure at least 1
end
processed_pixels = ProcessedPixel[]
sizehint!(processed_pixels, size(timing_matrix, 1))
empty_pixel_count = 0
# DEBUG: Scan coverage analysis
println("DEBUG: Scan coverage analysis:")
covered_pixels_count = 0
for i in 1:size(timing_matrix, 1)
first_scan = timing_matrix[i, 4]
last_scan = timing_matrix[i, 5]
if first_scan <= last_scan && first_scan >= 1 && last_scan <= size(scans, 1)
covered_pixels_count += 1
end
end
println(" Pixels with potential scans: $covered_pixels_count/$(size(timing_matrix, 1))")
println(" Total scans available: $(size(scans, 1))")
for i in 1:size(timing_matrix, 1)
pixel_info = timing_matrix[i, :]
first_scan = pixel_info[3]
last_scan = pixel_info[4]
x = pixel_info[1]
y = pixel_info[2]
first_scan = pixel_info[4]
last_scan = pixel_info[5]
# Skip if there are no scans for this pixel
(first_scan > last_scan) && continue
if first_scan > last_scan || first_scan < 1 || last_scan > size(scans, 1)
empty_pixel_count += 1
push!(processed_pixels, ProcessedPixel((x, y), Float64[], Float32[]))
continue
end
# Calculate coordinates
x = ((pixel_info[1] - 1) % width) + 1
y = fld(pixel_info[1] - 1, width) + 1
mz, intensity = RenderPixel(pixel_info, scans, spectra, scan_time_deltas, pixel_time_deltas)
mz, intensity = RenderPixel(pixel_info, scans, msi_data, scan_time_deltas, pixel_time_deltas)
push!(processed_pixels, ProcessedPixel((x, y), mz, intensity))
end
@info "Found and created $empty_pixel_count empty pixels out of $(size(timing_matrix, 1)) total."
return processed_pixels, (width, height)
end
function ExportImzml(target_file::String, pixels::Vector{ProcessedPixel}, dims::Tuple{Int, Int})
@warn "ExportImzml is not yet implemented. The processed pixel data has been generated but not saved to .imzML/.ibd files."
# TODO: Implement the logic to write the .imzML (XML) and .ibd (binary) files.
# 1. Open target_file.imzML and target_file.ibd for writing.
# 2. Write all m/z and intensity arrays sequentially to the .ibd file, tracking offsets and lengths.
# 3. Write the .imzML XML structure, including:
# - Boilerplate headers
# - <fileDescription>, <referenceableParamGroupList>, <scanSettingsList>
# - A <spectrumList> with a <spectrum> for each pixel.
# - Each <spectrum> must contain cvParams for x/y coords and external data pointers to the .ibd file.
ibd_file = replace(target_file, r"\.imzML$"i => ".ibd")
# Ensure we have valid pixels
if isempty(pixels)
@error "No pixels to export"
return false
end
binary_meta = BinaryMetadata[]
sizehint!(binary_meta, length(pixels))
try
# Step 1: Write .ibd file with proper binary format
open(ibd_file, "w") do ibd_stream
# Write UUID as first 16 bytes (placeholder)
write(ibd_stream, zeros(UInt8, 16))
# Write in little-endian format (standard for imzML)
for pixel in pixels
# Write m/z array as 64-bit little-endian floats
mz_offset = position(ibd_stream)
for val in pixel.mz
write(ibd_stream, htol(Float64(val))) # Host to little-endian
end
mz_length = position(ibd_stream) - mz_offset
# Write intensity array as 32-bit little-endian floats
int_offset = position(ibd_stream)
for val in pixel.intensity
write(ibd_stream, htol(Float32(val))) # Host to little-endian
end
int_length = position(ibd_stream) - int_offset
push!(binary_meta, BinaryMetadata(mz_offset, mz_length, int_offset, int_length))
end
end
# Step 2: Write proper .imzML XML file
open(target_file, "w") do imzml_stream
# Use indexedmzML wrapper
write(imzml_stream, """<?xml version="1.0" encoding="ISO-8859-1"?>
<indexedmzML 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_idx.xsd">
<mzML version="1.1" id="$(splitext(basename(target_file))[1])">
""")
# CV List
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>
""")
# File Description
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"/>
</fileContent>
</fileDescription>
""")
# Referenceable Param Groups
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="IMS" accession="IMS:1000101" name="external data" value="true"/>
<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="IMS" accession="IMS:1000101" name="external data" value="true"/>
<cvParam cvRef="MS" accession="MS:1000521" name="32-bit float"/>
</referenceableParamGroup>
</referenceableParamGroupList>
""")
# Sample List
write(imzml_stream, """ <sampleList count="1">
<sample id="sample1" name="ImagingSample">
<cvParam cvRef="MS" accession="MS:1000001" name="sample number" value="1"/>
</sample>
</sampleList>
""")
# Software List - OPEN SOURCE
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>
""")
# Scan Settings IMAGING
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"/>
</scanSettings>
</scanSettingsList>
""")
# Instrument Configuration MALDI/TOF
write(imzml_stream, """ <instrumentConfigurationList count="1">
<instrumentConfiguration id="instrument1">
<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"/>
</analyzer>
<detector order="3">
<cvParam cvRef="MS" accession="MS:1000253" name="electron multiplier"/>
</detector>
</componentList>
<softwareRef ref="MSIConverter"/>
</instrumentConfiguration>
</instrumentConfigurationList>
""")
# Data Processing
write(imzml_stream, """ <dataProcessingList count="1">
<dataProcessing id="conversionProcessing">
<processingMethod order="1" softwareRef="MSIConverter">
<cvParam cvRef="MS" accession="MS:1000544" name="Conversion to imzML"/>
</processingMethod>
</dataProcessing>
</dataProcessingList>
""")
# Run and Spectrum List
spectrum_offsets = UInt64[]
write(imzml_stream, """ <run defaultInstrumentConfigurationRef="instrument1" id="run1" sampleRef="sample1">
<spectrumList count="$(length(pixels))" defaultDataProcessingRef="conversionProcessing">
""")
# Write each spectrum
for (i, pixel) in enumerate(pixels)
meta = binary_meta[i]
x, y = pixel.coords
spectrum_start = position(imzml_stream)
push!(spectrum_offsets, spectrum_start)
write(imzml_stream, """ <spectrum id="Scan=$(i)" defaultArrayLength="$(length(pixel.mz))" index="$(i-1)">
<cvParam cvRef="MS" accession="MS:1000511" name="ms level" value="1"/>
<cvParam cvRef="MS" accession="MS:1000128" name="profile spectrum"/>
<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:1000102" name="external offset" value="$(meta.mz_offset)"/>
<cvParam cvRef="IMS" accession="IMS:1000103" name="external array length" value="$(length(pixel.mz))"/>
<cvParam cvRef="IMS" accession="IMS:1000104" name="external encoded length" value="$(meta.mz_length)"/>
<binary/>
</binaryDataArray>
<binaryDataArray encodedLength="0">
<referenceableParamGroupRef ref="intensityArray"/>
<cvParam cvRef="IMS" accession="IMS:1000102" name="external offset" value="$(meta.int_offset)"/>
<cvParam cvRef="IMS" accession="IMS:1000103" name="external array length" value="$(length(pixel.intensity))"/>
<cvParam cvRef="IMS" accession="IMS:1000104" name="external encoded length" value="$(meta.int_length)"/>
<binary/>
</binaryDataArray>
</binaryDataArrayList>
</spectrum>
""")
end
write(imzml_stream, """ </spectrumList>
</run>
</mzML>
""")
# Index List
index_list_start = position(imzml_stream)
write(imzml_stream, """ <indexList count="1">
<index name="spectrum">
""")
for (i, offset) in enumerate(spectrum_offsets)
write(imzml_stream, " <offset idRef=\"Scan=$(i)\">$offset</offset>\n")
end
write(imzml_stream, """ </index>
</indexList>
<indexListOffset>$index_list_start</indexListOffset>
</indexedmzML>
""")
end
println("Successfully created: $target_file")
println("Successfully created: $ibd_file")
return true
catch e
@error "Failed to export imzML files" exception=(e, catch_backtrace())
# Clean up partial files
isfile(ibd_file) && rm(ibd_file, force=true)
isfile(target_file) && rm(target_file, force=true)
return false
end
end
"""
ImportMzmlFile(source_file::String, sync_file::String, target_file::String)
@ -241,13 +620,15 @@ Main workflow function to convert a .mzML file to an .imzML file.
* `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)
function ImportMzmlFile(source_file::String, sync_file::String, target_file::String; img_width::Int=0, img_height::Int=0)
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)
timing_matrix = MatchAcquireTime(sync_file, scans; img_width=img_width, img_height=img_height)
println("Step 3: Converting spectra...")
processed_pixels, (width, height) = ConvertMzmlToImzml(source_file, timing_matrix, scans)
@ -264,7 +645,7 @@ function ImportMzmlFile(source_file::String, sync_file::String, target_file::Str
if success
println("Conversion successful: $target_file")
else
println("Conversion failed. Exporting is not yet implemented.")
println("Conversion failed.")
end
return success
end

15
src/ParserHelpers.jl
View File

@ -51,7 +51,7 @@ Retrieves an attribute's value from an XML tag string.
"""
function get_attribute(source::AbstractString, tag::String = "([^=]+)")
# Construct the regex pattern string
pattern_str = "\s" * tag * "=\"([^"]*)\""
pattern_str = "\\s" * tag * "=\"([^\"]*)\""
regStr = Regex(pattern_str)
return match(regStr, source)
end
@ -107,16 +107,3 @@ mutable struct CVParams
is_compressed::Bool
axis_type::Symbol
end
function update_cv_params!(params::CVParams, acc::String)
if acc == "MS:1000514"; params.axis_type = :mz;
elseif acc == "MS:1000515"; params.axis_type = :intensity;
elseif acc == "MS:1000519"; params.format = Int32;
elseif acc == "MS:1000521"; params.format = Float32;
elseif acc == "MS:1000522"; params.format = Int64;
elseif acc == "MS:1000523"; params.format = Float64;
elseif acc == "MS:1000574"; params.is_compressed = true;
end
end

281
src/imzML.jl
View File

@ -1,40 +1,4 @@
using Images, ImageFiltering, StatsBase, Statistics, CairoMakie, ColorSchemes, DataFrames, CSV, Printf, Dates
include("ParserHelpers.jl")
# --- Lazy Loading Data Structures ---
# Struct to hold metadata for a single spectrum, for lazy loading.
struct SpectrumMetadata
x::Int32
y::Int32
mz_offset::Int64
intensity_offset::Int64
mz_count::Int32
intensity_count::Int32
end
# Struct to hold the parsed imzML metadata and the handle to the .ibd file.
mutable struct ImzMLData
ibd_handle::IO
mz_format::Type
intensity_format::Type
spectra_metadata::Vector{SpectrumMetadata}
width::Int
height::Int
mz_is_first::Bool # To know the order in the ibd file
function ImzMLData(ibd_handle::IO, mz_format::Type, intensity_format::Type, spectra_metadata::Vector{SpectrumMetadata}, width::Int, height::Int, mz_is_first::Bool)
obj = new(ibd_handle, mz_format, intensity_format, spectra_metadata, width, height, mz_is_first)
finalizer(obj) do o
if isopen(o.ibd_handle)
close(o.ibd_handle)
end
end
return obj
end
end
using Images, Statistics, CairoMakie, DataFrames, Printf, ColorSchemes
# --- Extracted from imzML.jl ---
@ -43,7 +7,7 @@ This file provides a library for parsing `.imzML` and `.ibd` files in pure Julia
It is intended to be included by a parent script.
Core Functions:
- `load_imzml`: The main function that orchestrates the parsing.
- `load_imzml_lazy`: The main function that orchestrates the parsing.
- Helper functions for reading XML metadata and binary spectral data.
"""
@ -174,70 +138,66 @@ end
"""
load_imzml(file_path::String)
Main function to parse an `.imzML`/.ibd file pair and prepare for lazy loading.
It parses the metadata from the `.imzML` file, including spectrum coordinates and
byte offsets for data in the `.ibd` file, but does not load the spectral data itself.
load_imzml_lazy(file_path::String; cache_size=100)
# Returns
- An `ImzMLData` struct containing the open `.ibd` file handle and all necessary metadata
for on-demand data loading.
Main function to parse an `.imzML`/.ibd file pair and prepare for lazy loading.
It now returns a unified MSIData object.
"""
function load_imzml(file_path::String)
function load_imzml_lazy(file_path::String; cache_size=100)
println("DEBUG: Checking for .imzML file at $file_path")
if !isfile(file_path)
error("Provided path is not a file: $(file_path)")
end
ibd_path = replace(file_path, ".imzML" => ".ibd")
ibd_path = replace(file_path, r"\.(imzML|mzML)"i => ".ibd")
println("DEBUG: Checking for .ibd file at $ibd_path")
if !isfile(ibd_path)
error("Corresponding .ibd file not found for: $(file_path)")
end
stream = open(file_path)
hIbd = open(ibd_path)
# We don't use a try/finally block because we need to return the open hIbd handle.
# The finalizer on the ImzMLData struct will be responsible for closing it.
println("DEBUG: Opening file streams for .imzML and .ibd")
stream = open(file_path, "r")
hIbd = open(ibd_path, "r")
try
println("DEBUG: Configuring axes...")
axis = axes_config_img(stream)
println("DEBUG: Getting image dimensions...")
imgDim = get_img_dimensions(stream)
width, height, num_spectra = imgDim
println("DEBUG: Image dimensions: $(width)x$(height), $num_spectra spectra.")
# Determine which format is for m/z and which for intensity from the referenceableParamGroups
mz_config_idx = findfirst(a -> a.Axis == 1, axis)
int_config_idx = findfirst(a -> a.Axis == 2, axis)
mz_format = axis[mz_config_idx].Format
intensity_format = axis[int_config_idx].Format
println("DEBUG: m/z format: $mz_format, Intensity format: $intensity_format")
# The original parsing logic is highly optimized and relies on fixed offsets within the XML structure.
# We reuse it to get the initial .ibd file offset and the structural attributes of the <spectrum> tags.
# --- NEW PARSING LOGIC based on the old, working code ---
println("DEBUG: Learning file structure from first spectrum...")
start_of_spectra_xml = position(stream)
attr = get_spectrum_attributes(stream, hIbd)
# get_spectrum_attributes moves the hIbd pointer to the start of the first spectrum's binary data.
current_ibd_offset = position(hIbd)
# We need to re-scan the spectra list from the beginning to get metadata for each spectrum.
seek(stream, start_of_spectra_xml)
println("DEBUG: Initial IBD offset: $current_ibd_offset")
spectra_metadata = Vector{SpectrumMetadata}(undef, num_spectra)
# attr[3] is 4 if the first binary data array is intensity, 3 if it's m/z.
mz_is_first = attr[3] == 3
println("DEBUG: Parsing metadata for $num_spectra spectra using skip-based method...")
for k in 1:num_spectra
# These skips are based on the structure of the first spectrum tag, assuming all are identical.
# Use skip values learned from the first spectrum, assuming all are identical.
skip(stream, attr[5]) # Skip to X coordinate value
val_tag = find_tag(stream, r"value=\"(\d+)\"")
x = parse(Int32, val_tag.captures[1])
val_tag_x = find_tag(stream, r"value=\"(\d+)\"")
x = parse(Int32, val_tag_x.captures[1])
skip(stream, attr[6]) # Skip to Y coordinate value
val_tag = find_tag(stream, r"value=\"(\d+)\"")
y = parse(Int32, val_tag.captures[1])
val_tag_y = find_tag(stream, r"value=\"(\d+)\"")
y = parse(Int32, val_tag_y.captures[1])
skip(stream, attr[7]) # Skip to array length value
val_tag = find_tag(stream, r"value=\"(\d+)\"")
nPoints = parse(Int32, val_tag.captures[1])
val_tag_len = find_tag(stream, r"value=\"(\d+)\"")
nPoints = parse(Int32, val_tag_len.captures[1])
mz_len_bytes = nPoints * sizeof(mz_format)
int_len_bytes = nPoints * sizeof(intensity_format)
@ -251,18 +211,42 @@ function load_imzml(file_path::String)
mz_offset = int_offset + int_len_bytes
end
spectra_metadata[k] = SpectrumMetadata(x, y, mz_offset, int_offset, nPoints, nPoints)
# Create modern SpectrumAsset objects
mz_asset = SpectrumAsset(mz_format, false, mz_offset, nPoints, :mz)
int_asset = SpectrumAsset(intensity_format, false, int_offset, nPoints, :intensity)
spectra_metadata[k] = SpectrumMetadata(x, y, "", mz_asset, int_asset)
# Advance the offset for the next spectrum's data.
current_ibd_offset += mz_len_bytes + int_len_bytes
skip(stream, attr[8]) # Skip to the end of the spectrum tag
end
# --- END OF NEW PARSING LOGIC ---
println("DEBUG: Metadata parsing complete.")
# Build coordinate map for imzML files
println("DEBUG: Building coordinate map...")
coordinate_map = zeros(Int, width, height)
for (idx, meta) in enumerate(spectra_metadata)
if 1 <= meta.x <= width && 1 <= meta.y <= height
coordinate_map[meta.x, meta.y] = idx
end
end
println("DEBUG: Coordinate map built.")
# Close the .imzML file stream, but leave the .ibd stream open for lazy reading.
close(stream)
return ImzMLData(hIbd, mz_format, intensity_format, spectra_metadata, width, height, mz_is_first)
source = ImzMLSource(hIbd, mz_format, intensity_format)
println("DEBUG: Creating MSIData object.")
return MSIData(source, spectra_metadata, (width, height), coordinate_map, cache_size)
catch e
close(stream)
close(hIbd)
rethrow(e)
end
end
# --- End of content from imzML.jl ---
@ -310,9 +294,8 @@ end
load_slices(folder, masses, tolerance)
Loads image slices for multiple masses from all `.imzML` files in a directory.
This function is optimized to read the binary data file (`.ibd`) only once per file.
It iterates through each spectrum, reads its data, and finds all target masses
before proceeding to the next spectrum, minimizing disk I/O.
This function is now refactored to use the new MSIData architecture and its
caching capabilities.
"""
function load_slices(folder, masses, tolerance)
files = filter(f -> endswith(f, ".imzML"), readdir(folder, join=true))
@ -331,43 +314,27 @@ function load_slices(folder, masses, tolerance)
push!(names, name)
@info "Processing $(i)/$(n_files): $(name)"
imzML_data = @time load_imzml(file) # This is fast, it's fine.
# Load data using the new lazy loader, returning an MSIData object
msi_data = @time load_imzml_lazy(file)
# Create empty images for all slices for the current file
current_file_slices = [zeros(Float64, imzML_data.width, imzML_data.height) for _ in 1:n_slices]
width, height = msi_data.image_dims
current_file_slices = [zeros(Float64, width, height) for _ in 1:n_slices]
hIbd = imzML_data.ibd_handle
mz_format = imzML_data.mz_format
int_format = imzML_data.intensity_format
# Pre-allocate arrays to be reused for each spectrum
mz_array = Array{mz_format}(undef, 0)
intensity_array = Array{int_format}(undef, 0)
# Iterate through each spectrum ONCE
for meta in imzML_data.spectra_metadata
# Resize arrays if the number of points in the spectrum has changed
if length(mz_array) != meta.mz_count
resize!(mz_array, meta.mz_count)
resize!(intensity_array, meta.intensity_count)
end
# Read spectrum data ONCE from the .ibd file
seek(hIbd, meta.mz_offset)
read!(hIbd, mz_array)
seek(hIbd, meta.intensity_offset)
read!(hIbd, intensity_array)
# Use the high-performance iterator to process all spectra
_iterate_spectra_fast(msi_data) do spec_idx, mz_array, intensity_array
meta = msi_data.spectra_metadata[spec_idx]
# Now, check for all masses of interest in this single spectrum
for (j, mass) in enumerate(masses)
intensity = find_mass(mz_array, intensity_array, mass, tolerance)
if intensity > 0.0
if 1 <= meta.x <= imzML_data.width && 1 <= meta.y <= imzML_data.height
current_file_slices[j][meta.x, meta.y] = intensity
if 1 <= meta.x <= width && 1 <= meta.y <= height
current_file_slices[j][meta.y, meta.x] = intensity
end
end
end
end # end of spectra loop
end # end of fast iterator
# Assign the generated images to the main list
for j in 1:n_slices
@ -379,6 +346,70 @@ function load_slices(folder, masses, tolerance)
return (img_list, names)
end
"""
plot_slice(msi_data::MSIData, mass::Float64, tolerance::Float64, output_dir::String; stage_name="slice", bins=256)
Generates and saves a plot of a single image slice for a given m/z value.
This function closely imitates the logic of the original `GetSlice` but uses
the modern `MSIData` access patterns and robust peak finding.
"""
function plot_slice(msi_data::MSIData, mass::Real, tolerance::Real, output_dir::String; stage_name="slice_mz_$(mass)", bins=256)
# 1. Create an empty image for the slice, with dimensions matching plotting expectations
width, height = msi_data.image_dims
slice_matrix = zeros(Float64, height, width)
# 2. Iterate through each spectrum to build the slice
println("Generating slice for m/z $mass...")
_iterate_spectra_fast(msi_data) do spec_idx, mz_array, intensity_array
meta = msi_data.spectra_metadata[spec_idx]
# Find the peak intensity using the modern, robust find_mass
intensity = find_mass(mz_array, intensity_array, mass, tolerance)
if intensity > 0.0
# Populate the matrix using (y, x) indexing
if 1 <= meta.x <= width && 1 <= meta.y <= height
slice_matrix[meta.y, meta.x] = intensity
end
end
end
println("Slice generation complete.")
# 3. Plot the resulting slice matrix using CairoMakie
println("Plotting slice...")
fig = Figure(size = (600, 500))
ax = CairoMakie.Axis(fig[1, 1],
aspect=DataAspect(),
title=@sprintf("Slice for m/z: %.2f", mass)
)
hidedecorations!(ax)
# Use mass-specific bounds for colorrange, ensuring a valid range
min_val, max_val = extrema(slice_matrix)
if min_val == max_val
max_val = min_val + 1.0 # Ensure the color range has a non-zero width
end
hm = heatmap!(ax, slice_matrix,
colormap=cgrad(ColorSchemes.viridis, bins),
colorrange=(min_val, max_val)
)
Colorbar(fig[1, 2], hm, label="Intensity")
colgap!(fig.layout, 5)
# 4. Save the plot
mkpath(output_dir)
filename = "$(stage_name).png"
save_path = joinpath(output_dir, filename)
save(save_path, fig)
@info "Saved slice plot to $save_path"
return fig
end
# ============================================================================
#
#
@ -647,53 +678,3 @@ function plot_slices(slices, names, masses, output_dir; stage_name, bins=256, dp
end
return fig
end
# =============================================================================
#
#
# Main Workflow
#
# ============================================================================
function load_config(path)
config = Dict{String, Any}()
if !isfile(path)
@warn "Config file not found at $(path). Using default parameters."
# Define defaults here in case the file is missing
config["masses"] = [100]
config["tolerance"] = 0.1
config["color_depth"] = 255
config["outlier_prob"] = 0.98
return config
end
for line in eachline(path)
line = strip(line)
if isempty(line) || startswith(line, "#")
continue # Skip empty lines and comments
end
parts = split(line, '=', limit=2)
if length(parts) != 2
@warn "Skipping malformed line in config: $(line)"
continue
end
key = strip(parts[1])
value_str = strip(parts[2])
try
if key == "masses"
values = [parse(Float64, s) for s in split(value_str, ',')]
config[key] = values
elseif key == "tolerance"
config[key] = parse(Float64, value_str)
elseif key == "color_depth"
config[key] = parse(Int, value_str)
elseif key == "outlier_prob"
config[key] = parse(Float64, value_str)
end
catch e
@error "Could not parse value for key '$(key)': $(value_str)"
end
end
return config
end

246
src/imzML_old.jl
View File

@ -1,246 +0,0 @@
# include( "Common.jl" );
# *******************************************************************
# Load Spectra and return a matrix
# fileName: Full name path
# *******************************************************************
"""
LoadImzml( fileName )
Load an imzML file as a matrix. Each column stores x-pixel position,
y-pixel position, x-axis data and y-axis data.
# Arguments
* `fileName`: Full path name of the imzML file
# Examples
```julia
# Load DESI MSI Carcinoma image data
spectra = LoadImzml( "80TopL, 50TopR, 70BottomL, 60BottomR-centroid.imzML" )
size( spectra )
(4, 18632)
```
"""
function LoadImzml( fileName )
# Open file handles
if !isfile(fileName)
error("provided path is not a file")
end
if endswith(fileName, ".imzML")
stream = open(fileName)
hIbd = open(replace(fileName, ".imzML" => ".ibd"))
else
stream = open(fileName * ".imzML")
hIbd = open(fileName * ".ibd")
end
# Get axes types and image dimensions
axis = AxesConfigImg( stream )
imgDim = GetImgDimensions( stream )
format = [ axis[1].Format, axis[2].Format ]
# Locate spectrum attributes
start = position( stream )
attr = GetSpectrumAttributes( stream, hIbd )
# Load spectra
seek( stream, start )
spectra = LoadImgData( stream, hIbd, attr, imgDim[3], format )
close( stream )
close( hIbd )
return spectra
end
# *******************************************************************
# Get axes value type
# *******************************************************************
function AxesConfigImg( stream )
# Locate which axis is defined at first
tag = FindTag( stream, r"^\s*<(referenceableParamGroup )" )
value = GetAttribute( tag.captures[1], "intensityArray" )
order = 1 + ( value !== nothing )
# Read first axis configuration
axis = Array{ SpecDim, 1 }( undef, 2 )
axis[ order ] = ConfigureSpecDim( stream )
# Read second axis configuration
FindTag( stream, r"^\s*<(referenceableParamGroup )" )
axis[ xor( order,3 ) ] = ConfigureSpecDim( stream )
return axis
end
# *******************************************************************
# Get vector's storage options and image dimensions
# *******************************************************************
function GetImgDimensions( stream )
# Looks for "scanSettings" tag
# FindImgTag( stream, "scanSettings" )
FindTag( stream, r"^\s*<(scanSettings )" )
# Initial values for dimension retrieve
n = 2
dim = [ 0, 0, 0 ]
currLine = ""
matchInfo = RegexMatch
while( n > 0 )
# Next field
currLine = readline( stream )
matchInfo = match( r"^\s*<(cvParam)", currLine )
index = length( matchInfo.captures[1] )
matchInfo = GetAttribute( currLine[index:end], "accession" )
# Get axis identity
if( matchInfo.captures[1] == "IMS:1000042" # max X
|| matchInfo.captures[1] == "IMS:1000043" ) # max Y
# Read dimension's pixels
axis = matchInfo.captures[1][end] - '1'
index += matchInfo.offsets[1] + length( matchInfo.captures[1] )
matchInfo = GetAttribute( currLine[index:end], "value" )
dim[axis] = parse( Int32, matchInfo.captures[1] )
n -= 1
end
end
# Load stored spectra counter
matchInfo = FindTag( stream, r"^\s*<spectrumList(.+)" )
matchInfo = GetAttribute( matchInfo.captures[1], "count" )
dim[ 3 ] = parse( Int32, matchInfo.captures[1] )
return dim
end
# *******************************************************************
# Length of "spectrum" tag without attribute values
# *******************************************************************
function GetSpectrumTag( stream )
offset = position( stream )
tag = FindTag( stream, r"^\s*<spectrum (.+)" )
first = 1
while true
value = match( r"[^=]+=\"([^\"]+)\"", tag.captures[1][first:end] )
if value === nothing
break
end
first += value.offsets[1] + length( value.captures[1] )
offset += length( value.captures[1] )
end
return offset
end
# seek( stream, 7141 )
# *******************************************************************
# Locate spectrum attribute's
# [1] 1: x position stored first 2: y position stored first
# [2] Second stored position
# [3] 1: mzArray stored first 2: intensityArray stored first
# [4] Second stored axis
# [5] skip chars to find first dimension
# [6] skip chars to find second dimension
# [7] skip chars to find vector length
# [8] skip chars to find next spectrum
# *******************************************************************
function GetSpectrumAttributes( stream, hIbd )
# Reserve memory for vector configuration
skip = Vector{ UInt32 }( undef, 8 )
offset = GetSpectrumTag( stream )
# Get axis order and position of pixel coordinate values
tag = FindTag( stream, r" accession=\"IMS:100005(\d)\"(.+)" )
skip[1] = tag.captures[1][1] + 1 - '0'
skip[2] = xor( skip[1], 3 )
value = match( r" value=\"(\d+)\".+", tag.captures[2] )
skip[5] = position( stream ) - offset - length( value.match ) - 2
value = match( r" value=\"\d+\".+", readline( stream ) )
skip[6] = value.offset
# Get axis order
offset = position( stream )
tag = FindTag( stream, r"^\s*<referenceableParamGroupRef(.+)" )
value = GetAttribute( tag.captures[1], "ref" )
skip[3] = ( value.captures[1] == "intensityArray" ) + 3
skip[4] = xor( skip[3], 7 )
k = 2
while k != 0
tag = FindTag( stream, r" accession=\"IMS:100010(\d)\"(.+)" )
# Set IBD offset for first spectrum
if tag.captures[1][1] == '2'
value = GetAttribute( tag.match, "value" )
seek( hIbd, parse( Int32, value.captures[1] ) )
k -= 1
continue
elseif tag.captures[1][1] == '3'
skip[7] = position( stream ) - offset - length( tag.match )
offset = position( stream )
k -= 1
end
end
# Compute characters to skip
FindTag( stream, r"^\s*</spectrum" )
skip[8] = position( stream ) - offset
return skip
end
# *******************************************************************
# Read spectral data
# *******************************************************************
function LoadImgData( stream, hIbd, attr, imgDim, format )
# Reserve spectra memory & update IBD seek offset
spectra = Array{Any}( undef, ( 4, imgDim ) )
contador = 0;
for k in 1:imgDim
# Load image coordinates
skip( stream, attr[5] )
value = FindTag( stream, r"value=\"(\d+)\"" )
spectra[ attr[1],k ] = parse( Int32, value.captures[1] )
skip( stream, attr[6] )
value = FindTag( stream, r"value=\"(\d+)\"" )
spectra[ attr[2],k ] = parse( Int32, value.captures[1] )
# Get vector length
skip( stream, attr[7] )
value = FindTag( stream, r"value=\"(\d+)\"" )
nPoints = parse( Int32, value.captures[1] )
# Reserve memory and read spectrum
spectra[ 3,k ] = Array{ format[1] }( undef, nPoints )
spectra[ 4,k ] = Array{ format[2] }( undef, nPoints )
read!( hIbd, spectra[ attr[3],k ] )
read!( hIbd, spectra[ attr[4],k ] )
skip( stream, attr[8] )
end
return spectra
end

217
src/mzML.jl
View File

@ -1,78 +1,89 @@
# /home/pixel/Documents/Cinvestav_2025/JuliaMSI/src/mzML_lazy.jl
# /home/pixel/Documents/Cinvestav_2025/JuliaMSI/src/mzML.jl
using Base64, Libz
include("ParserHelpers.jl")
# ============================================================================
# Data Structures for Lazy Loading
# ============================================================================
struct SpectrumAsset
format::Type
is_compressed::Bool
offset::Int64
encoded_length::Int32
axis_type::Symbol # :mz or :intensity
end
struct MzMLSpectrumMetadata
id::String
mz_asset::SpectrumAsset
intensity_asset::SpectrumAsset
end
mutable struct MzMLData
file_handle::IO
spectra_metadata::Vector{MzMLSpectrumMetadata}
function MzMLData(file_handle::IO, spectra_metadata::Vector{MzMLSpectrumMetadata})
obj = new(file_handle, spectra_metadata)
finalizer(obj) do o
if isopen(o.file_handle)
close(o.file_handle)
end
end
return obj
end
end
# ============================================================================
# mzML Parser Implementation
# ============================================================================
# This file is responsible for parsing metadata from .mzML files.
# It has been refactored to produce a unified MSIData object.
# This function is kept internal to the mzML parsing process
function get_spectrum_asset_metadata(stream)
start_pos = position(stream)
bda_tag = find_tag(stream, r"<binaryDataArray\s+encodedLength=\"(\d+)\"" )
if bda_tag === nothing; error("Cannot find binaryDataArray"); end
if bda_tag === nothing
error("Cannot find binaryDataArray")
end
encoded_length = parse(Int32, bda_tag.captures[1])
params = CVParams(Float64, false, :mz)
# Initialize parameters as separate variables
data_format = Float64
compression_flag = false
axis = :mz
while true
line = readline(stream)
if occursin("</binaryDataArray>", line); break; end
if occursin("</binaryDataArray>", line)
break
end
if occursin("<cvParam", line)
accession = get_attribute(line, "accession")
if accession === nothing; continue; end
update_cv_params!(params, accession.captures[1])
if accession === nothing
continue
end
acc_str = accession.captures[1]
# println("DEBUG: Processing accession: $acc_str")
# Direct inline logic - NO FUNCTION CALLS
if acc_str == "MS:1000514"
axis = :mz
# println("DEBUG: Set axis_type to :mz")
elseif acc_str == "MS:1000515"
axis = :intensity
# println("DEBUG: Set axis_type to :intensity")
elseif acc_str == "MS:1000518"
data_format = Int16
# println("DEBUG: Set format to Int16")
elseif acc_str == "MS:1000519"
data_format = Int32
# println("DEBUG: Set format to Int32")
elseif acc_str == "MS:1000520"
data_format = Float64
# println("DEBUG: Set format to Float64")
elseif acc_str == "MS:1000521"
data_format = Float32
# println("DEBUG: Set format to Float32")
elseif acc_str == "MS:1000522"
data_format = Int64
# println("DEBUG: Set format to Int64")
elseif acc_str == "MS:1000523"
data_format = Float64
# println("DEBUG: Set format to Float64")
elseif acc_str == "MS:1000574"
compression_flag = true
# println("DEBUG: Set is_compressed to true")
elseif acc_str == "MS:1000576"
compression_flag = false
# println("DEBUG: Set is_compressed to false")
else
# println("DEBUG: Unknown accession: $acc_str")
end
end
end
seek(stream, start_pos)
binary_tag = find_tag(stream, r"<binary>")
if binary_tag === nothing; error("Cannot find binary tag"); end
readuntil(stream, "<binary>")
binary_offset = position(stream)
find_tag(stream, r"</binaryDataArray>") # Position for next call
# Move stream to the end of the binary data array for the next iteration
readuntil(stream, "</binaryDataArray>")
return SpectrumAsset(params.format, params.is_compressed, binary_offset, encoded_length, params.axis_type)
# 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
function parse_spectrum_metadata(stream, offset::Int64)
seek(stream, offset)
id_match = find_tag(stream, r"<spectrum\s+index=\"\d+\"\s+id=\"([^"]+)\"")
id_match = find_tag(stream, r"<spectrum\s+index=\"\d+\"\s+id=\"([^\"]+)" )
id = id_match === nothing ? "" : id_match.captures[1]
asset1 = get_spectrum_asset_metadata(stream)
@ -81,111 +92,113 @@ function parse_spectrum_metadata(stream, offset::Int64)
mz_asset = asset1.axis_type == :mz ? asset1 : asset2
int_asset = asset1.axis_type == :intensity ? asset1 : asset2
return MzMLSpectrumMetadata(id, mz_asset, int_asset)
# Create the new unified metadata object
# For mzML, x and y coordinates are not applicable, so we use 0.
return SpectrumMetadata(0, 0, id, mz_asset, int_asset)
end
function parse_offset_list(stream)
offsets = Int64[]
while true
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)
m = match(r"<offset\s+idRef=\"\S+\">(\d+)</offset>", line)
if m !== nothing
push!(offsets, parse(Int64, m.captures[1]))
elseif occursin("</index>", line) || occursin("</indexedmzML>", line)
# 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
function load_mzml_lazy(file_path::String)
stream = open(file_path)
# This is the main lazy-loading function for mzML, now returning an MSIData object.
function load_mzml_lazy(file_path::String; cache_size::Int=100)
println("DEBUG: Opening file stream for $file_path")
stream = open(file_path, "r")
try
println("DEBUG: Seeking to end of file.")
seekend(stream)
skip(stream, -4096) # Read last 4KB
println("DEBUG: Skipping back to read footer.")
skip(stream, -4096)
footer = read(stream, String)
println("DEBUG: Footer read successfully.")
index_offset_match = match(r"<indexListOffset>(\d+)</indexListOffset>", footer)
if index_offset_match === nothing
close(stream)
error("Could not find <indexListOffset>. File may not be an indexed mzML.")
end
println("DEBUG: Found indexListOffset.")
index_offset = parse(Int64, index_offset_match.captures[1])
println("DEBUG: Seeking to index list at offset $index_offset.")
seek(stream, index_offset)
println("DEBUG: Searching for '<index name=\"spectrum\">'.")
if find_tag(stream, r"<index\s+name=\"spectrum\"") === nothing
close(stream)
error("Could not find spectrum index.")
end
println("DEBUG: Found spectrum index tag. ")
println("DEBUG: Parsing spectrum offsets...")
spectrum_offsets = parse_offset_list(stream)
if isempty(spectrum_offsets)
close(stream)
error("No spectrum offsets found.")
end
println("DEBUG: Found $(length(spectrum_offsets)) spectrum offsets.")
println("DEBUG: Parsing metadata for each spectrum...")
spectra_metadata = [parse_spectrum_metadata(stream, offset) for offset in spectrum_offsets]
println("DEBUG: Metadata parsing complete.")
# We return the stream inside MzMLData, so we don't close it here.
# The finalizer on MzMLData will handle it.
return MzMLData(stream, spectra_metadata)
# 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
source = MzMLSource(stream, mz_format, intensity_format)
println("DEBUG: Creating MSIData object.")
return MSIData(source, spectra_metadata, (0, 0), nothing, cache_size)
catch e
close(stream) # Close file on error
close(stream) # Ensure stream is closed on error
rethrow(e)
end
end
# ============================================================================
# Data Access Functions
# ============================================================================
function read_binary_vector(io::IO, asset::SpectrumAsset)
seek(io, asset.offset)
# Read the raw base64 data. The regex is slow. Reading bytes is better.
# The old code did `Vector{UInt8}(base64Vec.captures[1])`. This implies ASCII.
raw_b64 = read(io, asset.encoded_length)
decoded_bytes = base64decode(String(raw_b64))
bytes_io = IOBuffer(asset.is_compressed ? Libz.inflate(decoded_bytes) : decoded_bytes)
n_elements = Int(bytes_io.size / sizeof(asset.format))
out_array = Array{asset.format}(undef, n_elements)
read!(bytes_io, out_array)
return out_array
end
function get_spectrum(data::MzMLData, index::Int)
meta = data.spectra_metadata[index]
mz_array = read_binary_vector(data.file_handle, meta.mz_asset)
intensity_array = read_binary_vector(data.file_handle, meta.intensity_asset)
return (id=meta.id, mz=mz_array, intensity=intensity_array)
end
"""
LoadMzml(fileName::String)
Eagerly loads all spectra from a .mzML file.
Provided for backward compatibility.
Provided for backward compatibility. Now uses the new MSIData architecture.
"""
function LoadMzml(fileName::String)
mzml_data = load_mzml_lazy(fileName)
# Use the lazy loader to get the MSIData object
msi_data = load_mzml_lazy(fileName, cache_size=0) # No need to cache if we load all
num_spectra = length(mzml_data.spectra_metadata)
num_spectra = length(msi_data.spectra_metadata)
spectra_matrix = Array{Any}(undef, (2, num_spectra))
for i in 1:num_spectra
spectrum = get_spectrum(mzml_data, i)
spectra_matrix[1, i] = spectrum.mz
spectra_matrix[2, i] = spectrum.intensity
# Use the new GetSpectrum API
mz, intensity = GetSpectrum(msi_data, i)
spectra_matrix[1, i] = mz
spectra_matrix[2, i] = intensity
end
# Eager load, so we can close the handle now.
close(mzml_data.file_handle)
# The finalizer on MSIData will close the handle, so no need to close it here.
return spectra_matrix
end

350
src/mzML_old.jl
View File

@ -1,350 +0,0 @@
# *******************************************************************
# Load scans from mzML file
# fileName: Full name path
# *******************************************************************
"""
`LoadMzml( fileName )`
Load the mzML file as a Matrix. Each column stores a single scan,
x-axis and y-axis scan data are stored in the first and second
row respectively.
# Arguments
* `fileName`: Full path name of the mzML file
# Examples
```julia
spectra = LoadMzml( "T9_A1.mzML" )
size( spectra )
(2, 10)
using Plots
plot( spectra[1,4], spectra[2,4] )
```
The following image will appear in an auxiliary window
![](./Plots.bmp)
"""
function LoadMzml( fileName )
# Open file hFile = open( joinpath( data_dir, "Col_1.mzML" ), "r+" )
stream = open( fileName, "r+" )
intensity = []
# Load footer
seekend( stream )
skip( stream, -160 )
raw = read( stream, 160 )
footer = String( raw )
# Search "indexListOffset" tag
index = match( r"<indexListOffset>(?<Ofst>\d+)</indexListOffset>", footer )
if index !== nothing
# Jump to index table
seek( stream, parse( Int64, index.captures[1] ) )
# Find first "index" tag
index = FindTag( stream, r"<index +name=\"spectrum" )
# Load offset list and load spectra
if index !== nothing
intensity = LoadSpectra( stream, ParseOffsetList( stream ) )
end
end
# Close file handle
close( stream )
# Return spectral matrix
return intensity
end
# *******************************************************************
# Parse indexList offsets
# stream: Source file stream
# *******************************************************************
function ParseOffsetList( stream )
# List of 8 elements
entry = zeros( Int64, 8 )
used = 0
unused = 8
while true
# Get offset string
raw = readline( stream )
index = match( r">(?<Ofst>\d+)<", String( raw ) )
# Exit condition
if index === nothing
break
end
# Save index
used += 1
unused -= 1
entry[used] = parse( Int64, index.captures[1] )
# Grow list when needed
if unused == 0
append!( entry, zeros( Int64, 8 ) )
unused = 8
end
end
return entry[1:used]
end
# *******************************************************************
# Read vector data
# stream: Source file stream
# *******************************************************************
function ReadVector( stream, axis )
# Read base64 vector
base64Vec = FindTag( stream, r"<binary>([^<]+)<" )
base64Vec = Vector{ UInt8 }( base64Vec.captures[1] )
# Fill a IO stream with decoded data
if axis.Packed == 1
io = IOBuffer( Libz.inflate( Decode64( base64Vec ) ) )
else
io = IOBuffer( Decode64( base64Vec ) )
end
# Convert IO stream content as dataType vector
nElem = Int32( io.size / sizeof( axis.Format ) )
out = Array{ axis.Format }( undef, nElem )
read!( io, out )
return out
end
# *******************************************************************
# Get the decoding options for the first vector in spectrum
# stream: Valid file stream
# lastTag: Last tag to skip i.e. "binaryDataArray"
# *******************************************************************
function AxisConfig( stream )
offset = position( stream )
skip = 0
while true
# Load line
raw = readline( stream )
tag = match( r"^\s*<([^\s]+)", raw )
start = length( tag.match ) + 1
# Is a "cvParam" tag?
if tag.captures[1] == "cvParam"
# Get accesion field value
attribute = GetAttribute( raw[ start:end ], "accession" )
# Subtract file name length from skip field
if attribute.captures[1] != "MS:1000796"
continue
else
start += attribute.offset + length( attribute.match )
value = GetAttribute( raw[ start:end ], "value" )
skip += length( value.captures[1] )
continue
end
# Exit condition
elseif tag.captures[1] == "binaryDataArray"
value = GetAttribute( raw[ start:end ], "encodedLength" )
offset = position( stream ) - offset
front = ConfigureSpecDim( stream )
front.Skip += offset - length( value.captures[1] ) - skip
return ( front )
end
end
end
# *******************************************************************
# Load Spectra and return a matrix
# stream: Source file stream
# *******************************************************************
function LoadSpectra( stream, offset )
# Move file pointer towards first spectrum tag
seek( stream, offset[1] )
# Decodes axis information and skip counter
firstAxis = AxisConfig( stream )
secondAxis = AxisConfig( stream )
# Allocates memory for spectra storage
specCount = length( offset )
spectrum = Array{ Any }( undef, ( 2, specCount ) )
index = 1
order = 1 + ( firstAxis.Axis == 2 )
for i in 1:specCount
# Load first vector
seek( stream, offset[ i ] )
skip( stream, firstAxis.Skip )
spectrum[ order, i ] = ReadVector( stream, firstAxis )
# Load second vector
skip( stream, secondAxis.Skip )
spectrum[ xor( order,3 ), i ] = ReadVector( stream, secondAxis )
end
return( spectrum )
end
# *******************************************************************
# Decoding base64 blocks
# *******************************************************************
# Decoding exchange table
base64dec = Vector{UInt8}( [
62, # (01) + -> 43d, 0x2B
0, # (02) ,
0, # (03) -
0, # (04) .
63, # (05) /
52, # (06) 0
53, # (07) 1
54, # (08) 2
55, # (09) 3
56, # (10) 4
57, # (11) 5
58, # (12) 6
59, # (13) 7
60, # (14) 8
61, # (15) 9
0, # (16) :
0, # (17) ;
0, # (18) <
0, # (19) =
0, # (20) >
0, # (21) ?
0, # (22) @
0, # (23) A
1, # (24) B
2, # (25) C
3, # (26) D
4, # (27) E
5, # (28) F
6, # (29) G
7, # (40) H
8, # (41) I
9, # (42) J
10, # (43) K
11, # (44) L
12, # (45) M
13, # (46) N
14, # (47) O
15, # (48) P
16, # (49) Q
17, # (40) R
18, # (41) S
19, # (42) T
20, # (43) U
21, # (44) V
22, # (45) W
23, # (46) X
24, # (47) Y
25, # (48) Z
0, # (49) [
0, # (50) \
0, # (51) ]
0, # (52) ^
0, # (53) _
0, # (54) `
26, # (55) a
27, # (56) b
28, # (57) c
29, # (58) d
30, # (59) e
31, # (60) f
32, # (61) g
33, # (62) h
34, # (63) i
35, # (64) j
36, # (65) k
37, # (66) l
38, # (67) m
39, # (68) n
40, # (69) o
41, # (70) p
42, # (71) q
43, # (72) r
44, # (73) s
45, # (74) t
46, # (75) u
47, # (76) v
48, # (77) w
49, # (78) x
50, # (79) y
51, # (80) z
] )
# base64 decoding
function DecodeTriplet( data, base64 )
# Character to binary equivalent
data = data .- 0x2A
data[1] = base64[ data[1] ]
data[2] = base64[ data[2] ]
data[3] = base64[ data[3] ]
data[4] = base64[ data[4] ]
# Get Triplet
decode = Array{UInt8}( undef,3 )
decode[1] = data[1] << 2 + data[2] >> 4
decode[2] = data[2] << 4 + data[3] >> 2
decode[3] = data[3] << 6 + data[4]
return decode
end
function Decode64( data )
# Reserve uninitialized memory
block = divrem( sizeof( data ), 4 )
decode = Array{UInt8}( undef, 3*block[1] )
# Recover binary from Chars
index = 1
first = 1
last = 3
for i in 1:block[1]
decode[first:last] = DecodeTriplet( data[index:index+3], base64dec )
index += 4
first += 3
last += 3
end
# Adjust block size
if data[end-1] == 0x3D
decode = decode[1:end-2]
elseif data[end] == 0x3D
decode = decode[1:end-1]
end
return decode
end

21
test/readme.md Normal file
View File

@ -0,0 +1,21 @@
to run the tests:
time julia --project=. test/run_tests.jl
to enable the different test cases, there are boolean variables
test1 = true
test2 = true
test3 = true
set them to true to perform that specific test.
test 1 consists in mzml validation and spectrum plotting
test 2 consists in imzml conversion and validation
test 3 consists in imzml validation and spectrum and slice plotting
use the global variables:
TEST_MZML_FILE, SPECTRUM_TO_PLOT, CONVERSION_SOURCE_MZML, CONVERSION_SYNC_FILE, CONVERSION_TARGET_IMZML, TEST_IMZML_FILE, MZ_VALUE_FOR_SLICE, MZ_TOLERANCE, COORDS_TO_PLOT, RESULTS_DIR
to modify input that the test script will recieve.

299
test/run_tests.jl Normal file
View File

@ -0,0 +1,299 @@
# 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 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/CE4_BF_R1/CE4_BF_R1.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"
# The corresponding synchronization text file.
const CONVERSION_SYNC_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.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"
# 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_TOLERANCE = 0.1
const MZ_TOLERANCE = 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 = true
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... ")
mz, intensity = @time GetSpectrum(msi_data, idx)
@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
# 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
# ===================================================================
# 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
# Get the msi data from the mzml
println("Plotting a sample spectrum from $TEST_MZML_FILE...")
msi_data = OpenMSIData(TEST_MZML_FILE)
mz, intensity = GetSpectrum(msi_data, SPECTRUM_TO_PLOT)
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")
# 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)
# 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...")
msi_data = OpenMSIData(TEST_IMZML_FILE)
x_coord, y_coord = COORDS_TO_PLOT
# Get the x y coordinate spectrum data
mz, intensity = GetSpectrum(msi_data, x_coord, y_coord)
# 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")
# 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")
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...")
msi_data = OpenMSIData(TEST_IMZML_FILE)
plot_slice(msi_data, MZ_VALUE_FOR_SLICE, 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()