From bdd931579695664a4fc3f94a5a2b93c3db5cf854 Mon Sep 17 00:00:00 2001 From: Pixelguy14 Date: Tue, 25 Nov 2025 12:53:57 -0600 Subject: [PATCH] Fixed preprocessing module, implemented precalculations, this code gives a prediction of the values for the parameters in every step for preprocessing, added test suites for both the parameter generation and the preprocessing pipeline, included a descriptive processing that returns an image with what every parameter does in every step, added instrument metadata structure with the parser for each type of file, to retrieve several important parameters to automate precalculations and preprocessing --- Manifest.toml | 2 +- Project.toml | 2 + src/MSIData.jl | 171 ++- src/MSI_src.jl | 69 +- src/MzmlConverter.jl | 54 +- src/Precalculations.jl | 1537 +++++++++++++++++++++++ src/Preprocessing.jl | 1605 +++++++++++++++---------- src/PreprocessingVersioning.jl | 446 ------- src/imzML.jl | 330 ++--- src/mzML.jl | 148 ++- test/run_descriptive_preprocessing.jl | 222 ++++ test/run_precalculation_dict.jl | 157 +++ test/run_precalculation_example.jl | 281 +++++ test/run_preprocessing.jl | 813 ++++++++----- test/run_tests.jl | 6 +- 15 files changed, 4167 insertions(+), 1676 deletions(-) create mode 100644 src/Precalculations.jl delete mode 100644 src/PreprocessingVersioning.jl create mode 100644 test/run_descriptive_preprocessing.jl create mode 100644 test/run_precalculation_dict.jl create mode 100644 test/run_precalculation_example.jl 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