From 4e907dee7caea4d409caa07a9cfb25fd14ee4cf6 Mon Sep 17 00:00:00 2001 From: meibujun Date: Fri, 21 Nov 2025 17:56:02 +0800 Subject: [PATCH] Implement GenomicPro multi-omics pipeline --- Project.toml | 15 ++++ config/example_pipeline.json | 42 +++++++++++ readme.md | 80 ++++++++++++++++++++- scripts/genomicpro_cli.jl | 18 +++++ src/GenomicPro.jl | 27 +++++++ src/analysis.jl | 108 ++++++++++++++++++++++++++++ src/integration.jl | 44 ++++++++++++ src/io.jl | 111 +++++++++++++++++++++++++++++ src/pipeline.jl | 17 +++++ src/preprocessing.jl | 126 +++++++++++++++++++++++++++++++++ src/reporting.jl | 102 ++++++++++++++++++++++++++ src/types.jl | 96 +++++++++++++++++++++++++ src/utils.jl | 119 +++++++++++++++++++++++++++++++ src/validation.jl | 38 ++++++++++ test/config/test_pipeline.json | 43 +++++++++++ test/data/proteomics.csv | 5 ++ test/data/transcriptomics.csv | 5 ++ test/runtests.jl | 38 ++++++++++ 18 files changed, 1033 insertions(+), 1 deletion(-) create mode 100644 Project.toml create mode 100644 config/example_pipeline.json create mode 100755 scripts/genomicpro_cli.jl create mode 100644 src/GenomicPro.jl create mode 100644 src/analysis.jl create mode 100644 src/integration.jl create mode 100644 src/io.jl create mode 100644 src/pipeline.jl create mode 100644 src/preprocessing.jl create mode 100644 src/reporting.jl create mode 100644 src/types.jl create mode 100644 src/utils.jl create mode 100644 src/validation.jl create mode 100644 test/config/test_pipeline.json create mode 100644 test/data/proteomics.csv create mode 100644 test/data/transcriptomics.csv create mode 100644 test/runtests.jl diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..dc3ff8a --- /dev/null +++ b/Project.toml @@ -0,0 +1,15 @@ +name = "GenomicPro" +uuid = "0c0bde55-41a4-47af-8c3e-91f915c1c200" +authors = ["AI Developer "] +version = "0.1.0" + +[deps] +CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +JSON3 = "0f8b85d8-493c-5b1a-9879-9ecda8c0e6a9" + +[compat] +CSV = "0.10" +DataFrames = "1" +JSON3 = "1" +julia = "1.12" diff --git a/config/example_pipeline.json b/config/example_pipeline.json new file mode 100644 index 0000000..03b5b7b --- /dev/null +++ b/config/example_pipeline.json @@ -0,0 +1,42 @@ +{ + "datasets": [ + { + "name": "transcriptomics", + "path": "data/transcriptomics.csv", + "delimiter": ",", + "feature_column": "Feature", + "normalization": "zscore", + "imputation": "mean", + "metadata": { + "missing_strings": ["NA", ""] + } + }, + { + "name": "proteomics", + "path": "data/proteomics.csv", + "delimiter": ",", + "feature_column": "Feature", + "normalization": "zscore", + "imputation": "median" + } + ], + "integration": { + "strategy": "concatenate" + }, + "analysis": { + "run_pca": true, + "pca_components": 3, + "clustering": "kmeans", + "cluster_count": 3 + }, + "report": { + "output_path": "reports/pipeline_report.json", + "text_output_path": "reports/pipeline_report.txt" + }, + "quality_control": { + "min_samples": 2, + "min_features": 2, + "min_nonmissing_ratio": 0.5, + "min_variance": 1e-8 + } +} diff --git a/readme.md b/readme.md index e173354..e6ae164 100644 --- a/readme.md +++ b/readme.md @@ -1 +1,79 @@ -meibujun julia +# GenomicPro + +GenomicPro 是一个使用 Julia 语言实现的多组学（multi-omics）分析流水线框架，提供从数据加载、质量控制、归一化、组学数据对齐、集成分析到报告生成的一站式能力。该实现基于 Julia v1.12.1，侧重于易用性与可扩展性，支持 CLI 运行以及通过配置文件管理复杂工作流。 + +## 关键特性 + +- **多组学数据加载**：支持从 CSV/TSV 表格读取不同组学层的数据，灵活配置分隔符、特征列、缺失值标记等。 +- **缺失值填补与质量控制**：内置均值/中位数/零值填补策略，支持对低方差、缺失率过高的特征进行过滤。 +- **标准化与变换**：提供 z-score、min-max、robust scaling 以及可选对数变换，适配不同测序平台的特性。 +- **数据对齐与集成**：自动对齐样本，支持特征拼接（concatenate）与加权求和（weighted sum）两种集成策略。 +- **统计/机器学习分析**：内置 PCA 降维与 KMeans 聚类，可根据配置决定是否执行并调整参数。 +- **报告生成**：输出 JSON 与文本摘要，涵盖数据概览、质量指标、分析结果，方便集成到下游系统。 + +## 项目结构 + +``` +Project.toml +src/ + GenomicPro.jl # 模块入口 + types.jl # 核心类型定义 + io.jl # 数据读取与配置解析 + preprocessing.jl # 缺失值填补、归一化、对齐等预处理 + integration.jl # 多组学集成算法 + analysis.jl # PCA、KMeans 等分析 + reporting.jl # 报告与摘要生成 + validation.jl # 数据校验与质量过滤 + utils.jl # 通用辅助函数 + pipeline.jl # 流水线编排 +scripts/ + genomicpro_cli.jl # 命令行入口 +config/ + example_pipeline.json# 配置示例 +``` + +## 安装依赖 + +```bash +julia --project -e 'using Pkg; Pkg.instantiate()' +``` + +## 运行示例流水线 + +准备示例数据（可参考 `test/data` 目录中的 CSV），并在根目录创建输出目录： + +```bash +mkdir -p data reports +cp test/data/*.csv data/ +``` + +然后运行命令行脚本： + +```bash +julia scripts/genomicpro_cli.jl config/example_pipeline.json +``` + +执行完成后，JSON 与文本报告分别保存在 `reports/pipeline_report.json` 与 `reports/pipeline_report.txt`。 + +## 在代码中使用 + +```julia +using GenomicPro + +config = read_pipeline_config("config/example_pipeline.json") +result = run_pipeline(config) +println(render_summary(result)) +``` + +## 测试 + +```bash +julia --project -e 'using Pkg; Pkg.test()' +``` + +## 扩展建议 + +- 接入更多集成策略（如基于网络的融合、协同矩阵分解等）。 +- 支持批量效应校正与更丰富的归一化方法。 +- 引入更高级的聚类与分类算法，或与 MLJ 生态衔接。 +- 构建交互式可视化报告，支持在浏览器中探索多组学结果。 diff --git a/scripts/genomicpro_cli.jl b/scripts/genomicpro_cli.jl new file mode 100755 index 0000000..993132c --- /dev/null +++ b/scripts/genomicpro_cli.jl @@ -0,0 +1,18 @@ +#!/usr/bin/env julia +using Pkg +Pkg.activate(joinpath(@__DIR__, "..")) +using GenomicPro + +function print_usage() + println("Usage: genomicpro_cli.jl ") +end + +function main() + isempty(ARGS) && return (print_usage(); exit(1)) + config_path = ARGS[1] + config = read_pipeline_config(config_path) + result = run_pipeline(config) + println(render_summary(result)) +end + +main() diff --git a/src/GenomicPro.jl b/src/GenomicPro.jl new file mode 100644 index 0000000..43651a7 --- /dev/null +++ b/src/GenomicPro.jl @@ -0,0 +1,27 @@ +module GenomicPro + +export DatasetConfig, IntegrationConfig, AnalysisConfig, ReportConfig, QualityControlConfig, + PipelineConfig, PipelineResult, OmicsDataset, IntegratedDataset, + load_omics_dataset, read_pipeline_config, run_pipeline, align_datasets, + integrate_datasets, normalize!, impute_missing!, filter_low_quality!, + run_pca, run_kmeans, build_report, save_report, render_summary + +using CSV +using DataFrames +using JSON3 +using LinearAlgebra +using Statistics +using Random +using Dates + +include("types.jl") +include("utils.jl") +include("validation.jl") +include("io.jl") +include("preprocessing.jl") +include("integration.jl") +include("analysis.jl") +include("reporting.jl") +include("pipeline.jl") + +end diff --git a/src/analysis.jl b/src/analysis.jl new file mode 100644 index 0000000..3247126 --- /dev/null +++ b/src/analysis.jl @@ -0,0 +1,108 @@ +function run_pca(dataset::IntegratedDataset; components::Int=3) + components > 0 || throw(ArgumentError("Number of components must be positive")) + n_samples = length(dataset.samples) + n_samples > 1 || throw(ArgumentError("PCA requires at least two samples")) + mat = dataset.matrix + centered = mat .- mean(mat, dims=2) + data = transpose(centered) + svd_res = svd(data; full=false) + k = min(components, length(svd_res.S)) + scores = svd_res.U[:, 1:k] * Diagonal(svd_res.S[1:k]) + loadings = svd_res.V[:, 1:k] + denom = max(n_samples - 1, 1) + explained_variance = (svd_res.S .^ 2) / denom + explained_ratio = explained_variance ./ sum(explained_variance) + return PCAResult(scores, loadings, explained_variance[1:k], explained_ratio[1:k]) +end + +function _initialise_centroids(data::Matrix{Float64}, k::Int, seed::Int) + Random.seed!(seed) + indices = randperm(size(data, 2))[1:k] + return data[:, indices] +end + +function _assign_clusters(data::Matrix{Float64}, centroids::Matrix{Float64}) + assignments = Vector{Int}(undef, size(data, 2)) + for j in 1:size(data, 2) + column = view(data, :, j) + best_idx = 1 + best_dist = Inf + for c in 1:size(centroids, 2) + dist = sum((column .- view(centroids, :, c)).^2) + if dist < best_dist + best_dist = dist + best_idx = c + end + end + assignments[j] = best_idx + end + return assignments +end + +function _update_centroids(data::Matrix{Float64}, assignments::Vector{Int}, k::Int) + dim = size(data, 1) + centroids = zeros(Float64, dim, k) + counts = zeros(Int, k) + for j in 1:length(assignments) + cluster = assignments[j] + centroids[:, cluster] .+= view(data, :, j) + counts[cluster] += 1 + end + for c in 1:k + if counts[c] > 0 + centroids[:, c] ./= counts[c] + end + end + return centroids, counts +end + +function _compute_inertia(data::Matrix{Float64}, assignments::Vector{Int}, centroids::Matrix{Float64}) + total = 0.0 + for j in 1:size(data, 2) + cluster = assignments[j] + diff = view(data, :, j) .- view(centroids, :, cluster) + total += sum(diff .^ 2) + end + return total +end + +function run_kmeans(dataset::IntegratedDataset; k::Int=3, maxiter::Int=300, tol::Float64=1e-4, seed::Int=42) + k > 0 || throw(ArgumentError("Number of clusters must be positive")) + data = dataset.matrix + size(data, 2) >= k || throw(ArgumentError("Number of clusters cannot exceed sample count")) + centroids = _initialise_centroids(data, k, seed) + assignments = _assign_clusters(data, centroids) + previous_inertia = Inf + for iter in 1:maxiter + centroids, counts = _update_centroids(data, assignments, k) + for c in 1:k + if counts[c] == 0 + centroids[:, c] = data[:, rand(1:size(data, 2))] + end + end + assignments = _assign_clusters(data, centroids) + inertia = _compute_inertia(data, assignments, centroids) + if abs(previous_inertia - inertia) < tol + previous_inertia = inertia + break + end + previous_inertia = inertia + end + return ClusteringResult(assignments, centroids, _compute_inertia(data, assignments, centroids)) +end + +function run_analysis(dataset::IntegratedDataset, config::AnalysisConfig) + pca_result = nothing + if config.run_pca + pca_result = run_pca(dataset; components=config.pca_components) + end + clustering_result = nothing + if config.clustering === :kmeans + clustering_result = run_kmeans(dataset; k=config.cluster_count, seed=config.random_seed) + elseif config.clustering === nothing || config.clustering === :none + clustering_result = nothing + else + throw(ArgumentError("Unsupported clustering strategy: $(config.clustering)")) + end + return AnalysisResult(pca=pca_result, clustering=clustering_result) +end diff --git a/src/integration.jl b/src/integration.jl new file mode 100644 index 0000000..9baad32 --- /dev/null +++ b/src/integration.jl @@ -0,0 +1,44 @@ +function integrate_concatenate(datasets::Vector{OmicsDataset}) + combined_matrix = vcat([ds.matrix for ds in datasets]...) + combined_features = String[] + for ds in datasets + append!(combined_features, [string(ds.name, "::", feature) for feature in ds.features]) + end + metadata = Dict{String, Any}( + "strategy" => "concatenate", + "source_datasets" => [ds.name for ds in datasets] + ) + return IntegratedDataset("Integrated", combined_matrix, combined_features, copy(datasets[1].samples), metadata) +end + +function integrate_weighted_sum(datasets::Vector{OmicsDataset}, weights::Dict{String, Float64}) + reference_features = datasets[1].features + for ds in datasets[2:end] + ds.features == reference_features || throw(ArgumentError("Weighted sum integration requires identical feature ordering")) + end + weight_vec = [weights[ds.name] for ds in datasets] + weight_sum = sum(weight_vec) + weight_sum == 0 && throw(ArgumentError("Weights must not sum to zero")) + normalised = weight_vec ./ weight_sum + combined_matrix = zeros(Float64, size(datasets[1].matrix)) + for (w, ds) in zip(normalised, datasets) + combined_matrix .+= w .* ds.matrix + end + metadata = Dict{String, Any}( + "strategy" => "weighted_sum", + "weights" => Dict(ds.name => weights[ds.name] for ds in datasets) + ) + return IntegratedDataset("Integrated", combined_matrix, copy(reference_features), copy(datasets[1].samples), metadata) +end + +function integrate_datasets(datasets::Vector{OmicsDataset}, config::IntegrationConfig) + isempty(datasets) && throw(ArgumentError("No datasets provided for integration")) + if config.strategy == :concatenate + return integrate_concatenate(datasets) + elseif config.strategy == :weighted_sum + weights = _ensure_weights(datasets, config.weights) + return integrate_weighted_sum(datasets, weights) + else + throw(ArgumentError("Unsupported integration strategy: $(config.strategy)")) + end +end diff --git a/src/io.jl b/src/io.jl new file mode 100644 index 0000000..1492ac7 --- /dev/null +++ b/src/io.jl @@ -0,0 +1,111 @@ +function load_omics_dataset(config::DatasetConfig) + raw_missing = get(config.metadata, "missing_strings", String[]) + missing_strings = raw_missing isa AbstractVector ? [String(v) for v in raw_missing] : String[] + delim = _to_char(config.delimiter) + df = CSV.read(config.path, DataFrame; delim=delim, missingstring=missing_strings, ignorerepeated=true) + feature_idx = _feature_column_index(df, config.feature_column) + features = String.(df[!, feature_idx]) + sample_cols = Symbol[] + for (idx, col) in enumerate(names(df)) + idx == feature_idx && continue + push!(sample_cols, col) + end + samples = String.(sample_cols) + nfeatures = length(features) + nsamples = length(samples) + matrix = Matrix{Float64}(undef, nfeatures, nsamples) + for (j, colname) in enumerate(sample_cols) + matrix[:, j] = _coerce_to_float_column(df[!, colname]) + end + metadata = Dict{String, Any}( + "source_path" => config.path, + "delimiter" => String(delim), + "missing_strings" => missing_strings, + "sample_columns" => samples + ) + for (k, v) in config.metadata + metadata[k] = v + end + return OmicsDataset(config.name, matrix, features, samples, metadata) +end + +function read_pipeline_config(path::AbstractString) + config_json = JSON3.read(read(path, String)) + haskey(config_json, "datasets") || throw(ArgumentError("Pipeline configuration missing 'datasets' section")) + dataset_configs = DatasetConfig[] + for entry in config_json["datasets"] + metadata = Dict{String, Any}() + if haskey(entry, "metadata") && entry["metadata"] !== nothing + for (k, v) in entry["metadata"] + key = String(k) + if v isa AbstractVector + metadata[key] = [String(x) for x in v] + elseif v isa AbstractString + metadata[key] = String(v) + else + metadata[key] = v + end + end + end + feature_column = 1 + if haskey(entry, "feature_column") && entry["feature_column"] !== nothing + value = entry["feature_column"] + if value isa Integer + feature_column = Int(value) + else + feature_column = Symbol(String(value)) + end + end + push!(dataset_configs, DatasetConfig( + name = String(entry["name"]), + path = String(entry["path"]), + delimiter = haskey(entry, "delimiter") && entry["delimiter"] !== nothing ? _to_char(String(entry["delimiter"])) : ',', + feature_column = feature_column, + normalization = haskey(entry, "normalization") && entry["normalization"] !== nothing ? Symbol(String(entry["normalization"])) : :zscore, + log_base = haskey(entry, "log_base") && entry["log_base"] !== nothing ? Float64(entry["log_base"]) : nothing, + imputation = haskey(entry, "imputation") && entry["imputation"] !== nothing ? Symbol(String(entry["imputation"])) : :mean, + min_nonmissing_ratio = haskey(entry, "min_nonmissing_ratio") && entry["min_nonmissing_ratio"] !== nothing ? Float64(entry["min_nonmissing_ratio"]) : 0.7, + min_variance = haskey(entry, "min_variance") && entry["min_variance"] !== nothing ? Float64(entry["min_variance"]) : 1e-8, + metadata = metadata + )) + end + integration_section = haskey(config_json, "integration") ? config_json["integration"] : Dict() + weights_dict = Dict{String, Float64}() + if haskey(integration_section, "weights") && integration_section["weights"] !== nothing + for (k, v) in integration_section["weights"] + weights_dict[String(k)] = Float64(v) + end + end + integration_config = IntegrationConfig( + strategy = haskey(integration_section, "strategy") && integration_section["strategy"] !== nothing ? Symbol(String(integration_section["strategy"])) : :concatenate, + weights = weights_dict + ) + analysis_section = haskey(config_json, "analysis") ? config_json["analysis"] : Dict() + analysis_config = AnalysisConfig( + run_pca = Bool(get(analysis_section, "run_pca", true)), + pca_components = Int(get(analysis_section, "pca_components", 3)), + clustering = haskey(analysis_section, "clustering") && analysis_section["clustering"] !== nothing ? Symbol(String(analysis_section["clustering"])) : :kmeans, + cluster_count = Int(get(analysis_section, "cluster_count", 3)), + random_seed = Int(get(analysis_section, "random_seed", 42)) + ) + report_section = haskey(config_json, "report") ? config_json["report"] : Dict() + report_config = ReportConfig( + output_path = haskey(report_section, "output_path") && report_section["output_path"] !== nothing ? String(report_section["output_path"]) : nothing, + format = haskey(report_section, "format") && report_section["format"] !== nothing ? Symbol(String(report_section["format"])) : :json, + text_output_path = haskey(report_section, "text_output_path") && report_section["text_output_path"] !== nothing ? String(report_section["text_output_path"]) : nothing, + include_datasets = Bool(get(report_section, "include_datasets", true)), + include_analysis = Bool(get(report_section, "include_analysis", true)) + ) + qc_section = haskey(config_json, "quality_control") ? config_json["quality_control"] : Dict() + quality_control = QualityControlConfig( + min_samples = Int(get(qc_section, "min_samples", 1)), + min_features = Int(get(qc_section, "min_features", 1)), + min_nonmissing_ratio = Float64(get(qc_section, "min_nonmissing_ratio", 0.5)), + min_variance = Float64(get(qc_section, "min_variance", 1e-8)) + ) + return PipelineConfig(dataset_configs; + integration_config=integration_config, + analysis_config=analysis_config, + report_config=report_config, + quality_control=quality_control) +end diff --git a/src/pipeline.jl b/src/pipeline.jl new file mode 100644 index 0000000..14b4d56 --- /dev/null +++ b/src/pipeline.jl @@ -0,0 +1,17 @@ +function run_pipeline(config::PipelineConfig) + datasets = OmicsDataset[] + for ds_config in config.dataset_configs + dataset = load_omics_dataset(ds_config) + validate_dataset(dataset, config.quality_control) + impute_missing!(dataset, ds_config.imputation) + normalize!(dataset, ds_config.normalization; log_base=ds_config.log_base) + filter_low_quality!(dataset, ds_config, config.quality_control) + push!(datasets, dataset) + end + aligned = align_datasets(datasets) + integrated = integrate_datasets(aligned, config.integration_config) + analysis = run_analysis(integrated, config.analysis_config) + report = build_report(aligned, integrated, analysis, config) + save_report(report, config.report_config) + return PipelineResult(aligned, integrated, analysis, report) +end diff --git a/src/preprocessing.jl b/src/preprocessing.jl new file mode 100644 index 0000000..9be7226 --- /dev/null +++ b/src/preprocessing.jl @@ -0,0 +1,126 @@ +function impute_missing!(dataset::OmicsDataset, strategy::Symbol) + for row in eachrow(dataset.matrix) + if strategy == :mean + value = _nanmean(row) + isnan(value) && (value = 0.0) + _replace_nan!(row, value) + elseif strategy == :median + value = _nanmedian(row) + isnan(value) && (value = 0.0) + _replace_nan!(row, value) + elseif strategy == :zero + _replace_nan!(row, 0.0) + else + throw(ArgumentError("Unsupported imputation strategy: $(strategy)")) + end + end + dataset.metadata["imputation"] = String(strategy) + return dataset +end + +function log_transform!(dataset::OmicsDataset, base::Float64) + base > 0 || throw(ArgumentError("Log base must be positive")) + offset = 0.0 + minimum_value = minimum(dataset.matrix) + if minimum_value <= 0 + offset = abs(minimum_value) + 1e-6 + end + dataset.matrix .= log.(dataset.matrix .+ offset .+ 1e-9) ./ log(base) + dataset.metadata["log_base"] = base + dataset.metadata["log_offset"] = offset + return dataset +end + +function zscore_normalize!(dataset::OmicsDataset) + for row in eachrow(dataset.matrix) + μ = mean(row) + σ = std(row) + σ ≈ 0 && (σ = 1.0) + row .-= μ + row ./= σ + end + dataset.metadata["normalization"] = "zscore" + return dataset +end + +function minmax_scale!(dataset::OmicsDataset) + for row in eachrow(dataset.matrix) + minv = minimum(row) + maxv = maximum(row) + range = maxv - minv + range ≈ 0 && (range = 1.0) + row .-= minv + row ./= range + end + dataset.metadata["normalization"] = "minmax" + return dataset +end + +function robust_scale!(dataset::OmicsDataset) + for row in eachrow(dataset.matrix) + med = median(row) + q1 = quantile(row, 0.25) + q3 = quantile(row, 0.75) + iqr = q3 - q1 + iqr ≈ 0 && (iqr = 1.0) + row .-= med + row ./= iqr + end + dataset.metadata["normalization"] = "robust" + return dataset +end + +function normalize!(dataset::OmicsDataset, method::Symbol; log_base::Union{Nothing, Float64}=nothing) + if log_base !== nothing + log_transform!(dataset, log_base) + end + if method == :zscore + zscore_normalize!(dataset) + elseif method == :minmax + minmax_scale!(dataset) + elseif method == :robust + robust_scale!(dataset) + elseif method == :none + dataset.metadata["normalization"] = "none" + else + throw(ArgumentError("Unsupported normalization method: $(method)")) + end + return dataset +end + +function filter_low_quality!(dataset::OmicsDataset, config::DatasetConfig, qc::QualityControlConfig) + apply_quality_filters!(dataset, config, qc) + dataset.metadata["post_filter_features"] = length(dataset.features) + return dataset +end + +function align_datasets(datasets::Vector{OmicsDataset}) + isempty(datasets) && throw(ArgumentError("No datasets provided for alignment")) + sample_sets = map(ds -> Set(ds.samples), datasets) + common_samples = reduce(intersect, sample_sets) + isempty(common_samples) && throw(ArgumentError("Datasets do not share common samples")) + ordered_samples = [sample for sample in datasets[1].samples if sample in common_samples] + aligned = OmicsDataset[] + for ds in datasets + sample_index = Dict(sample => idx for (idx, sample) in enumerate(ds.samples)) + matrix = Matrix{Float64}(undef, size(ds.matrix, 1), length(ordered_samples)) + for (j, sample) in enumerate(ordered_samples) + idx = sample_index[sample] + matrix[:, j] = ds.matrix[:, idx] + end + metadata = copy(ds.metadata) + metadata["aligned_samples"] = ordered_samples + push!(aligned, OmicsDataset(ds.name, matrix, copy(ds.features), ordered_samples, metadata)) + end + return aligned +end + +function summarise_dataset(dataset::OmicsDataset) + Dict( + :name => dataset.name, + :samples => length(dataset.samples), + :features => length(dataset.features), + :missing_ratio => dataset_missing_ratio(dataset), + :metadata => dataset.metadata + ) +end diff --git a/src/reporting.jl b/src/reporting.jl new file mode 100644 index 0000000..38ceab0 --- /dev/null +++ b/src/reporting.jl @@ -0,0 +1,102 @@ +function dataset_statistics(dataset::OmicsDataset) + variances = dataset_variances(dataset) + variance_summary = isempty(variances) ? Dict(:mean => 0.0, :min => 0.0, :max => 0.0) : Dict( + :mean => mean(variances), + :min => minimum(variances), + :max => maximum(variances) + ) + Dict( + :name => dataset.name, + :samples => length(dataset.samples), + :features => length(dataset.features), + :missing_ratio => dataset_missing_ratio(dataset), + :variance_summary => variance_summary, + :metadata => dataset.metadata + ) +end + +function analysis_statistics(result::AnalysisResult) + stats = Dict{Symbol, Any}() + if result.pca !== nothing + stats[:pca] = Dict( + :explained_ratio => result.pca.explained_ratio, + :explained_variance => result.pca.explained_variance + ) + end + if result.clustering !== nothing + stats[:clustering] = Dict( + :assignments => result.clustering.assignments, + :inertia => result.clustering.inertia + ) + end + return stats +end + +function build_report(datasets::Vector{OmicsDataset}, integrated::IntegratedDataset, analysis::AnalysisResult, config::PipelineConfig) + report = Dict{Symbol, Any}( + :timestamp => Dates.format(Dates.now(), Dates.RFC3339), + :integrated => Dict( + :samples => length(integrated.samples), + :features => length(integrated.features), + :metadata => integrated.metadata + ), + :quality_control => Dict( + :min_samples => config.quality_control.min_samples, + :min_features => config.quality_control.min_features, + :min_nonmissing_ratio => config.quality_control.min_nonmissing_ratio, + :min_variance => config.quality_control.min_variance + ) + ) + if config.report_config.include_datasets + report[:datasets] = [dataset_statistics(ds) for ds in datasets] + end + if config.report_config.include_analysis + report[:analysis] = analysis_statistics(analysis) + end + return report +end + +function save_report(report::Dict{Symbol, Any}, config::ReportConfig) + if config.output_path !== nothing + mkpath(dirname(config.output_path)) + open(config.output_path, "w") do io + JSON3.write(io, report) + end + end + if config.text_output_path !== nothing + mkpath(dirname(config.text_output_path)) + open(config.text_output_path, "w") do io + println(io, render_summary(report)) + end + end + return report +end + +function render_summary(report::Dict{Symbol, Any}) + io = IOBuffer() + println(io, "GenomicPro Multi-Omics Report") + println(io, "Generated at: ", report[:timestamp]) + integrated = report[:integrated] + println(io, "Integrated dataset: $(integrated[:samples]) samples x $(integrated[:features]) features") + if haskey(report, :datasets) + println(io, "Datasets:") + for ds in report[:datasets] + println(io, " - $(ds[:name]): $(ds[:samples]) samples, $(ds[:features]) features, missing=$(round(ds[:missing_ratio], digits=3))") + end + end + if haskey(report, :analysis) + if haskey(report[:analysis], :pca) + ratios = report[:analysis][:pca][:explained_ratio] + println(io, "PCA explained variance ratios: ", join(round.(ratios, digits=3), ", ")) + end + if haskey(report[:analysis], :clustering) + inertia = report[:analysis][:clustering][:inertia] + println(io, "Clustering inertia: ", round(inertia, digits=3)) + end + end + return String(take!(io)) +end + +function render_summary(result::PipelineResult) + return render_summary(result.report) +end diff --git a/src/types.jl b/src/types.jl new file mode 100644 index 0000000..9e3196f --- /dev/null +++ b/src/types.jl @@ -0,0 +1,96 @@ +Base.@kwdef struct DatasetConfig + name::String + path::String + delimiter::Char = ',' + feature_column::Union{Int, Symbol, String} = 1 + normalization::Symbol = :zscore + log_base::Union{Nothing, Float64} = nothing + imputation::Symbol = :mean + min_nonmissing_ratio::Float64 = 0.7 + min_variance::Float64 = 1e-8 + metadata::Dict{String, Any} = Dict{String, Any}() +end + +Base.@kwdef struct IntegrationConfig + strategy::Symbol = :concatenate + weights::Dict{String, Float64} = Dict{String, Float64}() +end + +Base.@kwdef struct AnalysisConfig + run_pca::Bool = true + pca_components::Int = 3 + clustering::Union{Nothing, Symbol} = :kmeans + cluster_count::Int = 3 + random_seed::Int = 42 +end + +Base.@kwdef struct ReportConfig + output_path::Union{Nothing, String} = nothing + format::Symbol = :json + text_output_path::Union{Nothing, String} = nothing + include_datasets::Bool = true + include_analysis::Bool = true +end + +Base.@kwdef struct QualityControlConfig + min_samples::Int = 1 + min_features::Int = 1 + min_nonmissing_ratio::Float64 = 0.5 + min_variance::Float64 = 1e-8 +end + +struct PipelineConfig + dataset_configs::Vector{DatasetConfig} + integration_config::IntegrationConfig + analysis_config::AnalysisConfig + report_config::ReportConfig + quality_control::QualityControlConfig +end + +Base.@kwdef mutable struct OmicsDataset + name::String + matrix::Matrix{Float64} + features::Vector{String} + samples::Vector{String} + metadata::Dict{String, Any} = Dict{String, Any}() +end + +Base.@kwdef mutable struct IntegratedDataset + name::String + matrix::Matrix{Float64} + features::Vector{String} + samples::Vector{String} + metadata::Dict{String, Any} = Dict{String, Any}() +end + +struct PCAResult + scores::Matrix{Float64} + loadings::Matrix{Float64} + explained_variance::Vector{Float64} + explained_ratio::Vector{Float64} +end + +struct ClusteringResult + assignments::Vector{Int} + centroids::Matrix{Float64} + inertia::Float64 +end + +Base.@kwdef struct AnalysisResult + pca::Union{Nothing, PCAResult} = nothing + clustering::Union{Nothing, ClusteringResult} = nothing +end + +struct PipelineResult + datasets::Vector{OmicsDataset} + integrated::IntegratedDataset + analysis::AnalysisResult + report::Dict{Symbol, Any} +end + +PipelineConfig(dataset_configs::Vector{DatasetConfig}; + integration_config::IntegrationConfig=IntegrationConfig(), + analysis_config::AnalysisConfig=AnalysisConfig(), + report_config::ReportConfig=ReportConfig(), + quality_control::QualityControlConfig=QualityControlConfig()) = + PipelineConfig(dataset_configs, integration_config, analysis_config, report_config, quality_control) diff --git a/src/utils.jl b/src/utils.jl new file mode 100644 index 0000000..357ecea --- /dev/null +++ b/src/utils.jl @@ -0,0 +1,119 @@ +function _stringify(v) + v === nothing && return "" + v isa AbstractString && return String(v) + return string(v) +end + +function _to_char(delim) + delim isa Char && return delim + delim isa AbstractString && !isempty(delim) && return delim[1] + throw(ArgumentError("Delimiter must be a single character, got $(delim)")) +end + +function _feature_column_index(df::DataFrame, feature_column) + if feature_column isa Int + 1 <= feature_column <= ncol(df) || throw(ArgumentError("feature_column index out of range")) + return feature_column + elseif feature_column isa Symbol + return findfirst(==(feature_column), propertynames(df)) || throw(ArgumentError("feature column $(feature_column) not found")) + elseif feature_column isa AbstractString + return findfirst(==(Symbol(feature_column)), propertynames(df)) || throw(ArgumentError("feature column $(feature_column) not found")) + else + throw(ArgumentError("Unsupported feature_column type: $(typeof(feature_column))")) + end +end + +function _parse_float(value) + if value === missing + return NaN + elseif value isa Number + return Float64(value) + elseif value isa AbstractString + stripped = strip(value) + isempty(stripped) && return NaN + try + return parse(Float64, stripped) + catch err + throw(ArgumentError("Cannot parse numeric value from '$(value)': $(err)")) + end + else + return Float64(value) + end +end + +function _coerce_to_float_column(column) + result = Vector{Float64}(undef, length(column)) + for (i, v) in enumerate(column) + result[i] = _parse_float(v) + end + return result +end + +function _nanmask(vec::AbstractVector{Float64}) + mask = BitVector(undef, length(vec)) + @inbounds for i in eachindex(vec) + mask[i] = isnan(vec[i]) + end + return mask +end + +function _nanmean(vec::AbstractVector{Float64}) + total = 0.0 + count = 0 + @inbounds for v in vec + if !isnan(v) + total += v + count += 1 + end + end + return count == 0 ? NaN : total / count +end + +function _nanmedian(vec::AbstractVector{Float64}) + filtered = filter(!isnan, vec) + isempty(filtered) && return NaN + sorted = sort(filtered) + mid = length(sorted) ÷ 2 + if isodd(length(sorted)) + return sorted[mid + 1] + else + return (sorted[mid] + sorted[mid + 1]) / 2 + end +end + +function _replace_nan!(vec::AbstractVector{Float64}, value::Float64) + @inbounds for i in eachindex(vec) + if isnan(vec[i]) + vec[i] = value + end + end + return vec +end + +function _variance(vec::AbstractVector{Float64}) + clean = filter(!isnan, vec) + length(clean) <= 1 && return 0.0 + μ = mean(clean) + return sum((x - μ)^2 for x in clean) / (length(clean) - 1) +end + +function _center_rows!(matrix::Matrix{Float64}) + for i in axes(matrix, 1) + row = view(matrix, i, :) + μ = mean(row) + row .-= μ + end + return matrix +end + +function _nan_ratio(vec::AbstractVector{Float64}) + count = count(isnan, vec) + return length(vec) == 0 ? 0.0 : count / length(vec) +end + +function _ensure_weights(datasets::Vector{OmicsDataset}, weights::Dict{String, Float64}) + isempty(weights) && return Dict(ds.name => 1.0 for ds in datasets) + missing = filter(name -> !haskey(weights, name), getfield.(datasets, :name)) + isempty(missing) || throw(ArgumentError("Missing weights for datasets: $(join(missing, ", "))")) + return weights +end diff --git a/src/validation.jl b/src/validation.jl new file mode 100644 index 0000000..ba4ed90 --- /dev/null +++ b/src/validation.jl @@ -0,0 +1,38 @@ +function validate_dataset(dataset::OmicsDataset, qc::QualityControlConfig) + nsamples = length(dataset.samples) + nfeatures = length(dataset.features) + nsamples >= qc.min_samples || throw(ArgumentError("Dataset $(dataset.name) does not have enough samples")) + nfeatures >= qc.min_features || throw(ArgumentError("Dataset $(dataset.name) does not have enough features")) + length(unique(dataset.samples)) == nsamples || throw(ArgumentError("Dataset $(dataset.name) has duplicate sample identifiers")) + length(unique(dataset.features)) == nfeatures || throw(ArgumentError("Dataset $(dataset.name) has duplicate feature identifiers")) + return dataset +end + +function dataset_missing_ratio(dataset::OmicsDataset) + isempty(dataset.features) && return 0.0 + ratios = map(row -> _nan_ratio(row), eachrow(dataset.matrix)) + return mean(ratios) +end + +function dataset_variances(dataset::OmicsDataset) + return map(row -> _variance(row), eachrow(dataset.matrix)) +end + +function apply_quality_filters!(dataset::OmicsDataset, config::DatasetConfig, qc::QualityControlConfig) + mask = trues(length(dataset.features)) + filtered = 0 + for (i, row) in enumerate(eachrow(dataset.matrix)) + missing_ratio = _nan_ratio(row) + variance = _variance(row) + if missing_ratio > (1 - config.min_nonmissing_ratio) || variance < max(config.min_variance, qc.min_variance) + mask[i] = false + filtered += 1 + end + end + if filtered > 0 + dataset.matrix = dataset.matrix[mask, :] + dataset.features = dataset.features[mask] + end + dataset.metadata["filtered_features"] = filtered + return dataset +end diff --git a/test/config/test_pipeline.json b/test/config/test_pipeline.json new file mode 100644 index 0000000..3c4c0b0 --- /dev/null +++ b/test/config/test_pipeline.json @@ -0,0 +1,43 @@ +{ + "datasets": [ + { + "name": "transcriptomics", + "path": "test/data/transcriptomics.csv", + "delimiter": ",", + "feature_column": "Feature", + "normalization": "zscore", + "imputation": "mean", + "metadata": { + "missing_strings": ["NA"] + } + }, + { + "name": "proteomics", + "path": "test/data/proteomics.csv", + "delimiter": ",", + "feature_column": "Feature", + "normalization": "zscore", + "imputation": "median" + } + ], + "integration": { + "strategy": "concatenate" + }, + "analysis": { + "run_pca": true, + "pca_components": 2, + "clustering": "kmeans", + "cluster_count": 2, + "random_seed": 7 + }, + "report": { + "include_datasets": true, + "include_analysis": true + }, + "quality_control": { + "min_samples": 3, + "min_features": 2, + "min_nonmissing_ratio": 0.5, + "min_variance": 1e-8 + } +} diff --git a/test/data/proteomics.csv b/test/data/proteomics.csv new file mode 100644 index 0000000..55aadce --- /dev/null +++ b/test/data/proteomics.csv @@ -0,0 +1,5 @@ +Feature,SampleA,SampleB,SampleC +Prot1,2.1,2.3,2.0 +Prot2,1.2,1.5,1.1 +Prot3,0.9,1.0,NA +Prot4,1.8,1.9,1.7 diff --git a/test/data/transcriptomics.csv b/test/data/transcriptomics.csv new file mode 100644 index 0000000..08b1aef --- /dev/null +++ b/test/data/transcriptomics.csv @@ -0,0 +1,5 @@ +Feature,SampleA,SampleB,SampleC +Gene1,10,12,8 +Gene2,5,7,6 +Gene3,3,NA,4 +Gene4,8,9,7 diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..801e941 --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,38 @@ +using Test +using GenomicPro + +@testset "Configuration" begin + config = read_pipeline_config("test/config/test_pipeline.json") + @test length(config.dataset_configs) == 2 + @test config.integration_config.strategy == :concatenate + @test config.analysis_config.cluster_count == 2 +end + +@testset "Dataset loading" begin + config = read_pipeline_config("test/config/test_pipeline.json") + ds = load_omics_dataset(config.dataset_configs[1]) + @test size(ds.matrix, 1) == 4 + @test size(ds.matrix, 2) == 3 + @test ds.features[1] == "Gene1" +end + +@testset "Preprocessing" begin + config = read_pipeline_config("test/config/test_pipeline.json") + ds = load_omics_dataset(config.dataset_configs[1]) + impute_missing!(ds, :mean) + @test !any(isnan, ds.matrix) + normalize!(ds, :zscore) + @test isapprox(mean(ds.matrix[1, :]), 0.0; atol=1e-8) +end + +@testset "Integration and analysis" begin + config = read_pipeline_config("test/config/test_pipeline.json") + result = run_pipeline(config) + @test length(result.integrated.samples) == 3 + @test size(result.integrated.matrix, 1) == 8 + @test result.analysis.pca !== nothing + @test result.analysis.clustering !== nothing + @test length(result.analysis.clustering.assignments) == 3 + @test haskey(result.report, :datasets) + @test haskey(result.report, :analysis) +end