added docstrings to most functions in the src folder, altered msi data calculations for fast iterators, improved total spectrumtimes, added closer colorbar bin calculations, plus more small changes
This commit is contained in:
Pixelguy14
parent
22c8d9637b
commit
7a0360ef68
@ -2,7 +2,7 @@
|
||||
|
||||
julia_version = "1.11.7"
|
||||
manifest_format = "2.0"
|
||||
project_hash = "0f3959248d174e05d23b9e9267b00c6d2ead256d"
|
||||
project_hash = "d4aecbadf6a54893101b3b89b57e8ae0ab391000"
|
||||
|
||||
[[deps.ATK_jll]]
|
||||
deps = ["Artifacts", "Glib_jll", "JLLWrappers", "Libdl"]
|
||||
|
||||
@ -15,6 +15,7 @@ GenieFramework = "a59fdf5c-6bf0-4f5d-949c-a137c9e2f353"
|
||||
Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
|
||||
Libz = "2ec943e9-cfe8-584d-b93d-64dcb6d567b7"
|
||||
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
|
||||
Mmap = "a63ad114-7e13-5084-954f-fe012c677804"
|
||||
NativeFileDialog = "e1fe445b-aa65-4df4-81c1-2041507f0fd4"
|
||||
NaturalSort = "c020b1a1-e9b0-503a-9c33-f039bfc54a85"
|
||||
PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5"
|
||||
|
||||
29
app.jl
29
app.jl
@ -15,6 +15,7 @@ using LinearAlgebra
|
||||
using NativeFileDialog # Opens the file explorer depending on the OS
|
||||
using StipplePlotly
|
||||
using Base.Filesystem: mv # To rename files in the system
|
||||
using Printf # Required for @sprintf macro in colorbar generation
|
||||
|
||||
# Bring MSIData into App module's scope
|
||||
using .MSI_src: MSIData, OpenMSIData, GetSpectrum, IterateSpectra, ImzMLSource, _iterate_spectra_fast, MzMLSource, find_mass, ViridisPalette, get_mz_slice, quantize_intensity, save_bitmap, median_filter, save_bitmap, downsample_spectrum, TrIQ
|
||||
@ -295,16 +296,16 @@ include("./julia_imzML_visual.jl")
|
||||
|
||||
fTime = time()
|
||||
eTime = round(fTime - sTime, digits=3)
|
||||
msg = "File loaded in $(eTime) seconds. Calculating total spectrum..."
|
||||
# msg = "File loaded in $(eTime) seconds. Calculating total spectrum..."
|
||||
msg = "File loaded in $(eTime) seconds."
|
||||
|
||||
# Enable UI controls
|
||||
btnStartDisable = !(msi_data.source isa ImzMLSource)
|
||||
btnPlotDisable = false
|
||||
btnSpectraDisable = false
|
||||
SpectraEnabled = true
|
||||
|
||||
"""
|
||||
# --- Automatically generate and display the sum spectrum plot ---
|
||||
# This is still part of the same async task
|
||||
progressSpectraPlot = true
|
||||
sTime = time()
|
||||
|
||||
@ -314,7 +315,7 @@ include("./julia_imzML_visual.jl")
|
||||
fTime = time()
|
||||
eTime = round(fTime - sTime, digits=3)
|
||||
msg = "Total spectrum plot loaded in $(eTime) seconds."
|
||||
|
||||
"""
|
||||
catch e
|
||||
msi_data = nothing
|
||||
msg = "Error loading file: $e"
|
||||
@ -324,7 +325,11 @@ include("./julia_imzML_visual.jl")
|
||||
SpectraEnabled = false
|
||||
@error "File loading failed" exception=(e, catch_backtrace())
|
||||
finally
|
||||
# This now correctly runs after everything is finished
|
||||
# This block will always run at the end of the async task
|
||||
GC.gc()
|
||||
if Sys.islinux()
|
||||
ccall(:malloc_trim, Int32, (Int32,), 0)
|
||||
end
|
||||
progress = false
|
||||
progressSpectraPlot = false
|
||||
end
|
||||
@ -347,7 +352,7 @@ include("./julia_imzML_visual.jl")
|
||||
msg = "No .imzML file loaded or selected file is not an .imzML. Please select a valid file."
|
||||
warning_msg = true
|
||||
elseif Nmass > 0 && Tol > 0 && Tol <= 1 && colorLevel > 1 && colorLevel < 257
|
||||
msg = "Creating image for Nmass=$(Nmass) Tol=$(Tol). Please be patient."
|
||||
msg = "Creating image for m/z=$(Nmass) Tol=$(Tol). Please be patient."
|
||||
try
|
||||
# Use the new get_mz_slice with the centralized MSIData object
|
||||
slice = get_mz_slice(msi_data, Nmass, Tol)
|
||||
@ -371,11 +376,8 @@ include("./julia_imzML_visual.jl")
|
||||
current_triq = "TrIQ_$(text_nmass).bmp"
|
||||
msgtriq = "TrIQ image with the Nmass of $(replace(text_nmass, "_" => "."))"
|
||||
|
||||
bound = MSI_src.get_outlier_thres(slice, triqProb)
|
||||
levels = range(bound[1], stop=bound[2], length=8)
|
||||
levels = vcat(levels, 2 * levels[end] - levels[end-1])
|
||||
Colorbar(fig[1, 1], colormap=cgrad(:viridis, colorLevel, categorical=true), limits=(0, bound[2]), ticks=levels, tickformat=log_tick_formatter, label="Intensity", size=25)
|
||||
save("public/colorbar_TrIQ_$(text_nmass).png", fig)
|
||||
colorbar_path = joinpath("public", "colorbar_TrIQ_$(text_nmass).png")
|
||||
generate_colorbar_image(slice, colorLevel, colorbar_path, use_triq=true, triq_prob=triqProb)
|
||||
colorbarT = "/colorbar_TrIQ_$(text_nmass).png?t=$(timestamp)"
|
||||
current_col_triq = "colorbar_TrIQ_$(text_nmass).png"
|
||||
|
||||
@ -400,9 +402,8 @@ include("./julia_imzML_visual.jl")
|
||||
current_msi = "MSI_$(text_nmass).bmp"
|
||||
msgimg = "Image with the Nmass of $(replace(text_nmass, "_" => "."))"
|
||||
|
||||
levels = range(0, maximum(slice), length=8)
|
||||
Colorbar(fig[1, 1], colormap=cgrad(:viridis, colorLevel, categorical=true), limits=(0, maximum(slice)), ticks=levels, tickformat=log_tick_formatter, label="Intensity", size=25)
|
||||
save("public/colorbar_MSI_$(text_nmass).png", fig)
|
||||
colorbar_path = joinpath("public", "colorbar_MSI_$(text_nmass).png")
|
||||
generate_colorbar_image(slice, colorLevel, colorbar_path)
|
||||
colorbar = "/colorbar_MSI_$(text_nmass).png?t=$(timestamp)"
|
||||
current_col_msi = "colorbar_MSI_$(text_nmass).png"
|
||||
|
||||
|
||||
@ -325,7 +325,66 @@ function log_tick_formatter(values::Vector{Float64})
|
||||
formValues[i]=value
|
||||
end
|
||||
return map((v, e) -> e == 0 ? "$(round(v, sigdigits=2))" : "$(round(v, sigdigits=2))x10" * Makie.UnicodeFun.to_superscript(e), formValues, exponents)
|
||||
end
|
||||
|
||||
function generate_colorbar_image(slice_data::AbstractMatrix, color_levels::Int, output_path::String; use_triq::Bool=false, triq_prob::Float64=0.98)
|
||||
# 1. Determine bounds based on whether TrIQ is used
|
||||
min_val, max_val = if use_triq
|
||||
MSI_src.get_outlier_thres(slice_data, triq_prob)
|
||||
else
|
||||
extrema(slice_data)
|
||||
end
|
||||
|
||||
# 2. Replicate the tick calculation logic from plot_slices
|
||||
bins = color_levels
|
||||
levels = range(min_val, stop=max_val, length=bins + 1)
|
||||
level_range = levels[end] - levels[1]
|
||||
|
||||
if level_range == 0
|
||||
levels = range(min_val - 0.1, stop=max_val + 0.1, length=bins + 1)
|
||||
level_range = 0.2
|
||||
end
|
||||
|
||||
exponent = level_range > 0 ? floor(log10(level_range)) / 3 : 0
|
||||
scale = 10^(3 * exponent)
|
||||
scaled_levels = levels ./ scale
|
||||
|
||||
format_num = level_range > 0 ? floor(log10(level_range)) % 3 : 0
|
||||
labels = if format_num == 0
|
||||
[ @sprintf("%3.2f", lvl) for lvl in scaled_levels]
|
||||
elseif format_num == 1
|
||||
[ @sprintf("%3.2f", lvl) for lvl in scaled_levels]
|
||||
else
|
||||
[ @sprintf("%3.2f", lvl) for lvl in scaled_levels]
|
||||
end
|
||||
|
||||
divisors = 2:7
|
||||
remainders = (bins - 1) .% divisors
|
||||
best_divisor = divisors[findlast(x -> x == minimum(remainders), remainders)]
|
||||
tick_indices = round.(Int, range(1, stop=bins + 1, length=best_divisor + 1))
|
||||
|
||||
if !(1 in tick_indices)
|
||||
pushfirst!(tick_indices, 1)
|
||||
end
|
||||
if !((bins + 1) in tick_indices)
|
||||
push!(tick_indices, bins + 1)
|
||||
end
|
||||
unique!(sort!(tick_indices))
|
||||
|
||||
tick_positions = levels[tick_indices]
|
||||
tick_labels = labels[tick_indices]
|
||||
|
||||
# 3. Create and save the colorbar image
|
||||
fig = Figure(size=(150, 250))
|
||||
Colorbar(fig[1, 1],
|
||||
colormap=cgrad(:viridis, bins),
|
||||
limits=(min_val, max_val),
|
||||
label=(scale == 1 ? "Intensity" : "Intensity ×10^$(round(Int, 3 * exponent))"),
|
||||
ticks=(tick_positions, tick_labels),
|
||||
labelsize=20,
|
||||
ticklabelsize=16
|
||||
)
|
||||
save(output_path, fig)
|
||||
end
|
||||
|
||||
# meanSpectrumPlot recieves the local directory of the image as a string,
|
||||
@ -408,7 +467,6 @@ function xySpectrumPlot(data::MSIData, xCoord::Int, yCoord::Int, imgWidth::Int,
|
||||
plotdata = [trace]
|
||||
plotlayout = layout
|
||||
|
||||
# Return the full data for other uses, and the plot data
|
||||
return plotdata, plotlayout, mz, intensity
|
||||
end
|
||||
|
||||
|
||||
299
src/MSIData.jl
299
src/MSIData.jl
@ -6,27 +6,50 @@ including caching and iteration logic, for handling large mzML and imzML dataset
|
||||
efficiently.
|
||||
"""
|
||||
|
||||
using Base64, Libz, Serialization # For reading binary data
|
||||
using Base64, Libz, Serialization, Printf # 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
|
||||
"""
|
||||
ImzMLSource <: MSDataSource
|
||||
|
||||
A data source for `.imzML` files, holding a handle to the binary `.ibd` file
|
||||
and the expected format for m/z and intensity arrays.
|
||||
"""
|
||||
struct ImzMLSource <: MSDataSource
|
||||
ibd_handle::IO
|
||||
mz_format::Type
|
||||
intensity_format::Type
|
||||
end
|
||||
|
||||
# Concrete type for mzML data sources
|
||||
"""
|
||||
MzMLSource <: MSDataSource
|
||||
|
||||
A data source for `.mzML` files, holding a handle to the `.mzML` file itself
|
||||
(which contains the binary data encoded in Base64) and the expected data formats.
|
||||
"""
|
||||
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)
|
||||
"""
|
||||
SpectrumAsset
|
||||
|
||||
Contains metadata for a single binary data array (m/z or intensity) within a spectrum.
|
||||
|
||||
# Fields
|
||||
- `format`: The data type of the elements (e.g., `Float32`, `Int64`).
|
||||
- `is_compressed`: A boolean flag indicating if the data is compressed (e.g., with zlib).
|
||||
- `offset`: The byte offset of the data within the file (`.ibd` for imzML, `.mzML` for mzML).
|
||||
- `encoded_length`: The length of the data. For mzML, this is the Base64 encoded length.
|
||||
For imzML, this is the number of elements in the array.
|
||||
- `axis_type`: A symbol (`:mz` or `:intensity`) indicating the type of data.
|
||||
"""
|
||||
struct SpectrumAsset
|
||||
format::Type
|
||||
is_compressed::Bool
|
||||
@ -37,7 +60,17 @@ struct SpectrumAsset
|
||||
axis_type::Symbol
|
||||
end
|
||||
|
||||
# Metadata for a single spectrum, common to both formats
|
||||
"""
|
||||
SpectrumMetadata
|
||||
|
||||
Contains all metadata for a single spectrum, common to both imzML and mzML formats.
|
||||
|
||||
# Fields
|
||||
- `x`, `y`: The spatial coordinates of the spectrum (for imzML only).
|
||||
- `id`: The unique identifier string for the spectrum (for mzML only).
|
||||
- `mz_asset`: A `SpectrumAsset` for the m/z array.
|
||||
- `int_asset`: A `SpectrumAsset` for the intensity array.
|
||||
"""
|
||||
struct SpectrumMetadata
|
||||
# For imzML
|
||||
x::Int32
|
||||
@ -51,7 +84,22 @@ struct SpectrumMetadata
|
||||
int_asset::SpectrumAsset
|
||||
end
|
||||
|
||||
# The main unified data object
|
||||
"""
|
||||
MSIData
|
||||
|
||||
The primary object for interacting with mass spectrometry data. It provides a unified
|
||||
interface for both `.mzML` and `.imzML` files and includes an LRU cache for
|
||||
efficient repeated access to spectra.
|
||||
|
||||
# Fields
|
||||
- `source`: The underlying `MSDataSource` (`ImzMLSource` or `MzMLSource`).
|
||||
- `spectra_metadata`: A vector of `SpectrumMetadata` for all spectra in the file.
|
||||
- `image_dims`: A tuple `(width, height)` of the spatial dimensions (for imzML).
|
||||
- `coordinate_map`: A matrix mapping `(x, y)` coordinates to a linear spectrum index (for imzML).
|
||||
- `cache`: A dictionary holding cached spectra.
|
||||
- `cache_order`: A vector tracking the usage order for the LRU cache.
|
||||
- `cache_size`: The maximum number of spectra to store in the cache.
|
||||
"""
|
||||
mutable struct MSIData
|
||||
source::MSDataSource
|
||||
spectra_metadata::Vector{SpectrumMetadata}
|
||||
@ -81,6 +129,26 @@ end
|
||||
|
||||
# --- Internal function for reading binary data --- #
|
||||
|
||||
"""
|
||||
read_binary_vector(io::IO, asset::SpectrumAsset)
|
||||
|
||||
Reads and decodes a single binary data vector (like m/z or intensity array)
|
||||
from a `.mzML` file. The data is expected to be Base64-encoded and may be
|
||||
compressed.
|
||||
|
||||
This internal function handles:
|
||||
1. Reading the raw Base64 string.
|
||||
2. Decoding from Base64.
|
||||
3. Decompressing the data if `asset.is_compressed` is true.
|
||||
4. Converting the byte order from network (big-endian) to host order.
|
||||
|
||||
# Arguments
|
||||
- `io`: The IO stream of the `.mzML` file.
|
||||
- `asset`: The `SpectrumAsset` containing metadata for the array.
|
||||
|
||||
# Returns
|
||||
- A `Vector` of the appropriate type containing the decoded data.
|
||||
"""
|
||||
function read_binary_vector(io::IO, asset::SpectrumAsset)
|
||||
seek(io, asset.offset)
|
||||
raw_b64 = read(io, asset.encoded_length)
|
||||
@ -94,6 +162,20 @@ function read_binary_vector(io::IO, asset::SpectrumAsset)
|
||||
end
|
||||
|
||||
# Overload for different source types
|
||||
"""
|
||||
read_spectrum_from_disk(source::ImzMLSource, meta::SpectrumMetadata)
|
||||
|
||||
Reads a single spectrum's m/z and intensity arrays directly from the `.ibd`
|
||||
binary file for an `.imzML` dataset. It handles reading the raw binary data
|
||||
and converting it from little-endian to the host's native byte order.
|
||||
|
||||
# Arguments
|
||||
- `source`: An `ImzMLSource` containing the file handle and data formats.
|
||||
- `meta`: The `SpectrumMetadata` for the spectrum to be read.
|
||||
|
||||
# Returns
|
||||
- A tuple `(mz, intensity)` containing the two requested data arrays.
|
||||
"""
|
||||
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.
|
||||
@ -106,9 +188,27 @@ function read_spectrum_from_disk(source::ImzMLSource, meta::SpectrumMetadata)
|
||||
seek(source.ibd_handle, meta.int_asset.offset)
|
||||
read!(source.ibd_handle, intensity)
|
||||
|
||||
# imzML data is little-endian. Convert to host byte order.
|
||||
mz .= ltoh.(mz)
|
||||
intensity .= ltoh.(intensity)
|
||||
|
||||
return mz, intensity
|
||||
end
|
||||
|
||||
"""
|
||||
read_spectrum_from_disk(source::MzMLSource, meta::SpectrumMetadata)
|
||||
|
||||
Reads a single spectrum's m/z and intensity arrays from a `.mzML` file.
|
||||
This is a high-level function that orchestrates the reading process by calling
|
||||
`read_binary_vector` for each of the m/z and intensity arrays.
|
||||
|
||||
# Arguments
|
||||
- `source`: A `MzMLSource` containing the file handle.
|
||||
- `meta`: The `SpectrumMetadata` for the spectrum to be read.
|
||||
|
||||
# Returns
|
||||
- A tuple `(mz, intensity)` containing the two requested data arrays.
|
||||
"""
|
||||
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)
|
||||
@ -182,13 +282,14 @@ end
|
||||
using Serialization
|
||||
|
||||
function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000)
|
||||
println("Calculating total spectrum (2-pass method for imzML)...")
|
||||
println("Calculating total spectrum for imzML (2-pass method)...")
|
||||
total_start_time = time_ns()
|
||||
|
||||
# 1. First Pass: Find the global m/z range
|
||||
pass1_start_time = time_ns()
|
||||
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)
|
||||
@ -196,36 +297,66 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000)
|
||||
global_max_mz = max(global_max_mz, local_max)
|
||||
end
|
||||
end
|
||||
pass1_duration = (time_ns() - pass1_start_time) / 1e9
|
||||
|
||||
if !isfinite(global_min_mz)
|
||||
@warn "Could not determine a valid m/z range. All spectra might be empty."
|
||||
@warn "Could not determine a valid m/z range for imzML. All spectra might be empty."
|
||||
return (Float64[], Float64[])
|
||||
end
|
||||
println(" Global m/z range found: [$(global_min_mz), $(global_max_mz)]")
|
||||
|
||||
# 2. Define Bins
|
||||
# 2. Define Bins and precompute constants
|
||||
mz_bins = range(global_min_mz, stop=global_max_mz, length=num_bins)
|
||||
intensity_sum = zeros(Float64, num_bins)
|
||||
bin_step = step(mz_bins)
|
||||
inv_bin_step = 1.0 / bin_step # Precompute reciprocal to avoid division
|
||||
min_mz = global_min_mz
|
||||
|
||||
# 3. Second Pass: Sum intensities into bins
|
||||
# 3. Second Pass: Optimized binning
|
||||
pass2_start_time = time_ns()
|
||||
println(" Pass 2: Summing intensities into $num_bins bins...")
|
||||
_iterate_spectra_fast(msi_data) do idx, mz, intensity
|
||||
if isempty(mz)
|
||||
return # continue equivalent
|
||||
return
|
||||
end
|
||||
for i in eachindex(mz)
|
||||
bin_index = clamp(round(Int, (mz[i] - global_min_mz) / bin_step) + 1, 1, num_bins)
|
||||
intensity_sum[bin_index] += intensity[i]
|
||||
|
||||
# Use views to avoid bounds checks on the input arrays
|
||||
mz_view = mz
|
||||
intensity_view = intensity
|
||||
|
||||
# Manual binning without clamp/round for SIMD
|
||||
@inbounds for i in eachindex(mz_view)
|
||||
# Calculate raw bin index (much faster than round+clamp)
|
||||
raw_index = (mz_view[i] - min_mz) * inv_bin_step + 1.0
|
||||
bin_index = trunc(Int, raw_index)
|
||||
|
||||
# Manual bounds checking (faster than clamp)
|
||||
if 1 <= bin_index <= num_bins
|
||||
intensity_sum[bin_index] += intensity_view[i]
|
||||
elseif bin_index < 1
|
||||
intensity_sum[1] += intensity_view[i]
|
||||
else # bin_index > num_bins
|
||||
intensity_sum[num_bins] += intensity_view[i]
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
println("Total spectrum calculation complete.")
|
||||
pass2_duration = (time_ns() - pass2_start_time) / 1e9
|
||||
|
||||
total_duration = (time_ns() - total_start_time) / 1e9
|
||||
println("\n--- imzML Profiling ---")
|
||||
@printf " Pass 1 (I/O only): %.2f seconds\n" pass1_duration
|
||||
@printf " Pass 2 (I/O + Binning): %.2f seconds\n" pass2_duration
|
||||
@printf " Est. Binning Overhead: %.2f seconds\n" (pass2_duration - pass1_duration)
|
||||
@printf " Total Function Time: %.2f seconds\n" total_duration
|
||||
println("-------------------------\n")
|
||||
|
||||
println("Total spectrum calculation complete for imzML.")
|
||||
return (collect(mz_bins), intensity_sum)
|
||||
end
|
||||
|
||||
function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000)
|
||||
println("Calculating total spectrum (single-pass optimization for mzML)...")
|
||||
println("Calculating total spectrum for mzML (optimized single-pass method)...")
|
||||
total_start_time = time_ns()
|
||||
num_spectra = length(msi_data.spectra_metadata)
|
||||
if num_spectra == 0
|
||||
return (Float64[], Float64[])
|
||||
@ -233,31 +364,33 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000)
|
||||
|
||||
temp_path, temp_io = mktemp()
|
||||
try
|
||||
# --- Pass 1: Write spectra to temp file and find min/max m/z ---
|
||||
println(" Pass 1: Caching spectra and finding global m/z range...")
|
||||
# --- Pass 1: Read from source, find m/z range, and write decoded spectra to temp file ---
|
||||
pass1_start_time = time_ns()
|
||||
println(" Pass 1: Caching decoded spectra and finding global m/z range...")
|
||||
global_min_mz = Inf
|
||||
global_max_mz = -Inf
|
||||
|
||||
|
||||
_iterate_spectra_fast(msi_data) do idx, mz, intensity
|
||||
Serialization.serialize(temp_io, (mz, intensity))
|
||||
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
|
||||
Serialization.serialize(temp_io, (mz, intensity))
|
||||
end
|
||||
|
||||
flush(temp_io)
|
||||
pass1_duration = (time_ns() - pass1_start_time) / 1e9
|
||||
|
||||
if !isfinite(global_min_mz)
|
||||
@warn "Could not determine a valid m/z range. All spectra might be empty."
|
||||
@warn "Could not determine a valid m/z range for mzML. All spectra might be empty."
|
||||
return (Float64[], Float64[])
|
||||
end
|
||||
println(" Global m/z range found: [$(global_min_mz), $(global_max_mz)]")
|
||||
|
||||
# --- Pass 2: Read from temp file and bin intensities ---
|
||||
# --- Pass 2: Read from fast temp file and bin intensities ---
|
||||
pass2_start_time = time_ns()
|
||||
println(" Pass 2: Reading from cache and summing intensities into $num_bins bins...")
|
||||
|
||||
seekstart(temp_io)
|
||||
mz_bins = range(global_min_mz, stop=global_max_mz, length=num_bins)
|
||||
intensity_sum = zeros(Float64, num_bins)
|
||||
@ -265,7 +398,6 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000)
|
||||
|
||||
while !eof(temp_io)
|
||||
mz, intensity = Serialization.deserialize(temp_io)::Tuple{AbstractVector, AbstractVector}
|
||||
|
||||
if isempty(mz)
|
||||
continue
|
||||
end
|
||||
@ -274,18 +406,27 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000)
|
||||
intensity_sum[bin_index] += intensity[i]
|
||||
end
|
||||
end
|
||||
|
||||
println("Total spectrum calculation complete.")
|
||||
pass2_duration = (time_ns() - pass2_start_time) / 1e9
|
||||
|
||||
total_duration = (time_ns() - total_start_time) / 1e9
|
||||
println("\n--- mzML Profiling ---")
|
||||
@printf " Pass 1 (Read+Decode+Cache): %.2f seconds\n" pass1_duration
|
||||
@printf " Pass 2 (Read Cache+Bin): %.2f seconds\n" pass2_duration
|
||||
@printf " Total Function Time: %.2f seconds\n" total_duration
|
||||
println("----------------------\n")
|
||||
|
||||
println("Total spectrum calculation complete for mzML.")
|
||||
return (collect(mz_bins), intensity_sum)
|
||||
|
||||
finally
|
||||
close(temp_io)
|
||||
rm(temp_path, force=true)
|
||||
println(" Temporary cache file removed.")
|
||||
end
|
||||
end
|
||||
|
||||
"""
|
||||
get_total_spectrum(msi_data::MSIData; num_bins::Int=20000) -> Tuple{Vector{Float64}, Vector{Float64}}
|
||||
get_total_spectrum(msi_data::MSIData; num_bins::Int=2000) -> Tuple{Vector{Float64}, Vector{Float64}}
|
||||
|
||||
Calculates the sum of all spectra in the dataset by binning.
|
||||
This function dispatches to a specialized implementation based on the file type
|
||||
@ -301,7 +442,6 @@ function get_total_spectrum(msi_data::MSIData; num_bins::Int=2000)
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
get_average_spectrum(msi_data::MSIData; num_bins::Int=20000) -> Tuple{Vector{Float64}, Vector{Float64}}
|
||||
|
||||
@ -328,6 +468,12 @@ end
|
||||
|
||||
# --- Iterator Implementation --- #
|
||||
|
||||
"""
|
||||
MSIDataIterator
|
||||
|
||||
A stateful iterator for sequentially accessing spectra from an `MSIData` object.
|
||||
This is created by the `IterateSpectra` function.
|
||||
"""
|
||||
struct MSIDataIterator
|
||||
data::MSIData
|
||||
end
|
||||
@ -336,7 +482,7 @@ end
|
||||
IterateSpectra(data::MSIData)
|
||||
|
||||
Returns an iterator that yields each spectrum, processing the file sequentially
|
||||
with minimal memory overhead.
|
||||
with minimal memory overhead. This iterator supports caching via `GetSpectrum`.
|
||||
|
||||
This function is the core of the "Event-driven" access pattern.
|
||||
"""
|
||||
@ -348,6 +494,20 @@ end
|
||||
Base.length(it::MSIDataIterator) = length(it.data.spectra_metadata)
|
||||
Base.eltype(it::MSIDataIterator) = Tuple{Int, Tuple{Vector, Vector}} # Index, (mz, intensity)
|
||||
|
||||
"""
|
||||
Base.iterate(it::MSIDataIterator, state=1)
|
||||
|
||||
Advances the iterator, returning the next spectrum in the sequence. It calls
|
||||
`GetSpectrum`, which means the iteration process is cached according to the
|
||||
`MSIData` object's cache settings.
|
||||
|
||||
# Arguments
|
||||
- `it`: The `MSIDataIterator` instance.
|
||||
- `state`: The index of the next spectrum to retrieve.
|
||||
|
||||
# Returns
|
||||
- A tuple `((index, spectrum), next_state)` or `nothing` if iteration is complete.
|
||||
"""
|
||||
function Base.iterate(it::MSIDataIterator, state=1)
|
||||
if state > length(it.data.spectra_metadata)
|
||||
return nothing # End of iteration
|
||||
@ -364,35 +524,14 @@ end
|
||||
# --- High-performance Internal Iterator --- #
|
||||
|
||||
"""
|
||||
_iterate_spectra_fast(f::Function, data::MSIData)
|
||||
_iterate_spectra_fast_impl(f::Function, data::MSIData, source::ImzMLSource)
|
||||
|
||||
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.
|
||||
Internal implementation of the fast iterator for `.imzML` files. It reads
|
||||
data directly from the `.ibd` file stream and reuses pre-allocated buffers
|
||||
to minimize memory allocations and overhead, making it ideal for bulk processing tasks.
|
||||
"""
|
||||
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)
|
||||
|
||||
# Optimized implementation: allocates arrays per spectrum and uses an efficient I/O pattern.
|
||||
for i in 1:length(data.spectra_metadata)
|
||||
meta = data.spectra_metadata[i]
|
||||
nPoints = meta.mz_asset.encoded_length
|
||||
@ -402,19 +541,43 @@ function _iterate_spectra_fast_impl(f::Function, data::MSIData, source::ImzMLSou
|
||||
continue
|
||||
end
|
||||
|
||||
mz_view = view(mz_buffer, 1:nPoints)
|
||||
int_view = view(int_buffer, 1:nPoints)
|
||||
# Allocate fresh arrays for each spectrum.
|
||||
mz = Vector{source.mz_format}(undef, nPoints)
|
||||
intensity = Vector{source.intensity_format}(undef, nPoints)
|
||||
|
||||
seek(source.ibd_handle, meta.mz_asset.offset)
|
||||
read!(source.ibd_handle, mz_view)
|
||||
# --- I/O OPTIMIZATION: Seek only once per spectrum ---
|
||||
# The mz and intensity data are contiguous, so after reading the first,
|
||||
# we can immediately read the second without a costly second seek.
|
||||
if meta.mz_asset.offset < meta.int_asset.offset
|
||||
# m/z is first, seek to it.
|
||||
seek(source.ibd_handle, meta.mz_asset.offset)
|
||||
# Read m/z, then immediately read intensity.
|
||||
read!(source.ibd_handle, mz)
|
||||
read!(source.ibd_handle, intensity)
|
||||
else
|
||||
# intensity is first, seek to it.
|
||||
seek(source.ibd_handle, meta.int_asset.offset)
|
||||
# Read intensity, then immediately read m/z.
|
||||
read!(source.ibd_handle, intensity)
|
||||
read!(source.ibd_handle, mz)
|
||||
end
|
||||
|
||||
seek(source.ibd_handle, meta.int_asset.offset)
|
||||
read!(source.ibd_handle, int_view)
|
||||
# Convert byte order.
|
||||
mz .= ltoh.(mz)
|
||||
intensity .= ltoh.(intensity)
|
||||
|
||||
f(i, mz_view, int_view)
|
||||
f(i, mz, intensity)
|
||||
end
|
||||
end
|
||||
|
||||
"""
|
||||
_iterate_spectra_fast_impl(f::Function, data::MSIData, source::MzMLSource)
|
||||
|
||||
Internal implementation of the fast iterator for `.mzML` files. It iterates
|
||||
through each spectrum, decodes the Base64 data on the fly, and calls the
|
||||
provided function. It bypasses the `GetSpectrum` cache to avoid storing
|
||||
all decoded spectra in memory.
|
||||
"""
|
||||
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
|
||||
@ -430,6 +593,16 @@ function _iterate_spectra_fast_impl(f::Function, data::MSIData, source::MzMLSour
|
||||
end
|
||||
end
|
||||
|
||||
"""
|
||||
_iterate_spectra_fast(f::Function, data::MSIData)
|
||||
|
||||
An internal, high-performance iterator for bulk processing that bypasses the cache.
|
||||
It dispatches to a specialized implementation based on the data source type
|
||||
(`ImzMLSource` or `MzMLSource`).
|
||||
|
||||
- `f`: A function to call for each spectrum, with signature `f(index, mz, intensity)`.
|
||||
- `data`: The `MSIData` object.
|
||||
"""
|
||||
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)
|
||||
|
||||
@ -20,7 +20,7 @@ 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)
|
||||
function OpenMSIData(filepath::String; cache_size=300)
|
||||
if endswith(lowercase(filepath), ".mzml")
|
||||
return load_mzml_lazy(filepath, cache_size=cache_size)
|
||||
elseif endswith(lowercase(filepath), ".imzml")
|
||||
|
||||
@ -11,13 +11,33 @@ using DataFrames, Printf, CSV
|
||||
|
||||
|
||||
# A struct to hold the processed pixel data before exporting
|
||||
"""
|
||||
ProcessedPixel
|
||||
|
||||
A temporary struct to hold the data for a single, fully rendered pixel before
|
||||
it is written to the final `.ibd` file.
|
||||
|
||||
# Fields
|
||||
- `coords`: A tuple `(x, y)` of the pixel's spatial coordinates.
|
||||
- `mz`: The m/z array for the pixel.
|
||||
- `intensity`: The calculated intensity array for the pixel.
|
||||
"""
|
||||
struct ProcessedPixel
|
||||
coords::Tuple{Int, Int}
|
||||
mz::Vector{Float64} # Assuming m/z is consistent, can be optimized later
|
||||
intensity::Vector{Float32}
|
||||
end
|
||||
|
||||
# This struct will hold the metadata needed for the XML file
|
||||
"""
|
||||
BinaryMetadata
|
||||
|
||||
A struct to hold the file offset and length for a spectrum's binary data arrays
|
||||
(m/z and intensity) after they have been written to the `.ibd` file.
|
||||
|
||||
# Fields
|
||||
- `mz_offset`, `mz_length`: Byte offset and length for the m/z array.
|
||||
- `int_offset`, `int_length`: Byte offset and length for the intensity array.
|
||||
"""
|
||||
struct BinaryMetadata
|
||||
mz_offset::UInt64
|
||||
mz_length::UInt64
|
||||
@ -28,8 +48,15 @@ end
|
||||
"""
|
||||
GetMzmlScanTime_linebyline(fileName::String)
|
||||
|
||||
(Internal) Slow, line-by-line parser for scan times. Used as a fallback if
|
||||
the .mzML file index is missing or corrupt.
|
||||
Parses a `.mzML` file line-by-line to extract the scan start time for each
|
||||
spectrum. This is a slow and memory-intensive fallback method used only when
|
||||
the faster, index-based `GetMzmlScanTime` fails.
|
||||
|
||||
# Arguments
|
||||
- `fileName`: Path to the `.mzML` file.
|
||||
|
||||
# Returns
|
||||
- A `Matrix{Int64}` where each row is `[spectrum_index, time_in_milliseconds]`.
|
||||
"""
|
||||
function GetMzmlScanTime_linebyline(fileName::String)
|
||||
times = Tuple{Int64, Int64}[]
|
||||
@ -350,6 +377,28 @@ function MatchAcquireTime(sync_file_path::String, scans::Matrix{Int64}; img_widt
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
RenderPixel(pixel_info, scans, msi_data, scan_time_deltas, pixel_time_deltas)
|
||||
|
||||
Reconstructs the spectrum for a single pixel by combining intensities from the
|
||||
raw MS scans that occurred during the pixel's acquisition time. It uses a
|
||||
weighted interpolation scheme based on the relative timing of scans and pixels.
|
||||
|
||||
This function handles three cases:
|
||||
1. A single scan falls entirely within the pixel's time window.
|
||||
2. The pixel's time window is covered by two partial scans.
|
||||
3. The pixel's time window covers one or more full scans plus two partial scans.
|
||||
|
||||
# Arguments
|
||||
- `pixel_info`: A row from the timing matrix containing the pixel's time and scan indices.
|
||||
- `scans`: The matrix of scan times.
|
||||
- `msi_data`: The `MSIData` object for the source `.mzML` file.
|
||||
- `scan_time_deltas`: Pre-calculated time durations for each scan.
|
||||
- `pixel_time_deltas`: Pre-calculated time durations for each pixel.
|
||||
|
||||
# Returns
|
||||
- A tuple `(mz_array, intensity_array)` for the rendered pixel spectrum.
|
||||
"""
|
||||
function RenderPixel(pixel_info, scans, msi_data::MSIData, scan_time_deltas, pixel_time_deltas)
|
||||
pixel_time = pixel_info[3]
|
||||
first_scan = pixel_info[4]
|
||||
@ -408,6 +457,26 @@ function RenderPixel(pixel_info, scans, msi_data::MSIData, scan_time_deltas, pix
|
||||
return (mz_array, new_intensity)
|
||||
end
|
||||
|
||||
"""
|
||||
ConvertMzmlToImzml(source_file, target_ibd_file, timing_matrix, scans)
|
||||
|
||||
Orchestrates the conversion of spectra from a `.mzML` file into a binary `.ibd`
|
||||
file. It iterates through each pixel defined in the `timing_matrix`, calls
|
||||
`RenderPixel` to reconstruct the pixel's spectrum, and writes the resulting
|
||||
m/z and intensity arrays to the `.ibd` file in little-endian byte order.
|
||||
|
||||
# Arguments
|
||||
- `source_file`: Path to the source `.mzML` file.
|
||||
- `target_ibd_file`: Path for the output `.ibd` binary file.
|
||||
- `timing_matrix`: The output from `MatchAcquireTime`, mapping pixels to scans.
|
||||
- `scans`: The matrix of scan times from `GetMzmlScanTime`.
|
||||
|
||||
# Returns
|
||||
- A tuple `(binary_meta_vec, coords_vec, (width, height))` containing:
|
||||
- A vector of `BinaryMetadata` for each spectrum.
|
||||
- A vector of `(x, y)` coordinate tuples.
|
||||
- A tuple of the final image dimensions.
|
||||
"""
|
||||
function ConvertMzmlToImzml(source_file::String, target_ibd_file::String, timing_matrix::Matrix{Int64}, scans::Matrix{Int64})
|
||||
if size(timing_matrix, 1) == 0
|
||||
# Create an empty .ibd file if there's nothing to process
|
||||
@ -489,6 +558,23 @@ function ConvertMzmlToImzml(source_file::String, target_ibd_file::String, timing
|
||||
return binary_meta_vec, coords_vec, (width, height)
|
||||
end
|
||||
|
||||
"""
|
||||
ExportImzml(target_file, binary_meta, coords, dims)
|
||||
|
||||
Generates the `.imzML` metadata file. This XML file contains all the necessary
|
||||
metadata to interpret the corresponding `.ibd` binary file, including references
|
||||
to external data offsets, image dimensions, and CV parameters describing the
|
||||
experiment and data format.
|
||||
|
||||
# Arguments
|
||||
- `target_file`: The path for the output `.imzML` file.
|
||||
- `binary_meta`: A vector of `BinaryMetadata` structs with offset and length info.
|
||||
- `coords`: A vector of `(x, y)` coordinates for each spectrum.
|
||||
- `dims`: A tuple `(width, height)` of the final image dimensions.
|
||||
|
||||
# Returns
|
||||
- `true` on success, `false` on failure.
|
||||
"""
|
||||
function ExportImzml(target_file::String, binary_meta::Vector{BinaryMetadata}, coords::Vector{Tuple{Int, Int}}, dims::Tuple{Int, Int})
|
||||
ibd_file = replace(target_file, r"\.imzML$"i => ".ibd")
|
||||
|
||||
|
||||
@ -59,7 +59,15 @@ end
|
||||
"""
|
||||
configure_spec_dim(stream)
|
||||
|
||||
Fills a `SpecDim` struct by parsing `cvParam` tags from the stream.
|
||||
Parses a block of `<cvParam>` tags from an XML stream to configure a `SpecDim`
|
||||
struct. It reads accessions to determine the data format (e.g., `Float32`),
|
||||
compression status (`zlib`), and axis type (m/z vs. intensity).
|
||||
|
||||
# Arguments
|
||||
- `stream`: An IO stream positioned at the start of the `cvParam` block.
|
||||
|
||||
# Returns
|
||||
- A `SpecDim` struct populated with the parsed configuration.
|
||||
"""
|
||||
function configure_spec_dim(stream)
|
||||
axis = SpecDim(Float64, false, 1, 0)
|
||||
@ -102,6 +110,18 @@ end
|
||||
#
|
||||
# ============================================================================
|
||||
|
||||
"""
|
||||
CVParams
|
||||
|
||||
A mutable struct used during the parsing of mzML files to temporarily hold
|
||||
CV parameter information for a data array before it is converted into a
|
||||
`SpectrumAsset`.
|
||||
|
||||
# Fields
|
||||
- `format`: The data type (e.g., `Float64`).
|
||||
- `is_compressed`: A boolean indicating if the data is compressed.
|
||||
- `axis_type`: A symbol indicating the array type (`:mz` or `:intensity`).
|
||||
"""
|
||||
mutable struct CVParams
|
||||
format::Type
|
||||
is_compressed::Bool
|
||||
|
||||
73
src/imzML.jl
73
src/imzML.jl
@ -1,3 +1,4 @@
|
||||
# src/imzML.jl
|
||||
using Images, Statistics, CairoMakie, DataFrames, Printf, ColorSchemes, StatsBase
|
||||
|
||||
# --- Extracted from imzML.jl ---
|
||||
@ -171,7 +172,10 @@ function load_imzml_lazy(file_path::String; cache_size=100)
|
||||
int_config_idx = findfirst(a -> a.Axis == 2, axis)
|
||||
mz_format = axis[mz_config_idx].Format
|
||||
intensity_format = axis[int_config_idx].Format
|
||||
mz_is_compressed = axis[mz_config_idx].Packed
|
||||
int_is_compressed = axis[int_config_idx].Packed
|
||||
println("DEBUG: m/z format: $mz_format, Intensity format: $intensity_format")
|
||||
println("DEBUG: m/z compressed: $mz_is_compressed, Intensity compressed: $int_is_compressed")
|
||||
|
||||
# --- NEW PARSING LOGIC based on the old, working code ---
|
||||
println("DEBUG: Learning file structure from first spectrum...")
|
||||
@ -212,8 +216,8 @@ function load_imzml_lazy(file_path::String; cache_size=100)
|
||||
end
|
||||
|
||||
# Create modern SpectrumAsset objects
|
||||
mz_asset = SpectrumAsset(mz_format, false, mz_offset, nPoints, :mz)
|
||||
int_asset = SpectrumAsset(intensity_format, false, int_offset, nPoints, :intensity)
|
||||
mz_asset = SpectrumAsset(mz_format, mz_is_compressed, mz_offset, nPoints, :mz)
|
||||
int_asset = SpectrumAsset(intensity_format, int_is_compressed, int_offset, nPoints, :intensity)
|
||||
|
||||
spectra_metadata[k] = SpectrumMetadata(x, y, "", mz_asset, int_asset)
|
||||
|
||||
@ -449,6 +453,21 @@ end
|
||||
#
|
||||
# ============================================================================
|
||||
|
||||
"""
|
||||
get_outlier_thres(img, prob=0.98)
|
||||
|
||||
Calculates dynamic intensity range bounds for an image based on a cumulative
|
||||
probability histogram. This function replicates an R algorithm to find an
|
||||
intensity threshold that corresponds to a given cumulative probability, which
|
||||
is used to exclude outliers before normalization.
|
||||
|
||||
# Arguments
|
||||
- `img`: The input image matrix.
|
||||
- `prob`: The cumulative probability threshold (default: 0.98) for outlier detection.
|
||||
|
||||
# Returns
|
||||
- A tuple `(low, high)` representing the calculated lower and upper intensity bounds.
|
||||
"""
|
||||
function get_outlier_thres(img, prob=0.98)
|
||||
# DO NOT filter zeros. Use all pixel values like R does.
|
||||
int_values = vec(img)
|
||||
@ -493,6 +512,21 @@ function get_outlier_thres(img, prob=0.98)
|
||||
return (low, actual_threshold)
|
||||
end
|
||||
|
||||
"""
|
||||
set_pixel_depth(img, bounds, depth)
|
||||
|
||||
Quantizes the intensity values of an image into a specified number of bins
|
||||
(`depth`) within a given intensity range (`bounds`). Pixels outside the bounds
|
||||
are clipped.
|
||||
|
||||
# Arguments
|
||||
- `img`: The input image matrix.
|
||||
- `bounds`: A tuple `(min, max)` specifying the intensity range for quantization.
|
||||
- `depth`: The number of quantization levels (bins).
|
||||
|
||||
# Returns
|
||||
- A `Matrix{UInt8}` with pixel values quantized to the specified depth.
|
||||
"""
|
||||
function set_pixel_depth(img, bounds, depth)
|
||||
min_val, max_val = bounds
|
||||
bins = depth - 1
|
||||
@ -595,6 +629,19 @@ function quantize_intensity(slice::AbstractMatrix{<:Real}, levels::Integer=256)
|
||||
return image
|
||||
end
|
||||
|
||||
"""
|
||||
median_filter(img)
|
||||
|
||||
Applies a 3x3 median filter to the input image. This is a simple noise
|
||||
reduction technique that replaces each pixel's value with the median value of
|
||||
its 3x3 neighborhood.
|
||||
|
||||
# Arguments
|
||||
- `img`: The input image matrix.
|
||||
|
||||
# Returns
|
||||
- A new matrix containing the filtered image.
|
||||
"""
|
||||
function median_filter(img)
|
||||
# 3x3 median filter implementation
|
||||
return mapwindow(median, img, (3, 3))
|
||||
@ -685,6 +732,28 @@ function display_statistics(slices, names, masses)
|
||||
return all_dfs
|
||||
end
|
||||
|
||||
"""
|
||||
plot_slices(slices, names, masses, output_dir; stage_name, bins=256, dpi=150, global_bounds=nothing)
|
||||
|
||||
Generates and saves a grid of image slice plots. Each row corresponds to a
|
||||
file and each column to a mass, creating a comprehensive overview.
|
||||
|
||||
# Arguments
|
||||
- `slices`: A 2D array of image slice matrices (`n_files` x `n_masses`).
|
||||
- `names`: A vector of file names, used for titling rows.
|
||||
- `masses`: A vector of m/z values, used for titling columns.
|
||||
- `output_dir`: The directory where the output plots will be saved.
|
||||
|
||||
# Keyword Arguments
|
||||
- `stage_name`: A string used to name the output files (e.g., "raw", "normalized").
|
||||
- `bins`: The number of color levels in the heatmap palette.
|
||||
- `dpi`: The resolution for the saved image files.
|
||||
- `global_bounds`: A vector of `(min, max)` tuples, one for each mass, to ensure a consistent
|
||||
color scale across all files for a given mass. If `nothing`, bounds are calculated automatically.
|
||||
|
||||
# Returns
|
||||
- The generated `Figure` object from Makie.
|
||||
"""
|
||||
function plot_slices(slices, names, masses, output_dir; stage_name, bins=256, dpi=150, global_bounds=nothing)
|
||||
n_files, n_masses = size(slices)
|
||||
|
||||
|
||||
73
src/mzML.jl
73
src/mzML.jl
@ -1,9 +1,22 @@
|
||||
# /home/pixel/Documents/Cinvestav_2025/JuliaMSI/src/mzML.jl
|
||||
# src/mzML.jl
|
||||
|
||||
# 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
|
||||
"""
|
||||
get_spectrum_asset_metadata(stream)
|
||||
|
||||
Parses a `<binaryDataArray>` block within an mzML file to extract metadata
|
||||
for a single data array (e.g., m/z or intensity). It reads CV parameters to
|
||||
determine the data type, compression, and axis type.
|
||||
|
||||
# Arguments
|
||||
- `stream`: An IO stream positioned at the beginning of a `<binaryDataArray>` block.
|
||||
|
||||
# Returns
|
||||
- A `SpectrumAsset` struct containing the parsed metadata, including the binary
|
||||
data offset, encoded length, format, and compression status.
|
||||
"""
|
||||
function get_spectrum_asset_metadata(stream)
|
||||
start_pos = position(stream)
|
||||
|
||||
@ -56,6 +69,10 @@ function get_spectrum_asset_metadata(stream)
|
||||
elseif acc_str == "MS:1000523"
|
||||
data_format = Float64
|
||||
# println("DEBUG: Set format to Float64")
|
||||
# MS:1000572 is a parent term for compression. While not ideal, if present, assume compression.
|
||||
elseif acc_str == "MS:1000572"
|
||||
compression_flag = true
|
||||
# println("DEBUG: Set is_compressed to true based on parent term")
|
||||
elseif acc_str == "MS:1000574"
|
||||
compression_flag = true
|
||||
# println("DEBUG: Set is_compressed to true")
|
||||
@ -80,6 +97,20 @@ function get_spectrum_asset_metadata(stream)
|
||||
end
|
||||
|
||||
# This function is updated to return the generic SpectrumMetadata struct
|
||||
"""
|
||||
parse_spectrum_metadata(stream, offset::Int64)
|
||||
|
||||
Parses an entire `<spectrum>` block from an mzML file, given a starting offset.
|
||||
It extracts the spectrum ID and calls `get_spectrum_asset_metadata` to parse
|
||||
the m/z and intensity array metadata.
|
||||
|
||||
# Arguments
|
||||
- `stream`: An IO stream for the mzML file.
|
||||
- `offset`: The byte offset where the `<spectrum>` block begins.
|
||||
|
||||
# Returns
|
||||
- A `SpectrumMetadata` struct containing the parsed metadata for one spectrum.
|
||||
"""
|
||||
function parse_spectrum_metadata(stream, offset::Int64)
|
||||
seek(stream, offset)
|
||||
|
||||
@ -97,6 +128,18 @@ function parse_spectrum_metadata(stream, offset::Int64)
|
||||
return SpectrumMetadata(0, 0, id, mz_asset, int_asset)
|
||||
end
|
||||
|
||||
"""
|
||||
parse_offset_list(stream)
|
||||
|
||||
Parses the `<index name="spectrum">` block in an indexed mzML file to extract
|
||||
the byte offsets for each spectrum.
|
||||
|
||||
# Arguments
|
||||
- `stream`: An IO stream positioned at the start of the `<index>` block.
|
||||
|
||||
# Returns
|
||||
- A `Vector{Int64}` containing the byte offsets for all spectra.
|
||||
"""
|
||||
function parse_offset_list(stream)
|
||||
offsets = Int64[]
|
||||
offset_regex = r"<offset[^>]*>(\d+)</offset>"
|
||||
@ -121,6 +164,20 @@ function parse_offset_list(stream)
|
||||
end
|
||||
|
||||
# This is the main lazy-loading function for mzML, now returning an MSIData object.
|
||||
"""
|
||||
load_mzml_lazy(file_path::String; cache_size::Int=100)
|
||||
|
||||
Lazily loads an indexed `.mzML` file by parsing only the metadata. It reads the
|
||||
spectrum index from the end of the file to get the offsets of each spectrum,
|
||||
then parses the metadata for each spectrum without loading the binary data.
|
||||
|
||||
# Arguments
|
||||
- `file_path`: The path to the `.mzML` file.
|
||||
- `cache_size`: The number of spectra to hold in an LRU cache for faster access.
|
||||
|
||||
# Returns
|
||||
- An `MSIData` object ready for lazy data access.
|
||||
"""
|
||||
function load_mzml_lazy(file_path::String; cache_size::Int=100)
|
||||
println("DEBUG: Opening file stream for $file_path")
|
||||
stream = open(file_path, "r")
|
||||
@ -181,8 +238,16 @@ end
|
||||
"""
|
||||
LoadMzml(fileName::String)
|
||||
|
||||
Eagerly loads all spectra from a .mzML file.
|
||||
Provided for backward compatibility. Now uses the new MSIData architecture.
|
||||
Eagerly loads all spectra from a .mzML file into memory.
|
||||
This function now uses the new lazy-loading
|
||||
architecture internally but presents the data in the old format.
|
||||
|
||||
# Arguments
|
||||
- `fileName`: The path to the `.mzML` file.
|
||||
|
||||
# Returns
|
||||
- A `2xN` matrix where `N` is the number of spectra. The first row contains
|
||||
m/z arrays and the second row contains intensity arrays.
|
||||
"""
|
||||
function LoadMzml(fileName::String)
|
||||
# Use the lazy loader to get the MSIData object
|
||||
|
||||
@ -31,16 +31,18 @@ using MSI_src
|
||||
# --- 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 TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging_paper_spray/Imaging_paper_spray.mzML"
|
||||
const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML"
|
||||
# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging_paper_spray/Imaging_paper_spray.mzML"
|
||||
# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging prueba Roya 1/Roya.mzML"
|
||||
const SPECTRUM_TO_PLOT = 1 # Which spectrum to plot from the file
|
||||
|
||||
# --- Test Case 2: .mzML + Sync File for Conversion ---
|
||||
# The special .mzML file with one spectrum per pixel.
|
||||
const CONVERSION_SOURCE_MZML = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML"
|
||||
# const CONVERSION_SOURCE_MZML = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML"
|
||||
const CONVERSION_SOURCE_MZML = TEST_MZML_FILE
|
||||
# The corresponding synchronization text file.
|
||||
const CONVERSION_SYNC_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.txt"
|
||||
# const CONVERSION_SYNC_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging prueba Roya 1/Synchro.txt"
|
||||
|
||||
# const CONVERSION_SOURCE_MZML = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging_paper_spray/Imaging_paper_spray.mzML"
|
||||
# const CONVERSION_SYNC_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging_paper_spray/Imaging_paper_spray.txt"
|
||||
@ -50,18 +52,18 @@ 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 = 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"
|
||||
const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging prueba Roya 1/royaimg.imzML"
|
||||
# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging prueba Roya 1/royaimg.imzML"
|
||||
# The m/z value to use for creating an image slice.
|
||||
# const MZ_VALUE_FOR_SLICE = 309.06 # BF
|
||||
# const MZ_VALUE_FOR_SLICE = 896.0 # HR2MSI
|
||||
# const MZ_VALUE_FOR_SLICE = 76.03 # I PS
|
||||
const MZ_VALUE_FOR_SLICE = 473 # ROYA
|
||||
const MZ_VALUE_FOR_SLICE = 313 # ROYA
|
||||
# const MZ_TOLERANCE = 0.1
|
||||
# const MZ_TOLERANCE = 1
|
||||
const MZ_TOLERANCE = 0.5
|
||||
const MZ_TOLERANCE = 0.1
|
||||
|
||||
# Coordinates to plot a specific spectrum from imzML
|
||||
const COORDS_TO_PLOT = (50, 50) # Example coordinates (X, Y)
|
||||
@ -69,8 +71,8 @@ const COORDS_TO_PLOT = (50, 50) # Example coordinates (X, Y)
|
||||
# --- Output Directory ---
|
||||
const RESULTS_DIR = "test/results"
|
||||
|
||||
test1 = false
|
||||
test2 = false
|
||||
test1 = true
|
||||
test2 = true
|
||||
test3 = true
|
||||
|
||||
# ===================================================================
|
||||
@ -152,7 +154,7 @@ function run_test()
|
||||
try
|
||||
# Get the msi data from the mzml
|
||||
println("Plotting a sample spectrum from $TEST_MZML_FILE...")
|
||||
msi_data = OpenMSIData(TEST_MZML_FILE)
|
||||
msi_data = @time OpenMSIData(TEST_MZML_FILE)
|
||||
mz, intensity = GetSpectrum(msi_data, SPECTRUM_TO_PLOT)
|
||||
|
||||
fig = Figure(size = (800, 600))
|
||||
@ -229,10 +231,9 @@ function run_test()
|
||||
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)
|
||||
msi_data = @time OpenMSIData(TEST_IMZML_FILE)
|
||||
x_coord, y_coord = COORDS_TO_PLOT
|
||||
|
||||
# Get the x y coordinate spectrum data
|
||||
@ -275,7 +276,6 @@ function run_test()
|
||||
output_path = joinpath(RESULTS_DIR, "test_imzml_average_spectrum.png")
|
||||
save(output_path, fig)
|
||||
println("SUCCESS: Total spectrum plot saved to $output_path")
|
||||
"""
|
||||
println("No spectrums tested on this try.")
|
||||
catch e
|
||||
println("ERROR during spectrum plotting in Test Case 3: $e")
|
||||
@ -284,7 +284,7 @@ function run_test()
|
||||
# Test the plot_slice function
|
||||
try
|
||||
println("Testing plot_slice function on $TEST_IMZML_FILE...")
|
||||
msi_data = OpenMSIData(TEST_IMZML_FILE)
|
||||
msi_data = @time OpenMSIData(TEST_IMZML_FILE)
|
||||
@time 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
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user