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

integrated MSI_src API into the UI, fixed interface elements, made some processes asynchronous, modified spectrum plots, added sum spectrum, removed use of obsolete julia_mzML_imzML library

Browse Source
This commit is contained in:
Pixelguy14 2025-10-02 13:28:15 -06:00
parent cd1115368f
commit 22c8d9637b
10 changed files with 1076 additions and 708 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Manifest.toml
View File

@ -2,7 +2,7 @@
julia_version = "1.11.7" julia_version = "1.11.7"
manifest_format = "2.0" manifest_format = "2.0"
project_hash = "896c16630591d2732da8e3315c4b64c618fd5ebb" project_hash = "0f3959248d174e05d23b9e9267b00c6d2ead256d"
[[deps.ATK_jll]] [[deps.ATK_jll]]
deps = ["Artifacts", "Glib_jll", "JLLWrappers", "Libdl"] deps = ["Artifacts", "Glib_jll", "JLLWrappers", "Libdl"]

2
Project.toml
View File

@ -19,6 +19,8 @@ NativeFileDialog = "e1fe445b-aa65-4df4-81c1-2041507f0fd4"
NaturalSort = "c020b1a1-e9b0-503a-9c33-f039bfc54a85" NaturalSort = "c020b1a1-e9b0-503a-9c33-f039bfc54a85"
PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5" PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
StipplePlotly = "ec984513-233d-481d-95b0-a3b58b97af2b" StipplePlotly = "ec984513-233d-481d-95b0-a3b58b97af2b"
julia_mzML_imzML = "38eb50d3-2fb6-4afa-992a-964ed8562ed9" julia_mzML_imzML = "38eb50d3-2fb6-4afa-992a-964ed8562ed9"

805
app.jl
View File

@ -6,7 +6,8 @@ using Libz
using PlotlyBase using PlotlyBase
using CairoMakie using CairoMakie
using Colors using Colors
using julia_mzML_imzML # using julia_mzML_imzML
using MSI_src # Import the new MSIData library
using Statistics using Statistics
using NaturalSort using NaturalSort
using Images using Images
@ -14,6 +15,10 @@ using LinearAlgebra
using NativeFileDialog # Opens the file explorer depending on the OS using NativeFileDialog # Opens the file explorer depending on the OS
using StipplePlotly using StipplePlotly
using Base.Filesystem: mv # To rename files in the system using Base.Filesystem: mv # To rename files in the system
# 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
include("./julia_imzML_visual.jl") include("./julia_imzML_visual.jl")
@genietools @genietools
@ -57,6 +62,7 @@ include("./julia_imzML_visual.jl")
@in compareBtn=false # To open dialog @in compareBtn=false # To open dialog
@in createMeanPlot=false # To generate mean spectrum plot @in createMeanPlot=false # To generate mean spectrum plot
@in createXYPlot=false # To generate an spectrum plot according to the xy values inputed @in createXYPlot=false # To generate an spectrum plot according to the xy values inputed
@in createSumPlot=false # To generate a sum of all the spectrum plots
@in image3dPlot=false # To generate 3d plot based on current image @in image3dPlot=false # To generate 3d plot based on current image
@in triq3dPlot=false # To generate 3d plot based on current triq image @in triq3dPlot=false # To generate 3d plot based on current triq image
@in imageCPlot=false # To generate contour plots of current image @in imageCPlot=false # To generate contour plots of current image
@ -110,10 +116,12 @@ include("./julia_imzML_visual.jl")
@out msgimgComp="" @out msgimgComp=""
@out msgtriqComp="" @out msgtriqComp=""
# Centralized MSIData object
@out msi_data::Union{MSIData, Nothing} = nothing
# Saves the route where imzML and mzML files are located # Saves the route where imzML and mzML files are located
@out full_route="" @out full_route=""
@out full_routeMz=""
@out full_routeMz2=""
# For the creation of images with a more specific mass charge # For the creation of images with a more specific mass charge
@out text_nmass="" @out text_nmass=""
@ -145,10 +153,12 @@ include("./julia_imzML_visual.jl")
layoutImg=PlotlyBase.Layout( layoutImg=PlotlyBase.Layout(
xaxis=PlotlyBase.attr( xaxis=PlotlyBase.attr(
visible=false, visible=false,
scaleanchor="y" scaleanchor="y",
range=[0, 0]
), ),
yaxis=PlotlyBase.attr( yaxis=PlotlyBase.attr(
visible=false visible=false,
range=[0, 0]
), ),
margin=attr(l=0,r=0,t=0,b=0,pad=0) margin=attr(l=0,r=0,t=0,b=0,pad=0)
) )
@ -168,10 +178,10 @@ include("./julia_imzML_visual.jl")
# Interface Plot Spectrum # Interface Plot Spectrum
layoutSpectra=PlotlyBase.Layout( layoutSpectra=PlotlyBase.Layout(
title="Spectrum plot", title="Spectrum plot",
hovermode="closest",
xaxis=PlotlyBase.attr( xaxis=PlotlyBase.attr(
title="<i>m/z</i>", title="<i>m/z</i>",
showgrid=true, showgrid=true
tickformat = ".3g"
), ),
yaxis=PlotlyBase.attr( yaxis=PlotlyBase.attr(
title="Intensity", title="Intensity",
@ -249,319 +259,305 @@ include("./julia_imzML_visual.jl")
# Reactive handlers watch a variable and execute a block of code when its value changes # Reactive handlers watch a variable and execute a block of code when its value changes
# The onbutton handler will set the variable to false after the block is executed # The onbutton handler will set the variable to false after the block is executed
@onbutton btnSearch begin @onbutton btnSearch @time begin
full_route=pick_file(; filterlist="imzML,imzml,mzML,mzml") # This part is synchronous and blocking, which is unavoidable
msg="" picked_route = pick_file(; filterlist="imzML,imzml,mzML,mzml")
if full_route != ""
if endswith(full_route, "imzml") if isempty(picked_route)
alt_route=replace(full_route, r"\.[^.]*$" => ".imzML") msg = "No file selected."
mv(full_route, alt_route) warning_msg = true
full_route=alt_route return
# to detect if there's an mzml named wrong
alt_routeMz=replace(full_route, r"\.[^.]*$" => ".mzml")
if isfile(alt_routeMz)
alt_routeMz2=replace(alt_routeMz, r"\.[^.]*$" => ".mzML")
mv(alt_routeMz, alt_routeMz2)
alt_routeMz=alt_routeMz2
end
end
if endswith(full_route, "mzml")
alt_route=replace(full_route, r"\.[^.]*$" => ".mzML")
mv(full_route, alt_route)
full_route = alt_route
# to detect if there's an imzml named wrong
alt_routeMz=replace(full_route, r"\.[^.]*$" => ".imzml")
if isfile(alt_routeMz)
alt_routeMz2=replace(alt_routeMz, r"\.[^.]*$" => ".imzML")
mv(alt_routeMz, alt_routeMz2)
alt_routeMz=alt_routeMz2
end
end
end end
if full_route==""
msg="No file selected"
warning_msg=true
btnStartDisable=true
btnSpectraDisable=true
SpectraEnabled=false
else
if endswith(full_route, "imzML") # Case if the file loaded is imzML
btnStartDisable=false
btnPlotDisable=false
# Splitting the route with regex from imzml to mzml so the plotting can work
full_routeMz=replace(full_route, r"\.[^.]*$" => ".mzML")
if isfile(full_routeMz)
# Start the mean spectrum creation on loading
progressSpectraPlot=true
btnPlotDisable=true
btnStartDisable=true
msg="Loading mean spectrum plot..."
sTime=time()
plotdata, plotlayout, xSpectraMz, ySpectraMz=meanSpectrumPlot(full_routeMz)
selectedTab="tab2"
progressSpectraPlot=false
btnPlotDisable=false
# We enable coord search and spectra plot creation, also re-enable main process # UI updates immediately
btnSpectraDisable=false progress = true
SpectraEnabled=true msg = "Opening file: $(basename(picked_route))..."
btnStartDisable=false
@async begin
fTime=time() try
eTime=round(fTime-sTime,digits=3) # --- Normalize file extension and path ---
msg="Plot loaded in $(eTime) seconds" if endswith(picked_route, "imzml")
full_route = replace(picked_route, r"\.imzml$"i => ".imzML")
mv(picked_route, full_route, force=true)
elseif endswith(picked_route, "mzml")
full_route = replace(picked_route, r"\.mzml$"i => ".mzML")
mv(picked_route, full_route, force=true)
else else
# If there's no mzML file, we deny access again full_route = picked_route
btnSpectraDisable=true
SpectraEnabled=false
end end
else # Case if the file loaded is mzML
full_routeMz=full_route
btnSpectraDisable=false
SpectraEnabled=true
if isfile(full_routeMz)
# Start the mean spectrum creation on loading
progressSpectraPlot=true
btnPlotDisable=true
btnStartDisable=true
msg="Loading mean spectrum plot..."
sTime=time()
plotdata, plotlayout, xSpectraMz, ySpectraMz=meanSpectrumPlot(full_routeMz)
selectedTab="tab2"
progressSpectraPlot=false
btnPlotDisable=false
# We enable coord search and spectra plot creation # --- Load data using the new MSIData library ---
btnSpectraDisable=false sTime = time()
SpectraEnabled=true msi_data = OpenMSIData(full_route)
w, h = msi_data.image_dims
imgWidth, imgHeight = w > 0 ? (w, h) : (500, 500)
fTime=time() fTime = time()
eTime=round(fTime-sTime,digits=3) eTime = round(fTime - sTime, digits=3)
msg="Plot loaded in $(eTime) seconds" msg = "File loaded in $(eTime) seconds. Calculating total spectrum..."
end
# Splitting the route the same way # Enable UI controls
full_route=replace(full_route, r"\.[^.]*$" => ".imzML") btnStartDisable = !(msi_data.source isa ImzMLSource)
if isfile(full_route) btnPlotDisable = false
btnStartDisable=false btnSpectraDisable = false
else SpectraEnabled = true
btnStartDisable=true
full_route=full_routeMz # --- Automatically generate and display the sum spectrum plot ---
end # This is still part of the same async task
progressSpectraPlot = true
sTime = time()
plotdata, plotlayout, xSpectraMz, ySpectraMz = sumSpectrumPlot(msi_data)
selectedTab = "tab2"
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"
warning_msg = true
btnStartDisable = true
btnSpectraDisable = true
SpectraEnabled = false
@error "File loading failed" exception=(e, catch_backtrace())
finally
# This now correctly runs after everything is finished
progress = false
progressSpectraPlot = false
end end
xCoord=0
yCoord=0
end end
end end
@onbutton mainProcess begin @onbutton mainProcess @time begin
progress=true # Start progress button animation # UI updates immediately
btnStartDisable=true # We disable the button to avoid multiple requests progress = true
btnPlotDisable=true btnStartDisable = true
btnSpectraDisable=true btnPlotDisable = true
text_nmass=replace(string(Nmass), "." => "_") btnSpectraDisable = true
sTime=time()
if isfile(full_route) && Nmass > 0 && Tol > 0 && Tol <=1 && colorLevel > 1 && colorLevel < 257 @async begin
msg="File exists, Nmass=$(Nmass) Tol=$(Tol). Loading file will begin, please be patient."
try try
spectra=LoadImzml(full_route) text_nmass = replace(string(Nmass), "." => "_")
msg="File loaded. Creating spectra with the specific mass and tolerance, please be patient." sTime = time()
slice=GetMzSliceJl(spectra,Nmass,Tol)
fig=CairoMakie.Figure(size=(150, 250)) # Container if msi_data === nothing || !(msi_data.source isa ImzMLSource)
# Append a query string to force the image to refresh msg = "No .imzML file loaded or selected file is not an .imzML. Please select a valid file."
timestamp=string(time_ns()) warning_msg = true
if triqEnabled # If we have TrIQ elseif Nmass > 0 && Tol > 0 && Tol <= 1 && colorLevel > 1 && colorLevel < 257
if triqProb < 0.8 || triqProb > 1 msg = "Creating image for Nmass=$(Nmass) Tol=$(Tol). Please be patient."
msg="Incorrect TrIQ values, please adjust accordingly and try again."
warning_msg=true
else
image_path=joinpath("./public", "TrIQ_$(text_nmass).bmp")
valid_slice=false
while Tol <= 1.0 && !valid_slice
try
slice=GetMzSliceJl(spectra, Nmass, Tol)
sliceTriq=TrIQ(slice, colorLevel, triqProb)
if MFilterEnabled # If the Median filter is ON
sliceTriq=medianFilterjl(sliceTriq)
end
valid_slice=true
catch e
msg="Warning: insufficient tolerance, inputs modified to allow the creation of an image regardless=$Tol: $e"
Tol += 0.1
end
end
sliceTriq=reverse(sliceTriq, dims=2)
SaveBitmapCl(joinpath("public", "TrIQ_$(text_nmass).bmp"),sliceTriq,ViridisPalette)
# Use timestamp to refresh image interface container
imgIntT="/TrIQ_$(text_nmass).bmp?t=$(timestamp)"
plotdataImgT, plotlayoutImgT, imgWidth, imgHeight=loadImgPlot(imgIntT)
# Get current image
current_triq="TrIQ_$(text_nmass).bmp"
msgtriq="TrIQ image with the Nmass of $(replace(text_nmass, "_" => "."))"
# Create colorbar
bound=julia_mzML_imzML.GetOutlierThres(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)
colorbarT="/colorbar_TrIQ_$(text_nmass).png?t=$(timestamp)"
# Get current colorbar
current_col_triq="colorbar_TrIQ_$(text_nmass).png"
# We update the directory to include the new placed images.
triq_bmp=sort(filter(filename -> startswith(filename, "TrIQ_") && endswith(filename, ".bmp"), readdir("public")),lt=natural)
col_triq_png=sort(filter(filename -> startswith(filename, "colorbar_TrIQ_") && endswith(filename, ".png"), readdir("public")),lt=natural)
fTime=time()
eTime=round(fTime-sTime,digits=3)
msg="The file has been created in $(eTime) seconds successfully inside the 'public' folder of the app"
selectedTab="tab1"
end
else # If we don't use TrIQ
image_path=joinpath("./public", "MSI_$(text_nmass).bmp")
try try
sliceQuant=IntQuantCl(slice,Int(colorLevel-1)) # Use the new get_mz_slice with the centralized MSIData object
if MFilterEnabled # If the Median filter is ON slice = get_mz_slice(msi_data, Nmass, Tol)
sliceQuant=medianFilterjl(sliceQuant) fig = CairoMakie.Figure(size=(150, 250)) # Container
timestamp = string(time_ns())
if triqEnabled # If we have TrIQ
if triqProb < 0.8 || triqProb > 1
msg = "Incorrect TrIQ values, please adjust accordingly and try again."
warning_msg = true
else
sliceTriq = TrIQ(slice, colorLevel, triqProb)
if MFilterEnabled
sliceTriq = round.(UInt8, median_filter(sliceTriq))
end
sliceTriq = reverse(sliceTriq, dims=2)
save_bitmap(joinpath("public", "TrIQ_$(text_nmass).bmp"), sliceTriq, ViridisPalette)
imgIntT = "/TrIQ_$(text_nmass).bmp?t=$(timestamp)"
plotdataImgT, plotlayoutImgT, imgWidth, imgHeight = loadImgPlot(imgIntT)
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)
colorbarT = "/colorbar_TrIQ_$(text_nmass).png?t=$(timestamp)"
current_col_triq = "colorbar_TrIQ_$(text_nmass).png"
triq_bmp = sort(filter(filename -> startswith(filename, "TrIQ_") && endswith(filename, ".bmp"), readdir("public")), lt=natural)
col_triq_png = sort(filter(filename -> startswith(filename, "colorbar_TrIQ_") && endswith(filename, ".png"), readdir("public")), lt=natural)
fTime = time()
eTime = round(fTime - sTime, digits=3)
msg = "The TrIQ image has been created in $(eTime) seconds successfully inside the 'public' folder of the app"
selectedTab = "tab1"
end
else # If we don't use TrIQ
sliceQuant = quantize_intensity(slice, colorLevel)
if MFilterEnabled
sliceQuant = round.(UInt8, median_filter(sliceQuant))
end
sliceQuant = reverse(sliceQuant, dims=2)
save_bitmap(joinpath("public", "MSI_$(text_nmass).bmp"), sliceQuant, ViridisPalette)
imgInt = "/MSI_$(text_nmass).bmp?t=$(timestamp)"
plotdataImg, plotlayoutImg, imgWidth, imgHeight = loadImgPlot(imgInt)
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 = "/colorbar_MSI_$(text_nmass).png?t=$(timestamp)"
current_col_msi = "colorbar_MSI_$(text_nmass).png"
msi_bmp = sort(filter(filename -> startswith(filename, "MSI_") && endswith(filename, ".bmp"), readdir("public")), lt=natural)
col_msi_png = sort(filter(filename -> startswith(filename, "colorbar_MSI_") && endswith(filename, ".png"), readdir("public")), lt=natural)
selectedTab = "tab0"
fTime = time()
eTime = round(fTime - sTime, digits=3)
msg = "The image has been created in $(eTime) seconds successfully inside the 'public' folder of the app"
end end
catch e catch e
msg="Warning: $e" msg = "There was an error creating the image: $e"
warning_msg = true
@error "Image creation failed" exception=(e, catch_backtrace())
end end
sliceQuant=reverse(sliceQuant, dims=2) else
SaveBitmapCl(joinpath("public", "MSI_$(text_nmass).bmp"),sliceQuant,ViridisPalette) msg = "Invalid parameters. Nmass, Tol, or colorLevel are incorrect."
# Use timestamp to refresh image interface container warning_msg = true
imgInt="/MSI_$(text_nmass).bmp?t=$(timestamp)" @error msg
plotdataImg, plotlayoutImg, imgWidth, imgHeight=loadImgPlot(imgInt)
# Get current image
current_msi="MSI_$(text_nmass).bmp"
msgimg="Image with the Nmass of $(replace(text_nmass, "_" => "."))"
# Create colorbar
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="/colorbar_MSI_$(text_nmass).png?t=$(timestamp)"
# Get current colorbar
current_col_msi="colorbar_MSI_$(text_nmass).png"
# We update the directory to include the new placed images.
msi_bmp=sort(filter(filename -> startswith(filename, "MSI_") && endswith(filename, ".bmp"), readdir("public")),lt=natural)
col_msi_png=sort(filter(filename -> startswith(filename, "colorbar_MSI_") && endswith(filename, ".png"), readdir("public")),lt=natural)
selectedTab="tab0"
fTime=time()
eTime=round(fTime-sTime,digits=3)
msg="The file has been created in $(eTime) seconds successfully inside the 'public' folder of the app"
end end
catch e finally
msg="There was an error loading the ImzML file, please verify the file accordingly and try again. $(e)" # This block will always run at the end of the async task
warning_msg=true GC.gc()
if Sys.islinux()
ccall(:malloc_trim, Int32, (Int32,), 0)
end
btnStartDisable = false
btnPlotDisable = false
btnOpticalDisable = false
progress = false
btnSpectraDisable = false
SpectraEnabled = true
end end
else
msg="File does not exist or a parameter is incorrect, please try again."
warning_msg=true
end
GC.gc() # Trigger garbage collection
if Sys.islinux()
ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure julia returns the freed memory to OS
end
btnStartDisable=false
btnPlotDisable=false
btnOpticalDisable=false
progress=false
if isfile(full_routeMz)
# We enable coord search and spectra plot creation
btnSpectraDisable=false
SpectraEnabled=true
end end
end end
@onbutton createMeanPlot begin @onbutton createMeanPlot begin
msg="Mean spectrum plot selected" if msi_data === nothing
sTime=time() msg = "No data loaded. Please select a file first."
if isfile(full_routeMz) # Check if the file exists warning_msg = true
progressSpectraPlot=true return
btnPlotDisable=true end
btnStartDisable=true
msg="Loading plot..." # UI updates immediately
plotdata, plotlayout, xSpectraMz, ySpectraMz=meanSpectrumPlot(full_routeMz) progressSpectraPlot = true
GC.gc() # Trigger garbage collection btnPlotDisable = true
if Sys.islinux() btnStartDisable = true
ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure julia returns the freed memory to OS msg = "Loading plot..."
@async begin
try
sTime = time()
plotdata, plotlayout, xSpectraMz, ySpectraMz = meanSpectrumPlot(msi_data)
selectedTab = "tab2"
fTime = time()
eTime = round(fTime - sTime, digits=3)
msg = "Plot loaded in $(eTime) seconds"
catch e
msg = "Could not generate mean spectrum plot: $e"
warning_msg = true
@error "Mean spectrum plotting failed" exception=(e, catch_backtrace())
finally
# This runs after the async task is finished
progressSpectraPlot = false
btnPlotDisable = false
btnSpectraDisable = false
if msi_data !== nothing && msi_data.source isa ImzMLSource
btnStartDisable = false
end
end end
selectedTab="tab2"
fTime=time()
eTime=round(fTime-sTime,digits=3)
msg="Plot loaded in $(eTime) seconds"
else
msg="there was an error with the mzML, please try again"
warning_msg=true
end end
progressSpectraPlot=false end
btnPlotDisable=false
if endswith(full_route, "imzML") @onbutton createSumPlot begin
btnStartDisable=false if msi_data === nothing
msg = "No data loaded. Please select a file first."
warning_msg = true
return
end end
if isfile(full_routeMz)
# We enable coord search and spectra plot creation # UI updates immediately
btnSpectraDisable=false progressSpectraPlot = true
SpectraEnabled=true btnPlotDisable = true
btnStartDisable = true
msg = "Loading total spectrum plot..."
@async begin
try
sTime = time()
plotdata, plotlayout, xSpectraMz, ySpectraMz = sumSpectrumPlot(msi_data)
selectedTab = "tab2"
fTime = time()
eTime = round(fTime - sTime, digits=3)
msg = "Total plot loaded in $(eTime) seconds"
catch e
msg = "Could not generate total spectrum plot: $e"
warning_msg = true
@error "Total spectrum plotting failed" exception=(e, catch_backtrace())
finally
# This runs after the async task is finished
progressSpectraPlot = false
btnPlotDisable = false
btnSpectraDisable = false
if msi_data !== nothing && msi_data.source isa ImzMLSource
btnStartDisable = false
end
end
end end
end end
@onbutton createXYPlot begin @onbutton createXYPlot begin
msg="XY spectrum plot selected" if msi_data === nothing
sTime=time() msg = "No data loaded. Please select a file first."
if isfile(full_routeMz) # Check if the file exists warning_msg = true
progressSpectraPlot=true return
btnStartDisable=true
btnPlotDisable=true
btnSpectraDisable=true
msg="Loading plot..."
spectraMz=LoadMzml(full_routeMz)
layoutSpectra=PlotlyBase.Layout(
title="($xCoord, $yCoord) Specific spectrum plot",
xaxis=PlotlyBase.attr(
title="<i>m/z</i>",
showgrid=true
),
yaxis=PlotlyBase.attr(
title="Intensity",
showgrid=true
),
autosize=false,
margin=attr(l=0,r=0,t=120,b=0,pad=0)
)
if xCoord < 1
xCoord=1
elseif xCoord > imgWidth
xCoord=imgWidth
end
if yCoord > -1
yCoord=-1
elseif yCoord < -imgHeight
yCoord=-imgHeight
end
xSpectraMz=spectraMz[1,abs(xCoord)]
ySpectraMz=spectraMz[2,abs(yCoord)]
traceSpectra=PlotlyBase.scatter(x=xSpectraMz, y=ySpectraMz, mode="lines")
plotdata=[traceSpectra] # We add the data from spectra to the plot
plotlayout=layoutSpectra
GC.gc() # Trigger garbage collection
if Sys.islinux()
ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure julia returns the freed memory to OS
end
selectedTab="tab2"
fTime=time()
eTime=round(fTime-sTime,digits=3)
msg="Plot loaded in $(eTime) seconds"
else
msg="there was an error with the mzML or the coordenates, please try again"
warning_msg=true
end end
progressSpectraPlot=false
btnPlotDisable=false # UI updates immediately
if endswith(full_route, "imzML") progressSpectraPlot = true
btnStartDisable=false btnStartDisable = true
end btnPlotDisable = true
if isfile(full_routeMz) btnSpectraDisable = true
# We enable coord search and spectra plot creation msg = "Loading plot..."
btnSpectraDisable=false
SpectraEnabled=true @async begin
try
sTime = time()
# The UI uses negative Y values, so we adjust before calling the plot function
y = yCoord < 0 ? abs(yCoord) : yCoord
plotdata, plotlayout, xSpectraMz, ySpectraMz = xySpectrumPlot(msi_data, xCoord, y, imgWidth, imgHeight)
# Update UI coordinates
xCoord = plotlayout.title == "Spectrum #$(xCoord)" ? xCoord : clamp(xCoord, 1, imgWidth)
yCoord = plotlayout.title == "Spectrum #$(xCoord)" ? 0 : -clamp(y, 1, imgHeight)
selectedTab = "tab2"
fTime = time()
eTime = round(fTime - sTime, digits=3)
msg = "Plot loaded in $(eTime) seconds"
catch e
msg = "Could not retrieve spectrum: $e"
warning_msg = true
@error "Spectrum plotting failed" exception=(e, catch_backtrace())
finally
# This runs after the async task is finished
progressSpectraPlot = false
btnPlotDisable = false
btnSpectraDisable = false
if msi_data !== nothing && msi_data.source isa ImzMLSource
btnStartDisable = false
end
end
end end
end end
@ -796,14 +792,21 @@ include("./julia_imzML_visual.jl")
cleaned_imgInt=replace(imgInt, r"\?.*" => "") cleaned_imgInt=replace(imgInt, r"\?.*" => "")
cleaned_imgInt=lstrip(cleaned_imgInt, '/') cleaned_imgInt=lstrip(cleaned_imgInt, '/')
var=joinpath( "./public", cleaned_imgInt ) var=joinpath( "./public", cleaned_imgInt )
sTime=time()
if isfile(var) if !isfile(var)
progressPlot=true msg="Image could not be 3d plotted"
btnPlotDisable=true warning_msg=true
btnStartDisable=true return
btnSpectraDisable=true end
progressPlot=true
btnPlotDisable=true
btnStartDisable=true
btnSpectraDisable=true
@async begin
try try
sTime=time()
plotdata3d, plotlayout3d=loadSurfacePlot(imgInt) plotdata3d, plotlayout3d=loadSurfacePlot(imgInt)
GC.gc() # Trigger garbage collection GC.gc() # Trigger garbage collection
if Sys.islinux() if Sys.islinux()
@ -816,34 +819,39 @@ include("./julia_imzML_visual.jl")
catch e catch e
msg="Failed to load and process image: $e" msg="Failed to load and process image: $e"
warning_msg=true warning_msg=true
finally
progressPlot=false
btnPlotDisable=false
btnStartDisable=false
if msi_data !== nothing
# We enable coord search and spectra plot creation
btnSpectraDisable=false
SpectraEnabled=true
end
end end
else
msg="Image could not be 3d plotted"
warning_msg=true
end end
progressPlot=false end # 3d plot for TrIQ
btnPlotDisable=false
btnStartDisable=false
if isfile(full_routeMz)
# We enable coord search and spectra plot creation
btnSpectraDisable=false
SpectraEnabled=true
end
end
# 3d plot for TrIQ
@onbutton triq3dPlot begin @onbutton triq3dPlot begin
msg="TrIQ 3D plot selected" msg="TrIQ 3D plot selected"
cleaned_imgIntT=replace(imgIntT, r"\?.*" => "") cleaned_imgIntT=replace(imgIntT, r"\?.*" => "")
cleaned_imgIntT=lstrip(cleaned_imgIntT, '/') cleaned_imgIntT=lstrip(cleaned_imgIntT, '/')
var=joinpath( "./public", cleaned_imgIntT ) var=joinpath( "./public", cleaned_imgIntT )
sTime=time()
if isfile(var) if !isfile(var)
progressPlot=true msg="Image could not be 3d plotted"
btnPlotDisable=true warning_msg=true
btnStartDisable=true return
btnSpectraDisable=true end
progressPlot=true
btnPlotDisable=true
btnStartDisable=true
btnSpectraDisable=true
@async begin
try try
sTime=time()
plotdata3d, plotlayout3d=loadSurfacePlot(imgIntT) plotdata3d, plotlayout3d=loadSurfacePlot(imgIntT)
GC.gc() # Trigger garbage collection GC.gc() # Trigger garbage collection
if Sys.islinux() if Sys.islinux()
@ -856,18 +864,16 @@ include("./julia_imzML_visual.jl")
catch e catch e
msg="Failed to load and process image: $e" msg="Failed to load and process image: $e"
warning_msg=true warning_msg=true
finally
progressPlot=false
btnPlotDisable=false
btnStartDisable=false
if msi_data !== nothing
# We enable coord search and spectra plot creation
btnSpectraDisable=false
SpectraEnabled=true
end
end end
else
msg="Image could not be 3d plotted"
warning_msg=true
end
progressPlot=false
btnPlotDisable=false
btnStartDisable=false
if isfile(full_routeMz)
# We enable coord search and spectra plot creation
btnSpectraDisable=false
SpectraEnabled=true
end end
end end
@ -877,15 +883,21 @@ include("./julia_imzML_visual.jl")
cleaned_imgInt=replace(imgInt, r"\?.*" => "") cleaned_imgInt=replace(imgInt, r"\?.*" => "")
cleaned_imgInt=lstrip(cleaned_imgInt, '/') cleaned_imgInt=lstrip(cleaned_imgInt, '/')
var=joinpath("./public", cleaned_imgInt) var=joinpath("./public", cleaned_imgInt)
sTime=time()
if !isfile(var)
if isfile(var) msg="Image could not be 2D plotted"
progressPlot=true warning_msg=true
btnPlotDisable=true return
btnStartDisable=true end
btnSpectraDisable=true
progressPlot=true
btnPlotDisable=true
btnStartDisable=true
btnSpectraDisable=true
@async begin
try try
img=load(var) sTime=time()
plotdataC,plotlayoutC=loadContourPlot(imgInt) plotdataC,plotlayoutC=loadContourPlot(imgInt)
GC.gc() # Trigger garbage collection GC.gc() # Trigger garbage collection
if Sys.islinux() if Sys.islinux()
@ -898,18 +910,16 @@ include("./julia_imzML_visual.jl")
catch e catch e
msg="Failed to load and process image: $e" msg="Failed to load and process image: $e"
warning_msg=true warning_msg=true
finally
progressPlot=false
btnPlotDisable=false
btnStartDisable=false
if msi_data !== nothing
# We enable coord search and spectra plot creation
btnSpectraDisable=false
SpectraEnabled=true
end
end end
else
msg="Image could not be 2D plotted"
warning_msg=true
end
progressPlot=false
btnPlotDisable=false
btnStartDisable=false
if isfile(full_routeMz)
# We enable coord search and spectra plot creation
btnSpectraDisable=false
SpectraEnabled=true
end end
end end
# Contour 2d plot for TrIQ # Contour 2d plot for TrIQ
@ -918,15 +928,21 @@ include("./julia_imzML_visual.jl")
cleaned_imgIntT=replace(imgIntT, r"\?.*" => "") cleaned_imgIntT=replace(imgIntT, r"\?.*" => "")
cleaned_imgIntT=lstrip(cleaned_imgIntT, '/') cleaned_imgIntT=lstrip(cleaned_imgIntT, '/')
var=joinpath("./public", cleaned_imgIntT) var=joinpath("./public", cleaned_imgIntT)
sTime=time()
if !isfile(var)
if isfile(var) msg="Image could not be 2D plotted"
progressPlot=true warning_msg=true
btnPlotDisable=true return
btnStartDisable=true end
btnSpectraDisable=true
progressPlot=true
btnPlotDisable=true
btnStartDisable=true
btnSpectraDisable=true
@async begin
try try
img=load(var) sTime=time()
plotdataC,plotlayoutC=loadContourPlot(imgIntT) plotdataC,plotlayoutC=loadContourPlot(imgIntT)
GC.gc() # Trigger garbage collection GC.gc() # Trigger garbage collection
if Sys.islinux() if Sys.islinux()
@ -939,18 +955,16 @@ include("./julia_imzML_visual.jl")
catch e catch e
msg="Failed to load and process image: $e" msg="Failed to load and process image: $e"
warning_msg=true warning_msg=true
finally
progressPlot=false
btnPlotDisable=false
btnStartDisable=false
if msi_data !== nothing
# We enable coord search and spectra plot creation
btnSpectraDisable=false
SpectraEnabled=true
end
end end
else
msg="Image could not be 2D plotted"
warning_msg=true
end
progressPlot=false
btnPlotDisable=false
btnStartDisable=false
if isfile(full_routeMz)
# We enable coord search and spectra plot creation
btnSpectraDisable=false
SpectraEnabled=true
end end
end end
@ -961,48 +975,39 @@ include("./julia_imzML_visual.jl")
# To include a visualization in the spectrum plot indicating where is the selected mass # To include a visualization in the spectrum plot indicating where is the selected mass
@onchange Nmass begin @onchange Nmass begin
if !isempty(xSpectraMz) if !isempty(xSpectraMz)
traceSpectra=PlotlyBase.scatter(x=xSpectraMz, y=ySpectraMz, mode="lines",name="Spectra",showlegend=false) # Use a stem plot for the main spectrum for consistency
trace2=PlotlyBase.scatter(x=[Nmass, Nmass],y=[0, maximum(ySpectraMz)],mode="lines",line=attr(color="red", width=0.5),name="<i>m/z</i> selected",showlegend=false) traceSpectra = PlotlyBase.stem(x=xSpectraMz, y=ySpectraMz, marker=attr(size=1, color="blue", opacity=0.5), name="Spectrum", hoverinfo="x", hovertemplate="<b>m/z</b>: %{x:.4f}<extra></extra>", showlegend=false)
plotdata=[traceSpectra,trace2] # We add the data from spectra and the red line to the plot
# Keep this as a scatter plot to draw the vertical line
trace2 = PlotlyBase.scatter(x=[Nmass, Nmass], y=[0, maximum(ySpectraMz)], mode="lines", line=attr(color="red", width=0.5), name="<i>m/z</i> selected", showlegend=false)
plotdata = [traceSpectra, trace2]
end end
end end
# Event detection for clicking on the images # Event detection for clicking on the images
@onchange data_click begin @onchange data_click begin
if selectedTab == "tab1" || selectedTab == "tab0"
# This is for the image heatmaps
cursor_data = data_click["cursor"]
x = Int32(round(cursor_data["x"]))
y = Int32(round(cursor_data["y"])) # y is negative in the UI
# Update the reactive coordinates, which will trigger the crosshair update
xCoord = clamp(x, 1, imgWidth)
yCoord = clamp(y, -imgHeight, -1)
end
end
@onchange xCoord, yCoord begin
if selectedTab == "tab1" if selectedTab == "tab1"
cursor_data=data_click["cursor"] plotdataImgT = filter(trace -> !(get(trace, :name, "") in ["Line X", "Line Y"]), plotdataImgT)
xCoord=Int32(round(cursor_data["x"])) trace1, trace2 = crossLinesPlot(xCoord, yCoord, imgWidth, -imgHeight)
yCoord=Int32(round(cursor_data["y"])) plotdataImgT = append!(plotdataImgT, [trace1, trace2])
if xCoord < 1
xCoord=1
elseif xCoord > imgWidth
xCoord=imgWidth
end
if yCoord > -1
yCoord=-1
elseif yCoord < -imgHeight
yCoord=-imgHeight
end # Get the x and y values from the click of the cursor and make sure they don't exceed image proportions
plotdataImgT=filter(trace -> !(get(trace, :name, "") in ["Line X", "Line Y"]), plotdataImgT)
trace1, trace2=crossLinesPlot(xCoord, yCoord, imgWidth, -imgHeight)
plotdataImgT=append!(plotdataImgT, [trace1, trace2])
elseif selectedTab == "tab0" elseif selectedTab == "tab0"
cursor_data=data_click["cursor"] plotdataImg = filter(trace -> !(get(trace, :name, "") in ["Line X", "Line Y", "Optical"]), plotdataImg)
xCoord=Int32(round(cursor_data["x"])) trace1, trace2 = crossLinesPlot(xCoord, yCoord, imgWidth, -imgHeight)
yCoord=Int32(round(cursor_data["y"])) plotdataImg = append!(plotdataImg, [trace1, trace2])
if xCoord < 1
xCoord=1
elseif xCoord > imgWidth
xCoord=imgWidth
end
if yCoord > -1
yCoord=-1
elseif yCoord < -imgHeight
yCoord=-imgHeight
end # Get the x and y values from the click of the cursor and make sure they don't exceed image proportions
plotdataImg=filter(trace -> !(get(trace, :name, "") in ["Line X", "Line Y","Optical"]), plotdataImg)
trace1, trace2=crossLinesPlot(xCoord, yCoord, imgWidth, -imgHeight)
plotdataImg=append!(plotdataImg, [trace1, trace2])
end end
end end

5
app.jl.html
View File

@ -69,6 +69,11 @@
<q-item-label>Mean spectrum plot</q-item-label> <q-item-label>Mean spectrum plot</q-item-label>
</q-item-section> </q-item-section>
</q-item> </q-item>
<q-item clickable v-close-popup v-on:click="createSumPlot=true">
<q-item-section>
<q-item-label>Sum Spectrum plot</q-item-label>
</q-item-section>
</q-item>
<q-item clickable v-close-popup v-on:click="createXYPlot=true"> <q-item clickable v-close-popup v-on:click="createXYPlot=true">
<q-item-section> <q-item-section>
<q-item-label>Spectrum plot (X,Y)</q-item-label> <q-item-label>Spectrum plot (X,Y)</q-item-label>

248
julia_imzML_visual.jl
View File

@ -1,84 +1,3 @@
# IntQuantCl is originally a function of julia mzMl imzML with the adition
# of altering the scale of the colors according to colorlevel.
function IntQuantCl( slice , colorLevel)
# Compute scale factor for amplitude discretization
lower = minimum( slice )
scale = colorLevel / 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
# SaveBitmap originally a function of the mzML imzML library in julia,
# This function dinamically adjust the color palete adjusting to the ammount of colors
# available in pixmap
function SaveBitmapCl( name, pixMap::Array{UInt8,2}, colorTable::Array{UInt32,1} )
# Get image dimensions
dim = size( pixMap )
if length( dim ) != 2
return 0
end
# Normalize pixel values to get a more accurate reading of the image
minVal = minimum(pixMap)
maxVal = maximum(pixMap)
pixMap = round.(UInt8, 255 * (pixMap .- minVal) ./ (maxVal - minVal))
# 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
# SaveBitmap originally a function of the mzML imzML library in julia,
# now has an adjustment for NaN values in case they exist to mantain data integrity
function GetMzSliceJl(imzML, mass, tolerance)
# Alloc space for slice
width = maximum(imzML[1, :])
height = maximum(imzML[2, :])
image = fill(0.0, width, height)
for i in 1:size(imzML)[2]
index = julia_mzML_imzML.FindMass(imzML[3, i], mass, tolerance)
if index != 0
image[imzML[1, i], imzML[2, i]] = imzML[4, i][index]
end
end
# Adjustment for NaN values with 0
replace!(image, NaN => 0.0)
return image
end
# == Search functions == # == Search functions ==
# Functions that recieve a list to update, and the current direction both as string for # Functions that recieve a list to update, and the current direction both as string for
# searching in the directory the position the list is going # searching in the directory the position the list is going
@ -129,10 +48,12 @@ function loadImgPlot(interfaceImg::String)
layout=PlotlyBase.Layout( layout=PlotlyBase.Layout(
xaxis=PlotlyBase.attr( xaxis=PlotlyBase.attr(
visible=false, visible=false,
scaleanchor="y" scaleanchor="y",
range=[0, width]
), ),
yaxis=PlotlyBase.attr( yaxis=PlotlyBase.attr(
visible=false visible=false,
range=[-height, 0]
), ),
margin=attr(l=0,r=0,t=0,b=0,pad=0) margin=attr(l=0,r=0,t=0,b=0,pad=0)
) )
@ -167,6 +88,7 @@ function loadImgPlot(interfaceImg::String)
plotlayout=layout plotlayout=layout
return plotdata, plotlayout, width, height return plotdata, plotlayout, width, height
end end
# loadImgPlot recieves the local directory of the image as a string, the local directory o the overlay image # loadImgPlot recieves the local directory of the image as a string, the local directory o the overlay image
# and the transparency its required to have. Returns the layout and data for the heatmap plotly plot # and the transparency its required to have. Returns the layout and data for the heatmap plotly plot
# this function loads the image into a plot # this function loads the image into a plot
@ -286,6 +208,7 @@ function loadContourPlot(interfaceImg::String)
plotlayout=layout plotlayout=layout
return plotdata, plotlayout return plotdata, plotlayout
end end
# loadSurfacePlot recieves the local directory of the image as a string, # loadSurfacePlot recieves the local directory of the image as a string,
# returns the layout and data for the surface plotly plot # returns the layout and data for the surface plotly plot
# this function loads the image and applies a gaussian filter # this function loads the image and applies a gaussian filter
@ -363,6 +286,7 @@ function loadSurfacePlot(interfaceImg::String)
plotlayout=layout3D plotlayout=layout3D
return plotdata, plotlayout return plotdata, plotlayout
end end
# This function recieves the x and y coords currently selected, and the dimentions of # This function recieves the x and y coords currently selected, and the dimentions of
# the image to create two traces that will display in a cross section # the image to create two traces that will display in a cross section
function crossLinesPlot(x, y, maxwidth, maxheight) function crossLinesPlot(x, y, maxwidth, maxheight)
@ -404,77 +328,117 @@ function log_tick_formatter(values::Vector{Float64})
end end
# Median filter: an adaptation of the R medianfilter, which averages the matrix
# with the close pixels just from the sides to reduce noise.
# This one in particular is a midpoint fiter from a 3x3 neighbour area
function medianFilterjl(pixMap)
height, width = size(pixMap)
padded_pixMap = padarray(pixMap, 1)
target = zeros(eltype(pixMap), height, width)
for j in 1:width
for i in 1:height
neighbors = []
for dj in j:j+2
for di in i:i+2
push!(neighbors, padded_pixMap[di, dj])
end
end
target[i, j] = median(neighbors)
end
end
return target
end
function padarray(A, padsize)
h, w = size(A)
padded = zeros(eltype(A), h + 2*padsize, w + 2*padsize)
padded[padsize+1:end-padsize, padsize+1:end-padsize] .= A
padded[1:padsize, padsize+1:end-padsize] .= A[1:padsize, :]
padded[end-padsize+1:end, padsize+1:end-padsize] .= A[end-padsize+1:end, :]
padded[padsize+1:end-padsize, 1:padsize] .= A[:, 1:padsize]
padded[padsize+1:end-padsize, end-padsize+1:end] .= A[:, end-padsize+1:end]
padded[1:padsize, 1:padsize] .= A[1, 1]
padded[1:padsize, end-padsize+1:end] .= A[1, end]
padded[end-padsize+1:end, 1:padsize] .= A[end, 1]
padded[end-padsize+1:end, end-padsize+1:end] .= A[end, end]
return padded
end
# meanSpectrumPlot recieves the local directory of the image as a string, # meanSpectrumPlot recieves the local directory of the image as a string,
# returns the layout and data for the surface plotly plot # returns the layout and data for the surface plotly plot
# this function loads the spectra data and makes a mean to display # this function loads the spectra data and makes a mean to display
# its values in the spectrum plot # its values in the spectrum plot
function meanSpectrumPlot(mzmlRoute::String) function meanSpectrumPlot(data::MSIData)
spectraMz=LoadMzml(mzmlRoute) layout = PlotlyBase.Layout(
xSpectraMz=Float64[] title="Average Spectrum Plot",
ySpectraMz=Float64[] hovermode="closest",
layout=PlotlyBase.Layout(
title="Mean spectrum plot",
xaxis=PlotlyBase.attr( xaxis=PlotlyBase.attr(
title="<i>m/z</i>", title="<i>m/z</i>",
showgrid=true
),
yaxis=PlotlyBase.attr(
title="Average Intensity",
showgrid=true, showgrid=true,
tickformat = ".3g" tickformat=".3g"
),
margin=attr(l=0, r=0, t=120, b=0, pad=0)
)
# Use the new, efficient function from the backend
xSpectraMz, ySpectraMz = get_average_spectrum(data)
if isempty(xSpectraMz)
@warn "Average spectrum is empty."
trace = PlotlyBase.stem(x=Float64[], y=Float64[])
else
trace = PlotlyBase.stem(x=xSpectraMz, y=ySpectraMz, marker=attr(size=1, color="blue", opacity=0.5), name="Average", hoverinfo="x",hovertemplate="<b>m/z</b>: %{x:.4f}<extra></extra>")
end
plotdata = [trace]
plotlayout = layout
return plotdata, plotlayout, xSpectraMz, ySpectraMz
end
function xySpectrumPlot(data::MSIData, xCoord::Int, yCoord::Int, imgWidth::Int, imgHeight::Int)
local mz::AbstractVector, intensity::AbstractVector
local plot_title::String
is_imaging = data.source isa ImzMLSource
if is_imaging
# For imaging data, use (X, Y) coordinates
x = clamp(xCoord, 1, imgWidth)
y = clamp(yCoord, 1, imgHeight)
mz, intensity = GetSpectrum(data, x, y)
plot_title = "Spectrum at ($x, $y)"
else
# For non-imaging data, treat xCoord as the spectrum index
index = clamp(xCoord, 1, length(data.spectra_metadata))
mz, intensity = GetSpectrum(data, index)
plot_title = "Spectrum #$index"
end
layout = PlotlyBase.Layout(
title=plot_title,
hovermode="closest",
xaxis=PlotlyBase.attr(
title="<i>m/z</i>",
showgrid=true
), ),
yaxis=PlotlyBase.attr( yaxis=PlotlyBase.attr(
title="Intensity", title="Intensity",
showgrid=true, showgrid=true,
tickformat = ".3g" tickformat=".3g"
), ),
autosize=false, margin=attr(l=0, r=0, t=120, b=0, pad=0)
margin=attr(l=0,r=0,t=120,b=0,pad=0)
) )
try
xSpectraMz=mean(spectraMz[1,:]) # Downsample for plotting performance
ySpectraMz=mean(spectraMz[2,:]) mz_down, int_down = MSI_src.downsample_spectrum(mz, intensity)
catch e
xSpectraMz=spectraMz[1,1] trace = PlotlyBase.stem(x=mz_down, y=int_down, marker=attr(size=1, color="blue", opacity=0.5), name="Spectrum", hoverinfo="x", hovertemplate="<b>m/z</b>: %{x:.4f}<extra></extra>")
ySpectraMz=spectraMz[2,1]
plotdata = [trace]
plotlayout = layout
# Return the full data for other uses, and the plot data
return plotdata, plotlayout, mz, intensity
end
function sumSpectrumPlot(data::MSIData)
layout = PlotlyBase.Layout(
title="Total Spectrum Plot",
hovermode="closest",
xaxis=PlotlyBase.attr(
title="<i>m/z</i>",
showgrid=true
),
yaxis=PlotlyBase.attr(
title="Total Intensity",
showgrid=true,
tickformat=".3g"
),
margin=attr(l=0, r=0, t=120, b=0, pad=0)
)
# Use the get_total_spectrum function from the backend
xSpectraMz, ySpectraMz = get_total_spectrum(data)
if isempty(xSpectraMz)
@warn "Total spectrum is empty."
trace = PlotlyBase.stem(x=Float64[], y=Float64[])
else
trace = PlotlyBase.stem(x=xSpectraMz, y=ySpectraMz, marker=attr(size=1, color="blue", opacity=0.5), name="Total", hoverinfo="x",hovertemplate="<b>m/z</b>: %{x:.4f}<extra></extra>")
end end
trace=PlotlyBase.stem(x=xSpectraMz, y=ySpectraMz,marker=attr(size=1, color="blue", opacity=0.1))
plotdata=[trace] # We add the data from spectra to the plot plotdata = [trace]
plotlayout=layout plotlayout = layout
return plotdata, plotlayout, xSpectraMz, ySpectraMz return plotdata, plotlayout, xSpectraMz, ySpectraMz
end end

97
src/MSIData.jl
View File

@ -6,7 +6,7 @@ including caching and iteration logic, for handling large mzML and imzML dataset
efficiently. efficiently.
""" """
using Base64, Libz # For reading binary data using Base64, Libz, Serialization # For reading binary data
# Abstract type for different data sources (e.g., mzML, imzML) # Abstract type for different data sources (e.g., mzML, imzML)
# This allows dispatching to the correct binary reading logic. # This allows dispatching to the correct binary reading logic.
@ -179,17 +179,10 @@ function GetSpectrum(data::MSIData, x::Int, y::Int)
return GetSpectrum(data, index) # Call the existing method return GetSpectrum(data, index) # Call the existing method
end end
""" using Serialization
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. function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000)
The m/z range is determined dynamically by finding the global min/max m/z values println("Calculating total spectrum (2-pass method for imzML)...")
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 # 1. First Pass: Find the global m/z range
println(" Pass 1: Finding global m/z range...") println(" Pass 1: Finding global m/z range...")
@ -222,7 +215,6 @@ function get_total_spectrum(msi_data::MSIData; num_bins::Int=20000)
return # continue equivalent return # continue equivalent
end end
for i in eachindex(mz) 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) bin_index = clamp(round(Int, (mz[i] - global_min_mz) / bin_step) + 1, 1, num_bins)
intensity_sum[bin_index] += intensity[i] intensity_sum[bin_index] += intensity[i]
end end
@ -232,6 +224,83 @@ function get_total_spectrum(msi_data::MSIData; num_bins::Int=20000)
return (collect(mz_bins), intensity_sum) return (collect(mz_bins), intensity_sum)
end end
function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000)
println("Calculating total spectrum (single-pass optimization for mzML)...")
num_spectra = length(msi_data.spectra_metadata)
if num_spectra == 0
return (Float64[], Float64[])
end
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...")
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
end
flush(temp_io)
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)]")
# --- Pass 2: Read from temp file and bin intensities ---
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)
bin_step = step(mz_bins)
while !eof(temp_io)
mz, intensity = Serialization.deserialize(temp_io)::Tuple{AbstractVector, AbstractVector}
if isempty(mz)
continue
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]
end
end
println("Total spectrum calculation complete.")
return (collect(mz_bins), intensity_sum)
finally
close(temp_io)
rm(temp_path, force=true)
end
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.
This function dispatches to a specialized implementation based on the file type
(.imzML or .mzML) for optimal performance.
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=2000)
if msi_data.source isa ImzMLSource
return get_total_spectrum_imzml(msi_data, num_bins=num_bins)
else # MzMLSource
return get_total_spectrum_mzml(msi_data, num_bins=num_bins)
end
end
""" """
get_average_spectrum(msi_data::MSIData; num_bins::Int=20000) -> Tuple{Vector{Float64}, Vector{Float64}} get_average_spectrum(msi_data::MSIData; num_bins::Int=20000) -> Tuple{Vector{Float64}, Vector{Float64}}
@ -241,7 +310,7 @@ This is effectively the total ion chromatogram (TIC) divided by the number of sp
Returns a tuple containing two vectors: the binned m/z axis and the averaged intensities. 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) function get_average_spectrum(msi_data::MSIData; num_bins::Int=2000)
# This function uses the exact same logic as get_total_spectrum... # This function uses the exact same logic as get_total_spectrum...
mz_bins, intensity_sum = get_total_spectrum(msi_data, num_bins=num_bins) mz_bins, intensity_sum = get_total_spectrum(msi_data, num_bins=num_bins)
@ -365,5 +434,3 @@ function _iterate_spectra_fast(f::Function, data::MSIData)
# Dispatch to the correct implementation based on the source type # Dispatch to the correct implementation based on the source type
_iterate_spectra_fast_impl(f, data, data.source) _iterate_spectra_fast_impl(f, data, data.source)
end end

304
src/MzmlConverter.jl
View File

@ -26,26 +26,19 @@ struct BinaryMetadata
end end
""" """
GetMzmlScanTime(fileName::String) GetMzmlScanTime_linebyline(fileName::String)
Parses a .mzML file to extract the scan start time for each spectrum. (Internal) Slow, line-by-line parser for scan times. Used as a fallback if
the .mzML file index is missing or corrupt.
# Arguments
* `fileName`: Path to the .mzML file.
# Returns
- A `Matrix{Int64}` where each row is `[spectrum_index, time_in_milliseconds]`.
""" """
function GetMzmlScanTime(fileName::String) function GetMzmlScanTime_linebyline(fileName::String)
times = Tuple{Int64, Int64}[] times = Tuple{Int64, Int64}[]
# Pre-allocate based on an estimate to reduce re-allocations
try try
file_size = filesize(fileName) file_size = filesize(fileName)
estimated_spectra = max(1000, file_size ÷ 10000) # Heuristic estimated_spectra = max(1000, file_size ÷ 10000)
sizehint!(times, estimated_spectra) sizehint!(times, estimated_spectra)
catch catch
# Ignore if filesize fails, sizehint! is just an optimization
end end
open(fileName, "r") do stream open(fileName, "r") do stream
@ -60,27 +53,23 @@ function GetMzmlScanTime(fileName::String)
current_index = 0 current_index = 0
current_time = nothing current_time = nothing
# Process index from the start tag itself
idx_match = match(r"index=\"(\d+)\"", line) idx_match = match(r"index=\"(\d+)\"", line)
if idx_match !== nothing if idx_match !== nothing
current_index = parse(Int, idx_match.captures[1]) + 1 current_index = parse(Int, idx_match.captures[1]) + 1
end end
end end
elseif state == :in_spectrum elseif state == :in_spectrum
# Look for time cvParam
if occursin("<cvParam", line) if occursin("<cvParam", line)
if occursin("MS:1000016", line) || occursin("MS:1000015", line) if occursin("MS:1000016", line) || occursin("MS:1000015", line)
time_match = match(r"value=\"([\d\.]+)\"", line) time_match = match(r"value=\"([\d\.]+)\"", line)
if time_match !== nothing if time_match !== nothing
time_val = parse(Float64, time_match.captures[1]) time_val = parse(Float64, time_match.captures[1])
# MS:1000016 is minutes, MS:1000015 is seconds
unit_scale = occursin("MS:1000016", line) ? 60000 : 1000 unit_scale = occursin("MS:1000016", line) ? 60000 : 1000
current_time = round(Int64, time_val * unit_scale) current_time = round(Int64, time_val * unit_scale)
end end
end end
end end
# Check for end of spectrum
if occursin("</spectrum>", line) if occursin("</spectrum>", line)
state = :outside_spectrum state = :outside_spectrum
if current_index > 0 && current_time !== nothing if current_index > 0 && current_time !== nothing
@ -91,9 +80,100 @@ function GetMzmlScanTime(fileName::String)
end end
end end
# Sort by spectrum index to ensure correct order before normalization
sort!(times, by = x -> x[1]) sort!(times, by = x -> x[1])
if !isempty(times)
first_time = times[1][2]
for i in eachindex(times)
times[i] = (times[i][1], times[i][2] - first_time)
end
end
if isempty(times)
return Matrix{Int64}(undef, 0, 2)
end
return permutedims(hcat(collect.(times)...))
end
"""
GetMzmlScanTime(fileName::String)
Parses a .mzML file to extract the scan start time for each spectrum.
Uses the indexed part of the .mzML file for fast access, falling back to
a slower line-by-line parse if the index is not present.
# Arguments
* `fileName`: Path to the .mzML file.
# Returns
- A `Matrix{Int64}` where each row is `[spectrum_index, time_in_milliseconds]`.
"""
function GetMzmlScanTime(fileName::String)
times = Tuple{Int64, Int64}[]
try
open(fileName, "r") do stream
# 1. Find and parse the spectrum index offsets
seekend(stream)
end_chunk_size = min(filesize(stream), 8192)
seek(stream, filesize(stream) - end_chunk_size)
footer = read(stream, String)
index_offset_match = match(r"<indexListOffset>(\d+)</indexListOffset>", footer)
if index_offset_match === nothing
@warn "No <indexListOffset> found. Falling back to slow line-by-line parsing for scan times. This may be memory intensive."
return GetMzmlScanTime_linebyline(fileName)
end
index_offset = parse(Int64, index_offset_match.captures[1])
seek(stream, index_offset)
# The find_tag function is defined in ParserHelpers.jl
if find_tag(stream, r"<index\s+name=\"spectrum\"") === nothing
@warn "Could not find spectrum index. Falling back to slow line-by-line parsing."
return GetMzmlScanTime_linebyline(fileName)
end
# The parse_offset_list function is defined in mzML.jl
spectrum_offsets = parse_offset_list(stream)
# 2. Iterate through offsets and parse time for each spectrum
for (idx, offset) in enumerate(spectrum_offsets)
seek(stream, offset)
current_time = nothing
# Read a limited number of lines to find the time
for _ in 1:100
if eof(stream) break end
line = readline(stream)
if occursin("MS:1000016", line) || occursin("MS:1000015", line) # scan start time
time_match = match(r"value=\"([\d\.]+)\"", line)
if time_match !== nothing
time_val = parse(Float64, time_match.captures[1])
unit_scale = occursin("MS:1000016", line) ? 60000 : 1000 # minutes vs seconds
current_time = round(Int64, time_val * unit_scale)
break
end
end
if occursin("<binary>", line) || occursin("</spectrum>", line)
break
end
end
if current_time !== nothing
push!(times, (idx, current_time))
end
end
end
catch e
@error "Failed to parse scan times with indexed method. Falling back to line-by-line." exception=(e, catch_backtrace())
return GetMzmlScanTime_linebyline(fileName)
end
# The index is already sorted by spectrum index, so no need to sort `times`.
# Normalize times relative to the first scan # Normalize times relative to the first scan
if !isempty(times) if !isempty(times)
first_time = times[1][2] first_time = times[1][2]
@ -102,7 +182,6 @@ function GetMzmlScanTime(fileName::String)
end end
end end
# Convert vector of tuples to a matrix
if isempty(times) if isempty(times)
return Matrix{Int64}(undef, 0, 2) return Matrix{Int64}(undef, 0, 2)
end end
@ -329,131 +408,113 @@ function RenderPixel(pixel_info, scans, msi_data::MSIData, scan_time_deltas, pix
return (mz_array, new_intensity) return (mz_array, new_intensity)
end end
function ConvertMzmlToImzml(source_file::String, timing_matrix::Matrix{Int64}, scans::Matrix{Int64}) function ConvertMzmlToImzml(source_file::String, target_ibd_file::String, timing_matrix::Matrix{Int64}, scans::Matrix{Int64})
if size(timing_matrix, 1) == 0 if size(timing_matrix, 1) == 0
return ProcessedPixel[], (0, 0) # Create an empty .ibd file if there's nothing to process
open(target_ibd_file, "w") do ibd_stream
write(ibd_stream, zeros(UInt8, 16)) # UUID placeholder
end
return BinaryMetadata[], Tuple{Int, Int}[], (0, 0)
end end
# Get image dimensions from timing_matrix (assuming columns 1 and 2 are X and Y)
width = maximum(timing_matrix[:, 1]) width = maximum(timing_matrix[:, 1])
height = maximum(timing_matrix[:, 2]) height = maximum(timing_matrix[:, 2])
# Open the mzML file using the new memory-efficient API
msi_data = OpenMSIData(source_file) msi_data = OpenMSIData(source_file)
# Pre-calculate time deltas robustly
scan_time_deltas = zeros(Int64, size(scans, 1)) scan_time_deltas = zeros(Int64, size(scans, 1))
if size(scans, 1) > 1 if size(scans, 1) > 1
for i in 1:(size(scans, 1) - 1) for i in 1:(size(scans, 1) - 1)
delta = scans[i+1, 2] - scans[i, 2] delta = scans[i+1, 2] - scans[i, 2]
scan_time_deltas[i] = max(1, delta) # Ensure delta is at least 1 scan_time_deltas[i] = max(1, delta)
end end
scan_time_deltas[end] = max(1, scan_time_deltas[end-1]) # Use previous delta for last scan, ensure at least 1 scan_time_deltas[end] = max(1, scan_time_deltas[end-1])
end end
pixel_time_deltas = zeros(Int64, size(timing_matrix, 1)) pixel_time_deltas = zeros(Int64, size(timing_matrix, 1))
if size(timing_matrix, 1) > 1 if size(timing_matrix, 1) > 1
for i in 1:(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] delta = timing_matrix[i+1, 3] - timing_matrix[i, 3]
pixel_time_deltas[i] = max(1, delta) # Ensure delta is at least 1 pixel_time_deltas[i] = max(1, delta)
end end
pixel_time_deltas[end] = max(1, pixel_time_deltas[end-1]) # Use previous delta for last pixel, ensure at least 1 pixel_time_deltas[end] = max(1, pixel_time_deltas[end-1])
end end
processed_pixels = ProcessedPixel[] binary_meta_vec = BinaryMetadata[]
sizehint!(processed_pixels, size(timing_matrix, 1)) sizehint!(binary_meta_vec, size(timing_matrix, 1))
coords_vec = Tuple{Int, Int}[]
sizehint!(coords_vec, size(timing_matrix, 1))
empty_pixel_count = 0 empty_pixel_count = 0
# DEBUG: Scan coverage analysis open(target_ibd_file, "w") do ibd_stream
println("DEBUG: Scan coverage analysis:") write(ibd_stream, zeros(UInt8, 16)) # UUID placeholder
covered_pixels_count = 0
for i in 1:size(timing_matrix, 1) for i in 1:size(timing_matrix, 1)
first_scan = timing_matrix[i, 4] pixel_info = timing_matrix[i, :]
last_scan = timing_matrix[i, 5] x, y = pixel_info[1], pixel_info[2]
if first_scan <= last_scan && first_scan >= 1 && last_scan <= size(scans, 1) push!(coords_vec, (x, y))
covered_pixels_count += 1
first_scan = pixel_info[4]
last_scan = pixel_info[5]
if first_scan > last_scan || first_scan < 1 || last_scan > size(scans, 1)
empty_pixel_count += 1
# For empty pixels, offsets point to the current end of file, with zero length
current_pos = position(ibd_stream)
push!(binary_meta_vec, BinaryMetadata(current_pos, 0, current_pos, 0))
continue
end
mz, intensity = RenderPixel(pixel_info, scans, msi_data, scan_time_deltas, pixel_time_deltas)
# Write m/z array
mz_offset = position(ibd_stream)
for val in mz
write(ibd_stream, htol(Float64(val)))
end
mz_length = position(ibd_stream) - mz_offset
# Write intensity array
int_offset = position(ibd_stream)
for val in intensity
write(ibd_stream, htol(Float32(val)))
end
int_length = position(ibd_stream) - int_offset
push!(binary_meta_vec, BinaryMetadata(mz_offset, mz_length, int_offset, int_length))
end end
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) @info "Found and processed $empty_pixel_count empty pixels out of $(size(timing_matrix, 1)) total."
pixel_info = timing_matrix[i, :] return binary_meta_vec, coords_vec, (width, height)
x = pixel_info[1]
y = pixel_info[2]
first_scan = pixel_info[4]
last_scan = pixel_info[5]
if first_scan > last_scan || first_scan < 1 || last_scan > size(scans, 1)
empty_pixel_count += 1
push!(processed_pixels, ProcessedPixel((x, y), Float64[], Float32[]))
continue
end
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 end
function ExportImzml(target_file::String, pixels::Vector{ProcessedPixel}, dims::Tuple{Int, Int}) 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") ibd_file = replace(target_file, r"\.imzML$"i => ".ibd")
# Ensure we have valid pixels if isempty(binary_meta)
if isempty(pixels) @warn "No binary metadata to export; creating empty imzML file."
@error "No pixels to export" # Still create a valid, empty imzML file
return false
end end
binary_meta = BinaryMetadata[]
sizehint!(binary_meta, length(pixels))
try try
# Step 1: Write .ibd file with proper binary format # The .ibd file is now written by ConvertMzmlToImzml.
open(ibd_file, "w") do ibd_stream # This function is only responsible for the .imzML XML metadata file.
# 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 open(target_file, "w") do imzml_stream
# Use indexedmzML wrapper # XML Header
write(imzml_stream, """<?xml version="1.0" encoding="ISO-8859-1"?> 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"> <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])"> <mzML version="1.1" id="$(splitext(basename(target_file))[1])">
""") """)
# CV List # CV List, File Description, etc. (static parts)
write(imzml_stream, """ <cvList count="3"> 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="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="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"/> <cv id="IMS" fullName="Imaging MS Ontology" version="0.9.1" URI="http://www.maldi-msi.org/download/imzml/imagingMS.obo"/>
</cvList> </cvList>
""") """)
# File Description
write(imzml_stream, """ <fileDescription> write(imzml_stream, """ <fileDescription>
<fileContent> <fileContent>
<cvParam cvRef="MS" accession="MS:1000579" name="MS1 spectrum"/> <cvParam cvRef="MS" accession="MS:1000579" name="MS1 spectrum"/>
@ -462,8 +523,6 @@ function ExportImzml(target_file::String, pixels::Vector{ProcessedPixel}, dims::
</fileContent> </fileContent>
</fileDescription> </fileDescription>
""") """)
# Referenceable Param Groups
write(imzml_stream, """ <referenceableParamGroupList count="2"> write(imzml_stream, """ <referenceableParamGroupList count="2">
<referenceableParamGroup id="mzArray"> <referenceableParamGroup id="mzArray">
<cvParam cvRef="MS" accession="MS:1000576" name="no compression"/> <cvParam cvRef="MS" accession="MS:1000576" name="no compression"/>
@ -479,24 +538,18 @@ function ExportImzml(target_file::String, pixels::Vector{ProcessedPixel}, dims::
</referenceableParamGroup> </referenceableParamGroup>
</referenceableParamGroupList> </referenceableParamGroupList>
""") """)
# Sample List
write(imzml_stream, """ <sampleList count="1"> write(imzml_stream, """ <sampleList count="1">
<sample id="sample1" name="ImagingSample"> <sample id="sample1" name="ImagingSample">
<cvParam cvRef="MS" accession="MS:1000001" name="sample number" value="1"/> <cvParam cvRef="MS" accession="MS:1000001" name="sample number" value="1"/>
</sample> </sample>
</sampleList> </sampleList>
""") """)
# Software List - OPEN SOURCE
write(imzml_stream, """ <softwareList count="1"> write(imzml_stream, """ <softwareList count="1">
<software id="MSIConverter" version="1.0"> <software id="MSIConverter" version="1.0">
<cvParam cvRef="MS" accession="MS:1000799" name="custom unreleased software tool" value="mzML to imzML converter"/> <cvParam cvRef="MS" accession="MS:1000799" name="custom unreleased software tool" value="mzML to imzML converter"/>
</software> </software>
</softwareList> </softwareList>
""") """)
# Scan Settings IMAGING
write(imzml_stream, """ <scanSettingsList count="1"> write(imzml_stream, """ <scanSettingsList count="1">
<scanSettings id="scanSettings1"> <scanSettings id="scanSettings1">
<cvParam cvRef="IMS" accession="IMS:1000042" name="max count of pixel x" value="$(dims[1])"/> <cvParam cvRef="IMS" accession="IMS:1000042" name="max count of pixel x" value="$(dims[1])"/>
@ -506,8 +559,6 @@ function ExportImzml(target_file::String, pixels::Vector{ProcessedPixel}, dims::
</scanSettings> </scanSettings>
</scanSettingsList> </scanSettingsList>
""") """)
# Instrument Configuration MALDI/TOF
write(imzml_stream, """ <instrumentConfigurationList count="1"> write(imzml_stream, """ <instrumentConfigurationList count="1">
<instrumentConfiguration id="instrument1"> <instrumentConfiguration id="instrument1">
<componentList count="3"> <componentList count="3">
@ -525,8 +576,6 @@ function ExportImzml(target_file::String, pixels::Vector{ProcessedPixel}, dims::
</instrumentConfiguration> </instrumentConfiguration>
</instrumentConfigurationList> </instrumentConfigurationList>
""") """)
# Data Processing
write(imzml_stream, """ <dataProcessingList count="1"> write(imzml_stream, """ <dataProcessingList count="1">
<dataProcessing id="conversionProcessing"> <dataProcessing id="conversionProcessing">
<processingMethod order="1" softwareRef="MSIConverter"> <processingMethod order="1" softwareRef="MSIConverter">
@ -539,18 +588,21 @@ function ExportImzml(target_file::String, pixels::Vector{ProcessedPixel}, dims::
# Run and Spectrum List # Run and Spectrum List
spectrum_offsets = UInt64[] spectrum_offsets = UInt64[]
write(imzml_stream, """ <run defaultInstrumentConfigurationRef="instrument1" id="run1" sampleRef="sample1"> write(imzml_stream, """ <run defaultInstrumentConfigurationRef="instrument1" id="run1" sampleRef="sample1">
<spectrumList count="$(length(pixels))" defaultDataProcessingRef="conversionProcessing"> <spectrumList count="$(length(binary_meta))" defaultDataProcessingRef="conversionProcessing">
""") """)
# Write each spectrum # Write each spectrum's metadata
for (i, pixel) in enumerate(pixels) for (i, meta) in enumerate(binary_meta)
meta = binary_meta[i] x, y = coords[i]
x, y = pixel.coords
spectrum_start = position(imzml_stream) spectrum_start = position(imzml_stream)
push!(spectrum_offsets, spectrum_start) push!(spectrum_offsets, spectrum_start)
write(imzml_stream, """ <spectrum id="Scan=$(i)" defaultArrayLength="$(length(pixel.mz))" index="$(i-1)"> # Calculate number of points from byte length
mz_points = meta.mz_length ÷ sizeof(Float64)
int_points = meta.int_length ÷ sizeof(Float32)
write(imzml_stream, """ <spectrum id="Scan=$(i)" defaultArrayLength="$(mz_points)" index="$(i-1)">
<cvParam cvRef="MS" accession="MS:1000511" name="ms level" value="1"/> <cvParam cvRef="MS" accession="MS:1000511" name="ms level" value="1"/>
<cvParam cvRef="MS" accession="MS:1000128" name="profile spectrum"/> <cvParam cvRef="MS" accession="MS:1000128" name="profile spectrum"/>
<scanList count="1"> <scanList count="1">
@ -563,14 +615,14 @@ function ExportImzml(target_file::String, pixels::Vector{ProcessedPixel}, dims::
<binaryDataArray encodedLength="0"> <binaryDataArray encodedLength="0">
<referenceableParamGroupRef ref="mzArray"/> <referenceableParamGroupRef ref="mzArray"/>
<cvParam cvRef="IMS" accession="IMS:1000102" name="external offset" value="$(meta.mz_offset)"/> <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:1000103" name="external array length" value="$(mz_points)"/>
<cvParam cvRef="IMS" accession="IMS:1000104" name="external encoded length" value="$(meta.mz_length)"/> <cvParam cvRef="IMS" accession="IMS:1000104" name="external encoded length" value="$(meta.mz_length)"/>
<binary/> <binary/>
</binaryDataArray> </binaryDataArray>
<binaryDataArray encodedLength="0"> <binaryDataArray encodedLength="0">
<referenceableParamGroupRef ref="intensityArray"/> <referenceableParamGroupRef ref="intensityArray"/>
<cvParam cvRef="IMS" accession="IMS:1000102" name="external offset" value="$(meta.int_offset)"/> <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:1000103" name="external array length" value="$(int_points)"/>
<cvParam cvRef="IMS" accession="IMS:1000104" name="external encoded length" value="$(meta.int_length)"/> <cvParam cvRef="IMS" accession="IMS:1000104" name="external encoded length" value="$(meta.int_length)"/>
<binary/> <binary/>
</binaryDataArray> </binaryDataArray>
@ -603,10 +655,10 @@ function ExportImzml(target_file::String, pixels::Vector{ProcessedPixel}, dims::
return true return true
catch e catch e
@error "Failed to export imzML files" exception=(e, catch_backtrace()) @error "Failed to export imzML metadata file" exception=(e, catch_backtrace())
# Clean up partial files # Clean up partial .imzML file
isfile(ibd_file) && rm(ibd_file, force=true)
isfile(target_file) && rm(target_file, force=true) isfile(target_file) && rm(target_file, force=true)
# Do not delete the .ibd file as it might be useful for debugging
return false return false
end end
end end
@ -630,17 +682,15 @@ function ImportMzmlFile(source_file::String, sync_file::String, target_file::Str
println("Step 2: Matching acquisition times...") println("Step 2: Matching acquisition times...")
timing_matrix = MatchAcquireTime(sync_file, scans; img_width=img_width, img_height=img_height) timing_matrix = MatchAcquireTime(sync_file, scans; img_width=img_width, img_height=img_height)
println("Step 3: Converting spectra...") println("Step 3: Converting spectra and writing .ibd file...")
processed_pixels, (width, height) = ConvertMzmlToImzml(source_file, timing_matrix, scans) ibd_file = replace(target_file, r"\.imzML$"i => ".ibd")
binary_meta, coords, (width, height) = ConvertMzmlToImzml(source_file, ibd_file, timing_matrix, scans)
# Flip image vertically to match R script output # Flip image vertically to match R script output
for i in eachindex(processed_pixels) flipped_coords = [(x, height - y + 1) for (x, y) in coords]
x, y = processed_pixels[i].coords
processed_pixels[i] = ProcessedPixel((x, height - y + 1), processed_pixels[i].mz, processed_pixels[i].intensity)
end
println("Step 4: Exporting to .imzML/.ibd format...") println("Step 4: Exporting .imzML metadata file...")
success = ExportImzml(target_file, processed_pixels, (width, height)) success = ExportImzml(target_file, binary_meta, flipped_coords, (width, height))
if success if success
println("Conversion successful: $target_file") println("Conversion successful: $target_file")

214
src/imzML.jl
View File

@ -1,4 +1,4 @@
using Images, Statistics, CairoMakie, DataFrames, Printf, ColorSchemes using Images, Statistics, CairoMakie, DataFrames, Printf, ColorSchemes, StatsBase
# --- Extracted from imzML.jl --- # --- Extracted from imzML.jl ---
@ -347,6 +347,38 @@ function load_slices(folder, masses, tolerance)
end end
"""
get_mz_slice(data::MSIData, mass::Real, tolerance::Real)
Extracts an image slice for a given m/z value without plotting.
This is a performant function that iterates through spectra once.
# Returns
- A `Matrix{Float64}` representing the intensity slice.
"""
function get_mz_slice(data::MSIData, mass::Real, tolerance::Real)
width, height = data.image_dims
slice_matrix = zeros(Float64, height, width) # Note: height, width for (y,x) indexing
_iterate_spectra_fast(data) do spec_idx, mz_array, intensity_array
meta = data.spectra_metadata[spec_idx]
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
# The original function had a NaN replacement, which is good practice to keep.
replace!(slice_matrix, NaN => 0.0)
return slice_matrix
end
""" """
plot_slice(msi_data::MSIData, mass::Float64, tolerance::Float64, output_dir::String; stage_name="slice", bins=256) plot_slice(msi_data::MSIData, mass::Float64, tolerance::Float64, output_dir::String; stage_name="slice", bins=256)
@ -486,6 +518,30 @@ function set_pixel_depth(img, bounds, depth)
return result return result
end end
"""
TrIQ(pixMap, depth, prob=0.98)
Applies TrIQ (Treshold Intensity Quantization) normalization to an image.
This function first computes dynamic intensity range bounds by identifying outliers
based on a cumulative probability, then sets the pixel depth (quantizes intensities)
within these bounds.
# Arguments
- `pixMap`: The input image matrix (e.g., a slice from `GetMzSliceJl`).
- `depth`: The number of grey levels (bins) to quantize the intensities into.
- `prob`: The target cumulative probability (e.g., 0.98) used to determine outlier thresholds.
# Returns
- A new image matrix with intensities quantized to the specified depth within the TrIQ bounds.
"""
function TrIQ(pixMap, depth, prob=0.98)
# Compute new dynamic range
bounds = get_outlier_thres(pixMap, prob)
# Set intensity dynamic range
return set_pixel_depth(pixMap, bounds, depth)
end
""" """
norm_slices_hist(slices, bins; prob=0.98) norm_slices_hist(slices, bins; prob=0.98)
@ -515,11 +571,73 @@ function norm_slices_hist(slices, bins; prob=0.98)
return (norm_img=norm_img, bounds=mass_bounds) # bounds is now a vector, one per mass return (norm_img=norm_img, bounds=mass_bounds) # bounds is now a vector, one per mass
end end
"""
quantize_intensity(slice::AbstractMatrix{<:Real}, levels::Integer=256)
Linearly scales the intensity values in a slice to a specified number of levels.
The output is an array of `UInt8` values.
This is a modernized version of `IntQuantCl`.
"""
function quantize_intensity(slice::AbstractMatrix{<:Real}, levels::Integer=256)
max_val = maximum(slice)
if max_val <= 0
return zeros(UInt8, size(slice))
end
# Scale relative to the absolute maximum value to preserve the zero point.
# The original logic used 'colorLevel' which was the max value (e.g., 255).
scale = (levels - 1) / max_val
# round is equivalent to floor(x+0.5) for positive numbers.
# clamp is used for robustness against floating point inaccuracies.
image = round.(UInt8, clamp.(slice .* scale, 0, levels - 1))
return image
end
function median_filter(img) function median_filter(img)
# 3x3 median filter implementation # 3x3 median filter implementation
return mapwindow(median, img, (3, 3)) return mapwindow(median, img, (3, 3))
end end
"""
downsample_spectrum(mz, intensity, n_points=2000)
Reduces the number of points in a spectrum for faster plotting, while preserving peaks.
It divides the m/z range into `n_points` bins and keeps only the most intense point from each bin.
"""
function downsample_spectrum(mz::AbstractVector, intensity::AbstractVector, n_points::Integer=2000)
if isempty(mz) || length(mz) <= n_points
return mz, intensity
end
mz_min, mz_max = extrema(mz)
bin_width = (mz_max - mz_min) / n_points
# We use a vector of tuples to store the max intensity and its corresponding mz for each bin
# (max_intensity, mz_value)
bins = fill((0.0, 0.0), n_points)
for i in eachindex(mz)
# Determine the bin for the current point
# Bin indices are 1-based
bin_index = min(n_points, floor(Int, (mz[i] - mz_min) / bin_width) + 1)
# If the current point's intensity is higher than what's in the bin, replace it
if intensity[i] > bins[bin_index][1]
bins[bin_index] = (intensity[i], mz[i])
end
end
# Filter out empty bins and separate the mz and intensity values
final_mz = [b[2] for b in bins if b[1] > 0.0]
final_intensity = [b[1] for b in bins if b[1] > 0.0]
return final_mz, final_intensity
end
# ============================================================================ # ============================================================================
# #
# #
@ -678,3 +796,97 @@ function plot_slices(slices, names, masses, output_dir; stage_name, bins=256, dp
end end
return fig return fig
end end
"""
save_bitmap(name::String, pixMap::Matrix{UInt8}, colorTable::Vector{UInt32})
Saves an 8-bit indexed image as a BMP file.
This is a modernized version of `SaveBitmapCl`.
"""
function save_bitmap(name::String, pixMap::Matrix{UInt8}, colorTable::Vector{UInt32})
# Get image dimensions
height, width = size(pixMap)
# Normalize pixel values to stretch contrast, as in the original SaveBitmapCl
minVal, maxVal = extrema(pixMap)
if maxVal > minVal
pixMap = round.(UInt8, 255 * (pixMap .- minVal) ./ (maxVal - minVal))
end
# Compute row padding (each row must be a multiple of 4 bytes)
padding = (4 - (width % 4)) % 4
# BMP color table must have 256 entries for 8-bit images
fullColorTable = Vector{UInt32}(undef, 256)
if length(colorTable) <= 256
fullColorTable[1:length(colorTable)] .= colorTable
fullColorTable[length(colorTable)+1:end] .= 0
else
fullColorTable .= colorTable[1:256]
end
# Compute file dimensions
offset = 14 + 40 + (256 * 4) # 14(file) + 40(info) + 1024(palette) = 1078
imgBytes = height * (width + padding)
fileSize = offset + imgBytes
open(name, "w") do stream
# === File Header (14 bytes) ===
write(stream, UInt16(0x4D42)) # "BM"
write(stream, UInt32(fileSize))
write(stream, UInt16(0)) # Reserved
write(stream, UInt16(0)) # Reserved
write(stream, UInt32(offset))
# === Info Header (40 bytes) ===
write(stream, UInt32(40)) # Info header size
write(stream, Int32(width))
write(stream, Int32(height)) # Positive for bottom-up storage
write(stream, UInt16(1)) # Number of color planes
write(stream, UInt16(8)) # Bits per pixel
write(stream, UInt32(0)) # Compression (BI_RGB)
write(stream, UInt32(imgBytes))
write(stream, Int32(2835)) # Pels per meter X (~72 DPI)
write(stream, Int32(2835)) # Pels per meter Y (~72 DPI)
write(stream, UInt32(256)) # Colors in color table
write(stream, UInt32(0)) # Important colors (0 = all)
# === Color Table ===
write(stream, fullColorTable)
# === Image Pixels (written bottom-up) ===
row_buffer = Vector{UInt8}(undef, width + padding)
row_buffer[width+1:end] .= 0 # pre-fill padding bytes
for i in height:-1:1
row_data = @view pixMap[i, :]
row_buffer[1:width] = row_data
write(stream, row_buffer)
end
end
end
# ********************************************************************
# Viridis color palette (256 colors)
# ********************************************************************
"""
generate_palette(colorscheme, n_colors=256)
Generates a BMP-compatible UInt32 color palette from a ColorScheme.
"""
function generate_palette(colorscheme, n_colors=256)
palette = Vector{UInt32}(undef, n_colors)
colors = get(colorscheme, range(0, 1, length=n_colors))
for i in 1:n_colors
c = colors[i]
r = round(UInt8, c.r * 255)
g = round(UInt8, c.g * 255)
b = round(UInt8, c.b * 255)
# BMP color table format is 0x00RRGGBB, written little-endian becomes BB GG RR 00
palette[i] = (UInt32(r) << 16) | (UInt32(g) << 8) | UInt32(b)
end
return palette
end
const ViridisPalette = generate_palette(ColorSchemes.viridis)

83
test/readme.md
View File

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

24
test/run_tests.jl
View File

@ -32,7 +32,8 @@ using MSI_src
# A regular, non-imaging mzML file. # A regular, non-imaging mzML file.
# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/mzML/T9_A1.mzML" # const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/mzML/T9_A1.mzML"
# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML" # const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML"
const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/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 SPECTRUM_TO_PLOT = 1 # Which spectrum to plot from the file const SPECTRUM_TO_PLOT = 1 # Which spectrum to plot from the file
# --- Test Case 2: .mzML + Sync File for Conversion --- # --- Test Case 2: .mzML + Sync File for Conversion ---
@ -49,15 +50,18 @@ const CONVERSION_TARGET_IMZML = "test/results/converted_mzml.imzML"
# --- Test Case 3: Standard .imzML file --- # --- Test Case 3: Standard .imzML file ---
# An existing imzML file (can be the one generated from Case 2). # 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/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_paper_spray/Imaging_paper_spray.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. # The m/z value to use for creating an image slice.
const MZ_VALUE_FOR_SLICE = 309.06 # BF # const MZ_VALUE_FOR_SLICE = 309.06 # BF
# const MZ_VALUE_FOR_SLICE = 896.0 # HR2MSI # const MZ_VALUE_FOR_SLICE = 896.0 # HR2MSI
# const MZ_VALUE_FOR_SLICE = 76.03 # I PS # const MZ_VALUE_FOR_SLICE = 76.03 # I PS
const MZ_VALUE_FOR_SLICE = 473 # ROYA
# const MZ_TOLERANCE = 0.1 # const MZ_TOLERANCE = 0.1
const MZ_TOLERANCE = 1 # const MZ_TOLERANCE = 1
const MZ_TOLERANCE = 0.5
# Coordinates to plot a specific spectrum from imzML # Coordinates to plot a specific spectrum from imzML
const COORDS_TO_PLOT = (50, 50) # Example coordinates (X, Y) const COORDS_TO_PLOT = (50, 50) # Example coordinates (X, Y)
@ -65,8 +69,8 @@ const COORDS_TO_PLOT = (50, 50) # Example coordinates (X, Y)
# --- Output Directory --- # --- Output Directory ---
const RESULTS_DIR = "test/results" const RESULTS_DIR = "test/results"
test1 = true test1 = false
test2 = true test2 = false
test3 = true test3 = true
# =================================================================== # ===================================================================
@ -216,8 +220,7 @@ function run_test()
end end
# --- Test Case 3: Process an existing .imzML file --- # --- Test Case 3: Process an existing .imzML file ---
println(" println("" * "="^20 * " Test Case 3: Processing .imzML " * "="^20)
" * "="^20 * " Test Case 3: Processing .imzML " * "="^20)
if test3 == true if test3 == true
# Run new, stronger validation for imzML # Run new, stronger validation for imzML
validate_msi_data(TEST_IMZML_FILE) validate_msi_data(TEST_IMZML_FILE)
@ -226,6 +229,7 @@ function run_test()
if isfile(TEST_IMZML_FILE) if isfile(TEST_IMZML_FILE)
# Add spectrum plotting for imzML to match Test Case 1 # Add spectrum plotting for imzML to match Test Case 1
try try
"""
# Get the msi data from the imzml # Get the msi data from the imzml
println("Plotting a sample spectrum from $TEST_IMZML_FILE...") println("Plotting a sample spectrum from $TEST_IMZML_FILE...")
msi_data = OpenMSIData(TEST_IMZML_FILE) msi_data = OpenMSIData(TEST_IMZML_FILE)
@ -271,6 +275,8 @@ function run_test()
output_path = joinpath(RESULTS_DIR, "test_imzml_average_spectrum.png") output_path = joinpath(RESULTS_DIR, "test_imzml_average_spectrum.png")
save(output_path, fig) save(output_path, fig)
println("SUCCESS: Total spectrum plot saved to $output_path") println("SUCCESS: Total spectrum plot saved to $output_path")
"""
println("No spectrums tested on this try.")
catch e catch e
println("ERROR during spectrum plotting in Test Case 3: $e") println("ERROR during spectrum plotting in Test Case 3: $e")
end end
@ -279,7 +285,7 @@ function run_test()
try try
println("Testing plot_slice function on $TEST_IMZML_FILE...") println("Testing plot_slice function on $TEST_IMZML_FILE...")
msi_data = OpenMSIData(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") @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 # The success message is now inside plot_slice
catch e catch e
println("ERROR during plot_slice test in Test Case 3: $e") println("ERROR during plot_slice test in Test Case 3: $e")