From e308292a34225eee152bde8b35148f65496366b9 Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Thu, 13 Aug 2026 21:56:46 +0800 Subject: [PATCH 01/29] fix: move container option to process level config, and remove stale temp var option --- modules/caret_multimodal/predict/main.nf | 5 +---- modules/caret_multimodal/preprocess/main.nf | 7 +----- .../caret_multimodal/select_feature/main.nf | 5 +---- modules/caret_multimodal/train/main.nf | 8 ++----- modules/prepare_data/parse_metadata/main.nf | 7 ++---- modules/prepare_data/prepare_mae_data/main.nf | 6 +---- modules/prepare_data/prepare_mu_data/main.nf | 5 +---- .../prepare_data/uncompress_record/main.nf | 6 ++--- nextflow.config | 22 +++++++++++++++---- 9 files changed, 29 insertions(+), 42 deletions(-) diff --git a/modules/caret_multimodal/predict/main.nf b/modules/caret_multimodal/predict/main.nf index ffada37..e842c78 100755 --- a/modules/caret_multimodal/predict/main.nf +++ b/modules/caret_multimodal/predict/main.nf @@ -3,12 +3,9 @@ include { getPublishPath } from "${modulesDir}/functions" process CARET_MULTIMODAL_PREDICT { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${dataset_name}-${fold_name}" debug "${params.debug}" - container "${ onSockeye ? - 'caret_multimodal.sif' : - 'tonyliang19/caret_multimodal:latest' }" + label 'caret_multimodal' publishDir ( path: "${params.outdir}/${getPublishPath(task.process)}/${dataset_name}/${fold_name}", diff --git a/modules/caret_multimodal/preprocess/main.nf b/modules/caret_multimodal/preprocess/main.nf index 6ae03c6..e1f7390 100755 --- a/modules/caret_multimodal/preprocess/main.nf +++ b/modules/caret_multimodal/preprocess/main.nf @@ -21,17 +21,12 @@ include { getPublishPath } from "${modulesDir}/functions" // TODO: Replace name here process CARET_MULTIMODAL_PREPROCESS { - // temp variables to use - def onSockeye = workflow.projectDir.toString().contains('/scratch') // process level configuration debug "${params.debug}" // debugs true or false by param in MESSI.config tag "${dataset_name}" // identifier of process when ran in parallel // By var before to determine what container to use // Uses apptainer if true otherwise docker - container "${ onSockeye ? - 'caret_multimodal.sif' : - 'tonyliang19/caret_multimodal:latest' }" - + label "caret_multimodal" label "process_low" /* outputs store to outdir for saving and inspecting purpose */ diff --git a/modules/caret_multimodal/select_feature/main.nf b/modules/caret_multimodal/select_feature/main.nf index acd3eb2..85171bb 100755 --- a/modules/caret_multimodal/select_feature/main.nf +++ b/modules/caret_multimodal/select_feature/main.nf @@ -4,13 +4,10 @@ include { getPublishPath } from "${modulesDir}/functions" process CARET_MULTIMODAL_SELECT_FEATURE { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${dataset_name}" debug true label 'process_high' - container "${ onSockeye ? - 'caret_multimodal.sif' : - 'tonyliang19/caret_multimodal:latest' }" + label 'caret_multimodal' publishDir ( path: "${params.outdir}/${task.process.tokenize(':').join('/').toLowerCase()}/${dataset_name}", diff --git a/modules/caret_multimodal/train/main.nf b/modules/caret_multimodal/train/main.nf index 40fb4eb..d184360 100755 --- a/modules/caret_multimodal/train/main.nf +++ b/modules/caret_multimodal/train/main.nf @@ -17,14 +17,9 @@ include { getPublishPath } from "${modulesDir}/functions" process CARET_MULTIMODAL_TRAIN { - // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') // Identifier for each dataset and fold combination tag "${dataset_name}-${fold_path.name}" - // TODO: rename to the actual image name used - container "${ onSockeye ? - 'caret_multimodal.sif' : - 'tonyliang19/caret_multimodal:latest' }" + // Parse the output directory to migrate results to publishDir ( path: "${params.outdir}/${getPublishPath(task.process)}/${dataset_name}/${fold_path.name}", @@ -34,6 +29,7 @@ process CARET_MULTIMODAL_TRAIN { // TODO: add your custom labels to tell if process // consumes large RAM/ROM and gpu access or not label 'process_medium' + label 'caret_multimodal' // TODO: Change arg name to mae_path or mu_path diff --git a/modules/prepare_data/parse_metadata/main.nf b/modules/prepare_data/parse_metadata/main.nf index 05114a5..2222ee1 100644 --- a/modules/prepare_data/parse_metadata/main.nf +++ b/modules/prepare_data/parse_metadata/main.nf @@ -11,16 +11,13 @@ //include { getPublishPath } from "${modulesDir}/functions" // TODO: Replace name here process PARSE_METADATA { - // temp variables to use - def onSockeye = workflow.projectDir.toString().contains('/scratch') // process level configuration debug "${params.debug}" // debugs true or false by param in MESSI.config // tag "${dataset_name}" // identifier of process when ran in parallel // By var before to determine what container to use // Uses apptainer if true otherwise docker - container "${ onSockeye ? - 'save_simulate.sif' : - 'tonyliang19/save_simulate:latest'}" + label "generic" + /* outputs store to outdir for saving and inspecting purpose */ publishDir ( diff --git a/modules/prepare_data/prepare_mae_data/main.nf b/modules/prepare_data/prepare_mae_data/main.nf index b97b1cc..3e0b452 100644 --- a/modules/prepare_data/prepare_mae_data/main.nf +++ b/modules/prepare_data/prepare_mae_data/main.nf @@ -8,14 +8,10 @@ */ process PREPARE_MAE_DATA { - def onSockeye = workflow.projectDir.toString().contains('/scratch') /* process metadata and configs */ tag "${dataset_name}" - // TODO: move this containter to somewhere else? label 'process_low' - container "${ onSockeye ? - 'save_simulate.sif' : - 'tonyliang19/save_simulate:latest' }" + label 'generic' publishDir ( path: "${params.outdir}/${task.process.tokenize(':').join('/').toLowerCase()}/${dataset_name}", saveAS: { file }, diff --git a/modules/prepare_data/prepare_mu_data/main.nf b/modules/prepare_data/prepare_mu_data/main.nf index 8b0ebff..47a1c9f 100644 --- a/modules/prepare_data/prepare_mu_data/main.nf +++ b/modules/prepare_data/prepare_mu_data/main.nf @@ -10,10 +10,7 @@ process PREPARE_MU_DATA { /* process metadata and configs */ tag "${dataset_name}" label 'process_low' - // TODO: move this containter to somewhere else? - container "${ onSockeye ? - 'save_simulate.sif' : - 'tonyliang19/save_simulate:latest' }" + label 'generic' publishDir ( path: "${params.outdir}/${task.process.tokenize(':').join('/').toLowerCase()}/${dataset_name}", saveAS: { file }, diff --git a/modules/prepare_data/uncompress_record/main.nf b/modules/prepare_data/uncompress_record/main.nf index de3fbd2..dbb2c43 100644 --- a/modules/prepare_data/uncompress_record/main.nf +++ b/modules/prepare_data/uncompress_record/main.nf @@ -3,10 +3,8 @@ process UNCOMPRESS_RECORD { /* process metadata and configs */ tag "${row_map.dataset_name}" label 'process_low' - // TODO: move this containter to somewhere else? - container "${ onSockeye ? - 'save_simulate.sif' : - 'tonyliang19/save_simulate:latest' }" + label 'generic' + publishDir ( path: "${params.outdir}/${task.process.tokenize(':').join('/').toLowerCase()}/${row_map.dataset_name}", saveAS: { file }, diff --git a/nextflow.config b/nextflow.config index bc88ede..3228614 100644 --- a/nextflow.config +++ b/nextflow.config @@ -136,6 +136,9 @@ params { cpu_training_array_size = 1 // Training array size cpu_prediction_array_size = 1 // Prediction array size feature_selection_array_size = 1 // Feature selection array size + + + } // Widely-used environment variables @@ -162,7 +165,7 @@ profiles { docker.enabled = true docker.temp = 'auto' // This is temporal fix - docker.runOptions = "-u \$(id -u):\$(id -g) -v ${projectDir}:/home/rstudio/${projectDir} --entrypoint ''" + docker.runOptions = "-u \$(id -u):\$(id -g) --entrypoint '' -v ${projectDir}:/home/rstudio/${projectDir}" apptainer.enabled = false singularity.enabled = false } @@ -224,9 +227,20 @@ profiles { // Process option cache and their shell option process { - cache = 'lenient' - // Capture exit codes from upstream processes when piping - shell = ['/bin/bash', '-euo', 'pipefail'] + cache = 'lenient' + // Capture exit codes from upstream processes when piping + shell = ['/bin/bash', '-euo', 'pipefail'] + + /* + MODULE CONTAINER CONFIG + */ + + // Assigning module containers here + + withLabel: caret_multimodal { container = 'tonyliang19/caret_multimodal:latest' } + withLabel: generic { container = 'tonyliang19/save_simulate:latest' } + withLabel: mogonet { container = 'tonyliang19/mogonet:latest' } + } // Pipeline information after exectuion From b0c6ba5550e7ceab94a3f63201704f2d01a64ebd Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Tue, 18 Aug 2026 00:21:53 +0800 Subject: [PATCH 02/29] fix: moved all container definition into config --- modules/calculate_metrics/main.nf | 8 ++------ modules/cooperative_learning/predict/main.nf | 5 +---- modules/cooperative_learning/preprocess/main.nf | 5 +---- modules/cooperative_learning/select_feature/main.nf | 7 ++----- modules/cooperative_learning/train/main.nf | 5 +---- modules/diablo/downstream/main.nf | 5 +---- modules/diablo/predict/main.nf | 5 +---- modules/diablo/preprocess/main.nf | 5 +---- modules/diablo/select_feature/main.nf | 5 +---- modules/diablo/train/main.nf | 5 +---- modules/integrao/predict/main.nf | 5 +---- modules/integrao/preprocess/main.nf | 5 +---- modules/integrao/select_feature/main.nf | 5 +---- modules/integrao/train/main.nf | 5 +---- modules/local/samplesheet_check/main.nf | 5 +---- modules/merge_result_table/main.nf | 5 +---- modules/merge_selected_features/main.nf | 5 +---- modules/mofa/predict/main.nf | 5 +---- modules/mofa/preprocess/main.nf | 5 +---- modules/mofa/select_feature/main.nf | 5 +---- modules/mofa/train/main.nf | 5 +---- modules/mogonet/predict/main.nf | 7 ++----- modules/mogonet/preprocess/main.nf | 5 +---- modules/mogonet/select_feature/main.nf | 5 +---- modules/mogonet/train/main.nf | 5 +---- modules/rgcca/predict/main.nf | 9 ++------- modules/rgcca/preprocess/main.nf | 5 +---- modules/rgcca/select_feature/main.nf | 5 +---- modules/rgcca/train/main.nf | 5 +---- modules/sklearn/predict/main.nf | 5 +---- modules/sklearn/preprocess/main.nf | 5 +---- modules/sklearn/select_feature/main.nf | 5 +---- modules/sklearn/train/main.nf | 5 +---- modules/split_train_test/main.nf | 7 ++----- nextflow.config | 7 ++++++- 35 files changed, 45 insertions(+), 145 deletions(-) diff --git a/modules/calculate_metrics/main.nf b/modules/calculate_metrics/main.nf index fb2575e..dfc71c6 100644 --- a/modules/calculate_metrics/main.nf +++ b/modules/calculate_metrics/main.nf @@ -1,11 +1,7 @@ process CALCULATE_METRICS { - def onSockeye = workflow.projectDir.toString().contains('/scratch') debug true label 'process_single' - container "${ onSockeye ? - 'mogonet.sif' : - 'tonyliang19/mogonet:latest' }" - + label 'mogonet' // Should be a generic python container? publishDir ( path: "${params.outdir}/${task.process.tokenize(':').join('/').toLowerCase()}", mode: 'copy', @@ -28,4 +24,4 @@ process CALCULATE_METRICS { --threshold=${threshold} > \ ${task.process.tokenize(':')[-1].toLowerCase()}.log """ -} \ No newline at end of file +} diff --git a/modules/cooperative_learning/predict/main.nf b/modules/cooperative_learning/predict/main.nf index 674d0f6..dd8724c 100644 --- a/modules/cooperative_learning/predict/main.nf +++ b/modules/cooperative_learning/predict/main.nf @@ -2,13 +2,10 @@ include { getPublishPath } from "${modulesDir}/functions" process COOPERATIVE_LEARNING_PREDICT { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${dataset_name}-${fold_name}" debug true label 'process_single' - container "${ onSockeye ? - 'codia.sif' : - 'tonyliang19/codia:latest' }" + label 'codia' publishDir ( path: "${params.outdir}/${getPublishPath(task.process)}/${dataset_name}/${fold_name}", diff --git a/modules/cooperative_learning/preprocess/main.nf b/modules/cooperative_learning/preprocess/main.nf index d146759..b405711 100644 --- a/modules/cooperative_learning/preprocess/main.nf +++ b/modules/cooperative_learning/preprocess/main.nf @@ -5,14 +5,11 @@ include { getPublishPath } from "${modulesDir}/functions" process COOPERATIVE_LEARNING_PREPROCESS { // Temp variables - def onSockeye = workflow.projectDir.toString().contains('/scratch') /* Directives for process */ debug "${params.debug}" // default is true tag "${dataset_name}" label 'process_single' - container "${ onSockeye ? - 'codia.sif' : - 'tonyliang19/codia:latest'}" + label 'codia' // More special publish dir by getting method/method_abc to method/abc publishDir ( diff --git a/modules/cooperative_learning/select_feature/main.nf b/modules/cooperative_learning/select_feature/main.nf index d7836f2..cddec60 100644 --- a/modules/cooperative_learning/select_feature/main.nf +++ b/modules/cooperative_learning/select_feature/main.nf @@ -4,13 +4,10 @@ include { getPublishPath } from "${modulesDir}/functions" process COOPERATIVE_LEARNING_SELECT_FEATURE { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${dataset_name}" debug true label 'process_high' - container "${ onSockeye ? - 'codia.sif' : - 'tonyliang19/codia:latest' }" + label 'codia' publishDir ( path: "${params.outdir}/${task.process.tokenize(':').join('/').toLowerCase()}/${dataset_name}", @@ -21,7 +18,7 @@ process COOPERATIVE_LEARNING_SELECT_FEATURE { // Labels label 'low_mem' label 'cpu' - label 'codia' + label 'codia' /* Input and output blocks*/ diff --git a/modules/cooperative_learning/train/main.nf b/modules/cooperative_learning/train/main.nf index bfccb08..6fc9dff 100644 --- a/modules/cooperative_learning/train/main.nf +++ b/modules/cooperative_learning/train/main.nf @@ -7,13 +7,10 @@ include { getPublishPath } from "${modulesDir}/functions" process COOPERATIVE_LEARNING_TRAIN { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${dataset_name}-${fold_path.name}" debug "${params.debug}" label 'process_high' - container "${ onSockeye ? - 'codia.sif' : - 'tonyliang19/codia:latest' }" + label 'codia' // More special publish dir by getting method/method_abc to method/abc publishDir ( path: "${params.outdir}/${getPublishPath(task.process)}/${dataset_name}/${fold_path.name}", diff --git a/modules/diablo/downstream/main.nf b/modules/diablo/downstream/main.nf index 8600457..ce7ddea 100644 --- a/modules/diablo/downstream/main.nf +++ b/modules/diablo/downstream/main.nf @@ -22,16 +22,13 @@ include { getPublishPath } from "${modulesDir}/functions" // TODO: Replace name here process DIABLO_DOWNSTREAM { // temp variables to use - def onSockeye = workflow.projectDir.toString().contains('/scratch') // process level configuration debug "${params.debug}" // debugs true or false by param in MESSI.config label 'process_single' tag "${dataset_name}" // identifier of process when ran in parallel // By var before to determine what container to use // Uses apptainer if true otherwise docker - container "${ onSockeye ? - 'codia.sif' : - 'tonyliang19/codia:latest'}" + label 'codia' /* outputs store to outdir for saving and inspecting purpose */ publishDir ( diff --git a/modules/diablo/predict/main.nf b/modules/diablo/predict/main.nf index 5e563e7..b11a39d 100644 --- a/modules/diablo/predict/main.nf +++ b/modules/diablo/predict/main.nf @@ -3,13 +3,10 @@ include { getPublishPath } from "${modulesDir}/functions" process DIABLO_PREDICT { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${dataset_name}-${fold_name}-${method}" debug true label 'process_single' - container "${ onSockeye ? - 'codia.sif' : - 'tonyliang19/codia:latest' }" + label 'codia' publishDir ( //path: "${params.outdir}/${task.process.tokenize(':')[-1].toLowerCase()}/${dataset_name}/${fold_name}", diff --git a/modules/diablo/preprocess/main.nf b/modules/diablo/preprocess/main.nf index b053580..15b3452 100644 --- a/modules/diablo/preprocess/main.nf +++ b/modules/diablo/preprocess/main.nf @@ -27,16 +27,13 @@ include { getPublishPath } from "${modulesDir}/functions" process DIABLO_PREPROCESS { // temp variables to use - def onSockeye = workflow.projectDir.toString().contains('/scratch') // process level configuration debug "${params.debug}" // debugs true or false by param in MESSI.config label 'process_single' tag "${dataset_name}" // identifier of process when ran in parallel // By var before to determine what container to use // Uses apptainer if true otherwise docker - container "${ onSockeye ? - 'codia.sif' : - 'tonyliang19/codia:latest'}" + label 'codia' /* outputs store to outdir for saving and inspecting purpose */ publishDir ( diff --git a/modules/diablo/select_feature/main.nf b/modules/diablo/select_feature/main.nf index a465a8d..099a734 100644 --- a/modules/diablo/select_feature/main.nf +++ b/modules/diablo/select_feature/main.nf @@ -4,13 +4,10 @@ include { getPublishPath } from "${modulesDir}/functions" process DIABLO_SELECT_FEATURE { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${dataset_name}-design_${design}-ncomp_${ncomp}" debug true label 'process_high' - container "${ onSockeye ? - 'codia.sif' : - 'tonyliang19/codia:latest' }" + label 'codia' publishDir ( path: "${params.outdir}/${task.process.tokenize(':').join('/').toLowerCase()}/${dataset_name}", diff --git a/modules/diablo/train/main.nf b/modules/diablo/train/main.nf index 49f5970..e42c039 100644 --- a/modules/diablo/train/main.nf +++ b/modules/diablo/train/main.nf @@ -6,15 +6,12 @@ include { getPublishPath } from "${modulesDir}/functions" process DIABLO_TRAIN { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') def method_name = "diablo" tag "${dataset_name}-${fold_path.name}-design_${design}-ncomp_${ncomp}" //debug "${params.debug}" debug true label 'process_medium' - container "${ onSockeye ? - 'codia.sif' : - 'tonyliang19/codia:latest' }" + label 'codia' // Parse this path publishDir ( //path: "${params.outdir}/${task.process.tokenize(':')[-1].toLowerCase()}/${dataset_name}/${fold_path.name}", diff --git a/modules/integrao/predict/main.nf b/modules/integrao/predict/main.nf index e3ebaeb..c281125 100755 --- a/modules/integrao/predict/main.nf +++ b/modules/integrao/predict/main.nf @@ -3,12 +3,9 @@ include { getPublishPath } from "${modulesDir}/functions" process INTEGRAO_PREDICT { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${dataset_name}-${fold_name}" debug "${params.debug}" - container "${ onSockeye ? - 'integrao.sif' : - 'tonyliang19/integrao:latest' }" + label 'integrao' publishDir ( path: "${params.outdir}/${getPublishPath(task.process)}/${dataset_name}/${fold_name}", diff --git a/modules/integrao/preprocess/main.nf b/modules/integrao/preprocess/main.nf index a5d4224..671b558 100755 --- a/modules/integrao/preprocess/main.nf +++ b/modules/integrao/preprocess/main.nf @@ -22,15 +22,12 @@ include { getPublishPath } from "${modulesDir}/functions" // TODO: Replace name here process INTEGRAO_PREPROCESS { // temp variables to use - def onSockeye = workflow.projectDir.toString().contains('/scratch') // process level configuration debug "${params.debug}" // debugs true or false by param in MESSI.config tag "${dataset_name}" // identifier of process when ran in parallel // By var before to determine what container to use // Uses apptainer if true otherwise docker - container "${ onSockeye ? - 'integrao.sif' : - 'tonyliang19/integrao:latest' }" + label 'integrao' label "process_low" diff --git a/modules/integrao/select_feature/main.nf b/modules/integrao/select_feature/main.nf index 8bf3316..734bb7a 100755 --- a/modules/integrao/select_feature/main.nf +++ b/modules/integrao/select_feature/main.nf @@ -4,14 +4,11 @@ include { getPublishPath } from "${modulesDir}/functions" process INTEGRAO_SELECT_FEATURE { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${dataset_name}" debug true label 'process_high' label 'gpu' - container "${ onSockeye ? - 'integrao.sif' : - 'tonyliang19/integrao:latest' }" + label 'integrao' publishDir ( path: "${params.outdir}/${task.process.tokenize(':').join('/').toLowerCase()}/${dataset_name}", diff --git a/modules/integrao/train/main.nf b/modules/integrao/train/main.nf index 95ea531..f8c2110 100755 --- a/modules/integrao/train/main.nf +++ b/modules/integrao/train/main.nf @@ -18,13 +18,10 @@ include { getPublishPath } from "${modulesDir}/functions" process INTEGRAO_TRAIN { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') // Identifier for each dataset and fold combination tag "${dataset_name}-${fold_path.name}" // TODO: rename to the actual image name used - container "${ onSockeye ? - 'integrao.sif' : - 'tonyliang19/integrao:latest' }" + label 'integrao' // Parse the output directory to migrate results to publishDir ( path: "${params.outdir}/${getPublishPath(task.process)}/${dataset_name}/${fold_path.name}", diff --git a/modules/local/samplesheet_check/main.nf b/modules/local/samplesheet_check/main.nf index c4739b0..4e04f0f 100644 --- a/modules/local/samplesheet_check/main.nf +++ b/modules/local/samplesheet_check/main.nf @@ -2,10 +2,7 @@ process SAMPLESHEET_CHECK { tag "$samplesheet" label 'process_single_low' - def onSockeye = workflow.projectDir.toString().contains('/scratch') - container "${ onSockeye ? - 'mogonet.sif' : - 'tonyliang19/mogonet:latest' }" + label 'mogonet' input: path(samplesheet) diff --git a/modules/merge_result_table/main.nf b/modules/merge_result_table/main.nf index 0227cb8..6794d42 100644 --- a/modules/merge_result_table/main.nf +++ b/modules/merge_result_table/main.nf @@ -1,11 +1,8 @@ process MERGE_RESULT_TABLE { - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${method_name}" debug true label 'process_single' - container "${ onSockeye ? - 'codia.sif' : - 'tonyliang19/codia:latest' }" + label 'codia' publishDir ( //path: "${params.outdir}/${task.process.tokenize(':')[-1].toLowerCase()}/${method_name}", diff --git a/modules/merge_selected_features/main.nf b/modules/merge_selected_features/main.nf index 362fd45..86904fa 100644 --- a/modules/merge_selected_features/main.nf +++ b/modules/merge_selected_features/main.nf @@ -1,11 +1,8 @@ process MERGE_SELECTED_FEATURES { - def onSockeye = workflow.projectDir.toString().contains('/scratch') //tag "${method_name}" debug true label 'process_single' - container "${ onSockeye ? - 'codia.sif' : - 'tonyliang19/codia:latest' }" + label 'codia' publishDir ( //path: "${params.outdir}/${task.process.tokenize(':')[-1].toLowerCase()}/${method_name}", diff --git a/modules/mofa/predict/main.nf b/modules/mofa/predict/main.nf index 0c31e3a..b1d1b26 100644 --- a/modules/mofa/predict/main.nf +++ b/modules/mofa/predict/main.nf @@ -3,13 +3,10 @@ include { getPublishPath } from "${modulesDir}/functions" process MOFA_PREDICT { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${dataset_name}-${fold_name}" debug true label 'process_single' - container "${ onSockeye ? - 'mofa.sif' : - 'tonyliang19/mofa:latest' }" + label 'mofa' publishDir ( path: "${params.outdir}/${getPublishPath(task.process)}/${dataset_name}/${fold_name}", diff --git a/modules/mofa/preprocess/main.nf b/modules/mofa/preprocess/main.nf index 15c4b80..ff7efbe 100644 --- a/modules/mofa/preprocess/main.nf +++ b/modules/mofa/preprocess/main.nf @@ -27,16 +27,13 @@ include { getPublishPath } from "${modulesDir}/functions" process MOFA_PREPROCESS { // temp variables to use - def onSockeye = workflow.projectDir.toString().contains('/scratch') // process level configuration debug "${params.debug}" // debugs true or false by param in MESSI.config label 'process_single' tag "${dataset_name}-num_factor_${num_factors}" // identifier of process when ran in parallel // By var before to determine what container to use // Uses apptainer if true otherwise docker - container "${ onSockeye ? - 'mofa.sif' : - 'tonyliang19/mofa:latest'}" + label 'mofa' /* outputs store to outdir for saving and inspecting purpose */ publishDir ( diff --git a/modules/mofa/select_feature/main.nf b/modules/mofa/select_feature/main.nf index 41054e5..699bab0 100644 --- a/modules/mofa/select_feature/main.nf +++ b/modules/mofa/select_feature/main.nf @@ -4,13 +4,10 @@ include { getPublishPath } from "${modulesDir}/functions" process MOFA_SELECT_FEATURE { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${dataset_name}" debug true label 'process_high' - container "${ onSockeye ? - 'mofa.sif' : - 'tonyliang19/mofa:latest' }" + label 'mofa' publishDir ( path: "${params.outdir}/${task.process.tokenize(':').join('/').toLowerCase()}/${dataset_name}", diff --git a/modules/mofa/train/main.nf b/modules/mofa/train/main.nf index fbba820..603ddfd 100644 --- a/modules/mofa/train/main.nf +++ b/modules/mofa/train/main.nf @@ -6,14 +6,11 @@ include { getPublishPath } from "${modulesDir}/functions" process MOFA_TRAIN { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${dataset_name}-${fold_path.name}" //debug "${params.debug}" debug true label 'process_medium' - container "${ onSockeye ? - 'mofa.sif' : - 'tonyliang19/mofa:latest' }" + label 'mofa' // Parse this path publishDir ( path: "${params.outdir}/${getPublishPath(task.process)}/${dataset_name}/${fold_path.name}", diff --git a/modules/mogonet/predict/main.nf b/modules/mogonet/predict/main.nf index 6578977..63821f6 100644 --- a/modules/mogonet/predict/main.nf +++ b/modules/mogonet/predict/main.nf @@ -2,13 +2,10 @@ include { getPublishPath } from "${modulesDir}/functions" process MOGONET_PREDICT { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${dataset_name}-${fold_name}" debug true label 'process_low' - container "${ onSockeye ? - 'mogonet.sif' : - 'tonyliang19/mogonet:latest' }" + label 'mogonet' publishDir ( path: "${params.outdir}/${getPublishPath(task.process)}/${dataset_name}/${fold_name}", @@ -19,7 +16,7 @@ process MOGONET_PREDICT { // Labels label 'low_mem' label 'gpu' - label 'mogonet' + label 'mogonet' input: tuple val(dataset_name), val(fold_name), path(model) diff --git a/modules/mogonet/preprocess/main.nf b/modules/mogonet/preprocess/main.nf index 571212a..8fa0e12 100644 --- a/modules/mogonet/preprocess/main.nf +++ b/modules/mogonet/preprocess/main.nf @@ -31,16 +31,13 @@ include { getPublishPath } from "${modulesDir}/functions" process MOGONET_PREPROCESS { // Temp variables - def onSockeye = workflow.projectDir.toString().contains('/scratch') def output_dir = "fold" // process configurations label 'process_low' debug "${params.debug}" tag "${dataset_name}" // If on sockeye, then use sif file otherwise assuming local (use docker) - container "${ onSockeye ? - 'mogonet.sif' : - 'tonyliang19/save_simulate:latest' }" + label 'mogonet' /* Outputs are stored to outdir (fold), for saving purpose, passing input/output should be done with channels */ publishDir ( diff --git a/modules/mogonet/select_feature/main.nf b/modules/mogonet/select_feature/main.nf index 5833224..3c16934 100644 --- a/modules/mogonet/select_feature/main.nf +++ b/modules/mogonet/select_feature/main.nf @@ -4,13 +4,10 @@ include { getPublishPath } from "${modulesDir}/functions" process MOGONET_SELECT_FEATURE { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${dataset_name}-he${he_base_dim}" debug true label 'process_high' - container "${ onSockeye ? - 'mogonet.sif' : - 'tonyliang19/mogonet:latest' }" + label 'mogonet' publishDir ( path: "${params.outdir}/${task.process.tokenize(':').join('/').toLowerCase()}/${dataset_name}", diff --git a/modules/mogonet/train/main.nf b/modules/mogonet/train/main.nf index 4a8c842..c4d746a 100644 --- a/modules/mogonet/train/main.nf +++ b/modules/mogonet/train/main.nf @@ -11,16 +11,13 @@ include { getPublishPath } from "${modulesDir}/functions" process MOGONET_TRAIN { // This label is defined in nextflow.config // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') // Process configurations tag "${dataset_name}-${fold_path.name}-he${he_base_dim}" label 'python' label 'process_medium' label 'gpu' debug "${params.debug}" - container "${ onSockeye ? - 'mogonet.sif' : - 'tonyliang19/mogonet:latest' }" + label 'mogonet' // Saving outputs publishDir ( path: "${params.outdir}/${getPublishPath(task.process)}/${dataset_name}/${fold_path.name}", diff --git a/modules/rgcca/predict/main.nf b/modules/rgcca/predict/main.nf index c6dcf0e..431dd65 100644 --- a/modules/rgcca/predict/main.nf +++ b/modules/rgcca/predict/main.nf @@ -3,14 +3,10 @@ include { getPublishPath } from "${modulesDir}/functions" process RGCCA_PREDICT { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${dataset_name}-${fold_name}-${method}" debug "${params.debug}" - label 'process_single' - container "${ onSockeye ? - 'rgcca.sif' : - 'tonyliang19/rgcca:latest' }" - + label 'process_single' + label 'rgcca' publishDir ( path: "${params.outdir}/${getPublishPath(task.process)}/${dataset_name}/${fold_name}", mode: 'copy', @@ -23,7 +19,6 @@ process RGCCA_PREDICT { */ label 'low_mem' label 'cpu' - label 'codia' /* Given RGCCA itself has many methods available, so we have an extra input here diff --git a/modules/rgcca/preprocess/main.nf b/modules/rgcca/preprocess/main.nf index 4e5dab9..37126ee 100644 --- a/modules/rgcca/preprocess/main.nf +++ b/modules/rgcca/preprocess/main.nf @@ -22,14 +22,11 @@ include { getPublishPath } from "${modulesDir}/functions" // TODO: Replace name here process RGCCA_PREPROCESS { // temp variables to use - def onSockeye = workflow.projectDir.toString().contains('/scratch') // process level configuration debug "${params.debug}" // default is true tag "${dataset_name}" label 'process_single' - container "${ onSockeye ? - 'rgcca.sif' : - 'tonyliang19/rgcca:latest'}" + label 'rgcca' /* outputs store to outdir for saving and inspecting purpose */ publishDir ( diff --git a/modules/rgcca/select_feature/main.nf b/modules/rgcca/select_feature/main.nf index 951d11d..963211b 100644 --- a/modules/rgcca/select_feature/main.nf +++ b/modules/rgcca/select_feature/main.nf @@ -4,13 +4,10 @@ include { getPublishPath } from "${modulesDir}/functions" process RGCCA_SELECT_FEATURE { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${dataset_name}-design_${design}-ncomp_${ncomp}" debug true label 'process_low' - container "${ onSockeye ? - 'rgcca.sif' : - 'tonyliang19/rgcca:latest' }" + label 'rgcca' publishDir ( path: "${params.outdir}/${task.process.tokenize(':').join('/').toLowerCase()}/${dataset_name}", diff --git a/modules/rgcca/train/main.nf b/modules/rgcca/train/main.nf index 306f977..714232d 100644 --- a/modules/rgcca/train/main.nf +++ b/modules/rgcca/train/main.nf @@ -18,14 +18,11 @@ include { getPublishPath } from "${modulesDir}/functions" process RGCCA_TRAIN { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') // Identifier for each dataset and fold combination tag "${dataset_name}-${fold_path.name}-${method}-design_${design}-ncomp_${ncomp}" label 'process_single' // TODO: rename to the actual image name used - container "${ onSockeye ? - 'rgcca.sif' : - 'tonyliang19/rgcca:latest' }" + label 'rgcca' // Parse the output directory to migrate results to publishDir ( path: "${params.outdir}/${getPublishPath(task.process)}/${dataset_name}/${fold_path.name}", diff --git a/modules/sklearn/predict/main.nf b/modules/sklearn/predict/main.nf index 44343d8..eb8fe18 100755 --- a/modules/sklearn/predict/main.nf +++ b/modules/sklearn/predict/main.nf @@ -3,12 +3,9 @@ include { getPublishPath } from "${modulesDir}/functions" process SKLEARN_PREDICT { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${dataset_name}-${fold_name}-${model_name}" debug "${params.debug}" - container "${ onSockeye ? - 'sklearn.sif' : - 'tonyliang19/sklearn:latest' }" + label 'sklearn' publishDir ( path: "${params.outdir}/${getPublishPath(task.process)}/${dataset_name}/${fold_name}/${model_name}", diff --git a/modules/sklearn/preprocess/main.nf b/modules/sklearn/preprocess/main.nf index f27a4b3..c339ced 100755 --- a/modules/sklearn/preprocess/main.nf +++ b/modules/sklearn/preprocess/main.nf @@ -22,15 +22,12 @@ include { getPublishPath } from "${modulesDir}/functions" // TODO: Replace name here process SKLEARN_PREPROCESS { // temp variables to use - def onSockeye = workflow.projectDir.toString().contains('/scratch') // process level configuration debug "${params.debug}" // debugs true or false by param in MESSI.config tag "${dataset_name}" // identifier of process when ran in parallel // By var before to determine what container to use // Uses apptainer if true otherwise docker - container "${ onSockeye ? - 'sklearn.sif' : - 'tonyliang19/sklearn:latest' }" + label 'sklearn' label "process_low" diff --git a/modules/sklearn/select_feature/main.nf b/modules/sklearn/select_feature/main.nf index 34ecc3f..52a0b0e 100755 --- a/modules/sklearn/select_feature/main.nf +++ b/modules/sklearn/select_feature/main.nf @@ -4,13 +4,10 @@ include { getPublishPath } from "${modulesDir}/functions" process SKLEARN_SELECT_FEATURE { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${model_name}-${dataset_name}" debug true label 'process_high' - container "${ onSockeye ? - 'sklearn.sif' : - 'tonyliang19/sklearn:latest' }" + label 'sklearn' publishDir ( path: "${params.outdir}/${task.process.tokenize(':').join('/').toLowerCase()}/${dataset_name}", diff --git a/modules/sklearn/train/main.nf b/modules/sklearn/train/main.nf index 0250cb8..e3b3c59 100755 --- a/modules/sklearn/train/main.nf +++ b/modules/sklearn/train/main.nf @@ -18,13 +18,10 @@ include { getPublishPath } from "${modulesDir}/functions" process SKLEARN_TRAIN { // Vars stuff - def onSockeye = workflow.projectDir.toString().contains('/scratch') // Identifier for each dataset and fold combination tag "${dataset_name}-${fold_path.name}-${model_name}" // TODO: rename to the actual image name used - container "${ onSockeye ? - 'sklearn.sif' : - 'tonyliang19/sklearn:latest' }" + label 'sklearn' // Parse the output directory to migrate results to publishDir ( path: "${params.outdir}/${getPublishPath(task.process)}/${dataset_name}/${fold_path.name}", diff --git a/modules/split_train_test/main.nf b/modules/split_train_test/main.nf index 6d74028..658cd37 100644 --- a/modules/split_train_test/main.nf +++ b/modules/split_train_test/main.nf @@ -2,13 +2,10 @@ process SPLIT_TRAIN_TEST { // process metadata and configs // Vars stuff //def isRemote = workflow.containerEngine == 'apptainer' && !workflow.profile == 'standard' - def onSockeye = workflow.projectDir.toString().contains('/scratch') tag "${dataset_name}" label 'process_single' + label 'generic' debug "${params.debug}" - container "${ onSockeye ? - 'save_simulate.sif' : - 'tonyliang19/save_simulate:latest' }" publishDir ( path: "${params.outdir}/${task.process.tokenize(':').join('/').toLowerCase()}", //saveAS: { fn -> fn.endsWith(".log") ? "${id}/$fn" : fn }, @@ -44,4 +41,4 @@ process SPLIT_TRAIN_TEST { touch ${output_dir}/b.txt touch test.log """ -} \ No newline at end of file +} diff --git a/nextflow.config b/nextflow.config index 3228614..7f5b6c9 100644 --- a/nextflow.config +++ b/nextflow.config @@ -240,7 +240,12 @@ process { withLabel: caret_multimodal { container = 'tonyliang19/caret_multimodal:latest' } withLabel: generic { container = 'tonyliang19/save_simulate:latest' } withLabel: mogonet { container = 'tonyliang19/mogonet:latest' } - + withLabel: mofa { container = 'tonyliang19/mofa:latest' } + withLabel: rgcca { container = 'tonyliang19/rgcca:latest' } + withLabel: codia { container = 'tonyliang19/codia:latest' } + withLabel: integrao { container = 'tonyliang19/integrao:latest' } + withLabel: sklearn { container = 'tonyliang19/sklearn:latest' } + } // Pipeline information after exectuion From 09c4905dce366c2499e81b1db0602b7c62224a92 Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Mon, 31 Aug 2026 13:17:15 +0800 Subject: [PATCH 03/29] fix: changed wrong syntax, adding args to pass in splitting modules, cli and config --- modules/split_train_test/main.nf | 2 ++ .../resources/usr/bin/split_tr_te.py | 28 +++++++++++++------ nextflow.config | 1 + subworkflows/splitting/main.nf | 3 +- 4 files changed, 25 insertions(+), 9 deletions(-) diff --git a/modules/split_train_test/main.nf b/modules/split_train_test/main.nf index 658cd37..4d7b27c 100644 --- a/modules/split_train_test/main.nf +++ b/modules/split_train_test/main.nf @@ -18,6 +18,7 @@ process SPLIT_TRAIN_TEST { tuple val(dataset_name), path(mu_path) val(split_type) val(num_splits) + val(outcome_type) val(output_dir) output: @@ -31,6 +32,7 @@ process SPLIT_TRAIN_TEST { split_tr_te.py ${mu_path} \ --split_type=${split_type} \ --num_splits=${num_splits} \ + --outcome_type=${outcome_type} \ --output_dir=${dataset_name}/${output_dir} > \ ${dataset_name}/${dataset_name}-${task.process.tokenize(':')[-1].toLowerCase()}.log """ 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 224e3c9..faa070c 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 @@ -12,11 +12,13 @@ MDATA path to find MuData format Options: - --split_type=SPLIT_TYPE Type of split: skf, sgkf, logo [default: skf] - --num_splits=NUM_SPLITS Number of splits to generate [default: 10] - --seed=SEED Random number seed to reproduce [default: 329] - --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] + --split_type=SPLIT_TYPE Type of split: skf, sgkf, logo [default: skf] + --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] + --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] """ from docopt import docopt @@ -37,6 +39,8 @@ class SplitConfig: k: int seed: int identifier_col: str = "sample_name" + outcome_type: str = "classification" + event_col: str = "event" output_dir: str = "splits" split_txt_name: str = "fold" @@ -44,14 +48,19 @@ class SplitConfig: # ------------------------------------------------------------------------------ # Utility functions (functional, stateless) # ------------------------------------------------------------------------------ -def load_mudata(path, identifier_col): +def load_mudata(path, identifier_col, outcome_type="classification", event_col="os_event"): mdata = mudata.read(path) block_key = list(mdata.mod.keys())[0] block = mdata.mod[block_key] X = block.to_df() - y = block.obs["response"].values groups = block.obs[identifier_col].values + + if outcome_type == "survival": + y = block.obs[event_col].values + else: + y = block.obs["response"].values + return X, y, groups def create_splitter(split_type, k, seed): @@ -111,13 +120,16 @@ 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"), + outcome_type=args.get("--outcome_type", "classification"), output_dir=args["--output_dir"], split_txt_name=args["--split_txt_name"], ) # Load the MuData and extract X, y, groups # Groups default is to sample_name - X, y, groups = load_mudata(args["MDATA"], cfg.identifier_col) + X, y, groups = load_mudata(args["MDATA"], cfg.identifier_col, outcome_type=cfg.outcome_type, + event_col=cfg.event_col) engine = SplitEngine(cfg) folds = engine.run(X, y, groups) diff --git a/nextflow.config b/nextflow.config index 7f5b6c9..3458e90 100644 --- a/nextflow.config +++ b/nextflow.config @@ -104,6 +104,7 @@ 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" /* Threshold to convert the predicted probability to class */ threshold = 0.5 /* For feature selction workflow */ diff --git a/subworkflows/splitting/main.nf b/subworkflows/splitting/main.nf index 3329c10..4b84ca2 100644 --- a/subworkflows/splitting/main.nf +++ b/subworkflows/splitting/main.nf @@ -7,6 +7,7 @@ workflow SPLITTING { // Load params num_splits = params.k_fold_number // value of K fold output_dir = params.split_dir // output directory to store the splits , def is "splits" + outcome_type = params.outcome_type // Outcome type one of classification or survival, default is classification // Workflow of splitting starts here take: // ch_datasets // tuple of idendifier, path of mae data, path of mu data @@ -28,7 +29,7 @@ workflow SPLITTING { log.info "Using ${split_type} splitting strategy based on 'sample name' column" - SPLIT_TRAIN_TEST ( mu_data , split_type, num_splits, output_dir ) + SPLIT_TRAIN_TEST ( mu_data , split_type, num_splits, outcome_type, output_dir ) emit: splits_indices = SPLIT_TRAIN_TEST.out.splits_indices ch_logs = SPLIT_TRAIN_TEST.out.split_log From 0819287bb4d08c4091b983bd74e5144d20734dcb Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Mon, 31 Aug 2026 13:29:09 +0800 Subject: [PATCH 04/29] fix: lint conf syntax, included back sklearn --- nextflow.config | 14 +++++++------- subworkflows/cross_validation/main.nf | 3 ++- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/nextflow.config b/nextflow.config index 3458e90..6e0c26d 100644 --- a/nextflow.config +++ b/nextflow.config @@ -239,13 +239,13 @@ process { // Assigning module containers here withLabel: caret_multimodal { container = 'tonyliang19/caret_multimodal:latest' } - withLabel: generic { container = 'tonyliang19/save_simulate:latest' } - withLabel: mogonet { container = 'tonyliang19/mogonet:latest' } - withLabel: mofa { container = 'tonyliang19/mofa:latest' } - withLabel: rgcca { container = 'tonyliang19/rgcca:latest' } - withLabel: codia { container = 'tonyliang19/codia:latest' } - withLabel: integrao { container = 'tonyliang19/integrao:latest' } - withLabel: sklearn { container = 'tonyliang19/sklearn:latest' } + withLabel: generic { container = 'tonyliang19/save_simulate:latest' } + withLabel: mogonet { container = 'tonyliang19/mogonet:latest' } + withLabel: mofa { container = 'tonyliang19/mofa:latest' } + 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 } diff --git a/subworkflows/cross_validation/main.nf b/subworkflows/cross_validation/main.nf index 054d272..a7fcb44 100644 --- a/subworkflows/cross_validation/main.nf +++ b/subworkflows/cross_validation/main.nf @@ -29,7 +29,8 @@ include { printBanner } from "${modulesDir}/functions" def shouldRunPython() { return !( params.skip_mogonet && - params.skip_integrao + params.skip_integrao && + params.skip_sklearn ) } From 5f7cd2e902916070e9471935f01150a2ade2d123 Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Mon, 31 Aug 2026 13:30:10 +0800 Subject: [PATCH 05/29] rm: remove old files --- modules/goat/main.nf | 38 -- modules/goat/resource/usr/bin/run_goat.py | 0 modules/simulation/simulate_intersim/main.nf | 71 ---- .../resources/usr/bin/save_mudata.py | 29 -- .../resources/usr/bin/simulate_InterSIM.R | 355 ------------------ modules/simulation/simulate_mvn_data/main.nf | 75 ---- .../resources/usr/bin/debug_simulate_data.R | 202 ---------- .../resources/usr/bin/gen_simul_metadata.R | 41 -- .../resources/usr/bin/simulate_data.R | 160 -------- .../resources/usr/bin/unique_matrices.R | 39 -- 10 files changed, 1010 deletions(-) delete mode 100644 modules/goat/main.nf delete mode 100755 modules/goat/resource/usr/bin/run_goat.py delete mode 100644 modules/simulation/simulate_intersim/main.nf delete mode 100644 modules/simulation/simulate_intersim/resources/usr/bin/save_mudata.py delete mode 100755 modules/simulation/simulate_intersim/resources/usr/bin/simulate_InterSIM.R delete mode 100644 modules/simulation/simulate_mvn_data/main.nf delete mode 100755 modules/simulation/simulate_mvn_data/resources/usr/bin/debug_simulate_data.R delete mode 100755 modules/simulation/simulate_mvn_data/resources/usr/bin/gen_simul_metadata.R delete mode 100755 modules/simulation/simulate_mvn_data/resources/usr/bin/simulate_data.R delete mode 100755 modules/simulation/simulate_mvn_data/resources/usr/bin/unique_matrices.R diff --git a/modules/goat/main.nf b/modules/goat/main.nf deleted file mode 100644 index 1502d23..0000000 --- a/modules/goat/main.nf +++ /dev/null @@ -1,38 +0,0 @@ -// Nextflow process -// Author: Tony Liang - - -// Handles the GOAT method listed the Reference section of top-level README - -process TRAIN_GOAT { - // This label is defined in nextflow.config - label 'python' - label 'gpu' - debug true - - tag "${id}" - - input: - tuple val(id), path(mu_path) - each test_splits - output: - path('*mod*'), emit: mod - path('*prediction*'), emit: pred - path('*log*'), emit: log - script: - // Note this Python bin is directly from the container defined above ^ - // Some of the variables here that you don see defined explitcily in the input - // is likely defined in the top-level nextflow.config - """ - run_goat.py - """ - stub: - """ - echo ${id} - echo ${mu_path} - echo 'some text' > text.log - touch prediction.csv - touch model - """ - -} \ No newline at end of file diff --git a/modules/goat/resource/usr/bin/run_goat.py b/modules/goat/resource/usr/bin/run_goat.py deleted file mode 100755 index e69de29..0000000 diff --git a/modules/simulation/simulate_intersim/main.nf b/modules/simulation/simulate_intersim/main.nf deleted file mode 100644 index a0d2b69..0000000 --- a/modules/simulation/simulate_intersim/main.nf +++ /dev/null @@ -1,71 +0,0 @@ -/* - This process aims to simulate the data, takes in a map/dictionary of - parameters to control the data-generating process. - - Input: val(grid), combination of parameters - output_format, output format to generate - Whereas you could acess its parameters by grid., i.e. grid.number - - Possible parameters are: - - dataset_name # Unique identifier of each combination - - number # of observations - - latent_predictors # of latent predictors - - num_predictors # of variables) - - sigma Controls noise of data - - Output is then organized by two elements: - - ~~1 tuple of paths to x, y, z matrices in the sim_data channel~~ - - ~~1 path to the log file of the process~~ - -*/ - -process SIMULATE_INTERSIM { - // Vars stuff - //def isRemote = workflow.containerEngine == 'apptainer' && !workflow.profile == 'standard' - debug false - def onSockeye = workflow.projectDir.toString().contains('/scratch') - // process metadata and configs - tag "${grid.dataset_name}-${output_format}" - label 'process_low' - container "${ onSockeye ? - 'intersim.sif' : - 'tonyliang19/intersim:latest' }" - publishDir ( - path: "${params.outdir}/${task.process.tokenize(':')[-1].toLowerCase()}/${grid.dataset_name}", - saveAS: { file }, - mode: 'copy', - overwrite: true - ) - // Input goes here - input: - val grid // Combinations of parameters, key-value pair, access element by grid.xxx - // Possible output for downstream - output: - // The optional MUST be true here, since it should output one of MAE or MuData at a time - // This make sure it is unique and matching glob - tuple val(grid.dataset_name), path("${grid.dataset_name}_mae_data"), emit: sim_mae - tuple val(grid.dataset_name), path("${grid.dataset_name}.h5mu"), emit: sim_mu - tuple val(grid.dataset_name), path('*.log'), emit: sim_log - - script: - """ - simulate_InterSIM.R --dataset_name=${grid.dataset_name} \ - --number_obs=${grid.num_obs} \ - --effect=${grid.effect} \ - --sigma=${grid.sigma} \ - --corr=${grid.corr} \ - --noise=${grid.noise} \ - --number_noise_vars=${grid.num_noise_vars} > \ - ${grid.dataset_name}_${task.process.tokenize(':')[-1].toLowerCase()}.log - """ - // You need extra param to run stub like -stub - stub: - """ - mkdir ${grid.dataset_name}_mae_data - touch ${grid.dataset_name}_mae_data/experiments.hdf5 - touch ${grid.dataset_name}_mae_data/mae.rds - touch ${grid.dataset_name}.h5mu - touch ${grid.dataset_name}.log - """ - -} \ No newline at end of file diff --git a/modules/simulation/simulate_intersim/resources/usr/bin/save_mudata.py b/modules/simulation/simulate_intersim/resources/usr/bin/save_mudata.py deleted file mode 100644 index 6e32068..0000000 --- a/modules/simulation/simulate_intersim/resources/usr/bin/save_mudata.py +++ /dev/null @@ -1,29 +0,0 @@ -from mudata import MuData -from anndata import AnnData -import numpy as np -import pandas as pd - -def save_mudata(X_dict, meta_df, feat_names_dict, dataset_name): - # X_dict is dict of all the omics count data in array - # meta_df is dataframe of relevant metadata information like subject - # identifier, response - # feat_names_dict is the feature names of X_dict - - # Combine to mudata - mu_dict = {} - for b, mat in X_dict.items(): - var_names = feat_names_dict[b] - var = pd.DataFrame(var_names, columns=["feature"]) - ann = AnnData(X = mat, - obs = meta_df, - var = var) - ann.obs_names = meta_df['sample_name'] - ann.var_names = var_names - mu_dict[b] = ann - - # store to mudata - mdata = MuData(mu_dict) - # Write to file - output_path = f"{dataset_name}.h5mu" - mdata.write(output_path) - return None diff --git a/modules/simulation/simulate_intersim/resources/usr/bin/simulate_InterSIM.R b/modules/simulation/simulate_intersim/resources/usr/bin/simulate_InterSIM.R deleted file mode 100755 index 7235f40..0000000 --- a/modules/simulation/simulate_intersim/resources/usr/bin/simulate_InterSIM.R +++ /dev/null @@ -1,355 +0,0 @@ -#!/usr/bin/env Rscript - -# ============================================================================== -doc <- "This script is to use the InterSIM package to simulate synthetic data -based on TCGA ovarian cancer study. It creates dna methylation, gene expression -and protein data. We then apply a specific transformation to generate a -dummy bernoulli distribution response of mimicking patient has cancer or not - -Usage: - simulate_InterSIM.R [options] - -Options: - --help Display this help message - --dataset_name=DNAME Name of the dataset [default: empty] - --number_obs=N Number of observations to generate [default: 30] - --number_noise_vars=H Number of noise variables to create [default: 200] - --effect=EFFECT Cluster mean shift on each view [default: 2] - --sigma=SIGMA Covariance structure in the each omics's data, one of def, indep [default: indep] - --noise=NOISE Gaussian noise standard deviation [default: 1] - --corr=CORR Correlation between each omics, one of: 0, 0.5, 1 . [default: 0] - --transformation=TR Transformation to apply of the generated X to derive the response variable [default: rev_logit] -" -# ============================================================================== -# Parse cli arguments -opt <- docopt::docopt(doc) -# Load library -library(InterSIM) -library(dplyr) -library(magrittr) - -# Helper to convert mae to h5mu on mudata -save_h5mu <- function(mae, dataset_name) { - # Takes the MAE experiment and save this as MuData - exps <- mae@ExperimentList |> lapply(t) - col_data <- mae@colData |> as.data.frame() - feat_names <- mae@ExperimentList |> lapply(rownames) - # Then calls python code here - reticulate::use_python("/usr/bin/python") - # Gather the pipeline dir (THIS IS VERY UGGLY FIX) - bin_dir <- Sys.getenv("PATH") |> - strsplit(":") |> - unlist() |> - tail(1) - is_scratch <- stringr::str_detect(bin_dir, pattern = "scratch") - if (is_scratch) { - pipeline_dir <- gsub("/bin", "", bin_dir) - } else { - pipeline_dir <- "" - } - reticulate::source_python(here::here(pipeline_dir, "modules/simulation/simulate_intersim/resources/usr/bin/save_mudata.py")) - # Save it to mudata - save_mudata(exps, col_data, feat_names, dataset_name) -} - - -# Helper fun to generate the bernoulli distributed response variable from the -# counts X using certain criteria -# TODO: Need to better describe this criteria - -# Args: -# X: List of matrices denoting the count matrix of omics -# Each matrix should have dimension of n x P_j, where n is number of -# observations, and P_j is number of features in j-th omics -# criteria: Custom criteria of how to transform from X to get bernoulli -# distributed Y -# Options are: -# - composite_score (rowSums >= median(rowSums)) -# - rev_logit: logit^-1 on the X to get [0, 1] values on each -# entry of the cbinded matrices M = [m1 m2 m3, ...] then take -# Bern( 1 - rowMeans(M) ) -# tol: tolerance of comparing difference of sample var and population of the response -# generate_Y <- function(X, transformation=c("rev_logit", "composite_score"), tol=0.01, id_name="sample_name") { -# # SHOULD not create the y by clustering it, otherwise method will catch it 100%? -# transformation <- match.arg(transformation) -# message("Using transformation: ", transformation) -# # NOTE: Always using row here, since InterSim output row as subject, col as -# # feature - - -# # ======================= -# # THIS IS FOR DEBUG -# #X <- dd@ExperimentList |> lapply(t) -# # ======================= - -# if (transformation == "composite_score") { -# z <- rowSums(do.call(cbind, lapply(X, rowSums))) -# # Calculate a avg score of either mean or median -# median_score <- median(z) -# # Named vector of 1 and 0s -# response <- ifelse(z >= median_score, 1, 0) -# } - -# if (transformation == "rev_logit") { -# # TODO: the rowmeans could be bad, since its almost always > 0 -# # resulting in a not so fair bernoulli random variable -# z <- InterSIM::rev.logit(do.call(cbind, X)) |> -# rowMeans() -# # TODO: try z (more positives) or 1 - z (more negative) in prob -# response <- rbinom(n = length(z) , size = 1 , prob = 1 - z) -# # Named vector of 1 and 0s -# names(response) <- names(z) -# } -# # Check if the response follows bernoulli distribution with some tolerance -# sample_var <- var(response) -# phat <- mean(response) -# pop_var <- phat * (1 - phat) -# diff <- abs(pop_var - sample_var) -# is_bernoulli <- diff <= tol -# if (!is_bernoulli) warning("Difference of population and sample variance: ", round(diff, 3), " which exceeded tolerance of ", tol) - -# # Lastly assign the rownames to df as well -# meta_df <- response |> -# as.data.frame() |> -# tibble::rownames_to_column(var={{ id_name }}) |> -# select({{ id_name }}, response) -# # Manually assign the rownames back, since MAE uses rownames of colData to -# # match those colnames of the X matrix (P_j x n) -# rownames(meta_df) <- meta_df |> pull( {{ id_name }} ) -# return(meta_df) -# } - -generate_Y <- function(response_df, id_name="sample_name", response_name="response") { - # Maybe should not create the y by clustering it, otherwise method will catch it 100%? - if(!("subjects" %in% colnames(response_df) && "cluster.id" %in% colnames(response_df))) { - stop("The columns 'subjects' and 'cluster.id' must exist in the response dataset.") - } - meta_df <- response_df %>% - rename( - {{ id_name }} := subjects, - {{ response_name }} := cluster.id - ) %>% - # The response vector is numeric already - # INTERSIM gives 3 clusters, so we only keep the first two 1 and 2 - filter( !!sym(response_name) != 3) %>% - mutate( {{ response_name }} := !!sym(response_name) - 1) - # Manually assign the rownames back, since MAE uses rownames of colData to - # match those colnames of the X matrix (P_j x n) - rownames(meta_df) <- meta_df |> pull( {{ id_name }} ) - return(meta_df) -} - - -# ============================================================================ -# HANDLING THE X LIST OF OMICS MATRICES -# ============================================================================ -# TODO: many functions here should be moved to other files - -# Helper to apply gaussian noise to each col or row of a matrix -add_gaussian_noise <- function(x, noise_mean = 0, noise_sd = 1) { - # X could be a column vector or row vector - noise <- rnorm(length(x), mean=noise_mean, sd=noise_sd) - new_x <- x + noise - return(new_x) -} - -# Helper to generate H noise variables given matrix of n x p , where -# n is number of subjects, p is number of variables. -# By design, these noisy variables have higher sd than existing ones - -generate_gaussian_noise_vars <- function(omic_matrix, H=100, sd_multiplier = 1.5) { - # Convenient var - n <- nrow(omic_matrix) - # Calculate standard deviation of existing variables - existing_sd <- apply(omic_matrix, 2, sd, na.rm = TRUE) # This is a vector - # Purposedly let sd of noise vars higher - noise_data <- rnorm(n * H, mean = 0, sd = sd_multiplier * mean(existing_sd) ) - noise_vars_matrix <- matrix(noise_data, nrow=n, ncol=H) - # Generate H noise variables - # noise_vars_matrix <- replicate(H, { - # # For each new noise variable, use a standard deviation larger than existing ones - # noise_sd <- sd_multiplier * mean(existing_sd, na.rm = TRUE) - - # # Generate random values with increased standard deviation - # runif(n = nrow(matrix), min = min(matrix, na.rm = TRUE), max = max(matrix, na.rm = TRUE)) - # }, simplify = TRUE) - return(noise_vars_matrix) -} - -generate_bimodal_noise_vars <- function(n, H=100, - beta_shape1 = 2, beta_shape2 = 2, - scale_1 = 2, shift_1 = 1, - scale_2 = 2, shift_2 = 3, - bimodal = TRUE) { - - - # Generate beta-distributed data for n rows and H columns - noise_data <- rbeta(n*H, shape1 = beta_shape1, shape2 = beta_shape2) - # Stored to matrix - noise_vars_matrix <- matrix(noise_data, nrow = n, ncol = H) - if (!bimodal) { - message("Returning beta distributed noise vars") - return(noise_vars_matrix) - } - # Create bimodal effect by applying different scaling and shifting to the two halves - first_half <- 1:(n / 2) - second_half <- (n / 2 + 1):n - # Then modify these two halves - noise_vars_matrix[first_half, ] <- noise_vars_matrix[first_half, ] * scale_1 + shift_1 - noise_vars_matrix[second_half, ] <- noise_vars_matrix[second_half, ] * scale_2 + shift_2 - return(noise_vars_matrix) -} - - -# Get the X omics matrices and add numbers of noise variable in each matrix -# And add additional noise to existing vars (including those noise variables) -generate_X <- function(X_raw, meta_df, dataset_name, H, noise_mean=0, noise_sd=1, sd_multiplier=1.5) { - # Rename its prefix of dat. - X_names <- gsub("dat.", "", names(X_raw)) - # The sample names in clusters 1 and 2 - keep_samples <- rownames(meta_df) - X <- lapply(names(X_raw), function(omic_name, H, noise_mean, noise_sd, sd_multiplier) { - omic <- X_raw[[omic_name]] - # Keep relevant observations, since some belong to cluster 3 which is not included - if (!all(keep_samples %in% rownames(omic))) { - stop("Some keep_samples are not found in the omic matrix.") - } - x_mat <- omic[keep_samples, ] - # Some conveninent vars - n <- nrow(x_mat) - # Check for methylation data, since they are beta distributed of [0, 1] - is_beta_distributed <- all(x_mat >= 0 & x_mat <= 1) - if (is_beta_distributed) { - # This specifically handles methyl data - beta_shape = 2 - beta_scale = 2 - shift_1 = 1 - shift_2 = 3 - # Uses bimodal distribute noise - noise_vars_matrix <- generate_bimodal_noise_vars( - n=n, H=H, - beta_shape1=beta_shape, beta_shape2=beta_shape, - scale_1=beta_scale, shift_1=shift_1, - scale_2=beta_scale, shift_2=shift_2 - ) - } - - # Otherwise always use gaussian noise variables - noise_vars_matrix <- generate_gaussian_noise_vars(omic_matrix=x_mat, H=H, sd_multiplier=sd_multiplier) - # Append dummy name to columns - # TODO: this bit could be redundant of removing prefix? - #omic_name_no_prefix <- gsub("dat.", "", omic_name) - colnames(noise_vars_matrix) <- paste(omic_name, "noise_var", seq_len(H), sep="_") - # Then combine the noise variables to the var matrix - x_mat_full <- cbind(x_mat, noise_vars_matrix) - # https://stats.stackexchange.com/questions/144410/how-to-add-noise-to-a-random-variable-whose-range-is-the-unit-interval - # TODO: might need to check if this doing right here - # Then for each column add gaussian noise - x_mat_full <- apply(x_mat_full, 2, function(col) add_gaussian_noise(col, noise_mean=noise_mean, noise_sd=noise_sd)) - # Also append the dataset name in front of each column name - # And transpose the X to MultiAssayExperiment format - return(t(x_mat_full)) - }, H=H, noise_mean=noise_mean, noise_sd=noise_sd, sd_multiplier=sd_multiplier) - names(X) <- X_names - return(X) -} - - - -# Main entrace of the function -main <- function(dataset_name, - n, effect, noise, H=200, - cluster.sample.prop = c(0.45,0.45,0.1), - p.DMP=0.2, p.DEG=NULL, p.DEP=NULL, - sigma=c("indep", "def"), - corr=0, - transformation="rev_logit" -) { - - # Stop when no custom dataset name is provided - if (dataset_name == "empty") stop("Did not provided a custom dataset for simulation of intersim") - # Match args of sigma - sigma <- match.arg(sigma) - if (sigma == "def") { - sigma <- NULL - } - - - # Internal calculate a seed to make sure the data is reproducible - seed <- sum(base::utf8ToInt(dataset_name)) - set.seed(seed) - # First generate the count data from InterSIM - # TODO: The interSIM pkg doesnt have a way to change number of features in each omic - # fixed to their defaults .... - - # Given cluster propotions are c(0.45, 0.45, 0.1), where last cluster is always dropped after creation - # So need to adjust that raw n to cancel this effect and having enough obsercations as stated. - # Using this formula: n* = ceiling(raw_n / 0.9) - # For example, if one want to simulate n = 50, then n* need to be ceiling(50 / 0.9) = 56 - # Then, 0.45 * 56 = 25.2 , 0.1 * 56 = 5.6 - # We can then only keep floor(25.2 + 25.2 ) = 50 which yields original n required - adjusted_n <- ceiling(n / 0.9) - dat <- InterSIM(n.sample=adjusted_n, - cluster.sample.prop=cluster.sample.prop, - delta.methyl = effect, delta.expr = effect, delta.protein = effect, - p.DMP=p.DMP,p.DEG=p.DEG, p.DEP=p.DEP, - sigma.methyl=sigma, sigma.expr=sigma, sigma.protein=sigma, - cor.methyl.expr=corr, cor.expr.protein=corr) - - - # Ignore its cluster assignment for now - n_list <- length(dat) - # We retaining only cluster 1 and 2 subjects, so need to first process meta then on X - Y_df <- generate_Y(response_df=dat[[n_list]]) - # Process the X as well with suitable parameters to control noise generation - # Let X be a J length list of n (row) x p (column) matrix - # 1. Generate H noise variables of either: - # a) Normal distributed with mean 0 , sd = sd_multiplier * mean(existing_sd of each column) of current omic and iterate for J omics - # This way noise has a higher spread than actual signal variables - # - # b) Multimodal distributed derived from Beta Distribution shapes are fixed to be same - # - # - # 2. Add only gaussian noise of mean = 0, sd = noise_sd , where noise sd is cli param to vary - # Note: this get added to those previous noise variables, so could be double source of noise - # And it also gets added to non normally distributed omics like the ones of Methylation - # which is stricly Beta distributed. - X <- generate_X(X_raw = dat[1:n_list - 1], meta_df=Y_df, dataset_name=dataset_name, H=H, noise_sd=noise) - # Construct the X and Y here - mae <- MultiAssayExperiment::MultiAssayExperiment( - experiments = X, - colData = Y_df - ) - - - # Then should convert mae to mudata h5mu - #if (tolower(output_format) == "mudata") { - - # Directly saves h5mu - save_h5mu(mae, dataset_name) - #} - - #if (tolower(output_format) == "mae") { - # Also saving it as mae - # Directly saves mae - MultiAssayExperiment::saveHDF5MultiAssayExperiment(mae, prefix="", - dir=paste0(dataset_name, "_", "mae_data"), - replace = T) - #} - return(mae) -} - - -# Call the main function -dat <- main( - dataset_name = opt$dataset_name, - n = as.numeric(opt$number_obs), - H = as.numeric(opt$number_noise_vars), - effect = as.numeric(opt$effect), - sigma = opt$sigma, - corr = as.numeric(opt$corr), - transformation = opt$transformation, - noise = as.numeric(opt$noise) -) -# View the output of this -print(dat) \ No newline at end of file diff --git a/modules/simulation/simulate_mvn_data/main.nf b/modules/simulation/simulate_mvn_data/main.nf deleted file mode 100644 index e803a9c..0000000 --- a/modules/simulation/simulate_mvn_data/main.nf +++ /dev/null @@ -1,75 +0,0 @@ -/* - This process aims to simulate the data, takes in a map/dictionary of - parameters to control the data-generating process. - - Input: val(grid), combination of parameters - Whereas you could acess its parameters by grid., i.e. grid.number - - Possible parameters are: - - dataset_name # Unique identifier of each combination - - number # of observations - - latent_predictors # of latent predictors - - num_predictors # of variables) - - sigma Controls noise of data - - Output is then organized by two elements: - - ~~1 tuple of paths to x, y, z matrices in the sim_data channel~~ - - ~~1 path to the log file of the process~~ - -*/ - -process SIMULATE_MVN_DATA { - // Vars stuff - //def isRemote = workflow.containerEngine == 'apptainer' && !workflow.profile == 'standard' - debug false - def onSockeye = workflow.projectDir.toString().contains('/scratch') - // process metadata and configs - tag "${grid.dataset_name}" - label 'process_medium' - container "${ onSockeye ? - 'save_simulate.sif' : - 'tonyliang19/save_simulate:latest' }" - label 'mae_mu' - publishDir ( - path: "${params.outdir}/${task.process.tokenize(':')[-1].toLowerCase()}/${grid.dataset_name}", - saveAS: { file }, - mode: 'copy', - overwrite: true - ) - // Input goes here - input: - val grid // Combinations of parameters, key-value pair, access element by grid.xxx - // Possible output for downstream - output: - // The optional MUST be true here, since it should output one of MAE or MuData at a time - tuple val(grid.dataset_name), path("${grid.dataset_name}*mae*"), emit: sim_mae - tuple val(grid.dataset_name), path("${grid.dataset_name}*.h5mu"), emit: sim_mu - tuple val(grid.dataset_name), path('*.log'), emit: sim_log - - script: - """ - echo -e "\nThis is try ${grid.dataset_name}" - simulate_data.R --dataset_name=${grid.dataset_name} \ - --number=${grid.num_obs} \ - --num_predictors=${grid.num_predictors} \ - --block_num=${grid.block_num} \ - --latent_predictors=${grid.latent_predictors} \ - --sigma=${grid.sigma} \ - --sy=${grid.sy} \ - --sp=${grid.sp} \ - --u_std=${grid.u_std} \ - --fct_str=${grid.fct_str} > \ - ${grid.dataset_name}_${task.process.tokenize(':')[-1].toLowerCase()}.log - echo -e "\nDone with ${grid.dataset_name}" - """ - // You need extra param to run stub like -stub - stub: - """ - mkdir ${grid.dataset_name}_mae_data - touch ${grid.dataset_name}_mae_data/experiments.hdf5 - touch ${grid.dataset_name}_mae_data/mae.rds - touch ${grid.dataset_name}.h5mu - touch ${grid.dataset_name}.log - """ - -} \ No newline at end of file diff --git a/modules/simulation/simulate_mvn_data/resources/usr/bin/debug_simulate_data.R b/modules/simulation/simulate_mvn_data/resources/usr/bin/debug_simulate_data.R deleted file mode 100755 index 20387b7..0000000 --- a/modules/simulation/simulate_mvn_data/resources/usr/bin/debug_simulate_data.R +++ /dev/null @@ -1,202 +0,0 @@ -#!/usr/bin/env Rscript - -# Doc section -------------------------------------------------------------- -'This is the script to simulate data with X block , Z block, and a response -vector R. The output is a list in R. By default it writes the rds file -to current directory. - -Author: Tony Liang - -Usage: - simulate_data.R [options] - -Options: - -n N --number=N Number of observations [default: 200] - --num_predictors=P Number of predictors to use [default: 500] - -m M --block_num=M Number of blocks to generate [default: 3] - --latent_predictors=IMP Number of latent predictors [default: 30] - -s SIGMA --sigma=SIGMA Noise strength to add into response [default: 39] - --sy=SY Standard deviation of Y [default: 1] - --sp=SP Standard deviation of block components [default: 3] - --u_std=U_STD Standard deviation of normal distributed U [default: 1] - --ft_str=FCT_STR Factor strength [default: 7] - --task=TASK Type of response, continous/binary/categorical [default: continous] - --tr=TRANSFORM Transformation on response, one of sigmoid or softmax [default: sigmoid] - --output_format=OUT_FMT Type of output data to write. One of MAE or MuData [default: MAE] - --name=NAME Name of data to save [default: sim_data] -' -> doc -cat("\nThis is the chr ver\n") -opt_chr <- docopt::docopt(doc) -print(opt_chr) -cat("\nThis is the num ver\n") -opt <- lapply(opt_chr, function(x) ifelse(grepl("^\\d+\\.?\\d*$", x), - as.numeric(x), x)) -print(opt) -#============================================================================== - -# number <- n <- 200 -# num_predictors <- p <- 300 -# m <- blocks_num <- 3 -# p_imp <- latent_predictors <- 30 -# sigma <- 39 -# sy <- 1 -# sp <- 3 -# u_std <- 1 -# ft_str <- factor_strength <- 7 -# #task <- c("continous", "binary", "categorical") -# task <- "binary" -# tr <- "sigmoid" -# #tr <- c("sigmoid", "softmax") -# name <- prefix <- "sim_data" -# output_format <- "rds" - - - -# library(magrittr) -# # Parse above doc -# opt_chr <- docopt::docopt(doc) -# # Source common helpers -# source(here::here("modules/R/generic_helpers.R")) -# load_helpers(helper_path = "modules/R/simulate_data/helpers/") -# # Convert all options to numeric -# opt <- lapply(opt_chr, function(x) ifelse(grepl("^\\d+\\.?\\d*$", x), -# as.numeric(x), x)) - -# # Generating process of data - -# # TODO: consider add another block of "omics", so another matrix -# # TODO: consider combining the px , pz and possibly the extra block to same "base" number of features -# # such that px = base + c1 , pz = base + c2 , pw = base + c3 -# # TODO: Now same thing happens with sy, sx, sz, so all needs to be fixed - -# number <- n <- 200 -# num_predictors <- p <- 500 -# p_imp <- 30 -# sigma <- 39 -# sy <- 1 -# sx <- sz <- sw <- 3 -# u_std <- 1 -# ft_str <- factor_strength <- 7 -# task <- c("continous", "binary", "categorical") -# tr <- c("sigmoid", "softmax") -# name <- prefix <- "sim_data" - -# simulate_data <- function(n, p, p_imp, sigma, -# sy, sx, sz, sw,u_std, ft_str, -# task = c("continous","binary", "categorical"), -# tr=c("sigmoid", "softmax"), -# output_format=c("mae", "mudata", "csv", "rds"), -# threshold=0.5) { -# # Match arguments -# task <- match.arg(task) -# tr <- match.arg(tr) -# output_format <- match.arg(output_format) -# # get predictors -# num_ps <- getNumPredictors(p = p, ft_str = ft_str, sigma=sigma) -# px <- num_ps$px -# pz <- num_ps$pz -# pw <- num_ps$pw -# args_used <- c(as.list(environment())) -# logging_params(args_used) # custom function to format and log -# # Record time to track execution time -# start_time <- Sys.time() - -# # Check p_imp not greater than any of px or pz -# if ( p_imp > max(px, pz, pw)) { -# print("p_imp cannot be > px or pz or pz, using default: 30") -# p_imp <- 30 -# } -# # Simulate data based on the factor model -# x = matrix(rnorm(n*px), n, px) -# z = matrix(rnorm(n*pz), n, pz) -# w = matrix(rnorm(n*pw), n, pw) -# U = matrix(rep(0, n*p_imp), n, p_imp) - -# # Relate matrices by u -# for (m in seq(p_imp)){ -# u = rnorm(n, sd = u_std) -# x[, m] = x[, m] + sx*u -# z[, m] = z[, m] + sz*u -# w[, m] = w[, m] + sw*u -# U[, m] = U[, m] + sy*u - -# } -# # Center and not scale these accordingly -# x = scale(x, center=TRUE, scale=FALSE) -# z = scale(z, center=TRUE, scale=FALSE) -# w = scale(w, center=TRUE, scale=FALSE) - -# # Assign column names to use 1 ... n -# colnames(x) = paste0("x", 1:ncol(x)) -# colnames(z) = paste0("z", 1:ncol(z)) -# colnames(w) = paste0("w", 1:ncol(w)) - -# # Assign rownames to be some string -# common_rows = 1:nrow(x) -# rownames(x) <- rownames(z) <- rownames(w) <- paste0("pat-", common_rows) - -# # Create beta matrix and Y -# beta_U = c(rep(ft_str, p_imp)) -# mu_all = U %*% beta_U -# y <- mu_all + sigma * rnorm(n) -# # Convert y to suitable outcome with user-transformation -# #tr <- match.arg(tr) -# threshold = 0.5 -# y_temp <- y %>% -# transformation(tr=tr) %>% -# (function(y) ifelse(y < threshold, "no", "yes")) %>% -# factor(levels = c("no", "yes")) -# # task <- match.arg(task) -# if (tolower(task) == "binary") { -# y <- y_temp %>% -# as.integer() - 1 -# } - -# if (tolower(task) == "categorical") { -# y <- y_temp -# } - -# elapsed <- Sys.time() - start_time -# cat("\nTime taken to simulate data:", round(elapsed, 6), "seconds\n") -# cat("\nData simulated!\n") -# return(dat=list(blocks = list(x=x, z=z, w=w), response=y)) -# } - -# main <- function(number, num_predictors, p_imp, -# sigma, sy, sx, sz, sw, u_std, -# factor_strength, task, tr, output_format, -# name="sim_data", prefix="sim_data") { - -# cat("Generating simulation... \n") -# # Invoke simulate data - -# dat <- simulate_data(n=number, p=num_predictors, -# p_imp=p_imp, sigma=sigma, sy = sy, sx=sx, sz=sz, sw=sw, -# u_std=u_std, ft_str=factor_strength, task=task,tr=tr, -# output_format=output_format) -# # Time to write data -# write_start <- Sys.time() -# #cat("\nWriting to disk\n") -# #lapply(names(dat), function(name) saveFile(object=dat[[name]], -# # name=name)) -# # ------------------------------------- -# # To MAE or MuData -# saveFile(dat, name = name, output_format=output_format, prefix=prefix) - -# logging_write_disk(write_start = write_start) -# # Write metadata to file as well -# return(dat) -# } - -# # Set seed to guarantee reproducible result (DELETE Later) -# set.seed(329) - -# dat <- main(number = opt$number, num_predictors = opt$num_predictors, -# p_imp = opt$latent_predictors, sigma = opt$sigma, -# sy = opt$sy, sx = opt$sx, sz = opt$sz, sw = opt$sw, -# u_std = opt$u_std, factor_strength = opt$ft_str, -# task = opt_chr$task, tr=opt_chr$tr, output_format=opt$output_format -# ) - - - diff --git a/modules/simulation/simulate_mvn_data/resources/usr/bin/gen_simul_metadata.R b/modules/simulation/simulate_mvn_data/resources/usr/bin/gen_simul_metadata.R deleted file mode 100755 index 19cc11d..0000000 --- a/modules/simulation/simulate_mvn_data/resources/usr/bin/gen_simul_metadata.R +++ /dev/null @@ -1,41 +0,0 @@ -library(magrittr) -# Helper to convert z to y -convert_z2y <- function(z, tr, task, response_name) { - # This makes it binary factor of 1 and 0s - y_bin_num <- z %>% - transformation(tr=tr) %>% - rbinom(length(z), 1, .) %>% # Bernoulli random var - as.factor() - # Determine by task and return suitable response y - if (tolower(task) == "binary") return(y_bin_numeric) - if (tolower(task) == "categorical") { - y <- y_bin_num %>% - tibble::as_tibble_col(column_name = response_name) %>% - dplyr::mutate({{ response_name }} := case_when( - !!sym(response_name) == 1 ~ "yes", - !!sym(response_name) == 0 ~ "no") - ) %>% - dplyr::mutate(!!response_name := factor(!!sym(response_name))) - return(y) - } -} - -# Generate dummy metadata, note response should be contained here -# and sample_names to be common rownames -gen_simul_metadata <- function(blocks, z, tr, task, id_name="sample_names", - response_name="response") { - # Convert this z to y response depending on task liked - # one of binary (num) or categorical (chr) - # TODO: need to put this elsewhere and clearer to denote what - # to transformation applied to Z to y - y <- convert_z2y(z, tr, task, response_name) - # Create dummy metadata - age <- sample(18:80, size = length(z), replace = TRUE) - sample_names <- intersect |> - Reduce(lapply(blocks, rownames)) - - # Combine these together in a df - metadata <- data.frame(sample_names, y, age) - names(metadata) <- c(id_name, response_name, "age") - return(metadata) -} diff --git a/modules/simulation/simulate_mvn_data/resources/usr/bin/simulate_data.R b/modules/simulation/simulate_mvn_data/resources/usr/bin/simulate_data.R deleted file mode 100755 index ad20fc8..0000000 --- a/modules/simulation/simulate_mvn_data/resources/usr/bin/simulate_data.R +++ /dev/null @@ -1,160 +0,0 @@ -#!/usr/bin/env Rscript - -# Doc section -------------------------------------------------------------- -'This is the script to simulate data with X block , Z block, and a response -vector R. The output is a list in R. By default it writes the rds file -to current directory. - -Author: Tony Liang - -Usage: - simulate_data.R [options] - -Options: - -n N --number=N Number of observations [default: 200] - --num_predictors=P Number of predictors to use [default: 500] - -m M --block_num=M Number of blocks to generate [default: 3] - --latent_predictors=IMP Number of latent predictors [default: 30] - -s SIGMA --sigma=SIGMA Noise strength to add into response [default: 39] - --sy=SY Standard deviation of Y [default: 1] - --sp=SP Standard deviation of block components [default: 3] - --u_std=U_STD Standard deviation of normal distributed U [default: 1] - --fct_str=FCT_STR Factor strength [default: 7] - --task=TASK Type of response, binary/categorical [default: categorical] - --tr=TRANSFORM Transformation on response, one of sigmoid or softmax [default: sigmoid] - --dataset_name=DATASET_NAME Name of data to save [default: sim_data] - --seed=SEED Seed to reproduce [default: 1] - --y_name=Y_NAME Column name for the response var [default: response] -' -> doc - -# Load libraries -library(dplyr) -library(magrittr) -library(here) - -# TODO: very uggly fix -bin_dir <- Sys.getenv("PATH") |> - strsplit(":") |> - unlist() |> - tail(1) -pipeline_dir <- gsub("/bin", "", bin_dir) -print(pipeline_dir) -# Loading scripts -source(here(pipeline_dir, "bin/rhelpers.R")) # This is included in nextflow bin path -# Load utils specific to simulation data -rp <- resource_helper_path(here(pipeline_dir, "modules/simulation/simulate_mvn_data")) -source(here(rp, "gen_simul_metadata.R")) -source(here(rp, "unique_matrices.R")) -# Loading generic utils -load_utils(here(pipeline_dir, "bin/logging")) -load_utils(here(pipeline_dir, "bin/preprocessing")) -source(here(pipeline_dir, "bin/savers/saveFile.R")) -# ============================================================================ -# Parse above doc -opt_chr <- docopt::docopt(doc) -opt <- opt2num(opt_chr) -# Functions ----------------------------------------------------------------- -# Generating process of data -# Source from https://www.pnas.org/doi/full/10.1073/pnas.2202113119 -simulate_data <- function( - n, p, m, p_imp, sigma, sy, - sp, u_std, fct_str, - task = c("continuous","binary", "categorical"), - tr=c("sigmoid", "softmax"), - y_name - ) { - # Match arguments ---------------------------------------------------------- - task <- match.arg(task) - tr <- match.arg(tr) - # Logging stuff ------------------------------------------------------------ - args_used <- c(as.list(environment())) - logging_params(args_used) # custom function to format and log - # Record time to track execution time - start_time <- Sys.time() - # Check p_imp not greater than any of px or pz - if ( p_imp > p) { - print("p_imp cannot be > p, using default: 30") - p_imp <- 30 - } - - # Instantiation ----------------------------------------------------------- - blocks <- unique_matrices(m=m, n=n, p=p, mu1=2, mu2=8, shift=TRUE) - U = matrix(rep(0, n*p_imp), n, p_imp) - - # Relate U until p_imp columns in each mat of the mat_list - for (j in seq(p_imp)) { - # Random noise with sd of u_std (default 1) - u = rnorm(n, sd = u_std) - # Add noise to each j-th column that are "latent predictor" in block - blocks <- lapply(blocks, function(block) { - block[, j] <- block[, j] + sp * u - return(block) - }) - # Do same thing to the U matrix - U[, j] = U[, j] + sy * u - } - # Center and not scale these accordingly - blocks <- lapply(blocks, scale, center=TRUE, scale=FALSE) - # Create beta matrix and Y - beta_U = c(rep(fct_str, p_imp)) - mu_all = U %*% beta_U - # Continuous way, transform it later - z <- mu_all + sigma * rnorm(n) - # Metadata generation ----------------------------------------------------- - metadata <- gen_simul_metadata(blocks=blocks, z=z, tr=tr, task=task, response_name=y_name) - # More logging to exit - elapsed <- Sys.time() - start_time - cat("\nTime taken to simulate data:", round(elapsed, 6), "seconds\n") - cat("\nData simulated!\n") - # TODO: This seems a bit hardcoded, but should have blocks that have all the omics - # while metadata should contain at least two columns, sample_names and response - return(dat=list(blocks=blocks, metadata=metadata)) -} - -# Main entrance of the scripts -main <- function(number, num_predictors, blocks_num, - p_imp, sigma, sy, sp, u_std, - factor_strength, task, tr, - dataset_name, y_name, prefix="") { - - cat("Generating simulation... \n") - - # Internal calculate a seed to make sure the data is reproducible - seed <- sum(base::utf8ToInt(dataset_name)) - set.seed(seed) - # Invoke simulate data - dat <- simulate_data(n=number, p=num_predictors, m=blocks_num, - p_imp=p_imp, sigma=sigma, sy = sy, sp=sp, - u_std=u_std, fct_str=factor_strength, task=task, - tr=tr, y_name=y_name) - # Time to write data - write_start <- Sys.time() - # To MAE and MuData - saveFile(dat, name=dataset_name, prefix=prefix, output_format="MAE") - saveFile(dat, name=dataset_name, prefix=prefix, output_format="MuData") - # Log to end - logging_write_disk(write_start = write_start) - # Write metadata to file as well - return(dat) -} - -# Set seed to guarantee reproducible result (DELETE Later) -#SEED = opt$seed -#set.seed(SEED) -dat <- main( - number = opt$number, - num_predictors = opt$num_predictors, - blocks_num = opt$block_num, - p_imp = opt$latent_predictors, - sigma = opt$sigma, - sy = opt$sy, - sp = opt$sp, - u_std = opt$u_std, - factor_strength = opt$fct_str, - task = opt_chr$task, - tr = opt_chr$tr, - dataset_name = opt_chr$dataset_name, - y_name = opt_chr$y_name - ) - - diff --git a/modules/simulation/simulate_mvn_data/resources/usr/bin/unique_matrices.R b/modules/simulation/simulate_mvn_data/resources/usr/bin/unique_matrices.R deleted file mode 100755 index dc1379d..0000000 --- a/modules/simulation/simulate_mvn_data/resources/usr/bin/unique_matrices.R +++ /dev/null @@ -1,39 +0,0 @@ -# Helpers for unique matrices -gen_p <- function(m, n , p, size=1) { - mu <- p / n - noise <- floor(rnorm(1, mean=mu, sd=p%/%m)) + m - unique_p <- p + noise - return(unique_p) -} - -# Could shift means by nth observations now -gen_rand_data <- function(n, unique_p, mu1, mu2, shift) { - # index of nth number to have different mu - j <- floor(n / 2) - if (shift == FALSE) { - # If want not shift then we have same means - mu1 <- mu2 <- 5 # randomly chose 5 - } - means <- c(rep(mu1, j), rep(mu2, n-j)) - dat <- rnorm(n * unique_p ,mean=rep(means, unique_p)) - return(dat) -} - -# Main function ---------------------------------------------- -unique_matrices <- function(m, n, p, mu1=2, mu2=8, shift=TRUE) { - mat_list <- list() - # m should not be greater than 26 though - for (i in 1:m) { - label <- LETTERS[i] # Use letters as unique labels (a, b, c, ...) - bname <- paste0(label, "_BLOCK") - #unique_p <- getNumPredictors(p, ft_str, sigma) - unique_p <- gen_p(m, n, p) - dat <- gen_rand_data(n, unique_p, mu1=mu1, mu2=mu2, shift=shift) - mat <- matrix(dat, n, unique_p) - colnames(mat) <- paste0(label, 1:ncol(mat)) - rownames(mat) <- paste0("pat-", 1:nrow(mat)) - mat_list[[bname]] <- mat - } - return(mat_list) -} - From 0ac61acdd23d4472779ba1bab32cf0213a3b249a Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Mon, 31 Aug 2026 13:31:27 +0800 Subject: [PATCH 06/29] rm: remove redundant staled configs --- nextflow.config | 37 ++----------------------------------- 1 file changed, 2 insertions(+), 35 deletions(-) diff --git a/nextflow.config b/nextflow.config index 6e0c26d..a7d56db 100644 --- a/nextflow.config +++ b/nextflow.config @@ -27,9 +27,10 @@ params { max_memory = '128.GB' max_cpus = 8 max_time = '240.h' + + // Input parameters // Options to control behavior of workflows - // simulate = true // Default to not run simulations, more specific params at subworkflow/simulation/params.config train = true // To false if just want to test the splitting and not go further at train stage // This setting is higher than the one inside a .nf file // Param to control if minimal setting or some more thorough test @@ -38,8 +39,6 @@ params { // Options to run part of workflows // TODO: redundant parameters here - runSimulation = false // Run simulation of data or not (default: false) - runSimple = false // This is used in creating simulation combination runInnerCV = false // Default to run inner cv in all possible methods runAllMethods = true // Run all multiomics integration methods, overrides runPython or runR runPython = true // Run methods implemented in Python only @@ -67,38 +66,6 @@ params { he_base_dim = 100 // Base dimension for the HE layer, default to 100 for real data, 1 for simulated data /* Prepare data related */ - // For simulation use, these params makes a grid of values to simulate on - // For those singleton list, could add more items if want to test it, but - // will make more possible combinations out (so more jobs to launch) - skip_simulation = true - skip_sim_intersim = true - skip_sim_MVN = true - /* Common simulation params */ - num_obs = [50, 100, 500] // Numbers of observations to create - dataset_base_name = "sim-data" // This is appended to a number from 1 to ... like sim-data - output_formats = ["MAE", "MuData"] // Formats of data copies to write to - y_name = "response" // Column name of the response variable - /* InterSIM related simulation params */ - intersim_sigma = ["def"] - //intersim_sigma = ["def", "indep"] // Within block covariance - intersim_corr = [0, 0.5, 1] // Between block covariance - intersim_effect = [0, 0.5, 1] // Cluster mean of each omic - intersim_noise = [1] - //intersim_noise = [1, 10, 100] // Standard deviation of normal distribution noise with mean 0 - intersim_num_noise_vars = [200] - //intersim_num_noise_vars = [100 , 200 , 500] // Number of noise variables to create in each omics - /* Cooperative learning mvn simulation params */ - // runSimple = true // True use only number, num_pred and output_format, otherwise use all params below - num_predictors = [500, 1000, 2500] // Base number of predictors in each omics - block_num = [3, 4, 5] // Numbers of omics to produce - latent_predictors = [30] // Numbers of latent predictors - mvn_sigma = [39] // Noise in the response variable - sy = [1] // Std of Y - sp = [3] // Std of each block (omics) - u_std = [1] // Std of U vector - fct_str = [7] // Strength of betas with y - task = ["categorical"] // Type of response to generate, another option in "binary" - tr = ["sigmoid"] // Transformation to apply on z to y, another option is "softmax" /* Splitting related */ split_type = "skf" // Type of splitting to perform, one of "skf" (stratified k fold) or // or "sgkf" (stratidied group k fold ) "logo" (leave one group out) From 2da5e23ee90d685a3953ab4881cd205d2b13186d Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Mon, 31 Aug 2026 15:42:12 +0800 Subject: [PATCH 07/29] fix: rename metric labels --- .../resources/usr/bin/calculate_metrics.py | 32 +++++++++++-------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/modules/calculate_metrics/resources/usr/bin/calculate_metrics.py b/modules/calculate_metrics/resources/usr/bin/calculate_metrics.py index 559ea3a..b8b6188 100755 --- a/modules/calculate_metrics/resources/usr/bin/calculate_metrics.py +++ b/modules/calculate_metrics/resources/usr/bin/calculate_metrics.py @@ -42,29 +42,35 @@ def calculate_metrics(group, threshold, average='binary', round_digit=3): y_pred = np.where(phat >= threshold, 1, 0) # Initialize dict for metrics metrics = { - 'auc': 0.0, + 'roc_auc': 0.0, + 'pr_auc': 0.0, + 'f1_score': 0.0, + 'log_loss': 0.0, + 'matthews_corrcoef': 0.0, 'accuracy': 0.0, 'balanced_accuracy': 0.0, 'precision': 0.0, - 'average_precision_score': 0.0, 'recall': 0.0, - 'f1_score': 0.0, - 'log_loss': 0.0, + } try: # Check if both classes are present in the current group if len(set(y_true)) > 1: # More than one unique value in y_true # AUC requires predicted probability (phat) not y_pred - metrics['auc'] = roc_auc_score(y_true, phat) - metrics['accuracy'] = accuracy_score(y_true, y_pred) - metrics['balanced_accuracy'] = balanced_accuracy_score(y_true, y_pred) - metrics['precision'] = precision_score(y_true, y_pred, average=average) - metrics['average_precision_score'] = average_precision_score(y_true, y_pred) - metrics['recall'] = recall_score(y_true, y_pred, average=average) - metrics['f1_score'] = f1_score(y_true, y_pred, average=average) - # For the losses as well - metrics['log_loss'] = log_loss(y_true, y_pred) + metrics['roc_auc'] = roc_auc_score(y_true, phat) + # Same thing with precision-recall auc + metrics['pr_auc'] = average_precision_score(y_true, phat) + # Rest could use the predicted class (y_pred) to calculate metrics + metrics['f1_score'] = f1_score(y_true, y_pred, average=average) + metrics['log_loss'] = log_loss(y_true, phat) + metrics['matthews_corrcoef'] = matthews_corrcoef(y_true, y_pred) + metrics['accuracy'] = accuracy_score(y_true, y_pred) + metrics['balanced_accuracy'] = balanced_accuracy_score(y_true, y_pred) + metrics['precision'] = precision_score(y_true, y_pred, average=average) + metrics['recall'] = recall_score(y_true, y_pred, average=average) + + else: raise ValueError("Only one class present") From 939ac4c042e0f9cfcadaed4e03f71a1f4fdf5704 Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Mon, 31 Aug 2026 15:42:41 +0800 Subject: [PATCH 08/29] rm: remove deprecated workflows --- subworkflows/methods/sgmr/main.nf | 12 -- subworkflows/prepare_data/main.nf | 3 - subworkflows/simulation/main.nf | 193 ------------------------------ 3 files changed, 208 deletions(-) delete mode 100644 subworkflows/methods/sgmr/main.nf delete mode 100644 subworkflows/simulation/main.nf diff --git a/subworkflows/methods/sgmr/main.nf b/subworkflows/methods/sgmr/main.nf deleted file mode 100644 index baf83de..0000000 --- a/subworkflows/methods/sgmr/main.nf +++ /dev/null @@ -1,12 +0,0 @@ -// This looks very hard-coded, TODO: fix this later and improve it -def method_name = "sgmr" - -workflow SGMR { - take: - mae_copy - main: - // TODO: This method is still boilerplate - log.info "This is SGMR method name: ${method_name}" - emit: - csv_results = Channel.empty() -} \ No newline at end of file diff --git a/subworkflows/prepare_data/main.nf b/subworkflows/prepare_data/main.nf index c8d84fc..60c89b0 100644 --- a/subworkflows/prepare_data/main.nf +++ b/subworkflows/prepare_data/main.nf @@ -1,6 +1,3 @@ -// Include simulation workflow -include { SIMULATION } from "${subworkflowDir}/simulation" - // Include modules here def prepare_dir = "${modulesDir}/prepare_data" include { PREPARE_MAE_DATA } from "${prepare_dir}/prepare_mae_data" diff --git a/subworkflows/simulation/main.nf b/subworkflows/simulation/main.nf deleted file mode 100644 index 4c1596d..0000000 --- a/subworkflows/simulation/main.nf +++ /dev/null @@ -1,193 +0,0 @@ -/* - This workflow is related for all simulations of data, to generate data with known ground-truth, relationships, - and properties like noise and signal to access method's robustness. - - Currently available strategies: - - Multivariate normal from CPLR - - InterSIM - - Mosim (TODO) -*/ - -// Include modules -def simulate_dir = "${modulesDir}/simulation" -include { SIMULATE_MVN_DATA } from "${simulate_dir}/simulate_mvn_data" -include { SIMULATE_INTERSIM } from "${simulate_dir}/simulate_intersim" -// Include functions -// include { createSimCombination } from "${modulesDir}/functions" - - -process COMPRESS_DATA_GZ { - debug false - def onSockeye = workflow.projectDir.toString().contains('/scratch') - // process metadata and configs - tag "${dataset_name}" - label 'process_low' - publishDir ( - path: "${params.outdir}/${task.process.tokenize(':')[-1].toLowerCase()}", - saveAS: { file }, - mode: 'copy', - overwrite: true - ) - - // These are required to handle to symlink from MAE data - stageInMode 'copy' - stageOutMode 'copy' - - input: - tuple val(dataset_name), path(mae_path), path(mu_path) - output: - path("${dataset_name}.tar.gz"), emit: ch_sim_datasets - - script: - """ - mkdir -p ${dataset_name} - cp -r ${mae_path} ${dataset_name} - cp ${mu_path} ${dataset_name} - - tar -czf ${dataset_name}.tar.gz ${dataset_name} - """ -} - - -workflow SIMULATION { - // WORKFLOW PARAMS - dataset_base_name = params.dataset_base_name - // TODO: output format need to be removed here and inside the - // scripts, just left here for compatibility - output_format = params.output_formats - num_obs = params.num_obs - y_name = params.y_name - // Strategy of simulation - skip_sim_MVN = params.skip_sim_MVN - skip_sim_intersim = params.skip_sim_intersim - - // Main code entrance - main: - // The simulation relies on a grid of parameters to simulate represented a groovy map - /* - For example: - grid = [ [ param1: ... , param2: ... , paramn: ...] , ... ] - - This way for each full param set in the grid, process could access those values by - its key names like param1, param2, and so on. - - NOTE: Some grid have less keys since certain simulation strategy do no support various - parameter, i.e. some only support changing number of observations and not number of variables - - big_grid = [ n: n1, p: p1, ... , s = 1] - small_grid = [ n: n1 , ... , s = 1] - */ - // Create simulation grids from default parameters - // mvn_sim_grid = createSimCombination(params) // Will determine if giving large grid or small grid - // Combination of parameters of simulation - - //ch_sim_params_comb.view() - - - /* - ======================================================================== - Setup the grid of exploitable params of simulation - ======================================================================== - - */ - - // Common params - ch_num_obs = Channel.fromList(num_obs) - ch_y_name = Channel.of(y_name) - - // ===================================================================== - // Make up the grid for InterSIM - // ===================================================================== - ch_intersim_num_noise_vars = Channel.fromList(params.intersim_num_noise_vars) - ch_intersim_effect = Channel.fromList(params.intersim_effect) - ch_intersim_noise = Channel.fromList(params.intersim_noise) - ch_intersim_sigma = Channel.fromList(params.intersim_sigma) - ch_intersim_corr = Channel.fromList(params.intersim_corr) - - - - // Assign to map for easy access later - Channel.of(dataset_base_name) - .combine(ch_num_obs) - .combine(ch_intersim_num_noise_vars) - .combine(ch_intersim_effect) - // TODO: might need better naming vs noise and sigma, since sigma here is more like covariance of omics - .combine(ch_intersim_noise) - .combine(ch_intersim_sigma) - .combine(ch_intersim_corr) - .map { m -> - [ dataset_name: "${m[0]}_strategy-intersim_n-${m[1]}_H-${m[2]}_effect-${m[3]}_e-${m[4]}_sigma-${m[5]}_corr-${m[6]}", - num_obs: m[1], num_noise_vars: m[2], effect: m[3], noise: m[4], sigma: m[5], corr: m[6] ] - } - .set { intersim_grid } - - // ===================================================================== - // Make up the grid for cplr mvn - // ===================================================================== - ch_mvn_num_predictors = Channel.fromList(params.num_predictors) - ch_mvn_block_num = Channel.fromList(params.block_num) - ch_mvn_latent_pred = Channel.fromList(params.latent_predictors) - ch_mvn_sigma = Channel.fromList(params.mvn_sigma) - ch_mvn_sy = Channel.fromList(params.sy) - ch_mvn_sp = Channel.fromList(params.sp) - ch_mvn_u_std = Channel.fromList(params.u_std) - ch_mvn_fct_str = Channel.fromList(params.fct_str) - - // Assign to map for easy access later - // - Channel.of(dataset_base_name) - .combine(ch_num_obs) - .combine(ch_mvn_num_predictors) - .combine(ch_mvn_sigma) - .combine(ch_mvn_block_num) - .combine(ch_mvn_latent_pred) - .combine(ch_mvn_sy) - .combine(ch_mvn_sp) - .combine(ch_mvn_u_std) - .combine(ch_mvn_fct_str) - .map { m -> - [ dataset_name: "${m[0]}_strategy-mvn_n-${m[1]}_p-${m[2]}_sigma-${m[3]}_j-${m[4]}_latp-${m[5]}_sy-${m[6]}_sp-${m[7]}_ustd-${m[8]}_fctstr-${m[9]}", - num_obs: m[1], num_predictors: m[2], sigma: m[3], block_num: m[4], latent_predictors: m[5], sy: m[6], sp: m[7], u_std: m[8], fct_str: m[9] - ] - } - .set { mvn_sim_grid } - - // Assign to map for easy access later - - /* - ======================================================================== - ACTUAL RUNNERS - ======================================================================== - */ - - // 1. When strategy is intersim pkg of simulation - ch_intersim = Channel.empty() - if (!skip_sim_intersim) { - // TODO: need to add their right params - SIMULATE_INTERSIM ( intersim_grid ) - ch_intersim = SIMULATE_INTERSIM.out.sim_mae.join(SIMULATE_INTERSIM.out.sim_mu) - } - // 2. When strategy is cplr's mvn simulation - ch_mvn = Channel.empty() - if (!skip_sim_MVN) { - // TODO: need to add their right params - SIMULATE_MVN_DATA ( mvn_sim_grid ) - ch_mvn = SIMULATE_MVN_DATA.out.sim_mae.join(SIMULATE_MVN_DATA.out.sim_mu) - } - - // Lastly mix these results - Channel.empty() - // TODO: Need to make sure they're mixable and have similar cardinality of arg in each channel - // prefarably have something to say what is the strategy, is simulated or not - // Mix each strategy together - .mix( ch_intersim ) - .mix( ch_mvn ) - .set { ch_sim_datasets } - - COMPRESS_DATA_GZ ( ch_sim_datasets ) - - // For each of these channels, we then compress them as tar.gz - // COMPRESS_DATA_GZ ( ch_sim_datasets ) - emit: - ch_sim_datasets = COMPRESS_DATA_GZ.out.ch_sim_datasets -} From 7944abd338f0febb807c4196eae96e032cf55f31 Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Mon, 31 Aug 2026 15:43:03 +0800 Subject: [PATCH 09/29] feat: sklearn to use pca50 or no reduction for more models --- modules/sklearn/train/main.nf | 15 ++-- .../resources/usr/bin/combine_mdata2df.py | 2 +- .../usr/bin/load_classifier_class.py | 74 ++++++++++++------- .../train/resources/usr/bin/sklearn_train.py | 45 +++++++++-- nextflow.config | 11 +-- 5 files changed, 99 insertions(+), 48 deletions(-) diff --git a/modules/sklearn/train/main.nf b/modules/sklearn/train/main.nf index e3b3c59..d83f00e 100755 --- a/modules/sklearn/train/main.nf +++ b/modules/sklearn/train/main.nf @@ -19,7 +19,7 @@ include { getPublishPath } from "${modulesDir}/functions" process SKLEARN_TRAIN { // Vars stuff // Identifier for each dataset and fold combination - tag "${dataset_name}-${fold_path.name}-${model_name}" + tag "${dataset_name}-${fold_path.name}-${model_name}-${reduction}" // TODO: rename to the actual image name used label 'sklearn' // Parse the output directory to migrate results to @@ -36,7 +36,9 @@ process SKLEARN_TRAIN { // TODO: Change arg name to mae_path or mu_path input: tuple val(dataset_name), path(data_path), path(fold_path) + each(reduction) each(model_name) + /* TODO: this part needs to be defined by user, mostly required is the template @@ -49,9 +51,9 @@ process SKLEARN_TRAIN { sample_script.py outputs my_model.pkl, then define path('my_model.pkl') in nextflow */ 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 + tuple val(dataset_name), val(fold_path.name), val("${model_name}-${reduction}"), path('*model*'), emit: model + tuple val(dataset_name), val(fold_path.name), val("${model_name}-${reduction}"), path('*test_data*'), emit: test_data + tuple val(dataset_name), val(fold_path.name), val("${model_name}-${reduction}"), path('*log*'), emit: log script: /* The script here should be found under method/resources/usr/bin/ , @@ -66,8 +68,9 @@ process SKLEARN_TRAIN { sklearn_train.py \ --fold_path=${fold_path} \ --label=${data_label} \ - --model_name=${model_name} > \ - ${data_label}-${model_name}-${getPublishPath(task.process).tokenize('/')[-1].toLowerCase()}.log + --model_name=${model_name} \ + --reduction=${reduction} > \ + ${data_label}-${model_name}-${reduction}-${getPublishPath(task.process).tokenize('/')[-1].toLowerCase()}.log echo ${data_label} > ${data_label} """ diff --git a/modules/sklearn/train/resources/usr/bin/combine_mdata2df.py b/modules/sklearn/train/resources/usr/bin/combine_mdata2df.py index 842f8ec..8e22d87 100644 --- a/modules/sklearn/train/resources/usr/bin/combine_mdata2df.py +++ b/modules/sklearn/train/resources/usr/bin/combine_mdata2df.py @@ -20,4 +20,4 @@ def combine_mdata2df(mdata, concat=True): 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 + return merged_df, mod_names \ No newline at end of file diff --git a/modules/sklearn/train/resources/usr/bin/load_classifier_class.py b/modules/sklearn/train/resources/usr/bin/load_classifier_class.py index 08f6fd4..b9d91eb 100644 --- a/modules/sklearn/train/resources/usr/bin/load_classifier_class.py +++ b/modules/sklearn/train/resources/usr/bin/load_classifier_class.py @@ -15,21 +15,21 @@ def load_classifier_class(model_name, random_state=42, probability=True): learning_rate = stats.uniform(0.01, 1.1) max_features = ["sqrt", "log2", 100, 500, 1000, None] # Dict to store relevant information of sklearn classifiers - model_info = { - # Logistic regression has built-in predict proba and coef - "Logit": { + + # 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": { + } + # 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 have risk of overfitting, hence require pruning of trees - "Decision_Tree": { + } + # 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": { @@ -38,9 +38,9 @@ def load_classifier_class(model_name, random_state=42, probability=True): "min_samples_leaf": min_sample_leaf, "max_leaf_nodes": [10, 100, 1000, None] } - }, - # RandomForest shuold in general work better than single decision tree - "Random_Forest": { + } + # 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": { @@ -48,19 +48,10 @@ def load_classifier_class(model_name, random_state=42, probability=True): "max_leaf_nodes": [10, 100, 1000, None], "min_samples_leaf": min_sample_leaf } - }, - # AdaBoost is a simpler boosting algorithm - "AdaBoost": { - "class_path": "sklearn.ensemble.AdaBoostClassifier", - "default_params": {"algorithm": "SAMME", "random_state": random_state}, - "params_dist": { - "n_estimators": n_estimators, - "learning_rate": learning_rate, - "algorithm": ['SAMME', 'SAMME.R'] - } - }, - # GradientBoost should tune large number of estimators with slow learning rate - "GradientBoost": { + } + + # 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": { @@ -69,7 +60,38 @@ def load_classifier_class(model_name, random_state=42, probability=True): "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 diff --git a/modules/sklearn/train/resources/usr/bin/sklearn_train.py b/modules/sklearn/train/resources/usr/bin/sklearn_train.py index a09e2ab..9d0215c 100755 --- a/modules/sklearn/train/resources/usr/bin/sklearn_train.py +++ b/modules/sklearn/train/resources/usr/bin/sklearn_train.py @@ -11,11 +11,14 @@ -h --help Show this message --fold_path=FOLD_PATH Directory containing one split directory of relevant input [default: empty] --label=LABEL Label of dataset and fold iteration [default: empty] + --reduction=REDUCTION Name of the reducer to use for dimensionality reduction [default: empty] --model_name=MOD Name of the classifier to run from sklearn [default: empty] """ from docopt import docopt -from sklearn.pipeline import make_pipeline +from sklearn.decomposition import PCA +from sklearn.compose import ColumnTransformer +from sklearn.pipeline import make_pipeline, Pipeline from sklearn.preprocessing import StandardScaler import mudata @@ -29,12 +32,28 @@ from combine_mdata2df import combine_mdata2df + +# Now also including a reducer model, with setting ncomps = 50 +def build_reducer(X_df, mod_names, n_comp=50, random_state=42): + transformers = [] + for m in mod_names: + prefix = f"{m}_" + cols = [c for c in X_df.columns if c.startswith(prefix)] + k = min(n_comp, len(cols), X_df.shape[0] - 1) + transformers.append((m, Pipeline([ + ("scale", StandardScaler()), + ("pca", PCA(n_components=k, random_state=random_state)), + ]), cols)) + return ColumnTransformer(transformers, remainder="drop") + + + # Train data is MuData # 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"): +def train(train_data, model_name, target_col="response", reduction="empty"): # Convert the mdata to merged dataframe column wise - merged_df = combine_mdata2df(train_data) + merged_df, mod_names = combine_mdata2df(train_data) # Transform the mudata into X and y for sklearn # X_df contains all count data from the views # y_df contains response @@ -46,25 +65,34 @@ def train(train_data, model_name, target_col="response"): clf = classifier_class(**params) print("\n", clf) # Tells to first apply standard scaling, then the model - clf = make_pipeline(StandardScaler(), clf) + if reduction != "empty": + # Build reducer + reducer = build_reducer(X_df, mod_names, n_comp=50) + # Create a pipeline with reducer and classifier + clf = make_pipeline(reducer, clf) + else: + clf = make_pipeline(StandardScaler(), clf) # Fitting model clf.fit(X_df, y_df["response"]) return clf # from upstream process -def main(fold_path, label, model_name): +def main(fold_path, label, model_name, reduction="empty"): # Load data first from fold_path train_data, test_data = load_tr_te(fold_path) # Train sklearn classifier available options are: # 'AdaBoost', 'Decision Tree', 'Gaussian Process', 'Linear SVM', 'Naive Bayes', # 'Nearest Neighbors', 'Neural Net', 'QDA', 'RBF SVM', 'Random Forest' # Case sensitive - model = train(train_data, model_name=model_name) + model = train(train_data, model_name=model_name, reduction=reduction) # Parse to more machine readable label model_label = model_name.lower().replace(" ", "_") # Parse label and choose output file to write - model_file = f"{label}-{model_label}-model.pkl" + if reduction != "empty": + model_file = f"{label}-{model_label}-{reduction}_model.pkl" + else: + model_file = f"{label}-{model_label}-model.pkl" # TODO: Decide to use pickle or joblib to write the trained model? joblib.dump(model, model_file) # Also write the test mudata to file so that it is passed to downstream @@ -81,6 +109,7 @@ def main(fold_path, label, model_name): main( fold_path = args["--fold_path"], label = args["--label"], - model_name = args["--model_name"] + model_name = args["--model_name"], + reduction = args["--reduction"] ) diff --git a/nextflow.config b/nextflow.config index a7d56db..e5d556c 100644 --- a/nextflow.config +++ b/nextflow.config @@ -39,6 +39,7 @@ params { // Options to run part of workflows // TODO: redundant parameters here + runSimple = false // Legacy param runInnerCV = false // Default to run inner cv in all possible methods runAllMethods = true // Run all multiomics integration methods, overrides runPython or runR runPython = true // Run methods implemented in Python only @@ -83,13 +84,9 @@ params { // SKLEARN // Available classifiers - /* - sklearn_classifier_names = [ - 'AdaBoost', 'Decision_Tree', 'GradientBoost', 'Linear_SVM', 'Logit', 'Random_Forest' - ] - */ - sklearn_classifier_names = ['Linear_SVM', 'Logit'] - + //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"] // For preprocessing data filter_low_var = "0" // This treats as false in R, one of "1" or "0" // Resource parameters From ec7e648e5d364cb82c9a24889c8397c1cdeb056d Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Mon, 31 Aug 2026 15:43:23 +0800 Subject: [PATCH 10/29] fix: update method options --- subworkflows/cross_validation/python/main.nf | 19 +++++-------------- subworkflows/cross_validation/r/main.nf | 5 ++--- subworkflows/methods/sklearn/main.nf | 7 ++++--- 3 files changed, 11 insertions(+), 20 deletions(-) diff --git a/subworkflows/cross_validation/python/main.nf b/subworkflows/cross_validation/python/main.nf index df6f105..2516d49 100644 --- a/subworkflows/cross_validation/python/main.nf +++ b/subworkflows/cross_validation/python/main.nf @@ -1,12 +1,11 @@ // Methods to include -include { INTEGRAO } from "${subworkflowDir}/methods/integrao" -include { SKLEARN } from "${subworkflowDir}/methods/sklearn" -include { MOGONET } from "${subworkflowDir}/methods/mogonet" -include { GOAT } from "${subworkflowDir}/methods/goat" +include { INTEGRAO } from "${subworkflowDir}/methods/integrao" +include { SKLEARN } from "${subworkflowDir}/methods/sklearn" +include { MOGONET } from "${subworkflowDir}/methods/mogonet" // 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 @@ -54,13 +53,6 @@ workflow CV_PYTHON { mogonet_results = MOGONET.out.csv_results } - // GOAT - goat_results = Channel.empty() - if (!skip_goat) { - GOAT ( mu_copy ) - goat_results = GOAT.out.csv_results - } - // ======================================================================== // Collect all result and mix it to merge it more Channel.empty() @@ -68,7 +60,6 @@ workflow CV_PYTHON { .mix( integrao_results ) .mix( sklearn_results ) .mix( mogonet_results ) - .mix( goat_results ) .map { it -> [ language_name, it[0], it[1] ] // Ch [R, method name, path of summary csv of method] } diff --git a/subworkflows/cross_validation/r/main.nf b/subworkflows/cross_validation/r/main.nf index 6b6001b..41a20f1 100644 --- a/subworkflows/cross_validation/r/main.nf +++ b/subworkflows/cross_validation/r/main.nf @@ -1,11 +1,10 @@ // Methods to include -include { CARET_MULTIMODAL } from "${subworkflowDir}/methods/caret_multimodal" -include { DEMO_LOGIT } from "${subworkflowDir}/methods/demo_logit" +include { CARET_MULTIMODAL } from "${subworkflowDir}/methods/caret_multimodal" +include { DEMO_LOGIT } from "${subworkflowDir}/methods/demo_logit" include { COOPERATIVE_LEARNING } from "${subworkflowDir}/methods/cooperative_learning" include { DIABLO } from "${subworkflowDir}/methods/diablo" include { MOFA } from "${subworkflowDir}/methods/mofa" include { RGCCA } from "${subworkflowDir}/methods/rgcca" -include { SGMR } from "${subworkflowDir}/methods/sgmr" // This module to collect results include { MERGE_RESULT_TABLE } from "${modulesDir}/merge_result_table" include { printBanner } from "${modulesDir}/functions" diff --git a/subworkflows/methods/sklearn/main.nf b/subworkflows/methods/sklearn/main.nf index 98e0eb1..71a778c 100755 --- a/subworkflows/methods/sklearn/main.nf +++ b/subworkflows/methods/sklearn/main.nf @@ -31,8 +31,9 @@ def saveMode = "method" workflow SKLEARN { // Classifier to train for sklearn - model_name = Channel.fromList(params.sklearn_classifier_names) - + model_name = Channel.fromList(params.sklearn_classifier_names) + // Reduction method (pca or empty) + reduction = Channel.value(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, @@ -72,7 +73,7 @@ workflow SKLEARN { */ // TODO: implement inner fold cross-validation in each method's train - SKLEARN_TRAIN ( train_input, model_name ) + SKLEARN_TRAIN ( train_input, reduction, model_name ) // Do some transformation to make a multiMap that has two branches for predict From e5ef195cc3e31aec1d9b42b111f9dc09278a03ed Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Thu, 3 Sep 2026 14:50:06 +0800 Subject: [PATCH 11/29] fix: include select feature for sklearn models --- .../usr/bin/load_classifier_class.py | 74 ++++++++++++------- .../resources/usr/bin/run_random_search_cv.py | 7 +- .../usr/bin/sklearn_select_features.py | 2 +- 3 files changed, 55 insertions(+), 28 deletions(-) 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 index 08f6fd4..b9d91eb 100644 --- a/modules/sklearn/select_feature/resources/usr/bin/load_classifier_class.py +++ b/modules/sklearn/select_feature/resources/usr/bin/load_classifier_class.py @@ -15,21 +15,21 @@ def load_classifier_class(model_name, random_state=42, probability=True): learning_rate = stats.uniform(0.01, 1.1) max_features = ["sqrt", "log2", 100, 500, 1000, None] # Dict to store relevant information of sklearn classifiers - model_info = { - # Logistic regression has built-in predict proba and coef - "Logit": { + + # 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": { + } + # 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 have risk of overfitting, hence require pruning of trees - "Decision_Tree": { + } + # 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": { @@ -38,9 +38,9 @@ def load_classifier_class(model_name, random_state=42, probability=True): "min_samples_leaf": min_sample_leaf, "max_leaf_nodes": [10, 100, 1000, None] } - }, - # RandomForest shuold in general work better than single decision tree - "Random_Forest": { + } + # 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": { @@ -48,19 +48,10 @@ def load_classifier_class(model_name, random_state=42, probability=True): "max_leaf_nodes": [10, 100, 1000, None], "min_samples_leaf": min_sample_leaf } - }, - # AdaBoost is a simpler boosting algorithm - "AdaBoost": { - "class_path": "sklearn.ensemble.AdaBoostClassifier", - "default_params": {"algorithm": "SAMME", "random_state": random_state}, - "params_dist": { - "n_estimators": n_estimators, - "learning_rate": learning_rate, - "algorithm": ['SAMME', 'SAMME.R'] - } - }, - # GradientBoost should tune large number of estimators with slow learning rate - "GradientBoost": { + } + + # 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": { @@ -69,7 +60,38 @@ def load_classifier_class(model_name, random_state=42, probability=True): "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 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 8e94adf..4ec6f2b 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 @@ -1,11 +1,16 @@ from sklearn.model_selection import RandomizedSearchCV +from sklearn.preprocessing import StandardScaler +from sklearn.pipeline import make_pipeline + # Run a randomized search cv to find optimal params def run_random_search_cv(clf_instance, X, Y, param_distributions, n_iter=10, random_state=42, n_jobs=-1): # 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) # Default to use all processors to speed up process clf_cv = RandomizedSearchCV( - clf_instance, param_distributions=param_distributions, + pipeline, param_distributions=param_distributions, n_iter=n_iter, random_state=random_state, n_jobs=n_jobs ) 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 97aae4b..62ee2ef 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 @@ -79,7 +79,7 @@ def main(mu_path, dataset_name, model_name, block_num=0, n_iter=10, random_state 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"]) + opt_clf.fit(X_df, y_df["response"].ravel()) # Extract the classifier from the pipeline classifier = opt_clf.steps[-1][1] # Adjust this based on your pipeline's step name # Then could either extract their weights or feature importance From 7c99113a31fa66df3698f7cdae701dbdbdf878fa Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Thu, 3 Sep 2026 14:50:39 +0800 Subject: [PATCH 12/29] fix: correct sklearn reduction to list param rather than singleton --- subworkflows/methods/sklearn/main.nf | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/subworkflows/methods/sklearn/main.nf b/subworkflows/methods/sklearn/main.nf index 71a778c..2bc2b3a 100755 --- a/subworkflows/methods/sklearn/main.nf +++ b/subworkflows/methods/sklearn/main.nf @@ -33,7 +33,7 @@ workflow SKLEARN { // Classifier to train for sklearn model_name = Channel.fromList(params.sklearn_classifier_names) // Reduction method (pca or empty) - reduction = Channel.value(params.sklearn_reduction) + 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, From 6f0ee70b81d32ae5478ff2e09d78db003ab047f3 Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Thu, 3 Sep 2026 14:51:07 +0800 Subject: [PATCH 13/29] fix: correct naming of parameter opt --- conf/real_data.config | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/conf/real_data.config b/conf/real_data.config index a57a822..6b881b7 100644 --- a/conf/real_data.config +++ b/conf/real_data.config @@ -25,7 +25,7 @@ params { skip_mogonet = false skip_mofa = false // Run feature selection - select_feature = true + selectFeature = true k_fold_number = 5 filter_low_var = "1" // This get casted as bool in both Python and R // number of component From 0ce67f6caa00cbbeed7d4ddb4c2447c04f6f656c Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Fri, 4 Sep 2026 17:53:42 +0800 Subject: [PATCH 14/29] fix: make sure sklearn cv search works even with make_pipeline --- .../resources/usr/bin/run_random_search_cv.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) 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..d3ca13b 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,27 @@ 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 ) # Could then fit this cv object of the X and Y of data search = clf_cv.fit(X, Y) - optimal_params = search.best_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 From b66e821fb4c6659c4b61541b0abccf5be6446e95 Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Fri, 11 Sep 2026 23:44:15 +0800 Subject: [PATCH 15/29] fix: use eipy container for sklearn as it contains xgboost, and use right classifiers --- nextflow.config | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/nextflow.config b/nextflow.config index e5d556c..1d9e33f 100644 --- a/nextflow.config +++ b/nextflow.config @@ -84,9 +84,12 @@ params { // SKLEARN // 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_classifier_names = ['Random_Forest', 'MLP', 'XGBoost'] + //sklearn_classifier_names = ["XGBoost"] + //sklearn_classifier_names = ['Logit', 'MLP', 'Decision_Tree'] sklearn_reduction = ["empty", "pca50"] + // If sklearn wants to do single input mode, default: false + sklearn_single_mode = true // For preprocessing data filter_low_var = "0" // This treats as false in R, one of "1" or "0" // Resource parameters @@ -209,7 +212,7 @@ 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 } From e190ae8e3abede7fe238b650315ee850ea7f4849 Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Sat, 12 Sep 2026 00:43:03 +0800 Subject: [PATCH 16/29] fix: moved common sklearn utils to bin/python_utils, and standardize funs --- bin/python_utils/combine_mdata2df.py | 46 +++++++ .../python_utils}/load_classifier_class.py | 19 ++- .../resources/usr/bin/combine_mdata2df.py | 23 ---- .../resources/usr/bin/sklearn_predict.py | 2 +- .../resources/usr/bin/combine_mdata2df.py | 23 ---- .../usr/bin/load_classifier_class.py | 113 ------------------ .../usr/bin/sklearn_select_features.py | 16 ++- .../resources/usr/bin/combine_mdata2df.py | 23 ---- .../train/resources/usr/bin/sklearn_train.py | 6 +- 9 files changed, 74 insertions(+), 197 deletions(-) create mode 100644 bin/python_utils/combine_mdata2df.py rename {modules/sklearn/train/resources/usr/bin => bin/python_utils}/load_classifier_class.py (91%) mode change 100644 => 100755 delete mode 100644 modules/sklearn/predict/resources/usr/bin/combine_mdata2df.py delete mode 100644 modules/sklearn/select_feature/resources/usr/bin/combine_mdata2df.py delete mode 100644 modules/sklearn/select_feature/resources/usr/bin/load_classifier_class.py delete mode 100644 modules/sklearn/train/resources/usr/bin/combine_mdata2df.py 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 91% rename from modules/sklearn/train/resources/usr/bin/load_classifier_class.py rename to bin/python_utils/load_classifier_class.py index b9d91eb..7bc36e7 --- 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 @@ -77,6 +82,16 @@ def load_classifier_class(model_name, random_state=42, probability=True): "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": 10, + "objective": "binary:logistic"}, + "params_dist": { "learning_rate": learning_rate, + "max_depth": stats.randint(2, 9) + } + } # ======================== # Lastly merge all together into a single dictionary @@ -91,7 +106,9 @@ def load_classifier_class(model_name, random_state=42, probability=True): "Gradient_Boost": gradient_boost_dict, # MLPClassifier is a neural network classifier "MLP": mlp_dict, - "Hist_Gradient_Boost": hist_gradient_boost_dict + "Hist_Gradient_Boost": hist_gradient_boost_dict, + # XGBoost goest here + "XGBoost": xgboost_dict } # Check if valid name of model was input 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/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/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/sklearn_select_features.py b/modules/sklearn/select_feature/resources/usr/bin/sklearn_select_features.py index 62ee2ef..8d80f25 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 @@ -60,26 +60,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) + optimal_params_dict = run_random_search_cv(clf_instance, X=X_df, Y=y, 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()) + opt_clf.fit(X_df, y) # Extract the classifier from the pipeline classifier = opt_clf.steps[-1][1] # Adjust this based on your pipeline's step name # Then could either extract their weights or feature importance 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) From e0956f7336d50b891716694de95a08c014da0b0b Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Sat, 12 Sep 2026 00:49:03 +0800 Subject: [PATCH 17/29] fix: refactor code for sklearn select feature --- .../resources/usr/bin/run_random_search_cv.py | 17 +++++++++-------- .../usr/bin/sklearn_select_features.py | 11 ++++------- 2 files changed, 13 insertions(+), 15 deletions(-) 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 d3ca13b..8a1eb16 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 @@ -24,11 +24,12 @@ def run_random_search_cv(clf_instance, X, Y, param_distributions, n_iter=10, ran ) # Could then fit this cv object of the X and Y of data search = clf_cv.fit(X, Y) - # 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 + # # 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 8d80f25..7ad7fcc 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 @@ -71,13 +71,10 @@ def main(mu_path, dataset_name, model_name, block_num=0, n_iter=10, random_state 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, 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) + search = run_random_search_cv(clf_instance, X=X_df, Y=y, param_distributions=param_dist, n_iter=n_iter, random_state=random_state) + opt_clf = search.best_estimator_ + print("Best parameters:", search.best_params_) + print("Best pipeline:", opt_clf) # Extract the classifier from the pipeline classifier = opt_clf.steps[-1][1] # Adjust this based on your pipeline's step name # Then could either extract their weights or feature importance From 44547e91f89d052a654c08b435f6c48a12304cb4 Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Wed, 16 Sep 2026 01:39:42 +0800 Subject: [PATCH 18/29] feat: update sklearn random search cv params as roc --- .../select_feature/resources/usr/bin/run_random_search_cv.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) 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 8a1eb16..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 @@ -20,7 +20,10 @@ def run_random_search_cv(clf_instance, X, Y, param_distributions, n_iter=10, ran clf_cv = RandomizedSearchCV( 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) From 36a59639d59570ce51113e9e1644bd18bce343ab Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Wed, 16 Sep 2026 02:35:21 +0800 Subject: [PATCH 19/29] fix: finalize models in sklearn --- bin/python_utils/load_classifier_class.py | 91 +++++++++-------------- 1 file changed, 36 insertions(+), 55 deletions(-) diff --git a/bin/python_utils/load_classifier_class.py b/bin/python_utils/load_classifier_class.py index 7bc36e7..ce5ad86 100755 --- a/bin/python_utils/load_classifier_class.py +++ b/bin/python_utils/load_classifier_class.py @@ -14,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 @@ -33,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], @@ -55,42 +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)} - } - # 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)} + "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) + } } + # Real xgboost xgboost_dict = { - "class_path": "xgboost.XGBClassifier", - "default_params": {"random_state": random_state , - "n_estimators": 10, - "objective": "binary:logistic"}, - "params_dist": { "learning_rate": learning_rate, - "max_depth": stats.randint(2, 9) - } + "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) + } } # ======================== @@ -98,16 +84,11 @@ 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 goest here + # XGBoost goes here "XGBoost": xgboost_dict } From 09fbe25a402901e2fecbc08ba664610f5793562c Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Wed, 16 Sep 2026 02:35:39 +0800 Subject: [PATCH 20/29] feat: add new option of sklearn single modality --- nextflow.config | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/nextflow.config b/nextflow.config index 1d9e33f..2640515 100644 --- a/nextflow.config +++ b/nextflow.config @@ -82,14 +82,15 @@ params { // DIABLO diablo_design_connection = ['full', 'null'] // connection for design matrix // SKLEARN + single_modality_mode = false // Available classifiers - sklearn_classifier_names = ['Random_Forest', 'MLP', 'XGBoost'] + //sklearn_classifier_names = ['Random_Forest', 'MLP', 'XGBoost'] //sklearn_classifier_names = ["XGBoost"] - //sklearn_classifier_names = ['Logit', 'MLP', 'Decision_Tree'] + sklearn_classifier_names = ['Logit'] sklearn_reduction = ["empty", "pca50"] // If sklearn wants to do single input mode, default: false - sklearn_single_mode = true + sklearn_single_mode = false // For preprocessing data filter_low_var = "0" // This treats as false in R, one of "1" or "0" // Resource parameters @@ -114,6 +115,7 @@ env { subworkflowDir = "$projectDir/subworkflows" modulesDir = "$projectDir/modules" configDir = "$projectDir/conf" + PYTHONPATH = "${projectDir}/bin/python_utils" } // Load base.config by default for all pipelines From eaefa010d386b7edc6221f7cef987e9123b7d2ba Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Wed, 16 Sep 2026 02:38:36 +0800 Subject: [PATCH 21/29] feat: include functionality to handle single modality running, limited to sklearn logistic regressions only --- subworkflows/cross_validation/main.nf | 19 ++++++-- subworkflows/cross_validation/python/main.nf | 10 +++- subworkflows/methods/sklearn/main.nf | 7 ++- subworkflows/prepare_data/main.nf | 48 ++++++++++++-------- workflows/messi_benchmark.nf | 21 ++++++++- 5 files changed, 78 insertions(+), 27 deletions(-) diff --git a/subworkflows/cross_validation/main.nf b/subworkflows/cross_validation/main.nf index a7fcb44..1da106f 100644 --- a/subworkflows/cross_validation/main.nf +++ b/subworkflows/cross_validation/main.nf @@ -73,7 +73,11 @@ workflow CROSS_VALIDATION { 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 +114,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,7 +136,7 @@ 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 { @@ -132,7 +145,7 @@ workflow CROSS_VALIDATION { } 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() diff --git a/subworkflows/cross_validation/python/main.nf b/subworkflows/cross_validation/python/main.nf index 2516d49..2eb910e 100644 --- a/subworkflows/cross_validation/python/main.nf +++ b/subworkflows/cross_validation/python/main.nf @@ -25,6 +25,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,7 +49,8 @@ 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 } diff --git a/subworkflows/methods/sklearn/main.nf b/subworkflows/methods/sklearn/main.nf index 2bc2b3a..af7ce59 100755 --- a/subworkflows/methods/sklearn/main.nf +++ b/subworkflows/methods/sklearn/main.nf @@ -33,7 +33,12 @@ workflow SKLEARN { // Classifier to train for sklearn 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, diff --git a/subworkflows/prepare_data/main.nf b/subworkflows/prepare_data/main.nf index 60c89b0..8694866 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,9 @@ 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 /* Workflow starts here */ // Workflow required input take: @@ -46,22 +48,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 } @@ -99,6 +86,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 +117,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..de0a107 100644 --- a/workflows/messi_benchmark.nf +++ b/workflows/messi_benchmark.nf @@ -60,11 +60,28 @@ 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 // From 43753b93b217def24114ec529c8b6ae2fb582779 Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Wed, 16 Sep 2026 03:07:11 +0800 Subject: [PATCH 22/29] feat: include missing scripts of splitting modalities --- modules/prepare_data/split_modality/main.nf | 36 +++++++++++++ .../resources/usr/bin/split_modalities.py | 52 +++++++++++++++++++ 2 files changed, 88 insertions(+) create mode 100644 modules/prepare_data/split_modality/main.nf create mode 100755 modules/prepare_data/split_modality/resources/usr/bin/split_modalities.py 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..7d03d66 --- /dev/null +++ b/modules/prepare_data/split_modality/resources/usr/bin/split_modalities.py @@ -0,0 +1,52 @@ +#!/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 + + +def split_modality(mdata, mod_names, dataset_name=""): + out_files = [] + 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}'") + sub = mudata.MuData({mod: mdata[mod].copy()}) + out_file = f"{dataset_name}-{mod}_processed.h5mu" + sub.write(out_file) + out_files.append(out_file) + print(f"[INFO] Wrote {out_file} ({sub.n_obs} obs x {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)) From 9245afc746ed97d7fbb6e979b41d2815db4794e3 Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Sat, 26 Sep 2026 18:14:26 +0800 Subject: [PATCH 23/29] fix: coerce scheme of rgcca and diablo to 'horst' --- .../resources/usr/bin/diablo_select_features.R | 2 +- modules/diablo/train/resources/usr/bin/run_diablo.R | 8 +++++--- modules/diablo/train/resources/usr/bin/tune_diablo.R | 3 ++- .../resources/usr/bin/rgcca_select_features.R | 5 +++++ modules/rgcca/train/resources/usr/bin/run_rgcca.R | 6 +++++- 5 files changed, 18 insertions(+), 6 deletions(-) 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..07d932f 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 @@ -97,7 +97,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 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/rgcca/select_feature/resources/usr/bin/rgcca_select_features.R b/modules/rgcca/select_feature/resources/usr/bin/rgcca_select_features.R index df7cf46..9175d6e 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 @@ -83,6 +83,10 @@ main <- function(mae_path, dataset_name, ncomp=2, design="full", prediction_mode # 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,6 +114,7 @@ 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 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") } From 1fbaaeb60a0ebd83bb1c4abb214ee6ba504db25f Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Sat, 26 Sep 2026 18:15:31 +0800 Subject: [PATCH 24/29] fix: make sure modality names are valid when exporting to filenames --- .../resources/usr/bin/split_modalities.py | 31 +++++++++++++++++-- 1 file changed, 28 insertions(+), 3 deletions(-) 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 index 7d03d66..16d2d38 100755 --- a/modules/prepare_data/split_modality/resources/usr/bin/split_modalities.py +++ b/modules/prepare_data/split_modality/resources/usr/bin/split_modalities.py @@ -19,20 +19,45 @@ 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}'") + 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()}) - out_file = f"{dataset_name}-{mod}_processed.h5mu" + + 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} ({sub.n_obs} obs x {mdata[mod].n_vars} vars)") + + print( + f"[INFO] Wrote {out_file} " + f"({sub.n_obs} obs × {mdata[mod].n_vars} vars)" + ) + return out_files From 06ede3afd90c5065b4539bd37271cc7713b47db4 Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Sat, 26 Sep 2026 18:16:07 +0800 Subject: [PATCH 25/29] feat: add new boolean to run classification or survival --- nextflow.config | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/nextflow.config b/nextflow.config index 2640515..837fda5 100644 --- a/nextflow.config +++ b/nextflow.config @@ -72,9 +72,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 @@ -84,11 +90,8 @@ params { // SKLEARN single_modality_mode = false // Available classifiers - - //sklearn_classifier_names = ['Random_Forest', 'MLP', 'XGBoost'] - //sklearn_classifier_names = ["XGBoost"] - sklearn_classifier_names = ['Logit'] - sklearn_reduction = ["empty", "pca50"] + sklearn_classifier_names = ['Logit', 'Random_Forest', 'MLP', 'XGBoost'] + sklearn_reduction = ["empty", "pca50"] // Preprocessing done to data // If sklearn wants to do single input mode, default: false sklearn_single_mode = false // For preprocessing data From 8f4a27e986bd34d079f498f6fd09f20d6d2347fc Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Mon, 28 Sep 2026 17:41:40 +0800 Subject: [PATCH 26/29] fix: remove old skipper of methods --- nextflow.config | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/nextflow.config b/nextflow.config index 837fda5..94e9902 100644 --- a/nextflow.config +++ b/nextflow.config @@ -56,8 +56,9 @@ params { skip_rgcca = false skip_mogonet = false skip_mofa = false - skip_goat = true - skip_sgmr = true + //skip_goat = true + //skip_sgmr = true + skip_sksurv = true skip_sklearn = true /*Parameters for methods*/ From b598a38e2576270079273ae5d8a50fb8c090d69e Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Mon, 28 Sep 2026 17:42:39 +0800 Subject: [PATCH 27/29] feat: add sksurv modules (not tested, need to add container def too --- modules/sksurv/predict/main.nf | 43 +++++ .../usr/bin/combine_mdata2df_survival.py | 74 ++++++++ .../usr/bin/generate_result_table_survival.py | 38 ++++ .../usr/bin/load_survival_model_class.py | 165 ++++++++++++++++++ .../resources/usr/bin/manual_set_seed.py | 25 +++ .../resources/usr/bin/sksurv_predict.py | 76 ++++++++ modules/sksurv/preprocess/main.nf | 46 +++++ .../resources/usr/bin/load_test_splits.py | 27 +++ .../resources/usr/bin/sksurv_preprocess.py | 53 ++++++ .../resources/usr/bin/tr_te_split_mdata.py | 12 ++ modules/sksurv/train/main.nf | 45 +++++ .../usr/bin/combine_mdata2df_survival.py | 74 ++++++++ .../usr/bin/load_survival_model_class.py | 165 ++++++++++++++++++ .../train/resources/usr/bin/load_tr_te.py | 15 ++ .../resources/usr/bin/manual_set_seed.py | 25 +++ .../train/resources/usr/bin/sksurv_train.py | 59 +++++++ subworkflows/cross_validation/python/main.nf | 15 +- subworkflows/methods/sklearn/main.nf | 4 +- subworkflows/methods/sksurv/main.nf | 82 +++++++++ 19 files changed, 1039 insertions(+), 4 deletions(-) create mode 100644 modules/sksurv/predict/main.nf create mode 100644 modules/sksurv/predict/resources/usr/bin/combine_mdata2df_survival.py create mode 100644 modules/sksurv/predict/resources/usr/bin/generate_result_table_survival.py create mode 100644 modules/sksurv/predict/resources/usr/bin/load_survival_model_class.py create mode 100644 modules/sksurv/predict/resources/usr/bin/manual_set_seed.py create mode 100755 modules/sksurv/predict/resources/usr/bin/sksurv_predict.py create mode 100644 modules/sksurv/preprocess/main.nf create mode 100644 modules/sksurv/preprocess/resources/usr/bin/load_test_splits.py create mode 100755 modules/sksurv/preprocess/resources/usr/bin/sksurv_preprocess.py create mode 100644 modules/sksurv/preprocess/resources/usr/bin/tr_te_split_mdata.py create mode 100644 modules/sksurv/train/main.nf create mode 100644 modules/sksurv/train/resources/usr/bin/combine_mdata2df_survival.py create mode 100644 modules/sksurv/train/resources/usr/bin/load_survival_model_class.py create mode 100644 modules/sksurv/train/resources/usr/bin/load_tr_te.py create mode 100644 modules/sksurv/train/resources/usr/bin/manual_set_seed.py create mode 100755 modules/sksurv/train/resources/usr/bin/sksurv_train.py create mode 100644 subworkflows/methods/sksurv/main.nf diff --git a/modules/sksurv/predict/main.nf b/modules/sksurv/predict/main.nf new file mode 100644 index 0000000..8bc4c4e --- /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 'sklearn' + + 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..9ad4c7c --- /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 'sklearn' + 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..9c3933a --- /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 'sklearn' + 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/subworkflows/cross_validation/python/main.nf b/subworkflows/cross_validation/python/main.nf index 2eb910e..e6d726b 100644 --- a/subworkflows/cross_validation/python/main.nf +++ b/subworkflows/cross_validation/python/main.nf @@ -18,6 +18,7 @@ workflow CV_PYTHON { skip_sklearn = params.skip_sklearn // boolean: true/false skip_mogonet = params.skip_mogonet // boolean: true/false skip_goat = params.skip_goat // boolean: true/false + skip_sksurv = params.skip_sksurv // boolean: true/false // Method specific parameters he_base_dim = params.he_base_dim // Inputs of workflow @@ -53,6 +54,13 @@ workflow CV_PYTHON { 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() @@ -65,9 +73,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] } diff --git a/subworkflows/methods/sklearn/main.nf b/subworkflows/methods/sklearn/main.nf index af7ce59..b0595c1 100755 --- a/subworkflows/methods/sklearn/main.nf +++ b/subworkflows/methods/sklearn/main.nf @@ -31,7 +31,9 @@ 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) // Set reduction to empty if single_mode is true if ( params.single_modality_mode ) { diff --git a/subworkflows/methods/sksurv/main.nf b/subworkflows/methods/sksurv/main.nf new file mode 100644 index 0000000..3fb6a44 --- /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 ) + // ===================================================================== + + emit: + csv_results = MERGE_RESULT_TABLE.out.csv_results +} From 6fbc86ff1e33637fd4ed76ba85983dad37701708 Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Mon, 28 Sep 2026 17:43:10 +0800 Subject: [PATCH 28/29] feat: mofa try survival, need to run later --- .../resources/usr/bin/preprocess_mofa.R | 9 +- modules/mofa/train/main.nf | 6 +- .../resources/usr/bin/run_mofa_survival.R | 129 ++++++++++++++++++ 3 files changed, 141 insertions(+), 3 deletions(-) create mode 100644 modules/mofa/train/resources/usr/bin/run_mofa_survival.R diff --git a/modules/mofa/preprocess/resources/usr/bin/preprocess_mofa.R b/modules/mofa/preprocess/resources/usr/bin/preprocess_mofa.R index 8b44516..5271a61 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: @@ -138,6 +138,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 diff --git a/modules/mofa/train/main.nf b/modules/mofa/train/main.nf index 603ddfd..d331482 100644 --- a/modules/mofa/train/main.nf +++ b/modules/mofa/train/main.nf @@ -31,9 +31,11 @@ process MOFA_TRAIN { tuple val(dataset_name), val(fold_path.name), path('*log*'), emit: log script: def data_label = "${dataset_name}-${fold_path.name}" + // Dynamically determine script based on classification or survival task + def script_names = params.outcome_type == "classification" ? "run_mofa.R" : "run_mofa_survival.R" if (run_inner_cv) """ - run_mofa.R \ + ${script_names} \ --mae_path=${mae_path} \ --label=${data_label} \ --fold_path=${fold_path} \ @@ -43,7 +45,7 @@ process MOFA_TRAIN { """ else """ - run_mofa.R \ + ${script_names} \ --mae_path=${mae_path} \ --label=${data_label} \ --fold_path=${fold_path} > \ 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)) +} From bb719ea7f0ae6f0c83ea0f7efea80ee0750c8c89 Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Mon, 28 Sep 2026 17:44:00 +0800 Subject: [PATCH 29/29] feat: this THU specific launcher local script wrapper, del later --- launcher_local.sh | 76 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 76 insertions(+) create mode 100755 launcher_local.sh diff --git a/launcher_local.sh b/launcher_local.sh new file mode 100755 index 0000000..21abfb4 --- /dev/null +++ b/launcher_local.sh @@ -0,0 +1,76 @@ +#!/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 +# ============================================================================ + +# 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_TEMP="/data1/tliang19/tools/nextflow/nf-tmp" +#SELECT_FEAT='true' # Or use 'false' +SELECT_FEAT='false' + +#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" + +F="false" +T="true" + +FILTER="1" # True in nextflow +#FILTER="0" + +#K_FOLD_NUMBER=2 +K_FOLD_NUMBER=5 +#SK_MODS="Logit,XGBoost" +#SK_MODS="Logit,MLP" + +#MAX_MEM="4G" +MAX_MEM="6G" + +# Specify number of CPUs and memory +nextflow run main.nf \ + -profile standard,docker,test \ + --single_modality_mode $T \ + --samplesheet ${CSV_FILE} \ + --filter_low_var $FILTER \ + --outdir results \ + --max_memory $MAX_MEM \ + --skip_rgcca $F \ + --skip_mofa $T \ + --skip_sklearn $T \ + --skip_caret_multimodal $T \ + --skip_diablo $F \ + --skip_cplr $T \ + --skip_mogonet $T \ + --skip_integrao $T \ + --pipeline_dir ./ \ + --k_fold_number $K_FOLD_NUMBER \ + --selectFeature $F \ + --publish_relevant $T \ + -resume