diff --git a/app.jl b/app.jl index 8dd47c6..79525d7 100644 --- a/app.jl +++ b/app.jl @@ -2,7 +2,7 @@ module App # ==Packages == -using GenieFramework # Set up Genie development environment. +using GenieFramework using Pkg using Libz using PlotlyBase @@ -22,7 +22,7 @@ using Dates using Base.Threads # Bring MSIData into App module's scope -using .MSI_src: MSIData, OpenMSIData, process_spectrum, IterateSpectra, ImzMLSource, _iterate_spectra_fast, MzMLSource, find_mass, ViridisPalette, get_mz_slice, get_multiple_mz_slices, quantize_intensity, save_bitmap, median_filter, save_bitmap, downsample_spectrum, TrIQ, precompute_analytics, ImportMzmlFile, generate_colorbar_image, load_and_prepare_mask, set_global_mz_range! +using .MSI_src: MSIData, OpenMSIData, process_spectrum, IterateSpectra, ImzMLSource, _iterate_spectra_fast, MzMLSource, find_mass, ViridisPalette, get_mz_slice, get_multiple_mz_slices, quantize_intensity, save_bitmap, median_filter, save_bitmap, downsample_spectrum, TrIQ, precompute_analytics, ImportMzmlFile, generate_colorbar_image, load_and_prepare_mask, set_global_mz_range!, main_precalculation, MutableSpectrum if !@isdefined(increment_image) include("./julia_imzML_visual.jl") @@ -72,6 +72,10 @@ end # == Reactive code == # Reactive code to make the UI interactive @app begin + # == Loading Screen Variables == + @in is_initializing = true + @in initialization_message = "Initializing..." + # == Reactive variables == # reactive variables exist in both the Julia backend and the browser with two-way synchronization # @out variables can only be modified by the backend @@ -207,6 +211,139 @@ end # == Pre Processing Variables == @in pre_tab = "stabilization" + ## Preprocessing Parameters + @in progressPrep=false + @in stabilization_method="sqrt" + @in smoothing_method="sg" + @in smoothing_window = "" + @in smoothing_order = "" + @in baseline_method="snip" + @in baseline_iterations = "" + @in baseline_window = "" + @in normalization_method="tic" + @in alignment_method="lowess" + @in alignment_span = "" + @in alignment_tolerance = "" + @in alignment_tolerance_unit="mz" + @in alignment_max_shift_ppm = "" + @in alignment_min_matched_peaks = "" + @in peak_picking_method="profile" + @in peak_picking_snr_threshold = "" + @in peak_picking_half_window = "" + @in peak_picking_min_peak_prominence = "" + @in peak_picking_merge_peaks_tolerance = "" + @in peak_picking_min_peak_width_ppm = "" + @in peak_picking_max_peak_width_ppm = "" + @in peak_picking_min_peak_shape_r2 = "" + @in binning_method="adaptive" + @in binning_tolerance = "" + @in binning_tolerance_unit="ppm" + @in binning_frequency_threshold = "" + @in binning_min_peak_per_bin = "" + @in binning_max_bin_width_ppm = "" + @in binning_intensity_weighted_centers=true + @in binning_num_uniform_bins = "" + @in calibration_fit_order = "" + @in calibration_ppm_tolerance = "" + @in peak_selection_min_snr = "" + @in peak_selection_min_fwhm_ppm = "" + @in peak_selection_max_fwhm_ppm = "" + @in peak_selection_min_shape_r2 = "" + @in peak_selection_frequency_threshold = "" + @in peak_selection_correlation_threshold = "" + + # == Pipeline Step Order and Mask Route Variables == + @in pipeline_step_order = [ + Dict("name" => "stabilization", "label" => "Stabilization", "enabled" => true), + Dict("name" => "smoothing", "label" => "Smoothing", "enabled" => true), + Dict("name" => "baseline_correction", "label" => "Baseline Correction", "enabled" => true), + Dict("name" => "peak_picking", "label" => "Peak Picking", "enabled" => true), + Dict("name" => "peak_selection", "label" => "Peak Selection", "enabled" => true), + Dict("name" => "calibration", "label" => "Calibration", "enabled" => true), + Dict("name" => "peak_alignment", "label" => "Peak Alignment", "enabled" => true), + Dict("name" => "normalization", "label" => "Normalization", "enabled" => true), + Dict("name" => "peak_binning", "label" => "Peak Binning", "enabled" => true) + ] + @in preprocessing_mask_route = "" + @in selected_spectrum_id_for_plot = 1 + @in feature_matrix_result = nothing + @in bin_info_result = nothing + + @in reference_peaks_list = [ + Dict("mz" => 137.0244, "label" => "DHB_fragment"), + Dict("mz" => 155.0349, "label" => "DHB_M+H"), + ] + + # --- Methods for Reference Peaks List --- + function addReferencePeak() + push!(reference_peaks_list, Dict("mz" => 0.0, "label" => "")) + reference_peaks_list = deepcopy(reference_peaks_list) # Force reactivity + end + + function removeReferencePeak(index::Int) + deleteat!(reference_peaks_list, index) + reference_peaks_list = deepcopy(reference_peaks_list) # Force reactivity + end + + # Step reordering functions + function moveStepUp(index::Int) + if index > 1 + pipeline_step_order = deepcopy(pipeline_step_order) + temp = pipeline_step_order[index] + pipeline_step_order[index] = pipeline_step_order[index-1] + pipeline_step_order[index-1] = temp + pipeline_step_order = pipeline_step_order # Force reactivity + end + end + + function moveStepDown(index::Int) + if index < length(pipeline_step_order) + pipeline_step_order = deepcopy(pipeline_step_order) + temp = pipeline_step_order[index] + pipeline_step_order[index] = pipeline_step_order[index+1] + pipeline_step_order[index+1] = temp + pipeline_step_order = pipeline_step_order # Force reactivity + end + end + @in save_feature_matrix_btn = false + + # Trigger for running the full pipeline + @in run_full_pipeline = false + @out current_pipeline_step = "" # To indicate which step is currently running in the full pipeline + + @in export_params_btn = false + @in import_params_btn = false + @in imported_params_file = nothing + + @in suggested_smoothing_window = "" + @in suggested_smoothing_order = "" + @in suggested_baseline_iterations = "" + @in suggested_baseline_window = "" + @in suggested_alignment_span = "" + @in suggested_alignment_tolerance = "" + @in suggested_alignment_max_shift_ppm = "" + @in suggested_alignment_min_matched_peaks = "" + @in suggested_peak_picking_snr_threshold = "" + @in suggested_peak_picking_half_window = "" + @in suggested_peak_picking_min_peak_prominence = "" + @in suggested_peak_picking_merge_peaks_tolerance = "" + @in suggested_peak_picking_min_peak_width_ppm = "" + @in suggested_peak_picking_max_peak_width_ppm = "" + @in suggested_peak_picking_min_peak_shape_r2 = "" + @in suggested_binning_tolerance = "" + @in suggested_binning_frequency_threshold = "" + @in suggested_binning_min_peak_per_bin = "" + @in suggested_binning_max_bin_width_ppm = "" + @in suggested_binning_num_uniform_bins = "" + @in suggested_calibration_fit_order = "" + @in suggested_calibration_ppm_tolerance = "" + @in suggested_peak_selection_min_snr = "" + @in suggested_peak_selection_min_fwhm_ppm = "" + @in suggested_peak_selection_max_fwhm_ppm = "" + @in suggested_peak_selection_min_shape_r2 = "" + @in suggested_peak_selection_frequency_threshold = "" + @in suggested_peak_selection_correlation_threshold = "" + # == Batch Summary Dialog == @in showBatchSummary = false @out batch_summary = "" @@ -323,11 +460,19 @@ end # Create conection to frontend @out plotdata=[traceSpectra] @out plotlayout=layoutSpectra + + @in idSpectrum=0 @in xCoord=0 @in yCoord=0 @out xSpectraMz = Vector{Float64}() @out ySpectraMz = Vector{Float64}() + # UI plot data for Preprocessing + @out plotdata_before = [traceSpectra] + @out plotlayout_before = layoutSpectra + @out plotdata_after = [traceSpectra] + @out plotlayout_after = layoutSpectra + # Interactive plot reactions @in data_click=Dict{String,Any}() #@in data_selected=Dict{String,Any}() # Selected is for areas, this can work for the masks @@ -413,7 +558,7 @@ end msg = "Opening file: $(basename(picked_route))..." try - dataset_name = replace(basename(picked_route), r"(\.(imzML|imzml|mzML|mzml))$"i => "") + dataset_name = replace(basename(picked_route), r"(\.(imzML|imzml|mzML|mzml))$ "i => "") registry = load_registry(registry_path) existing_entry = get(registry, dataset_name, nothing) @@ -469,6 +614,188 @@ end precompute_analytics(loaded_data) + # Auto-suggest parameters + try + println("Calling main_precalculation to get recommended parameters...") + recommended_params = main_precalculation(loaded_data) + + for (step_name, params) in recommended_params + for (param_key, value) in params + # Convert value to appropriate type before assignment + processed_value = if value === nothing + nothing + elseif value isa Tuple + @warn "Skipping invalid parameter suggestion (tuple): $value for $param_key" + "" # Set to empty string for safety + elseif value isa Number + value + else + string(value) + end + + if processed_value !== nothing + if step_name == :Smoothing + if param_key == :window + suggested_smoothing_window = string(processed_value) + smoothing_window = string(processed_value) + println(" suggested_smoothing_window set to $(suggested_smoothing_window)") + elseif param_key == :order + suggested_smoothing_order = string(processed_value) + smoothing_order = string(processed_value) + println(" suggested_smoothing_order set to $(suggested_smoothing_order)") + end + elseif step_name == :BaselineCorrection + if param_key == :iterations + suggested_baseline_iterations = string(processed_value) + baseline_iterations = string(processed_value) + println(" suggested_baseline_iterations set to $(suggested_baseline_iterations)") + elseif param_key == :window + suggested_baseline_window = string(processed_value) + baseline_window = string(processed_value) + println(" suggested_baseline_window set to $(suggested_baseline_window)") + end + elseif step_name == :PeakAlignment + if param_key == :span + suggested_alignment_span = string(processed_value) + alignment_span = string(processed_value) + println(" suggested_alignment_span set to $(suggested_alignment_span)") + elseif param_key == :tolerance + suggested_alignment_tolerance = string(processed_value) + alignment_tolerance = string(processed_value) + println(" suggested_alignment_tolerance set to $(suggested_alignment_tolerance)") + elseif param_key == :max_shift_ppm + suggested_alignment_max_shift_ppm = string(processed_value) + alignment_max_shift_ppm = string(processed_value) + println(" suggested_alignment_max_shift_ppm set to $(suggested_alignment_max_shift_ppm)") + elseif param_key == :min_matched_peaks + suggested_alignment_min_matched_peaks = string(processed_value) + alignment_min_matched_peaks = string(processed_value) + println(" suggested_alignment_min_matched_peaks set to $(suggested_alignment_min_matched_peaks)") + end + elseif step_name == :Calibration + if param_key == :fit_order + suggested_calibration_fit_order = string(processed_value) + calibration_fit_order = string(processed_value) + println(" suggested_calibration_fit_order set to $(suggested_calibration_fit_order)") + elseif param_key == :ppm_tolerance + suggested_calibration_ppm_tolerance = string(processed_value) + calibration_ppm_tolerance = string(processed_value) + println(" suggested_calibration_ppm_tolerance set to $(suggested_calibration_ppm_tolerance)") + end + elseif step_name == :PeakPicking + if param_key == :snr_threshold + suggested_peak_picking_snr_threshold = string(processed_value) + peak_picking_snr_threshold = string(processed_value) + println(" suggested_peak_picking_snr_threshold set to $(suggested_peak_picking_snr_threshold)") + elseif param_key == :half_window + suggested_peak_picking_half_window = string(processed_value) + peak_picking_half_window = string(processed_value) + println(" suggested_peak_picking_half_window set to $(suggested_peak_picking_half_window)") + elseif param_key == :min_peak_prominence + suggested_peak_picking_min_peak_prominence = string(processed_value) + peak_picking_min_peak_prominence = string(processed_value) + println(" suggested_peak_picking_min_peak_prominence set to $(suggested_peak_picking_min_peak_prominence)") + elseif param_key == :merge_peaks_tolerance + suggested_peak_picking_merge_peaks_tolerance = string(processed_value) + peak_picking_merge_peaks_tolerance = string(processed_value) + println(" suggested_peak_picking_merge_peaks_tolerance set to $(suggested_peak_picking_merge_peaks_tolerance)") + elseif param_key == :min_peak_width_ppm + suggested_peak_picking_min_peak_width_ppm = string(processed_value) + peak_picking_min_peak_width_ppm = string(processed_value) + println(" suggested_peak_picking_min_peak_width_ppm set to $(suggested_peak_picking_min_peak_width_ppm)") + elseif param_key == :max_peak_width_ppm + suggested_peak_picking_max_peak_width_ppm = string(processed_value) + peak_picking_max_peak_width_ppm = string(processed_value) + println(" suggested_peak_picking_max_peak_width_ppm set to $(suggested_peak_picking_max_peak_width_ppm)") + elseif param_key == :min_peak_shape_r2 + suggested_peak_picking_min_peak_shape_r2 = string(processed_value) + peak_picking_min_peak_shape_r2 = string(processed_value) + println(" suggested_peak_picking_min_peak_shape_r2 set to $(suggested_peak_picking_min_peak_shape_r2)") + end + elseif step_name == :PeakSelection + if param_key == :min_snr + suggested_peak_selection_min_snr = string(processed_value) + peak_selection_min_snr = string(processed_value) + println(" suggested_peak_selection_min_snr set to $(suggested_peak_selection_min_snr)") + elseif param_key == :min_fwhm_ppm + suggested_peak_selection_min_fwhm_ppm = string(processed_value) + peak_selection_min_fwhm_ppm = string(processed_value) + println(" suggested_peak_selection_min_fwhm_ppm set to $(suggested_peak_selection_min_fwhm_ppm)") + elseif param_key == :max_fwhm_ppm + suggested_peak_selection_max_fwhm_ppm = string(processed_value) + peak_selection_max_fwhm_ppm = string(processed_value) + println(" suggested_peak_selection_max_fwhm_ppm set to $(suggested_peak_selection_max_fwhm_ppm)") + elseif param_key == :min_shape_r2 + suggested_peak_selection_min_shape_r2 = string(processed_value) + peak_selection_min_shape_r2 = string(processed_value) + println(" suggested_peak_selection_min_shape_r2 set to $(suggested_peak_selection_min_shape_r2)") + elseif param_key == :frequency_threshold + suggested_peak_selection_frequency_threshold = string(processed_value) + peak_selection_frequency_threshold = string(processed_value) + println(" suggested_peak_selection_frequency_threshold set to $(suggested_peak_selection_frequency_threshold)") + elseif param_key == :correlation_threshold + suggested_peak_selection_correlation_threshold = string(processed_value) + peak_selection_correlation_threshold = string(processed_value) + println(" suggested_peak_selection_correlation_threshold set to $(suggested_peak_selection_correlation_threshold)") + end + elseif step_name == :PeakBinning + if param_key == :tolerance + suggested_binning_tolerance = string(processed_value) + binning_tolerance = string(processed_value) + println(" suggested_binning_tolerance set to $(suggested_binning_tolerance)") + elseif param_key == :frequency_threshold + suggested_binning_frequency_threshold = string(processed_value) + binning_frequency_threshold = string(processed_value) + println(" suggested_binning_frequency_threshold set to $(suggested_binning_frequency_threshold)") + elseif param_key == :min_peak_per_bin + suggested_binning_min_peak_per_bin = string(processed_value) + binning_min_peak_per_bin = string(processed_value) + println(" suggested_binning_min_peak_per_bin set to $(suggested_binning_min_peak_per_bin)") + elseif param_key == :max_bin_width_ppm + suggested_binning_max_bin_width_ppm = string(processed_value) + binning_max_bin_width_ppm = string(processed_value) + println(" suggested_binning_max_bin_width_ppm set to $(suggested_binning_max_bin_width_ppm)") + elseif param_key == :num_uniform_bins + suggested_binning_num_uniform_bins = string(processed_value) + binning_num_uniform_bins = string(processed_value) + println(" suggested_binning_num_uniform_bins set to $(suggested_binning_num_uniform_bins)") + end + end + end + end + end + + # Also set method types for steps + if haskey(recommended_params, :Smoothing) && haskey(recommended_params[:Smoothing], :method) + smoothing_method = string(recommended_params[:Smoothing][:method]) + end + if haskey(recommended_params, :BaselineCorrection) && haskey(recommended_params[:BaselineCorrection], :method) + baseline_method = string(recommended_params[:BaselineCorrection][:method]) + end + if haskey(recommended_params, :Normalization) && haskey(recommended_params[:Normalization], :method) + normalization_method = string(recommended_params[:Normalization][:method]) + end + if haskey(recommended_params, :PeakAlignment) && haskey(recommended_params[:PeakAlignment], :method) + alignment_method = string(recommended_params[:PeakAlignment][:method]) + end + if haskey(recommended_params, :PeakPicking) && haskey(recommended_params[:PeakPicking], :method) + peak_picking_method = string(recommended_params[:PeakPicking][:method]) + end + if haskey(recommended_params, :PeakBinning) && haskey(recommended_params[:PeakBinning], :method) + binning_method = string(recommended_params[:PeakBinning][:method]) + end + if haskey(recommended_params, :PeakAlignment) && haskey(recommended_params[:PeakAlignment], :tolerance_unit) + alignment_tolerance_unit = string(recommended_params[:PeakAlignment][:tolerance_unit]) + end + if haskey(recommended_params, :PeakBinning) && haskey(recommended_params[:PeakBinning], :tolerance_unit) + binning_tolerance_unit = string(recommended_params[:PeakBinning][:tolerance_unit]) + end + + msg = "File loaded and parameters suggested." + catch e + @warn "Could not suggest parameters. Using defaults. Error: $e" + end + metadata_columns = [ Dict("name" => "parameter", "label" => "Parameter", "field" => "parameter", "align" => "left"), Dict("name" => "value", "label" => "Value", "field" => "value", "align" => "left"), @@ -511,15 +838,271 @@ end btnMetadataDisable = true @error "File loading failed" exception=(e, catch_backtrace()) finally - GC.gc() + GC.gc() # Trigger garbage collection if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) + ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS end progress = false progressSpectraPlot = false end end + @onbutton export_params_btn begin + params_to_export = Dict( + "pipeline_step_order" => pipeline_step_order, + "enable_standards" => enable_standards, # Export global flag + "stabilization_method" => stabilization_method, + "smoothing_method" => smoothing_method, + "smoothing_window" => smoothing_window, + "smoothing_order" => smoothing_order, + "baseline_method" => baseline_method, + "baseline_iterations" => baseline_iterations, + "baseline_window" => baseline_window, + "normalization_method" => normalization_method, + "alignment_method" => alignment_method, + "alignment_span" => alignment_span, + "alignment_tolerance" => alignment_tolerance, + "alignment_tolerance_unit" => alignment_tolerance_unit, + "alignment_max_shift_ppm" => alignment_max_shift_ppm, + "alignment_min_matched_peaks" => alignment_min_matched_peaks, + "peak_picking_method" => peak_picking_method, + "peak_picking_snr_threshold" => peak_picking_snr_threshold, + "peak_picking_half_window" => peak_picking_half_window, + "peak_picking_min_peak_prominence" => peak_picking_min_peak_prominence, + "peak_picking_merge_peaks_tolerance" => peak_picking_merge_peaks_tolerance, + "peak_picking_min_peak_width_ppm" => peak_picking_min_peak_width_ppm, + "peak_picking_max_peak_width_ppm" => peak_picking_max_peak_width_ppm, + "peak_picking_min_peak_shape_r2" => peak_picking_min_peak_shape_r2, + "binning_method" => binning_method, + "binning_tolerance" => binning_tolerance, + "binning_tolerance_unit" => binning_tolerance_unit, + "binning_frequency_threshold" => binning_frequency_threshold, + "binning_min_peak_per_bin" => binning_min_peak_per_bin, + "binning_max_bin_width_ppm" => binning_max_bin_width_ppm, + "binning_intensity_weighted_centers" => binning_intensity_weighted_centers, + "binning_num_uniform_bins" => binning_num_uniform_bins, + "calibration_fit_order" => calibration_fit_order, + "calibration_ppm_tolerance" => calibration_ppm_tolerance, + "peak_selection_min_snr" => peak_selection_min_snr, + "peak_selection_min_fwhm_ppm" => peak_selection_min_fwhm_ppm, + "peak_selection_max_fwhm_ppm" => peak_selection_max_fwhm_ppm, + "peak_selection_min_shape_r2" => peak_selection_min_shape_r2, + "peak_selection_frequency_threshold" => peak_selection_frequency_threshold, + "peak_selection_correlation_threshold" => peak_selection_correlation_threshold, + "reference_peaks_list" => reference_peaks_list + ) + json_string = JSON.json(params_to_export) + js_script = """ + var element = document.createElement('a'); + element.setAttribute('href', 'data:text/json;charset=utf-8,' + encodeURIComponent(`$json_string`)); + element.setAttribute('download', 'preprocessing_params.json'); + element.style.display = 'none'; + document.body.appendChild(element); + element.click(); + document.body.removeChild(element); + """ + run_js(js_script) + msg = "Parameters exported." + end + + @onchange import_params_btn begin + if import_params_btn + try + json_string = String(imported_params_file.data) + params = JSON.parse(json_string) + + for (key, value) in params + if key == "reference_peaks_list" + reference_peaks_list = value + elseif key == "pipeline_step_order" + pipeline_step_order = value + elseif key == "enable_standards" # Import global flag + enable_standards = value + else + # Use getfield and setproperty! to update reactive variables by name + if hasfield(typeof(@__MODULE__), Symbol(key)) + getfield(@__MODULE__, Symbol(key))[] = value + end + end + end + msg = "Parameters imported successfully." + catch e + msg = "Failed to import parameters: $e" + warning_msg = true + finally + import_params_btn = false # Reset button state + end + end + end + + @onbutton run_full_pipeline begin + progressPrep = true + current_pipeline_step = "Initializing..." + + try + # 1. Data Preparation + if isempty(selected_folder_main) + msg = "No dataset selected. Please process a file first." + warning_msg = true + return + end + + registry = load_registry(registry_path) + entry = registry[selected_folder_main] + target_path = entry["source_path"] + + if msi_data === nothing || full_route != target_path + if msi_data !== nothing + close(msi_data) + end + full_route = target_path + msi_data = OpenMSIData(target_path) + end + + # Apply mask if enabled + spectrum_indices_to_process = collect(1:length(msi_data.spectra_metadata)) + if maskEnabled && !isempty(preprocessing_mask_route) && isfile(preprocessing_mask_route) + current_pipeline_step = "Loading mask..." + mask_matrix = load_and_prepare_mask(preprocessing_mask_route, msi_data.image_dims) + masked_indices_set = get_masked_spectrum_indices(msi_data, mask_matrix) + spectrum_indices_to_process = collect(masked_indices_set) + end + + # Initialize spectra data structure + current_pipeline_step = "Loading spectra..." + current_spectra = Vector{MutableSpectrum}(undef, length(spectrum_indices_to_process)) + + Threads. @threads for i in 1:length(spectrum_indices_to_process) + original_idx = spectrum_indices_to_process[i] + mz, intensity = GetSpectrum(msi_data, original_idx) + current_spectra[i] = MutableSpectrum(original_idx, mz, intensity, []) + end + + # 2. Parameter Assembly + current_pipeline_step = "Configuring parameters..." + final_params = Dict( + :Stabilization => Dict(:method => Symbol(stabilization_method)), + :Smoothing => Dict( + :method => Symbol(smoothing_method), + :window => parse(Int, smoothing_window), + :order => parse(Int, smoothing_order) + ), + :BaselineCorrection => Dict( + :method => Symbol(baseline_method), + :iterations => parse(Int, baseline_iterations), + :window => parse(Int, baseline_window) + ), + :Normalization => Dict(:method => Symbol(normalization_method)), + :PeakPicking => Dict( + :method => Symbol(peak_picking_method), + :snr_threshold => parse(Float64, peak_picking_snr_threshold), + :half_window => parse(Int, peak_picking_half_window), + :min_peak_prominence => parse(Float64, peak_picking_min_peak_prominence), + :merge_peaks_tolerance => parse(Float64, peak_picking_merge_peaks_tolerance) + ), + :PeakSelection => Dict( + :min_snr => parse(Float64, peak_selection_min_snr), + :min_fwhm_ppm => parse(Float64, peak_selection_min_fwhm_ppm), + :max_fwhm_ppm => parse(Float64, peak_selection_max_fwhm_ppm), + :min_shape_r2 => parse(Float64, peak_selection_min_shape_r2) + ), + :Calibration => Dict( + :method => :internal_standards, + :ppm_tolerance => parse(Float64, calibration_ppm_tolerance), + :fit_order => parse(Int, calibration_fit_order) + ), + :PeakAlignment => Dict( + :method => Symbol(alignment_method), + :tolerance => parse(Float64, alignment_tolerance), + :tolerance_unit => Symbol(alignment_tolerance_unit) + ), + :PeakBinning => Dict( + :method => Symbol(binning_method), + :tolerance => parse(Float64, binning_tolerance), + :tolerance_unit => Symbol(binning_tolerance_unit), + :min_peak_per_bin => parse(Int, binning_min_peak_per_bin) + ) + ) + + # Build pipeline steps from enabled steps in order + pipeline_stp = [step["name"] for step in pipeline_step_order if step["enabled"]] + + # Convert reference peaks + ref_peaks = Dict(p["mz"] => p["label"] for p in reference_peaks_list) + + # 3. Execute Pipeline + current_pipeline_step = "Running preprocessing pipeline..." + feature_matrix_result, bin_info_result = execute_full_preprocessing( + current_spectra, + final_params, + pipeline_stp, + ref_peaks, + maskEnabled ? preprocessing_mask_route : nothing + ) do step + current_pipeline_step = "Processing: $step" + end + + # 4. Update Results Display + current_pipeline_step = "Updating results..." + + # Find the spectrum to display in "after" plot + display_spectrum_idx = findfirst(s -> s.id == selected_spectrum_id_for_plot, current_spectra) + if display_spectrum_idx !== nothing + processed_spectrum = current_spectra[display_spectrum_idx] + + # Create "after" plot data + after_trace = PlotlyBase.scatter( + x=processed_spectrum.mz, + y=processed_spectrum.intensity, + mode="lines", + name="Processed Spectrum" + ) + + traces_after = [after_trace] + + # Add peaks if they exist + if !isempty(processed_spectrum.peaks) + peak_mzs = [p.mz for p in processed_spectrum.peaks] + peak_intensities = [p.intensity for p in processed_spectrum.peaks] # Fixed: Should be peak.intensity + peak_trace = PlotlyBase.scatter( + x=peak_mzs, + y=peak_intensities, + mode="markers", + name="Picked Peaks", + marker=attr(color="red", size=8) + ) + push!(traces_after, peak_trace) + end + + plotdata_after = traces_after + plotlayout_after = PlotlyBase.Layout( + title="After Preprocessing (Spectrum $selected_spectrum_id_for_plot)", + xaxis_title="m/z", + yaxis_title="Intensity" + ) + end + + # Save feature matrix if binning was performed + if feature_matrix_result !== nothing + output_dir = joinpath("public", selected_folder_main, "preprocessing_results") + mkpath(output_dir) + save_feature_matrix(feature_matrix_result, bin_info_result, output_dir) + msg = "Pipeline completed successfully. Feature matrix saved." + else + msg = "Pipeline completed successfully." + end + + catch e + msg = "Error during pipeline execution: $e" + warning_msg = true + @error "Pipeline failed" exception=(e, catch_backtrace()) + finally + progressPrep = false + current_pipeline_step = "" + GC.gc() + end + end + # This new handler correctly adds the file from full_route to the batch list. @onbutton btnAddBatch begin if isempty(full_route) || full_route == "unknown (manually added)" @@ -704,7 +1287,7 @@ end files_without_mask = 0 for (file_idx, file_path) in enumerate(current_selected_files) - progress_message = "Processing file $file_idx/$num_files: $(basename(file_path))" + progress_message = "Processing file $(file_idx)/$(num_files): $(basename(file_path))" overall_progress = current_step / total_steps all_params = ( @@ -759,7 +1342,7 @@ end Errors by category: • Load failures: $(length(errors["load_errors"])) - • Slice generation: $(length(errors["slice_errors"])) + • Slice generation: $(length(errors["slice_errors"])) • I/O issues: $(length(errors["io_errors"])) Detailed errors: @@ -832,9 +1415,9 @@ end SpectraEnabled = true overall_progress = 0.0 println("Done") - GC.gc() + GC.gc() # Trigger garbage collection if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) + ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS end end end @@ -891,6 +1474,8 @@ end end plotdata, plotlayout, xSpectraMz, ySpectraMz = meanSpectrumPlot(msi_data, selected_folder_main, mask_path=mask_path_for_plot) + plotdata_before = plotdata + plotlayout_before = plotlayout selectedTab = "tab2" fTime = time() eTime = round(fTime - sTime, digits=3) @@ -905,9 +1490,9 @@ end btnPlotDisable = false btnSpectraDisable = false btnStartDisable = false - GC.gc() + GC.gc() # Trigger garbage collection if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) + ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS end end end @@ -922,37 +1507,37 @@ end progressSpectraPlot = true btnPlotDisable = true btnStartDisable = true - msg = "Loading total spectrum plot for $(selected_folder_main)..." + msg = "Loading total spectrum plot for $(selected_folder_main)..." - try - sTime = time() - registry = load_registry(registry_path) - entry = registry[selected_folder_main] - target_path = entry["source_path"] - - if target_path == "unknown (manually added)" - msg = "Dataset selected contained no route." - warning_msg = true - return - end - - if msi_data === nothing || full_route != target_path - if msi_data !== nothing - close(msi_data) - end - msg = "Reloading $(basename(target_path)) for analysis..." - full_route = target_path - msi_data = OpenMSIData(target_path) - if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing - raw_min = entry["metadata"]["global_min_mz"] - raw_max = entry["metadata"]["global_max_mz"] - min_val = isa(raw_min, Dict) ? get(raw_min, "value", raw_min) : raw_min - max_val = isa(raw_max, Dict) ? get(raw_max, "value", raw_max) : raw_max - set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val)) - else - precompute_analytics(msi_data) - end - end + try + sTime = time() + registry = load_registry(registry_path) + entry = registry[selected_folder_main] + target_path = entry["source_path"] + + if target_path == "unknown (manually added)" + msg = "Dataset selected contained no route." + warning_msg = true + return + end + + if msi_data === nothing || full_route != target_path + if msi_data !== nothing + close(msi_data) + end + msg = "Reloading $(basename(target_path)) for analysis..." + full_route = target_path + msi_data = OpenMSIData(target_path) + if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing + raw_min = entry["metadata"]["global_min_mz"] + raw_max = entry["metadata"]["global_max_mz"] + min_val = isa(raw_min, Dict) ? get(raw_min, "value", raw_min) : raw_min + max_val = isa(raw_max, Dict) ? get(raw_max, "value", raw_max) : raw_max + set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val)) + else + precompute_analytics(msi_data) + end + end local mask_path_for_plot::Union{String, Nothing} = nothing if maskEnabled && get(entry, "has_mask", false) mask_path_for_plot = get(entry, "mask_path", "") @@ -963,6 +1548,8 @@ end end plotdata, plotlayout, xSpectraMz, ySpectraMz = sumSpectrumPlot(msi_data, selected_folder_main, mask_path=mask_path_for_plot) + plotdata_before = plotdata + plotlayout_before = plotlayout selectedTab = "tab2" fTime = time() eTime = round(fTime - sTime, digits=3) @@ -977,9 +1564,9 @@ end btnPlotDisable = false btnSpectraDisable = false btnStartDisable = false - GC.gc() + GC.gc() # Trigger garbage collection if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) + ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS end end end @@ -1003,7 +1590,7 @@ end # Add error handling for registry access if !haskey(registry, selected_folder_main) - msg = "Dataset '$selected_folder_main' not found in registry." + msg = "Dataset '$(selected_folder_main)' not found in registry." warning_msg = true return end @@ -1047,6 +1634,8 @@ end # Convert to positive coordinates for processing y_positive = yCoord < 0 ? abs(yCoord) : yCoord plotdata, plotlayout, xSpectraMz, ySpectraMz = xySpectrumPlot(msi_data, xCoord, y_positive, imgWidth, imgHeight, selected_folder_main, mask_path=mask_path_for_plot) + plotdata_before = plotdata + plotlayout_before = plotlayout # Update coordinates based on actual plot title # Extract title text from the Dict safely @@ -1092,9 +1681,9 @@ end btnPlotDisable = false btnSpectraDisable = false btnStartDisable = false - GC.gc() + GC.gc() # Trigger garbage collection if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) + ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS end end end @@ -1375,9 +1964,9 @@ end end msi_data = nothing log_memory_usage("Folder Changed (msi_data cleared)", msi_data) - GC.gc() + GC.gc() # Trigger garbage collection if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) + ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS end if !isempty(selected_folder_main) @@ -1623,9 +2212,9 @@ end btnStartDisable=false btnSpectraDisable=false SpectraEnabled=true - GC.gc() + GC.gc() # Trigger garbage collection if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) + ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS end end end @@ -1650,7 +2239,7 @@ end try # --- Get Mask Path --- local mask_path_for_plot::Union{String, Nothing} = nothing - if maskEnabled[] && !isempty(selected_folder_main) + if maskEnabled && !isempty(selected_folder_main) registry = load_registry(registry_path) entry = get(registry, selected_folder_main, nothing) if entry !== nothing && get(entry, "has_mask", false) @@ -1686,9 +2275,9 @@ end btnStartDisable=false btnSpectraDisable=false SpectraEnabled=true - GC.gc() + GC.gc() # Trigger garbage collection if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) + ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS end end end @@ -1789,16 +2378,35 @@ end # To include a visualization in the spectrum plot indicating where is the selected mass @onchange Nmass begin if !isempty(xSpectraMz) - # Main spectrum trace - traceSpectra = PlotlyBase.stem( - x=xSpectraMz, - y=ySpectraMz, - marker=attr(size=1, color="blue", opacity=0.5), - name="Spectrum", - hoverinfo="x", - hovertemplate="m/z: %{x:.4f}", - showlegend=false - ) + df = msi_data.spectrum_stats_df + profile_count = 0 + if df !== nothing && hasproperty(df, :Mode) + profile_count = count(==(MSI_src.PROFILE), df.Mode) + end + + if profile_count > 0 + # Main spectrum trace + traceSpectra = PlotlyBase.scatter( + x=xSpectraMz, + y=ySpectraMz, + marker=attr(size=1, color="blue", opacity=0.5), + name="Spectrum", + hoverinfo="x", + hovertemplate="m/z: %{x:.4f}", + showlegend=false + ) + else + # Main spectrum trace + traceSpectra = PlotlyBase.stem( + x=xSpectraMz, + y=ySpectraMz, + marker=attr(size=1, color="blue", opacity=0.5), + name="Spectrum", + hoverinfo="x", + hovertemplate="m/z: %{x:.4f}", + showlegend=false + ) + end # Parse all valid masses from the comma-separated string mass_strs = split(Nmass, ',', keepempty=false) @@ -1955,9 +2563,11 @@ end @onchange isready @time begin if isready && !registry_init_done + initialization_message = "Pre-compiling functions at startup..." warmup_init() + initialization_message = "Pre-compilation finished." try - println("Synchronizing registry with filesystem on backend init...") + initialization_message = "Synchronizing registry with filesystem on backend init..." reg_path = abspath(joinpath(@__DIR__, "public", "registry.json")) registry = isfile(reg_path) ? load_registry(reg_path) : Dict{String, Any}() @@ -1986,10 +2596,8 @@ end end if !isempty(new_folders) || !isempty(removed_folders) - println("Registry changed, saving...") - open(reg_path, "w") do f - JSON.print(f, registry, 4) - end + initialization_message = "Registry changed, saving..." + save_registry(reg_path, registry) end all_folders = sort(collect(keys(registry)), lt=natural) @@ -2006,8 +2614,11 @@ end selected_files = String[] finally registry_init_done = true + is_initializing = false # Hide loading screen when initialization is complete end end + is_initializing = false # Hide loading screen when initialization is complete (current code is hidden due to incompatibility) + initialization_message = "Done." log_memory_usage("App Ready", msi_data) end @@ -2019,4 +2630,4 @@ end # == Pages == # Register a new route and the page that will be loaded on access @page("/", "app.jl.html") -end +end \ No newline at end of file diff --git a/app.jl.html b/app.jl.html index 47a3673..ae5a7f3 100644 --- a/app.jl.html +++ b/app.jl.html @@ -4,6 +4,16 @@

JuliaMSI 

+ + +
@@ -16,219 +26,408 @@ -
imzML & mzML Data Pre-Treatment
-
- - - - - +
imzML & mzML Data Pre-Treatment
+ + +
+ + + + + +
+ + +
+ +
+ + +
+
Preview Spectrum:
+ + + + Mean spectrum plot + + + Sum Spectrum plot + + + Spectrum plot (X,Y) + + + +
+
+ +
+
+ +
+
+ +
+
+
+ + + + +
Internal Standards
+
Reference peaks for calibration and alignment
+
+ + + + + + + + +
+
+ +
+
+ +
+
+
+
+ + + + + +
+
+
+ + +
Preprocessing Pipeline
+ + + + + + + + + + + + + + + + + + + + + +
+
+ +
+
+ +
+
+
+
+ + + +
Method
+
The algorithm to use for baseline correction.
+
+ +
+
+
+
+ +
Parameters
+
+ +
+
+ +
+
+ +
+
+
+
+ + + +
Method
+
The normalization method to apply.
+
+ +
+
+
+
+
+
+ + + +
Method
+
The alignment algorithm.
+
+ +
+
+
+
+ +
Parameters
+
+ +
+
+ +
+
+ +
+
+ +
+
+ +
+
+ +
+
+
+
+ + + +
Parameters
+
+ +
+
+ +
+
+ +
+
+
+
+ + + +
Method
+
The peak detection algorithm.
+
+ +
+
+
+
+ +
Parameters
+
+ +
+
+ +
+
+ +
+
+
+
+ +
+
+ +
+
+
+
+ +
+
+ +
+
+ +
+
+ + + +
Peak Quality Filters
+
+ +
+
+ +
+
+ +
+
+
+
+ +
+
+ +
+
+
+
+ +
+
+ +
+
+
+
+ + + +
Method
+
The binning strategy.
+
+ +
+
+
+ +
Parameters
+
+ +
+
+ +
+
+ +
+
+
+
+ +
+
+ +
+
+
+
+ +
+
+ +
+
+ +
+
+ +
+
+ + +
+ + + +
+ + + + + + +
+ + + +
+
imzML & mzML Data Processor

Please make sure the ibd and imzML file are located in the same directory and have the same name. @@ -330,6 +529,10 @@

+
+ +
Before Preprocessing
- + @@ -790,4 +993,4 @@ - \ No newline at end of file + diff --git a/app_snippet.html b/app_snippet.html new file mode 100644 index 0000000..398a10d --- /dev/null +++ b/app_snippet.html @@ -0,0 +1,665 @@ + + + + +
+
+ +
+ + + + + + + + +
imzML & mzML Data Pre-Treatment
+
+ + + + + +
+ + + + {{ file }} + + + + + + + +
+
+ + + + + + + + + + + + +
+ + + +
+ + + + + + +
Method
+
Applies a variance-stabilizing transformation to the intensity vector.
+
+ +
+
+
+
+
+
+
+
+ + + + +
Method
+
The smoothing algorithm.
+
+ +
+
+
+
+ + +
Parameters
+
+ + + + +
+
+ + + + +
Method
+
The algorithm to use for baseline correction.
+
+ +
+
+
+
+
+ + +
Parameters
+
+ + + + +
+
+ + + + +
Method
+
The normalization method to apply.
+
+ +
+
+
+
+
+
+
+ + + + +
Method
+
The alignment algorithm.
+
+ +
+
+
+
+
+ + +
Parameters
+
+ + + + + + + +
+
+ + +
Define Internal Standards / Reference Peaks
+

Provide m/z values and optional labels for internal standards or reference peaks. These are used for mass calibration and alignment.

+ + + + + + +
+
+ +
+
+ +
+
+
+
+ + + + + +
+
+ + +

This step uses the peaks defined in the 'Internal Standards' tab to correct the m/z axis.

+ + +
Parameters
+
+ + + + +
+
+ + +

Select the appropriate peak picking method and set its parameters.

+ + +
Method
+
The peak detection algorithm.
+
+ +
+
+
+
+
+ + +
Parameters
+
+ + + + + + + + + +
+
+ + + + +
Peak Quality Filters
+
+ + + + + + + + +
+
+ + + + +
Method
+
The binning strategy.
+
+ +
+
+
+
+ + +
Parameters
+
+ + + + + + + + + +
+
+
+
+ +
imzML & mzML Data Processor
+

Please make sure the ibd and imzML file are located in the same directory and have the same name. +
It may take a while to generate the slice / spectrum, please be patient. +
To generate the contour or surface plots, you have to select the desired slice first using the + interface. +

+
+ + + + + +
+ + + + {{ file }} + + + + + + + + +
+
+ + +
+
+ +
+
+ +
+
+
+ +
+
+ + + +
+
+
+ +
+
+
+ +
+
+ + + + + + + Mean spectrum plot + + + + + Sum Spectrum plot + + + + + Spectrum plot (X,Y) + + + + +
+
+
+ +
+
+ +
+
+
+
+
+ + + + + + +
+ +
{{ progress_message }}
+
+
+

{{msg}}

+
+ + + + + + + Image topography Plot + + + + + + TrIQ topography Plot + + + + + + Image surface Plot + + + + + + TrIQ Surface Plot + + + + + + + + Over normal image + + + + + Over TrIQ image + + + +
+ + Transparency: {{ imgTrans }} +
+
+
+ +
mzML to imzML Converter
+

Select the .mzML file and the corresponding .txt synchronization file to convert them into an .imzML/.ibd + pair.

+ + + + + + + + + + + +

{{msg_conversion}}

+
+
+
+
+
+ +
+
+
Spectrum View
+ + +
Before Preprocessing
+
+ + + +
+ + +
After Preprocessing
+
+ + + + +
+
+
+ + + + +
Image visualizer
+
+ + + + +
+ +
+
+ +
+
+ +
+
+

+
+ + + +
TrIQ visualizer
+
+ + + + +
+ +
+
+ +
+
+ +
+
+

+
+ + +
+ + + + + + + + Mean spectrum plot + + + + + Sum Spectrum plot + + + + + Spectrum plot (X,Y) + + + + +
+ +
+ + + + + + + + + +
+
+
+
+
diff --git a/app_snippet.jl b/app_snippet.jl new file mode 100644 index 0000000..118af40 --- /dev/null +++ b/app_snippet.jl @@ -0,0 +1,1524 @@ +# app.jl + +module App +# ==Packages == +using GenieFramework +using Pkg +using Libz +using PlotlyBase +using CairoMakie +using Colors +using MSI_src # Import the new MSIData library +using Statistics +using NaturalSort +using Images +using LinearAlgebra +using NativeFileDialog # Opens the file explorer depending on the OS +using StipplePlotly +using Base.Filesystem: mv # To rename files in the system +using Printf # Required for @sprintf macro in colorbar generation +using JSON +using Dates +using Base.Threads + +# Bring MSIData into App module's scope +using .MSI_src: MSIData, OpenMSIData, process_spectrum, IterateSpectra, ImzMLSource, _iterate_spectra_fast, MzMLSource, find_mass, ViridisPalette, get_mz_slice, get_multiple_mz_slices, quantize_intensity, save_bitmap, median_filter, save_bitmap, downsample_spectrum, TrIQ, precompute_analytics, ImportMzmlFile, generate_colorbar_image, load_and_prepare_mask, set_global_mz_range!, main_precalculation, MutableSpectrum + +if !@isdefined(increment_image) + include("./julia_imzML_visual.jl") +end + +# --- Memory Validation Logging --- +if get(ENV, "GENIE_ENV", "dev") != "prod" + function get_rss_mb() + if !Sys.islinux() + return 0.0 + end + try + pid = getpid() + cmd = `ps -p $pid -o rss=` + rss_kb_str = read(cmd, String) + rss_kb = parse(Int, strip(rss_kb_str)) + return round(rss_kb / 1024, digits=2) + catch e + @warn "Could not get RSS via `ps` command. Error: $e" + return 0.0 + end + end + + function log_memory_usage(context::String, msi_data_val) + rss_mb = get_rss_mb() + + msi_data_size_mb = 0 + if msi_data_val !== nothing + msi_data_size_mb = round(Base.summarysize(msi_data_val) / (1024^2), digits=2) + end + + gc_time_s = round(GC.time(), digits=3) + + println("--- MEMORY LOG [$(context)] ---") + println(" Timestamp: $(now())") + println(" Process RSS: $(rss_mb) MB") + println(" msi_data size: $(msi_data_size_mb) MB") + println(" Cumulative GC time: $(gc_time_s) s") + println("--------------------------") + end +else + log_memory_usage(context::String, msi_data_val) = nothing # No-op for production +end + +@genietools + +# == Reactive code == +# Reactive code to make the UI interactive +@app begin + # == Loading Screen Variables == + @in is_initializing = true + @in initialization_message = "Initializing..." + + # == Reactive variables == + # reactive variables exist in both the Julia backend and the browser with two-way synchronization + # @out variables can only be modified by the backend + # @in variables can be modified by both the backend and the browser + # variables must be initialized with constant values, or variables defined outside of the @app block + + ## Interface non Variables + @out btnStartDisable=true + @out btnPlotDisable=false + @out btnSpectraDisable=false + # Loading animations + @in progress=false + @in progressPlot=false + @in progressSpectraPlot=false + # Text field validations + @in triqEnabled=false + @in SpectraEnabled=false + @in MFilterEnabled=false + @in maskEnabled=false + # Dialogs + @in warning_msg=false + @in CompareDialog=false + + ## Interface Variables + @in file_route="" + @in file_name="" + @in Nmass="0.0" + @in Tol=0.1 + @in triqProb=0.98 + @in colorLevel=20 + + ## Interface Buttons + @in btnSearch=false # To search for files in your device + @in btnAddBatch = false + @in clear_batch_btn = false + @out batch_file_count = 0 + @in mainProcess=false # To generate images + @in compareBtn=false # To open dialog + @in createMeanPlot=false # To generate mean spectrum plot + @in createXYPlot=false # To generate an spectrum plot according to the xy values inputed + @in createSumPlot=false # To generate a sum of all the spectrum plots + @in image3dPlot=false # To generate 3d plot based on current image + @in triq3dPlot=false # To generate 3d plot based on current triq image + @in imageCPlot=false # To generate contour plots of current image + @in triqCPlot=false # To generate contour plots of current triq image + # Image change buttons + @in imgPlus=false + @in imgMinus=false + @in imgPlusT=false + @in imgMinusT=false + # Image change comparative buttons + @in imgPlusCompLeft=false + @in imgMinusCompLeft=false + @in imgPlusTCompLeft=false + @in imgMinusTCompLeft=false + @in imgPlusCompRight=false + @in imgMinusCompRight=false + @in imgPlusTCompRight=false + @in imgMinusTCompRight=false + + ## Tabulation variables + @out tabIDs=["tab0","tab1","tab2","tab3","tab4"] + @out tabLabels=["Image", "TrIQ", "Spectrum Plot", "Topography Plot","Surface Plot"] + @in selectedTab="tab0" + @out CompTabIDsLeft=["tab0","tab1","tab2","tab3","tab4"] + @out CompTabLabelsLeft=["Image", "TrIQ", "Spectrum Plot", "Topography Plot","Surface Plot"] + @in CompSelectedTabLeft="tab0" + @out CompTabIDsRight=["tab0","tab1","tab2","tab3","tab4"] + @out CompTabLabelsRight=["Image", "TrIQ", "Spectrum Plot", "Topography Plot","Surface Plot"] + @in CompSelectedTabRight="tab0" + + # Interface Images + @out imgInt="/.bmp" # image Interface + @out imgIntT="/.bmp" # image Interface TrIQ + @out colorbar="/.png" + @out colorbarT="/.png" + # Interface controlling for the comparative view + @out imgIntCompLeft="/.bmp" + @out imgIntTCompLeft="/.bmp" + @out colorbarCompLeft="/.png" + @out colorbarTCompLeft="/.png" + @out imgIntCompRight="/.bmp" + @out imgIntTCompRight="/.bmp" + @out colorbarCompRight="/.png" + @out colorbarTCompRight="/.png" + @out imgWidth=0 + @out imgHeight=0 + + # Optical Image Overlay & Transparency + @in imgTrans=1.0 + @in progressOptical=false + @out btnOpticalDisable=true + @in btnOptical=false + @in btnOpticalT=false + @in opticalOverTriq=false + @out imgRoute="" + + # Messages to interface variables + @out msg="" + @out msgimg="" + @out msgtriq="" + # Reiteration of the messages under the image to know which spectra is being visualized + @out msgimgCompLeft="" + @out msgtriqCompLeft="" + @out msgimgCompRight="" + @out msgtriqCompRight="" + + # Centralized MSIData object + @out msi_data::Union{MSIData, Nothing} = nothing + + # Metadata table variables + @in showMetadataDialog = false + @in showMetadataBtn = false + @out metadata_columns = [] + @out metadata_rows = [] + @out btnMetadataDisable = false + @in selected_folder_metadata = "" + + # Saves the route where imzML and mzML files are located + @out full_route="" + + # == Converter Tab Variables == + @in left_tab = "generator" + @out mzml_full_route = "" + @out sync_full_route = "" + @in btnSearchMzml = false + @in btnSearchSync = false + @in convert_process = false + @out progress_conversion = false + @out msg_conversion = "" + @out btnConvertDisable = true + + # == Pre Processing Variables == + @in pre_tab = "stabilization" + + ## Preprocessing Parameters + @in progressPrep=false + @in stabilization_method="sqrt" + @in smoothing_method="sg" + @in smoothing_window = "" + @in smoothing_order = "" + @in baseline_method="snip" + @in baseline_iterations = "" + @in baseline_window = "" + @in normalization_method="tic" + @in alignment_method="lowess" + @in alignment_span = "" + @in alignment_tolerance = "" + @in alignment_tolerance_unit="mz" + @in alignment_max_shift_ppm = "" + @in alignment_min_matched_peaks = "" + @in peak_picking_method="profile" + @in peak_picking_snr_threshold = "" + @in peak_picking_half_window = "" + @in peak_picking_min_peak_prominence = "" + @in peak_picking_merge_peaks_tolerance = "" + @in peak_picking_min_peak_width_ppm = "" + @in peak_picking_max_peak_width_ppm = "" + @in peak_picking_min_peak_shape_r2 = "" + @in binning_method="adaptive" + @in binning_tolerance = "" + @in binning_tolerance_unit="ppm" + @in binning_frequency_threshold = "" + @in binning_min_peak_per_bin = "" + @in binning_max_bin_width_ppm = "" + @in binning_intensity_weighted_centers=true + @in binning_num_uniform_bins = "" + @in calibration_fit_order = "" + @in calibration_ppm_tolerance = "" + @in peak_selection_min_snr = "" + @in peak_selection_min_fwhm_ppm = "" + @in peak_selection_max_fwhm_ppm = "" + @in peak_selection_min_shape_r2 = "" + @in peak_selection_frequency_threshold = "" + @in peak_selection_correlation_threshold = "" + + @in reference_peaks_list = [ + Dict("mz" => 137.0244, "label" => "DHB_fragment"), + Dict("mz" => 155.0349, "label" => "DHB_M+H"), + ] + + # --- Methods for Reference Peaks List --- + function addReferencePeak() + push!(reference_peaks_list, Dict("mz" => 0.0, "label" => "")) + reference_peaks_list = deepcopy(reference_peaks_list) # Force reactivity + end + + function removeReferencePeak(index::Int) + deleteat!(reference_peaks_list, index) + reference_peaks_list = deepcopy(reference_peaks_list) # Force reactivity + end + + # Individual step enable/disable flags + @in enable_stabilization = true + @in enable_smoothing = true + @in enable_baseline = true + @in enable_normalization = true + @in enable_standards = true + @in enable_alignment = true + @in enable_peak_picking = true + @in enable_binning = true + @in enable_calibration = true + @in enable_peak_selection = true + + # Trigger for running the full pipeline + @in run_full_pipeline = false + @out current_pipeline_step = "" # To indicate which step is currently running in the full pipeline + + @in export_params_btn = false + @in import_params_btn = false + @in imported_params_file = nothing + + @in suggested_smoothing_window = "" + @in suggested_smoothing_order = "" + @in suggested_baseline_iterations = "" + @in suggested_baseline_window = "" + @in suggested_alignment_span = "" + @in suggested_alignment_tolerance = "" + @in suggested_alignment_max_shift_ppm = "" + @in suggested_alignment_min_matched_peaks = "" + @in suggested_peak_picking_snr_threshold = "" + @in suggested_peak_picking_half_window = "" + @in suggested_peak_picking_min_peak_prominence = "" + @in suggested_peak_picking_merge_peaks_tolerance = "" + @in suggested_peak_picking_min_peak_width_ppm = "" + @in suggested_peak_picking_max_peak_width_ppm = "" + @in suggested_peak_picking_min_peak_shape_r2 = "" + @in suggested_binning_tolerance = "" + @in suggested_binning_frequency_threshold = "" + @in suggested_binning_min_peak_per_bin = "" + @in suggested_binning_max_bin_width_ppm = "" + @in suggested_binning_num_uniform_bins = "" + @in suggested_calibration_fit_order = "" + @in suggested_calibration_ppm_tolerance = "" + @in suggested_peak_selection_min_snr = "" + @in suggested_peak_selection_min_fwhm_ppm = "" + @in suggested_peak_selection_max_fwhm_ppm = "" + @in suggested_peak_selection_min_shape_r2 = "" + @in suggested_peak_selection_frequency_threshold = "" + @in suggested_peak_selection_correlation_threshold = "" + + # == Batch Summary Dialog == + @in showBatchSummary = false + @out batch_summary = "" + + # == Batch Processing & Registry Variables == + @private registry_init_done = false + @in refetch_folders = false + @in selected_files = String[] + @in available_folders = String[] + @in image_available_folders = String[] + @out registry_path = abspath(joinpath(@__DIR__, "public", "registry.json")) + # Progress reporting + @out overall_progress = 0.0 + @out progress_message = "" + + # == Folder-based UI State == + @in selected_folder_main = "" + @in selected_folder_compare_left = "" + @in selected_folder_compare_right = "" + + + # For the creation of images with a more specific mass charge + @out text_nmass="" + + # For image search image lists we apply a filter that searches specific type of images into our public folder, then we sort it in a "numerical" order + @in msi_bmp=sort(filter(filename -> startswith(filename, "MSI_") && endswith(filename, ".bmp"), readdir("public")),lt=natural) + @in col_msi_png=sort(filter(filename -> startswith(filename, "colorbar_MSI_") && endswith(filename, ".png"), readdir("public")),lt=natural) + @in triq_bmp=sort(filter(filename -> startswith(filename, "TrIQ_") && endswith(filename, ".bmp"), readdir("public")),lt=natural) + @in col_triq_png=sort(filter(filename -> startswith(filename, "colorbar_TrIQ_") && endswith(filename, ".png"), readdir("public")),lt=natural) + + # Set current image for the list to display + @out current_msi="" + @out current_col_msi="" + @out current_triq="" + @out current_col_triq="" + # We reiterate the process to display in the comparative view + @out current_msiCompLeft="" + @out current_col_msiCompLeft="" + @out current_triqCompLeft="" + @out current_col_triqCompLeft="" + @out current_msiCompRight="" + @out current_col_msiCompRight="" + @out current_triqCompRight="" + @out current_col_triqCompRight="" + + ## Time measurement variables + @out sTime=time() + @out fTime=time() + @out eTime=time() + + ## Plots + # Local image to plot + layoutImg=PlotlyBase.Layout( + title=PlotlyBase.attr( + text="", + font=PlotlyBase.attr( + family="Roboto, Lato, sans-serif", + size=14, + color="black" + ) + ), + xaxis=PlotlyBase.attr( + visible=false, + scaleanchor="y", + range=[0, 0] + ), + yaxis=PlotlyBase.attr( + visible=false, + range=[0, 0] + ), + margin=attr(l=0,r=0,t=0,b=0,pad=0) + ) + traceImg=PlotlyBase.heatmap(x=Vector{Float64}(), y=Vector{Float64}()) + @out plotdataImg=[traceImg] + @out plotlayoutImg=layoutImg + # For the image in the comparative view + @out plotdataImgCompLeft=[traceImg] + @out plotlayoutImgCompLeft=layoutImg + @out plotdataImgCompRight=[traceImg] + @out plotlayoutImgCompRight=layoutImg + + # For triq image + @out plotdataImgT=[traceImg] + @out plotlayoutImgT=layoutImg + # For the triq image in the comparative view + @out plotdataImgTCompLeft=[traceImg] + @out plotlayoutImgTCompLeft=layoutImg + @out plotdataImgTCompRight=[traceImg] + @out plotlayoutImgTCompRight=layoutImg + # Interface Plot Spectrum + layoutSpectra=PlotlyBase.Layout( + title=PlotlyBase.attr( + text="Spectrum plot", + font=PlotlyBase.attr( + family="Roboto, Lato, sans-serif", + size=18, + color="black" + ) + ), + hovermode="closest", + xaxis=PlotlyBase.attr( + title="m/z", + showgrid=true + ), + yaxis=PlotlyBase.attr( + title="Intensity", + showgrid=true, + tickformat = ".3g" + ), + margin=attr(l=0,r=0,t=120,b=0,pad=0) + ) + # Dummy 2D scatter plot + traceSpectra=PlotlyBase.scatter(x=Vector{Float64}(), y=Vector{Float64}(), mode="lines", marker=attr(size=1, color="blue", opacity=0.1)) + # Create conection to frontend + @out plotdata=[traceSpectra] + @out plotlayout=layoutSpectra + @in xCoord=0 + @in yCoord=0 + @out xSpectraMz = Vector{Float64}() + @out ySpectraMz = Vector{Float64}() + + # UI plot data for Preprocessing + @out plotdata_before = [traceSpectra] + @out plotlayout_before = layoutSpectra + @out plotdata_after = [traceSpectra] + @out plotlayout_after = layoutSpectra + + # Interactive plot reactions + @in data_click=Dict{String,Any}() + #@in data_selected=Dict{String,Any}() # Selected is for areas, this can work for the masks + # + + # Interface Plot Surface + layoutContour=PlotlyBase.Layout( + title=PlotlyBase.attr( + text="2D Topographic map", + font=PlotlyBase.attr( + family="Roboto, Lato, sans-serif", + size=18, + color="black" + ) + ), + xaxis=PlotlyBase.attr( + visible=false, + scaleanchor="y" + ), + yaxis=PlotlyBase.attr( + visible=false + ), + margin=attr(l=0,r=0,t=100,b=0,pad=0) + ) + # Dummy 2D surface plot + traceContour=PlotlyBase.contour(x=Vector{Float64}(), y=Vector{Float64}(), mode="lines") + # Create conection to frontend + @out plotdataC=[traceContour] + @out plotlayoutC=layoutContour + + # Interface Plot 3d + # Define the layout for the 3D plot + layout3D=PlotlyBase.Layout( + title=PlotlyBase.attr( + text="3D Surface plot", + font=PlotlyBase.attr( + family="Roboto, Lato, sans-serif", + size=18, + color="black" + ) + ), + scene=attr( + xaxis_title="X", + yaxis_title="Y", + zaxis_title="Z", + xaxis_nticks=20, + yaxis_nticks=20, + zaxis_nticks=4, + camera=attr(eye=attr(x=0, y=-1, z=0.5)), + aspectratio=attr(x=1, y=1, z=0.2) + ), + margin=attr(l=0,r=0,t=120,b=0,pad=0) + ) + + # Dummy 3D surface plot + x=1:10 + y=1:10 + z=[sin(i * j / 10) for i in x, j in y] + trace3D=PlotlyBase.surface(x=Vector{Float64}(), y=Vector{Float64}(), z=Matrix{Float64}(undef, 0, 0), + contours_z=attr( + show=true, + usecolormap=true, + highlightcolor="limegreen", + project_z=true + ), colorscale="Viridis") + # Create conection to frontend + @out plotdata3d=[trace3D] + @out plotlayout3d=layout3D + + # == Reactive handlers == + # Reactive handlers watch a variable and execute a block of code when its value changes + # The onbutton handler will set the variable to false after the block is executed + + # This handler correctly uses pick_file and loads the selected file + # as the active dataset for the UI. + @onbutton btnSearch begin + picked_route = pick_file(; filterlist="imzML,imzml,mzML,mzml") + if isempty(picked_route) + return + end + + progress = true + msg = "Opening file: $(basename(picked_route))..." + + try + dataset_name = replace(basename(picked_route), r"(\.(imzML|imzml|mzML|mzml))$ "i => "") + registry = load_registry(registry_path) + existing_entry = get(registry, dataset_name, nothing) + + # --- Fast Load Path --- + is_same_file = (existing_entry !== nothing && existing_entry["source_path"] == picked_route) + if is_same_file && !isempty(get(existing_entry, "metadata", Dict())) + msg = "Fast loading pre-processed file: $(dataset_name)" + println(msg) + + full_route = existing_entry["source_path"] + metadata_rows = existing_entry["metadata"]["summary"] + + dims_str = first(filter(r -> r["parameter"] == "Image Dimensions", metadata_rows))["value"] + dims = parse.(Int, split(dims_str, " x ")) + imgWidth, imgHeight = dims[1], dims[2] + + msi_data = nothing # Ensure data is not held in memory + log_memory_usage("Fast Load (msi_data cleared)", msi_data) + btnMetadataDisable = false + btnStartDisable = false + btnPlotDisable = false + btnSpectraDisable = false + SpectraEnabled = true + selected_folder_main = dataset_name + + # Update folder lists in UI + all_folders = sort(collect(keys(registry)), lt=natural) + img_folders = filter(folder -> get(get(registry, folder, Dict()), "is_imzML", false), all_folders) + available_folders = deepcopy(all_folders) + image_available_folders = deepcopy(img_folders) + + msg = "Successfully loaded pre-processed dataset: $(dataset_name)" + progress = false + return + end + + # --- Full Load Path --- + msg = "Performing first-time analysis for: $(basename(picked_route))..." + local local_full_route + if endswith(picked_route, r"imzml"i) + local_full_route = replace(picked_route, r"\.imzml$"i => ".imzML") + if picked_route != local_full_route + mv(picked_route, local_full_route, force=true) + end + else + local_full_route = picked_route + end + full_route = local_full_route + + sTime = time() + loaded_data = OpenMSIData(local_full_route) + is_imzML = loaded_data.source isa ImzMLSource + + precompute_analytics(loaded_data) + + # Auto-suggest parameters + try + println("Calling main_precalculation to get recommended parameters...") + recommended_params = main_precalculation(loaded_data) + + for (step_name, params) in recommended_params + for (param_key, value) in params + # Convert value to appropriate type before assignment + processed_value = if value === nothing + nothing + elseif value isa Tuple + @warn "Skipping invalid parameter suggestion (tuple): $value for $param_key" + "" # Set to empty string for safety + elseif value isa Number + value + else + string(value) + end + + if processed_value !== nothing + if step_name == :Smoothing + if param_key == :window + suggested_smoothing_window = string(processed_value) + smoothing_window = string(processed_value) + println(" suggested_smoothing_window set to $(suggested_smoothing_window)") + elseif param_key == :order + suggested_smoothing_order = string(processed_value) + smoothing_order = string(processed_value) + println(" suggested_smoothing_order set to $(suggested_smoothing_order)") + end + elseif step_name == :BaselineCorrection + if param_key == :iterations + suggested_baseline_iterations = string(processed_value) + baseline_iterations = string(processed_value) + println(" suggested_baseline_iterations set to $(suggested_baseline_iterations)") + elseif param_key == :window + suggested_baseline_window = string(processed_value) + baseline_window = string(processed_value) + println(" suggested_baseline_window set to $(suggested_baseline_window)") + end + elseif step_name == :PeakAlignment + if param_key == :span + suggested_alignment_span = string(processed_value) + alignment_span = string(processed_value) + println(" suggested_alignment_span set to $(suggested_alignment_span)") + elseif param_key == :tolerance + suggested_alignment_tolerance = string(processed_value) + alignment_tolerance = string(processed_value) + println(" suggested_alignment_tolerance set to $(suggested_alignment_tolerance)") + elseif param_key == :max_shift_ppm + suggested_alignment_max_shift_ppm = string(processed_value) + alignment_max_shift_ppm = string(processed_value) + println(" suggested_alignment_max_shift_ppm set to $(suggested_alignment_max_shift_ppm)") + elseif param_key == :min_matched_peaks + suggested_alignment_min_matched_peaks = string(processed_value) + alignment_min_matched_peaks = string(processed_value) + println(" suggested_alignment_min_matched_peaks set to $(suggested_alignment_min_matched_peaks)") + end + elseif step_name == :Calibration + if param_key == :fit_order + suggested_calibration_fit_order = string(processed_value) + calibration_fit_order = string(processed_value) + println(" suggested_calibration_fit_order set to $(suggested_calibration_fit_order)") + elseif param_key == :ppm_tolerance + suggested_calibration_ppm_tolerance = string(processed_value) + calibration_ppm_tolerance = string(processed_value) + println(" suggested_calibration_ppm_tolerance set to $(suggested_calibration_ppm_tolerance)") + end + elseif step_name == :PeakPicking + if param_key == :snr_threshold + suggested_peak_picking_snr_threshold = string(processed_value) + peak_picking_snr_threshold = string(processed_value) + println(" suggested_peak_picking_snr_threshold set to $(suggested_peak_picking_snr_threshold)") + elseif param_key == :half_window + suggested_peak_picking_half_window = string(processed_value) + peak_picking_half_window = string(processed_value) + println(" suggested_peak_picking_half_window set to $(suggested_peak_picking_half_window)") + elseif param_key == :min_peak_prominence + suggested_peak_picking_min_peak_prominence = string(processed_value) + peak_picking_min_peak_prominence = string(processed_value) + println(" suggested_peak_picking_min_peak_prominence set to $(suggested_peak_picking_min_peak_prominence)") + elseif param_key == :merge_peaks_tolerance + suggested_peak_picking_merge_peaks_tolerance = string(processed_value) + peak_picking_merge_peaks_tolerance = string(processed_value) + println(" suggested_peak_picking_merge_peaks_tolerance set to $(suggested_peak_picking_merge_peaks_tolerance)") + elseif param_key == :min_peak_width_ppm + suggested_peak_picking_min_peak_width_ppm = string(processed_value) + peak_picking_min_peak_width_ppm = string(processed_value) + println(" suggested_peak_picking_min_peak_width_ppm set to $(suggested_peak_picking_min_peak_width_ppm)") + elseif param_key == :max_peak_width_ppm + suggested_peak_picking_max_peak_width_ppm = string(processed_value) + peak_picking_max_peak_width_ppm = string(processed_value) + println(" suggested_peak_picking_max_peak_width_ppm set to $(suggested_peak_picking_max_peak_width_ppm)") + elseif param_key == :min_peak_shape_r2 + suggested_peak_picking_min_peak_shape_r2 = string(processed_value) + peak_picking_min_peak_shape_r2 = string(processed_value) + println(" suggested_peak_picking_min_peak_shape_r2 set to $(suggested_peak_picking_min_peak_shape_r2)") + end + elseif step_name == :PeakSelection + if param_key == :min_snr + suggested_peak_selection_min_snr = string(processed_value) + peak_selection_min_snr = string(processed_value) + println(" suggested_peak_selection_min_snr set to $(suggested_peak_selection_min_snr)") + elseif param_key == :min_fwhm_ppm + suggested_peak_selection_min_fwhm_ppm = string(processed_value) + peak_selection_min_fwhm_ppm = string(processed_value) + println(" suggested_peak_selection_min_fwhm_ppm set to $(suggested_peak_selection_min_fwhm_ppm)") + elseif param_key == :max_fwhm_ppm + suggested_peak_selection_max_fwhm_ppm = string(processed_value) + peak_selection_max_fwhm_ppm = string(processed_value) + println(" suggested_peak_selection_max_fwhm_ppm set to $(suggested_peak_selection_max_fwhm_ppm)") + elseif param_key == :min_shape_r2 + suggested_peak_selection_min_shape_r2 = string(processed_value) + peak_selection_min_shape_r2 = string(processed_value) + println(" suggested_peak_selection_min_shape_r2 set to $(suggested_peak_selection_min_shape_r2)") + elseif param_key == :frequency_threshold + suggested_peak_selection_frequency_threshold = string(processed_value) + peak_selection_frequency_threshold = string(processed_value) + println(" suggested_peak_selection_frequency_threshold set to $(suggested_peak_selection_frequency_threshold)") + elseif param_key == :correlation_threshold + suggested_peak_selection_correlation_threshold = string(processed_value) + peak_selection_correlation_threshold = string(processed_value) + println(" suggested_peak_selection_correlation_threshold set to $(suggested_peak_selection_correlation_threshold)") + end + elseif step_name == :PeakBinning + if param_key == :tolerance + suggested_binning_tolerance = string(processed_value) + binning_tolerance = string(processed_value) + println(" suggested_binning_tolerance set to $(suggested_binning_tolerance)") + elseif param_key == :frequency_threshold + suggested_binning_frequency_threshold = string(processed_value) + binning_frequency_threshold = string(processed_value) + println(" suggested_binning_frequency_threshold set to $(suggested_binning_frequency_threshold)") + elseif param_key == :min_peak_per_bin + suggested_binning_min_peak_per_bin = string(processed_value) + binning_min_peak_per_bin = string(processed_value) + println(" suggested_binning_min_peak_per_bin set to $(suggested_binning_min_peak_per_bin)") + elseif param_key == :max_bin_width_ppm + suggested_binning_max_bin_width_ppm = string(processed_value) + binning_max_bin_width_ppm = string(processed_value) + println(" suggested_binning_max_bin_width_ppm set to $(suggested_binning_max_bin_width_ppm)") + elseif param_key == :num_uniform_bins + suggested_binning_num_uniform_bins = string(processed_value) + binning_num_uniform_bins = string(processed_value) + println(" suggested_binning_num_uniform_bins set to $(suggested_binning_num_uniform_bins)") + end + end + end + end + end + + # Also set method types for steps + if haskey(recommended_params, :Smoothing) && haskey(recommended_params[:Smoothing], :method) + smoothing_method = string(recommended_params[:Smoothing][:method]) + end + if haskey(recommended_params, :BaselineCorrection) && haskey(recommended_params[:BaselineCorrection], :method) + baseline_method = string(recommended_params[:BaselineCorrection][:method]) + end + if haskey(recommended_params, :Normalization) && haskey(recommended_params[:Normalization], :method) + normalization_method = string(recommended_params[:Normalization][:method]) + end + if haskey(recommended_params, :PeakAlignment) && haskey(recommended_params[:PeakAlignment], :method) + alignment_method = string(recommended_params[:PeakAlignment][:method]) + end + if haskey(recommended_params, :PeakPicking) && haskey(recommended_params[:PeakPicking], :method) + peak_picking_method = string(recommended_params[:PeakPicking][:method]) + end + if haskey(recommended_params, :PeakBinning) && haskey(recommended_params[:PeakBinning], :method) + binning_method = string(recommended_params[:PeakBinning][:method]) + end + if haskey(recommended_params, :PeakAlignment) && haskey(recommended_params[:PeakAlignment], :tolerance_unit) + alignment_tolerance_unit = string(recommended_params[:PeakAlignment][:tolerance_unit]) + end + if haskey(recommended_params, :PeakBinning) && haskey(recommended_params[:PeakBinning], :tolerance_unit) + binning_tolerance_unit = string(recommended_params[:PeakBinning][:tolerance_unit]) + end + + msg = "File loaded and parameters suggested." + catch e + @warn "Could not suggest parameters. Using defaults. Error: $e" + end + + metadata_columns = [ + Dict("name" => "parameter", "label" => "Parameter", "field" => "parameter", "align" => "left"), + Dict("name" => "value", "label" => "Value", "field" => "value", "align" => "left"), + ] + summary_stats = extract_metadata(loaded_data, local_full_route) + metadata_rows = summary_stats["summary"] + btnMetadataDisable = isempty(metadata_rows) + + w, h = loaded_data.image_dims + imgWidth, imgHeight = w > 0 ? (w, h) : (500, 500) + + update_registry(registry_path, dataset_name, local_full_route, summary_stats, is_imzML) + + # Update folder lists in UI + registry = load_registry(registry_path) + all_folders = sort(collect(keys(registry)), lt=natural) + img_folders = filter(folder -> get(get(registry, folder, Dict()), "is_imzML", false), all_folders) + available_folders = deepcopy(all_folders) + image_available_folders = deepcopy(img_folders) + + selected_folder_main = dataset_name + msi_data = loaded_data + log_memory_usage("Full Load", msi_data) + + eTime = round(time() - sTime, digits=3) + msg = "Active file loaded in $(eTime) seconds. Dataset '$(dataset_name)' is ready for analysis." + + btnStartDisable = false + btnPlotDisable = false + btnSpectraDisable = false + SpectraEnabled = true + + catch e + msi_data = nothing + msg = "Error loading active file: $e" + warning_msg = true + btnStartDisable = true + btnSpectraDisable = true + SpectraEnabled = false + btnMetadataDisable = true + @error "File loading failed" exception=(e, catch_backtrace()) + finally + GC.gc() # Trigger garbage collection + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS + end + progress = false + progressSpectraPlot = false + end + end + + @onbutton export_params_btn begin + params_to_export = Dict( + "stabilization_method" => stabilization_method, + "smoothing_method" => smoothing_method, + "smoothing_window" => smoothing_window, + "smoothing_order" => smoothing_order, + "baseline_method" => baseline_method, + "baseline_iterations" => baseline_iterations, + "baseline_window" => baseline_window, + "normalization_method" => normalization_method, + "alignment_method" => alignment_method, + "alignment_span" => alignment_span, + "alignment_tolerance" => alignment_tolerance, + "alignment_tolerance_unit" => alignment_tolerance_unit, + "alignment_max_shift_ppm" => alignment_max_shift_ppm, + "alignment_min_matched_peaks" => alignment_min_matched_peaks, + "peak_picking_method" => peak_picking_method, + "peak_picking_snr_threshold" => peak_picking_snr_threshold, + "peak_picking_half_window" => peak_picking_half_window, + "peak_picking_min_peak_prominence" => peak_picking_min_peak_prominence, + "peak_picking_merge_peaks_tolerance" => peak_picking_merge_peaks_tolerance, + "peak_picking_min_peak_width_ppm" => peak_picking_min_peak_width_ppm, + "peak_picking_max_peak_width_ppm" => peak_picking_max_peak_width_ppm, + "peak_picking_min_peak_shape_r2" => peak_picking_min_peak_shape_r2, + "binning_method" => binning_method, + "binning_tolerance" => binning_tolerance, + "binning_tolerance_unit" => binning_tolerance_unit, + "binning_frequency_threshold" => binning_frequency_threshold, + "binning_min_peak_per_bin" => binning_min_peak_per_bin, + "binning_max_bin_width_ppm" => binning_max_bin_width_ppm, + "binning_intensity_weighted_centers" => binning_intensity_weighted_centers, + "binning_num_uniform_bins" => binning_num_uniform_bins, + "calibration_fit_order" => calibration_fit_order, + "calibration_ppm_tolerance" => calibration_ppm_tolerance, + "peak_selection_min_snr" => peak_selection_min_snr, + "peak_selection_min_fwhm_ppm" => peak_selection_min_fwhm_ppm, + "peak_selection_max_fwhm_ppm" => peak_selection_max_fwhm_ppm, + "peak_selection_min_shape_r2" => peak_selection_min_shape_r2, + "peak_selection_frequency_threshold" => peak_selection_frequency_threshold, + "peak_selection_correlation_threshold" => peak_selection_correlation_threshold, + "reference_peaks_list" => reference_peaks_list + ) + json_string = JSON.json(params_to_export) + js_script = """ + var element = document.createElement('a'); + element.setAttribute('href', 'data:text/json;charset=utf-8,' + encodeURIComponent(`$json_string`)); + element.setAttribute('download', 'preprocessing_params.json'); + element.style.display = 'none'; + document.body.appendChild(element); + element.click(); + document.body.removeChild(element); + """ + run_js(js_script) + msg = "Parameters exported." + end + + @onchange import_params_btn begin + if import_params_btn + try + json_string = String(imported_params_file.data) + params = JSON.parse(json_string) + + for (key, value) in params + if key == "reference_peaks_list" + reference_peaks_list = value + else + # Use getfield and setproperty! to update reactive variables by name + if hasfield(typeof(@__MODULE__), Symbol(key)) + getfield(@__MODULE__, Symbol(key))[] = value + end + end + end + msg = "Parameters imported successfully." + catch e + msg = "Failed to import parameters: $e" + warning_msg = true + finally + import_params_btn = false # Reset button state + end + end + end + + @onbutton run_full_pipeline begin + progressPrep = true + msg = "Running preprocessing pipeline..." + try + # 1. Get data from "before" plot + if isempty(plotdata_before.traces) || isempty(plotdata_before.traces[1].x) + msg = "Please generate a spectrum plot first (e.g., Mean, Sum, or X,Y)." + warning_msg = true + return + end + + mz_data = plotdata_before.traces[1].x + intensity_data = plotdata_before.traces[1].y + + # 2. Create a temporary MutableSpectrum + temp_spectrum = MutableSpectrum(mz_data, intensity_data, [], 1) # id=1 is arbitrary + + # 3. Run the pipeline + pipeline_spectra = [temp_spectrum] + + if enable_stabilization + current_pipeline_step = "Stabilizing..." + params = Dict(:method => Symbol(stabilization_method)) + apply_intensity_transformation(pipeline_spectra, params) + end + if enable_smoothing + current_pipeline_step = "Smoothing..." + params = Dict(:method => Symbol(smoothing_method), :window => parse(Int, smoothing_window), :order => parse(Int, smoothing_order)) + apply_smoothing(pipeline_spectra, params) + end + if enable_baseline + current_pipeline_step = "Correcting baseline..." + params = Dict(:method => Symbol(baseline_method), :iterations => parse(Int, baseline_iterations), :window => parse(Int, baseline_window)) + apply_baseline_correction(pipeline_spectra, params) + end + if enable_normalization + current_pipeline_step = "Normalizing..." + params = Dict(:method => Symbol(normalization_method)) + apply_normalization(pipeline_spectra, params) + end + if enable_standards && enable_calibration + current_pipeline_step = "Calibrating..." + ref_peaks = Dict(p["mz"] => p["label"] for p in reference_peaks_list) + params = Dict(:ppm_tolerance => parse(Float64, calibration_ppm_tolerance), :fit_order => parse(Int, calibration_fit_order)) + apply_calibration(pipeline_spectra, params) + end + if enable_alignment + # Alignment is a no-op for a single spectrum, but we call it for completeness + current_pipeline_step = "Aligning (skipped for single spectrum)..." + end + if enable_peak_picking + current_pipeline_step = "Picking peaks..." + params = Dict( + :method => Symbol(peak_picking_method), + :snr_threshold => parse(Float64, peak_picking_snr_threshold), + :half_window => parse(Int, peak_picking_half_window), + :min_peak_prominence => parse(Float64, peak_picking_min_peak_prominence), + :merge_peaks_tolerance => parse(Float64, peak_picking_merge_peaks_tolerance) + ) + apply_peak_picking(pipeline_spectra, params) + end + + # 4. Update "After" plot + processed_spectrum = pipeline_spectra[1] + + # Main spectrum trace + after_trace = PlotlyBase.scatter(x=processed_spectrum.mz, y=processed_spectrum.intensity, mode="lines", name="Processed Spectrum") + + traces_after = [after_trace] + + # Add peaks if they exist + if !isempty(processed_spectrum.peaks) + peak_mzs = [p.mz for p in processed_spectrum.peaks] + peak_intensities = [p.intensity for p in processed_spectrum.peaks] + peak_trace = PlotlyBase.scatter(x=peak_mzs, y=peak_intensities, mode="markers", name="Picked Peaks", marker=attr(color="red", size=8)) + push!(traces_after, peak_trace) + end + + plotdata_after = traces_after + plotlayout_after = PlotlyBase.Layout(title="After Preprocessing") + msg = "Pipeline finished." + + catch e + msg = "Error during pipeline execution: $e" + warning_msg = true + @error "Pipeline failed" exception=(e, catch_backtrace()) + finally + progressPrep = false + current_pipeline_step = "" + end + end + + # This new handler correctly adds the file from full_route to the batch list. + @onbutton btnAddBatch begin + if isempty(full_route) || full_route == "unknown (manually added)" + msg = "No active file selected to add to batch." + warning_msg = true + return + end + + if !(full_route in selected_files) + push!(selected_files, full_route) + selected_files = deepcopy(selected_files) # Force reactivity + batch_file_count = length(selected_files) + msg = "File added to batch." + else + msg = "File is already in the batch list." + warning_msg = true + end + end + + @onbutton clear_batch_btn begin + selected_files = String[] + batch_file_count = 0 + msg = "Batch cleared" + end + + @onchange selected_files begin + batch_file_count = length(selected_files) + end + + @onchange full_route begin + if !isempty(full_route) && !(full_route in selected_files) + push!(selected_files, full_route) + selected_files = deepcopy(selected_files) # Force reactivity + batch_file_count = length(selected_files) + msg = "File automatically added to batch" + end + end + + @onbutton showMetadataBtn begin + if !isempty(available_folders) + if !isempty(selected_folder_main) + selected_folder_metadata = selected_folder_main + elseif !isempty(available_folders) + selected_folder_metadata = first(available_folders) + end + showMetadataDialog = true + else + msg = "No processed datasets available." + warning_msg = true + end + end + + @onchange selected_folder_metadata begin + if !isempty(selected_folder_metadata) + registry = load_registry(registry_path) + dataset_info = get(registry, selected_folder_metadata, nothing) + + if dataset_info !== nothing && haskey(dataset_info, "metadata") && !isempty(get(dataset_info["metadata"], "summary", [])) + metadata_rows = dataset_info["metadata"]["summary"] + btnMetadataDisable = false + else + metadata_rows = [] + btnMetadataDisable = true + msg = "Metadata not found in registry for $(selected_folder_metadata)." + end + end + end + + + @onbutton mainProcess @time begin + # --- UI State Update --- + progress = true + btnStartDisable = true + btnPlotDisable = true + btnSpectraDisable = true + overall_progress = 0.0 + progress_message = "Preparing batch process..." + + # --- CAPTURE CURRENT VALUES HERE --- + current_selected_files = selected_files + current_nmass = Nmass + current_tol = Tol + current_color_level = colorLevel + current_triq_enabled = triqEnabled + current_triq_prob = triqProb + current_mfilter_enabled = MFilterEnabled + current_mask_enabled = maskEnabled + current_registry_path = registry_path + + println("starting main process with $(length(current_selected_files)) files") + total_time_start = time() + try + # --- 1. Parameter Validation --- + if isempty(current_selected_files) + progress_message = "No .imzML files in batch. Please add files first." + warning_msg = true + println(progress_message) + return + end + + masses = Float64[] + try + masses = [parse(Float64, strip(m)) for m in split(current_nmass, ',', keepempty=false)] + catch e + progress_message = "Invalid m/z value(s). Please provide a comma-separated list of numbers. Error: $e" + warning_msg = true + return + end + + if isempty(masses) + progress_message = "No valid m/z values found. Please provide comma-separated positive numbers." + warning_msg = true + return + end + + # --- 2. Batch Processing Loop --- + num_files = length(current_selected_files) + total_steps = num_files + current_step = 0 + errors = Dict("load_errors" => String[], "slice_errors" => String[], "io_errors" => String[]) + newly_created_folders = String[] + files_without_mask = 0 + + for (file_idx, file_path) in enumerate(current_selected_files) + progress_message = "Processing file $(file_idx)/$(num_files): $(basename(file_path))" + overall_progress = current_step / total_steps + + all_params = ( + tolerance = current_tol, + colorL = current_color_level, + triqE = current_triq_enabled, + triqP = current_triq_prob, + medianF = current_mfilter_enabled, + registry = current_registry_path, + fileIdx = file_idx, + nFiles = num_files + ) + + success, error_msg = process_file_safely(file_path, masses, all_params, progress_message, overall_progress, use_mask=current_mask_enabled) + + if !success + push!(errors["load_errors"], error_msg) + else + push!(newly_created_folders, replace(basename(file_path), r"\.imzML$"i => "")) + end + current_step += 1 + end + + # --- 3. Final Report --- + total_time_end = round(time() - total_time_start, digits=3) + + registry = load_registry(current_registry_path) + all_folders = sort(collect(keys(registry)), lt=natural) + img_folders = filter(folder -> get(get(registry, folder, Dict()), "is_imzML", false), all_folders) + available_folders = deepcopy(all_folders) + image_available_folders = deepcopy(img_folders) + + if !isempty(newly_created_folders) + selected_folder_main = first(newly_created_folders) + end + + successful_files = length(newly_created_folders) + total_errors = sum(length, values(errors)) + + if total_errors == 0 + msg = "Successfully processed all $(successful_files) file(s) in $(total_time_end) seconds." + else + msg = "Batch completed in $(total_time_end) seconds with $(total_errors) error(s)." + warning_msg = true + end + + mask_summary = current_mask_enabled ? "\nFiles processed without a mask: $(files_without_mask)" : "" + + batch_summary = """ + Processed $(successful_files)/$(num_files) files successfully. + $(mask_summary) + + Errors by category: + • Load failures: $(length(errors["load_errors"])) + • Slice generation: $(length(errors["slice_errors"])) + • I/O issues: $(length(errors["io_errors"])) + + Detailed errors: + $(join(vcat(values(errors)...), "\n")) + """ + showBatchSummary = true + + # Update UI to display the last generated image + if !isempty(newly_created_folders) + timestamp = string(time_ns()) + folder_path = joinpath("public", selected_folder_main) + + if current_triq_enabled + triq_files = filter(filename -> startswith(filename, "TrIQ_") && endswith(filename, ".bmp"), readdir(folder_path)) + col_triq_files = filter(filename -> startswith(filename, "colorbar_TrIQ_") && endswith(filename, ".png"), readdir(folder_path)) + + if !isempty(triq_files) + latest_triq = triq_files[argmax([mtime(joinpath(folder_path, f)) for f in triq_files])] + current_triq = latest_triq + imgIntT = "/$(selected_folder_main)/$(current_triq)?t=$(timestamp)" + plotdataImgT, plotlayoutImgT, _, _ = loadImgPlot(imgIntT) + text_nmass = replace(current_triq, r"TrIQ_|.bmp" => "") + msgtriq = "TrIQ m/z: $(replace(text_nmass, "_" => "."))" + + if !isempty(col_triq_files) + latest_col_triq = col_triq_files[argmax([mtime(joinpath(folder_path, f)) for f in col_triq_files])] + current_col_triq = latest_col_triq + colorbarT = "/$(selected_folder_main)/$(current_col_triq)?t=$(timestamp)" + else + colorbarT = "" + end + selectedTab = "tab1" + end + else # Not TrIQ enabled, display regular MSI image + msi_files = filter(filename -> startswith(filename, "MSI_") && endswith(filename, ".bmp"), readdir(folder_path)) + col_msi_files = filter(filename -> startswith(filename, "colorbar_MSI_") && endswith(filename, ".png"), readdir(folder_path)) + + if !isempty(msi_files) + latest_msi = msi_files[argmax([mtime(joinpath(folder_path, f)) for f in msi_files])] + current_msi = latest_msi + imgInt = "/$(selected_folder_main)/$(current_msi)?t=$(timestamp)" + plotdataImg, plotlayoutImg, _, _ = loadImgPlot(imgInt) + text_nmass = replace(current_msi, r"MSI_|.bmp" => "") + msgimg = "m/z: $(replace(text_nmass, "_" => "."))" + + if !isempty(col_msi_files) + latest_col_msi = col_msi_files[argmax([mtime(joinpath(folder_path, f)) for f in col_msi_files])] + current_col_msi = latest_col_msi + colorbar = "/$(selected_folder_main)/$(current_col_msi)?t=$(timestamp)" + else + colorbar = "" + end + selectedTab = "tab0" + end + end + end + + catch e + println("Error in main process: $e") + msg = "Batch processing failed: $e" + warning_msg = true + @error "Main process failed" exception=(e, catch_backtrace()) + finally + # --- UI State Reset --- + progress = false + btnStartDisable = false + btnPlotDisable = false + btnOpticalDisable = false + btnSpectraDisable = false + SpectraEnabled = true + overall_progress = 0.0 + println("Done") + GC.gc() # Trigger garbage collection + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS + end + end + end + + @onbutton createMeanPlot @time begin + if isempty(selected_folder_main) + msg = "No dataset selected. Please process a file and select a folder first." + warning_msg = true + return + end + + progressSpectraPlot = true + btnPlotDisable = true + btnStartDisable = true + msg = "Loading plot for $(selected_folder_main)..." + + try + sTime = time() + registry = load_registry(registry_path) + entry = registry[selected_folder_main] + target_path = entry["source_path"] + + if target_path == "unknown (manually added)" + msg = "Dataset selected contained no route." + warning_msg = true + return + end + + if msi_data === nothing || full_route != target_path + if msi_data !== nothing + close(msi_data) + end + msg = "Reloading $(basename(target_path)) for analysis..." + full_route = target_path + msi_data = OpenMSIData(target_path) + if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing + raw_min = entry["metadata"]["global_min_mz"] + raw_max = entry["metadata"]["global_max_mz"] + min_val = isa(raw_min, Dict) ? get(raw_min, "value", raw_min) : raw_min + max_val = isa(raw_max, Dict) ? get(raw_max, "value", raw_max) : raw_max + set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val)) + else + precompute_analytics(msi_data) + end + end + + local mask_path_for_plot::Union{String, Nothing} = nothing + if maskEnabled && get(entry, "has_mask", false) + mask_path_for_plot = get(entry, "mask_path", "") + if !isfile(mask_path_for_plot) + @warn "Mask not found for plotting: $(mask_path_for_plot). Plotting without mask." + mask_path_for_plot = nothing + end + end + + plotdata, plotlayout, xSpectraMz, ySpectraMz = meanSpectrumPlot(msi_data, selected_folder_main, mask_path=mask_path_for_plot) + plotdata_before = plotdata + plotlayout_before = plotlayout + selectedTab = "tab2" + fTime = time() + eTime = round(fTime - sTime, digits=3) + msg = "Plot loaded in $(eTime) seconds" + log_memory_usage("Mean Plot Generated", msi_data) + catch e + msg = "Could not generate mean spectrum plot: $e" + warning_msg = true + @error "Mean spectrum plotting failed" exception=(e, catch_backtrace()) + finally + progressSpectraPlot = false + btnPlotDisable = false + btnSpectraDisable = false + btnStartDisable = false + GC.gc() # Trigger garbage collection + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS + end + end + end + + @onbutton createSumPlot @time begin + if isempty(selected_folder_main) + msg = "No dataset selected. Please process a file and select a folder first." + warning_msg = true + return + end + + progressSpectraPlot = true + btnPlotDisable = true + btnStartDisable = true + msg = "Loading total spectrum plot for $(selected_folder_main)..." + + try + sTime = time() + registry = load_registry(registry_path) + entry = registry[selected_folder_main] + target_path = entry["source_path"] + + if target_path == "unknown (manually added)" + msg = "Dataset selected contained no route." + warning_msg = true + return + end + + if msi_data === nothing || full_route != target_path + if msi_data !== nothing + close(msi_data) + end + msg = "Reloading $(basename(target_path)) for analysis..." + full_route = target_path + msi_data = OpenMSIData(target_path) + if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing + raw_min = entry["metadata"]["global_min_mz"] + raw_max = entry["metadata"]["global_max_mz"] + min_val = isa(raw_min, Dict) ? get(raw_min, "value", raw_min) : raw_min + max_val = isa(raw_max, Dict) ? get(raw_max, "value", raw_max) : raw_max + set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val)) + else + precompute_analytics(msi_data) + end + end + local mask_path_for_plot::Union{String, Nothing} = nothing + if maskEnabled && get(entry, "has_mask", false) + mask_path_for_plot = get(entry, "mask_path", "") + if !isfile(mask_path_for_plot) + @warn "Mask not found for plotting: $(mask_path_for_plot). Plotting without mask." + mask_path_for_plot = nothing + end + end + + plotdata, plotlayout, xSpectraMz, ySpectraMz = sumSpectrumPlot(msi_data, selected_folder_main, mask_path=mask_path_for_plot) + plotdata_before = plotdata + plotlayout_before = plotlayout + selectedTab = "tab2" + fTime = time() + eTime = round(fTime - sTime, digits=3) + msg = "Total plot loaded in $(eTime) seconds" + log_memory_usage("Sum Plot Generated", msi_data) + catch e + msg = "Could not generate total spectrum plot: $e" + warning_msg = true + @error "Total spectrum plotting failed" exception=(e, catch_backtrace()) + finally + progressSpectraPlot = false + btnPlotDisable = false + btnSpectraDisable = false + btnStartDisable = false + GC.gc() # Trigger garbage collection + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS + end + end + end + + @onbutton createXYPlot @time begin + if isempty(selected_folder_main) + msg = "No dataset selected. Please process a file and select a folder first." + warning_msg = true + return + end + + progressSpectraPlot = true + btnStartDisable = true + btnPlotDisable = true + btnSpectraDisable = true + msg = "Loading plot for $(selected_folder_main)..." + + try + sTime = time() + registry = load_registry(registry_path) + + if !haskey(registry, selected_folder_main) + msg = "Dataset '$(selected_folder_main)' not found in registry." + warning_msg = true + return + end + + entry = registry[selected_folder_main] + target_path = entry["source_path"] + + if target_path == "unknown (manually added)" + msg = "Dataset selected contained no route." + warning_msg = true + return + end + + if msi_data === nothing || full_route != target_path + if msi_data !== nothing + close(msi_data) + end + msg = "Reloading $(basename(target_path)) for analysis..." + full_route = target_path + msi_data = OpenMSIData(target_path) + if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing + raw_min = entry["metadata"]["global_min_mz"] + raw_max = entry["metadata"]["global_max_mz"] + min_val = isa(raw_min, Dict) ? get(raw_min, "value", raw_min) : raw_min + max_val = isa(raw_max, Dict) ? get(raw_max, "value", raw_max) : raw_max + set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val)) + else + precompute_analytics(msi_data) + end + end + + local mask_path_for_plot::Union{String, Nothing} = nothing + if maskEnabled && get(entry, "has_mask", false) + mask_path_for_plot = get(entry, "mask_path", "") + if !isfile(mask_path_for_plot) + @warn "Mask not found for plotting: $(mask_path_for_plot). Plotting without mask." + mask_path_for_plot = nothing + end + end + + y_positive = yCoord < 0 ? abs(yCoord) : yCoord + plotdata, plotlayout, xSpectraMz, ySpectraMz = xySpectrumPlot(msi_data, xCoord, y_positive, imgWidth, imgHeight, selected_folder_main, mask_path=mask_path_for_plot) + plotdata_before = plotdata + plotlayout_before = plotlayout + + actual_title = if plotlayout.title isa Dict && haskey(plotlayout.title, :text) + plotlayout.title[:text] + elseif plotlayout.title isa Dict && haskey(plotlayout.title, "text") + plotlayout.title["text"] + else + string(plotlayout.title) # Fallback + end + + if occursin("Masked Spectrum at", actual_title) + # Extract coordinates from masked spectrum title + coords_match = match(r"Masked Spectrum at \((\d+), (\d+)\)", actual_title) + if coords_match !== nothing + xCoord = parse(Int, coords_match.captures[1]) + yCoord = -parse(Int, coords_match.captures[2]) # Negative for display + end + elseif occursin("Spectrum at", actual_title) + # Extract coordinates from regular spectrum title + coords_match = match(r"Spectrum at \((\d+), (\d+)\)", actual_title) + if coords_match !== nothing + xCoord = parse(Int, coords_match.captures[1]) + yCoord = -parse(Int, coords_match.captures[2]) # Negative for display + end + else + # For non-imaging data or fallback, just clamp the coordinates + xCoord = clamp(xCoord, 1, imgWidth) + yCoord = yCoord < 0 ? yCoord : -clamp(yCoord, 1, imgHeight) + end + + selectedTab = "tab2" + fTime = time() + eTime = round(fTime - sTime, digits=3) + msg = "Plot loaded in $(eTime) seconds" + log_memory_usage("XY Plot Generated", msi_data) + catch e + msg = "Could not retrieve spectrum: $e" + warning_msg = true + @error "Spectrum plotting failed" exception=(e, catch_backtrace()) + finally + progressSpectraPlot = false + btnPlotDisable = false + btnSpectraDisable = false + btnStartDisable = false + GC.gc() # Trigger garbage collection + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS + end + end + end + + GC.gc() # Trigger garbage collection + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure julia returns the freed memory to OS + end +end +# == Pages == +# Register a new route and the page that will be loaded on access +@page("/", "app.jl.html") +end \ No newline at end of file diff --git a/julia_imzML_visual.jl b/julia_imzML_visual.jl index cedb8d9..0195551 100644 --- a/julia_imzML_visual.jl +++ b/julia_imzML_visual.jl @@ -926,6 +926,27 @@ function update_registry(registry_path, dataset_name, source_path, metadata=noth end end +""" + save_registry(registry_path, registry_data) + +Saves the dataset registry to a JSON file, ensuring thread-safe access. + +# Arguments +- `registry_path`: Path to the `registry.json` file. +- `registry_data`: The dictionary containing the registry data to save. +""" +function save_registry(registry_path, registry_data) + lock(REGISTRY_LOCK) do + try + open(registry_path, "w") do f + JSON.print(f, registry_data, 4) + end + catch e + @error "Failed to write to registry.json: $e" + end + end +end + """ process_file_safely(file_path, masses, params, progress_message_ref, overall_progress_ref) diff --git a/public/css/genieapp.css b/public/css/genieapp.css index 9aa8eed..19bc18c 100644 --- a/public/css/genieapp.css +++ b/public/css/genieapp.css @@ -41,6 +41,38 @@ color: rgb(0, 0, 0) !important; } +/* Placeholder color for custom-standout q-inputs */ +.custom-standout .q-field__native::placeholder { + color: #CFD8DC !important; /* Lighter grey for placeholder text */ + opacity: 1 !important; /* Ensure full visibility */ +} + +.custom-standout .q-field__control::before { + border-color: #7f8389 !important; /* Keep the original border color */ +} + +/* Placeholder color for various browsers */ +.custom-standout input::placeholder { + color: #CFD8DC !important; + opacity: 1 !important; +} +.custom-standout input::-webkit-input-placeholder { /* WebKit, Blink, Edge */ + color: #CFD8DC !important; + opacity: 1 !important; +} +.custom-standout input::-moz-placeholder { /* Mozilla Firefox 19+ */ + color: #CFD8DC !important; + opacity: 1 !important; +} +.custom-standout input:-ms-input-placeholder { /* Internet Explorer 10-11 */ + color: #CFD8DC !important; + opacity: 1 !important; +} +.custom-standout input::-ms-input-placeholder { /* Microsoft Edge */ + color: #CFD8DC !important; + opacity: 1 !important; +} + #tabHeader { color: #009f90; } @@ -82,10 +114,6 @@ font-size: 0.9rem; } -#intDivStyle-left .q-tab-panels, #intDivStyle-right .q-tab-panels { - /* Removed fixed height */ -} - #intDivStyle-left .q-tab-panel, #intDivStyle-right .q-tab-panel { height: auto; /* Let content define height */ overflow-y: hidden; /* Remove scrollbar */ @@ -118,3 +146,24 @@ .q-option-group > div { margin-bottom: 8px; } + +/* Loading Overlay Styles */ +.loading-overlay { + position: fixed; + top: 0; + left: 0; + width: 100%; + height: 100%; + background-color: rgba(0, 0, 0, 0.7); + display: flex; + justify-content: center; + align-items: center; + z-index: 9999; + text-align: center; +} +.loading-content .q-spinner { + color: white; +} +.loading-content .text-h6 { + color: white; +} \ No newline at end of file diff --git a/src/MSI_src.jl b/src/MSI_src.jl index bc3c7ee..3dd7243 100644 --- a/src/MSI_src.jl +++ b/src/MSI_src.jl @@ -30,19 +30,31 @@ export run_preprocessing_analysis, BaselineCorrection, Normalization, PeakPicking, - PeakBinningParams, + PeakBinning, get_masked_spectrum_indices, - detect_peaks_profile, - detect_peaks_centroid, - smooth_spectrum, - apply_baseline_correction, - apply_normalization, - bin_peaks, + detect_peaks_profile_core, + detect_peaks_centroid_core, + smooth_spectrum_core, + apply_baseline_correction_core, + apply_normalization_core, + bin_peaks_core, PeakSelection, PeakAlignment, - find_calibration_peaks, - align_peaks_lowess, - MutableSpectrum + find_calibration_peaks_core, + align_peaks_lowess_core, + MutableSpectrum, + transform_intensity_core, + calibrate_spectra_core + +export apply_baseline_correction, + apply_smoothing, + apply_peak_picking, + apply_calibration, + apply_peak_alignment, + apply_normalization, + apply_peak_binning, + apply_intensity_transformation, + save_feature_matrix # Include all source files directly into the main module include("BloomFilters.jl") @@ -55,6 +67,7 @@ include("MzmlConverter.jl") include("Preprocessing.jl") include("ImageProcessing.jl") include("Precalculations.jl") +include("PreprocessingPipeline.jl") using Setfield # For immutable struct updates diff --git a/src/Precalculations.jl b/src/Precalculations.jl index 849a036..70e28f7 100644 --- a/src/Precalculations.jl +++ b/src/Precalculations.jl @@ -5,10 +5,21 @@ using StatsBase # For mean, std, median, quantile, mad # ============================================================================= """ - find_ppm_error_by_region(msi_data, region_masks, reference_peaks)::Dict + find_ppm_error_by_region(msi_data::MSIData, region_masks::Dict, reference_peaks::Dict) -> Dict -Calculates PPM error statistics for different spatial regions. -`region_masks` is a Dict mapping region names (e.g., :tumor) to BitMatrix masks. +Calculates and reports mass accuracy (PPM error) statistics for different spatial regions +defined by masks. This is useful for identifying spatial variations in calibration. + +# Arguments +- `msi_data::MSIData`: The main MSI data object. +- `region_masks::Dict{Symbol, BitMatrix}`: A dictionary mapping region names (e.g., `:tumor`, `:stroma`) + to `BitMatrix` masks. The dimensions of each mask must match `msi_data.image_dims`. +- `reference_peaks::Dict{Float64, String}`: A dictionary of known reference peaks, mapping + theoretical m/z to a name. + +# Returns +- `Dict{Symbol, NamedTuple}`: A dictionary where keys are region names and values are `NamedTuple`s + containing the mass accuracy report for that region, as generated by `analyze_mass_accuracy`. """ function find_ppm_error_by_region(msi_data::MSIData, region_masks::Dict, reference_peaks::Dict) regional_reports = Dict{Symbol, NamedTuple}() @@ -37,6 +48,28 @@ function find_ppm_error_by_region(msi_data::MSIData, region_masks::Dict, referen return regional_reports end +""" + analyze_mass_accuracy(msi_data, reference_peaks; ...) -> NamedTuple + +Analyzes the mass accuracy for a given subset of spectra by comparing detected peaks +against a list of known reference masses. + +# Arguments +- `msi_data::MSIData`: The main MSI data object. +- `reference_peaks::Dict{Float64, String}`: A dictionary of known reference peaks, mapping + theoretical m/z to a name. +- `spectrum_indices::AbstractVector{Int}`: A vector of indices for the spectra to be analyzed. +- `peak_detection_snr_threshold::Float64`: The Signal-to-Noise ratio threshold to use for + detecting peaks within the spectra. +- `ppm_tolerance_for_matching::Float64`: The tolerance in Parts Per Million (PPM) used to match + a detected peak to a reference peak. + +# Returns +- `NamedTuple`: A report containing summary statistics of the PPM errors found, including: + - `mean_ppm_error`, `median_ppm_error`, `std_ppm_error`, `min_ppm_error`, `max_ppm_error` + - `total_matched_peaks`: The total count of successful matches between detected and reference peaks. + - `total_spectra_analyzed`: The number of spectra processed. +""" function analyze_mass_accuracy( msi_data::MSIData, reference_peaks::Dict{Float64, String}; # m/z => name @@ -57,7 +90,7 @@ function analyze_mass_accuracy( return end - detected_peaks = detect_peaks_profile(mz, intensity; snr_threshold=peak_detection_snr_threshold) + detected_peaks = detect_peaks_profile_core(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 @@ -88,7 +121,8 @@ function analyze_mass_accuracy( min_ppm_error = NaN, max_ppm_error = NaN, total_matched_peaks = 0, - total_spectra_analyzed = total_spectra_processed + total_spectra_analyzed = total_spectra_processed, + ppm_error_distribution = Float64[] ) end @@ -107,7 +141,8 @@ function analyze_mass_accuracy( min_ppm_error = min_err, max_ppm_error = max_err, total_matched_peaks = total_matched_peaks, - total_spectra_analyzed = total_spectra_processed + total_spectra_analyzed = total_spectra_processed, + ppm_error_distribution = all_ppm_errors ) end @@ -118,7 +153,16 @@ end """ calculate_adaptive_bin_tolerance(ppm_error_distribution) -> Float64 -Calculates an appropriate binning tolerance based on observed mass accuracy. +Calculates an appropriate binning tolerance in PPM based on the observed mass accuracy +distribution. The strategy is to set the tolerance to capture the vast majority of +peaks from the same analyte, typically using `mean + 3 * standard_deviation`. + +# Arguments +- `ppm_error_distribution::Vector{Float64}`: A vector of PPM error values from a mass + accuracy analysis. + +# Returns +- `Float64`: The suggested binning tolerance in PPM, capped between 10.0 and 100.0. """ function calculate_adaptive_bin_tolerance(ppm_error_distribution::Vector{Float64}) if isempty(ppm_error_distribution) @@ -148,10 +192,20 @@ function calculate_adaptive_bin_tolerance(ppm_error_distribution::Vector{Float64 end """ - calculate_preprocessing_hints(data::MSIData; sample_size::Int=100)::Dict{Symbol, Any} + calculate_preprocessing_hints(data::MSIData; sample_indices)::Dict{Symbol, Any} -Analyzes a sample of spectra to determine optimal default parameters for -preprocessing steps, returning a dictionary of hints. +Analyzes a sample of spectra to determine initial "hints" for preprocessing parameters. +This function provides quick, data-driven defaults for noise level, SNR, and smoothing. + +# Arguments +- `data::MSIData`: The main MSI data object. +- `sample_indices::AbstractVector{Int}`: The indices of spectra to sample for the analysis. + +# Returns +- `Dict{Symbol, Any}`: A dictionary of hints, including: + - `:estimated_noise`: The mean noise level estimated using Median Absolute Deviation (MAD). + - `:suggested_snr`: A default SNR threshold (typically 3.0). + - `:suggested_smoothing_window`: A suggested window size for smoothing, based on instrument resolution if available. """ 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...") @@ -218,9 +272,22 @@ function calculate_preprocessing_hints(data::MSIData; sample_indices::AbstractVe end """ - analyze_instrument_characteristics(msi_data::MSIData)::Dict + analyze_instrument_characteristics(msi_data::MSIData; sample_indices)::Dict -Analyzes instrument metadata and data characteristics to determine acquisition properties. +Analyzes instrument metadata and spectral data to infer key acquisition properties. +It combines information from the `msi_data.instrument_metadata` with direct analysis +of the spectra. + +# Arguments +- `msi_data::MSIData`: The main MSI data object. +- `sample_indices::AbstractVector{Int}`: The indices of spectra to sample for the analysis. + +# Returns +- `Dict{Symbol, Any}`: A dictionary summarizing instrument characteristics: + - `:acquisition_mode`: Inferred as `:profile`, `:centroid`, or `:mixed`. + - `:mz_axis_type`: Inferred as `:regular` or `:irregular` based on m/z step consistency. + - `:dynamic_range`: An estimate of the intensity dynamic range in orders of magnitude. + - Other fields from `instrument_metadata` like `:resolution`, `:polarity`, etc. """ function analyze_instrument_characteristics(msi_data::MSIData; sample_indices::AbstractVector{Int})::Dict results = Dict{Symbol, Any}() @@ -337,9 +404,20 @@ function analyze_instrument_characteristics(msi_data::MSIData; sample_indices::A end """ - analyze_signal_quality(msi_data::MSIData; sample_size::Int=100)::Dict + analyze_signal_quality(msi_data::MSIData; sample_indices)::Dict -Analyzes noise characteristics, signal-to-noise ratios, and overall signal quality. +Analyzes a sample of spectra to assess signal quality, including noise levels, +Signal-to-Noise Ratio (SNR), and Total Ion Current (TIC) variation. + +# Arguments +- `msi_data::MSIData`: The main MSI data object. +- `sample_indices::AbstractVector{Int}`: The indices of spectra to sample for the analysis. + +# Returns +- `Dict{Symbol, Any}`: A dictionary of signal quality metrics: + - `:noise_mean`, `:noise_std`, `:noise_cv`: Statistics of the noise level. + - `:snr_mean`, `:snr_median`, `:snr_95th`: Distribution statistics of the estimated SNR. + - `:tic_mean`, `:tic_std`, `:tic_cv`: Statistics of the Total Ion Current. """ function analyze_signal_quality(msi_data::MSIData; sample_indices::AbstractVector{Int})::Dict results = Dict{Symbol, Any}() @@ -415,9 +493,21 @@ function analyze_signal_quality(msi_data::MSIData; sample_indices::AbstractVecto end """ - analyze_mass_accuracy_global(msi_data::MSIData, reference_peaks::Dict, sample_size::Int)::Dict + analyze_mass_accuracy_global(msi_data, reference_peaks; spectrum_indices)::Dict -Analyzes mass accuracy across the dataset using reference peaks. +Performs a global mass accuracy analysis across a sample of spectra and suggests an +adaptive binning tolerance. + +# Arguments +- `msi_data::MSIData`: The main MSI data object. +- `reference_peaks::Dict`: A dictionary of known reference peaks. +- `spectrum_indices::AbstractVector{Int}`: The indices of spectra to sample for the analysis. + +# Returns +- `Dict{Symbol, Any}`: A dictionary containing: + - `:global_accuracy`: The `NamedTuple` report from `analyze_mass_accuracy`. + - `:suggested_bin_tolerance`: An adaptive tolerance in PPM for peak binning, derived + from the mass accuracy results. """ function analyze_mass_accuracy_global(msi_data::MSIData, reference_peaks::Dict; spectrum_indices::AbstractVector{Int})::Dict @@ -431,7 +521,7 @@ function analyze_mass_accuracy_global(msi_data::MSIData, reference_peaks::Dict; 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 + total_spectra_analyzed = 0, ppm_error_distribution = Float64[] ) results[:global_accuracy] = empty_report results[:suggested_bin_tolerance] = 20.0 # Default @@ -444,7 +534,7 @@ function analyze_mass_accuracy_global(msi_data::MSIData, reference_peaks::Dict; 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 + accuracy_report.ppm_error_distribution ) println(" - Mean PPM error: $(round(accuracy_report.mean_ppm_error, digits=2))") @@ -454,9 +544,21 @@ function analyze_mass_accuracy_global(msi_data::MSIData, reference_peaks::Dict; end """ - analyze_spatial_regions(msi_data::MSIData, region_masks::Dict, reference_peaks::Dict)::Dict + analyze_spatial_regions(msi_data, region_masks, reference_peaks)::Dict -Analyzes different spatial regions for variations in signal quality and mass accuracy. +Analyzes different spatial regions for variations in mass accuracy. This function is a +wrapper around `find_ppm_error_by_region` and summarizes the results. + +# Arguments +- `msi_data::MSIData`: The main MSI data object. +- `region_masks::Dict`: A dictionary of named `BitMatrix` masks for each region. +- `reference_peaks::Dict`: A dictionary of known reference peaks. + +# Returns +- `Dict{Symbol, Any}`: A dictionary containing: + - `:regional_ppm_errors`: A dictionary mapping each region name to its mass accuracy report. + - `:max_regional_ppm_difference`: The difference between the highest and lowest mean PPM + error across all analyzed regions. """ function analyze_spatial_regions(msi_data::MSIData, region_masks::Dict, reference_peaks::Dict)::Dict results = Dict{Symbol, Any}() @@ -484,9 +586,24 @@ function analyze_spatial_regions(msi_data::MSIData, region_masks::Dict, referenc end """ - analyze_peak_characteristics(msi_data::MSIData, sample_size::Int)::Dict + analyze_peak_characteristics(msi_data, instrument_analysis, mass_accuracy_results; spectrum_indices)::Dict -Analyzes peak shape, width, and quality characteristics. +Analyzes peak shape, width, and quality from a sample of spectra. The behavior +adapts based on whether the data is in `:profile` or `:centroid` mode. + +# Arguments +- `msi_data::MSIData`: The main MSI data object. +- `instrument_analysis::Dict`: The output from `analyze_instrument_characteristics`. +- `mass_accuracy_results`: The output from `analyze_mass_accuracy_global`. +- `spectrum_indices::AbstractVector{Int}`: The indices of spectra to sample. + +# Returns +- `Dict{Symbol, Any}`: A dictionary of peak metrics: + - `:mean_fwhm_ppm`, `:median_fwhm_ppm`: Statistics of Full Width at Half Maximum (FWHM). + For centroid data, this is estimated from mass accuracy. + - `:mean_gaussian_r2`: The average goodness-of-fit to a Gaussian shape (profile data only). + - `:peak_resolution_estimate`: An estimate of instrument resolution based on FWHM. + - `:mean_peaks_per_spectrum`: The average number of peaks detected per spectrum. """ function analyze_peak_characteristics(msi_data::MSIData, instrument_analysis::Dict, mass_accuracy_results; spectrum_indices::AbstractVector{Int})::Dict @@ -508,12 +625,21 @@ function analyze_peak_characteristics(msi_data::MSIData, instrument_analysis::Di 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 + if acquisition_mode == :centroid + # Much more permissive defaults for centroid data + results[:mean_fwhm_ppm] = 50.0 + results[:median_fwhm_ppm] = 50.0 + results[:mean_gaussian_r2] = 0.0 # Disable shape filtering for centroids + results[:peak_resolution_estimate] = 20000.0 + results[:mean_peaks_per_spectrum] = 1000 + else + 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 + end return results end @@ -537,7 +663,7 @@ function analyze_peak_characteristics(msi_data::MSIData, instrument_analysis::Di 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) + peaks = detect_peaks_profile_core(mz, intensity; snr_threshold=2.0) push!(peak_counts, length(peaks)) if !isempty(peaks) @@ -600,12 +726,10 @@ function analyze_peak_characteristics(msi_data::MSIData, instrument_analysis::Di 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 @@ -619,16 +743,16 @@ function analyze_peak_characteristics(msi_data::MSIData, instrument_analysis::Di 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 + # For centroid data, use much more permissive parameters + # Don't estimate FWHM from mass accuracy - use reasonable defaults + results[:mean_fwhm_ppm] = 50.0 # Reasonable default for centroid data + results[:median_fwhm_ppm] = 50.0 + results[:mean_gaussian_r2] = 0.0 # Disable shape filtering for centroids + results[:peak_resolution_estimate] = 20000.0 # Reasonable estimate 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))") + println(" - Using permissive FWHM for centroid data: 50.0 ppm") + println(" - Shape filtering disabled for centroid data") end # Common prints @@ -641,7 +765,16 @@ end """ generate_preprocessing_recommendations(analysis_results::Dict)::Dict{Symbol, Any} -Generates intelligent preprocessing recommendations based on the analysis results. +Generates intelligent preprocessing recommendations by synthesizing the results from +various analysis functions (`analyze_instrument_characteristics`, `analyze_signal_quality`, etc.). + +# Arguments +- `analysis_results::Dict`: A dictionary containing the comprehensive analysis results from + the `run_preprocessing_analysis` pipeline. + +# Returns +- `Dict{Symbol, Any}`: A dictionary where keys are preprocessing step names (e.g., `:smoothing`, + `:peak_picking`) and values are dictionaries of recommended parameters for that step. """ function generate_preprocessing_recommendations(analysis_results::Dict)::Dict{Symbol, Any} recommendations = Dict{Symbol, Any}() @@ -651,14 +784,17 @@ function generate_preprocessing_recommendations(analysis_results::Dict)::Dict{Sy mass_accuracy = get(analysis_results, :mass_accuracy, Dict()) peak_analysis = get(analysis_results, :peak_analysis, Dict()) + # Stabilization Recommendations + recommendations[:stabilization] = generate_stabilization_recommendations(signal_analysis) + # 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) + # Peak Picking Recommendations - pass instrument analysis + recommendations[:peak_picking] = generate_peak_picking_recommendations(signal_analysis, peak_analysis, inst_analysis) # Normalization Recommendations recommendations[:normalization] = generate_normalization_recommendations(signal_analysis) @@ -672,6 +808,7 @@ function generate_preprocessing_recommendations(analysis_results::Dict)::Dict{Sy return recommendations end +"""Generate baseline correction recommendations based on data properties.""" function generate_baseline_recommendations(inst_analysis, signal_analysis) recommendations = Dict{Symbol, Any}() @@ -691,6 +828,7 @@ function generate_baseline_recommendations(inst_analysis, signal_analysis) return recommendations end +"""Generate smoothing recommendations based on peak width and m/z step.""" function generate_smoothing_recommendations(inst_analysis, peak_analysis) recommendations = Dict{Symbol, Any}() @@ -713,19 +851,33 @@ function generate_smoothing_recommendations(inst_analysis, peak_analysis) return recommendations end -function generate_peak_picking_recommendations(signal_analysis, peak_analysis) +"""Generate peak picking recommendations from signal and peak analyses.""" +function generate_peak_picking_recommendations(signal_analysis, peak_analysis, inst_analysis) recommendations = Dict{Symbol, Any}() - snr_threshold = get(signal_analysis, :suggested_snr, 3.0) - fwhm_ppm = get(peak_analysis, :mean_fwhm_ppm, 20.0) + acquisition_mode = get(inst_analysis, :acquisition_mode, :profile) - 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 + if acquisition_mode == :centroid + # Much more permissive parameters for centroid data + recommendations[:snr_threshold] = 2.0 # Lower threshold for centroid + recommendations[:min_peak_width_ppm] = 0.0 # No minimum width for centroids + recommendations[:max_peak_width_ppm] = 200.0 # Very wide maximum for centroids + recommendations[:reason] = "Centroid data: using permissive parameters" + else + # Existing profile mode logic + 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 + recommendations[:max_peak_width_ppm] = fwhm_ppm * 3.0 + recommendations[:reason] = "Profile data: using standard parameters" + end return recommendations end +"""Generate normalization recommendations based on TIC variation.""" function generate_normalization_recommendations(signal_analysis) recommendations = Dict{Symbol, Any}() @@ -742,6 +894,7 @@ function generate_normalization_recommendations(signal_analysis) return recommendations end +"""Generate alignment recommendations based on calibration status and mass error.""" function generate_alignment_recommendations(mass_accuracy, inst_analysis) recommendations = Dict{Symbol, Any}() @@ -765,6 +918,7 @@ function generate_alignment_recommendations(mass_accuracy, inst_analysis) return recommendations end +"""Generate binning recommendations based on peak width and mass accuracy.""" function generate_binning_recommendations(peak_analysis, mass_accuracy) recommendations = Dict{Symbol, Any}() @@ -794,18 +948,47 @@ function generate_binning_recommendations(peak_analysis, mass_accuracy) return recommendations end +"""Generate intensity stabilization recommendations.""" +function generate_stabilization_recommendations(signal_analysis) + recommendations = Dict{Symbol, Any}() + # Default to sqrt, as it's a common and generally robust transformation. + # More advanced logic could analyze intensity distribution skewness if needed. + recommendations[:method] = :sqrt + 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} + run_preprocessing_analysis(msi_data; ...) Runs a comprehensive pre-analysis pipeline to determine optimal preprocessing parameters. -This function analyzes the dataset to provide intelligent defaults for preprocessing steps. +This function orchestrates a series of analysis steps on a sample of the dataset to +provide intelligent defaults for a full preprocessing workflow. + +The pipeline consists of several phases: +1. **Instrument & Data Characteristics**: Infers acquisition mode, m/z axis type, etc. +2. **Noise & Signal Quality**: Estimates noise, SNR, and TIC variation. +3. **Mass Accuracy**: Calculates PPM error against reference peaks (if provided). +4. **Spatial Regions**: Analyzes regional variations (if masks are provided). +5. **Peak Characteristics**: Measures peak width, shape, and density. +6. **Recommendations**: Synthesizes all analysis results into actionable parameter suggestions. + +# Arguments +- `msi_data::MSIData`: The main MSI data object. +- `reference_peaks::Dict`: Optional. Known m/z values for mass accuracy analysis. +- `region_masks::Dict`: Optional. Named `BitMatrix` masks for regional analysis. +- `sample_size::Int`: The number of spectra to sample for the analysis. +- `mask_path::String`: Optional path to a PNG mask to restrict analysis to a specific ROI. +- `spectrum_indices::AbstractVector{Int}`: Optional vector of indices to restrict analysis to, + overriding `mask_path` and `sample_size` for selection. + +# Returns +- `Dict{Symbol, Any}`: A nested dictionary containing the results of each analysis phase + and a final `:recommendations` dictionary. The recommendations are also stored in + `msi_data.preprocessing_hints`. """ function run_preprocessing_analysis(msi_data::MSIData; reference_peaks::Dict{Float64, String}=Dict{Float64, String}(), @@ -1019,7 +1202,7 @@ function calculate_resolution_fwhm(mz::Real, profile_mz::AbstractVector{<:Real}, return fwhm > 0 ? Float64(mz) / fwhm : NaN end -# Helper functions for FWHM calculation +"""Helper to find the last index in a vector with a value below a threshold, used for FWHM.""" function find_last_below(v::AbstractVector{<:Real}, threshold::Real) for i in length(v):-1:2 if v[i] >= threshold && v[i-1] < threshold @@ -1029,6 +1212,7 @@ function find_last_below(v::AbstractVector{<:Real}, threshold::Real) return 0 end +"""Helper to find the first index in a vector with a value below a threshold, used for FWHM.""" function find_first_below(v::AbstractVector{<:Real}, threshold::Real) for i in 1:(length(v)-1) if v[i] >= threshold && v[i+1] < threshold @@ -1270,10 +1454,11 @@ function main_precalculation(msi_data::MSIData; 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()) + mass_accuracy = get(analysis_results, :mass_accuracy, nothing) # Can be nothing inst_analysis = get(analysis_results, :instrument_analysis, Dict()) # Initialize parameter dictionaries for each step + stab_params = Dict{Symbol, Any}() cal_params = Dict{Symbol, Any}() sm_params = Dict{Symbol, Any}() bc_params = Dict{Symbol, Any}() @@ -1284,11 +1469,22 @@ function main_precalculation(msi_data::MSIData; pb_params = Dict{Symbol, Any}() # --- Populate Parameters for each step --- + # Stabilization + if !isempty(recs) && haskey(recs, :stabilization) + stab_rec = recs[:stabilization] + stab_params[:method] = get(stab_rec, :method, :sqrt) + else + stab_params[:method] = :sqrt # Default + end # 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) + mean_ppm_error = NaN + suggested_bin_tol = NaN + if mass_accuracy !== nothing + mean_ppm_error = get(get(mass_accuracy, :global_accuracy, Dict()), :mean_ppm_error, NaN) + suggested_bin_tol = get(mass_accuracy, :suggested_bin_tolerance, NaN) + end if !isempty(recs) && haskey(recs, :alignment) && get(recs[:alignment], :required, false) calibration_required = true @@ -1400,89 +1596,114 @@ function main_precalculation(msi_data::MSIData; end # Peak Picking - # Robust method selection with fallback - acquisition_mode = get(inst_analysis, :acquisition_mode, :profile) # Default to profile + # Robust method selection based on acquisition mode + acquisition_mode = get(inst_analysis, :acquisition_mode, :unknown) 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 acquisition_mode == :centroid + # Much more permissive parameters for centroid data + pp_params[:snr_threshold] = 2.0 # Lower for centroid + pp_params[:min_peak_width_ppm] = 0.0 # No minimum width + pp_params[:max_peak_width_ppm] = 200.0 # Very wide maximum + pp_params[:min_peak_shape_r2] = 0.0 # Disable shape filtering + + # Much lower prominence threshold for centroid + estimated_noise = get(signal_analysis, :noise_mean, NaN) + if isfinite(estimated_noise) + pp_params[:min_peak_prominence] = max(estimated_noise * 0.5, 0.001) # Much lower + else + pp_params[:min_peak_prominence] = 0.001 # Very permissive + end + + pp_params[:half_window] = 2 # Smaller window for centroid data + pp_params[:merge_peaks_tolerance] = 10.0 # More permissive merging + else # Profile mode logic + if !isempty(recs) && haskey(recs, :peak_picking) + pk_rec = recs[:peak_picking] + pp_params[:snr_threshold] = get(pk_rec, :snr_threshold, 2.0) # Default to 2.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] = 2.0 # Default to 2.0 + pp_params[:min_peak_width_ppm] = nothing + pp_params[:max_peak_width_ppm] = nothing + end + + # Robust prominence calculation with safety cap for profile + estimated_noise = get(signal_analysis, :noise_mean, NaN) + if isfinite(estimated_noise) + calculated_prominence = round(estimated_noise * 2, digits=4) + pp_params[:min_peak_prominence] = min(calculated_prominence, 0.005) + else + pp_params[:min_peak_prominence] = 0.005 + 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 + 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 + # Robust half_window calculation with safety floor for 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 + window_points = fwhm_mz / avg_mz_step + calculated_half_window = ceil(Int, window_points / 2) + pp_params[:half_window] = max(calculated_half_window, 3) + else + pp_params[:half_window] = 5 + end + + # Robust R^2 calculation with floor for profile + mean_r2 = get(peak_analysis, :mean_gaussian_r2, NaN) + if isfinite(mean_r2) + pp_params[:min_peak_shape_r2] = round(max(0.0, mean_r2 * 0.8), digits=2) + else + pp_params[:min_peak_shape_r2] = 0.0 + end 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) + if acquisition_mode == :centroid + # Much more permissive parameters for centroid data + ps_params[:min_snr] = 1.5 # Even lower than picking threshold + ps_params[:min_fwhm_ppm] = 0.0 # No minimum + ps_params[:max_fwhm_ppm] = 500.0 # Very wide maximum + ps_params[:min_shape_r2] = 0.0 # Disable shape filtering + ps_params[:frequency_threshold] = nothing # User input dependent / Hard to determine + ps_params[:correlation_threshold] = nothing # Hard to determine + else # Profile mode logic + 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 + mean_r2 = get(peak_analysis, :mean_gaussian_r2, NaN) + if isfinite(mean_r2) + ps_params[:min_shape_r2] = round(max(0.0, mean_r2 * 0.8), digits=2) + else + ps_params[:min_shape_r2] = 0.0 # Adjusted fallback + end 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 + ps_params[:frequency_threshold] = nothing + ps_params[:correlation_threshold] = nothing end - ps_params[:frequency_threshold] = nothing # User input dependent / Hard to determine - ps_params[:correlation_threshold] = nothing # Hard to determine # Peak Binning @@ -1525,6 +1746,7 @@ function main_precalculation(msi_data::MSIData; end return Dict( + :Stabilization => stab_params, :Calibration => cal_params, :Smoothing => sm_params, :BaselineCorrection => bc_params, @@ -1532,6 +1754,6 @@ function main_precalculation(msi_data::MSIData; :PeakPicking => pp_params, :PeakAlignment => pa_params, :PeakSelection => ps_params, - :PeakBinningParams => pb_params + :PeakBinning => pb_params ) end diff --git a/src/Preprocessing.jl b/src/Preprocessing.jl index 9200870..95cb15d 100644 --- a/src/Preprocessing.jl +++ b/src/Preprocessing.jl @@ -21,6 +21,7 @@ using ContinuousWavelets # For CWT peak detection using ImageFiltering # For localmaxima in detect_peaks_wavelet using Interpolations # For calibration using Loess # For robust peak alignment +using Base.Threads # For multithreading in apply functions # ============================================================================= # Data Structures @@ -30,6 +31,13 @@ using Loess # For robust peak alignment FeatureMatrix A struct to hold the final feature matrix generated from the preprocessing pipeline. + +# Fields +- `matrix::Array{Float64,2}`: The feature matrix where rows correspond to samples (spectra) + and columns correspond to features (m/z bins). +- `mz_bins::Vector{Tuple{Float64,Float64}}`: A vector of tuples defining the start and + end m/z for each bin (column) in the `matrix`. +- `sample_ids::Vector{Int}`: A vector of identifiers for each sample (row) in the `matrix`. """ struct FeatureMatrix matrix::Array{Float64,2} @@ -38,9 +46,18 @@ struct FeatureMatrix end """ + MutableSpectrum + 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. + +# Fields +- `id::Int`: A unique identifier for the spectrum. +- `mz::AbstractVector{Float64}`: The m/z values of the spectrum. +- `intensity::AbstractVector{Float64}`: The intensity values corresponding to the m/z values. +- `peaks::Vector{NamedTuple}`: A vector of detected peaks, each a `NamedTuple` with fields + like `:mz`, `:intensity`, `:fwhm`, etc. """ mutable struct MutableSpectrum id::Int @@ -62,8 +79,18 @@ abstract type AbstractPreprocessingStep end """ 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. +A preprocessing step for mass calibration. Corrects systematic mass errors in the m/z axis. +Set a parameter to `nothing` to use an auto-determined value from the data where applicable. + +# Arguments +- `method::Symbol`: The calibration method. Currently supports `:internal_standards`. +- `internal_standards::Union{Dict{Float64, String}, Nothing}`: A dictionary mapping theoretical + m/z values of internal standards to their names. +- `base_peak_mz_references::Union{Vector{Float64}, Nothing}`: A vector of reference m/z values + for base peak calibration (not yet implemented). +- `ppm_tolerance::Union{Float64, Nothing}`: The tolerance in parts-per-million (ppm) for matching + peaks to internal standards. +- `fit_order::Int`: The polynomial order for the calibration fit (e.g., 1 for linear, 2 for quadratic). """ struct Calibration <: AbstractPreprocessingStep method::Symbol @@ -80,7 +107,16 @@ end """ BaselineCorrection(; method=:snip, ...) -A preprocessing step for baseline correction. +A preprocessing step for baseline correction. This step estimates and subtracts the +background noise (baseline) from the spectral intensities. + +# Arguments +- `method::Symbol`: The algorithm to use. Options include `:snip` (Sensitive Nonlinear + Iterative Peak clipping), `:convex_hull`, and `:median`. +- `iterations::Union{Int, Nothing}`: The number of iterations for the SNIP algorithm. + A higher number results in a more aggressive baseline. +- `window::Union{Int, Nothing}`: The window size for the `:median` method, determining + the local region for median calculation. """ struct BaselineCorrection <: AbstractPreprocessingStep method::Symbol @@ -95,7 +131,15 @@ end """ Smoothing(; method=:savitzky_golay, ...) -A preprocessing step for spectral smoothing. +A preprocessing step for spectral smoothing. This helps to reduce high-frequency noise +in the intensity data. + +# Arguments +- `method::Symbol`: The smoothing algorithm. Options are `:savitzky_golay` and `:moving_average`. +- `window::Union{Int, Nothing}`: The size of the smoothing window. For Savitzky-Golay, + this must be an odd integer. +- `order::Union{Int, Nothing}`: The polynomial order for the Savitzky-Golay filter. Must be + less than the window size. """ struct Smoothing <: AbstractPreprocessingStep method::Symbol @@ -112,7 +156,15 @@ end """ Normalization(; method=:tic) -A preprocessing step for intensity normalization. +A preprocessing step for intensity normalization. This corrects for variations in total +ion current between different spectra, making them more comparable. + +# Arguments +- `method::Symbol`: The normalization method. Options include: + - `:tic`: Total Ion Current normalization (divides by the sum of intensities). + - `:median`: Divides by the median intensity. + - `:rms`: Root Mean Square normalization. + - `:none`: No normalization is applied. """ struct Normalization <: AbstractPreprocessingStep method::Symbol @@ -125,7 +177,25 @@ end """ PeakPicking(; method=nothing, ...) -A preprocessing step for peak detection. +A preprocessing step for peak detection. This step identifies peaks (signals of interest) +in the profile or centroided spectra. + +# Arguments +- `method::Union{Symbol, Nothing}`: The peak detection algorithm. + - `:profile`: For profile-mode data, using local maxima and quality filters. + - `:wavelet`: Continuous Wavelet Transform (CWT) based peak detection. + - `:centroid`: For centroid-mode data, essentially a filtering step. +- `snr_threshold::Union{Float64, Nothing}`: Signal-to-Noise Ratio threshold. Peaks with SNR + below this value are discarded. +- `half_window::Union{Int, Nothing}`: The number of data points to the left and right of a + potential peak to consider for local maximum detection (for `:profile`). +- `min_peak_prominence::Union{Float64, Nothing}`: The minimum required prominence of a peak, + expressed as a fraction of its height. +- `merge_peaks_tolerance::Union{Float64, Nothing}`: The m/z tolerance within which to merge + adjacent peaks, keeping the more intense one. +- `min_peak_width_ppm, max_peak_width_ppm`: Minimum and maximum acceptable peak width (FWHM) in ppm. +- `min_peak_shape_r2`: Minimum R-squared value from a Gaussian fit to the peak, used as a + quality measure for peak shape. """ struct PeakPicking <: AbstractPreprocessingStep method::Union{Symbol, Nothing} # :profile, :wavelet, :centroid @@ -145,7 +215,20 @@ end """ PeakAlignment(; method=:lowess, ...) -A preprocessing step for peak alignment. +A preprocessing step for peak alignment. This corrects for m/z shifts between spectra, +ensuring that the same analyte peak appears at the same m/z across all samples. + +# Arguments +- `method::Symbol`: The alignment algorithm. Options: `:lowess`, `:linear`, `:ransac`. +- `span::Union{Float64, Nothing}`: The span parameter for LOWESS regression, controlling smoothness. +- `tolerance::Union{Float64, Nothing}`: The tolerance for matching peaks between the target + and reference spectrum. +- `tolerance_unit::Union{Symbol, Nothing}`: The unit for `tolerance`, either `:mz` (absolute) + or `:ppm` (relative). +- `max_shift_ppm::Union{Float64, Nothing}`: The maximum allowed m/z shift in ppm to prevent + spurious peak matches. +- `min_matched_peaks::Union{Int, Nothing}`: The minimum number of matching peaks required + to perform the alignment. """ struct PeakAlignment <: AbstractPreprocessingStep method::Symbol @@ -163,7 +246,19 @@ end """ PeakSelection(; frequency_threshold=nothing, ...) -A preprocessing step for peak selection (filtering). +A preprocessing step for peak selection (filtering). After peak detection, this step +filters the detected peaks based on various quality criteria to remove noise and +irrelevant signals. + +# Arguments +- `frequency_threshold::Union{Float64, Nothing}`: The minimum fraction of spectra in which a + peak must be present to be kept. +- `min_snr::Union{Float64, Nothing}`: Minimum Signal-to-Noise Ratio. +- `min_fwhm_ppm, max_fwhm_ppm`: Minimum and maximum Full Width at Half Maximum in ppm. +- `min_shape_r2::Union{Float64, Nothing}`: Minimum R-squared value from a Gaussian fit, + filtering for good peak shape. +- `correlation_threshold::Union{Float64, Nothing}`: Minimum correlation with neighboring + peaks (not yet implemented). """ struct PeakSelection <: AbstractPreprocessingStep frequency_threshold::Union{Float64, Nothing} @@ -179,11 +274,24 @@ struct PeakSelection <: AbstractPreprocessingStep end """ - PeakBinningParams(; method=:adaptive, ...) + PeakBinning(; method=:adaptive, ...) -A preprocessing step for peak binning. +A preprocessing step for peak binning. This step groups peaks from all spectra into +common m/z bins to generate a feature matrix. + +# Arguments +- `method::Symbol`: The binning strategy. + - `:adaptive`: Creates bins based on the density of detected peaks. + - `:uniform`: Creates a fixed number of equally spaced bins over the m/z range. +- `tolerance, tolerance_unit`: Tolerance for grouping peaks into a bin in `:adaptive` mode. +- `frequency_threshold`: The minimum fraction of spectra a bin must contain a peak in to be kept. +- `min_peak_per_bin`: The minimum number of individual peaks required to form a bin in `:adaptive` mode. +- `max_bin_width_ppm`: Maximum width of a bin in ppm for `:adaptive` mode. +- `intensity_weighted_centers`: If `true`, calculates bin centers as an intensity-weighted + average of the peaks within it. +- `num_uniform_bins`: The number of bins to create for the `:uniform` method. """ -struct PeakBinningParams <: AbstractPreprocessingStep +struct PeakBinning <: AbstractPreprocessingStep method::Symbol tolerance::Union{Float64, Nothing} tolerance_unit::Union{Symbol, Nothing} @@ -193,7 +301,7 @@ struct PeakBinningParams <: AbstractPreprocessingStep 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) + function PeakBinning(; 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 @@ -237,7 +345,7 @@ Checks include: - `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). +- `m/z` values are monotonically non-decreasing. """ function validate_spectrum(mz::AbstractVector{<:Real}, intensity::AbstractVector{<:Real})::Bool # 1. Check for empty vectors @@ -278,11 +386,22 @@ end # ============================================================================= """ - transform_intensity(intensity; method=:sqrt) -> Vector + transform_intensity_core(intensity; method=:sqrt) -> Vector -Applies a variance-stabilizing transformation to the intensity vector. +Applies a variance-stabilizing transformation to the intensity vector. This can +help to make the variance of the signal more constant across the intensity range, +which is often an assumption of downstream statistical methods. + +# Arguments +- `intensity::AbstractVector{<:Real}`: The input intensity values. +- `method::Symbol`: The transformation to apply. Options are: + - `:sqrt`: Square root transformation. + - `:log`: Natural log transformation. + - `:log2`: Base-2 log transformation. + - `:log10`: Base-10 log transformation. + - `:log1p`: Natural log of `1 + x`, useful for data with zeros. """ -function transform_intensity(intensity::AbstractVector{<:Real}; method::Symbol=:sqrt) +function transform_intensity_core(intensity::AbstractVector{<:Real}; method::Symbol=:sqrt) if method === :sqrt return sqrt.(max.(zero(eltype(intensity)), intensity)) elseif method === :log1p @@ -299,7 +418,7 @@ function transform_intensity(intensity::AbstractVector{<:Real}; method::Symbol=: end """ - smooth_spectrum(y::AbstractVector{<:Real}; method::Symbol=:savitzky_golay, window::Int=9, order::Int=2) -> Vector + smooth_spectrum_core(y::AbstractVector{<:Real}; method::Symbol=:savitzky_golay, window::Int=9, order::Int=2) -> Vector Applies a smoothing filter to the intensity data. @@ -309,7 +428,7 @@ Applies a smoothing filter to the intensity data. - `window`: The window size for the filter. - `order`: The polynomial order for Savitzky-Golay """ -function smooth_spectrum(y::AbstractVector{<:Real}; method::Symbol=:savitzky_golay, window::Int=9, order::Int=2) +function smooth_spectrum_core(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 @@ -399,7 +518,9 @@ end """ convex_hull_baseline(y) -> Vector -Estimates the baseline of a spectrum using the convex hull algorithm. +Estimates the baseline of a spectrum using the convex hull algorithm. This method +finds the lower convex hull of the spectrum, which is then used as the baseline. +It is generally faster than SNIP but can be less flexible. """ function convex_hull_baseline(y::AbstractVector{<:Real}) n = length(y) @@ -455,7 +576,7 @@ function median_baseline(y::AbstractVector{<:Real}; window::Int=20) end """ - apply_baseline_correction(y::AbstractVector{<:Real}; method::Symbol=:snip, iterations::Int=100, window::Int=20) -> Vector + apply_baseline_correction_core(y::AbstractVector{<:Real}; method::Symbol=:snip, iterations::Int=100, window::Int=20) -> Vector Applies a baseline correction algorithm to the intensity data. @@ -465,7 +586,7 @@ Applies a baseline correction algorithm to the intensity data. - `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) +function apply_baseline_correction_core(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 @@ -485,7 +606,9 @@ end """ tic_normalize(y) -> Vector -Normalizes spectrum intensities to the Total Ion Current (TIC). +Normalizes spectrum intensities to the Total Ion Current (TIC). Each intensity value +is divided by the sum of all intensities in the spectrum. This method assumes that the +total number of ions produced is similar for all samples. """ function tic_normalize(y::AbstractVector{<:Real}) s = sum(y) @@ -496,6 +619,18 @@ end pqn_normalize(M) -> Matrix Performs Probabilistic Quotient Normalization (PQN) on a matrix of spectra. +This is a more robust normalization method that is less sensitive to a small +number of highly abundant, variable peaks compared to TIC. + +# Steps: +1. A reference spectrum is calculated (typically the median spectrum across all samples). +2. For each spectrum, the quotients of its intensities and the reference spectrum's + intensities are calculated. +3. The median of these quotients is found for each spectrum. +4. Each spectrum is divided by its median quotient. + +# Arguments +- `M::AbstractMatrix{<:Real}`: A matrix where columns are spectra and rows are m/z bins. """ function pqn_normalize(M::AbstractMatrix{<:Real}) M_float = collect(float.(M)) @@ -532,7 +667,7 @@ function rms_normalize(y::AbstractVector{<:Real}) end """ - apply_normalization(y::AbstractVector{<:Real}; method::Symbol=:tic) -> Vector + apply_normalization_core(y::AbstractVector{<:Real}; method::Symbol=:tic) -> Vector Applies a per-spectrum normalization algorithm to the intensity data. @@ -540,7 +675,7 @@ Applies a per-spectrum normalization algorithm to the intensity data. - `y`: The intensity data. - `method`: The normalization method (:tic, :median, :rms, or :none). """ -function apply_normalization(y::AbstractVector{<:Real}; method::Symbol=:tic) +function apply_normalization_core(y::AbstractVector{<:Real}; method::Symbol=:tic)::Vector if method === :tic return tic_normalize(y) elseif method === :median @@ -555,188 +690,12 @@ function apply_normalization(y::AbstractVector{<:Real}; method::Symbol=:tic) 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; ...) -> Vector{NamedTuple} + detect_peaks_profile_core(mz, y; ...) -> Vector{NamedTuple} Enhanced peak detection for profile-mode spectra with advanced filtering and quality metrics. @@ -748,7 +707,7 @@ Returns a vector of NamedTuples, each representing a detected peak with: - `snr`: Signal-to-Noise Ratio - `prominence`: Peak prominence """ -function detect_peaks_profile(mz::AbstractVector{<:Real}, y::AbstractVector{<:Real}; +function detect_peaks_profile_core(mz::AbstractVector{<:Real}, y::AbstractVector{<:Real}; half_window::Int=10, snr_threshold::Float64=2.0, min_peak_prominence::Float64=0.1, @@ -762,7 +721,7 @@ function detect_peaks_profile(mz::AbstractVector{<:Real}, y::AbstractVector{<:Re 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; method=:savitzky_golay, window=max(5, 2*half_window+1), order=2) # Use smoothed data for detection + ys = smooth_spectrum_core(y; method=:savitzky_golay, window=max(5, 2*half_window+1), order=2) # Use smoothed data for detection candidate_peak_indices = Int[] for i in 2:n-1 @@ -818,16 +777,30 @@ function detect_peaks_profile(mz::AbstractVector{<:Real}, y::AbstractVector{<:Re end """ - detect_peaks_wavelet(mz, intensity; ...) -> Vector{NamedTuple} + detect_peaks_wavelet_core(mz, intensity; ...) -> Vector{NamedTuple} -Detects peaks using Continuous Wavelet Transform (CWT). +Detects peaks using Continuous Wavelet Transform (CWT). CWT is effective at +identifying peaks at different scales (widths), making it robust for complex spectra. -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) +# Arguments +- `mz::AbstractVector`: The m/z values of the spectrum. +- `intensity::AbstractVector`: The intensity values of the spectrum. +- `scales`: A range of scales to use for the CWT. Corresponds to the widths of + the features to be detected. +- `snr_threshold`: The minimum Signal-to-Noise Ratio for a CWT-detected local maximum + in the original spectrum to be considered a peak. +- `half_window`: Used for calculating peak quality metrics like FWHM and shape R^2. + +# Returns +A vector of `NamedTuple`s, each representing a detected peak with: +- `mz`: m/z value of the peak. +- `intensity`: Intensity of the peak from the original spectrum. +- `fwhm`: Full Width at Half Maximum (in ppm). +- `shape_r2`: Goodness-of-fit to a Gaussian shape. +- `snr`: Signal-to-Noise Ratio. +- `prominence`: Peak prominence (estimated as peak intensity for this method). """ -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}}} +function detect_peaks_wavelet_core(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}}[] @@ -876,7 +849,7 @@ function detect_peaks_wavelet(mz::AbstractVector, intensity::AbstractVector; sca end """ - detect_peaks_centroid(mz, y; ...) -> Vector{NamedTuple} + detect_peaks_centroid_core(mz, y; ...) -> Vector{NamedTuple} Filters peaks in centroid-mode data based on intensity threshold. @@ -884,7 +857,7 @@ 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}; snr_threshold::Float64=0.0) +function detect_peaks_centroid_core(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}}[] @@ -903,11 +876,11 @@ end # ============================================================================= """ - align_peaks_lowess(ref_mz, tgt_mz; ...) + align_peaks_lowess_core(ref_mz, tgt_mz; ...) Enhanced peak alignment with PPM tolerance and other constraints. """ -function align_peaks_lowess(ref_mz::Vector{<:Real}, tgt_mz::Vector{<:Real}; +function align_peaks_lowess_core(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, @@ -1030,15 +1003,15 @@ function align_peaks_lowess(ref_mz::Vector{<:Real}, tgt_mz::Vector{<:Real}; end """ - find_calibration_peaks(mz, intensity, reference_masses; ...) + find_calibration_peaks_core(mz, intensity, reference_masses; ...) Finds peaks that match a list of reference masses. """ -function find_calibration_peaks(mz::AbstractVector, intensity::AbstractVector, reference_masses::AbstractVector; ppm_tolerance=20.0) +function find_calibration_peaks_core(mz::AbstractVector, intensity::AbstractVector, reference_masses::AbstractVector; ppm_tolerance=20.0) matched_peaks = Dict{Float64, Float64}() - # detect_peaks_profile returns Vector{NamedTuple}, so we need to extract mz values - detected_peaks_list = detect_peaks_profile(mz, intensity) + # detect_peaks_profile_core returns Vector{NamedTuple}, so we need to extract mz values + detected_peaks_list = detect_peaks_profile_core(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] @@ -1057,16 +1030,16 @@ function find_calibration_peaks(mz::AbstractVector, intensity::AbstractVector, r end """ - calibrate_spectra(spectra, internal_standards; ...) + calibrate_spectra_core(spectra, internal_standards; ...) Calibrates spectra using internal standards. """ -function calibrate_spectra(spectra::Vector, internal_standards::Vector; ppm_tolerance=20.0) +function calibrate_spectra_core(spectra::Vector, internal_standards::Vector; ppm_tolerance=20.0) calibrated_spectra = similar(spectra) for (i, spec) in enumerate(spectra) mz, intensity = spec[1], spec[2] - matched_peaks = find_calibration_peaks(mz, intensity, internal_standards; ppm_tolerance=ppm_tolerance) + matched_peaks = find_calibration_peaks_core(mz, intensity, internal_standards; ppm_tolerance=ppm_tolerance) if length(matched_peaks) < 2 @warn "Spectrum $i: Not enough calibration peaks found. Skipping." calibrated_spectra[i] = spec @@ -1092,16 +1065,16 @@ end # ============================================================================= """ - bin_peaks(all_pk_mz::Vector{Vector{Float64}}, + bin_peaks_core(all_pk_mz::Vector{Vector{Float64}}, all_pk_int::Vector{Vector{Float64}}, - params::PeakBinningParams) -> Tuple{FeatureMatrix, Vector{Tuple{Float64,Float64}}} + params::PeakBinning) -> Tuple{FeatureMatrix, Vector{Tuple{Float64,Float64}}} 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. +- `params`: A `PeakBinning` struct. # Returns - `Tuple{FeatureMatrix, Vector{Tuple{Float64,Float64}}}`: A tuple containing the generated FeatureMatrix and the bin definitions. @@ -1111,8 +1084,8 @@ 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(spectra::Vector{MutableSpectrum}, - params::PeakBinningParams) +function bin_peaks_core(spectra::Vector{MutableSpectrum}, + params::PeakBinning) ns = length(spectra) # Number of spectra ns == 0 && return FeatureMatrix(zeros(0,0), Tuple{Float64,Float64}[], Int[]), Tuple{Float64,Float64}[] diff --git a/src/PreprocessingPipeline.jl b/src/PreprocessingPipeline.jl new file mode 100644 index 0000000..012f0fb --- /dev/null +++ b/src/PreprocessingPipeline.jl @@ -0,0 +1,455 @@ +# src/PreprocessingPipeline.jl + +using Base.Threads # For multithreading +using Printf # For @sprintf +using Interpolations # For linear_interpolation +using DataFrames # For saving feature matrix +using CSV # For saving feature matrix + +# This file provides a set of functions that apply preprocessing steps to a vector of +# `MutableSpectrum` objects. Each function takes the vector of spectra and a dictionary +# of parameters, modifying the spectra in-place where appropriate. This mirrors the +# logic from `test/run_preprocessing.jl` but is intended for use in the main application. + +# =================================================================== +# PREPROCESSING PIPELINE FUNCTIONS (IN-PLACE) +# =================================================================== + +""" + apply_baseline_correction(spectra::Vector{MutableSpectrum}, params::Dict) + +Applies baseline correction to the intensity data of each spectrum. This function +modifies the `.intensity` field of each `MutableSpectrum` object in-place. + +# Parameters from `params` Dict: +- `:method` (Symbol): The algorithm to use. Supports `:snip`, `:convex_hull`, `:median`. Defaults to `:snip`. +- `:iterations` (Int): The number of iterations for the SNIP algorithm. Defaults to 100. +- `:window` (Int): The window size for the Median algorithm. Defaults to 20. +""" +function apply_baseline_correction(spectra::Vector{MutableSpectrum}, params::Dict) + method = get(params, :method, :snip) + iterations = get(params, :iterations, 100) + window = get(params, :window, 20) + + Threads.@threads for s in spectra + if validate_spectrum(s.mz, s.intensity) + baseline = apply_baseline_correction_core(s.intensity; method=method, iterations=iterations, window=window) + s.intensity = max.(0.0, s.intensity .- baseline) + end + end +end + +""" + apply_intensity_transformation(spectra::Vector{MutableSpectrum}, params::Dict) + +Applies an intensity transformation to the intensity data of each spectrum. This function +modifies the `.intensity` field of each `MutableSpectrum` object in-place. + +# Parameters from `params` Dict: +- `:method` (Symbol): The transformation to apply. Supports `:sqrt`, `:log`, `:log2`, `:log10`, `:log1p`. Defaults to `:sqrt`. +""" +function apply_intensity_transformation(spectra::Vector{MutableSpectrum}, params::Dict) + method = get(params, :method, :sqrt) + + Threads.@threads for s in spectra + if validate_spectrum(s.mz, s.intensity) + s.intensity = transform_intensity_core(s.intensity; method=method) + end + end +end + +""" + apply_smoothing(spectra::Vector{MutableSpectrum}, params::Dict) + +Applies a smoothing filter to the intensity data of each spectrum. This function +modifies the `.intensity` field of each `MutableSpectrum` object in-place. + +# Parameters from `params` Dict: +- `:method` (Symbol): The smoothing algorithm. Supports `:savitzky_golay`, `:moving_average`. Defaults to `:savitzky_golay`. +- `:window` (Int): The size of the smoothing window. Defaults to 9. +- `:order` (Int): The polynomial order for the Savitzky-Golay filter. Defaults to 2. +""" +function apply_smoothing(spectra::Vector{MutableSpectrum}, params::Dict) + method = get(params, :method, :savitzky_golay) + window = get(params, :window, 9) + order = get(params, :order, 2) + + Threads.@threads for s in spectra + if validate_spectrum(s.mz, s.intensity) + smoothed_intensity = max.(0.0, smooth_spectrum_core(s.intensity; method=method, window=window, order=order)) + s.intensity = smoothed_intensity + end + end +end + +""" + apply_peak_picking(spectra::Vector{MutableSpectrum}, params::Dict) + +Detects peaks in each spectrum and stores them in the `.peaks` field of each +`MutableSpectrum` object, modifying it in-place. + +# Parameters from `params` Dict: +- `:method` (Symbol): The peak detection algorithm. Supports `:profile`, `:wavelet`, `:centroid`. Defaults to `:profile`. +- `:snr_threshold` (Float64): Signal-to-Noise Ratio threshold. +- `:half_window` (Int): Half-window size for local maxima detection. +- `:min_peak_prominence` (Float64): Minimum required prominence for a peak. +- `:merge_peaks_tolerance` (Float64): m/z tolerance to merge adjacent peaks. +""" +function apply_peak_picking(spectra::Vector{MutableSpectrum}, params::Dict) + 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) + + Threads.@threads for s in spectra + if validate_spectrum(s.mz, s.intensity) + if method == :profile + s.peaks = detect_peaks_profile_core(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_core(s.mz, s.intensity; snr_threshold=snr_threshold, half_window=half_window) + elseif method == :centroid + s.peaks = detect_peaks_centroid_core(s.mz, s.intensity; snr_threshold=snr_threshold) + else + s.peaks = detect_peaks_profile_core(s.mz, s.intensity; snr_threshold=snr_threshold, half_window=half_window) + end + else + s.peaks = [] + end + end +end + +""" + apply_peak_selection(spectra::Vector{MutableSpectrum}, params::Dict) + +Filters peaks within each spectrum based on quality criteria. This step removes peaks +that do not meet the specified thresholds for signal-to-noise ratio (SNR), +full width at half maximum (FWHM), and peak shape. + +This function modifies the `.peaks` field of each `MutableSpectrum` object in the `spectra` vector in-place. + +# Parameters from `params` Dict: +- `:min_snr` (Float64): The minimum Signal-to-Noise Ratio required for a peak to be kept. +- `:min_fwhm_ppm` (Float64): The minimum FWHM (in ppm) for a peak. +- `:max_fwhm_ppm` (Float64): The maximum FWHM (in ppm) for a peak. +- `:min_shape_r2` (Float64): The minimum R² value from a Gaussian fit, measuring peak shape quality. +""" +function apply_peak_selection(spectra::Vector{MutableSpectrum}, params::Dict) + min_snr = get(params, :min_snr, 0.0) + min_fwhm = get(params, :min_fwhm_ppm, 0.0) + max_fwhm = get(params, :max_fwhm_ppm, Inf) + min_r2 = get(params, :min_shape_r2, 0.0) + + # Handle `nothing` from params, which can happen if pre-calculation fails. + min_snr = isnothing(min_snr) ? 0.0 : min_snr + min_fwhm = isnothing(min_fwhm) ? 0.0 : min_fwhm + max_fwhm = isnothing(max_fwhm) ? Inf : max_fwhm + min_r2 = isnothing(min_r2) ? 0.0 : min_r2 + + Threads.@threads for s in spectra + if !isempty(s.peaks) + filter!(p -> + p.snr >= min_snr && + (min_fwhm <= p.fwhm <= max_fwhm) && + p.shape_r2 >= min_r2, + s.peaks + ) + end + end +end + +""" + apply_calibration(spectra::Vector{MutableSpectrum}, params::Dict, reference_peaks::Dict) + +Performs mass calibration on each spectrum using a list of internal standards. +This function modifies the `.mz` axis of each `MutableSpectrum` object in-place. + +# Parameters from `params` Dict: +- `:method` (Symbol): The calibration method. Only `:internal_standards` is currently meaningful. +- `:ppm_tolerance` (Float64): The tolerance in PPM for matching detected peaks to reference masses. +- `:fit_order` (Int): The polynomial order for the calibration fit (not yet used in this implementation, defaults to linear). +""" +function apply_calibration(spectra::Vector{MutableSpectrum}, params::Dict, reference_peaks::Dict) + method = get(params, :method, :none) + ppm_tolerance = get(params, :ppm_tolerance, 20.0) + + if method == :none || isempty(reference_peaks) + return + end + + reference_masses = collect(keys(reference_peaks)) + + Threads.@threads for i in 1:length(spectra) + s = spectra[i] + if validate_spectrum(s.mz, s.intensity) + matched_peaks = find_calibration_peaks_core(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 + else + @warn "Spectrum $(s.id): insufficient reference peaks ($(length(matched_peaks)) found), skipping calibration." + end + end + end +end + +""" + apply_peak_alignment(spectra::Vector{MutableSpectrum}, params::Dict) + +Aligns the m/z axis of all spectra to a chosen reference spectrum. This function +modifies both the `.mz` axis and the m/z values within the `.peaks` field of each +`MutableSpectrum` object in-place. + +# Parameters from `params` Dict: +- `:method` (Symbol): The alignment algorithm. Supports `:lowess`, `:linear`, `:ransac`. +- `:tolerance` (Float64): The tolerance for matching peaks between spectra. +- `:tolerance_unit` (Symbol): The unit for tolerance, `:mz` or `:ppm`. +""" +function apply_peak_alignment(spectra::Vector{MutableSpectrum}, params::Dict) + method = get(params, :method, :none) + tolerance = get(params, :tolerance, 0.002) + tolerance_unit = get(params, :tolerance_unit, :mz) + + if method == :none + return + end + + ref_find_idx = findfirst(s -> !isempty(s.peaks), spectra) + if ref_find_idx === nothing + @warn "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] + + 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_core(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 + +""" + apply_normalization(spectra::Vector{MutableSpectrum}, params::Dict) + +Applies intensity normalization to each spectrum. This function modifies the +`.intensity` field of each `MutableSpectrum` object in-place. + +# Parameters from `params` Dict: +- `:method` (Symbol): The normalization method. Supports `:tic`, `:median`, `:rms`, `:none`. +""" +function apply_normalization(spectra::Vector{MutableSpectrum}, params::Dict) + method = get(params, :method, :tic) + + Threads.@threads for s in spectra + if validate_spectrum(s.mz, s.intensity) + s.intensity = apply_normalization_core(s.intensity; method=method) + end + end +end +function apply_peak_binning(spectra::Vector{MutableSpectrum}, params::Dict) + tolerance = get(params, :tolerance, 20.0) + tolerance_unit = get(params, :tolerance_unit, :ppm) + min_peak_per_bin = get(params, :min_peak_per_bin, 3) + + if isempty(spectra) || all(s -> isempty(s.peaks), spectra) + @warn "No peaks found for binning. Returning empty feature matrix." + return nothing, nothing + end + + all_peaks = Vector{Tuple{Float64, Float64}}() + for s in spectra + for p in s.peaks + push!(all_peaks, (p.mz, p.intensity)) + end + end + + if isempty(all_peaks) + @warn "No peaks collected for binning." + return nothing, nothing + end + + sort!(all_peaks, by=x->x[1]) + + bin_centers = Float64[] + bin_intensities = Float64[] + + i = 1 + while i <= length(all_peaks) + current_bin_start = i + current_peak = all_peaks[i] + + j = i + 1 + while j <= length(all_peaks) + next_peak = all_peaks[j] + tol = (tolerance_unit == :ppm) ? (current_peak[1] * tolerance / 1e6) : tolerance + + if (next_peak[1] - current_peak[1]) <= tol + j += 1 + else + break + end + end + + current_bin_end = j - 1 + bin_size = current_bin_end - current_bin_start + 1 + + if bin_size >= min_peak_per_bin + bin_peaks = all_peaks[current_bin_start:current_bin_end] + + mz_sum = sum(p[1] for p in bin_peaks) + intensity_sum = sum(p[2] for p in bin_peaks) + + mz_center = mz_sum / bin_size + avg_intensity = intensity_sum / bin_size + + push!(bin_centers, mz_center) + push!(bin_intensities, avg_intensity) + end + + i = j + end + + if !isempty(bin_centers) + n_bins = length(bin_centers) + feature_matrix = Matrix{Float64}(undef, 2, n_bins) + + for i in 1:n_bins + feature_matrix[1, i] = bin_centers[i] + feature_matrix[2, i] = bin_intensities[i] + end + + bin_info = [(bin_centers[i], bin_intensities[i]) for i in 1:n_bins] + return feature_matrix, bin_info + else + @warn "No bins created after filtering" + return nothing, nothing + end +end + +""" + save_feature_matrix(feature_matrix::Matrix{Float64}, bin_info, output_dir::String) -> Tuple{String, String} + +Saves the aggregated `2 x n_bins` feature matrix into two different CSV formats. + +1. **Simple Format (`feature_matrix_simple.csv`):** A two-column CSV with "mz" and "intensity". +2. **Standard Format (`feature_matrix_standard.csv`):** A row-based format where m/z values are headers and there is a single data row for the aggregated spectrum. + +# Arguments +- `feature_matrix::Matrix{Float64}`: The `2 x n_bins` matrix from `apply_peak_binning`. +- `bin_info`: The associated bin information (currently unused but kept for compatibility). +- `output_dir::String`: The directory where the output CSV files will be saved. + +# Returns +- A tuple containing the paths to the two saved files. +""" +function save_feature_matrix(feature_matrix::Matrix{Float64}, bin_info, output_dir::String) + # Save as simple CSV with m/z and intensity rows + csv_path = joinpath(output_dir, "feature_matrix_simple.csv") + + open(csv_path, "w") do io + write(io, "mz,intensity\n") + for i in 1:size(feature_matrix, 2) + mz = feature_matrix[1, i] + intensity = feature_matrix[2, i] + write(io, "$mz,$intensity\n") + end + end + @info "Saved simple feature matrix: $csv_path" + + # Also save in a more standard format for MSI + csv_path_standard = joinpath(output_dir, "feature_matrix_standard.csv") + + open(csv_path_standard, "w") do io + write(io, "sample_type,") + mz_headers = [@sprintf("mz_%.4f", feature_matrix[1, i]) for i in 1:size(feature_matrix, 2)] + write(io, join(mz_headers, ",") * "\n") + + write(io, "aggregated_spectrum,") + intensity_values = [feature_matrix[2, i] for i in 1:size(feature_matrix, 2)] + write(io, join(string.(intensity_values), ",") * "\n") + end + @info "Saved standard format matrix: $csv_path_standard" + + return csv_path, csv_path_standard +end + +function execute_full_preprocessing(spectra::Vector{MutableSpectrum}, params::Dict, + pipeline_steps::Vector{String}, reference_peaks::Dict, + mask_path::Union{String, Nothing}=nothing; + progress_callback::Function=(step -> nothing)) + + println("Starting preprocessing pipeline with $(length(spectra)) spectra") + println("Steps: $(join(pipeline_steps, " -> "))") + + # These variables will be populated by the pipeline steps + feature_matrix = nothing + bin_definitions = nothing + + # Apply pipeline steps, modifying `spectra` in-place + for step in pipeline_steps + progress_callback(step) + println("\n" * "-"^60) + println("PROCESSING STEP: $step") + println("-"^60) + + if step == "stabilization" + println(" Applying intensity transformation (stabilization)") + apply_intensity_transformation(spectra, get(params, :Stabilization, Dict())) + + elseif step == "baseline_correction" + println(" Applying baseline correction") + apply_baseline_correction(spectra, get(params, :BaselineCorrection, Dict())) + + elseif step == "smoothing" + println(" Applying smoothing") + apply_smoothing(spectra, get(params, :Smoothing, Dict())) + + elseif step == "peak_picking" + println(" Applying peak picking") + apply_peak_picking(spectra, get(params, :PeakPicking, Dict())) + + elseif step == "peak_selection" + println(" Applying peak selection") + apply_peak_selection(spectra, get(params, :PeakSelection, Dict())) + + elseif step == "calibration" + println(" Applying calibration") + apply_calibration(spectra, get(params, :Calibration, Dict()), reference_peaks) + + elseif step == "peak_alignment" + println(" Applying peak alignment") + apply_peak_alignment(spectra, get(params, :PeakAlignment, Dict())) + + elseif step == "normalization" + println(" Applying normalization") + apply_normalization(spectra, get(params, :Normalization, Dict())) + + elseif step == "peak_binning" + println(" Applying peak binning") + feature_matrix, bin_definitions = apply_peak_binning(spectra, get(params, :PeakBinning, Dict())) + + else + @warn "Unknown step: $step, skipping" + end + + println("✓ Completed step: $step") + end + + return feature_matrix, bin_definitions +end \ No newline at end of file diff --git a/src/mzML.jl b/src/mzML.jl index f795f9b..bcf1f17 100644 --- a/src/mzML.jl +++ b/src/mzML.jl @@ -20,6 +20,7 @@ const DATA_FORMAT_ACCESSIONS = Dict{String, DataType}( ) function parse_instrument_metadata_mzml(stream::IO) + println("DEBUG: Starting mzML instrument metadata parsing...") # Initialize with default values from the InstrumentMetadata constructor instrument_meta = InstrumentMetadata() @@ -67,37 +68,53 @@ function parse_instrument_metadata_mzml(stream::IO) if acc == "MS:1000031" # instrument model instrument_model = val + #println("DEBUG: Instrument model: $instrument_model") elseif acc == "MS:1001496" # mass resolving power (more specific) resolution = tryparse(Float64, val) + #println("DEBUG: Resolution (specific): $resolution") elseif acc == "MS:1000011" && resolution === nothing # resolution (less specific) resolution = tryparse(Float64, val) + #println("DEBUG: Resolution (less specific): $resolution") elseif acc == "MS:1000016" # mass accuracy (ppm) mass_accuracy_ppm = tryparse(Float64, val) + #println("DEBUG: Mass accuracy (ppm): $mass_accuracy_ppm") elseif acc == "MS:1000130" # positive scan polarity = :positive + #println("DEBUG: Polarity: positive") elseif acc == "MS:1000129" # negative scan polarity = :negative + #println("DEBUG: Polarity: negative") elseif acc == "MS:1000592" # external calibration calibration_status = :external + #println("DEBUG: Calibration status: external") elseif acc == "MS:1000593" # internal calibration calibration_status = :internal + #println("DEBUG: Calibration status: internal") elseif acc == "MS:1000747" && calibration_status == :uncalibrated # instrument specific calibration calibration_status = :internal # Assume as a form of internal calibration + #println("DEBUG: Calibration status: instrument specific (internal)") elseif acc == "MS:1000867" # laser wavelength laser_settings["wavelength_nm"] = tryparse(Float64, val) + #println("DEBUG: Laser wavelength: $(laser_settings["wavelength_nm"]) nm") elseif acc == "MS:1000868" # laser fluence laser_settings["fluence"] = tryparse(Float64, val) + #println("DEBUG: Laser fluence: $(laser_settings["fluence"])") elseif acc == "MS:1000869" # laser repetition rate laser_settings["repetition_rate_hz"] = tryparse(Float64, val) + #println("DEBUG: Laser repetition rate: $(laser_settings["repetition_rate_hz"]) Hz") # NEW: Vendor Preprocessing terms elseif acc == "MS:1000579" # baseline correction push!(vendor_preprocessing_steps, "Baseline Correction") + #println("DEBUG: Vendor preprocessing step: Baseline Correction") elseif acc == "MS:1000580" # smoothing push!(vendor_preprocessing_steps, "Smoothing") + #println("DEBUG: Vendor preprocessing step: Smoothing") elseif acc == "MS:1000578" # data transformation (e.g., centroiding) push!(vendor_preprocessing_steps, "Data Transformation: $(name)") + #println("DEBUG: Vendor preprocessing step: Data Transformation: $(name)") elseif acc == "MS:1000800" # deisotoping push!(vendor_preprocessing_steps, "Deisotoping") + #println("DEBUG: Vendor preprocessing step: Deisotoping") end end end @@ -106,6 +123,8 @@ function parse_instrument_metadata_mzml(stream::IO) @warn "Could not fully parse instrument metadata from mzML header. Using defaults. Error: $e" end + println("DEBUG: Finished mzML instrument metadata parsing.") + # Always return a valid object return InstrumentMetadata( resolution, @@ -137,14 +156,15 @@ determine the data type, compression, and axis type. function get_spectrum_asset_metadata(stream::IO) start_pos = position(stream) - bda_tag = find_tag(stream, r"") binary_offset = position(stream) + #println("DEBUG: Binary data offset: $binary_offset") # Move stream to the end of the binary data array for the next iteration readuntil(stream, "") + #println("DEBUG: Exiting get_spectrum_asset_metadata.") # Create SpectrumAsset directly from the variables return SpectrumAsset(data_format, compression_flag, binary_offset, encoded_length, axis) @@ -223,13 +250,16 @@ function parse_spectrum_metadata(stream::IO, offset::Int64) id_match = match(r"= 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 + return b 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 + baseline_iter_20 = simple_snip(raw_signal, 20) + baseline_iter_200 = simple_snip(raw_signal, 200) + corrected_signal = raw_signal .- baseline_iter_200 + + lines!(ax, mz_range, raw_signal, color=(:grey, 0.6), label="Raw Signal") + lines!(ax, mz_range, corrected_signal, color=:green, linewidth=2.5, label="Corrected Signal") + l1 = lines!(ax, mz_range, baseline_iter_20, color=(:red, 0.7), linestyle=:dash, linewidth=2, label="Baseline (iterations: 20)") + l2 = lines!(ax, mz_range, baseline_iter_200, color=:red, linewidth=2.5, label="Baseline (iterations: 200)") + + text!(ax, "Parameter: `iterations`\nMore iterations result in a more aggressive baseline that follows the signal floor more closely.", + position=Point2f(0.05, 0.95), space=:relative, align=(:left, :top), fontsize=12) + axislegend(ax, position=:rt) + save("descriptive_plot_baseline_correction.png", fig) + println("Saved: descriptive_plot_baseline_correction.png") end -# Create and display the plot -fig = create_msi_parameter_plot() +""" +Plots the effect of smoothing, comparing different window sizes. +""" +function plot_smoothing_details() + fig = Figure(size=(1200, 700), fontsize=14) + ax = Axis(fig[1, 1], title="Detailed View: Smoothing (Savitzky-Golay)", xlabel="m/z", ylabel="Intensity") -# Save the plot -save("msi_preprocessing_parameters.png", fig) -println("Plot saved as 'msi_preprocessing_parameters.png'") + mz_range = range(400, 700, length=1000) + true_signal = 0.8 .* exp.(-0.002 .* (mz_range .- 500).^2) .+ 0.6 .* exp.(-0.001 .* (mz_range .- 600).^2) + noisy_signal = true_signal .+ 0.1 .* randn(length(mz_range)) + + function simple_moving_average(y, window) + smoothed = similar(y) + for i in 1:length(y) + start_idx = max(1, i - window ÷ 2) + end_idx = min(length(y), i + window ÷ 2) + smoothed[i] = mean(y[start_idx:end_idx]) + end + return smoothed + end + + smoothed_small_window = simple_moving_average(noisy_signal, 5) + smoothed_large_window = simple_moving_average(noisy_signal, 21) + + lines!(ax, mz_range, noisy_signal, color=(:red, 0.4), label="Noisy Signal") + lines!(ax, mz_range, true_signal, color=:black, linestyle=:dash, linewidth=2, label="True Signal") + lines!(ax, mz_range, smoothed_small_window, color=:blue, linewidth=2, label="Smoothed (window: 5)") + lines!(ax, mz_range, smoothed_large_window, color=:purple, linewidth=2, label="Smoothed (window: 21)") + + text!(ax, "Parameter: `window`\nA larger window increases smoothing but may broaden peaks.", + position=Point2f(0.05, 0.95), space=:relative, align=(:left, :top), fontsize=12) + axislegend(ax, position=:rt) + save("descriptive_plot_smoothing.png", fig) + println("Saved: descriptive_plot_smoothing.png") +end + +""" +Visualizes the peak picking process for profile-mode data, illustrating the +effects of SNR threshold and peak prominence. +""" +function plot_peak_picking_profile_details() + fig = Figure(size=(1200, 700), fontsize=14) + ax = Axis(fig[1, 1], title="Detailed View: Peak Picking (Profile Mode)", xlabel="m/z", ylabel="Intensity") + + mz = 1:200 + base_signal = 10 .* exp.(-((mz .- 50).^2) ./ (2*3^2)) .+ 7 .* exp.(-((mz .- 120).^2) ./ (2*5^2)) .+ 2 + small_peak_signal = zeros(200) + for i in 80:90 + small_peak_signal[i] = 3 * exp(-((i - 85)^2) / 2.0) + end + noise = 0.5 .* randn(200) + intensity = base_signal .+ small_peak_signal .+ noise + + noise_level = median(abs.(intensity .- median(intensity))) * 1.4826 # MAD + snr_threshold_val = 3.0 + intensity_threshold = noise_level * snr_threshold_val + + picked_peaks_mz = [50, 85, 120] + picked_peaks_intensity = intensity[picked_peaks_mz] + + hlines!(ax, [noise_level], color=:gray, linestyle=:dot, label="Est. Noise Level") + hlines!(ax, [intensity_threshold], color=:orange, linestyle=:dash, label="SNR Threshold (snr_threshold = 3.0)") + lines!(ax, mz, intensity, color=:blue, label="Profile Spectrum") + scatter!(ax, picked_peaks_mz, picked_peaks_intensity, color=:green, markersize=15, strokewidth=2, label="Peaks passing SNR") + scatter!(ax, [25], [intensity[25]], color=:red, marker=:x, markersize=15, label="Local max below SNR") + + # Illustrate prominence + arrows!(ax, [120, 120], [intensity[135], intensity[120]], [0, 0], [intensity[120]-intensity[135], 0], color=:purple) + text!(ax, 125, (intensity[120]+intensity[135])/2, text="Prominence", color=:purple) + + axislegend(ax, position=:rt) + save("descriptive_plot_peakpicking_profile.png", fig) + println("Saved: descriptive_plot_peakpicking_profile.png") +end + +""" +Visualizes peak picking (filtering) for centroid-mode data based on an SNR +(intensity) threshold. +""" +function plot_peak_picking_centroid_details() + fig = Figure(size=(1200, 700), fontsize=14) + ax = Axis(fig[1, 1], title="Detailed View: Peak Picking (Centroid Mode)", xlabel="m/z", ylabel="Intensity") + + mz = [100, 150, 200, 250, 300, 350, 400] + intensity = [10, 5, 25, 8, 3, 18, 12] + snr_threshold_val = 10.0 + + stem!(ax, mz, intensity, color=:gray, label="Input Centroids") + + selected_mask = intensity .>= snr_threshold_val + stem!(ax, mz[selected_mask], intensity[selected_mask], color=:green, trunkwidth=3, label="Selected Peaks (intensity >= 10)") + stem!(ax, mz[.!selected_mask], intensity[.!selected_mask], color=:red, trunkwidth=3, label="Rejected Peaks (intensity < 10)") + + hlines!(ax, [snr_threshold_val], color=:orange, linestyle=:dash, label="Intensity Threshold (snr_threshold)") + + text!(ax, "Parameter: `snr_threshold`\nIn centroid mode, this acts as a direct intensity filter.", + position=Point2f(0.05, 0.95), space=:relative, align=(:left, :top), fontsize=12) + axislegend(ax, position=:rt) + save("descriptive_plot_peakpicking_centroid.png", fig) + println("Saved: descriptive_plot_peakpicking_centroid.png") +end + +""" +Plots the effect of different normalization methods on a set of spectra. +""" +function plot_normalization_details() + fig = Figure(size=(1200, 700), fontsize=14) + + mz = range(300, 400, length=500) + spec1 = 1.5 .* exp.(-((mz .- 350).^2) ./ 50) + spec2 = 0.8 .* exp.(-((mz .- 350).^2) ./ 50) + + ax1 = Axis(fig[1, 1], title="Before Normalization", ylabel="Absolute Intensity") + lines!(ax1, mz, spec1, label="Spectrum A (High TIC)") + lines!(ax1, mz, spec2, label="Spectrum B (Low TIC)") + axislegend(ax1) + + ax2 = Axis(fig[1, 2], title="After Normalization (TIC)", ylabel="Relative Intensity") + lines!(ax2, mz, spec1 ./ sum(spec1), label="Spectrum A (Normalized)") + lines!(ax2, mz, spec2 ./ sum(spec2), label="Spectrum B (Normalized)") + + text!(ax2, "Effect: Spectra are scaled to have the same total area, making their intensities comparable.", + position=Point2f(0.05, 0.95), space=:relative, align=(:left, :top), fontsize=12, justification=:left) + + save("descriptive_plot_normalization.png", fig) + println("Saved: descriptive_plot_normalization.png") +end + +""" +Illustrates mass calibration, showing how a calibration curve corrects measured +m/z values based on reference peaks. +""" +function plot_calibration_details() + fig = Figure(size=(1200, 700), fontsize=14) + ax = Axis(fig[1, 1], title="Detailed View: Calibration", xlabel="Measured m/z", ylabel="m/z Error (Measured - Reference)") + + ref_mz = [200, 400, 600, 800, 1000] + measured_mz = ref_mz .+ [0.1, 0.15, 0.2, 0.25, 0.3] .+ 0.02 .* randn(5) + errors = measured_mz .- ref_mz + + # Fit a linear model to the error + A = [ones(5) measured_mz] + coeffs = A \ errors + correction_func(m) = m - (coeffs[1] .+ coeffs[2] .* m) + + fit_line = coeffs[1] .+ coeffs[2] .* measured_mz + + scatter!(ax, measured_mz, errors, color=:red, markersize=15, label="Measured Error") + lines!(ax, measured_mz, fit_line, color=:blue, label="Calibration Curve (fit_order=1)") + + text!(ax, "Parameter: `fit_order`\nA curve is fit to the error of known reference peaks.\nThis curve is then used to correct all m/z values.", + position=Point2f(0.05, 0.95), space=:relative, align=(:left, :top), fontsize=12) + axislegend(ax, position=:rb) + save("descriptive_plot_calibration.png", fig) + println("Saved: descriptive_plot_calibration.png") +end + +""" +Visualizes peak alignment by showing multiple spectra with misaligned peaks +before and after the alignment process. +""" +function plot_alignment_details() + fig = Figure(size=(1600, 600), fontsize=14) + ax1 = Axis(fig[1, 1], title="Before Alignment", xlabel="m/z", yticklabelsvisible=false, ygridvisible=false) + ax2 = Axis(fig[1, 2], title="After Alignment", xlabel="m/z", yticklabelsvisible=false, ygridvisible=false) + + mz = range(490, 510, length=1000) + shifts = [-0.5, 0.0, 0.8] + colors = [:blue, :green, :purple] + + for (i, shift) in enumerate(shifts) + peak_center = 500 + shift + spectrum = exp.(-((mz .- peak_center).^2) ./ 0.1) + lines!(ax1, mz, spectrum .+ i, color=colors[i]) + vlines!(ax1, [peak_center], color=(colors[i], 0.5), linestyle=:dash) + + # After alignment, all peaks are at 500 + aligned_spectrum = exp.(-((mz .- 500).^2) ./ 0.1) + lines!(ax2, mz, aligned_spectrum .+ i, color=colors[i]) + end + vlines!(ax2, [500], color=:red, linestyle=:dash, label="Reference m/z") + axislegend(ax2) + + save("descriptive_plot_alignment.png", fig) + println("Saved: descriptive_plot_alignment.png") +end + +""" +Illustrates peak selection by filtering a population of peaks based on +FWHM and SNR criteria. +""" +function plot_peak_selection_details() + fig = Figure(size=(1200, 700), fontsize=14) + ax = Axis(fig[1, 1], title="Detailed View: Peak Selection", xlabel="FWHM (ppm)", ylabel="Signal-to-Noise Ratio (SNR)") + + n_peaks = 100 + fwhm = rand(n_peaks) .* 150 + snr = rand(n_peaks) .* 20 + + min_fwhm_ppm = 20.0 + max_fwhm_ppm = 100.0 + min_snr = 5.0 + + selected_mask = (fwhm .>= min_fwhm_ppm) .& (fwhm .<= max_fwhm_ppm) .& (snr .>= min_snr) + + scatter!(ax, fwhm[.!selected_mask], snr[.!selected_mask], color=(:red, 0.5), label="Rejected Peaks") + scatter!(ax, fwhm[selected_mask], snr[selected_mask], color=:green, label="Selected Peaks") + + vlines!(ax, [min_fwhm_ppm, max_fwhm_ppm], color=:blue, linestyle=:dash, label="FWHM bounds") + hlines!(ax, [min_snr], color=:orange, linestyle=:dash, label="SNR bound") + + poly!(ax, BBox(min_fwhm_ppm, max_fwhm_ppm, min_snr, 22), color=(:green, 0.1)) + text!(ax, "Selection Region", position=(60, 12), color=:green, fontsize=14) + + axislegend(ax) + save("descriptive_plot_peak_selection.png", fig) + println("Saved: descriptive_plot_peak_selection.png") +end + +""" +Visualizes the adaptive peak binning process, showing how peaks from different +spectra are grouped into a common bin based on a PPM tolerance. +""" +function plot_peak_binning_details() + fig = Figure(size=(1200, 700), fontsize=14) + ax = Axis(fig[1, 1], title="Detailed View: Adaptive Peak Binning", xlabel="m/z", yticklabelsvisible=false) + + ref_mz = 500.0 + tolerance_ppm = 50.0 + tol_mz = ref_mz * tolerance_ppm / 1e6 + + bin_start = ref_mz - tol_mz/2 + bin_end = ref_mz + tol_mz/2 + + peaks_mz = [ref_mz - 0.01, ref_mz + 0.005, ref_mz + 0.02, ref_mz - 0.015] + peak_intensities = [0.8, 1.0, 0.9, 0.7] + peak_colors = [:blue, :green, :purple, :orange] + + vspan!(ax, bin_start, bin_end, color=(:gray, 0.2), label="Bin (tolerance: 50 ppm)") + stem!(ax, peaks_mz, peak_intensities, color=peak_colors, markersize=15) + + # Show bin center + bin_center = mean(peaks_mz) + vlines!(ax, [bin_center], color=:red, linestyle=:dash, label="Calculated Bin Center") + + text!(ax, "Parameter: `tolerance`\nPeaks from different spectra within the tolerance window are grouped into a single feature.", + position=Point2f(0.05, 0.95), space=:relative, align=(:left, :top), fontsize=12) + axislegend(ax, position=:rt) + xlims!(ax, ref_mz - tol_mz*2, ref_mz + tol_mz*2) + save("descriptive_plot_peak_binning.png", fig) + println("Saved: descriptive_plot_peak_binning.png") +end + + +""" +Main function to generate and save all descriptive plots. +""" +function create_and_save_all_plots() + println("Generating detailed descriptive plots for preprocessing steps...") + + plot_baseline_correction_details() + plot_smoothing_details() + plot_peak_picking_profile_details() + plot_peak_picking_centroid_details() + plot_normalization_details() + plot_calibration_details() + plot_alignment_details() + plot_peak_selection_details() + plot_peak_binning_details() + + println("\nAll descriptive plots have been saved in the current directory.") +end + +# Execute the plot generation +if abspath(PROGRAM_FILE) == @__FILE__ + create_and_save_all_plots() +end -# Display the plot (if in an interactive environment) -fig diff --git a/test/run_preprocessing.jl b/test/run_preprocessing.jl index 2e0e329..e244f8a 100644 --- a/test/run_preprocessing.jl +++ b/test/run_preprocessing.jl @@ -5,6 +5,8 @@ import Pkg using CairoMakie using DataFrames # For creating dataframes using CSV +using Statistics +using Interpolations # --- Load the MSI_src Module --- Pkg.activate(joinpath(@__DIR__, "..")) @@ -15,12 +17,13 @@ using MSI_src # CONFIG: PLEASE FILL IN YOUR FILE PATHS HERE # =================================================================== -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/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" +const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML" +# const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML" -# const MASK_ROUTE = "/home/pixel/Documents/Cinvestav_2025/JuliaMSI/public/css/masks/Stomach_DHB_uncompressed.png" -const MASK_ROUTE = "" +const MASK_ROUTE = "/home/pixel/Documents/Cinvestav_2025/JuliaMSI/public/css/masks/Stomach_DHB_uncompressed.png" +# const MASK_ROUTE = "" const OUTPUT_DIR = "./test/results/preprocessing_results" @@ -41,6 +44,55 @@ const reference_peaks = Dict( 124.0393 => "Tropine [M+H]+", ) +const PIPELINE_STP = [ + "stabilization", + #"baseline_correction", + "smoothing", + "peak_picking", + "peak_selection", + #"calibration", + "peak_alignment", + "normalization", + "peak_binning" +] + +# =================================================================== +# 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( + :Stabilization => Dict( + :method => :sqrt # Default stabilization method + ), + :PeakPicking => Dict( + #:snr_threshold => 1.5, + #:half_window => 5, + :snr_threshold => 8.0, + #:merge_peaks_tolerance => 2.5, + #:half_window => 2 + #:half_window => 3 + ), + # :Smoothing => Dict( + # :window => 11 + # ) + #:PeakAlignment => Dict( + # :tolerance => 0.01 + #) + #:PeakSelection => Dict( + #:frequency_threshold => 0, + #:min_shape_r2 => 0.5, + #:max_fwhm_ppm => 30.0, + #:min_fwhm_ppm => 2.0 + #), + #:PeakBinning => Dict( + #:tolerance => 30.0, + #:frequency_threshold => 0 + #) +) + # =================================================================== # HELPER FUNCTIONS # =================================================================== @@ -79,7 +131,7 @@ end # PREPROCESSING PIPELINE FUNCTIONS (IN-PLACE) # =================================================================== -function apply_baseline_correction(spectra::Vector{MutableSpectrum}, params::Dict, msi_data::MSIData) +function apply_baseline_correction_core(spectra::Vector{MutableSpectrum}, params::Dict, msi_data::MSIData) print_step_header("Baseline Correction") method = get(params, :method, :snip) @@ -91,7 +143,7 @@ function apply_baseline_correction(spectra::Vector{MutableSpectrum}, params::Dic 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) + baseline = MSI_src.apply_baseline_correction_core(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 @@ -99,7 +151,7 @@ function apply_baseline_correction(spectra::Vector{MutableSpectrum}, params::Dic 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) + baseline = MSI_src.apply_baseline_correction_core(s.intensity; method=method, iterations=iterations, window=window) s.intensity = max.(0.0, s.intensity .- baseline) end end @@ -117,14 +169,14 @@ function apply_smoothing(spectra::Vector{MutableSpectrum}, params::Dict) 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)) + smoothed_intensity = max.(0.0, smooth_spectrum_core(s.intensity; method=method, window=window, order=order)) plot_spectrum_step(s.mz, smoothed_intensity, "smoothing", spectrum_index=s.id) end end 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)) + smoothed_intensity = max.(0.0, smooth_spectrum_core(s.intensity; method=method, window=window, order=order)) s.intensity = smoothed_intensity end end @@ -153,13 +205,13 @@ function apply_peak_picking(spectra::Vector{MutableSpectrum}, params::Dict) 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) + s.peaks = detect_peaks_profile_core(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) + s.peaks = detect_peaks_centroid_core(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) + s.peaks = detect_peaks_profile_core(s.mz, s.intensity; snr_threshold=snr_threshold, half_window=half_window) end else s.peaks = [] @@ -176,6 +228,51 @@ function apply_peak_picking(spectra::Vector{MutableSpectrum}, params::Dict) end end +function apply_peak_selection(spectra::Vector{MutableSpectrum}, params::Dict) + print_step_header("Peak Selection") + + min_snr = get(params, :min_snr, 0.0) + min_fwhm = get(params, :min_fwhm_ppm, 0.0) + max_fwhm = get(params, :max_fwhm_ppm, Inf) + min_r2 = get(params, :min_shape_r2, 0.0) + + # Handle `nothing` values from params, default to non-filtering values + min_snr = isnothing(min_snr) ? 0.0 : min_snr + min_fwhm = isnothing(min_fwhm) ? 0.0 : min_fwhm + max_fwhm = isnothing(max_fwhm) ? Inf : max_fwhm + min_r2 = isnothing(min_r2) ? 0.0 : min_r2 + + println(" - Min SNR: $min_snr") + println(" - FWHM Range (ppm): [$min_fwhm, $max_fwhm]") + println(" - Min Shape R²: $min_r2") + + total_peaks_before = sum(s -> length(s.peaks), spectra) + + Threads.@threads for s in spectra + if !isempty(s.peaks) + s.peaks = filter(p -> + p.snr >= min_snr && + (min_fwhm <= p.fwhm <= max_fwhm) && + p.shape_r2 >= min_r2, + s.peaks + ) + end + end + + total_peaks_after = sum(s -> length(s.peaks), spectra) + println(" - Peaks before: $total_peaks_before, Peaks after: $total_peaks_after") + + # Safely plot the first spectrum to show effect of filtering + if !isempty(spectra) + s1_idx = findfirst(s -> s.id == 1, spectra) + if s1_idx !== nothing + s1 = spectra[s1_idx] + plot_spectrum_step(s1.mz, s1.intensity, "peak_selection", peaks=s1.peaks, spectrum_index=s1.id) + println(" - Spectrum 1 now has $(length(s1.peaks)) peaks after selection.") + end + end +end + function apply_calibration(spectra::Vector{MutableSpectrum}, params::Dict, reference_peaks::Dict) print_step_header("Calibration") @@ -195,7 +292,7 @@ function apply_calibration(spectra::Vector{MutableSpectrum}, params::Dict, refer 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) + matched_peaks = find_calibration_peaks_core(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))) @@ -248,7 +345,7 @@ function apply_peak_alignment(spectra::Vector{MutableSpectrum}, params::Dict) 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) + alignment_func = align_peaks_lowess_core(ref_peaks_mz, current_peaks_mz; method=method, tolerance=tolerance, tolerance_unit=tolerance_unit) s.mz = alignment_func.(s.mz) # Update m/z axis @@ -261,7 +358,7 @@ function apply_peak_alignment(spectra::Vector{MutableSpectrum}, params::Dict) end end -function apply_normalization(spectra::Vector{MutableSpectrum}, params::Dict) +function apply_normalization_core(spectra::Vector{MutableSpectrum}, params::Dict) print_step_header("Normalization") method = get(params, :method, :tic) @@ -273,7 +370,7 @@ function apply_normalization(spectra::Vector{MutableSpectrum}, params::Dict) if s1 !== nothing s = spectra[s1] if validate_spectrum(s.mz, s.intensity) - normalized_intensity = MSI_src.apply_normalization(s.intensity; method=method) + normalized_intensity = MSI_src.apply_normalization_core(s.intensity; method=method) plot_spectrum_step(s.mz, normalized_intensity, "normalization") end end @@ -281,7 +378,7 @@ function apply_normalization(spectra::Vector{MutableSpectrum}, params::Dict) Threads.@threads for s in spectra if validate_spectrum(s.mz, s.intensity) - s.intensity = MSI_src.apply_normalization(s.intensity; method=method) + s.intensity = MSI_src.apply_normalization_core(s.intensity; method=method) end end end @@ -293,8 +390,6 @@ function apply_peak_binning(spectra::Vector{MutableSpectrum}, params::Dict) 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) @@ -304,92 +399,129 @@ function apply_peak_binning(spectra::Vector{MutableSpectrum}, params::Dict) 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") + # Collect all peaks with their intensities + all_peaks = Vector{Tuple{Float64, Float64}}() # (mz, intensity) + for s in spectra + for p in s.peaks + push!(all_peaks, (p.mz, p.intensity)) 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") + if isempty(all_peaks) + @warn "No peaks collected for binning." + return nothing, nothing + end - return csv_path, bins_path + sort!(all_peaks, by=x->x[1]) + + # Create bins - just store mz_center and intensity + bin_centers = Float64[] + bin_intensities = Float64[] + + i = 1 + while i <= length(all_peaks) + current_bin_start = i + current_peak = all_peaks[i] + + # Find all peaks in this bin + j = i + 1 + while j <= length(all_peaks) + next_peak = all_peaks[j] + # Calculate tolerance + tol = (tolerance_unit == :ppm) ? (current_peak[1] * tolerance / 1e6) : tolerance + + if (next_peak[1] - current_peak[1]) <= tol + j += 1 + else + break + end + end + + current_bin_end = j - 1 + bin_size = current_bin_end - current_bin_start + 1 + + # Check if we have enough peaks in this bin + if bin_size >= min_peak_per_bin + bin_peaks_core = all_peaks[current_bin_start:current_bin_end] + + # Calculate m/z center and average intensity + mz_sum = 0.0 + intensity_sum = 0.0 + for peak in bin_peaks_core + mz_sum += peak[1] + intensity_sum += peak[2] + end + + mz_center = mz_sum / bin_size + avg_intensity = intensity_sum / bin_size + + push!(bin_centers, mz_center) + push!(bin_intensities, avg_intensity) + end + + i = j # Move to next potential bin + end + + # Create the 2-row matrix + if !isempty(bin_centers) + n_bins = length(bin_centers) + feature_matrix = Matrix{Float64}(undef, 2, n_bins) + + for i in 1:n_bins + feature_matrix[1, i] = bin_centers[i] + feature_matrix[2, i] = bin_intensities[i] + end + + println(" - Created feature matrix: 2 × $n_bins") + println(" - Number of bins created: $n_bins") + println(" - m/z range: $(round(bin_centers[1], digits=4)) - $(round(bin_centers[end], digits=4))") + + # Return bin info as a vector of tuples for compatibility + bin_info = [(bin_centers[i], bin_intensities[i]) for i in 1:n_bins] + return feature_matrix, bin_info + else + @warn "No bins created after filtering" + return nothing, nothing + end 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 - ) -) +function save_feature_matrix(feature_matrix::Matrix{Float64}, bin_info) + print_step_header("Saving Results") + + # Save as simple CSV with m/z and intensity rows + csv_path = joinpath(OUTPUT_DIR, "feature_matrix_simple.csv") + + open(csv_path, "w") do io + # Write header + write(io, "mz,intensity\n") + + # Write data: m/z values in first column, intensities in second + for i in 1:size(feature_matrix, 2) + mz = feature_matrix[1, i] + intensity = feature_matrix[2, i] + write(io, "$mz,$intensity\n") + end + end + println(" - Saved simple feature matrix: $csv_path") + + # Also save in a more standard format for MSI + csv_path_standard = joinpath(OUTPUT_DIR, "feature_matrix_standard.csv") + + open(csv_path_standard, "w") do io + # Write header with m/z values as column names + write(io, "sample_type,") + mz_headers = [@sprintf("mz_%.4f", feature_matrix[1, i]) for i in 1:size(feature_matrix, 2)] + write(io, join(mz_headers, ",") * "\n") + + # Write the aggregated intensity values + write(io, "aggregated_spectrum,") + intensity_values = [feature_matrix[2, i] for i in 1:size(feature_matrix, 2)] + write(io, join(string.(intensity_values), ",") * "\n") + end + println(" - Saved standard format matrix: $csv_path_standard") + + return csv_path, csv_path_standard +end # =================================================================== # MAIN PREPROCESSING PIPELINE @@ -446,28 +578,45 @@ function run_preprocessing_pipeline() end # Define pipeline steps - pipeline_steps = [ - "baseline_correction", - "smoothing", - "peak_picking", - "calibration", - "peak_alignment", - "normalization", - "peak_binning" - ] + pipeline_steps = PIPELINE_STP 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, []) + local spectrum_indices_to_process::AbstractVector{Int} + if !isempty(MASK_ROUTE) + println("Applying mask from: $(MASK_ROUTE)") + try + mask_matrix = MSI_src.load_and_prepare_mask(MASK_ROUTE, msi_data.image_dims) + masked_indices_set = MSI_src.get_masked_spectrum_indices(msi_data, mask_matrix) + spectrum_indices_to_process = collect(masked_indices_set) + println("Mask applied. $(length(spectrum_indices_to_process)) spectra are within the masked region.") + catch e + @error "Failed to load or apply mask: $e. Proceeding without mask." + spectrum_indices_to_process = 1:length(msi_data.spectra_metadata) + end + else + spectrum_indices_to_process = 1:length(msi_data.spectra_metadata) + end + + if isempty(spectrum_indices_to_process) + @warn "No spectra available for processing after applying mask/filter. Exiting pipeline." + close(msi_data) + return + end + + num_spectra_to_process = length(spectrum_indices_to_process) + current_spectra = Vector{MutableSpectrum}(undef, num_spectra_to_process) + + Threads.@threads for i in 1:num_spectra_to_process + original_idx = spectrum_indices_to_process[i] + mz, intensity = MSI_src.GetSpectrum(msi_data, original_idx) + current_spectra[i] = MSI_src.MutableSpectrum(original_idx, mz, intensity, []) - # Plot raw spectrum without masks - if idx == 1 - plot_spectrum_step(mz, intensity, "raw_unmasked_spectrum") + # Plot raw spectrum without masks (only the first processed one) + if i == 1 + plot_spectrum_step(mz, intensity, "raw_unmasked_spectrum", spectrum_index=original_idx) end end @@ -481,14 +630,42 @@ function run_preprocessing_pipeline() println("PROCESSING STEP: $step") println("-"^60) - if step == "baseline_correction" - @time apply_baseline_correction(current_spectra, auto_params[:BaselineCorrection], msi_data) + if step == "stabilization" + print_step_header("Intensity Transformation (Stabilization)") + method = get(auto_params[:Stabilization], :method, :sqrt) + println(" Method: $method") + + # Safely plot the first spectrum before transformation + if !isempty(current_spectra) + s = current_spectra[1] + if validate_spectrum(s.mz, s.intensity) + # Use a temporary spectrum to plot before and after + initial_intensity = deepcopy(s.intensity) + plot_spectrum_step(s.mz, initial_intensity, "stabilization_before", spectrum_index=s.id) + end + end + + @time apply_intensity_transformation(current_spectra, auto_params[:Stabilization]) + + # Safely plot the first spectrum after transformation + if !isempty(current_spectra) + s = current_spectra[1] + if validate_spectrum(s.mz, s.intensity) + plot_spectrum_step(s.mz, s.intensity, "stabilization_after", spectrum_index=s.id) + end + end + + elseif step == "baseline_correction" + @time apply_baseline_correction_core(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 == "peak_selection" + @time apply_peak_selection(current_spectra, auto_params[:PeakSelection]) elseif step == "calibration" @time apply_calibration(current_spectra, auto_params[:Calibration], reference_peaks) @@ -497,14 +674,14 @@ function run_preprocessing_pipeline() @time apply_peak_alignment(current_spectra, auto_params[:PeakAlignment]) elseif step == "normalization" - @time apply_normalization(current_spectra, auto_params[:Normalization]) + @time apply_normalization_core(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]) + feature_matrix, bin_info = @time apply_peak_binning(current_spectra, auto_params[:PeakBinning]) if feature_matrix !== nothing - @time save_feature_matrix(feature_matrix, bin_definitions) + @time save_feature_matrix(feature_matrix, bin_info) end else