+ label="Select Right Dataset" class="q-ma-sm" style="min-width: 200px;"
+ v-on:focus="refetch_folders = true">
@@ -561,7 +765,7 @@
Dataset Summary
+ style="min-width: 250px;" standout="custom-standout" v-on:focus="refetch_folders = true">
diff --git a/julia_imzML_visual.jl b/julia_imzML_visual.jl
index b9b40c4..cedb8d9 100644
--- a/julia_imzML_visual.jl
+++ b/julia_imzML_visual.jl
@@ -666,10 +666,10 @@ function xySpectrumPlot(data::MSIData, xCoord::Int, yCoord::Int, imgWidth::Int,
# Downsample for plotting performance
mz_down, int_down = MSI_src.downsample_spectrum(mz, intensity)
- trace = if spectrum_mode == MSI_src.PROFILE
- PlotlyBase.scatter(x=mz_down, y=int_down, mode="lines", marker=attr(size=1, color="blue", opacity=0.5), name="Spectrum", hoverinfo="x", hovertemplate="m/z: %{x:.4f}")
- else
+ trace = if spectrum_mode == MSI_src.CENTROID
PlotlyBase.stem(x=mz_down, y=int_down, marker=attr(size=1, color="blue", opacity=0.5), name="Spectrum", hoverinfo="x", hovertemplate="m/z: %{x:.4f}")
+ else
+ PlotlyBase.scatter(x=mz_down, y=int_down, mode="lines", marker=attr(size=1, color="blue", opacity=0.5), name="Spectrum", hoverinfo="x", hovertemplate="m/z: %{x:.4f}")
end
plotdata = [trace]
diff --git a/mask.jl b/mask.jl
index c426bf4..45365a3 100644
--- a/mask.jl
+++ b/mask.jl
@@ -151,6 +151,7 @@ end
@out available_folders = String[]
@out image_available_folders = String[]
@private registry_init_done = false
+ @in refetch_folders = false
@out imgInt = "" # Path to current slice
@out current_msi = ""
@@ -711,6 +712,23 @@ end
end
end
+ @onbutton refetch_folders begin
+ # First, re-run the logic to populate the folder list
+ 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)
+
+ # Apply the logic to select the first item if none is selected
+ if !isempty(image_available_folders) && isempty(selected_folder_main)
+ # This ensures the list is up-to-date
+ image_available_folders = deepcopy(img_folders)
+ selected_folder_main = first(image_available_folders)
+ end
+ end
+
@onchange isready begin
if isready && !registry_init_done
sleep(1.0) # Give frontend time to initialize
diff --git a/mask.jl.html b/mask.jl.html
index 60d7011..2ae26f3 100644
--- a/mask.jl.html
+++ b/mask.jl.html
@@ -15,7 +15,7 @@
Step 1: Select Slice
+ class="q-ma-sm col" standout="custom-standout" :disable="is_editing_mask" v-on:focus="refetch_folders = true">
diff --git a/public/css/genieapp.css b/public/css/genieapp.css
index 98ddaf2..9aa8eed 100644
--- a/public/css/genieapp.css
+++ b/public/css/genieapp.css
@@ -45,9 +45,10 @@
color: #009f90;
}
-#intDivStyle-left {
+#intDivStyle-left, #intDivStyle-right {
border-radius: 10px;
padding: 10px;
+ background-color: #fdfdfd; /* Add a background color to distinguish from the outer div */
}
/* Tab styling to match your theme */
@@ -81,13 +82,13 @@
font-size: 0.9rem;
}
-#intDivStyle-left .q-tab-panels {
- height: 700px; /* Set this to accommodate your tallest content */
+#intDivStyle-left .q-tab-panels, #intDivStyle-right .q-tab-panels {
+ /* Removed fixed height */
}
-#intDivStyle-left .q-tab-panel {
- height: 100%;
- overflow-y: auto; /* Add scroll if content overflows */
+#intDivStyle-left .q-tab-panel, #intDivStyle-right .q-tab-panel {
+ height: auto; /* Let content define height */
+ overflow-y: hidden; /* Remove scrollbar */
}
.pixelated-plot svg image {
@@ -98,4 +99,22 @@
*{
font-family: 'Roboto', 'Lato', sans-serif;
-}
\ No newline at end of file
+}
+
+.q-card {
+ border-radius: 8px !important;
+ box-shadow: 0 2px 4px rgba(0,0,0,0.1) !important;
+ margin-bottom: 16px !important;
+}
+
+.q-card__section--dark {
+ background: #f7f7f7 !important;
+}
+
+.q-radio__label {
+ font-size: 1rem !important;
+}
+
+.q-option-group > div {
+ margin-bottom: 8px;
+}
diff --git a/src/BloomFilters.jl b/src/BloomFilters.jl
index 2d95ed2..e5a6790 100644
--- a/src/BloomFilters.jl
+++ b/src/BloomFilters.jl
@@ -37,9 +37,20 @@ mutable struct StreamingBloomFilter
end
"""
- BloomFilter(expected_elements::Int, false_positive_rate::Float64=0.01; seed::UInt64=0x12345678)
+ BloomFilter{T}(expected_elements::Int, false_positive_rate::Float64=0.01; kwargs...) -> Return type
Creates a Bloom filter optimized for the expected number of elements and desired false positive rate.
+
+# Arguments
+
+- `expected_elements::Int`: Argument description
+- `false_positive_rate::Float64`: Argument description
+ (**Default**: `0.01`)
+
+# Keywords
+
+- `seed::Union{UInt32,UInt64}`: Keyword description
+ (**Default**: `0x12345678`)
"""
function BloomFilter{T}(expected_elements::Int, false_positive_rate::Float64=0.01; seed::Union{UInt32,UInt64}=0x12345678) where T
# Convert seed to UInt64 for consistency
@@ -56,9 +67,12 @@ end
"""
optimal_bit_size(n::Int, p::Float64) -> Int
-Calculates the optimal number of bits for a Bloom filter given:
-- n: expected number of elements
-- p: desired false positive rate
+Calculates the optimal number of bits for a Bloom filter
+
+# Arguments
+
+- `n::Int`: Expected number of elements
+- `p::Float64`: Desired false positive rate
"""
function optimal_bit_size(n::Int, p::Float64)::Int
if p <= 0.0 || p >= 1.0
@@ -73,6 +87,11 @@ end
optimal_hash_count(n::Int, m::Int) -> Int
Calculates the optimal number of hash functions for a Bloom filter.
+
+# Arguments
+
+- `n::Int`: Expected number of elements
+- `p::Float64`: Desired false positive rate
"""
function optimal_hash_count(n::Int, m::Int)::Int
if n <= 0 || m <= 0
@@ -109,6 +128,11 @@ end
Base.push!(bf::BloomFilter{T}, item::T)
Adds an element to the Bloom filter.
+
+# Arguments
+
+- `bf::BloomFilter{T}`: Argument description
+- `item::T`: Argument description
"""
function Base.push!(bf::BloomFilter{T}, item::T) where T
hashes = hash_functions(item, bf.hash_count, bf.size, bf.seed)
@@ -169,27 +193,39 @@ function false_positive_rate(bf::BloomFilter)::Float64
end
"""
- fill_ratio(bf::BloomFilter) -> Float64
+ fill_ratio(bf::BloomFilter)::Float64 return count(bf.bits) / length(bf.bits) end -> Return type
Returns the fraction of bits that are set to 1.
+
+# Arguments
+
+- `bf::BloomFilter`: Argument description
"""
function fill_ratio(bf::BloomFilter)::Float64
return count(bf.bits) / length(bf.bits)
end
"""
- is_empty(bf::BloomFilter) -> Bool
+ is_empty(bf::BloomFilter)::Bool return bf.count == 0 end -> Return type
Checks if the Bloom filter is empty (no elements added).
+
+# Arguments
+
+- `bf::BloomFilter`: Argument description
"""
function is_empty(bf::BloomFilter)::Bool
return bf.count == 0
end
"""
- reset!(bf::BloomFilter)
+ reset!(bf::BloomFilter) fill!(bf.bits, false) bf.count = 0 return bf end -> Return type
Clears the Bloom filter, removing all elements.
+
+# Arguments
+
+- `bf::BloomFilter`: Argument description
"""
function reset!(bf::BloomFilter)
fill!(bf.bits, false)
@@ -198,6 +234,20 @@ function reset!(bf::BloomFilter)
end
# Specialized constructor for empty Bloom filters
+"""
+ BloomFilter{T}(; kwargs...) -> Return type
+
+Description of the function
+
+# Keywords
+
+- `size::Int`: Keyword description
+ (**Default**: `100`)
+- `hash_count::Int`: Keyword description
+ (**Default**: `3`)
+- `seed::Union{UInt32,UInt64}`: Keyword description
+ (**Default**: `0x12345678`)
+"""
function BloomFilter{T}(;size::Int=100, hash_count::Int=3, seed::Union{UInt32,UInt64}=0x12345678) where T
seed_uint64 = UInt64(seed)
bits = falses(size)
diff --git a/src/MSIData.jl b/src/MSIData.jl
index ed47a8b..6b7c8d7 100644
--- a/src/MSIData.jl
+++ b/src/MSIData.jl
@@ -420,6 +420,12 @@ function read_binary_vector(data::MSIData, io::IO, asset::SpectrumAsset)
return out_array
end
+function read_binary_vector(data::MSIData, ts_handle::ThreadSafeFileHandle, asset::SpectrumAsset)
+ lock(ts_handle.lock) do
+ return read_binary_vector(data, ts_handle.handle, asset)
+ end
+end
+
# Overload for different source types
"""
read_spectrum_from_disk(source::ImzMLSource, meta::SpectrumMetadata)
diff --git a/src/MSI_src.jl b/src/MSI_src.jl
index bf5df59..d878bd2 100644
--- a/src/MSI_src.jl
+++ b/src/MSI_src.jl
@@ -28,8 +28,12 @@ export FeatureMatrix,
snip_baseline,
tic_normalize,
pqn_normalize,
+ median_normalize,
detect_peaks_profile,
+ detect_peaks_wavelet,
+ detect_peaks_centroid,
align_peaks_lowess,
+ find_calibration_peaks,
bin_peaks,
plot_stage_spectrum,
calculate_ppm_error,
diff --git a/src/Preprocessing.jl b/src/Preprocessing.jl
index 62e1e30..a1a8cff 100644
--- a/src/Preprocessing.jl
+++ b/src/Preprocessing.jl
@@ -16,6 +16,8 @@ using SavitzkyGolay # For SavitzkyGolay filtering
using Dates # For now()
using CSV # For writing CSV files
using DataFrames # For creating dataframes
+using ContinuousWavelets # For CWT peak detection
+using Interpolations # For calibration
# =============================================================================
# Data Structures
@@ -24,19 +26,29 @@ using DataFrames # For creating dataframes
"""
FeatureMatrix
-Structure to hold the final feature matrix, including m/z bin boundaries and sample indices.
-
-# Fields
-- `matrix`: The numerical matrix where rows are samples and columns are features (m/z bins).
-- `mz_bins`: A vector of tuples `(low_mz, high_mz)` for each feature column.
-- `sample_ids`: A vector of indices corresponding to the original spectra.
+A struct to hold the final feature matrix generated from the preprocessing pipeline.
"""
struct FeatureMatrix
- matrix::Array{Float64,2} # samples × features
- mz_bins::Vector{Tuple{Float64,Float64}} # [(low, high), ...]
+ matrix::Array{Float64,2}
+ mz_bins::Vector{Tuple{Float64,Float64}}
sample_ids::Vector{Int}
end
+"""
+ QCParameters
+
+A struct to hold comprehensive quality control parameters for validating spectra.
+"""
+struct QCParameters
+ min_tic_threshold::Float64 # Minimum total ion current
+ max_tic_threshold::Float64 # Maximum TIC (saturation check)
+ snr_threshold::Float64 # Minimum signal-to-noise for a spectrum to be considered valid
+ peak_count_range::Tuple{Int,Int} # Acceptable peak count range (min, max)
+ spatial_consistency_threshold::Float64 # For MSI spatial coherence (stub)
+ calibration_accuracy_ppm::Float64 # Maximum allowed PPM error for calibration lock masses
+end
+
+
# =============================================================================
# 0) Quality Control (QC)
# =============================================================================
@@ -52,7 +64,7 @@ qc_is_empty(mz::AbstractVector, intensity::AbstractVector)::Bool =
"""
qc_is_regular(mz) -> Bool
-Checks that the m/z axis is monotonically non-decreasing, as expected in a profile spectrum.
+Checks that the m/z axis is monotonically non-decreasing.
"""
function qc_is_regular(mz::AbstractVector)
n = length(mz)
@@ -73,26 +85,33 @@ end
transform_intensity(intensity; method=:sqrt) -> Vector
Applies a variance-stabilizing transformation to the intensity vector.
-Supported methods: `:sqrt` (default) and `:log1p`.
"""
function transform_intensity(intensity::AbstractVector{<:Real}; method::Symbol=:sqrt)
if method === :sqrt
return sqrt.(max.(zero(eltype(intensity)), intensity))
elseif method === :log1p
return log1p.(max.(zero(eltype(intensity)), intensity))
+ elseif method === :log
+ return log.(max.(eps(eltype(intensity)), intensity))
+ elseif method === :log2
+ return log2.(max.(eps(eltype(intensity)), intensity))
+ elseif method === :log10
+ return log10.(max.(eps(eltype(intensity)), intensity))
else
return collect(float.(intensity))
end
end
"""
- smooth_spectrum(y; window=21, order=2) -> Vector
+ smooth_spectrum(y; window=9, order=2) -> Vector
-Applies a Savitzky–Golay filter if `SavitzkyGolay.jl` is available.
-Otherwise, falls back to a simple (non-phase-correct) moving average.
+Applies a Savitzky–Golay filter to smooth the intensity data.
"""
function smooth_spectrum(y::AbstractVector{<:Real}; window::Int=9, order::Int=2)
win = isodd(window) ? window : window + 1
+ if length(y) < win
+ return y # Cannot smooth if data is smaller than window
+ end
res = SavitzkyGolay.savitzky_golay(collect(float.(y)), win, order)
return res.y
end
@@ -104,22 +123,20 @@ end
"""
snip_baseline(y, iterations=100) -> Vector
-Estimates the baseline using a simple 1D SNIP (Statistics-sensitive Non-linear
-Iterative Peak-clipping) algorithm. `iterations` controls the aggressiveness.
+Estimates the baseline of a spectrum using the SNIP algorithm.
"""
-function snip_baseline(y::AbstractVector{<:Real}, iterations::Int=100)
+function snip_baseline(y::AbstractVector{<:Real}; iterations::Int=100)
n = length(y)
- b = collect(float.(y)) # work copy
+ b = collect(float.(y))
buf = similar(b)
for k in 1:iterations
copyto!(buf, b)
@inbounds for i in 2:n-1
buf[i] = min(b[i], 0.5 * (b[i-1] + b[i+1]))
end
- # Handle endpoints
buf[1] = min(b[1], b[2])
buf[end] = min(b[end], b[end-1])
- b, buf = buf, b # Swap buffers
+ b, buf = buf, b
end
return b
end
@@ -131,7 +148,7 @@ end
"""
tic_normalize(y) -> Vector
-Normalizes intensities to the Total Ion Current (TIC). If the sum is zero, returns a copy.
+Normalizes spectrum intensities to the Total Ion Current (TIC).
"""
function tic_normalize(y::AbstractVector{<:Real})
s = sum(y)
@@ -141,108 +158,154 @@ end
"""
pqn_normalize(M) -> Matrix
-Performs Probabilistic Quotient Normalization on a matrix `M` where columns are spectra.
+Performs Probabilistic Quotient Normalization (PQN) on a matrix of spectra.
"""
function pqn_normalize(M::AbstractMatrix{<:Real})
M_float = collect(float.(M))
- # Calculate reference spectrum (median across all spectra)
ref = mapslices(median, M_float; dims=2)[:,1]
-
- # Calculate quotients for each spectrum relative to the reference
Q = similar(M_float)
@inbounds for j in axes(M_float, 2)
Q[:, j] = M_float[:, j] ./ (ref .+ eps(eltype(M_float)))
end
-
- # Find the median quotient for each spectrum (scaling factor)
- s = [median( @view Q[:, j]) for j in axes(Q, 2)]
-
- # Normalize the original matrix
+ s = [median(@view Q[:, j]) for j in axes(Q, 2)]
@inbounds for j in axes(M_float, 2)
M_float[:, j] ./= (s[j] + eps(eltype(M_float)))
end
return M_float
end
+"""
+ median_normalize(y) -> Vector
+
+Normalizes spectrum intensities by dividing by the median intensity.
+"""
+function median_normalize(y::AbstractVector{<:Real})
+ m = median(y)
+ return m <= 0 ? collect(float.(y)) : collect(float.(y)) ./ m
+end
+
# =============================================================================
-# 4) Peak Detection (for Profile Data)
+# 4) Peak Detection
# =============================================================================
"""
- detect_peaks_profile(mz, y; half_window=10, snr_threshold=2.0)
+ detect_peaks_profile(mz, y; ...)
-Detects local maxima with a signal-to-noise threshold (using MAD for noise estimation).
-Assumes profile-mode data and a monotonic m/z axis.
+Enhanced peak detection for profile-mode spectra with advanced filtering.
"""
-function detect_peaks_profile(mz::AbstractVector{<:Real},
- y::AbstractVector{<:Real};
+function detect_peaks_profile(mz::AbstractVector{<:Real}, y::AbstractVector{<:Real};
half_window::Int=10,
- snr_threshold::Float64=2.0)
+ snr_threshold::Float64=2.0,
+ min_peak_width::Real=0.01,
+ max_peak_width::Real=2.0,
+ peak_shape_threshold::Float64=0.7, # Stub
+ merge_peaks_tolerance::Float64=0.002,
+ min_peak_prominence::Float64=0.1)
n = length(y)
n < 3 && return (Float64[], Float64[])
- # Noise estimation using Median Absolute Deviation (robust to peaks)
noise_level = mad(y, normalize=true) + eps(Float64)
-
- # Smooth the spectrum to make peak detection more robust
ys = smooth_spectrum(y; window=max(5, 2*half_window+1), order=2)
- peak_idx = Int[]
+ peak_indices = Int[]
@inbounds for i in 2:n-1
left = max(1, i - half_window)
right = min(n, i + half_window)
- local_max = ys[i]
- # A point is a peak if it's the maximum in its neighborhood and above the SNR threshold
- if local_max >= maximum( @view ys[left:right]) && (local_max > snr_threshold * noise_level)
- # Ensure we only record one point for flat-topped peaks
- if isempty(peak_idx) || (i - last(peak_idx) > half_window)
- push!(peak_idx, i)
- end
+ # Prominence check
+ prominence = ys[i] - max(minimum(@view ys[left:i]), minimum(@view ys[i:right]))
+
+ if ys[i] >= maximum(@view ys[left:right]) &&
+ (ys[i] > snr_threshold * noise_level) &&
+ (prominence > min_peak_prominence * ys[i])
+ push!(peak_indices, i)
end
end
+
+ # (STUB) Peak shape validation would be applied here
+ # e.g., fitting a Gaussian and checking R^2 > peak_shape_threshold
- # Return original intensities at peak locations
- pk_mz = [float(mz[i]) for i in peak_idx]
- pk_int = [float(y[i]) for i in peak_idx]
+ # Merge close peaks
+ if !isempty(peak_indices) && merge_peaks_tolerance > 0
+ merged_indices = [peak_indices[1]]
+ for i in 2:length(peak_indices)
+ if (mz[peak_indices[i]] - mz[last(merged_indices)]) > merge_peaks_tolerance
+ push!(merged_indices, peak_indices[i])
+ elseif y[peak_indices[i]] > y[last(merged_indices)]
+ merged_indices[end] = peak_indices[i] # Replace with more intense peak
+ end
+ end
+ peak_indices = merged_indices
+ end
+
+ # (STUB) Peak width filtering would be applied here
+ # This would require calculating FWHM for each peak, which is computationally intensive.
+
+ pk_mz = [float(mz[i]) for i in peak_indices]
+ pk_int = [float(y[i]) for i in peak_indices]
return (pk_mz, pk_int)
end
"""
- detect_peaks_centroid(mz, y; intensity_threshold=0.0)
+ detect_peaks_wavelet(mz, intensity; ...)
-Filters centroided data based on a minimum intensity threshold.
+Detects peaks using Continuous Wavelet Transform (CWT).
"""
-function detect_peaks_centroid(mz::AbstractVector{<:Real},
- y::AbstractVector{<:Real};
- intensity_threshold::Float64=0.0)
+function detect_peaks_wavelet(mz::AbstractVector, intensity::AbstractVector; scales=1:10, snr_threshold=3.0)
+ n = length(intensity)
+ n < 10 && return (Float64[], Float64[])
+ cwt_res = ContinuousWavelets.cwt(intensity, ContinuousWavelets.morl)
+ peak_indices = Int[]
+ noise_level = mad(intensity, normalize=true) + eps(Float64)
+
+ for i in 2:n-1
+ if intensity[i] > intensity[i-1] && intensity[i] > intensity[i+1] && intensity[i] > snr_threshold * noise_level
+ if abs(cwt_res[i, end]) > 0
+ push!(peak_indices, i)
+ end
+ end
+ end
+
+ return (mz[peak_indices], intensity[peak_indices])
+end
+
+"""
+ detect_peaks_centroid(mz, y; ...)
+
+Filters peaks in centroid-mode data.
+"""
+function detect_peaks_centroid(mz::AbstractVector{<:Real}, y::AbstractVector{<:Real}; intensity_threshold::Float64=0.0)
keep_indices = findall(y .>= intensity_threshold)
-
return (mz[keep_indices], y[keep_indices])
end
# =============================================================================
-# 5) Peak Alignment
+# 5) Peak Alignment & Calibration
# =============================================================================
"""
- align_peaks_lowess(ref_mz, tgt_mz; tolerance=0.002) -> warp::Function
+ align_peaks_lowess(ref_mz, tgt_mz; ...)
-Generates a warping function `warp(x)` to map target m/z values to reference m/z values.
-Uses a lightweight LOWESS-like approach with linear interpolation.
+Enhanced peak alignment with PPM tolerance and other constraints.
"""
-function align_peaks_lowess(ref_mz::Vector{<:Real},
- tgt_mz::Vector{<:Real};
- tolerance::Float64=0.002)
- # Efficiently match peaks between sorted lists
- pairs = Tuple{Float64,Float64}[] # (target_mz, reference_mz)
+function align_peaks_lowess(ref_mz::Vector{<:Real}, tgt_mz::Vector{<:Real};
+ tolerance::Float64=0.002,
+ tolerance_unit::Symbol=:mz,
+ max_shift_ppm::Float64=50.0,
+ min_matched_peaks::Int=5)
+ pairs = Tuple{Float64,Float64}[]
i = 1; j = 1
while i <= length(tgt_mz) && j <= length(ref_mz)
+ tol = (tolerance_unit == :ppm) ? (ref_mz[j] * tolerance / 1e6) : tolerance
dt = tgt_mz[i] - ref_mz[j]
- if abs(dt) <= tolerance
- push!(pairs, (float(tgt_mz[i]), float(ref_mz[j])))
+
+ if abs(dt) <= tol
+ # Max shift check
+ if abs(dt) * 1e6 / ref_mz[j] <= max_shift_ppm
+ push!(pairs, (float(tgt_mz[i]), float(ref_mz[j])))
+ end
i += 1; j += 1
elseif dt < 0
i += 1
@@ -251,42 +314,68 @@ function align_peaks_lowess(ref_mz::Vector{<:Real},
end
end
- if length(pairs) < 3
- @warn "Too few matching peaks for alignment. Returning identity function."
+ if length(pairs) < min_matched_peaks
+ @warn "Too few matching peaks ($(length(pairs)) < $min_matched_peaks). Returning identity function."
return x -> float.(x)
end
t = [p[1] for p in pairs]
r = [p[2] for p in pairs]
+
+ # (STUB) A robust regression model (e.g., RANSAC) would be better here.
+ itp = linear_interpolation(t, r, extrapolation_bc=Line())
+ return itp
+end
- # Lightly smooth the mapping to reduce noise
- t_s = smooth_spectrum(t; window=5, order=2)
- r_s = smooth_spectrum(r; window=5, order=2)
+"""
+ find_calibration_peaks(mz, intensity, reference_masses; ...)
- # Return a function that performs linear interpolation for warping
- function warp(x::AbstractVector{<:Real})
- out = similar(collect(float.(x)))
- for (k, xv) in enumerate(x)
- if xv <= t_s[1]
- # Linear extrapolation at the start
- m = (r_s[2]-r_s[1]) / (t_s[2]-t_s[1] + eps())
- out[k] = r_s[1] + m*(xv - t_s[1])
- elseif xv >= t_s[end]
- # Linear extrapolation at the end
- m = (r_s[end]-r_s[end-1]) / (t_s[end]-t_s[end-1] + eps())
- out[k] = r_s[end-1] + m*(xv - t_s[end-1])
- else
- # Linear interpolation for points in the middle
- lo = searchsortedlast(t_s, xv)
- hi = lo + 1
- α = (xv - t_s[lo]) / (t_s[hi] - t_s[lo] + eps())
- out[k] = (1-α)*r_s[lo] + α*r_s[hi]
- end
+Finds peaks that match a list of reference masses.
+"""
+function find_calibration_peaks(mz::AbstractVector, intensity::AbstractVector, reference_masses::AbstractVector; ppm_tolerance=20.0)
+ matched_peaks = Dict{Float64, Float64}()
+ detected_mz, _ = detect_peaks_profile(mz, intensity)
+
+ for ref_mass in reference_masses
+ tol = ref_mass * ppm_tolerance / 1e6
+ candidates = findall(m -> abs(m - ref_mass) <= tol, detected_mz)
+ if !isempty(candidates)
+ closest_peak_idx = argmin(abs.(detected_mz[candidates] .- ref_mass))
+ matched_peaks[ref_mass] = detected_mz[candidates[closest_peak_idx]]
end
- return out
end
+ return matched_peaks
+end
- return warp
+"""
+ calibrate_spectra(spectra, internal_standards; ...)
+
+Calibrates spectra using internal standards.
+"""
+function calibrate_spectra(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)
+ if length(matched_peaks) < 2
+ @warn "Spectrum $i: Not enough calibration peaks found. Skipping."
+ calibrated_spectra[i] = spec
+ continue
+ end
+ measured = sort(collect(values(matched_peaks)))
+ theoretical = sort(collect(keys(matched_peaks)))
+ itp = linear_interpolation(measured, theoretical, extrapolation_bc=Line())
+
+ new_mz = itp(mz)
+
+ if length(spec) == 3
+ calibrated_spectra[i] = (new_mz, intensity, spec[3])
+ else
+ calibrated_spectra[i] = (new_mz, intensity)
+ end
+ end
+ return calibrated_spectra
end
# =============================================================================
@@ -294,92 +383,199 @@ end
# =============================================================================
"""
- _find_bin_index(x, bins) -> Int
+ bin_peaks(all_pk_mz, all_pk_int, tolerance; ...)
-Efficiently finds the index of the bin `(low, high)` that contains `x` using binary search.
-Returns 0 if not found.
-"""
-function _find_bin_index(x::Float64, bins::Vector{Tuple{Float64,Float64}})
- lo, hi = 1, length(bins)
- while lo <= hi
- mid = (lo + hi) >>> 1
- b = bins[mid]
- if x < b[1]
- hi = mid - 1
- elseif x > b[2]
- lo = mid + 1
- else
- return mid
- end
- end
- return 0
-end
-
-"""
- bin_peaks(all_pk_mz, all_pk_int, tolerance; frequency_threshold=0.25)
-
-Groups peaks from all spectra into consensus m/z bins and creates a feature matrix.
-Filters out features that do not appear in a minimum fraction of spectra.
+Enhanced peak binning with adaptive and PPM-based parameters.
"""
function bin_peaks(all_pk_mz::Vector{<:AbstractVector{<:Real}},
all_pk_int::Vector{<:AbstractVector{<:Real}},
- tolerance::Float64; frequency_threshold::Float64=0.25)
+ tolerance::Float64;
+ frequency_threshold::Float64=0.25,
+ tolerance_unit::Symbol=:mz,
+ adaptive_tolerance::Bool=false, # Stub
+ min_peak_per_bin::Int=2,
+ max_bin_width_ppm::Float64=100.0,
+ intensity_weighted_centers::Bool=true)
ns = length(all_pk_mz)
ns == 0 && return (zeros(0,0), Tuple{Float64,Float64}[])
- # 1) Collect all unique peak m/z values and sort them
- flat_mz = Float64[]
- for v in all_pk_mz
- append!(flat_mz, float.(v))
- end
- sort!(flat_mz)
- isempty(flat_mz) && return (zeros(ns, 0), Tuple{Float64,Float64}[])
+ # Collect all peaks with their intensities and original spectrum index
+ all_peaks = [(float(m), float(i), s_idx) for s_idx in 1:ns for (m, i) in zip(all_pk_mz[s_idx], all_pk_int[s_idx])]
+ sort!(all_peaks, by=p->p[1])
+
+ isempty(all_peaks) && return (zeros(ns, 0), Tuple{Float64,Float64}[])
- # 2) Create contiguous m/z bins based on tolerance
- bins = Tuple{Float64,Float64}[]
- cur_lo = flat_mz[1]
- cur_hi = flat_mz[1]
- for x in @view flat_mz[2:end]
- if x - cur_hi <= tolerance
- cur_hi = x # Extend the current bin
+ # Create bins
+ bins = []
+ current_bin_peaks = [all_peaks[1]]
+ for p in @view all_peaks[2:end]
+ bin_center = mean(first.(current_bin_peaks))
+ tol = (tolerance_unit == :ppm) ? (bin_center * tolerance / 1e6) : tolerance
+
+ if p[1] - last(current_bin_peaks)[1] <= tol
+ push!(current_bin_peaks, p)
else
- push!(bins, (cur_lo, cur_hi)) # Finalize old bin
- cur_lo = x; cur_hi = x # Start a new one
- end
- end
- push!(bins, (cur_lo, cur_hi))
-
- # 3) Create the feature matrix (samples x features) using max intensity per bin
- X = zeros(Float64, ns, length(bins))
- for i in 1:ns
- for (mzv, iv) in zip(all_pk_mz[i], all_pk_int[i])
- bidx = _find_bin_index(float(mzv), bins)
- if bidx > 0
- X[i, bidx] = max(X[i, bidx], float(iv))
+ if length(current_bin_peaks) >= min_peak_per_bin
+ push!(bins, current_bin_peaks)
end
+ current_bin_peaks = [p]
+ end
+ end
+ length(current_bin_peaks) >= min_peak_per_bin && push!(bins, current_bin_peaks)
+
+ # Filter bins by max width
+ filter!(b -> (last(b)[1] - first(b)[1]) * 1e6 / mean(p[1] for p in b) <= max_bin_width_ppm, bins)
+
+ # Create feature matrix
+ X = zeros(Float64, ns, length(bins))
+ final_bins_boundaries = Vector{Tuple{Float64,Float64}}(undef, length(bins))
+
+ for (j, bin_peaks) in enumerate(bins)
+ # Calculate bin center
+ local bin_center
+ if intensity_weighted_centers
+ weights = [p[2] for p in bin_peaks]
+ bin_center = sum(p[1]*p[2] for p in bin_peaks) / sum(weights)
+ else
+ bin_center = mean(p[1] for p in bin_peaks)
+ end
+
+ final_bins_boundaries[j] = (first(bin_peaks)[1], last(bin_peaks)[1])
+
+ for p in bin_peaks
+ s_idx = p[3]
+ X[s_idx, j] = max(X[s_idx, j], p[2])
end
end
- # 4) Filter features by minimum frequency
+ # Filter by frequency
if frequency_threshold > 0
present_count = vec(sum(X .> 0, dims=1))
min_count = ceil(Int, frequency_threshold * ns)
keep_mask = findall(present_count .>= min_count)
X = X[:, keep_mask]
- bins = bins[keep_mask]
+ final_bins_boundaries = final_bins_boundaries[keep_mask]
end
- return (X, bins)
+ return (X, final_bins_boundaries)
end
+# =============================================================================
+# 7) Spatial & Advanced Processing (Stubs & New Functions)
+# =============================================================================
+
+"""
+ find_ppm_error_by_region(msi_data, region_masks, reference_peaks) -> Dict
+
+Calculates PPM error statistics for different spatial regions.
+`region_masks` is a Dict mapping region names (e.g., :tumor) to BitMatrix masks.
+"""
+function find_ppm_error_by_region(msi_data::MSIData, region_masks::Dict, reference_peaks::Dict)
+ regional_reports = Dict{Symbol, NamedTuple}()
+
+ width, height = msi_data.image_dims
+
+ for (region_name, mask) in region_masks
+ mask_height, mask_width = size(mask)
+ if mask_width != width || mask_height != height
+ @warn "Mask dimensions ($(mask_width)x$(mask_height)) for region '$region_name' do not match image dimensions ($(width)x$(height)). Skipping."
+ continue
+ end
+
+ indices = [
+ i for i in 1:length(msi_data.spectra_metadata)
+ if msi_data.spectra_metadata[i].x > 0 &&
+ msi_data.spectra_metadata[i].y > 0 &&
+ mask[msi_data.spectra_metadata[i].y, msi_data.spectra_metadata[i].x]
+ ]
+
+ if isempty(indices) continue end
+
+ # Call analyze_mass_accuracy with the specific indices for the region
+ regional_reports[region_name] = analyze_mass_accuracy(msi_data, reference_peaks; spectrum_indices=indices)
+ end
+ return regional_reports
+end
+
+"""
+ regional_calibration(msi_data, region_masks, reference_peaks)
+
+(STUB) Applies different calibration models to different spatial regions.
+"""
+function regional_calibration(msi_data::MSIData, region_masks::Dict, reference_peaks::Dict)
+ @warn "regional_calibration is a stub and not fully implemented."
+ # 1. For each region in region_masks:
+ # 2. Create a region-specific calibration model using `calibrate_spectra`.
+ # 3. Apply the model to all spectra within that region.
+ # 4. Return a new MSIData object or modified spectra vector.
+ return msi_data # Return unmodified for now
+end
+
+# =============================================================================
+# 8) Advanced Peak Quality & Adaptive Parameters
+# =============================================================================
+
+"""
+ calculate_peak_quality_metrics(peak_mz, peak_intensity, local_spectrum) -> Dict
+
+(STUB) Calculates advanced quality metrics for a single peak.
+"""
+function calculate_peak_quality_metrics(peak_mz, peak_intensity, local_spectrum)
+ # A full implementation would calculate:
+ # - Sharpness (e.g., ratio of height to FWHM)
+ # - Symmetry (e.g., ratio of left/right half-widths)
+ # - Signal-to-Noise (using local noise estimation)
+ # - Isolation score (how close are other peaks)
+ return Dict(:sharpness => 1.0, :symmetry => 1.0, :snr => 10.0)
+end
+
+
+
+
+"""
+ calculate_adaptive_bin_tolerance(ppm_error_distribution) -> Float64
+
+Calculates an appropriate binning tolerance based on observed mass accuracy.
+"""
+function calculate_adaptive_bin_tolerance(ppm_error_distribution::Vector{Float64})
+ if isempty(ppm_error_distribution)
+ return 20.0 # Default if no data
+ end
+ # A robust strategy: mean + 3 * std deviation to capture ~99.7% of peaks
+ return mean(ppm_error_distribution) + 3 * std(ppm_error_distribution)
+end
+
+
# =============================================================================
# 7) Plotting Helper
# =============================================================================
"""
- plot_stage_spectrum(mz, intensity; title, ...)
+ plot_stage_spectrum(mz, intensity; title, xlabel, ylabel) -> Figure
-Returns a `CairoMakie.Figure` for a single spectrum trace. The caller is responsible for saving.
+A helper function to generate a plot of a single spectrum using `CairoMakie`. This is
+useful for visualizing the output of different preprocessing steps.
+
+# Arguments
+- `mz::AbstractVector`: The m/z vector for the x-axis.
+- `intensity::AbstractVector`: The intensity vector for the y-axis.
+
+# Keyword Arguments
+- `title::AbstractString`: The title of the plot.
+- `xlabel::AbstractString`: The label for the x-axis. Defaults to "m/z".
+- `ylabel::AbstractString`: The label for the y-axis. Defaults to "Intensity".
+
+# Returns
+- `CairoMakie.Figure`: A `Figure` object containing the plot. The caller is responsible for displaying or saving it.
+
+# Example
+```julia
+using CairoMakie
+mz = 100:0.1:110;
+intensity = rand(length(mz));
+fig = plot_stage_spectrum(mz, intensity, title="My Spectrum");
+save("my_spectrum.png", fig);
+```
"""
function plot_stage_spectrum(mz::AbstractVector, intensity::AbstractVector;
title::AbstractString, xlabel::AbstractString="m/z",
@@ -395,255 +591,57 @@ function plot_stage_spectrum(mz::AbstractVector, intensity::AbstractVector;
return fig
end
-# =============================================================================
-# 8) Pipeline Orchestrator
-# =============================================================================
-
-"""
- run_preprocessing_pipeline(spectra; steps, params, on_stage)
-
-Executes a flexible preprocessing pipeline on a vector of spectra.
-
-# Arguments
-- `spectra`: A vector of `(mz, intensity)` tuples.
-- `steps`: A vector of symbols defining the pipeline order (e.g., `[:qc, :smooth, :baseline, :peaks, :bin]`).
-- `params`: A dictionary of parameters for each step.
-- `on_stage`: An optional callback function `on_stage(stage_symbol; idx, mz, intensity)` executed after each step for logging or visualization.
-
-# Returns
-- A `FeatureMatrix` if `:bin` is in the steps, otherwise the vector of processed spectra.
-"""
-function run_preprocessing_pipeline(spectra::Vector;
- steps::Vector{Symbol},
- params::Dict=Dict(),
- on_stage::Function=(;kwargs...)->nothing)
-
- processed = deepcopy(spectra) # Don't mutate the original input
- reference_peaks = nothing # For alignment
-
- _emit(stage::Symbol, idx::Int, mz, y) = on_stage(stage; idx=idx, mz=mz, intensity=y)
-
- for step in steps
- @info "Running step: $step"
-
- if step === :qc
- for (i, (mz, y)) in enumerate(processed)
- (isempty(mz) || isempty(y)) && continue
- _emit(:qc_raw, i, mz, y)
- qc_is_empty(mz, y) && @warn "Spectrum at index $i is empty."
- !qc_is_regular(mz) && @warn "m/z axis at index $i is not monotonic."
- end
-
- elseif step === :transform
- meth = get(params, :transform_method, :sqrt)
- for i in eachindex(processed)
- mz, y = processed[i]
- y_new = transform_intensity(y; method=meth)
- processed[i] = (mz, y_new)
- _emit(:transform, i, mz, y_new)
- end
-
- elseif step === :smooth
- win = get(params, :sg_window, 21)
- ord = get(params, :sg_order, 2)
- for i in eachindex(processed)
- mz, y = processed[i]
- y_smooth = smooth_spectrum(y; window=win, order=ord)
- processed[i] = (mz, y_smooth)
- _emit(:smooth, i, mz, y_smooth)
- end
-
- elseif step === :baseline
- iters = get(params, :snip_iterations, 100)
- for i in eachindex(processed)
- mz, y = processed[i]
- baseline = snip_baseline(y, iters)
- y_corrected = max.(0.0, y .- baseline)
- processed[i] = (mz, y_corrected)
- _emit(:baseline, i, mz, y_corrected)
- end
-
- elseif step === :normalize
- mode = get(params, :normalize_method, :tic)
- if mode === :tic
- for i in eachindex(processed)
- mz, y = processed[i]
- y_norm = tic_normalize(y)
- processed[i] = (mz, y_norm)
- _emit(:normalize, i, mz, y_norm)
- end
- elseif mode === :pqn
- # Note: PQN assumes spectra are on a common m/z grid.
- matrix = hcat([float.(p[2]) for p in processed]...)
- matrix_norm = pqn_normalize(matrix)
- for i in eachindex(processed)
- mz, _ = processed[i]
- processed[i] = (mz, view(matrix_norm, :, i))
- _emit(:normalize, i, mz, processed[i][2])
- end
- end
-
- elseif step === :peaks
- peak_results = Vector{Tuple{Vector{Float64},Vector{Float64}}}(undef, length(processed))
- hw = get(params, :peak_half_window, 10)
- snr = get(params, :peak_snr, 2.0)
- for (i, (mz, y)) in enumerate(processed)
- pk_mz, pk_int = detect_peaks_profile(mz, y; half_window=hw, snr_threshold=snr)
- peak_results[i] = (pk_mz, pk_int)
- # Emit with original mz axis but maybe stem plot of peaks?
- _emit(:peaks, i, pk_mz, pk_int)
- end
- processed = peak_results
- # Set reference for alignment
- reference_peaks = isempty(processed) ? nothing : processed[1][1]
-
- elseif step === :align
- reference_peaks === nothing && (@error "Alignment requires a :peaks step first."; continue)
- tol = get(params, :align_tolerance, 0.002)
- for i in 2:length(processed)
- tgt_peaks, intens = processed[i]
- warp_func = align_peaks_lowess(reference_peaks, tgt_peaks; tolerance=tol)
- processed[i] = (warp_func(tgt_peaks), intens)
- _emit(:align, i, processed[i][1], processed[i][2])
- end
-
- elseif step === :bin
- all_pks = [s[1] for s in processed]
- all_ints = [s[2] for s in processed]
- tol = get(params, :bin_tolerance, 0.002)
- freq = get(params, :bin_min_frequency, 0.25)
- mat, mz_bins = bin_peaks(all_pks, all_ints, tol; frequency_threshold=freq)
- return FeatureMatrix(mat, mz_bins, collect(1:length(processed)))
- end
- end
-
- @warn "Pipeline finished without a :bin step; returning processed spectra."
- return processed
-end
-
-"""
- run_preprocessing_pipeline(msi_data::MSIData, indices::Vector{Int}; steps, params, on_stage)
-
-Executes a flexible preprocessing pipeline on a subset of spectra from an MSIData object,
-with mode-aware logic for centroid and profile data.
-"""
-function run_preprocessing_pipeline(msi_data::MSIData, indices::Vector{Int};
- steps::Vector{Symbol},
- params::Dict=Dict(),
- on_stage::Function=(;kwargs...)->nothing)
-
- # This version of the pipeline is mode-aware.
- # It processes spectra directly from the MSIData object.
-
- processed_spectra = Vector{Tuple}(undef, length(indices))
-
- # First, load all spectra and apply initial steps that run on individual spectra
- for (i, spec_idx) in enumerate(indices)
- mz, intensity = GetSpectrum(msi_data, spec_idx)
- mode = msi_data.spectra_metadata[spec_idx].mode
-
- on_stage(:qc_raw; idx=spec_idx, mz=mz, intensity=intensity)
-
- for step in steps
- if step === :transform
- meth = get(params, :transform_method, :sqrt)
- intensity = transform_intensity(intensity; method=meth)
- on_stage(:transform; idx=spec_idx, mz=mz, intensity=intensity)
- elseif step === :smooth && mode == PROFILE
- win = get(params, :sg_window, 21)
- ord = get(params, :sg_order, 2)
- intensity = smooth_spectrum(intensity; window=win, order=ord)
- on_stage(:smooth; idx=spec_idx, mz=mz, intensity=intensity)
- elseif step === :baseline && mode == PROFILE
- iters = get(params, :snip_iterations, 100)
- baseline = snip_baseline(intensity, iters)
- intensity = max.(0.0, intensity .- baseline)
- on_stage(:baseline; idx=spec_idx, mz=mz, intensity=intensity)
- elseif step === :normalize
- norm_mode = get(params, :normalize_method, :tic)
- if norm_mode === :tic
- intensity = tic_normalize(intensity)
- on_stage(:normalize; idx=spec_idx, mz=mz, intensity=intensity)
- end
- end
- end
- processed_spectra[i] = (mz, intensity, mode) # Store mode for peak detection
- end
-
- # Now, handle steps that require all spectra (like PQN) or are the final steps
- final_result = nothing
- for step in steps
- if step === :normalize && get(params, :normalize_method, :tic) === :pqn
- # Note: PQN assumes spectra are on a common m/z grid.
- matrix = hcat([float.(p[2]) for p in processed_spectra]...)
- matrix_norm = pqn_normalize(matrix)
- for i in eachindex(processed_spectra)
- mz, _, mode = processed_spectra[i]
- processed_spectra[i] = (mz, view(matrix_norm, :, i), mode)
- on_stage(:normalize; idx=indices[i], mz=mz, intensity=processed_spectra[i][2])
- end
- elseif step === :peaks
- peak_results = Vector{Tuple{Vector{Float64},Vector{Float64}}}(undef, length(processed_spectra))
- hw = get(params, :peak_half_window, 10)
- snr = get(params, :peak_snr, 2.0)
- intensity_thresh = get(params, :peak_intensity_threshold, 0.0)
-
- for (i, (mz, y, mode)) in enumerate(processed_spectra)
- spec_idx = indices[i]
- if mode == PROFILE
- pk_mz, pk_int = detect_peaks_profile(mz, y; half_window=hw, snr_threshold=snr)
- else # CENTROID
- pk_mz, pk_int = detect_peaks_centroid(mz, y; intensity_threshold=intensity_thresh)
- end
- peak_results[i] = (pk_mz, pk_int)
- on_stage(:peaks; idx=spec_idx, mz=pk_mz, intensity=pk_int)
- end
- processed_spectra = peak_results # Now contains peak lists
-
- elseif step === :align
- # Alignment requires a reference peak list, typically from the first spectrum
- reference_peaks = isempty(processed_spectra) ? nothing : processed_spectra[1][1]
- if reference_peaks === nothing
- @error "Alignment requires a :peaks step first."; continue
- end
- tol = get(params, :align_tolerance, 0.002)
- for i in 2:length(processed_spectra)
- tgt_peaks, intens = processed_spectra[i]
- warp_func = align_peaks_lowess(reference_peaks, tgt_peaks; tolerance=tol)
- processed_spectra[i] = (warp_func(tgt_peaks), intens)
- on_stage(:align; idx=indices[i], mz=processed_spectra[i][1], intensity=processed_spectra[i][2])
- end
-
- elseif step === :bin
- all_pks = [s[1] for s in processed_spectra]
- all_ints = [s[2] for s in processed_spectra]
- tol = get(params, :bin_tolerance, 0.002)
- freq = get(params, :bin_min_frequency, 0.25)
- mat, mz_bins = bin_peaks(all_pks, all_ints, tol; frequency_threshold=freq)
- final_result = FeatureMatrix(mat, mz_bins, indices)
- break # Binning is the last step
- end
- end
-
- if final_result !== nothing
- return final_result
- else
- @warn "Pipeline finished without a :bin step; returning processed spectra."
- return processed_spectra
- end
-end
# =============================================================================
# 9) Quality Control Metrics
# =============================================================================
"""
- calculate_ppm_error(measured_mz::Float64, theoretical_mz::Float64) -> Float64
+ get_common_calibration_standards(type::Symbol) -> Dict{Float64, String}
-Calculates mass accuracy in parts-per-million (PPM).
+Returns a dictionary of common m/z values for calibration standards.
+
+# Arguments
+- `type::Symbol`: The type of standards to return. Currently supports `:maldi_pos`.
+
+# Returns
+- `Dict{Float64, String}`: A dictionary mapping theoretical m/z to compound names.
+"""
+function get_common_calibration_standards(type::Symbol)
+ if type == :maldi_pos
+ # Common peptide/protein standards for positive mode MALDI
+ return Dict{Float64, String}(
+ 1046.5420 => "Angiotensin II",
+ 1060.5690 => "Bradykinin",
+ 1296.6853 => "Angiotensin I",
+ 2465.1989 => "ACTH clip (1-24)"
+ )
+ else
+ @warn "Unsupported calibration standard type: $type. Returning empty dictionary."
+ return Dict{Float64, String}()
+ end
+end
+
+"""
+ calculate_ppm_error(measured_mz::Real, theoretical_mz::Real) -> Float64
+
+Calculates the mass accuracy error in parts-per-million (PPM) between a measured
+and a theoretical m/z value.
+
+# Arguments
+- `measured_mz::Real`: The experimentally measured m/z value.
+- `theoretical_mz::Real`: The known, theoretical m/z value of a compound.
+
+# Returns
+- `Float64`: The calculated PPM error. Returns `Inf` if `theoretical_mz` is zero.
# Formula
-PPM = 10⁶ × |measured_mz - theoretical_mz| / theoretical_mz
+`PPM = 10^6 * |measured_mz - theoretical_mz| / theoretical_mz`
+
+# Example
+```julia
+calculate_ppm_error(100.005, 100.0) # returns 50.0
+```
"""
function calculate_ppm_error(measured_mz::Real, theoretical_mz::Real)
if theoretical_mz == 0
@@ -653,30 +651,38 @@ function calculate_ppm_error(measured_mz::Real, theoretical_mz::Real)
end
"""
- calculate_ppm_error_bulk(measured_mz::Vector{Float64}, theoretical_mz::Vector{Float64}) -> Vector{Float64}
+ calculate_ppm_error_bulk(measured_mz::Vector{<:Real}, theoretical_mz::Vector{<:Real}) -> Vector{Float64}
-Calculates PPM errors for multiple mass values.
+Calculates PPM errors for multiple pairs of measured and theoretical mass values.
+
+# Arguments
+- `measured_mz::Vector{<:Real}`: A vector of experimentally measured m/z values.
+- `theoretical_mz::Vector{<:Real}`: A vector of known, theoretical m/z values.
+
+# Returns
+- `Vector{Float64}`: A vector containing the calculated PPM error for each pair.
"""
function calculate_ppm_error_bulk(measured_mz::Vector{Real}, theoretical_mz::Vector{Real})
return [calculate_ppm_error(m, t) for (m, t) in zip(measured_mz, theoretical_mz)]
end
"""
- calculate_resolution_fwhm(mz::Float64, profile_mz::Vector{Float64},
- profile_intensity::Vector{Float64}) -> Float64
+ calculate_resolution_fwhm(mz::Real, profile_mz::AbstractVector{<:Real}, profile_intensity::AbstractVector{<:Real}) -> Float64
-Calculates mass resolution using Full Width at Half Maximum (FWHM).
-
-# Formula
-Resolution = m / Δm, where Δm is FWHM
+Calculates the mass resolution of a peak in profile-mode data using the Full Width at
+Half Maximum (FWHM) method. Resolution is a measure of an instrument's ability to
+distinguish between two peaks of slightly different mass-to-charge ratios.
# Arguments
-- `mz`: Peak centroid m/z
-- `profile_mz`: Full m/z array from profile data
-- `profile_intensity`: Full intensity array from profile data
+- `mz::Real`: The m/z value of the peak's centroid.
+- `profile_mz::AbstractVector{<:Real}`: The full m/z array from the profile-mode spectrum.
+- `profile_intensity::AbstractVector{<:Real}`: The full intensity array from the profile-mode spectrum.
# Returns
-Resolution or NaN if cannot be calculated
+- `Float64`: The calculated resolution (`m / Δm`). Returns `NaN` if the FWHM cannot be determined (e.g., peak is at the edge of the spectrum).
+
+# Formula
+`Resolution = m / Δm`, where `Δm` is the FWHM.
"""
function calculate_resolution_fwhm(mz::Real, profile_mz::AbstractVector{<:Real},
profile_intensity::AbstractVector{<:Real})
@@ -733,36 +739,48 @@ function find_first_below(v::AbstractVector{<:Real}, threshold::Real)
end
"""
- analyze_mass_accuracy(msi_data, reference_peaks; ppm_tolerance=20.0)
+ analyze_mass_accuracy(msi_data, reference_peaks; ppm_tolerance=5.0, sample_spectra=100) -> NamedTuple
-Analyzes mass accuracy across the dataset using known reference peaks.
+Analyzes the mass accuracy across a sample of spectra from an MSI dataset by comparing
+detected peaks against a list of known reference m/z values. It calculates PPM error
+statistics and suggests an optimal PPM tolerance for future analyses.
# Arguments
-- `msi_data`: Your MSI dataset
-- `reference_peaks`: Dict of theoretical m/z values -> compound names
-- `ppm_tolerance`: Initial tolerance for peak matching
+- `msi_data`: An `MSIData` object containing the dataset.
+- `reference_peaks::Dict{Float64,String}`: A dictionary mapping theoretical m/z values to compound names.
+
+# Keyword Arguments
+- `ppm_tolerance::Float64`: The initial PPM tolerance used to match detected peaks to reference peaks. Defaults to 5.0.
+- `sample_spectra::Int`: The number of spectra to sample from the dataset for the analysis. Defaults to 100.
# Returns
-Comprehensive mass accuracy report
+- `NamedTuple`: A comprehensive report containing statistics like `mean_ppm`, `std_ppm`, `optimal_ppm`, `n_matches`, and lists of matched peaks and all PPM errors.
"""
function analyze_mass_accuracy(msi_data, reference_peaks::Dict{Float64,String};
- ppm_tolerance::Float64=5.0, sample_spectra=100)
+ ppm_tolerance::Float64=5.0, sample_spectra=100,
+ spectrum_indices=nothing)
println("\n[ MASS ACCURACY ANALYSIS ]")
println("PPM tolerance: ", ppm_tolerance, " ppm")
- println("Spectra to sample: ", sample_spectra)
theoretical_mz = sort(collect(keys(reference_peaks)))
ppm_errors = Float64[]
matched_peaks = Tuple{Float64,Float64,String}[] # (theoretical, measured, compound)
- # Sample spectra across the dataset
- if sample_spectra >= length(msi_data.spectra_metadata)
- spectrum_indices = 1:length(msi_data.spectra_metadata)
+ # Determine which spectra to process
+ local indices_to_process
+ if spectrum_indices === nothing
+ println("Spectra to sample: ", sample_spectra)
+ if sample_spectra >= length(msi_data.spectra_metadata)
+ indices_to_process = 1:length(msi_data.spectra_metadata)
+ else
+ indices_to_process = round.(Int, range(1, length(msi_data.spectra_metadata), length=sample_spectra))
+ end
else
- spectrum_indices = round.(Int, range(1, length(msi_data.spectra_metadata), length=sample_spectra))
+ indices_to_process = spectrum_indices
+ println("Analyzing $(length(indices_to_process)) specified spectra.")
end
- for idx in spectrum_indices
+ for idx in indices_to_process
mz, intensity = GetSpectrum(msi_data, idx)
# Detect peaks in this spectrum
@@ -788,7 +806,7 @@ function analyze_mass_accuracy(msi_data, reference_peaks::Dict{Float64,String};
if isempty(ppm_errors)
@warn "No peaks matched within $ppm_tolerance ppm tolerance"
- return (mean_ppm=NaN, std_ppm=NaN, min_ppm=NaN, max_ppm=NaN, optimal_ppm=NaN, n_matches=0, matched_peaks=[], all_ppm_errors=[])
+ return (mean_ppm=NaN, std_ppm=NaN, min_ppm=NaN, max_ppm=NaN, optimal_ppm=NaN, n_matches=0, matched_peaks=[], all_ppm_errors=Float64[])
end
# Calculate statistics
@@ -813,62 +831,28 @@ function analyze_mass_accuracy(msi_data, reference_peaks::Dict{Float64,String};
end
"""
- get_common_calibration_standards(standard_type::Symbol)
+ generate_qc_report(msi_data, filename; reference_peaks, output_dir, sample_spectra) -> Tuple
-Returns common calibration masses for different instrument types.
+Generates a comprehensive Quality Control (QC) report for an MSI dataset. It assesses
+mass accuracy and resolution, saving the results to CSV files and a summary text file.
-# Supported standards
-- `:maldi_pos`: Common MALDI-TOF positive mode calibrants
-- `:maldi_neg`: Common MALDI-TOF negative mode calibrants
-- `:esi_pos`: ESI positive mode calibrants
-- `:lcms`: LC-MS commonly used standards
+# Arguments
+- `msi_data`: The `MSIData` object to be analyzed.
+- `filename::String`: The name of the input data file, used for reporting.
+
+# Keyword Arguments
+- `reference_peaks::Dict`: A dictionary of reference peaks for mass accuracy analysis. If `nothing`, defaults are loaded using `get_common_calibration_standards`.
+- `output_dir::String`: The directory where the report files will be saved. Defaults to `"qc_results"`.
+- `sample_spectra::Int`: The number of spectra to sample for the analysis. Defaults to 100.
+
+# Returns
+- `Tuple`: A tuple containing the detailed `accuracy_report` (NamedTuple) and `resolution_results` (Vector).
+
+# Side Effects
+- Creates a directory specified by `output_dir`.
+- Writes `mass_accuracy_results.csv`, `resolution_results.csv`, and `qc_summary.txt` into the output directory.
"""
-function get_common_calibration_standards(standard_type::Symbol=:maldi_pos)
- standards = Dict{Float64,String}()
-
- if standard_type == :maldi_pos
- standards = Dict(
- 104.10754 => "C5H4N2 (Imidazole)",
- 175.11995 => "C6H15O4P (Glycerophosphocholine fragment)",
- 226.15687 => "C10H20NO4P (Phosphocholine)",
- 322.04810 => "[Glu1]-Fibrinopeptide B fragment",
- 379.09247 => "C12H22O11 (Sucrose)",
- 515.32539 => "C26H52NO7P (PC(16:0/0:0))",
- 622.02896 => "C20H12O5S2 (1-Hydroxypyrene-3,6,8-trisulfate)",
- 757.39917 => "C37H74NO8P (PC(34:1))",
- 1046.54198 => "Angiotensin I",
- 1296.68477 => "ACTH clip 1-17",
- 1570.67744 => "ACTH clip 18-39",
- 2465.19829 => "ACTH clip 7-38"
- )
- elseif standard_type == :maldi_neg
- standards = Dict(
- 112.98563 => "C2F3O2 (Trifluoroacetate)",
- 152.99568 => "C2F6S (Perfluoroethylsulfonate)",
- 214.00166 => "C4F7O2 (Heptafluorobutyrate)",
- 264.93278 => "C6F6 (Hexafluorobenzene)",
- 362.96198 => "C8F15O2 (Perfluorooctanoate)",
- 466.96714 => "C10F17O2S (Perfluorooctanesulfonate)"
- )
- elseif standard_type == :esi_pos
- standards = Dict(
- 118.08626 => "C5H12NO2 (Valine)",
- 175.11900 => "C6H15O4P (Phosphocholine fragment)",
- 524.26496 => "C23H48NO7P (LysoPC(16:0))",
- 622.02896 => "C20H12O5S2 (Standard)",
- 922.00980 => "C18H18O6N3S3 (Ultramark 1621)"
- )
- end
-
- return standards
-end
-
-"""
- generate_qc_report(msi_data; reference_peaks, output_dir)
-
-Generates a comprehensive QC report including mass accuracy and resolution.
-"""
-function generate_qc_report(msi_data, filename::String; reference_peaks=nothing, output_dir="qc_results", sample_spectra=100)
+function generate_qc_report(msi_data, filename::String; reference_peaks=nothing, output_dir="qc_results", sample_spectra=100, spectrum_indices=nothing)
println("\n[ QC REPORT GENERATION ]")
println("Input file: ", filename)
println("Output directory: ", output_dir)
@@ -883,7 +867,7 @@ function generate_qc_report(msi_data, filename::String; reference_peaks=nothing,
println("Using $(length(reference_peaks)) reference masses")
# 1. Analyze mass accuracy
- accuracy_report = analyze_mass_accuracy(msi_data, reference_peaks, sample_spectra=sample_spectra)
+ accuracy_report = analyze_mass_accuracy(msi_data, reference_peaks; sample_spectra=sample_spectra, spectrum_indices=spectrum_indices)
println("\n" * "="^50)
println("MASS ACCURACY REPORT")
@@ -976,4 +960,3 @@ function generate_qc_report(msi_data, filename::String; reference_peaks=nothing,
println("\nQC report saved to: $output_dir")
return accuracy_report, resolution_results
end
-
diff --git a/src/PreprocessingVersioning.jl b/src/PreprocessingVersioning.jl
new file mode 100644
index 0000000..aa26537
--- /dev/null
+++ b/src/PreprocessingVersioning.jl
@@ -0,0 +1,446 @@
+# src/PreprocessingVersioning.jl
+
+"""
+This module provides a versioning helper functions in our data preprocessing, module
+inspired by the functionality of the R package MALDIquant. It includes
+functions for quality control, intensity transformation, smoothing, baseline correction,
+normalization, peak picking, alignment, and feature matrix generation.
+
+This module is unfinished and may not end up in the release version.
+"""
+
+# =============================================================================
+# Dependencies
+# =============================================================================
+
+using UUIDs, JLD2, CodecBase # For preprocessing versioning
+
+# =============================================================================
+# Data Structures
+# =============================================================================
+
+"""
+ VersionedSpectralData
+
+A struct to hold a snapshot of spectral data at a specific point in the
+preprocessing workflow. Each transformation creates a new instance of this
+struct, forming a history of all operations.
+
+# Fields
+- `version_id::String`: A unique identifier for this version of the data.
+- `parent_id::String`: The identifier of the version from which this one was derived.
+- `spectra::Vector`: The spectral data itself, typically a vector of `(mz, intensity)` tuples.
+- `processing_step::String`: A description of the operation that created this version (e.g., "smooth", "baseline").
+- `timestamp::DateTime`: The time at which this version was created.
+- `parameters::Dict`: A dictionary of the parameters used in the processing step.
+"""
+struct VersionedSpectralData
+ version_id::String
+ parent_id::String
+ spectra::Vector
+ processing_step::String
+ timestamp::DateTime
+ parameters::Dict
+end
+
+"""
+ MSISession
+
+Manages the state of a preprocessing session, including the history of all
+data versions and the current working version.
+
+# Fields
+- `session_id::String`: A unique identifier for the entire session.
+- `original_data::VersionedSpectralData`: The initial, unprocessed spectral data.
+- `history::Vector{VersionedSpectralData}`: A chronological list of all data versions created during the session.
+- `current_version::VersionedSpectralData`: The version of the data currently being worked on or displayed.
+- `processing_steps::Vector{String}`: A human-readable log of the processing steps applied.
+"""
+mutable struct MSISession
+ session_id::String
+ original_data::VersionedSpectralData
+ history::Vector{VersionedSpectralData}
+ current_version::VersionedSpectralData
+ processing_steps::Vector{String}
+end
+
+# =============================================================================
+# Session Management
+# =============================================================================
+
+"""
+ initialize_session(raw_spectra::Vector, session_name::String) -> MSISession
+
+Creates a new preprocessing session from raw spectral data.
+
+# Arguments
+- `raw_spectra::Vector`: A vector of `(mz, intensity)` tuples representing the initial dataset.
+- `session_name::String`: A name for the session.
+
+# Returns
+- `MSISession`: A new session object initialized with the provided data.
+"""
+function initialize_session(raw_spectra::Vector, session_name::String)::MSISession
+ original = VersionedSpectralData(
+ string(uuid4()),
+ "root",
+ raw_spectra,
+ "Original Data",
+ now(),
+ Dict()
+ )
+ return MSISession(session_name, original, [original], original, [])
+end
+
+"""
+ get_processing_history(session::MSISession) -> Vector{String}
+
+Returns a list of the processing steps applied during the session.
+"""
+get_processing_history(session::MSISession) = [v.processing_step for v in session.history]
+
+"""
+ undo_last_step!(session::MSISession)
+
+Reverts the session to the state before the last processing step was applied.
+This modifies the session in place.
+"""
+function undo_last_step!(session::MSISession)
+ if length(session.history) > 1
+ pop!(session.history)
+ pop!(session.processing_steps)
+ session.current_version = last(session.history)
+ end
+ return session
+end
+
+"""
+ revert_to_version!(session::MSISession, version_id::String)
+
+Reverts the session to a specific version in its history, discarding all
+subsequent changes. This modifies the session in place.
+"""
+function revert_to_version!(session::MSISession, version_id::String)
+ target_idx = findfirst(v -> v.version_id == version_id, session.history)
+ if !isnothing(target_idx)
+ session.history = session.history[1:target_idx]
+ session.current_version = session.history[end]
+ # Also truncate the descriptive processing steps
+ num_steps_to_keep = max(0, target_idx - 1)
+ session.processing_steps = session.processing_steps[1:num_steps_to_keep]
+ end
+ return session
+end
+
+# =============================================================================
+# Versioned Preprocessing
+# =============================================================================
+
+"""
+ apply_processing_step(spectra::Vector, step::Symbol, params::Dict) -> Vector
+
+A dispatcher that applies a single, specified preprocessing function to a set of spectra.
+This is the core function called by the versioned workflow.
+
+# Arguments
+- `spectra::Vector`: The input spectral data.
+- `step::Symbol`: The symbol representing the processing step (e.g., `:smooth`, `:baseline`).
+- `params::Dict`: A dictionary of parameters for the step.
+
+# Returns
+- `Vector`: The processed spectral data.
+"""
+function apply_processing_step(spectra::Vector, step::Symbol, params::Dict)
+ processed = deepcopy(spectra)
+ params_sym = Dict(Symbol(k) => v for (k,v) in params) # Ensure keys are symbols
+
+ if step === :transform
+ for i in eachindex(processed)
+ mz, y = processed[i]
+ processed[i] = (mz, transform_intensity(y; params_sym...))
+ end
+ elseif step === :smooth
+ for i in eachindex(processed)
+ mz, y = processed[i]
+ processed[i] = (mz, smooth_spectrum(y; params_sym...))
+ end
+ elseif step === :baseline
+ for i in eachindex(processed)
+ mz, y = processed[i]
+ baseline = snip_baseline(y; params_sym...)
+ processed[i] = (mz, max.(0.0, y .- baseline))
+ end
+ elseif step === :normalize
+ mode = get(params, :normalize_method, :tic)
+ if mode === :pqn
+ matrix = hcat([float.(p[2]) for p in processed]...)
+ matrix_norm = pqn_normalize(matrix)
+ for i in eachindex(processed)
+ mz, _ = processed[i]
+ processed[i] = (mz, view(matrix_norm, :, i))
+ end
+ else # :tic or :median
+ norm_func = (mode === :tic) ? tic_normalize : median_normalize
+ for i in eachindex(processed)
+ mz, y = processed[i]
+ processed[i] = (mz, norm_func(y))
+ end
+ end
+ elseif step === :peaks
+ peak_results = Vector{Tuple{Vector{Float64},Vector{Float64}}}(undef, length(processed))
+ for (i, (mz, y)) in enumerate(processed)
+ pk_mz, pk_int = detect_peaks_profile(mz, y; params_sym...)
+ peak_results[i] = (pk_mz, pk_int)
+ end
+ return peak_results
+ elseif step === :peaks_wavelet
+ peak_results = Vector{Tuple{Vector{Float64},Vector{Float64}}}(undef, length(processed))
+ for (i, (mz, y)) in enumerate(processed)
+ pk_mz, pk_int = detect_peaks_wavelet(mz, y; params_sym...)
+ peak_results[i] = (pk_mz, pk_int)
+ end
+ return peak_results
+ elseif step === :calibrate
+ return calibrate_spectra(processed; params_sym...)
+ else
+ @warn "Unsupported processing step: $step"
+ end
+ return processed
+end
+
+"""
+ run_versioned_preprocessing(session::MSISession, step::Symbol, params::Dict) -> VersionedSpectralData
+
+Creates a new version of spectral data by applying a single processing step to the
+current version in the session.
+
+# Arguments
+- `session::MSISession`: The current processing session.
+- `step::Symbol`: The processing step to apply.
+- `params::Dict`: Parameters for the processing step.
+
+# Returns
+- `VersionedSpectralData`: A new data version with the transformation applied.
+"""
+function run_versioned_preprocessing(session::MSISession, step::Symbol, params::Dict)
+ # Create a new version based on the current one
+ parent_version = session.current_version
+ new_version_id = string(uuid4())
+
+ # Apply the processing step
+ new_spectra = apply_processing_step(parent_version.spectra, step, params)
+
+ # Create the new versioned data object
+ new_version = VersionedSpectralData(
+ new_version_id,
+ parent_version.version_id,
+ new_spectra,
+ string(step),
+ now(),
+ params
+ )
+
+ return new_version
+end
+
+"""
+ apply_processing_with_version!(session::MSISession, step::Symbol, ui_params::Dict)
+
+A high-level function to apply a processing step, create a new version, and
+update the session state in place.
+
+# Arguments
+- `session::MSISession`: The session to modify.
+- `step::Symbol`: The processing step to apply.
+- `ui_params::Dict`: A dictionary of parameters, typically from a UI.
+"""
+function apply_processing_with_version!(session::MSISession, step::Symbol, ui_params::Dict)
+ # In a real app, you would validate and convert UI params here.
+ processing_params = ui_params
+
+ # Create the new version
+ new_version = run_versioned_preprocessing(session, step, processing_params)
+
+ # Update the session
+ push!(session.history, new_version)
+ session.current_version = new_version
+
+ # Describe the step for the history log
+ step_description = "$(string(step)) with params: " * join(["$k=$v" for (k,v) in processing_params], ", ")
+ push!(session.processing_steps, step_description)
+
+ return session
+end
+
+# =============================================================================
+# Data Serialization & Export
+# =============================================================================
+
+"""
+ save_spectral_version(data::VersionedSpectralData, filepath::String)
+
+Saves a `VersionedSpectralData` object to a file using the JLD2 format.
+"""
+function save_spectral_version(data::VersionedSpectralData, filepath::String)
+ JLD2.save_object(filepath, data)
+end
+
+"""
+ load_spectral_version(filepath::String) -> VersionedSpectralData
+
+Loads a `VersionedSpectralData` object from a JLD2 file.
+"""
+function load_spectral_version(filepath::String)::VersionedSpectralData
+ return JLD2.load_object(filepath)
+end
+
+"""
+ _encode_base64(data::AbstractVector{<:Real}) -> String
+
+Helper function to convert a numeric vector into a Base64 encoded string.
+"""
+function _encode_base64(data::AbstractVector{T}) where T <: Real
+ bytes = reinterpret(UInt8, data)
+ return base64encode(bytes)
+end
+
+"""
+ export_to_mzml(session::MSISession, filepath::String)
+
+Exports the current version of the spectral data in a session to a standard
+.mzML file.
+
+# Arguments
+- `session::MSISession`: The current session.
+- `filepath::String`: The path for the output .mzML file.
+"""
+function export_to_mzml(session::MSISession, filepath::String)
+ spectra_to_export = session.current_version.spectra
+
+ open(filepath, "w") do f
+ # XML Header
+ write(f, "\n")
+ write(f, "\n")
+
+ # CV List
+ write(f, " \n")
+ write(f, " \n")
+ write(f, " \n")
+ write(f, " \n")
+
+ # File Description
+ write(f, " \n \n")
+ write(f, " \n")
+ write(f, " \n \n")
+
+ # Run and Spectrum List
+ write(f, " \n")
+ write(f, " \n")
+
+ for (i, (mz, intensity)) in enumerate(spectra_to_export)
+ # Ensure data is in the correct format for encoding
+ mz_64 = convert(Vector{Float64}, mz)
+ int_32 = convert(Vector{Float32}, intensity)
+
+ # Base64 encode the binary data
+ mz_b64 = _encode_base64(mz_64)
+ int_b64 = _encode_base64(int_32)
+
+ write(f, " \n")
+ write(f, " \n")
+ # Assuming profile mode for processed data, could be made dynamic
+ write(f, " \n")
+
+ write(f, " \n")
+
+ # m/z array
+ write(f, " \n")
+ write(f, " \n")
+ write(f, " \n")
+ write(f, " \n")
+ write(f, " $(mz_b64)\n")
+ write(f, " \n")
+
+ # Intensity array
+ write(f, " \n")
+ write(f, " \n")
+ write(f, " \n")
+ write(f, " \n")
+ write(f, " $(int_b64)\n")
+ write(f, " \n")
+
+ write(f, " \n")
+ write(f, " \n")
+ end
+
+ write(f, " \n")
+ write(f, " \n")
+
+ # Dummy instrument and data processing info
+ write(f, " \n")
+ write(f, " \n")
+ write(f, " \n")
+ write(f, " \n")
+ write(f, " \n")
+ write(f, " \n")
+
+ write(f, "\n")
+ end
+ println("Successfully exported current data to $filepath")
+end
+
+
+# =============================================================================
+# Stubs for UI Integration
+# =============================================================================
+
+"""
+ validate_ui_parameters(step::Symbol, ui_params::Dict) -> Tuple{Bool, String}
+
+(Stub) Validates parameters from a UI before they are used in a processing step.
+"""
+function validate_ui_parameters(step::Symbol, ui_params::Dict)::Tuple{Bool, String}
+ # In a real implementation, this would check types, ranges, etc.
+ # For example, for :smooth, ensure 'sg_window' is an odd integer.
+ println("Validating parameters for step: $step")
+ return (true, "Parameters are valid.")
+end
+
+"""
+ run_processing_with_progress(session, step, params, progress_callback)
+
+(Stub) A wrapper for running a processing step that includes a progress reporting callback.
+"""
+function run_processing_with_progress(session, step, params, progress_callback)
+ progress_callback(0.0, "Starting $step...")
+
+ # This is a simplified example. Real implementation would need to
+ # hook into the loops inside apply_processing_step.
+ new_version = run_versioned_preprocessing(session, step, params)
+
+ progress_callback(1.0, "Finished $step.")
+ return new_version
+end
+
+"""
+ safe_processing_application!(session::MSISession, step::Symbol, params::Dict)
+
+(Stub) A safe wrapper to apply a processing step that includes error handling and rollback.
+"""
+function safe_processing_application!(session::MSISession, step::Symbol, params::Dict)
+ num_history = length(session.history)
+ try
+ println("Safely applying step: $step")
+ apply_processing_with_version!(session, step, params)
+ catch e
+ @error "Processing step '$step' failed!" exception=(e, catch_backtrace())
+ # Roll back to the previous state
+ while length(session.history) > num_history
+ undo_last_step!(session)
+ end
+ println("Session has been rolled back to the previous state.")
+ return session # Return the rolled-back session
+ end
+ return session # Return the updated session
+end
+
diff --git a/test/readme.md b/test/readme.md
index 197af0d..1965e10 100644
--- a/test/readme.md
+++ b/test/readme.md
@@ -40,10 +40,10 @@ Before running the tests, you must edit the `test/run_tests.jl` file to point to
2. **Execute the Different Test Scripts**:
Run the following command from the project's root directory. This will install the necessary dependencies and run the tests.
```bash
- julia --project=. test/run_tests.jl
+ julia --threads auto --project=. test/run_tests.jl
```
```bash
- julia --project=. test/run_preprocessing.jl
+ julia --threads auto --project=. test/run_preprocessing.jl
```
3. **Check the Results**:
diff --git a/test/run_preprocessing.jl b/test/run_preprocessing.jl
index 4ea526a..bc70023 100644
--- a/test/run_preprocessing.jl
+++ b/test/run_preprocessing.jl
@@ -1,468 +1,371 @@
# test/run_preprocessing.jl
# ===================================================================
-# Test Environment for the Preprocessing.jl Module
+# Preprocessing Test Suite for JuliaMSI
# ===================================================================
-# This script tests the full preprocessing pipeline on single spectra
-# and total spectra from both .mzML and .imzML files.
-# It generates an overlay plot showing all preprocessing stages and
-# saves the resulting feature matrix to a CSV file.
+# This script provides a customizable workflow to test and visualize
+# the effects of different preprocessing steps on individual spectra
+# from .mzML or .imzML files.
#
# Instructions:
-# 1. Ensure the file paths in the "CONFIG" section are correct.
+# 1. Configure the file paths and parameters in the "CONFIG" section.
# 2. Run the script from the project's root directory:
-# julia test/run_preprocessing.jl
-# 3. Check the `test/results/` folder for output plots and CSVs.
+# julia --threads auto --project=. test/run_preprocessing.jl
+# 3. Check the configured `RESULTS_DIR` for output plots and CSV files.
# ===================================================================
using Printf
using CairoMakie
+using DataFrames
+using CSV
+using Statistics
+using Interpolations
import Pkg
-using DataFrames # For saving FeatureMatrix to CSV
-using CSV # For saving FeatureMatrix to CSV
-using Statistics # For mean()
-# --- Load Modules ---
-# Activate the project environment to access dependencies
+# --- Load the MSI_src Module ---
+# Activate the project environment at the parent directory of this test script
Pkg.activate(joinpath(@__DIR__, ".."))
-using MSI_src # This brings in Preprocessing.jl functions via export
+using MSI_src
+
# ===================================================================
-# CONFIG: Test files and parameters
+# CONFIG: CUSTOMIZE YOUR PREPROCESSING WORKFLOW HERE
# ===================================================================
-# --- Test Files ---
-# An mzML file for testing spectrum-based processing
-# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML"
-# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/set de datos MS/Leaf_profile_LD_LTP_MS.mzML"
-const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/set de datos MS/Escopolamina_tuneo_fraq_20ev.mzML"
-#const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/set de datos MS/Atropina_tuneo_fraq_20ev.mzML"
+# --- Input and Output ---
+const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML" # Can be .mzML or .imzML
+const RESULTS_DIR = "test/results/preprocessing"
+const NUM_SPECTRA_TO_PROCESS = nothing # Set to `nothing` to process all spectra
-const MZML_SPECTRUM_ID = 1
+# --- Internal Standards for Calibration and QC ---
+# replace these with m/z values of known compounds present in your dataset.
+const INTERNAL_STANDARDS = Dict(
+ "P13" => 432.6584,
+ "P15" => 464.6059,
+ "P17" => 526.5534,
+ "P21" => 650.4485,
+ "P25" => 774.3435,
+ "P29" => 898.2385,
+ "P31" => 950.1861,
+ "P33" => 1022.1336,
+ "P37" => 1146.0286,
+ "P45" => 1593.8187,
+ "Unknown 1" => 772.433,
+ "Unknown 2" => 772.5253
+)
-# An imzML file for testing
-# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.imzML"
-# const IMZML_COORDS = (50, 50)
-const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML"
-const IMZML_COORDS = (1997, 639)
+# --- Preprocessing Step Parameters ---
-# --- Output Directory ---
-const RESULTS_DIR = "test/results"
+# Step 0: Quality Control (QC)
+const QC_PARAMS = (
+ ppm_tolerance = 5.0, # PPM tolerance for matching internal standards
+)
+
+# Step 1: Calibration
+const CALIBRATION_PARAMS = (
+ enabled = true,
+ ppm_tolerance = 5.0, # PPM tolerance for finding calibration peaks
+)
+
+# Step 2: Smoothing
+# Note: Smoothing is less effective and often unnecessary for centroid data.
+const SMOOTHING_PARAMS = (
+ enabled = false,
+ window = 9,
+ order = 2,
+)
+
+# Step 3: Baseline Correction
+# Note: Baseline correction is less effective and often unnecessary for centroid data.
+const BASELINE_CORRECTION_PARAMS = (
+ enabled = false,
+ iterations = 100,
+)
+
+# Step 4: Normalization
+const NORMALIZATION_PARAMS = (
+ enabled = true,
+ method = :tic, # :tic, :median, or :none
+)
+
+# Step 5: Peak Detection
+const PEAK_DETECTION_PARAMS = (
+ enabled = true,
+ method = :centroid, # Use :centroid for centroided data, :profile for profile data
+ # --- Parameters for :centroid method ---
+ intensity_threshold = 0.0, # Filters out peaks below this absolute intensity
+ # --- Parameters for :profile method ---
+ half_window = 10,
+ snr_threshold = 3.0,
+ min_peak_prominence = 0.05,
+ merge_peaks_tolerance = 0.002,
+)
+
+# Step 6: Peak Alignment (Warping)
+# This is performed after collecting peaks from all spectra.
+const PEAK_ALIGNMENT_PARAMS = (
+ enabled = true,
+ tolerance = 0.002,
+ tolerance_unit = :mz, # :mz or :ppm
+ min_matched_peaks = 3, # Lowered for sparse centroid data
+)
+
+# Step 7: Peak Binning
+const PEAK_BINNING_PARAMS = (
+ enabled = true,
+ tolerance = 0.1,
+ tolerance_unit = :mz, # :mz or :ppm
+ frequency_threshold = 0.1, # Min fraction of spectra a peak must be in
+)
# ===================================================================
-# HELPER FUNCTIONS FOR PLOTTING
+# HELPER FUNCTIONS
# ===================================================================
"""
- plot_overlay_stages(collected_data, output_path, title)
+ plot_spectrum_on_fig(fig_pos, mz, intensity, title)
-Creates a single plot overlaying spectra from different preprocessing stages.
+Helper function to plot a spectrum on a specific position of a CairoMakie figure.
"""
-function plot_overlay_stages(collected_data, output_path, title)
- fig = Figure(size = (1400, 800))
- ax = Axis(fig[1, 1], title=title, xlabel="m/z", ylabel="Intensity")
-
- colors = Makie.wong_colors() # A good set of distinct colors
-
- for (i, (stage, mz, intensity)) in enumerate(collected_data)
- color = colors[mod1(i, length(colors))] # Cycle through colors
-
- # Plot the spectrum as a line
- lines!(ax, mz, intensity, color=color, label=string(stage))
-
- # If it's the peaks stage, also mark the peak tops
- if stage == :peaks
- scatter!(ax, mz, intensity, color=color, marker=:circle, markersize=8, label="$(string(stage)) (tops)")
- end
- end
- axislegend(ax, position=:rt) # Right top position
- save(output_path, fig)
- println("SUCCESS: Overlay plot saved to $output_path")
+function plot_spectrum_on_fig(fig_pos, mz, intensity, title)
+ ax = Axis(fig_pos, title=title, xlabel="m/z", ylabel="Intensity")
+ lines!(ax, mz, intensity)
end
# ===================================================================
-# TEST DEFINITIONS
+# MAIN PROCESSING SCRIPT
# ===================================================================
-"""
- test_full_pipeline(msi_data, spectrum_id; output_dir, file_type_prefix)
+function run_preprocessing_suite()
+ println("Starting Preprocessing Test Suite...")
+ mkpath(RESULTS_DIR)
-Tests the full preprocessing pipeline on a single spectrum and saves a plot
-for each intermediate step using the `on_stage` callback.
-`spectrum_id` can be an `Int` (for mzML) or a `Tuple{Int, Int}` (for imzML).
-"""
-function test_full_pipeline(msi_data, spectrum_id; output_dir, file_type_prefix, mz_tolerance=0.002)
- println("\n--- Testing Full Preprocessing Pipeline on Spectrum: $spectrum_id (File Type: $file_type_prefix) ---")
+ # --- Load Data ---
+ println("Loading data from: $TEST_FILE")
+ if !isfile(TEST_FILE)
+ println("ERROR: Test file not found. Please check the path in the CONFIG section.")
+ return
+ end
+ msi_data = @time OpenMSIData(TEST_FILE)
+ println("Data loaded successfully. Found $(length(msi_data.spectra_metadata)) spectra.")
- # 1. Determine the spectrum index
- local spec_idx
- if spectrum_id isa Int
- spec_idx = spectrum_id
- else # Tuple for imzML
- spec_idx = msi_data.coordinate_map[spectrum_id...]
+ # --- Determine which spectra to process ---
+ total_spectra = length(msi_data.spectra_metadata)
+ indices_to_process = if NUM_SPECTRA_TO_PROCESS === nothing
+ 1:total_spectra
+ else
+ unique(round.(Int, range(1, total_spectra, length=min(NUM_SPECTRA_TO_PROCESS, total_spectra))))
+ end
+ println("Will process $(length(indices_to_process)) spectra.")
+
+ # --- Data storage for results ---
+ qc_results = DataFrame(spectrum_idx=Int[], metric=String[], value=Float64[], compound=String[])
+ all_processed_peaks_mz = Vector{Vector{Float64}}()
+ all_processed_peaks_int = Vector{Vector{Float64}}()
+ spectrum_indices_with_peaks = Int[]
+
+ # --- Main Loop: Process each spectrum individually ---
+ println("\n" * "="^20 * " Processing Individual Spectra " * "="^20)
+ for (i, idx) in enumerate(indices_to_process)
+ print("\rProcessing spectrum #$idx ($(i)/$(length(indices_to_process)))...")
+
+ process_spectrum(msi_data, idx) do mz, intensity
+ if qc_is_empty(mz, intensity) || !qc_is_regular(mz)
+ # @warn "Skipping empty or irregular spectrum #$idx"
+ return
+ end
+
+ original_mz, original_intensity = copy(mz), copy(intensity)
+ processed_mz, processed_intensity = copy(mz), copy(intensity)
+
+ # --- Step 0: QC (PPM and Resolution) ---
+ # This happens before any modification to the m/z axis
+ let ref_masses = collect(values(INTERNAL_STANDARDS))
+ matched_peaks = find_calibration_peaks(processed_mz, processed_intensity, ref_masses, ppm_tolerance=QC_PARAMS.ppm_tolerance)
+
+ for (compound, theoretical_mz) in INTERNAL_STANDARDS
+ # Find the measured m/z that corresponds to this theoretical_mz
+ measured_mz = 0.0
+ for (theo, meas) in matched_peaks
+ if theo == theoretical_mz
+ measured_mz = meas
+ break
+ end
+ end
+
+ if measured_mz > 0
+ # PPM Error
+ ppm_error = calculate_ppm_error(measured_mz, theoretical_mz)
+ push!(qc_results, (idx, "ppm_error", ppm_error, compound))
+
+ # Resolution
+ resolution = calculate_resolution_fwhm(measured_mz, processed_mz, processed_intensity)
+ if !isnan(resolution)
+ push!(qc_results, (idx, "resolution", resolution, compound))
+ end
+ end
+ end
+ end
+
+ # --- Step 1: Calibration ---
+ if CALIBRATION_PARAMS.enabled
+ ref_masses = collect(values(INTERNAL_STANDARDS))
+ matched_peaks = find_calibration_peaks(processed_mz, processed_intensity, ref_masses, ppm_tolerance=CALIBRATION_PARAMS.ppm_tolerance)
+ if length(matched_peaks) >= 2
+ measured = sort(collect(values(matched_peaks)))
+ theoretical = sort(collect(keys(matched_peaks)))
+ itp = linear_interpolation(measured, theoretical, extrapolation_bc=Line())
+ processed_mz = itp(processed_mz)
+ end
+ end
+
+ # --- Step 2: Smoothing ---
+ if SMOOTHING_PARAMS.enabled
+ processed_intensity = smooth_spectrum(processed_intensity, window=SMOOTHING_PARAMS.window, order=SMOOTHING_PARAMS.order)
+ end
+
+ # --- Step 3: Baseline Correction ---
+ if BASELINE_CORRECTION_PARAMS.enabled
+ baseline = snip_baseline(processed_intensity, iterations=BASELINE_CORRECTION_PARAMS.iterations)
+ processed_intensity .-= baseline
+ processed_intensity = max.(0, processed_intensity) # Ensure non-negativity
+ end
+
+ # --- Step 4: Normalization ---
+ if NORMALIZATION_PARAMS.enabled && NORMALIZATION_PARAMS.method != :none
+ if NORMALIZATION_PARAMS.method == :tic
+ processed_intensity = tic_normalize(processed_intensity)
+ elseif NORMALIZATION_PARAMS.method == :median
+ processed_intensity = median_normalize(processed_intensity)
+ end
+ end
+
+ # --- Step 5: Peak Detection ---
+ local pk_mz, pk_int
+ pk_mz, pk_int = Float64[], Float64[] # Initialize empty
+ if PEAK_DETECTION_PARAMS.enabled
+ if PEAK_DETECTION_PARAMS.method == :profile
+ pk_mz, pk_int = detect_peaks_profile(processed_mz, processed_intensity,
+ half_window=PEAK_DETECTION_PARAMS.half_window,
+ snr_threshold=PEAK_DETECTION_PARAMS.snr_threshold,
+ min_peak_prominence=PEAK_DETECTION_PARAMS.min_peak_prominence,
+ merge_peaks_tolerance=PEAK_DETECTION_PARAMS.merge_peaks_tolerance)
+ elseif PEAK_DETECTION_PARAMS.method == :wavelet
+ pk_mz, pk_int = detect_peaks_wavelet(processed_mz, processed_intensity)
+ else # :centroid
+ pk_mz, pk_int = detect_peaks_centroid(processed_mz, processed_intensity)
+ end
+
+ if !isempty(pk_mz)
+ push!(all_processed_peaks_mz, pk_mz)
+ push!(all_processed_peaks_int, pk_int)
+ push!(spectrum_indices_with_peaks, idx)
+ end
+ end
+
+ # --- Visualization of a sample spectrum ---
+ if i == 1 # Only plot the first processed spectrum
+ println("\nGenerating example plots for spectrum #$idx...")
+ fig = Figure(size=(1200, 800))
+
+ plot_spectrum_on_fig(fig[1,1], original_mz, original_intensity, "1. Original Spectrum")
+ plot_spectrum_on_fig(fig[2,1], processed_mz, processed_intensity, "2. After All Steps (Before Peak Picking)")
+
+ # Plot detected peaks
+ ax = Axis(fig[3,1], title="3. Detected Peaks")
+ if !isempty(pk_mz)
+ stem!(ax, pk_mz, pk_int)
+ end
+
+ save(joinpath(RESULTS_DIR, "example_spectrum_processing.png"), fig)
+ println("Saved example processing plots.")
+ end
+ end
+ end
+ println("\nIndividual spectrum processing complete.")
+
+ # --- Save QC Results ---
+ if !isempty(qc_results)
+ println("\n" * "="^20 * " Generating QC Report " * "="^20)
+ CSV.write(joinpath(RESULTS_DIR, "qc_results.csv"), qc_results)
+ println("QC results saved to qc_results.csv")
+
+ # Generate summary plots for QC
+ fig = Figure(size=(1200, 600))
+
+ # PPM Error Histogram
+ ppm_errors = filter(row -> row.metric == "ppm_error", qc_results).value
+ if !isempty(ppm_errors)
+ ax1 = Axis(fig[1,1], title="PPM Error Distribution", xlabel="PPM Error")
+ hist!(ax1, ppm_errors, bins=30)
+ end
+
+ # Resolution Histogram
+ resolutions = filter(row -> row.metric == "resolution", qc_results).value
+ if !isempty(resolutions)
+ ax2 = Axis(fig[1,2], title="Resolution Distribution", xlabel="Resolution (FWHM)")
+ hist!(ax2, resolutions, bins=30)
+ end
+
+ save(joinpath(RESULTS_DIR, "qc_summary_plots.png"), fig)
+ println("QC summary plots saved.")
end
- if spec_idx == 0
- println("SKIPPED: No spectrum found at coordinates $spectrum_id.")
+ if isempty(all_processed_peaks_mz)
+ @warn "No peaks were detected in any of the processed spectra. Skipping alignment and binning."
return
end
- # 2. Define the pipeline steps in the desired order
- pipeline_steps = [
- :qc,
- :transform,
- :smooth,
- :baseline,
- :normalize,
- :peaks,
- :align, # Align requires multiple spectra, but we'll run it on a single one for now (will warn)
- :bin
- ]
+ # --- Step 6: Peak Alignment ---
+ println("\n" * "="^20 * " Aligning Peaks " * "="^20)
+ aligned_peaks_mz = copy(all_processed_peaks_mz)
+ if PEAK_ALIGNMENT_PARAMS.enabled
+ # Create a reference peak list (e.g., from the spectrum with the most peaks)
+ ref_idx = argmax(length.(all_processed_peaks_mz))
+ reference_peaks = all_processed_peaks_mz[ref_idx]
+ println("Using spectrum $(spectrum_indices_with_peaks[ref_idx]) as alignment reference.")
- # Define parameters for each step
- params = Dict(
- :transform_method => :sqrt,
- :sg_window => 15,
- :sg_order => 2,
- :snip_iterations => 100,
- :normalize_method => :tic,
- :peak_half_window => 10,
- :peak_snr => 3.0,
- :peak_intensity_threshold => 0.0, # For centroid peak detection
- :align_tolerance => mz_tolerance,
- :bin_tolerance => mz_tolerance,
- :bin_min_frequency => 0.0 # Keep all bins for a single spectrum
- )
+ for i in 1:length(aligned_peaks_mz)
+ if i == ref_idx continue end
+ alignment_func = align_peaks_lowess(reference_peaks, aligned_peaks_mz[i],
+ tolerance=PEAK_ALIGNMENT_PARAMS.tolerance,
+ tolerance_unit=PEAK_ALIGNMENT_PARAMS.tolerance_unit,
+ min_matched_peaks=PEAK_ALIGNMENT_PARAMS.min_matched_peaks)
+
+ aligned_peaks_mz[i] = alignment_func(aligned_peaks_mz[i])
+ end
+ println("Peak alignment complete.")
+ end
- # 3. Define the on_stage callback to collect data for overlay plot and save separate plots
- collected_stage_data = []
- stage_counter = Ref(0) # Initialize counter for sequential naming
- normalized_spectrum = nothing # Variable to hold the normalized spectrum
-
- function stage_callback(stage; idx, mz, intensity)
- stage_counter[] += 1 # Increment counter
- println(" -> Generating plot for stage: $stage")
+ # --- Step 7: Peak Binning & Feature Matrix Generation ---
+ println("\n" * "="^20 * " Binning Peaks " * "="^20)
+ if PEAK_BINNING_PARAMS.enabled
+ feature_matrix, mz_bins = bin_peaks(aligned_peaks_mz, all_processed_peaks_int,
+ PEAK_BINNING_PARAMS.tolerance,
+ tolerance_unit=PEAK_BINNING_PARAMS.tolerance_unit,
+ frequency_threshold=PEAK_BINNING_PARAMS.frequency_threshold)
- local fig # Make fig available in the whole function scope
-
- if stage == :normalize
- normalized_spectrum = (mz, intensity)
- fig = plot_stage_spectrum(mz, intensity, title="Stage: $stage (Spectrum $spectrum_id)")
- elseif stage == :peaks && normalized_spectrum !== nothing
- # For the peaks stage, plot the normalized spectrum as a base layer
- fig = Figure(size = (1400, 500))
- ax = Axis(fig[1, 1], title="Stage: Peaks (Spectrum $spectrum_id)", xlabel="m/z", ylabel="Intensity")
- lines!(ax, normalized_spectrum[1], normalized_spectrum[2], color=:gray, label="Normalized Spectrum")
- scatter!(ax, mz, intensity, color=:red, marker=:circle, markersize=8, label="Detected Peaks")
- axislegend(ax)
+ if !isempty(feature_matrix)
+ df = DataFrame(feature_matrix, :auto)
+ rename!(df, ["bin_$(i)" for i in 1:size(df, 2)])
+ insertcols!(df, 1, :spectrum_idx => spectrum_indices_with_peaks)
+
+ CSV.write(joinpath(RESULTS_DIR, "feature_matrix.csv"), df)
+ println("Feature matrix saved to feature_matrix.csv")
+
+ # Save bin m/z ranges
+ bin_df = DataFrame(bin_index=1:length(mz_bins), mz_start=[b[1] for b in mz_bins], mz_end=[b[2] for b in mz_bins])
+ CSV.write(joinpath(RESULTS_DIR, "bin_definitions.csv"), bin_df)
+ println("Bin definitions saved to bin_definitions.csv")
else
- # Default plotting for all other stages
- fig = plot_stage_spectrum(mz, intensity, title="Stage: $stage (Spectrum $spectrum_id)")
+ @warn "Feature matrix was empty after binning."
end
-
- # Save the figure
- stage_output_path = joinpath(output_dir, "$(file_type_prefix)_$(spectrum_id)_$(stage_counter[])_$(stage).png")
- save(stage_output_path, fig)
-
- # Collect data for overlay plot
- push!(collected_stage_data, (stage, mz, intensity))
end
- # 4. Run the pipeline on the single spectrum
- println("Running pipeline with steps: $pipeline_steps")
- processed_result = run_preprocessing_pipeline(
- msi_data,
- [spec_idx], # The pipeline expects a vector of indices
- steps=pipeline_steps,
- params=params,
- on_stage=stage_callback
- )
-
- # 5. Generate and save the overlay plot
- overlay_output_path = joinpath(output_dir, "$(file_type_prefix)_$(spectrum_id)_all_stages_overlay.png")
- plot_overlay_stages(collected_stage_data, overlay_output_path, "Preprocessing Stages Overlay (Spectrum $spectrum_id)")
-
- # 6. Save feature matrix if generated
- if processed_result isa FeatureMatrix
- feature_matrix_output_path = joinpath(output_dir, "$(file_type_prefix)_$(spectrum_id)_feature_matrix.csv")
- # Convert mz_bins to a more readable format for CSV
- mz_labels = ["$(round(b[1], digits=4))_$(round(b[2], digits=4))" for b in processed_result.mz_bins]
- df = DataFrame(processed_result.matrix, Symbol.(mz_labels))
- CSV.write(feature_matrix_output_path, df)
- println("SUCCESS: Feature matrix saved to $feature_matrix_output_path")
- else
- @warn "Pipeline did not return a FeatureMatrix for Spectrum $spectrum_id."
- processed_result
- end
-
- println("--- Pipeline test finished for Spectrum: $spectrum_id (File Type: $file_type_prefix) ---")
- println("Check the '$(output_dir)' directory for output plots and CSVs.")
+ println("\nPreprocessing Test Suite finished successfully!")
end
-"""
- test_full_pipeline_on_total_spectrum(msi_data; output_dir, file_type_prefix)
-
-Tests the full preprocessing pipeline on the *total spectrum* (sum of all spectra)
-and saves a plot for each intermediate step.
-"""
-function test_full_pipeline_on_total_spectrum(msi_data; output_dir, file_type_prefix, mz_tolerance=0.002)
- println("\n--- Testing Full Preprocessing Pipeline on TOTAL Spectrum (File Type: $file_type_prefix) ---")
-
- # 1. Get the total spectrum
- total_mz, total_intensity = get_total_spectrum(msi_data)
- total_spectrum = (total_mz, total_intensity)
-
- if qc_is_empty(total_mz, total_intensity)
- println("SKIPPED: Total spectrum is empty.")
- return
- end
-
- # 2. Define the pipeline steps and parameters (same as for single spectrum)
- pipeline_steps = [
- :qc,
- :transform,
- :smooth,
- :baseline,
- :normalize,
- :peaks,
- :align, # Align requires multiple spectra, but we'll run it on a single one for now (will warn)
- :bin
- ]
-
- params = Dict(
- :transform_method => :sqrt,
- :sg_window => 15,
- :sg_order => 2,
- :snip_iterations => 100,
- :normalize_method => :tic,
- :peak_half_window => 10,
- :peak_snr => 3.0,
- :align_tolerance => mz_tolerance,
- :bin_tolerance => mz_tolerance,
- :bin_min_frequency => 0.0 # Keep all bins for a single spectrum
- )
-
- # 3. Define the on_stage callback
- collected_stage_data = []
- stage_counter = Ref(0) # Initialize counter for sequential naming
- normalized_spectrum_total = nothing # Variable to hold the normalized spectrum
-
- function stage_callback_total(stage; idx, mz, intensity)
- stage_counter[] += 1 # Increment counter
- println(" -> Generating plot for stage: $stage (Total Spectrum)")
-
- local fig
-
- if stage == :normalize
- normalized_spectrum_total = (mz, intensity)
- fig = plot_stage_spectrum(mz, intensity, title="Stage: $stage (Total Spectrum)")
- elseif stage == :peaks && normalized_spectrum_total !== nothing
- fig = Figure(size = (1400, 500))
- ax = Axis(fig[1, 1], title="Stage: Peaks (Total Spectrum)", xlabel="m/z", ylabel="Intensity")
- lines!(ax, normalized_spectrum_total[1], normalized_spectrum_total[2], color=:gray, label="Normalized Spectrum")
- scatter!(ax, mz, intensity, color=:red, marker=:circle, markersize=8, label="Detected Peaks")
- axislegend(ax)
- else
- fig = plot_stage_spectrum(mz, intensity, title="Stage: $stage (Total Spectrum)")
- end
-
- # Save separate plot with sequential name
- stage_output_path = joinpath(output_dir, "$(file_type_prefix)_total_$(stage_counter[])_$(stage).png")
- save(stage_output_path, fig)
-
- # Collect data for overlay plot
- push!(collected_stage_data, (stage, mz, intensity))
- end
-
- # 4. Run the pipeline on the single total spectrum
- println("Running pipeline with steps: $pipeline_steps")
- processed_result = run_preprocessing_pipeline(
- [total_spectrum], # Pass the total spectrum as a vector of one spectrum
- steps=pipeline_steps,
- params=params,
- on_stage=stage_callback_total
- )
-
- # 5. Generate and save the overlay plot
- overlay_output_path = joinpath(output_dir, "$(file_type_prefix)_total_all_stages_overlay.png")
- plot_overlay_stages(collected_stage_data, overlay_output_path, "Preprocessing Stages Overlay (Total Spectrum)")
-
- # 6. Save feature matrix if generated
- if processed_result isa FeatureMatrix
- feature_matrix_output_path = joinpath(output_dir, "$(file_type_prefix)_total_feature_matrix.csv")
- # Convert mz_bins to a more readable format for CSV
- mz_labels = ["$(round(b[1], digits=4))_$(round(b[2], digits=4))" for b in processed_result.mz_bins]
- df = DataFrame(processed_result.matrix, Symbol.(mz_labels))
- CSV.write(feature_matrix_output_path, df)
- println("SUCCESS: Feature matrix saved to $feature_matrix_output_path")
- else
- @warn "Pipeline did not return a FeatureMatrix for Total Spectrum."
- processed_result
- end
-
- println("--- Pipeline test finished for TOTAL Spectrum (File Type: $file_type_prefix) ---")
- println("Check the '$(output_dir)' directory for output plots and CSVs.")
-end
-
-
-# ===================================================================
-# TEST RUNNER
-# ===================================================================
-
-function run_preprocessing_tests()
- println("="^80)
- println("STARTING PREPROCESSING TEST SUITE")
- println("="^80)
-
- # --- Test Case 1: Run full pipeline on a single mzML spectrum ---
- println("\n" * "="^20 * " Test Case 1: Full Pipeline on .mzML Spectrum " * "="^20)
- println("FILE: ", TEST_MZML_FILE)
- if isfile(TEST_MZML_FILE)
- try
- msi_data_mzml = OpenMSIData(TEST_MZML_FILE)
-
- # Dynamically determine tolerance
- println("\n--- Calculating optimal tolerance for .mzML data ---")
- report_mzml = analyze_mass_accuracy(msi_data_mzml, get_common_calibration_standards(:maldi_pos))
- mz_tolerance_mzml = 0.002 # Default
- if haskey(report_mzml, :optimal_ppm) && !isnan(report_mzml.optimal_ppm) && !isempty(report_mzml.matched_peaks)
- avg_mz = mean([p[1] for p in report_mzml.matched_peaks])
- mz_tolerance_mzml = avg_mz * report_mzml.optimal_ppm / 1e6
- println("Optimal PPM: $(round(report_mzml.optimal_ppm, digits=2)), Average m/z: $(round(avg_mz, digits=2))")
- println("Calculated m/z tolerance: $(round(mz_tolerance_mzml, digits=5))")
- else
- println("Could not determine optimal tolerance, using default: $mz_tolerance_mzml")
- end
-
- # Create a dedicated subdirectory for the output plots
- mzml_output_dir = joinpath(RESULTS_DIR, "mzml_pipeline_stages")
- mkpath(mzml_output_dir)
-
- test_full_pipeline(msi_data_mzml, MZML_SPECTRUM_ID, output_dir=mzml_output_dir, file_type_prefix="mzml", mz_tolerance=mz_tolerance_mzml)
- test_full_pipeline_on_total_spectrum(msi_data_mzml, output_dir=mzml_output_dir, file_type_prefix="mzml", mz_tolerance=mz_tolerance_mzml)
- catch e
- println("ERROR in .mzML pipeline test: $e")
- showerror(stdout, e, catch_backtrace())
- end
- else
- println("SKIPPED: File not found: $TEST_MZML_FILE")
- end
-
- # --- Test Case 2: Run full pipeline on a single imzML spectrum ---
- println("\n" * "="^20 * " Test Case 2: Full Pipeline on .imzML Spectrum " * "="^20)
- println("FILE: ", TEST_IMZML_FILE)
- if isfile(TEST_IMZML_FILE)
- try
- msi_data_imzml = OpenMSIData(TEST_IMZML_FILE)
-
- # Dynamically determine tolerance
- println("\n--- Calculating optimal tolerance for .imzML data ---")
- report_imzml = analyze_mass_accuracy(msi_data_imzml, get_common_calibration_standards(:maldi_pos))
- mz_tolerance_imzml = 0.002 # Default
- if haskey(report_imzml, :optimal_ppm) && !isnan(report_imzml.optimal_ppm) && !isempty(report_imzml.matched_peaks)
- avg_mz = mean([p[1] for p in report_imzml.matched_peaks])
- mz_tolerance_imzml = avg_mz * report_imzml.optimal_ppm / 1e6
- println("Optimal PPM: $(round(report_imzml.optimal_ppm, digits=2)), Average m/z: $(round(avg_mz, digits=2))")
- println("Calculated m/z tolerance: $(round(mz_tolerance_imzml, digits=5))")
- else
- println("Could not determine optimal tolerance, using default: $mz_tolerance_imzml")
- end
-
- # Create a dedicated subdirectory for the output plots
- imzml_output_dir = joinpath(RESULTS_DIR, "imzml_pipeline_stages")
- mkpath(imzml_output_dir)
-
- test_full_pipeline(msi_data_imzml, IMZML_COORDS, output_dir=imzml_output_dir, file_type_prefix="imzml", mz_tolerance=mz_tolerance_imzml)
- test_full_pipeline_on_total_spectrum(msi_data_imzml, output_dir=imzml_output_dir, file_type_prefix="imzml", mz_tolerance=mz_tolerance_imzml)
- # generate_qc_report(msi_data_imzml, TEST_IMZML_FILE, output_dir=imzml_output_dir)
- custom_reference_peaks = Dict(
- 31.974 => "Red Phosphorus",
- 432.6584 => "P13",
- 464.6059 => "P15",
- 526.5534 => "P17",
- 650.4485 => "P21",
- 774.3435 => "P25",
- 898.2385 => "P29",
- 950.1861 => "P31",
- 1022.1336 => "P33",
- 1146.0286 => "P37",
- 1593.8187 => "P45",
- 772.433 => "Unknown 1",
- 772.5253 => "Unknown 2"
- )
- n_samples = length(msi_data_imzml.spectra_metadata)
- generate_qc_report(msi_data_imzml, TEST_IMZML_FILE, reference_peaks=custom_reference_peaks, output_dir=imzml_output_dir, sample_spectra=n_samples)
- catch e
- println("ERROR in .imzML pipeline test: $e")
- showerror(stdout, e, catch_backtrace())
- end
- else
- println("SKIPPED: File not found: $TEST_IMZML_FILE")
- end
-
- # --- Test Case 3: Generate QC Report for .imzML data ---
- println("\n" * "="^20 * " Test Case 3: QC Report Generation for .imzML " * "="^20)
- println("FILE: ", TEST_IMZML_FILE)
- if isfile(TEST_IMZML_FILE)
- try
- msi_data_imzml = OpenMSIData(TEST_IMZML_FILE)
-
- # Create a dedicated subdirectory for the QC report
- qc_output_dir = joinpath(RESULTS_DIR, "qc_report")
- mkpath(qc_output_dir)
-
- println("\n--- Generating comprehensive QC report ---")
- custom_reference_peaks = Dict(
- 31.974 => "Red Phosphorus",
- 432.6584 => "P13",
- 464.6059 => "P15",
- 526.5534 => "P17",
- 650.4485 => "P21",
- 774.3435 => "P25",
- 898.2385 => "P29",
- 950.1861 => "P31",
- 1022.1336 => "P33",
- 1146.0286 => "P37",
- 1593.8187 => "P45",
- 772.433 => "Unknown 1",
- 772.5253 => "Unknown 2"
- )
- # You can control the number of spectra sampled for the QC report.
- # For the most accurate results, you can sample all spectra, but it will take longer.
- # To sample all, use: n_samples = length(msi_data_imzml.spectra_metadata)
- n_samples = length(msi_data_imzml.spectra_metadata)
-
- #generate_qc_report(msi_data_imzml, TEST_IMZML_FILE, output_dir=qc_output_dir)
- generate_qc_report(msi_data_imzml, TEST_IMZML_FILE, reference_peaks=custom_reference_peaks, output_dir=qc_output_dir, sample_spectra=n_samples)
-
- println("\n--- Analyzing specific reference peaks ---")
- #=
- reference_peaks = Dict(
- 104.10754 => "Imidazole",
- 175.11995 => "GPC fragment",
- 226.15687 => "Phosphocholine"
- )
- =#
- report = analyze_mass_accuracy(msi_data_imzml, custom_reference_peaks)
- if haskey(report, :optimal_ppm)
- println("Optimal PPM tolerance with specific peaks: $(round(report.optimal_ppm, digits=2)) ppm")
- else
- println("Could not determine optimal PPM with specific peaks.")
- end
-
- catch e
- println("ERROR in QC report generation test: $e")
- showerror(stdout, e, catch_backtrace())
- end
- else
- println("SKIPPED: File not found: $TEST_IMZML_FILE")
- end
-
- println("\nPreprocessing tests finished.")
-end
# --- Execute ---
-# Ensure the results directory exists
-mkpath(RESULTS_DIR)
-@time run_preprocessing_tests()
\ No newline at end of file
+run_preprocessing_suite()