diff --git a/Manifest.toml b/Manifest.toml
index c6e5873..d9db884 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -2,7 +2,7 @@
julia_version = "1.11.7"
manifest_format = "2.0"
-project_hash = "7843f141174fe17a3abc4f25c3819d3195459b5a"
+project_hash = "4955dfdbb7a233ddefa69e3cce2833f416ffe1bc"
[[deps.ANSIColoredPrinters]]
git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
diff --git a/Project.toml b/Project.toml
index d21383f..94a4fe2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -35,11 +35,13 @@ Loess = "4345ca2d-374a-55d4-8d30-97f9976e7612"
Mmap = "a63ad114-7e13-5084-954f-fe012c677804"
NativeFileDialog = "e1fe445b-aa65-4df4-81c1-2041507f0fd4"
NaturalSort = "c020b1a1-e9b0-503a-9c33-f039bfc54a85"
+Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
SavitzkyGolay = "c4bf5708-b6a6-4fbe-bcd0-6850ed671584"
Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
StipplePlotly = "ec984513-233d-481d-95b0-a3b58b97af2b"
diff --git a/src/MSIData.jl b/src/MSIData.jl
index 6b7c8d7..7789a89 100644
--- a/src/MSIData.jl
+++ b/src/MSIData.jl
@@ -6,7 +6,7 @@ including caching and iteration logic, for handling large mzML and imzML dataset
efficiently.
"""
-using Base64, Libz, Serialization, Printf, DataFrames, Base.Threads # For reading binary data
+using Base64, Libz, Serialization, Printf, DataFrames, Base.Threads, StatsBase # For reading binary data
const FILE_HANDLE_LOCK = ReentrantLock()
@@ -28,9 +28,16 @@ function ThreadSafeFileHandle(path::String, mode::String="r")
return ThreadSafeFileHandle(path, handle, ReentrantLock())
end
-function Base.seek(tsfh::ThreadSafeFileHandle, pos::Integer)
+"""
+ read_at!(tsfh::ThreadSafeFileHandle, a::AbstractArray, pos::Integer)
+
+Atomically seeks to a position and reads data into an array. This is the thread-safe
+way to read from a specific offset in the file.
+"""
+function read_at!(tsfh::ThreadSafeFileHandle, a::AbstractArray, pos::Integer)
lock(tsfh.lock) do
seek(tsfh.handle, pos)
+ read!(tsfh.handle, a)
end
end
@@ -40,12 +47,6 @@ function Base.read(tsfh::ThreadSafeFileHandle, n::Integer)
end
end
-function Base.read!(tsfh::ThreadSafeFileHandle, a::AbstractArray)
- lock(tsfh.lock) do
- return read!(tsfh.handle, a)
- end
-end
-
function Base.close(tsfh::ThreadSafeFileHandle)
close(tsfh.handle)
end
@@ -66,6 +67,20 @@ function Base.filesize(tsfh::ThreadSafeFileHandle)
end
end
+# Overload Base.seek for ThreadSafeFileHandle
+function Base.seek(tsfh::ThreadSafeFileHandle, pos::Integer)
+ lock(tsfh.lock) do
+ seek(tsfh.handle, pos)
+ end
+end
+
+# Overload Base.read! for ThreadSafeFileHandle
+function Base.read!(tsfh::ThreadSafeFileHandle, a::AbstractArray)
+ lock(tsfh.lock) do
+ read!(tsfh.handle, a)
+ end
+end
+
# Abstract type for different data sources (e.g., mzML, imzML)
# This allows dispatching to the correct binary reading logic.
abstract type MSDataSource end
@@ -131,6 +146,50 @@ An enum representing the type of a spectrum's data.
"""
@enum SpectrumMode CENTROID=1 PROFILE=2 UNKNOWN=3
+"""
+ InstrumentMetadata
+
+Contains metadata about the instrument configuration and acquisition parameters.
+This information is parsed from the file header and is crucial for guiding
+preprocessing steps.
+
+# Fields
+- `resolution`: Instrument resolution at a specific m/z, if available.
+- `acquisition_mode`: The overall data acquisition mode (`:profile`, `:centroid`, or `:mixed`).
+- `mz_axis_type`: The nature of the m/z axis (`:regular`, `:irregular`, `:mixed`).
+- `calibration_status`: The calibration state of the data (`:internal`, `:external`, `:uncalibrated`).
+- `instrument_model`: The vendor and model of the instrument.
+- `mass_accuracy_ppm`: Manufacturer-specified mass accuracy in parts-per-million.
+- `laser_settings`: A dictionary of MSI-specific laser parameters.
+- `polarity`: The polarity of the acquisition (`:positive`, `:negative`, `:unknown`).
+"""
+struct InstrumentMetadata
+ resolution::Union{Float64, Nothing}
+ acquisition_mode::Symbol
+ mz_axis_type::Symbol
+ calibration_status::Symbol
+ instrument_model::String
+ mass_accuracy_ppm::Union{Float64, Nothing}
+ laser_settings::Union{Dict, Nothing}
+ polarity::Symbol
+ vendor_preprocessing_steps::Union{Vector{String}, Nothing} # New field for provenance
+end
+
+# Default constructor for initializing with unknown values
+function InstrumentMetadata()
+ return InstrumentMetadata(
+ nothing, # resolution
+ :unknown, # acquisition_mode
+ :unknown, # mz_axis_type
+ :uncalibrated, # calibration_status
+ "Not Available",# instrument_model
+ nothing, # mass_accuracy_ppm
+ nothing, # laser_settings
+ :unknown, # polarity
+ nothing # vendor_preprocessing
+ )
+end
+
"""
SpectrumMetadata
@@ -150,6 +209,7 @@ struct SpectrumMetadata
# For mzML
id::String
+ type::Symbol # NEW: e.g., :sample, :blank, :qc
mode::SpectrumMode
# Common binary data info
@@ -180,6 +240,7 @@ efficient repeated access to spectra.
mutable struct MSIData
source::MSDataSource
spectra_metadata::Vector{SpectrumMetadata}
+ instrument_metadata::Union{InstrumentMetadata, Nothing}
image_dims::Tuple{Int, Int}
coordinate_map::Union{Matrix{Int}, Nothing}
@@ -198,13 +259,14 @@ mutable struct MSIData
spectrum_stats_df::Union{DataFrame, Nothing}
bloom_filters::Union{Vector{<:BloomFilter}, Nothing}
analytics_ready::AtomicFlag
+ preprocessing_hints::Union{Dict{Symbol, Any}, Nothing} # For auto-determined parameters
- function MSIData(source, metadata, dims, coordinate_map, cache_size)
- obj = new(source, metadata, dims, coordinate_map,
+ function MSIData(source, metadata, instrument_meta, dims, coordinate_map, cache_size)
+ obj = new(source, metadata, instrument_meta, dims, coordinate_map,
Dict(), [], cache_size, ReentrantLock(),
SimpleBufferPool(),
Threads.Atomic{Float64}(Inf), Threads.Atomic{Float64}(-Inf),
- nothing, nothing, AtomicFlag())
+ nothing, nothing, AtomicFlag(), nothing)
# Ensure file handles are closed when the object is garbage collected
finalizer(obj) do o
@@ -407,16 +469,17 @@ function read_binary_vector(data::MSIData, io::IO, asset::SpectrumAsset)
# Calculate number of elements
n_elements = length(decoded_bytes) ÷ sizeof(asset.format)
- out_array = Vector{asset.format}(undef, n_elements)
-
- # Use unsafe_wrap to avoid copy - this is the critical optimization
- temp_array = unsafe_wrap(Vector{asset.format}, pointer(decoded_bytes), n_elements)
-
- # Convert byte order in-place
- @inbounds for i in 1:n_elements
- out_array[i] = ltoh(temp_array[i])
+
+ if n_elements * sizeof(asset.format) != length(decoded_bytes)
+ throw(FileFormatError("Size of decoded byte array is not a multiple of the element size."))
end
+ # Reinterpret the byte array as an array of the target type. This does not copy.
+ reinterpreted_array = reinterpret(asset.format, decoded_bytes)
+
+ # Allocate the final output array and convert byte order while copying.
+ out_array = [ltoh(x) for x in reinterpreted_array]
+
return out_array
end
@@ -442,28 +505,18 @@ and converting it from little-endian to the host's native byte order.
- A tuple `(mz, intensity)` containing the two requested data arrays.
"""
function read_spectrum_from_disk(source::ImzMLSource, meta::SpectrumMetadata)
- # For imzML, the binary data is raw, not base64 encoded.
+ # For imzML, the binary data is raw, not base64 encoded.
# The `encoded_length` field in this case holds the number of points.
mz = Array{source.mz_format}(undef, meta.mz_asset.encoded_length)
intensity = Array{source.intensity_format}(undef, meta.int_asset.encoded_length)
- # FIX: Read arrays in the order they appear in the file for efficiency,
- # but ensure we seek to the correct offset for both.
- if meta.mz_asset.offset < meta.int_asset.offset
- seek(source.ibd_handle, meta.mz_asset.offset)
- read!(source.ibd_handle, mz)
- seek(source.ibd_handle, meta.int_asset.offset)
- read!(source.ibd_handle, intensity)
- else
- seek(source.ibd_handle, meta.int_asset.offset)
- read!(source.ibd_handle, intensity)
- seek(source.ibd_handle, meta.mz_asset.offset)
- read!(source.ibd_handle, mz)
- end
+ # Use the new atomic read_at! method for thread-safety
+ read_at!(source.ibd_handle, mz, meta.mz_asset.offset)
+ read_at!(source.ibd_handle, intensity, meta.int_asset.offset)
# imzML data is little-endian. Convert to host byte order.
mz .= ltoh.(mz)
- intensity .= ltoh.(intensity) # FIX: Added missing conversion for intensity
+ intensity .= ltoh.(intensity)
validate_spectrum_data(mz, intensity, meta.id)
@@ -733,6 +786,25 @@ function precompute_analytics(msi_data::MSIData)
g_min_mz = minimum(thread_local_min_mz)
g_max_mz = maximum(thread_local_max_mz)
+ # Only update if we found valid ranges
+ if isfinite(g_min_mz) && isfinite(g_max_mz) && g_min_mz < g_max_mz
+ Threads.atomic_xchg!(msi_data.global_min_mz, g_min_mz)
+ Threads.atomic_xchg!(msi_data.global_max_mz, g_max_mz)
+ else
+ # Fallback: calculate from first spectrum
+ try
+ if num_spectra > 0
+ mz, _ = GetSpectrum(msi_data, 1)
+ if !isempty(mz)
+ Threads.atomic_xchg!(msi_data.global_min_mz, minimum(mz))
+ Threads.atomic_xchg!(msi_data.global_max_mz, maximum(mz))
+ end
+ end
+ catch e
+ @warn "Could not determine global m/z range from first spectrum: $e"
+ end
+ end
+
# Create results
computed_stats_df = DataFrame(
SpectrumID = 1:num_spectra,
@@ -745,10 +817,6 @@ function precompute_analytics(msi_data::MSIData)
Mode = modes
)
- # Set atomic values
- Threads.atomic_xchg!(msi_data.global_min_mz, g_min_mz)
- Threads.atomic_xchg!(msi_data.global_max_mz, g_max_mz)
-
# Set other data
msi_data.spectrum_stats_df = computed_stats_df
msi_data.bloom_filters = bloom_filters # Will be created on-demand
@@ -1480,30 +1548,3 @@ function get_bloom_filter(data::MSIData, index::Int)::BloomFilter{Int}
return data.bloom_filters[index]
end
end
-
-"""
- monitor_memory_usage(phase::String)
-
-Utility function to print current memory usage (allocated bytes and RSS)
-at different phases of execution.
-
-# Arguments
-- `phase`: A string describing the current phase (e.g., "Initialization", "Processing chunk 1").
-"""
-function monitor_memory_usage(phase::String)
- # Use Julia's internal memory tracking
- bytes_allocated = Base.gc_bytes()
- mb_allocated = bytes_allocated / 1024^2
-
- # Get process memory from system
- if Sys.islinux()
- proc_status = read("/proc/self/status", String)
- vm_rss_match = match(r"VmRSS:\s*(\d+)\s*kB", proc_status)
- if vm_rss_match !== nothing
- vm_rss_mb = parse(Int, vm_rss_match[1]) / 1024
- @printf "[%s] Memory: %.2f MB allocated, %.2f MB RSS\n" phase mb_allocated vm_rss_mb
- end
- else
- @printf "[%s] Memory: %.2f MB allocated\n" phase mb_allocated
- end
-end
diff --git a/src/MSI_src.jl b/src/MSI_src.jl
index d878bd2..bc3c7ee 100644
--- a/src/MSI_src.jl
+++ b/src/MSI_src.jl
@@ -7,7 +7,6 @@ export OpenMSIData,
GetSpectrum,
IterateSpectra,
ImportMzmlFile,
- plot_slices,
plot_slice,
get_total_spectrum,
get_average_spectrum,
@@ -16,31 +15,34 @@ export OpenMSIData,
generate_colorbar_image,
process_image_pipeline,
load_and_prepare_mask,
- set_global_mz_range!
+ set_global_mz_range!,
+ get_global_mz_range,
+ MSIData,
+ _iterate_spectra_fast,
+ validate_spectrum
-# Export the public Preprocessing API
-export FeatureMatrix,
- run_preprocessing_pipeline,
- qc_is_empty,
- qc_is_regular,
- transform_intensity,
- smooth_spectrum,
- snip_baseline,
- tic_normalize,
- pqn_normalize,
- median_normalize,
+# Export the public preprocessing & precalculations API
+export run_preprocessing_analysis,
+ main_precalculation,
+ FeatureMatrix,
+ Calibration,
+ Smoothing,
+ BaselineCorrection,
+ Normalization,
+ PeakPicking,
+ PeakBinningParams,
+ get_masked_spectrum_indices,
detect_peaks_profile,
- detect_peaks_wavelet,
detect_peaks_centroid,
- align_peaks_lowess,
- find_calibration_peaks,
+ smooth_spectrum,
+ apply_baseline_correction,
+ apply_normalization,
bin_peaks,
- plot_stage_spectrum,
- calculate_ppm_error,
- calculate_resolution_fwhm,
- analyze_mass_accuracy,
- generate_qc_report,
- get_common_calibration_standards
+ PeakSelection,
+ PeakAlignment,
+ find_calibration_peaks,
+ align_peaks_lowess,
+ MutableSpectrum
# Include all source files directly into the main module
include("BloomFilters.jl")
@@ -52,6 +54,9 @@ include("imzML.jl")
include("MzmlConverter.jl")
include("Preprocessing.jl")
include("ImageProcessing.jl")
+include("Precalculations.jl")
+
+using Setfield # For immutable struct updates
# --- Main Entry Point --- #
@@ -63,14 +68,28 @@ Opens a .mzML or .imzML file and prepares it for data access.
This is the main entry point for the new data access API.
"""
-function OpenMSIData(filepath::String; cache_size=300)
+function OpenMSIData(filepath::String; cache_size=300, spectrum_type_map::Union{Dict{Int, Symbol}, Nothing}=nothing)
+ local msi_data
if endswith(lowercase(filepath), ".mzml")
- return load_mzml_lazy(filepath, cache_size=cache_size)
+ msi_data = load_mzml_lazy(filepath, cache_size=cache_size)
elseif endswith(lowercase(filepath), ".imzml")
- return load_imzml_lazy(filepath, cache_size=cache_size)
+ msi_data = load_imzml_lazy(filepath, cache_size=cache_size)
else
error("Unsupported file type: $filepath. Please provide a .mzML or .imzML file.")
end
+
+ # Apply spectrum type map if provided
+ if spectrum_type_map !== nothing
+ for (idx, type_symbol) in spectrum_type_map
+ if 1 <= idx <= length(msi_data.spectra_metadata)
+ msi_data.spectra_metadata[idx] = @set msi_data.spectra_metadata[idx].type = type_symbol
+ else
+ @warn "Spectrum index $idx out of bounds for spectrum_type_map. Skipping."
+ end
+ end
+ end
+
+ return msi_data
end
end # module MSI_src
diff --git a/src/MzmlConverter.jl b/src/MzmlConverter.jl
index f908e5b..64db42a 100644
--- a/src/MzmlConverter.jl
+++ b/src/MzmlConverter.jl
@@ -601,7 +601,7 @@ function ConvertMzmlToImzml(source_file::String, target_ibd_file::String, timing
close(msi_data)
end
- return binary_meta_vec, coords_vec, (width, height), pixel_modes, ibd_uuid
+ return binary_meta_vec, coords_vec, (width, height), pixel_modes, ibd_uuid, msi_data.instrument_metadata
end
"""
@@ -621,7 +621,7 @@ experiment and data format.
# Returns
- `true` on success, `false` on failure.
"""
-function ExportImzml(target_file::String, binary_meta::Vector{BinaryMetadata}, coords::Vector{Tuple{Int, Int}}, dims::Tuple{Int, Int}, modes::Vector{SpectrumMode}, ibd_uuid::UUID)
+function ExportImzml(target_file::String, binary_meta::Vector{BinaryMetadata}, coords::Vector{Tuple{Int, Int}}, dims::Tuple{Int, Int}, modes::Vector{SpectrumMode}, ibd_uuid::UUID, instrument_meta::InstrumentMetadata)
ibd_file = replace(target_file, r"\.imzML$"i => ".ibd")
if isempty(binary_meta)
@@ -652,6 +652,15 @@ function ExportImzml(target_file::String, binary_meta::Vector{BinaryMetadata}, c
+$(
+ if instrument_meta.polarity == :positive
+ """ """
+ elseif instrument_meta.polarity == :negative
+ """ """
+ else
+ ""
+ end
+)
""")
@@ -693,12 +702,35 @@ function ExportImzml(target_file::String, binary_meta::Vector{BinaryMetadata}, c
""")
write(imzml_stream, """
+$(
+ if instrument_meta.instrument_model != "Not Available"
+ """ """
+ else
+ ""
+ end
+)
+$(
+ if instrument_meta.resolution !== nothing
+ """ """
+ else
+ ""
+ end
+)
+$(
+ if instrument_meta.acquisition_mode == :centroid
+ """ """
+ elseif instrument_meta.acquisition_mode == :profile
+ """ """
+ else
+ ""
+ end
+)
@@ -708,11 +740,23 @@ function ExportImzml(target_file::String, binary_meta::Vector{BinaryMetadata}, c
""")
- write(imzml_stream, """
+ write(imzml_stream, """
+$(
+ if instrument_meta.vendor_preprocessing_steps !== nothing
+ join([
+ """
+
+ """
+ for (i, step) in enumerate(instrument_meta.vendor_preprocessing_steps)
+ ], "\n")
+ else
+ ""
+ end
+)
""")
@@ -812,13 +856,13 @@ function ImportMzmlFile(source_file::String, sync_file::String, target_file::Str
println("Step 3: Converting spectra and writing .ibd file...")
ibd_file = replace(target_file, r"\.imzML$"i => ".ibd")
- binary_meta, coords, (width, height), pixel_modes, ibd_uuid = ConvertMzmlToImzml(source_file, ibd_file, timing_matrix, scans)
+ binary_meta, coords, (width, height), pixel_modes, ibd_uuid, instrument_meta = ConvertMzmlToImzml(source_file, ibd_file, timing_matrix, scans)
# Flip image vertically to match R script output
flipped_coords = [(x, height - y + 1) for (x, y) in coords]
println("Step 4: Exporting .imzML metadata file...")
- success = ExportImzml(target_file, binary_meta, flipped_coords, (width, height), pixel_modes, ibd_uuid)
+ success = ExportImzml(target_file, binary_meta, flipped_coords, (width, height), pixel_modes, ibd_uuid, instrument_meta)
if success
println("Conversion successful: $target_file")
diff --git a/src/Precalculations.jl b/src/Precalculations.jl
new file mode 100644
index 0000000..849a036
--- /dev/null
+++ b/src/Precalculations.jl
@@ -0,0 +1,1537 @@
+using StatsBase # For mean, std, median, quantile, mad
+
+# =============================================================================
+# 7) Spatial & Advanced Processing (Stubs & New Functions)
+# =============================================================================
+
+"""
+ find_ppm_error_by_region(msi_data, region_masks, reference_peaks)::Dict
+
+Calculates PPM error statistics for different spatial regions.
+`region_masks` is a Dict mapping region names (e.g., :tumor) to BitMatrix masks.
+"""
+function find_ppm_error_by_region(msi_data::MSIData, region_masks::Dict, reference_peaks::Dict)
+ regional_reports = Dict{Symbol, NamedTuple}()
+
+ width, height = msi_data.image_dims
+
+ for (region_name, mask) in region_masks
+ mask_height, mask_width = size(mask)
+ if mask_width != width || mask_height != height
+ @warn "Mask dimensions ($(mask_width)x$(mask_height)) for region '$region_name' do not match image dimensions ($(width)x$(height)). Skipping."
+ continue
+ end
+
+ indices = [
+ i for i in 1:length(msi_data.spectra_metadata)
+ if msi_data.spectra_metadata[i].x > 0 &&
+ msi_data.spectra_metadata[i].y > 0 &&
+ mask[msi_data.spectra_metadata[i].y, msi_data.spectra_metadata[i].x]
+ ]
+
+ if isempty(indices) continue end
+
+ # Call analyze_mass_accuracy with the specific indices for the region
+ regional_reports[region_name] = analyze_mass_accuracy(msi_data, reference_peaks; spectrum_indices=indices)
+ end
+ return regional_reports
+end
+
+function analyze_mass_accuracy(
+ msi_data::MSIData,
+ reference_peaks::Dict{Float64, String}; # m/z => name
+ spectrum_indices::AbstractVector{Int},
+ peak_detection_snr_threshold::Float64 = 2.0,
+ ppm_tolerance_for_matching::Float64 = 50.0
+ )::NamedTuple
+ println("Analyzing mass accuracy for $(length(spectrum_indices)) spectra...")
+
+ all_ppm_errors = Float64[]
+ total_matched_peaks = 0
+ total_spectra_processed = 0
+
+ _iterate_spectra_fast(msi_data, spectrum_indices) do idx, mz, intensity
+ total_spectra_processed += 1
+ if !validate_spectrum(mz, intensity)
+ @warn "Spectrum $idx is invalid, skipping mass accuracy analysis for it."
+ return
+ end
+
+ detected_peaks = detect_peaks_profile(mz, intensity; snr_threshold=peak_detection_snr_threshold)
+
+ for ref_mz in keys(reference_peaks)
+ # Find the closest detected peak to this reference m/z within tolerance
+ min_ppm_error = Inf
+ best_matched_peak_mz = nothing
+
+ for p in detected_peaks
+ ppm_error = calculate_ppm_error(p.mz, ref_mz)
+ if ppm_error <= ppm_tolerance_for_matching && ppm_error < min_ppm_error
+ min_ppm_error = ppm_error
+ best_matched_peak_mz = p.mz
+ end
+ end
+
+ if best_matched_peak_mz !== nothing
+ push!(all_ppm_errors, min_ppm_error)
+ total_matched_peaks += 1
+ end
+ end
+ end
+
+ if isempty(all_ppm_errors)
+ @warn "No reference peaks matched in any of the $(total_spectra_processed) processed spectra."
+ return (
+ mean_ppm_error = NaN,
+ median_ppm_error = NaN,
+ std_ppm_error = NaN,
+ min_ppm_error = NaN,
+ max_ppm_error = NaN,
+ total_matched_peaks = 0,
+ total_spectra_analyzed = total_spectra_processed
+ )
+ end
+
+ # Calculate summary statistics
+ mean_err = mean(all_ppm_errors)
+ median_err = median(all_ppm_errors)
+ std_err = std(all_ppm_errors)
+ min_err = minimum(all_ppm_errors)
+ max_err = maximum(all_ppm_errors)
+
+ println("Mass accuracy analysis complete for $(total_spectra_processed) spectra.")
+ return (
+ mean_ppm_error = mean_err,
+ median_ppm_error = median_err,
+ std_ppm_error = std_err,
+ min_ppm_error = min_err,
+ max_ppm_error = max_err,
+ total_matched_peaks = total_matched_peaks,
+ total_spectra_analyzed = total_spectra_processed
+ )
+end
+
+# =============================================================================
+# 8) Advanced Peak Quality & Adaptive Parameters
+# =============================================================================
+
+"""
+ calculate_adaptive_bin_tolerance(ppm_error_distribution) -> Float64
+
+Calculates an appropriate binning tolerance based on observed mass accuracy.
+"""
+function calculate_adaptive_bin_tolerance(ppm_error_distribution::Vector{Float64})
+ if isempty(ppm_error_distribution)
+ return 20.0 # Default if no data
+ end
+
+ # Filter out NaN and infinite values
+ valid_errors = filter(x -> isfinite(x), ppm_error_distribution)
+
+ if isempty(valid_errors)
+ return 20.0 # Default if no valid data
+ end
+
+ # A robust strategy: mean + 3 * std deviation to capture ~99.7% of peaks
+ mean_err = mean(valid_errors)
+ std_err = std(valid_errors)
+
+ # Ensure we don't get NaN
+ if !isfinite(mean_err) || !isfinite(std_err)
+ return 20.0
+ end
+
+ tolerance = mean_err + 3 * std_err
+
+ # Cap at reasonable value and ensure finite
+ return min(max(tolerance, 10.0), 100.0)
+end
+
+"""
+ calculate_preprocessing_hints(data::MSIData; sample_size::Int=100)::Dict{Symbol, Any}
+
+Analyzes a sample of spectra to determine optimal default parameters for
+preprocessing steps, returning a dictionary of hints.
+"""
+function calculate_preprocessing_hints(data::MSIData; sample_indices::AbstractVector{Int})::Dict{Symbol, Any}
+ println("Calculating preprocessing hints from a sample of $(length(sample_indices)) spectra...")
+
+ num_spectra = length(data.spectra_metadata)
+ if num_spectra == 0 || isempty(sample_indices)
+ @warn "No spectra in dataset/sample to calculate hints from."
+ return Dict{Symbol, Any}(
+ :estimated_noise => 1.0,
+ :suggested_snr => 3.0,
+ :suggested_smoothing_window => 9
+ )
+ end
+
+ # Initialize hints with defaults immediately
+ hints = Dict{Symbol, Any}(
+ :estimated_noise => 1.0,
+ :suggested_snr => 3.0,
+ :suggested_smoothing_window => 9
+ )
+
+ all_noise_levels = Float64[]
+
+ for idx in sample_indices
+ try
+ mz, intensity = GetSpectrum(data, idx)
+ if !isempty(intensity)
+ noise = mad(intensity, normalize=true)
+ if isfinite(noise) && noise > 0
+ push!(all_noise_levels, noise)
+ end
+ end
+ catch e
+ @warn "Could not access spectrum #$idx to calculate hints: $e. Skipping."
+ end
+ end
+
+ if !isempty(all_noise_levels)
+ estimated_noise = mean(all_noise_levels)
+ hints[:estimated_noise] = estimated_noise
+ println(" - Estimated Noise Level: $(round(estimated_noise, digits=4))")
+ println(" - Suggested SNR Threshold: 3.0")
+ else
+ @warn "Could not estimate noise from sample. Using default hints."
+ # Defaults already set in `hints` initialization
+ end
+
+ # Suggest smoothing window based on resolution (if available)
+ if data.instrument_metadata !== nothing && data.instrument_metadata.resolution !== nothing
+ res = data.instrument_metadata.resolution
+ # Update suggested smoothing window if resolution is known
+ if res > 40000
+ hints[:suggested_smoothing_window] = 5
+ elseif res > 10000
+ hints[:suggested_smoothing_window] = 7
+ else
+ hints[:suggested_smoothing_window] = 9
+ end
+ end
+ println(" - Suggested Smoothing Window: $(hints[:suggested_smoothing_window])")
+
+ println("Preprocessing hints calculated.")
+ return hints
+end
+
+"""
+ analyze_instrument_characteristics(msi_data::MSIData)::Dict
+
+Analyzes instrument metadata and data characteristics to determine acquisition properties.
+"""
+function analyze_instrument_characteristics(msi_data::MSIData; sample_indices::AbstractVector{Int})::Dict
+ results = Dict{Symbol, Any}()
+
+ # Use instrument metadata if available
+ if msi_data.instrument_metadata !== nothing
+ inst = msi_data.instrument_metadata
+ results[:resolution] = inst.resolution
+ results[:mass_accuracy_ppm] = inst.mass_accuracy_ppm
+ results[:instrument_model] = inst.instrument_model
+ results[:polarity] = inst.polarity
+ results[:calibration_status] = inst.calibration_status
+ results[:vendor_preprocessing] = inst.vendor_preprocessing_steps
+ end
+
+ # Analyze data characteristics from spectra
+ println(" Analyzing data characteristics from spectra...")
+
+ spectrum_modes = Set{SpectrumMode}()
+ mz_step_sizes = Float64[]
+ intensity_ranges = Tuple{Float64, Float64}[] # (min, max) per spectrum
+
+ if isempty(sample_indices)
+ @warn "No indices to sample for instrument characteristics analysis."
+ results[:acquisition_mode] = :unknown
+ results[:mz_axis_type] = :unknown
+ results[:dynamic_range] = 3.0
+ return results
+ end
+
+ _iterate_spectra_fast(msi_data, sample_indices) do idx, mz, intensity
+ # Record spectrum mode
+ push!(spectrum_modes, msi_data.spectra_metadata[idx].mode)
+
+ # Calculate m/z step statistics (for profile data)
+ if length(mz) > 1 && msi_data.spectra_metadata[idx].mode == PROFILE
+ steps = diff(mz)
+ if !isempty(steps)
+ push!(mz_step_sizes, mean(steps))
+ end
+ end
+
+ # Record intensity range
+ if !isempty(intensity)
+ push!(intensity_ranges, (minimum(intensity), maximum(intensity)))
+ end
+ end
+
+ # Determine acquisition mode
+ if length(spectrum_modes) == 1
+ results[:acquisition_mode] = first(spectrum_modes) == CENTROID ? :centroid : :profile
+ else
+ results[:acquisition_mode] = :mixed
+ end
+
+ # Determine m/z axis regularity
+ if !isempty(mz_step_sizes)
+ avg_step = mean(mz_step_sizes)
+ step_std = std(mz_step_sizes)
+ results[:mz_axis_type] = step_std / avg_step < 0.01 ? :regular : :irregular
+ results[:average_mz_step] = avg_step
+ else
+ results[:mz_axis_type] = :unknown
+ end
+
+ # Enhanced dynamic range calculation
+ if !isempty(intensity_ranges)
+ println("DEBUG: Found $(length(intensity_ranges)) intensity ranges")
+
+ # FIX: Handle cases where min intensities are zero
+ max_intensities = Float64[]
+ min_positive_intensities = Float64[]
+
+ for (min_val, max_val) in intensity_ranges
+ if max_val > 1e-6 # Valid maximum
+ push!(max_intensities, max_val)
+
+ # Find the smallest positive intensity in the spectrum
+ # For now, use a reasonable estimate: 1% of the noise level
+ # In practice, you'd want to sample the actual spectrum
+ if max_val > 0
+ # Estimate minimum detectable signal as ~3x noise level
+ estimated_min_signal = max_val * 1e-4 # Conservative estimate
+ push!(min_positive_intensities, estimated_min_signal)
+ end
+ end
+ end
+
+ println("DEBUG: Valid max intensities: $(length(max_intensities)), estimated min intensities: $(length(min_positive_intensities))")
+
+ if !isempty(max_intensities) && !isempty(min_positive_intensities)
+ avg_max = mean(max_intensities)
+ avg_min = mean(min_positive_intensities)
+ if avg_min > 0
+ dynamic_range = log10(avg_max / avg_min)
+ results[:dynamic_range] = dynamic_range
+ println("DEBUG: Dynamic range calculated: $dynamic_range (avg_max=$avg_max, avg_min=$avg_min)")
+ else
+ results[:dynamic_range] = 0.0
+ end
+ else
+ # Estimate based on typical values
+ results[:dynamic_range] = 3.0 # Typical for MS data
+ end
+ else
+ results[:dynamic_range] = 3.0 # Default estimate
+ end
+
+ println(" - Acquisition mode: $(results[:acquisition_mode])")
+ println(" - m/z axis type: $(results[:mz_axis_type])")
+ println(" - Dynamic range: $(round(get(results, :dynamic_range, 0), digits=2))")
+
+ return results
+end
+
+"""
+ analyze_signal_quality(msi_data::MSIData; sample_size::Int=100)::Dict
+
+Analyzes noise characteristics, signal-to-noise ratios, and overall signal quality.
+"""
+function analyze_signal_quality(msi_data::MSIData; sample_indices::AbstractVector{Int})::Dict
+ results = Dict{Symbol, Any}()
+
+ println(" Analyzing signal quality from $(length(sample_indices)) sample spectra...")
+
+ # Use the existing function for basic noise estimation
+ hints_from_calc = calculate_preprocessing_hints(msi_data, sample_indices=sample_indices)
+
+ # Copy relevant hints to results, or handle cases where they might be missing
+ results[:estimated_noise] = get(hints_from_calc, :estimated_noise, 1.0)
+ results[:suggested_snr] = get(hints_from_calc, :suggested_snr, 3.0)
+ results[:suggested_smoothing_window] = get(hints_from_calc, :suggested_smoothing_window, 9)
+
+ # Enhanced noise analysis
+ noise_levels = Float64[]
+ snr_distribution = Float64[]
+ tic_values = Float64[]
+
+ if isempty(sample_indices)
+ @warn "No indices to sample for signal quality analysis."
+ return results
+ end
+
+ _iterate_spectra_fast(msi_data, sample_indices) do idx, mz, intensity
+ if !isempty(intensity)
+ # Noise estimation using MAD
+ noise = mad(intensity, normalize=true)
+ push!(noise_levels, noise)
+
+ # FIX: More robust SNR calculation
+ valid_intensity = intensity[intensity .> 0] # Remove zeros
+ if !isempty(valid_intensity)
+ # Use robust signal estimate (95th percentile instead of max)
+ signal_estimate = quantile(valid_intensity, 0.95)
+ noise_robust = max(noise, 1e-6) # Avoid division by zero
+
+ if noise_robust > 0 && isfinite(signal_estimate)
+ snr_val = signal_estimate / noise_robust
+ # Cap unrealistic SNR values
+ push!(snr_distribution, min(snr_val, 1e6))
+ end
+ end
+
+ # Total ion count
+ push!(tic_values, sum(intensity))
+ end
+ end
+
+ if !isempty(noise_levels)
+ results[:noise_mean] = mean(noise_levels)
+ results[:noise_std] = std(noise_levels)
+ results[:noise_cv] = results[:noise_std] / results[:noise_mean] # Coefficient of variation
+ end
+
+ if !isempty(snr_distribution)
+ results[:snr_mean] = mean(snr_distribution)
+ results[:snr_median] = median(snr_distribution)
+ results[:snr_95th] = quantile(snr_distribution, 0.95)
+ end
+
+ if !isempty(tic_values)
+ results[:tic_mean] = mean(tic_values)
+ results[:tic_std] = std(tic_values)
+ results[:tic_cv] = results[:tic_std] / results[:tic_mean]
+ end
+
+ println(" - Estimated noise level: $(round(get(results, :noise_mean, 0), digits=4))")
+ println(" - Average SNR: $(round(get(results, :snr_mean, 0), digits=2))")
+ println(" - TIC CV: $(round(get(results, :tic_cv, 0) * 100, digits=1))%")
+
+ return results
+end
+
+"""
+ analyze_mass_accuracy_global(msi_data::MSIData, reference_peaks::Dict, sample_size::Int)::Dict
+
+Analyzes mass accuracy across the dataset using reference peaks.
+"""
+function analyze_mass_accuracy_global(msi_data::MSIData, reference_peaks::Dict;
+ spectrum_indices::AbstractVector{Int})::Dict
+ results = Dict{Symbol, Any}()
+
+ println(" Analyzing mass accuracy using $(length(reference_peaks)) reference peaks on $(length(spectrum_indices)) spectra...")
+
+ if isempty(spectrum_indices)
+ @warn "No indices to sample for mass accuracy analysis."
+ # Return a structure indicating no analysis was performed
+ empty_report = (
+ mean_ppm_error = NaN, median_ppm_error = NaN, std_ppm_error = NaN,
+ min_ppm_error = NaN, max_ppm_error = NaN, total_matched_peaks = 0,
+ total_spectra_analyzed = 0
+ )
+ results[:global_accuracy] = empty_report
+ results[:suggested_bin_tolerance] = 20.0 # Default
+ return results
+ end
+
+ # Use the existing analyze_mass_accuracy function with the provided indices
+ accuracy_report = analyze_mass_accuracy(msi_data, reference_peaks;
+ spectrum_indices=spectrum_indices)
+
+ results[:global_accuracy] = accuracy_report
+ results[:suggested_bin_tolerance] = calculate_adaptive_bin_tolerance(
+ collect(Iterators.flatten([accuracy_report.mean_ppm_error])) # Simplified - would need actual distribution
+ )
+
+ println(" - Mean PPM error: $(round(accuracy_report.mean_ppm_error, digits=2))")
+ println(" - Suggested bin tolerance: $(round(results[:suggested_bin_tolerance], digits=2)) ppm")
+
+ return results
+end
+
+"""
+ analyze_spatial_regions(msi_data::MSIData, region_masks::Dict, reference_peaks::Dict)::Dict
+
+Analyzes different spatial regions for variations in signal quality and mass accuracy.
+"""
+function analyze_spatial_regions(msi_data::MSIData, region_masks::Dict, reference_peaks::Dict)::Dict
+ results = Dict{Symbol, Any}()
+
+ println(" Analyzing $(length(region_masks)) spatial regions...")
+
+ # Use the existing function
+ regional_reports = find_ppm_error_by_region(msi_data, region_masks, reference_peaks)
+
+ results[:regional_ppm_errors] = regional_reports
+
+ # Calculate regional variations
+ if !isempty(regional_reports)
+ mean_errors = [report.mean_ppm_error for report in values(regional_reports) if isfinite(report.mean_ppm_error)]
+ if !isempty(mean_errors)
+ results[:max_regional_ppm_difference] = maximum(mean_errors) - minimum(mean_errors)
+ end
+ end
+
+ for (region, report) in regional_reports
+ println(" - $region: $(round(report.mean_ppm_error, digits=2)) ppm (n=$(report.total_spectra_analyzed))")
+ end
+
+ return results
+end
+
+"""
+ analyze_peak_characteristics(msi_data::MSIData, sample_size::Int)::Dict
+
+Analyzes peak shape, width, and quality characteristics.
+"""
+function analyze_peak_characteristics(msi_data::MSIData, instrument_analysis::Dict, mass_accuracy_results;
+ spectrum_indices::AbstractVector{Int})::Dict
+ results = Dict{Symbol, Any}()
+
+ # Determine acquisition mode from the analysis results, not from metadata
+ acquisition_mode = get(instrument_analysis, :acquisition_mode, :unknown)
+
+ # Calculate estimated_mean_ppm_error for fallback logic
+ estimated_mean_ppm_error = 30.0 # Default if mass_accuracy_results is nothing or invalid
+ if mass_accuracy_results !== nothing && haskey(mass_accuracy_results, :global_accuracy) && isa(mass_accuracy_results[:global_accuracy], NamedTuple)
+ ppm_error_report = mass_accuracy_results[:global_accuracy]
+ mean_ppm_error_val = get(ppm_error_report, :mean_ppm_error, 30.0)
+ if isfinite(mean_ppm_error_val)
+ estimated_mean_ppm_error = mean_ppm_error_val
+ end
+ end
+
+ if isempty(spectrum_indices)
+ @warn "No indices to sample for peak characteristics analysis."
+ # Provide sensible defaults if no analysis can be run
+ estimated_fwhm = estimated_mean_ppm_error * (acquisition_mode == :profile ? 3 : 2)
+ results[:mean_fwhm_ppm] = estimated_fwhm
+ results[:median_fwhm_ppm] = estimated_fwhm
+ results[:mean_gaussian_r2] = acquisition_mode == :profile ? 0.7 : 0.9
+ results[:peak_resolution_estimate] = 1e6 / estimated_fwhm
+ results[:mean_peaks_per_spectrum] = 0
+ return results
+ end
+
+ if acquisition_mode == PROFILE
+ println(" Analyzing peak characteristics for PROFILE mode from $(length(spectrum_indices)) sample spectra...")
+
+ peak_widths_ppm = Float64[]
+ r_squared_values = Float64[]
+ peak_counts = Int[]
+ fwhm_values = Float64[]
+
+ spectra_analyzed = 0
+ peaks_analyzed = 0
+
+ _iterate_spectra_fast(msi_data, spectrum_indices) do idx, mz, intensity
+ if length(mz) < 10 # Skip spectra with too few points
+ return
+ end
+
+ spectra_analyzed += 1
+ meta = msi_data.spectra_metadata[idx]
+
+ # Detect peaks with lower SNR threshold to find more peaks
+ peaks = detect_peaks_profile(mz, intensity; snr_threshold=2.0)
+ push!(peak_counts, length(peaks))
+
+ if !isempty(peaks)
+ # Analyze the strongest 3 peaks per spectrum
+ sorted_peaks = sort(peaks, by=p->p.intensity, rev=true)
+ for peak in sorted_peaks[1:min(3, length(sorted_peaks))]
+ try
+ # Find the index of the peak in the original mz array
+ peak_idx = argmin(abs.(mz .- peak.mz))
+
+ fwhm_delta_m = calculate_robust_fwhm(mz, intensity, peak_idx)
+
+ if !isnan(fwhm_delta_m) && fwhm_delta_m > 0.001 && fwhm_delta_m < 0.5 # Reasonable range in Da
+ fwhm_ppm = 1e6 * fwhm_delta_m / peak.mz
+ if 5.0 < fwhm_ppm < 500.0 # Reasonable ppm range
+ push!(peak_widths_ppm, fwhm_ppm)
+ push!(fwhm_values, fwhm_delta_m)
+
+ r2 = _fit_gaussian_and_r2(mz, intensity, peak_idx, 5)
+ push!(r_squared_values, r2)
+ peaks_analyzed += 1
+
+ if peaks_analyzed <= 3
+ println("DEBUG: Peak at m/z $(peak.mz), FWHM = $(fwhm_ppm) ppm, R² = $r2")
+ end
+ end
+ end
+ catch e
+ continue
+ end
+ end
+ end
+ end
+
+ println("DEBUG: Analyzed $peaks_analyzed peaks from $spectra_analyzed spectra")
+
+ if !isempty(peak_widths_ppm)
+ results[:mean_fwhm_ppm] = mean(peak_widths_ppm)
+ results[:median_fwhm_ppm] = median(peak_widths_ppm)
+ results[:mean_gaussian_r2] = mean(r_squared_values)
+ results[:peak_resolution_estimate] = 1e6 / results[:mean_fwhm_ppm]
+ println(" - Actual FWHM measurement: $(round(results[:mean_fwhm_ppm], digits=2)) ppm")
+ else
+ # Fallback if no valid peaks found in PROFILE mode
+ estimated_fwhm = estimated_mean_ppm_error * 3 # FWHM typically wider than mass error
+ results[:mean_fwhm_ppm] = estimated_fwhm
+ results[:median_fwhm_ppm] = estimated_fwhm
+ results[:mean_gaussian_r2] = 0.7
+ results[:peak_resolution_estimate] = 1e6 / estimated_fwhm
+ println(" - Estimated FWHM (from mass accuracy): $estimated_fwhm ppm")
+ end
+
+ if !isempty(peak_counts)
+ results[:mean_peaks_per_spectrum] = mean(peak_counts)
+ else
+ results[:mean_peaks_per_spectrum] = 0
+ end
+
+ else # CENTROID mode
+ println(" Analyzing peak characteristics for CENTROID mode from $(length(spectrum_indices)) sample spectra...")
+
+ peak_counts = Int[]
+ # peak_intensities = Float64[] # Not strictly needed for these metrics
+
+ _iterate_spectra_fast(msi_data, spectrum_indices) do idx, mz, intensity
+ if !isempty(mz)
+ push!(peak_counts, length(mz))
+ # append!(peak_intensities, intensity) # If needed for other metrics
+ end
+ end
+
+ if !isempty(peak_counts)
+ results[:mean_peaks_per_spectrum] = mean(peak_counts)
+ results[:median_peaks_per_spectrum] = median(peak_counts)
+ results[:total_peaks_detected] = sum(peak_counts)
+ else
+ results[:mean_peaks_per_spectrum] = 0
+ results[:median_peaks_per_spectrum] = 0
+ results[:total_peaks_detected] = 0
+ end
+
+ # For centroid data, FWHM is estimated from mass accuracy
+ estimated_fwhm = estimated_mean_ppm_error * 2 # Centroid peaks are narrower, factor of 2-3 is common
+ results[:mean_fwhm_ppm] = estimated_fwhm
+ results[:median_fwhm_ppm] = estimated_fwhm # Defaulting median to mean
+ results[:mean_gaussian_r2] = 0.9 # Default for centroid, assuming good peak shapes
+ results[:peak_resolution_estimate] = 1e6 / estimated_fwhm
+
+ println(" - Mean peaks per spectrum: $(round(results[:mean_peaks_per_spectrum], digits=1))")
+ println(" - Estimated peak width: $(round(estimated_fwhm, digits=2)) ppm")
+ println(" - Estimated resolution: $(round(results[:peak_resolution_estimate], digits=0))")
+ end
+
+ # Common prints
+ println(" - Mean Gaussian R²: $(round(get(results, :mean_gaussian_r2, 0.0), digits=3))")
+ println(" - Mean peaks per spectrum: $(round(get(results, :mean_peaks_per_spectrum, 0.0), digits=1))")
+
+ return results
+end
+
+"""
+ generate_preprocessing_recommendations(analysis_results::Dict)::Dict{Symbol, Any}
+
+Generates intelligent preprocessing recommendations based on the analysis results.
+"""
+function generate_preprocessing_recommendations(analysis_results::Dict)::Dict{Symbol, Any}
+ recommendations = Dict{Symbol, Any}()
+
+ inst_analysis = get(analysis_results, :instrument_analysis, Dict())
+ signal_analysis = get(analysis_results, :signal_analysis, Dict())
+ mass_accuracy = get(analysis_results, :mass_accuracy, Dict())
+ peak_analysis = get(analysis_results, :peak_analysis, Dict())
+
+ # Baseline Correction Recommendations
+ recommendations[:baseline_correction] = generate_baseline_recommendations(inst_analysis, signal_analysis)
+
+ # Smoothing Recommendations
+ recommendations[:smoothing] = generate_smoothing_recommendations(inst_analysis, peak_analysis)
+
+ # Peak Picking Recommendations
+ recommendations[:peak_picking] = generate_peak_picking_recommendations(signal_analysis, peak_analysis)
+
+ # Normalization Recommendations
+ recommendations[:normalization] = generate_normalization_recommendations(signal_analysis)
+
+ # Alignment Recommendations
+ recommendations[:alignment] = generate_alignment_recommendations(mass_accuracy, inst_analysis)
+
+ # Binning Recommendations
+ recommendations[:binning] = generate_binning_recommendations(peak_analysis, mass_accuracy)
+
+ return recommendations
+end
+
+function generate_baseline_recommendations(inst_analysis, signal_analysis)
+ recommendations = Dict{Symbol, Any}()
+
+ # Determine baseline correction method based on data characteristics
+ acquisition_mode = get(inst_analysis, :acquisition_mode, :unknown)
+ noise_level = get(signal_analysis, :noise_mean, 1.0)
+
+ if acquisition_mode == :profile
+ recommendations[:method] = "SNIP"
+ recommendations[:window_size] = 200 # Default, can be optimized
+ recommendations[:iterations] = 100
+ else
+ recommendations[:method] = "linear"
+ recommendations[:noise_threshold] = noise_level * 3
+ end
+
+ return recommendations
+end
+
+function generate_smoothing_recommendations(inst_analysis, peak_analysis)
+ recommendations = Dict{Symbol, Any}()
+
+ mz_step = get(inst_analysis, :average_mz_step, 0.01)
+ fwhm_ppm = get(peak_analysis, :mean_fwhm_ppm, 20.0)
+
+ # Convert FWHM from ppm to m/z units for typical m/z
+ typical_mz = 500.0
+ fwhm_mz = typical_mz * fwhm_ppm / 1e6
+
+ # Savitzky-Golay window should be ~FWHM in points
+ window_points = max(5, min(21, round(Int, fwhm_mz / mz_step)))
+ # Ensure odd number
+ window_points = isodd(window_points) ? window_points : window_points + 1
+
+ recommendations[:method] = "Savitzky-Golay"
+ recommendations[:window_size] = window_points
+ recommendations[:polynomial_order] = 3
+
+ return recommendations
+end
+
+function generate_peak_picking_recommendations(signal_analysis, peak_analysis)
+ recommendations = Dict{Symbol, Any}()
+
+ snr_threshold = get(signal_analysis, :suggested_snr, 3.0)
+ fwhm_ppm = get(peak_analysis, :mean_fwhm_ppm, 20.0)
+
+ recommendations[:snr_threshold] = snr_threshold
+ recommendations[:min_peak_width_ppm] = fwhm_ppm * 0.5 # Avoid detecting noise as peaks
+ recommendations[:max_peak_width_ppm] = fwhm_ppm * 3.0 # Avoid merging distinct peaks
+
+ return recommendations
+end
+
+function generate_normalization_recommendations(signal_analysis)
+ recommendations = Dict{Symbol, Any}()
+
+ tic_cv = get(signal_analysis, :tic_cv, 0.5)
+
+ if tic_cv < 0.3 # Low TIC variation
+ recommendations[:method] = "TIC"
+ recommendations[:reason] = "Low TIC variation across spectra"
+ else # High TIC variation
+ recommendations[:method] = "RMS"
+ recommendations[:reason] = "High TIC variation, using robust normalization"
+ end
+
+ return recommendations
+end
+
+function generate_alignment_recommendations(mass_accuracy, inst_analysis)
+ recommendations = Dict{Symbol, Any}()
+
+ calibration_status = get(inst_analysis, :calibration_status, :uncalibrated)
+
+ mean_ppm_error_val = 10.0 # Default value
+ if mass_accuracy !== nothing && haskey(mass_accuracy, :global_accuracy) && isa(mass_accuracy[:global_accuracy], NamedTuple)
+ ppm_error_report = mass_accuracy[:global_accuracy]
+ mean_ppm_error_val = get(ppm_error_report, :mean_ppm_error, 10.0)
+ end
+
+ if calibration_status == :uncalibrated || mean_ppm_error_val > 20.0
+ recommendations[:method] = "reference_based"
+ recommendations[:max_ppm_shift] = 50.0
+ recommendations[:required] = true
+ else
+ recommendations[:method] = "none"
+ recommendations[:required] = false
+ end
+
+ return recommendations
+end
+
+function generate_binning_recommendations(peak_analysis, mass_accuracy)
+ recommendations = Dict{Symbol, Any}()
+
+ fwhm_ppm = get(peak_analysis, :mean_fwhm_ppm, 20.0)
+
+ # Handle the case when mass_accuracy is nothing or has NaN values
+ suggested_tolerance = 20.0 # Default value
+ if mass_accuracy !== nothing
+ tol = get(mass_accuracy, :suggested_bin_tolerance, 20.0)
+ if !isnan(tol) && isfinite(tol)
+ suggested_tolerance = tol
+ end
+ end
+
+ # Ensure fwhm_ppm is valid
+ if isnan(fwhm_ppm) || !isfinite(fwhm_ppm)
+ fwhm_ppm = 20.0
+ end
+
+ # Bin width should be ~FWHM/2 to preserve resolution while reducing data size
+ bin_width_ppm = max(fwhm_ppm * 0.5, suggested_tolerance)
+
+ recommendations[:method] = "adaptive"
+ recommendations[:bin_width_ppm] = bin_width_ppm
+ recommendations[:min_peaks_per_bin] = 3
+
+ return recommendations
+end
+
+# =============================================================================
+# Pre-Analysis Pipeline for Auto Parameter Determination
+# =============================================================================
+
+"""
+ run_preprocessing_analysis(msi_data::MSIData;
+ reference_peaks::Dict{Float64, String}=Dict(),
+ region_masks::Dict=Dict(),
+ sample_size::Int=100)::Dict{Symbol, Any}
+
+Runs a comprehensive pre-analysis pipeline to determine optimal preprocessing parameters.
+This function analyzes the dataset to provide intelligent defaults for preprocessing steps.
+"""
+function run_preprocessing_analysis(msi_data::MSIData;
+ reference_peaks::Dict{Float64, String}=Dict{Float64, String}(),
+ region_masks::Dict{Symbol, BitMatrix}=Dict{Symbol, BitMatrix}(),
+ sample_size::Int=100,
+ mask_path::Union{String, Nothing}=nothing,
+ spectrum_indices::Union{AbstractVector{Int}, Nothing}=nothing)::Dict{Symbol, Any}
+
+ println("="^60)
+ println("RUNNING PRE-ANALYSIS PIPELINE")
+ println("="^60)
+
+ analysis_results = Dict{Symbol, Any}()
+
+ local all_available_indices::AbstractVector{Int}
+
+ if spectrum_indices !== nothing
+ println("Using provided list of $(length(spectrum_indices)) spectrum indices.")
+ all_available_indices = spectrum_indices
+ elseif mask_path !== nothing
+ println("Applying mask from: $(mask_path)")
+ try
+ mask_matrix = load_and_prepare_mask(mask_path, msi_data.image_dims)
+ masked_indices_set = get_masked_spectrum_indices(msi_data, mask_matrix)
+ all_available_indices = collect(masked_indices_set)
+ println("Mask applied. $(length(all_available_indices)) spectra are within the masked region.")
+ catch e
+ @error "Failed to load or apply mask: $e. Proceeding without mask."
+ all_available_indices = 1:length(msi_data.spectra_metadata)
+ end
+ else
+ all_available_indices = 1:length(msi_data.spectra_metadata)
+ end
+
+ if isempty(all_available_indices)
+ @warn "No spectra available for analysis (after applying mask/filter). Returning empty results."
+ return analysis_results
+ end
+
+ # Create a single sample set from the available indices
+ num_available = length(all_available_indices)
+ indices_to_sample = if num_available > sample_size
+ # Use StatsBase.sample for sampling without replacement
+ sample(all_available_indices, sample_size, replace=false)
+ else
+ # Use all available indices if they are fewer than the sample size
+ collect(all_available_indices)
+ end
+
+ # Ensure analytics are precomputed
+ if !is_set(msi_data.analytics_ready)
+ println("Pre-computing basic analytics...")
+ precompute_analytics(msi_data)
+ end
+
+ # Phase 1: Instrument and Data Characteristics
+ println("\n--- Phase 1: Instrument & Data Characteristics ---")
+ instrument_analysis = analyze_instrument_characteristics(msi_data, sample_indices=indices_to_sample)
+ analysis_results[:instrument_analysis] = instrument_analysis
+
+ # Phase 2: Noise and Signal Quality Analysis
+ println("\n--- Phase 2: Noise & Signal Quality Analysis ---")
+ signal_analysis = analyze_signal_quality(msi_data, sample_indices=indices_to_sample)
+ analysis_results[:signal_analysis] = signal_analysis
+
+ # Phase 3: Mass Accuracy Analysis
+ println("\n--- Phase 3: Mass Accuracy Analysis ---")
+ if !isempty(reference_peaks)
+ mass_accuracy_analysis = analyze_mass_accuracy_global(msi_data, reference_peaks, spectrum_indices=indices_to_sample)
+ analysis_results[:mass_accuracy] = mass_accuracy_analysis
+ else
+ println("No reference peaks provided - skipping mass accuracy analysis")
+ analysis_results[:mass_accuracy] = nothing
+ end
+
+ # Phase 4: Spatial Region Analysis (if masks provided)
+ # This function is not affected by the global mask, as it analyzes specific, named regions.
+ println("\n--- Phase 4: Spatial Region Analysis ---")
+ if !isempty(region_masks)
+ regional_analysis = analyze_spatial_regions(msi_data, region_masks, reference_peaks)
+ analysis_results[:regional_analysis] = regional_analysis
+ else
+ println("No region masks provided - skipping regional analysis")
+ analysis_results[:regional_analysis] = nothing
+ end
+
+ # Phase 5: Peak Characteristics Analysis
+ println("\n--- Phase 5: Peak Characteristics Analysis ---")
+ peak_analysis = analyze_peak_characteristics(msi_data, instrument_analysis, analysis_results[:mass_accuracy], spectrum_indices=indices_to_sample)
+ analysis_results[:peak_analysis] = peak_analysis
+
+ # Phase 6: Generate Preprocessing Recommendations
+ println("\n--- Phase 6: Generating Preprocessing Recommendations ---")
+ recommendations = generate_preprocessing_recommendations(analysis_results)
+ analysis_results[:recommendations] = recommendations
+
+ # Store in MSIData object
+ msi_data.preprocessing_hints = recommendations
+
+ println("\n" * "="^60)
+ println("PRE-ANALYSIS COMPLETE")
+ println("="^60)
+
+ return analysis_results
+end
+
+"""
+ calculate_ppm_error(measured_mz::Real, theoretical_mz::Real) -> Float64
+
+Calculates the mass accuracy error in parts-per-million (PPM) between a measured
+and a theoretical m/z value.
+
+# Arguments
+- `measured_mz::Real`: The experimentally measured m/z value.
+- `theoretical_mz::Real`: The known, theoretical m/z value of a compound.
+
+# Returns
+- `Float64`: The calculated PPM error. Returns `Inf` if `theoretical_mz` is zero.
+
+# Formula
+`PPM = 10^6 * |measured_mz - theoretical_mz| / theoretical_mz`
+
+# Example
+```julia
+calculate_ppm_error(100.005, 100.0) # returns 50.0
+```
+"""
+function calculate_ppm_error(measured_mz::Real, theoretical_mz::Real)
+ if theoretical_mz == 0
+ return Inf
+ end
+ return 1e6 * abs(Float64(measured_mz) - Float64(theoretical_mz)) / Float64(theoretical_mz)
+end
+
+"""
+ calculate_ppm_error_bulk(measured_mz::Vector{<:Real}, theoretical_mz::Vector{<:Real}) -> Vector{Float64}
+
+Calculates PPM errors for multiple pairs of measured and theoretical mass values.
+
+# Arguments
+- `measured_mz::Vector{<:Real}`: A vector of experimentally measured m/z values.
+- `theoretical_mz::Vector{<:Real}`: A vector of known, theoretical m/z values.
+
+# Returns
+- `Vector{Float64}`: A vector containing the calculated PPM error for each pair.
+"""
+function calculate_ppm_error_bulk(measured_mz::Vector{Real}, theoretical_mz::Vector{Real})
+ return [calculate_ppm_error(m, t) for (m, t) in zip(measured_mz, theoretical_mz)]
+end
+
+"""
+ calculate_resolution_fwhm(mz::Real, profile_mz::AbstractVector{<:Real}, profile_intensity::AbstractVector{<:Real}) -> Float64
+
+Calculates the mass resolution of a peak in profile-mode data using the Full Width at
+Half Maximum (FWHM) method. Resolution is a measure of an instrument's ability to
+distinguish between two peaks of slightly different mass-to-charge ratios.
+
+# Arguments
+- `mz::Real`: The m/z value of the peak's centroid.
+- `profile_mz::AbstractVector{<:Real}`: The full m/z array from the profile-mode spectrum.
+- `profile_intensity::AbstractVector{<:Real}`: The full intensity array from the profile-mode spectrum.
+
+# Returns
+- `Float64`: The calculated resolution (`m / Δm`). Returns `NaN` if the FWHM cannot be determined (e.g., peak is at the edge of the spectrum).
+
+# Formula
+`Resolution = m / Δm`, where `Δm` is the FWHM.
+"""
+function calculate_resolution_fwhm(mz::Real, profile_mz::AbstractVector{<:Real},
+ profile_intensity::AbstractVector{<:Real})
+
+ # Find peak center index
+ peak_idx = argmin(abs.(profile_mz .- mz))
+ peak_height = Float64(profile_intensity[peak_idx])
+ half_max = peak_height / 2
+
+ # Find left half-maximum point (interpolate for accuracy)
+ left_idx = find_last_below(profile_intensity[1:peak_idx], half_max)
+ if left_idx == 0 || left_idx == length(profile_intensity[1:peak_idx])
+ return NaN
+ end
+
+ # Linear interpolation for left FWHM
+ x1, x2 = Float64(profile_mz[left_idx]), Float64(profile_mz[left_idx+1])
+ y1, y2 = Float64(profile_intensity[left_idx]), Float64(profile_intensity[left_idx+1])
+ denominator = y2 - y1
+ if denominator == 0
+ left_fwhm = NaN
+ else
+ left_fwhm = x1 + (x2 - x1) * (half_max - y1) / denominator
+ end
+
+ # Find right half-maximum point
+ right_slice = profile_intensity[peak_idx:end]
+ right_offset = find_first_below(right_slice, half_max)
+ if right_offset == 0 || right_offset == length(right_slice)
+ return NaN
+ end
+
+ right_idx = peak_idx + right_offset - 1
+ x1, x2 = Float64(profile_mz[right_idx-1]), Float64(profile_mz[right_idx])
+ y1, y2 = Float64(profile_intensity[right_idx-1]), Float64(profile_intensity[right_idx])
+ denominator = y2 - y1
+ if denominator == 0
+ right_fwhm = NaN
+ else
+ right_fwhm = x1 + (x2 - x1) * (half_max - y1) / denominator
+ end
+
+ fwhm = right_fwhm - left_fwhm
+ return fwhm > 0 ? Float64(mz) / fwhm : NaN
+end
+
+# Helper functions for FWHM calculation
+function find_last_below(v::AbstractVector{<:Real}, threshold::Real)
+ for i in length(v):-1:2
+ if v[i] >= threshold && v[i-1] < threshold
+ return i-1
+ end
+ end
+ return 0
+end
+
+function find_first_below(v::AbstractVector{<:Real}, threshold::Real)
+ for i in 1:(length(v)-1)
+ if v[i] >= threshold && v[i+1] < threshold
+ return i+1
+ end
+ end
+ return 0
+end
+
+"""
+ _calculate_fwhm_delta_m(mz::AbstractVector{<:Real}, intensity::AbstractVector{<:Real}, peak_idx::Int) -> Float64
+
+Calculates the Full Width at Half Maximum (FWHM) in m/z units (Δm) for a peak.
+Returns `NaN` if FWHM cannot be determined.
+"""
+function _calculate_fwhm_delta_m(mz::AbstractVector{<:Real}, intensity::AbstractVector{<:Real}, peak_idx::Int)
+ peak_height = Float64(intensity[peak_idx])
+ half_max = peak_height / 2
+
+ # Find left half-maximum point
+ left_idx = find_last_below(intensity[1:peak_idx], half_max)
+ if left_idx == 0 || left_idx == length(intensity[1:peak_idx])
+ return NaN
+ end
+
+ # Linear interpolation for left FWHM m/z
+ x1, x2 = Float64(mz[left_idx]), Float64(mz[left_idx+1])
+ y1, y2 = Float64(intensity[left_idx]), Float64(intensity[left_idx+1])
+ denominator = y2 - y1
+ if denominator == 0
+ left_fwhm_mz = NaN
+ else
+ left_fwhm_mz = x1 + (x2 - x1) * (half_max - y1) / denominator
+ end
+
+ # Find right half-maximum point
+ right_slice = intensity[peak_idx:end]
+ right_offset = find_first_below(right_slice, half_max)
+ if right_offset == 0 || right_offset == length(right_slice)
+ return NaN
+ end
+
+ right_idx = peak_idx + right_offset - 1
+ x1, x2 = Float64(mz[right_idx-1]), Float64(mz[right_idx])
+ y1, y2 = Float64(intensity[right_idx-1]), Float64(intensity[right_idx])
+ denominator = y2 - y1
+ if denominator == 0
+ right_fwhm_mz = NaN
+ else
+ right_fwhm_mz = x1 + (x2 - x1) * (half_max - y1) / denominator
+ end
+
+ return right_fwhm_mz - left_fwhm_mz
+end
+
+"""
+ _fit_gaussian_and_r2(mz::AbstractVector{<:Real}, intensity::AbstractVector{<:Real}, peak_idx::Int, half_window::Int) -> Float64
+
+Estimates Gaussian parameters for a peak and returns a pseudo R^2 value.
+This is an approximation and not a full non-linear least squares fit.
+"""
+function _fit_gaussian_and_r2(mz::AbstractVector{<:Real}, intensity::AbstractVector{<:Real}, peak_idx::Int, half_window::Int)
+ n = length(mz)
+ if n < 3 || peak_idx <= 0 || peak_idx > n
+ return 0.0
+ end
+
+ # Define the region around the peak
+ start_idx = max(1, peak_idx - half_window)
+ end_idx = min(n, peak_idx + half_window)
+
+ # Ensure there's enough data to fit
+ if (end_idx - start_idx + 1) < 3
+ return 0.0
+ end
+
+ x_data = mz[start_idx:end_idx]
+ y_data = intensity[start_idx:end_idx]
+
+ # Estimate Gaussian parameters
+ # Amplitude (A): peak intensity
+ A_est = intensity[peak_idx]
+ # Mean (μ): m/z at peak intensity
+ mu_est = mz[peak_idx]
+ # Standard deviation (σ): related to FWHM. FWHM = 2 * sqrt(2 * ln(2)) * σ ≈ 2.355 * σ
+ # So, σ ≈ FWHM / 2.355
+ fwhm_delta_m = _calculate_fwhm_delta_m(mz, intensity, peak_idx)
+ if isnan(fwhm_delta_m) || fwhm_delta_m <= 0
+ return 0.0 # Cannot estimate sigma without a valid FWHM
+ end
+ sigma_est = fwhm_delta_m / 2.355
+
+ # If sigma is too small, it might lead to division by zero or very sharp peaks
+ if sigma_est < eps(Float64)
+ return 0.0
+ end
+
+ # Gaussian function
+ gaussian(x, A, mu, sigma) = A * exp.(-(x .- mu).^2 ./ (2 * sigma^2))
+
+ # Generate estimated Gaussian curve
+ y_est = gaussian(x_data, A_est, mu_est, sigma_est)
+
+ # Calculate pseudo R-squared
+ # R^2 = 1 - (SS_res / SS_tot)
+ # SS_res = sum((y_data - y_est).^2)
+ # SS_tot = sum((y_data - mean(y_data)).^2)
+
+ SS_res = sum((y_data .- y_est).^2)
+ SS_tot = sum((y_data .- mean(y_data)).^2)
+
+ if SS_tot == 0
+ return 1.0 # Perfect fit if all y_data are the same
+ end
+
+ r_squared = 1.0 - (SS_res / SS_tot)
+ return max(0.0, r_squared) # R^2 can be negative if fit is worse than mean, cap at 0
+end
+
+"""
+ calculate_robust_fwhm(mz::AbstractVector{<:Real}, intensity::AbstractVector{<:Real}, peak_idx::Int) -> Float64
+
+A more robust FWHM calculation that handles edge cases better.
+"""
+function calculate_robust_fwhm(mz::AbstractVector{<:Real}, intensity::AbstractVector{<:Real}, peak_idx::Int)
+ n = length(mz)
+ if n < 5 || peak_idx < 3 || peak_idx > n-2
+ return NaN
+ end
+
+ peak_height = Float64(intensity[peak_idx])
+ half_max = peak_height / 2.0
+
+ # Find left half-maximum with bounds checking
+ left_idx = peak_idx
+ while left_idx > 1 && intensity[left_idx] >= half_max
+ left_idx -= 1
+ end
+
+ if left_idx == 1 || left_idx >= n-1
+ return NaN
+ end
+
+ # Linear interpolation for left FWHM
+ x1, x2 = Float64(mz[left_idx]), Float64(mz[left_idx+1])
+ y1, y2 = Float64(intensity[left_idx]), Float64(intensity[left_idx+1])
+
+ if y2 == y1 # Avoid division by zero
+ left_fwhm_mz = x1
+ else
+ left_fwhm_mz = x1 + (x2 - x1) * (half_max - y1) / (y2 - y1)
+ end
+
+ # Find right half-maximum
+ right_idx = peak_idx
+ while right_idx < n && intensity[right_idx] >= half_max
+ right_idx += 1
+ end
+
+ if right_idx == n || right_idx <= 2
+ return NaN
+ end
+
+ # Linear interpolation for right FWHM
+ x1, x2 = Float64(mz[right_idx-1]), Float64(mz[right_idx])
+ y1, y2 = Float64(intensity[right_idx-1]), Float64(intensity[right_idx])
+
+ if y2 == y1 # Avoid division by zero
+ right_fwhm_mz = x1
+ else
+ right_fwhm_mz = x1 + (x2 - x1) * (half_max - y1) / (y2 - y1)
+ end
+
+ fwhm = right_fwhm_mz - left_fwhm_mz
+
+ # Validate result
+ if fwhm <= 0 || !isfinite(fwhm) || fwhm > 1.0 # Unreasonably large
+ return NaN
+ end
+
+ return fwhm
+end
+
+"""
+ main_precalculation(msi_data::MSIData; ...)
+
+Runs the non-verbose pre-analysis pipeline and returns a dictionary of recommended
+preprocessing parameters based on data characteristics and heuristics.
+
+This function serves as a quiet entry point to the analysis engine, translating the
+analytical results into a concrete set of parameters for a preprocessing pipeline.
+It also categorizes parameters that cannot be automatically determined.
+
+# Arguments
+- `msi_data::MSIData`: The main MSI data object.
+- `reference_peaks::Dict`: Optional dictionary of reference m/z values for mass accuracy analysis.
+- `region_masks::Dict`: Optional dictionary of named `BitMatrix` masks for regional analysis.
+- `sample_size::Int`: The number of spectra to sample for statistical analysis.
+- `mask_path::String`: Optional path to a PNG mask file to restrict the analysis to a specific ROI.
+- `spectrum_indices::AbstractVector{Int}`: Optional vector of spectrum indices to restrict analysis to.
+
+# Returns
+- A `Dict` with two keys:
+ - `"recommended_parameters"`: A `Dict{Symbol, Any}` of suggested parameter values.
+ - `"unsupported_parameters"`: A `Dict` categorizing parameters that could not be determined.
+"""
+function main_precalculation(msi_data::MSIData;
+ reference_peaks::Dict{Float64, String}=Dict{Float64, String}(),
+ region_masks::Dict{Symbol, BitMatrix}=Dict{Symbol, BitMatrix}(),
+ sample_size::Int=100,
+ mask_path::Union{String, Nothing}=nothing,
+ spectrum_indices::Union{AbstractVector{Int}, Nothing}=nothing)::Dict
+
+ # --- Run analysis pipeline quietly ---
+ local analysis_results
+ original_stdout = stdout
+ # Redirect stdout to the system's null device to robustly silence output
+ null_stream = open(Sys.iswindows() ? "nul" : "/dev/null", "w")
+ redirect_stdout(null_stream)
+ try
+ analysis_results = run_preprocessing_analysis(msi_data,
+ reference_peaks=reference_peaks,
+ region_masks=region_masks,
+ sample_size=sample_size,
+ mask_path=mask_path,
+ spectrum_indices=spectrum_indices
+ )
+ finally
+ redirect_stdout(original_stdout)
+ close(null_stream)
+ end
+
+ if isempty(analysis_results)
+ @warn "Preprocessing analysis returned no results. Cannot generate recommendations."
+ return Dict() # Return an empty dictionary if no results
+ end
+
+ # --- Safely get nested dictionaries ---
+ recs = get(analysis_results, :recommendations, Dict())
+ signal_analysis = get(analysis_results, :signal_analysis, Dict())
+ peak_analysis = get(analysis_results, :peak_analysis, Dict())
+ mass_accuracy = get(analysis_results, :mass_accuracy, Dict())
+ inst_analysis = get(analysis_results, :instrument_analysis, Dict())
+
+ # Initialize parameter dictionaries for each step
+ cal_params = Dict{Symbol, Any}()
+ sm_params = Dict{Symbol, Any}()
+ bc_params = Dict{Symbol, Any}()
+ norm_params = Dict{Symbol, Any}()
+ pp_params = Dict{Symbol, Any}()
+ pa_params = Dict{Symbol, Any}()
+ ps_params = Dict{Symbol, Any}()
+ pb_params = Dict{Symbol, Any}()
+
+ # --- Populate Parameters for each step ---
+
+ # Calibration & Alignment (Note: These are intertwined in the current logic)
+ calibration_required = false
+ mean_ppm_error = get(get(mass_accuracy, :global_accuracy, Dict()), :mean_ppm_error, NaN)
+ suggested_bin_tol = get(mass_accuracy, :suggested_bin_tolerance, NaN)
+
+ if !isempty(recs) && haskey(recs, :alignment) && get(recs[:alignment], :required, false)
+ calibration_required = true
+ cal_params[:method] = :internal_standards
+ cal_params[:fit_order] = 2 # Default from struct
+ if isfinite(suggested_bin_tol)
+ cal_params[:ppm_tolerance] = suggested_bin_tol
+ else
+ cal_params[:ppm_tolerance] = nothing
+ end
+ cal_params[:internal_standards] = nothing # User input dependent
+ cal_params[:base_peak_mz_references] = nothing # User input dependent
+
+ if isfinite(mean_ppm_error)
+ pa_params[:method] = mean_ppm_error > 30.0 ? :lowess : :linear
+ else
+ pa_params[:method] = :lowess # Default
+ end
+
+ tic_cv = get(signal_analysis, :tic_cv, NaN)
+ if isfinite(tic_cv) && get(pa_params, :method, :none) == :lowess
+ pa_params[:span] = round(max(0.3, min(0.8, 1.0 - tic_cv / 2)), digits=2)
+ else
+ pa_params[:span] = nothing
+ end
+
+ if isfinite(suggested_bin_tol)
+ pa_params[:tolerance] = suggested_bin_tol
+ pa_params[:tolerance_unit] = :ppm
+ else
+ pa_params[:tolerance] = nothing
+ pa_params[:tolerance_unit] = :ppm # Default
+ end
+ pa_params[:max_shift_ppm] = get(recs[:alignment], :max_ppm_shift, 50.0)
+ pa_params[:min_matched_peaks] = nothing # User input dependent
+ else
+ cal_params[:method] = :none
+ cal_params[:ppm_tolerance] = nothing
+ cal_params[:fit_order] = nothing
+ cal_params[:internal_standards] = nothing
+ cal_params[:base_peak_mz_references] = nothing
+
+ pa_params[:method] = :none
+ pa_params[:span] = nothing
+ pa_params[:tolerance] = nothing
+ pa_params[:tolerance_unit] = :ppm
+ pa_params[:max_shift_ppm] = nothing
+ pa_params[:min_matched_peaks] = nothing
+ end
+
+ # Smoothing
+ if !isempty(recs) && haskey(recs, :smoothing)
+ sm_rec = recs[:smoothing]
+ if get(sm_rec, :method, "") == "Savitzky-Golay"
+ sm_params[:method] = :savitzky_golay
+ sm_params[:window] = get(sm_rec, :window_size, nothing)
+ sm_params[:order] = get(sm_rec, :polynomial_order, nothing)
+ else
+ sm_params[:method] = :none
+ sm_params[:window] = nothing
+ sm_params[:order] = nothing
+ end
+ else
+ sm_params[:method] = :none
+ sm_params[:window] = nothing
+ sm_params[:order] = nothing
+ end
+
+ # Baseline Correction
+ if !isempty(recs) && haskey(recs, :baseline_correction)
+ bl_rec = recs[:baseline_correction]
+ if get(bl_rec, :method, "") == "SNIP"
+ bc_params[:method] = :snip
+ bc_params[:iterations] = get(bl_rec, :iterations, nothing)
+ else
+ bc_params[:method] = :none
+ bc_params[:iterations] = nothing
+ end
+ else
+ bc_params[:method] = :none
+ bc_params[:iterations] = nothing
+ end
+
+ if get(inst_analysis, :acquisition_mode, :unknown) == :profile
+ mean_fwhm_ppm = get(peak_analysis, :mean_fwhm_ppm, NaN)
+ avg_mz_step = get(inst_analysis, :average_mz_step, NaN)
+ if isfinite(mean_fwhm_ppm) && isfinite(avg_mz_step) && avg_mz_step > 0
+ fwhm_mz = 500.0 * mean_fwhm_ppm / 1e6 # At typical m/z 500
+ bc_params[:window] = ceil(Int, fwhm_mz / avg_mz_step * 2)
+ else
+ bc_params[:window] = nothing
+ end
+ else
+ bc_params[:window] = nothing
+ end
+
+ # Normalization
+ if !isempty(recs) && haskey(recs, :normalization)
+ method = get(recs[:normalization], :method, "")
+ if method == "TIC"
+ norm_params[:method] = :tic
+ elseif method == "RMS"
+ norm_params[:method] = :rms
+ else
+ norm_params[:method] = :none
+ end
+ else
+ norm_params[:method] = :none
+ end
+
+ # Peak Picking
+ # Robust method selection with fallback
+ acquisition_mode = get(inst_analysis, :acquisition_mode, :profile) # Default to profile
+ pp_params[:method] = acquisition_mode == :profile ? :profile : :centroid
+
+ if !isempty(recs) && haskey(recs, :peak_picking)
+ pk_rec = recs[:peak_picking]
+ pp_params[:snr_threshold] = get(pk_rec, :snr_threshold, 3.0) # Default to 3.0
+ pp_params[:min_peak_width_ppm] = get(pk_rec, :min_peak_width_ppm, nothing)
+ pp_params[:max_peak_width_ppm] = get(pk_rec, :max_peak_width_ppm, nothing)
+ else
+ pp_params[:snr_threshold] = 3.0 # Default to 3.0
+ pp_params[:min_peak_width_ppm] = nothing
+ pp_params[:max_peak_width_ppm] = nothing
+ end
+
+ # Robust prominence calculation with safety cap
+ estimated_noise = get(signal_analysis, :noise_mean, NaN)
+ if isfinite(estimated_noise)
+ # Never be more aggressive than 0.05, a reasonable upper limit
+ calculated_prominence = round(estimated_noise * 2, digits=4)
+ pp_params[:min_peak_prominence] = min(calculated_prominence, 0.05)
+ else
+ # Fallback to a safe, non-aggressive value if noise couldn't be estimated
+ pp_params[:min_peak_prominence] = 0.01
+ end
+
+ if isfinite(suggested_bin_tol)
+ pp_params[:merge_peaks_tolerance] = round(suggested_bin_tol / 2, digits=4)
+ else
+ pp_params[:merge_peaks_tolerance] = nothing
+ end
+
+ # Robust half_window calculation with safety floor
+ mean_fwhm_ppm = get(peak_analysis, :mean_fwhm_ppm, NaN)
+ avg_mz_step = get(inst_analysis, :average_mz_step, NaN)
+ if isfinite(mean_fwhm_ppm) && isfinite(avg_mz_step) && avg_mz_step > 0
+ fwhm_mz = 500.0 * mean_fwhm_ppm / 1e6 # At typical m/z 500
+ window_points = fwhm_mz / avg_mz_step
+ calculated_half_window = ceil(Int, window_points / 2)
+ # Ensure half_window is at least 3, a safe absolute minimum
+ pp_params[:half_window] = max(calculated_half_window, 3)
+ else
+ # Fallback if calculation is not possible
+ pp_params[:half_window] = 5
+ end
+
+ # Robust R^2 calculation with floor
+ mean_r2 = get(peak_analysis, :mean_gaussian_r2, NaN)
+ if isfinite(mean_r2)
+ pp_params[:min_peak_shape_r2] = round(max(0.5, mean_r2 * 0.8), digits=2)
+ else
+ # Fallback to a reasonable default
+ pp_params[:min_peak_shape_r2] = 0.6
+ end
+
+
+ # Peak Selection
+ if haskey(pp_params, :snr_threshold) && pp_params[:snr_threshold] !== nothing
+ ps_params[:min_snr] = pp_params[:snr_threshold]
+ else
+ ps_params[:min_snr] = nothing
+ end
+ if !isempty(peak_analysis)
+ mean_fwhm = get(peak_analysis, :mean_fwhm_ppm, NaN)
+ if isfinite(mean_fwhm)
+ ps_params[:min_fwhm_ppm] = round(mean_fwhm * 0.5, digits=2)
+ ps_params[:max_fwhm_ppm] = round(mean_fwhm * 2.5, digits=2)
+ else
+ ps_params[:min_fwhm_ppm] = nothing
+ ps_params[:max_fwhm_ppm] = nothing
+ end
+ if isfinite(mean_r2)
+ ps_params[:min_shape_r2] = round(max(0.5, mean_r2 * 0.8), digits=2)
+ else
+ ps_params[:min_shape_r2] = nothing
+ end
+ else
+ ps_params[:min_fwhm_ppm] = nothing
+ ps_params[:max_fwhm_ppm] = nothing
+ ps_params[:min_shape_r2] = nothing
+ end
+ ps_params[:frequency_threshold] = nothing # User input dependent / Hard to determine
+ ps_params[:correlation_threshold] = nothing # Hard to determine
+
+
+ # Peak Binning
+ if !isempty(recs) && haskey(recs, :binning)
+ bin_rec = recs[:binning]
+ if get(bin_rec, :method, "") == "adaptive"
+ pb_params[:method] = :adaptive
+ bin_width = get(bin_rec, :bin_width_ppm, NaN)
+ if isfinite(bin_width)
+ pb_params[:tolerance] = bin_width
+ pb_params[:max_bin_width_ppm] = round(bin_width * 3, digits=2)
+ else
+ pb_params[:tolerance] = nothing
+ pb_params[:max_bin_width_ppm] = nothing
+ end
+ pb_params[:tolerance_unit] = :ppm
+ pb_params[:min_peak_per_bin] = get(bin_rec, :min_peaks_per_bin, nothing)
+ pb_params[:intensity_weighted_centers] = true # Default from struct
+ pb_params[:num_uniform_bins] = nothing # User input dependent
+ pb_params[:frequency_threshold] = nothing # Hard to determine
+ else
+ pb_params[:method] = :none
+ pb_params[:tolerance] = nothing
+ pb_params[:max_bin_width_ppm] = nothing
+ pb_params[:tolerance_unit] = nothing
+ pb_params[:min_peak_per_bin] = nothing
+ pb_params[:intensity_weighted_centers] = true
+ pb_params[:num_uniform_bins] = nothing
+ pb_params[:frequency_threshold] = nothing
+ end
+ else
+ pb_params[:method] = :none
+ pb_params[:tolerance] = nothing
+ pb_params[:max_bin_width_ppm] = nothing
+ pb_params[:tolerance_unit] = nothing
+ pb_params[:min_peak_per_bin] = nothing
+ pb_params[:intensity_weighted_centers] = true
+ pb_params[:num_uniform_bins] = nothing
+ pb_params[:frequency_threshold] = nothing
+ end
+
+ return Dict(
+ :Calibration => cal_params,
+ :Smoothing => sm_params,
+ :BaselineCorrection => bc_params,
+ :Normalization => norm_params,
+ :PeakPicking => pp_params,
+ :PeakAlignment => pa_params,
+ :PeakSelection => ps_params,
+ :PeakBinningParams => pb_params
+ )
+end
diff --git a/src/Preprocessing.jl b/src/Preprocessing.jl
index a1a8cff..9200870 100644
--- a/src/Preprocessing.jl
+++ b/src/Preprocessing.jl
@@ -1,9 +1,10 @@
# src/Preprocessing.jl
"""
This module provides a comprehensive workflow for mass spectrometry imaging (MSI) data
-preprocessing, inspired by the functionality of the R package MALDIquant. It includes
-functions for quality control, intensity transformation, smoothing, baseline correction,
-normalization, peak picking, alignment, and feature matrix generation.
+preprocessing, inspired by the functionality of the R package MALDIquant and Cardinal.
+It includes functions for quality control, intensity transformation, smoothing,
+baseline correction, normalization, peak picking, alignment, and feature matrix
+generation.
"""
# =============================================================================
@@ -17,7 +18,9 @@ using Dates # For now()
using CSV # For writing CSV files
using DataFrames # For creating dataframes
using ContinuousWavelets # For CWT peak detection
+using ImageFiltering # For localmaxima in detect_peaks_wavelet
using Interpolations # For calibration
+using Loess # For robust peak alignment
# =============================================================================
# Data Structures
@@ -35,19 +38,165 @@ struct FeatureMatrix
end
"""
- QCParameters
-
-A struct to hold comprehensive quality control parameters for validating spectra.
+A mutable struct to hold spectrum data. Using a mutable struct allows
+in-place modification of fields (like intensity or m/z), which dramatically
+reduces memory allocations compared to creating new immutable tuples at each step.
"""
-struct QCParameters
- min_tic_threshold::Float64 # Minimum total ion current
- max_tic_threshold::Float64 # Maximum TIC (saturation check)
- snr_threshold::Float64 # Minimum signal-to-noise for a spectrum to be considered valid
- peak_count_range::Tuple{Int,Int} # Acceptable peak count range (min, max)
- spatial_consistency_threshold::Float64 # For MSI spatial coherence (stub)
- calibration_accuracy_ppm::Float64 # Maximum allowed PPM error for calibration lock masses
+mutable struct MutableSpectrum
+ id::Int
+ mz::AbstractVector{Float64}
+ intensity::AbstractVector{Float64}
+ peaks::Vector{NamedTuple{(:mz, :intensity, :fwhm, :shape_r2, :snr, :prominence), NTuple{6, Float64}}}
end
+# --- Core Pipeline Structs ---
+
+"""
+An abstract type for all preprocessing steps. Allows for a modular and extensible
+pipeline where users can define a sequence of operations.
+"""
+abstract type AbstractPreprocessingStep end
+
+# --- Preprocessing Step Definitions ---
+
+"""
+ Calibration(; method=:internal_standards, ...)
+
+A preprocessing step for mass calibration. Set a parameter to `nothing` to use an
+auto-determined value from the data where applicable.
+"""
+struct Calibration <: AbstractPreprocessingStep
+ method::Symbol
+ internal_standards::Union{Dict{Float64, String}, Nothing}
+ base_peak_mz_references::Union{Vector{Float64}, Nothing}
+ ppm_tolerance::Union{Float64, Nothing}
+ fit_order::Int
+
+ function Calibration(; method=:internal_standards, internal_standards=nothing, base_peak_mz_references=nothing, ppm_tolerance=nothing, fit_order=2)
+ new(method, internal_standards, base_peak_mz_references, ppm_tolerance, fit_order)
+ end
+end
+
+"""
+ BaselineCorrection(; method=:snip, ...)
+
+A preprocessing step for baseline correction.
+"""
+struct BaselineCorrection <: AbstractPreprocessingStep
+ method::Symbol
+ iterations::Union{Int, Nothing} # For SNIP
+ window::Union{Int, Nothing} # For median
+
+ function BaselineCorrection(; method=:snip, iterations=nothing, window=nothing)
+ new(method, iterations, window)
+ end
+end
+
+"""
+ Smoothing(; method=:savitzky_golay, ...)
+
+A preprocessing step for spectral smoothing.
+"""
+struct Smoothing <: AbstractPreprocessingStep
+ method::Symbol
+ window::Union{Int, Nothing}
+ order::Union{Int, Nothing} # For Savitzky-Golay
+
+ function Smoothing(; method=:savitzky_golay, window=nothing, order=nothing)
+ new(method, window, order)
+ end
+end
+
+
+
+"""
+ Normalization(; method=:tic)
+
+A preprocessing step for intensity normalization.
+"""
+struct Normalization <: AbstractPreprocessingStep
+ method::Symbol
+
+ function Normalization(; method=:tic)
+ new(method)
+ end
+end
+
+"""
+ PeakPicking(; method=nothing, ...)
+
+A preprocessing step for peak detection.
+"""
+struct PeakPicking <: AbstractPreprocessingStep
+ method::Union{Symbol, Nothing} # :profile, :wavelet, :centroid
+ snr_threshold::Union{Float64, Nothing}
+ half_window::Union{Int, Nothing}
+ min_peak_prominence::Union{Float64, Nothing}
+ merge_peaks_tolerance::Union{Float64, Nothing}
+ min_peak_width_ppm::Union{Float64, Nothing}
+ max_peak_width_ppm::Union{Float64, Nothing}
+ min_peak_shape_r2::Union{Float64, Nothing}
+
+ function PeakPicking(; method=nothing, snr_threshold=nothing, half_window=nothing, min_peak_prominence=nothing, merge_peaks_tolerance=nothing, min_peak_width_ppm=nothing, max_peak_width_ppm=nothing, min_peak_shape_r2=nothing)
+ new(method, snr_threshold, half_window, min_peak_prominence, merge_peaks_tolerance, min_peak_width_ppm, max_peak_width_ppm, min_peak_shape_r2)
+ end
+end
+
+"""
+ PeakAlignment(; method=:lowess, ...)
+
+A preprocessing step for peak alignment.
+"""
+struct PeakAlignment <: AbstractPreprocessingStep
+ method::Symbol
+ span::Union{Float64, Nothing}
+ tolerance::Union{Float64, Nothing}
+ tolerance_unit::Union{Symbol, Nothing}
+ max_shift_ppm::Union{Float64, Nothing}
+ min_matched_peaks::Union{Int, Nothing}
+
+ function PeakAlignment(; method=:lowess, span=nothing, tolerance=nothing, tolerance_unit=nothing, max_shift_ppm=nothing, min_matched_peaks=nothing)
+ new(method, span, tolerance, tolerance_unit, max_shift_ppm, min_matched_peaks)
+ end
+end
+
+"""
+ PeakSelection(; frequency_threshold=nothing, ...)
+
+A preprocessing step for peak selection (filtering).
+"""
+struct PeakSelection <: AbstractPreprocessingStep
+ frequency_threshold::Union{Float64, Nothing}
+ min_snr::Union{Float64, Nothing}
+ min_fwhm_ppm::Union{Float64, Nothing}
+ max_fwhm_ppm::Union{Float64, Nothing}
+ min_shape_r2::Union{Float64, Nothing}
+ correlation_threshold::Union{Float64, Nothing}
+
+ function PeakSelection(; frequency_threshold=nothing, min_snr=nothing, min_fwhm_ppm=nothing, max_fwhm_ppm=nothing, min_shape_r2=nothing, correlation_threshold=nothing)
+ new(frequency_threshold, min_snr, min_fwhm_ppm, max_fwhm_ppm, min_shape_r2, correlation_threshold)
+ end
+end
+
+"""
+ PeakBinningParams(; method=:adaptive, ...)
+
+A preprocessing step for peak binning.
+"""
+struct PeakBinningParams <: AbstractPreprocessingStep
+ method::Symbol
+ tolerance::Union{Float64, Nothing}
+ tolerance_unit::Union{Symbol, Nothing}
+ frequency_threshold::Union{Float64, Nothing}
+ min_peak_per_bin::Union{Int, Nothing}
+ max_bin_width_ppm::Union{Float64, Nothing}
+ intensity_weighted_centers::Bool
+ num_uniform_bins::Union{Int, Nothing}
+
+ function PeakBinningParams(; method=:adaptive, tolerance=nothing, tolerance_unit=nothing, frequency_threshold=nothing, min_peak_per_bin=nothing, max_bin_width_ppm=nothing, intensity_weighted_centers=true, num_uniform_bins=nothing)
+ new(method, tolerance, tolerance_unit, frequency_threshold, min_peak_per_bin, max_bin_width_ppm, intensity_weighted_centers, num_uniform_bins)
+ end
+end
# =============================================================================
# 0) Quality Control (QC)
@@ -77,6 +226,53 @@ function qc_is_regular(mz::AbstractVector)
return true
end
+"""
+ validate_spectrum(mz, intensity) -> Bool
+
+Performs comprehensive validation on a mass spectrum (m/z and intensity arrays).
+Returns `true` if the spectrum is valid, `false` otherwise.
+
+Checks include:
+- Both `mz` and `intensity` are non-empty.
+- `mz` and `intensity` have the same length.
+- All `m/z` values are finite and non-negative.
+- All `intensity` values are finite and non-negative.
+- `m/z` values are strictly increasing (no duplicates or decreasing values).
+"""
+function validate_spectrum(mz::AbstractVector{<:Real}, intensity::AbstractVector{<:Real})::Bool
+ # 1. Check for empty vectors
+ if isempty(mz) || isempty(intensity)
+ @warn "Spectrum validation failed: m/z or intensity vector is empty."
+ return false
+ end
+
+ # 2. Check for length mismatch
+ if length(mz) != length(intensity)
+ @warn "Spectrum validation failed: m/z and intensity vectors have different lengths ($(length(mz)) vs $(length(intensity)))"
+ return false
+ end
+
+ # 3. Check for NaN/Inf in mz and non-negativity
+ if any(!isfinite, mz) || any(<(0), mz)
+ @warn "Spectrum validation failed: m/z vector contains NaN, Inf, or negative values."
+ return false
+ end
+
+ # 4. Check for NaN/Inf in intensity and non-negativity
+ if any(!isfinite, intensity) || any(<(0), intensity)
+ #@warn "Spectrum validation failed: intensity vector contains NaN, Inf, or negative values."
+ return false
+ end
+
+ # 5. Check for strictly increasing m/z values (and thus no duplicates)
+ if !qc_is_regular(mz)
+ #@warn "Spectrum validation failed: m/z vector is not strictly increasing (contains duplicates or decreasing values)."
+ return false
+ end
+
+ return true
+end
+
# =============================================================================
# 1) Intensity Transformation & Smoothing
# =============================================================================
@@ -103,17 +299,61 @@ function transform_intensity(intensity::AbstractVector{<:Real}; method::Symbol=:
end
"""
- smooth_spectrum(y; window=9, order=2) -> Vector
+ smooth_spectrum(y::AbstractVector{<:Real}; method::Symbol=:savitzky_golay, window::Int=9, order::Int=2) -> Vector
-Applies a Savitzky–Golay filter to smooth the intensity data.
+Applies a smoothing filter to the intensity data.
+
+# Arguments
+- `y`: The intensity data.
+- `method`: The smoothing method (:savitzky_golay or :moving_average).
+- `window`: The window size for the filter.
+- `order`: The polynomial order for Savitzky-Golay
"""
-function smooth_spectrum(y::AbstractVector{<:Real}; window::Int=9, order::Int=2)
- win = isodd(window) ? window : window + 1
- if length(y) < win
- return y # Cannot smooth if data is smaller than window
+function smooth_spectrum(y::AbstractVector{<:Real}; method::Symbol=:savitzky_golay, window::Int=9, order::Int=2)
+ if window < 3
+ throw(ArgumentError("Window size must be at least 3"))
end
- res = SavitzkyGolay.savitzky_golay(collect(float.(y)), win, order)
- return res.y
+ if order >= window
+ throw(ArgumentError("Polynomial order must be less than window size"))
+ end
+ if order < 0
+ throw(ArgumentError("Polynomial order cannot be negative"))
+ end
+ if method === :savitzky_golay
+ win = isodd(window) ? window : window + 1
+ if length(y) < win
+ return y # Cannot smooth if data is smaller than window
+ end
+ res = SavitzkyGolay.savitzky_golay(collect(float.(y)), win, order)
+ return res.y
+ elseif method === :moving_average
+ return moving_average_smooth(y, window)
+ else
+ @warn "Unsupported smoothing method: $method. Returning original intensity."
+ return collect(float.(y))
+ end
+end
+
+"""
+ moving_average_smooth(y::AbstractVector{<:Real}, window::Int) -> Vector
+
+Applies a simple moving average filter to the intensity data.
+"""
+function moving_average_smooth(y::AbstractVector{<:Real}, window::Int)
+ n = length(y)
+ if n < window
+ return collect(float.(y)) # Cannot smooth if data is smaller than window
+ end
+
+ smoothed_y = zeros(Float64, n)
+ half_window = div(window, 2)
+
+ for i in 1:n
+ start_idx = max(1, i - half_window)
+ end_idx = min(n, i + half_window)
+ smoothed_y[i] = mean(@view y[start_idx:end_idx])
+ end
+ return smoothed_y
end
# =============================================================================
@@ -121,24 +361,121 @@ end
# =============================================================================
"""
- snip_baseline(y, iterations=100) -> Vector
+ _snip_baseline_impl(y, iterations=100) -> Vector
-Estimates the baseline of a spectrum using the SNIP algorithm.
+Estimates the baseline of a spectrum using the SNIP algorithm (internal implementation).
"""
-function snip_baseline(y::AbstractVector{<:Real}; iterations::Int=100)
+function _snip_baseline_impl(y::AbstractVector{<:Real}; iterations::Int=100)
n = length(y)
- b = collect(float.(y))
- buf = similar(b)
+
+ # Initialize two buffers. b1 holds the current baseline estimate, b2 for the next.
+ # Always convert to Float64 to ensure type stability and avoid copying if already correct type
+ b1 = collect(float.(y))
+ b2 = similar(b1)
+
+ current_b = b1
+ next_b = b2
+
for k in 1:iterations
- copyto!(buf, b)
- @inbounds for i in 2:n-1
- buf[i] = min(b[i], 0.5 * (b[i-1] + b[i+1]))
+ # Calculate next baseline estimate into `next_b` based on `current_b`
+ # Boundary conditions
+ if n > 1
+ next_b[1] = min(current_b[1], current_b[2])
+ next_b[n] = min(current_b[n], current_b[n-1])
end
- buf[1] = min(b[1], b[2])
- buf[end] = min(b[end], b[end-1])
- b, buf = buf, b
+
+ @inbounds for i in 2:n-1
+ next_b[i] = min(current_b[i], 0.5 * (current_b[i-1] + current_b[i+1]))
+ end
+
+ # Swap references for the next iteration (no data copy here)
+ current_b, next_b = next_b, current_b
+ end
+
+ # Return the final baseline estimate (which is in current_b after the last swap)
+ return current_b
+end
+
+"""
+ convex_hull_baseline(y) -> Vector
+
+Estimates the baseline of a spectrum using the convex hull algorithm.
+"""
+function convex_hull_baseline(y::AbstractVector{<:Real})
+ n = length(y)
+ if n < 3
+ return zeros(Float64, n)
+ end
+
+ # Find upper convex hull points
+ hull_indices = Int[1]
+ for i in 2:n
+ while length(hull_indices) >= 2 &&
+ (y[hull_indices[end]] - y[hull_indices[end-1]]) * (i - hull_indices[end]) <=
+ (y[i] - y[hull_indices[end]]) * (hull_indices[end] - hull_indices[end-1])
+ pop!(hull_indices)
+ end
+ push!(hull_indices, i)
+ end
+
+ # Interpolate between hull points
+ baseline = zeros(Float64, n)
+ for i in 1:(length(hull_indices)-1)
+ idx1 = hull_indices[i]
+ idx2 = hull_indices[i+1]
+
+ # Linear interpolation
+ for j in idx1:idx2
+ baseline[j] = y[idx1] + (y[idx2] - y[idx1]) * (j - idx1) / (idx2 - idx1)
+ end
+ end
+ return baseline
+end
+
+"""
+ median_baseline(y; window=20) -> Vector
+
+Estimates the baseline of a spectrum using a moving median filter.
+"""
+function median_baseline(y::AbstractVector{<:Real}; window::Int=20)
+ n = length(y)
+ if n < window
+ return zeros(Float64, n)
+ end
+
+ baseline = zeros(Float64, n)
+ half_window = div(window, 2)
+
+ for i in 1:n
+ start_idx = max(1, i - half_window)
+ end_idx = min(n, i + half_window)
+ baseline[i] = median(@view y[start_idx:end_idx])
+ end
+ return baseline
+end
+
+"""
+ apply_baseline_correction(y::AbstractVector{<:Real}; method::Symbol=:snip, iterations::Int=100, window::Int=20) -> Vector
+
+Applies a baseline correction algorithm to the intensity data.
+
+# Arguments
+- `y`: The intensity data.
+- `method`: The baseline correction method (:snip, :convex_hull, or :median).
+- `iterations`: Iterations for SNIP method.
+- `window`: Window size for median method.
+"""
+function apply_baseline_correction(y::AbstractVector{<:Real}; method::Symbol=:snip, iterations::Int=100, window::Int=20)
+ if method === :snip
+ return _snip_baseline_impl(y, iterations=iterations)
+ elseif method === :convex_hull
+ return convex_hull_baseline(y)
+ elseif method === :median
+ return median_baseline(y, window=window)
+ else
+ @warn "Unsupported baseline correction method: $method. Returning zero baseline."
+ return zeros(Float64, length(y))
end
- return b
end
# =============================================================================
@@ -184,31 +521,251 @@ function median_normalize(y::AbstractVector{<:Real})
return m <= 0 ? collect(float.(y)) : collect(float.(y)) ./ m
end
+"""
+ rms_normalize(y) -> Vector
+
+Normalizes spectrum intensities to the Root Mean Square (RMS).
+"""
+function rms_normalize(y::AbstractVector{<:Real})
+ s = sqrt(sum(abs2, y) / length(y))
+ return s <= 0 ? collect(float.(y)) : collect(float.(y)) ./ s
+end
+
+"""
+ apply_normalization(y::AbstractVector{<:Real}; method::Symbol=:tic) -> Vector
+
+Applies a per-spectrum normalization algorithm to the intensity data.
+
+# Arguments
+- `y`: The intensity data.
+- `method`: The normalization method (:tic, :median, :rms, or :none).
+"""
+function apply_normalization(y::AbstractVector{<:Real}; method::Symbol=:tic)
+ if method === :tic
+ return tic_normalize(y)
+ elseif method === :median
+ return median_normalize(y)
+ elseif method === :rms
+ return rms_normalize(y)
+ elseif method === :none
+ return collect(float.(y))
+ else
+ @warn "Unsupported normalization method: $method. Returning original intensity."
+ return collect(float.(y))
+ end
+end
+
+"""
+ remove_matrix_peaks_from_spectrum(mz::AbstractVector{<:Real}, intensity::AbstractVector{<:Real},
+ matrix_peak_mzs::AbstractVector{<:Real};
+ tolerance::Float64=0.002, tolerance_unit::Symbol=:mz,
+ removal_method::Symbol=:zero_out) -> Vector{Float64}
+
+Removes or reduces intensity around specified matrix peaks in a single spectrum.
+
+# Arguments
+- `mz`: The m/z vector of the spectrum.
+- `intensity`: The intensity vector of the spectrum.
+- `matrix_peak_mzs`: A list of m/z values identified as matrix peaks.
+- `tolerance`: The m/z tolerance for matching matrix peaks.
+- `tolerance_unit`: Unit of tolerance (:mz or :ppm).
+- `removal_method`: How to remove the peaks (:zero_out or :subtract).
+
+# Returns
+- `Vector{Float64}`: The modified intensity vector.
+"""
+function remove_matrix_peaks_from_spectrum(mz::AbstractVector{<:Real}, intensity::AbstractVector{<:Real},
+ matrix_peak_mzs::AbstractVector{<:Real};
+ tolerance::Float64=0.002, tolerance_unit::Symbol=:mz,
+ removal_method::Symbol=:zero_out)
+ modified_intensity = copy(intensity)
+
+ for matrix_mz in matrix_peak_mzs
+ # Calculate dynamic tolerance if in PPM
+ current_tolerance = (tolerance_unit == :ppm) ? (matrix_mz * tolerance / 1e6) : tolerance
+
+ # Find indices within the tolerance window
+ indices_to_modify = findall(m -> abs(m - matrix_mz) <= current_tolerance, mz)
+
+ if !isempty(indices_to_modify)
+ if removal_method == :zero_out
+ modified_intensity[indices_to_modify] .= 0.0
+ elseif removal_method == :subtract
+ # Baseline-aware subtraction: replace peak with a line connecting its "feet"
+ idx_start = indices_to_modify[1]
+ idx_end = indices_to_modify[end]
+
+ # Ensure we are not at the very edge of the spectrum
+ if idx_start > 1 && idx_end < length(modified_intensity)
+ y1 = modified_intensity[idx_start - 1]
+ y2 = modified_intensity[idx_end + 1]
+ x1 = idx_start - 1
+ x2 = idx_end + 1
+
+ # Linearly interpolate the baseline under the peak
+ for i in idx_start:idx_end
+ # y = y1 + (y2 - y1) * (x - x1) / (x2 - x1)
+ baseline_val = y1 + (y2 - y1) * (i - x1) / (x2 - x1)
+ # Set the intensity to the baseline, but don't increase it (e.g., if baseline is above signal)
+ modified_intensity[i] = min(modified_intensity[i], baseline_val)
+ end
+ else
+ # If peak is at the edge, we can't interpolate, so just zero it out
+ modified_intensity[indices_to_modify] .= 0.0
+ end
+ else
+ @warn "Unsupported matrix peak removal method: $removal_method. Skipping."
+ end
+ end
+ end
+ return modified_intensity
+end
+
+function identify_matrix_peaks_from_blanks(
+ msi_data::MSIData,
+ blank_spectrum_tag::Symbol,
+ snr_threshold::Float64;
+ frequency_threshold::Float64=0.5, # e.g., peak must be in 50% of blanks
+ bin_tolerance::Float64=0.005, # m/z tolerance for binning blank peaks
+ bin_tolerance_unit::Symbol=:mz
+ )::Vector{Float64}
+ if !(0 < frequency_threshold <= 1.0)
+ throw(ArgumentError("`frequency_threshold` must be between 0 and 1 (exclusive of 0). Got: $frequency_threshold"))
+ end
+ println("Identifying matrix peaks from blank spectra (tag: $blank_spectrum_tag)...")
+ blank_indices = Int[]
+ for (i, meta) in enumerate(msi_data.spectra_metadata)
+ if meta.type == blank_spectrum_tag
+ push!(blank_indices, i)
+ end
+ end
+
+ if isempty(blank_indices)
+ @warn "No blank spectra found with tag: $blank_spectrum_tag. Cannot identify matrix peaks."
+ return Float64[]
+ end
+
+ all_blank_peaks = Vector{NamedTuple}[]
+ num_blank_spectra = 0
+
+ # Collect peaks from each blank spectrum
+ _iterate_spectra_fast(msi_data, blank_indices) do idx, mz, intensity
+ num_blank_spectra += 1
+ if !validate_spectrum(mz, intensity)
+ @warn "Blank spectrum $idx is invalid, skipping peak detection for it."
+ push!(all_blank_peaks, NamedTuple[]) # Add empty list to maintain count
+ return
+ end
+ # Assuming profile mode for matrix peak detection
+ peaks = detect_peaks_profile(mz, intensity; snr_threshold=snr_threshold)
+ push!(all_blank_peaks, peaks)
+ end
+
+ if isempty(all_blank_peaks) || all(isempty, all_blank_peaks)
+ @warn "No peaks detected in any blank spectra. Cannot identify matrix peaks."
+ return Float64[]
+ end
+
+ # Flatten all peaks from blank spectra for initial binning
+ flat_blank_peaks = NamedTuple[]
+ for peaks_in_spec in all_blank_peaks
+ append!(flat_blank_peaks, peaks_in_spec)
+ end
+ sort!(flat_blank_peaks, by=p->p.mz)
+
+ # Bin the detected peaks from all blanks to find common m/z features
+ # This is a simplified binning for matrix peak identification
+ binned_matrix_features = Dict{Float64, Int}() # mz_center => count of spectra it appeared in
+
+ i = 1
+ while i <= length(flat_blank_peaks)
+ current_bin_start_idx = i
+ current_peak = flat_blank_peaks[i]
+
+ # Calculate dynamic tolerance if in PPM
+ current_bin_mz = current_peak.mz
+ tol_val = (bin_tolerance_unit == :ppm) ? (current_bin_mz * bin_tolerance / 1e6) : bin_tolerance
+
+ j = i + 1
+ while j <= length(flat_blank_peaks) && (flat_blank_peaks[j].mz - current_peak.mz) <= tol_val
+ j += 1
+ end
+ current_bin_end_idx = j - 1
+
+ # Calculate a representative m/z for the bin (e.g., intensity-weighted average)
+ peaks_in_bin = flat_blank_peaks[current_bin_start_idx:current_bin_end_idx]
+ if !isempty(peaks_in_bin)
+ sum_intensity = sum(p.intensity for p in peaks_in_bin)
+ if sum_intensity > 0
+ bin_center_mz = sum(p.mz * p.intensity for p in peaks_in_bin) / sum(sum_intensity)
+ else
+ bin_center_mz = mean(p.mz for p in peaks_in_bin)
+ end
+
+ # Check how many blank spectra this feature appeared in
+ spectra_count = 0
+ # For each blank spectrum, check if it contains a peak within the current bin's tolerance
+ for spec_peaks in all_blank_peaks
+ if any(p -> abs(p.mz - bin_center_mz) <= tol_val, spec_peaks)
+ spectra_count += 1
+ end
+ end
+ binned_matrix_features[bin_center_mz] = spectra_count
+ end
+ i = j
+ end
+
+ # Filter based on frequency_threshold
+ matrix_peak_mzs = Float64[]
+ min_spectra_count = ceil(Int, num_blank_spectra * frequency_threshold)
+
+ for (mz_center, count) in binned_matrix_features
+ if count >= min_spectra_count
+ push!(matrix_peak_mzs, mz_center)
+ end
+ end
+ sort!(matrix_peak_mzs)
+
+ println("Identified $(length(matrix_peak_mzs)) potential matrix peaks from $(num_blank_spectra) blank spectra.")
+ return matrix_peak_mzs
+end
+
+
# =============================================================================
# 4) Peak Detection
# =============================================================================
"""
- detect_peaks_profile(mz, y; ...)
+ detect_peaks_profile(mz, y; ...) -> Vector{NamedTuple}
-Enhanced peak detection for profile-mode spectra with advanced filtering.
+Enhanced peak detection for profile-mode spectra with advanced filtering and quality metrics.
+
+Returns a vector of NamedTuples, each representing a detected peak with:
+- `mz`: m/z value of the peak
+- `intensity`: Intensity of the peak
+- `fwhm`: Full Width at Half Maximum (Δm)
+- `shape_r2`: Pseudo R^2 for peak shape (simplified)
+- `snr`: Signal-to-Noise Ratio
+- `prominence`: Peak prominence
"""
-function detect_peaks_profile(mz::AbstractVector{<:Real}, y::AbstractVector{<:Real};
+function detect_peaks_profile(mz::AbstractVector{<:Real}, y::AbstractVector{<:Real};
half_window::Int=10,
snr_threshold::Float64=2.0,
- min_peak_width::Real=0.01,
- max_peak_width::Real=2.0,
- peak_shape_threshold::Float64=0.7, # Stub
+ min_peak_prominence::Float64=0.1,
merge_peaks_tolerance::Float64=0.002,
- min_peak_prominence::Float64=0.1)
+ # New parameters for quality filtering (used for calculation, not filtering here)
+ min_peak_width_ppm::Float64=0.0, # Not used for filtering in this function
+ max_peak_width_ppm::Float64=Inf, # Not used for filtering in this function
+ min_peak_shape_r2::Float64=0.0 # Not used for filtering in this function
+ )
n = length(y)
- n < 3 && return (Float64[], Float64[])
+ n < 3 && return NamedTuple{(:mz, :intensity, :fwhm, :shape_r2, :snr, :prominence), Tuple{Float64, Float64, Float64, Float64, Float64, Float64}}[]
noise_level = mad(y, normalize=true) + eps(Float64)
- ys = smooth_spectrum(y; window=max(5, 2*half_window+1), order=2)
+ ys = smooth_spectrum(y; method=:savitzky_golay, window=max(5, 2*half_window+1), order=2) # Use smoothed data for detection
- peak_indices = Int[]
- @inbounds for i in 2:n-1
+ candidate_peak_indices = Int[]
+ for i in 2:n-1
left = max(1, i - half_window)
right = min(n, i + half_window)
@@ -218,67 +775,127 @@ function detect_peaks_profile(mz::AbstractVector{<:Real}, y::AbstractVector{<:Re
if ys[i] >= maximum(@view ys[left:right]) &&
(ys[i] > snr_threshold * noise_level) &&
(prominence > min_peak_prominence * ys[i])
- push!(peak_indices, i)
+ push!(candidate_peak_indices, i)
end
end
- # (STUB) Peak shape validation would be applied here
- # e.g., fitting a Gaussian and checking R^2 > peak_shape_threshold
-
# Merge close peaks
- if !isempty(peak_indices) && merge_peaks_tolerance > 0
- merged_indices = [peak_indices[1]]
- for i in 2:length(peak_indices)
- if (mz[peak_indices[i]] - mz[last(merged_indices)]) > merge_peaks_tolerance
- push!(merged_indices, peak_indices[i])
- elseif y[peak_indices[i]] > y[last(merged_indices)]
- merged_indices[end] = peak_indices[i] # Replace with more intense peak
+ if !isempty(candidate_peak_indices) && merge_peaks_tolerance > 0
+ merged_indices = [candidate_peak_indices[1]]
+ for i in 2:length(candidate_peak_indices)
+ if (mz[candidate_peak_indices[i]] - mz[last(merged_indices)]) > merge_peaks_tolerance
+ push!(merged_indices, candidate_peak_indices[i])
+ elseif y[candidate_peak_indices[i]] > y[last(merged_indices)]
+ merged_indices[end] = candidate_peak_indices[i] # Replace with more intense peak
end
end
- peak_indices = merged_indices
+ candidate_peak_indices = merged_indices
end
- # (STUB) Peak width filtering would be applied here
- # This would require calculating FWHM for each peak, which is computationally intensive.
+ detected_peaks = NamedTuple{(:mz, :intensity, :fwhm, :shape_r2, :snr, :prominence), Tuple{Float64, Float64, Float64, Float64, Float64, Float64}}[]
+ for p_idx in candidate_peak_indices
+ peak_mz = float(mz[p_idx])
+ peak_int = float(y[p_idx])
+
+ fwhm_delta_m = calculate_robust_fwhm(mz, y, p_idx)
+ fwhm_ppm = if isnan(fwhm_delta_m) || fwhm_delta_m <= 0
+ 0.0
+ else
+ 1e6 * fwhm_delta_m / peak_mz
+ end
+ shape_r2 = _fit_gaussian_and_r2(mz, y, p_idx, half_window)
+
+ peak_snr = peak_int / noise_level # Recalculate SNR based on final peak_int
+
+ left = max(1, p_idx - half_window)
+ right = min(n, p_idx + half_window)
+ prominence = ys[p_idx] - max(minimum(@view ys[left:p_idx]), minimum(@view ys[p_idx:right]))
- pk_mz = [float(mz[i]) for i in peak_indices]
- pk_int = [float(y[i]) for i in peak_indices]
+ push!(detected_peaks, (mz=peak_mz, intensity=peak_int, fwhm=fwhm_ppm, shape_r2=shape_r2, snr=peak_snr, prominence=prominence))
+ end
- return (pk_mz, pk_int)
+ return detected_peaks
end
"""
- detect_peaks_wavelet(mz, intensity; ...)
+ detect_peaks_wavelet(mz, intensity; ...) -> Vector{NamedTuple}
Detects peaks using Continuous Wavelet Transform (CWT).
-"""
-function detect_peaks_wavelet(mz::AbstractVector, intensity::AbstractVector; scales=1:10, snr_threshold=3.0)
- n = length(intensity)
- n < 10 && return (Float64[], Float64[])
- cwt_res = ContinuousWavelets.cwt(intensity, ContinuousWavelets.morl)
-
- peak_indices = Int[]
- noise_level = mad(intensity, normalize=true) + eps(Float64)
- for i in 2:n-1
- if intensity[i] > intensity[i-1] && intensity[i] > intensity[i+1] && intensity[i] > snr_threshold * noise_level
- if abs(cwt_res[i, end]) > 0
- push!(peak_indices, i)
+Returns a vector of NamedTuples, each representing a detected peak with:
+- `mz`: m/z value of the peak
+- `intensity`: Intensity of the peak
+- `snr`: Signal-to-Noise Ratio (simplified)
+"""
+function detect_peaks_wavelet(mz::AbstractVector, intensity::AbstractVector; scales=1:10, snr_threshold=3.0, half_window=10)::Vector{NamedTuple{(:mz, :intensity, :fwhm, :shape_r2, :snr, :prominence), Tuple{Float64, Float64, Float64, Float64, Float64, Float64}}}
+ n = length(intensity)
+ n < 10 && return NamedTuple{(:mz, :intensity, :fwhm, :shape_r2, :snr, :prominence), Tuple{Float64, Float64, Float64, Float64, Float64, Float64}}[]
+
+ # Compute CWT
+ cwt_res = ContinuousWavelets.cwt(intensity, ContinuousWavelets.morl) # Morlet is good for peaks
+
+ noise_level = mad(intensity, normalize=true) + eps(Float64)
+
+ candidate_indices = Set{Int}() # Use a Set to store unique peak indices
+
+ # Find local maxima in the CWT coefficients (magnitude)
+ # This assumes `localmaxima` is available (e.g., from ImageFiltering, often a dependency of ContinuousWavelets).
+ abs_cwt_res = abs.(cwt_res)
+ # Iterate over all (m/z index, scale index) pairs that are local maxima in the CWT matrix
+ for (m_idx, scale_idx) in Tuple.(localmaxima(abs_cwt_res))
+ # Ensure m_idx is not at the very edges of the intensity array to avoid index out of bounds
+ if m_idx > 1 && m_idx < n
+ # Check if this CWT maximum corresponds to a local maximum in the original intensity
+ # AND if the original intensity is above the SNR threshold
+ if intensity[m_idx] > intensity[m_idx-1] && intensity[m_idx] > intensity[m_idx+1] && intensity[m_idx] > snr_threshold * noise_level
+ push!(candidate_indices, m_idx)
end
end
end
+ peak_indices = collect(candidate_indices)
+ sort!(peak_indices) # Ensure order
+
+ detected_peaks = NamedTuple{(:mz, :intensity, :fwhm, :shape_r2, :snr, :prominence), Tuple{Float64, Float64, Float64, Float64, Float64, Float64}}[]
+ for p_idx in peak_indices
+ peak_mz = float(mz[p_idx])
+ peak_int = float(intensity[p_idx])
+ peak_snr = peak_int / noise_level
+
+ fwhm_delta_m = calculate_robust_fwhm(mz, intensity, p_idx)
+ fwhm_ppm = if isnan(fwhm_delta_m) || fwhm_delta_m <= 0
+ 0.0
+ else
+ 1e6 * fwhm_delta_m / peak_mz
+ end
+ shape_r2 = _fit_gaussian_and_r2(mz, intensity, p_idx, half_window)
+
+ push!(detected_peaks, (mz=peak_mz, intensity=peak_int, fwhm=fwhm_ppm, shape_r2=shape_r2, snr=peak_snr, prominence=peak_int))
+ end
- return (mz[peak_indices], intensity[peak_indices])
+ return detected_peaks
end
"""
- detect_peaks_centroid(mz, y; ...)
+ detect_peaks_centroid(mz, y; ...) -> Vector{NamedTuple}
-Filters peaks in centroid-mode data.
+Filters peaks in centroid-mode data based on intensity threshold.
+
+Returns a vector of NamedTuples, each representing a detected peak with:
+- `mz`: m/z value of the peak
+- `intensity`: Intensity of the peak
"""
-function detect_peaks_centroid(mz::AbstractVector{<:Real}, y::AbstractVector{<:Real}; intensity_threshold::Float64=0.0)
- keep_indices = findall(y .>= intensity_threshold)
- return (mz[keep_indices], y[keep_indices])
+function detect_peaks_centroid(mz::AbstractVector{<:Real}, y::AbstractVector{<:Real}; snr_threshold::Float64=0.0)
+ noise_level = mad(y, normalize=true) + eps(Float64)
+ detected_peaks = NamedTuple{(:mz, :intensity, :fwhm, :shape_r2, :snr, :prominence), Tuple{Float64, Float64, Float64, Float64, Float64, Float64}}[]
+
+ for i in eachindex(mz)
+ snr = noise_level > 0 ? y[i] / noise_level : y[i] > 0 ? Inf : 0.0
+ if snr >= snr_threshold
+ # For centroid data, FWHM and shape are not applicable. Return placeholders.
+ push!(detected_peaks, (mz=float(mz[i]), intensity=float(y[i]), fwhm=0.0, shape_r2=1.0, snr=snr, prominence=y[i]))
+ end
+ end
+ return detected_peaks
end
# =============================================================================
@@ -291,6 +908,8 @@ end
Enhanced peak alignment with PPM tolerance and other constraints.
"""
function align_peaks_lowess(ref_mz::Vector{<:Real}, tgt_mz::Vector{<:Real};
+ method::Symbol=:linear, # :linear, :lowess, or :ransac
+ span::Float64=0.75, # Span for LOWESS
tolerance::Float64=0.002,
tolerance_unit::Symbol=:mz,
max_shift_ppm::Float64=50.0,
@@ -298,11 +917,12 @@ function align_peaks_lowess(ref_mz::Vector{<:Real}, tgt_mz::Vector{<:Real};
pairs = Tuple{Float64,Float64}[]
i = 1; j = 1
while i <= length(tgt_mz) && j <= length(ref_mz)
+ # Dynamic tolerance for PPM
tol = (tolerance_unit == :ppm) ? (ref_mz[j] * tolerance / 1e6) : tolerance
dt = tgt_mz[i] - ref_mz[j]
if abs(dt) <= tol
- # Max shift check
+ # Max shift check to avoid spurious matches
if abs(dt) * 1e6 / ref_mz[j] <= max_shift_ppm
push!(pairs, (float(tgt_mz[i]), float(ref_mz[j])))
end
@@ -315,16 +935,98 @@ function align_peaks_lowess(ref_mz::Vector{<:Real}, tgt_mz::Vector{<:Real};
end
if length(pairs) < min_matched_peaks
- @warn "Too few matching peaks ($(length(pairs)) < $min_matched_peaks). Returning identity function."
+ #@warn "Too few matching peaks ($(length(pairs)) < $min_matched_peaks). Returning identity function."
return x -> float.(x)
end
- t = [p[1] for p in pairs]
- r = [p[2] for p in pairs]
+ t = [p[1] for p in pairs] # Target m/z (x-axis)
+ r = [p[2] for p in pairs] # Reference m/z (y-axis)
- # (STUB) A robust regression model (e.g., RANSAC) would be better here.
- itp = linear_interpolation(t, r, extrapolation_bc=Line())
- return itp
+ if method == :linear
+ itp = linear_interpolation(t, r, extrapolation_bc=Line())
+ return itp
+ elseif method == :lowess
+ try
+ model = loess(t, r; span=span)
+ # Predict over the original target m/z values (t)
+ predicted_r = Loess.predict(model, t)
+ # Create a linear interpolation from original t and Loess-predicted r,
+ # allowing linear extrapolation.
+ itp = linear_interpolation(t, predicted_r, extrapolation_bc=Line())
+ return itp
+ catch e
+ @warn "LOWESS fitting failed: $e. Falling back to linear interpolation."
+ itp = linear_interpolation(t, r, extrapolation_bc=Line())
+ return itp
+ end
+ elseif method == :ransac
+ # RANSAC implementation for linear model
+ num_iterations = 100
+ best_model_inliers_count = -1
+ best_inliers_t = Float64[]
+ best_inliers_r = Float64[]
+
+ n_pairs = length(pairs)
+ if n_pairs < 2 # Need at least 2 points to fit a line
+ @warn "Not enough matched peaks ($n_pairs) for RANSAC. Falling back to linear interpolation."
+ itp = linear_interpolation(t, r, extrapolation_bc=Line())
+ return itp
+ end
+
+ for iter in 1:num_iterations
+ # 1. Randomly select 2 points
+ sample_indices = StatsBase.sample(1:n_pairs, 2, replace=false)
+ p1_x, p1_y = pairs[sample_indices[1]]
+ p2_x, p2_y = pairs[sample_indices[2]]
+
+ # Avoid vertical line for simplicity in this basic implementation
+ if isapprox(p1_x, p2_x, atol=1e-9)
+ continue
+ end
+
+ # 2. Fit a line y = mx + b
+ m = (p2_y - p1_y) / (p2_x - p1_x)
+ b = p1_y - m * p1_x
+
+ current_inliers_t = Float64[]
+ current_inliers_r = Float64[]
+ current_inliers_count = 0
+
+ # 3. Find inliers
+ for (px, py) in pairs
+ predicted_y = m * px + b
+ residual = abs(py - predicted_y)
+
+ # Determine inlier threshold based on tolerance_unit
+ current_inlier_threshold = (tolerance_unit == :ppm) ? (px * tolerance / 1e6) : tolerance
+
+ if residual <= current_inlier_threshold
+ current_inliers_count += 1
+ push!(current_inliers_t, px)
+ push!(current_inliers_r, py)
+ end
+ end
+
+ # 4. Evaluate model
+ if current_inliers_count > best_model_inliers_count
+ best_model_inliers_count = current_inliers_count
+ best_inliers_t = current_inliers_t
+ best_inliers_r = current_inliers_r
+ end
+ end
+
+ if best_model_inliers_count >= 2 # At least 2 inliers to form a line
+ return linear_interpolation(best_inliers_t, best_inliers_r, extrapolation_bc=Line())
+ else
+ @warn "RANSAC failed to find a robust model (found $(best_model_inliers_count) inliers). Falling back to linear interpolation."
+ itp = linear_interpolation(t, r, extrapolation_bc=Line())
+ return itp
+ end
+ else
+ @warn "Unsupported alignment method: $method. Falling back to linear interpolation."
+ itp = linear_interpolation(t, r, extrapolation_bc=Line())
+ return itp
+ end
end
"""
@@ -334,14 +1036,21 @@ Finds peaks that match a list of reference masses.
"""
function find_calibration_peaks(mz::AbstractVector, intensity::AbstractVector, reference_masses::AbstractVector; ppm_tolerance=20.0)
matched_peaks = Dict{Float64, Float64}()
- detected_mz, _ = detect_peaks_profile(mz, intensity)
+
+ # detect_peaks_profile returns Vector{NamedTuple}, so we need to extract mz values
+ detected_peaks_list = detect_peaks_profile(mz, intensity)
+
+ # Extract only the m/z values into a new vector for easier processing
+ detected_mz_values = [p.mz for p in detected_peaks_list]
for ref_mass in reference_masses
tol = ref_mass * ppm_tolerance / 1e6
- candidates = findall(m -> abs(m - ref_mass) <= tol, detected_mz)
+ # Find candidates within the detected m/z values
+ candidates = findall(m -> abs(m - ref_mass) <= tol, detected_mz_values)
if !isempty(candidates)
- closest_peak_idx = argmin(abs.(detected_mz[candidates] .- ref_mass))
- matched_peaks[ref_mass] = detected_mz[candidates[closest_peak_idx]]
+ # Find the closest detected peak to the reference mass among candidates
+ closest_peak_idx = argmin(abs.(detected_mz_values[candidates] .- ref_mass))
+ matched_peaks[ref_mass] = detected_mz_values[candidates[closest_peak_idx]]
end
end
return matched_peaks
@@ -383,580 +1092,182 @@ end
# =============================================================================
"""
- bin_peaks(all_pk_mz, all_pk_int, tolerance; ...)
+ bin_peaks(all_pk_mz::Vector{Vector{Float64}},
+ all_pk_int::Vector{Vector{Float64}},
+ params::PeakBinningParams) -> Tuple{FeatureMatrix, Vector{Tuple{Float64,Float64}}}
-Enhanced peak binning with adaptive and PPM-based parameters.
+Enhanced peak binning with adaptive and PPM-based parameters, or uniform binning.
+
+# Arguments
+- `all_pk_mz`: A vector of m/z vectors for all spectra.
+- `all_pk_int`: A vector of intensity vectors for all spectra.
+- `params`: A `PeakBinningParams` struct.
+
+# Returns
+- `Tuple{FeatureMatrix, Vector{Tuple{Float64,Float64}}}`: A tuple containing the generated FeatureMatrix and the bin definitions.
+
+# Thread Safety
+The use of `Threads.@threads` in the `:adaptive` and `:uniform` methods is safe.
+- In the `:adaptive` method, the loop is over the bins (`j` index). Each thread writes only to its assigned column `X[:, j]`, so there are no write conflicts between threads.
+- In the `:uniform` method, the loop is over the spectra (`s_idx`). Writes to `X[s_idx, bin_idx]` could theoretically conflict if different peaks from the same spectrum (`s_idx`) are processed by different threads. However, the loop is over `s_idx`, meaning each thread handles a distinct spectrum, making writes to `X[s_idx, :]` exclusive to that thread and thus safe.
"""
-function bin_peaks(all_pk_mz::Vector{<:AbstractVector{<:Real}},
- all_pk_int::Vector{<:AbstractVector{<:Real}},
- tolerance::Float64;
- frequency_threshold::Float64=0.25,
- tolerance_unit::Symbol=:mz,
- adaptive_tolerance::Bool=false, # Stub
- min_peak_per_bin::Int=2,
- max_bin_width_ppm::Float64=100.0,
- intensity_weighted_centers::Bool=true)
- ns = length(all_pk_mz)
- ns == 0 && return (zeros(0,0), Tuple{Float64,Float64}[])
+function bin_peaks(spectra::Vector{MutableSpectrum},
+ params::PeakBinningParams)
+ ns = length(spectra) # Number of spectra
+ ns == 0 && return FeatureMatrix(zeros(0,0), Tuple{Float64,Float64}[], Int[]), Tuple{Float64,Float64}[]
- # Collect all peaks with their intensities and original spectrum index
- all_peaks = [(float(m), float(i), s_idx) for s_idx in 1:ns for (m, i) in zip(all_pk_mz[s_idx], all_pk_int[s_idx])]
- sort!(all_peaks, by=p->p[1])
-
- isempty(all_peaks) && return (zeros(ns, 0), Tuple{Float64,Float64}[])
-
- # Create bins
- bins = []
- current_bin_peaks = [all_peaks[1]]
- for p in @view all_peaks[2:end]
- bin_center = mean(first.(current_bin_peaks))
- tol = (tolerance_unit == :ppm) ? (bin_center * tolerance / 1e6) : tolerance
-
- if p[1] - last(current_bin_peaks)[1] <= tol
- push!(current_bin_peaks, p)
- else
- if length(current_bin_peaks) >= min_peak_per_bin
- push!(bins, current_bin_peaks)
+ if params.method == :uniform
+ println("Performing uniform binning with $(params.num_uniform_bins) bins.")
+ # Determine global m/z range from all peaks
+ min_mz = Inf
+ max_mz = -Inf
+ for s in spectra
+ if !isempty(s.peaks)
+ # Extract m/z values from peaks NamedTuple
+ pk_mz_vec = [p.mz for p in s.peaks]
+ min_mz = min(min_mz, minimum(pk_mz_vec))
+ max_mz = max(max_mz, maximum(pk_mz_vec))
end
- current_bin_peaks = [p]
end
- end
- length(current_bin_peaks) >= min_peak_per_bin && push!(bins, current_bin_peaks)
- # Filter bins by max width
- filter!(b -> (last(b)[1] - first(b)[1]) * 1e6 / mean(p[1] for p in b) <= max_bin_width_ppm, bins)
+ if !isfinite(min_mz) || !isfinite(max_mz) || min_mz == max_mz
+ @warn "Could not determine valid m/z range for uniform binning. Returning empty FeatureMatrix."
+ return FeatureMatrix(zeros(0,0), Tuple{Float64,Float64}[], Int[]), Tuple{Float64,Float64}[]
+ end
- # Create feature matrix
- X = zeros(Float64, ns, length(bins))
- final_bins_boundaries = Vector{Tuple{Float64,Float64}}(undef, length(bins))
-
- for (j, bin_peaks) in enumerate(bins)
- # Calculate bin center
- local bin_center
- if intensity_weighted_centers
- weights = [p[2] for p in bin_peaks]
- bin_center = sum(p[1]*p[2] for p in bin_peaks) / sum(weights)
- else
- bin_center = mean(p[1] for p in bin_peaks)
+ # Create uniform bins
+ mz_edges = collect(range(min_mz, stop=max_mz, length=params.num_uniform_bins + 1))
+ bin_definitions = Vector{Tuple{Float64,Float64}}(undef, params.num_uniform_bins)
+ for i in 1:params.num_uniform_bins
+ bin_definitions[i] = (mz_edges[i], mz_edges[i+1])
end
- final_bins_boundaries[j] = (first(bin_peaks)[1], last(bin_peaks)[1])
-
- for p in bin_peaks
- s_idx = p[3]
- X[s_idx, j] = max(X[s_idx, j], p[2])
- end
- end
-
- # Filter by frequency
- if frequency_threshold > 0
- present_count = vec(sum(X .> 0, dims=1))
- min_count = ceil(Int, frequency_threshold * ns)
- keep_mask = findall(present_count .>= min_count)
- X = X[:, keep_mask]
- final_bins_boundaries = final_bins_boundaries[keep_mask]
- end
-
- return (X, final_bins_boundaries)
-end
-
-# =============================================================================
-# 7) Spatial & Advanced Processing (Stubs & New Functions)
-# =============================================================================
-
-"""
- find_ppm_error_by_region(msi_data, region_masks, reference_peaks) -> Dict
-
-Calculates PPM error statistics for different spatial regions.
-`region_masks` is a Dict mapping region names (e.g., :tumor) to BitMatrix masks.
-"""
-function find_ppm_error_by_region(msi_data::MSIData, region_masks::Dict, reference_peaks::Dict)
- regional_reports = Dict{Symbol, NamedTuple}()
-
- width, height = msi_data.image_dims
-
- for (region_name, mask) in region_masks
- mask_height, mask_width = size(mask)
- if mask_width != width || mask_height != height
- @warn "Mask dimensions ($(mask_width)x$(mask_height)) for region '$region_name' do not match image dimensions ($(width)x$(height)). Skipping."
- continue
- end
-
- indices = [
- i for i in 1:length(msi_data.spectra_metadata)
- if msi_data.spectra_metadata[i].x > 0 &&
- msi_data.spectra_metadata[i].y > 0 &&
- mask[msi_data.spectra_metadata[i].y, msi_data.spectra_metadata[i].x]
- ]
+ X = zeros(Float64, ns, params.num_uniform_bins)
- if isempty(indices) continue end
-
- # Call analyze_mass_accuracy with the specific indices for the region
- regional_reports[region_name] = analyze_mass_accuracy(msi_data, reference_peaks; spectrum_indices=indices)
- end
- return regional_reports
-end
-
-"""
- regional_calibration(msi_data, region_masks, reference_peaks)
-
-(STUB) Applies different calibration models to different spatial regions.
-"""
-function regional_calibration(msi_data::MSIData, region_masks::Dict, reference_peaks::Dict)
- @warn "regional_calibration is a stub and not fully implemented."
- # 1. For each region in region_masks:
- # 2. Create a region-specific calibration model using `calibrate_spectra`.
- # 3. Apply the model to all spectra within that region.
- # 4. Return a new MSIData object or modified spectra vector.
- return msi_data # Return unmodified for now
-end
-
-# =============================================================================
-# 8) Advanced Peak Quality & Adaptive Parameters
-# =============================================================================
-
-"""
- calculate_peak_quality_metrics(peak_mz, peak_intensity, local_spectrum) -> Dict
-
-(STUB) Calculates advanced quality metrics for a single peak.
-"""
-function calculate_peak_quality_metrics(peak_mz, peak_intensity, local_spectrum)
- # A full implementation would calculate:
- # - Sharpness (e.g., ratio of height to FWHM)
- # - Symmetry (e.g., ratio of left/right half-widths)
- # - Signal-to-Noise (using local noise estimation)
- # - Isolation score (how close are other peaks)
- return Dict(:sharpness => 1.0, :symmetry => 1.0, :snr => 10.0)
-end
-
-
-
-
-"""
- calculate_adaptive_bin_tolerance(ppm_error_distribution) -> Float64
-
-Calculates an appropriate binning tolerance based on observed mass accuracy.
-"""
-function calculate_adaptive_bin_tolerance(ppm_error_distribution::Vector{Float64})
- if isempty(ppm_error_distribution)
- return 20.0 # Default if no data
- end
- # A robust strategy: mean + 3 * std deviation to capture ~99.7% of peaks
- return mean(ppm_error_distribution) + 3 * std(ppm_error_distribution)
-end
-
-
-# =============================================================================
-# 7) Plotting Helper
-# =============================================================================
-
-"""
- plot_stage_spectrum(mz, intensity; title, xlabel, ylabel) -> Figure
-
-A helper function to generate a plot of a single spectrum using `CairoMakie`. This is
-useful for visualizing the output of different preprocessing steps.
-
-# Arguments
-- `mz::AbstractVector`: The m/z vector for the x-axis.
-- `intensity::AbstractVector`: The intensity vector for the y-axis.
-
-# Keyword Arguments
-- `title::AbstractString`: The title of the plot.
-- `xlabel::AbstractString`: The label for the x-axis. Defaults to "m/z".
-- `ylabel::AbstractString`: The label for the y-axis. Defaults to "Intensity".
-
-# Returns
-- `CairoMakie.Figure`: A `Figure` object containing the plot. The caller is responsible for displaying or saving it.
-
-# Example
-```julia
-using CairoMakie
-mz = 100:0.1:110;
-intensity = rand(length(mz));
-fig = plot_stage_spectrum(mz, intensity, title="My Spectrum");
-save("my_spectrum.png", fig);
-```
-"""
-function plot_stage_spectrum(mz::AbstractVector, intensity::AbstractVector;
- title::AbstractString, xlabel::AbstractString="m/z",
- ylabel::AbstractString="Intensity")
- # This dynamic import is for script-like use; in a package, Makie would be a full dependency.
- @eval begin
- import CairoMakie
- using CairoMakie
- end
- fig = CairoMakie.Figure(size = (1400, 500))
- ax = CairoMakie.Axis(fig[1, 1], title=title, xlabel=xlabel, ylabel=ylabel)
- CairoMakie.lines!(ax, mz, intensity)
- return fig
-end
-
-
-# =============================================================================
-# 9) Quality Control Metrics
-# =============================================================================
-
-"""
- get_common_calibration_standards(type::Symbol) -> Dict{Float64, String}
-
-Returns a dictionary of common m/z values for calibration standards.
-
-# Arguments
-- `type::Symbol`: The type of standards to return. Currently supports `:maldi_pos`.
-
-# Returns
-- `Dict{Float64, String}`: A dictionary mapping theoretical m/z to compound names.
-"""
-function get_common_calibration_standards(type::Symbol)
- if type == :maldi_pos
- # Common peptide/protein standards for positive mode MALDI
- return Dict{Float64, String}(
- 1046.5420 => "Angiotensin II",
- 1060.5690 => "Bradykinin",
- 1296.6853 => "Angiotensin I",
- 2465.1989 => "ACTH clip (1-24)"
- )
- else
- @warn "Unsupported calibration standard type: $type. Returning empty dictionary."
- return Dict{Float64, String}()
- end
-end
-
-"""
- calculate_ppm_error(measured_mz::Real, theoretical_mz::Real) -> Float64
-
-Calculates the mass accuracy error in parts-per-million (PPM) between a measured
-and a theoretical m/z value.
-
-# Arguments
-- `measured_mz::Real`: The experimentally measured m/z value.
-- `theoretical_mz::Real`: The known, theoretical m/z value of a compound.
-
-# Returns
-- `Float64`: The calculated PPM error. Returns `Inf` if `theoretical_mz` is zero.
-
-# Formula
-`PPM = 10^6 * |measured_mz - theoretical_mz| / theoretical_mz`
-
-# Example
-```julia
-calculate_ppm_error(100.005, 100.0) # returns 50.0
-```
-"""
-function calculate_ppm_error(measured_mz::Real, theoretical_mz::Real)
- if theoretical_mz == 0
- return Inf
- end
- return 1e6 * abs(Float64(measured_mz) - Float64(theoretical_mz)) / Float64(theoretical_mz)
-end
-
-"""
- calculate_ppm_error_bulk(measured_mz::Vector{<:Real}, theoretical_mz::Vector{<:Real}) -> Vector{Float64}
-
-Calculates PPM errors for multiple pairs of measured and theoretical mass values.
-
-# Arguments
-- `measured_mz::Vector{<:Real}`: A vector of experimentally measured m/z values.
-- `theoretical_mz::Vector{<:Real}`: A vector of known, theoretical m/z values.
-
-# Returns
-- `Vector{Float64}`: A vector containing the calculated PPM error for each pair.
-"""
-function calculate_ppm_error_bulk(measured_mz::Vector{Real}, theoretical_mz::Vector{Real})
- return [calculate_ppm_error(m, t) for (m, t) in zip(measured_mz, theoretical_mz)]
-end
-
-"""
- calculate_resolution_fwhm(mz::Real, profile_mz::AbstractVector{<:Real}, profile_intensity::AbstractVector{<:Real}) -> Float64
-
-Calculates the mass resolution of a peak in profile-mode data using the Full Width at
-Half Maximum (FWHM) method. Resolution is a measure of an instrument's ability to
-distinguish between two peaks of slightly different mass-to-charge ratios.
-
-# Arguments
-- `mz::Real`: The m/z value of the peak's centroid.
-- `profile_mz::AbstractVector{<:Real}`: The full m/z array from the profile-mode spectrum.
-- `profile_intensity::AbstractVector{<:Real}`: The full intensity array from the profile-mode spectrum.
-
-# Returns
-- `Float64`: The calculated resolution (`m / Δm`). Returns `NaN` if the FWHM cannot be determined (e.g., peak is at the edge of the spectrum).
-
-# Formula
-`Resolution = m / Δm`, where `Δm` is the FWHM.
-"""
-function calculate_resolution_fwhm(mz::Real, profile_mz::AbstractVector{<:Real},
- profile_intensity::AbstractVector{<:Real})
-
- # Find peak center index
- peak_idx = argmin(abs.(profile_mz .- mz))
- peak_height = Float64(profile_intensity[peak_idx])
- half_max = peak_height / 2
-
- # Find left half-maximum point (interpolate for accuracy)
- left_idx = find_last_below(profile_intensity[1:peak_idx], half_max)
- if left_idx == 0 || left_idx == length(profile_intensity[1:peak_idx])
- return NaN
- end
-
- # Linear interpolation for left FWHM
- x1, x2 = Float64(profile_mz[left_idx]), Float64(profile_mz[left_idx+1])
- y1, y2 = Float64(profile_intensity[left_idx]), Float64(profile_intensity[left_idx+1])
- left_fwhm = x1 + (x2 - x1) * (half_max - y1) / (y2 - y1)
-
- # Find right half-maximum point
- right_slice = profile_intensity[peak_idx:end]
- right_offset = find_first_below(right_slice, half_max)
- if right_offset == 0 || right_offset == length(right_slice)
- return NaN
- end
-
- right_idx = peak_idx + right_offset - 1
- x1, x2 = Float64(profile_mz[right_idx-1]), Float64(profile_mz[right_idx])
- y1, y2 = Float64(profile_intensity[right_idx-1]), Float64(profile_intensity[right_idx])
- right_fwhm = x1 + (x2 - x1) * (half_max - y1) / (y2 - y1)
-
- fwhm = right_fwhm - left_fwhm
- return fwhm > 0 ? Float64(mz) / fwhm : NaN
-end
-
-# Helper functions for FWHM calculation
-function find_last_below(v::AbstractVector{<:Real}, threshold::Real)
- for i in length(v):-1:2
- if v[i] >= threshold && v[i-1] < threshold
- return i-1
- end
- end
- return 0
-end
-
-function find_first_below(v::AbstractVector{<:Real}, threshold::Real)
- for i in 1:(length(v)-1)
- if v[i] >= threshold && v[i+1] < threshold
- return i+1
- end
- end
- return 0
-end
-
-"""
- analyze_mass_accuracy(msi_data, reference_peaks; ppm_tolerance=5.0, sample_spectra=100) -> NamedTuple
-
-Analyzes the mass accuracy across a sample of spectra from an MSI dataset by comparing
-detected peaks against a list of known reference m/z values. It calculates PPM error
-statistics and suggests an optimal PPM tolerance for future analyses.
-
-# Arguments
-- `msi_data`: An `MSIData` object containing the dataset.
-- `reference_peaks::Dict{Float64,String}`: A dictionary mapping theoretical m/z values to compound names.
-
-# Keyword Arguments
-- `ppm_tolerance::Float64`: The initial PPM tolerance used to match detected peaks to reference peaks. Defaults to 5.0.
-- `sample_spectra::Int`: The number of spectra to sample from the dataset for the analysis. Defaults to 100.
-
-# Returns
-- `NamedTuple`: A comprehensive report containing statistics like `mean_ppm`, `std_ppm`, `optimal_ppm`, `n_matches`, and lists of matched peaks and all PPM errors.
-"""
-function analyze_mass_accuracy(msi_data, reference_peaks::Dict{Float64,String};
- ppm_tolerance::Float64=5.0, sample_spectra=100,
- spectrum_indices=nothing)
- println("\n[ MASS ACCURACY ANALYSIS ]")
- println("PPM tolerance: ", ppm_tolerance, " ppm")
-
- theoretical_mz = sort(collect(keys(reference_peaks)))
- ppm_errors = Float64[]
- matched_peaks = Tuple{Float64,Float64,String}[] # (theoretical, measured, compound)
-
- # Determine which spectra to process
- local indices_to_process
- if spectrum_indices === nothing
- println("Spectra to sample: ", sample_spectra)
- if sample_spectra >= length(msi_data.spectra_metadata)
- indices_to_process = 1:length(msi_data.spectra_metadata)
- else
- indices_to_process = round.(Int, range(1, length(msi_data.spectra_metadata), length=sample_spectra))
- end
- else
- indices_to_process = spectrum_indices
- println("Analyzing $(length(indices_to_process)) specified spectra.")
- end
-
- for idx in indices_to_process
- mz, intensity = GetSpectrum(msi_data, idx)
-
- # Detect peaks in this spectrum
- detected_peaks, _ = detect_peaks_profile(mz, intensity, snr_threshold=3.0)
-
- # Match detected peaks to reference peaks
- for (i, theoretical) in enumerate(theoretical_mz)
- # Find closest detected peak within tolerance
- distances = abs.(detected_peaks .- theoretical)
- if !isempty(distances)
- min_idx = argmin(distances)
- min_distance = distances[min_idx]
-
- ppm_error = calculate_ppm_error(detected_peaks[min_idx], theoretical)
-
- if ppm_error <= ppm_tolerance
- push!(ppm_errors, ppm_error)
- push!(matched_peaks, (theoretical, detected_peaks[min_idx], reference_peaks[theoretical]))
+ Threads.@threads for s_idx in 1:ns
+ s = spectra[s_idx] # Get the current MutableSpectrum
+ for p in s.peaks # Iterate over peaks directly
+ m = p.mz
+ i = p.intensity
+ # Find which bin this peak belongs to
+ bin_idx = searchsortedlast(mz_edges, m)
+ if bin_idx > 0 && bin_idx <= params.num_uniform_bins
+ X[s_idx, bin_idx] = max(X[s_idx, bin_idx], i) # Take max intensity in bin
end
end
end
- end
-
- if isempty(ppm_errors)
- @warn "No peaks matched within $ppm_tolerance ppm tolerance"
- return (mean_ppm=NaN, std_ppm=NaN, min_ppm=NaN, max_ppm=NaN, optimal_ppm=NaN, n_matches=0, matched_peaks=[], all_ppm_errors=Float64[])
- end
-
- # Calculate statistics
- mean_ppm = mean(ppm_errors)
- std_ppm = std(ppm_errors)
- min_ppm = minimum(ppm_errors)
- max_ppm = maximum(ppm_errors)
-
- # Determine optimal ppm tolerance (mean + 3σ covers ~99.7% of peaks for normal distribution)
- optimal_ppm = mean_ppm + 3 * std_ppm
-
- return (
- mean_ppm = mean_ppm,
- std_ppm = std_ppm,
- min_ppm = min_ppm,
- max_ppm = max_ppm,
- optimal_ppm = optimal_ppm,
- n_matches = length(ppm_errors),
- matched_peaks = matched_peaks,
- all_ppm_errors = ppm_errors
- )
-end
-"""
- generate_qc_report(msi_data, filename; reference_peaks, output_dir, sample_spectra) -> Tuple
+ if params.frequency_threshold !== nothing && params.frequency_threshold > 0
+ present_count = vec(sum(X .> 0, dims=1))
+ min_count = ceil(Int, params.frequency_threshold * ns)
+ keep_mask = findall(present_count .>= min_count)
+ X = X[:, keep_mask]
+ bin_definitions = bin_definitions[keep_mask]
+ end
-Generates a comprehensive Quality Control (QC) report for an MSI dataset. It assesses
-mass accuracy and resolution, saving the results to CSV files and a summary text file.
+ return FeatureMatrix(X, bin_definitions, collect(1:ns)), bin_definitions
-# Arguments
-- `msi_data`: The `MSIData` object to be analyzed.
-- `filename::String`: The name of the input data file, used for reporting.
+ elseif params.method == :adaptive
+ println("Performing adaptive binning with tolerance $(params.tolerance) $(params.tolerance_unit).")
+ # Collect all peaks with their intensities and original spectrum ID directly from MutableSpectrum objects
+ all_peaks = Vector{Tuple{Float64, Float64, Int}}()
+ sizehint!(all_peaks, sum(length(s.peaks) for s in spectra)) # Pre-allocate memory
-# Keyword Arguments
-- `reference_peaks::Dict`: A dictionary of reference peaks for mass accuracy analysis. If `nothing`, defaults are loaded using `get_common_calibration_standards`.
-- `output_dir::String`: The directory where the report files will be saved. Defaults to `"qc_results"`.
-- `sample_spectra::Int`: The number of spectra to sample for the analysis. Defaults to 100.
-
-# Returns
-- `Tuple`: A tuple containing the detailed `accuracy_report` (NamedTuple) and `resolution_results` (Vector).
-
-# Side Effects
-- Creates a directory specified by `output_dir`.
-- Writes `mass_accuracy_results.csv`, `resolution_results.csv`, and `qc_summary.txt` into the output directory.
-"""
-function generate_qc_report(msi_data, filename::String; reference_peaks=nothing, output_dir="qc_results", sample_spectra=100, spectrum_indices=nothing)
- println("\n[ QC REPORT GENERATION ]")
- println("Input file: ", filename)
- println("Output directory: ", output_dir)
- mkpath(output_dir)
-
- # Use default calibrants if none provided
- if reference_peaks === nothing
- reference_peaks = get_common_calibration_standards(:maldi_pos)
- end
-
- println("Generating QC Report...")
- println("Using $(length(reference_peaks)) reference masses")
-
- # 1. Analyze mass accuracy
- accuracy_report = analyze_mass_accuracy(msi_data, reference_peaks; sample_spectra=sample_spectra, spectrum_indices=spectrum_indices)
-
- println("\n" * "="^50)
- println("MASS ACCURACY REPORT")
- println("="^50)
- println("Mean PPM error: $(round(accuracy_report.mean_ppm, digits=2)) ppm")
- println("Std PPM error: $(round(accuracy_report.std_ppm, digits=2)) ppm")
- println("Min PPM error: $(round(accuracy_report.min_ppm, digits=2)) ppm")
- println("Max PPM error: $(round(accuracy_report.max_ppm, digits=2)) ppm")
- if haskey(accuracy_report, :optimal_ppm)
- println("Optimal PPM tolerance: $(round(accuracy_report.optimal_ppm, digits=2)) ppm")
- else
- println("Optimal PPM tolerance: Not available")
- end
- println("Number of matches: $(accuracy_report.n_matches)")
-
- # 2. Calculate resolution for a few representative peaks
- println("\n" * "="^50)
- println("RESOLUTION ANALYSIS")
- println("="^50)
-
- resolution_results = []
- sample_spectra = min(10, length(msi_data.spectra_metadata))
-
- for (i, idx) in enumerate(round.(Int, range(1, length(msi_data.spectra_metadata), length=sample_spectra)))
- process_spectrum(msi_data, idx) do mz, intensity
- if !qc_is_empty(mz, intensity)
- # Test resolution on the most intense peak
- max_intensity_idx = argmax(intensity)
- test_mz = mz[max_intensity_idx]
-
- resolution = calculate_resolution_fwhm(test_mz, mz, intensity)
- if !isnan(resolution)
- push!(resolution_results, resolution)
- println("Spectrum $idx: Resolution = $(round(resolution))")
- end
+ for s in spectra
+ for p in s.peaks
+ push!(all_peaks, (p.mz, p.intensity, s.id))
end
end
+ sort!(all_peaks, by=p->p[1])
+
+ isempty(all_peaks) && return FeatureMatrix(zeros(0,0), Tuple{Float64,Float64}[], Int[]), Tuple{Float64,Float64}[]
+
+ # Create bins using indices into all_peaks to avoid copying large vectors
+ bins_indices = Vector{UnitRange{Int}}() # Stores UnitRange (start_idx:end_idx) for each bin
+ if isempty(all_peaks)
+ # Handle empty all_peaks case if it was not caught earlier
+ return FeatureMatrix(zeros(0,0), Tuple{Float64,Float64}[], Int[]), Tuple{Float64,Float64}[]
+ end
+
+ current_bin_start_idx = 1
+ for i in 2:length(all_peaks)
+ p_prev = all_peaks[i-1]
+ p_current = all_peaks[i]
+
+ # Approximate bin center for tolerance calculation. Using current_bin_start_idx or p_current.mz is sufficient.
+ # A more precise bin_center would require iterating over current_bin_peaks, which we are trying to avoid for memory.
+ # For tolerance calculation, a single mz value (e.g., p_current.mz or all_peaks[current_bin_start_idx].mz) is often sufficient
+ # because the m/z range within a small bin is typically very narrow.
+ approx_bin_mz = p_prev[1] # Use the m/z of the previous peak in the bin as a proxy for bin center
+ tol = (params.tolerance_unit == :ppm) ? (approx_bin_mz * params.tolerance / 1e6) : params.tolerance
+
+ if (p_current[1] - p_prev[1]) <= tol
+ # Peak is within tolerance, continue current bin
+ else
+ # Peak is outside tolerance, close current bin and start new one
+ current_bin_end_idx = i - 1
+ # Apply min_peak_per_bin filter when closing the bin
+ if params.min_peak_per_bin === nothing || (current_bin_end_idx - current_bin_start_idx + 1) >= params.min_peak_per_bin
+ push!(bins_indices, current_bin_start_idx:current_bin_end_idx)
+ end
+ current_bin_start_idx = i
+ end
+ end
+ # Add the last bin
+ current_bin_end_idx = length(all_peaks)
+ if params.min_peak_per_bin === nothing || (current_bin_end_idx - current_bin_start_idx + 1) >= params.min_peak_per_bin
+ push!(bins_indices, current_bin_start_idx:current_bin_end_idx)
+ end
+
+ # Filter bins by max width
+ if params.max_bin_width_ppm !== nothing
+ filter!(range_idx -> begin
+ bin_peaks_view = @view all_peaks[range_idx]
+ min_mz_bin = first(bin_peaks_view)[1]
+ max_mz_bin = last(bin_peaks_view)[1]
+ bin_center_approx = (min_mz_bin + max_mz_bin) / 2
+ (max_mz_bin - min_mz_bin) * 1e6 / bin_center_approx <= params.max_bin_width_ppm
+ end, bins_indices)
+ end
+
+ # Create feature matrix
+ X = zeros(Float64, ns, length(bins_indices))
+ final_bins_boundaries = Vector{Tuple{Float64,Float64}}(undef, length(bins_indices))
+
+ Threads.@threads for j in 1:length(bins_indices)
+ range_idx = bins_indices[j]
+ bin_peaks_view = @view all_peaks[range_idx] # Create a view into all_peaks for the current bin
+
+ # Calculate bin center
+ local bin_center
+ if params.intensity_weighted_centers && sum(p[2] for p in bin_peaks_view) > 0
+ weights = [p[2] for p in bin_peaks_view]
+ bin_center = sum(p[1] * p[2] for p in bin_peaks_view) / sum(weights)
+ else
+ bin_center = mean(p[1] for p in bin_peaks_view)
+ end
+
+ final_bins_boundaries[j] = (first(bin_peaks_view)[1], last(bin_peaks_view)[1])
+
+ for p in bin_peaks_view
+ s_idx = p[3]
+ # This is thread-safe because each thread writes to a unique column j
+ X[s_idx, j] = max(X[s_idx, j], p[2])
+ end
+ end
+
+ # Filter by frequency
+ if params.frequency_threshold !== nothing && params.frequency_threshold > 0
+ present_count = vec(sum(X .> 0, dims=1))
+ min_count = ceil(Int, params.frequency_threshold * ns)
+ keep_mask = findall(present_count .>= min_count)
+ X = X[:, keep_mask]
+ final_bins_boundaries = final_bins_boundaries[keep_mask]
+ end
+
+ return FeatureMatrix(X, final_bins_boundaries, collect(1:ns)), final_bins_boundaries
+ else
+ @warn "Unsupported binning method: $(params.method). Returning empty FeatureMatrix."
+ return FeatureMatrix(zeros(0,0), Tuple{Float64,Float64}[], Int[]), Tuple{Float64,Float64}[]
end
-
- if !isempty(resolution_results)
- avg_resolution = mean(resolution_results)
- println("\nAverage resolution: $(round(avg_resolution))")
- println("Resolution range: $(round(minimum(resolution_results))) - $(round(maximum(resolution_results)))")
- end
-
- # 3. Save detailed results
-
- # Save PPM error distribution
- ppm_df = DataFrame(
- theoretical_mz = [p[1] for p in accuracy_report.matched_peaks],
- measured_mz = [p[2] for p in accuracy_report.matched_peaks],
- compound = [p[3] for p in accuracy_report.matched_peaks],
- ppm_error = accuracy_report.all_ppm_errors
- )
-
- CSV.write(joinpath(output_dir, "mass_accuracy_results.csv"), ppm_df)
-
- # Save resolution results
- if !isempty(resolution_results)
- res_df = DataFrame(resolution = resolution_results)
- CSV.write(joinpath(output_dir, "resolution_results.csv"), res_df)
- end
-
- # 4. Create summary
- summary = """
- QC REPORT SUMMARY
- =================
- Date: $(now())
- File: $(filename)
- Spectra analyzed: $(length(msi_data.spectra_metadata))
-
- MASS ACCURACY:
- - Mean PPM: $(round(accuracy_report.mean_ppm, digits=2)) ppm
- - Std PPM: $(round(accuracy_report.std_ppm, digits=2)) ppm
- - Recommended tolerance: $(round(accuracy_report.optimal_ppm, digits=2)) ppm (mean + 3 * std)
-
- RESOLUTION:
- - Average: $(isempty(resolution_results) ? "N/A" : string(round(mean(resolution_results))))
- - Range: $(isempty(resolution_results) ? "N/A" : "$(round(minimum(resolution_results))) - $(round(maximum(resolution_results)))")
-
- RECOMMENDATIONS:
- - Use $(round(accuracy_report.optimal_ppm, digits=2)) ppm for peak matching
- - Instrument performance: $(accuracy_report.mean_ppm < 5 ? "Excellent" : accuracy_report.mean_ppm < 10 ? "Good" : "Needs calibration")
- """
-
- open(joinpath(output_dir, "qc_summary.txt"), "w") do f
- write(f, summary)
- end
-
- println("\nQC report saved to: $output_dir")
- return accuracy_report, resolution_results
end
diff --git a/src/PreprocessingVersioning.jl b/src/PreprocessingVersioning.jl
deleted file mode 100644
index aa26537..0000000
--- a/src/PreprocessingVersioning.jl
+++ /dev/null
@@ -1,446 +0,0 @@
-# src/PreprocessingVersioning.jl
-
-"""
-This module provides a versioning helper functions in our data preprocessing, module
-inspired by the functionality of the R package MALDIquant. It includes
-functions for quality control, intensity transformation, smoothing, baseline correction,
-normalization, peak picking, alignment, and feature matrix generation.
-
-This module is unfinished and may not end up in the release version.
-"""
-
-# =============================================================================
-# Dependencies
-# =============================================================================
-
-using UUIDs, JLD2, CodecBase # For preprocessing versioning
-
-# =============================================================================
-# Data Structures
-# =============================================================================
-
-"""
- VersionedSpectralData
-
-A struct to hold a snapshot of spectral data at a specific point in the
-preprocessing workflow. Each transformation creates a new instance of this
-struct, forming a history of all operations.
-
-# Fields
-- `version_id::String`: A unique identifier for this version of the data.
-- `parent_id::String`: The identifier of the version from which this one was derived.
-- `spectra::Vector`: The spectral data itself, typically a vector of `(mz, intensity)` tuples.
-- `processing_step::String`: A description of the operation that created this version (e.g., "smooth", "baseline").
-- `timestamp::DateTime`: The time at which this version was created.
-- `parameters::Dict`: A dictionary of the parameters used in the processing step.
-"""
-struct VersionedSpectralData
- version_id::String
- parent_id::String
- spectra::Vector
- processing_step::String
- timestamp::DateTime
- parameters::Dict
-end
-
-"""
- MSISession
-
-Manages the state of a preprocessing session, including the history of all
-data versions and the current working version.
-
-# Fields
-- `session_id::String`: A unique identifier for the entire session.
-- `original_data::VersionedSpectralData`: The initial, unprocessed spectral data.
-- `history::Vector{VersionedSpectralData}`: A chronological list of all data versions created during the session.
-- `current_version::VersionedSpectralData`: The version of the data currently being worked on or displayed.
-- `processing_steps::Vector{String}`: A human-readable log of the processing steps applied.
-"""
-mutable struct MSISession
- session_id::String
- original_data::VersionedSpectralData
- history::Vector{VersionedSpectralData}
- current_version::VersionedSpectralData
- processing_steps::Vector{String}
-end
-
-# =============================================================================
-# Session Management
-# =============================================================================
-
-"""
- initialize_session(raw_spectra::Vector, session_name::String) -> MSISession
-
-Creates a new preprocessing session from raw spectral data.
-
-# Arguments
-- `raw_spectra::Vector`: A vector of `(mz, intensity)` tuples representing the initial dataset.
-- `session_name::String`: A name for the session.
-
-# Returns
-- `MSISession`: A new session object initialized with the provided data.
-"""
-function initialize_session(raw_spectra::Vector, session_name::String)::MSISession
- original = VersionedSpectralData(
- string(uuid4()),
- "root",
- raw_spectra,
- "Original Data",
- now(),
- Dict()
- )
- return MSISession(session_name, original, [original], original, [])
-end
-
-"""
- get_processing_history(session::MSISession) -> Vector{String}
-
-Returns a list of the processing steps applied during the session.
-"""
-get_processing_history(session::MSISession) = [v.processing_step for v in session.history]
-
-"""
- undo_last_step!(session::MSISession)
-
-Reverts the session to the state before the last processing step was applied.
-This modifies the session in place.
-"""
-function undo_last_step!(session::MSISession)
- if length(session.history) > 1
- pop!(session.history)
- pop!(session.processing_steps)
- session.current_version = last(session.history)
- end
- return session
-end
-
-"""
- revert_to_version!(session::MSISession, version_id::String)
-
-Reverts the session to a specific version in its history, discarding all
-subsequent changes. This modifies the session in place.
-"""
-function revert_to_version!(session::MSISession, version_id::String)
- target_idx = findfirst(v -> v.version_id == version_id, session.history)
- if !isnothing(target_idx)
- session.history = session.history[1:target_idx]
- session.current_version = session.history[end]
- # Also truncate the descriptive processing steps
- num_steps_to_keep = max(0, target_idx - 1)
- session.processing_steps = session.processing_steps[1:num_steps_to_keep]
- end
- return session
-end
-
-# =============================================================================
-# Versioned Preprocessing
-# =============================================================================
-
-"""
- apply_processing_step(spectra::Vector, step::Symbol, params::Dict) -> Vector
-
-A dispatcher that applies a single, specified preprocessing function to a set of spectra.
-This is the core function called by the versioned workflow.
-
-# Arguments
-- `spectra::Vector`: The input spectral data.
-- `step::Symbol`: The symbol representing the processing step (e.g., `:smooth`, `:baseline`).
-- `params::Dict`: A dictionary of parameters for the step.
-
-# Returns
-- `Vector`: The processed spectral data.
-"""
-function apply_processing_step(spectra::Vector, step::Symbol, params::Dict)
- processed = deepcopy(spectra)
- params_sym = Dict(Symbol(k) => v for (k,v) in params) # Ensure keys are symbols
-
- if step === :transform
- for i in eachindex(processed)
- mz, y = processed[i]
- processed[i] = (mz, transform_intensity(y; params_sym...))
- end
- elseif step === :smooth
- for i in eachindex(processed)
- mz, y = processed[i]
- processed[i] = (mz, smooth_spectrum(y; params_sym...))
- end
- elseif step === :baseline
- for i in eachindex(processed)
- mz, y = processed[i]
- baseline = snip_baseline(y; params_sym...)
- processed[i] = (mz, max.(0.0, y .- baseline))
- end
- elseif step === :normalize
- mode = get(params, :normalize_method, :tic)
- if mode === :pqn
- matrix = hcat([float.(p[2]) for p in processed]...)
- matrix_norm = pqn_normalize(matrix)
- for i in eachindex(processed)
- mz, _ = processed[i]
- processed[i] = (mz, view(matrix_norm, :, i))
- end
- else # :tic or :median
- norm_func = (mode === :tic) ? tic_normalize : median_normalize
- for i in eachindex(processed)
- mz, y = processed[i]
- processed[i] = (mz, norm_func(y))
- end
- end
- elseif step === :peaks
- peak_results = Vector{Tuple{Vector{Float64},Vector{Float64}}}(undef, length(processed))
- for (i, (mz, y)) in enumerate(processed)
- pk_mz, pk_int = detect_peaks_profile(mz, y; params_sym...)
- peak_results[i] = (pk_mz, pk_int)
- end
- return peak_results
- elseif step === :peaks_wavelet
- peak_results = Vector{Tuple{Vector{Float64},Vector{Float64}}}(undef, length(processed))
- for (i, (mz, y)) in enumerate(processed)
- pk_mz, pk_int = detect_peaks_wavelet(mz, y; params_sym...)
- peak_results[i] = (pk_mz, pk_int)
- end
- return peak_results
- elseif step === :calibrate
- return calibrate_spectra(processed; params_sym...)
- else
- @warn "Unsupported processing step: $step"
- end
- return processed
-end
-
-"""
- run_versioned_preprocessing(session::MSISession, step::Symbol, params::Dict) -> VersionedSpectralData
-
-Creates a new version of spectral data by applying a single processing step to the
-current version in the session.
-
-# Arguments
-- `session::MSISession`: The current processing session.
-- `step::Symbol`: The processing step to apply.
-- `params::Dict`: Parameters for the processing step.
-
-# Returns
-- `VersionedSpectralData`: A new data version with the transformation applied.
-"""
-function run_versioned_preprocessing(session::MSISession, step::Symbol, params::Dict)
- # Create a new version based on the current one
- parent_version = session.current_version
- new_version_id = string(uuid4())
-
- # Apply the processing step
- new_spectra = apply_processing_step(parent_version.spectra, step, params)
-
- # Create the new versioned data object
- new_version = VersionedSpectralData(
- new_version_id,
- parent_version.version_id,
- new_spectra,
- string(step),
- now(),
- params
- )
-
- return new_version
-end
-
-"""
- apply_processing_with_version!(session::MSISession, step::Symbol, ui_params::Dict)
-
-A high-level function to apply a processing step, create a new version, and
-update the session state in place.
-
-# Arguments
-- `session::MSISession`: The session to modify.
-- `step::Symbol`: The processing step to apply.
-- `ui_params::Dict`: A dictionary of parameters, typically from a UI.
-"""
-function apply_processing_with_version!(session::MSISession, step::Symbol, ui_params::Dict)
- # In a real app, you would validate and convert UI params here.
- processing_params = ui_params
-
- # Create the new version
- new_version = run_versioned_preprocessing(session, step, processing_params)
-
- # Update the session
- push!(session.history, new_version)
- session.current_version = new_version
-
- # Describe the step for the history log
- step_description = "$(string(step)) with params: " * join(["$k=$v" for (k,v) in processing_params], ", ")
- push!(session.processing_steps, step_description)
-
- return session
-end
-
-# =============================================================================
-# Data Serialization & Export
-# =============================================================================
-
-"""
- save_spectral_version(data::VersionedSpectralData, filepath::String)
-
-Saves a `VersionedSpectralData` object to a file using the JLD2 format.
-"""
-function save_spectral_version(data::VersionedSpectralData, filepath::String)
- JLD2.save_object(filepath, data)
-end
-
-"""
- load_spectral_version(filepath::String) -> VersionedSpectralData
-
-Loads a `VersionedSpectralData` object from a JLD2 file.
-"""
-function load_spectral_version(filepath::String)::VersionedSpectralData
- return JLD2.load_object(filepath)
-end
-
-"""
- _encode_base64(data::AbstractVector{<:Real}) -> String
-
-Helper function to convert a numeric vector into a Base64 encoded string.
-"""
-function _encode_base64(data::AbstractVector{T}) where T <: Real
- bytes = reinterpret(UInt8, data)
- return base64encode(bytes)
-end
-
-"""
- export_to_mzml(session::MSISession, filepath::String)
-
-Exports the current version of the spectral data in a session to a standard
-.mzML file.
-
-# Arguments
-- `session::MSISession`: The current session.
-- `filepath::String`: The path for the output .mzML file.
-"""
-function export_to_mzml(session::MSISession, filepath::String)
- spectra_to_export = session.current_version.spectra
-
- open(filepath, "w") do f
- # XML Header
- write(f, "\n")
- write(f, "\n")
-
- # CV List
- write(f, " \n")
- write(f, " \n")
- write(f, " \n")
- write(f, " \n")
-
- # File Description
- write(f, " \n \n")
- write(f, " \n")
- write(f, " \n \n")
-
- # Run and Spectrum List
- write(f, " \n")
- write(f, " \n")
-
- for (i, (mz, intensity)) in enumerate(spectra_to_export)
- # Ensure data is in the correct format for encoding
- mz_64 = convert(Vector{Float64}, mz)
- int_32 = convert(Vector{Float32}, intensity)
-
- # Base64 encode the binary data
- mz_b64 = _encode_base64(mz_64)
- int_b64 = _encode_base64(int_32)
-
- write(f, " \n")
- write(f, " \n")
- # Assuming profile mode for processed data, could be made dynamic
- write(f, " \n")
-
- write(f, " \n")
-
- # m/z array
- write(f, " \n")
- write(f, " \n")
- write(f, " \n")
- write(f, " \n")
- write(f, " $(mz_b64)\n")
- write(f, " \n")
-
- # Intensity array
- write(f, " \n")
- write(f, " \n")
- write(f, " \n")
- write(f, " \n")
- write(f, " $(int_b64)\n")
- write(f, " \n")
-
- write(f, " \n")
- write(f, " \n")
- end
-
- write(f, " \n")
- write(f, " \n")
-
- # Dummy instrument and data processing info
- write(f, " \n")
- write(f, " \n")
- write(f, " \n")
- write(f, " \n")
- write(f, " \n")
- write(f, " \n")
-
- write(f, "\n")
- end
- println("Successfully exported current data to $filepath")
-end
-
-
-# =============================================================================
-# Stubs for UI Integration
-# =============================================================================
-
-"""
- validate_ui_parameters(step::Symbol, ui_params::Dict) -> Tuple{Bool, String}
-
-(Stub) Validates parameters from a UI before they are used in a processing step.
-"""
-function validate_ui_parameters(step::Symbol, ui_params::Dict)::Tuple{Bool, String}
- # In a real implementation, this would check types, ranges, etc.
- # For example, for :smooth, ensure 'sg_window' is an odd integer.
- println("Validating parameters for step: $step")
- return (true, "Parameters are valid.")
-end
-
-"""
- run_processing_with_progress(session, step, params, progress_callback)
-
-(Stub) A wrapper for running a processing step that includes a progress reporting callback.
-"""
-function run_processing_with_progress(session, step, params, progress_callback)
- progress_callback(0.0, "Starting $step...")
-
- # This is a simplified example. Real implementation would need to
- # hook into the loops inside apply_processing_step.
- new_version = run_versioned_preprocessing(session, step, params)
-
- progress_callback(1.0, "Finished $step.")
- return new_version
-end
-
-"""
- safe_processing_application!(session::MSISession, step::Symbol, params::Dict)
-
-(Stub) A safe wrapper to apply a processing step that includes error handling and rollback.
-"""
-function safe_processing_application!(session::MSISession, step::Symbol, params::Dict)
- num_history = length(session.history)
- try
- println("Safely applying step: $step")
- apply_processing_with_version!(session, step, params)
- catch e
- @error "Processing step '$step' failed!" exception=(e, catch_backtrace())
- # Roll back to the previous state
- while length(session.history) > num_history
- undo_last_step!(session)
- end
- println("Session has been rolled back to the previous state.")
- return session # Return the rolled-back session
- end
- return session # Return the updated session
-end
-
diff --git a/src/imzML.jl b/src/imzML.jl
index 6d66485..51dd401 100644
--- a/src/imzML.jl
+++ b/src/imzML.jl
@@ -17,6 +17,134 @@ Core Functions:
#
# ============================================================================
+function parse_imzml_header(stream::IO)
+ # Initialize all return values and temporary variables
+ instrument_meta = InstrumentMetadata()
+ param_groups = Dict{String, SpecDim}()
+ img_dims = zeros(Int32, 3)
+ vendor_preprocessing_steps = String[] # Initialize vendor_preprocessing_steps
+
+ # Temp vars for parsing
+ resolution = instrument_meta.resolution
+ instrument_model = instrument_meta.instrument_model
+ mass_accuracy_ppm = instrument_meta.mass_accuracy_ppm
+ polarity = instrument_meta.polarity
+ calibration_status = instrument_meta.calibration_status
+ laser_settings = Dict{String, Any}()
+
+ try
+ seekstart(stream)
+
+ in_scan_settings = false
+ dim_axes_found = 0
+ found_spectrum_count = false
+
+ # This loop performs a single pass to get most metadata
+ while !eof(stream)
+ line = readline(stream)
+
+ # State-machine like logic to enter/exit sections
+ if occursin("", line)
+ in_scan_settings = false
+ end
+
+ # --- Parse Scan Settings for Image Dimensions ---
+ if in_scan_settings && occursin(" 0 && !eof(stream)
- currLine = readline(stream)
- if occursin("", line)
- break
- end
- end
- spectrum_xml = String(take!(spectrum_buffer))
-
- spectrum_mode = global_mode
- if occursin("MS:1000127", spectrum_xml)
- spectrum_mode = CENTROID
- elseif occursin("MS:1000128", spectrum_xml)
- spectrum_mode = PROFILE
- end
- seek(stream, current_pos)
-
- # CRITICAL FIX: Use correct data types that match the converter output
- mz_asset = SpectrumAsset(Float64, mz_is_compressed, mz_offset, nPoints, :mz)
- int_asset = SpectrumAsset(Float32, int_is_compressed, int_offset, nPoints, :intensity)
-
- spectra_metadata[k] = SpectrumMetadata(x, y, "", spectrum_mode, mz_asset, int_asset)
-
- current_ibd_offset += mz_len_bytes + int_len_bytes
-
- # Skip to next spectrum
- line = readuntil(stream, "")
- end
-
- return spectra_metadata
-end
-=#
-
function parse_compressed(stream::IO, hIbd::Union{IO, ThreadSafeFileHandle}, param_groups::Dict{String, SpecDim},
width::Int32, height::Int32, num_spectra::Int32,
default_mz_format::DataType, default_intensity_format::DataType,
@@ -532,7 +548,7 @@ function parse_compressed(stream::IO, hIbd::Union{IO, ThreadSafeFileHandle}, par
int_is_compressed ? int_info.encoded_length : int_info.array_length, :intensity)
end
- spectra_metadata[k] = SpectrumMetadata(x, y, "", spectrum_mode, mz_asset, int_asset)
+ spectra_metadata[k] = SpectrumMetadata(x, y, "", :sample, spectrum_mode, mz_asset, int_asset)
end
return spectra_metadata
@@ -675,7 +691,7 @@ function parse_neofx(stream::IO, hIbd::Union{IO, ThreadSafeFileHandle}, param_gr
int_is_compressed ? int_info.encoded_length : int_info.array_length, :intensity)
end
- spectra_metadata[k] = SpectrumMetadata(x, y, "", spectrum_mode, mz_asset, int_asset)
+ spectra_metadata[k] = SpectrumMetadata(x, y, "", :sample, spectrum_mode, mz_asset, int_asset)
end
return spectra_metadata
diff --git a/src/mzML.jl b/src/mzML.jl
index 1550a08..f795f9b 100644
--- a/src/mzML.jl
+++ b/src/mzML.jl
@@ -19,6 +19,107 @@ const DATA_FORMAT_ACCESSIONS = Dict{String, DataType}(
"MS:1000523" => Float64
)
+function parse_instrument_metadata_mzml(stream::IO)
+ # Initialize with default values from the InstrumentMetadata constructor
+ instrument_meta = InstrumentMetadata()
+
+ # Create temporary variables to hold parsed values
+ resolution = instrument_meta.resolution
+ instrument_model = instrument_meta.instrument_model
+ mass_accuracy_ppm = instrument_meta.mass_accuracy_ppm
+ polarity = instrument_meta.polarity
+ calibration_status = instrument_meta.calibration_status
+ laser_settings = Dict{String, Any}() # Use Any for heterogeneous values
+ vendor_preprocessing_steps = String[] # Initialize vendor_preprocessing_steps
+
+ try
+ # Reset stream and read a sufficiently large header block to find metadata.
+ seekstart(stream)
+ header_block = ""
+ header_read_limit = 50000 # Read up to 50KB of header
+ bytes_read = 0
+
+ while !eof(stream) && bytes_read < header_read_limit
+ line_pos = position(stream)
+ line = readline(stream)
+ bytes_read += (position(stream) - line_pos)
+
+ # Stop at the start of the main data section
+ if occursin("", line) || occursin("", line)
+ break
+ end
+ header_block *= line * "\n"
+ end
+
+ header_io = IOBuffer(header_block)
+ while !eof(header_io)
+ line = readline(header_io)
+
+ if occursin(" m == CENTROID, modes)
+ num_profile = count(m -> m == PROFILE, modes)
+
+ acq_mode_symbol = if num_centroid > 0 && num_profile == 0
+ :centroid
+ elseif num_profile > 0 && num_centroid == 0
+ :profile
+ elseif num_centroid > 0 && num_profile > 0
+ :mixed
+ else
+ :unknown
+ end
+
+ final_instrument_meta = InstrumentMetadata(
+ instrument_meta.resolution,
+ acq_mode_symbol, # Update with parsed mode
+ instrument_meta.mz_axis_type,
+ instrument_meta.calibration_status,
+ instrument_meta.instrument_model,
+ instrument_meta.mass_accuracy_ppm,
+ instrument_meta.laser_settings,
+ instrument_meta.polarity,
+ instrument_meta.vendor_preprocessing_steps # Add this new field
+ )
+
source = MzMLSource(ts_stream, mz_format, intensity_format)
println("DEBUG: Creating MSIData object.")
- return MSIData(source, spectra_metadata, (0, 0), nothing, cache_size)
+ return MSIData(source, spectra_metadata, final_instrument_meta, (0, 0), nothing, cache_size)
catch e
close(ts_stream) # Ensure stream is closed on error
diff --git a/test/run_descriptive_preprocessing.jl b/test/run_descriptive_preprocessing.jl
new file mode 100644
index 0000000..a5085da
--- /dev/null
+++ b/test/run_descriptive_preprocessing.jl
@@ -0,0 +1,222 @@
+using CairoMakie
+using Statistics
+using Colors
+
+function create_msi_parameter_plot()
+ # Create the figure with subplots for each preprocessing step
+ fig = Figure(size=(1600, 1200), fontsize=12)
+
+ # Define the preprocessing steps and their parameters
+ steps = [
+ ("BaselineCorrection", ["iterations: 10", "method: SNIP", "window: 3.0"]),
+ ("Calibration", ["fit_order: 2", "method: internal_standards", "ppm_tolerance: 20.0"]),
+ ("Normalization", ["method: rms"]),
+ ("PeakAlignment", ["max_shift_ppm: 50.0", "method: linear", "tolerance: 20.0 ppm"]),
+ ("PeakBinningParams", ["max_bin_width_ppm: 60.0", "method: adaptive", "min_peak_per_bin: 3"]),
+ ("PeakPicking", ["half_window: 5", "method: centroid", "snr_threshold: 15.0"]),
+ ("PeakSelection", ["max_fwhm_ppm: 95.39", "min_fwhm_ppm: 19.08", "min_snr: 3.0"]),
+ ("Smoothing", ["method: savitzky_golay", "order: 3", "window: 5"])
+ ]
+
+ # Create a grid of subplots
+ g = fig[1, 1] = GridLayout()
+
+ # Plot each preprocessing step
+ for (idx, (step_name, params)) in enumerate(steps)
+ row, col = fldmod1(idx, 2)
+ ax = Axis(g[row, col], title=step_name, titlesize=14)
+
+ # Generate simulated m/z values and intensities
+ mz_range = range(100, 1000, length=500)
+
+ if step_name == "BaselineCorrection"
+ # Simulate spectrum with baseline
+ true_signal = 0.5 .* exp.(-0.001 .* (mz_range .- 400).^2) .+
+ 0.3 .* exp.(-0.002 .* (mz_range .- 600).^2)
+ baseline = 0.1 .+ 0.05 .* sin.(mz_range ./ 50)
+ noisy_signal = true_signal .+ baseline .+ 0.02 .* randn(length(mz_range))
+ corrected = noisy_signal .- baseline
+
+ lines!(ax, mz_range, noisy_signal, color=:blue, linewidth=2, label="Raw")
+ lines!(ax, mz_range, baseline, color=:red, linewidth=2, linestyle=:dash, label="Baseline")
+ lines!(ax, mz_range, corrected, color=:green, linewidth=2, label="Corrected")
+
+ elseif step_name == "Calibration"
+ # Simulate calibration shift
+ reference_peaks = [200, 400, 600, 800]
+ measured_peaks = reference_peaks .+ 2.0 .* randn(length(reference_peaks))
+
+ scatter!(ax, reference_peaks, fill(0.5, length(reference_peaks)),
+ color=:red, markersize=15, label="Reference")
+ scatter!(ax, measured_peaks, fill(0.3, length(measured_peaks)),
+ color=:blue, markersize=10, label="Measured")
+
+ # Add calibration lines
+ for i in 1:length(reference_peaks)
+ lines!(ax, [measured_peaks[i], reference_peaks[i]], [0.3, 0.5],
+ color=:black, linewidth=1, linestyle=:dash)
+ end
+
+ elseif step_name == "Normalization"
+ # Simulate normalization effect
+ spectra = [
+ 0.8 .* exp.(-0.001 .* (mz_range .- 300).^2) .+ 0.2 .* randn(length(mz_range)),
+ 1.2 .* exp.(-0.001 .* (mz_range .- 300).^2) .+ 0.2 .* randn(length(mz_range)),
+ 0.9 .* exp.(-0.001 .* (mz_range .- 300).^2) .+ 0.2 .* randn(length(mz_range))
+ ]
+
+ normalized_spectra = [spec ./ std(spec) for spec in spectra]
+
+ for (i, spec) in enumerate(spectra)
+ lines!(ax, mz_range, spec .+ i*0.3, color=RGBA(1, 0, 0, 0.6), linewidth=2,
+ label=i==1 ? "Before Norm" : "")
+ end
+ for (i, spec) in enumerate(normalized_spectra)
+ lines!(ax, mz_range, spec .+ i*0.3, color=RGBA(0, 0, 1, 0.6), linewidth=2,
+ label=i==1 ? "After Norm" : "")
+ end
+
+ elseif step_name == "PeakAlignment"
+ # Simulate peak alignment
+ base_peaks = [300, 500, 700]
+ shifts = [-15, 5, 10]
+
+ for (i, shift) in enumerate(shifts)
+ shifted_peaks = base_peaks .+ shift
+ aligned_peaks = base_peaks
+
+ scatter!(ax, shifted_peaks, fill(i, length(shifted_peaks)),
+ color=:red, markersize=12, label=i==1 ? "Before Align" : "")
+ scatter!(ax, aligned_peaks, fill(i+0.3, length(aligned_peaks)),
+ color=:green, markersize=12, label=i==1 ? "After Align" : "")
+
+ # Show alignment lines
+ for j in 1:length(base_peaks)
+ lines!(ax, [shifted_peaks[j], aligned_peaks[j]], [i, i+0.3],
+ color=:black, linewidth=1, linestyle=:dash)
+ end
+ end
+
+ elseif step_name == "PeakBinningParams"
+ # Simulate peak binning with centroids
+ raw_peaks_mz = 100:25:900
+ raw_peaks_intensity = rand(length(raw_peaks_mz))
+
+ # Create binned peaks (wider bins)
+ bin_centers = 150:60:850
+ bin_intensities = [sum(raw_peaks_intensity[abs.(raw_peaks_mz .- center) .< 30])
+ for center in bin_centers] .* 0.8
+
+ # Profile mode (continuous)
+ profile_signal = zeros(length(mz_range))
+ for (mz, int) in zip(raw_peaks_mz, raw_peaks_intensity)
+ profile_signal .+= int .* exp.(-0.001 .* (mz_range .- mz).^2)
+ end
+
+ lines!(ax, mz_range, profile_signal, color=:blue, linewidth=2, label="Profile")
+ barplot!(ax, bin_centers, bin_intensities, color=RGBA(1, 0, 0, 0.7),
+ width=50, label="Binned Centroids")
+
+ elseif step_name == "PeakPicking"
+ # Simulate peak picking from profile to centroids
+ profile_signal = 0.6 .* exp.(-0.0005 .* (mz_range .- 400).^2) .+
+ 0.4 .* exp.(-0.0008 .* (mz_range .- 650).^2) .+
+ 0.1 .* randn(length(mz_range))
+
+ # Simulate picked peaks (centroids)
+ peak_positions = [380, 405, 640, 660]
+ peak_intensities = [0.5, 0.6, 0.35, 0.4]
+
+ lines!(ax, mz_range, profile_signal, color=:blue, linewidth=2, label="Profile Spectrum")
+ scatter!(ax, peak_positions, peak_intensities, color=:red, markersize=20,
+ label="Picked Centroids", strokewidth=2)
+
+ elseif step_name == "PeakSelection"
+ # Simulate peak selection based on criteria
+ all_peaks_mz = 200:50:800
+ all_peaks_fwhm = rand(length(all_peaks_mz)) .* 100 .+ 10
+ all_peaks_snr = rand(length(all_peaks_mz)) .* 10
+
+ # Selection criteria
+ selected = (all_peaks_fwhm .>= 19.08) .& (all_peaks_fwhm .<= 95.39) .& (all_peaks_snr .>= 3.0)
+
+ scatter!(ax, all_peaks_mz[.!selected], all_peaks_fwhm[.!selected],
+ color=:red, markersize=15, label="Rejected")
+ scatter!(ax, all_peaks_mz[selected], all_peaks_fwhm[selected],
+ color=:green, markersize=15, label="Selected")
+
+ # Add selection criteria lines
+ hlines!(ax, [19.08, 95.39], color=:black, linestyle=:dash, linewidth=2)
+ text!(ax, 850, 50; text="FWHM bounds", color=:black, fontsize=10)
+
+ elseif step_name == "Smoothing"
+ # Simulate smoothing effect
+ true_signal = 0.7 .* exp.(-0.001 .* (mz_range .- 450).^2) .+
+ 0.5 .* exp.(-0.0008 .* (mz_range .- 650).^2)
+ noisy_signal = true_signal .+ 0.1 .* randn(length(mz_range))
+
+ # Simple smoothing simulation
+ smoothed = similar(noisy_signal)
+ window = 5
+ for i in 1:length(noisy_signal)
+ start_idx = max(1, i - window ÷ 2)
+ end_idx = min(length(noisy_signal), i + window ÷ 2)
+ smoothed[i] = mean(noisy_signal[start_idx:end_idx])
+ end
+
+ lines!(ax, mz_range, noisy_signal, color=:red, linewidth=1, label="Noisy")
+ lines!(ax, mz_range, smoothed, color=:blue, linewidth=2, label="Smoothed")
+ lines!(ax, mz_range, true_signal, color=:green, linewidth=1, linestyle=:dash, label="True")
+ end
+
+ # Add parameter text using a more reliable approach
+ param_text = join(params, "\n")
+
+ text!(ax, param_text, position=Point2f(0.05, 0.95),
+ space=:relative, align=(:left, :top), color=:black,
+ fontsize=10, font=:regular)
+
+ # Add legend for selected plots
+ if idx <= 4
+ axislegend(ax, position=:rt, framevisible=true, backgroundcolor=RGBA(1,1,1,0.8))
+ end
+
+ # Customize axes
+ ax.xlabel = "m/z"
+ ax.ylabel = idx in [1,3,5,7] ? "Intensity" : ""
+ ax.xgridvisible = false
+ ax.ygridvisible = false
+ end
+
+ # Add overall title
+ Label(fig[0, :], "MSI Preprocessing Pipeline Parameters and Simulations",
+ fontsize=18, font=:bold, padding=(0, 0, 10, 0))
+
+ # Add explanation
+ explanation = """
+ Simulation of MSI preprocessing parameters showing:
+ • Blue lines: Profile/continuous spectra
+ • Red bars/points: Centroid data
+ • Dashed lines: Reference/true signals
+ • Green: Processed/corrected data
+ • Each subplot demonstrates key parameters for the preprocessing step
+ """
+
+ Label(fig[2, :], explanation, fontsize=12, tellwidth=false, padding=(10, 10, 10, 10))
+
+ # Adjust layout
+ colgap!(g, 20)
+ rowgap!(g, 20)
+
+ fig
+end
+
+# Create and display the plot
+fig = create_msi_parameter_plot()
+
+# Save the plot
+save("msi_preprocessing_parameters.png", fig)
+println("Plot saved as 'msi_preprocessing_parameters.png'")
+
+# Display the plot (if in an interactive environment)
+fig
diff --git a/test/run_precalculation_dict.jl b/test/run_precalculation_dict.jl
new file mode 100644
index 0000000..2726407
--- /dev/null
+++ b/test/run_precalculation_dict.jl
@@ -0,0 +1,157 @@
+# test/run_precalculation_example.jl
+
+using Printf
+import Pkg
+
+# --- Load the MSI_src Module ---
+Pkg.activate(joinpath(@__DIR__, ".."))
+using MSI_src
+
+# ===================================================================
+# CONFIG: PLEASE FILL IN YOUR FILE PATHS HERE
+# ===================================================================
+
+# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML"
+#const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/mzML/Col_1.mzML"
+const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Thricoderma_etc/Imaging_interaccion_trichoderma_vs_streptomyces.mzML"
+
+# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.imzML"
+#const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML"
+const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Thricoderma_etc/Imaging_interaccion_trichoderma_vs_streptomyces.imzML"
+
+#const MASK_ROUTE = "/home/pixel/Documents/Cinvestav_2025/JuliaMSI/public/css/masks/Stomach_DHB_uncompressed.png"
+const MASK_ROUTE = ""
+
+const reference_peaks = Dict(
+ # DHB Matrix peaks (should be present)
+ 137.0244 => "DHB_fragment",
+ 155.0349 => "DHB_M+H",
+ 177.0168 => "DHB_M+Na",
+
+ # Common lipids in your mass range
+ 496.3398 => "PC_16:0_16:0",
+ 520.3398 => "PC_16:0_18:1",
+ 760.5851 => "PC_16:0_18:1_Na",
+
+ # Common contaminants
+ 391.2843 => "PDMS",
+ 413.2662 => "PDMS_Na",
+
+ # Add some high mass peaks
+ 842.5092 => "Protein_standard",
+ 1045.532 => "Protein_standard",
+
+ 290.1747 => "Atropine [M+H]+",
+ 304.1903 => "Scopolamine [M+H]+",
+ 124.0393 => "Tropine [M+H]+",
+)
+
+# ===================================================================
+# HELPER FUNCTIONS FOR PRINTING
+# ===================================================================
+
+function print_header(title::String)
+ println("\n" * "="^80)
+ println("$(title)")
+ println("="^80)
+end
+
+function check_data_range(msi_data::MSIData)
+ println("\n--- Data Range Analysis ---")
+
+ min_mz, max_mz = get_global_mz_range(msi_data)
+ if isfinite(min_mz) && isfinite(max_mz) && min_mz < max_mz
+ println("Global m/z range: [$(min_mz), $(max_mz)]")
+ else
+ println("Global m/z range: Not yet determined or invalid (initial: [$(min_mz), $(max_mz)])")
+ end
+
+ # Check a few spectra to see actual m/z values
+ println("\nChecking first few spectra for actual m/z values:")
+ for i in 1:min(3, length(msi_data.spectra_metadata))
+ try
+ mz, intensity = GetSpectrum(msi_data, i)
+ if !isempty(mz)
+ println("Spectrum $i: m/z range [$(minimum(mz)), $(maximum(mz))], length=$(length(mz))")
+ # Print first and last few m/z values
+ if length(mz) > 10
+ println(" First 5 m/z: $(mz[1:5])")
+ println(" Last 5 m/z: $(mz[end-4:end])")
+ end
+ end
+ catch e
+ println("Spectrum $i: Error - $e")
+ end
+ end
+end
+
+# ===================================================================
+# MAIN EXAMPLE RUNNER
+# ===================================================================
+
+function run_precalculation_example()
+ # --- Process mzML file ---
+ print_header("Processing mzML File: $(basename(TEST_MZML_FILE))")
+ if !isfile(TEST_MZML_FILE)
+ println("SKIPPING: mzML file not found at $(TEST_MZML_FILE)")
+ else
+ try
+ msi_data_mzml = @time OpenMSIData(TEST_MZML_FILE)
+ check_data_range(msi_data_mzml)
+ analysis_results_mzml = main_precalculation(msi_data_mzml, reference_peaks=reference_peaks)
+
+ println("\n" * "*"^80)
+ println("MZML PREPROCESSING ANALYSIS RESULTS")
+ println("*"^80)
+ for (step_name, params) in analysis_results_mzml
+ println("\n--- $(uppercase(string(step_name))) Parameters ---")
+ if isempty(params)
+ println(" No recommended parameters.")
+ else
+ for (param_name, param_value) in params
+ println(" $param_name: $param_value")
+ end
+ end
+ end
+
+ close(msi_data_mzml) # Close file handles
+ catch e
+ println("ERROR processing mzML file: $e")
+ showerror(stdout, e, catch_backtrace())
+ end
+ end
+
+ # --- Process imzML file ---
+ print_header("Processing imzML File: $(basename(TEST_IMZML_FILE))")
+ if !isfile(TEST_IMZML_FILE)
+ println("SKIPPING: imzML file not found at $(TEST_IMZML_FILE)")
+ else
+ try
+ msi_data_imzml = @time OpenMSIData(TEST_IMZML_FILE)
+ check_data_range(msi_data_imzml)
+ analysis_results_imzml = main_precalculation(msi_data_imzml, reference_peaks=reference_peaks, mask_path=MASK_ROUTE)
+
+ println("\n" * "*"^80)
+ println("IMZML PREPROCESSING ANALYSIS RESULTS")
+ println("*"^80)
+ for (step_name, params) in analysis_results_imzml
+ println("\n--- $(uppercase(string(step_name))) Parameters ---")
+ if isempty(params)
+ println(" No recommended parameters.")
+ else
+ for (param_name, param_value) in params
+ println(" $param_name: $param_value")
+ end
+ end
+ end
+
+ close(msi_data_imzml) # Close file handles
+ catch e
+ println("ERROR processing imzML file: $e")
+ showerror(stdout, e, catch_backtrace())
+ end
+ end
+end
+
+# --- Execute ---
+@time run_precalculation_example()
diff --git a/test/run_precalculation_example.jl b/test/run_precalculation_example.jl
new file mode 100644
index 0000000..05f7cec
--- /dev/null
+++ b/test/run_precalculation_example.jl
@@ -0,0 +1,281 @@
+# test/run_precalculation_example.jl
+
+using Printf
+import Pkg
+
+# --- Load the MSI_src Module ---
+Pkg.activate(joinpath(@__DIR__, ".."))
+using MSI_src
+
+# ===================================================================
+# CONFIG: PLEASE FILL IN YOUR FILE PATHS HERE
+# ===================================================================
+
+# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML"
+# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/mzML/Col_1.mzML"
+const TEST_MZML_FILE = ""
+# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.imzML"
+#const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML"
+const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Thricoderma_etc/Imaging_interaccion_trichoderma_vs_streptomyces.imzML"
+#const MASK_ROUTE = "/home/pixel/Documents/Cinvestav_2025/JuliaMSI/public/css/masks/Stomach_DHB_uncompressed.png"
+const MASK_ROUTE = ""
+
+#=
+const reference_peaks = Dict(
+ # Common ESI positive mode reference compounds
+ 121.0509 => "Purine",
+ 149.0233 => "HP-921",
+ 322.0481 => "Hexakis(1H,1H,3H-tetrafluoropropoxy)phosphazine",
+ 622.0290 => "Hexakis(2,2-difluoroethoxy)phosphazine",
+
+ # Atropine and related compounds in positive mode
+ 290.1747 => "Atropine [M+H]+",
+ 304.1903 => "Scopolamine [M+H]+",
+ 124.0393 => "Tropine [M+H]+",
+
+ # Common contaminants and lock masses
+ 391.2843 => "Polydimethylcyclosiloxane [M+H]+",
+ 413.2662 => "Polydimethylcyclosiloxane [M+Na]+",
+ 429.2402 => "Polydimethylcyclosiloxane [M+K]+"
+)
+=#
+
+const reference_peaks = Dict(
+ # DHB Matrix peaks (should be present)
+ 137.0244 => "DHB_fragment",
+ 155.0349 => "DHB_M+H",
+ 177.0168 => "DHB_M+Na",
+
+ # Common lipids in your mass range
+ 496.3398 => "PC_16:0_16:0",
+ 520.3398 => "PC_16:0_18:1",
+ 760.5851 => "PC_16:0_18:1_Na",
+
+ # Common contaminants
+ 391.2843 => "PDMS",
+ 413.2662 => "PDMS_Na",
+
+ # Add some high mass peaks
+ 842.5092 => "Protein_standard",
+ 1045.532 => "Protein_standard",
+)
+
+# ===================================================================
+# HELPER FUNCTIONS FOR PRINTING
+# ===================================================================
+
+function print_header(title::String)
+ println("\n" * "="^80)
+ println("$(title)")
+ println("="^80)
+end
+
+function print_subheader(title::String)
+ println("\n" * "-"^80)
+ println("$(title)")
+ println("-"^80)
+end
+
+function print_param(key, value)
+ if value === nothing || (isa(value, Number) && isnan(value))
+ println(" " * "" * rpad(key, 30) * ": unable to get, user needs to input manually")
+ elseif isa(value, Symbol) && startswith(string(key), "method") # Heuristic for method selection
+ println(" " * "" * rpad(key, 30) * ": automatically detected this as the most optimal method: $(value)")
+ else
+ println(" " * "" * rpad(key, 30) * ": $(value)")
+ end
+end
+
+function print_recommendations(recommendations::Dict)
+ for (step, params) in recommendations
+ print_subheader("Recommendations for $(String(step))")
+ for (key, value) in params
+ print_param(key, value)
+ end
+ end
+end
+
+function check_data_range(msi_data::MSIData)
+ println("\n--- Data Range Analysis ---")
+
+ min_mz, max_mz = get_global_mz_range(msi_data)
+ if isfinite(min_mz) && isfinite(max_mz) && min_mz < max_mz
+ println("Global m/z range: [$(min_mz), $(max_mz)]")
+ else
+ println("Global m/z range: Not yet determined or invalid (initial: [$(min_mz), $(max_mz)])")
+ end
+
+ # Check a few spectra to see actual m/z values
+ println("\nChecking first few spectra for actual m/z values:")
+ for i in 1:min(3, length(msi_data.spectra_metadata))
+ try
+ mz, intensity = GetSpectrum(msi_data, i)
+ if !isempty(mz)
+ println("Spectrum $i: m/z range [$(minimum(mz)), $(maximum(mz))], length=$(length(mz))")
+ # Print first and last few m/z values
+ if length(mz) > 10
+ println(" First 5 m/z: $(mz[1:5])")
+ println(" Last 5 m/z: $(mz[end-4:end])")
+ end
+ end
+ catch e
+ println("Spectrum $i: Error - $e")
+ end
+ end
+end
+
+# ===================================================================
+# MAIN EXAMPLE RUNNER
+# ===================================================================
+
+function run_precalculation_example()
+ # --- Process mzML file ---
+ print_header("Processing mzML File: $(basename(TEST_MZML_FILE))")
+ if !isfile(TEST_MZML_FILE)
+ println("SKIPPING: mzML file not found at $(TEST_MZML_FILE)")
+ else
+ try
+ msi_data_mzml = @time OpenMSIData(TEST_MZML_FILE)
+ check_data_range(msi_data_mzml)
+ analysis_results_mzml = run_preprocessing_analysis(msi_data_mzml, reference_peaks=reference_peaks)
+
+ println("\n" * "*"^80)
+ println("MZML PREPROCESSING ANALYSIS RESULTS")
+ println("*"^80)
+
+ for (phase, results) in analysis_results_mzml
+ if phase == :recommendations
+ print_subheader("Generated Preprocessing Recommendations")
+ print_recommendations(results)
+ else
+ print_subheader("Phase: $(String(phase))")
+ if isa(results, Dict)
+ for (key, value) in results
+ print_param(key, value)
+ end
+ elseif isa(results, NamedTuple)
+ for field in fieldnames(typeof(results))
+ print_param(field, getfield(results, field))
+ end
+ else
+ println(" $(String(phase)) results: $(results)")
+ end
+ end
+ end
+
+ close(msi_data_mzml) # Close file handles
+ catch e
+ println("ERROR processing mzML file: $e")
+ showerror(stdout, e, catch_backtrace())
+ end
+ end
+
+ # --- Process imzML file ---
+ print_header("Processing imzML File: $(basename(TEST_IMZML_FILE))")
+ if !isfile(TEST_IMZML_FILE)
+ println("SKIPPING: imzML file not found at $(TEST_IMZML_FILE)")
+ else
+ try
+ msi_data_imzml = @time OpenMSIData(TEST_IMZML_FILE)
+ check_data_range(msi_data_imzml)
+ analysis_results_imzml = run_preprocessing_analysis(msi_data_imzml, reference_peaks=reference_peaks, mask_path=MASK_ROUTE)
+
+ println("\n" * "*"^80)
+ println("IMZML PREPROCESSING ANALYSIS RESULTS")
+ println("*"^80)
+
+ for (phase, results) in analysis_results_imzml
+ if phase == :recommendations
+ print_subheader("Generated Preprocessing Recommendations")
+ print_recommendations(results)
+ else
+ print_subheader("Phase: $(String(phase))")
+ if isa(results, Dict)
+ for (key, value) in results
+ print_param(key, value)
+ end
+ elseif isa(results, NamedTuple)
+ for field in fieldnames(typeof(results))
+ print_param(field, getfield(results, field))
+ end
+ else
+ println(" $(String(phase)) results: $(results)")
+ end
+ end
+ end
+
+ close(msi_data_imzml) # Close file handles
+ catch e
+ println("ERROR processing imzML file: $e")
+ showerror(stdout, e, catch_backtrace())
+ end
+ end
+end
+
+# --- Execute ---
+@time run_precalculation_example()
+
+# =============================================================================
+# Example Usage
+# =============================================================================
+
+#=
+"""
+ example_preanalysis_workflow(msi_data::MSIData)
+
+Demonstrates how to run the pre-analysis pipeline.
+"""
+function example_preanalysis_workflow(msi_data::MSIData)
+ # Load your MSIData object (this would come from your actual data loading)
+ # msi_data = load_imzml_dataset("path/to/your/data.imzML") # Assuming msi_data is passed
+
+ # Define reference peaks for mass accuracy analysis
+ reference_peaks = Dict(
+ 89.04767 => "Alanin",
+ 147.07642 => "Lysin",
+ 189.12392 => "Unknown",
+ 524.26496 => "PC(34:1) [M+H]+"
+ )
+
+ # Define region masks if you have spatial annotations
+ region_masks = Dict{Symbol, BitMatrix}()
+ # region_masks[:tumor] = load_mask("tumor_mask.png") # Not defined, comment out
+ # region_masks[:stroma] = load_mask("stroma_mask.png") # Not defined, comment out
+
+ println("Starting comprehensive pre-analysis...")
+
+ # Run the complete pre-analysis pipeline
+ analysis_results = run_preprocessing_analysis(
+ msi_data, # Your MSIData object
+ reference_peaks=reference_peaks,
+ region_masks=region_masks,
+ sample_size=200 # Adjust based on dataset size
+ )
+
+ # Access the recommendations
+ recommendations = analysis_results[:recommendations]
+
+ println("\n" * "="^60)
+ println("PREPROCESSING RECOMMENDATIONS")
+ println("="^60)
+
+ for (step, params) in recommendations
+ println("\n$step:")
+ for (key, value) in params
+ println(" - $key: $value")
+ end
+ end
+
+ return analysis_results
+end
+
+# You can also run individual analysis steps:
+function run_targeted_analysis(msi_data::MSIData)
+ # Just analyze signal quality and peak characteristics
+ signal_analysis = analyze_signal_quality(msi_data, sample_size=100)
+ peak_analysis = analyze_peak_characteristics(msi_data, sample_size=50)
+
+ return (signal_analysis, peak_analysis)
+end
+
+=#
diff --git a/test/run_preprocessing.jl b/test/run_preprocessing.jl
index bc70023..2e0e329 100644
--- a/test/run_preprocessing.jl
+++ b/test/run_preprocessing.jl
@@ -1,371 +1,534 @@
-# test/run_preprocessing.jl
-
-# ===================================================================
-# Preprocessing Test Suite for JuliaMSI
-# ===================================================================
-# This script provides a customizable workflow to test and visualize
-# the effects of different preprocessing steps on individual spectra
-# from .mzML or .imzML files.
-#
-# Instructions:
-# 1. Configure the file paths and parameters in the "CONFIG" section.
-# 2. Run the script from the project's root directory:
-# julia --threads auto --project=. test/run_preprocessing.jl
-# 3. Check the configured `RESULTS_DIR` for output plots and CSV files.
-# ===================================================================
+# run_preprocessing.jl
using Printf
-using CairoMakie
-using DataFrames
-using CSV
-using Statistics
-using Interpolations
import Pkg
+using CairoMakie
+using DataFrames # For creating dataframes
+using CSV
# --- Load the MSI_src Module ---
-# Activate the project environment at the parent directory of this test script
Pkg.activate(joinpath(@__DIR__, ".."))
using MSI_src
# ===================================================================
-# CONFIG: CUSTOMIZE YOUR PREPROCESSING WORKFLOW HERE
+# CONFIG: PLEASE FILL IN YOUR FILE PATHS HERE
# ===================================================================
-# --- Input and Output ---
-const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML" # Can be .mzML or .imzML
-const RESULTS_DIR = "test/results/preprocessing"
-const NUM_SPECTRA_TO_PROCESS = nothing # Set to `nothing` to process all spectra
+const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Thricoderma_etc/Imaging_interaccion_trichoderma_vs_streptomyces.imzML"
+# const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/set de datos MS/Atropina_tuneo_fraq_20ev.mzML"
+# const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML"
-# --- Internal Standards for Calibration and QC ---
-# replace these with m/z values of known compounds present in your dataset.
-const INTERNAL_STANDARDS = Dict(
- "P13" => 432.6584,
- "P15" => 464.6059,
- "P17" => 526.5534,
- "P21" => 650.4485,
- "P25" => 774.3435,
- "P29" => 898.2385,
- "P31" => 950.1861,
- "P33" => 1022.1336,
- "P37" => 1146.0286,
- "P45" => 1593.8187,
- "Unknown 1" => 772.433,
- "Unknown 2" => 772.5253
-)
+# const MASK_ROUTE = "/home/pixel/Documents/Cinvestav_2025/JuliaMSI/public/css/masks/Stomach_DHB_uncompressed.png"
+const MASK_ROUTE = ""
-# --- Preprocessing Step Parameters ---
+const OUTPUT_DIR = "./test/results/preprocessing_results"
-# Step 0: Quality Control (QC)
-const QC_PARAMS = (
- ppm_tolerance = 5.0, # PPM tolerance for matching internal standards
-)
-
-# Step 1: Calibration
-const CALIBRATION_PARAMS = (
- enabled = true,
- ppm_tolerance = 5.0, # PPM tolerance for finding calibration peaks
-)
-
-# Step 2: Smoothing
-# Note: Smoothing is less effective and often unnecessary for centroid data.
-const SMOOTHING_PARAMS = (
- enabled = false,
- window = 9,
- order = 2,
-)
-
-# Step 3: Baseline Correction
-# Note: Baseline correction is less effective and often unnecessary for centroid data.
-const BASELINE_CORRECTION_PARAMS = (
- enabled = false,
- iterations = 100,
-)
-
-# Step 4: Normalization
-const NORMALIZATION_PARAMS = (
- enabled = true,
- method = :tic, # :tic, :median, or :none
-)
-
-# Step 5: Peak Detection
-const PEAK_DETECTION_PARAMS = (
- enabled = true,
- method = :centroid, # Use :centroid for centroided data, :profile for profile data
- # --- Parameters for :centroid method ---
- intensity_threshold = 0.0, # Filters out peaks below this absolute intensity
- # --- Parameters for :profile method ---
- half_window = 10,
- snr_threshold = 3.0,
- min_peak_prominence = 0.05,
- merge_peaks_tolerance = 0.002,
-)
-
-# Step 6: Peak Alignment (Warping)
-# This is performed after collecting peaks from all spectra.
-const PEAK_ALIGNMENT_PARAMS = (
- enabled = true,
- tolerance = 0.002,
- tolerance_unit = :mz, # :mz or :ppm
- min_matched_peaks = 3, # Lowered for sparse centroid data
-)
-
-# Step 7: Peak Binning
-const PEAK_BINNING_PARAMS = (
- enabled = true,
- tolerance = 0.1,
- tolerance_unit = :mz, # :mz or :ppm
- frequency_threshold = 0.1, # Min fraction of spectra a peak must be in
+# Reference peaks for calibration and alignment
+const reference_peaks = Dict(
+ 137.0244 => "DHB_fragment",
+ 155.0349 => "DHB_M+H",
+ 177.0168 => "DHB_M+Na",
+ 496.3398 => "PC_16:0_16:0",
+ 520.3398 => "PC_16:0_18:1",
+ 760.5851 => "PC_16:0_18:1_Na",
+ 391.2843 => "PDMS",
+ 413.2662 => "PDMS_Na",
+ 842.5092 => "Protein_standard",
+ 1045.532 => "Protein_standard",
+ 290.1747 => "Atropine [M+H]+",
+ 304.1903 => "Scopolamine [M+H]+",
+ 124.0393 => "Tropine [M+H]+",
)
# ===================================================================
# HELPER FUNCTIONS
# ===================================================================
-"""
- plot_spectrum_on_fig(fig_pos, mz, intensity, title)
-
-Helper function to plot a spectrum on a specific position of a CairoMakie figure.
-"""
-function plot_spectrum_on_fig(fig_pos, mz, intensity, title)
- ax = Axis(fig_pos, title=title, xlabel="m/z", ylabel="Intensity")
- lines!(ax, mz, intensity)
+function ensure_output_dir()
+ if !isdir(OUTPUT_DIR)
+ mkdir(OUTPUT_DIR)
+ println("Created output directory: $OUTPUT_DIR")
+ end
end
+function plot_spectrum_step(mz, intensity, step_name; peaks=nothing, spectrum_index=1)
+ fig = Figure(size=(1200, 600))
+ ax = Axis(fig[1, 1], xlabel="m/z", ylabel="Intensity", title="Step: $step_name (Spectrum $spectrum_index)")
+
+ lines!(ax, mz, intensity, color=:blue, linewidth=1)
+
+ if peaks !== nothing && !isempty(peaks)
+ peak_mz = [p.mz for p in peaks]
+ peak_intensity = [p.intensity for p in peaks]
+ scatter!(ax, peak_mz, peak_intensity, color=:red, markersize=5)
+ end
+
+ filename = joinpath(OUTPUT_DIR, "spectrum_$(step_name)_index_$spectrum_index.png")
+ save(filename, fig)
+ println(" - Saved plot: $(basename(filename))")
+end
+
+function print_step_header(step_name::String)
+ println("\n" * ">"^50)
+ println("APPLYING STEP: $step_name")
+ println(">"^50)
+end
# ===================================================================
-# MAIN PROCESSING SCRIPT
+# PREPROCESSING PIPELINE FUNCTIONS (IN-PLACE)
# ===================================================================
-function run_preprocessing_suite()
- println("Starting Preprocessing Test Suite...")
- mkpath(RESULTS_DIR)
+function apply_baseline_correction(spectra::Vector{MutableSpectrum}, params::Dict, msi_data::MSIData)
+ print_step_header("Baseline Correction")
+
+ method = get(params, :method, :snip)
+ iterations = get(params, :iterations, 100)
+ window = get(params, :window, 20)
+ println(" Method: $method, Iterations: $iterations, Window: $window")
- # --- Load Data ---
- println("Loading data from: $TEST_FILE")
- if !isfile(TEST_FILE)
- println("ERROR: Test file not found. Please check the path in the CONFIG section.")
- return
- end
- msi_data = @time OpenMSIData(TEST_FILE)
- println("Data loaded successfully. Found $(length(msi_data.spectra_metadata)) spectra.")
-
- # --- Determine which spectra to process ---
- total_spectra = length(msi_data.spectra_metadata)
- indices_to_process = if NUM_SPECTRA_TO_PROCESS === nothing
- 1:total_spectra
- else
- unique(round.(Int, range(1, total_spectra, length=min(NUM_SPECTRA_TO_PROCESS, total_spectra))))
- end
- println("Will process $(length(indices_to_process)) spectra.")
-
- # --- Data storage for results ---
- qc_results = DataFrame(spectrum_idx=Int[], metric=String[], value=Float64[], compound=String[])
- all_processed_peaks_mz = Vector{Vector{Float64}}()
- all_processed_peaks_int = Vector{Vector{Float64}}()
- spectrum_indices_with_peaks = Int[]
-
- # --- Main Loop: Process each spectrum individually ---
- println("\n" * "="^20 * " Processing Individual Spectra " * "="^20)
- for (i, idx) in enumerate(indices_to_process)
- print("\rProcessing spectrum #$idx ($(i)/$(length(indices_to_process)))...")
-
- process_spectrum(msi_data, idx) do mz, intensity
- if qc_is_empty(mz, intensity) || !qc_is_regular(mz)
- # @warn "Skipping empty or irregular spectrum #$idx"
- return
- end
-
- original_mz, original_intensity = copy(mz), copy(intensity)
- processed_mz, processed_intensity = copy(mz), copy(intensity)
-
- # --- Step 0: QC (PPM and Resolution) ---
- # This happens before any modification to the m/z axis
- let ref_masses = collect(values(INTERNAL_STANDARDS))
- matched_peaks = find_calibration_peaks(processed_mz, processed_intensity, ref_masses, ppm_tolerance=QC_PARAMS.ppm_tolerance)
-
- for (compound, theoretical_mz) in INTERNAL_STANDARDS
- # Find the measured m/z that corresponds to this theoretical_mz
- measured_mz = 0.0
- for (theo, meas) in matched_peaks
- if theo == theoretical_mz
- measured_mz = meas
- break
- end
- end
-
- if measured_mz > 0
- # PPM Error
- ppm_error = calculate_ppm_error(measured_mz, theoretical_mz)
- push!(qc_results, (idx, "ppm_error", ppm_error, compound))
-
- # Resolution
- resolution = calculate_resolution_fwhm(measured_mz, processed_mz, processed_intensity)
- if !isnan(resolution)
- push!(qc_results, (idx, "resolution", resolution, compound))
- end
- end
- end
- end
-
- # --- Step 1: Calibration ---
- if CALIBRATION_PARAMS.enabled
- ref_masses = collect(values(INTERNAL_STANDARDS))
- matched_peaks = find_calibration_peaks(processed_mz, processed_intensity, ref_masses, ppm_tolerance=CALIBRATION_PARAMS.ppm_tolerance)
- if length(matched_peaks) >= 2
- measured = sort(collect(values(matched_peaks)))
- theoretical = sort(collect(keys(matched_peaks)))
- itp = linear_interpolation(measured, theoretical, extrapolation_bc=Line())
- processed_mz = itp(processed_mz)
- end
- end
-
- # --- Step 2: Smoothing ---
- if SMOOTHING_PARAMS.enabled
- processed_intensity = smooth_spectrum(processed_intensity, window=SMOOTHING_PARAMS.window, order=SMOOTHING_PARAMS.order)
- end
-
- # --- Step 3: Baseline Correction ---
- if BASELINE_CORRECTION_PARAMS.enabled
- baseline = snip_baseline(processed_intensity, iterations=BASELINE_CORRECTION_PARAMS.iterations)
- processed_intensity .-= baseline
- processed_intensity = max.(0, processed_intensity) # Ensure non-negativity
- end
-
- # --- Step 4: Normalization ---
- if NORMALIZATION_PARAMS.enabled && NORMALIZATION_PARAMS.method != :none
- if NORMALIZATION_PARAMS.method == :tic
- processed_intensity = tic_normalize(processed_intensity)
- elseif NORMALIZATION_PARAMS.method == :median
- processed_intensity = median_normalize(processed_intensity)
- end
- end
-
- # --- Step 5: Peak Detection ---
- local pk_mz, pk_int
- pk_mz, pk_int = Float64[], Float64[] # Initialize empty
- if PEAK_DETECTION_PARAMS.enabled
- if PEAK_DETECTION_PARAMS.method == :profile
- pk_mz, pk_int = detect_peaks_profile(processed_mz, processed_intensity,
- half_window=PEAK_DETECTION_PARAMS.half_window,
- snr_threshold=PEAK_DETECTION_PARAMS.snr_threshold,
- min_peak_prominence=PEAK_DETECTION_PARAMS.min_peak_prominence,
- merge_peaks_tolerance=PEAK_DETECTION_PARAMS.merge_peaks_tolerance)
- elseif PEAK_DETECTION_PARAMS.method == :wavelet
- pk_mz, pk_int = detect_peaks_wavelet(processed_mz, processed_intensity)
- else # :centroid
- pk_mz, pk_int = detect_peaks_centroid(processed_mz, processed_intensity)
- end
-
- if !isempty(pk_mz)
- push!(all_processed_peaks_mz, pk_mz)
- push!(all_processed_peaks_int, pk_int)
- push!(spectrum_indices_with_peaks, idx)
- end
- end
-
- # --- Visualization of a sample spectrum ---
- if i == 1 # Only plot the first processed spectrum
- println("\nGenerating example plots for spectrum #$idx...")
- fig = Figure(size=(1200, 800))
-
- plot_spectrum_on_fig(fig[1,1], original_mz, original_intensity, "1. Original Spectrum")
- plot_spectrum_on_fig(fig[2,1], processed_mz, processed_intensity, "2. After All Steps (Before Peak Picking)")
-
- # Plot detected peaks
- ax = Axis(fig[3,1], title="3. Detected Peaks")
- if !isempty(pk_mz)
- stem!(ax, pk_mz, pk_int)
- end
-
- save(joinpath(RESULTS_DIR, "example_spectrum_processing.png"), fig)
- println("Saved example processing plots.")
- end
+ # Safely plot the first spectrum
+ if !isempty(spectra)
+ s = spectra[1]
+ if validate_spectrum(s.mz, s.intensity)
+ baseline = MSI_src.apply_baseline_correction(s.intensity; method=method, iterations=iterations, window=window)
+ corrected_intensity = max.(0.0, s.intensity .- baseline)
+ plot_spectrum_step(s.mz, corrected_intensity, "baseline_correction")
end
end
- println("\nIndividual spectrum processing complete.")
- # --- Save QC Results ---
- if !isempty(qc_results)
- println("\n" * "="^20 * " Generating QC Report " * "="^20)
- CSV.write(joinpath(RESULTS_DIR, "qc_results.csv"), qc_results)
- println("QC results saved to qc_results.csv")
-
- # Generate summary plots for QC
- fig = Figure(size=(1200, 600))
-
- # PPM Error Histogram
- ppm_errors = filter(row -> row.metric == "ppm_error", qc_results).value
- if !isempty(ppm_errors)
- ax1 = Axis(fig[1,1], title="PPM Error Distribution", xlabel="PPM Error")
- hist!(ax1, ppm_errors, bins=30)
+ Threads.@threads for s in spectra
+ if validate_spectrum(s.mz, s.intensity)
+ baseline = MSI_src.apply_baseline_correction(s.intensity; method=method, iterations=iterations, window=window)
+ s.intensity = max.(0.0, s.intensity .- baseline)
end
+ end
+end
- # Resolution Histogram
- resolutions = filter(row -> row.metric == "resolution", qc_results).value
- if !isempty(resolutions)
- ax2 = Axis(fig[1,2], title="Resolution Distribution", xlabel="Resolution (FWHM)")
- hist!(ax2, resolutions, bins=30)
+function apply_smoothing(spectra::Vector{MutableSpectrum}, params::Dict)
+ print_step_header("Smoothing")
+
+ method = get(params, :method, :savitzky_golay)
+ window = get(params, :window, 9)
+ order = get(params, :order, 2)
+ println(" Method: $method, Window: $window, Order: $order")
+
+ # Safely plot the first spectrum
+ if !isempty(spectra)
+ s = spectra[1]
+ if validate_spectrum(s.mz, s.intensity)
+ smoothed_intensity = max.(0.0, smooth_spectrum(s.intensity; method=method, window=window, order=order))
+ plot_spectrum_step(s.mz, smoothed_intensity, "smoothing", spectrum_index=s.id)
end
-
- save(joinpath(RESULTS_DIR, "qc_summary_plots.png"), fig)
- println("QC summary plots saved.")
end
- if isempty(all_processed_peaks_mz)
- @warn "No peaks were detected in any of the processed spectra. Skipping alignment and binning."
- return
- end
-
- # --- Step 6: Peak Alignment ---
- println("\n" * "="^20 * " Aligning Peaks " * "="^20)
- aligned_peaks_mz = copy(all_processed_peaks_mz)
- if PEAK_ALIGNMENT_PARAMS.enabled
- # Create a reference peak list (e.g., from the spectrum with the most peaks)
- ref_idx = argmax(length.(all_processed_peaks_mz))
- reference_peaks = all_processed_peaks_mz[ref_idx]
- println("Using spectrum $(spectrum_indices_with_peaks[ref_idx]) as alignment reference.")
-
- for i in 1:length(aligned_peaks_mz)
- if i == ref_idx continue end
- alignment_func = align_peaks_lowess(reference_peaks, aligned_peaks_mz[i],
- tolerance=PEAK_ALIGNMENT_PARAMS.tolerance,
- tolerance_unit=PEAK_ALIGNMENT_PARAMS.tolerance_unit,
- min_matched_peaks=PEAK_ALIGNMENT_PARAMS.min_matched_peaks)
-
- aligned_peaks_mz[i] = alignment_func(aligned_peaks_mz[i])
+ Threads.@threads for s in spectra
+ if validate_spectrum(s.mz, s.intensity)
+ smoothed_intensity = max.(0.0, smooth_spectrum(s.intensity; method=method, window=window, order=order))
+ s.intensity = smoothed_intensity
end
- println("Peak alignment complete.")
end
+end
- # --- Step 7: Peak Binning & Feature Matrix Generation ---
- println("\n" * "="^20 * " Binning Peaks " * "="^20)
- if PEAK_BINNING_PARAMS.enabled
- feature_matrix, mz_bins = bin_peaks(aligned_peaks_mz, all_processed_peaks_int,
- PEAK_BINNING_PARAMS.tolerance,
- tolerance_unit=PEAK_BINNING_PARAMS.tolerance_unit,
- frequency_threshold=PEAK_BINNING_PARAMS.frequency_threshold)
-
- if !isempty(feature_matrix)
- df = DataFrame(feature_matrix, :auto)
- rename!(df, ["bin_$(i)" for i in 1:size(df, 2)])
- insertcols!(df, 1, :spectrum_idx => spectrum_indices_with_peaks)
-
- CSV.write(joinpath(RESULTS_DIR, "feature_matrix.csv"), df)
- println("Feature matrix saved to feature_matrix.csv")
-
- # Save bin m/z ranges
- bin_df = DataFrame(bin_index=1:length(mz_bins), mz_start=[b[1] for b in mz_bins], mz_end=[b[2] for b in mz_bins])
- CSV.write(joinpath(RESULTS_DIR, "bin_definitions.csv"), bin_df)
- println("Bin definitions saved to bin_definitions.csv")
+function apply_peak_picking(spectra::Vector{MutableSpectrum}, params::Dict)
+ print_step_header("Peak Picking")
+
+ method = get(params, :method, :profile)
+ snr_threshold = get(params, :snr_threshold, 3.0)
+ half_window = get(params, :half_window, 10)
+ min_peak_prominence = get(params, :min_peak_prominence, 0.1)
+ merge_peaks_tolerance = get(params, :merge_peaks_tolerance, 0.002)
+ println(" Method: $method, SNR: $snr_threshold, Half Window: $half_window")
+
+ found_invalid_spectrum = Threads.Atomic{Bool}(false)
+
+ Threads.@threads for s in spectra
+ is_valid = validate_spectrum(s.mz, s.intensity)
+
+ if !is_valid && !found_invalid_spectrum[]
+ if Threads.atomic_xchg!(found_invalid_spectrum, true) == false
+ # Debug printing for the first invalid spectrum found
+ end
+ end
+
+ if is_valid
+ if method == :profile
+ s.peaks = detect_peaks_profile(s.mz, s.intensity; snr_threshold=snr_threshold, half_window=half_window, min_peak_prominence=min_peak_prominence, merge_peaks_tolerance=merge_peaks_tolerance)
+ elseif method == :wavelet
+ s.peaks = detect_peaks_wavelet(s.mz, s.intensity; snr_threshold=snr_threshold, half_window=half_window)
+ elseif method == :centroid
+ s.peaks = detect_peaks_centroid(s.mz, s.intensity; snr_threshold=snr_threshold)
+ else
+ s.peaks = detect_peaks_profile(s.mz, s.intensity; snr_threshold=snr_threshold, half_window=half_window)
+ end
else
- @warn "Feature matrix was empty after binning."
+ s.peaks = []
end
end
- println("\nPreprocessing Test Suite finished successfully!")
+ # Safely plot the first spectrum
+ if !isempty(spectra)
+ s1 = findfirst(s -> s.id == 1, spectra)
+ if s1 !== nothing
+ plot_spectrum_step(spectra[s1].mz, spectra[s1].intensity, "peak_picking", peaks=spectra[s1].peaks)
+ println(" - Detected $(length(spectra[s1].peaks)) peaks in spectrum 1")
+ end
+ end
end
+function apply_calibration(spectra::Vector{MutableSpectrum}, params::Dict, reference_peaks::Dict)
+ print_step_header("Calibration")
+
+ method = get(params, :method, :none)
+ ppm_tolerance = get(params, :ppm_tolerance, 20.0)
+ println(" Method: $method, PPM Tolerance: $ppm_tolerance")
-# --- Execute ---
-run_preprocessing_suite()
+ if method == :none || isempty(reference_peaks)
+ println(" - Skipping calibration (no method or reference peaks)")
+ return
+ end
+
+ reference_masses = collect(keys(reference_peaks))
+ calibration_info = Vector{String}(undef, length(spectra))
+
+ Threads.@threads for i in 1:length(spectra)
+ s = spectra[i]
+ info_message = ""
+ if validate_spectrum(s.mz, s.intensity)
+ matched_peaks = find_calibration_peaks(s.mz, s.intensity, reference_masses; ppm_tolerance=ppm_tolerance)
+ if length(matched_peaks) >= 2
+ measured = sort(collect(values(matched_peaks)))
+ theoretical = sort(collect(keys(matched_peaks)))
+ itp = linear_interpolation(measured, theoretical, extrapolation_bc=Line())
+ s.mz = itp(s.mz) # Modify mz-axis in-place
+ info_message = " - Spectrum $(s.id): calibrated using $(length(matched_peaks)) reference peaks"
+ else
+ info_message = " - Spectrum $(s.id): insufficient reference peaks ($(length(matched_peaks)) found), skipping"
+ end
+ end
+ calibration_info[i] = info_message
+ end
+
+ # Print summary of results
+ printed_count = 0
+ for info in calibration_info
+ if !isempty(info) && printed_count < 5
+ println(info)
+ printed_count += 1
+ end
+ end
+end
+
+function apply_peak_alignment(spectra::Vector{MutableSpectrum}, params::Dict)
+ print_step_header("Peak Alignment")
+
+ method = get(params, :method, :none)
+ tolerance = get(params, :tolerance, 0.002)
+ tolerance_unit = get(params, :tolerance_unit, :mz)
+ println(" Method: $method, Tolerance: $tolerance $tolerance_unit")
+
+ if method == :none
+ println(" - Skipping peak alignment")
+ return
+ end
+
+ ref_find_idx = findfirst(s -> !isempty(s.peaks), spectra)
+ if ref_find_idx === nothing
+ println(" - Insufficient spectra with peaks for alignment. Skipping.")
+ return
+ end
+
+ ref_spectrum = spectra[ref_find_idx]
+ ref_peaks_mz = [p.mz for p in ref_spectrum.peaks]
+ println(" - Using spectrum $(ref_spectrum.id) as reference with $(length(ref_peaks_mz)) peaks")
+
+ Threads.@threads for s in spectra
+ if s.id == ref_spectrum.id || isempty(s.peaks)
+ continue
+ end
+
+ current_peaks_mz = [p.mz for p in s.peaks]
+ alignment_func = align_peaks_lowess(ref_peaks_mz, current_peaks_mz; method=method, tolerance=tolerance, tolerance_unit=tolerance_unit)
+
+ s.mz = alignment_func.(s.mz) # Update m/z axis
+
+ # Update peak m/z values
+ for i in 1:length(s.peaks)
+ old_peak = s.peaks[i]
+ aligned_peak_mz = alignment_func(old_peak.mz)
+ s.peaks[i] = (mz=aligned_peak_mz, intensity=old_peak.intensity, fwhm=old_peak.fwhm, shape_r2=old_peak.shape_r2, snr=old_peak.snr, prominence=old_peak.prominence)
+ end
+ end
+end
+
+function apply_normalization(spectra::Vector{MutableSpectrum}, params::Dict)
+ print_step_header("Normalization")
+
+ method = get(params, :method, :tic)
+ println(" Method: $method")
+
+ # Safely plot the first spectrum
+ if !isempty(spectra)
+ s1 = findfirst(s -> s.id == 1, spectra)
+ if s1 !== nothing
+ s = spectra[s1]
+ if validate_spectrum(s.mz, s.intensity)
+ normalized_intensity = MSI_src.apply_normalization(s.intensity; method=method)
+ plot_spectrum_step(s.mz, normalized_intensity, "normalization")
+ end
+ end
+ end
+
+ Threads.@threads for s in spectra
+ if validate_spectrum(s.mz, s.intensity)
+ s.intensity = MSI_src.apply_normalization(s.intensity; method=method)
+ end
+ end
+end
+
+function apply_peak_binning(spectra::Vector{MutableSpectrum}, params::Dict)
+ print_step_header("Peak Binning")
+
+ method = get(params, :method, :adaptive)
+ tolerance = get(params, :tolerance, 20.0)
+ tolerance_unit = get(params, :tolerance_unit, :ppm)
+ min_peak_per_bin = get(params, :min_peak_per_bin, 3)
+ max_bin_width_ppm = get(params, :max_bin_width_ppm, 150.0)
+ intensity_weighted_centers = get(params, :intensity_weighted_centers, true)
+ println(" Method: $method, Tolerance: $tolerance $tolerance_unit")
+
+ if isempty(spectra) || all(s -> isempty(s.peaks), spectra)
+ @warn "No peaks found for binning. Returning empty feature matrix."
+ return nothing, nothing
+ end
+
+ println(" - Binning peaks from $(length(spectra)) spectra")
+
+ binning_params = PeakBinningParams(
+ method=method,
+ tolerance=tolerance,
+ tolerance_unit=tolerance_unit,
+ min_peak_per_bin=min_peak_per_bin,
+ max_bin_width_ppm=max_bin_width_ppm,
+ intensity_weighted_centers=intensity_weighted_centers
+ )
+
+ feature_matrix, bin_definitions = bin_peaks(spectra, binning_params)
+
+ if feature_matrix !== nothing && bin_definitions !== nothing
+ println(" - Generated feature matrix: $(size(feature_matrix.matrix))")
+ println(" - Number of bins: $(length(bin_definitions))")
+ end
+
+ return feature_matrix, bin_definitions
+end
+
+function save_feature_matrix(feature_matrix, bin_definitions)
+ print_step_header("Saving Results")
+
+ # Save feature matrix as CSV
+ csv_path = joinpath(OUTPUT_DIR, "feature_matrix.csv")
+
+ # Create column headers
+ bin_headers = ["bin_$(i)_$(round(def[1], digits=4))-$(round(def[2], digits=4))"
+ for (i, def) in enumerate(bin_definitions)]
+
+ # Open file for writing
+ open(csv_path, "w") do io
+ # Write header row
+ write(io, "spectrum_index," * join(bin_headers, ",") * "\n")
+
+ # Write data rows
+ for r_idx in 1:size(feature_matrix.matrix, 1)
+ write(io, "$(feature_matrix.sample_ids[r_idx]),")
+ for c_idx in 1:size(feature_matrix.matrix, 2)
+ write(io, "$(feature_matrix.matrix[r_idx, c_idx])")
+ if c_idx < size(feature_matrix.matrix, 2)
+ write(io, ",")
+ end
+ end
+ write(io, "\n")
+ end
+ end
+ println(" - Saved feature matrix: $csv_path")
+
+ # Save bin definitions (still using DataFrame as it's small)
+ bins_path = joinpath(OUTPUT_DIR, "bin_definitions.csv")
+ bins_df = DataFrame(
+ bin_index = 1:length(bin_definitions),
+ mz_start = [def[1] for def in bin_definitions],
+ mz_end = [def[2] for def in bin_definitions],
+ mz_center = [(def[1] + def[2])/2 for def in bin_definitions]
+ )
+ CSV.write(bins_path, bins_df)
+ println(" - Saved bin definitions: $bins_path")
+
+ return csv_path, bins_path
+end
+
+# ===================================================================
+# USER OVERRIDES: Manually specify parameters here
+# ===================================================================
+# This dictionary allows you to override any auto-detected parameters.
+# The structure should match the output of `main_precalculation`.
+# Example: Force a less aggressive peak prominence threshold.
+const USER_OVERRIDES = Dict(
+ :PeakPicking => Dict(
+ #:snr_threshold => 1.5,
+ :half_window => 5,
+ :snr_threshold => 15.0,
+ #:merge_peaks_tolerance => 2.5,
+ #:half_window => 2
+ ),
+ # :Smoothing => Dict(
+ # :window => 11
+ # )
+ #:PeakAlignment => Dict(
+ # :tolerance => 0.01
+ #)
+ :PeakBinningParams => Dict(
+ :tolerance => 30.0
+ )
+)
+
+# ===================================================================
+# MAIN PREPROCESSING PIPELINE
+# ===================================================================
+
+function run_preprocessing_pipeline()
+ println("="^80)
+ println("MSI PREPROCESSING PIPELINE")
+ println("="^80)
+
+ # Ensure output directory exists
+ ensure_output_dir()
+
+ # Load data
+ println("\nLoading data: $(basename(TEST_FILE))")
+ msi_data = OpenMSIData(TEST_FILE)
+
+ # Precompute analytics
+ println("\nPrecomputing analytics...")
+ precompute_analytics(msi_data)
+
+ # Run precalculation to get automatic parameters
+ println("\nRunning precalculation for automatic parameter determination...")
+ auto_params = main_precalculation(msi_data, reference_peaks=reference_peaks, mask_path=MASK_ROUTE)
+
+ # --- Apply User Overrides ---
+ if !isempty(USER_OVERRIDES)
+ println("\nApplying user overrides...")
+ for (step, params) in USER_OVERRIDES
+ if haskey(auto_params, step)
+ for (param, value) in params
+ @info "OVERRIDE: For step `:$step`, setting `:$param` to `$value` (was `$(get(auto_params[step], param, "not set"))`)"
+ end
+ merge!(auto_params[step], params)
+ else
+ @warn "User override for non-existent step `:$step` ignored."
+ end
+ end
+ end
+
+ # --- Final Parameter Summary ---
+ println("\n" * "-"^60)
+ println("FINAL PARAMETERS FOR PIPELINE RUN")
+ println("-"^60)
+ for (step, params) in sort(collect(pairs(auto_params)), by=p -> p.first)
+ println(" - Step: $step")
+ if isempty(params)
+ println(" (No parameters)")
+ continue
+ end
+ for (param, value) in sort(collect(pairs(params)), by=p -> p.first)
+ println(" - $(rpad(param, 25)): $value")
+ end
+ end
+
+ # Define pipeline steps
+ pipeline_steps = [
+ "baseline_correction",
+ "smoothing",
+ "peak_picking",
+ "calibration",
+ "peak_alignment",
+ "normalization",
+ "peak_binning"
+ ]
+
+ println("\nPipeline steps: $(join(pipeline_steps, " -> "))")
+
+ # Initialize `current_spectra` as a Vector of mutable structs for in-place modification
+ println("\nInitializing spectra data structure...")
+ num_spectra = length(msi_data.spectra_metadata)
+ current_spectra = Vector{MutableSpectrum}(undef, num_spectra)
+ _iterate_spectra_fast(msi_data) do idx, mz, intensity
+ current_spectra[idx] = MutableSpectrum(idx, mz, intensity, [])
+
+ # Plot raw spectrum without masks
+ if idx == 1
+ plot_spectrum_step(mz, intensity, "raw_unmasked_spectrum")
+ end
+ end
+
+ # These variables will be populated by the pipeline steps
+ feature_matrix = nothing
+ bin_definitions = nothing
+
+ # Apply pipeline steps, modifying `current_spectra` in-place
+ for step in pipeline_steps
+ println("\n" * "-"^60)
+ println("PROCESSING STEP: $step")
+ println("-"^60)
+
+ if step == "baseline_correction"
+ @time apply_baseline_correction(current_spectra, auto_params[:BaselineCorrection], msi_data)
+
+ elseif step == "smoothing"
+ apply_smoothing(current_spectra, auto_params[:Smoothing])
+
+ elseif step == "peak_picking"
+ @time apply_peak_picking(current_spectra, auto_params[:PeakPicking])
+
+ elseif step == "calibration"
+ @time apply_calibration(current_spectra, auto_params[:Calibration], reference_peaks)
+
+ elseif step == "peak_alignment"
+ @time apply_peak_alignment(current_spectra, auto_params[:PeakAlignment])
+
+ elseif step == "normalization"
+ @time apply_normalization(current_spectra, auto_params[:Normalization])
+
+ elseif step == "peak_binning"
+ # This step is different as it generates the final matrix, not modifying spectra in-place
+ feature_matrix, bin_definitions = @time apply_peak_binning(current_spectra, auto_params[:PeakBinningParams])
+
+ if feature_matrix !== nothing
+ @time save_feature_matrix(feature_matrix, bin_definitions)
+ end
+
+ else
+ @warn "Unknown step: $step, skipping"
+ end
+
+ # Print progress
+ println("✓ Completed step: $step")
+ end
+
+ # Clean up
+ close(msi_data)
+ GC.gc()
+
+ println("\n" * "="^80)
+ println("PREPROCESSING PIPELINE COMPLETED SUCCESSFULLY!")
+ println("Results saved to: $OUTPUT_DIR")
+ println("="^80)
+end
+
+# ===================================================================
+# EXECUTE PIPELINE
+# ===================================================================
+
+if abspath(PROGRAM_FILE) == @__FILE__
+ @time run_preprocessing_pipeline()
+end
\ No newline at end of file
diff --git a/test/run_tests.jl b/test/run_tests.jl
index 9dba6b8..6de7a38 100644
--- a/test/run_tests.jl
+++ b/test/run_tests.jl
@@ -75,9 +75,9 @@ const COORDS_TO_PLOT = (50, 50) # Example coordinates (X, Y)
# --- Output Directory ---
const RESULTS_DIR = "test/results"
-test1 = false
-test2 = true
-test3 = false
+test1 = true
+test2 = false
+test3 = true
# ===================================================================
# DATA VALIDATION UTILITY