diff --git a/Makefile b/Makefile index 4aaa3ee..ab30f14 100755 --- a/Makefile +++ b/Makefile @@ -72,6 +72,7 @@ clean: @rm -rf .nextflow rm -rf tmp find work/ -type f -print0 | xargs -0 -P 8 -n 100 rm -f + @rm -rf work @rm -rf plugins @rm -rf .config @rm -f .bash_history diff --git a/bin/misc_utils/checkers.R b/bin/misc_utils/checkers.R index 84c36ef..65ea28e 100755 --- a/bin/misc_utils/checkers.R +++ b/bin/misc_utils/checkers.R @@ -1,10 +1,54 @@ -# Check if matched sample names, then return those check_common_samples <- function(data) { - sample_names <- lapply(data$X, rownames) |> unlist() |> unique() - matched <- length(sample_names) == length(data$Y) - if (matched) { - return(sample_names) + if (length(data$X) == 0L) { + stop("data$X is empty.") + } + + sample_names_list <- lapply(data$X, rownames) + + valid_names <- vapply( + sample_names_list, + function(ids) { + !is.null(ids) && + length(ids) > 0L && + !anyNA(ids) && + all(nzchar(ids)) && + anyDuplicated(ids) == 0L + }, + logical(1) + ) + + if (!all(valid_names)) { + stop("Each modality must have non-empty, unique sample row names.") + } + + sample_names <- sample_names_list[[1]] + + matched <- vapply( + sample_names_list, + function(ids) identical(ids, sample_names), + logical(1) + ) + + if (!all(matched)) { + stop("Sample names or their order differ between modalities.") + } + + if (is.null(data$Y)) { + stop("data$Y is missing.") + } + + y_count <- if (is.null(dim(data$Y))) { + length(data$Y) } else { - stop("Failed to find match samples, check!") + nrow(data$Y) + } + + message("Found ", length(sample_names), " samples in each modality.") + message("Found ", y_count, " observations in Y.") + + if (length(sample_names) != y_count) { + stop("Sample count in X does not match observation count in Y.") } -} + + return(sample_names) +} \ No newline at end of file diff --git a/bin/misc_utils/extract_Xy.R b/bin/misc_utils/extract_Xy.R index be69562..0977813 100755 --- a/bin/misc_utils/extract_Xy.R +++ b/bin/misc_utils/extract_Xy.R @@ -29,7 +29,7 @@ parseY <- function(Y, verbose = FALSE) { } # Convert MAE to list of X and Y -extract_Xy <- function(mae, verbose_target=FALSE) { +extract_Xy <- function(mae, outcome_type="classification", verbose_target=FALSE) { # Note need to transpose back to p * n # TODO: add check for dimension match # COOP LR do not like delayed matrix, so transform it to S3 matrix @@ -43,11 +43,22 @@ extract_Xy <- function(mae, verbose_target=FALSE) { # MAE response would be a dataframe? # In simulated data this is always atomic, but in real, they are # transformed to dataframe first, so need to pull it - y_temp <- mae$response - if (is.data.frame(y_temp)) { - Y <- parseY(y_temp |> dplyr::pull(response), verbose = verbose_target) - } else { - Y <- parseY(y_temp, verbose = verbose_target) + if (outcome_type == "survival") { + #stop("Survival outcome type is not implemented yet") + y_df <- SummarizedExperiment::colData(mae) |> as.data.frame() + # And get the time and status columns + out <- y_df[, c("time", "status")] + return(list(X=X, Y=out)) + } + if (outcome_type == "classification") { + y_temp <- mae$response + if (is.data.frame(y_temp)) { + Y <- parseY(y_temp |> dplyr::pull(response), verbose = verbose_target) + } else { + Y <- parseY(y_temp, verbose = verbose_target) + } + return(list(X=X, Y=Y)) } - return(list(X=X, Y=Y)) + # Otherwise stop + stop("Unsupported outcome_type: ", outcome_type) } diff --git a/bin/python_utils/combine_mdata2df.py b/bin/python_utils/combine_mdata2df.py new file mode 100644 index 0000000..844e5c2 --- /dev/null +++ b/bin/python_utils/combine_mdata2df.py @@ -0,0 +1,46 @@ +import mudata +import pandas as pd +import numpy as np + + +# def combine_mdata2df(mdata, concat=True): + +# mod_names = list(mdata.mod.keys()) + +# X_df = pd.concat( [mdata[k].to_df().add_prefix(f"{k}_") for k in mod_names], axis=1 ) +# # Also extract the observation df +# y_df = mdata[mod_names[0]].obs +# # Should contain response column and is of string yes or no +# assert y_df["response"].isin(['yes', 'no']).all(), "Column contains values other than 'yes' or 'no'" +# # Converting to numeric binary +# y_df.loc[:, "response"] = np.where(y_df["response"] == "yes", 1, 0) +# # When supplied not to concat meaning y requires other informations more than response +# if not concat: +# return X_df, y_df +# # Merge all DataFrames together, and drop irrevelant meta information +# merged_df = pd.concat([ X_df, y_df[["response"]] ] , axis=1) +# return merged_df, mod_names + + +def combine_mdata2df(mdata, target_col="response"): + # Takes in a mudata and extract X and y components + # Get the Xs as a dataframe of combining all modality together columnwise + mod_names = list(mdata.mod.keys()) + # We also add the modality in front of every feature just like "epigenomics_some_feature_name" + X_df = pd.concat( + [ + mdata[mod].to_df().add_prefix(f"{mod}_") + for mod in mod_names + ], + axis=1, + ) + # Also extract the observation df + y_df = mdata[mod_names[0]].obs.copy() + + assert y_df[target_col].isin(["yes", "no"]).all(), ( + "Column contains values other than 'yes' or 'no'" + ) + # Converting to numeric binary + y_df[target_col] = np.where(y_df[target_col] == "yes", 1, 0) + + return X_df, y_df, mod_names diff --git a/modules/sklearn/train/resources/usr/bin/load_classifier_class.py b/bin/python_utils/load_classifier_class.py old mode 100644 new mode 100755 similarity index 55% rename from modules/sklearn/train/resources/usr/bin/load_classifier_class.py rename to bin/python_utils/load_classifier_class.py index b9d91eb..ce5ad86 --- a/modules/sklearn/train/resources/usr/bin/load_classifier_class.py +++ b/bin/python_utils/load_classifier_class.py @@ -1,6 +1,11 @@ +# This script is imported for sklearn classifiers usage + + import scipy.stats as stats import importlib +print("This date should be shown as : 2026") + def load_classifier_class(model_name, random_state=42, probability=True): # Define models with their respective classes, default parameters, and distributions @@ -9,11 +14,16 @@ def load_classifier_class(model_name, random_state=42, probability=True): # https://stackoverflow.com/questions/33843981/under-what-parameters-are-svc-and-linearsvc-in-scikit-learn-equivalent # Common params - C = stats.loguniform(1e-4, 1e4) - min_sample_leaf = stats.randint(1,6) - n_estimators = stats.randint(50, 501) - learning_rate = stats.uniform(0.01, 1.1) - max_features = ["sqrt", "log2", 100, 500, 1000, None] + C = stats.loguniform(1e-4, 1e4) + min_sample_leaf = stats.randint(1,6) + base_n_estimators = 10 + n_estimators = stats.randint(50, 501) + learning_rate = stats.uniform(0.01, 1.1) + max_features = ["sqrt", "log2", 100, 500, 1000, None] + max_depth = 10 + + + # Dict to store relevant information of sklearn classifiers # For logistic regression, it has built-in predict proba and coef @@ -28,21 +38,14 @@ def load_classifier_class(model_name, random_state=42, probability=True): "default_params": {"C": 1.0, "kernel": "linear", "random_state": random_state, "probability": probability}, "params_dist": {"C": C, "kernel": ["linear"] } } - # Decision Tree classifier, risk of overfitting, hence require pruning of trees - decision_tree_dict = { - "class_path": "sklearn.tree.DecisionTreeClassifier", - "default_params": {"max_depth": 10, "random_state": random_state}, - "params_dist": { - "criterion": ['gini', 'entropy', 'log_loss'], - "max_depth": stats.randint(5, 41), - "min_samples_leaf": min_sample_leaf, - "max_leaf_nodes": [10, 100, 1000, None] - } - } # Random Forest classifier, should in general work better than single decision tree random_forest_dict = { "class_path": "sklearn.ensemble.RandomForestClassifier", - "default_params": {"n_estimators": 10, "max_features": "sqrt", "max_depth": 10, "random_state": random_state, "n_jobs": -1}, + "default_params": { + "n_estimators": base_n_estimators, "max_features": "sqrt", + "max_depth": max_depth, + "random_state": random_state, "n_jobs": -1 + }, "params_dist": { "max_features": max_features, "max_leaf_nodes": [10, 100, 1000, None], @@ -50,32 +53,30 @@ def load_classifier_class(model_name, random_state=42, probability=True): } } - # Gradient Boost classifier, should tune large number of estimators with slow learning rate - gradient_boost_dict = { - "class_path": "sklearn.ensemble.GradientBoostingClassifier", - "default_params": {"loss":'log_loss', "learning_rate":0.1, "n_estimators":10}, - "params_dist": { - "n_estimators": n_estimators, - "learning_rate": learning_rate, - "min_samples_leaf": min_sample_leaf, - "max_features": max_features - } - } # MLPClassifier is a neural network classifier mlp_dict = { - "class_path": "sklearn.neural_network.MLPClassifier", - "default_params": {"hidden_layer_sizes": (64,), "max_iter": 500, - "random_state": random_state}, - "params_dist": {"hidden_layer_sizes": [(32,), (64,), (128,), (64, 32)], - "alpha": stats.loguniform(1e-5, 1e-1)} + "class_path": "sklearn.neural_network.MLPClassifier", + "default_params": { + "hidden_layer_sizes": (64,), "max_iter": 500, + "random_state": random_state + }, + "params_dist": { + "hidden_layer_sizes": [(32,), (64,), (128,), (64, 32)], + "alpha": stats.loguniform(1e-5, 1e-1) + } } - # Mimic xgboost - hist_gradient_boost_dict = { - "class_path": "sklearn.ensemble.HistGradientBoostingClassifier", - "default_params": {"random_state": random_state}, - "params_dist": {"learning_rate": learning_rate, - "max_leaf_nodes": [15, 31, 63], - "min_samples_leaf": stats.randint(5, 30)} + + # Real xgboost + xgboost_dict = { + "class_path": "xgboost.XGBClassifier", + "default_params": { + "random_state": random_state , + "n_estimators": base_n_estimators, + "objective": "binary:logistic"}, + "params_dist": { + "learning_rate": learning_rate, + "max_depth": stats.randint(2, 9) + } } # ======================== @@ -83,15 +84,12 @@ def load_classifier_class(model_name, random_state=42, probability=True): model_info = { "Logit": logit_dict, "Linear_SVM": linear_svm_dict, - # Decision Tree have risk of overfitting, hence require pruning of trees - "Decision_Tree": decision_tree_dict , - # RandomForest shuold in general work better than single decision tree + # RandomForest should in general work better than single decision tree "Random_Forest": random_forest_dict, - # GradientBoost should tune large number of estimators with slow learning rate - "Gradient_Boost": gradient_boost_dict, # MLPClassifier is a neural network classifier "MLP": mlp_dict, - "Hist_Gradient_Boost": hist_gradient_boost_dict + # XGBoost goes here + "XGBoost": xgboost_dict } # Check if valid name of model was input diff --git a/bin/rhelpers.R b/bin/rhelpers.R index cdeb8da..a4b7aeb 100755 --- a/bin/rhelpers.R +++ b/bin/rhelpers.R @@ -22,3 +22,71 @@ opt2num <- function(opt_chr) { as.numeric(x), x)) return(opt) } + + +make_parameter_record <- function(value, treatment) { + allowed_treatments <- c( + "tuned", + "fixed", + "default", + "data-derived" + ) + + if (!treatment %in% allowed_treatments) { + stop( + "Unknown parameter treatment: ", + treatment + ) + } + + list( + value = value, + treatment = treatment + ) +} + + +write_selected_hyperparameters <- function( + dataset_name, + method_name, + parameters, + selection = NULL, + output_path = NULL +) { + if (!requireNamespace("jsonlite", quietly = TRUE)) { + stop("The jsonlite package is required") + } + + if (is.null(output_path)) { + safe_method_name <- method_name |> + tolower() |> + gsub("[^a-z0-9_-]+", "_", x = _) + + output_path <- paste0( + safe_method_name, + "-", + dataset_name, + "_selected_hyperparameters.json" + ) + } + + result <- list( + dataset = dataset_name, + method = method_name, + analysis_stage = "model_selection", + selection = selection, + parameters = parameters + ) + + jsonlite::write_json( + result, + path = output_path, + pretty = TRUE, + auto_unbox = TRUE, + null = "null", + na = "null", + digits = NA + ) + + invisible(output_path) +} \ No newline at end of file diff --git a/bin/split_mae.R b/bin/split_mae.R index bcdc993..dc14c27 100755 --- a/bin/split_mae.R +++ b/bin/split_mae.R @@ -13,6 +13,9 @@ Options: --split_dir=SPLIT_DIR Directory containing list of txt file [default: empty] --dataset_name=NAME Name of dataset that is splitting [default: empty] --transpose Transpose the data as method requires [default: False] + --outcome_type=TYPE classification or survival [default: classification] + --time_col=TIME_COL colData column of survival time [default: time] + --status_col=STAT_COL colData column of event indicator [default: status] " # Parase docopt @@ -63,21 +66,60 @@ load_test_splits <- function(split_dir, pattern=".txt", ...) { # } -# Use this function to reconstruct mae -reconstruct_mae <- function(mae) { - # Given an mae with delayed matrices, we could load it into - # memory and make it of HDF5 arrays instead - X <- mae@ExperimentList |> lapply(as.matrix) - y <- mae$response - # Construct MAE - new_mae <- MultiAssayExperiment::MultiAssayExperiment(experiments = X) - new_mae$response <- y - return(new_mae) -} +reconstruct_mae <- function( + mae, + outcome_type = "classification", + response_col = "response", + time_col = "time", + status_col = "status" +) { + if (!outcome_type %in% c("classification", "survival")) { + stop("outcome_type must be 'classification' or 'survival'") + } + + cd <- SummarizedExperiment::colData(mae) |> as.data.frame() + + required_cols <- if (outcome_type == "classification") { + response_col + } else { + c(time_col, status_col) + } + + missing_cols <- setdiff(required_cols, colnames(cd)) + if (length(missing_cols) > 0L) { + stop("Missing outcome columns: ", paste(missing_cols, collapse = ", ")) + } + # Materialize each experiment as an in-memory matrix. + X <- lapply( + MultiAssayExperiment::experiments(mae), + function(x) { + # For SummarizedExperiment, use its first assay. + if (methods::is(x, "SummarizedExperiment")) { + x <- SummarizedExperiment::assay(x) + } + as.matrix(x) + } + ) + + if (outcome_type == "classification") { + cd$response <- cd[[response_col]] + } else { + # Extract both before assignment in case source names overlap. + time <- cd[[time_col]] + status <- cd[[status_col]] + cd$time <- time + cd$status <- status + } + + MultiAssayExperiment::MultiAssayExperiment( + experiments = X, + colData = cd + ) +} # Actual fun to split each MAE to train and test portion -split_mae <- function(mae_path, split_dir, dataset_name) { +split_mae <- function(mae_path, split_dir, dataset_name, outcome_type="classification") { # Read in the MAE # Note the prefix "" is required here? mae <- MultiAssayExperiment::loadHDF5MultiAssayExperiment(dir=mae_path, prefix="") @@ -92,8 +134,8 @@ split_mae <- function(mae_path, split_dir, dataset_name) { # First subset both split <- test_splits[[fold_name]] # TODO: Transpose data only when method requires it to - tr_mae <- mae[, -split, drop=TRUE] |> reconstruct_mae() - te_mae <- mae[, split, drop=TRUE] |> reconstruct_mae() + tr_mae <- mae[, -split, drop=TRUE] |> reconstruct_mae(outcome_type = outcome_type) + te_mae <- mae[, split, drop=TRUE] |> reconstruct_mae(outcome_type = outcome_type) # Then save each fold's train and test portion as subdirectory of fold name cat("\nSaving for", fold_name, "\n") if (!dir.exists(fold_name)) { @@ -111,4 +153,4 @@ split_mae <- function(mae_path, split_dir, dataset_name) { } } -split_mae(mae_path=opt$mae_path, split_dir=opt$split_dir, dataset_name=opt$dataset_name) +split_mae(mae_path=opt$mae_path, split_dir=opt$split_dir, dataset_name=opt$dataset_name, outcome_type=opt$outcome_type) diff --git a/conf/base.config b/conf/base.config index 44db74b..d34b0b9 100644 --- a/conf/base.config +++ b/conf/base.config @@ -16,7 +16,7 @@ process { time = { check_max( 4.h * task.attempt, 'time' ) } errorStrategy = { task.exitStatus in ((130..145) + 104) ? 'retry' : 'finish' } - maxRetries = 1 + maxRetries = 2 maxErrors = '-1' // Resources allocation group by usage of resources diff --git a/conf/dev-sockeye.config b/conf/dev-sockeye.config index e56581b..054603e 100644 --- a/conf/dev-sockeye.config +++ b/conf/dev-sockeye.config @@ -25,50 +25,27 @@ params { project_data_dir = allocation_code ? "/arc/project/${allocation_code}/datasets/messi_demo_data" : "invalid project data dir" } +apptainer { + enabled = true + autoMounts = true + cacheDir = "${params.apptainer_cache_dir}" + libraryDir = "${params.apptainer_cache_dir}" + runOptions = "--bind ${params.project_data_dir}" + docker.enabled = false + conda.enabled = false +} + executor { - name = 'local' +// name = 'local' cpus = params.max_cpus memory = params.max_memory } process { scratch = true - // NOTE this overrides the base config resource requirement - // withLabel:process_single { - // cpus = { check_max( 1 , 'cpus' ) } - // memory = { check_max( 256.MB * task.attempt, 'memory' ) } - // time = { check_max( 15.min * task.attempt, 'time' ) } - // } - // withLabel:process_low { - // cpus = { check_max( 2 * task.attempt, 'cpus' ) } - // memory = { check_max( 512.MB * task.attempt, 'memory' ) } - // time = { check_max( 30.min * task.attempt, 'time' ) } - // } - // withLabel:process_medium { - // cpus = { check_max( 4 * task.attempt, 'cpus' ) } - // memory = { check_max( 1.GB * task.attempt, 'memory' ) } - // time = { check_max( 1.h * task.attempt, 'time' ) } - // } - // withLabel:process_high { - // cpus = { check_max( 6 * task.attempt, 'cpus' ) } - // memory = { check_max( 2.GB * task.attempt, 'memory' ) } - // time = { check_max( 2.h * task.attempt, 'time' ) } - // } - - withLabel: 'mae_mu' { + withLabel: 'mae_mu' { container = "save_simulate.sif" - } + } } -// replace to apptainer later -// after supplying cacheDir, you could use container with image name only like above, -// no need to specify file:// ... -apptainer { - enabled = true - autoMounts = true - // Add the bind path to project space demo data - runOptions = "--bind ${params.project_data_dir} --bind '/scratch/st-singha53-1/tliang19/multiomics_debug_space/data/all_tar_gz'" - cacheDir = params.apptainer_cache_dir - docker.enabled = false -} diff --git a/conf/modules.config b/conf/modules.config index bafd192..11a0eae 100644 --- a/conf/modules.config +++ b/conf/modules.config @@ -8,11 +8,11 @@ process { if (params.publish_relevant) { process { // Suppress all publish output that are not these 3 processes - withName: '^(?!.*(MERGE_SELECTED_FEATURES|MERGE_RESULT_TABLE|PARSE_METADATA|CALCULATE_METRICS)).*$' { + withName: '^(?!.*(MERGE_SELECTED_FEATURES|MERGE_RESULT_TABLE|MERGE_SELECTED_HYPERPARAMETERS|PARSE_METADATA|CALCULATE_METRICS)).*$' { publishDir = [ enabled: false ] } } } - \ No newline at end of file + diff --git a/conf/real_data.config b/conf/real_data.config index 6b881b7..69b92f3 100644 --- a/conf/real_data.config +++ b/conf/real_data.config @@ -13,10 +13,14 @@ params { config_profile_name = 'Real Data Profile' config_profile_description = 'Evaluation on real datasets only' - + outcome_type = 'classification' // Full test mode, run all methods and all data //samplesheet = "${params.data_dir}/samplesheet_real_data.csv" + // And running single modality ones + single_modality_mode = true + // trigger to skip method or not + skip_sklearn = false skip_caret_multimodal = false skip_cplr = false skip_diablo = false diff --git a/conf/run_survival.config b/conf/run_survival.config new file mode 100644 index 0000000..1d73695 --- /dev/null +++ b/conf/run_survival.config @@ -0,0 +1,38 @@ +/* +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + Nextflow config file for running on survival tasks +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + Defines input files and everything required to run on real datasets + + Use as follows: + nextflow run main.nf -profile survival, + +---------------------------------------------------------------------------------------- +*/ + +params { + config_profile_name = 'Survival only' + config_profile_description = 'Evaluation on survival task' + outcome_type = 'survival' + + // Full test mode, run all methods and all data + single_modality_mode = true + // trigger to skip method or not + skip_sklearn = true + skip_caret_multimodal = true + skip_cplr = false + skip_diablo = true + skip_integrao = true + skip_rgcca = true + skip_mogonet = true + skip_mofa = false + skip_sksurv = false + // Run feature selection + selectFeature = false + k_fold_number = 5 + filter_low_var = "1" // This get casted as bool in both Python and R + // number of component + num_comps = 2 + // HE for mogonet + he_base_dim = 100 // Base dimension for the HE layer, default to 100 for real data +} diff --git a/conf/simulated_data.config b/conf/simulated_data.config index d78f9d0..8100dfc 100644 --- a/conf/simulated_data.config +++ b/conf/simulated_data.config @@ -13,11 +13,17 @@ params { config_profile_name = 'Simluated Data Profile' config_profile_description = 'Evaluation on simulated datasets only' + outcome_type = 'classification' + runDownstreamAnalysis = false // Skip downstream analysis // Full test mode, run all methods and all data //samplesheet = "${params.data_dir}/samplesheet_simulated_data.csv" + // Run single mod + single_modality_mode = true + // skippers skip_caret_multimodal = false + skip_sklearn = false skip_cplr = false skip_diablo = false skip_integrao = false diff --git a/launch_MESSI_pipeline.sh b/launch_MESSI_pipeline.sh index 74b5630..c3726d5 100755 --- a/launch_MESSI_pipeline.sh +++ b/launch_MESSI_pipeline.sh @@ -36,9 +36,13 @@ module load CVMFS_CC # module load java/11.0.16_8 # module load nextflow/23.04.3 # THESE ARE NEWER VERSIONS -module load apptainer/1.3.4 -module load java/17.0.6 -module load nextflow/24.04.4 +#module load apptainer/1.3.4 +#module load java/17.0.6 +#module load nextflow/24.04.4 +# Even newer +module load apptainer/1.3.5 +module load java/21.0.1 +module load nextflow/24.10.2 # Source the load script with env vars setup #source bin/helper.sh # ============================================================================= @@ -58,8 +62,9 @@ export NXF_OFFLINE='true' # The NXF script to run, located on the repo root directory NXF_SRC_MAIN=$PIPELINE_DIR/main.nf # Profile order matters, since the later one overrides the prior ones -PROFILE=sockeye,simulated_data +#PROFILE=sockeye,simulated_data #PROFILE=sockeye,real_data +PROFILE=sockeye,run_survival # Or use this one for development usage #PROFILE=sockeye,test,debug #PARAMS_FILE=remote_params.yaml @@ -68,26 +73,43 @@ PROFILE=sockeye,simulated_data timestamp=$(date +"%Y%m%d_%H%M%S") PROFILE_SUFFIX=${PROFILE#sockeye,} PROFILE_SUFFIX="${PROFILE_SUFFIX//,/–}" -OUTDIR=${timestamp}-job${SLURM_JOB_ID}-MESSI_results-${PROFILE_SUFFIX} +#OUTDIR=${timestamp}-job${SLURM_JOB_ID}-MESSI_results-${PROFILE_SUFFIX} # Samplesheet should maybe have it in profile only and not with CLI as it overrides it +#SAMPLESHEET="data/samplesheet_simulated_full.csv" +#SAMPLESHEET=data/samplesheet_bulk.csv +#SAMPLESHEET="data/samplesheet_multimodal.csv" +#SAMPLESHEET="data/samplesheet_covid_only.csv" +#SAMPLESHEET="data/samplesheet_htx_only.csv" +#SAMPLESHEET="data/samplesheet_all.csv" +SAMPLESHEET="data/samplesheet_momix_full.csv" +#SAMPLESHEET=data/samplesheet_momix_3.csv #SAMPLESHEET=data/samplesheet_simulated_data.csv #SAMPLESHEET=data/samplesheet_simulated_part1.csv #SAMPLESHEET=data/samplesheet_simulated_part2.csv -SAMPLESHEET=data/samplesheet_simulated_part3.csv +#SAMPLESHEET=data/samplesheet_simulated_part3.csv #SAMPLESHEET=data/samplesheet_feat_selection.csv #SAMPLESHEET=data/samplesheet_test_full.csv #SAMPLESHEET=data/samplesheet_test_small.csv #SAMPLESHEET=data/samplesheet_325-405.csv #echo "Running pipeline with ${NXF_SRC_MAIN}" -#echo "Running data under '${SAMPLESHEET}'" +echo "Running data under '${SAMPLESHEET}'" +SAMPLESHEET_NAME="${SAMPLESHEET##*/}" +SAMPLESHEET_NAME="${SAMPLESHEET_NAME%.*}" +OUTDIR="${timestamp}-job${SLURM_JOB_ID}-MESSI_results-${PROFILE_SUFFIX}-${SAMPLESHEET_NAME}" # ============================================================================= # 4. Run the pipeline on the work dir +T=true +F=false + +OUTCOME="survival" + # The ansi-log option is used for redirecting output nextflow run ${NXF_SRC_MAIN} \ -profile ${PROFILE} \ --outdir ${OUTDIR} \ --samplesheet ${SAMPLESHEET} \ - -ansi-log false + -ansi-log false \ + -resume # ============================================================================= # 5. Compress output results and move back to submitted directory #cd ${TMPDIR} diff --git a/launcher_local.sh b/launcher_local.sh new file mode 100755 index 0000000..5aed1b6 --- /dev/null +++ b/launcher_local.sh @@ -0,0 +1,119 @@ +#!/bin/bash + +# ============================================================================ +# This is a wrapper script that triggers a SLURM BATCH script, hence most +# logics are in the other script. Here is more for parsing command line args +# +# Author: Tony Liang +# ============================================================================ + + +# Load modules? +#module load apptainer/1.3.4 +#module load java/17.0.6 +#module load nextflow/24.10.2 +module load apptainer/1.3.5 +module load java/21.0.1 +module load nextflow/24.10.2 + +# 1) gcccore/.12.3 (H) 5) libfabric/1.18.0 9) flexiblas/3.3.1 13) java/21 -> java/21.0.1 +# 2) gcc/12.3 (t) 6) pmix/4.2.4 10) imkl/2023.2.0 14) nextflow/24.10.2 +# 3) hwloc/2.9.1 7) ucc/1.2.0 11) CVMFS_CC/2023 +# 4) ucx/1.14.1 8) openmpi/4.1.5 (m) 12) apptainer/1.3.5 + + + +# pick a writable dir — adjust to your HPC's scratch/work path +export NXF_HOME=$SCRATCH_PATH/.nextflow + +# keep these off /home too +export NXF_WORK=./work +export NXF_TEMP=$SCRATCH_PATH/nxf_tmp +export NXF_PLUGINS_DIR=$NXF_HOME/nxf_plugins +export CAPSULE_CACHE_DIR=$SCRATCH_PATH/.capsule # Java/Capsule launcher cache + +# Java itself may also try to write to $HOME; give it a writable tmp +export JAVA_TOOL_OPTIONS="-Djava.io.tmpdir=/scratch/tliang19/tmp" + +# Locates dir of this script +SCRIPT_DIR=$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd ) +LAUNCHER_SCRIPT=${SCRIPT_DIR}/launch_MESSI_pipeline.sh +# The hidden env file +ENV_FILE=${SCRIPT_DIR}/.env +# Check if file exists or not +if [ ! -f ${ENV_FILE} ]; then + echo "missing .env file" + exit 1 +else + set -o allexport + source .env # This sources the .env file + set +o allexport + if [ "${ALLOCATION_CODE}" = "REPLACE" ] || [ "${MAIL_USER}" = "REPLACE" ]; then + echo -e "\nERROR: Did not change ALLOCATION_CODE or MAIL_USER\n" + exit 1 + fi +fi + + + + +# Then call the pipeline here +export NXF_OFFLINE='true' +#export NXF_PLUGINS_DIR=./nxf_plugins +#export NXF_OFFLINE=true +export NXF_DISABLE_CHECK_LATEST='true' + + +#export NXF_TEMP="/data1/tliang19/tools/nextflow/nf-tmp" +#SELECT_FEAT='true' # Or use 'false' +SELECT_FEAT='false' + +#CSV_FILE="data/samplesheet_momix.csv" +CSV_FILE="data/samplesheet_momix_3.csv" +#CSV_FILE="data/local_samplesheet.csv" +#CSV_FILE="data/samplesheet_bulk.csv" +#CSV_FILE="data/samplesheet_multimodal.csv" +#CSV_FILE="data/samplesheet_htx_only.csv" +#CSV_FILE="data/samplesheet_covid_only.csv" +#CSV_FILE="data/all_datasets.csv" +#CSV_FILE="data/samplesheet_test_small.csv" +F="false" +T="true" + +FILTER="1" # True in nextflow +#FILTER="0" +OUTCOME_TYPE="classification" +#OUTCOME_TYPE="survival" +K_FOLD_NUMBER=5 +#K_FOLD_NUMBER=5 +#SK_MODS="Logit,XGBoost" +#SK_MODS="Logit,MLP" + +#MAX_MEM="4G" +MAX_MEM="6G" + + +#PROFILE="standard,docker,test" +PROFILE="arc_local,apptainer" +# Specify number of CPUs and memory +nextflow run main.nf \ + -profile $PROFILE \ + -resume \ + --single_modality_mode $F \ + --samplesheet ${CSV_FILE} \ + --filter_low_var $FILTER \ + --outdir results \ + --skip_sksurv $T \ + --skip_rgcca $T \ + --skip_mofa $T \ + --skip_sklearn $T \ + --skip_caret_multimodal $T \ + --skip_diablo $T \ + --skip_cplr $F \ + --skip_mogonet $T \ + --skip_integrao $T \ + --pipeline_dir ./ \ + --k_fold_number $K_FOLD_NUMBER \ + --selectFeature $F \ + --outcome_type $OUTCOME_TYPE \ + --publish_relevant $F diff --git a/modules/calculate_metrics/main.nf b/modules/calculate_metrics/main.nf index dfc71c6..fbf1eb0 100644 --- a/modules/calculate_metrics/main.nf +++ b/modules/calculate_metrics/main.nf @@ -1,7 +1,7 @@ process CALCULATE_METRICS { debug true label 'process_single' - label 'mogonet' // Should be a generic python container? + label 'sksurv' // This has generic python , sklearn and sksurv modules publishDir ( path: "${params.outdir}/${task.process.tokenize(':').join('/').toLowerCase()}", mode: 'copy', @@ -12,16 +12,26 @@ process CALCULATE_METRICS { // Input output blocks input: tuple val(method_name), path(result_table) - val(threshold) + val(threshold) // This gets ignored for classification + val(outcome_type) // One of classification or survival output: tuple val(method_name), path('*.csv'), emit: metric_table path('*log*'), emit: log script: - def script_name = "calculate_metrics.py" // Execute this script in resources/usr/bin + //def script_name = "calculate_metrics.py" // Execute this script in resources/usr/bin + if (outcome_type == "classification") """ - ${script_name} \ + calculate_classification_metrics.py \ --result_path=${result_table} \ --threshold=${threshold} > \ ${task.process.tokenize(':')[-1].toLowerCase()}.log """ + else if (outcome_type == "survival") + """ + calculate_survival_metrics.py \ + --result_path=${result_table} > \ + ${task.process.tokenize(':')[-1].toLowerCase()}.log + """ + else + error "Did not implement other outcome types than 'classification' or 'survival'" } diff --git a/modules/calculate_metrics/resources/usr/bin/calculate_metrics.py b/modules/calculate_metrics/resources/usr/bin/calculate_classification_metrics.py similarity index 100% rename from modules/calculate_metrics/resources/usr/bin/calculate_metrics.py rename to modules/calculate_metrics/resources/usr/bin/calculate_classification_metrics.py diff --git a/modules/calculate_metrics/resources/usr/bin/calculate_survival_metrics.py b/modules/calculate_metrics/resources/usr/bin/calculate_survival_metrics.py new file mode 100755 index 0000000..02ddd98 --- /dev/null +++ b/modules/calculate_metrics/resources/usr/bin/calculate_survival_metrics.py @@ -0,0 +1,232 @@ +#!/usr/bin/env python +""" +This script is used to calculate survival metrics from merged out-of-fold (OOF) +results of all compared methods. Each row is one sample predicted while its fold +was the test set, with columns: + + method_name, dataset, fold identifiers + time, status observed follow-up time and event indicator (1 = event) + lp risk score (higher = higher risk, e.g. Cox linear predictor) + surv_ predicted survival probability S(t) at horizon t + +Metrics are computed for each fold (training set = the other folds of the same +method + dataset, used for the IPCW censoring weights) and then summarized as +mean / sd across folds (`_mean`, `_std`). Probability-based +metrics are also computed on the pooled OOF predictions (`_pooled`). +The output is one row per method + dataset. C-index on `lp` is NOT pooled because risk scores from different fold +models are on different scales and should not be ranked against each other. +Instead, `c_harrell_s` / `c_uno_s` use 1 - S(t) at --c_horizon as the risk +score, which is a probability and therefore comparable across fold models. + +`ipa` (index of prediction accuracy) = 1 - IBS / IBS of a null model that gives +everyone the Kaplan-Meier curve of the training set. 0 = no better than the null +model, so unlike the raw IBS it is comparable across datasets. + +Author: Tony Liang + +Usage: + calculate_survival_metrics.py [options] + +Options: + --result_path=RES_PATH Path to the merged result table [default: empty] + --out_path=OUT_PATH Path of the output summary table [default: metrics_summary.csv] + --round_digit=DIGITS Number of digits to round metrics [default: 3] + --c_horizon=T Horizon t of the 1 - S(t) based C-index [default: 730] +""" + +# Import libraries +import re + +from docopt import docopt +import pandas as pd +import numpy as np +from sksurv.util import Surv +from sksurv.nonparametric import kaplan_meier_estimator +from sksurv.metrics import ( + concordance_index_censored, concordance_index_ipcw, + cumulative_dynamic_auc, brier_score, integrated_brier_score +) + + +ID_COLS = ['method_name', 'dataset', 'fold'] + + +def get_horizons(df): + """Parse evaluation horizons from the `surv_` column names.""" + cols = [c for c in df.columns if re.fullmatch(r'surv_\d+(\.\d+)?', c)] + cols = sorted(cols, key=lambda c: float(c.split('_')[1])) + return cols, np.array([float(c.split('_')[1]) for c in cols]) + + +def safe(fun, *args, **kwargs): + """Run a metric, return NaN (with a message) if it is not computable.""" + try: + return fun(*args, **kwargs) + except ValueError as e: + print(f" skipped {fun.__name__}: {e}") + return None + + + +# Helper to calculate metrics of one evaluation (one fold, or pooled) +def calculate_metrics(train, test, surv_cols, horizons, c_horizon, use_lp=True): + """ + Arguments: + + train Rows used to estimate the censoring distribution for IPCW + (the other folds; for pooled evaluation, all rows) + ng `surv_cols` + c_horizon Horizon t whose 1 - S(t) is used as a fold-comparable risk score + use_lp Whether to compute C-indices from `lp` + """ + y_train = Surv.from_arrays(event=train['status'].astype(bool), time=train['time']) + # IPCW needs the censoring KM of train at each test time, which sksurv cannot + # extrapolate past the last train time. Test rows after it are still at risk at + # every evaluated time (< au) + y_test = Surv.from_arrays(event=test['status'].astype(bool), + time=test['time'].clip(upper=train['time'].max())) + + metrics = {'n': len(test), 'n_events': int(test['status'].sum())} + c_s = f's{c_horizon:g}' + for name in ['c_harrell', 'c_uno', f'c_harrell_{c_s}', f'c_uno_{c_s}']: + metrics[name] = np.nan + for t in horizons: + metrics[f'auc_{t:g}'] = np.nan + for t in horizons: + metrics[f'brier_{t:g}'] = np.nan + metrics['mean_auc'] = np.nan + for name in ['ibs', 'ibs_null', 'ipa']: + metrics[name] = np.nan + # S(t) is NaN after the last training event of the fold model (no information) + surv_all = test[surv_cols].to_numpy() + predicted = ~np.isnan(surv_all).any(axis=0) + if not predicted.all(): + print(f" skipped horizons without S(t) prediction: {horizons[~predicted]}") + # IPCW metrics only work inside the follow-up range of both train and test + valid = predicted & (horizons >= test['time'].min()) \ + & (horizons < test['time'].max()) & (horizons < train['time'].max()) + times = horizons[valid] + surv = surv_all[:, valid] + tau = times.max() if len(times) else None + + # Discrimination of the overall risk score: `lp` (within one fold model only) + # and 1 - S(t) at c_horizon (comparable across fold models, so also pooled) + risk_scores = {} + if c_horizon in horizons[predicted]: + risk_scores[f'_{c_s}'] = 1 - test[f'surv_{c_horizon:g}'].to_numpy() + else: + print(f" skipped C-index on 1 - S({c_horizon:g}): no S(t) prediction") + if use_lp: + risk_scores = {'': test['lp'].to_numpy(), **risk_scores} + for suffix, risk in risk_scores.items(): + res = safe(concordance_index_censored, + test['status'].astype(bool), test['time'], risk) + if res is not None: + metrics[f'c_harrell{suffix}'] = res[0] + # Uno's C is truncated at tau to keep the IPCW weights stable + res = safe(concordance_index_ipcw, y_train, y_test, risk, tau=tau) + if res is not None: + metrics[f'c_uno{suffix}'] = res[0] + + if len(times) == 0: + print(" no valid horizon inside the follow-up range") + return metrics + + # Time-dependent (cumulative/dynamic) AUC, risk at t is 1 - S(t) + res = safe(cumulative_dynamic_auc, y_train, y_test, 1 - surv, times) + if res is not None: + for t, a in zip(times, res[0]): + metrics[f'auc_{t:g}'] = a + metrics['mean_auc'] = res[1] + # Integrated Brier score over the horizons (needs at least 2) + if len(times) >= 2: + res = safe(integrated_brier_score, y_train, y_test, surv, times) + if res is not None: + metrics['ibs'] = res + # Null model: Kaplan-Meier curve of the training set for everyone + km_t, km_s = kaplan_meier_estimator(train['status'].astype(bool), train['time']) + s_null = km_s[np.searchsorted(km_t, times, side='right') - 1] + res = safe(integrated_brier_score, y_train, y_test, + np.tile(s_null, (len(test), 1)), times) + if res is not None: + metrics['ibs_null'] = res + metrics['ipa'] = 1 - metrics['ibs'] / res + + return metrics + + + + # Main entrance of the script +def main(result_path, out_path, round_digit, c_horizon): + df = pd.read_csv(result_path) + surv_cols, horizons = get_horizons(df) + print(f"Horizons found: {horizons}") + if c_horizon not in horizons: + raise ValueError(f"--c_horizon={c_horizon:g} is not one of the horizons {horizons}") + + rows = [] + for (method, dataset), group in df.groupby(['method_name', 'dataset']): + folds = sorted(group['fold'].unique()) + print(f"\n{method} / {dataset}: {len(folds)} folds, {len(group)} samples") + # Keep only horizons predicted in every fold, so per-fold metrics (e.g. IBS + # range, Uno's tau) are comparable to each other and to the pooled ones + common = group[surv_cols].notna().all(axis=0).to_numpy() + if not common.all(): + print(f" horizons dropped (no S(t) in some fold): {horizons[~common]}") + group = group.drop(columns=np.array(surv_cols)[~common]) + g_cols, g_horizons = list(np.array(surv_cols)[common]), horizons[common] + if c_horizon not in g_horizons: + print(f" --c_horizon={c_horizon:g} dropped, no 1 - S(t) C-index") + # Per fold: censoring distribution from the training part of that fold + for fold in folds: + print(f" {fold}") + test = group[group['fold'] == fold] + train = group[group['fold'] != fold] + rows.append({'method_name': method, 'dataset': dataset, 'fold': fold, + **calculate_metrics(train, test, g_cols, g_horizons, c_horizon)}) + # Pooled OOF predictions, probability-based metrics only + print(" pooled") + rows.append({'method_name': method, 'dataset': dataset, 'fold': 'pooled', + **calculate_metrics(group, group, g_cols, g_horizons, c_horizon, + use_lp=False)}) + # Same columns for every group, even if some horizons were dropped + per_fold = pd.DataFrame(rows) + c_s = f's{c_horizon:g}' + per_fold = per_fold.reindex(columns=ID_COLS + [ + 'n', 'n_events', 'c_harrell', 'c_uno', f'c_harrell_{c_s}', f'c_uno_{c_s}', + *[f'auc_{t:g}' for t in horizons], *[f'brier_{t:g}' for t in horizons], + 'mean_auc', 'ibs', 'ibs_null', 'ipa']) + metric_cols = [c for c in per_fold.columns if c not in ID_COLS + ['n', 'n_events']] + # Mean and sd across folds, plus the pooled values + is_pooled = per_fold['fold'] == 'pooled' + grouped = per_fold[~is_pooled].groupby(['method_name', 'dataset']) + summary = grouped[metric_cols].agg(['mean', 'std']) + summary.columns = [f'{m}_{s}' for m, s in summary.columns] + summary.insert(0, 'n_folds', grouped.size()) + pooled = per_fold[is_pooled].set_index(['method_name', 'dataset']) + summary.insert(1, 'n', pooled['n']) + summary.insert(2, 'n_events', pooled['n_events']) + pooled = pooled[metric_cols] + pooled = pooled.dropna(axis=1, how='all').add_suffix('_pooled') + summary = summary.join(pooled).reset_index() + + # ===================================================== + # Write out the summary dataframe + print(f"\nSaving metrics to '{out_path}'") + summary.round(round_digit).to_csv(out_path, index=False) + + return summary + + + +# Execute the fun here +if __name__ == '__main__': + # Parse docopt + args = docopt(__doc__) + # Execute runner + main( + result_path=args['--result_path'], + out_path=args['--out_path'], + round_digit=int(args['--round_digit']), + c_horizon=float(args['--c_horizon']) + ) diff --git a/modules/caret_multimodal/select_feature/main.nf b/modules/caret_multimodal/select_feature/main.nf index 85171bb..7850c85 100755 --- a/modules/caret_multimodal/select_feature/main.nf +++ b/modules/caret_multimodal/select_feature/main.nf @@ -20,6 +20,7 @@ process CARET_MULTIMODAL_SELECT_FEATURE { tuple val(dataset_name), path(data_path) output: path("*.csv"), emit: features + path("*.json"), emit: hyperparameters path("*plot*"), optional: true, emit: plot path("*.log"), emit: log diff --git a/modules/caret_multimodal/select_feature/resources/usr/bin/caret_multimodal_select_feature.R b/modules/caret_multimodal/select_feature/resources/usr/bin/caret_multimodal_select_feature.R index 03d22d8..376581c 100755 --- a/modules/caret_multimodal/select_feature/resources/usr/bin/caret_multimodal_select_feature.R +++ b/modules/caret_multimodal/select_feature/resources/usr/bin/caret_multimodal_select_feature.R @@ -25,6 +25,7 @@ library(dplyr) library(magrittr) library(tibble) library(caret) +library(jsonlite) bin_dir <- Sys.getenv("PATH") |> strsplit(":") |> @@ -119,6 +120,8 @@ main <- function(mae_path, dataset_name, nfolds, criteria_order="coef") { # TODO: Run a cv on X and Y to get hyperparameters # Then fit the model cv_model <- run_caret_multimodal_cv(X, Y, nfolds = nfolds) + # Get tuned summary + tuning_summary <- as.data.frame(summary(cv_model)) # TODO: Extract the features out from your final model and wrangle to df for downstream usage feature_weights <- caretMultimodal:::compute_feature_contributions.caret_stack(cv_model, n_features = Inf) feats_df <- feature_weights %>% @@ -143,7 +146,64 @@ main <- function(mae_path, dataset_name, nfolds, criteria_order="coef") { output_format <- "csv" feats_file <- paste0(comb_name , "_", "features_selected", ".", output_format) write.csv(x=feats_df, file=feats_file, row.names=FALSE) + # And the hyperparams section as well + tuning_summary <- as.data.frame(summary(cv_model)) + + lambda_grid <- sort(unique(10^seq(-4, 3, length = 20))) + + hyperparameter_record <- list( + run_id = paste(method, dataset_name, sep = "-"), + dataset = dataset_name, + method = method, + analysis_stage = "model_selection", + + selection = list( + strategy = "cross_validated_late_fusion_stacking", + folds = nfolds, + + selected_model_summary = tuning_summary + ), + + parameters = list( + base_learner = list( + value = "glmnet", + treatment = "fixed" + ), + stack_learner = list( + value = "glmnet", + treatment = "fixed" + ), + alpha = list( + value = 0, + treatment = "fixed" + ), + lambda = list( + value = list( + candidate_values = lambda_grid, + selected_models = tuning_summary + ), + treatment = "tuned" + ), + seed = list( + value = seed, + treatment = "data-derived" + ) + ) + ) + + hyperparameter_file <- paste0( + method, "-", dataset_name, + "_selected_hyperparameters.json" + ) + jsonlite::write_json( + hyperparameter_record, + path = hyperparameter_file, + pretty = TRUE, + auto_unbox = TRUE, + dataframe = "rows", + na = "null" + ) return(feats_df) } diff --git a/modules/caret_multimodal/train/resources/usr/bin/caret_multimodal_train.R b/modules/caret_multimodal/train/resources/usr/bin/caret_multimodal_train.R index 042b94d..2f3f298 100755 --- a/modules/caret_multimodal/train/resources/usr/bin/caret_multimodal_train.R +++ b/modules/caret_multimodal/train/resources/usr/bin/caret_multimodal_train.R @@ -5,18 +5,16 @@ doc <- "This script is to run caret_multimodal method, run it only against the t processed data. Output is the path to this trained model and the test MAE. Usage: - caret_multimodal_train.R [options] + caret_multimodal_train.R [options] Options: - --fold_path=FOLD_PATH Path to read current test fold - --label=LABEL Label of id and fold of data [default: data] - --prefix=PREFIX Prefix to read HDF5 [default: pre] - --run_inner_cv Run inner cv with train data or not [default: false] - --method_name=METHOD_NAME Name of the method [default: caret_multimodal] + --fold_path=FOLD_PATH Path to read current test fold + --label=LABEL Label of id and fold of data [default: data] + --prefix=PREFIX Prefix to read HDF5 [default: pre] + --run_inner_cv Run inner cv with train data or not [default: false] + --method_name=METHOD_NAME Name of the method [default: caret_multimodal] " -# TODO: You could possibly add more args for the cli - # Load libraries library(dplyr) @@ -26,17 +24,17 @@ library(stringr) library(caret) # Load scripts ======================================================== # Gather the pipeline dir (THIS IS VERY UGGLY FIX) -bin_dir <- Sys.getenv("PATH") |> - strsplit(":") |> - unlist() |> - tail(1) +bin_dir <- Sys.getenv("PATH") |> + strsplit(":") |> + unlist() |> + tail(1) # Determine if running on cluster deploy mode or local mode is_scratch <- stringr::str_detect(bin_dir, pattern = "scratch") if (is_scratch) { - pipeline_dir <- gsub("/bin", "", bin_dir) + pipeline_dir <- gsub("/bin", "", bin_dir) } else { - pipeline_dir <- "" + pipeline_dir <- "" } @@ -49,107 +47,236 @@ load_utils(here(pipeline_dir, "bin/plotting")) # Parse docopt opt <- docopt::docopt(doc) -# TODO: implement your training logic here -train_model <- function(train_data, inner_cv=FALSE) { - #alphas <- c(0.7, 0.775, 0.850, 0.925, 1) - #lambdas <- seq(0.001, 0.1, by = 0.01) - alphas <- c(0) # Ridge regression only, no Lasso or elastic-net, since we have many features and want to keep them all - lambdas <- 10^seq(-4, 3, length = 20) # Use log space lambda values - tuneGrid <- expand.grid(alpha = alphas, lambda = lambdas) - # Caret relies on tuneGrid and trainControl for hyperparameter tuning - if (inner_cv) { - trControl <- trainControl( - method = "repeatedcv", - number = 5, - repeats = 3, - classProbs = TRUE, - summaryFunction = twoClassSummary, - savePredictions = "final" - ) - - } else { - # Default with no inner cv, simple train-test - trControl <- trainControl( - method = "cv", - number = 5, # This is internal cv for tuning the hyperparameters, not the same as the outer cv fold split, which is done by nextflow - classProbs = TRUE, - summaryFunction = twoClassSummary, - savePredictions = "final" - ) - } - # ========================================================= - # First fit individual models for each modality - base_models <- caretMultimodal::caret_list( - target = train_data$Y, - data_list = train_data$X, - method = "glmnet", - tuneGrid = tuneGrid, - trControl = trControl - ) - message("\nFinished fitting base models for each modality, now stacking them together...") - # Then fit the ensemble model - stack_model <- caretMultimodal::caret_stack( - caret_list = base_models, - method = "glmnet", - tuneGrid = tuneGrid, - trControl = trControl - ) - message("\nFinished fitting the stacked model.") - # Return the stacked model - return(stack_model) +# Stacking helpers ==================================================== + +# Save the objects needed to reproduce/debug a caret_stack failure. +# Returns the normalized path of the debug file. +.save_stack_debug <- function(base_models, target, tuneGrid, trControl, + error_msg, label) { + debug_file <- paste0(label, "-caret_stack_debug.rds") + saveRDS( + list( + base_models = base_models, + target = target, + tuneGrid = tuneGrid, + trControl = trControl, + rng_state = .Random.seed, + error = error_msg + ), + file = debug_file + ) + normalizePath(debug_file, mustWork = FALSE) +} + +# glmnet refuses exactly zero-variance predictors (C++ error 7777), which is +# what happens when the base models carry no signal: caret selects a +# degenerate intercept-only lambda, the OOF predictions become constant, and +# the stack glmnet crashes -> caret stops with "Stopping". +# +# caret_stack() recomputes its OOF features from the base models' stored +# `pred` data frames, so we nudge constant class-probability columns there by +# a tiny deterministic offset (+/- epsilon, no RNG involved). Only the +# probability columns are touched -- never alpha/lambda, which caret uses as +# join keys when extracting the best-tune predictions. +# +# Returns list(base_models, nudged_columns). +.nudge_constant_oof <- function(base_models, epsilon = 1e-6) { + nudged_cols <- character(0) + for (i in seq_along(base_models)) { + m <- base_models[[i]] + p <- m$pred + prob_cols <- intersect(m$levels, names(p)) + const_cols <- prob_cols[vapply( + prob_cols, + function(cl) length(unique(p[[cl]])) == 1, + logical(1) + )] + if (length(const_cols)) { + nudge <- rep(c(-1, 1), length.out = nrow(p)) * epsilon + for (cl in const_cols) { + p[[cl]] <- p[[cl]] + nudge + } + base_models[[i]]$pred <- p + nudged_cols <- c(nudged_cols, paste0(names(base_models)[i], ".", const_cols)) + } + } + list(base_models = base_models, nudged_columns = nudged_cols) +} + +# Fit the stacked model with a safety net for no-signal data: +# 1st attempt: plain caret_stack (the normal path, unchanged). +# On error: save debug objects, nudge constant OOF columns, retry +# caret_stack itself, and tag the result with a +# `stack_fallback` attribute so retried folds can be +# identified downstream. +# If the retry also fails, stop with a message pointing to the debug file. +.stack_with_fallback <- function(base_models, tuneGrid, trControl, target, label) { + stack_model <- tryCatch( + { + caretMultimodal::caret_stack( + caret_list = base_models, + method = "glmnet", + tuneGrid = tuneGrid, + trControl = trControl + ) + }, + error = function(e) { + debug_file <- .save_stack_debug( + base_models = base_models, + target = target, + tuneGrid = tuneGrid, + trControl = trControl, + error_msg = conditionMessage(e), + label = label + ) + + nudged <- .nudge_constant_oof(base_models) + + retried <- tryCatch( + { + caretMultimodal::caret_stack( + caret_list = nudged$base_models, + method = "glmnet", + tuneGrid = tuneGrid, + trControl = trControl + ) + }, + error = function(e2) { + stop( + "Stacking failed again after nudging constant OOF predictions: ", + conditionMessage(e2), + "\nDebug objects saved to: ", debug_file, + call. = FALSE + ) + } + ) + + attr(retried, "stack_fallback") <- list( + reason = conditionMessage(e), + nudged_columns = nudged$nudged_columns, + debug_file = debug_file + ) + + warning( + "\nStacking failed: ", conditionMessage(e), + "\nRetried caret_stack on nudged (near-constant) OOF predictions.", + "\nDebug objects saved to: ", debug_file + ) + + retried + } + ) + stack_model +} + +# Training ============================================================ + +train_model <- function(train_data, label = "data", inner_cv = FALSE) { + #alphas <- c(0.7, 0.775, 0.850, 0.925, 1) + #lambdas <- seq(0.001, 0.1, by = 0.01) + alphas <- c(0) # Ridge regression only, no Lasso or elastic-net, since we have many features and want to keep them all + lambdas <- 10^seq(-4, 3, length = 20) # Use log space lambda values + tuneGrid <- expand.grid(alpha = alphas, lambda = lambdas) + # Caret relies on tuneGrid and trainControl for hyperparameter tuning + if (inner_cv) { + trControl <- trainControl( + method = "repeatedcv", + number = 5, + repeats = 3, + classProbs = TRUE, + summaryFunction = twoClassSummary, + savePredictions = "final" + ) + + } else { + # Default with no inner cv, simple train-test + trControl <- trainControl( + method = "cv", + number = 5, # This is internal cv for tuning the hyperparameters, not the same as the outer cv fold split, which is done by nextflow + classProbs = TRUE, + summaryFunction = twoClassSummary, + savePredictions = "final" + ) + } + + + + + # ========================================================= + # First fit individual models for each modality + base_models <- caretMultimodal::caret_list( + target = train_data$Y, + data_list = train_data$X, + method = "glmnet", + tuneGrid = tuneGrid, + trControl = trControl + ) + message("\nFinished fitting base models for each modality, now stacking them together...") + + # Then fit the ensemble model (with a retry safety net for no-signal data) + stack_model <- .stack_with_fallback( + base_models = base_models, + tuneGrid = tuneGrid, + trControl = trControl, + target = train_data$Y, + label = label + ) + + message("\nFinished fitting the stacked model.") + return(stack_model) } get_seed <- function(dataset_name) { - d_int <- utf8ToInt(dataset_name) # Convert dataset name to integer - # For caretMultimodal only, add a "hacky" constant to the seed to make it different from other methods - # As it could go into problem with internal cv splitting when parallized run - seed <- sum(d_int) + 1 - message("\nSeed used: ", seed) - return(seed) + d_int <- utf8ToInt(dataset_name) # Convert dataset name to integer + # For caretMultimodal only, add a "hacky" constant to the seed to make it different from other methods + # As it could go into problem with internal cv splitting when parallized run + seed <- sum(d_int) + 1 + message("\nSeed used: ", seed) + return(seed) } # Main function to run main <- function(fold_path, label, prefix, method_name, inner_cv=FALSE) { - seed <- get_seed(label) # Set seed based on dataset name for reproducibility - set.seed(seed) - # Log the params used - args_used <- c(as.list(environment())) - logging_params(args_used) - cat("\nLooking at this fold: ", fold_path, "\n") - # Finds the subdirectories containing tr and te MAE - d <- list.files(path=fold_path, full.names = TRUE) - train_path <- d[str_detect(d, pattern = "_tr")] - test_path <- d[str_detect(d, pattern = "_te")] - # Then should read in the MAE and convert it to list of X and Y - # This converts usual MAE format of p_i x N to N x p_i - train_data <- load_MAE(train_path, prefix="train") |> extract_Xy(verbose_target=TRUE) - test_data <- load_MAE(test_path, prefix="test") |> extract_Xy(verbose_target=TRUE) - sample_names <- check_common_samples(train_data) - # Sanity check for debugging - cat("\nTotal of", length(sample_names), "samples:\n", sample_names) - # Train the model here - model <- train_model(train_data, inner_cv=FALSE) - if (is.null(model)) { - stop("Model did not implement yet, model is NULL") - } - - # Filenames to write out - model_file <- paste(label, paste(method_name, "model.rds", sep="_"), sep="-") - test_file <- paste(label, "test_data.rds", sep="-") - cat("\nSaving files to", label, "\n") - # Write out to disk - saveRDS(object = model, model_file) - saveRDS(object = test_data, test_file) - return(model) + seed <- get_seed(label) # Set seed based on dataset name for reproducibility + set.seed(seed) + # Log the params used + args_used <- c(as.list(environment())) + logging_params(args_used) + cat("\nLooking at this fold: ", fold_path, "\n") + # Finds the subdirectories containing tr and te MAE + d <- list.files(path=fold_path, full.names = TRUE) + train_path <- d[str_detect(d, pattern = "_tr")] + test_path <- d[str_detect(d, pattern = "_te")] + # Then should read in the MAE and convert it to list of X and Y + # This converts usual MAE format of p_i x N to N x p_i + train_data <- load_MAE(train_path, prefix="train") |> extract_Xy(verbose_target=TRUE) + test_data <- load_MAE(test_path, prefix="test") |> extract_Xy(verbose_target=TRUE) + sample_names <- check_common_samples(train_data) + # Sanity check for debugging + cat("\nTotal of", length(sample_names), "samples:\n", sample_names) + # Train the model here + model <- train_model(train_data, label = label, inner_cv = inner_cv) + if (is.null(model)) { + stop("Model did not implement yet, model is NULL") + } + + # Filenames to write out + model_file <- paste(label, paste(method_name, "model.rds", sep = "_"), sep = "-") + test_file <- paste(label, "test_data.rds", sep = "-") + cat("\nSaving files to", label, "\n") + # Write out to disk + saveRDS(object = model, model_file) + saveRDS(object = test_data, test_file) + return(model) } # Call the function here main( - fold_path = opt$fold_path, - label = opt$label, - prefix = opt$prefix, - method_name = opt$method_name, - inner_cv = as.logical(opt$run_inner_cv) + fold_path = opt$fold_path, + label = opt$label, + prefix = opt$prefix, + method_name = opt$method_name, + inner_cv = as.logical(opt$run_inner_cv) ) -cat("Done") \ No newline at end of file +cat("Done") diff --git a/modules/cooperative_learning/predict/main.nf b/modules/cooperative_learning/predict/main.nf index dd8724c..d5bd858 100644 --- a/modules/cooperative_learning/predict/main.nf +++ b/modules/cooperative_learning/predict/main.nf @@ -22,6 +22,7 @@ process COOPERATIVE_LEARNING_PREDICT { tuple val(dataset_name), val(fold_name), path(model) tuple val(dataset_name), val(fold_name), path(test_path) val(method_name) + val(outcome_type) output: tuple val(dataset_name), val(fold_name), val(method_name), path("*result*"), emit: result_table path('*log*'), optional: true, emit: log @@ -32,6 +33,7 @@ process COOPERATIVE_LEARNING_PREDICT { --model=${model} \ --test_path=${test_path} \ --label=${data_label} \ + --outcome_type=${outcome_type} \ --method_name=${method_name} > \ ${data_label}-${getPublishPath(task.process).tokenize('/')[-1]}.log diff --git a/modules/cooperative_learning/predict/resources/usr/bin/predict_cooperative_learning.R b/modules/cooperative_learning/predict/resources/usr/bin/predict_cooperative_learning.R index 4b5265b..9086068 100755 --- a/modules/cooperative_learning/predict/resources/usr/bin/predict_cooperative_learning.R +++ b/modules/cooperative_learning/predict/resources/usr/bin/predict_cooperative_learning.R @@ -1,10 +1,12 @@ #!/usr/bin/env Rscript -# Script to predict cooperative learning (simulate now) -doc <- "This script is to make predictions on test data of particular fold, +# Script to predict cooperative learning +doc <- "This script is to make predictions on test data of particular fold, using a model trained with cooperative learning method from multiview package -Output type is a path containing the predicted probabilities +classification: output table has the predicted probabilities (phat) +survival: output table has the risk score (lp) and S(t) at the horizons, + same format as the Python sksurv methods Usage: predict_cooperative_learning.R [options] @@ -15,6 +17,7 @@ Options: --label=LABEL Label of id and fold of data [default: data-fold_i] --method_name=METHOD Method name input from upstream [default: empty] --output_ext=EXT Extension of output table to save [default: csv] + --outcome_type=TYPE classification or survival [default: classification] " library(multiview) library(magrittr) @@ -29,24 +32,19 @@ pipeline_dir <- gsub("/bin", "", bin_dir) source(here::here(pipeline_dir, "bin/wrangling/get_result_table.R")) # Parse docopt opt <- docopt::docopt(doc) -# Default to use AveragedPredict and max.dist -main <- function(model_path, test_path, label, output_ext, method_name, - s=0.005, type="response", digit=3) { - if (method_name == "empty") { - stop("You did not provide method name") - } - # Load test data (from same fold test portion) - test_data <- readRDS(test_path) - cat("\nRead model from", model_path, "\n") - cat("\nRead test data from", test_path, "\n") + + +# Classification: predicted probability of the positive class +predict_classification <- function(model_path, test_path, label, method_name, type="response", digit=3,s=0.005) { # Load model (from same fold train portion) model <- readRDS(model_path) + test_data <- readRDS(test_path) # TODO: NEED a better way to handle this # Check if model is of cv object or not # When its cv - if ("cv.multivew" %in% class(model)) { - cat("\nReceived internal CV model, using lambda.1se") + if ("cv.multiview" %in% class(model)) { + message("\nReceived internal CV model, using lambda.1se") s <- "lambda.1se" } # Predict and get result @@ -59,13 +57,88 @@ main <- function(model_path, test_path, label, output_ext, method_name, method_name=method_name, test_data=test_data, digit=digit ) + return(result_table) +} + +# Survival: risk score and S(t) at the horizons from the multiview cox model. +# To be comparable with the sksurv methods (sksurv_predict.py), the output: +# - lp is the linear predictor, higher = riskier +# - S(t) at the same horizons (days), NA after the last train event +# (both handled by baseline_surv from run_cooperative_learning.R) +# - has the same columns: sample_name, lp, surv_, time, status, +# method_name, dataset, fold + +predict_survival <- function(model_path, test_path, label, method_name, digit=3, type="link") { + # Load model (from same fold train portion) + model <- readRDS(model_path) + test_data <- readRDS(test_path) + message("\nRead model from", model_path, "\n") + message("\nRead test data from", test_path, "\n") + + # Risk score, lambda chosen by inner CV in run_cooperative_learning.R + lp <- predict(model$cvfit, newx = test_data$X, s = "lambda.min", + type = type) |> as.numeric() + # Cox model: S(t | x) = S0(t) ^ exp(lp), S0 estimated on the train fold + surv <- outer(exp(lp), model$baseline_surv, function(r, s0) s0 ^ r) + # Result table, dataset and fold parsed from the label (dataset-..fold_i..) + result_table <- data.frame(sample_name = rownames(test_data$X[[1]]), + lp = lp) + for (h in names(model$baseline_surv)) { + result_table[[paste0("surv_", h)]] <- surv[, h] + } + result_table$time <- test_data$Y$time + result_table$status <- test_data$Y$status + result_table$method_name <- method_name + result_table$dataset <- sub("-[^-]*fold_[0-9]+.*$", "", label) + result_table$fold <- regmatches(label, regexpr("fold_[0-9]+", label)) + return(result_table) + +} +# Default to use AveragedPredict and max.dist +main <- function(model_path, test_path, label, output_ext, method_name, + s=0.005, outcome_type="classification", type="response", digit=3) { + + if (method_name == "empty") { + stop("You did not provide method name") + } + # Load test data (from same fold test portion) + + message("\nRead model from", model_path, "\n") + message("\nRead test data from", test_path, "\n") + if (outcome_type == "classification") { + message("\nPredicting classification\n") + + result_table <- predict_classification( + model_path = model_path, + test_path = test_path, + label = label, + method_name = method_name, + type = type, + digit = digit, + s = s + ) + } else if (outcome_type == "survival") { + message("\nPredicting survival\n") + + # For survival let s = lambda.min + result_table <- predict_survival( + model_path = model_path, + test_path = test_path, + label = label, + method_name = method_name, + digit = digit + ) + + } else { + stop("Unknown outcome type: ", outcome_type) + } # Write to files - cat("\nSaving as", output_ext, "format\n") + message("\nSaving as", output_ext, "format\n") result_file <- paste(label, paste0("result_table", ".", output_ext), sep="-") # Save to disk write.csv(result_table, result_file, row.names = FALSE) - return(pred_probs) + return(result_table) } @@ -74,7 +147,8 @@ main(model_path=opt$model_path, test_path=opt$test_path, label=opt$label, output_ext=opt$output_ext, - method_name=opt$method_name + method_name=opt$method_name, + outcome_type=opt$outcome_type ) -cat("Done") \ No newline at end of file +message("Done") diff --git a/modules/cooperative_learning/preprocess/main.nf b/modules/cooperative_learning/preprocess/main.nf index b405711..52ee889 100644 --- a/modules/cooperative_learning/preprocess/main.nf +++ b/modules/cooperative_learning/preprocess/main.nf @@ -20,6 +20,7 @@ process COOPERATIVE_LEARNING_PREPROCESS { /*Input and output blocks*/ input: tuple val(dataset_name), path(mae_path), path(split_dir) + val(outcome_type) output: /* Series of folder, each represent a fold, whereas within fold there are @@ -38,6 +39,7 @@ process COOPERATIVE_LEARNING_PREPROCESS { --mae_path=${mae_path} \ --split_dir=${split_dir} \ --dataset_name=${dataset_name} \ + --outcome_type=${outcome_type} \ --transpose > \ ${dataset_name}-${getPublishPath(task.process).tokenize('/')[-1].toLowerCase()}.log """ diff --git a/modules/cooperative_learning/select_feature/main.nf b/modules/cooperative_learning/select_feature/main.nf index cddec60..a6e047f 100644 --- a/modules/cooperative_learning/select_feature/main.nf +++ b/modules/cooperative_learning/select_feature/main.nf @@ -27,6 +27,7 @@ process COOPERATIVE_LEARNING_SELECT_FEATURE { tuple val(dataset_name), path(mae_path) output: path("*.csv"), emit: features + path("*.json"), emit: hyperparameters path("*plot*"), emit: plot path("*.log"), emit: log diff --git a/modules/cooperative_learning/select_feature/resources/usr/bin/cplr_select_features.R b/modules/cooperative_learning/select_feature/resources/usr/bin/cplr_select_features.R index 98b045d..d8df81d 100755 --- a/modules/cooperative_learning/select_feature/resources/usr/bin/cplr_select_features.R +++ b/modules/cooperative_learning/select_feature/resources/usr/bin/cplr_select_features.R @@ -62,6 +62,9 @@ nfolds=5, criteria_order="standardized_coef") { set.seed(seed) # PARAMS method <- "cooperative_learning" + model_family <- binomial() + lambda_rule <- "lambda.min" + # Log the params used args_used <- c(as.list(environment())) logging_params(args_used) @@ -81,12 +84,12 @@ nfolds=5, criteria_order="standardized_coef") { } # fit the model (but with cv) , note: their default nfolds is actual 10 but that takes long time cv_model <- cv.multiview(x_list = X, y = Y, - family = binomial(), type.measure=type.measure, + family = model_family, type.measure=type.measure, rho=rho, alpha=alpha, nfolds=nfolds) # Get the lambda to use - s <- cv_model$lambda.min + selected_lambda <- cv_model$lambda.min # Plot the loss vs different lambdas, the model under the hood is using glmnet par(mar=c(1,1,1,1)) @@ -96,7 +99,7 @@ nfolds=5, criteria_order="standardized_coef") { dev.off() # Now get those features out by taking top n percent of features in each view - feats_df <- coef_ordered(cv_model, s=s) %>% + feats_df <- coef_ordered(cv_model, s=selected_lambda) %>% as_tibble() %>% rename(feature=view_col) %>% mutate( @@ -119,6 +122,57 @@ nfolds=5, criteria_order="standardized_coef") { # write it to disk feats_file <- paste0(method, "-", dataset_name, "_", "features_selected", ".csv") write.csv(x=feats_df, file=feats_file, row.names=FALSE) + # =========================== + # And the hyperparameters part + hyperparameters <- list( + lambda = make_parameter_record( + value = as.numeric(selected_lambda), + treatment = "tuned" + ), + + rho = make_parameter_record( + value = as.numeric(rho), + treatment = "fixed" + ), + + alpha = make_parameter_record( + value = as.numeric(alpha), + treatment = "fixed" + ), + + family = make_parameter_record( + value = list( + family = model_family$family, + link = model_family$link + ), + treatment = "fixed" + ) + ) + + selection <- list( + strategy = "cv.multiview", + tuned_parameter = "lambda", + selected_rule = lambda_rule, + lambda_sequence = "automatically generated", + candidate_count = length(cv_model$lambda), + metric = type.measure, + validation = "kfold", + folds = as.integer(nfolds), + random_seed = as.integer(seed), + best_cv_error = min(cv_model$cvm, na.rm = TRUE) + ) + + hyperparameter_file <- write_selected_hyperparameters( + dataset_name = dataset_name, + method_name = method, + parameters = hyperparameters, + selection = selection + ) + + message( + "Selected hyperparameters written to: ", + hyperparameter_file + ) return(feats_df) } diff --git a/modules/cooperative_learning/train/main.nf b/modules/cooperative_learning/train/main.nf index 6fc9dff..2cdf2be 100644 --- a/modules/cooperative_learning/train/main.nf +++ b/modules/cooperative_learning/train/main.nf @@ -26,6 +26,7 @@ process COOPERATIVE_LEARNING_TRAIN { // Input, output blocks input: tuple val(dataset_name), path(mae_path), path(fold_path) + val(outcome_type) output: tuple val(dataset_name), val(fold_path.name), path('*model.rds'), emit: model tuple val(dataset_name), val(fold_path.name), path('*test_data.rds'), emit: test_data @@ -38,6 +39,7 @@ process COOPERATIVE_LEARNING_TRAIN { run_cooperative_learning.R \ --mae_path=${mae_path} \ --label=${data_label} \ + --outcome_type=${outcome_type} \ --fold_path=${fold_path} > \ ${data_label}-${getPublishPath(task.process).tokenize('/')[-1].toLowerCase()}.log """ diff --git a/modules/cooperative_learning/train/resources/usr/bin/run_cooperative_learning.R b/modules/cooperative_learning/train/resources/usr/bin/run_cooperative_learning.R index 5c014a0..531ca6a 100755 --- a/modules/cooperative_learning/train/resources/usr/bin/run_cooperative_learning.R +++ b/modules/cooperative_learning/train/resources/usr/bin/run_cooperative_learning.R @@ -1,10 +1,14 @@ #!/usr/bin/env Rscript # Script to run cooperative learning (simulate now) -doc <- "This script is to run cooperative learning method from multiview package, -train only. Tt could possibly be ran on a inner CV model, output is a model for +doc <- "This script is to run cooperative learning method from multiview package, +train only. Tt could possibly be ran on a inner CV model, output is a model for prediction usage in downstream. +classification: multiview binomial with default lambda path. +survival: multiview cox, lambda chosen by inner CV (cv.multiview) on the + train fold, plus the baseline survival S0(t) at the horizons. + Usage: run_cooperative_learning.R [options] @@ -16,6 +20,10 @@ Options: --prefix=PREFIX Prefix to read HDF5 [default: pre] --rho=RHO Weight on the agreement penalty, rho=0 is a form of early fusion, and rho=1 is a form of late fusion. [default: 0.5] --alpha=ALPHA Elastic net mixing param with 0 <= alpha <=1, when 1 = lasso, when 0 ridge. [default: 1] + --outcome_type=TYPE classification or survival [default: classification] + --time_col=TIME_COL colData column of survival time [default: time] + --status_col=STAT_COL colData column of event indicator [default: status] + --horizons=HORIZONS Comma separated days to predict S(t) at [default: 365,548,730,1095,1461,1826] " # Parase docopt opt <- docopt::docopt(doc) @@ -47,8 +55,72 @@ get_seed <- function(dataset_name) { } +run_cplr_classification <- function(train_data, rho, alpha) { + # use default settings + model <- multiview(x_list=train_data$X, + y=train_data$Y, + rho=rho, + family=binomial(), + alpha=alpha + ) + return(model) +} + +# Breslow estimate of the baseline survival S0(t) = S(t | lp = 0) at the horizons: +# H0(t) = sum over event times t_i <= t of d_i / sum_{j: t_j >= t_i} exp(lp_j) +# S0(t) = exp(-H0(t)) +# NA after the last event: the curve is flat there, as in sksurv_predict.py +breslow_baseline <- function(lp, time, status, horizons) { + event_times <- sort(unique(time[status == 1])) + risk <- exp(lp) + dH <- sapply(event_times, function(t) sum(status[time == t]) / sum(risk[time >= t])) + H0 <- cumsum(dH) + idx <- findInterval(horizons, event_times) + baseline_surv <- exp(-ifelse(idx == 0, 0, H0[pmax(idx, 1)])) + baseline_surv[horizons > max(event_times)] <- NA + names(baseline_surv) <- horizons + return(baseline_surv) +} + + + + +run_cplr_survival <- function(train_data, rho, alpha, horizons) { + train_x <- train_data$X + train_y <- train_data$Y + # glmnet cox rejects time <= 0 (e.g. TCGA patients with 0 days of follow-up): + # for fitting only, move them to half of the smallest positive time + train_time <- train_y$time + if (any(train_time <= 0)) { + eps <- min(train_time[train_time > 0]) / 2 + message("\n", sum(train_time <= 0), " train samples with time <= 0 set to ", eps) + train_time <- pmax(train_time, eps) + } + # lambda chosen by CV C-index: with rho > 0 the CV partial likelihood deviance + # of multiview cox increases from the first lambda on, so lambda.min would + # always be the all-zero model. S(t) stays calibrated since S0 is refitted below + model <- cv.multiview(x_list = train_x, + y = as.matrix(train_y), + family = "cox", rho = rho, alpha = alpha, + type.measure = "C") + coefs <- coef(model, s = "lambda.min") + cat("\nFitted model,", sum(coefs != 0), "non-zero coefficients at lambda.min\n") + + # Baseline survival from the train fold, predict then only needs + # S(t | x) = S0(t) ^ exp(lp) + train_lp <- as.numeric(predict(model, newx = train_x, s = "lambda.min", type = "link")) + baseline_surv <- breslow_baseline(train_lp, train_time, train_y$status, horizons) + message("\nBaseline S0(t) at ", paste(horizons, collapse = ", "), ": ", + paste(round(baseline_surv, 3), collapse = ", ")) + + model <- list(outcome_type = "survival", cvfit = model, baseline_surv = baseline_surv) + return(model) +} + + + # This the main function to execute -main <- function(mae_path, label, fold_path, inner_cv, prefix, rho, alpha) { +main <- function(mae_path, label, fold_path, inner_cv, prefix, rho, alpha, outcome_type="classification") { seed <- get_seed(label) # Set seed based on dataset name set.seed(seed) args_used <- c(as.list(environment())) @@ -58,8 +130,8 @@ main <- function(mae_path, label, fold_path, inner_cv, prefix, rho, alpha) { train_path <- d[str_detect(d, pattern = "_tr")] test_path <- d[str_detect(d, pattern = "_te")] # Then should read in the MAE and convert it to list of X and Y - train_data <- load_MAE(train_path, prefix="train") |> extract_Xy() - test_data <- load_MAE(test_path, prefix="test") |> extract_Xy() + train_data <- load_MAE(train_path, prefix="train") |> extract_Xy(outcome_type=outcome_type) + test_data <- load_MAE(test_path, prefix="test") |> extract_Xy(outcome_type=outcome_type) sample_names <- check_common_samples(train_data) cat("\nTotal of", length(sample_names), "samples:\n", sample_names) # Train a model to run inner cv or not @@ -69,13 +141,15 @@ main <- function(mae_path, label, fold_path, inner_cv, prefix, rho, alpha) { model <- "A" } else { cat("\nNot running inner cv per fold\n") - # use default settings - model <- multiview(x_list=train_data$X, - y=train_data$Y, - rho=rho, - family=binomial(), - alpha=alpha - ) + if (outcome_type == "survival") { + horizons <- c(365,548,730,1095,1461,1826) # Fixed for now + model <- run_cplr_survival(train_data, rho=rho, alpha=alpha, horizons=horizons) + } else if (outcome_type == "classification") { + model <- run_cplr_classification(train_data, rho=rho, alpha=alpha) + } else { + stop("Outcome type must be either 'classification' or 'survival'") + } + cat("\nFitted model\n") } # Files names to write out @@ -91,11 +165,12 @@ main <- function(mae_path, label, fold_path, inner_cv, prefix, rho, alpha) { } # Call the function here -main(mae_path = opt$mae_path, - label = opt$label, - fold_path = opt$fold_path, - inner_cv = opt$inner_cv, - prefix = opt$prefix, - rho = as.numeric(opt$rho), - alpha = as.numeric(opt$alpha) +main(mae_path = opt$mae_path, + label = opt$label, + fold_path = opt$fold_path, + inner_cv = opt$inner_cv, + prefix = opt$prefix, + rho = as.numeric(opt$rho), + alpha = as.numeric(opt$alpha), + outcome_type = opt$outcome_type ) \ No newline at end of file diff --git a/modules/diablo/select_feature/main.nf b/modules/diablo/select_feature/main.nf index 099a734..0e2750c 100644 --- a/modules/diablo/select_feature/main.nf +++ b/modules/diablo/select_feature/main.nf @@ -23,6 +23,7 @@ process DIABLO_SELECT_FEATURE { output: path("*.csv"), emit: features path("*plot*"), emit: plot, optional: true + path("*.json"), emit: hyperparameters path("*.log"), emit: log script: diff --git a/modules/diablo/select_feature/resources/usr/bin/diablo_select_features.R b/modules/diablo/select_feature/resources/usr/bin/diablo_select_features.R index 7866bcc..4290875 100755 --- a/modules/diablo/select_feature/resources/usr/bin/diablo_select_features.R +++ b/modules/diablo/select_feature/resources/usr/bin/diablo_select_features.R @@ -83,6 +83,15 @@ main <- function(mae_path, dataset_name, n_percent, design, ncomp) { logging_head_names(X=X, n = 10) # Then make the list of keepX, now defaults to take 100 percent of it keepX <- createKeepX(X, n_percent=n_percent, ncomp=ncomp) + + # Convert keepX into a JSON-friendly named list. + # Each list element represents one modality, and each value represents + # the number of retained features for one component. + keepX_json <- lapply( + keepX, + function(x) as.integer(unname(x)) + ) + cat("\nKeep X is the following:", "\n", unlist(keepX), "\n") # Then fit the model @@ -97,7 +106,7 @@ main <- function(mae_path, dataset_name, n_percent, design, ncomp) { } design_mat <- getDesign(X, corr = corr) - model <- block.splsda(X, Y, keepX=keepX, design=design_mat, ncomp=ncomp) + model <- block.splsda(X, Y, keepX=keepX, design=design_mat, ncomp=ncomp, scheme="horst") # Loop over components and extract features # Extract the features out from var_list and wrangle to df for downstream usage @@ -129,17 +138,64 @@ main <- function(mae_path, dataset_name, n_percent, design, ncomp) { feats_file <- paste0(comb_name , "_", "features_selected", ".", output_format) write.csv(x=feats_df, file=feats_file, row.names=FALSE) - # Plot the selected variables out - if (ncomp < 2) { - warning("Selected variables plot can only be shown when ncomp >= 1") - } else { - par(mar=c(1,1,1,1)) - getPlotDevice(name = "diablo_feature_selection_plot", dataset_name=comb_name, - height=8, width=8, device="svg") - plotVar(model) - dev.off() - - } + # ========================== + method_name <- paste(method, design, sep = "-") + # And the hyperparameters part + hyperparameters <- list( + ncomp = make_parameter_record( + value = as.integer(ncomp), + treatment = "fixed" + ), + + + keepX = make_parameter_record( + value = keepX_json, + treatment = "data-derived" + ), + + design = make_parameter_record( + value = design, + treatment = "fixed" + ), + + design_correlation = make_parameter_record( + value = as.numeric(corr), + treatment = "data-derived" + ), + + scheme = make_parameter_record( + value = "horst", + treatment = "fixed" + ) + ) + + selection <- list( + strategy = "none" + ) + + hyperparameter_file <- write_selected_hyperparameters( + dataset_name = dataset_name, + method_name = method_name, + parameters = hyperparameters, + selection = selection + ) + + message( + "Selected hyperparameters written to: ", + hyperparameter_file + ) + + # # Plot the selected variables out + # if (ncomp < 2) { + # warning("Selected variables plot can only be shown when ncomp >= 1") + # } else { + # par(mar=c(1,1,1,1)) + # getPlotDevice(name = "diablo_feature_selection_plot", dataset_name=comb_name, + # height=8, width=8, device="svg") + # plotVar(model) + # dev.off() + + # } return(feats_df) } diff --git a/modules/diablo/train/resources/usr/bin/run_diablo.R b/modules/diablo/train/resources/usr/bin/run_diablo.R index 5c00f61..3af1ab9 100755 --- a/modules/diablo/train/resources/usr/bin/run_diablo.R +++ b/modules/diablo/train/resources/usr/bin/run_diablo.R @@ -76,7 +76,8 @@ main <- function(mae_path, label, fold_path, design, ncomp, run_inner_cv, prefix sample_names <- check_common_samples(train_data) cat("\nTotal of", length(sample_names), "samples:\n", sample_names) - + # Set scheme to horst + scheme <- "horst" # Get the design matrix here if (design == "full") { corr <- 1 @@ -105,13 +106,14 @@ main <- function(mae_path, label, fold_path, design, ncomp, run_inner_cv, prefix model <- mixOmics::block.splsda(X = X, Y = Y, ncomp = tuned_output$ncomp, keepX = tuned_output$keepX, - design = design_mat) + design = design_mat, + scheme = scheme) #--------------------------------------------------------------------------- } else { cat("\nNot running inner cv per fold\n") # use default settings # Use a fully connected design on def , could also use null - model <- mixOmics::block.splsda(X = train_data$X, Y = train_data$Y, design=design_mat, ncomp=ncomp) + model <- mixOmics::block.splsda(X = train_data$X, Y = train_data$Y, design=design_mat, ncomp=ncomp, scheme=scheme) cat("\nFitted model\n") } diff --git a/modules/diablo/train/resources/usr/bin/tune_diablo.R b/modules/diablo/train/resources/usr/bin/tune_diablo.R index 12e6ea7..c2c1473 100755 --- a/modules/diablo/train/resources/usr/bin/tune_diablo.R +++ b/modules/diablo/train/resources/usr/bin/tune_diablo.R @@ -77,7 +77,8 @@ tune_diablo <- function(base_model, design, dist="centroids.dist", validation='M design = design, validation = validation, folds = folds, nrepeat = nrepeat, - dist = dist) + dist = dist, + scheme="horst") cat("\nFinished tuning keepX\n") list_keepX <- tune_features$choice.keepX diff --git a/modules/integrao/select_feature/main.nf b/modules/integrao/select_feature/main.nf index 734bb7a..a491d6d 100755 --- a/modules/integrao/select_feature/main.nf +++ b/modules/integrao/select_feature/main.nf @@ -21,6 +21,7 @@ process INTEGRAO_SELECT_FEATURE { tuple val(dataset_name), path(data_path) output: path("*.csv"), emit: features + path("*.json"), emit: hyperparameters path("*plot*"), optional: true, emit: plot path("*.log"), emit: log diff --git a/modules/integrao/select_feature/resources/usr/bin/integrao_select_feature.py b/modules/integrao/select_feature/resources/usr/bin/integrao_select_feature.py index 73b5d55..50cc40a 100755 --- a/modules/integrao/select_feature/resources/usr/bin/integrao_select_feature.py +++ b/modules/integrao/select_feature/resources/usr/bin/integrao_select_feature.py @@ -18,7 +18,7 @@ --embedding_dims=E_DIMS Latent space dimensionality [default: 64]. --fusing_iteration=FUSE_IT SNF diffusion iterations [default: 30]. --normalization_factor=NF Normalization factor [default: 1.0]. - --alignment_epochs=ALI_EPO Unsupervised alignment epochs [default: 300]. + --alignment_epochs=ALI_EPO Unsupervised alignment epochs [default: 500]. --finetune_epochs=FINT_EPO Supervised fine-tuning epochs [default: 800]. --beta=BETA Beta loss weight [default: 1.0]. --mu=MU Mu loss weight [default: 0.5]. @@ -33,6 +33,102 @@ import numpy as np import torch import os +import json +from pathlib import Path + + + + + + + +def make_json_serializable(value): + if isinstance(value, np.generic): + return value.item() + + if isinstance(value, np.ndarray): + return value.tolist() + + if isinstance(value, dict): + return { + str(key): make_json_serializable(item) + for key, item in value.items() + } + + if isinstance(value, (list, tuple)): + return [ + make_json_serializable(item) + for item in value + ] + + return value + + +def make_parameter_record(value, treatment): + allowed_treatments = { + "tuned", + "fixed", + "default", + "data-derived" + } + + if treatment not in allowed_treatments: + raise ValueError( + f"Unknown parameter treatment: {treatment}" + ) + + return { + "value": make_json_serializable(value), + "treatment": treatment + } + + +def write_integrao_hyperparameters( + dataset_name, + method, + hparams, + n_times +): + seeds = list(range(1, n_times + 1)) + + parameters = { + parameter: make_parameter_record( + value=value, + treatment="fixed" + ) + for parameter, value in hparams.items() + } + + selection = { + "strategy": "none", + "feature_scoring": "integrated_gradients", + "training_repetitions": n_times, + "seeds": seeds, + "aggregation": "mean_across_seeds" + } + + result = { + "schema_version": "1.0", + "run_id": f"{dataset_name}__{method}", + "dataset": dataset_name, + "method": method, + "analysis_stage": "model_selection", + "selection": selection, + "parameters": parameters + } + + output_path = f"{method}-{dataset_name}_selected_hyperparameters.json" + + + with open(output_path, "w", encoding="utf-8") as handle: + json.dump( + result, + handle, + indent=2, + ensure_ascii=False + ) + + return output_path # ------------------------- @@ -49,6 +145,14 @@ def get_default_hyperparameters(): num_classes=2, ) return DEFAULT_HPARAMS + +def resolve_hyperparameters(hparams=None): + return { + **get_default_hyperparameters(), + **(hparams or {}) + } + + # ------------------------- @@ -244,16 +348,30 @@ def main( # Common objects to be used through the place full_dfs, modality_names, truelabel = load_data(mdata_path) + + resolved_hparams = resolve_hyperparameters(hparams) # Run the feature selection loop all_feat_dfs = run_feat_selection_loop( outdir=outdir, full_dfs=full_dfs, modality_names=modality_names, - truelabel=truelabel, hparams=hparams, dataset_name=dataset_name, n_times=n_times) + truelabel=truelabel, hparams=resolved_hparams, dataset_name=dataset_name, n_times=n_times) # --- Aggregate feature importances --- agg_feat_df = summarize_feature_importance(all_feat_dfs) # --- Now add metadata info --- feats_df = wrangle_result_table(agg_feat_df, method, dataset_name) # Lastly write out to file filename = f"{method}-{dataset_name}_features_selected.csv" + # ANd write the hyperparams to file + hyperparameter_path = write_integrao_hyperparameters( + dataset_name=dataset_name, + method=method, + hparams=resolved_hparams, + n_times=n_times + ) + print( + "Selected hyperparameters written to: " + f"{hyperparameter_path}" + ) + feats_df.to_csv(filename, index=False, header=True) return feats_df diff --git a/modules/merge_result_table/main.nf b/modules/merge_result_table/main.nf index 6794d42..9ed9025 100644 --- a/modules/merge_result_table/main.nf +++ b/modules/merge_result_table/main.nf @@ -1,43 +1,48 @@ process MERGE_RESULT_TABLE { - tag "${method_name}" - debug true - label 'process_single' - label 'codia' + tag "${method_name}" + debug true + label 'process_single' + label 'codia' - publishDir ( - //path: "${params.outdir}/${task.process.tokenize(':')[-1].toLowerCase()}/${method_name}", - path: "${params.outdir}/${task.process.tokenize(':').join('/').toLowerCase()}", - mode: 'copy', - overwrite: true - ) + publishDir ( + //path: "${params.outdir}/${task.process.tokenize(':')[-1].toLowerCase()}/${method_name}", + path: "${params.outdir}/${task.process.tokenize(':').join('/').toLowerCase()}", + mode: 'copy', + overwrite: true + ) // Input output blocks - input: - tuple val(method_name), path(list_result_tables) + input: + tuple val(method_name), path(list_result_tables) val(saveMode) + val(outcome_type) output: tuple val(method_name), path('*.csv'), optional: true, emit: csv_results - tuple val(method_name), path('*.txt'), optional: true, emit: txt_results path('*log*'), optional: true, emit: log + script: - // This script assumes to read in list of files - // The "${variable}" or ${variable.join(' ')} does the trick to make it one whole string - // Different branch to call different script under resources/usr/bin to merge and collect results - if (saveMode == "method") + def file_list = list_result_tables.collect { it.toString().replace("'", "'\"'\"'") }.join('\n') + + // Different branch to call different script under resources/usr/bin to merge and collect results + if (saveMode == "method") """ + printf '%s\\n' '${file_list}' > input_files.txt combine_tables.R \ - --tables="${list_result_tables}" \ + --input_list=input_files.txt \ + --outcome_type=${outcome_type} \ --method_name=${method_name} \ --methodMode > \ ${method_name}-${task.process.tokenize(':')[-1].toLowerCase()}.log """ - else if (saveMode == "language") + else if (saveMode == "language") """ + printf '%s\\n' '${file_list}' > input_files.txt combine_tables.R \ - --tables="${list_result_tables}" \ + --input_list=input_files.txt \ + --outcome_type=${outcome_type} \ --method_name=${method_name} > \ ${method_name}-${task.process.tokenize(':')[-1].toLowerCase()}.log """ - else + else error "Did not provide save mode to be language or method" - } \ No newline at end of file + } diff --git a/modules/merge_result_table/resources/usr/bin/combine_tables.R b/modules/merge_result_table/resources/usr/bin/combine_tables.R index 8147609..e5893bf 100755 --- a/modules/merge_result_table/resources/usr/bin/combine_tables.R +++ b/modules/merge_result_table/resources/usr/bin/combine_tables.R @@ -6,50 +6,170 @@ Author: Tony Liang Usage: combine_tables.R [options] - + Options: - --tables=TABLES Joined path of list of the result table of particular fold or dataset of method [default: empty] + --input_list=FILE Text file containing one input path per line --method_name=MNAME Name of method run on [default: empty] --methodMode Collecting results for method specific [default: false] + --outcome_type=OUTCOME Outcome type to merge results from. One of 'classification' or 'survival'. [default: classification] " -# Parase docopt +# Parse docopt opt <- docopt::docopt(doc) -# Helper to check format of table and transform it -convert_table_format <- function(table) { - # Make sure first col is sample name - cols <- colnames(table) - # First col is always sample_name - match_first_col <- cols[1] == "sample_name" - if (!match_first_col) { - stop("First column is not sample_name, wrong naming or missed somewhere") + + +# ====================================================================== + +# Check required columns and return them in a consistent order +select_relevant_cols <- function(table, relevant_cols) { + missing_cols <- setdiff(relevant_cols, colnames(table)) + + if (length(missing_cols) > 0L) { + stop( + "Result table is missing required columns: ", + paste(missing_cols, collapse = ", "), + call. = FALSE + ) } - # TODO: Should at least contain sample name, phat, method_name, dataset - # TODO: allow this extra column of fold - relevant_cols <- c("sample_name", "y", "phat", "method_name", "dataset", "fold") - contains_relevant_cols <- cols %in% relevant_cols - if(!all(contains_relevant_cols)) { - warning("\nResult table might not contain all relevant columns:\n ", - paste(relevant_cols, collapse=", "), "\n", "current is: ", - paste(cols, collapse=", "), "\n") + + return(table[, relevant_cols, drop = FALSE]) +} + + + + +#clean_classification_table <- function(table, cols) { +# relevant_cols <- c("sample_name", "y", "phat", "method_name", "dataset", "fold") +# contains_relevant_cols <- cols %in% relevant_cols +# if(!all(contains_relevant_cols)) { +# warning("\nResult table might not contain all relevant columns:\n ", +# paste(relevant_cols, collapse=", "), "\n", "current is: ", +# paste(cols, collapse=", "), "\n") # TODO: Find a better way for this? - table <- table[, relevant_cols] - } +# table <- table[, relevant_cols] +# } # Also check if has right column types - bv <- c(1, 0) - is_binary_y <- all(is.element(table$y, bv)) - if (!is_binary_y) { - message("\nConverting y to binary output\n") - table$y <- ifelse(table$y == "yes", 1, 0) +# bv <- c(1, 0) +# is_binary_y <- all(is.element(table$y, bv)) +# if (!is_binary_y) { +# message("\nConverting y to binary output\n") +# table$y <- ifelse(table$y == "yes", 1, 0) +# } +# return(table) +#} + + +#clean_survival_table <- function(table, cols) { +# relevant_cols <- c("sample_name", "lp", "surv_365", "surv_548", "surv_730", "surv_1095", "surv_1461", "surv_1826", "time", "status", "method_name", "dataset", "fold") +# contains_relevant_cols <- cols %in% relevant_cols +# if (!all(contains_relevant_cols)) { +# warning("\nResult table might not contain all relevant columns:\n ", +# paste(relevant_cols, collapse=", "), "\n", "current is: ", +# paste(cols, collapse=", "), "\n") +# # And reorder it +# table <- table[, relevant_cols] +# } +# return(table) +#} + + +# Helper to check format of table and transform it +#convert_table_format <- function(table, outcome_type) { +# # Make sure first col is sample name +# cols <- colnames(table) +# # First col is always sample_name +# match_first_col <- cols[1] == "sample_name" +# if (!match_first_col) { +# stop("First column is not sample_name, wrong naming or missed somewhere") +# } + # TODO: Should at least contain sample name, phat, method_name, dataset + # TODO: allow this extra column of fold +# if (outcome_type == "classification") { +# table <- clean_classification_table(table, cols) +# +# } +# if (outcome_type == "survival") { +# table <- clean_survival_table(table, cols) +# } + # Otherwise return nothing and prompt to error if outcome type is incorrect +# return(NULL) +#} + + + +clean_classification_table <- function(table) { + relevant_cols <- c( + "sample_name", "y", "phat", + "method_name", "dataset", "fold" + ) + + table <- select_relevant_cols(table, relevant_cols) + + # Normalize labels before validating and converting + y <- tolower(trimws(as.character(table$y))) + valid_labels <- c("0", "1", "no", "yes") + + invalid <- is.na(y) | !(y %in% valid_labels) + + if (any(invalid)) { + stop( + "Invalid or missing values in y: ", + paste(unique(y[invalid]), collapse = ", "), + ". Expected 0/1 or yes/no.", + call. = FALSE + ) } + + table$y <- as.integer(y %in% c("1", "yes")) + + return(table) +} + + +clean_survival_table <- function(table) { + relevant_cols <- c( + "sample_name", "lp", + "surv_365", "surv_548", "surv_730", + "surv_1095", "surv_1461", "surv_1826", + "time", "status", + "method_name", "dataset", "fold" + ) + + table <- select_relevant_cols(table, relevant_cols) + return(table) } -main <- function(tables, method_name, methodMode, readMode="csv", pattern="-result.*") { + + + + +convert_table_format <- function(table, outcome_type) { + # sample_name can be anywhere in the input; + # the cleaning functions move it to the first column. + if (outcome_type == "classification") { + return(clean_classification_table(table)) + } else if (outcome_type == "survival") { + return(clean_survival_table(table)) + } else { + stop( + "Unsupported outcome_type: ", outcome_type, + ". Expected 'classification' or 'survival'.", + call. = FALSE + ) + } +} + + + + +main <- function(input_list, method_name, methodMode, outcome_type, readMode="csv", pattern="-result.*") { # Special script to handle here - tables <- strsplit(tables, " ") |> unlist() + tables <- readLines(input_list, warn=FALSE) + tables <- tables[nzchar(trimws(tables))] +# tables <- strsplit(tables, " ") |> unlist() # Store to list and bind by rows laters to_bind <- list() # Check which string to replace instead @@ -62,9 +182,10 @@ main <- function(tables, method_name, methodMode, readMode="csv", pattern="-resu # readMode <- "table" #} #} - for (table_path in tables) { + for (i in seq_along(tables)) { # TODO: need to make this label and identifier better # Get everything before last hypen - to retrieve unique label + table_path <- tables[[i]] label <- gsub(pattern, "", table_path) message("\nThis is label: ", label, "\n") #table <- switch( @@ -75,7 +196,7 @@ main <- function(tables, method_name, methodMode, readMode="csv", pattern="-resu # Force sample name to be character, as it could come in numbers as well as id names table <- read.csv(table_path, header=TRUE, colClasses=c("sample_name"="character")) # Check the format of each table aligns before adding into list - to_bind[[label]] <- convert_table_format(table) # Add it to list + to_bind[[i]] <- convert_table_format(table, outcome_type) # Add it to list } message("\nMerging", length(to_bind), "tables\n") # Flatten these tables by merging rows @@ -90,4 +211,4 @@ main <- function(tables, method_name, methodMode, readMode="csv", pattern="-resu return(merged_table) } -main(tables=opt$tables, method_name=opt$method_name, methodMode=opt$methodMode) +main(input_list=opt$input_list, method_name=opt$method_name, methodMode=opt$methodMode, outcome_type=opt$outcome_type) diff --git a/modules/merge_selected_hyperparameters/main.nf b/modules/merge_selected_hyperparameters/main.nf new file mode 100644 index 0000000..290171c --- /dev/null +++ b/modules/merge_selected_hyperparameters/main.nf @@ -0,0 +1,25 @@ +process MERGE_SELECTED_HYPERPARAMETERS { + label "generic" + tag "Merge selected hyperparameters" + + publishDir ( + path: "${params.outdir}/${task.process.tokenize(':')[-1].toLowerCase()}", + mode: 'copy', + overwrite: true + ) + + input: + path(input_files) + + output: + path("selected_hyperparameters.json"), emit: merged_json + + path("selected_hyperparameters.csv"), emit: merged_csv + + script: + """ + merge_selected_hyperparameters.py \ + --inputs ${input_files} > \ + ${task.process.tokenize(':')[-1].toLowerCase()}.log + """ +} \ No newline at end of file diff --git a/modules/merge_selected_hyperparameters/resources/usr/bin/merge_selected_hyperparameters.py b/modules/merge_selected_hyperparameters/resources/usr/bin/merge_selected_hyperparameters.py new file mode 100755 index 0000000..96010b5 --- /dev/null +++ b/modules/merge_selected_hyperparameters/resources/usr/bin/merge_selected_hyperparameters.py @@ -0,0 +1,293 @@ +#!/usr/bin/env python + +""" +Merge model-selection hyperparameter JSON files and generate a long-format CSV. +""" + +import argparse +import csv +import json +from pathlib import Path + + +REQUIRED_FIELDS = { + "dataset", + "method", + "analysis_stage", + "parameters", +} + + +def parse_args(): + parser = argparse.ArgumentParser( + description=( + "Merge selected-hyperparameter JSON files and convert the " + "parameters into a long-format CSV." + ) + ) + + parser.add_argument( + "--inputs", + nargs="+", + required=True, + help="Input selected-hyperparameter JSON files." + ) + + parser.add_argument( + "--json-output", + default="selected_hyperparameters.json", + help="Path for the merged JSON output." + ) + + parser.add_argument( + "--csv-output", + default="selected_hyperparameters.csv", + help="Path for the flattened CSV output." + ) + + return parser.parse_args() + + +def read_hyperparameter_json(path): + path = Path(path) + + with path.open("r", encoding="utf-8") as handle: + result = json.load(handle) + + missing_fields = REQUIRED_FIELDS.difference(result) + + if missing_fields: + missing = ", ".join(sorted(missing_fields)) + raise ValueError( + f"{path}: missing required fields: {missing}" + ) + + if not isinstance(result["parameters"], dict): + raise TypeError( + f"{path}: 'parameters' must be a JSON object" + ) + + return result + + +def value_type(value): + """ + Return a simple cross-language description of a JSON value. + """ + + if value is None: + return "null" + + if isinstance(value, bool): + return "boolean" + + if isinstance(value, int): + return "integer" + + if isinstance(value, float): + return "numeric" + + if isinstance(value, str): + return "string" + + if isinstance(value, list): + return "array" + + if isinstance(value, dict): + return "object" + + return type(value).__name__ + + +def csv_value(value): + """ + Convert a scalar JSON value into a CSV-safe representation. + """ + + if value is None: + return "" + + if isinstance(value, bool): + return str(value).lower() + + if isinstance(value, (dict, list)): + return json.dumps(value, ensure_ascii=False) + + return value + + +def flatten_value(parameter_path, value): + """ + Recursively flatten a nested parameter value. + + Examples + -------- + keepX = { + "genomics": [50, 25], + "proteomics": [20, 20] + } + + becomes: + + keepX.genomics[1] = 50 + keepX.genomics[2] = 25 + keepX.proteomics[1] = 20 + keepX.proteomics[2] = 20 + """ + + if isinstance(value, dict): + if not value: + yield parameter_path, value + return + + for key in sorted(value): + child_path = f"{parameter_path}.{key}" + yield from flatten_value(child_path, value[key]) + + return + + if isinstance(value, list): + if not value: + yield parameter_path, value + return + + # Use one-based indexing because these usually represent + # latent components or ordered model layers. + for index, item in enumerate(value, start=1): + child_path = f"{parameter_path}[{index}]" + yield from flatten_value(child_path, item) + + return + + yield parameter_path, value + + +def make_parameter_rows(run): + selection = run.get("selection", {}) + rows = [] + + for parameter, record in run["parameters"].items(): + if not isinstance(record, dict): + raise TypeError( + f"{run['dataset']}/{run['method']}/{parameter}: " + "parameter record must be a JSON object" + ) + + if "value" not in record: + raise ValueError( + f"{run['dataset']}/{run['method']}/{parameter}: " + "parameter record is missing 'value'" + ) + + treatment = record.get("treatment", "unknown") + + for parameter_path, selected_value in flatten_value( + parameter, + record["value"] + ): + rows.append({ + "dataset": run["dataset"], + "method": run["method"], + "analysis_stage": run["analysis_stage"], + "parameter_path": parameter_path, + "selected_value": csv_value(selected_value), + "treatment": treatment, + # "selection_metric": selection.get("metric", ""), + # "cv_type": selection.get("cv_type", ""), + # "n_iter": selection.get("n_iter", ""), + # "random_state": selection.get("random_state", ""), + # "best_cv_score": selection.get("best_cv_score", ""), + }) + + return rows + + +def write_merged_json(runs, output_path): + merged_result = { + "schema_version": "1.0", + "analysis_stage": "model_selection", + "n_runs": len(runs), + "runs": runs, + } + + with Path(output_path).open("w", encoding="utf-8") as handle: + json.dump( + merged_result, + handle, + indent=2, + ensure_ascii=False + ) + + +def write_long_csv(rows, output_path): + columns = [ + "dataset", + "method", + "analysis_stage", + "parameter_path", + "selected_value", + "treatment", + # "selection_metric", + # "cv_type", + # "n_iter", + # "random_state", + # "best_cv_score", + ] + + with Path(output_path).open( + "w", + encoding="utf-8", + newline="" + ) as handle: + writer = csv.DictWriter(handle, fieldnames=columns) + writer.writeheader() + writer.writerows(rows) + + +def main(): + args = parse_args() + + runs = [ + read_hyperparameter_json(path) + for path in args.inputs + ] + + # Deterministic ordering makes outputs easier to compare with git diff. + runs.sort( + key=lambda run: ( + str(run.get("dataset", "")), + str(run.get("method", "")), + ) + ) + + rows = [] + + for run in runs: + rows.extend(make_parameter_rows(run)) + + rows.sort( + key=lambda row: ( + row["dataset"], + row["method"], + row["parameter_path"], + ) + ) + + write_merged_json( + runs=runs, + output_path=args.json_output + ) + + write_long_csv( + rows=rows, + output_path=args.csv_output + ) + + print( + f"Merged {len(runs)} model-selection runs " + f"into {len(rows)} parameter rows." + ) + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/modules/mofa/predict/main.nf b/modules/mofa/predict/main.nf index b1d1b26..a32b119 100644 --- a/modules/mofa/predict/main.nf +++ b/modules/mofa/predict/main.nf @@ -21,6 +21,7 @@ process MOFA_PREDICT { input: tuple val(dataset_name), val(fold_name), path(model) tuple val(dataset_name), val(fold_name), path(test_path) + val(outcome_type) output: tuple val(dataset_name), val(fold_name), val("mofa"), path("*result*"), emit: result_table //tuple val(dataset_name), val(fold_name), val("mofa"), path("*weight*"), emit: weight @@ -31,7 +32,8 @@ process MOFA_PREDICT { predict_mofa.R \ --model=${model} \ --test_path=${test_path} \ - --label=${data_label} > \ + --label=${data_label} \ + --outcome_type=${outcome_type} > \ ${data_label}-${getPublishPath(task.process).tokenize('/')[-1]}.log echo mofa > method_name diff --git a/modules/mofa/predict/resources/usr/bin/predict_mofa.R b/modules/mofa/predict/resources/usr/bin/predict_mofa.R index 135839c..26600d2 100755 --- a/modules/mofa/predict/resources/usr/bin/predict_mofa.R +++ b/modules/mofa/predict/resources/usr/bin/predict_mofa.R @@ -14,6 +14,7 @@ Options: --test_path=TEST_PATH Path containing test data [default: null] --label=LABEL Label of id and fold of data [default: data-fold_i] --output_ext=EXT Extension of output table to save [default: csv] + --outcome_type=TYPE Outcome type to PREDICT from. One of 'classification' or 'survival'. [default: classification] " library(here) @@ -34,9 +35,10 @@ source(here(pipeline_dir, "bin/wrangling/get_result_table.R")) # Parse docopt opt <- docopt::docopt(doc) -main <- function(model_path, test_path, label, output_ext, method_name="mofa", - digit = 3, s = 0) { - # Load model (from same fold train portion) + +# Fun for classification +predict_classification <- function(model_path, test_path, label, method_name, digit, s=0) { + # Load model (from same fold train portion) model <- readRDS(model_path) # Load test data (from same fold test portion) test_data <- readRDS(test_path) # NOTE here data is n x p @@ -58,12 +60,85 @@ main <- function(model_path, test_path, label, output_ext, method_name="mofa", method_name=method_name, test_data=test_data, digit=digit ) + return(result_table) +} + +# Survival: risk score and S(t) at the horizons from the coxnet model. +# To be comparable with the sksurv methods (sksurv_predict.py), the output: +# - lp is the linear predictor, higher = riskier +# - S(t) at the same horizons (days), NA after the last train event +# (both handled by baseline_surv from run_mofa.R) +# - has the same columns: sample_name, lp, surv_, time, status, +# method_name, dataset, fold +predict_survival <- function(model_path, test_path, label, method_name, digit, s="lambda.min") { + # Load model (from same fold train portion) + model <- readRDS(model_path) + # Load test data (from same fold test portion) + test_data <- readRDS(test_path) # NOTE here data is n x p + test_x <- as.matrix(test_data$X$embeddings) + # Risk score, lambda chosen by inner CV in run_mofa.R + beta <- coef(model$cvfit, s = s) + + message("test_x: ", paste(dim(test_x), collapse = " x ")) + message("beta: ", paste(dim(beta), collapse = " x ")) + + stopifnot(ncol(test_x) == nrow(beta)) + + message("Predicting survival on test data with lambda = ", s) + lp <- stats::predict(model$cvfit, newx = test_x, s = s, + type = "link") |> as.numeric() + + # Cox model: S(t | x) = S0(t) ^ exp(lp), S0 estimated on the train fold + surv <- outer(exp(lp), model$baseline_surv, function(r, s0) s0 ^ r) + + # Result table, dataset and fold parsed from the label (dataset-..fold_i..) + result_table <- data.frame(sample_name = rownames(test_x), lp = round(lp, digit)) + for (h in names(model$baseline_surv)) { + result_table[[paste0("surv_", h)]] <- round(surv[, h], digit) + } + # And wrangle here + result_table$time <- test_data$Y$time + result_table$status <- test_data$Y$status + result_table$method_name <- method_name + result_table$dataset <- sub("-[^-]*fold_[0-9]+.*$", "", label) + result_table$fold <- regmatches(label, regexpr("fold_[0-9]+", label)) + # Return it out + return(result_table) +} + + + +main <- function(model_path, test_path, label, output_ext, outcome_type="classification", method_name="mofa", + digit = 3, s = 0) { + + # Init result table to be empty + result_table <- data.frame() + + if (outcome_type == "classification") { + result_table <- predict_classification( + model_path=model_path, test_path=test_path, + label=label, method_name=method_name, digit=digit, s=s + ) + } else if (outcome_type == "survival") { + s <- "lambda.min" # Coerce it to use lambda min, survival always runs on cv glmnet + result_table <- predict_survival( + model_path=model_path, test_path=test_path, + label=label, method_name=method_name, digit=digit, s=s + ) + } else { + stop("Unsupported outcome_type: ", outcome_type, + ". Expected 'classification' or 'survival'.", + call. = FALSE) + } + + + # Write to files cat("\nSaving as", output_ext, "format\n") result_file <- paste(label, paste0("result_table", ".", output_ext), sep="-") # Save to disk write.csv(result_table, result_file, row.names = FALSE) - return(pred_probs) + return(result_table) } @@ -71,6 +146,7 @@ main <- function(model_path, test_path, label, output_ext, method_name="mofa", main(model_path=opt$model_path, test_path=opt$test_path, label=opt$label, + outcome_type=opt$outcome_type, output_ext=opt$output_ext ) diff --git a/modules/mofa/preprocess/main.nf b/modules/mofa/preprocess/main.nf index ff7efbe..81edae4 100644 --- a/modules/mofa/preprocess/main.nf +++ b/modules/mofa/preprocess/main.nf @@ -47,6 +47,7 @@ process MOFA_PREPROCESS { /* data name identifier, MAE portion of data , and dir that contains txt in it" */ tuple val(dataset_name), path(mae_path), path(split_dir) val(num_factors) + val(outcome_type) // TODO: This part could be different dependent on method output: /* @@ -65,7 +66,8 @@ process MOFA_PREPROCESS { --mae_path=${mae_path} \ --split_dir=${split_dir} \ --dataset_name=${dataset_name} \ - --num_factors=${num_factors} > \ + --num_factors=${num_factors} \ + --outcome_type=${outcome_type} > \ ${dataset_name}-${getPublishPath(task.process).tokenize('/')[-1].toLowerCase()}.log echo ${dataset_name} > dataset_name diff --git a/modules/mofa/preprocess/resources/usr/bin/preprocess_mofa.R b/modules/mofa/preprocess/resources/usr/bin/preprocess_mofa.R index 8b44516..f1f8f1a 100755 --- a/modules/mofa/preprocess/resources/usr/bin/preprocess_mofa.R +++ b/modules/mofa/preprocess/resources/usr/bin/preprocess_mofa.R @@ -2,7 +2,7 @@ # Script to prepare mofa input doc <- "This script is to run MOFA method from MOFA2 package, train only -it could possibly be ran on a inner CV model, output is a modelel for prediction +it could possibly be ran on a inner CV model, output is a model for prediction usage in downstream. Usage: @@ -13,6 +13,7 @@ Options: --split_dir=SPLIT_DIR Directory containing list of txt file [default: empty] --dataset_name=NAME Name of dataset that is splitting [default: empty] --num_factors=NUM_FACTOR Number of factors to supply into MOFA [default: 1] + --outcome_type=OUTCOME Outcome type to merge results from. One of 'classification' or 'survival'. [default: classification] " # Parse cli args @@ -20,9 +21,9 @@ opt <- docopt::docopt(doc) # Load libraries library(MOFA2) -library(MultiAssayExperiment) +suppressPackageStartupMessages(library(MultiAssayExperiment)) library(here) -library(dplyr) +suppressPackageStartupMessages(library(dplyr)) # Python related default_python <- "/usr/bin/python" @@ -59,18 +60,57 @@ load_test_splits <- function(split_dir, pattern=".txt", ...) { } -# Use this function to reconstruct mae -reconstruct_mae <- function(mae) { - # Given an mae with delayed matrices, we could load it into - # memory and make it of HDF5 arrays instead - X <- mae@ExperimentList |> lapply(as.matrix) - y <- mae$response - # Construct MAE - new_mae <- MultiAssayExperiment::MultiAssayExperiment(experiments = X) - new_mae$response <- y - return(new_mae) -} +reconstruct_mae <- function( + mae, + outcome_type = "classification", + response_col = "response", + time_col = "time", + status_col = "status" +) { + if (!outcome_type %in% c("classification", "survival")) { + stop("outcome_type must be 'classification' or 'survival'") + } + + required_cols <- if (outcome_type == "classification") { + response_col + } else { + c(time_col, status_col) + } + + cd <- as.data.frame(mae@colData) + missing_cols <- setdiff(required_cols, colnames(cd)) + if (length(missing_cols) > 0L) { + stop("Missing outcome columns: ", paste(missing_cols, collapse = ", ")) + } + + # Materialize each experiment as an in-memory matrix. + X <- lapply( + MultiAssayExperiment::experiments(mae), + function(x) { + # For SummarizedExperiment, use its first assay. + if (methods::is(x, "SummarizedExperiment")) { + x <- SummarizedExperiment::assay(x) + } + as.matrix(x) + } + ) + + if (outcome_type == "classification") { + cd$response <- cd[[response_col]] + } else { + # Extract both before assignment in case source names overlap. + time <- cd[[time_col]] + status <- cd[[status_col]] + cd$time <- time + cd$status <- status + } + + MultiAssayExperiment::MultiAssayExperiment( + experiments = X, + colData = cd + ) +} get_seed <- function(dataset_name) { d_int <- utf8ToInt(dataset_name) # Convert dataset name to integer @@ -82,7 +122,7 @@ get_seed <- function(dataset_name) { # ================================================================================= # MAIN entrance point -main <- function(mae_path, split_dir, dataset_name, num_factors) { +main <- function(mae_path, split_dir, dataset_name, num_factors, outcome_type="classification") { # Seed for reproducibility seed <- get_seed(dataset_name) # Convert dataset name to integer and sum it to get a seed set.seed(seed) @@ -138,6 +178,13 @@ main <- function(mae_path, split_dir, dataset_name, num_factors) { mofa_emb <- load_model(file=mofa_emb_file, remove_inactive_factors = FALSE) # Although to make predictions, need its embeddings and use a glmnet on prediction factors <- get_factors(mofa_emb, factors="all") + + # ================================================================================= + # NOTE: The MOFA factors are then extracted and construct as MAEs as input objects + # for downstream model fitting and prediction + # For classification it uses glmnent family binomial + # For survival it uses glmnet family cox + # ================================================================================= # Then use this new mae mae <- MultiAssayExperiment(experiments = list(embeddings=factors$group1 |> t() ), colData = raw_col_data @@ -149,8 +196,8 @@ main <- function(mae_path, split_dir, dataset_name, num_factors) { # First subset both split <- test_splits[[fold_name]] # TODO: Transpose data only when method requires it to - tr_mae <- mae[, -split, drop=TRUE] |> reconstruct_mae() - te_mae <- mae[, split, drop=TRUE] |> reconstruct_mae() + tr_mae <- mae[, -split, drop=TRUE] |> reconstruct_mae(outcome_type=outcome_type) + te_mae <- mae[, split, drop=TRUE] |> reconstruct_mae(outcome_type=outcome_type) # Then save each fold's train and test portion as subdirectory of fold name cat("\nSaving for", fold_name, "\n") if (!dir.exists(fold_name)) { @@ -173,5 +220,6 @@ main( mae_path=opt$mae_path, split_dir=opt$split_dir, dataset_name=opt$dataset_name, - num_factors=as.numeric(opt$num_factors) + num_factors=as.numeric(opt$num_factors), + outcome_type=opt$outcome_type ) diff --git a/modules/mofa/select_feature/main.nf b/modules/mofa/select_feature/main.nf index 699bab0..a2de993 100644 --- a/modules/mofa/select_feature/main.nf +++ b/modules/mofa/select_feature/main.nf @@ -21,6 +21,7 @@ process MOFA_SELECT_FEATURE { val(num_factors) output: path("*.csv"), emit: features + path("*.json"), emit: hyperparameters path("*plot*"), optional:true, emit: plot path("*.log"), emit: log diff --git a/modules/mofa/select_feature/resources/usr/bin/mofa_select_features.R b/modules/mofa/select_feature/resources/usr/bin/mofa_select_features.R index f3c5f2e..1054d11 100755 --- a/modules/mofa/select_feature/resources/usr/bin/mofa_select_features.R +++ b/modules/mofa/select_feature/resources/usr/bin/mofa_select_features.R @@ -26,27 +26,39 @@ default_python <- "/usr/bin/python" reticulate::use_python(default_python) # Fun to prepare mofa object for training -mofa_pre <- function(mae, num_factors, scale_views=FALSE) { - # TODO: Could add more options here - # See https://du-bii.github.io/module-6-Integrative-Bioinformatics/2019/Session2-3/practical_MOFA.html +mofa_pre <- function( + mae, + num_factors, + seed, + scale_views = FALSE +) { mofa_obj <- create_mofa(mae) + + # Data options data_opts <- get_default_data_options(mofa_obj) - # Options of data data_opts$scale_views <- scale_views - + + # Model options model_opts <- get_default_model_options(mofa_obj) - # Also add number of factors from nextflow params model_opts$num_factors <- num_factors + # Training options train_opts <- get_default_training_options(mofa_obj) - # Notice this overrides the previous object created - mofa_obj <- prepare_mofa( + train_opts$seed <- seed + + prepared_obj <- prepare_mofa( object = mofa_obj, data_options = data_opts, model_options = model_opts, training_options = train_opts ) - return(mofa_obj) + + return(list( + object = prepared_obj, + data_options = data_opts, + model_options = model_opts, + training_options = train_opts + )) } @@ -57,6 +69,99 @@ get_seed <- function(dataset_name) { return(seed) } +tag_mofa_options <- function(options, prefix, treatment_overrides = NULL) { + tagged <- lapply(names(options), function(option_name) { + treatment <- "default" + + if ( + !is.null(treatment_overrides) && + option_name %in% names(treatment_overrides) + ) { + treatment <- unname(treatment_overrides[[option_name]]) + } + + list( + value = options[[option_name]], + treatment = treatment + ) + }) + + names(tagged) <- paste(prefix, names(options), sep = ".") + return(tagged) +} + + +write_selected_hyperparameters <- function( + dataset_name, + data_options, + model_options, + training_options, + method = "mofa" +) { + data_parameters <- tag_mofa_options( + options = data_options, + prefix = "data", + treatment_overrides = c( + scale_views = "fixed", + views = "data-derived", + groups = "data-derived" + ) + ) + + model_parameters <- tag_mofa_options( + options = model_options, + prefix = "model", + treatment_overrides = c( + num_factors = "fixed" + ) + ) + + training_parameters <- tag_mofa_options( + options = training_options, + prefix = "training", + treatment_overrides = c( + seed = "data-derived" + ) + ) + + result <- list( + schema_version = "1.0", + run_id = paste(method, dataset_name, sep = "-"), + dataset = dataset_name, + method = method, + analysis_stage = "model_selection", + selection = list( + strategy = "fixed_configuration", + hyperparameter_tuning = FALSE + ), + parameters = c( + data_parameters, + model_parameters, + training_parameters + ) + ) + + output_file <- paste0( + method, + "-", + dataset_name, + "_selected_hyperparameters.json" + ) + + jsonlite::write_json( + x = result, + path = output_file, + pretty = TRUE, + auto_unbox = TRUE, + null = "null", + na = "null" + ) + + message("\nSelected hyperparameters written to: ", output_file) + + return(output_file) +} + # Main entrypoint here @@ -68,7 +173,14 @@ main <- function(mae_path, dataset_name, num_factors, raw_mae <- loadHDF5MultiAssayExperiment(mae_path) # Then starting to use mofa here # TODO: add options into scale view or scale anything - mofa_obj <- mofa_pre(mae=raw_mae, num_factors=num_factors, scale_views=scale_views) + mofa_setup <- mofa_pre( + mae = raw_mae, + num_factors = num_factors, + seed = seed, + scale_views = scale_views + ) + + mofa_obj <- mofa_setup$object # Then train on the mofa obj mofa_emb_file <- paste( paste(method, dataset_name, sep="-"), "emb.hdf5", sep="_" @@ -105,6 +217,14 @@ main <- function(mae_path, dataset_name, num_factors, # write it to disk feats_file <- paste0(method, "-", dataset_name, "_", "features_selected", ".csv") write.csv(x=feats_df, file=feats_file, row.names=FALSE) + # And the hyperparams part + write_selected_hyperparameters( + dataset_name = dataset_name, + data_options = mofa_setup$data_options, + model_options = mofa_setup$model_options, + training_options = mofa_setup$training_options, + method = method + ) return(feats_df) } diff --git a/modules/mofa/train/main.nf b/modules/mofa/train/main.nf index 603ddfd..6987de1 100644 --- a/modules/mofa/train/main.nf +++ b/modules/mofa/train/main.nf @@ -25,6 +25,7 @@ process MOFA_TRAIN { input: tuple val(dataset_name), path(mae_path), path(fold_path) val(run_inner_cv) + val(outcome_type) output: tuple val(dataset_name), val(fold_path.name), path('*model.rds'), emit: model tuple val(dataset_name), val(fold_path.name), path('*test_data.rds'), emit: test_data @@ -37,6 +38,7 @@ process MOFA_TRAIN { --mae_path=${mae_path} \ --label=${data_label} \ --fold_path=${fold_path} \ + --outcome_type=${outcome_type} \ --run_inner_cv > \ ${data_label}-${getPublishPath(task.process).tokenize('/')[-1].toLowerCase()}.log echo ${dataset_name} > dataset_name @@ -46,7 +48,8 @@ process MOFA_TRAIN { run_mofa.R \ --mae_path=${mae_path} \ --label=${data_label} \ - --fold_path=${fold_path} > \ + --fold_path=${fold_path} \ + --outcome_type=${outcome_type} > \ ${data_label}-${getPublishPath(task.process).tokenize('/')[-1].toLowerCase()}.log echo ${dataset_name} > dataset_name """ diff --git a/modules/mofa/train/resources/usr/bin/run_mofa.R b/modules/mofa/train/resources/usr/bin/run_mofa.R index 99d2167..18db099 100755 --- a/modules/mofa/train/resources/usr/bin/run_mofa.R +++ b/modules/mofa/train/resources/usr/bin/run_mofa.R @@ -14,14 +14,17 @@ Options: --fold_path=FOLD_PATH Path to read current test fold --prefix=PREFIX Prefix to read HDF5 [default: pre] --run_inner_cv Run inner cv with train data or not [default: false] + --outcome_type=TYPE classification or survival [default: classification] + --time_col=TIME_COL colData column of survival time [default: time] + --status_col=STAT_COL colData column of event indicator [default: status] " # Load libraries -library(dplyr) -library(glmnet) +suppressPackageStartupMessages(library(dplyr)) +suppressPackageStartupMessages(library(glmnet)) library(here) -library(MOFA2) -library(MultiAssayExperiment) +suppressPackageStartupMessages(library(MOFA2)) +suppressPackageStartupMessages(library(MultiAssayExperiment)) library(stringr) # Gather the pipeline dir (THIS IS VERY UGGLY FIX) @@ -49,12 +52,47 @@ get_seed <- function(dataset_name) { } + +run_survival <- function(train_x, train_y, horizons, alpha=1) { + # glmnet cox rejects time <= 0 (e.g. TCGA patients with 0 days of follow-up): + # for fitting only, move them to half of the smallest positive time + train_time <- train_y$time + train_status <- train_y$status + if (any(train_time <= 0)) { + eps <- min(train_time[train_time > 0]) / 2 + message("\n", sum(train_time <= 0), " train samples with time <= 0 set to ", eps) + train_time <- pmax(train_time, eps) + } + cvfit <- cv.glmnet(train_x, survival::Surv(train_time, train_status), + family = "cox", alpha = alpha) + message("\ncoxnet: ", sum(coef(cvfit, s = "lambda.min") != 0), + " non-zero coefficients at lambda.min") + # Baseline survival S0(t) = S(t | lp = 0) at the horizons, estimated on the + # train fold, so predict_mofa.R only needs S(t | x) = S0(t) ^ exp(lp). + # NA after the last train event: the curve is flat there, as in sksurv_predict.py + s0_fit <- survival::survfit(cvfit, s = "lambda.min", x = train_x, + y = survival::Surv(train_time, train_y$status), + newx = matrix(0, 1, ncol(train_x))) + idx <- findInterval(horizons, s0_fit$time) + baseline_surv <- ifelse(idx == 0, 1, s0_fit$surv[pmax(idx, 1)]) + t_max_event <- max(train_time[train_y$status == 1]) + baseline_surv[horizons > t_max_event] <- NA + names(baseline_surv) <- horizons + + message("\nBaseline S0(t) at ", paste(horizons, collapse = ", "), ": ", + paste(round(baseline_surv, 3), collapse = ", ")) + + model <- list(outcome_type = "survival", cvfit = cvfit, baseline_surv = baseline_surv) + return(model) +} + + # Need to force use python default_python <- "/usr/bin/python" reticulate::use_python(default_python) # Main function to run -main <- function(mae_path, label, fold_path, run_inner_cv, prefix) { +main <- function(mae_path, label, fold_path, run_inner_cv, prefix, outcome_type="classification") { seed <- get_seed(label) # Get seed for reproducibility, this is important for mofa training and cv set.seed(seed) # Log the params used @@ -65,8 +103,13 @@ main <- function(mae_path, label, fold_path, run_inner_cv, prefix) { train_path <- d[str_detect(d, pattern = "_tr")] test_path <- d[str_detect(d, pattern = "_te")] # Then should read in the MAE - train_data <- load_MAE(train_path, prefix="train") |> extract_Xy() - test_data <- load_MAE(test_path, prefix="test") |> extract_Xy() + # When outcome is survival, then the response is a dataframe with time and status columns + # When outcome is classification, then the response is a factor with 2 levels + train_data <- load_MAE(train_path, prefix="train") |> extract_Xy(outcome_type=outcome_type) + message("\nTrain data has ", nrow(train_data$X$embeddings), " samples and ", ncol(train_data$X$embeddings), " features") + message("\nTrain data has ", length(train_data$Y), " samples in response variable") + message("\nHead of y is ", train_data$Y |> head()) + test_data <- load_MAE(test_path, prefix="test") |> extract_Xy(outcome_type=outcome_type) sample_names <- check_common_samples(train_data) cat("\nTotal of", length(sample_names), "samples:\n", sample_names) @@ -96,12 +139,20 @@ main <- function(mae_path, label, fold_path, run_inner_cv, prefix) { test_data$X$embeddings <- test_emb } - # Glmnet model (ACTUALLY using this to predict) - model <- glmnet( - x = train_x, - y = train_y, - family = "binomial" - ) + if (outcome_type == "classification") { + # Glmnet model (ACTUALLY using this to predict) + model <- glmnet( + x = train_x, + y = train_y, + family = "binomial" + ) + } + + if (outcome_type == "survival") { + horizons <- c(365,548,730,1095,1461,1826) + model <- run_survival(train_x=train_x, train_y=train_y, horizons=horizons) + } + } # ========================================= message("\nSaving files to", label, "\n") @@ -117,6 +168,7 @@ main(mae_path = opt$mae_path, label = opt$label, fold_path = opt$fold_path, prefix = opt$prefix, - run_inner_cv = opt$run_inner_cv + run_inner_cv = opt$run_inner_cv, + outcome_type = opt$outcome_type ) cat("Done") \ No newline at end of file diff --git a/modules/mofa/train/resources/usr/bin/run_mofa_survival.R b/modules/mofa/train/resources/usr/bin/run_mofa_survival.R new file mode 100644 index 0000000..e28cffe --- /dev/null +++ b/modules/mofa/train/resources/usr/bin/run_mofa_survival.R @@ -0,0 +1,129 @@ +#!/usr/bin/env Rscript +# Survival variant of MOFA training (task = "survival"). +# Reads the train fold of the MAE, fits the mofa survival model, saves +# model RDS + test_data RDS (same conventions as the classification scripts). +doc <- "Train mofa survival model. +Usage: + run_mofa_survival.R [options] +Options: + --mae_path=MAE_PATH Path to the train-fold MAE directory + --label=LABEL Dataset-fold label + --fold_path=FOLD_PATH Path to the fold directory (contains train/test prefixes) + --prefix=PREFIX HDF5 prefix [default: pre] +" +opt <- docopt::docopt(doc) +pipeline_dir <- gsub("/bin", "", (Sys.getenv("PATH") |> strsplit(":") |> unlist() |> tail(1))) +# Load scripts ======================================================== +source(here(pipeline_dir, "bin/rhelpers.R")) # This is included in nextflow bin path +# Loading generic utils +load_utils(here(pipeline_dir, "bin/logging")) +load_utils(here(pipeline_dir, "bin/preprocessing")) +load_utils(here(pipeline_dir, "bin/misc_utils")) + + + +get_seed <- function(dataset_name) { + d_int <- utf8ToInt(dataset_name) # Convert dataset name to integer + seed <- sum(d_int) + message("\nSeed used: ", seed) + return(seed) +} + +# Main function to run +main <- function(mae_path, label, fold_path, run_inner_cv, prefix) { + seed <- get_seed(label) # Get seed for reproducibility, this is important for mofa training and cv + set.seed(seed) + # Log the params used + args_used <- c(as.list(environment())) + logging_params(args_used) + cat("\nLooking at this fold:", fold_path, "\n") + d <- list.files(path=fold_path, full.names = TRUE) + train_path <- d[str_detect(d, pattern = "_tr")] + test_path <- d[str_detect(d, pattern = "_te")] + # Then should read in the MAE + train_data <- load_MAE(train_path, prefix="train") |> extract_Xy() + test_data <- load_MAE(test_path, prefix="test") |> extract_Xy() + sample_names <- check_common_samples(train_data) + cat("\nTotal of", length(sample_names), "samples:\n", sample_names) + + # Filenames to write out + model_file <- paste(label, "mofa_cox_survival_model.rds", sep="-") + test_file <- paste(label, "mofa_cox_survival_test_data.rds", sep="-") + + # Train a model to run inner cv or not + if (run_inner_cv) { + message("\nThe inner CV option in each fold is not implemented for MOFA survival yet\n") + #--------------------------------------------------------------------------- + } else { + message("\nNot running inner cv per fold\n") + # The input data are already MOFA embeddings, so we can directly use them to fit a glmnet model for survival + train_x <- train_data$X$embeddings + # Print head of embeddings + cat( train_x |> head() ) + train_y <- train_data$Y + + # Do an additional check if train_x only has 1 column then add a dummy 0 column + # to it, since glmnet expect X to be at least N x 2 + if (ncol(train_x) == 1) { + train_x <- cbind(train_x, dummy_zero=0) + # Now since adding dummy col, test data needs to be updated too + test_emb <- as.matrix(test_data$X$embeddings) + test_emb <- cbind(test_emb, dummy_zero=0) + test_data$X$embeddings <- test_emb + } + + # Glmnet model (ACTUALLY using this to predict) + model <- glmnet( + x = train_x, + y = train_y, + family = "cox" + ) + } + # ========================================= + message("\nSaving files to", label, "\n") + # Write out to disk + # Saving the test data for later use + # The mofa model hdf5 is saved once its training is done + saveRDS(object = test_data, file=test_file) + saveRDS(object = model, file=model_file) + + model_file <- paste(opt$label, "mofa_cox_survival_model.rds", sep = "-") + test_file <- paste(opt$label, paste0("mofa_cox_survival_test_data.rds"), sep = "-") + saveRDS(model, model_file) + saveRDS(list(blocks = d_tr$blocks, time = d_tr$time, status = d_tr$status, + test_mae_path = test_path), test_file) + cat("Saved", model_file, "and", test_file, "\n") + + return(model) +} +# Call the function here +main(mae_path = opt$mae_path, + label = opt$label, + fold_path = opt$fold_path, + prefix = opt$prefix, + run_inner_cv = opt$run_inner_cv +) +cat("Done") + + +train_survival_mofa <- function(blocks, y, n_factors = 10) { + reticulate::use_virtualenv("/workspace/.venv", required = FALSE) + blocks_t <- lapply(blocks, t) # MOFA2 expects features x samples + mofa <- MOFA2::create_mofa(blocks_t) + data_opts <- MOFA2::get_default_data_options(mofa) + model_opts <- MOFA2::get_default_model_options(mofa) + model_opts$num_factors <- n_factors + train_opts <- MOFA2::get_default_training_options(mofa) + train_opts$maxiter <- 1000 + mofa <- MOFA2::prepare_mofa(mofa, data_options = data_opts, + model_options = model_opts, training_options = train_opts) + model <- MOFA2::run_mofa(mofa, outfile = tempfile(fileext = ".hdf5")) + Z <- MOFA2::get_factors(model, factors = "all")[[1]] + fit <- glmnet::cv.glmnet(as.matrix(Z), Surv(y$time, y$status), family = "cox", alpha = 0.5) + # Per-factor linear calibration for manual projection of new data + calib <- lapply(seq_len(ncol(Z)), function(j) lm(Z[, j] ~ 0)$coef) + feat_means <- lapply(blocks, colMeans) + W <- MOFA2::get_weights(model) + list(method = "mofa_cox", model = model, fit = fit, calib = calib, + feat_means = feat_means, W = W, factor_names = colnames(Z)) +} diff --git a/modules/mogonet/select_feature/main.nf b/modules/mogonet/select_feature/main.nf index 3c16934..0ee189d 100644 --- a/modules/mogonet/select_feature/main.nf +++ b/modules/mogonet/select_feature/main.nf @@ -23,8 +23,9 @@ process MOGONET_SELECT_FEATURE { tuple val(dataset_name), path(mu_path) val(he_base_dim) output: - path("*.csv"), emit: features - path("*.log"), emit: log + path("*.csv"), emit: features + path("*.json"), emit: hyperparameters + path("*.log"), emit: log script: // The selected features, note have to run preprocess inside this step diff --git a/modules/mogonet/select_feature/resources/usr/bin/mogonet_select_features.py b/modules/mogonet/select_feature/resources/usr/bin/mogonet_select_features.py index d3a88fa..76c2639 100755 --- a/modules/mogonet/select_feature/resources/usr/bin/mogonet_select_features.py +++ b/modules/mogonet/select_feature/resources/usr/bin/mogonet_select_features.py @@ -21,11 +21,111 @@ import mudata import hashlib import copy +import json # Custom functions import from prepare_mogonet_single_split import prepare_mogonet_feat_select from feature_importance import cal_feat_imp, summarize_imp_feat +def make_json_serializable(value): + """Convert common NumPy values into JSON-compatible Python values.""" + + if isinstance(value, np.generic): + return value.item() + + if isinstance(value, np.ndarray): + return value.tolist() + + if isinstance(value, dict): + return { + str(key): make_json_serializable(item) + for key, item in value.items() + } + + if isinstance(value, (list, tuple)): + return [make_json_serializable(item) for item in value] + + return value + + +def write_selected_hyperparameters( + dataset_name, + he_base_dim, + adj_parameter, + lr_e_pretrain, + lr_e, + lr_c, + num_epoch_pretrain, + num_epoch, + num_class, + random_state, +): + method = "mogonet" + + result = { + "schema_version": "1.0", + "run_id": f"{method}-{dataset_name}", + "dataset": dataset_name, + "method": method, + "analysis_stage": "model_selection", + "selection": { + "strategy": "fixed_configuration", + "hyperparameter_tuning": False, + }, + "parameters": { + "he_base_dim": { + "value": he_base_dim, + "treatment": "fixed", + }, + "adj_parameter": { + "value": adj_parameter, + "treatment": "fixed", + }, + "lr_e_pretrain": { + "value": lr_e_pretrain, + "treatment": "fixed", + }, + "lr_e": { + "value": lr_e, + "treatment": "fixed", + }, + "lr_c": { + "value": lr_c, + "treatment": "fixed", + }, + "num_epoch_pretrain": { + "value": num_epoch_pretrain, + "treatment": "fixed", + }, + "num_epoch": { + "value": num_epoch, + "treatment": "fixed", + }, + "num_class": { + "value": num_class, + "treatment": "fixed", + }, + "random_state": { + "value": random_state, + "treatment": "fixed", + }, + }, + } + + output_path = f"{method}-{dataset_name}_selected_hyperparameters.json" + + with open(output_path, "w", encoding="utf-8") as handle: + json.dump( + result, + handle, + indent=2, + ensure_ascii=False, + ) + + print(f"Selected hyperparameters written to: {output_path}") + + return output_path + # def get_view_list(data_folder): # # First list all *_featname.csv in current data folder # pattern = "_featname.csv" @@ -60,7 +160,15 @@ def seed_everything(seed: int): # TODO: The feat importance is always 0? # and make sure to have topn to be big number # like 10% of datasize -def main(mu_path, dataset_name, n_percent, random_state=123, block_num=0, test_size=0.25, reps=5, num_class=2, he_base_dim=2, adj_parameter=5, num_epoch=200): +def main(mu_path, dataset_name, n_percent, random_state=123, block_num=0, test_size=0.25, reps=5, + num_class=2, + he_base_dim=2, + adj_parameter=5, + lr_e_pretrain=1e-3, + lr_e=5e-4, + lr_c=1e-3, + num_epoch_pretrain=50, + num_epoch=200): """ Parameters ---------- @@ -103,6 +211,19 @@ def main(mu_path, dataset_name, n_percent, random_state=123, block_num=0, test_s filename = f"{method}-{dataset_name}_features_selected.csv" # And write it to file feats_df.to_csv(filename, index=False) + # Also write out the hyperparams + write_selected_hyperparameters( + dataset_name=dataset_name, + he_base_dim=he_base_dim, + adj_parameter=adj_parameter, + lr_e_pretrain=lr_e_pretrain, + lr_e=lr_e, + lr_c=lr_c, + num_epoch_pretrain=num_epoch_pretrain, + num_epoch=num_epoch, + num_class=num_class, + random_state=random_state, + ) return(feats_df) diff --git a/modules/prepare_data/prepare_mae_data/main.nf b/modules/prepare_data/prepare_mae_data/main.nf index 3e0b452..46ac1f5 100644 --- a/modules/prepare_data/prepare_mae_data/main.nf +++ b/modules/prepare_data/prepare_mae_data/main.nf @@ -23,6 +23,7 @@ process PREPARE_MAE_DATA { input: tuple val(dataset_name), path(mae_path) val(filter_low_var) + val(outcome_type) output: tuple val(dataset_name), path("${dataset_name}*processed*mae_data"), emit: mae_data // Directory containing MultiAssayExperiment tuple val(dataset_name), path("*.log"), emit:log @@ -30,7 +31,8 @@ process PREPARE_MAE_DATA { """ transform_mae_format.R --mae_path=${mae_path} \ --dataset_name=${dataset_name} \ - --filter_low_var=${filter_low_var} > \ + --filter_low_var=${filter_low_var} \ + --outcome_type=${outcome_type} > \ ${dataset_name}.log """ } diff --git a/modules/prepare_data/prepare_mae_data/resources/usr/bin/save_mae.R b/modules/prepare_data/prepare_mae_data/resources/usr/bin/save_mae.R index e53687d..dacc3ff 100755 --- a/modules/prepare_data/prepare_mae_data/resources/usr/bin/save_mae.R +++ b/modules/prepare_data/prepare_mae_data/resources/usr/bin/save_mae.R @@ -1,19 +1,41 @@ # Use this function instead to save MAE # We have checked formats already -save_mae <- function(object, dataset_name, prefix) { +save_mae <- function(object, dataset_name, prefix, outcome_type="classification") { # Extract from the list with blocks and metadata blocks <- object$blocks - metadata <- data.frame(response = object$response) + + if (outcome_type == "classification") { + # For classification, we only need the response column + metadata <- data.frame(response = object$y) + } else if (outcome_type == "survival") { + # For survival, we need both time and status columns + # y here would be df for survival + # So it already contains the df + metadata <- data.frame(response = object$y$response, + time = object$y$time, + status = object$y$status) + message("\nHead of time: ", head(metadata$time)) + message("\nHead of status: ", head(metadata$status)) + + } else { + stop("Unsupported outcome_type: ", outcome_type, + ". Expected 'classification' or 'survival'.", + call. = FALSE) + } + + # Lastly assign the sample names to rownames of colData + rownames(metadata) <- colnames(object$blocks[[1]]) # Note, this assumes all blocks have the same sample names + # Construct new mae # TODO: Fix this or make it more robust # NOTE: this metadata is solely the response, others are discard mae <- MultiAssayExperiment::MultiAssayExperiment( experiments = blocks, - metadata = metadata + colData = metadata ) # Note, the delayed matrix is affecting the subset of metadata # so manually add response here - mae$response <- metadata + mae$response <- metadata$response # Save to HDF5 format # Hardcode this prefix MultiAssayExperiment::saveHDF5MultiAssayExperiment( diff --git a/modules/prepare_data/prepare_mae_data/resources/usr/bin/transform_mae_format.R b/modules/prepare_data/prepare_mae_data/resources/usr/bin/transform_mae_format.R index 239be01..e0562b9 100755 --- a/modules/prepare_data/prepare_mae_data/resources/usr/bin/transform_mae_format.R +++ b/modules/prepare_data/prepare_mae_data/resources/usr/bin/transform_mae_format.R @@ -13,6 +13,7 @@ Options: --dataset_name=DNAME Name of dataset to provide as id [default: empty] --replace_na_val=NA_VAL Value to replace NAs inside the data [default: 0] --filter_low_var=FIL_LOW_VAR Filter low variance or not [default: 0] + --outcome_type=OUTCOME_TYPE Outcome type . One of 'classification' or 'survival'. [default: classification] " library(here) library(mixOmics) @@ -42,7 +43,7 @@ opt <- docopt::docopt(doc) # 3. Need to have a common observations names set as sample_names # 4. Also requires to supply a dataset_name -main <- function(mae_path, dataset_name, prefix="", replace_na_val=0, center=TRUE, scale=FALSE, filter_low_var=FALSE) { +main <- function(mae_path, dataset_name, prefix="", outcome_type="classification", replace_na_val=0, center=TRUE, scale=FALSE, filter_low_var=FALSE) { # Fail quickly if (mae_path == "empty") stop ("Need to provide a path to directory containing MAE") if (dataset_name == "empty") stop("Need to provide a dataset name") @@ -56,13 +57,23 @@ main <- function(mae_path, dataset_name, prefix="", replace_na_val=0, center=TRU # Another preprocessing step X <- preprocess_view(X, replace_na_val=replace_na_val, scale=scale, filter_low_var=filter_low_var) # This would transform response to factor chr - y <- check_response(y=mae$response) + + message("\nOutcome type is: ", outcome_type) + if (outcome_type == "classification") { + y <- check_response(y=mae$response) + } else if (outcome_type == "survival") { + # TODO: uggly fix here .... + y <- mae@colData |> as.data.frame() |> dplyr::select(response, time, status) |> + dplyr::mutate(response = check_response(response)) + } else { + stop("Outcome type is not supported, please use classification or survival") + } # Put together these inputs and resave - dat <- list(blocks=X, response=y) + dat <- list(blocks=X, y=y) # Save each of this to both MAE and Mu # Add prefix of processed in the dataset name dname <- paste0(dataset_name, "_", "processed") - new_mae <- save_mae(dat, dataset_name=dname, prefix=prefix) + new_mae <- save_mae(dat, dataset_name=dname, prefix=prefix, outcome_type=outcome_type) return(new_mae) } @@ -70,4 +81,5 @@ main <- function(mae_path, dataset_name, prefix="", replace_na_val=0, center=TRU # Execute the main function here main(mae_path=opt$mae_path, dataset_name=opt$dataset_name, replace_na_val=as.numeric(opt$replace_na_val), +outcome_type=opt$outcome_type, filter_low_var=as.logical(as.numeric(opt$filter_low_var))) diff --git a/modules/prepare_data/prepare_mu_data/main.nf b/modules/prepare_data/prepare_mu_data/main.nf index 47a1c9f..be42f27 100644 --- a/modules/prepare_data/prepare_mu_data/main.nf +++ b/modules/prepare_data/prepare_mu_data/main.nf @@ -22,13 +22,15 @@ process PREPARE_MU_DATA { input: tuple val(dataset_name), path(mu_path) val(filter_low_var) + val(outcome_type) output: tuple val(dataset_name), path("${dataset_name}*processed*.h5mu"), emit: mu_data // Path to processed h5mu file tuple val(dataset_name), path("*.log"), emit: log """ transform_mudata_format.py --mu_path=${mu_path} \ --dataset_name=${dataset_name} \ - --filter_low_var=${filter_low_var} > \ + --filter_low_var=${filter_low_var} \ + --outcome_type=${outcome_type} > \ ${dataset_name}.log """ } diff --git a/modules/prepare_data/prepare_mu_data/resources/usr/bin/transform_mudata_format.py b/modules/prepare_data/prepare_mu_data/resources/usr/bin/transform_mudata_format.py index c4bd6bd..52d119f 100755 --- a/modules/prepare_data/prepare_mu_data/resources/usr/bin/transform_mudata_format.py +++ b/modules/prepare_data/prepare_mu_data/resources/usr/bin/transform_mudata_format.py @@ -14,6 +14,7 @@ --var_threshold=VAR_THRES Threhold for variance to filter features from [default: 0.16] --replace_na_val=NA_VAL Value to replace NANs in omics [default: 0] --filter_low_var=FIL_LOW_VAR Filter low variance or not [default: 0] + --outcome_type=OUTCOME_TYPE Outcome type, one of 'classification' or 'survival' [default: classification] """ # import libraries @@ -37,7 +38,7 @@ def transform_mudata_format( - mu_path, dataset_name, identifier_col="sample_name", + mu_path, dataset_name, identifier_col="sample_name", outcome_type="classification", var_threshold=0.16, replace_na_val=0, scale=False, convert_to="categorical", filter_low_var=False @@ -48,7 +49,7 @@ def transform_mudata_format( data = raw_data.copy() # Mudata has stricter format, so only need to check if each obs contain the # required columns - accepted_cols = set(["response", "sample_names", "sample_name"]) + accepted_cols = set(["response", "sample_names", "sample_name", "sample"]) new_mu_dict = {} for modality in data.mod: @@ -109,6 +110,7 @@ def transform_mudata_format( # TODO: Remove the replace na val as its not doing anything here transform_mudata_format( mu_path=args['--mu_path'] , dataset_name=args['--dataset_name'], + outcome_type=args['--outcome_type'], var_threshold=float(args['--var_threshold']), replace_na_val=float(args['--replace_na_val']), filter_low_var=bool(int(args['--filter_low_var'])) - ) \ No newline at end of file + ) diff --git a/modules/prepare_data/split_modality/main.nf b/modules/prepare_data/split_modality/main.nf new file mode 100644 index 0000000..c13e162 --- /dev/null +++ b/modules/prepare_data/split_modality/main.nf @@ -0,0 +1,36 @@ +/* + Process to expand a MuData into one single-modality MuData per + omics. Part of the single-modality mode (params.single_modality_mode). + + Input : tuple val(dataset_name), path("_processed.h5mu") + Output: tuple val(dataset_name), path("-_processed.h5mu" ...) + NOTE: the val is still the PARENT dataset name; the subworkflow + derives the child (unimodal) dataset names from the file names via + flatMap, because a process output tuple cannot fan out per file. +*/ + +process SPLIT_MODALITY { + tag "${dataset_name}" + label 'process_low' + label 'generic' + publishDir ( + path: "${params.outdir}/${task.process.tokenize(':').join('/').toLowerCase()}/${dataset_name}", + mode: 'copy', + overwrite: true + ) + + input: + tuple val(dataset_name), path(mu_path) + + output: + tuple val(dataset_name), path("${dataset_name}-*_processed.h5mu"), emit: unimodal_files + tuple val(dataset_name), path("*.log"), emit: log + + script: + """ + split_modalities.py \ + --mu_path=${mu_path} \ + --dataset_name=${dataset_name} > \ + ${dataset_name}-split_modality.log + """ +} diff --git a/modules/prepare_data/split_modality/resources/usr/bin/split_modalities.py b/modules/prepare_data/split_modality/resources/usr/bin/split_modalities.py new file mode 100755 index 0000000..16d2d38 --- /dev/null +++ b/modules/prepare_data/split_modality/resources/usr/bin/split_modalities.py @@ -0,0 +1,77 @@ +#!/usr/bin/env python +""" +Split a MuData into one single-modality MuData per view. + +Used by the single-modality mode (params.single_modality_mode): each output is a +first-class dataset named '{dataset_name}-{modality}' so that downstream +splitting / cross-validation / methods run UNCHANGED. + +Requirements on the input (already enforced by transform_mudata_format.py): + - each modality's obs contains 'response' (yes/no) and 'sample_name' + +Usage: + split_modalities.py [options] + +Options: + --mu_path=MU_PATH Path to the processed MuData [default: empty] + --dataset_name=DNAME Parent dataset name [default: empty] +""" + +from docopt import docopt +import mudata +import re + +def make_safe_id(value): + """Safe for filenames, Nextflow work outputs, and shell paths.""" + value = str(value).strip() + value = re.sub(r"\s+", "_", value) # spaces → _ + value = re.sub(r"[^A-Za-z0-9._-]", "_", value) # other unsafe chars → _ + value = re.sub(r"_+", "_", value) # collapse ___ + return value.strip("._-") + + +def split_modality(mdata, mod_names, dataset_name=""): + out_files = [] + + safe_dataset_name = make_safe_id(dataset_name) + + for mod in mod_names: + mod_obs = mdata[mod].obs + + for col in ("response", "sample_name"): + if col not in mod_obs.columns: + raise KeyError( + f"{dataset_name}/{mod}: mod-level obs missing '{col}'" + ) + + # Keep the original modality key inside the MuData object. + sub = mudata.MuData({mod: mdata[mod].copy()}) + + safe_mod = make_safe_id(mod) + out_file = f"{safe_dataset_name}-{safe_mod}_processed.h5mu" + + sub.write(out_file) + out_files.append(out_file) + + print( + f"[INFO] Wrote {out_file} " + f"({sub.n_obs} obs × {mdata[mod].n_vars} vars)" + ) + + return out_files + + +def main(mu_path, dataset_name): + mdata = mudata.read(mu_path) + mod_names = list(mdata.mod.keys()) + if len(mod_names) < 2: + print(f"[WARN] {dataset_name} has a single modality ({mod_names}); nothing to split") + out_files = split_modality(mdata, mod_names, dataset_name) + return out_files + + +if __name__ == "__main__": + args = docopt(__doc__) + files = main(args["--mu_path"], args["--dataset_name"]) + # Print file names one per line so the wrapping process can glob them + print("\n".join(files)) diff --git a/modules/rgcca/select_feature/main.nf b/modules/rgcca/select_feature/main.nf index 963211b..3e8959b 100644 --- a/modules/rgcca/select_feature/main.nf +++ b/modules/rgcca/select_feature/main.nf @@ -22,6 +22,7 @@ process RGCCA_SELECT_FEATURE { each(design) output: path("*.csv"), emit: features + path("*.json"), emit: hyperparameters path("*plot*"), emit: plot path("*.log"), emit: log diff --git a/modules/rgcca/select_feature/resources/usr/bin/rgcca_select_features.R b/modules/rgcca/select_feature/resources/usr/bin/rgcca_select_features.R index df7cf46..562ee4e 100755 --- a/modules/rgcca/select_feature/resources/usr/bin/rgcca_select_features.R +++ b/modules/rgcca/select_feature/resources/usr/bin/rgcca_select_features.R @@ -80,9 +80,13 @@ main <- function(mae_path, dataset_name, ncomp=2, design="full", prediction_mode # It should now look like list(X1=X1, X2=X2, ... , XN=XN, response=Y) rgcca_input <- X rgcca_input[["response"]] <- as.factor(data_list$Y) - + block_names <- names(rgcca_input) # This is number of omics including the response block, so H + 1 J <- length(rgcca_input) + + # Set scheme to horst for same comparison with diablo + scheme <- "horst" + # Set up the connection matrix if (design == "full") { # Full means 1 everywhere not of diagonal, meaning every omics @@ -110,8 +114,8 @@ main <- function(mae_path, dataset_name, ncomp=2, design="full", prediction_mode blocks = rgcca_input, response = length(rgcca_input), connection = connection, method = "rgcca", + scheme = scheme, sparsity = 1, # Use the default value - #par_type = par_type, tau = 1, # Fix tau 1 for all components so gets non-zero weights for all features ncomp = ncomp, prediction_model = prediction_model, @@ -140,7 +144,7 @@ main <- function(mae_path, dataset_name, ncomp=2, design="full", prediction_mode filter(view != "response") |> as_tibble() - + # Now wrangle this to dataframe for downstream usage feats_df <- weights_df %>% @@ -163,6 +167,85 @@ main <- function(mae_path, dataset_name, ncomp=2, design="full", prediction_mode comb_name <- paste(design, dataset_name, sep="-") feats_file <- paste0(comb_name, "_", "features_selected", ".", "csv") write.csv(x=feats_df, file=feats_file, row.names=FALSE) + # ------------------------------------ + # And the hyperparameters section + # rgcca_cv stores the CV-selected parameter combination here. + selected_tau <- as.numeric(cv_out$best_params) + selected_tau_json <- as.list(stats::setNames(selected_tau, block_names)) + connection_json <- lapply( + seq_len(nrow(connection)), + function(i) { + as.list( + stats::setNames( + as.numeric(connection[i, ]), + block_names + ) + ) + } + ) + + names(connection_json) <- block_names + + method_name <- paste("rgcca", design, sep = "-") + + hyperparameters <- list( + tau = make_parameter_record( + value = selected_tau_json, + treatment = "tuned" + ), + + ncomp = make_parameter_record( + value = as.integer(ncomp), + treatment = "fixed" + ), + + sparsity = make_parameter_record( + value = 1, + treatment = "fixed" + ), + + scheme = make_parameter_record( + value = scheme, + treatment = "fixed" + ), + + design = make_parameter_record( + value = design, + treatment = "fixed" + ), + + connection = make_parameter_record( + value = connection_json, + treatment = "data-derived" + ) + ) + + selection <- list( + strategy = "rgcca_cv", + tuned_parameter = cv_out$par_type, + candidate_generation = "automatic", + candidate_sets = nrow(cv_out$params), + prediction_model = prediction_model, + validation = validation, + folds = as.integer(nfolds), + repeats = as.integer(reps), + metric = metric, + random_seed = as.integer(seed) + ) + + hyperparameter_file <- write_selected_hyperparameters( + dataset_name = dataset_name, + method_name = method_name, + parameters = hyperparameters, + selection = selection + ) + + message( + "Selected hyperparameters written to: ", + hyperparameter_file + ) + + return(feats_df) } diff --git a/modules/rgcca/train/resources/usr/bin/run_rgcca.R b/modules/rgcca/train/resources/usr/bin/run_rgcca.R index 89006d4..0261966 100755 --- a/modules/rgcca/train/resources/usr/bin/run_rgcca.R +++ b/modules/rgcca/train/resources/usr/bin/run_rgcca.R @@ -79,6 +79,9 @@ main <- function(mae_path, label, fold_path, inner_cv, prefix, method, design, n sample_names <- check_common_samples(train_data) cat("\nTotal of", length(sample_names), "samples:\n", sample_names) + # Set scheme to horst for same comparison with diablo + scheme <- "horst" + # Also make up the connection matrix based on the design chosen # one of full or null # This is number of omics including the response block, so H + 1 @@ -109,7 +112,8 @@ main <- function(mae_path, label, fold_path, inner_cv, prefix, method, design, n # use default settings # The response block is always set at the end of the list of data model <- rgcca(train_data, tau=tau, connection=connection, - method=method, response=length(train_data), ncomp=ncomp + method=method, response=length(train_data), ncomp=ncomp, + scheme=scheme ) message("\nFitted model\n") } diff --git a/modules/sklearn/predict/resources/usr/bin/combine_mdata2df.py b/modules/sklearn/predict/resources/usr/bin/combine_mdata2df.py deleted file mode 100644 index fe451cc..0000000 --- a/modules/sklearn/predict/resources/usr/bin/combine_mdata2df.py +++ /dev/null @@ -1,23 +0,0 @@ -import mudata -import pandas as pd -import numpy as np - - -def combine_mdata2df(mdata, concat=True): - # Takes in a mudata and extract X and y components - # Get the Xs as a dataframe of combining all modality together columnwise - mod_names = list(mdata.mod.keys()) - # We also add the modality in front of every feature just like "epigenomics_some_feature_name" - X_df = pd.concat( [mdata[k].to_df().add_prefix(f"{k}_") for k in mod_names], axis=1 ) - # Also extract the observation df - y_df = mdata[mod_names[0]].obs - # Should contain response column and is of string yes or no - assert y_df["response"].isin(['yes', 'no']).all(), "Column contains values other than 'yes' or 'no'" - # Converting to numeric binary - y_df.loc[:, "response"] = np.where(y_df["response"] == "yes", 1, 0) - # When supplied not to concat meaning y requires other informations more than response - if not concat: - return X_df, y_df - # Merge all DataFrames together, and drop irrevelant meta information - merged_df = pd.concat([ X_df, y_df[["response"]] ] , axis=1) - return merged_df diff --git a/modules/sklearn/predict/resources/usr/bin/sklearn_predict.py b/modules/sklearn/predict/resources/usr/bin/sklearn_predict.py index 049dbf0..15ac418 100755 --- a/modules/sklearn/predict/resources/usr/bin/sklearn_predict.py +++ b/modules/sklearn/predict/resources/usr/bin/sklearn_predict.py @@ -35,7 +35,7 @@ def main(model_path, test_path, label, method_name, target_col="response"): test_data = mudata.read(test_path) # Partion the mudata to x df and y df (which contains response and other metadata) # Y contains other metadata information, hence not doing in the form of merged way - test_X_df, test_y_df = combine_mdata2df(test_data, concat=False) + test_X_df, test_y_df, _ = combine_mdata2df(test_data) # Get predicted probabilities predicted_df = pd.DataFrame(model.predict_proba(test_X_df), columns=model.classes_) # Retrieve the class equals to 1 only (no string required) diff --git a/modules/sklearn/select_feature/main.nf b/modules/sklearn/select_feature/main.nf index 52a0b0e..ae7fa65 100755 --- a/modules/sklearn/select_feature/main.nf +++ b/modules/sklearn/select_feature/main.nf @@ -21,6 +21,7 @@ process SKLEARN_SELECT_FEATURE { each (model_name) output: path("*.csv"), emit: features + path("*.json"), emit: hyperparameters path("*.log"), emit: log script: diff --git a/modules/sklearn/select_feature/resources/usr/bin/combine_mdata2df.py b/modules/sklearn/select_feature/resources/usr/bin/combine_mdata2df.py deleted file mode 100644 index 842f8ec..0000000 --- a/modules/sklearn/select_feature/resources/usr/bin/combine_mdata2df.py +++ /dev/null @@ -1,23 +0,0 @@ -import mudata -import pandas as pd -import numpy as np - - -def combine_mdata2df(mdata, concat=True): - # Takes in a mudata and extract X and y components - # Get the Xs as a dataframe of combining all modality together columnwise - mod_names = list(mdata.mod.keys()) - # We also add the modality in front of every feature just like "epigenomics_some_feature_name" - X_df = pd.concat( [mdata[k].to_df().add_prefix(f"{k}_") for k in mod_names], axis=1 ) - # Also extract the observation df - y_df = mdata[mod_names[0]].obs - # Should contain response column and is of string yes or no - assert y_df["response"].isin(['yes', 'no']).all(), "Column contains values other than 'yes' or 'no'" - # Converting to numeric binary - y_df.loc[:, "response"] = np.where(y_df["response"] == "yes", 1, 0) - # When supplied not to concat meaning y requires other informations more than response - if not concat: - return X_df, y_df - # Merge all DataFrames together, and drop irrevelant meta information - merged_df = pd.concat([ X_df, y_df[["response"]] ] , axis=1) - return merged_df \ No newline at end of file diff --git a/modules/sklearn/select_feature/resources/usr/bin/get_feats_df.py b/modules/sklearn/select_feature/resources/usr/bin/get_feats_df.py index c42c951..0503f76 100644 --- a/modules/sklearn/select_feature/resources/usr/bin/get_feats_df.py +++ b/modules/sklearn/select_feature/resources/usr/bin/get_feats_df.py @@ -1,9 +1,21 @@ import pandas as pd +import warnings # Use this fun to extract coefficients or feature importance # from sklearn classifiers def get_feats_df(classifier, feat_names, model_name, dataset_name, modality_names): + columns = ["feature", "view", "coef", "method", "dataset_name"] + + # MLP has no native per-feature importance + if model_name.lower() == "mlp": + warnings.warn( + f"{dataset_name}/{model_name}: native feature importance " + "is unavailable; writing a header-only feature CSV." + ) + feats_df = pd.DataFrame(columns=columns) + return feats_df + # Init importance value importance = None # Check the type of classifier and extract feature importance accordingly @@ -35,9 +47,9 @@ def get_feats_df(classifier, feat_names, model_name, dataset_name, modality_name model_lower = model_name.lower().replace(" ", "_") feats_df['method'] = f"sklearn-{model_lower}" # Lastly alter it to right order and matching columns - right_order = ['feature', 'view', 'coef', 'method', 'dataset_name'] + # Arrange columns in the expected output order try: - feats_df = feats_df[right_order] + feats_df = feats_df[columns] except KeyError as e: - print(f"Sklearn select feature for '{dataset_name}', '{model_lower}' column not found: {e}") + print(f"Sklearn select feature for '{dataset_name}', '{model_lower}' column not found: {e}") return feats_df \ No newline at end of file diff --git a/modules/sklearn/select_feature/resources/usr/bin/load_classifier_class.py b/modules/sklearn/select_feature/resources/usr/bin/load_classifier_class.py deleted file mode 100644 index b9d91eb..0000000 --- a/modules/sklearn/select_feature/resources/usr/bin/load_classifier_class.py +++ /dev/null @@ -1,113 +0,0 @@ -import scipy.stats as stats -import importlib - - -def load_classifier_class(model_name, random_state=42, probability=True): - # Define models with their respective classes, default parameters, and distributions - # Adopted from https://scikit-learn.org/stable/auto_examples/classification/plot_classifier_comparison.html - # NOTE: keep to use SVC rather than LinearSVC class since the latter do not have predict_prob - # https://stackoverflow.com/questions/33843981/under-what-parameters-are-svc-and-linearsvc-in-scikit-learn-equivalent - - # Common params - C = stats.loguniform(1e-4, 1e4) - min_sample_leaf = stats.randint(1,6) - n_estimators = stats.randint(50, 501) - learning_rate = stats.uniform(0.01, 1.1) - max_features = ["sqrt", "log2", 100, 500, 1000, None] - # Dict to store relevant information of sklearn classifiers - - # For logistic regression, it has built-in predict proba and coef - logit_dict = { - "class_path": "sklearn.linear_model.LogisticRegression", - "default_params": {"C": 1.0, "penalty": "l2", "solver": "liblinear"}, - "params_dist": {"C": C, "penalty": ["l2"], "solver": ["liblinear"]} - } - # Use SVC with kernel linear and not LinearSVC, since the latter do not have predict_proba - linear_svm_dict ={ - "class_path": "sklearn.svm.SVC", - "default_params": {"C": 1.0, "kernel": "linear", "random_state": random_state, "probability": probability}, - "params_dist": {"C": C, "kernel": ["linear"] } - } - # Decision Tree classifier, risk of overfitting, hence require pruning of trees - decision_tree_dict = { - "class_path": "sklearn.tree.DecisionTreeClassifier", - "default_params": {"max_depth": 10, "random_state": random_state}, - "params_dist": { - "criterion": ['gini', 'entropy', 'log_loss'], - "max_depth": stats.randint(5, 41), - "min_samples_leaf": min_sample_leaf, - "max_leaf_nodes": [10, 100, 1000, None] - } - } - # Random Forest classifier, should in general work better than single decision tree - random_forest_dict = { - "class_path": "sklearn.ensemble.RandomForestClassifier", - "default_params": {"n_estimators": 10, "max_features": "sqrt", "max_depth": 10, "random_state": random_state, "n_jobs": -1}, - "params_dist": { - "max_features": max_features, - "max_leaf_nodes": [10, 100, 1000, None], - "min_samples_leaf": min_sample_leaf - } - } - - # Gradient Boost classifier, should tune large number of estimators with slow learning rate - gradient_boost_dict = { - "class_path": "sklearn.ensemble.GradientBoostingClassifier", - "default_params": {"loss":'log_loss', "learning_rate":0.1, "n_estimators":10}, - "params_dist": { - "n_estimators": n_estimators, - "learning_rate": learning_rate, - "min_samples_leaf": min_sample_leaf, - "max_features": max_features - } - } - # MLPClassifier is a neural network classifier - mlp_dict = { - "class_path": "sklearn.neural_network.MLPClassifier", - "default_params": {"hidden_layer_sizes": (64,), "max_iter": 500, - "random_state": random_state}, - "params_dist": {"hidden_layer_sizes": [(32,), (64,), (128,), (64, 32)], - "alpha": stats.loguniform(1e-5, 1e-1)} - } - # Mimic xgboost - hist_gradient_boost_dict = { - "class_path": "sklearn.ensemble.HistGradientBoostingClassifier", - "default_params": {"random_state": random_state}, - "params_dist": {"learning_rate": learning_rate, - "max_leaf_nodes": [15, 31, 63], - "min_samples_leaf": stats.randint(5, 30)} - } - - # ======================== - # Lastly merge all together into a single dictionary - model_info = { - "Logit": logit_dict, - "Linear_SVM": linear_svm_dict, - # Decision Tree have risk of overfitting, hence require pruning of trees - "Decision_Tree": decision_tree_dict , - # RandomForest shuold in general work better than single decision tree - "Random_Forest": random_forest_dict, - # GradientBoost should tune large number of estimators with slow learning rate - "Gradient_Boost": gradient_boost_dict, - # MLPClassifier is a neural network classifier - "MLP": mlp_dict, - "Hist_Gradient_Boost": hist_gradient_boost_dict - } - - # Check if valid name of model was input - if model_name not in model_info: - valid_models = list(model_info.keys()) - print(f"Valid models are: {valid_models}") - raise ValueError(f"Model '{model_name}' not valid, check spelling") - - # Retrieve model info - model = model_info[model_name] - class_path = model["class_path"] - default_params = model["default_params"] - params_dist = model["params_dist"] - - # Import and return the classifier class - module_path, class_name = class_path.rsplit(".", 1) - classifier_class = getattr(importlib.import_module(module_path), class_name) - print(f"class name is: {class_name}") - return classifier_class, default_params, params_dist diff --git a/modules/sklearn/select_feature/resources/usr/bin/run_random_search_cv.py b/modules/sklearn/select_feature/resources/usr/bin/run_random_search_cv.py index 4ec6f2b..9d5fa64 100644 --- a/modules/sklearn/select_feature/resources/usr/bin/run_random_search_cv.py +++ b/modules/sklearn/select_feature/resources/usr/bin/run_random_search_cv.py @@ -8,13 +8,31 @@ def run_random_search_cv(clf_instance, X, Y, param_distributions, n_iter=10, ran # Apply random search cross validation to find optimal hyperparameters # Make sure to use a pipeline to standardize the data before fitting the model pipeline = make_pipeline(StandardScaler(), clf_instance) + # When searching over a pipeline, param keys must be prefixed with the step + # name (e.g. 'logisticregression__C'), otherwise sklearn treats them as + # invalid Pipeline params. make_pipeline names the step by the lowercased + # class name, so grab it from the last step. + step_name = pipeline.steps[-1][0] + prefixed_param_distributions = { + f"{step_name}__{k}": v for k, v in param_distributions.items() + } # Default to use all processors to speed up process clf_cv = RandomizedSearchCV( - pipeline, param_distributions=param_distributions, + pipeline, param_distributions=prefixed_param_distributions, n_iter=n_iter, random_state=random_state, - n_jobs=n_jobs + n_jobs=n_jobs, + scoring="roc_auc", + # This gives us the best model with optimal hyperparameters after fitting + refit=True ) # Could then fit this cv object of the X and Y of data search = clf_cv.fit(X, Y) - optimal_params = search.best_params_ - return optimal_params + # # Strip the step prefix so callers can re-instantiate the bare classifier + # # with the returned params (classifier_class(**optimal_params)). + # prefix = f"{step_name}__" + # optimal_params = { + # k[len(prefix):] if k.startswith(prefix) else k: v + # for k, v in search.best_params_.items() + # } + # return optimal_params + return search diff --git a/modules/sklearn/select_feature/resources/usr/bin/sklearn_select_features.py b/modules/sklearn/select_feature/resources/usr/bin/sklearn_select_features.py index 62ee2ef..f073383 100755 --- a/modules/sklearn/select_feature/resources/usr/bin/sklearn_select_features.py +++ b/modules/sklearn/select_feature/resources/usr/bin/sklearn_select_features.py @@ -21,6 +21,8 @@ import mudata import os import copy +import json +import numpy as np from sklearn.pipeline import make_pipeline from sklearn.preprocessing import StandardScaler @@ -31,6 +33,103 @@ from load_classifier_class import load_classifier_class + + +def make_json_serializable(value): + """Convert common NumPy values into JSON-compatible Python values.""" + + if isinstance(value, np.generic): + return value.item() + + if isinstance(value, np.ndarray): + return value.tolist() + + if isinstance(value, dict): + return { + str(key): make_json_serializable(item) + for key, item in value.items() + } + + if isinstance(value, (list, tuple)): + return [make_json_serializable(item) for item in value] + + return value + + +def remove_pipeline_prefix(best_params, step_name): + """ + Convert parameters such as logisticregression__C into C. + """ + + prefix = f"{step_name}__" + + return { + key[len(prefix):] if key.startswith(prefix) else key: value + for key, value in best_params.items() + } + + +def write_selected_hyperparameters( + search, + init_params, + dataset_name, + model_name, + n_iter, + random_state, + scoring="roc_auc" + ): + """ + Write the selected and intentionally fixed model parameters to JSON. + """ + + step_name = search.best_estimator_.steps[-1][0] + + selected_params = remove_pipeline_prefix( + search.best_params_, + step_name=step_name + ) + + parameter_records = {} + + # Parameters intentionally fixed in load_classifier_class(). + # Do not label a parameter as fixed if it was subsequently tuned. + for parameter, value in init_params.items(): + if parameter not in selected_params: + parameter_records[parameter] = { + "value": make_json_serializable(value), + "treatment": "fixed" + } + + for parameter, value in selected_params.items(): + parameter_records[parameter] = { + "value": make_json_serializable(value), + "treatment": "tuned" + } + + method = f"sklearn-{model_name.lower().replace(' ', '_')}" + + result = { + "dataset": dataset_name, + "method": method, + "analysis_stage": "model_selection", + "selection": { + "strategy": "randomized_search", + "metric": scoring, + "n_iter": n_iter, + "random_state": random_state, + "best_cv_score": make_json_serializable(search.best_score_) + }, + "parameters": parameter_records + } + + output_path = f"{method}-{dataset_name}_selected_hyperparameters.json" + + with open(output_path, "w", encoding="utf-8") as handle: + json.dump(result, handle, indent=2, ensure_ascii=False) + + return output_path + + # This is the main entrance of the script # TODO: You need to re-implement the main logic # 1. Perform some kind of cv to find hyperparameters @@ -60,28 +159,24 @@ def main(mu_path, dataset_name, model_name, block_num=0, n_iter=10, random_state raw_mdata = mudata.read(mu_path) # Use a copy here to avoid mixing up stuff mdata = raw_mdata.copy() - # TODO: this is uggly solution now - # Take the modality out for usage later - modality_names = list(mdata.mod.keys()) - # Then convert the mdata to merged dataframe column wise - merged_df = combine_mdata2df(mdata) - # And split them to X and Y - X_df, y_df = merged_df.drop(columns=[target_col]), merged_df[[target_col]] + # Combine the MuData into a single dataframe for X and y + X_df, y_df, modality_names = combine_mdata2df(mdata, target_col=target_col) + # sklearn expects a one-dimensional target + y = y_df[target_col].to_numpy().ravel() + # For a model , apply a CV on full data to find optimal hyperparam # then instantiate new model with best param to get feature importance or weight print(f"Model name is '{model_name}'") classifier_class, init_params, param_dist = load_classifier_class(model_name=model_name) clf_instance = classifier_class(**init_params) # Instantiate object from class with model init params # Then apply random CV on the parameter distribution of given model - optimal_params_dict = run_random_search_cv(clf_instance, X=X_df, Y=y_df, param_distributions=param_dist, n_iter=n_iter, random_state=random_state) - # Now, instantiate new instance of the model with optimal params instead - opt_clf_instance = classifier_class(**optimal_params_dict) - print(opt_clf_instance) - # Apply scaling and fit final model - opt_clf = make_pipeline(StandardScaler(), opt_clf_instance) - opt_clf.fit(X_df, y_df["response"].ravel()) + search = run_random_search_cv(clf_instance, X=X_df, Y=y, param_distributions=param_dist, n_iter=n_iter, random_state=random_state) + print("Best parameters:", search.best_params_) + print("Best CV ROC-AUC:", search.best_score_) # Extract the classifier from the pipeline - classifier = opt_clf.steps[-1][1] # Adjust this based on your pipeline's step name + # Already refitted on the complete dataset because refit=True. + optimal_pipeline = search.best_estimator_ + classifier = optimal_pipeline.steps[-1][1] # Adjust this based on your pipeline's step name # Then could either extract their weights or feature importance feats_df = get_feats_df( classifier=classifier, feat_names=X_df.columns, @@ -89,9 +184,25 @@ def main(mu_path, dataset_name, model_name, block_num=0, n_iter=10, random_state modality_names=modality_names) # Fix naming here for output, specifically add sklearn and turn it to lower method = f"sklearn-{model_name.lower().replace(' ', '_')}" - filename = f"{method}-{dataset_name}_features_selected.csv" + feature_path = f"{method}-{dataset_name}_features_selected.csv" # And write it to file - feats_df.to_csv(filename, index=False) + feats_df.to_csv(feature_path, index=False) + + # Also writing out the hyperparms selected + hyperparameter_path = write_selected_hyperparameters( + search=search, + init_params=init_params, + dataset_name=dataset_name, + model_name=model_name, + n_iter=n_iter, + random_state=random_state, + scoring="roc_auc" + ) + + print(f"Selected features written to: {feature_path}") + print(f"Selected hyperparameters written to: {hyperparameter_path}") + + return(feats_df) diff --git a/modules/sklearn/train/resources/usr/bin/combine_mdata2df.py b/modules/sklearn/train/resources/usr/bin/combine_mdata2df.py deleted file mode 100644 index 8e22d87..0000000 --- a/modules/sklearn/train/resources/usr/bin/combine_mdata2df.py +++ /dev/null @@ -1,23 +0,0 @@ -import mudata -import pandas as pd -import numpy as np - - -def combine_mdata2df(mdata, concat=True): - # Takes in a mudata and extract X and y components - # Get the Xs as a dataframe of combining all modality together columnwise - mod_names = list(mdata.mod.keys()) - # We also add the modality in front of every feature just like "epigenomics_some_feature_name" - X_df = pd.concat( [mdata[k].to_df().add_prefix(f"{k}_") for k in mod_names], axis=1 ) - # Also extract the observation df - y_df = mdata[mod_names[0]].obs - # Should contain response column and is of string yes or no - assert y_df["response"].isin(['yes', 'no']).all(), "Column contains values other than 'yes' or 'no'" - # Converting to numeric binary - y_df.loc[:, "response"] = np.where(y_df["response"] == "yes", 1, 0) - # When supplied not to concat meaning y requires other informations more than response - if not concat: - return X_df, y_df - # Merge all DataFrames together, and drop irrevelant meta information - merged_df = pd.concat([ X_df, y_df[["response"]] ] , axis=1) - return merged_df, mod_names \ No newline at end of file diff --git a/modules/sklearn/train/resources/usr/bin/sklearn_train.py b/modules/sklearn/train/resources/usr/bin/sklearn_train.py index 9d0215c..2b68a47 100755 --- a/modules/sklearn/train/resources/usr/bin/sklearn_train.py +++ b/modules/sklearn/train/resources/usr/bin/sklearn_train.py @@ -52,12 +52,10 @@ def build_reducer(X_df, mod_names, n_comp=50, random_state=42): # model_name is name of classifier to use from sklearn # See main function of available options def train(train_data, model_name, target_col="response", reduction="empty"): - # Convert the mdata to merged dataframe column wise - merged_df, mod_names = combine_mdata2df(train_data) - # Transform the mudata into X and y for sklearn + # Combine the MuData into a single dataframe for X and y + X_df, y_df, mod_names = combine_mdata2df(train_data, target_col) # X_df contains all count data from the views # y_df contains response - X_df, y_df = merged_df.drop(columns=[target_col]), merged_df[[target_col]] # Dynamically load classifier class # The third is params_dist which is for cv tuning, so discard it classifier_class, params, _ = load_classifier_class(model_name) diff --git a/modules/sksurv/predict/main.nf b/modules/sksurv/predict/main.nf new file mode 100644 index 0000000..be54ee0 --- /dev/null +++ b/modules/sksurv/predict/main.nf @@ -0,0 +1,43 @@ +// Include the parse method process name output dir +include { getPublishPath } from "${modulesDir}/functions" + +process SKSURV_PREDICT { + tag "${dataset_name}-${fold_name}-${model_name}" + debug "${params.debug}" + label 'sksurv' + + publishDir ( + path: "${params.outdir}/${getPublishPath(task.process)}/${dataset_name}/${fold_name}/${model_name}", + mode: 'copy', + overwrite: true + ) + + label 'process_low' + + input: + tuple val(dataset_name), val(fold_name), val(model_name), path(model_path) + tuple val(dataset_name), val(fold_name), val(model_name), path(test_path) + val(method_name) + + output: + tuple val(dataset_name), val(fold_name), val(method_name), path("*result*"), emit: result_table + tuple val(dataset_name), val(fold_name), val(method_name), path('*log*'), emit: log + + script: + def data_label = "${dataset_name}-${fold_name}" + """ + sksurv_predict.py \ + --model_path=${model_path} \ + --test_path=${test_path} \ + --label=${data_label} \ + --method_name=${method_name}-${model_name} > \ + ${data_label}-${model_name}-${getPublishPath(task.process).tokenize('/')[-1]}.log + """ + stub: + """ + echo ${dataset_name} + echo ${fold_name} + echo ${model_path} + echo ${test_path} + """ +} diff --git a/modules/sksurv/predict/resources/usr/bin/combine_mdata2df_survival.py b/modules/sksurv/predict/resources/usr/bin/combine_mdata2df_survival.py new file mode 100644 index 0000000..0c44ca1 --- /dev/null +++ b/modules/sksurv/predict/resources/usr/bin/combine_mdata2df_survival.py @@ -0,0 +1,74 @@ +import mudata +import numpy as np +import pandas as pd + + +# Structured array dtype expected by scikit-survival +SURV_Y_DTYPE = [("status", "?"), ("time", " risk score (higher = riskier) + - predict_survival_function(X)-> list of StepFunction S(t) + """ + + def __init__(self, model, scaler=None): + self.model = model + self.scaler = scaler + + def _X(self, X): + return self.scaler.transform(X) if self.scaler is not None else X + + def predict(self, X): + return self.model.predict(self._X(X)) + + def predict_survival_function(self, X, **kwargs): + return self.model.predict_survival_function(self._X(X), **kwargs) + + def predict_cumulative_hazard_function(self, X, **kwargs): + return self.model.predict_cumulative_hazard_function(self._X(X), **kwargs) + + +def cindex_scorer(estimator, X, y): + """Harrell's C-index scorer usable by sklearn GridSearchCV with survival y.""" + lp = estimator.predict(X) + return concordance_index_censored(y["status"], y["time"], lp)[0] + + +def build_coxnet(X, y, l1_ratio=0.5, n_alphas=10, inner_folds=3): + """Elastic-net Cox with inner-CV alpha selection on the training fold.""" + set_seed(SEED) + # Alpha path from the training fold (same convention as cv.glmnet) + path_model = CoxnetSurvivalAnalysis(l1_ratio=l1_ratio) + path_model.fit(X, y) + alphas = np.asarray(path_model.alphas_) + # Keep a log-spaced subset of the path (down to max/100) + alphas = np.unique(np.logspace( + np.log10(alphas.max()), + np.log10(max(alphas.min(), alphas.max() / 100)), + n_alphas, + )) + # stratify inner CV on the event indicator; shrink folds if events are scarce + # (sklearn StratifiedKFold cannot stratify on structured arrays directly) + n_minority = int(min(y["status"].sum(), (~y["status"]).sum())) + inner_folds = max(2, min(inner_folds, n_minority)) + skf = StratifiedKFold(n_splits=inner_folds, shuffle=True, random_state=SEED) + cv_splits = list(skf.split(X, y["status"])) + + # Manual inner CV: per-fold try/except so numerically unstable alphas + # (small end of the path) get NaN instead of aborting the search + mean_scores = {} + for a in alphas: + scores = [] + for tr_i, va_i in cv_splits: + try: + m = _make_coxnet_pipe(l1_ratio).set_params(coxnet__alphas=[a]) + m.fit(X.iloc[tr_i], y[tr_i]) + scores.append(cindex_scorer(m, X.iloc[va_i], y[va_i])) + except (ArithmeticError, ValueError): + scores.append(np.nan) + mean_scores[a] = np.nanmean(scores) if np.any(~np.isnan(scores)) else np.nan + + # Refit on the full training fold, best CV score first; skip alphas whose + # refit fails numerically or yields a degenerate all-zero model + ranked = sorted(alphas, key=lambda a: -np.nan_to_num(mean_scores[a], nan=-np.inf)) + for a in ranked: + if np.isnan(mean_scores[a]): + continue + try: + final = _make_coxnet_pipe(l1_ratio).set_params(coxnet__alphas=[a]) + final.fit(X, y) + lp = final.predict(X) + if np.std(lp) < 1e-10: + continue # fully regularized (all coefficients zero) + return SurvivalPipeline(model=final.named_steps["coxnet"], + scaler=final.named_steps["scale"]) + except (ArithmeticError, ValueError): + continue + raise RuntimeError("coxnet refit failed for all candidate alphas") + + +from sklearn.pipeline import Pipeline as _SkPipeline + + +def _make_coxnet_pipe(l1_ratio): + return _SkPipeline([ + ("scale", StandardScaler()), + # fit_baseline_model=True enables predict_survival_function (Breslow) + ("coxnet", CoxnetSurvivalAnalysis( + l1_ratio=l1_ratio, alphas=[1.0], fit_baseline_model=True)), + ]) + + +def build_rsf(X, y): + set_seed(SEED) + model = RandomSurvivalForest( + n_estimators=200, + min_samples_leaf=10, + max_features=0.2, + n_jobs=-1, + random_state=SEED, + ) + model.fit(X, y) + return SurvivalPipeline(model=model) # trees need no scaling + + +def build_gbm(X, y): + set_seed(SEED) + model = GradientBoostingSurvivalAnalysis( + loss="coxph", + n_estimators=200, + learning_rate=0.1, + subsample=0.8, + max_features=0.2, + random_state=SEED, + ) + model.fit(X, y) + return SurvivalPipeline(model=model) # trees need no scaling + + +BUILDERS = { + "coxnet": build_coxnet, + "rsf": build_rsf, + "gbm": build_gbm, +} + + +def load_survival_model_class(model_name): + """Return the builder callable for a model name (raises on unknown names).""" + key = model_name.lower().replace(" ", "_").replace("-", "_") + if key not in BUILDERS: + raise ValueError( + f"Unknown survival model '{model_name}'. Available: {sorted(BUILDERS)}" + ) + return BUILDERS[key] diff --git a/modules/sksurv/predict/resources/usr/bin/manual_set_seed.py b/modules/sksurv/predict/resources/usr/bin/manual_set_seed.py new file mode 100644 index 0000000..e8a5120 --- /dev/null +++ b/modules/sksurv/predict/resources/usr/bin/manual_set_seed.py @@ -0,0 +1,25 @@ +# ============================================================================= +# Little utilities to use here +# ============================================================================= + +import random + + +def set_seed(seed: int): + """ + Helper function to set the seed in ``random`` and ``numpy`` (and torch + when installed) for reproducible behavior. + """ + random.seed(seed) + np.random.seed(seed) + try: + import torch + + torch.manual_seed(seed) + if torch.cuda.is_available(): + torch.cuda.manual_seed_all(seed) + except ImportError: + pass + + +import numpy as np diff --git a/modules/sksurv/predict/resources/usr/bin/sksurv_predict.py b/modules/sksurv/predict/resources/usr/bin/sksurv_predict.py new file mode 100755 index 0000000..deb1c3b --- /dev/null +++ b/modules/sksurv/predict/resources/usr/bin/sksurv_predict.py @@ -0,0 +1,76 @@ +#!/usr/bin/env python + +""" +Predict per-observation survival results with a trained scikit-survival model. +Outputs a result table with one row per test observation: risk score (lp) and +predicted survival probabilities S(t) at fixed horizons (days). + +Usage: + sksurv_predict.py [options] + +Options: + -h --help Show this message + --model_path=MODEL Path to trained model on specific fold [default: empty] + --test_path=TEST_PATH Path to the test MuData h5mu [default: empty] + --label=LABEL Label of dataset and fold iteration [default: empty] + --method_name=METHOD_NAME Method name ran [default: sksurv] +""" + +from docopt import docopt +import joblib +import numpy as np +import pandas as pd + +from combine_mdata2df_survival import combine_mdata2df_survival +from generate_result_table_survival import generate_result_table +from load_survival_model_class import SURV_HORIZONS + + +def survival_at_horizons(model, X_df, horizons, t_min_train, t_max_train): + """Evaluate predicted S(t) at fixed horizons, clipped to train follow-up.""" + sfs = model.predict_survival_function(X_df) + out = np.zeros((len(X_df), len(horizons))) + for j, h in enumerate(horizons): + if h < t_min_train: + out[:, j] = 1.0 # before any observed follow-up: S = 1 + else: + h_eff = min(h, t_max_train) + out[:, j] = [float(sf(h_eff)) for sf in sfs] + # defensive sanitization: Breslow baselines can diverge numerically + return np.clip(np.nan_to_num(out, nan=0.0, posinf=1.0, neginf=0.0), 0.0, 1.0) + + +def predict_model(model, test_data, horizons=None): + """Return (lp, surv_matrix, meta_df) for a test-fold MuData.""" + horizons = horizons or SURV_HORIZONS + X_df, meta_df, _ = combine_mdata2df_survival(test_data, concat=False) + lp = model.predict(X_df) + # training-fold follow-up range stored at fit time; fall back to test range + t_max_train = getattr(model, "_t_max_train", None) or float(meta_df["time"].max()) + t_min_train = getattr(model, "_t_min_train", None) or float(meta_df["time"].min()) + surv = survival_at_horizons(model, X_df, horizons, t_min_train, t_max_train) + return lp, surv, meta_df + + +def main(model_path, test_path, label, method_name): + model = joblib.load(model_path) + import mudata + test_data = mudata.read(test_path) + lp, surv, meta_df = predict_model(model, test_data) + result_table = generate_result_table( + lp=lp, surv=surv, meta_df=meta_df, method_name=method_name, label=label + ) + result_file = f"{label}-{method_name}-result.csv" + result_table.to_csv(result_file, index=False, header=True) + print(f"Wrote {result_file}") + return result_table + + +if __name__ == "__main__": + args = docopt(__doc__) + main( + model_path=args["--model_path"], + test_path=args["--test_path"], + label=args["--label"], + method_name=args["--method_name"], + ) diff --git a/modules/sksurv/preprocess/main.nf b/modules/sksurv/preprocess/main.nf new file mode 100644 index 0000000..358b453 --- /dev/null +++ b/modules/sksurv/preprocess/main.nf @@ -0,0 +1,46 @@ +// This should be template for preprocess_methods + +/* + Use this process prepare inputs, or any necessary + transformation/preprocessing steps to train for + a specific + + Survival variant: input MuData obs must contain 'time' and 'status' + (1 = event, 0 = censored). Fold txt files contain test indices. +*/ + +include { getPublishPath } from "${modulesDir}/functions" + +process SKSURV_PREPROCESS { + debug "${params.debug}" + tag "${dataset_name}" + label 'sksurv' + label "process_low" + + publishDir ( + path: "${params.outdir}/${getPublishPath(task.process)}/${dataset_name}", + mode: 'copy', + overwrite: true + ) + + input: + tuple val(dataset_name), path(data_path), path(split_dir) + + output: + tuple val(dataset_name), path("*fold*"), emit: fold_splits + tuple val(dataset_name), path('*.log*'), emit: log + + script: + """ + sksurv_preprocess.py \ + --data_path=${data_path} \ + --split_dir=${split_dir} \ + --dataset_name=${dataset_name} > \ + ${dataset_name}-${getPublishPath(task.process).tokenize('/')[-1].toLowerCase()}.log + """ + stub: + """ + echo ${dataset_name} + touch ${dataset_name}_fold_1 + """ +} diff --git a/modules/sksurv/preprocess/resources/usr/bin/load_test_splits.py b/modules/sksurv/preprocess/resources/usr/bin/load_test_splits.py new file mode 100644 index 0000000..1e1d41d --- /dev/null +++ b/modules/sksurv/preprocess/resources/usr/bin/load_test_splits.py @@ -0,0 +1,27 @@ +import pandas as pd +import numpy as np +import os +import glob + + +def load_test_splits(splits_dir, pattern="*fold*"): + """ + Parameters: + splits_dir: Directory containing test indices pre splitted earlier + + Return: + out_df: a dataframe, with index being test_fold_i for i to k + each row is an array of index indicating which original + row of data was "test" data + """ + # Requires pattern to match multiple files + files = glob.glob(os.path.join(splits_dir, pattern)) + split_dict = {} + for index, file in enumerate(files): + # Note we need to shift these back to 0 index + # idxs = [x - 1 for x in np.loadtxt(f).astype(int)] + test_idxs = np.loadtxt(file).astype(int) # Load them as integer from numpy + label = f"test_fold_split_{index+1}" + split_dict[label] = test_idxs + out_df = pd.DataFrame([split_dict]).T.rename(columns={0:"test_indices"}) + return out_df \ No newline at end of file diff --git a/modules/sksurv/preprocess/resources/usr/bin/sksurv_preprocess.py b/modules/sksurv/preprocess/resources/usr/bin/sksurv_preprocess.py new file mode 100755 index 0000000..e1f221a --- /dev/null +++ b/modules/sksurv/preprocess/resources/usr/bin/sksurv_preprocess.py @@ -0,0 +1,53 @@ +#!/usr/bin/env python + +""" +Preprocess survival MuData for training: partition the full MuData into +per-fold directories, each containing the train and test h5mu portion, +using pre-computed fold test indices (fold_*.txt). Same contract as +sklearn_preprocess.py of the classification method. + +Usage: + sksurv_preprocess.py [options] + +Options: + -h --help Show this message + --data_path=DATA_PATH Path to MuData h5mu [default: empty] + --split_dir=SPLIT_DIR Path to directory of fold txt files [default: splits] + --dataset_name=DNAME Name of the dataset [default: empty] +""" + +from docopt import docopt +import os +import mudata + +from load_test_splits import load_test_splits +from tr_te_split_mdata import tr_te_split_mdata + + +def main(mdata_path, split_dir, dataset_name, base_dir="fold", ext="h5mu"): + mdata = mudata.read(mdata_path) + # sanity: survival outcome columns must be present + if not {"time", "status"}.issubset(mdata.obs.columns): + first_mod = list(mdata.mod.keys())[0] + if not {"time", "status"}.issubset(mdata[first_mod].obs.columns): + raise ValueError("MuData obs must contain 'time' and 'status' columns") + + test_splits_df = load_test_splits(split_dir) + for i, split in enumerate(test_splits_df.index): + outdir = f"{base_dir}_{i+1}" + os.makedirs(outdir, exist_ok=True) + train_file = f"{outdir}/{dataset_name}-train_mu.{ext}" + test_file = f"{outdir}/{dataset_name}-test_mu.{ext}" + train_mu, test_mu = tr_te_split_mdata(mdata, splits_df=test_splits_df, split=split) + train_mu.write_h5mu(train_file) + test_mu.write_h5mu(test_file) + return None + + +if __name__ == "__main__": + args = docopt(__doc__) + main( + mdata_path=args["--data_path"], + split_dir=args["--split_dir"], + dataset_name=args["--dataset_name"], + ) diff --git a/modules/sksurv/preprocess/resources/usr/bin/tr_te_split_mdata.py b/modules/sksurv/preprocess/resources/usr/bin/tr_te_split_mdata.py new file mode 100644 index 0000000..6bd2c60 --- /dev/null +++ b/modules/sksurv/preprocess/resources/usr/bin/tr_te_split_mdata.py @@ -0,0 +1,12 @@ +import mudata +import pandas as pd + +def tr_te_split_mdata(mdata, splits_df, split): + # Given the original mdata and a splits_df containing indices of test folds, partion + # the original mudata to have a train and test copy + test_idx = splits_df.loc[split].iloc[0] + # Get the rest idxs except those in current split, sorted + train_idx = sorted(splits_df[splits_df.index != split].test_indices.explode().tolist()) + # Get train and test copy + train_mu, test_mu = mdata[train_idx].copy(), mdata[test_idx].copy() + return train_mu, test_mu \ No newline at end of file diff --git a/modules/sksurv/train/main.nf b/modules/sksurv/train/main.nf new file mode 100644 index 0000000..5b99db9 --- /dev/null +++ b/modules/sksurv/train/main.nf @@ -0,0 +1,45 @@ +// Nextflow process +// Survival training: fits a scikit-survival model on the train fold. + +include { getPublishPath } from "${modulesDir}/functions" + +process SKSURV_TRAIN { + tag "${dataset_name}-${fold_path.name}-${model_name}" + label 'sksurv' + label 'process_medium' + + publishDir ( + path: "${params.outdir}/${getPublishPath(task.process)}/${dataset_name}/${fold_path.name}", + mode: 'copy', + overwrite: true + ) + + input: + tuple val(dataset_name), path(data_path), path(fold_path) + each(model_name) + + output: + tuple val(dataset_name), val(fold_path.name), val("${model_name}"), path('*model*'), emit: model + tuple val(dataset_name), val(fold_path.name), val("${model_name}"), path('*test_data*'), emit: test_data + tuple val(dataset_name), val(fold_path.name), val("${model_name}"), path('*log*'), emit: log + + script: + def data_label = "${dataset_name}-${fold_path.name}" + """ + sksurv_train.py \ + --fold_path=${fold_path} \ + --label=${data_label} \ + --model_name=${model_name} > \ + ${data_label}-${model_name}-${getPublishPath(task.process).tokenize('/')[-1].toLowerCase()}.log + echo ${data_label} > ${data_label} + """ + + stub: + """ + echo ${dataset_name} + echo ${data_path} + echo 'some text' > text.log + touch prediction.csv + touch model + """ +} diff --git a/modules/sksurv/train/resources/usr/bin/combine_mdata2df_survival.py b/modules/sksurv/train/resources/usr/bin/combine_mdata2df_survival.py new file mode 100644 index 0000000..0c44ca1 --- /dev/null +++ b/modules/sksurv/train/resources/usr/bin/combine_mdata2df_survival.py @@ -0,0 +1,74 @@ +import mudata +import numpy as np +import pandas as pd + + +# Structured array dtype expected by scikit-survival +SURV_Y_DTYPE = [("status", "?"), ("time", " risk score (higher = riskier) + - predict_survival_function(X)-> list of StepFunction S(t) + """ + + def __init__(self, model, scaler=None): + self.model = model + self.scaler = scaler + + def _X(self, X): + return self.scaler.transform(X) if self.scaler is not None else X + + def predict(self, X): + return self.model.predict(self._X(X)) + + def predict_survival_function(self, X, **kwargs): + return self.model.predict_survival_function(self._X(X), **kwargs) + + def predict_cumulative_hazard_function(self, X, **kwargs): + return self.model.predict_cumulative_hazard_function(self._X(X), **kwargs) + + +def cindex_scorer(estimator, X, y): + """Harrell's C-index scorer usable by sklearn GridSearchCV with survival y.""" + lp = estimator.predict(X) + return concordance_index_censored(y["status"], y["time"], lp)[0] + + +def build_coxnet(X, y, l1_ratio=0.5, n_alphas=10, inner_folds=3): + """Elastic-net Cox with inner-CV alpha selection on the training fold.""" + set_seed(SEED) + # Alpha path from the training fold (same convention as cv.glmnet) + path_model = CoxnetSurvivalAnalysis(l1_ratio=l1_ratio) + path_model.fit(X, y) + alphas = np.asarray(path_model.alphas_) + # Keep a log-spaced subset of the path (down to max/100) + alphas = np.unique(np.logspace( + np.log10(alphas.max()), + np.log10(max(alphas.min(), alphas.max() / 100)), + n_alphas, + )) + # stratify inner CV on the event indicator; shrink folds if events are scarce + # (sklearn StratifiedKFold cannot stratify on structured arrays directly) + n_minority = int(min(y["status"].sum(), (~y["status"]).sum())) + inner_folds = max(2, min(inner_folds, n_minority)) + skf = StratifiedKFold(n_splits=inner_folds, shuffle=True, random_state=SEED) + cv_splits = list(skf.split(X, y["status"])) + + # Manual inner CV: per-fold try/except so numerically unstable alphas + # (small end of the path) get NaN instead of aborting the search + mean_scores = {} + for a in alphas: + scores = [] + for tr_i, va_i in cv_splits: + try: + m = _make_coxnet_pipe(l1_ratio).set_params(coxnet__alphas=[a]) + m.fit(X.iloc[tr_i], y[tr_i]) + scores.append(cindex_scorer(m, X.iloc[va_i], y[va_i])) + except (ArithmeticError, ValueError): + scores.append(np.nan) + mean_scores[a] = np.nanmean(scores) if np.any(~np.isnan(scores)) else np.nan + + # Refit on the full training fold, best CV score first; skip alphas whose + # refit fails numerically or yields a degenerate all-zero model + ranked = sorted(alphas, key=lambda a: -np.nan_to_num(mean_scores[a], nan=-np.inf)) + for a in ranked: + if np.isnan(mean_scores[a]): + continue + try: + final = _make_coxnet_pipe(l1_ratio).set_params(coxnet__alphas=[a]) + final.fit(X, y) + lp = final.predict(X) + if np.std(lp) < 1e-10: + continue # fully regularized (all coefficients zero) + return SurvivalPipeline(model=final.named_steps["coxnet"], + scaler=final.named_steps["scale"]) + except (ArithmeticError, ValueError): + continue + raise RuntimeError("coxnet refit failed for all candidate alphas") + + +from sklearn.pipeline import Pipeline as _SkPipeline + + +def _make_coxnet_pipe(l1_ratio): + return _SkPipeline([ + ("scale", StandardScaler()), + # fit_baseline_model=True enables predict_survival_function (Breslow) + ("coxnet", CoxnetSurvivalAnalysis( + l1_ratio=l1_ratio, alphas=[1.0], fit_baseline_model=True)), + ]) + + +def build_rsf(X, y): + set_seed(SEED) + model = RandomSurvivalForest( + n_estimators=200, + min_samples_leaf=10, + max_features=0.2, + n_jobs=-1, + random_state=SEED, + ) + model.fit(X, y) + return SurvivalPipeline(model=model) # trees need no scaling + + +def build_gbm(X, y): + set_seed(SEED) + model = GradientBoostingSurvivalAnalysis( + loss="coxph", + n_estimators=200, + learning_rate=0.1, + subsample=0.8, + max_features=0.2, + random_state=SEED, + ) + model.fit(X, y) + return SurvivalPipeline(model=model) # trees need no scaling + + +BUILDERS = { + "coxnet": build_coxnet, + "rsf": build_rsf, + "gbm": build_gbm, +} + + +def load_survival_model_class(model_name): + """Return the builder callable for a model name (raises on unknown names).""" + key = model_name.lower().replace(" ", "_").replace("-", "_") + if key not in BUILDERS: + raise ValueError( + f"Unknown survival model '{model_name}'. Available: {sorted(BUILDERS)}" + ) + return BUILDERS[key] diff --git a/modules/sksurv/train/resources/usr/bin/load_tr_te.py b/modules/sksurv/train/resources/usr/bin/load_tr_te.py new file mode 100644 index 0000000..174ba42 --- /dev/null +++ b/modules/sksurv/train/resources/usr/bin/load_tr_te.py @@ -0,0 +1,15 @@ +import os +import mudata + +# Given a directory to fold, that contains train_mu.h5mu +# and test_mu.h5mu, load these in memory +def load_tr_te(fold_dir): + # Find the train mu and test mu inside the fold path + d = os.listdir(fold_dir) + # Retrieve the paths to train data and test data + train_path = next((os.path.join(fold_dir, path) for path in d if 'train' in path), None) + test_path = next((os.path.join(fold_dir, path) for path in d if 'test' in path), None) + # Load as mudata then transform as dataframe + train_data = mudata.read(train_path) + test_data = mudata.read(test_path) + return train_data, test_data \ No newline at end of file diff --git a/modules/sksurv/train/resources/usr/bin/manual_set_seed.py b/modules/sksurv/train/resources/usr/bin/manual_set_seed.py new file mode 100644 index 0000000..e8a5120 --- /dev/null +++ b/modules/sksurv/train/resources/usr/bin/manual_set_seed.py @@ -0,0 +1,25 @@ +# ============================================================================= +# Little utilities to use here +# ============================================================================= + +import random + + +def set_seed(seed: int): + """ + Helper function to set the seed in ``random`` and ``numpy`` (and torch + when installed) for reproducible behavior. + """ + random.seed(seed) + np.random.seed(seed) + try: + import torch + + torch.manual_seed(seed) + if torch.cuda.is_available(): + torch.cuda.manual_seed_all(seed) + except ImportError: + pass + + +import numpy as np diff --git a/modules/sksurv/train/resources/usr/bin/sksurv_train.py b/modules/sksurv/train/resources/usr/bin/sksurv_train.py new file mode 100755 index 0000000..ad84b4c --- /dev/null +++ b/modules/sksurv/train/resources/usr/bin/sksurv_train.py @@ -0,0 +1,59 @@ +#!/usr/bin/env python + +""" +Train a scikit-survival model on the train fold of a MuData split. +Output is a fitted model (joblib) on the train portion, plus the test MuData +passed downstream — same contract as sklearn_train.py of the classification +method. + +Usage: + sksurv_train.py [options] + +Options: + -h --help Show this message + --fold_path=FOLD_PATH Directory containing one split (train/test h5mu) + --label=LABEL Label of dataset and fold iteration [default: empty] + --model_name=MOD Survival model to run: coxnet | rsf | gbm [default: coxnet] +""" + +from docopt import docopt +import joblib + +from load_tr_te import load_tr_te +from combine_mdata2df_survival import combine_mdata2df_survival +from load_survival_model_class import load_survival_model_class + + +def train_model(train_data, model_name): + """Fit a survival model on a train-fold MuData. Returns the fitted pipeline.""" + X_df, mod_names = None, None + merged_df, mod_names = combine_mdata2df_survival(train_data) + X_df = merged_df.drop(columns=["time", "status"]) + y = combine_mdata2df_survival(train_data, concat=False)[2] + builder = load_survival_model_class(model_name) + model = builder(X_df, y) + # store training-fold follow-up so prediction clips S(t) horizons correctly + model._t_min_train = float(y["time"].min()) + model._t_max_train = float(y["time"].max()) + print(f"\n[{model_name}] fitted on {X_df.shape[0]} samples x {X_df.shape[1]} features") + return model + + +def main(fold_path, label, model_name): + train_data, test_data = load_tr_te(fold_path) + model = train_model(train_data, model_name) + model_file = f"{label}-{model_name}-model.pkl" + joblib.dump(model, model_file) + test_file = f"{label}-test_data.h5mu" + test_data.write(test_file) + print(f"Saved {model_file} and {test_file}") + return model + + +if __name__ == "__main__": + args = docopt(__doc__) + main( + fold_path=args["--fold_path"], + label=args["--label"], + model_name=args["--model_name"], + ) diff --git a/modules/split_train_test/resources/usr/bin/split_tr_te.py b/modules/split_train_test/resources/usr/bin/split_tr_te.py index faa070c..7e5d04b 100755 --- a/modules/split_train_test/resources/usr/bin/split_tr_te.py +++ b/modules/split_train_test/resources/usr/bin/split_tr_te.py @@ -16,7 +16,7 @@ --num_splits=NUM_SPLITS Number of splits to generate [default: 10] --seed=SEED Random number seed to reproduce [default: 329] --outcome_type=OUTCOME_TYPE Type of outcome: classification or survival [default: classification] - --event_col=EVENT_COL Column name for survival event (survival only) [default: os_event] + --event_col=EVENT_COL Column name for survival event (survival only) [default: status] --output_dir=OUT_DIR Output folder to write fold ids txt [default: splits] --split_txt_name=SNAME Name of individual fold txt file [default: fold] """ @@ -40,7 +40,8 @@ class SplitConfig: seed: int identifier_col: str = "sample_name" outcome_type: str = "classification" - event_col: str = "event" + event_col: str = "status" + response_col: str = "response" output_dir: str = "splits" split_txt_name: str = "fold" @@ -48,7 +49,7 @@ class SplitConfig: # ------------------------------------------------------------------------------ # Utility functions (functional, stateless) # ------------------------------------------------------------------------------ -def load_mudata(path, identifier_col, outcome_type="classification", event_col="os_event"): +def load_mudata(path, identifier_col, outcome_type="classification", event_col="status", response_col="response"): mdata = mudata.read(path) block_key = list(mdata.mod.keys())[0] block = mdata.mod[block_key] @@ -59,7 +60,7 @@ def load_mudata(path, identifier_col, outcome_type="classification", event_col=" if outcome_type == "survival": y = block.obs[event_col].values else: - y = block.obs["response"].values + y = block.obs[response_col].values return X, y, groups @@ -120,7 +121,7 @@ def main(args): k=int(args["--num_splits"]), seed=int(args["--seed"]), identifier_col=args.get("--identifier_col", "sample_name"), - event_col=args.get("--event_col", "os_event"), + event_col=args.get("--event_col", "status"), outcome_type=args.get("--outcome_type", "classification"), output_dir=args["--output_dir"], split_txt_name=args["--split_txt_name"], diff --git a/nextflow.config b/nextflow.config index e5d556c..074a12d 100644 --- a/nextflow.config +++ b/nextflow.config @@ -58,6 +58,7 @@ params { skip_mofa = false skip_goat = true skip_sgmr = true + skip_sksurv = true skip_sklearn = true /*Parameters for methods*/ @@ -72,9 +73,15 @@ params { // or "sgkf" (stratidied group k fold ) "logo" (leave one group out) k_fold_number = 5 // Number of folds to generate (omitted if split type = logo) split_dir = "splits" // Directory to splitted indices - outcome_type = "classification" + outcome_type = "classification" // "classification" or "survival". Survival requires + // colData/metadata columns: time (numeric days) and status (0/1) + /* Survival task parameters (used when outcome_type = "survival") */ + time_col = "time" // colData column with follow-up time in days + event_col = "status" // colData column with event indicator (0/1) /* Threshold to convert the predicted probability to class */ threshold = 0.5 + + /* For feature selction workflow */ selectFeature = true // Default runs feature selection on all possible methods n_percent = 10 // N percent of features to be selected @@ -82,11 +89,14 @@ params { // DIABLO diablo_design_connection = ['full', 'null'] // connection for design matrix // SKLEARN + single_modality_mode = false // Available classifiers - - //sklearn_classifier_names = ['Linear_SVM', 'Logit', 'Decision_Tree', 'Gradient_Boost', 'Random_Forest', 'MLP', 'Hist_Gradient_Boost'] - sklearn_classifier_names = ['Logit', 'MLP', 'Decision_Tree'] - sklearn_reduction = ["empty", "pca50"] + sklearn_classifier_names = ['Logit', 'Random_Forest', 'MLP', 'XGBoost'] + sklearn_reduction = ["empty", "pca50"] // Preprocessing done to data + // Survival model stuff + survival_model_names = ["rsf"] + // If sklearn wants to do single input mode, default: false + sklearn_single_mode = false // For preprocessing data filter_low_var = "0" // This treats as false in R, one of "1" or "0" // Resource parameters @@ -111,6 +121,7 @@ env { subworkflowDir = "$projectDir/subworkflows" modulesDir = "$projectDir/modules" configDir = "$projectDir/conf" + PYTHONPATH = "${projectDir}/bin/python_utils" } // Load base.config by default for all pipelines @@ -136,9 +147,10 @@ profiles { } apptainer { - enabled = true - autoMounts = true - cacheDir = params.apptainer_cache_dir + apptainer.enabled = true + apptainer.autoMounts = true + apptainer.cacheDir = params.apptainer_cache_dir + apptainer.libraryDir = params.apptainer_cache_dir conda.enabled = false docker.enabled = false singularity.enabled = false @@ -162,6 +174,7 @@ profiles { includeConfig 'conf/sockeye.config' } + /* These two profiles CANNOT be used together simultaneously One for real data @@ -175,6 +188,13 @@ profiles { includeConfig 'conf/simulated_data.config' } + // For survival task + run_survival { + includeConfig 'conf/run_survival.config' + } + + // Legacy ones + run_mogonet_only { includeConfig 'conf/run_mogonet_only.config' } @@ -209,7 +229,8 @@ process { withLabel: rgcca { container = 'tonyliang19/rgcca:latest' } withLabel: codia { container = 'tonyliang19/codia:latest' } withLabel: integrao { container = 'tonyliang19/integrao:latest' } - withLabel: sklearn { container = 'tonyliang19/mogonet:latest' } // fix: should have a new sklearn container instead, for now use mogonet one + withLabel: sklearn { container = 'tonyliang19/eipy:latest' } // fix: should have a new sklearn container instead, for now use eipy bcz it contains XGBoost + withLabel: sksurv { container = 'tonyliang19/sksurv:latest' } } @@ -236,13 +257,13 @@ dag { // Meta data information manifest { - name = 'tonyliang19/multi-omics-pipeline' + name = 'CompBio-Lab/MESSI-pipeline' author = """Tony Liang""" version = '1.0.0dev' - homePage = 'https://github.com/tonyliang19/multi-omics-pipeline' + homePage = 'https://github.com/CompBio-Lab/MESSI-pipeline' description = """Benchmark Multiomics integration methods over single cell and spatial bulk data""" mainScript = 'main.nf' - nextflowVersion = '>=22.10.7' // version could be changed to 23.04.3 + nextflowVersion = '>=24.10.2' // Minimal version to work } // Load modules.config for DSL2 module specific options diff --git a/subworkflows/cross_validation/main.nf b/subworkflows/cross_validation/main.nf index a7fcb44..1299c88 100644 --- a/subworkflows/cross_validation/main.nf +++ b/subworkflows/cross_validation/main.nf @@ -30,7 +30,8 @@ def shouldRunPython() { return !( params.skip_mogonet && params.skip_integrao && - params.skip_sklearn + params.skip_sklearn && + params.skip_sksurv ) } @@ -42,7 +43,8 @@ def shouldRunR() { params.skip_diablo && params.skip_mofa && params.skip_rgcca && - params.skip_caret_multimodal + +params.skip_caret_multimodal ) } @@ -60,20 +62,25 @@ def lang = "all_langs" // ============================ workflow CROSS_VALIDATION { /* - This is a complex workflow, by training and cross validating at the same time + This is a complex workflow, by training and cross validating at the same t ime Moreover, this is ongoing for all methods, so it would crazily hard to fix due to its parallelization */ // runAllMethods = params.runAllMethods // runPython = params.runPython // runR = params.runR + outcome_type = params.outcome_type // inputs of workflow take: // datasets mae_data mu_data splits_indices - + // NEW (single-modality mode): channel of [unimodal_dataset_name, mu_path]. + // Routed to Python methods ONLY (currently sklearn), never to R methods + // or multi-block integration methods. Empty channel when + // params.single_modality_mode is false. + mu_data_unimodal main: // runAllMethods = params.runAllMethods // Determine if should run python or R or both @@ -110,8 +117,17 @@ workflow CROSS_VALIDATION { } .set { ch_full_exp } + + // ch_full_exp.mae_copy.count().view {"Total is ${it}"} // ch_full_exp.mae_copy.view() + + // The unimodal ones goes here + mu_data_unimodal + .join( splits_indices, by: 0 ) // The first element of the tuple is the dataset name, so join by that + .map { dname, mu_path, indices -> [ dname, mu_path, indices ] } + .set { unimodal_exp } + // Note, this cannot be viewed, you need ch_data.mae.view() or ch_data.mu.view() // Execute these two (they should be in parallel) // Execute language workflows or specific only @@ -123,16 +139,17 @@ workflow CROSS_VALIDATION { } else if ( runPython ) { log.info "Running methods in Python only" - CV_PYTHON ( ch_full_exp.mudata_copy ) + CV_PYTHON ( ch_full_exp.mudata_copy, unimodal_exp ) csv_results = CV_PYTHON.out.csv_results } else { log.info "None of the method could be run, check if parameter provided" + csv_results = Channel.empty() } } else { log.info "Running all CVs with both R and Python" - CV_PYTHON ( ch_full_exp.mudata_copy ) + CV_PYTHON ( ch_full_exp.mudata_copy, unimodal_exp ) CV_R ( ch_full_exp.mae_copy ) // Mix results together and set to csv_results Channel.empty() @@ -144,9 +161,9 @@ workflow CROSS_VALIDATION { .groupTuple(by: 0) .set { csv_results } // Collect all result and mix it to merge it more - MERGE_RESULT_TABLE ( csv_results, saveMode ) + MERGE_RESULT_TABLE ( csv_results, saveMode, outcome_type ) csv_results = MERGE_RESULT_TABLE.out.csv_results } emit: csv_results -} \ No newline at end of file +} diff --git a/subworkflows/cross_validation/python/main.nf b/subworkflows/cross_validation/python/main.nf index 2516d49..6fc903a 100644 --- a/subworkflows/cross_validation/python/main.nf +++ b/subworkflows/cross_validation/python/main.nf @@ -1,11 +1,12 @@ // Methods to include -include { INTEGRAO } from "${subworkflowDir}/methods/integrao" -include { SKLEARN } from "${subworkflowDir}/methods/sklearn" -include { MOGONET } from "${subworkflowDir}/methods/mogonet" +include { INTEGRAO } from "${subworkflowDir}/methods/integrao" +include { SKLEARN } from "${subworkflowDir}/methods/sklearn" +include { MOGONET } from "${subworkflowDir}/methods/mogonet" +include { SKSURV } from "${subworkflowDir}/methods/sksurv" // This module to collect results -include { MERGE_RESULT_TABLE } from "${modulesDir}/merge_result_table" +include { MERGE_RESULT_TABLE } from "${modulesDir}/merge_result_table" // Helper fun -include { printBanner } from "${modulesDir}/functions" +include { printBanner} from "${modulesDir}/functions" // Workflow specific params to use @@ -13,11 +14,17 @@ def language_name = "Python" def saveMode = "language" workflow CV_PYTHON { - // Skip or trigger method to run - skip_integrao = params.skip_integrao // boolean: true/false - skip_sklearn = params.skip_sklearn // boolean: true/false - skip_mogonet = params.skip_mogonet // boolean: true/false - skip_goat = params.skip_goat // boolean: true/false + + // Determine if should run classification or survival one at a time only + outcome_type = params.outcome_type + + // Skip if explicitly requested OR unsupported for this outcome + skip_integrao = params.skip_integrao || outcome_type != 'classification' + skip_sklearn = params.skip_sklearn || outcome_type != 'classification' + skip_mogonet = params.skip_mogonet || outcome_type != 'classification' + skip_goat = params.skip_goat || outcome_type != 'classification' + skip_sksurv = params.skip_sksurv || outcome_type != 'survival' + // Method specific parameters he_base_dim = params.he_base_dim // Inputs of workflow @@ -25,6 +32,13 @@ workflow CV_PYTHON { mu_copy // channel of (key, key/path_to_mu, split_indices), // where each split_indices/ contains // list of txt files. + + // NEW (single-modality mode): same tuple shape as mu_copy, but containing + // the expanded per-modality datasets ('-'). Only + // sklearn consumes them: single-block input is not meaningful for + // multi-block integration methods (INTEGRAO, MOGONET). Empty channel + // when params.single_modality_mode is false. + mu_copy_unimodal main: /* Need to first allocate empty output for each of the methods @@ -42,9 +56,17 @@ workflow CV_PYTHON { // SKLEARN sklearn_results = Channel.empty() if (!skip_sklearn) { - SKLEARN ( mu_copy ) + // Mix inputs of mudata and unimodal + SKLEARN ( mu_copy.mix( mu_copy_unimodal ) ) sklearn_results = SKLEARN.out.csv_results } + + // SKSURV + sksurv_results = Channel.empty() + if (!skip_sksurv) { + SKSURV ( mu_copy ) + sksurv_results = SKSURV.out.csv_results + } // MOGONET mogonet_results = Channel.empty() @@ -57,9 +79,10 @@ workflow CV_PYTHON { // Collect all result and mix it to merge it more Channel.empty() // Then these are outputs of methods - .mix( integrao_results ) - .mix( sklearn_results ) - .mix( mogonet_results ) + .mix( integrao_results ) + .mix( sklearn_results ) + .mix( mogonet_results ) + .mix( sksurv_results ) .map { it -> [ language_name, it[0], it[1] ] // Ch [R, method name, path of summary csv of method] } @@ -68,7 +91,7 @@ workflow CV_PYTHON { .set { csv_results } // ======================================================================== // Merge result tables together - MERGE_RESULT_TABLE ( csv_results, saveMode ) + MERGE_RESULT_TABLE ( csv_results, saveMode, outcome_type ) emit: csv_results = MERGE_RESULT_TABLE.out.csv_results -} \ No newline at end of file +} diff --git a/subworkflows/cross_validation/r/main.nf b/subworkflows/cross_validation/r/main.nf index 41a20f1..8514188 100644 --- a/subworkflows/cross_validation/r/main.nf +++ b/subworkflows/cross_validation/r/main.nf @@ -14,15 +14,19 @@ params.full_mode = false def language_name = "R" def saveMode = "language" workflow CV_R { + outcome_type = params.outcome_type // "classification" or "survival" + // Skip or trigger method to run - skip_caret_multimodal = params.skip_caret_multimodal // boolean: true/false - skip_demo_logit = params.skip_demo_logit // boolean: true/false - skip_cplr = params.skip_cplr // boolean: true/false - skip_diablo = params.skip_diablo // boolean: true/false - skip_rgcca = params.skip_rgcca // boolean: true/false - skip_sgmr = params.skip_sgmr // boolean: true/false - skip_mofa = params.skip_mofa // boolean: true/false + // Classification only + skip_caret_multimodal = params.skip_caret_multimodal || outcome_type != 'classification' + skip_demo_logit = params.skip_demo_logit || outcome_type != 'classification' + skip_diablo = params.skip_diablo || outcome_type != 'classification' + skip_rgcca = params.skip_rgcca || outcome_type != 'classification' + skip_sgmr = params.skip_sgmr || outcome_type != 'classification' + // Both classification and survival + skip_mofa = params.skip_mofa + skip_cplr = params.skip_cplr // Method specific parameters num_comps = params.num_comps take: @@ -112,7 +116,7 @@ workflow CV_R { .set { csv_results } // ======================================================================== // Merge result tables together - MERGE_RESULT_TABLE ( csv_results, saveMode ) + MERGE_RESULT_TABLE ( csv_results, saveMode, outcome_type ) emit: csv_results = MERGE_RESULT_TABLE.out.csv_results } diff --git a/subworkflows/feature_selection/main.nf b/subworkflows/feature_selection/main.nf index f35799d..9ac59db 100644 --- a/subworkflows/feature_selection/main.nf +++ b/subworkflows/feature_selection/main.nf @@ -11,7 +11,7 @@ include { RGCCA_SELECT_FEATURE } from "${modulesDir}/rgcca/sele include { SKLEARN_SELECT_FEATURE } from "${modulesDir}/sklearn/select_feature" include { MOFA_SELECT_FEATURE } from "${modulesDir}/mofa/select_feature" include { MERGE_SELECTED_FEATURES } from "${modulesDir}/merge_selected_features" - +include { MERGE_SELECTED_HYPERPARAMETERS } from "${modulesDir}/merge_selected_hyperparameters" // Define vars to use later //def saveMode = "language" //def lang = "all_langs" @@ -57,58 +57,74 @@ workflow FEATURE_SELECTION { //mae_data.view { "This is mae: $it"} /* The ones in R */ caret_multimodal_features = Channel.empty() + caret_multimodal_hyperparams = Channel.empty() if (!skip_caret_multimodal) { CARET_MULTIMODAL_SELECT_FEATURE ( mae_data ) caret_multimodal_features = CARET_MULTIMODAL_SELECT_FEATURE.out.features + caret_multimodal_hyperparams = CARET_MULTIMODAL_SELECT_FEATURE.out.hyperparameters } cooperative_learning_features = Channel.empty() + cooperative_learning_hyperparams = Channel.empty() if (!skip_cplr) { COOPERATIVE_LEARNING_SELECT_FEATURE (mae_data) cooperative_learning_features = COOPERATIVE_LEARNING_SELECT_FEATURE.out.features + cooperative_learning_hyperparams = COOPERATIVE_LEARNING_SELECT_FEATURE.out.hyperparameters } diablo_features = Channel.empty() + diablo_hyperparams = Channel.empty() if (!skip_diablo) { // Connection for its design matrix ch_design = Channel.fromList( diablo_design_connection ) DIABLO_SELECT_FEATURE ( mae_data, num_comps, ch_design) diablo_features = DIABLO_SELECT_FEATURE.out.features + diablo_hyperparams = DIABLO_SELECT_FEATURE.out.hyperparameters } integrao_features = Channel.empty() + integrao_hyperparams = Channel.empty() if (!skip_integrao) { INTEGRAO_SELECT_FEATURE ( mu_data ) integrao_features = INTEGRAO_SELECT_FEATURE.out.features + integrao_hyperparams = INTEGRAO_SELECT_FEATURE.out.hyperparameters } mofa_features = Channel.empty() + mofa_hyperparams = Channel.empty() if (!skip_mofa) { MOFA_SELECT_FEATURE ( mae_data, num_comps ) mofa_features = MOFA_SELECT_FEATURE.out.features + mofa_hyperparams = MOFA_SELECT_FEATURE.out.hyperparameters } rgcca_features = Channel.empty() + rgcca_hyperparams = Channel.empty() if (!skip_rgcca) { // RGCCA can use same design matrices like full or null as if in DIABLO ch_design = Channel.fromList ( diablo_design_connection ) RGCCA_SELECT_FEATURE (mae_data, num_comps, ch_design) rgcca_features = RGCCA_SELECT_FEATURE.out.features + rgcca_hyperparams = RGCCA_SELECT_FEATURE.out.hyperparameters } /* The ones in Python */ mogonet_features = Channel.empty() + mogonet_hyperparams = Channel.empty() if (!skip_mogonet) { MOGONET_SELECT_FEATURE ( mu_data, he_base_dim) mogonet_features = MOGONET_SELECT_FEATURE.out.features + mogonet_hyperparams = MOGONET_SELECT_FEATURE.out.hyperparameters } sklearn_features = Channel.empty() + sklearn_hyperparams = Channel.empty() if (!skip_sklearn) { // Classifier from sklearn ch_sk_classifiers = Channel.fromList( sklearn_classifier_names ) SKLEARN_SELECT_FEATURE (mu_data, ch_sk_classifiers) sklearn_features = SKLEARN_SELECT_FEATURE.out.features + sklearn_hyperparams = SKLEARN_SELECT_FEATURE.out.hyperparameters } @@ -129,7 +145,20 @@ workflow FEATURE_SELECTION { //.groupTuple(by: 0) //.map { lang, methods, list_csvs -> [lang, list_csvs] } // Ch [R, list of summary table only] .set { features_csv } + /* Similarly merge the hyperparam json as csv */ + Channel.empty() + .mix( caret_multimodal_hyperparams ) + .mix( cooperative_learning_hyperparams ) + .mix( diablo_hyperparams ) + .mix( integrao_hyperparams ) + .mix( mogonet_hyperparams ) + .mix( mofa_hyperparams ) + .mix( rgcca_hyperparams ) + .mix( sklearn_hyperparams ) + .collect() + .set { hyperparams_json } // ======================================================================== // Merge result tables together MERGE_SELECTED_FEATURES ( features_csv ) + MERGE_SELECTED_HYPERPARAMETERS ( hyperparams_json ) } diff --git a/subworkflows/methods/caret_multimodal/main.nf b/subworkflows/methods/caret_multimodal/main.nf index 3ea8f42..ee7e840 100755 --- a/subworkflows/methods/caret_multimodal/main.nf +++ b/subworkflows/methods/caret_multimodal/main.nf @@ -111,7 +111,7 @@ workflow CARET_MULTIMODAL { } .set { result_table } // Lastly merge it, this would be quite fast - MERGE_RESULT_TABLE ( result_table, saveMode ) + MERGE_RESULT_TABLE ( result_table, saveMode, params.outcome_type ) // ===================================================================== // And emit the result back to upstream (which is another merge of different method) diff --git a/subworkflows/methods/cooperative_learning/main.nf b/subworkflows/methods/cooperative_learning/main.nf index 9e7691d..fa6e6d7 100644 --- a/subworkflows/methods/cooperative_learning/main.nf +++ b/subworkflows/methods/cooperative_learning/main.nf @@ -14,6 +14,8 @@ def method_name = "cooperative_learning" def saveMode = "method" workflow COOPERATIVE_LEARNING { + + outcome_type = params.outcome_type take: mae_copy // ch of tuple dataset, path of mae data, directories of fold, containing all txts main: @@ -28,7 +30,7 @@ workflow COOPERATIVE_LEARNING { TODO: Need to turn output of this to 'train_input' */ - COOPERATIVE_LEARNING_PREPROCESS ( mae_copy ) + COOPERATIVE_LEARNING_PREPROCESS ( mae_copy, params.outcome_type ) mae_copy.join(COOPERATIVE_LEARNING_PREPROCESS.out.fold_splits, by: 0) .multiMap { it -> input_data: [ it[0], it[1] ] // [dataset_name, mae_data] @@ -44,7 +46,7 @@ workflow COOPERATIVE_LEARNING { 2. Training for each fold created through previous preprocessed data */ - COOPERATIVE_LEARNING_TRAIN ( train_input ) + COOPERATIVE_LEARNING_TRAIN ( train_input, outcome_type ) // Do some transformation to make a multiMap that has two branches for predict COOPERATIVE_LEARNING_TRAIN.out.model .join(COOPERATIVE_LEARNING_TRAIN.out.test_data, by: [0, 1]) @@ -59,7 +61,8 @@ workflow COOPERATIVE_LEARNING { COOPERATIVE_LEARNING_PREDICT ( predict_input.model, predict_input.test_data, - Channel.value(method_name) + Channel.value(method_name), + outcome_type ) /* 4. Collect results of predicted folds for each data and group by GSE dataset name @@ -74,7 +77,7 @@ workflow COOPERATIVE_LEARNING { .set { result_tables } // Run these by batch of result tables (K tables per data) // result_tables.view() - MERGE_RESULT_TABLE ( result_tables, saveMode ) + MERGE_RESULT_TABLE ( result_tables, saveMode, outcome_type ) emit: csv_results = MERGE_RESULT_TABLE.out.csv_results } diff --git a/subworkflows/methods/diablo/main.nf b/subworkflows/methods/diablo/main.nf index 255c2d7..df20717 100644 --- a/subworkflows/methods/diablo/main.nf +++ b/subworkflows/methods/diablo/main.nf @@ -81,8 +81,8 @@ workflow DIABLO { // Run these by batch of result tables (K tables per data) //result_tables.view() - MERGE_RESULT_TABLE ( result_tables, saveMode ) + MERGE_RESULT_TABLE ( result_tables, saveMode, params.outcome_type ) // After ran the training part, we should validate it emit: csv_results = MERGE_RESULT_TABLE.out.csv_results -} \ No newline at end of file +} diff --git a/subworkflows/methods/integrao/main.nf b/subworkflows/methods/integrao/main.nf index 32e9b7c..31f60fc 100755 --- a/subworkflows/methods/integrao/main.nf +++ b/subworkflows/methods/integrao/main.nf @@ -104,7 +104,7 @@ workflow INTEGRAO { .set { result_table } // Lastly merge it, this would be quite fast - MERGE_RESULT_TABLE ( result_table, saveMode ) + MERGE_RESULT_TABLE ( result_table, saveMode, params.outcome_type ) // ===================================================================== // And emit the result back to upstream (which is another merge of different method) diff --git a/subworkflows/methods/mofa/main.nf b/subworkflows/methods/mofa/main.nf index faeb340..0699ec4 100644 --- a/subworkflows/methods/mofa/main.nf +++ b/subworkflows/methods/mofa/main.nf @@ -9,6 +9,7 @@ include { MERGE_RESULT_TABLE } from "${modulesDir}/merge_result_table" def saveMode = "method" workflow MOFA { runInnerCV = params.runInnerCV + outcome_type = params.outcome_type take: mae_copy // ch of tuple dataset, path of mae data, directories of fold, containing all txts num_factors // Number of factors to tune in mofa model @@ -24,7 +25,7 @@ workflow MOFA { // // MOFA_DOWNSTREAM ( mae_copy ) // } - MOFA_PREPROCESS ( mae_copy, num_factors ) + MOFA_PREPROCESS ( mae_copy, num_factors, outcome_type ) /* ========================================================================== @@ -52,7 +53,8 @@ workflow MOFA { // As mentioned, there's a possible run inner cv option // runInnerCV: boolean, true or false - MOFA_TRAIN( train_input, runInnerCV ) + // outcome_type: string, "classification" or "survival" + MOFA_TRAIN( train_input, runInnerCV, outcome_type ) // Transform certain outputs here to use in prediction // Join outputs from trained models and prepare for prediction MOFA_TRAIN.out.model @@ -65,7 +67,8 @@ workflow MOFA { // Predict each fold with their corresponding model and test data within fold MOFA_PREDICT ( predict_input.model, - predict_input.test_data + predict_input.test_data, + outcome_type ) // Collect results of predicted folds for each data and group by GSE dataset name @@ -79,8 +82,8 @@ workflow MOFA { // Run these by batch of result tables (K tables per data) //result_tables.view() - MERGE_RESULT_TABLE ( result_tables, saveMode ) + MERGE_RESULT_TABLE ( result_tables, saveMode, params.outcome_type ) // After ran the training part, we should validate it emit: csv_results = MERGE_RESULT_TABLE.out.csv_results -} \ No newline at end of file +} diff --git a/subworkflows/methods/mogonet/main.nf b/subworkflows/methods/mogonet/main.nf index d2da7be..95f156b 100644 --- a/subworkflows/methods/mogonet/main.nf +++ b/subworkflows/methods/mogonet/main.nf @@ -79,7 +79,7 @@ workflow MOGONET { .set { result_tables } // Run these by batch of result tables (K tables per data) // result_tables.view() - MERGE_RESULT_TABLE ( result_tables, saveMode ) + MERGE_RESULT_TABLE ( result_tables, saveMode, params.outcome_type ) emit: csv_results = MERGE_RESULT_TABLE.out.csv_results // Uncomment this bit below and comment above to debug @@ -92,4 +92,4 @@ Unused code // .flatMap { sample, list -> list.collate(3).collect { [ sample, it ]}}.view() -*/ \ No newline at end of file +*/ diff --git a/subworkflows/methods/rgcca/main.nf b/subworkflows/methods/rgcca/main.nf index 861bd8f..fddf753 100644 --- a/subworkflows/methods/rgcca/main.nf +++ b/subworkflows/methods/rgcca/main.nf @@ -108,7 +108,7 @@ workflow RGCCA { } .set { result_table } // Lastly merge it, this would be quite fast - MERGE_RESULT_TABLE ( result_table, saveMode ) + MERGE_RESULT_TABLE ( result_table, saveMode, params.outcome_type ) // // ===================================================================== // } diff --git a/subworkflows/methods/sklearn/main.nf b/subworkflows/methods/sklearn/main.nf index 2bc2b3a..0f1a4c3 100755 --- a/subworkflows/methods/sklearn/main.nf +++ b/subworkflows/methods/sklearn/main.nf @@ -31,9 +31,16 @@ def saveMode = "method" workflow SKLEARN { // Classifier to train for sklearn - model_name = Channel.fromList(params.sklearn_classifier_names) + // Set classifier to Logit only if single_mode is true + model_name = Channel.fromList(params.sklearn_classifier_names) + // Reduction method (pca or empty) - reduction = Channel.fromList(params.sklearn_reduction) + // Set reduction to empty if single_mode is true + if ( params.single_modality_mode ) { + reduction = Channel.of("empty") + } else { + reduction = Channel.fromList(params.sklearn_reduction) + } take: // TODO: rename this data_copy to mae_copy or mu_copy depending on language data_copy // ch of tuple dataset, path of mae/mu data, @@ -114,7 +121,7 @@ workflow SKLEARN { } .set { result_table } // Lastly merge it, this would be quite fast - MERGE_RESULT_TABLE ( result_table, saveMode ) + MERGE_RESULT_TABLE ( result_table, saveMode, params.outcome_type ) // ===================================================================== // And emit the result back to upstream (which is another merge of different method) diff --git a/subworkflows/methods/sksurv/main.nf b/subworkflows/methods/sksurv/main.nf new file mode 100644 index 0000000..5af169e --- /dev/null +++ b/subworkflows/methods/sksurv/main.nf @@ -0,0 +1,82 @@ +/* =========================================================================== */ +/* + +Survival method workflow (scikit-survival), mirrors the sklearn classification +method workflow: preprocess -> train (per fold x model) -> predict -> merge. + +Author: Tony Liang +Date: 2026-09-18 +*/ +/* =========================================================================== */ + +def method_dir = "${modulesDir}/sksurv" +include { SKSURV_PREPROCESS } from "${method_dir}/preprocess" +include { SKSURV_TRAIN } from "${method_dir}/train" +include { SKSURV_PREDICT } from "${method_dir}/predict" +include { MERGE_RESULT_TABLE } from "${modulesDir}/merge_result_table" + +// Workflow related params +def method_name = "sksurv" +def saveMode = "method" + + +workflow SKSURV { + // Survival models to train: coxnet | rsf | gbm + model_name = Channel.fromList(params.survival_model_names) + take: + // dataset name + path of mu data + directory of fold txts + data_copy + + main: + // ====================================================================== + // 1. Preprocess: partition full MuData into per-fold train/test h5mu + SKSURV_PREPROCESS ( data_copy ) + data_copy.join( SKSURV_PREPROCESS.out.fold_splits, by:0 ) + .multiMap { it -> + input_data: [ it[0], it[1] ] // [dataset_name, mu_data] + data_folds: [ it[0], it[3].flatten() ] // [fold1, fold2, ... , foldk] + }.set{ interm } + interm.input_data + .combine(interm.data_folds.transpose(), by: 0) + .set { train_input } + // ====================================================================== + + /* + 2. Train for each fold; inner CV (alpha selection for coxnet) happens + inside the training fold only. + */ + SKSURV_TRAIN ( train_input, model_name ) + + SKSURV_TRAIN.out.model + .join(SKSURV_TRAIN.out.test_data, by: [0, 1, 2]) + .multiMap { it -> + model: [ it[0], it[1], it[2], it[3] ] // [ dataset_name, fold_name, model_name, model ] + test_data: [ it[0], it[1], it[2], it[4] ] // [ dataset_name, fold_name, model_name, test_data] + }.set { predict_input } + // ====================================================================== + + /* + 3. Predict per fold: per-observation risk score (lp) + S(t) horizons. + */ + SKSURV_PREDICT ( + predict_input.model, + predict_input.test_data, + Channel.value(method_name) + ) + // ====================================================================== + + /* + 4. Collect results of predicted folds and merge. + */ + SKSURV_PREDICT.out.result_table + .groupTuple(by: 2) + .map {it -> + [ it[2], it[3] ] + } + .set { result_table } + MERGE_RESULT_TABLE ( result_table, saveMode, params.outcome_type ) + // ===================================================================== + + emit: + csv_results = MERGE_RESULT_TABLE.out.csv_results +} diff --git a/subworkflows/prepare_data/main.nf b/subworkflows/prepare_data/main.nf index 60c89b0..b741b33 100644 --- a/subworkflows/prepare_data/main.nf +++ b/subworkflows/prepare_data/main.nf @@ -2,6 +2,7 @@ def prepare_dir = "${modulesDir}/prepare_data" include { PREPARE_MAE_DATA } from "${prepare_dir}/prepare_mae_data" include { PREPARE_MU_DATA } from "${prepare_dir}/prepare_mu_data" +include { SPLIT_MODALITY } from "${prepare_dir}/split_modality" include { PARSE_METADATA } from "${prepare_dir}/parse_metadata" include { UNCOMPRESS_RECORD } from "${prepare_dir}/uncompress_record" // More generic utitilies here @@ -21,8 +22,10 @@ workflow PREPARE_DATA { // ch_output = Channel.fromList( // Output format MAE and MuData // params.output_formats // ) - runSimple = params.runSimple - filter_low_var = params.filter_low_var // Filter the low variance features (default: "0") + filter_low_var = params.filter_low_var // Filter the low variance features (default: "0") + single_modality_mode = params.single_modality_mode // Single modality mode (default: "False") When true, each data is expanded to all modalities, + // and force to run in sklearn pipeline to test modality baselines + outcome_type = params.outcome_type // Outcome type (default: "classification") One of 'classification' or 'survival' /* Workflow starts here */ // Workflow required input take: @@ -46,22 +49,7 @@ workflow PREPARE_DATA { ch_datasets = Channel.empty() UNCOMPRESS_RECORD ( data_records ) - // // Join by common dataset_name and make it multimap - // UNCOMPRESS_RECORD.out - // .record - // // .map { dname, mae_path, mu_path -> - // // [ - // // [ dname, mae_path, mu_path] - // // // dataset_name: it[0], - // // // mae_path: it[1], - // // // mu_path: it[2] - // // ] - // // } - // .multiMap{ dname, mae_path, mu_path -> - // mae_pt: [dname, mae_path] - // mu_pt: [dname, mu_path] - // } - // .set { real_data } + // Expand to different parts of the channel for MAE and MuData UNCOMPRESS_RECORD.out.record.map { dname, mae_path, mu_path -> [ dname, mae_path ] }.set { mae_pt } UNCOMPRESS_RECORD.out.record.map { dname, mae_path, mu_path -> [ dname, mu_path ] }.set { mu_pt } @@ -69,10 +57,10 @@ workflow PREPARE_DATA { // /* ===================================================================== */ // // Have a process to check the right format for MAE // // output MAE back with suitable transformations? - PREPARE_MAE_DATA ( mae_pt, filter_low_var ) + PREPARE_MAE_DATA ( mae_pt, filter_low_var, outcome_type ) // // Have a process to check the right format for MuData // // output MuData back with suitable transformation? - PREPARE_MU_DATA ( mu_pt, filter_low_var ) + PREPARE_MU_DATA ( mu_pt, filter_low_var, outcome_type ) // // // TODO: Need to test this bit first // PREPARE_MAE_DATA.out // .mae_data @@ -99,6 +87,26 @@ workflow PREPARE_DATA { .map { dname, data_path -> data_path } .collect() ) + + // ===================================================================== + // NEW (single-modality mode): expand each processed MuData into one + // single-modality MuData per omics, named '-'. + // These are first-class datasets: SPLITTING runs on them independently + // and yields folds identical to the parent multimodal dataset (the + // splitter depends only on y and obs row order, both preserved here). + // ===================================================================== + mu_data_unimodal = Channel.empty() + if ( single_modality_mode ) { + SPLIT_MODALITY ( PREPARE_MU_DATA.out.mu_data ) + // Derive the child dataset name from the file name: + // 'rosmap-genomics_processed.h5mu' -> 'rosmap-genomics' + SPLIT_MODALITY.out.unimodal_files + .flatMap { dname, files -> + files.collect { f -> [ f.name.replace('_processed.h5mu', ''), f ] } + } + .set { mu_data_unimodal } + } + // TODO: Do the conversions later? assume two formats exists now // ch_datasets = Channel.empty() emit: @@ -110,7 +118,8 @@ workflow PREPARE_DATA { N = n(real_data) + n(simulated_data) N = n(real_data) */ - mae_data = PREPARE_MAE_DATA.out.mae_data - mu_data = PREPARE_MU_DATA.out.mu_data + mae_data = PREPARE_MAE_DATA.out.mae_data + mu_data = PREPARE_MU_DATA.out.mu_data + mu_data_unimodal = mu_data_unimodal } diff --git a/workflows/messi_benchmark.nf b/workflows/messi_benchmark.nf index 42c29ca..dc63853 100644 --- a/workflows/messi_benchmark.nf +++ b/workflows/messi_benchmark.nf @@ -60,17 +60,35 @@ workflow MESSI_BENCHMARK { // SUBWORKFLOW: Perform splitting with stratification to the response // variable on each dataset // - SPLITTING ( PREPARE_DATA.out.mu_data, params.split_type ) + + // NEW (single-modality mode): the unimodal datasets are mixed into the + // same SPLITTING call. When params.single_modality_mode is false, + // mu_data_unimodal is an empty channel and the mix is a no-op. + // Unimodal children get folds identical to their parent multimodal + // dataset (the splitter depends only on y and obs row order). + // + + SPLITTING ( + PREPARE_DATA.out.mu_data.mix ( PREPARE_DATA.out.mu_data_unimodal ), + params.split_type + ) // // SUBWORKFLOW: Perform cross validation for each dataset using these indices // - CROSS_VALIDATION ( PREPARE_DATA.out.mae_data, PREPARE_DATA.out.mu_data, SPLITTING.out.splits_indices ) + // NEW: 4th input routes the unimodal datasets; inside CROSS_VALIDATION + // they are consumed by python methods only (currently sklearn). + // + CROSS_VALIDATION ( + PREPARE_DATA.out.mae_data, PREPARE_DATA.out.mu_data, + SPLITTING.out.splits_indices, PREPARE_DATA.out.mu_data_unimodal + ) // // MODULE: Use the output of cross validation to calculate metrics // CALCULATE_METRICS ( CROSS_VALIDATION.out.csv_results, - params.threshold + params.threshold, + params.outcome_type ) } } @@ -120,4 +138,4 @@ workflow.onError = { ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ THE END ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -*/ \ No newline at end of file +*/