most stable version of JuliaMSI
This commit is contained in:
parent
0164e54d77
commit
977b416579
@ -37,7 +37,6 @@ NativeFileDialog = "e1fe445b-aa65-4df4-81c1-2041507f0fd4"
|
|||||||
NaturalSort = "c020b1a1-e9b0-503a-9c33-f039bfc54a85"
|
NaturalSort = "c020b1a1-e9b0-503a-9c33-f039bfc54a85"
|
||||||
Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
|
Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
|
||||||
PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5"
|
PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5"
|
||||||
Polyester = "f401cd3a-86c5-430c-a3c3-63b7194639e4"
|
|
||||||
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
|
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
|
||||||
ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
|
ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
|
||||||
SavitzkyGolay = "c4bf5708-b6a6-4fbe-bcd0-6850ed671584"
|
SavitzkyGolay = "c4bf5708-b6a6-4fbe-bcd0-6850ed671584"
|
||||||
|
|||||||
@ -7,43 +7,9 @@ efficiently.
|
|||||||
"""
|
"""
|
||||||
|
|
||||||
using Base64, Libz, Serialization, Printf, DataFrames, Base.Threads, StatsBase
|
using Base64, Libz, Serialization, Printf, DataFrames, Base.Threads, StatsBase
|
||||||
using .Platform # Import the new Platform module
|
|
||||||
|
|
||||||
const FILE_HANDLE_LOCK = ReentrantLock()
|
const FILE_HANDLE_LOCK = ReentrantLock()
|
||||||
|
|
||||||
"""
|
|
||||||
calculate_optimal_analytics_chunk_size(total_spectra::Int, system_caps::NamedTuple) -> Int
|
|
||||||
|
|
||||||
Calculates an optimal chunk size for analytics processing based on total number of spectra
|
|
||||||
and detected system capabilities (total memory, number of CPU threads).
|
|
||||||
|
|
||||||
The heuristic aims to balance memory usage and parallelism:
|
|
||||||
- Smaller chunks for less memory or fewer threads.
|
|
||||||
- Larger chunks for more memory or more threads, but with a practical upper limit.
|
|
||||||
"""
|
|
||||||
function calculate_optimal_analytics_chunk_size(total_spectra::Int, system_caps::NamedTuple)::Int
|
|
||||||
# Base chunk size - can be refined with more profiling
|
|
||||||
base_chunk_size = 5000
|
|
||||||
|
|
||||||
# Adjust based on available RAM (heuristic: more RAM, larger chunks)
|
|
||||||
# Assume 1GB RAM can handle a certain chunk size comfortably
|
|
||||||
# E.g., for 32GB RAM, allow 2x more than for 16GB
|
|
||||||
memory_factor = max(1.0, system_caps.total_memory_gb / 16.0) # Scale based on 16GB baseline
|
|
||||||
|
|
||||||
# Adjust based on CPU threads (heuristic: more threads, potentially more, smaller chunks to keep threads busy, or fewer large chunks)
|
|
||||||
# A simple approach for now: scale with sqrt of threads to avoid overly large chunks with many threads
|
|
||||||
thread_factor = max(1.0, sqrt(system_caps.num_cpu_threads / 4.0)) # Scale based on 4-thread baseline
|
|
||||||
|
|
||||||
# Combine factors, with a practical minimum and maximum
|
|
||||||
dynamic_chunk_size = round(Int, base_chunk_size * memory_factor * thread_factor)
|
|
||||||
|
|
||||||
# Ensure chunk size is within reasonable bounds
|
|
||||||
min_chunk_size = 500
|
|
||||||
max_chunk_size = 20000
|
|
||||||
|
|
||||||
return clamp(dynamic_chunk_size, min_chunk_size, max_chunk_size)
|
|
||||||
end
|
|
||||||
|
|
||||||
"""
|
"""
|
||||||
ThreadSafeFileHandle
|
ThreadSafeFileHandle
|
||||||
|
|
||||||
|
|||||||
@ -1,29 +0,0 @@
|
|||||||
# src/Platform.jl
|
|
||||||
|
|
||||||
module Platform
|
|
||||||
|
|
||||||
using Base.Threads # For Threads.nthreads()
|
|
||||||
|
|
||||||
# Export functions for system capability detection
|
|
||||||
export detect_system_capabilities
|
|
||||||
|
|
||||||
"""
|
|
||||||
detect_system_capabilities()
|
|
||||||
|
|
||||||
Detects and returns key system capabilities such as total memory and number of CPU threads.
|
|
||||||
|
|
||||||
# Returns
|
|
||||||
- A `NamedTuple` with fields:
|
|
||||||
- `total_memory_gb::Float64`: Total physical memory in gigabytes.
|
|
||||||
- `num_cpu_threads::Int`: Number of available CPU threads.
|
|
||||||
"""
|
|
||||||
function detect_system_capabilities()::NamedTuple
|
|
||||||
total_memory_bytes = Sys.total_memory()
|
|
||||||
total_memory_gb = total_memory_bytes / (1024^3) # Convert bytes to gigabytes
|
|
||||||
|
|
||||||
num_cpu_threads = Sys.CPU_THREADS # Total logical CPU cores
|
|
||||||
|
|
||||||
return (total_memory_gb=total_memory_gb, num_cpu_threads=num_cpu_threads)
|
|
||||||
end
|
|
||||||
|
|
||||||
end # module Platform
|
|
||||||
102
src/imzML.jl
102
src/imzML.jl
@ -721,65 +721,6 @@ function find_mass(mz_array::AbstractVector{<:Real}, intensity_array::AbstractVe
|
|||||||
return max_intensity
|
return max_intensity
|
||||||
end
|
end
|
||||||
|
|
||||||
#=
|
|
||||||
"""
|
|
||||||
load_slices(folder, masses, tolerance)
|
|
||||||
|
|
||||||
Loads image slices for multiple masses from all `.imzML` files in a directory.
|
|
||||||
This function is now refactored to use the new MSIData architecture and its
|
|
||||||
caching capabilities.
|
|
||||||
"""
|
|
||||||
function load_slices(folder::String, masses::AbstractVector{<:Real}, tolerance::Real)
|
|
||||||
files = filter(f -> endswith(f, ".imzML"), readdir(folder, join=true))
|
|
||||||
if isempty(files)
|
|
||||||
@warn "No .imzML files found in the specified directory: $folder"
|
|
||||||
return (Array{Matrix{Float64}}(undef, 0, 0), String[])
|
|
||||||
end
|
|
||||||
n_files = length(files)
|
|
||||||
n_slices = length(masses)
|
|
||||||
|
|
||||||
# FIXED: Use concrete type instead of Any
|
|
||||||
img_list = Array{Matrix{Float64}}(undef, n_files, n_slices)
|
|
||||||
names = String[]
|
|
||||||
|
|
||||||
for (i, file) in enumerate(files)
|
|
||||||
name = basename(file)
|
|
||||||
push!(names, name)
|
|
||||||
@info "Processing $(i)/$(n_files): $(name)"
|
|
||||||
|
|
||||||
# Load data using the new lazy loader, returning an MSIData object
|
|
||||||
msi_data = @time load_imzml_lazy(file)
|
|
||||||
|
|
||||||
# Create empty images for all slices for the current file
|
|
||||||
width, height = msi_data.image_dims
|
|
||||||
current_file_slices = [zeros(Float64, width, height) for _ in 1:n_slices]
|
|
||||||
|
|
||||||
# Use the high-performance iterator to process all spectra
|
|
||||||
_iterate_spectra_fast(msi_data) do spec_idx, mz_array, intensity_array
|
|
||||||
meta = msi_data.spectra_metadata[spec_idx]
|
|
||||||
|
|
||||||
# Now, check for all masses of interest in this single spectrum
|
|
||||||
for (j, mass) in enumerate(masses)
|
|
||||||
intensity = find_mass(mz_array, intensity_array, mass, tolerance)
|
|
||||||
if intensity > 0.0
|
|
||||||
if 1 <= meta.x <= width && 1 <= meta.y <= height
|
|
||||||
current_file_slices[j][meta.y, meta.x] = intensity
|
|
||||||
end
|
|
||||||
end
|
|
||||||
end
|
|
||||||
end # end of fast iterator
|
|
||||||
|
|
||||||
# Assign the generated images to the main list
|
|
||||||
for j in 1:n_slices
|
|
||||||
img_list[i, j] = current_file_slices[j]
|
|
||||||
end
|
|
||||||
|
|
||||||
end # end of files loop
|
|
||||||
|
|
||||||
return (img_list, names)
|
|
||||||
end
|
|
||||||
=#
|
|
||||||
|
|
||||||
"""
|
"""
|
||||||
get_mz_slice(data::MSIData, mass::Real, tolerance::Real)
|
get_mz_slice(data::MSIData, mass::Real, tolerance::Real)
|
||||||
|
|
||||||
@ -1294,49 +1235,6 @@ end
|
|||||||
# Analysis and Visualization
|
# Analysis and Visualization
|
||||||
#
|
#
|
||||||
# ============================================================================
|
# ============================================================================
|
||||||
#=
|
|
||||||
"""
|
|
||||||
display_statistics(slices::AbstractArray{<:AbstractMatrix{<:Real}, 2},
|
|
||||||
names::AbstractVector{String}, masses::AbstractVector{<:Real})
|
|
||||||
|
|
||||||
Calculates and prints key statistics for each slice.
|
|
||||||
"""
|
|
||||||
function display_statistics(slices::AbstractArray{<:AbstractMatrix{<:Real}, 2},
|
|
||||||
names::AbstractVector{String}, masses::AbstractVector{<:Real})
|
|
||||||
if isempty(slices)
|
|
||||||
@warn "Cannot display statistics for empty slice list."
|
|
||||||
return nothing
|
|
||||||
end
|
|
||||||
|
|
||||||
n_files, n_masses = size(slices)
|
|
||||||
stats_to_calc = Dict{String, Function}("Mean" => mean)
|
|
||||||
all_dfs = Dict{String, DataFrame}()
|
|
||||||
|
|
||||||
for (stat_name, stat_func) in stats_to_calc
|
|
||||||
# Pre-allocate matrix for better performance
|
|
||||||
stat_matrix = zeros(Float64, n_files, n_masses)
|
|
||||||
|
|
||||||
@inbounds for i in 1:n_files, j in 1:n_masses
|
|
||||||
flat_slice = vec(slices[i, j])
|
|
||||||
if !isempty(flat_slice)
|
|
||||||
stat_matrix[i, j] = stat_func(flat_slice)
|
|
||||||
end
|
|
||||||
end
|
|
||||||
|
|
||||||
# Create the DataFrame with masses as column headers
|
|
||||||
df = DataFrame(stat_matrix, Symbol.(masses))
|
|
||||||
# Insert the file names as the first column
|
|
||||||
insertcols!(df, 1, :Data => names)
|
|
||||||
mz_row = ["m/z"; masses...]
|
|
||||||
push!(df, mz_row)
|
|
||||||
|
|
||||||
println("\n--- Statistics: $(stat_name) ---")
|
|
||||||
println(df)
|
|
||||||
all_dfs[stat_name] = df
|
|
||||||
end
|
|
||||||
return all_dfs
|
|
||||||
end
|
|
||||||
=#
|
|
||||||
|
|
||||||
"""
|
"""
|
||||||
plot_slices(slices, names, masses, output_dir; stage_name, bins=256, dpi=150, global_bounds=nothing)
|
plot_slices(slices, names, masses, output_dir; stage_name, bins=256, dpi=150, global_bounds=nothing)
|
||||||
|
|||||||
124
src/mzML.jl
124
src/mzML.jl
@ -460,127 +460,3 @@ function load_mzml_lazy(file_path::String; cache_size::Int=100)
|
|||||||
rethrow(e)
|
rethrow(e)
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
||||||
#=
|
|
||||||
"""
|
|
||||||
LoadMzml(fileName::String)
|
|
||||||
|
|
||||||
Eagerly loads all spectra from a .mzML file into memory.
|
|
||||||
This function now uses the new lazy-loading
|
|
||||||
architecture internally but presents the data in the old format.
|
|
||||||
|
|
||||||
# Arguments
|
|
||||||
- `fileName`: The path to the `.mzML` file.
|
|
||||||
|
|
||||||
# Returns
|
|
||||||
- A `2xN` matrix where `N` is the number of spectra. The first row contains
|
|
||||||
m/z arrays and the second row contains intensity arrays.
|
|
||||||
"""
|
|
||||||
function LoadMzml(fileName::String)
|
|
||||||
# Use the lazy loader to get the MSIData object
|
|
||||||
msi_data = load_mzml_lazy(fileName, cache_size=0) # No need to cache if we load all
|
|
||||||
|
|
||||||
num_spectra = length(msi_data.spectra_metadata)
|
|
||||||
|
|
||||||
# FIXED: Use concrete typed arrays instead of Array{Any}
|
|
||||||
spectra_matrix = Vector{Tuple{Vector{Float64}, Vector{Float64}}}(undef, num_spectra)
|
|
||||||
|
|
||||||
# Pre-allocate and use bounds checking optimization
|
|
||||||
@inbounds for i in 1:num_spectra
|
|
||||||
# Use the new GetSpectrum API
|
|
||||||
mz, intensity = GetSpectrum(msi_data, i)
|
|
||||||
spectra_matrix[i] = (mz, intensity)
|
|
||||||
end
|
|
||||||
|
|
||||||
# Convert to the expected 2xN format if needed by downstream code
|
|
||||||
# Note: This maintains the original interface but with better typing
|
|
||||||
result_matrix = Array{Any}(undef, (2, num_spectra))
|
|
||||||
@inbounds for i in 1:num_spectra
|
|
||||||
result_matrix[1, i] = spectra_matrix[i][1]
|
|
||||||
result_matrix[2, i] = spectra_matrix[i][2]
|
|
||||||
end
|
|
||||||
|
|
||||||
return result_matrix
|
|
||||||
end
|
|
||||||
=#
|
|
||||||
|
|
||||||
#=
|
|
||||||
"""
|
|
||||||
load_mzml_batch(file_path::String, spectrum_indices::AbstractVector{Int})
|
|
||||||
|
|
||||||
Loads a specific batch of spectra from an mzML file efficiently.
|
|
||||||
Useful for parallel processing or when only specific spectra are needed.
|
|
||||||
|
|
||||||
# Arguments
|
|
||||||
- `file_path`: Path to the mzML file
|
|
||||||
- `spectrum_indices`: Indices of spectra to load
|
|
||||||
|
|
||||||
# Returns
|
|
||||||
- Vector of (mz_array, intensity_array) tuples
|
|
||||||
"""
|
|
||||||
function load_mzml_batch(file_path::String, spectrum_indices::AbstractVector{Int})
|
|
||||||
msi_data = load_mzml_lazy(file_path)
|
|
||||||
num_to_load = length(spectrum_indices)
|
|
||||||
|
|
||||||
# Pre-allocate result with concrete types
|
|
||||||
results = Vector{Tuple{Vector{Float64}, Vector{Float64}}}(undef, num_to_load)
|
|
||||||
|
|
||||||
@inbounds for (i, idx) in enumerate(spectrum_indices)
|
|
||||||
mz, intensity = GetSpectrum(msi_data, idx)
|
|
||||||
results[i] = (mz, intensity)
|
|
||||||
end
|
|
||||||
|
|
||||||
return results
|
|
||||||
end
|
|
||||||
=#
|
|
||||||
|
|
||||||
#=
|
|
||||||
"""
|
|
||||||
get_mzml_summary(file_path::String)::NamedTuple
|
|
||||||
|
|
||||||
Quickly extracts summary information from an mzML file without loading all metadata.
|
|
||||||
|
|
||||||
# Returns
|
|
||||||
- Named tuple with: num_spectra, mz_range, intensity_range, data_formats
|
|
||||||
"""
|
|
||||||
function get_mzml_summary(file_path::String)::NamedTuple
|
|
||||||
ts_stream = ThreadSafeFileHandle(file_path, "r")
|
|
||||||
try
|
|
||||||
index_offset = find_index_offset(ts_stream.handle)
|
|
||||||
seek(ts_stream.handle, index_offset)
|
|
||||||
|
|
||||||
if find_tag(ts_stream.handle, r"<index\s+name=\"spectrum\"") === nothing
|
|
||||||
throw(FileFormatError("Could not find spectrum index."))
|
|
||||||
end
|
|
||||||
|
|
||||||
spectrum_offsets = parse_offset_list(ts_stream.handle)
|
|
||||||
num_spectra = length(spectrum_offsets)
|
|
||||||
|
|
||||||
# Sample first few spectra to determine formats and ranges
|
|
||||||
sample_size = min(10, num_spectra)
|
|
||||||
sample_metadata = [parse_spectrum_metadata(ts_stream.handle, spectrum_offsets[i]) for i in 1:sample_size]
|
|
||||||
|
|
||||||
# Extract formats from sample
|
|
||||||
mz_format = sample_metadata[1].mz_asset.format
|
|
||||||
intensity_format = sample_metadata[1].int_asset.format
|
|
||||||
|
|
||||||
return (num_spectra=num_spectra,
|
|
||||||
mz_format=mz_format,
|
|
||||||
intensity_format=intensity_format,
|
|
||||||
sample_spectra=sample_size)
|
|
||||||
finally
|
|
||||||
close(ts_stream)
|
|
||||||
end
|
|
||||||
end
|
|
||||||
=#
|
|
||||||
|
|
||||||
#=
|
|
||||||
# Performance optimization: Cache frequently used regex patterns
|
|
||||||
const PRE_COMPILED_REGEX = (
|
|
||||||
encoded_length = r"<binaryDataArray\s+encodedLength=\"(\d+)\"",
|
|
||||||
spectrum_id = r"<spectrum\s+index=\"\d+\"\s+id=\"([^\"]+)",
|
|
||||||
offset = r"<offset[^>]*>(\d+)</offset>",
|
|
||||||
index_list = r"<index\s+name=\"spectrum\"",
|
|
||||||
index_offset = r"<indexListOffset>(\d+)</indexListOffset>"
|
|
||||||
)
|
|
||||||
=#
|
|
||||||
|
|||||||
Loading…
x
Reference in New Issue
Block a user