TCGA-LGG (Lower Grade Glioma) - Binary Survival Prediction

Cohort Summary:

Stage	Methylation	mRNA	miRNA	Clinical
Raw (features x samples)	20,115 x 685	18,328 x 701	548 x 531	14 x 1,110
Final aligned (samples x features)	511 x 20,114	511 x 18,328	511 x 548	511 x 15
After feature selection	511 x 400	511 x 200	511 x 100	511 x 15

Target Definition: Binary vital status, Alive (n=386) vs. Deceased (n=125).

Data Source: FireBrowse LGG

import pandas as pd
from pathlib import Path
root = Path("/home/vicente/Github/BioNeuralNet/LGG")

mirna_raw = pd.read_csv(root/"LGG.miRseq_RPKM_log2.txt", sep="\t",index_col=0,low_memory=False)                            
rna_raw = pd.read_csv(root / "LGG.uncv2.mRNAseq_RSEM_normalized_log2.txt", sep="\t",index_col=0,low_memory=False)
meth_raw = pd.read_csv(root/"LGG.meth.by_mean.data.txt", sep='\t',index_col=0,low_memory=False)
clinical_raw = pd.read_csv(root / "LGG.clin.merged.picked.txt",sep="\t", index_col=0, low_memory=False)

# display shapes and first few rows-columns of each file
display(mirna_raw.iloc[:3,:5])
display(mirna_raw.shape)

display(rna_raw.iloc[:3,:5])
display(rna_raw.shape)

display(meth_raw.iloc[:3,:5])
display(meth_raw.shape)

display(clinical_raw.iloc[:3,:5])
display(clinical_raw.shape)

Data Processing Summary

1. Transpose Data: All raw data (miRNA, RNA, etc.) is flipped so rows represent patients and columns represent features.

2. Standardize Patient IDs: Patient IDs in all tables are cleaned to the 12-character TCGA format (e.g., TCGA-AB-1234) for matching.

3. Handle Duplicates: Duplicate patient rows are averaged in the omics data. The first entry is kept for duplicate patients in the clinical data.

4. Find Common Patients: The script identifies the list of patients that exist in all datasets.

5. Sparse Feature Filter: Remove features where missingness exceeds the imputation limit (20%).

6. Final Data Aligment: All data tables are filtered down to only this common list of patients.

7. Impute Missing Values: Remaining missing data (NaNs) are estimated and filled using mean imputation.

8. Extract Target: The vital_status column is pulled from the processed clinical data to be used as the prediction target (y-variable).

from bioneuralnet.utils import m_transform, impute_simple
from bioneuralnet.utils import data_stats, sparse_filter

mirna = mirna_raw.T
rna = rna_raw.T
meth = meth_raw.T
clinical = clinical_raw.T

print(f"miRNA (samples, features): {mirna.shape}")
print(f"RNA (samples, features): {rna.shape}")
print(f"Methylation (samples, features): {meth.shape}")
print(f"Clinical (samples, features): {clinical.shape}")

def trim_barcode(idx):
    return idx.to_series().str.slice(0, 12)

# Standardized patient IDs across all files
meth.index = trim_barcode(meth.index)
rna.index = trim_barcode(rna.index)
mirna.index = trim_barcode(mirna.index)
clinical.index = clinical.index.str.upper()
clinical.index.name = "Patient_ID"

# Convert all data to numeric, coercing errors to NaN
meth = meth.apply(pd.to_numeric, errors='coerce')
rna = rna.apply(pd.to_numeric, errors='coerce')
mirna = mirna.apply(pd.to_numeric, errors='coerce')

# For any duplicate columns in the omics data -> average their values
meth = meth.groupby(meth.index).mean()
rna = rna.groupby(rna.index).mean()
mirna = mirna.groupby(mirna.index).mean()

# Drop unwanted columns
clinical.drop(columns=["Composite Element REF"], errors="ignore", inplace=True)
meth.drop(columns=["Composite Element REF"], errors="ignore", inplace=True)

# For any duplicate rows (patients) in the clinical data -> keep the first occurrence
clinical = clinical[~clinical.index.duplicated(keep='first')]

# Convert beta values to M-values using bioneuralnet utility with small epsilon to avoid log(0)
meth_m = m_transform(meth, eps=1e-7) 

print(f"RNA shape: {rna.shape}")
print(f"miRNA shape: {mirna.shape}")
print(f"Clinical shape: {clinical.shape}")

# Standardize column names
for df in [meth_m, rna, mirna]:
    df.columns = df.columns.str.replace(r"\?", "unknown_", regex=True)
    df.columns = df.columns.str.replace(r"\|", "_", regex=True)
    df.columns = df.columns.str.replace("-", "_", regex=False)
    df.columns = df.columns.str.replace(r"_+", "_", regex=True)
    df.columns = df.columns.str.strip("_")

# Patients common across all data files
common_patients = sorted(list(set(meth_m.index) & set(rna.index) & set(mirna.index) & set(clinical.index)))

# Subset to only common patients
X_meth = meth_m.loc[common_patients].copy()
X_rna = rna.loc[common_patients].copy()
X_mirna = mirna.loc[common_patients].copy()
clinical = clinical.loc[common_patients].copy()

# Calling data_stats on subset to prevent issues downstream
data_stats(X_mirna, "miRNA")
data_stats(X_rna, "RNA")
data_stats(X_meth, "Methylation")

# Optinal: this will drop sparse columns up to a user defined Threshold
# For this example we kept all see TCGA-BRCA where we use this sparse_filter
# X_mirna = sparse_filter(X_mirna, missing_fraction=0.20)
# X_rna = sparse_filter(X_rna, missing_fraction=0.20)
# X_meth = sparse_filter(X_meth, missing_fraction=0.20)

# Secondary aligment after sparse filter
final_patients = sorted(
    set(X_meth.index) & 
    set(X_rna.index) & 
    set(X_mirna.index) &
    set(clinical.index)
)

# subsetting based on secondary alignment
X_meth = X_meth.loc[final_patients]
X_rna = X_rna.loc[final_patients]
X_mirna = X_mirna.loc[final_patients]
clinical = clinical.loc[final_patients]

# impute the remaining missing values
# default impute method is set to: "mean"
X_meth = impute_simple(X_meth)
X_rna = impute_simple(X_rna)
X_mirna = impute_simple(X_mirna)

# target labels from clinical data
Y_label = clinical['vital_status']
Y_label = Y_label.to_frame(name="target")

# Inspect the first 3 rows and 5 colums.

display(X_meth.iloc[:3,:5])
display(X_meth.shape)

display(X_rna.iloc[:3,:5])
display(X_rna.shape)

display(X_mirna.iloc[:3,:5])
display(X_mirna.shape)

display(clinical.iloc[:3,:5])
display(clinical.shape)

display(Y_label.value_counts())

Feature Selection

Unsupervised feature selection was performed across all three omic modalities using Laplacian Score filtering.

from bioneuralnet.utils import laplacian_score

methylation_lap = laplacian_score(X_meth, n_keep=400)
rna_lap = laplacian_score(X_rna, n_keep=200)
mirna_lap = laplacian_score(X_mirna, n_keep=100)

import numpy as np
from bioneuralnet.utils import preprocess_clinical

columns_to_drop = [
    'Hybridization REF',
    'tumor_tissue_site',
    'date_of_initial_pathologic_diagnosis',
    'vital_status',
    'days_to_death',
    'days_to_last_followup',
    'days_to_last_known_alive'
]

continuous_cols = ['years_to_birth', 'karnofsky_performance_score']

clinical_numeric = preprocess_clinical(
    X=clinical,
    scale=True,
    impute=True,
    drop_columns=columns_to_drop,
    continuous_columns=continuous_cols)


print(f"\nNaN count: {clinical_numeric.isna().sum().sum()}")
print(f"\nInf count: {np.isinf(clinical_numeric.values).sum()}")
print(f"\nColumn dtypes:\n{clinical_numeric.dtypes}")
print(clinical_numeric.describe())

Data Availability

To facilitate rapid experimentation and reproduction of our results, the fully processed and feature-selected dataset used in this example has been made available directly within the package.

Users can load this dataset, bypassing all preceding data acquisition, preprocessing, and feature selection steps.

import pandas as pd
from bioneuralnet.datasets import DatasetLoader

tgca_lgg = DatasetLoader("lgg")
display(tgca_lgg.shape)

# The dataset is returned as a dictionary. We extract each file independently based on the name( Key).
X_methylation = tgca_lgg["methylation"]
X_rna = tgca_lgg["rna"]
X_mirna = tgca_lgg["mirna"]
clinical_numeric = tgca_lgg["clinical"]
Y_labels = tgca_lgg["target"]

display(X_methylation.iloc[:3,:5])
display(X_rna.iloc[:3,:5])
display(X_mirna.iloc[:3,:5])
display(clinical_numeric.iloc[:3,:5])
display(Y_labels.iloc[:3,:5])

Reproducibility and Seeding

This utility function propagates the seed to all sources of randomness, including random, numpy, and torch.

from bioneuralnet.utils import set_seed

SEED = 8183
set_seed(SEED)

from bioneuralnet.network import correlation_network, similarity_network

omics_lgg = pd.concat([X_methylation, X_rna, X_mirna], axis=1)

similarity_22 = similarity_network(omics_lgg, k=22)
spearman_12 = correlation_network(omics_lgg, method="spearman", k=12)

from bioneuralnet.network import NetworkAnalyzer

sim_22_analysis = NetworkAnalyzer(similarity_22, source_omics=[X_methylation, X_rna, X_mirna])
sim_22_analysis.basic_statistics(0.0001)
sim_22_analysis.hub_analysis(0.0001)
sim_22_analysis.find_strongest_edges(5)

spear_12_analysis = NetworkAnalyzer(spearman_12, source_omics=[X_methylation, X_rna, X_mirna])
spear_12_analysis.basic_statistics(0.0001)
spear_12_analysis.hub_analysis(0.0001)
spear_12_analysis.find_strongest_edges(5)

# A targeted subset of clinical variables is selected. This can be modified as needed for testing/experiments

clinical_numeric = clinical_numeric[[
    'years_to_birth',
    'radiation_therapy_yes',
    'histological_type_oligoastrocytoma',
    'histological_type_oligodendroglioma',
    'race_asian',
    'ethnicity_not hispanic or latino',
    'race_black or african american'
]].copy()

Classification using DPMON: Training and Evaluation

from bioneuralnet.downstream_task import DPMON

output_dir_base = Path("/home/vicente/Github/BioNeuralNet/dpmon_gcn_results/lgg")

current_output_dir = output_dir_base / "similarity_22"
current_output_dir.mkdir(parents=True, exist_ok=True)

dpmon_params_base = {
    "adjacency_matrix": similarity_22,
    "omics_list": omics_lgg,
    "phenotype_data": Y_labels,
    "phenotype_col": "target",
    "clinical_data": clinical_numeric,
    "tune_trials" : 20,
    "model": 'GCN',
    "tune": True, 
    "cv": True,   
    "n_folds": 5,
    "repeat_num": 5,
    "gpu": True,
    "seed": SEED,
    "output_dir": current_output_dir
}

gcn_dpmon = DPMON(**dpmon_params_base)
gcn_predictions, gcn_metrics, gcn_embeddings = gcn_dpmon.run()

from bioneuralnet.downstream_task import DPMON

output_dir_base = Path("/home/vicente/Github/BioNeuralNet/dpmon_gat_results/lgg")

current_output_dir = output_dir_base / "spearman_12"
current_output_dir.mkdir(parents=True, exist_ok=True)

dpmon_params_base = {
    "adjacency_matrix": spearman_12,
    "omics_list": omics_lgg,
    "phenotype_data": Y_labels,
    "phenotype_col": "target",
    "clinical_data": clinical_numeric,
    "tune_trials" : 20,
    "model": 'GAT',
    "tune": True, 
    "cv": True,   
    "n_folds": 5,
    "repeat_num": 5,
    "gpu": True,
    "seed": SEED,
    "output_dir": current_output_dir
}

gat_dpmon = DPMON(**dpmon_params_base)
gat_predictions, gat_metrics, gat_embeddings = gat_dpmon.run()

Latent Space Visulization

Project feature embeddings from the DPMON module into a 2D latent space. This allows users to inspect how different omics layers cluster and evaluate the distribution of learned representations within the model.

# for a more interpretable embedding plot we will map the node labels to omics types
node_labels = []
GAT_embeddings_array = gcn_embeddings.values
feature_names = gcn_embeddings.index

# using the sets from the orginal to serve as a loopup table
dna_feats = set(X_methylation.columns)
rna_feats = set(X_rna.columns)
mirna_feats = set(X_mirna.columns)

for feat in feature_names:
    if feat in dna_feats:
        # 1 for DNA methylation
        node_labels.append(1)
    elif feat in rna_feats:
        # 2 for RNA 
        node_labels.append(2)
    elif feat in mirna_feats:
        # and 3 for miRNA 
        node_labels.append(3)
    else:
        node_labels.append(0)

# our plotting function needs an array to squeeze the embedings into 2Dspace
node_labels = np.array(node_labels)

from bioneuralnet.metrics import plot_embeddings

gcn_embeddings_array = gcn_embeddings.values
feature_names = gcn_embeddings.index

# feature(omics) embeddings colored by omics type.
print(f"Plotting Feature Embeddings for GCN graph: GCN paired w/ Similarity 22")
plot_embeddings(gcn_embeddings_array, node_labels, legend_labels=["DNA_Methylation", "RNA", "miRNA"])

# Similar workflow for GAT embeddings.

gat_embeddings_array = gat_embeddings.values
feature_names = gat_embeddings_array.index

print(f"Plotting Feature Embeddings for GCN graph: GAT paired w/ Spearman 12")
plot_embeddings(gat_embeddings_array, node_labels, legend_labels=["DNA_Methylation", "RNA", "miRNA"])

from bioneuralnet.downstream_task import DPMON

output_dir_base = Path("/home/vicente/Github/BioNeuralNet/dpmon_gat_multiple_results/lgg")
all_results_gat = []

# To test several graphs sequentially. We can loop over the graphs dynamically like this:

comparison_runs = [
    {"name": "similarity_22", "graph": similarity_22, "omics": omics_lgg},
    {"name": "spearman_12", "graph": spearman_12, "omics": omics_lgg},
]

for run_config in comparison_runs:
    graph_name = run_config["name"]
    A_full = run_config["graph"]
    omics = run_config["omics"]
    
    current_output_dir = output_dir_base / graph_name
    current_output_dir.mkdir(parents=True, exist_ok=True)

    dpmon_params_base = {
        "adjacency_matrix": A_full,
        "omics_list": omics,
        "phenotype_data": Y_labels,
        "phenotype_col": "target",
        "clinical_data": clinical_numeric,
        "tune_trials" : 20,
        "model": 'GAT',
        "tune": True, 
        "cv": True,   
        "n_folds": 5,
        "repeat_num":5,
        "gpu": True,
        "seed": SEED,
        "output_dir": current_output_dir
    }
    
    dpmon_tunned = DPMON(**dpmon_params_base)
    predictions_df, metrics, embeddings = dpmon_tunned.run()

    all_results_gat.append({
        "graph": graph_name,
        "predictions": predictions_df,
        "metrics": metrics,
        "embeddings": embeddings
    })

# Once all the experiments are done, we can review the results for each graph configuration:
for res in all_results_gat:
    graph_name = res["graph"]
    graph_metrics = res["metrics"]
    
    acc_row = graph_metrics.loc[graph_metrics['Metric'] == 'Accuracy'].iloc[0]
    f1_macro_row = graph_metrics.loc[graph_metrics['Metric'] == 'F1 Macro'].iloc[0]
    f1_weighted_row = graph_metrics.loc[graph_metrics['Metric'] == 'F1 Weighted'].iloc[0]
    recall_row = graph_metrics.loc[graph_metrics['Metric'] == 'Recall'].iloc[0]
    precission_row = graph_metrics.loc[graph_metrics['Metric'] == 'Precision'].iloc[0]
    auc_row = graph_metrics.loc[graph_metrics['Metric'] == 'AUC'].iloc[0]
    aupr_row = graph_metrics.loc[graph_metrics['Metric'] == 'AUPR'].iloc[0]
    
    acc_avg, acc_std = acc_row['Average'], acc_row['StdDev']
    f1_macro_avg, f1_macro_std = f1_macro_row['Average'], f1_macro_row['StdDev']
    f1_weighted_avg, f1_weighted_std = f1_weighted_row['Average'], f1_weighted_row['StdDev']
    recall_avg, recall_std = recall_row['Average'], recall_row['StdDev']
    precision_avg, precision_std = precission_row['Average'], precission_row['StdDev']
    auc_avg, auc_std = auc_row['Average'], auc_row['StdDev']
    aupr_avg, aupr_std = aupr_row['Average'], aupr_row['StdDev']

    print(f"\nResults for: {graph_name}")
    print(f"Accuracy (Avg +/- Std): {acc_avg:.4f} +/- {acc_std:.4f}")
    print(f"F1 Macro (Avg +/- Std): {f1_macro_avg:.4f}  +/- {f1_macro_std:.4f}")
    print(f"F1 Weighted (Avg +/- Std): {f1_weighted_avg:.4f} +/- {f1_weighted_std:.4f}")
    print(f"Recall: {recall_avg:.4f} +/- {recall_std:.4f}")
    print(f"Precision: {precision_avg:.4f} +/- {precision_std:.4f}")
    print(f"AUC: {auc_avg:.4f} +/- {auc_std:.4f}")
    print(f"AUPR: {aupr_avg:.4f} +/- {aupr_std:.4f}")

from sklearn.model_selection import StratifiedKFold, ParameterSampler, RepeatedStratifiedKFold
from sklearn.metrics import accuracy_score, f1_score, recall_score, roc_auc_score, average_precision_score, precision_score
from sklearn.preprocessing import StandardScaler, label_binarize
from sklearn.linear_model import LogisticRegression
from sklearn.neural_network import MLPClassifier
from sklearn.pipeline import Pipeline
from sklearn.base import clone
from scipy.stats import loguniform, randint
from sklearn.ensemble import RandomForestClassifier
from sklearn.tree import DecisionTreeClassifier
from xgboost import XGBClassifier
from sklearn.svm import SVC

SEED = 8183
y = Y_labels["target"]
X = pd.concat([omics_lgg,clinical_numeric],axis=1)

print(f"Successfully created X vector with shape: {X.shape}")
print(f"Successfully created y vector with shape: {y.shape}")

pipe_lr = Pipeline([
    ('scaler', StandardScaler()),
    ('model', LogisticRegression(solver='lbfgs', max_iter=1000, penalty=None, random_state=SEED))
])

pipe_mlp = Pipeline([
    ('scaler', StandardScaler()),
    ('model', MLPClassifier(max_iter=500, early_stopping=True, n_iter_no_change=10, random_state=SEED))
])

pipe_xgb = Pipeline([
    ('scaler', StandardScaler()),
    ('model', XGBClassifier(eval_metric='logloss', tree_method='hist', max_bin=128, random_state=SEED))
])
pipe_rf = Pipeline([
    ('scaler', StandardScaler()),
    ('model', RandomForestClassifier(random_state=SEED))
])

pipe_svm = Pipeline([
    ('scaler', StandardScaler()),
    ('model', SVC(probability=True, random_state=SEED))
])

pipe_dt = Pipeline([
    ('scaler', StandardScaler()),
    ('model', DecisionTreeClassifier(random_state=SEED))
])

params_lr = {'model__penalty': ['l2'], 'model__C': loguniform(1e-4, 1e2)}

params_mlp = {
    'model__hidden_layer_sizes': [(100,), (100, 50), (50, 50)],
    'model__activation': ['relu', 'tanh'],
    'model__alpha': loguniform(1e-5, 1e-1),
    'model__learning_rate_init': loguniform(1e-4, 1e-2)
}
params_xgb = {
    'model__n_estimators': randint(50, 200),
    'model__learning_rate': loguniform(0.01, 0.3),
    'model__max_depth': randint(3, 7),
    'model__subsample': [0.8, 1.0], 
    'model__colsample_bytree': [0.8, 1.0]
}
params_rf = {
    'model__n_estimators': randint(100, 300),
    'model__max_depth': [10, 20, 30, None],
    'model__min_samples_split': randint(2, 10),
    'model__min_samples_leaf': randint(1, 5),
    'model__max_features': ['sqrt', 'log2']
}
params_svm = {
    'model__C': loguniform(1e-2, 1e3),
    'model__kernel': ['rbf', 'linear'],
    'model__gamma': ['scale', 'auto']
}

params_dt = {
    'model__max_depth': randint(3, 15),
    'model__min_samples_split': randint(2, 20),
    'model__criterion': ['gini', 'entropy']
}

models_to_tune = {
    "LogisticRegression": (pipe_lr, params_lr),
    "SVM": (pipe_svm, params_svm),
    "MLP": (pipe_mlp, params_mlp),
    "XGBoost": (pipe_xgb, params_xgb),
    "RandomForest": (pipe_rf, params_rf),
    "DecisionTree": (pipe_dt, params_dt),
}

all_results = {
    "LogisticRegression": {"acc": [], "f1_w": [], "f1_m": [], "recall": [], "auc": [], "precision": [], "aupr": []},
    "MLP": {"acc": [], "f1_w": [], "f1_m": [], "recall": [], "auc": [], "precision": [], "aupr": []},
    "XGBoost": {"acc": [], "f1_w": [], "f1_m": [], "recall": [], "auc": [], "precision": [], "aupr": []},
    "RandomForest": {"acc": [], "f1_w": [], "f1_m": [], "recall": [], "auc": [], "precision": [], "aupr": []},
    "SVM": {"acc": [], "f1_w": [], "f1_m": [], "recall": [], "auc": [], "precision": [], "aupr": []},
    "DecisionTree": {"acc": [], "f1_w": [], "f1_m": [], "recall": [], "auc": [], "precision": [], "aupr": []},
}



for model_name, (pipeline, param_dist) in models_to_tune.items():
    print(f"Evaluating model with nested CV: {model_name}")
    
    outer_cv = RepeatedStratifiedKFold(n_splits=5, n_repeats=5, random_state=SEED)

    # Inner fold loop is for finding the best hyperparameters per fold.
    for fold_idx, (train_idx, test_idx) in enumerate(outer_cv.split(X, y), start=1):
        X_train_outer, X_test_outer = X.iloc[train_idx], X.iloc[test_idx]
        y_train_outer, y_test_outer = y.iloc[train_idx], y.iloc[test_idx]
        inner_cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=SEED)
        
        best_score_fold = -np.inf
        best_params_fold = None

        # A fixed seed = same hyperparamters each fold.
        # No seed = different hyperparameters each fold. 
        # This adds more randomness, and may yield better generalization.
        param_sampler = list(ParameterSampler(param_dist, n_iter=20))
        
        for params in param_sampler:
            inner_scores = []
            
            for inner_train_idx, inner_val_idx in inner_cv.split(X_train_outer, y_train_outer):
                X_train_inner = X_train_outer.iloc[inner_train_idx]
                X_val_inner = X_train_outer.iloc[inner_val_idx]
                y_train_inner = y_train_outer.iloc[inner_train_idx]
                y_val_inner = y_train_outer.iloc[inner_val_idx]
                
                inner_pipeline = clone(pipeline)
                inner_pipeline.set_params(**params)
                inner_pipeline.fit(X_train_inner, y_train_inner)
                
                y_val_pred = inner_pipeline.predict(X_val_inner)
                score = f1_score(y_val_inner, y_val_pred, average='macro', zero_division=0)
                inner_scores.append(score)
            
            mean_score = np.mean(inner_scores)
            if mean_score > best_score_fold:
                best_score_fold = mean_score
                best_params_fold = params
        
        print(f"Outer fold {fold_idx}: best params (inner CV F1-M={best_score_fold:.4f})")
        print(f"{best_params_fold}")

        final_pipeline = clone(pipeline)
        final_pipeline.set_params(**best_params_fold)
        final_pipeline.fit(X_train_outer, y_train_outer)
        
        preds = final_pipeline.predict(X_test_outer)
        
        if hasattr(final_pipeline, "predict_proba"):
            proba = final_pipeline.predict_proba(X_test_outer)
        else:
            proba = None
        
        acc = accuracy_score(y_test_outer, preds)
        f1_w = f1_score(y_test_outer, preds, average='weighted', zero_division=0)
        f1_m = f1_score(y_test_outer, preds, average='macro', zero_division=0)
        recall = recall_score(y_test_outer, preds, average='macro', zero_division=0)
        precision = precision_score(y_test_outer, preds, average='macro', zero_division=0)
        
        auc = np.nan
        aupr = np.nan
        
        if proba is not None:
            try:
                if len(np.unique(y)) == 2:
                    auc = roc_auc_score(y_test_outer, proba[:, 1])
                    aupr = average_precision_score(y_test_outer, proba[:, 1])
                else:
                    auc = roc_auc_score(y_test_outer, proba, multi_class='ovr', average='macro')
                    y_test_bin = label_binarize(y_test_outer, classes=np.unique(y))
                    aupr = average_precision_score(y_test_bin, proba, average='macro')
            except Exception:
                auc = np.nan
                aupr = np.nan

        print(f"Fold {fold_idx} results: Acc={acc:.4f}, F1-W={f1_w:.4f}, "
            f"F1-M={f1_m:.4f}, Recall={recall:.4f}, Precision={precision:.4f}, AUC={auc:.4f}, AUPR={aupr:.4f}")
        
        all_results[model_name]["acc"].append(acc)
        all_results[model_name]["f1_w"].append(f1_w)
        all_results[model_name]["f1_m"].append(f1_m)
        all_results[model_name]["recall"].append(recall)
        all_results[model_name]["precision"].append(precision)
        all_results[model_name]["auc"].append(auc)
        all_results[model_name]["aupr"].append(aupr)

print("\nFINAL BASELINE RESULTS\n")
for model_name, metrics in all_results.items():
    avg_acc = np.mean(metrics["acc"])
    std_acc = np.std(metrics["acc"])
    avg_f1_w = np.mean(metrics["f1_w"])
    std_f1_w = np.std(metrics["f1_w"])
    avg_f1_m = np.mean(metrics["f1_m"])
    std_f1_m = np.std(metrics["f1_m"])
    avg_recall = np.mean(metrics["recall"])
    std_recall = np.std(metrics["recall"])
    avg_prec = np.mean(metrics["precision"])
    std_prec = np.std(metrics["precision"])
    avg_auc = np.nanmean(metrics["auc"])
    std_auc = np.nanstd(metrics["auc"])
    avg_aupr = np.nanmean(metrics["aupr"])
    std_aupr = np.nanstd(metrics["aupr"])
    
    print(f"\n{model_name}:")
    print(f"Accuracy: {avg_acc:.4f} +/- {std_acc:.4f}")
    print(f"Precision: {avg_prec:.4f} +/- {std_prec:.4f}")
    print(f"F1 Weighted: {avg_f1_w:.4f} +/- {std_f1_w:.4f}")
    print(f"F1 Macro: {avg_f1_m:.4f} +/- {std_f1_m:.4f}")
    print(f"Recall: {avg_recall:.4f} +/- {std_recall:.4f}")
    print(f"AUC: {avg_auc:.4f} +/- {std_auc:.4f}")
    print(f"AUPR: {avg_aupr:.4f} +/- {std_aupr:.4f}")

Phenotype-Driven Subgraph Detection:

Using HybridLouvain module, we iteratively pruned the global Spearman network. This process systematically isolates biologically meaningful subgraphs that maximize correlation with patient vital status.

from bioneuralnet.clustering import HybridLouvain

hl2 = HybridLouvain(
    G=spearman_12,
    B=omics_lgg,
    Y=Y_labels,
    k_L=0.7,
    k_P=0.8,
    max_iter=10,
    min_nodes=5,
    seed=SEED,
)

subnetworks2 = hl2.run(as_dfs=True)

from bioneuralnet.metrics import plot_network, cluster_correlation

driver_module2 = subnetworks2[4]
cluster_plot2 = cluster_correlation(driver_module2)

lgg_mapping2 = plot_network(
    cluster_plot2, 
    weight_threshold=0.35, 
    show_labels=True,
    show_edge_weights=True
)

print(lgg_mapping2.head(10))