From e308292a34225eee152bde8b35148f65496366b9 Mon Sep 17 00:00:00 2001 From: Tony Liang Date: Thu, 13 Aug 2026 21:56:46 +0800 Subject: [PATCH 01/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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))