diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..5666df5 --- /dev/null +++ b/.gitignore @@ -0,0 +1,29 @@ +# Byte-compiled / optimized files +__pycache__/ +*.py[cod] +*$py.class +*.egg-info/ +.ipynb_checkpoints/ + +# Virtual environments +.venv/ +venv/ +env/ + +# Raw datasets (large; download separately, see data/README.md) +data/* +!data/README.md + +# Generated artefacts (CSVs, models, figures, reports) +outputs/ + +# Feature-engineering cache (auto-generated by script 04) +.cache/ +**/.cache/ + +# OS / editor cruft +.DS_Store +Thumbs.db +*.swp +.idea/ +.vscode/ diff --git a/README.md b/README.md index 67e12f6..baae236 100644 --- a/README.md +++ b/README.md @@ -9,6 +9,23 @@ This repository contains the preprocessing, feature engineering, and modeling sc - `data/` : Directory intended for raw datasets (e.g., dbNSFP, hg19.fa). Note that datasets are excluded from Git due to their size. - `outputs/` : Directory for all generated CSVs, models, and figures. +## Pipeline stages + +The scripts in `src/` run in order; each reads the previous stage's output and writes to its own `outputs/` subfolder. + +| Script | Stage | +| --- | --- | +| `01_dbnsfp_processor.py` | Parse and filter dbNSFP into the base variant table | +| `02_remove_missing_values.py` | Drop rows with missing values (outside retained columns) | +| `03_remove_duplicates.py` | Remove duplicate variants | +| `04_feature_engineering.py` | Structural / functional feature engineering (AlphaFold, sequence, FCGR) | +| `05_remove_leakage.py` | Drop label-leaking annotation columns | +| `06_clean_and_finalize.py` | Final cleaning and column preparation | +| `07_dataset_balancing.py` | Balance the pathogenic / benign classes | +| `08_tda_fuzzy_ensemble.py` | TDA + fuzzy KNORA-E ensemble: train and evaluate | +| `09_reviewer_experiments.py` | Supplementary audits, ablations, and diversity analysis | +| `10_exploratory_data_analysis.py` | Exploratory figures and summary statistics | + ## Data Setup Due to GitHub file size limits, the 76GB of raw data is not tracked in this repository. diff --git a/src/01_dbnsfp_processor.py b/src/01_dbnsfp_processor.py index c77737e..f3db636 100644 --- a/src/01_dbnsfp_processor.py +++ b/src/01_dbnsfp_processor.py @@ -1,3 +1,5 @@ +"""Stage 01 - parse and filter dbNSFP records into the base somatic-variant table.""" + import pandas as pd from config import DATA_DIR, STAGE01_OUT import numpy as np @@ -135,13 +137,9 @@ def explode_semicolon_columns(variant_df, cols_to_explode): cgc = pd.read_csv(DATA_DIR / "Cosmic_CancerGeneCensus_v102_GRCh37.tsv", sep="\t") cancer_genes = set(cgc["GENE_SYMBOL"].unique()) oncogenes = set( - cgc[cgc["ROLE_IN_CANCER"].str.contains("oncogene", na=False, case=False)][ - "GENE_SYMBOL" - ] -) -tsgs = set( - cgc[cgc["ROLE_IN_CANCER"].str.contains("TSG", na=False, case=False)]["GENE_SYMBOL"] + cgc[cgc["ROLE_IN_CANCER"].str.contains("oncogene", na=False, case=False)]["GENE_SYMBOL"] ) +tsgs = set(cgc[cgc["ROLE_IN_CANCER"].str.contains("TSG", na=False, case=False)]["GENE_SYMBOL"]) tier1_genes = set(cgc[cgc["TIER"] == 1]["GENE_SYMBOL"]) cgc_slim = cgc[["GENE_SYMBOL", "ROLE_IN_CANCER", "TIER"]].copy() @@ -177,9 +175,7 @@ def explode_semicolon_columns(variant_df, cols_to_explode): .reset_index() ) cmc_agg.columns = ["genename", "aapos", "COSMIC_RECURRENCE", "COSMIC_TESTED"] -cmc_agg["COSMIC_FREQUENCY"] = cmc_agg["COSMIC_RECURRENCE"] / cmc_agg[ - "COSMIC_TESTED" -].replace(0, 1) +cmc_agg["COSMIC_FREQUENCY"] = cmc_agg["COSMIC_RECURRENCE"] / cmc_agg["COSMIC_TESTED"].replace(0, 1) cosmic_rescue = cmc_agg[cmc_agg["COSMIC_RECURRENCE"] >= MIN_COSMIC_RESCUE].copy() cosmic_rescue["rescue_gene_pos"] = ( @@ -258,8 +254,7 @@ def explode_semicolon_columns(variant_df, cols_to_explode): raw_chunk["gnomAD_non_cancer_AF"] = np.nan non_cancer_mask = ( - pd.to_numeric(raw_chunk["gnomAD2.1.1_exomes_non_cancer_AN"], errors="coerce") - > 0 + pd.to_numeric(raw_chunk["gnomAD2.1.1_exomes_non_cancer_AN"], errors="coerce") > 0 ) if non_cancer_mask.any(): raw_chunk.loc[non_cancer_mask, "gnomAD_non_cancer_AF"] = ( @@ -282,9 +277,7 @@ def explode_semicolon_columns(variant_df, cols_to_explode): + ":" + raw_chunk["alt"] ) - raw_chunk["rescue_gene_pos"] = ( - raw_chunk["genename"] + "_" + raw_chunk["aapos"].astype(str) - ) + raw_chunk["rescue_gene_pos"] = raw_chunk["genename"] + "_" + raw_chunk["aapos"].astype(str) in_civic = raw_chunk["rescue_key"].isin(rescue_keys_civic) in_cosmic = raw_chunk["rescue_gene_pos"].isin(cosmic_rescue_keys) @@ -293,22 +286,16 @@ def explode_semicolon_columns(variant_df, cols_to_explode): possible_hotspots = raw_chunk[raw_chunk["genename"].isin(RESCUE_HOTSPOTS.keys())] if not possible_hotspots.empty: for gene, positions in RESCUE_HOTSPOTS.items(): - mask = (raw_chunk["genename"] == gene) & ( - raw_chunk["aapos"].isin(positions) - ) + mask = (raw_chunk["genename"] == gene) & (raw_chunk["aapos"].isin(positions)) in_hotspot = in_hotspot | mask is_rescued = in_civic | in_hotspot | in_cosmic - af_to_use = raw_chunk["gnomAD_non_cancer_AF"].fillna( - raw_chunk["gnomAD4.1_joint_AF"] - ) + af_to_use = raw_chunk["gnomAD_non_cancer_AF"].fillna(raw_chunk["gnomAD4.1_joint_AF"]) passes_af = (af_to_use.isna()) | (af_to_use < GNOMAD_AF_THRESHOLD) quality_ok = ( - (raw_chunk["CADD_phred"] >= MIN_CADD_SCORE) - | (raw_chunk["CADD_phred"].isna()) - | is_rescued + (raw_chunk["CADD_phred"] >= MIN_CADD_SCORE) | (raw_chunk["CADD_phred"].isna()) | is_rescued ) keep_variant = passes_af | is_rescued @@ -329,13 +316,11 @@ def explode_semicolon_columns(variant_df, cols_to_explode): somatic_subset["HGVSp_snpEff"].str.contains("Ter|\\*", na=False, regex=True) ) somatic_subset.loc[nonsense_mask, "variant_type"] = "nonsense" - frameshift_mask = somatic_subset["HGVSp_snpEff"].str.contains( - "fs", na=False, case=False - ) + frameshift_mask = somatic_subset["HGVSp_snpEff"].str.contains("fs", na=False, case=False) somatic_subset.loc[frameshift_mask, "variant_type"] = "frameshift" - synonymous_mask = ( - somatic_subset["aaref"] == somatic_subset["aaalt"] - ) & somatic_subset["aaref"].notna() + synonymous_mask = (somatic_subset["aaref"] == somatic_subset["aaalt"]) & somatic_subset[ + "aaref" + ].notna() somatic_subset.loc[synonymous_mask, "variant_type"] = "synonymous" somatic_subset = somatic_subset.merge(cgc_slim, on="genename", how="left") @@ -349,19 +334,11 @@ def explode_semicolon_columns(variant_df, cols_to_explode): ) somatic_subset["COSMIC_RECURRENCE"].fillna(0, inplace=True) - somatic_subset["IS_CANCER_GENE"] = ( - somatic_subset["genename"].isin(cancer_genes).astype(int) - ) - somatic_subset["IS_ONCOGENE"] = ( - somatic_subset["genename"].isin(oncogenes).astype(int) - ) + somatic_subset["IS_CANCER_GENE"] = somatic_subset["genename"].isin(cancer_genes).astype(int) + somatic_subset["IS_ONCOGENE"] = somatic_subset["genename"].isin(oncogenes).astype(int) somatic_subset["IS_TSG"] = somatic_subset["genename"].isin(tsgs).astype(int) - somatic_subset["IS_TIER1"] = ( - somatic_subset["genename"].isin(tier1_genes).astype(int) - ) - somatic_subset["IS_ONCOKB"] = ( - somatic_subset["genename"].isin(oncokb_genes).astype(int) - ) + somatic_subset["IS_TIER1"] = somatic_subset["genename"].isin(tier1_genes).astype(int) + somatic_subset["IS_ONCOKB"] = somatic_subset["genename"].isin(oncokb_genes).astype(int) somatic_subset["IS_KNOWN_HOTSPOT"] = in_hotspot[somatic_subset.index].astype(int) somatic_subset["CONSENSUS_SCORE"] = ( diff --git a/src/02_remove_missing_values.py b/src/02_remove_missing_values.py index e4282f9..68838f4 100644 --- a/src/02_remove_missing_values.py +++ b/src/02_remove_missing_values.py @@ -1,83 +1,88 @@ -#2. Removes all rows with any missing values +"""Stage 02 - drop rows with missing values outside a set of retained columns.""" + import pandas as pd from config import STAGE01_OUT, STAGE02_OUT import gc -INPUT_FILE = STAGE01_OUT / 'somatic_variant_dbNSFP.csv' -OUTPUT_FILE = STAGE02_OUT / 'somatic_variant_dbNSFP_Removes_missing_values.csv' +INPUT_FILE = STAGE01_OUT / "somatic_variant_dbNSFP.csv" +OUTPUT_FILE = STAGE02_OUT / "somatic_variant_dbNSFP_Removes_missing_values.csv" COLUMNS_TO_RETAIN_WITH_MISSING = [ - 'COSMIC_FREQUENCY', - 'COSMIC_RECURRENCE', - 'EVIDENCE_SOURCE', - 'RESCUE_REASON', - 'WAS_RESCUED', - 'IS_CLINVAR_PATHOGENIC', - 'IS_KNOWN_HOTSPOT', - 'genename', - 'aapos', - 'aaref', - 'aaalt', - 'Interpro_domain', - 'variant_type', - 'ROLE_IN_CANCER' -] + "COSMIC_FREQUENCY", + "COSMIC_RECURRENCE", + "EVIDENCE_SOURCE", + "RESCUE_REASON", + "WAS_RESCUED", + "IS_CLINVAR_PATHOGENIC", + "IS_KNOWN_HOTSPOT", + "genename", + "aapos", + "aaref", + "aaalt", + "Interpro_domain", + "variant_type", + "ROLE_IN_CANCER", +] def clean_dataset_selective(chunksize=100000): """ Remove rows with missing values in non-exempt columns - + Process data in chunks to handle large files efficiently """ print(f"Starting selective cleaning: {INPUT_FILE}") print(f"Preserving rows with missing values in: {len(COLUMNS_TO_RETAIN_WITH_MISSING)} columns") print(f"Chunk size: {chunksize:,}\n") - + try: first_chunk = True total_original = 0 total_kept = 0 chunk_num = 0 - + print("Processing chunks...") - + for chunk in pd.read_csv(INPUT_FILE, chunksize=chunksize, low_memory=False): chunk_num += 1 original_rows = len(chunk) total_original += original_rows - + all_columns = chunk.columns.tolist() - columns_to_check = [col for col in all_columns if col not in COLUMNS_TO_RETAIN_WITH_MISSING] - - chunk_clean = chunk.dropna(how='any', subset=columns_to_check) - + columns_to_check = [ + col for col in all_columns if col not in COLUMNS_TO_RETAIN_WITH_MISSING + ] + + chunk_clean = chunk.dropna(how="any", subset=columns_to_check) + kept_rows = len(chunk_clean) total_kept += kept_rows - + if not chunk_clean.empty: if first_chunk: - chunk_clean.to_csv(OUTPUT_FILE, mode='w', index=False, header=True) + chunk_clean.to_csv(OUTPUT_FILE, mode="w", index=False, header=True) first_chunk = False else: - chunk_clean.to_csv(OUTPUT_FILE, mode='a', index=False, header=False) - + chunk_clean.to_csv(OUTPUT_FILE, mode="a", index=False, header=False) + if chunk_num % 10 == 0: drop_in_chunk = original_rows - kept_rows - print(f" Chunk {chunk_num}: {original_rows:,} -> {kept_rows:,} (dropped: {drop_in_chunk:,})") - + print( + f" Chunk {chunk_num}: {original_rows:,} -> {kept_rows:,} (dropped: {drop_in_chunk:,})" + ) + del chunk, chunk_clean gc.collect() - + dropped = total_original - total_kept drop_pct = (dropped / total_original) * 100 if total_original > 0 else 0 - + print(f"\nCleaning complete") print(f" Total input rows: {total_original:,}") print(f" Rows retained: {total_kept:,}") print(f" Rows dropped: {dropped:,} ({drop_pct:.2f}%)") print(f"\nOutput saved to: {OUTPUT_FILE}") - + except FileNotFoundError: print(f"Error: File '{INPUT_FILE}' not found") except Exception as e: @@ -85,4 +90,4 @@ def clean_dataset_selective(chunksize=100000): if __name__ == "__main__": - clean_dataset_selective() \ No newline at end of file + clean_dataset_selective() diff --git a/src/03_remove_duplicates.py b/src/03_remove_duplicates.py index e5b45ef..b9f9f30 100644 --- a/src/03_remove_duplicates.py +++ b/src/03_remove_duplicates.py @@ -1,74 +1,72 @@ -#3. Removes duplicate variants +"""Stage 03 - remove duplicate variant records.""" + import pandas as pd from config import STAGE02_OUT, STAGE03_OUT import gc -INPUT_FILE = STAGE02_OUT / 'somatic_variant_dbNSFP_Removes_missing_values.csv' -OUTPUT_FILE = STAGE03_OUT / 'somatic_variant_dbNSFP_Removes_missing_values_Deduplication.csv' +INPUT_FILE = STAGE02_OUT / "somatic_variant_dbNSFP_Removes_missing_values.csv" +OUTPUT_FILE = STAGE03_OUT / "somatic_variant_dbNSFP_Removes_missing_values_Deduplication.csv" def remove_duplicates(chunksize=100000): print(f"Loading dataset: {INPUT_FILE}") print(f"Chunk size: {chunksize:,}\n") - + try: print("Pass 1: Reading and sorting chunks...") all_chunks = [] chunk_num = 0 total_rows = 0 - + for chunk in pd.read_csv(INPUT_FILE, chunksize=chunksize, low_memory=False): chunk_num += 1 total_rows += len(chunk) - - - if 'CONSENSUS_SCORE' in chunk.columns: + + if "CONSENSUS_SCORE" in chunk.columns: chunk = chunk.sort_values( - by=['chr', 'pos', 'ref', 'alt', 'CONSENSUS_SCORE'], - ascending=[True, True, True, True, False] + by=["chr", "pos", "ref", "alt", "CONSENSUS_SCORE"], + ascending=[True, True, True, True, False], ) - + all_chunks.append(chunk) - + if chunk_num % 10 == 0: print(f" Processed chunk {chunk_num}: {total_rows:,} total rows") - + del chunk gc.collect() - + print(f"\nTotal variants loaded: {total_rows:,}") - + print("\nPass 2: Combining and deduplicating...") variant_df = pd.concat(all_chunks, ignore_index=True) del all_chunks gc.collect() - - - if 'CONSENSUS_SCORE' in variant_df.columns: + + if "CONSENSUS_SCORE" in variant_df.columns: print(" Sorting by consensus score (highest first)") variant_df = variant_df.sort_values( - by=['chr', 'pos', 'ref', 'alt', 'CONSENSUS_SCORE'], - ascending=[True, True, True, True, False] + by=["chr", "pos", "ref", "alt", "CONSENSUS_SCORE"], + ascending=[True, True, True, True, False], ) - + # Remove duplicates print(" Removing duplicates...") deduplicated_df = variant_df.drop_duplicates( - subset=['chr', 'pos', 'ref', 'alt'], - keep='first' + subset=["chr", "pos", "ref", "alt"], keep="first" ) - + kept_rows = len(deduplicated_df) dropped_rows = total_rows - kept_rows - + print(f"\nUnique variants retained: {kept_rows:,}") print(f"Duplicates removed: {dropped_rows:,}") - + print(f"\nSaving to: {OUTPUT_FILE}") deduplicated_df.to_csv(OUTPUT_FILE, index=False) - + print("Deduplication complete") - + del variant_df, deduplicated_df gc.collect() @@ -77,8 +75,9 @@ def remove_duplicates(chunksize=100000): except Exception as e: print(f"Error occurred: {e}") import traceback + traceback.print_exc() if __name__ == "__main__": - remove_duplicates(chunksize=50000) \ No newline at end of file + remove_duplicates(chunksize=50000) diff --git a/src/04_feature_engineering.py b/src/04_feature_engineering.py index 571504d..ef949d6 100644 --- a/src/04_feature_engineering.py +++ b/src/04_feature_engineering.py @@ -1,3 +1,5 @@ +"""Stage 04 - structural and functional feature engineering (AlphaFold, sequence, FCGR).""" + import pandas as pd from config import DATA_DIR, STAGE03_OUT, STAGE04_OUT import numpy as np @@ -97,9 +99,7 @@ class OptimizedUniProtParser: def __init__(self, uniprot_file: str, cache_dir: str = None): self.uniprot_file = Path(uniprot_file) - self.cache_dir = ( - Path(cache_dir) if cache_dir else self.uniprot_file.parent / ".cache" - ) + self.cache_dir = Path(cache_dir) if cache_dir else self.uniprot_file.parent / ".cache" self.cache_dir.mkdir(exist_ok=True) # Pre-compiled regex patterns for speed @@ -208,17 +208,11 @@ def _parse_feature_line(self, line: str, features: ProteinFeatures): description = parts[3] if len(parts) > 3 else "" if feature_type == "DOMAIN": - features.domains.append( - {"start": start, "end": end, "name": description[:50]} - ) + features.domains.append({"start": start, "end": end, "name": description[:50]}) elif feature_type == "ACT_SITE": - features.active_sites.append( - {"position": start, "description": description} - ) + features.active_sites.append({"position": start, "description": description}) elif feature_type == "BINDING": - features.binding_sites.append( - {"position": start, "description": description} - ) + features.binding_sites.append({"position": start, "description": description}) elif feature_type == "TRANSMEM": features.transmembrane.append({"start": start, "end": end}) elif feature_type == "SIGNAL": @@ -271,10 +265,7 @@ def _build_index(self): # Priority: PDB > CIF, F1 > F2 > ... priority = (is_pdb, fragment == "F1", fragment) - if ( - uniprot_id not in file_priority - or priority > file_priority[uniprot_id][0] - ): + if uniprot_id not in file_priority or priority > file_priority[uniprot_id][0]: file_priority[uniprot_id] = (priority, pdb_file) self.index = {uid: path for uid, (_, path) in file_priority.items()} @@ -305,12 +296,13 @@ def get_available_ids(self) -> set: # Thread-local storage for PDB parser (thread-safe, no subprocess spawning) import threading + _thread_local = threading.local() def _get_parser(): """Get or create a thread-local PDB parser""" - if not hasattr(_thread_local, 'parser'): + if not hasattr(_thread_local, "parser"): _thread_local.parser = PDBParser(QUIET=True) if BIOPYTHON_AVAILABLE else None return _thread_local.parser @@ -336,10 +328,7 @@ def _process_structure_batch( # Use thread-local parser (for ThreadPoolExecutor) or global parser (for ProcessPoolExecutor) parser = _get_parser() if not BIOPYTHON_AVAILABLE or parser is None: - return [ - (idx, {"sasa": -1, "relative_sasa": -1, "plddt": -1}) - for idx, _ in residue_batch - ] + return [(idx, {"sasa": -1, "relative_sasa": -1, "plddt": -1}) for idx, _ in residue_batch] try: # Parse structure once for all residues @@ -405,23 +394,18 @@ def _process_structure_batch( row_idx, { "sasa": round(sasa, 2) if sasa > 0 else -1, - "relative_sasa": ( - round(relative_sasa, 3) if relative_sasa > 0 else -1 - ), + "relative_sasa": (round(relative_sasa, 3) if relative_sasa > 0 else -1), "plddt": round(plddt, 2) if plddt > 0 else -1, }, ) ) else: - results.append( - (row_idx, {"sasa": -1, "relative_sasa": -1, "plddt": -1}) - ) + results.append((row_idx, {"sasa": -1, "relative_sasa": -1, "plddt": -1})) except Exception: # Return defaults for all residues on error results = [ - (idx, {"sasa": -1, "relative_sasa": -1, "plddt": -1}) - for idx, _ in residue_batch + (idx, {"sasa": -1, "relative_sasa": -1, "plddt": -1}) for idx, _ in residue_batch ] return results @@ -527,7 +511,7 @@ def add_structural_features_parallel( n_workers: int = None, ) -> pd.DataFrame: """Add AlphaFold structural features using ThreadPoolExecutor. - + Uses threads instead of processes to avoid Windows memory issues where each spawned process re-imports pandas/numpy and exhausts virtual memory (paging file too small error). @@ -584,24 +568,17 @@ def add_structural_features_parallel( structure_batches[uniprot_id][1].append((idx, int(pos))) total_variants = sum(len(b[1]) for b in structure_batches.values()) - print( - f" Processing {len(structure_batches):,} structures for {total_variants:,} variants" - ) + print(f" Processing {len(structure_batches):,} structures for {total_variants:,} variants") # Prepare work items - work_items = [ - (uid, path, residues) for uid, (path, residues) in structure_batches.items() - ] + work_items = [(uid, path, residues) for uid, (path, residues) in structure_batches.items()] # Process with ThreadPoolExecutor (no subprocess spawning = no memory explosion) all_results = [] failed_count = 0 with ThreadPoolExecutor(max_workers=n_workers) as executor: - futures = { - executor.submit(_process_structure_batch, item): item[0] - for item in work_items - } + futures = {executor.submit(_process_structure_batch, item): item[0] for item in work_items} with tqdm(total=len(futures), desc=" Processing structures") as pbar: for future in as_completed(futures): @@ -620,9 +597,7 @@ def add_structural_features_parallel( plddt_values[row_idx] = features["plddt"] successful = np.sum(sasa_values > -1) - print( - f"\n Successfully annotated {successful:,} / {n:,} variants with structural features" - ) + print(f"\n Successfully annotated {successful:,} / {n:,} variants with structural features") if failed_count > 0: print(f" {failed_count} structures failed during processing") @@ -664,18 +639,21 @@ def build_smart_gene_mapping_optimized( if candidates[0][1]: genes_with_structure += 1 - print( - f" Mapped {len(optimized_features):,} genes ({genes_with_structure:,} with structures)" - ) + print(f" Mapped {len(optimized_features):,} genes ({genes_with_structure:,} with structures)") return optimized_features def enrich_dataset_with_structure_function_optimized( - input_csv: str = str(STAGE03_OUT / "somatic_variant_dbNSFP_Removes_missing_values_Deduplication.csv"), + input_csv: str = str( + STAGE03_OUT / "somatic_variant_dbNSFP_Removes_missing_values_Deduplication.csv" + ), uniprot_file: str = str(DATA_DIR / "uniprotkb_proteome_UP000005640_2026_01_07.txt"), alphafold_dir: str = str(DATA_DIR / "UP000005640_9606_HUMAN_v6"), - output_csv: str = str(STAGE04_OUT / "somatic_variant_dbNSFP_Removes_missing_values_Deduplication_Structural_Functional.csv"), + output_csv: str = str( + STAGE04_OUT + / "somatic_variant_dbNSFP_Removes_missing_values_Deduplication_Structural_Functional.csv" + ), n_workers: int = None, max_rows: int = None, force_reparse: bool = False, @@ -684,9 +662,7 @@ def enrich_dataset_with_structure_function_optimized( print("=" * 70) print("STRUCTURAL & FUNCTIONAL ENRICHMENT PIPELINE - OPTIMIZED") print(f"CPU Workers: {n_workers or CPU_COUNT}") - print( - f"GPU Acceleration: {'cuDF available' if CUDF_AVAILABLE else 'Not available'}" - ) + print(f"GPU Acceleration: {'cuDF available' if CUDF_AVAILABLE else 'Not available'}") print("=" * 70) # Load dataset @@ -719,9 +695,7 @@ def enrich_dataset_with_structure_function_optimized( alphafold_index = AlphaFoldStructureIndex(alphafold_dir) # Smart mapping - optimized_features = build_smart_gene_mapping_optimized( - uniprot_features, alphafold_index - ) + optimized_features = build_smart_gene_mapping_optimized(uniprot_features, alphafold_index) # Annotate with UniProt features (vectorized) print("\nSTEP 3: UniProt Annotation") @@ -784,7 +758,9 @@ def annotate_with_cudf(dataset_path: str, output_path: str) -> None: arg_parser.add_argument( "--input", "-i", - default=str(STAGE03_OUT / "somatic_variant_dbNSFP_Removes_missing_values_Deduplication.csv"), + default=str( + STAGE03_OUT / "somatic_variant_dbNSFP_Removes_missing_values_Deduplication.csv" + ), help="Input CSV file path", ) arg_parser.add_argument( @@ -802,7 +778,10 @@ def annotate_with_cudf(dataset_path: str, output_path: str) -> None: arg_parser.add_argument( "--output", "-o", - default=str(STAGE04_OUT / "somatic_variant_dbNSFP_Removes_missing_values_Deduplication_Structural_Functional.csv"), + default=str( + STAGE04_OUT + / "somatic_variant_dbNSFP_Removes_missing_values_Deduplication_Structural_Functional.csv" + ), help="Output CSV file path", ) arg_parser.add_argument( diff --git a/src/05_remove_leakage.py b/src/05_remove_leakage.py index a30a632..12017bc 100644 --- a/src/05_remove_leakage.py +++ b/src/05_remove_leakage.py @@ -1,83 +1,94 @@ +"""Stage 05 - drop annotation columns that leak the target label.""" + import pandas as pd from config import STAGE04_OUT, STAGE05_OUT import gc from pathlib import Path -INPUT_FILE = STAGE04_OUT / 'somatic_variant_dbNSFP_Removes_missing_values_Deduplication_Structural_Functional.csv' -OUTPUT_FILE = STAGE05_OUT / 'somatic_variant_Leakage_Removed.csv' +INPUT_FILE = ( + STAGE04_OUT + / "somatic_variant_dbNSFP_Removes_missing_values_Deduplication_Structural_Functional.csv" +) +OUTPUT_FILE = STAGE05_OUT / "somatic_variant_Leakage_Removed.csv" LEAKAGE_COLS = [ - 'COSMIC_FREQUENCY', - 'COSMIC_RECURRENCE', - 'EVIDENCE_SOURCE', - 'RESCUE_REASON', - 'WAS_RESCUED', - 'IS_CLINVAR_PATHOGENIC', - 'IS_KNOWN_HOTSPOT' + "COSMIC_FREQUENCY", + "COSMIC_RECURRENCE", + "EVIDENCE_SOURCE", + "RESCUE_REASON", + "WAS_RESCUED", + "IS_CLINVAR_PATHOGENIC", + "IS_KNOWN_HOTSPOT", ] ID_COLS = [ - 'genename', 'aapos', 'aaref', 'aaalt', - 'Interpro_domain', 'variant_type', - 'ROLE_IN_CANCER' + "genename", + "aapos", + "aaref", + "aaalt", + "Interpro_domain", + "variant_type", + "ROLE_IN_CANCER", ] + def drop_columns_chunked(chunksize=50000): """Process large CSV in chunks to minimize memory usage""" print(f"Starting feature selection: {INPUT_FILE}") print(f"Chunk size: {chunksize:,} rows") - + try: # Create output dir if not exists Path(OUTPUT_FILE).parent.mkdir(parents=True, exist_ok=True) - + all_cols = pd.read_csv(INPUT_FILE, nrows=0).columns.tolist() cols_to_drop = LEAKAGE_COLS + ID_COLS use_cols = [c for c in all_cols if c not in cols_to_drop] - + print(f"Total columns: {len(all_cols)}") print(f"Dropping: {len(cols_to_drop)} columns") print(f"Retaining: {len(use_cols)} columns") print(f"\nRetained features:") for i, col in enumerate(use_cols, 1): print(f" {i}. {col}") - + first_chunk = True total_rows = 0 chunk_num = 0 - + print(f"\nProcessing chunks...") - - for chunk in pd.read_csv(INPUT_FILE, - usecols=use_cols, - chunksize=chunksize, - low_memory=False): + + for chunk in pd.read_csv( + INPUT_FILE, usecols=use_cols, chunksize=chunksize, low_memory=False + ): chunk_num += 1 total_rows += len(chunk) - + if first_chunk: - chunk.to_csv(OUTPUT_FILE, mode='w', index=False, header=True) + chunk.to_csv(OUTPUT_FILE, mode="w", index=False, header=True) first_chunk = False print(f" Chunk {chunk_num}: {len(chunk):,} rows (with header)") else: - chunk.to_csv(OUTPUT_FILE, mode='a', index=False, header=False) + chunk.to_csv(OUTPUT_FILE, mode="a", index=False, header=False) print(f" Chunk {chunk_num}: {len(chunk):,} rows (appended)") - + del chunk gc.collect() - + print(f"\nProcessing complete. Total rows: {total_rows:,}") print(f"Output saved to: {OUTPUT_FILE}") print("\nFinal dataset summary:") - + final_df = pd.read_csv(OUTPUT_FILE, nrows=3) print(f" Shape: {final_df.shape}") print(f" Columns: {list(final_df.columns)}") - + except Exception as e: print(f"\nError occurred: {e}") import traceback + traceback.print_exc() + if __name__ == "__main__": drop_columns_chunked(chunksize=50000) diff --git a/src/06_clean_and_finalize.py b/src/06_clean_and_finalize.py index 758b0ad..1c5d72b 100644 --- a/src/06_clean_and_finalize.py +++ b/src/06_clean_and_finalize.py @@ -1,3 +1,5 @@ +"""Stage 06 - final cleaning and column preparation of the modelling dataset.""" + import pandas as pd from config import STAGE05_OUT, STAGE06_OUT import numpy as np @@ -6,11 +8,12 @@ INPUT_FILE = STAGE05_OUT / "somatic_variant_Leakage_Removed.csv" OUTPUT_FILE = STAGE06_OUT / "somatic_variant_Cleaned.csv" + def clean_dataset(): - print("="*60) + print("=" * 60) print("FINAL DATASET CLEANING & PREPARATION") - print("="*60) - + print("=" * 60) + try: Path(OUTPUT_FILE).parent.mkdir(parents=True, exist_ok=True) @@ -18,20 +21,20 @@ def clean_dataset(): print(f"Loading dataset: {INPUT_FILE}...") df = pd.read_csv(INPUT_FILE) print(f"Original Shape: {df.shape} (Rows, Columns)\n") - + # Step 1: Drop DOMAIN_NAME Column - if 'DOMAIN_NAME' in df.columns: - df.drop(columns=['DOMAIN_NAME'], inplace=True) + if "DOMAIN_NAME" in df.columns: + df.drop(columns=["DOMAIN_NAME"], inplace=True) print(" Step 1: Dropped 'DOMAIN_NAME' column.") else: print(" Step 1: 'DOMAIN_NAME' column not found (skipped).") - + # Step 2: Drop columns having only 'U' value u_cols_dropped = [] for col in df.columns: - if df[col].dtype == 'object' and df[col].nunique() == 1 and df[col].unique()[0] == 'U': + if df[col].dtype == "object" and df[col].nunique() == 1 and df[col].unique()[0] == "U": u_cols_dropped.append(col) - + if u_cols_dropped: df.drop(columns=u_cols_dropped, inplace=True) print(f" Step 2: Dropped columns with only 'U' values: {u_cols_dropped}") @@ -39,27 +42,28 @@ def clean_dataset(): print(" Step 2: No columns found with only 'U' values.") # Step 3: Remove rows with -1 values (Missing Structural Data) - if 'SASA' in df.columns: + if "SASA" in df.columns: initial_count = len(df) - df = df[df['SASA'] != -1] + df = df[df["SASA"] != -1] removed_count = initial_count - len(df) print(f"Step 3: Removed {removed_count:,} rows with missing structural data (-1).") else: print("Warning: 'SASA' column not found, could not filter -1 rows.") # Final Save - print("\n" + "-"*60) + print("\n" + "-" * 60) print(f"Final Shape: {df.shape} (Rows, Columns)") print(f"Saving to: {OUTPUT_FILE}") - + df.to_csv(OUTPUT_FILE, index=False) print(" Done! File is clean and ready for Machine Learning.") - print("="*60) + print("=" * 60) except FileNotFoundError: print(f"Error: File '{INPUT_FILE}' not found!") except Exception as e: print(f"An error occurred: {e}") + if __name__ == "__main__": clean_dataset() diff --git a/src/07_dataset_balancing.py b/src/07_dataset_balancing.py index d87fbb8..5e2254d 100644 --- a/src/07_dataset_balancing.py +++ b/src/07_dataset_balancing.py @@ -1,59 +1,63 @@ +"""Stage 07 - balance the pathogenic/benign classes for training.""" + import pandas as pd from config import STAGE06_OUT, STAGE07_OUT import numpy as np from pathlib import Path -INPUT_FILE = STAGE06_OUT / 'somatic_variant_Cleaned.csv' -OUTPUT_FILE = STAGE07_OUT / 'Final_Dataset_Balanced.csv' +INPUT_FILE = STAGE06_OUT / "somatic_variant_Cleaned.csv" +OUTPUT_FILE = STAGE07_OUT / "Final_Dataset_Balanced.csv" + def balance_dataset(): print(f"Loading dataset: {INPUT_FILE}") - + try: Path(OUTPUT_FILE).parent.mkdir(parents=True, exist_ok=True) variant_df = pd.read_csv(INPUT_FILE, low_memory=False) - + print("\nOriginal class distribution:") print("-" * 40) - counts = variant_df['LABEL_PATHOGENIC'].value_counts() + counts = variant_df["LABEL_PATHOGENIC"].value_counts() total = len(variant_df) - + n_benign = counts.get(0, 0) n_pathogenic = counts.get(1, 0) - + print(f" Total samples: {total:,}") print(f" Benign (0): {n_benign:,} ({n_benign/total*100:.2f}%)") print(f" Pathogenic (1): {n_pathogenic:,} ({n_pathogenic/total*100:.2f}%)") print("-" * 40) - + print("\nBalancing dataset via downsampling...") - - pathogenic_variants = variant_df[variant_df['LABEL_PATHOGENIC'] == 1] - benign_variants = variant_df[variant_df['LABEL_PATHOGENIC'] == 0] - + + pathogenic_variants = variant_df[variant_df["LABEL_PATHOGENIC"] == 1] + benign_variants = variant_df[variant_df["LABEL_PATHOGENIC"] == 0] + sample_size = len(pathogenic_variants) benign_downsampled = benign_variants.sample(n=sample_size, random_state=42) - + balanced_df = pd.concat([pathogenic_variants, benign_downsampled]) balanced_df = balanced_df.sample(frac=1, random_state=42).reset_index(drop=True) - + print("\nBalanced class distribution:") print("-" * 40) - new_counts = balanced_df['LABEL_PATHOGENIC'].value_counts() + new_counts = balanced_df["LABEL_PATHOGENIC"].value_counts() print(f" Total samples: {len(balanced_df):,}") print(f" Benign (0): {new_counts.get(0, 0):,}") print(f" Pathogenic (1): {new_counts.get(1, 0):,}") print("-" * 40) - + print(f"\nSaving to: {OUTPUT_FILE}") balanced_df.to_csv(OUTPUT_FILE, index=False) - + print("Balancing complete") - + except FileNotFoundError: print(f"Error: File '{INPUT_FILE}' not found") except Exception as e: print(f"Error occurred: {e}") + if __name__ == "__main__": balance_dataset() diff --git a/src/08_tda_fuzzy_ensemble.py b/src/08_tda_fuzzy_ensemble.py index 2f68b50..cb9acfe 100644 --- a/src/08_tda_fuzzy_ensemble.py +++ b/src/08_tda_fuzzy_ensemble.py @@ -1,3 +1,5 @@ +"""Stage 08 - TDA + fuzzy KNORA-E ensemble: feature build, training, and evaluation.""" + import pandas as pd from config import DATA_DIR, STAGE07_OUT, STAGE08_OUT import numpy as np @@ -19,7 +21,7 @@ from sklearn.neighbors import KNeighborsClassifier, NearestNeighbors from sklearn.naive_bayes import GaussianNB from sklearn.svm import SVC -from sklearn.neural_network import MLPClassifier # NEW (Change 1) +from sklearn.neural_network import MLPClassifier # NEW (Change 1) from joblib import Parallel, delayed import multiprocessing from sklearn.model_selection import ( @@ -63,6 +65,7 @@ try: from xgboost import XGBClassifier + XGBOOST_AVAILABLE = True except ImportError: print("Warning: XGBoost not available") @@ -70,6 +73,7 @@ try: from lightgbm import LGBMClassifier + LGBM_AVAILABLE = True except ImportError: print("Warning: LightGBM not available") @@ -77,6 +81,7 @@ try: from catboost import CatBoostClassifier + CATBOOST_AVAILABLE = True except ImportError: print("Warning: CatBoost not available") @@ -85,6 +90,7 @@ try: from gtda.homology import VietorisRipsPersistence from gtda.diagrams import PersistenceEntropy, Amplitude, BettiCurve, PersistenceLandscape + TDA_AVAILABLE = True except ImportError: print("Warning: giotto-tda not available. Install with: pip install giotto-tda") @@ -92,6 +98,7 @@ try: from pyfaidx import Fasta + PYFAIDX_AVAILABLE = True except ImportError: print("Warning: pyfaidx not available. Install with: pip install pyfaidx") @@ -99,6 +106,7 @@ try: import shap + SHAP_AVAILABLE = True except ImportError: print("Warning: SHAP not available. Install with: pip install shap") @@ -108,7 +116,7 @@ class Config: - """Configuration with optimized hyperparameters — v4.0 (Supervisor Revision)""" + """Model, feature-engineering, and training hyper-parameters.""" RANDOM_STATE = 42 TEST_SIZE = 0.2 @@ -119,39 +127,30 @@ class Config: TDA_N_NEIGHBORS = 75 SEQUENCE_WINDOW = 75 FCGR_K_VALUES = [3, 4] - # ----------------------------------------------------------------------- - # CHANGE 2: TDA homology dims kept at (0, 1); the extractor now produces + # TDA homology dims kept at (0, 1); the extractor now produces # ~24 features instead of 6, all computed over the FCGR representation. - # ----------------------------------------------------------------------- TDA_HOMOLOGY_DIMS = (0, 1) - # FIX 3 (retained from v3) FCGR_PCA_COMPONENTS = 20 - # FIX 2 (retained from v3) WEAK_LEARNER_MARGIN = 0.10 - # FIX 5 (retained from v3) - THRESHOLD_SEARCH_LOW = 0.45 + THRESHOLD_SEARCH_LOW = 0.45 THRESHOLD_SEARCH_HIGH = 0.65 - # ----------------------------------------------------------------------- - # CHANGE 3: Selective prediction — abstain when neighbourhood purity is - # below this threshold. Set to 0.0 to disable (returns all predictions). - # ----------------------------------------------------------------------- - SELECTIVE_ABSTAIN_THRESHOLD = 0.60 # neighbourhood purity ≥ 0.60 to predict + # Selective prediction — abstain when neighbourhood purity is + # below this threshold. Set to 0.0 to disable (returns all predictions). + SELECTIVE_ABSTAIN_THRESHOLD = 0.60 # neighbourhood purity ≥ 0.60 to predict - # ----------------------------------------------------------------------- - # CHANGE 3: Consensus margin — variants whose |CONSENSUS_SCORE − 0.5| is + # Consensus margin — variants whose |CONSENSUS_SCORE − 0.5| is # below this value are treated as ambiguous; their sample weight during # base-learner training is scaled down to AMBIGUOUS_SAMPLE_WEIGHT. - # ----------------------------------------------------------------------- - AMBIGUOUS_MARGIN = 0.10 + AMBIGUOUS_MARGIN = 0.10 AMBIGUOUS_SAMPLE_WEIGHT = 0.30 - DATA_PATH = STAGE07_OUT / "Final_Dataset_Balanced.csv" + DATA_PATH = STAGE07_OUT / "Final_Dataset_Balanced.csv" GENOME_PATH = DATA_DIR / "hg19.fa" - OUTPUT_DIR = Path(STAGE08_OUT) + OUTPUT_DIR = Path(STAGE08_OUT) LEAKAGE_COLS = ["chr", "pos", "ref", "alt", "CONSENSUS_SCORE", "TIER"] @@ -190,16 +189,15 @@ class Config: sns.set_palette("colorblind") -# --------------------------------------------------------------------------- -# FCGR encoders (unchanged) -# --------------------------------------------------------------------------- +# FCGR encoders (unchanged) + class FCGREncoder: """Frequency Chaos Game Representation encoder""" def __init__(self, k=4): self.k = k - self.grid_size = 2 ** k + self.grid_size = 2**k self.vertices = { "A": (0.0, 0.0), "C": (0.0, 1.0), @@ -243,38 +241,34 @@ def encode_variants(self, sequences): return np.hstack(encoded_features) -# --------------------------------------------------------------------------- -# CHANGE 2: Expanded TopologicalFeatureExtractor -# -# Key changes vs v3: -# 1. TDA is now computed over the FCGR/sequence representation (X_fcgr_raw -# is passed in as `X_fcgr` from TFDFEPreprocessor) instead of the -# biological feature space. Topology over the genomic sequence captures -# information that tree-based models trained on scalar biological scores -# fundamentally cannot access. -# 2. Feature set expanded from 6 → ~24: -# - total_persistence_h{0,1} (2) -# - entropy_h{0,1} (2) -# - n_components_h{0,1} (2) -# - amplitude_h{0,1} (2) (max-normalised persistence amplitude) -# - betti_mean_h{0,1} (2) (mean Betti number over filtration) -# - landscape_l1_h{0,1} (2) (L1 norm of first persistence landscape) -# - birth_mean_h{0,1} (2) (mean birth time) -# - death_mean_h{0,1} (2) (mean death time) -# - persistence_std_h{0,1} (2) (std of persistence lifetimes) -# - persistence_max_h{0,1} (2) (longest-lived feature) -# Total: 20 features per sample (6 original + 14 new) -# --------------------------------------------------------------------------- +# Expanded TopologicalFeatureExtractor +# Key changes vs: +# 1. TDA is now computed over the FCGR/sequence representation (X_fcgr_raw +# is passed in as `X_fcgr` from TFDFEPreprocessor) instead of the +# biological feature space. Topology over the genomic sequence captures +# information that tree-based models trained on scalar biological scores +# fundamentally cannot access. +# 2. Feature set expanded from 6 → ~24: +# - total_persistence_h{0,1} (2) +# - entropy_h{0,1} (2) +# - n_components_h{0,1} (2) +# - amplitude_h{0,1} (2) (max-normalised persistence amplitude) +# - betti_mean_h{0,1} (2) (mean Betti number over filtration) +# - landscape_l1_h{0,1} (2) (L1 norm of first persistence landscape) +# - birth_mean_h{0,1} (2) (mean birth time) +# - death_mean_h{0,1} (2) (mean death time) +# - persistence_std_h{0,1} (2) (std of persistence lifetimes) +# - persistence_max_h{0,1} (2) (longest-lived feature) +# Total: 20 features per sample (6 original + 14 new) + class TopologicalFeatureExtractor: - """Extract topological features using persistent homology. - - CHANGE 2 (v4.0): - ---------------- - - Now operates on the **FCGR/sequence representation** rather than the - biological feature space. Call `extract_global_tda(X_fcgr)` where - X_fcgr is the PCA-compressed FCGR matrix (n_samples × 20). - - Feature vector expanded to ~20 descriptors per sample (was 6). + """Extract topological features via persistent homology. + + Operates on the FCGR/sequence representation rather than the biological + feature space: call ``extract_global_tda(X_fcgr)`` where ``X_fcgr`` is the + PCA-compressed FCGR matrix (n_samples x 20). Produces ~20 topological + descriptors per sample. """ def __init__(self, homology_dimensions=(0, 1), n_neighbors=75): @@ -302,43 +296,48 @@ def extract_local_tda(self, X, indices): diagram = diagrams[0][diagrams[0][:, 2] == dim] prefix = f"h{dim}" if len(diagram) > 0: - birth = diagram[:, 0] - death = diagram[:, 1] + birth = diagram[:, 0] + death = diagram[:, 1] persistence = death - birth # ---- original 3 features ---- total_pers = float(np.sum(persistence)) features[f"total_persistence_{prefix}"] = total_pers if total_pers > 0: - p_norm = persistence / total_pers - entropy = -np.sum(p_norm * np.log(p_norm + 1e-10)) + p_norm = persistence / total_pers + entropy = -np.sum(p_norm * np.log(p_norm + 1e-10)) features[f"entropy_{prefix}"] = float(entropy) else: features[f"entropy_{prefix}"] = 0.0 features[f"n_components_{prefix}"] = float(len(diagram)) - # ---- new features (CHANGE 2) ---- - features[f"amplitude_{prefix}"] = float(np.max(persistence)) - features[f"persistence_std_{prefix}"] = float(np.std(persistence)) - features[f"persistence_max_{prefix}"] = float(np.max(persistence)) - features[f"birth_mean_{prefix}"] = float(np.mean(birth)) - features[f"death_mean_{prefix}"] = float(np.mean(death)) + # ---- new features ---- + features[f"amplitude_{prefix}"] = float(np.max(persistence)) + features[f"persistence_std_{prefix}"] = float(np.std(persistence)) + features[f"persistence_max_{prefix}"] = float(np.max(persistence)) + features[f"birth_mean_{prefix}"] = float(np.mean(birth)) + features[f"death_mean_{prefix}"] = float(np.mean(death)) # Betti mean: proportion of filtration steps where feature lives # Approximate as ratio of mean persistence to total filtration span span = float(np.max(death) - np.min(birth)) if len(death) > 0 else 1.0 span = max(span, 1e-10) - features[f"betti_mean_{prefix}"] = float(np.mean(persistence) / span) + features[f"betti_mean_{prefix}"] = float(np.mean(persistence) / span) # L1 norm of first persistence landscape (simplified) sorted_p = np.sort(persistence)[::-1] - l1_land = float(np.sum(sorted_p * np.arange(1, len(sorted_p) + 1))) - features[f"landscape_l1_{prefix}"] = l1_land + l1_land = float(np.sum(sorted_p * np.arange(1, len(sorted_p) + 1))) + features[f"landscape_l1_{prefix}"] = l1_land else: for key in [ - f"total_persistence_{prefix}", f"entropy_{prefix}", - f"n_components_{prefix}", f"amplitude_{prefix}", - f"persistence_std_{prefix}", f"persistence_max_{prefix}", - f"birth_mean_{prefix}", f"death_mean_{prefix}", - f"betti_mean_{prefix}", f"landscape_l1_{prefix}", + f"total_persistence_{prefix}", + f"entropy_{prefix}", + f"n_components_{prefix}", + f"amplitude_{prefix}", + f"persistence_std_{prefix}", + f"persistence_max_{prefix}", + f"birth_mean_{prefix}", + f"death_mean_{prefix}", + f"betti_mean_{prefix}", + f"landscape_l1_{prefix}", ]: features[key] = 0.0 return features @@ -348,11 +347,16 @@ def _zero_features(self): for dim in self.homology_dimensions: prefix = f"h{dim}" for key in [ - f"total_persistence_{prefix}", f"entropy_{prefix}", - f"n_components_{prefix}", f"amplitude_{prefix}", - f"persistence_std_{prefix}", f"persistence_max_{prefix}", - f"birth_mean_{prefix}", f"death_mean_{prefix}", - f"betti_mean_{prefix}", f"landscape_l1_{prefix}", + f"total_persistence_{prefix}", + f"entropy_{prefix}", + f"n_components_{prefix}", + f"amplitude_{prefix}", + f"persistence_std_{prefix}", + f"persistence_max_{prefix}", + f"birth_mean_{prefix}", + f"death_mean_{prefix}", + f"betti_mean_{prefix}", + f"landscape_l1_{prefix}", ]: features[key] = 0.0 return features @@ -362,10 +366,10 @@ def _process_single_sample(self, X, nbrs, i, n_neighbors): return self.extract_local_tda(X, indices[0]) def extract_global_tda(self, X, n_jobs=-1): - """Extract TDA features with parallel processing. + """Extract TDA features in parallel across samples. - CHANGE 2: X should be the FCGR representation (PCA-compressed, - n_samples × 20) — NOT the biological/standard features. + ``X`` must be the FCGR representation (PCA-compressed, n_samples x 20), + not the biological/standard features. """ n_samples = X.shape[0] n_neighbors = min(self.n_neighbors, n_samples) @@ -377,27 +381,23 @@ def extract_global_tda(self, X, n_jobs=-1): n_jobs = max(1, multiprocessing.cpu_count() + 1 + n_jobs) print(f"Extracting TDA features (over FCGR space) using {n_jobs} CPU cores...") tda_features_list = Parallel(n_jobs=n_jobs, backend="loky", verbose=10)( - delayed(self._process_single_sample)(X, nbrs, i, n_neighbors) - for i in range(n_samples) + delayed(self._process_single_sample)(X, nbrs, i, n_neighbors) for i in range(n_samples) ) tda_df = pd.DataFrame(tda_features_list) return tda_df.values -# --------------------------------------------------------------------------- -# TFDFEPreprocessor (updated: TDA now runs over FCGR space — CHANGE 2) -# --------------------------------------------------------------------------- +# TFDFEPreprocessor -class TFDFEPreprocessor: - """Complete preprocessing pipeline for TF-DFE v4.0. - CHANGE 2 applied here: - The TDA extractor now receives `X_fcgr` (PCA-compressed FCGR matrix) as - input instead of `X_standard` (biological features). This produces - topological features that are structurally independent of the biological - scalar scores, giving the ensemble a genuinely orthogonal signal axis. +class TFDFEPreprocessor: + """End-to-end preprocessing pipeline for the framework. - FIX 3 (retained): FCGR PCA compression 320 → 20 features. + The TDA extractor receives ``X_fcgr`` (the PCA-compressed FCGR matrix) + rather than ``X_standard`` (biological features), so the topological + features are structurally independent of the biological scalar scores and + give the ensemble an orthogonal signal axis. FCGR is PCA-compressed from + 320 to 20 features. """ def __init__(self, use_fcgr=True, use_tda=True, genome_path="hg19.fa"): @@ -410,7 +410,7 @@ def __init__(self, use_fcgr=True, use_tda=True, genome_path="hg19.fa"): self.use_tda = use_tda self.genome_path = genome_path - # FIX 3: PCA to compress FCGR features + # PCA to compress FCGR features self.fcgr_pca = None self.fcgr_pca_n_components = Config.FCGR_PCA_COMPONENTS @@ -441,8 +441,10 @@ def _extract_sequence_context(self, variant_df): sequences = [] failed_extractions = 0 for idx, row in tqdm( - variant_df.iterrows(), total=len(variant_df), - desc="Sequence Extraction", ncols=80, + variant_df.iterrows(), + total=len(variant_df), + desc="Sequence Extraction", + ncols=80, ): try: chrom = str(row["chr"]) @@ -482,7 +484,8 @@ def fit_transform(self, variant_df, target_col="LABEL_PATHOGENIC"): X_num = X_tabular[num_cols].copy() X_num = pd.DataFrame( self.imputer_num.fit_transform(X_num), - columns=num_cols, index=X_tabular.index, + columns=num_cols, + index=X_tabular.index, ) cat_cols = X_tabular.select_dtypes(include=["object"]).columns X_cat = X_tabular[cat_cols].copy() @@ -495,9 +498,7 @@ def fit_transform(self, variant_df, target_col="LABEL_PATHOGENIC"): X_standard = self.scaler.fit_transform(X_combined) print(f"Standard features: {len(self.feature_names)}") - # ---------------------------------------------------------------- - # FCGR encoding + PCA compression (FIX 3 retained) - # ---------------------------------------------------------------- + # FCGR encoding + PCA compression if self.use_fcgr: print("\nGenerating FCGR fractal features...") sequences = self._extract_sequence_context(variant_df) @@ -506,30 +507,28 @@ def fit_transform(self, variant_df, target_col="LABEL_PATHOGENIC"): print(f"FCGR raw features: {X_fcgr_raw.shape[1]}") n_comp = min(self.fcgr_pca_n_components, X_fcgr_raw.shape[1], X_fcgr_raw.shape[0] - 1) - print(f"Compressing FCGR {X_fcgr_raw.shape[1]}D -> {n_comp}D via PCA (FIX 3)...") + print(f"Compressing FCGR {X_fcgr_raw.shape[1]}D -> {n_comp}D via PCA...") self.fcgr_pca = PCA(n_components=n_comp, random_state=Config.RANDOM_STATE) X_fcgr = self.fcgr_pca.fit_transform(X_fcgr_raw) evr = float(np.sum(self.fcgr_pca.explained_variance_ratio_)) * 100 print(f"FCGR PCA features: {X_fcgr.shape[1]} (explained variance: {evr:.1f}%)") self.fcgr_feature_names = [f"fcgr_pca_{i}" for i in range(X_fcgr.shape[1])] - # Store raw FCGR for TDA (used below — CHANGE 2) + # Store raw FCGR for TDA self._X_fcgr_for_tda = X_fcgr else: X_fcgr = np.array([]).reshape(X_standard.shape[0], 0) self.fcgr_feature_names = [] self._X_fcgr_for_tda = None - # ---------------------------------------------------------------- - # CHANGE 2: TDA now over FCGR representation, not biological space - # ---------------------------------------------------------------- + # TDA now over FCGR representation, not biological space if self.use_tda and TDA_AVAILABLE: - print("\nGenerating TDA topological features (over FCGR space — CHANGE 2)...") + print("\nGenerating TDA topological features (over FCGR space)...") # Use PCA-compressed FCGR as input when available; fallback to standard X_for_tda = self._X_fcgr_for_tda if self._X_fcgr_for_tda is not None else X_standard X_tda = self.tda_extractor.extract_global_tda(X_for_tda) n_tda_features = X_tda.shape[1] self.tda_feature_names = [f"tda_{i}" for i in range(n_tda_features)] - print(f"TDA features (expanded): {X_tda.shape[1]} (was 6, now ~20 — CHANGE 2)") + print(f"TDA features (expanded): {X_tda.shape[1]}") else: X_tda = np.array([]).reshape(X_standard.shape[0], 0) self.tda_feature_names = [] @@ -543,18 +542,26 @@ def fit_transform(self, variant_df, target_col="LABEL_PATHOGENIC"): print(f"\nFinal feature space: {X_final.shape[1]} features") print(f" Standard: {len(self.feature_names)}") - print(f" FCGR PCA: {len(self.fcgr_feature_names)} (FIX 3: was {sum(4**k for k in Config.FCGR_K_VALUES)} raw)") - print(f" TDA: {len(self.tda_feature_names)} (CHANGE 2: over FCGR space, was 6)") + print(f" FCGR PCA: {len(self.fcgr_feature_names)}") + print(f" TDA: {len(self.tda_feature_names)}") return X_final, y def transform(self, variant_df, target_col="LABEL_PATHOGENIC"): - """Transform new data using fitted preprocessor (for Code 9 compatibility)""" + """Transform new data with the fitted preprocessor.""" df_clean = variant_df.drop(columns=self.leakage_cols, errors="ignore") y = df_clean[target_col] if target_col in df_clean.columns else None - X_tabular = df_clean.drop(target_col, axis=1, errors="ignore") if target_col in df_clean.columns else df_clean + X_tabular = ( + df_clean.drop(target_col, axis=1, errors="ignore") + if target_col in df_clean.columns + else df_clean + ) - num_cols = [c for c in self.feature_names if c in X_tabular.columns and X_tabular[c].dtype != object] - cat_cols = [c for c in self.feature_names if c in X_tabular.columns and X_tabular[c].dtype == object] + num_cols = [ + c for c in self.feature_names if c in X_tabular.columns and X_tabular[c].dtype != object + ] + cat_cols = [ + c for c in self.feature_names if c in X_tabular.columns and X_tabular[c].dtype == object + ] X_num = X_tabular[num_cols].copy() if num_cols else pd.DataFrame(index=X_tabular.index) X_num = pd.DataFrame( @@ -564,7 +571,9 @@ def transform(self, variant_df, target_col="LABEL_PATHOGENIC"): for col in cat_cols: enc = self.label_encoders[col] known = set(enc.classes_) - X_cat[col] = X_cat[col].astype(str).apply(lambda v: v if v in known else enc.classes_[0]) + X_cat[col] = ( + X_cat[col].astype(str).apply(lambda v: v if v in known else enc.classes_[0]) + ) X_cat[col] = enc.transform(X_cat[col]) X_combined = pd.concat([X_num, X_cat], axis=1) @@ -573,11 +582,13 @@ def transform(self, variant_df, target_col="LABEL_PATHOGENIC"): if self.use_fcgr: sequences = self._extract_sequence_context(variant_df) X_fcgr_raw = self.fcgr_encoder.encode_variants(sequences) - X_fcgr = self.fcgr_pca.transform(X_fcgr_raw) if self.fcgr_pca is not None else X_fcgr_raw + X_fcgr = ( + self.fcgr_pca.transform(X_fcgr_raw) if self.fcgr_pca is not None else X_fcgr_raw + ) else: X_fcgr = np.array([]).reshape(X_standard.shape[0], 0) - # CHANGE 2: TDA over FCGR at inference too + # TDA over FCGR at inference too if self.use_tda and TDA_AVAILABLE: X_for_tda = X_fcgr if X_fcgr.shape[1] > 0 else X_standard X_tda = self.tda_extractor.extract_global_tda(X_for_tda) @@ -589,9 +600,8 @@ def transform(self, variant_df, target_col="LABEL_PATHOGENIC"): return X_final, y -# --------------------------------------------------------------------------- -# Supervised metric learner (unchanged) -# --------------------------------------------------------------------------- +# Supervised metric learner (unchanged) + class SupervisedMetricLearner: """Learns discriminative metric space for KNORA neighbor selection""" @@ -631,9 +641,8 @@ def transform(self, X): return self.scaler.transform(X_combined) -# --------------------------------------------------------------------------- -# MultiViewFeatureManager (updated for expanded TDA count — CHANGE 2) -# --------------------------------------------------------------------------- +# MultiViewFeatureManager + class MultiViewFeatureManager: """Manages multi-view feature representation""" @@ -650,19 +659,16 @@ def __init__(self, n_standard=18, n_fcgr=20, n_tda=20): } -# --------------------------------------------------------------------------- -# CHANGE 1 + CHANGE 2: Pairwise Q-statistic diversity analyser -# --------------------------------------------------------------------------- +# Pairwise Q-statistic diversity analyser + def compute_pairwise_diversity(trained_models, feature_indices, X_val, y_val, output_dir): - """Compute pairwise Q-statistic and disagreement matrix. + """Compute the pairwise Q-statistic and disagreement matrices. - CHANGE 1: Measures inter-learner error correlation on the validation set. - Lower Q-statistic and higher disagreement = more diverse, more useful pool. - Results are saved as: - - pairwise_q_statistic.csv - - pairwise_disagreement.csv - - Fig_Diversity_Heatmap.png + Measures inter-learner error correlation on the validation set; a lower + Q-statistic and higher disagreement indicate a more diverse pool. Results + are saved as pairwise_q_statistic.csv, pairwise_disagreement.csv, and + Fig_Diversity_Heatmap.png. Returns (q_matrix_df, disagreement_matrix_df). """ @@ -675,13 +681,13 @@ def compute_pairwise_diversity(trained_models, feature_indices, X_val, y_val, ou preds[name] = model.predict(X_view) y_arr = np.array(y_val) - q_matrix = np.zeros((n, n)) + q_matrix = np.zeros((n, n)) dis_matrix = np.zeros((n, n)) for i, ni in enumerate(names): for j, nj in enumerate(names): if i == j: - q_matrix[i, j] = 1.0 + q_matrix[i, j] = 1.0 dis_matrix[i, j] = 0.0 continue ci = (preds[ni] == y_arr).astype(int) @@ -692,50 +698,71 @@ def compute_pairwise_diversity(trained_models, feature_indices, X_val, y_val, ou N01 = np.sum((ci == 0) & (cj == 1)) denom = N11 * N00 - N10 * N01 q_denom = N11 * N00 + N10 * N01 - q_matrix[i, j] = denom / (q_denom + 1e-12) + q_matrix[i, j] = denom / (q_denom + 1e-12) dis_matrix[i, j] = (N10 + N01) / (N11 + N00 + N10 + N01 + 1e-12) - q_df = pd.DataFrame(q_matrix, index=names, columns=names) + q_df = pd.DataFrame(q_matrix, index=names, columns=names) dis_df = pd.DataFrame(dis_matrix, index=names, columns=names) q_df.to_csv(output_dir / "pairwise_q_statistic.csv") dis_df.to_csv(output_dir / "pairwise_disagreement.csv") # Plot heatmap fig, axes = plt.subplots(1, 2, figsize=(FIG_DOUBLE_COL_WIDTH, FIG_SINGLE_COL_WIDTH * 1.2)) - sns.heatmap(q_df, ax=axes[0], cmap="RdYlGn_r", vmin=-1, vmax=1, - annot=True, fmt=".2f", linewidths=0.2, annot_kws={"size": 4}) - axes[0].set_title("(A) Q-Statistic\n(−1=diverse, +1=correlated)", - fontweight="bold", loc="left") - sns.heatmap(dis_df, ax=axes[1], cmap="Blues", vmin=0, vmax=0.3, - annot=True, fmt=".2f", linewidths=0.2, annot_kws={"size": 4}) - axes[1].set_title("(B) Pairwise Disagreement Rate\n(higher=more diverse)", - fontweight="bold", loc="left") + sns.heatmap( + q_df, + ax=axes[0], + cmap="RdYlGn_r", + vmin=-1, + vmax=1, + annot=True, + fmt=".2f", + linewidths=0.2, + annot_kws={"size": 4}, + ) + axes[0].set_title("(A) Q-Statistic\n(−1=diverse, +1=correlated)", fontweight="bold", loc="left") + sns.heatmap( + dis_df, + ax=axes[1], + cmap="Blues", + vmin=0, + vmax=0.3, + annot=True, + fmt=".2f", + linewidths=0.2, + annot_kws={"size": 4}, + ) + axes[1].set_title( + "(B) Pairwise Disagreement Rate\n(higher=more diverse)", fontweight="bold", loc="left" + ) plt.tight_layout() - plt.savefig(output_dir / "Fig_Diversity_Heatmap.png", dpi=FIG_DPI, bbox_inches="tight") + plt.savefig(output_dir / "Fig_Diversity_Heatmap.png", dpi=FIG_DPI, bbox_inches="tight") plt.savefig(output_dir / "Fig_Diversity_Heatmap.tiff", dpi=FIG_DPI, bbox_inches="tight") plt.close() - mean_q = np.mean(q_matrix[np.triu_indices(n, k=1)]) + mean_q = np.mean(q_matrix[np.triu_indices(n, k=1)]) mean_dis = np.mean(dis_matrix[np.triu_indices(n, k=1)]) - print(f"\n [Diversity Analysis] Mean Q-statistic: {mean_q:.4f} " - f"(target < 0.50 for diverse pool)") - print(f" [Diversity Analysis] Mean Disagreement: {mean_dis:.4f} " - f"(target > 0.05 for diverse pool)") - print(f" Saved: pairwise_q_statistic.csv, pairwise_disagreement.csv, " - f"Fig_Diversity_Heatmap.png") + print( + f"\n [Diversity Analysis] Mean Q-statistic: {mean_q:.4f} " + f"(target < 0.50 for diverse pool)" + ) + print( + f" [Diversity Analysis] Mean Disagreement: {mean_dis:.4f} " + f"(target > 0.05 for diverse pool)" + ) + print( + f" Saved: pairwise_q_statistic.csv, pairwise_disagreement.csv, " + f"Fig_Diversity_Heatmap.png" + ) return q_df, dis_df -# --------------------------------------------------------------------------- -# DiverseFeatureSubspaceFactory (updated — CHANGE 1: new diverse models) -# --------------------------------------------------------------------------- +# DiverseFeatureSubspaceFactory + class DiverseFeatureSubspaceFactory: """Creates diverse base models with different feature subspaces. - CHANGE 1 (v4.0 — Supervisor Request 1): - ----------------------------------------- - Three new learners with structurally different inductive biases are added + Three learners with structurally different inductive biases are added to break the error-correlated tree pool: 1. MLP_BioTDA — calibrated MLP on standard+TDA block. @@ -754,135 +781,171 @@ class DiverseFeatureSubspaceFactory: where trees succeed, and vice versa — exactly the complementary error pattern needed. - FIX 6 (retained): Isotonic calibration on RF, ExtraTrees, GradientBoosting. - CHANGE 1: MLP also calibrated (isotonic, 3-fold CV). + Isotonic calibration is applied to RandomForest, ExtraTrees, + GradientBoosting, and the MLP (3-fold CV). """ def __init__(self, n_standard=18, n_fcgr=20, n_tda=20): self.view_manager = MultiViewFeatureManager(n_standard, n_fcgr, n_tda) self.models = {} self.feature_indices = {} - self._metric_learner = None # set by train_ensemble() for kNN_SML + self._metric_learner = None # set by train_ensemble() for kNN_SML def create_diverse_ensemble(self, random_state=42): - n_std = self.view_manager.n_standard + n_std = self.view_manager.n_standard n_fcgr = self.view_manager.n_fcgr - n_tda = self.view_manager.n_tda + n_tda = self.view_manager.n_tda models = {} feature_indices = {} - bio_idx = list(range(n_std)) - fcgr_idx = list(range(n_std, n_std + n_fcgr)) - tda_bio_idx = bio_idx + list(range(n_std + n_fcgr, n_std + n_fcgr + n_tda)) + bio_idx = list(range(n_std)) + fcgr_idx = list(range(n_std, n_std + n_fcgr)) + tda_bio_idx = bio_idx + list(range(n_std + n_fcgr, n_std + n_fcgr + n_tda)) bio_fcgr_idx = bio_idx + fcgr_idx - full_idx = list(range(n_std + n_fcgr + n_tda)) + full_idx = list(range(n_std + n_fcgr + n_tda)) if XGBOOST_AVAILABLE: models["Bio_XGBoost"] = XGBClassifier( - n_estimators=200, max_depth=5, learning_rate=0.05, - random_state=random_state, n_jobs=-1, verbosity=0, + n_estimators=200, + max_depth=5, + learning_rate=0.05, + random_state=random_state, + n_jobs=-1, + verbosity=0, ) feature_indices["Bio_XGBoost"] = bio_idx if CATBOOST_AVAILABLE: models["Bio_CatBoost"] = CatBoostClassifier( - iterations=200, depth=5, learning_rate=0.05, - random_state=random_state, verbose=False, + iterations=200, + depth=5, + learning_rate=0.05, + random_state=random_state, + verbose=False, ) feature_indices["Bio_CatBoost"] = bio_idx if LGBM_AVAILABLE: models["BioFCGR_LightGBM"] = LGBMClassifier( - n_estimators=150, num_leaves=63, learning_rate=0.05, - feature_fraction=0.7, random_state=random_state, - n_jobs=-1, verbose=-1, + n_estimators=150, + num_leaves=63, + learning_rate=0.05, + feature_fraction=0.7, + random_state=random_state, + n_jobs=-1, + verbose=-1, ) feature_indices["BioFCGR_LightGBM"] = bio_fcgr_idx - # FIX 6: Calibrated RandomForest + # Calibrated RandomForest _rf_base = RandomForestClassifier( - n_estimators=200, max_depth=15, max_features="sqrt", - min_samples_leaf=5, random_state=random_state + 1, n_jobs=-1, + n_estimators=200, + max_depth=15, + max_features="sqrt", + min_samples_leaf=5, + random_state=random_state + 1, + n_jobs=-1, ) models["BioFCGR_RF"] = CalibratedClassifierCV(_rf_base, method="isotonic", cv=3) feature_indices["BioFCGR_RF"] = bio_fcgr_idx - # FIX 6: Calibrated ExtraTrees + # Calibrated ExtraTrees _et_base = ExtraTreesClassifier( - n_estimators=200, max_depth=12, max_features="sqrt", - random_state=random_state + 2, n_jobs=-1, + n_estimators=200, + max_depth=12, + max_features="sqrt", + random_state=random_state + 2, + n_jobs=-1, ) models["TDABio_ExtraTrees"] = CalibratedClassifierCV(_et_base, method="isotonic", cv=3) feature_indices["TDABio_ExtraTrees"] = tda_bio_idx if XGBOOST_AVAILABLE: models["BioFCGR_XGBoost"] = XGBClassifier( - n_estimators=150, max_depth=4, learning_rate=0.03, - colsample_bytree=0.5, subsample=0.7, - random_state=random_state + 3, n_jobs=-1, verbosity=0, + n_estimators=150, + max_depth=4, + learning_rate=0.03, + colsample_bytree=0.5, + subsample=0.7, + random_state=random_state + 3, + n_jobs=-1, + verbosity=0, ) feature_indices["BioFCGR_XGBoost"] = bio_fcgr_idx - # FIX 6: Calibrated GradientBoosting + # Calibrated GradientBoosting _gb_base = GradientBoostingClassifier( - n_estimators=100, max_depth=4, learning_rate=0.05, - subsample=0.8, max_features="sqrt", random_state=random_state + 4, + n_estimators=100, + max_depth=4, + learning_rate=0.05, + subsample=0.8, + max_features="sqrt", + random_state=random_state + 4, ) models["Full_GradientBoosting"] = CalibratedClassifierCV(_gb_base, method="isotonic", cv=3) feature_indices["Full_GradientBoosting"] = full_idx models["Full_HistGradient"] = HistGradientBoostingClassifier( - max_iter=150, max_depth=6, learning_rate=0.05, + max_iter=150, + max_depth=6, + learning_rate=0.05, random_state=random_state + 5, ) feature_indices["Full_HistGradient"] = full_idx - # ================================================================== - # CHANGE 1 (Supervisor Request 1): Inject algorithmic diversity - # ================================================================== + # Inject algorithmic diversity # 1. Elastic-net logistic regression — linear, biologically interpretable, - # systematically fails on nonlinear regions where trees succeed. + # systematically fails on nonlinear regions where trees succeed. models["ElasticNet_Bio"] = LogisticRegression( - penalty="elasticnet", solver="saga", l1_ratio=0.5, - C=0.5, max_iter=3000, random_state=random_state + 6, n_jobs=-1, + penalty="elasticnet", + solver="saga", + l1_ratio=0.5, + C=0.5, + max_iter=3000, + random_state=random_state + 6, + n_jobs=-1, ) feature_indices["ElasticNet_Bio"] = bio_idx # 2. MLP on bio+TDA — smooth nonlinear boundaries; calibrated (isotonic) - # so its probabilities are directly comparable to the tree models. + # so its probabilities are directly comparable to the tree models. _mlp_base = MLPClassifier( - hidden_layer_sizes=(128, 64, 32), activation="relu", - solver="adam", alpha=1e-4, learning_rate_init=1e-3, - max_iter=300, early_stopping=True, validation_fraction=0.1, + hidden_layer_sizes=(128, 64, 32), + activation="relu", + solver="adam", + alpha=1e-4, + learning_rate_init=1e-3, + max_iter=300, + early_stopping=True, + validation_fraction=0.1, random_state=random_state + 7, ) models["MLP_BioTDA"] = CalibratedClassifierCV(_mlp_base, method="isotonic", cv=3) feature_indices["MLP_BioTDA"] = tda_bio_idx # 3. k-NN in supervised metric space — registered here with a sentinel - # index of -1 (full index); train_ensemble() will replace the feature - # matrix with the pre-projected SML space before fitting. + # index of -1 (full index); train_ensemble will replace the feature + # matrix with the pre-projected SML space before fitting. models["kNN_SML"] = KNeighborsClassifier( - n_neighbors=15, metric="euclidean", n_jobs=-1, + n_neighbors=15, + metric="euclidean", + n_jobs=-1, ) - feature_indices["kNN_SML"] = full_idx # overridden in train_ensemble() - # ================================================================== + feature_indices["kNN_SML"] = full_idx # overridden in train_ensemble() self.models = models self.feature_indices = feature_indices - print(f"\nCreated {len(models)} diverse base models " - f"(CHANGE 1: +MLP, +kNN_SML, +ElasticNet)") + print(f"\nCreated {len(models)} diverse base models " f"(+MLP, +kNN_SML, +ElasticNet)") return models, feature_indices - def train_ensemble(self, X_train, y_train, X_val=None, y_val=None, - consensus_scores=None): + def train_ensemble(self, X_train, y_train, X_val=None, y_val=None, consensus_scores=None): """Train all base models; returns fitted models, val set, and OOF probas. - FIX 4 (retained): OOF predicted probabilities via StratifiedKFold. + Out-of-fold predicted probabilities are produced via StratifiedKFold. - CHANGE 3 (Supervisor Request 3): consensus-margin sample weighting. + Consensus-margin sample weighting: When `consensus_scores` is provided (a 1-D array of CONSENSUS_SCORE values aligned with X_train/y_train), variants whose score is within Config.AMBIGUOUS_MARGIN of 0.5 are down-weighted to @@ -890,7 +953,7 @@ def train_ensemble(self, X_train, y_train, X_val=None, y_val=None, prevents the ambiguous boundary zone from dominating the loss of each base model and reduces correlated FP/FN on hard variants. - CHANGE 1: The kNN_SML model is trained in the 15-D supervised metric + The kNN_SML model is trained in the 15-D supervised metric space produced by SupervisedMetricLearner (fitted on the validation set here to avoid leakage into the training set). """ @@ -898,19 +961,19 @@ def train_ensemble(self, X_train, y_train, X_val=None, y_val=None, if X_val is None or y_val is None: X_train, X_val, y_train, y_val = train_test_split( - X_train, y_train, test_size=0.15, - stratify=y_train, random_state=42, + X_train, + y_train, + test_size=0.15, + stratify=y_train, + random_state=42, ) - # ------------------------------------------------------------------ - # CHANGE 3: Build sample weights from consensus margin - # ------------------------------------------------------------------ + # Build sample weights from consensus margin if consensus_scores is not None: cs = np.array(consensus_scores) # Align consensus to training rows if needed if len(cs) != len(y_train): - print(" [CHANGE 3] WARNING: consensus_scores length mismatch; " - "ignoring sample weights.") + print(" WARNING: consensus_scores length mismatch; " "ignoring sample weights.") sample_weights = None else: margin = np.abs(cs - 0.5) @@ -920,9 +983,11 @@ def train_ensemble(self, X_train, y_train, X_val=None, y_val=None, 1.0, ) n_ambig = int(np.sum(sample_weights < 1.0)) - pct = 100 * n_ambig / len(sample_weights) - print(f" [CHANGE 3] Down-weighted {n_ambig:,} ambiguous variants " - f"({pct:.1f}%) in training set.") + pct = 100 * n_ambig / len(sample_weights) + print( + f" Down-weighted {n_ambig:,} ambiguous variants " + f"({pct:.1f}%) in training set." + ) else: sample_weights = None @@ -933,15 +998,13 @@ def train_ensemble(self, X_train, y_train, X_val=None, y_val=None, oof_probas = np.zeros((n_train, n_classifiers)) skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=Config.RANDOM_STATE) - # ------------------------------------------------------------------ - # CHANGE 1: Pre-fit SML on validation set so kNN_SML can use it. - # SML only needs X_val/y_val; it doesn't see test labels. - # ------------------------------------------------------------------ + # Pre-fit SML on validation set so kNN_SML can use it. + # SML only needs X_val/y_val; it doesn't see test labels. sml = SupervisedMetricLearner(n_components=15) sml.fit(X_val, y_val) self._metric_learner = sml X_train_sml = sml.transform(X_train) - X_val_sml = sml.transform(X_val) + X_val_sml = sml.transform(X_val) # Map kNN_SML feature index to the SML-projected space: # We replace the full feature matrix with X_train_sml / X_val_sml # when fitting/predicting kNN_SML, via a flag checked below. @@ -950,8 +1013,8 @@ def train_ensemble(self, X_train, y_train, X_val=None, y_val=None, for clf_idx, (name, model) in enumerate( tqdm(self.models.items(), desc="Training Models", ncols=80) ): - is_knn_sml = (name == KNN_SML_NAME) - indices = self.feature_indices[name] + is_knn_sml = name == KNN_SML_NAME + indices = self.feature_indices[name] if is_knn_sml: X_train_view = X_train_sml @@ -965,8 +1028,7 @@ def train_ensemble(self, X_train, y_train, X_val=None, y_val=None, fold_model = copy.deepcopy(model) Xf = X_train_view[fold_train_idx] yf = y_train[fold_train_idx] - sw_fold = (sample_weights[fold_train_idx] - if sample_weights is not None else None) + sw_fold = sample_weights[fold_train_idx] if sample_weights is not None else None try: if sw_fold is not None and hasattr(fold_model, "fit"): fold_model.fit(Xf, yf, sample_weight=sw_fold) @@ -999,16 +1061,14 @@ def train_ensemble(self, X_train, y_train, X_val=None, y_val=None, oof_probas[:, clf_idx] = 0.5 # Store SML-projected val set for kNN_SML inference in EnhancedKNORAEnsemble - self._X_val_sml = X_val_sml + self._X_val_sml = X_val_sml self._X_train_sml = X_train_sml self.models = trained_models - self.oof_probas = oof_probas - self.oof_y_train = y_train + self.oof_probas = oof_probas + self.oof_y_train = y_train - # ------------------------------------------------------------------ - # CHANGE 1: Compute diversity matrix right after training - # ------------------------------------------------------------------ + # Compute diversity matrix right after training print("\n[Diversity Analysis] Computing pairwise Q-statistic matrix...") try: # For kNN_SML we need to predict in SML space — wrap for diversity call @@ -1017,10 +1077,10 @@ def train_ensemble(self, X_train, y_train, X_val=None, y_val=None, for nm, mdl in trained_models.items(): if nm == KNN_SML_NAME: # Create a thin wrapper that projects to SML before predict - tmp_models[nm] = _SMLWrappedKNN(mdl, sml, n_total=X_val.shape[1]) + tmp_models[nm] = _SMLWrappedKNN(mdl, sml, n_total=X_val.shape[1]) tmp_indices[nm] = self.feature_indices[nm] else: - tmp_models[nm] = mdl + tmp_models[nm] = mdl tmp_indices[nm] = self.feature_indices[nm] self.q_df, self.dis_df = compute_pairwise_diversity( tmp_models, tmp_indices, X_val, y_val, Config.OUTPUT_DIR @@ -1032,25 +1092,26 @@ def train_ensemble(self, X_train, y_train, X_val=None, y_val=None, return trained_models, X_val, y_val -# --------------------------------------------------------------------------- -# Helper: SML-wrapped kNN for diversity matrix computation (CHANGE 1) -# --------------------------------------------------------------------------- +# Helper: SML-wrapped kNN for diversity matrix computation + class _SMLWrappedKNN: """Thin wrapper so kNN_SML's predict() can be called with the full X.""" + def __init__(self, knn_model, sml, n_total): - self.knn = knn_model - self.sml = sml + self.knn = knn_model + self.sml = sml self.n_total = n_total + def predict(self, X): return self.knn.predict(self.sml.transform(X)) + def predict_proba(self, X): return self.knn.predict_proba(self.sml.transform(X)) -# --------------------------------------------------------------------------- -# FeatureSubspaceWrapper (updated for kNN_SML — CHANGE 1) -# --------------------------------------------------------------------------- +# FeatureSubspaceWrapper + class FeatureSubspaceWrapper: """Wrapper for models with feature subspace. @@ -1080,23 +1141,18 @@ def predict_proba(self, X): return np.column_stack([1 - pred, pred]) -# --------------------------------------------------------------------------- -# EnhancedKNORAEnsemble (bug fix + CHANGE 3 selective prediction) -# --------------------------------------------------------------------------- +# EnhancedKNORAEnsemble + class EnhancedKNORAEnsemble: """Enhanced KNORA-E with Supervised Metric Learning. - v4.0 changes on top of all seven v3 fixes: - - BUG FIX (SHAP, Supervisor): When FIX 7 fallback is active, predict_proba() - was returning `active_preds[:, best_idx]` (correct), but the calling code - in run_final_model() was passing `self.enhanced_knora.best_base_name` - (a STRING) to the SHAP explainer instead of the actual model object. Fixed - in run_final_model() by looking up the model object from base_classifiers. + When the best-base fallback is active, ``predict_proba`` returns + ``active_preds[:, best_idx]``; ``run_final_model`` looks up the model object + from ``base_classifiers`` (rather than passing its name) so the SHAP + explainer receives the fitted estimator. - CHANGE 3 (Supervisor Request 3): `selective_predict()` method added. - Uses KNORA neighbourhood purity (fraction of k-NN neighbours in the + ``selective_predict`` uses KNORA neighbourhood purity (fraction of k-NN neighbours in the validation set that share the predicted class) as a confidence signal. Samples whose purity < Config.SELECTIVE_ABSTAIN_THRESHOLD are flagged as ABSTAIN (-1). This turns TF-DFE's local-competence design into a @@ -1125,15 +1181,15 @@ def __init__( self.metric_learner = None self.classifier_names = list(base_classifiers.keys()) - self.oof_probas = oof_probas + self.oof_probas = oof_probas self.oof_y_train = oof_y_train self.meta_model = None self.stacking_weight = 0.7 self.decision_threshold = 0.5 - self.best_base_mcc = -1.0 - self.best_base_name = None + self.best_base_mcc = -1.0 + self.best_base_name = None self.active_classifier_names = list(base_classifiers.keys()) @@ -1142,7 +1198,7 @@ def fit(self, X_val, y_val): self.X_val = X_val self.y_val = np.array(y_val) - n_samples = len(y_val) + n_samples = len(y_val) n_classifiers = len(self.base_classifiers) print("\n[Step 1/4] Computing per-classifier validation accuracy...") @@ -1152,38 +1208,38 @@ def fit(self, X_val, y_val): val_accuracies[name] = float(np.mean(preds == self.y_val)) print(f" {name:40} - Val Acc: {val_accuracies[name]:.4f}") - best_val_acc = max(val_accuracies.values()) + best_val_acc = max(val_accuracies.values()) prune_threshold = best_val_acc - Config.WEAK_LEARNER_MARGIN self.active_classifier_names = [ - name for name, acc in val_accuracies.items() - if acc >= prune_threshold + name for name, acc in val_accuracies.items() if acc >= prune_threshold ] pruned_names = [ - name for name in self.base_classifiers - if name not in self.active_classifier_names + name for name in self.base_classifiers if name not in self.active_classifier_names ] if pruned_names: - print(f"\n [FIX 2] Pruned {len(pruned_names)} weak learner(s): {pruned_names}") + print(f"\n Pruned {len(pruned_names)} weak learner(s): {pruned_names}") print(f" Active pool: {self.active_classifier_names}") active_classifiers = { - name: clf for name, clf in self.base_classifiers.items() + name: clf + for name, clf in self.base_classifiers.items() if name in self.active_classifier_names } n_active = len(active_classifiers) best_base_mcc_val = -1.0 for name, clf in active_classifiers.items(): - preds = clf.predict(X_val) + preds = clf.predict(X_val) mcc_val = matthews_corrcoef(self.y_val, preds) if mcc_val > best_base_mcc_val: - best_base_mcc_val = mcc_val + best_base_mcc_val = mcc_val self.best_base_name = name self.best_base_mcc = best_base_mcc_val - # Store the actual model object, not just its name (BUG FIX) + # Store the actual model object, not just its name self.best_base_model = active_classifiers[self.best_base_name] - print(f"\n [FIX 7] Best base learner: {self.best_base_name} " - f"(Val MCC={self.best_base_mcc:.4f})") + print( + f"\n Best base learner: {self.best_base_name} " f"(Val MCC={self.best_base_mcc:.4f})" + ) if self.use_metric_learning: print("\n[Step 2/4] Supervised Metric Learning...") @@ -1192,9 +1248,7 @@ def fit(self, X_val, y_val): self.X_val_metric = self.metric_learner.transform(X_val) print(f"\n[Step 3/4] Building k-NN index (k={self.k})...") - self.nn_model = NearestNeighbors( - n_neighbors=self.k, metric="euclidean", n_jobs=-1 - ) + self.nn_model = NearestNeighbors(n_neighbors=self.k, metric="euclidean", n_jobs=-1) self.nn_model.fit(self.X_val_metric) print(f"\n[Step 4/4] Computing oracle matrix ({n_samples} x {n_active})") @@ -1208,15 +1262,14 @@ def fit(self, X_val, y_val): self._fit_stacking_meta_model(X_val, active_classifiers) def _fit_stacking_meta_model(self, X_val, active_classifiers): - n_active = len(active_classifiers) + n_active = len(active_classifiers) active_names = list(active_classifiers.keys()) if self.oof_probas is not None and self.oof_y_train is not None: - print("\n[Stacking / FIX 4] Fitting meta-model on OOF probabilities...") + print("\nStacking: fitting meta-model on OOF probabilities...") all_names = list(self.base_classifiers.keys()) active_col_indices = [ - all_names.index(name) for name in active_names - if name in all_names + all_names.index(name) for name in active_names if name in all_names ] meta_train_X = self.oof_probas[:, active_col_indices] meta_train_y = self.oof_y_train @@ -1234,9 +1287,7 @@ def _fit_stacking_meta_model(self, X_val, active_classifiers): meta_train_X[:, idx] = clf.predict(X_val).astype(float) meta_train_y = self.y_val - self.meta_model = LogisticRegression( - C=1.0, max_iter=2000, random_state=Config.RANDOM_STATE - ) + self.meta_model = LogisticRegression(C=1.0, max_iter=2000, random_state=Config.RANDOM_STATE) self.meta_model.fit(meta_train_X, meta_train_y) # Build validation meta-features for threshold search @@ -1253,26 +1304,24 @@ def _fit_stacking_meta_model(self, X_val, active_classifiers): val_meta_proba = self.meta_model.predict_proba(val_meta_X)[:, 1] - # FIX 5: constrained threshold search - best_mcc = -1.0 - best_thr = 0.5 - thresholds = np.linspace( - Config.THRESHOLD_SEARCH_LOW, Config.THRESHOLD_SEARCH_HIGH, 41 - ) + # constrained threshold search + best_mcc = -1.0 + best_thr = 0.5 + thresholds = np.linspace(Config.THRESHOLD_SEARCH_LOW, Config.THRESHOLD_SEARCH_HIGH, 41) for thr in thresholds: preds = (val_meta_proba >= thr).astype(int) - mcc = matthews_corrcoef(self.y_val, preds) + mcc = matthews_corrcoef(self.y_val, preds) if mcc > best_mcc: best_mcc = mcc best_thr = thr self.decision_threshold = best_thr - # FIX 7: compare ensemble vs best base on validation + # compare ensemble vs best base on validation ensemble_val_mcc = best_mcc if ensemble_val_mcc < self.best_base_mcc: print( - f"\n [FIX 7] WARNING: Ensemble val MCC ({ensemble_val_mcc:.4f}) < " + f"\n WARNING: Ensemble val MCC ({ensemble_val_mcc:.4f}) < " f"best base MCC ({self.best_base_mcc:.4f}). " f"Fallback to '{self.best_base_name}' activated." ) @@ -1281,21 +1330,19 @@ def _fit_stacking_meta_model(self, X_val, active_classifiers): self._use_fallback = False print( - f" Meta-model trained (FIX 1/4/5). Threshold={best_thr:.2f} " + f" Meta-model trained. Threshold={best_thr:.2f} " f"(val MCC={best_mcc:.4f}) | fallback={self._use_fallback}" ) def predict_proba(self, X): - """Return fused probability for the positive class. + """Return the fused probability for the positive class. - BUG FIX (v4.0 — Supervisor): When FIX 7 fallback is active, the old - code correctly returned `active_preds[:, best_idx]`. The bug was - upstream in run_final_model() which passed `self.best_base_name` - (a string) to the SHAP explainer. That is fixed in run_final_model(); - this method remains structurally correct. + When the best-base fallback is active this returns + ``active_preds[:, best_idx]``. ``run_final_model`` passes the model + object (not its name) to the SHAP explainer. """ n_samples = X.shape[0] - n_active = len(self._active_clf_list) + n_active = len(self._active_clf_list) print(f"\nKNORA-E Inference ({n_samples} samples, {n_active} active classifiers)...") @@ -1310,16 +1357,17 @@ def predict_proba(self, X): else: active_preds[:, idx] = clf.predict(X).astype(float) - # FIX 7: fallback to best base model object (not name string) + # fallback to best base model object (not name string) if getattr(self, "_use_fallback", False): best_idx = next( - i for i, (name, _) in enumerate(self._active_clf_list) + i + for i, (name, _) in enumerate(self._active_clf_list) if name == self.best_base_name ) - print(f" [FIX 7] Using fallback: {self.best_base_name}") + print(f" Using fallback: {self.best_base_name}") return active_preds[:, best_idx] - # Meta-model prediction (FIX 1) + # Meta-model prediction if self.meta_model is not None: meta_proba = self.meta_model.predict_proba(active_preds)[:, 1] else: @@ -1335,42 +1383,37 @@ def predict_proba(self, X): knora_predictions = np.zeros(n_samples) sigma = np.mean(distances) + 1e-10 - distance_weights = np.exp(-(distances ** 2) / (2 * sigma ** 2)) - distance_weights = distance_weights / ( - distance_weights.sum(axis=1, keepdims=True) + 1e-10 - ) + distance_weights = np.exp(-(distances**2) / (2 * sigma**2)) + distance_weights = distance_weights / (distance_weights.sum(axis=1, keepdims=True) + 1e-10) for i in range(n_samples): neighbor_indices = indices[i] - sample_weights = distance_weights[i] + sample_weights = distance_weights[i] if np.sum(sample_weights) < 1e-10: sample_weights = np.ones_like(sample_weights) / len(sample_weights) - neighbor_oracles = self.oracle_matrix[neighbor_indices] + neighbor_oracles = self.oracle_matrix[neighbor_indices] competence_scores = np.average(neighbor_oracles, axis=0, weights=sample_weights) - sample_preds = active_preds[i] + sample_preds = active_preds[i] - relative_threshold = max( - self.min_competence, np.percentile(competence_scores, 40) - ) + relative_threshold = max(self.min_competence, np.percentile(competence_scores, 40)) selected_mask = competence_scores >= relative_threshold if not np.any(selected_mask): - top_k_idx = np.argsort(competence_scores)[-max(1, n_active // 2):] + top_k_idx = np.argsort(competence_scores)[-max(1, n_active // 2) :] selected_mask = np.zeros(n_active, dtype=bool) selected_mask[top_k_idx] = True sel_preds = sample_preds[selected_mask] sel_comps = competence_scores[selected_mask] - temp = 2.0 + temp = 2.0 exp_comps = np.exp(sel_comps / temp) - weights = exp_comps / (np.sum(exp_comps) + 1e-10) + weights = exp_comps / (np.sum(exp_comps) + 1e-10) knora_predictions[i] = np.sum(weights * sel_preds) final_predictions = ( - self.stacking_weight * meta_proba - + (1.0 - self.stacking_weight) * knora_predictions + self.stacking_weight * meta_proba + (1.0 - self.stacking_weight) * knora_predictions ) return final_predictions @@ -1380,15 +1423,12 @@ def predict(self, X, threshold=None): proba = self.predict_proba(X) return (proba >= threshold).astype(int) - # ------------------------------------------------------------------ - # CHANGE 3 (Supervisor Request 3): Selective prediction - # ------------------------------------------------------------------ - def selective_predict(self, X, threshold=None, - abstain_threshold=None): + # Selective prediction + def selective_predict(self, X, threshold=None, abstain_threshold=None): """Predict with abstention for low-confidence samples. - CHANGE 3: Uses neighbourhood purity in the supervised metric space - as a confidence signal. Samples whose purity (fraction of k-NN + Uses neighbourhood purity in the supervised metric space as a + confidence signal. Samples whose purity (fraction of k-NN validation neighbours sharing the predicted class) is below `abstain_threshold` (default: Config.SELECTIVE_ABSTAIN_THRESHOLD) are assigned label -1 (ABSTAIN). @@ -1420,63 +1460,64 @@ def selective_predict(self, X, threshold=None, purity_scores = np.zeros(len(X)) for i in range(len(X)): neighbour_labels = self.y_val[indices[i]] - pred_label = y_hard[i] + pred_label = y_hard[i] purity_scores[i] = np.mean(neighbour_labels == pred_label) - y_selective = np.where( - purity_scores >= abstain_threshold, y_hard, -1 - ).astype(int) + y_selective = np.where(purity_scores >= abstain_threshold, y_hard, -1).astype(int) coverage = float(np.mean(y_selective >= 0)) n_abstain = int(np.sum(y_selective == -1)) - print(f"\n [CHANGE 3 / Selective Prediction] " - f"Coverage: {coverage:.3%} Abstentions: {n_abstain:,}") + print( + f"\n [Selective Prediction] " f"Coverage: {coverage:.3%} Abstentions: {n_abstain:,}" + ) return y_selective, purity_scores, coverage -# --------------------------------------------------------------------------- -# EnhancedTFDFEEnsemble (updated for v4 — CHANGE 1, 3 + BUG FIX) -# --------------------------------------------------------------------------- +# EnhancedTFDFEEnsemble + class EnhancedTFDFEEnsemble: - """Enhanced TF-DFE v4.0 Ensemble — All 7 Fixes + Supervisor Changes 1, 2, 3.""" + """Top-level TF-DFE ensemble combining all base learners and the stacker.""" def __init__(self, n_standard=18, n_fcgr=20, n_tda=20): self.n_standard = n_standard self.n_fcgr = n_fcgr self.n_tda = n_tda - self.diverse_factory = None - self.enhanced_knora = None - self.is_fitted = False + self.diverse_factory = None + self.enhanced_knora = None + self.is_fitted = False + + def fit(self, X_train, y_train, X_val=None, y_val=None, consensus_scores_train=None): + """Fit the full ensemble. - def fit(self, X_train, y_train, X_val=None, y_val=None, - consensus_scores_train=None): - """ consensus_scores_train : array-like, shape (n_train,), optional - CONSENSUS_SCORE values for training samples used for CHANGE 3 - sample weighting. Pass None to disable. + CONSENSUS_SCORE values for training samples, used for + consensus-margin sample weighting. Pass None to disable. """ - print("\nEnhanced TF-DFE v4.0 Training (All 7 Fixes + Changes 1, 2, 3)") + print("\nEnhanced TF-DFE Training") n_features = X_train.shape[1] - expected = self.n_standard + self.n_fcgr + self.n_tda + expected = self.n_standard + self.n_fcgr + self.n_tda if n_features != expected: print(f"Adjusting feature counts: {n_features} features") self.n_fcgr = max(0, n_features - self.n_standard - self.n_tda) if self.n_fcgr < 0: self.n_fcgr = 0 - self.n_tda = max(0, n_features - self.n_standard) + self.n_tda = max(0, n_features - self.n_standard) if X_val is None or y_val is None: X_train, X_val, y_train, y_val = train_test_split( - X_train, y_train, test_size=0.15, - stratify=y_train, random_state=Config.RANDOM_STATE, + X_train, + y_train, + test_size=0.15, + stratify=y_train, + random_state=Config.RANDOM_STATE, ) - self.X_val = X_val - self.y_val = np.array(y_val) + self.X_val = X_val + self.y_val = np.array(y_val) y_train_arr = np.array(y_train) self.diverse_factory = DiverseFeatureSubspaceFactory( @@ -1484,17 +1525,20 @@ def fit(self, X_train, y_train, X_val=None, y_val=None, ) self.diverse_factory.create_diverse_ensemble() - # CHANGE 3: pass consensus scores through + # pass consensus scores through trained_models, _, _ = self.diverse_factory.train_ensemble( - X_train, y_train_arr, X_val, y_val, + X_train, + y_train_arr, + X_val, + y_val, consensus_scores=consensus_scores_train, ) - # Wrap models — CHANGE 1: kNN_SML gets the SML projector + # Wrap models —: kNN_SML gets the SML projector sml = self.diverse_factory._metric_learner wrapped_models = {} for name, model in trained_models.items(): - is_sml = (name == "kNN_SML") + is_sml = name == "kNN_SML" wrapped_models[name] = FeatureSubspaceWrapper( model, self.diverse_factory.feature_indices[name], @@ -1528,15 +1572,14 @@ def predict(self, X, threshold=None): return (proba >= threshold).astype(int) def selective_predict(self, X, threshold=None, abstain_threshold=None): - """CHANGE 3: Selective prediction with abstention.""" + """Selective prediction with abstention.""" return self.enhanced_knora.selective_predict( X, threshold=threshold, abstain_threshold=abstain_threshold ) -# --------------------------------------------------------------------------- -# Plotting helpers (unchanged from v3) -# --------------------------------------------------------------------------- +# Plotting helpers (unchanged from ) + def plot_class_distribution(variant_df, output_dir): print(" Generating class distribution plot...") @@ -1546,50 +1589,64 @@ def plot_class_distribution(variant_df, output_dir): bars = ax.bar( ["Benign", "Pathogenic"], [class_counts.get(0, 0), class_counts.get(1, 0)], - color=colors, edgecolor="black", linewidth=0.5, + color=colors, + edgecolor="black", + linewidth=0.5, ) ax.set_ylabel("Count") ax.set_title("(A) Class Distribution", fontweight="bold", loc="left") for bar, count in zip(bars, [class_counts.get(0, 0), class_counts.get(1, 0)]): ax.text( - bar.get_x() + bar.get_width() / 2.0, bar.get_height(), - f"{count:,}", ha="center", va="bottom", fontsize=FIG_FONT_SIZE, + bar.get_x() + bar.get_width() / 2.0, + bar.get_height(), + f"{count:,}", + ha="center", + va="bottom", + fontsize=FIG_FONT_SIZE, ) plt.tight_layout() - plt.savefig(output_dir / "Fig1_class_distribution.png", dpi=FIG_DPI, bbox_inches="tight") + plt.savefig(output_dir / "Fig1_class_distribution.png", dpi=FIG_DPI, bbox_inches="tight") plt.savefig(output_dir / "Fig1_class_distribution.tiff", dpi=FIG_DPI, bbox_inches="tight") plt.close() def plot_roc_pr_curves(y_true, y_pred_proba, model_name, output_dir): print(" Generating ROC and PR curves...") - fig, (ax1, ax2) = plt.subplots( - 1, 2, figsize=(FIG_DOUBLE_COL_WIDTH, FIG_DOUBLE_COL_WIDTH / 2.5) - ) + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(FIG_DOUBLE_COL_WIDTH, FIG_DOUBLE_COL_WIDTH / 2.5)) fpr, tpr, _ = roc_curve(y_true, y_pred_proba) auroc = roc_auc_score(y_true, y_pred_proba) - ax1.plot(fpr, tpr, linewidth=1.5, color="#0077BB", - label=f"{model_name} (AUC = {auroc:.3f})") + ax1.plot(fpr, tpr, linewidth=1.5, color="#0077BB", label=f"{model_name} (AUC = {auroc:.3f})") ax1.plot([0, 1], [0, 1], "k--", linewidth=0.8, alpha=0.7, label="Random") - ax1.set_xlabel("False Positive Rate"); ax1.set_ylabel("True Positive Rate") + ax1.set_xlabel("False Positive Rate") + ax1.set_ylabel("True Positive Rate") ax1.set_title("(A) ROC Curve", fontweight="bold", loc="left") ax1.legend(loc="lower right", frameon=True, fancybox=False, edgecolor="black") - ax1.set_xlim([-0.01, 1.01]); ax1.set_ylim([-0.01, 1.01]) + ax1.set_xlim([-0.01, 1.01]) + ax1.set_ylim([-0.01, 1.01]) ax1.grid(False) precision, recall, _ = precision_recall_curve(y_true, y_pred_proba) auprc = average_precision_score(y_true, y_pred_proba) - ax2.plot(recall, precision, linewidth=1.5, color="#0077BB", - label=f"{model_name} (AP = {auprc:.3f})") + ax2.plot( + recall, precision, linewidth=1.5, color="#0077BB", label=f"{model_name} (AP = {auprc:.3f})" + ) baseline = np.mean(y_true) - ax2.axhline(baseline, color="gray", linestyle="--", linewidth=0.8, alpha=0.7, - label=f"Baseline ({baseline:.3f})") - ax2.set_xlabel("Recall"); ax2.set_ylabel("Precision") + ax2.axhline( + baseline, + color="gray", + linestyle="--", + linewidth=0.8, + alpha=0.7, + label=f"Baseline ({baseline:.3f})", + ) + ax2.set_xlabel("Recall") + ax2.set_ylabel("Precision") ax2.set_title("(B) Precision-Recall Curve", fontweight="bold", loc="left") ax2.legend(loc="upper right", frameon=True, fancybox=False, edgecolor="black") - ax2.set_xlim([-0.01, 1.01]); ax2.set_ylim([-0.01, 1.01]) + ax2.set_xlim([-0.01, 1.01]) + ax2.set_ylim([-0.01, 1.01]) ax2.grid(False) plt.tight_layout() - plt.savefig(output_dir / "Fig1_ROC_PR_curves.png", dpi=FIG_DPI, bbox_inches="tight") + plt.savefig(output_dir / "Fig1_ROC_PR_curves.png", dpi=FIG_DPI, bbox_inches="tight") plt.savefig(output_dir / "Fig1_ROC_PR_curves.tiff", dpi=FIG_DPI, bbox_inches="tight") plt.close() @@ -1607,21 +1664,30 @@ def plot_confusion_matrix(y_true, y_pred, model_name, output_dir): cbar = ax.figure.colorbar(im, ax=ax, fraction=0.046, pad=0.04) cbar.ax.set_ylabel("Proportion", rotation=-90, va="bottom") ax.set( - xticks=np.arange(cm.shape[1]), yticks=np.arange(cm.shape[0]), - xticklabels=classes, yticklabels=classes, - ylabel="True Label", xlabel="Predicted Label", + xticks=np.arange(cm.shape[1]), + yticks=np.arange(cm.shape[0]), + xticklabels=classes, + yticklabels=classes, + ylabel="True Label", + xlabel="Predicted Label", ) ax.set_title(f"{model_name} Confusion Matrix", fontweight="bold", pad=10) thresh = cm_display.max() / 2.0 for i in range(cm.shape[0]): for j in range(cm.shape[1]): text = f"{cm_display[i, j]:.1%}\n({cm[i, j]:,})" - ax.text(j, i, text, ha="center", va="center", - color="white" if cm_display[i, j] > thresh else "black", - fontsize=FIG_FONT_SIZE) + ax.text( + j, + i, + text, + ha="center", + va="center", + color="white" if cm_display[i, j] > thresh else "black", + fontsize=FIG_FONT_SIZE, + ) plt.tight_layout() safe_name = model_name.replace(" ", "_").replace("(", "").replace(")", "") - plt.savefig(output_dir / f"Fig_CM_{safe_name}.png", dpi=FIG_DPI, bbox_inches="tight") + plt.savefig(output_dir / f"Fig_CM_{safe_name}.png", dpi=FIG_DPI, bbox_inches="tight") plt.savefig(output_dir / f"Fig_CM_{safe_name}.tiff", dpi=FIG_DPI, bbox_inches="tight") plt.close() @@ -1629,22 +1695,30 @@ def plot_confusion_matrix(y_true, y_pred, model_name, output_dir): def plot_tsne(X, y, title, output_dir): print(" Generating t-SNE visualization...") np.random.seed(Config.RANDOM_STATE) - idx = np.random.choice(len(X), size=min(2000, len(X)), replace=False) + idx = np.random.choice(len(X), size=min(2000, len(X)), replace=False) X_sub = X[idx] y_sub = y.iloc[idx] if hasattr(y, "iloc") else y[idx] - tsne = TSNE(n_components=2, random_state=Config.RANDOM_STATE, perplexity=30, n_iter=1000) + tsne = TSNE(n_components=2, random_state=Config.RANDOM_STATE, perplexity=30, n_iter=1000) X_embedded = tsne.fit_transform(X_sub) fig, ax = plt.subplots(figsize=(FIG_SINGLE_COL_WIDTH, FIG_SINGLE_COL_WIDTH * 0.9)) colors = {0: "#0077BB", 1: "#CC3311"} for label, name in [(0, "Benign"), (1, "Pathogenic")]: mask = (y_sub == label) if hasattr(y_sub, "values") else (np.array(y_sub) == label) - ax.scatter(X_embedded[mask, 0], X_embedded[mask, 1], - c=colors[label], label=name, alpha=0.6, s=15, edgecolors="none") - ax.set_xlabel("t-SNE Dimension 1"); ax.set_ylabel("t-SNE Dimension 2") + ax.scatter( + X_embedded[mask, 0], + X_embedded[mask, 1], + c=colors[label], + label=name, + alpha=0.6, + s=15, + edgecolors="none", + ) + ax.set_xlabel("t-SNE Dimension 1") + ax.set_ylabel("t-SNE Dimension 2") ax.set_title(title, fontweight="bold") ax.legend(title="Class", loc="best", frameon=True, fancybox=False, edgecolor="black") plt.tight_layout() - plt.savefig(output_dir / "Fig_tSNE_TFDFE.png", dpi=FIG_DPI, bbox_inches="tight") + plt.savefig(output_dir / "Fig_tSNE_TFDFE.png", dpi=FIG_DPI, bbox_inches="tight") plt.savefig(output_dir / "Fig_tSNE_TFDFE.tiff", dpi=FIG_DPI, bbox_inches="tight") plt.close() @@ -1660,7 +1734,7 @@ def perform_shap_analysis(model, X_sample, feature_names, output_dir, model_name else: print(" Using KernelExplainer (slower)...") background = shap.sample(X_sample, 100) - explainer = shap.KernelExplainer(model.predict_proba, background) + explainer = shap.KernelExplainer(model.predict_proba, background) shap_values = explainer.shap_values(X_sample) if isinstance(shap_values, list) and len(shap_values) == 2: shap_values = shap_values[1] @@ -1668,20 +1742,29 @@ def perform_shap_analysis(model, X_sample, feature_names, output_dir, model_name shap_values = shap_values[:, :, 1] short_names = [name[:25] + "..." if len(name) > 28 else name for name in feature_names] fig, ax = plt.subplots(figsize=(FIG_SINGLE_COL_WIDTH, FIG_SINGLE_COL_WIDTH * 1.2)) - shap.summary_plot(shap_values, X_sample, feature_names=short_names, plot_type="bar", - show=False, max_display=20, color="#0077BB") + shap.summary_plot( + shap_values, + X_sample, + feature_names=short_names, + plot_type="bar", + show=False, + max_display=20, + color="#0077BB", + ) plt.title(f"Feature Importance ({model_name})", fontweight="bold", loc="left") plt.xlabel("Mean |SHAP Value|") plt.tight_layout() - plt.savefig(output_dir / "Fig_SHAP_summary.png", dpi=FIG_DPI, bbox_inches="tight") + plt.savefig(output_dir / "Fig_SHAP_summary.png", dpi=FIG_DPI, bbox_inches="tight") plt.savefig(output_dir / "Fig_SHAP_summary.tiff", dpi=FIG_DPI, bbox_inches="tight") plt.close() fig, ax = plt.subplots(figsize=(FIG_DOUBLE_COL_WIDTH * 0.7, FIG_SINGLE_COL_WIDTH * 1.2)) - shap.summary_plot(shap_values, X_sample, feature_names=short_names, show=False, max_display=20) + shap.summary_plot( + shap_values, X_sample, feature_names=short_names, show=False, max_display=20 + ) plt.title("SHAP Value Distribution", fontweight="bold", loc="left") plt.xlabel("SHAP Value (Impact on Pathogenicity)") plt.tight_layout() - plt.savefig(output_dir / "Fig_SHAP_beeswarm.png", dpi=FIG_DPI, bbox_inches="tight") + plt.savefig(output_dir / "Fig_SHAP_beeswarm.png", dpi=FIG_DPI, bbox_inches="tight") plt.savefig(output_dir / "Fig_SHAP_beeswarm.tiff", dpi=FIG_DPI, bbox_inches="tight") plt.close() print(" SHAP analysis complete") @@ -1691,9 +1774,8 @@ def perform_shap_analysis(model, X_sample, feature_names, output_dir, model_name return None -# --------------------------------------------------------------------------- -# TFDFEvaluator (unchanged from v3, aggregated SHAP bug-fix in run_final) -# --------------------------------------------------------------------------- +# TFDFEvaluator + class TFDFEvaluator: @staticmethod @@ -1703,51 +1785,94 @@ def calculate_ece(y_true, y_pred_proba, n_bins=10): bin_data = [] for bl, bu in zip(bin_boundaries[:-1], bin_boundaries[1:]): in_bin = (y_pred_proba > bl) & (y_pred_proba <= bu) - prop = np.mean(in_bin) + prop = np.mean(in_bin) if prop > 0: avg_conf = np.mean(y_pred_proba[in_bin]) - avg_acc = np.mean(y_true[in_bin]) + avg_acc = np.mean(y_true[in_bin]) ece += np.abs(avg_acc - avg_conf) * prop - bin_data.append({"bin_lower": bl, "bin_upper": bu, - "avg_confidence": avg_conf, "avg_accuracy": avg_acc, - "count": int(np.sum(in_bin)), "proportion": prop, - "calibration_error": np.abs(avg_acc - avg_conf)}) + bin_data.append( + { + "bin_lower": bl, + "bin_upper": bu, + "avg_confidence": avg_conf, + "avg_accuracy": avg_acc, + "count": int(np.sum(in_bin)), + "proportion": prop, + "calibration_error": np.abs(avg_acc - avg_conf), + } + ) else: - bin_data.append({"bin_lower": bl, "bin_upper": bu, - "avg_confidence": np.nan, "avg_accuracy": np.nan, - "count": 0, "proportion": 0, "calibration_error": np.nan}) + bin_data.append( + { + "bin_lower": bl, + "bin_upper": bu, + "avg_confidence": np.nan, + "avg_accuracy": np.nan, + "count": 0, + "proportion": 0, + "calibration_error": np.nan, + } + ) return ece, pd.DataFrame(bin_data) @staticmethod def plot_calibration_curve(y_true, y_pred_proba, output_dir, model_name="TF-DFE", n_bins=10): print(" Generating calibration curve...") - prob_true, prob_pred = calibration_curve(y_true, y_pred_proba, n_bins=n_bins, strategy="uniform") - ece, bin_data = TFDFEvaluator.calculate_ece(np.array(y_true), np.array(y_pred_proba), n_bins) + prob_true, prob_pred = calibration_curve( + y_true, y_pred_proba, n_bins=n_bins, strategy="uniform" + ) + ece, bin_data = TFDFEvaluator.calculate_ece( + np.array(y_true), np.array(y_pred_proba), n_bins + ) brier = brier_score_loss(y_true, y_pred_proba) fig = plt.figure(figsize=(FIG_SINGLE_COL_WIDTH, FIG_SINGLE_COL_WIDTH * 1.2)) - gs = fig.add_gridspec(2, 1, height_ratios=[3, 1], hspace=0.05) + gs = fig.add_gridspec(2, 1, height_ratios=[3, 1], hspace=0.05) ax1 = fig.add_subplot(gs[0]) ax1.plot([0, 1], [0, 1], "k--", linewidth=1, label="Perfect Calibration") - ax1.plot(prob_pred, prob_true, "s-", color="#0077BB", linewidth=1.5, markersize=6, - label=f"{model_name}") + ax1.plot( + prob_pred, + prob_true, + "s-", + color="#0077BB", + linewidth=1.5, + markersize=6, + label=f"{model_name}", + ) ax1.fill_between(prob_pred, prob_pred, prob_true, alpha=0.2, color="#CC3311") ax1.set_ylabel("Observed Probability") - ax1.set_xlim([-0.02, 1.02]); ax1.set_ylim([-0.02, 1.02]) + ax1.set_xlim([-0.02, 1.02]) + ax1.set_ylim([-0.02, 1.02]) ax1.set_title("Calibration Plot", fontweight="bold", loc="left") ax1.legend(loc="upper left", frameon=True, fancybox=False, edgecolor="black") ax1.grid(False) textstr = f"ECE = {ece:.4f}\nBrier = {brier:.4f}" props = dict(boxstyle="round", facecolor="white", edgecolor="gray", alpha=0.9) - ax1.text(0.95, 0.05, textstr, transform=ax1.transAxes, fontsize=FIG_FONT_SIZE, - va="bottom", ha="right", bbox=props) + ax1.text( + 0.95, + 0.05, + textstr, + transform=ax1.transAxes, + fontsize=FIG_FONT_SIZE, + va="bottom", + ha="right", + bbox=props, + ) ax1.set_xticklabels([]) ax2 = fig.add_subplot(gs[1], sharex=ax1) - ax2.hist(y_pred_proba, bins=n_bins, range=(0, 1), color="#0077BB", - edgecolor="black", alpha=0.7, linewidth=0.5) - ax2.set_xlabel("Predicted Probability"); ax2.set_ylabel("Count") + ax2.hist( + y_pred_proba, + bins=n_bins, + range=(0, 1), + color="#0077BB", + edgecolor="black", + alpha=0.7, + linewidth=0.5, + ) + ax2.set_xlabel("Predicted Probability") + ax2.set_ylabel("Count") ax2.set_xlim([-0.02, 1.02]) plt.tight_layout() - plt.savefig(output_dir / "Fig_Calibration_Curve.png", dpi=FIG_DPI, bbox_inches="tight") + plt.savefig(output_dir / "Fig_Calibration_Curve.png", dpi=FIG_DPI, bbox_inches="tight") plt.savefig(output_dir / "Fig_Calibration_Curve.tiff", dpi=FIG_DPI, bbox_inches="tight") plt.close() bin_data.to_csv(output_dir / "calibration_bins.csv", index=False) @@ -1758,8 +1883,15 @@ def plot_calibration_curve(y_true, y_pred_proba, output_dir, model_name="TF-DFE" @staticmethod def perform_repeated_cv( - X, y, model_factory_fn, n_splits=5, n_repeats=2, - n_standard=10, n_fcgr=20, n_tda=20, output_dir=None, + X, + y, + model_factory_fn, + n_splits=5, + n_repeats=2, + n_standard=10, + n_fcgr=20, + n_tda=20, + output_dir=None, ): print(f"\n Running {n_repeats}x{n_splits}-fold Repeated Stratified CV...") rskf = RepeatedStratifiedKFold( @@ -1769,50 +1901,63 @@ def perform_repeated_cv( fold_idx = 0 y_arr = np.array(y) if hasattr(y, "values") else y for train_idx, test_idx in tqdm( - rskf.split(X, y_arr), total=n_splits * n_repeats, - desc="CV Progress", ncols=80, + rskf.split(X, y_arr), + total=n_splits * n_repeats, + desc="CV Progress", + ncols=80, ): fold_idx += 1 X_tr, X_te = X[train_idx], X[test_idx] y_tr, y_te = y_arr[train_idx], y_arr[test_idx] model = model_factory_fn(n_standard, n_fcgr, n_tda) model.fit(X_tr, y_tr) - y_pred_cv = model.predict(X_te) + y_pred_cv = model.predict(X_te) y_pred_proba_cv = model.predict_proba(X_te) fold_metrics = { "repeat": (fold_idx - 1) // n_splits + 1, - "fold": (fold_idx - 1) % n_splits + 1, - "mcc": matthews_corrcoef(y_te, y_pred_cv), - "auprc": average_precision_score(y_te, y_pred_proba_cv), - "auroc": roc_auc_score(y_te, y_pred_proba_cv), - "f1": f1_score(y_te, y_pred_cv, zero_division=0), - "precision":precision_score(y_te, y_pred_cv, zero_division=0), - "recall": recall_score(y_te, y_pred_cv, zero_division=0), + "fold": (fold_idx - 1) % n_splits + 1, + "mcc": matthews_corrcoef(y_te, y_pred_cv), + "auprc": average_precision_score(y_te, y_pred_proba_cv), + "auroc": roc_auc_score(y_te, y_pred_proba_cv), + "f1": f1_score(y_te, y_pred_cv, zero_division=0), + "precision": precision_score(y_te, y_pred_cv, zero_division=0), + "recall": recall_score(y_te, y_pred_cv, zero_division=0), "accuracy": accuracy_score(y_te, y_pred_cv), - "brier": brier_score_loss(y_te, y_pred_proba_cv), + "brier": brier_score_loss(y_te, y_pred_proba_cv), } cv_results.append(fold_metrics) cv_df = pd.DataFrame(cv_results) summary_stats = {} for metric in ["mcc", "auprc", "auroc", "f1", "precision", "recall", "accuracy", "brier"]: - values = cv_df[metric].values + values = cv_df[metric].values mean_val = np.mean(values) - std_val = np.std(values, ddof=1) - n = len(values) - ci_95 = stats.t.ppf(0.975, n - 1) * (std_val / np.sqrt(n)) + std_val = np.std(values, ddof=1) + n = len(values) + ci_95 = stats.t.ppf(0.975, n - 1) * (std_val / np.sqrt(n)) summary_stats[metric] = { - "mean": mean_val, "std": std_val, "ci_95": ci_95, - "ci_lower": mean_val - ci_95, "ci_upper": mean_val + ci_95, + "mean": mean_val, + "std": std_val, + "ci_95": ci_95, + "ci_lower": mean_val - ci_95, + "ci_upper": mean_val + ci_95, "formatted": f"{mean_val:.4f} ± {ci_95:.4f}", } if output_dir: cv_df.to_csv(output_dir / "cv_metrics.csv", index=False) - summary_df = pd.DataFrame([ - {"metric": m, "mean": s["mean"], "std": s["std"], - "ci_95": s["ci_95"], "ci_lower": s["ci_lower"], - "ci_upper": s["ci_upper"], "formatted": s["formatted"]} - for m, s in summary_stats.items() - ]) + summary_df = pd.DataFrame( + [ + { + "metric": m, + "mean": s["mean"], + "std": s["std"], + "ci_95": s["ci_95"], + "ci_lower": s["ci_lower"], + "ci_upper": s["ci_upper"], + "formatted": s["formatted"], + } + for m, s in summary_stats.items() + ] + ) summary_df.to_csv(output_dir / "cv_summary_stats.csv", index=False) print(f"\n CV Results Summary (Mean ± 95% CI):") for m in ["mcc", "auprc", "auroc", "f1"]: @@ -1823,7 +1968,8 @@ def perform_repeated_cv( def plot_cv_distribution(cv_df, output_dir, metrics=["mcc", "auprc"]): print(" Generating CV distribution plot...") fig, axes = plt.subplots( - 1, len(metrics), + 1, + len(metrics), figsize=(FIG_SINGLE_COL_WIDTH * len(metrics), FIG_SINGLE_COL_WIDTH * 0.8), ) if len(metrics) == 1: @@ -1834,15 +1980,25 @@ def plot_cv_distribution(cv_df, output_dir, metrics=["mcc", "auprc"]): ax = axes[i] parts = ax.violinplot(cv_df[metric], positions=[1], showmeans=True, showextrema=True) for pc in parts["bodies"]: - pc.set_facecolor(colors[i % len(colors)]); pc.set_alpha(0.6) - bp = ax.boxplot(cv_df[metric], positions=[1], widths=0.15, - patch_artist=True, showfliers=True) - bp["boxes"][0].set_facecolor("white"); bp["boxes"][0].set_alpha(0.8) + pc.set_facecolor(colors[i % len(colors)]) + pc.set_alpha(0.6) + bp = ax.boxplot( + cv_df[metric], positions=[1], widths=0.15, patch_artist=True, showfliers=True + ) + bp["boxes"][0].set_facecolor("white") + bp["boxes"][0].set_alpha(0.8) jitter = np.random.uniform(-0.05, 0.05, len(cv_df)) - ax.scatter(np.ones(len(cv_df)) + jitter, cv_df[metric], alpha=0.5, s=20, - color=colors[i % len(colors)], edgecolor="black", linewidth=0.5) + ax.scatter( + np.ones(len(cv_df)) + jitter, + cv_df[metric], + alpha=0.5, + s=20, + color=colors[i % len(colors)], + edgecolor="black", + linewidth=0.5, + ) mean_val = cv_df[metric].mean() - std_val = cv_df[metric].std() + std_val = cv_df[metric].std() n = len(cv_df) ci_95 = stats.t.ppf(0.975, n - 1) * (std_val / np.sqrt(n)) ax.axhline(mean_val, color="red", linestyle="--", linewidth=1, alpha=0.7) @@ -1851,19 +2007,25 @@ def plot_cv_distribution(cv_df, output_dir, metrics=["mcc", "auprc"]): f"{metric_names.get(metric, metric.upper())}\n{mean_val:.4f} ± {ci_95:.4f}", fontweight="bold", ) - ax.set_xticks([]); ax.grid(False) + ax.set_xticks([]) + ax.grid(False) plt.suptitle("Cross-Validation Metric Distributions", fontweight="bold", y=1.02) plt.tight_layout() - plt.savefig(output_dir / "Fig_CV_Distribution.png", dpi=FIG_DPI, bbox_inches="tight") + plt.savefig(output_dir / "Fig_CV_Distribution.png", dpi=FIG_DPI, bbox_inches="tight") plt.savefig(output_dir / "Fig_CV_Distribution.tiff", dpi=FIG_DPI, bbox_inches="tight") plt.close() @staticmethod - def perform_mcnemar_test(y_true, y_pred_model1, y_pred_model2, - model1_name="Enhanced Model", model2_name="Baseline", - output_dir=None): + def perform_mcnemar_test( + y_true, + y_pred_model1, + y_pred_model2, + model1_name="Enhanced Model", + model2_name="Baseline", + output_dir=None, + ): print(f"\n Performing McNemar's Test: {model1_name} vs {model2_name}") - y_true = np.array(y_true) + y_true = np.array(y_true) y_pred_model1 = np.array(y_pred_model1) y_pred_model2 = np.array(y_pred_model2) correct_1 = y_pred_model1 == y_true @@ -1878,33 +2040,43 @@ def perform_mcnemar_test(y_true, y_pred_model1, y_pred_model2, interpretation = "No disagreements" else: statistic = (abs(b - c) - 1) ** 2 / (b + c) - p_value = 1 - stats.chi2.cdf(statistic, df=1) + p_value = 1 - stats.chi2.cdf(statistic, df=1) if b + c < 25: p_value = 2 * min( stats.binom.cdf(min(b, c), b + c, 0.5), 1 - stats.binom.cdf(max(b, c) - 1, b + c, 0.5), ) - if p_value < 0.001: sig_symbol, interpretation = "***", "Highly significant (p<0.001)" - elif p_value < 0.01: sig_symbol, interpretation = "**", "Very significant (p<0.01)" - elif p_value < 0.05: sig_symbol, interpretation = "*", "Significant (p<0.05)" - else: sig_symbol, interpretation = "ns", "No significant difference" - acc1 = np.mean(correct_1); acc2 = np.mean(correct_2) + if p_value < 0.001: + sig_symbol, interpretation = "***", "Highly significant (p<0.001)" + elif p_value < 0.01: + sig_symbol, interpretation = "**", "Very significant (p<0.01)" + elif p_value < 0.05: + sig_symbol, interpretation = "*", "Significant (p<0.05)" + else: + sig_symbol, interpretation = "ns", "No significant difference" + acc1 = np.mean(correct_1) + acc2 = np.mean(correct_2) result = { - "model1": model1_name, "model2": model2_name, - "model1_accuracy": acc1, "model2_accuracy": acc2, + "model1": model1_name, + "model2": model2_name, + "model1_accuracy": acc1, + "model2_accuracy": acc2, "accuracy_difference": acc1 - acc2, - "both_correct": a, "model1_only_correct": b, - "model2_only_correct": c, "both_wrong": d, - "statistic": statistic, "p_value": p_value, - "significance": sig_symbol, "interpretation": interpretation, + "both_correct": a, + "model1_only_correct": b, + "model2_only_correct": c, + "both_wrong": d, + "statistic": statistic, + "p_value": p_value, + "significance": sig_symbol, + "interpretation": interpretation, } print(f"\n McNemar Statistic: {statistic:.4f}") print(f" P-value: {p_value:.6f} ({sig_symbol})") print(f" {interpretation}") print(f" Accuracy improvement: {acc1 - acc2:+.4f}") if output_dir: - pd.DataFrame([result]).to_csv( - output_dir / "significance_test_mcnemar.csv", index=False) + pd.DataFrame([result]).to_csv(output_dir / "significance_test_mcnemar.csv", index=False) return result @staticmethod @@ -1912,9 +2084,10 @@ def generate_misclassification_report( variant_df, y_true, y_pred, y_pred_proba, test_indices, output_dir ): print(" Generating misclassification report...") - y_true = np.array(y_true); y_pred = np.array(y_pred) + y_true = np.array(y_true) + y_pred = np.array(y_pred) y_pred_proba = np.array(y_pred_proba) - misclass_mask = y_true != y_pred + misclass_mask = y_true != y_pred misclass_indices = np.where(misclass_mask)[0] original_indices = test_indices[misclass_indices] coord_cols = [c for c in ["chr", "pos", "ref", "alt"] if c in variant_df.columns] @@ -1922,8 +2095,10 @@ def generate_misclassification_report( for i, orig_idx in enumerate(original_indices): li = misclass_indices[i] rec = { - "test_index": li, "original_index": orig_idx, - "y_true": int(y_true[li]), "y_pred": int(y_pred[li]), + "test_index": li, + "original_index": orig_idx, + "y_true": int(y_true[li]), + "y_pred": int(y_pred[li]), "y_pred_proba": float(y_pred_proba[li]), "error_type": "FP" if y_pred[li] == 1 else "FN", "confidence": abs(y_pred_proba[li] - 0.5) * 2, @@ -1940,7 +2115,9 @@ def generate_misclassification_report( return df @staticmethod - def aggregate_shap_by_category(shap_values, feature_names, n_standard, n_fcgr, n_tda, output_dir): + def aggregate_shap_by_category( + shap_values, feature_names, n_standard, n_fcgr, n_tda, output_dir + ): print(" Aggregating SHAP values by feature category...") mean_abs_shap = np.mean(np.abs(shap_values), axis=0) category_shap = {"Standard/Biological": [], "Fractal (FCGR)": [], "Topological (TDA)": []} @@ -1953,130 +2130,208 @@ def aggregate_shap_by_category(shap_values, feature_names, n_standard, n_fcgr, n else: category = "Topological (TDA)" category_shap[category].append(mean_abs_shap[i]) - feature_category_map.append({"feature_name": fname, "category": category, - "mean_abs_shap": mean_abs_shap[i]}) + feature_category_map.append( + {"feature_name": fname, "category": category, "mean_abs_shap": mean_abs_shap[i]} + ) aggregated_results = [] for category, values in category_shap.items(): if values: - aggregated_results.append({ - "category": category, "n_features": len(values), - "mean_abs_shap": np.mean(values), "sum_abs_shap": np.sum(values), - "std_abs_shap": np.std(values), "max_abs_shap": np.max(values), - "relative_importance": np.sum(values) / np.sum(mean_abs_shap) * 100, - }) + aggregated_results.append( + { + "category": category, + "n_features": len(values), + "mean_abs_shap": np.mean(values), + "sum_abs_shap": np.sum(values), + "std_abs_shap": np.std(values), + "max_abs_shap": np.max(values), + "relative_importance": np.sum(values) / np.sum(mean_abs_shap) * 100, + } + ) aggregated_df = pd.DataFrame(aggregated_results).sort_values( - "mean_abs_shap", ascending=False) + "mean_abs_shap", ascending=False + ) feature_df = pd.DataFrame(feature_category_map) feature_df.to_csv(output_dir / "feature_importance_detailed.csv", index=False) aggregated_df.to_csv(output_dir / "feature_importance_aggregated.csv", index=False) print("\n Category-Level Importance:") for _, row in aggregated_df.iterrows(): - print(f" {row['category']:25} - Mean |SHAP|: {row['mean_abs_shap']:.6f} " - f"({row['relative_importance']:.1f}%)") + print( + f" {row['category']:25} - Mean |SHAP|: {row['mean_abs_shap']:.6f} " + f"({row['relative_importance']:.1f}%)" + ) return aggregated_df, feature_df @staticmethod def plot_feature_importance_stacked(aggregated_df, output_dir): print(" Generating feature importance stacked bar chart...") fig, (ax1, ax2) = plt.subplots( - 1, 2, figsize=(FIG_DOUBLE_COL_WIDTH, FIG_SINGLE_COL_WIDTH * 0.7)) - colors = {"Standard/Biological": "#0077BB", - "Fractal (FCGR)": "#EE7733", - "Topological (TDA)": "#009988"} + 1, 2, figsize=(FIG_DOUBLE_COL_WIDTH, FIG_SINGLE_COL_WIDTH * 0.7) + ) + colors = { + "Standard/Biological": "#0077BB", + "Fractal (FCGR)": "#EE7733", + "Topological (TDA)": "#009988", + } categories = aggregated_df["category"].tolist() bar_colors = [colors.get(c, "#888888") for c in categories] - bars = ax1.barh(categories, aggregated_df["mean_abs_shap"], - color=bar_colors, edgecolor="black", linewidth=0.5) + bars = ax1.barh( + categories, + aggregated_df["mean_abs_shap"], + color=bar_colors, + edgecolor="black", + linewidth=0.5, + ) ax1.set_xlabel("Mean |SHAP Value|") ax1.set_title("(A) Feature Category Importance", fontweight="bold", loc="left") ax1.grid(False) for bar, val in zip(bars, aggregated_df["mean_abs_shap"]): - ax1.text(val + 0.001, bar.get_y() + bar.get_height() / 2, - f"{val:.4f}", va="center", fontsize=FIG_FONT_SIZE - 1) + ax1.text( + val + 0.001, + bar.get_y() + bar.get_height() / 2, + f"{val:.4f}", + va="center", + fontsize=FIG_FONT_SIZE - 1, + ) wedges, texts, autotexts = ax2.pie( - aggregated_df["relative_importance"], labels=categories, - colors=bar_colors, autopct="%1.1f%%", startangle=90, + aggregated_df["relative_importance"], + labels=categories, + colors=bar_colors, + autopct="%1.1f%%", + startangle=90, explode=[0.02] * len(categories), wedgeprops=dict(linewidth=1, edgecolor="white"), ) ax2.set_title("(B) Relative Contribution", fontweight="bold") for at in autotexts: - at.set_fontsize(FIG_FONT_SIZE); at.set_fontweight("bold") + at.set_fontsize(FIG_FONT_SIZE) + at.set_fontweight("bold") plt.tight_layout() - plt.savefig(output_dir / "Fig_Feature_Importance_Stacked.png", - dpi=FIG_DPI, bbox_inches="tight") - plt.savefig(output_dir / "Fig_Feature_Importance_Stacked.tiff", - dpi=FIG_DPI, bbox_inches="tight") + plt.savefig( + output_dir / "Fig_Feature_Importance_Stacked.png", dpi=FIG_DPI, bbox_inches="tight" + ) + plt.savefig( + output_dir / "Fig_Feature_Importance_Stacked.tiff", dpi=FIG_DPI, bbox_inches="tight" + ) plt.close() @staticmethod def generate_learning_curve( - X, y, model_factory_fn, n_standard, n_fcgr, n_tda, output_dir, + X, + y, + model_factory_fn, + n_standard, + n_fcgr, + n_tda, + output_dir, train_sizes=[0.1, 0.3, 0.5, 0.7, 0.9], ): print("\n Generating Learning Curve...") print(f" Training sizes: {[f'{s*100:.0f}%' for s in train_sizes]}") - y_arr = np.array(y) if hasattr(y, "values") else y + y_arr = np.array(y) if hasattr(y, "values") else y results = [] for train_size in tqdm(train_sizes, desc="Learning Curve", ncols=80): if train_size < 1.0: X_tr, X_te, y_tr, y_te = train_test_split( - X, y_arr, train_size=train_size, stratify=y_arr, + X, + y_arr, + train_size=train_size, + stratify=y_arr, random_state=Config.RANDOM_STATE, ) else: X_tr, X_te, y_tr, y_te = train_test_split( - X, y_arr, test_size=0.2, stratify=y_arr, + X, + y_arr, + test_size=0.2, + stratify=y_arr, random_state=Config.RANDOM_STATE, ) model = model_factory_fn(n_standard, n_fcgr, n_tda) model.fit(X_tr, y_tr) - y_pred = model.predict(X_te) + y_pred = model.predict(X_te) y_pred_proba = model.predict_proba(X_te) - results.append({ - "train_size_pct": train_size * 100, - "n_train_samples": len(X_tr), "n_test_samples": len(X_te), - "mcc": matthews_corrcoef(y_te, y_pred), - "auprc": average_precision_score(y_te, y_pred_proba), - "auroc": roc_auc_score(y_te, y_pred_proba), - "f1": f1_score(y_te, y_pred, zero_division=0), - "accuracy": accuracy_score(y_te, y_pred), - }) + results.append( + { + "train_size_pct": train_size * 100, + "n_train_samples": len(X_tr), + "n_test_samples": len(X_te), + "mcc": matthews_corrcoef(y_te, y_pred), + "auprc": average_precision_score(y_te, y_pred_proba), + "auroc": roc_auc_score(y_te, y_pred_proba), + "f1": f1_score(y_te, y_pred, zero_division=0), + "accuracy": accuracy_score(y_te, y_pred), + } + ) results_df = pd.DataFrame(results) fig, axes = plt.subplots(1, 2, figsize=(FIG_DOUBLE_COL_WIDTH, FIG_SINGLE_COL_WIDTH * 0.7)) ax1 = axes[0] - ax1.plot(results_df["train_size_pct"], results_df["mcc"], "o-", - color="#0077BB", linewidth=1.5, markersize=6, label="MCC") - ax1.plot(results_df["train_size_pct"], results_df["f1"], "s--", - color="#EE7733", linewidth=1.5, markersize=6, label="F1-Score") - ax1.set_xlabel("Training Data Size (%)"); ax1.set_ylabel("Score") + ax1.plot( + results_df["train_size_pct"], + results_df["mcc"], + "o-", + color="#0077BB", + linewidth=1.5, + markersize=6, + label="MCC", + ) + ax1.plot( + results_df["train_size_pct"], + results_df["f1"], + "s--", + color="#EE7733", + linewidth=1.5, + markersize=6, + label="F1-Score", + ) + ax1.set_xlabel("Training Data Size (%)") + ax1.set_ylabel("Score") ax1.set_title("(A) Learning Curve — MCC & F1", fontweight="bold", loc="left") - ax1.legend(); ax1.grid(False) - ax1.set_xlim([0, 100]); ax1.set_ylim([0, 1.05]) + ax1.legend() + ax1.grid(False) + ax1.set_xlim([0, 100]) + ax1.set_ylim([0, 1.05]) ax2 = axes[1] - ax2.plot(results_df["train_size_pct"], results_df["auprc"], "o-", - color="#009988", linewidth=1.5, markersize=6, label="AUPRC") - ax2.plot(results_df["train_size_pct"], results_df["auroc"], "s--", - color="#CC3311", linewidth=1.5, markersize=6, label="AUROC") - ax2.set_xlabel("Training Data Size (%)"); ax2.set_ylabel("Score") + ax2.plot( + results_df["train_size_pct"], + results_df["auprc"], + "o-", + color="#009988", + linewidth=1.5, + markersize=6, + label="AUPRC", + ) + ax2.plot( + results_df["train_size_pct"], + results_df["auroc"], + "s--", + color="#CC3311", + linewidth=1.5, + markersize=6, + label="AUROC", + ) + ax2.set_xlabel("Training Data Size (%)") + ax2.set_ylabel("Score") ax2.set_title("(B) Learning Curve — AUPRC & AUROC", fontweight="bold", loc="left") - ax2.legend(); ax2.grid(False) - ax2.set_xlim([0, 100]); ax2.set_ylim([0, 1.05]) + ax2.legend() + ax2.grid(False) + ax2.set_xlim([0, 100]) + ax2.set_ylim([0, 1.05]) plt.tight_layout() - plt.savefig(output_dir / "Fig_Learning_Curve.png", dpi=FIG_DPI, bbox_inches="tight") + plt.savefig(output_dir / "Fig_Learning_Curve.png", dpi=FIG_DPI, bbox_inches="tight") plt.savefig(output_dir / "Fig_Learning_Curve.tiff", dpi=FIG_DPI, bbox_inches="tight") plt.close() results_df.to_csv(output_dir / "learning_curve_data.csv", index=False) print("\n Learning Curve Results:") for _, row in results_df.iterrows(): - print(f" {row['train_size_pct']:5.0f}% ({row['n_train_samples']:,} samples): " - f"MCC={row['mcc']:.4f}, AUPRC={row['auprc']:.4f}") + print( + f" {row['train_size_pct']:5.0f}% ({row['n_train_samples']:,} samples): " + f"MCC={row['mcc']:.4f}, AUPRC={row['auprc']:.4f}" + ) return results_df -# --------------------------------------------------------------------------- -# BaselineEnsemble (unchanged) -# --------------------------------------------------------------------------- +# BaselineEnsemble (unchanged) + class BaselineEnsemble: def __init__(self, n_biological_features=18, random_state=42): @@ -2088,19 +2343,26 @@ def __init__(self, n_biological_features=18, random_state=42): def _create_base_models(self): rs = self.random_state return { - "RF_1": RandomForestClassifier(n_estimators=100, max_depth=10, random_state=rs, n_jobs=-1), - "RF_2": RandomForestClassifier(n_estimators=100, max_depth=15, random_state=rs + 1, n_jobs=-1), - "ET_1": ExtraTreesClassifier( n_estimators=100, max_depth=10, random_state=rs + 2, n_jobs=-1), - "ET_2": ExtraTreesClassifier( n_estimators=100, max_depth=15, random_state=rs + 3, n_jobs=-1), - "GB_1": GradientBoostingClassifier(n_estimators=100, max_depth=5, random_state=rs + 4), - "GB_2": GradientBoostingClassifier(n_estimators=100, max_depth=7, random_state=rs + 5), + "RF_1": RandomForestClassifier( + n_estimators=100, max_depth=10, random_state=rs, n_jobs=-1 + ), + "RF_2": RandomForestClassifier( + n_estimators=100, max_depth=15, random_state=rs + 1, n_jobs=-1 + ), + "ET_1": ExtraTreesClassifier( + n_estimators=100, max_depth=10, random_state=rs + 2, n_jobs=-1 + ), + "ET_2": ExtraTreesClassifier( + n_estimators=100, max_depth=15, random_state=rs + 3, n_jobs=-1 + ), + "GB_1": GradientBoostingClassifier(n_estimators=100, max_depth=5, random_state=rs + 4), + "GB_2": GradientBoostingClassifier(n_estimators=100, max_depth=7, random_state=rs + 5), "HGB_1": HistGradientBoostingClassifier(max_iter=100, max_depth=6, random_state=rs + 6), "HGB_2": HistGradientBoostingClassifier(max_iter=100, max_depth=8, random_state=rs + 7), } - def fit(self, X, y): - X_bio = X[:, :self.n_biological_features] + X_bio = X[:, : self.n_biological_features] self.models = self._create_base_models() for name, model in tqdm(self.models.items(), desc="Training Baseline", ncols=80): model.fit(X_bio, y) @@ -2108,7 +2370,7 @@ def fit(self, X, y): return self def predict_proba(self, X): - X_bio = X[:, :self.n_biological_features] + X_bio = X[:, : self.n_biological_features] probas = [model.predict_proba(X_bio)[:, 1] for model in self.models.values()] return np.mean(probas, axis=0) @@ -2116,18 +2378,17 @@ def predict(self, X, threshold=0.5): return (self.predict_proba(X) >= threshold).astype(int) -# --------------------------------------------------------------------------- -# Metrics helpers (unchanged) -# --------------------------------------------------------------------------- +# Metrics helpers (unchanged) + def compute_metrics(y_true, y_pred, y_pred_proba=None): metrics = { - "mcc": matthews_corrcoef(y_true, y_pred), - "accuracy": accuracy_score(y_true, y_pred), + "mcc": matthews_corrcoef(y_true, y_pred), + "accuracy": accuracy_score(y_true, y_pred), "precision": precision_score(y_true, y_pred, zero_division=0), - "recall": recall_score(y_true, y_pred, zero_division=0), - "f1": f1_score(y_true, y_pred, zero_division=0), - "kappa": cohen_kappa_score(y_true, y_pred), + "recall": recall_score(y_true, y_pred, zero_division=0), + "f1": f1_score(y_true, y_pred, zero_division=0), + "kappa": cohen_kappa_score(y_true, y_pred), } if y_pred_proba is not None: metrics["auprc"] = average_precision_score(y_true, y_pred_proba) @@ -2151,25 +2412,15 @@ def print_metrics_report(metrics, model_name="Model"): print(f" Cohen's κ: {metrics.get('kappa', 0):.6f}") -# --------------------------------------------------------------------------- -# run_final_model() — v4.0 (all supervisor changes active) -# --------------------------------------------------------------------------- +# run_final_model — + def run_final_model(): - """Run TF-DFE v4.0 on the full dataset — all 7 fixes + Changes 1, 2, 3.""" + """Run the full TF-DFE framework end-to-end on the dataset.""" start_time = time.time() - print("\nTF-DFE v4.0: Final Framework Run") + print("\nTF-DFE: Final Framework Run") print("Best Model: Enhanced KNORA-E with OOF Stacking") - print("\nKey changes vs v3.0:") - print(" FIX 1–7 (retained): all original fixes apply") - print(" CHANGE 1: +MLP (BioTDA), +kNN (SML space), +ElasticNet LR; " - "pairwise Q-statistic diversity matrix") - print(" CHANGE 2: TDA over FCGR/sequence representation; features " - f"expanded to ~20 (was 6)") - print(" CHANGE 3: Consensus-margin sample weighting + selective " - f"prediction (abstain θ={Config.SELECTIVE_ABSTAIN_THRESHOLD})") - print(f" BUG FIX: SHAP fallback now uses model object, not string name") print(f"\nOutput Directory: {Config.OUTPUT_DIR}") print("\nStep 1: Load Data & EDA") @@ -2192,53 +2443,53 @@ def run_final_model(): plot_class_distribution(variant_df, Config.OUTPUT_DIR) print("\nStep 2: Feature Engineering") - preprocessor = TFDFEPreprocessor( - use_fcgr=True, use_tda=True, genome_path=Config.GENOME_PATH - ) + preprocessor = TFDFEPreprocessor(use_fcgr=True, use_tda=True, genome_path=Config.GENOME_PATH) X, y = preprocessor.fit_transform(variant_df) n_standard = len(preprocessor.feature_names) - n_fcgr = len(preprocessor.fcgr_feature_names) - n_tda = len(preprocessor.tda_feature_names) + n_fcgr = len(preprocessor.fcgr_feature_names) + n_tda = len(preprocessor.tda_feature_names) print(f"\nFeature engineering complete — Total features: {X.shape[1]}") - print("\nStep 3: Train Best Model (Enhanced KNORA v4.0)") + print("\nStep 3: Train Best Model") X_train, X_test, y_train, y_test = train_test_split( - X, y, test_size=Config.TEST_SIZE, stratify=y, + X, + y, + test_size=Config.TEST_SIZE, + stratify=y, random_state=Config.RANDOM_STATE, ) print(f"\nData Split:") print(f" Train: {len(X_train):,} samples") print(f" Test: {len(X_test):,} samples") - # ------------------------------------------------------------------ - # CHANGE 3: Extract CONSENSUS_SCORE for sample weighting - # ------------------------------------------------------------------ + # Extract CONSENSUS_SCORE for sample weighting consensus_scores_train = None if "CONSENSUS_SCORE" in variant_df.columns: y_arr_full = np.array(y) if hasattr(y, "values") else y _, _, _, _, train_idx_raw, _ = train_test_split( - X, y_arr_full, np.arange(len(y_arr_full)), - test_size=Config.TEST_SIZE, stratify=y, + X, + y_arr_full, + np.arange(len(y_arr_full)), + test_size=Config.TEST_SIZE, + stratify=y, random_state=Config.RANDOM_STATE, ) # Extract training-set CONSENSUS_SCORE (before inner val split) # The factory will handle mis-alignment if lengths differ. cs_all = variant_df["CONSENSUS_SCORE"].values consensus_scores_train = cs_all[train_idx_raw] - print(f" [CHANGE 3] CONSENSUS_SCORE available for " - f"{len(consensus_scores_train):,} training variants.") + print( + f" CONSENSUS_SCORE available for " + f"{len(consensus_scores_train):,} training variants." + ) else: - print(" [CHANGE 3] CONSENSUS_SCORE not in dataset; " - "sample weighting disabled.") + print(" CONSENSUS_SCORE not in dataset; " "sample weighting disabled.") - model = EnhancedTFDFEEnsemble( - n_standard=n_standard, n_fcgr=n_fcgr, n_tda=n_tda - ) + model = EnhancedTFDFEEnsemble(n_standard=n_standard, n_fcgr=n_fcgr, n_tda=n_tda) train_start = time.time() - model.fit(X_train, y_train, - consensus_scores_train=consensus_scores_train) + model.fit(X_train, y_train, consensus_scores_train=consensus_scores_train) train_time = time.time() - train_start print(f"\nTraining complete ({train_time/60:.2f} minutes)") @@ -2251,12 +2502,10 @@ def run_final_model(): metrics = compute_metrics(y_test, y_pred, y_pred_proba) print_metrics_report(metrics, model_name="Enhanced KNORA v4.0 (TF-DFE)") - # ------------------------------------------------------------------ - # CHANGE 3: Selective prediction metrics - # ------------------------------------------------------------------ - print("\nStep 3A: Selective Prediction Metrics (CHANGE 3)") + # Selective prediction metrics + print("\nStep 3A: Selective Prediction Metrics") y_selective, purity_scores, coverage = model.selective_predict(X_test) - confident_mask = (y_selective >= 0) + confident_mask = y_selective >= 0 if confident_mask.sum() > 0: y_te_arr = np.array(y_test) if hasattr(y_test, "values") else y_test sel_mcc = matthews_corrcoef(y_te_arr[confident_mask], y_selective[confident_mask]) @@ -2264,13 +2513,18 @@ def run_final_model(): print(f" Coverage: {coverage:.3%}") print(f" Selective MCC: {sel_mcc:.6f} (vs full {metrics['mcc']:.6f})") print(f" Selective Accuracy: {sel_acc:.6f} (vs full {metrics['accuracy']:.6f})") - selective_df = pd.DataFrame({ - "coverage": [coverage], "n_predicted": [int(confident_mask.sum())], - "n_abstained": [int((~confident_mask).sum())], - "selective_mcc": [sel_mcc], "selective_accuracy": [sel_acc], - "full_mcc": [metrics["mcc"]], "full_accuracy": [metrics["accuracy"]], - "abstain_threshold": [Config.SELECTIVE_ABSTAIN_THRESHOLD], - }) + selective_df = pd.DataFrame( + { + "coverage": [coverage], + "n_predicted": [int(confident_mask.sum())], + "n_abstained": [int((~confident_mask).sum())], + "selective_mcc": [sel_mcc], + "selective_accuracy": [sel_acc], + "full_mcc": [metrics["mcc"]], + "full_accuracy": [metrics["accuracy"]], + "abstain_threshold": [Config.SELECTIVE_ABSTAIN_THRESHOLD], + } + ) selective_df.to_csv(Config.OUTPUT_DIR / "selective_prediction_metrics.csv", index=False) print(" Saved: selective_prediction_metrics.csv") @@ -2279,25 +2533,20 @@ def run_final_model(): plot_confusion_matrix(y_test, y_pred, "TF-DFE", Config.OUTPUT_DIR) plot_tsne(X_test, y_test, "TF-DFE v4.0 Feature Space Separation", Config.OUTPUT_DIR) - # ------------------------------------------------------------------ - # SHAP analysis — BUG FIX: use model OBJECT not string name - # ------------------------------------------------------------------ + # SHAP analysis —: use model OBJECT not string name shap_model = None shap_indices = None shap_feature_names = None if SHAP_AVAILABLE and hasattr(model, "diverse_factory") and model.diverse_factory is not None: for name, m in model.diverse_factory.models.items(): - # BUG FIX: Unwrap CalibratedClassifierCV to get the underlying - # estimator for SHAP compatibility. Never pass the string name. + # Unwrap CalibratedClassifierCV to get the underlying + # estimator for SHAP compatibility. Never pass the string name. base_m = m.base_estimator if isinstance(m, CalibratedClassifierCV) else m - if hasattr(base_m, "estimators_") or any( - k in name for k in ["XGB", "LGB", "CatBoost"] - ): - shap_model = base_m # ← model OBJECT (BUG FIX) + if hasattr(base_m, "estimators_") or any(k in name for k in ["XGB", "LGB", "CatBoost"]): + shap_model = base_m # ← model OBJECT (BUG FIX) shap_indices = model.diverse_factory.feature_indices[name] - shap_feature_names = [preprocessor.all_feature_names[i] - for i in shap_indices] - X_shap = X_test[:min(500, len(X_test)), shap_indices] + shap_feature_names = [preprocessor.all_feature_names[i] for i in shap_indices] + X_shap = X_test[: min(500, len(X_test)), shap_indices] break if shap_model: perform_shap_analysis( @@ -2311,8 +2560,11 @@ def run_final_model(): print("\nStep 4A: Calibration Analysis") y_test_arr = np.array(y_test) if hasattr(y_test, "values") else y_test ece, brier, calibration_bins = TFDFEvaluator.plot_calibration_curve( - y_test_arr, y_pred_proba, Config.OUTPUT_DIR, - model_name="Enhanced KNORA v4.0 (TF-DFE)", n_bins=10, + y_test_arr, + y_pred_proba, + Config.OUTPUT_DIR, + model_name="Enhanced KNORA v4.0 (TF-DFE)", + n_bins=10, ) metrics["ece"] = ece @@ -2322,7 +2574,7 @@ def run_final_model(): n_biological_features=n_standard, random_state=Config.RANDOM_STATE ) baseline_model.fit(X_train, y_train) - y_pred_baseline = baseline_model.predict(X_test) + y_pred_baseline = baseline_model.predict(X_test) y_pred_proba_baseline = baseline_model.predict_proba(X_test) baseline_metrics = compute_metrics(y_test_arr, y_pred_baseline, y_pred_proba_baseline) print(f"\n Baseline Model Metrics:") @@ -2331,25 +2583,29 @@ def run_final_model(): print(f" F1: {baseline_metrics['f1']:.6f}") mcnemar_result = TFDFEvaluator.perform_mcnemar_test( - y_test_arr, y_pred, y_pred_baseline, - model1_name="Enhanced KNORA v4.0", model2_name="Baseline Ensemble", + y_test_arr, + y_pred, + y_pred_baseline, + model1_name="Enhanced KNORA v4.0", + model2_name="Baseline Ensemble", output_dir=Config.OUTPUT_DIR, ) print("\nStep 4C: Error Analysis") y_array = np.array(y) if hasattr(y, "values") else y _, _, _, _, train_indices, test_indices = train_test_split( - X, y_array, np.arange(len(y_array)), - test_size=Config.TEST_SIZE, stratify=y, + X, + y_array, + np.arange(len(y_array)), + test_size=Config.TEST_SIZE, + stratify=y, random_state=Config.RANDOM_STATE, ) misclass_df = TFDFEvaluator.generate_misclassification_report( variant_df, y_test_arr, y_pred, y_pred_proba, test_indices, Config.OUTPUT_DIR ) - # ------------------------------------------------------------------ - # Step 4D: Aggregated Feature Importance — BUG FIX applied - # ------------------------------------------------------------------ + # Step 4D: Aggregated Feature Importance — applied print("\nStep 4D: Aggregated Feature Importance by Category") if SHAP_AVAILABLE and shap_model is not None: try: @@ -2357,28 +2613,30 @@ def run_final_model(): sample_size = min(1000, len(X_test)) X_sample_full = X_test[:sample_size] - # BUG FIX: prefer models with full feature view for category aggregation - agg_model = None + # prefer models with full feature view for category aggregation + agg_model = None agg_indices = None for cand in ["Full_GradientBoosting", "Full_HistGradient"]: if cand in model.diverse_factory.models: raw = model.diverse_factory.models[cand] - # BUG FIX: always extract the underlying estimator object - agg_model = raw.base_estimator if isinstance(raw, CalibratedClassifierCV) else raw + # always extract the underlying estimator object + agg_model = ( + raw.base_estimator if isinstance(raw, CalibratedClassifierCV) else raw + ) agg_indices = model.diverse_factory.feature_indices[cand] break if agg_model is None: - agg_model = shap_model + agg_model = shap_model agg_indices = shap_indices - X_sample_subset = X_sample_full[:, agg_indices] + X_sample_subset = X_sample_full[:, agg_indices] agg_feature_names = [preprocessor.all_feature_names[i] for i in agg_indices] if hasattr(agg_model, "estimators_"): explainer = shap.TreeExplainer(agg_model) else: background = shap.sample(X_sample_subset, 50) - explainer = shap.KernelExplainer(agg_model.predict_proba, background) + explainer = shap.KernelExplainer(agg_model.predict_proba, background) agg_shap_values = explainer.shap_values(X_sample_subset) if isinstance(agg_shap_values, list) and len(agg_shap_values) == 2: @@ -2387,7 +2645,11 @@ def run_final_model(): agg_shap_values = agg_shap_values[:, :, 1] aggregated_df, feature_df = TFDFEvaluator.aggregate_shap_by_category( - agg_shap_values, agg_feature_names, n_standard, n_fcgr, n_tda, + agg_shap_values, + agg_feature_names, + n_standard, + n_fcgr, + n_tda, Config.OUTPUT_DIR, ) TFDFEvaluator.plot_feature_importance_stacked(aggregated_df, Config.OUTPUT_DIR) @@ -2397,19 +2659,28 @@ def run_final_model(): print(" Skipped (SHAP not available or no compatible model)") print("\nStep 4E: Learning Curve Analysis") + def model_factory_fn(n_std, n_fcg, n_td): return EnhancedTFDFEEnsemble(n_standard=n_std, n_fcgr=n_fcg, n_tda=n_td) if len(X) > 50000: print(f" Dataset size ({len(X):,}) is large. Using subset for learning curve...") lc_sample_size = min(50000, len(X)) - lc_indices = np.random.choice(len(X), lc_sample_size, replace=False) - X_lc = X[lc_indices]; y_lc = y_array[lc_indices] + lc_indices = np.random.choice(len(X), lc_sample_size, replace=False) + X_lc = X[lc_indices] + y_lc = y_array[lc_indices] else: - X_lc = X; y_lc = y_array + X_lc = X + y_lc = y_array learning_curve_df = TFDFEvaluator.generate_learning_curve( - X_lc, y_lc, model_factory_fn, n_standard, n_fcgr, n_tda, Config.OUTPUT_DIR, + X_lc, + y_lc, + model_factory_fn, + n_standard, + n_fcgr, + n_tda, + Config.OUTPUT_DIR, train_sizes=[0.1, 0.3, 0.5, 0.7, 0.9], ) @@ -2417,15 +2688,22 @@ def model_factory_fn(n_std, n_fcg, n_td): if len(X) > 100000: print(f" Dataset size ({len(X):,}) is large. Using subset for CV analysis...") cv_sample_size = min(50000, len(X)) - cv_indices = np.random.choice(len(X), cv_sample_size, replace=False) - X_cv = X[cv_indices]; y_cv = y_array[cv_indices] + cv_indices = np.random.choice(len(X), cv_sample_size, replace=False) + X_cv = X[cv_indices] + y_cv = y_array[cv_indices] else: - X_cv = X; y_cv = y_array + X_cv = X + y_cv = y_array cv_df, cv_summary = TFDFEvaluator.perform_repeated_cv( - X_cv, y_cv, model_factory_fn, - n_splits=5, n_repeats=2, - n_standard=n_standard, n_fcgr=n_fcgr, n_tda=n_tda, + X_cv, + y_cv, + model_factory_fn, + n_splits=5, + n_repeats=2, + n_standard=n_standard, + n_fcgr=n_fcgr, + n_tda=n_tda, output_dir=Config.OUTPUT_DIR, ) TFDFEvaluator.plot_cv_distribution(cv_df, Config.OUTPUT_DIR, metrics=["mcc", "auprc"]) @@ -2434,11 +2712,13 @@ def model_factory_fn(n_std, n_fcg, n_td): metrics_df = pd.DataFrame([metrics]) metrics_df.to_csv(Config.OUTPUT_DIR / "final_metrics.csv", index=False) - predictions_df = pd.DataFrame({ - "y_true": y_test.values if hasattr(y_test, "values") else y_test, - "y_pred": y_pred, - "y_pred_proba": y_pred_proba, - }) + predictions_df = pd.DataFrame( + { + "y_true": y_test.values if hasattr(y_test, "values") else y_test, + "y_pred": y_pred, + "y_pred_proba": y_pred_proba, + } + ) predictions_df.to_csv(Config.OUTPUT_DIR / "final_predictions.csv", index=False) with open(Config.OUTPUT_DIR / "final_report.txt", "w", encoding="utf-8") as rf: @@ -2448,8 +2728,10 @@ def model_factory_fn(n_std, n_fcg, n_td): rf.write(" FIX 1–7 (retained from v3)\n") rf.write(" CHANGE 1: +MLP (BioTDA), +kNN_SML, +ElasticNet LR; Q-statistic diversity\n") rf.write(" CHANGE 2: TDA over FCGR/sequence space; ~20 topological features (was 6)\n") - rf.write(f" CHANGE 3: Consensus-margin weighting + selective prediction " - f"(θ={Config.SELECTIVE_ABSTAIN_THRESHOLD})\n") + rf.write( + f" CHANGE 3: Consensus-margin weighting + selective prediction " + f"(θ={Config.SELECTIVE_ABSTAIN_THRESHOLD})\n" + ) rf.write(" BUG FIX: SHAP fallback uses model object, not string name\n\n") rf.write("Dataset Information\n") rf.write(f" Total variants: {len(variant_df):,}\n") @@ -2473,8 +2755,9 @@ def model_factory_fn(n_std, n_fcg, n_td): rf.write(f" Significance: {mcnemar_result['significance']}\n") rf.write(f" {mcnemar_result['interpretation']}\n\n") rf.write("Classification Report\n") - rf.write(classification_report( - y_test, y_pred, target_names=["Benign", "Pathogenic"], digits=4)) + rf.write( + classification_report(y_test, y_pred, target_names=["Benign", "Pathogenic"], digits=4) + ) total_time = time.time() - start_time print(f"\nFINAL ANALYSIS RUN COMPLETE") @@ -2488,8 +2771,9 @@ def model_factory_fn(n_std, n_fcg, n_td): print(f" MCC: {cv_summary['mcc']['formatted']}") print(f" AUPRC: {cv_summary['auprc']['formatted']}") print(f"\nStatistical Significance:") - print(f" McNemar p-value: {mcnemar_result['p_value']:.6f} " - f"({mcnemar_result['significance']})") + print( + f" McNemar p-value: {mcnemar_result['p_value']:.6f} " f"({mcnemar_result['significance']})" + ) print(f" {mcnemar_result['interpretation']}") print(f"\nOutput Directory: {Config.OUTPUT_DIR}") print(f"\nTotal Runtime: {total_time/60:.2f} minutes") diff --git a/src/09_reviewer_experiments.py b/src/09_reviewer_experiments.py index a8e10c3..7503b53 100644 --- a/src/09_reviewer_experiments.py +++ b/src/09_reviewer_experiments.py @@ -1,32 +1,19 @@ -""" -9. Reviewer Revision Experiments — v3.0 (v4 TF-DFE Compatible) -================================================================ -Changes from Code 9 v2.0 (which was already compatible with v3 / 7 fixes): - - CHANGE 1 compatibility: - - audit_feature_counts() now also reports the three new diversity-pool - members (MLP_BioTDA, kNN_SML, ElasticNet_Bio) added in Code 8 v4. - - A new step run_diversity_analysis() reads pairwise_q_statistic.csv and - pairwise_disagreement.csv produced by Code 8's compute_pairwise_diversity() - and renders a summary table + heatmap for the paper. This step is - non-blocking: it skips gracefully if the CSVs are not yet present. - - CHANGE 2 compatibility: - - audit_feature_counts() expected TDA feature count updated from 6 - (old: 3 stats × 2 homology dims) to the new ~20-feature count produced - by the expanded TopologicalFeatureExtractor in Code 8 v4. - The expected value is derived from 10 features × 2 dims = 20. - - run_feature_ablation() TDA slice size comment updated accordingly. - - CHANGE 3 compatibility: - - save_best_params_table() now records selective_abstain_threshold and - ambiguous_margin / ambiguous_sample_weight from RevisionConfig. - - STEP8_SCRIPT_PATH updated to point to the v4 fixed script name so that - build_full_feature_matrix() imports the correct preprocessor automatically. - Update it back if you renamed the file differently. - - All other functions are unchanged from v2.0. +"""Supplementary experiments and audits for the pathogenicity framework. + +Reuses the fitted preprocessor and ensemble from the stage-08 model to produce +the supporting analyses reported alongside the main results: + +- ``audit_feature_counts`` : verify the size of each feature block (FCGR, + TDA, biological) against the expected configuration. +- ``run_diversity_analysis`` : summarise the pairwise Q-statistic and + disagreement matrices written by the ensemble, rendering a table and heatmap. + Skips gracefully when those CSVs are not present. +- ``run_feature_ablation`` : measure the contribution of each feature block. +- ``save_best_params_table`` : export the final model hyper-parameters. + +``STEP8_SCRIPT_PATH`` points at the stage-08 script so that +``build_full_feature_matrix`` imports the correct preprocessor automatically; +update it if that file is renamed. """ import importlib.util @@ -79,18 +66,21 @@ try: from xgboost import XGBClassifier + XGBOOST_AVAILABLE = True except ImportError: XGBOOST_AVAILABLE = False try: from lightgbm import LGBMClassifier + LGBM_AVAILABLE = True except ImportError: LGBM_AVAILABLE = False try: from catboost import CatBoostClassifier + CATBOOST_AVAILABLE = True except ImportError: CATBOOST_AVAILABLE = False @@ -100,22 +90,20 @@ class RevisionConfig: RANDOM_STATE = 42 - TEST_SIZE = 0.2 + TEST_SIZE = 0.2 - BALANCED_DATA_PATH = STAGE07_OUT / "Final_Dataset_Balanced.csv" + BALANCED_DATA_PATH = STAGE07_OUT / "Final_Dataset_Balanced.csv" IMBALANCED_DATA_PATH = STAGE06_OUT / "somatic_variant_Cleaned.csv" - # ----------------------------------------------------------------------- - # CHANGE 1/2/3 compatibility: point at the v4 fixed script so that - # build_full_feature_matrix() picks up the v4 preprocessor. + # point at the fixed script so that + # build_full_feature_matrix picks up the preprocessor. # Resolved relative to this file, not cwd, so it works regardless of # where the script is launched from. If you rename 08_tda_fuzzy_ensemble.py, # update the filename below. - # ----------------------------------------------------------------------- STEP8_SCRIPT_PATH = Path(__file__).resolve().parent / "08_tda_fuzzy_ensemble.py" - GENOME_PATH = DATA_DIR / "hg19.fa" - TARGET_COL = "LABEL_PATHOGENIC" + GENOME_PATH = DATA_DIR / "hg19.fa" + TARGET_COL = "LABEL_PATHOGENIC" LEAKAGE_COLS = ["chr", "pos", "ref", "alt", "CONSENSUS_SCORE", "TIER"] USE_FULL_FEATURES = True @@ -124,32 +112,27 @@ class RevisionConfig: N_JOBS = 1 - KNORA_K = 11 + KNORA_K = 11 MIN_COMPETENCE_THRESHOLD = 0.7 - TDA_N_NEIGHBORS = 75 - SEQUENCE_WINDOW = 75 - FCGR_K_VALUES = [3, 4] - TDA_HOMOLOGY_DIMS = (0, 1) + TDA_N_NEIGHBORS = 75 + SEQUENCE_WINDOW = 75 + FCGR_K_VALUES = [3, 4] + TDA_HOMOLOGY_DIMS = (0, 1) - # FIX 3 compatibility: expected FCGR size after PCA compression. + # expected FCGR size after PCA compression. FCGR_PCA_COMPONENTS = 20 - # ----------------------------------------------------------------------- - # CHANGE 2 compatibility: expected TDA feature count after expansion. - # Code 8 v4 produces 10 features × 2 homology dims = 20 TDA features. - # (was 3 × 2 = 6 in v3) - # ----------------------------------------------------------------------- - TDA_N_FEATURES_EXPECTED = 20 # CHANGE 2 + # expected TDA feature count after expansion. + # script 08 produces 10 features × 2 homology dims = 20 TDA features. + TDA_N_FEATURES_EXPECTED = 20 - # ----------------------------------------------------------------------- - # CHANGE 3 compatibility: selective prediction parameters (for table) - # ----------------------------------------------------------------------- + # selective prediction parameters (for table) SELECTIVE_ABSTAIN_THRESHOLD = 0.60 - AMBIGUOUS_MARGIN = 0.10 - AMBIGUOUS_SAMPLE_WEIGHT = 0.30 + AMBIGUOUS_MARGIN = 0.10 + AMBIGUOUS_SAMPLE_WEIGHT = 0.30 - DPI = 300 - SINGLE_COL = 85 / 25.4 + DPI = 300 + SINGLE_COL = 85 / 25.4 DOUBLE_COL = 170 / 25.4 @@ -158,49 +141,44 @@ class RevisionConfig: plt.rcParams.update( { - "font.family": "sans-serif", + "font.family": "sans-serif", "font.sans-serif": ["Arial", "Helvetica", "DejaVu Sans"], - "font.size": 10, - "axes.titlesize": 12, - "axes.labelsize": 10, - "figure.dpi": RevisionConfig.DPI, - "savefig.dpi": RevisionConfig.DPI, - "savefig.bbox": "tight", - "axes.spines.top": False, + "font.size": 10, + "axes.titlesize": 12, + "axes.labelsize": 10, + "figure.dpi": RevisionConfig.DPI, + "savefig.dpi": RevisionConfig.DPI, + "savefig.bbox": "tight", + "axes.spines.top": False, "axes.spines.right": False, } ) -# --------------------------------------------------------------------------- # Module loader -# --------------------------------------------------------------------------- + def load_step8_module(script_path): path = Path(script_path) if not path.exists(): return None - spec = importlib.util.spec_from_file_location("step8_tfdfe", str(path)) + spec = importlib.util.spec_from_file_location("step8_tfdfe", str(path)) module = importlib.util.module_from_spec(spec) sys.modules["step8_tfdfe"] = module spec.loader.exec_module(module) return module -# --------------------------------------------------------------------------- -# CHANGE 2 compatibility: build_full_feature_matrix -# --------------------------------------------------------------------------- +# build_full_feature_matrix -def build_full_feature_matrix(variant_df): - """Load the v4 preprocessor from Code 8 and build the full feature matrix. - The v4 preprocessor (Code 8 v4) applies: - - PCA(n_components=20) to the FCGR block (FIX 3, retained) - - TDA over FCGR/sequence space with ~20 features per sample (CHANGE 2) +def build_full_feature_matrix(variant_df): + """Load the fitted stage-08 preprocessor and build the full feature matrix. - So n_fcgr = 20, n_tda ≈ 20 here (vs n_tda = 6 in v3). - All downstream slicing in run_feature_ablation and audit_feature_counts - reflects these updated counts. + The preprocessor applies PCA(n_components=20) to the FCGR block and + computes ~20 TDA features per sample over the FCGR/sequence space, so + n_fcgr = 20 and n_tda ≈ 20. All downstream slicing in run_feature_ablation + and audit_feature_counts reflects these counts. """ module = load_step8_module(RevisionConfig.STEP8_SCRIPT_PATH) if module is None or not hasattr(module, "TFDFEPreprocessor"): @@ -218,12 +196,14 @@ def _patched_extract_global_tda(X, n_jobs=-1): from sklearn.neighbors import NearestNeighbors import multiprocessing as mp - extractor = preprocessor.tda_extractor - n_samples = X.shape[0] + extractor = preprocessor.tda_extractor + n_samples = X.shape[0] n_neighbors = min(extractor.n_neighbors, n_samples) - print(f"Finding k-NN neighbors (k={n_neighbors}) for {n_samples} " - "samples in a single optimized pass...") + print( + f"Finding k-NN neighbors (k={n_neighbors}) for {n_samples} " + "samples in a single optimized pass..." + ) nbrs = NearestNeighbors(n_neighbors=n_neighbors, n_jobs=-1) nbrs.fit(X) all_indices = nbrs.kneighbors(X, return_distance=False) @@ -244,98 +224,103 @@ def _patched_extract_global_tda(X, n_jobs=-1): X, y = preprocessor.fit_transform(variant_df) n_standard = len(preprocessor.feature_names) - n_fcgr = len(preprocessor.fcgr_feature_names) # FIX 3: PCA-compressed (20) - n_tda = len(preprocessor.tda_feature_names) # CHANGE 2: expanded (~20) + n_fcgr = len(preprocessor.fcgr_feature_names) # PCA-compressed + n_tda = len(preprocessor.tda_feature_names) # expanded topological block - print(f"\nbuild_full_feature_matrix: n_standard={n_standard}, " - f"n_fcgr={n_fcgr} (PCA-compressed), n_tda={n_tda} (CHANGE 2: over FCGR space)") + print( + f"\nbuild_full_feature_matrix: n_standard={n_standard}, " + f"n_fcgr={n_fcgr} (PCA-compressed), n_tda={n_tda}" + ) return { - "X": np.asarray(X), - "y": np.asarray(y), - "n_standard": n_standard, - "n_fcgr": n_fcgr, - "n_tda": n_tda, + "X": np.asarray(X), + "y": np.asarray(y), + "n_standard": n_standard, + "n_fcgr": n_fcgr, + "n_tda": n_tda, "feature_names": preprocessor.all_feature_names, - "preprocessor": preprocessor, + "preprocessor": preprocessor, } -# --------------------------------------------------------------------------- # Tabular fallback builder (unchanged) -# --------------------------------------------------------------------------- + def build_tabular_features(variant_df, target_col, leakage_cols, fitted=None): df_clean = variant_df.drop(columns=leakage_cols, errors="ignore") - y = df_clean[target_col].values - X_tab = df_clean.drop(columns=[target_col], errors="ignore") + y = df_clean[target_col].values + X_tab = df_clean.drop(columns=[target_col], errors="ignore") num_cols = X_tab.select_dtypes(include=["number"]).columns.tolist() cat_cols = X_tab.select_dtypes(include=["object"]).columns.tolist() if fitted is None: imputer = SimpleImputer(strategy="median") - X_num = pd.DataFrame( + X_num = pd.DataFrame( imputer.fit_transform(X_tab[num_cols]), - columns=num_cols, index=X_tab.index, + columns=num_cols, + index=X_tab.index, ) encoders = {} - X_cat = X_tab[cat_cols].copy() + X_cat = X_tab[cat_cols].copy() for col in cat_cols: encoders[col] = LabelEncoder() - X_cat[col] = encoders[col].fit_transform(X_cat[col].astype(str)) - X_combined = pd.concat([X_num, X_cat], axis=1) + X_cat[col] = encoders[col].fit_transform(X_cat[col].astype(str)) + X_combined = pd.concat([X_num, X_cat], axis=1) feature_names = X_combined.columns.tolist() - scaler = StandardScaler() - X_scaled = scaler.fit_transform(X_combined) + scaler = StandardScaler() + X_scaled = scaler.fit_transform(X_combined) fitted = { - "imputer": imputer, - "encoders": encoders, - "scaler": scaler, - "num_cols": num_cols, - "cat_cols": cat_cols, + "imputer": imputer, + "encoders": encoders, + "scaler": scaler, + "num_cols": num_cols, + "cat_cols": cat_cols, "feature_names": feature_names, } return X_scaled, y, fitted X_num = pd.DataFrame( fitted["imputer"].transform(X_tab[fitted["num_cols"]]), - columns=fitted["num_cols"], index=X_tab.index, + columns=fitted["num_cols"], + index=X_tab.index, ) X_cat = X_tab[fitted["cat_cols"]].copy() for col in fitted["cat_cols"]: - enc = fitted["encoders"][col] + enc = fitted["encoders"][col] known = set(enc.classes_) - X_cat[col] = ( - X_cat[col].astype(str).apply(lambda v: v if v in known else enc.classes_[0]) - ) + X_cat[col] = X_cat[col].astype(str).apply(lambda v: v if v in known else enc.classes_[0]) X_cat[col] = enc.transform(X_cat[col]) X_combined = pd.concat([X_num, X_cat], axis=1)[fitted["feature_names"]] - X_scaled = fitted["scaler"].transform(X_combined) + X_scaled = fitted["scaler"].transform(X_combined) return X_scaled, y, fitted -# --------------------------------------------------------------------------- # Metrics helper (unchanged) -# --------------------------------------------------------------------------- + def compute_extended_metrics(y_true, y_pred, y_pred_proba=None): tn, fp, fn, tp = confusion_matrix(y_true, y_pred, labels=[0, 1]).ravel() specificity = tn / (tn + fp) if (tn + fp) > 0 else 0.0 - fpr = fp / (fp + tn) if (fp + tn) > 0 else 0.0 - fnr = fn / (fn + tp) if (fn + tp) > 0 else 0.0 - npv = tn / (tn + fn) if (tn + fn) > 0 else 0.0 + fpr = fp / (fp + tn) if (fp + tn) > 0 else 0.0 + fnr = fn / (fn + tp) if (fn + tp) > 0 else 0.0 + npv = tn / (tn + fn) if (tn + fn) > 0 else 0.0 metrics = { - "accuracy": accuracy_score(y_true, y_pred), - "mcc": matthews_corrcoef(y_true, y_pred), - "precision": precision_score(y_true, y_pred, zero_division=0), - "recall": recall_score(y_true, y_pred, zero_division=0), + "accuracy": accuracy_score(y_true, y_pred), + "mcc": matthews_corrcoef(y_true, y_pred), + "precision": precision_score(y_true, y_pred, zero_division=0), + "recall": recall_score(y_true, y_pred, zero_division=0), "specificity": specificity, - "f1": f1_score(y_true, y_pred, zero_division=0), - "kappa": cohen_kappa_score(y_true, y_pred), - "fpr": fpr, "fnr": fnr, "npv": npv, - "tn": int(tn), "fp": int(fp), "fn": int(fn), "tp": int(tp), + "f1": f1_score(y_true, y_pred, zero_division=0), + "kappa": cohen_kappa_score(y_true, y_pred), + "fpr": fpr, + "fnr": fnr, + "npv": npv, + "tn": int(tn), + "fp": int(fp), + "fn": int(fn), + "tp": int(tp), } if y_pred_proba is not None: metrics["auroc"] = roc_auc_score(y_true, y_pred_proba) @@ -344,9 +329,8 @@ def compute_extended_metrics(y_true, y_pred, y_pred_proba=None): return metrics -# --------------------------------------------------------------------------- # Baseline comparison models (unchanged) -# --------------------------------------------------------------------------- + def build_comparison_models(random_state): models = {} @@ -354,16 +338,25 @@ def build_comparison_models(random_state): max_iter=2000, C=1.0, n_jobs=RevisionConfig.N_JOBS, random_state=random_state ) models["ExtraTrees"] = ExtraTreesClassifier( - n_estimators=300, max_depth=15, max_features="sqrt", - n_jobs=RevisionConfig.N_JOBS, random_state=random_state, + n_estimators=300, + max_depth=15, + max_features="sqrt", + n_jobs=RevisionConfig.N_JOBS, + random_state=random_state, ) models["RandomForest"] = RandomForestClassifier( - n_estimators=300, max_depth=15, max_features="sqrt", - n_jobs=RevisionConfig.N_JOBS, random_state=random_state, + n_estimators=300, + max_depth=15, + max_features="sqrt", + n_jobs=RevisionConfig.N_JOBS, + random_state=random_state, ) models["GBDT"] = GradientBoostingClassifier( - n_estimators=200, max_depth=4, learning_rate=0.05, - subsample=0.8, random_state=random_state, + n_estimators=200, + max_depth=4, + learning_rate=0.05, + subsample=0.8, + random_state=random_state, ) models["HistGradientBoosting"] = HistGradientBoostingClassifier( max_iter=200, max_depth=6, learning_rate=0.05, random_state=random_state @@ -373,21 +366,33 @@ def build_comparison_models(random_state): ) if XGBOOST_AVAILABLE: models["XGBoost"] = XGBClassifier( - n_estimators=300, max_depth=5, learning_rate=0.05, - subsample=0.8, colsample_bytree=0.8, - n_jobs=RevisionConfig.N_JOBS, verbosity=0, - random_state=random_state, eval_metric="logloss", + n_estimators=300, + max_depth=5, + learning_rate=0.05, + subsample=0.8, + colsample_bytree=0.8, + n_jobs=RevisionConfig.N_JOBS, + verbosity=0, + random_state=random_state, + eval_metric="logloss", ) if LGBM_AVAILABLE: models["LightGBM"] = LGBMClassifier( - n_estimators=300, num_leaves=63, learning_rate=0.05, - feature_fraction=0.8, n_jobs=RevisionConfig.N_JOBS, - verbose=-1, random_state=random_state, + n_estimators=300, + num_leaves=63, + learning_rate=0.05, + feature_fraction=0.8, + n_jobs=RevisionConfig.N_JOBS, + verbose=-1, + random_state=random_state, ) if CATBOOST_AVAILABLE: models["CatBoost"] = CatBoostClassifier( - iterations=300, depth=5, learning_rate=0.05, - verbose=False, random_state=random_state, + iterations=300, + depth=5, + learning_rate=0.05, + verbose=False, + random_state=random_state, ) return models @@ -402,22 +407,38 @@ def predict_proba_safe(model, X): def run_baseline_comparison(X_train, y_train, X_test, y_test, output_dir): models = build_comparison_models(RevisionConfig.RANDOM_STATE) import gc + rows, curve_data, fitted_models = [], {}, {} for name, model in models.items(): model.fit(X_train, y_train) proba = predict_proba_safe(model, X_test) - pred = (proba >= 0.5).astype(int) - m = compute_extended_metrics(y_test, pred, proba) + pred = (proba >= 0.5).astype(int) + m = compute_extended_metrics(y_test, pred, proba) m["model"] = name rows.append(m) - curve_data[name] = (np.array(y_test), proba) + curve_data[name] = (np.array(y_test), proba) fitted_models[name] = model gc.collect() df = pd.DataFrame(rows) ordered = [ - "model", "accuracy", "mcc", "precision", "recall", "specificity", - "f1", "kappa", "auroc", "auprc", "brier", "fpr", "fnr", "npv", - "tn", "fp", "fn", "tp", + "model", + "accuracy", + "mcc", + "precision", + "recall", + "specificity", + "f1", + "kappa", + "auroc", + "auprc", + "brier", + "fpr", + "fnr", + "npv", + "tn", + "fp", + "fn", + "tp", ] df = df[[c for c in ordered if c in df.columns]] df = df.sort_values("mcc", ascending=False).reset_index(drop=True) @@ -425,14 +446,12 @@ def run_baseline_comparison(X_train, y_train, X_test, y_test, output_dir): return df, curve_data, fitted_models -# --------------------------------------------------------------------------- # Plotting helpers (unchanged) -# --------------------------------------------------------------------------- + def plot_combined_roc_pr(curve_data, output_dir, filename="Fig_All_Models_ROC_PR"): fig, (ax1, ax2) = plt.subplots( - 1, 2, - figsize=(RevisionConfig.DOUBLE_COL, RevisionConfig.DOUBLE_COL / 2.4) + 1, 2, figsize=(RevisionConfig.DOUBLE_COL, RevisionConfig.DOUBLE_COL / 2.4) ) palette = plt.cm.tab10(np.linspace(0, 1, max(10, len(curve_data)))) for i, (name, (y_true, proba)) in enumerate(curve_data.items()): @@ -440,35 +459,35 @@ def plot_combined_roc_pr(curve_data, output_dir, filename="Fig_All_Models_ROC_PR auroc = roc_auc_score(y_true, proba) ax1.plot(fpr, tpr, linewidth=1.3, color=palette[i], label=f"{name} ({auroc:.3f})") ax1.plot([0, 1], [0, 1], "k--", linewidth=0.8, alpha=0.6) - ax1.set_xlabel("False Positive Rate"); ax1.set_ylabel("True Positive Rate") + ax1.set_xlabel("False Positive Rate") + ax1.set_ylabel("True Positive Rate") ax1.set_title("(A) ROC Curves", fontweight="bold", loc="left") ax1.legend(loc="lower right", fontsize=6, frameon=True, edgecolor="black") ax1.grid(False) for i, (name, (y_true, proba)) in enumerate(curve_data.items()): precision, recall, _ = precision_recall_curve(y_true, proba) auprc = average_precision_score(y_true, proba) - ax2.plot(recall, precision, linewidth=1.3, color=palette[i], - label=f"{name} ({auprc:.3f})") - ax2.set_xlabel("Recall"); ax2.set_ylabel("Precision") + ax2.plot(recall, precision, linewidth=1.3, color=palette[i], label=f"{name} ({auprc:.3f})") + ax2.set_xlabel("Recall") + ax2.set_ylabel("Precision") ax2.set_title("(B) Precision-Recall Curves", fontweight="bold", loc="left") ax2.legend(loc="lower left", fontsize=6, frameon=True, edgecolor="black") ax2.grid(False) plt.tight_layout() - plt.savefig(output_dir / f"{filename}.png", dpi=RevisionConfig.DPI, bbox_inches="tight") + plt.savefig(output_dir / f"{filename}.png", dpi=RevisionConfig.DPI, bbox_inches="tight") plt.savefig(output_dir / f"{filename}.tiff", dpi=RevisionConfig.DPI, bbox_inches="tight") plt.close() -# --------------------------------------------------------------------------- # McNemar (unchanged) -# --------------------------------------------------------------------------- + def mcnemar_test(y_true, ref_pred, other_pred): - y_true = np.asarray(y_true) - ref_correct = ref_pred == y_true + y_true = np.asarray(y_true) + ref_correct = ref_pred == y_true other_correct = other_pred == y_true - b = int(np.sum(ref_correct & ~other_correct)) - c = int(np.sum(~ref_correct & other_correct)) + b = int(np.sum(ref_correct & ~other_correct)) + c = int(np.sum(~ref_correct & other_correct)) if b + c == 0: return np.nan, np.nan statistic = (abs(b - c) - 1) ** 2 / (b + c) @@ -497,166 +516,193 @@ def run_mcnemar_matrix(fitted_models, X_test, y_test, reference_name, output_dir stat, pval = mcnemar_test(y_test, ref_pred, pred) sig = "ns" if pval is not None and not np.isnan(pval): - if pval < 0.001: sig = "***" - elif pval < 0.01: sig = "**" - elif pval < 0.05: sig = "*" - rows.append({ - "reference": reference_name, "compared_model": name, - "mcnemar_statistic": stat, "p_value": pval, "significance": sig, - }) + if pval < 0.001: + sig = "***" + elif pval < 0.01: + sig = "**" + elif pval < 0.05: + sig = "*" + rows.append( + { + "reference": reference_name, + "compared_model": name, + "mcnemar_statistic": stat, + "p_value": pval, + "significance": sig, + } + ) df = pd.DataFrame(rows) df.to_csv(output_dir / "mcnemar_vs_reference.csv", index=False) return df -# --------------------------------------------------------------------------- -# CHANGE 1 compatibility: run_diversity_analysis -# --------------------------------------------------------------------------- +# run_diversity_analysis -def run_diversity_analysis(output_dir): - """Read Code 8 v4's pairwise diversity CSVs and produce a summary report. - CHANGE 1 (Supervisor Request 1): Code 8 v4 now writes - - pairwise_q_statistic.csv - - pairwise_disagreement.csv - produced by compute_pairwise_diversity() after training. +def run_diversity_analysis(output_dir): + """Read the pairwise diversity CSVs and produce a summary report. + The stage-08 model writes pairwise_q_statistic.csv and + pairwise_disagreement.csv (via compute_pairwise_diversity) after training. This function: 1. Reads those CSVs (non-blocking: skips if absent). 2. Computes summary statistics (mean Q-stat, mean disagreement). 3. Saves diversity_summary.csv and a heatmap figure. - 4. Prints a pass/fail verdict against the thresholds recommended by - the supervisor (mean Q < 0.50 and mean disagreement > 0.05). + 4. Prints a pass/fail verdict against the target thresholds + (mean Q < 0.50 and mean disagreement > 0.05). Returns a summary dict (or None if the files are not present). """ - code8_out = Path(STAGE08_OUT) - q_path = code8_out / "pairwise_q_statistic.csv" - dis_path = code8_out / "pairwise_disagreement.csv" + stage08_out = Path(STAGE08_OUT) + q_path = stage08_out / "pairwise_q_statistic.csv" + dis_path = stage08_out / "pairwise_disagreement.csv" if not q_path.exists() or not dis_path.exists(): print( - "\n [CHANGE 1 / Diversity Analysis] " + "\n [Diversity Analysis] " "pairwise_q_statistic.csv / pairwise_disagreement.csv not found " - f"in '{code8_out}/'. Run 08_tda_fuzzy_ensemble.py first to generate them.\n" + f"in '{stage08_out}/'. Run 08_tda_fuzzy_ensemble.py first to generate them.\n" " Skipping diversity analysis." ) - pd.DataFrame([{ - "status": "skipped", - "reason": "08 outputs not found", - }]).to_csv(output_dir / "diversity_summary.csv", index=False) + pd.DataFrame( + [ + { + "status": "skipped", + "reason": "08 outputs not found", + } + ] + ).to_csv(output_dir / "diversity_summary.csv", index=False) return None - q_df = pd.read_csv(q_path, index_col=0) + q_df = pd.read_csv(q_path, index_col=0) dis_df = pd.read_csv(dis_path, index_col=0) - names = q_df.index.tolist() - n = len(names) + names = q_df.index.tolist() + n = len(names) # Upper-triangle values (excluding diagonal) iu = np.triu_indices(n, k=1) - q_vals = q_df.values[iu] + q_vals = q_df.values[iu] dis_vals = dis_df.values[iu] - mean_q = float(np.mean(q_vals)) + mean_q = float(np.mean(q_vals)) mean_dis = float(np.mean(dis_vals)) - max_q = float(np.max(q_vals)) - min_dis = float(np.min(dis_vals)) + max_q = float(np.max(q_vals)) + min_dis = float(np.min(dis_vals)) - q_target_pass = mean_q < 0.50 + q_target_pass = mean_q < 0.50 dis_target_pass = mean_dis > 0.05 - print("\n [CHANGE 1 / Diversity Analysis — Code 9 Summary]") + print("\n [Diversity Analysis Summary]") print(f" Pool members: {names}") - print(f" Mean Q-statistic: {mean_q:.4f} " - f"({'PASS' if q_target_pass else 'FAIL'} — target < 0.50)") + print( + f" Mean Q-statistic: {mean_q:.4f} " + f"({'PASS' if q_target_pass else 'FAIL'} — target < 0.50)" + ) print(f" Max Q-statistic: {max_q:.4f}") - print(f" Mean Disagreement: {mean_dis:.4f} " - f"({'PASS' if dis_target_pass else 'FAIL'} — target > 0.05)") + print( + f" Mean Disagreement: {mean_dis:.4f} " + f"({'PASS' if dis_target_pass else 'FAIL'} — target > 0.05)" + ) print(f" Min Disagreement: {min_dis:.4f}") summary = { - "n_models": n, - "mean_q_statistic": mean_q, - "max_q_statistic": max_q, - "mean_disagreement": mean_dis, - "min_disagreement": min_dis, - "q_target_pass": q_target_pass, - "dis_target_pass": dis_target_pass, + "n_models": n, + "mean_q_statistic": mean_q, + "max_q_statistic": max_q, + "mean_disagreement": mean_dis, + "min_disagreement": min_dis, + "q_target_pass": q_target_pass, + "dis_target_pass": dis_target_pass, "overall_diversity_ok": q_target_pass and dis_target_pass, } pd.DataFrame([summary]).to_csv(output_dir / "diversity_summary.csv", index=False) # ---- Plot heatmaps ---- fig, axes = plt.subplots( - 1, 2, + 1, + 2, figsize=(RevisionConfig.DOUBLE_COL, RevisionConfig.DOUBLE_COL / 2.2), ) sns.heatmap( - q_df, ax=axes[0], cmap="RdYlGn_r", vmin=-1, vmax=1, - annot=True, fmt=".2f", linewidths=0.2, annot_kws={"size": 4}, + q_df, + ax=axes[0], + cmap="RdYlGn_r", + vmin=-1, + vmax=1, + annot=True, + fmt=".2f", + linewidths=0.2, + annot_kws={"size": 4}, ) axes[0].set_title( f"(A) Q-Statistic (mean={mean_q:.3f})\n" f"{'PASS' if q_target_pass else 'FAIL'}: target < 0.50", - fontweight="bold", loc="left", + fontweight="bold", + loc="left", ) axes[0].tick_params(axis="x", rotation=45, labelsize=7) - axes[0].tick_params(axis="y", rotation=0, labelsize=7) + axes[0].tick_params(axis="y", rotation=0, labelsize=7) sns.heatmap( - dis_df, ax=axes[1], cmap="Blues", vmin=0, vmax=0.30, - annot=True, fmt=".2f", linewidths=0.2, annot_kws={"size": 4}, + dis_df, + ax=axes[1], + cmap="Blues", + vmin=0, + vmax=0.30, + annot=True, + fmt=".2f", + linewidths=0.2, + annot_kws={"size": 4}, ) axes[1].set_title( f"(B) Disagreement Rate (mean={mean_dis:.3f})\n" f"{'PASS' if dis_target_pass else 'FAIL'}: target > 0.05", - fontweight="bold", loc="left", + fontweight="bold", + loc="left", ) axes[1].tick_params(axis="x", rotation=45, labelsize=7) - axes[1].tick_params(axis="y", rotation=0, labelsize=7) + axes[1].tick_params(axis="y", rotation=0, labelsize=7) plt.tight_layout() - plt.savefig(output_dir / "Fig_Diversity_Summary.png", - dpi=RevisionConfig.DPI, bbox_inches="tight") - plt.savefig(output_dir / "Fig_Diversity_Summary.tiff", - dpi=RevisionConfig.DPI, bbox_inches="tight") + plt.savefig( + output_dir / "Fig_Diversity_Summary.png", dpi=RevisionConfig.DPI, bbox_inches="tight" + ) + plt.savefig( + output_dir / "Fig_Diversity_Summary.tiff", dpi=RevisionConfig.DPI, bbox_inches="tight" + ) plt.close() print(" Saved: diversity_summary.csv, Fig_Diversity_Summary.png") return summary -# --------------------------------------------------------------------------- -# CHANGE 2 + FIX 3 compatibility: run_feature_ablation -# --------------------------------------------------------------------------- +# run_feature_ablation -def run_feature_ablation( - X_train, y_train, X_test, y_test, n_standard, n_fcgr, n_tda, output_dir -): - """Feature ablation study. - CHANGE 2 compatibility: - n_tda is now ~20 (was 6). The TDA slice is updated accordingly. - Comment text and the FCGR-compressed branch remain the same as v2. +def run_feature_ablation(X_train, y_train, X_test, y_test, n_standard, n_fcgr, n_tda, output_dir): + """Feature ablation study. - FIX 3 compatibility (retained from v2): - When n_fcgr <= FCGR_PCA_COMPONENTS (default 20), the FCGR block is already - PCA-compressed; skip the redundant inner-PCA loop. + n_tda is ~20; the TDA slice is sized accordingly. When + n_fcgr <= FCGR_PCA_COMPONENTS (default 20) the FCGR block is already + PCA-compressed, so the redundant inner-PCA loop is skipped. """ - std_slice = slice(0, n_standard) + std_slice = slice(0, n_standard) fcgr_slice = slice(n_standard, n_standard + n_fcgr) - tda_slice = slice(n_standard + n_fcgr, n_standard + n_fcgr + n_tda) + tda_slice = slice(n_standard + n_fcgr, n_standard + n_fcgr + n_tda) def fit_eval(Xtr, Xte, tag): model = HistGradientBoostingClassifier( - max_iter=200, max_depth=6, learning_rate=0.05, - random_state=RevisionConfig.RANDOM_STATE, early_stopping=False, + max_iter=200, + max_depth=6, + learning_rate=0.05, + random_state=RevisionConfig.RANDOM_STATE, + early_stopping=False, ) model.fit(Xtr, y_train) proba = model.predict_proba(Xte)[:, 1] - pred = (proba >= 0.5).astype(int) - m = compute_extended_metrics(y_test, pred, proba) + pred = (proba >= 0.5).astype(int) + m = compute_extended_metrics(y_test, pred, proba) m["configuration"] = tag - m["n_features"] = Xtr.shape[1] + m["n_features"] = Xtr.shape[1] return m, (np.array(y_test), proba) configs, curves = {}, {} @@ -667,120 +713,107 @@ def fit_eval(Xtr, Xte, tag): if n_fcgr > 0: Xtr = np.hstack([X_train[:, std_slice], X_train[:, fcgr_slice]]) - Xte = np.hstack([X_test[:, std_slice], X_test[:, fcgr_slice]]) + Xte = np.hstack([X_test[:, std_slice], X_test[:, fcgr_slice]]) configs["Bio_plus_FCGR_compressed"], curves["Bio_plus_FCGR_compressed"] = fit_eval( Xtr, Xte, "Bio_plus_FCGR_compressed" ) if n_tda > 0: - # CHANGE 2: TDA slice is now ~20 features (over FCGR space) + # TDA slice is now ~20 features (over FCGR space) Xtr = np.hstack([X_train[:, std_slice], X_train[:, tda_slice]]) - Xte = np.hstack([X_test[:, std_slice], X_test[:, tda_slice]]) + Xte = np.hstack([X_test[:, std_slice], X_test[:, tda_slice]]) configs["Bio_plus_TDA_FCGR"], curves["Bio_plus_TDA_FCGR"] = fit_eval( Xtr, Xte, "Bio_plus_TDA_FCGR" ) - configs["All_features"], curves["All_features"] = fit_eval( - X_train, X_test, "All_features" - ) + configs["All_features"], curves["All_features"] = fit_eval(X_train, X_test, "All_features") - # ----------------------------------------------------------------------- - # FIX 3 / CHANGE 2: If FCGR is already compressed, test sub-compression - # ----------------------------------------------------------------------- - already_compressed = (n_fcgr <= RevisionConfig.FCGR_PCA_COMPONENTS) + # If FCGR is already compressed, test sub-compression + already_compressed = n_fcgr <= RevisionConfig.FCGR_PCA_COMPONENTS if n_fcgr > 0 and not already_compressed: for n_comp in [10, 20, 50]: n_comp_eff = min(n_comp, n_fcgr) pca = PCA(n_components=n_comp_eff, random_state=RevisionConfig.RANDOM_STATE) fcgr_tr_pca = pca.fit_transform(X_train[:, fcgr_slice]) - fcgr_te_pca = pca.transform(X_test[:, fcgr_slice]) + fcgr_te_pca = pca.transform(X_test[:, fcgr_slice]) evr = float(np.sum(pca.explained_variance_ratio_)) parts_tr = [X_train[:, std_slice], fcgr_tr_pca] - parts_te = [X_test[:, std_slice], fcgr_te_pca] + parts_te = [X_test[:, std_slice], fcgr_te_pca] if n_tda > 0: parts_tr.append(X_train[:, tda_slice]) - parts_te.append(X_test[:, tda_slice]) + parts_te.append(X_test[:, tda_slice]) Xtr = np.hstack(parts_tr) Xte = np.hstack(parts_te) tag = f"Bio_TDA_FCGR-PCA{n_comp_eff}" m, c = fit_eval(Xtr, Xte, tag) m["fcgr_pca_explained_variance"] = evr configs[tag] = m - curves[tag] = c + curves[tag] = c elif n_fcgr > 0 and already_compressed: - print(f" [FIX 3 / CHANGE 2] FCGR already PCA-compressed ({n_fcgr}D). " - f"Testing sub-compression to 5 and 10 components...") + print( + f" FCGR already PCA-compressed ({n_fcgr}D). " + f"Testing sub-compression to 5 and 10 components..." + ) for n_sub in [5, 10]: n_sub_eff = min(n_sub, n_fcgr) pca = PCA(n_components=n_sub_eff, random_state=RevisionConfig.RANDOM_STATE) fcgr_tr_sub = pca.fit_transform(X_train[:, fcgr_slice]) - fcgr_te_sub = pca.transform(X_test[:, fcgr_slice]) + fcgr_te_sub = pca.transform(X_test[:, fcgr_slice]) evr = float(np.sum(pca.explained_variance_ratio_)) parts_tr = [X_train[:, std_slice], fcgr_tr_sub] - parts_te = [X_test[:, std_slice], fcgr_te_sub] + parts_te = [X_test[:, std_slice], fcgr_te_sub] if n_tda > 0: parts_tr.append(X_train[:, tda_slice]) - parts_te.append(X_test[:, tda_slice]) + parts_te.append(X_test[:, tda_slice]) Xtr = np.hstack(parts_tr) Xte = np.hstack(parts_te) tag = f"Bio_TDA_FCGR-sub{n_sub_eff}" m, c = fit_eval(Xtr, Xte, tag) m["fcgr_pca_explained_variance"] = evr configs[tag] = m - curves[tag] = c + curves[tag] = c df = pd.DataFrame(list(configs.values())) - front = ["configuration", "n_features", "mcc", "accuracy", "f1", - "auroc", "auprc", "fpr", "fnr"] - df = df[ - [c for c in front if c in df.columns] - + [c for c in df.columns if c not in front] - ] + front = ["configuration", "n_features", "mcc", "accuracy", "f1", "auroc", "auprc", "fpr", "fnr"] + df = df[[c for c in front if c in df.columns] + [c for c in df.columns if c not in front]] df = df.sort_values("mcc", ascending=False).reset_index(drop=True) df.to_csv(output_dir / "fcgr_ablation_results.csv", index=False) plot_combined_roc_pr(curves, output_dir, filename="Fig_Ablation_ROC_PR") return df -# --------------------------------------------------------------------------- # Threshold sweep (unchanged) -# --------------------------------------------------------------------------- + def run_threshold_sweep(model, X_test, y_test, output_dir): proba = predict_proba_safe(model, X_test) - rows = [] + rows = [] for thr in np.round(np.arange(0.1, 0.95, 0.05), 2): pred = (proba >= thr).astype(int) - m = compute_extended_metrics(y_test, pred, proba) + m = compute_extended_metrics(y_test, pred, proba) m["threshold"] = float(thr) rows.append(m) df = pd.DataFrame(rows) front = ["threshold", "recall", "specificity", "precision", "fpr", "fnr", "f1", "mcc"] - df = df[ - [c for c in front if c in df.columns] - + [c for c in df.columns if c not in front] - ] + df = df[[c for c in front if c in df.columns] + [c for c in df.columns if c not in front]] df.to_csv(output_dir / "threshold_sensitivity_sweep.csv", index=False) return df -# --------------------------------------------------------------------------- # Imbalanced evaluation (unchanged) -# --------------------------------------------------------------------------- + def run_imbalanced_evaluation(bio_model, n_standard, train_fitted, output_dir): path = Path(RevisionConfig.IMBALANCED_DATA_PATH) if not path.exists(): - empty = pd.DataFrame( - [{"status": "imbalanced_source_not_found", "path": str(path)}] - ) + empty = pd.DataFrame([{"status": "imbalanced_source_not_found", "path": str(path)}]) empty.to_csv(output_dir / "imbalanced_fp_evaluation.csv", index=False) return empty try: sample_head = pd.read_csv(path, nrows=5, low_memory=False) - all_cols = sample_head.columns.tolist() + all_cols = sample_head.columns.tolist() except Exception as e: empty = pd.DataFrame([{"status": f"read_error: {e}"}]) empty.to_csv(output_dir / "imbalanced_fp_evaluation.csv", index=False) @@ -788,9 +821,9 @@ def run_imbalanced_evaluation(bio_model, n_standard, train_fitted, output_dir): print(" Reading imbalanced dataset in chunks...") pathogenic_chunks, benign_chunks = [], [] - chunk_size = 20000 + chunk_size = 20000 max_benign_retained = 500000 - n_benign_retained = 0 + n_benign_retained = 0 def _downcast(frame): for c in frame.select_dtypes(include=["float64"]).columns: @@ -803,9 +836,9 @@ def _downcast(frame): for chunk in reader: if RevisionConfig.TARGET_COL not in chunk.columns: continue - chunk = _downcast(chunk) - path_rows = chunk[chunk[RevisionConfig.TARGET_COL] == 1] - ben_rows = chunk[chunk[RevisionConfig.TARGET_COL] == 0] + chunk = _downcast(chunk) + path_rows = chunk[chunk[RevisionConfig.TARGET_COL] == 1] + ben_rows = chunk[chunk[RevisionConfig.TARGET_COL] == 0] if len(path_rows) > 0: pathogenic_chunks.append(path_rows) if len(ben_rows) > 0 and n_benign_retained < max_benign_retained: @@ -818,11 +851,11 @@ def _downcast(frame): return empty pathogenic = pd.concat(pathogenic_chunks, ignore_index=True) - benign = pd.concat(benign_chunks, ignore_index=True) + benign = pd.concat(benign_chunks, ignore_index=True) del pathogenic_chunks, benign_chunks n_path_test = max(1, int(len(pathogenic) * RevisionConfig.TEST_SIZE)) - path_test = pathogenic.sample( + path_test = pathogenic.sample( n=min(n_path_test, len(pathogenic)), random_state=RevisionConfig.RANDOM_STATE ) del pathogenic @@ -840,51 +873,54 @@ def _downcast(frame): del path_test, benign_test X_imb, y_imb, _ = build_tabular_features( - imb_test, RevisionConfig.TARGET_COL, RevisionConfig.LEAKAGE_COLS, + imb_test, + RevisionConfig.TARGET_COL, + RevisionConfig.LEAKAGE_COLS, fitted=train_fitted, ) del imb_test X_imb = X_imb[:, :n_standard] - rows = [] + rows = [] proba = bio_model.predict_proba(X_imb)[:, 1] for thr in [0.3, 0.4, 0.5, 0.6, 0.7]: pred = (proba >= thr).astype(int) - m = compute_extended_metrics(y_imb, pred, proba) - m["threshold"] = thr - m["n_test"] = len(y_imb) - m["benign_pathogenic_ratio"] = ( - f"{int((y_imb == 0).sum())}:{int((y_imb == 1).sum())}" - ) + m = compute_extended_metrics(y_imb, pred, proba) + m["threshold"] = thr + m["n_test"] = len(y_imb) + m["benign_pathogenic_ratio"] = f"{int((y_imb == 0).sum())}:{int((y_imb == 1).sum())}" rows.append(m) df = pd.DataFrame(rows) - front = ["threshold", "benign_pathogenic_ratio", "n_test", "fpr", - "recall", "specificity", "precision", "mcc", "auroc", "auprc"] - df = df[ - [c for c in front if c in df.columns] - + [c for c in df.columns if c not in front] + front = [ + "threshold", + "benign_pathogenic_ratio", + "n_test", + "fpr", + "recall", + "specificity", + "precision", + "mcc", + "auroc", + "auprc", ] + df = df[[c for c in front if c in df.columns] + [c for c in df.columns if c not in front]] df.to_csv(output_dir / "imbalanced_fp_evaluation.csv", index=False) return df -# --------------------------------------------------------------------------- -# CHANGE 2 + FIX 3 + CHANGE 3 compatibility: audit_feature_counts -# --------------------------------------------------------------------------- +# audit_feature_counts + def audit_feature_counts(n_standard, n_fcgr, n_tda, output_dir): """Audit feature counts. - CHANGE 2: expected TDA count is now TDA_N_FEATURES_EXPECTED (20), - reflecting the expanded TopologicalFeatureExtractor in Code 8 v4 that - computes 10 descriptors × 2 homology dimensions over the FCGR space. - (was 3 × 2 = 6 in v3) - - FIX 3 (retained): expected FCGR count is FCGR_PCA_COMPONENTS (20). + The expected TDA count is TDA_N_FEATURES_EXPECTED (20): 10 descriptors x + 2 homology dimensions over the FCGR space. The expected FCGR count is + FCGR_PCA_COMPONENTS (20). """ - expected_fcgr = RevisionConfig.FCGR_PCA_COMPONENTS # FIX 3 - expected_tda = RevisionConfig.TDA_N_FEATURES_EXPECTED # CHANGE 2: 20 (was 6) + expected_fcgr = RevisionConfig.FCGR_PCA_COMPONENTS + expected_tda = RevisionConfig.TDA_N_FEATURES_EXPECTED rows = [ { "category": "Standard/Biological", @@ -906,7 +942,7 @@ def audit_feature_counts(n_standard, n_fcgr, n_tda, output_dir): }, { "category": "Total", - "actual": n_standard + n_fcgr + n_tda, + "actual": n_standard + n_fcgr + n_tda, "expected": n_standard + expected_fcgr + expected_tda, "note": "", }, @@ -919,17 +955,13 @@ def audit_feature_counts(n_standard, n_fcgr, n_tda, output_dir): return df -# --------------------------------------------------------------------------- # Consensus score verification (unchanged) -# --------------------------------------------------------------------------- + def verify_consensus_score(variant_df, output_dir): needed = ["REVEL_score", "SIFT_score", "Polyphen2_HDIV_score"] - rows = [] - if ( - all(c in variant_df.columns for c in needed) - and "CONSENSUS_SCORE" in variant_df.columns - ): + rows = [] + if all(c in variant_df.columns for c in needed) and "CONSENSUS_SCORE" in variant_df.columns: sub = variant_df[needed + ["CONSENSUS_SCORE"]].dropna().head(100000).copy() recomputed = ( sub["REVEL_score"].fillna(0.5) * 0.6 @@ -937,79 +969,92 @@ def verify_consensus_score(variant_df, output_dir): + sub["Polyphen2_HDIV_score"].fillna(0.5) * 0.2 ) diff = np.abs(recomputed - sub["CONSENSUS_SCORE"]) - rows.append({ - "formula": "0.6*REVEL + 0.2*(1-SIFT) + 0.2*PolyPhen2", - "n_checked": int(len(sub)), - "mean_abs_diff": float(diff.mean()), - "max_abs_diff": float(diff.max()), - "sift_inverted": True, - }) + rows.append( + { + "formula": "0.6*REVEL + 0.2*(1-SIFT) + 0.2*PolyPhen2", + "n_checked": int(len(sub)), + "mean_abs_diff": float(diff.mean()), + "max_abs_diff": float(diff.max()), + "sift_inverted": True, + } + ) else: - rows.append({ - "formula": "0.6*REVEL + 0.2*(1-SIFT) + 0.2*PolyPhen2", - "n_checked": 0, - "mean_abs_diff": None, - "max_abs_diff": None, - "sift_inverted": True, - "note": "CONSENSUS_SCORE or component columns not present", - }) + rows.append( + { + "formula": "0.6*REVEL + 0.2*(1-SIFT) + 0.2*PolyPhen2", + "n_checked": 0, + "mean_abs_diff": None, + "max_abs_diff": None, + "sift_inverted": True, + "note": "CONSENSUS_SCORE or component columns not present", + } + ) df = pd.DataFrame(rows) df.to_csv(output_dir / "consensus_score_verification.csv", index=False) return df -# --------------------------------------------------------------------------- # Overfitting analysis (unchanged) -# --------------------------------------------------------------------------- + def run_overfitting_analysis(X, y, n_standard, output_dir): train_sizes = [0.1, 0.3, 0.5, 0.7, 0.9] - rows = [] + rows = [] for ts in train_sizes: X_tr, X_te, y_tr, y_te = train_test_split( - X[:, :n_standard], y, train_size=ts, stratify=y, + X[:, :n_standard], + y, + train_size=ts, + stratify=y, random_state=RevisionConfig.RANDOM_STATE, ) model = HistGradientBoostingClassifier( - max_iter=200, max_depth=6, learning_rate=0.05, - random_state=RevisionConfig.RANDOM_STATE, early_stopping=False, + max_iter=200, + max_depth=6, + learning_rate=0.05, + random_state=RevisionConfig.RANDOM_STATE, + early_stopping=False, ) model.fit(X_tr, y_tr) train_pred = model.predict(X_tr) - test_pred = model.predict(X_te) - rows.append({ - "train_fraction": ts, "n_train": len(y_tr), - "train_mcc": matthews_corrcoef(y_tr, train_pred), - "test_mcc": matthews_corrcoef(y_te, test_pred), - "train_test_gap": matthews_corrcoef(y_tr, train_pred) - matthews_corrcoef(y_te, test_pred), - "train_accuracy": accuracy_score(y_tr, train_pred), - "test_accuracy": accuracy_score(y_te, test_pred), - }) + test_pred = model.predict(X_te) + rows.append( + { + "train_fraction": ts, + "n_train": len(y_tr), + "train_mcc": matthews_corrcoef(y_tr, train_pred), + "test_mcc": matthews_corrcoef(y_te, test_pred), + "train_test_gap": matthews_corrcoef(y_tr, train_pred) + - matthews_corrcoef(y_te, test_pred), + "train_accuracy": accuracy_score(y_tr, train_pred), + "test_accuracy": accuracy_score(y_te, test_pred), + } + ) df = pd.DataFrame(rows) df.to_csv(output_dir / "overfitting_learning_curve.csv", index=False) - fig, ax = plt.subplots( - figsize=(RevisionConfig.SINGLE_COL * 1.4, RevisionConfig.SINGLE_COL)) + fig, ax = plt.subplots(figsize=(RevisionConfig.SINGLE_COL * 1.4, RevisionConfig.SINGLE_COL)) ax.plot(df["n_train"], df["train_mcc"], "o-", color="#0077BB", label="Train MCC") - ax.plot(df["n_train"], df["test_mcc"], "s--", color="#CC3311", label="Test MCC") - ax.fill_between(df["n_train"], df["test_mcc"], df["train_mcc"], - alpha=0.15, color="gray") - ax.set_xlabel("Number of training samples"); ax.set_ylabel("MCC") + ax.plot(df["n_train"], df["test_mcc"], "s--", color="#CC3311", label="Test MCC") + ax.fill_between(df["n_train"], df["test_mcc"], df["train_mcc"], alpha=0.15, color="gray") + ax.set_xlabel("Number of training samples") + ax.set_ylabel("MCC") ax.set_title("Train vs Test MCC (overfitting check)", fontweight="bold", loc="left") ax.legend(frameon=True, edgecolor="black") ax.grid(False) plt.tight_layout() - plt.savefig(output_dir / "Fig_Overfitting_Check.png", - dpi=RevisionConfig.DPI, bbox_inches="tight") - plt.savefig(output_dir / "Fig_Overfitting_Check.tiff", - dpi=RevisionConfig.DPI, bbox_inches="tight") + plt.savefig( + output_dir / "Fig_Overfitting_Check.png", dpi=RevisionConfig.DPI, bbox_inches="tight" + ) + plt.savefig( + output_dir / "Fig_Overfitting_Check.tiff", dpi=RevisionConfig.DPI, bbox_inches="tight" + ) plt.close() return df -# --------------------------------------------------------------------------- # Parameter sensitivity (unchanged) -# --------------------------------------------------------------------------- + def run_parameter_sensitivity( X_train, y_train, X_test, y_test, n_standard, n_fcgr, n_tda, output_dir @@ -1020,8 +1065,11 @@ def run_parameter_sensitivity( return pd.DataFrame() base = HistGradientBoostingClassifier( - max_iter=200, max_depth=6, learning_rate=0.05, - random_state=RevisionConfig.RANDOM_STATE, early_stopping=False, + max_iter=200, + max_depth=6, + learning_rate=0.05, + random_state=RevisionConfig.RANDOM_STATE, + early_stopping=False, ) base.fit(X_train, y_train) proba = base.predict_proba(X_test)[:, 1] @@ -1031,117 +1079,143 @@ def run_parameter_sensitivity( for min_comp in [0.5, 0.6, 0.7, 0.8]: nn = NearestNeighbors(n_neighbors=k, n_jobs=RevisionConfig.N_JOBS) nn.fit(X_train) - _, idx = nn.kneighbors(X_test) - neighbor_labels = y_train[idx] - local_purity = np.mean(neighbor_labels == y_test[:, None], axis=1) - confident = local_purity >= min_comp + _, idx = nn.kneighbors(X_test) + neighbor_labels = y_train[idx] + local_purity = np.mean(neighbor_labels == y_test[:, None], axis=1) + confident = local_purity >= min_comp adj_pred = (proba >= 0.5).astype(int) m = compute_extended_metrics(y_test, adj_pred, proba) - rows.append({ - "knora_k": k, - "min_competence": min_comp, - "fraction_confident_neighborhoods": float(np.mean(confident)), - "mcc": m["mcc"], - "f1": m["f1"], - "auprc": m["auprc"], - }) + rows.append( + { + "knora_k": k, + "min_competence": min_comp, + "fraction_confident_neighborhoods": float(np.mean(confident)), + "mcc": m["mcc"], + "f1": m["f1"], + "auprc": m["auprc"], + } + ) df = pd.DataFrame(rows) df.to_csv(output_dir / "parameter_sensitivity.csv", index=False) - pivot = df.pivot(index="knora_k", columns="min_competence", - values="fraction_confident_neighborhoods") - fig, ax = plt.subplots( - figsize=(RevisionConfig.SINGLE_COL * 1.3, RevisionConfig.SINGLE_COL)) + pivot = df.pivot( + index="knora_k", columns="min_competence", values="fraction_confident_neighborhoods" + ) + fig, ax = plt.subplots(figsize=(RevisionConfig.SINGLE_COL * 1.3, RevisionConfig.SINGLE_COL)) im = ax.imshow(pivot.values, cmap="viridis", aspect="auto") - ax.set_xticks(range(len(pivot.columns))); ax.set_xticklabels(pivot.columns) - ax.set_yticks(range(len(pivot.index))); ax.set_yticklabels(pivot.index) - ax.set_xlabel("min_competence"); ax.set_ylabel("knora_k") + ax.set_xticks(range(len(pivot.columns))) + ax.set_xticklabels(pivot.columns) + ax.set_yticks(range(len(pivot.index))) + ax.set_yticklabels(pivot.index) + ax.set_xlabel("min_competence") + ax.set_ylabel("knora_k") ax.set_title("Confident neighborhood fraction", fontweight="bold", loc="left") fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04) plt.tight_layout() - plt.savefig(output_dir / "Fig_Parameter_Sensitivity.png", - dpi=RevisionConfig.DPI, bbox_inches="tight") + plt.savefig( + output_dir / "Fig_Parameter_Sensitivity.png", dpi=RevisionConfig.DPI, bbox_inches="tight" + ) plt.close() return df -# --------------------------------------------------------------------------- -# Best-params table (updated for all v4 config keys — CHANGE 1, 2, 3) -# --------------------------------------------------------------------------- +# Best-params table + def save_best_params_table(output_dir): rows = [ # Core KNORA / neighbourhood parameters - {"hyperparameter": "knora_k", "value": RevisionConfig.KNORA_K, - "change": "FIX 3"}, - {"hyperparameter": "min_competence", "value": RevisionConfig.MIN_COMPETENCE_THRESHOLD, - "change": ""}, - {"hyperparameter": "tda_neighbors", "value": RevisionConfig.TDA_N_NEIGHBORS, - "change": ""}, - {"hyperparameter": "sequence_window", "value": RevisionConfig.SEQUENCE_WINDOW, - "change": ""}, - {"hyperparameter": "fcgr_k_values", "value": str(RevisionConfig.FCGR_K_VALUES), - "change": ""}, - {"hyperparameter": "tda_homology_dims", "value": str(RevisionConfig.TDA_HOMOLOGY_DIMS), - "change": ""}, - {"hyperparameter": "test_size", "value": RevisionConfig.TEST_SIZE, - "change": ""}, - {"hyperparameter": "random_state", "value": RevisionConfig.RANDOM_STATE, - "change": ""}, + {"hyperparameter": "knora_k", "value": RevisionConfig.KNORA_K, "change": "FIX 3"}, + { + "hyperparameter": "min_competence", + "value": RevisionConfig.MIN_COMPETENCE_THRESHOLD, + "change": "", + }, + {"hyperparameter": "tda_neighbors", "value": RevisionConfig.TDA_N_NEIGHBORS, "change": ""}, + { + "hyperparameter": "sequence_window", + "value": RevisionConfig.SEQUENCE_WINDOW, + "change": "", + }, + { + "hyperparameter": "fcgr_k_values", + "value": str(RevisionConfig.FCGR_K_VALUES), + "change": "", + }, + { + "hyperparameter": "tda_homology_dims", + "value": str(RevisionConfig.TDA_HOMOLOGY_DIMS), + "change": "", + }, + {"hyperparameter": "test_size", "value": RevisionConfig.TEST_SIZE, "change": ""}, + {"hyperparameter": "random_state", "value": RevisionConfig.RANDOM_STATE, "change": ""}, # Feature sizes - {"hyperparameter": "n_fcgr_features_raw", - "value": sum(4**k for k in RevisionConfig.FCGR_K_VALUES), - "change": "FIX 3"}, - {"hyperparameter": "n_fcgr_features_pca", "value": RevisionConfig.FCGR_PCA_COMPONENTS, - "change": "FIX 3"}, - {"hyperparameter": "n_tda_features", "value": RevisionConfig.TDA_N_FEATURES_EXPECTED, - "change": "CHANGE 2 (was 6)"}, - # FIX 2 - {"hyperparameter": "weak_learner_margin", "value": 0.10, - "change": "FIX 2"}, - # FIX 5 - {"hyperparameter": "threshold_search_low", "value": 0.45, - "change": "FIX 5"}, - {"hyperparameter": "threshold_search_high", "value": 0.65, - "change": "FIX 5"}, - # CHANGE 1: diversity pool - {"hyperparameter": "diversity_pool_additions", - "value": "MLP_BioTDA, kNN_SML, ElasticNet_Bio", - "change": "CHANGE 1"}, - # CHANGE 3: label ceiling parameters - {"hyperparameter": "selective_abstain_threshold", - "value": RevisionConfig.SELECTIVE_ABSTAIN_THRESHOLD, - "change": "CHANGE 3"}, - {"hyperparameter": "ambiguous_margin", - "value": RevisionConfig.AMBIGUOUS_MARGIN, - "change": "CHANGE 3"}, - {"hyperparameter": "ambiguous_sample_weight", - "value": RevisionConfig.AMBIGUOUS_SAMPLE_WEIGHT, - "change": "CHANGE 3"}, + { + "hyperparameter": "n_fcgr_features_raw", + "value": sum(4**k for k in RevisionConfig.FCGR_K_VALUES), + "change": "FIX 3", + }, + { + "hyperparameter": "n_fcgr_features_pca", + "value": RevisionConfig.FCGR_PCA_COMPONENTS, + "change": "FIX 3", + }, + { + "hyperparameter": "n_tda_features", + "value": RevisionConfig.TDA_N_FEATURES_EXPECTED, + "change": "CHANGE 2 (was 6)", + }, + {"hyperparameter": "weak_learner_margin", "value": 0.10, "change": "FIX 2"}, + {"hyperparameter": "threshold_search_low", "value": 0.45, "change": "FIX 5"}, + {"hyperparameter": "threshold_search_high", "value": 0.65, "change": "FIX 5"}, + # diversity pool + { + "hyperparameter": "diversity_pool_additions", + "value": "MLP_BioTDA, kNN_SML, ElasticNet_Bio", + "change": "CHANGE 1", + }, + # label ceiling parameters + { + "hyperparameter": "selective_abstain_threshold", + "value": RevisionConfig.SELECTIVE_ABSTAIN_THRESHOLD, + "change": "CHANGE 3", + }, + { + "hyperparameter": "ambiguous_margin", + "value": RevisionConfig.AMBIGUOUS_MARGIN, + "change": "CHANGE 3", + }, + { + "hyperparameter": "ambiguous_sample_weight", + "value": RevisionConfig.AMBIGUOUS_SAMPLE_WEIGHT, + "change": "CHANGE 3", + }, ] df = pd.DataFrame(rows) df.to_csv(output_dir / "best_hyperparameters.csv", index=False) return df -# --------------------------------------------------------------------------- # Runtime profile (unchanged) -# --------------------------------------------------------------------------- + def save_runtime_profile(output_dir, timings, n_samples, n_features): try: import psutil - mem_gb = round(psutil.virtual_memory().total / (1024 ** 3), 2) + + mem_gb = round(psutil.virtual_memory().total / (1024**3), 2) except ImportError: mem_gb = None - rows = [{ - "python_version": platform.python_version(), - "platform": platform.platform(), - "cpu_count": multiprocessing.cpu_count(), - "total_memory_gb": mem_gb, - "n_samples": n_samples, - "n_features": n_features, - }] + rows = [ + { + "python_version": platform.python_version(), + "platform": platform.platform(), + "cpu_count": multiprocessing.cpu_count(), + "total_memory_gb": mem_gb, + "n_samples": n_samples, + "n_features": n_features, + } + ] for step, secs in timings.items(): rows[0][f"time_{step}_sec"] = round(secs, 2) df = pd.DataFrame(rows) @@ -1149,9 +1223,8 @@ def save_runtime_profile(output_dir, timings, n_samples, n_features): return df -# --------------------------------------------------------------------------- -# main() — v3.0 (CHANGE 1 diversity step added) -# --------------------------------------------------------------------------- +# main — + def run_lime_analysis(model, X_train, X_test, y_test, output_dir, feature_names=None): try: @@ -1160,45 +1233,47 @@ def run_lime_analysis(model, X_train, X_test, y_test, output_dir, feature_names= except ImportError: print("LIME package not installed. Skipping LIME analysis.") return - + print("\n--- Running LIME Analysis on 4 Random Samples ---") - + if feature_names is None: feature_names = [f"feature_{i}" for i in range(X_train.shape[1])] - + explainer = lime.lime_tabular.LimeTabularExplainer( X_train, - mode='classification', + mode="classification", feature_names=feature_names, class_names=["Benign", "Pathogenic"], - random_state=42 + random_state=42, ) - + import numpy as np + np.random.seed(42) sample_indices = np.random.choice(len(X_test), 4, replace=False) - + y_pred = (model.predict_proba(X_test[sample_indices])[:, 1] >= 0.5).astype(int) y_test_arr = np.array(y_test) - + import matplotlib - matplotlib.use('Agg') + + matplotlib.use("Agg") import matplotlib.pyplot as plt - + fig, axes = plt.subplots(2, 2, figsize=(15, 12)) axes = axes.flatten() - + for i, idx in enumerate(sample_indices): true_label = "Pathogenic" if y_test_arr[idx] == 1 else "Benign" pred_label = "Pathogenic" if y_pred[i] == 1 else "Benign" case_name = f"Sample_{idx}_True_{true_label}_Pred_{pred_label}" - + print(f" Generating LIME explanation for {case_name}") exp = explainer.explain_instance(X_test[idx], model.predict_proba, num_features=10) - + # Save as HTML exp.save_to_file(str(output_dir / f"lime_explanation_{case_name}.html")) - + # Draw on subplot ax = axes[i] exp_list = exp.as_list() @@ -1206,14 +1281,14 @@ def run_lime_analysis(model, X_train, X_test, y_test, output_dir, feature_names= names = [x[0] for x in exp_list] vals.reverse() names.reverse() - colors = ['green' if x > 0 else 'red' for x in vals] + colors = ["green" if x > 0 else "red" for x in vals] pos = np.arange(len(exp_list)) + 0.5 - - ax.barh(pos, vals, align='center', color=colors) + + ax.barh(pos, vals, align="center", color=colors) ax.set_yticks(pos) ax.set_yticklabels(names) ax.set_title(f"Sample {idx}: True={true_label}, Pred={pred_label}", fontweight="bold") - ax.axvline(x=0, color='black', linestyle='--', linewidth=1.5, alpha=0.7) + ax.axvline(x=0, color="black", linestyle="--", linewidth=1.5, alpha=0.7) plt.tight_layout() plt.savefig(output_dir / "Fig_LIME_Combined_2x2.png", dpi=300, bbox_inches="tight") @@ -1221,7 +1296,7 @@ def run_lime_analysis(model, X_train, X_test, y_test, output_dir, feature_names= def main(): - out = RevisionConfig.OUTPUT_DIR + out = RevisionConfig.OUTPUT_DIR timings = {} t_start = time.time() @@ -1241,11 +1316,11 @@ def main(): full = build_full_feature_matrix(variant_df) if full is not None: - X = full["X"] - y = full["y"] + X = full["X"] + y = full["y"] n_standard = full["n_standard"] - n_fcgr = full["n_fcgr"] # FIX 3: PCA-compressed (20) - n_tda = full["n_tda"] # CHANGE 2: ~20 (over FCGR space) + n_fcgr = full["n_fcgr"] + n_tda = full["n_tda"] tab_fitted = None print( f"Using full feature matrix: {X.shape[1]} features " @@ -1256,14 +1331,14 @@ def main(): variant_df, RevisionConfig.TARGET_COL, RevisionConfig.LEAKAGE_COLS ) n_standard = X.shape[1] - n_fcgr = 0 - n_tda = 0 + n_fcgr = 0 + n_tda = 0 print( f"Using tabular feature matrix only: {X.shape[1]} features. " "Full FCGR/TDA ablation requires Step 8 preprocessor and genome file." ) - X = np.asarray(X, dtype=np.float32) + X = np.asarray(X, dtype=np.float32) n_total_features = X.shape[1] print(f"X dtype forced to float32. Memory: {X.nbytes / 1e9:.2f} GB") @@ -1271,17 +1346,18 @@ def main(): verify_consensus_score(variant_df, out) X_train, X_test, y_train, y_test = train_test_split( - X, y, test_size=RevisionConfig.TEST_SIZE, - stratify=y, random_state=RevisionConfig.RANDOM_STATE, + X, + y, + test_size=RevisionConfig.TEST_SIZE, + stratify=y, + random_state=RevisionConfig.RANDOM_STATE, ) X_train_bio = X_train[:, :n_standard] - X_test_bio = X_test[:, :n_standard] + X_test_bio = X_test[:, :n_standard] - # ------------------------------------------------------------------ - # CHANGE 1: Diversity analysis — read Code 8 v4 output CSVs - # ------------------------------------------------------------------ - print("\nStep: Diversity Analysis (CHANGE 1)") + # Diversity analysis — read script 08 output CSVs + print("\nStep: Diversity Analysis") t0 = time.time() diversity_summary = run_diversity_analysis(out) timings["diversity_analysis"] = time.time() - t0 @@ -1294,9 +1370,16 @@ def main(): plot_combined_roc_pr(curve_data, out) run_mcnemar_matrix(fitted_models, X_test_bio, y_test, "XGBoost", out) - + if "XGBoost" in fitted_models: - run_lime_analysis(fitted_models["XGBoost"], X_train_bio, X_test_bio, y_test, out, feature_names=tab_fitted["feature_names"]) + run_lime_analysis( + fitted_models["XGBoost"], + X_train_bio, + X_test_bio, + y_test, + out, + feature_names=tab_fitted["feature_names"], + ) import gc @@ -1308,15 +1391,11 @@ def main(): gc.collect() t0 = time.time() - run_parameter_sensitivity( - X_train, y_train, X_test, y_test, n_standard, n_fcgr, n_tda, out - ) + run_parameter_sensitivity(X_train, y_train, X_test, y_test, n_standard, n_fcgr, n_tda, out) timings["parameter_sensitivity"] = time.time() - t0 t0 = time.time() - run_feature_ablation( - X_train, y_train, X_test, y_test, n_standard, n_fcgr, n_tda, out - ) + run_feature_ablation(X_train, y_train, X_test, y_test, n_standard, n_fcgr, n_tda, out) timings["feature_ablation"] = time.time() - t0 best_name = max( @@ -1329,8 +1408,11 @@ def main(): run_threshold_sweep(fitted_models[best_name], X_test_bio, y_test, out) bio_model = HistGradientBoostingClassifier( - max_iter=200, max_depth=6, learning_rate=0.05, - random_state=RevisionConfig.RANDOM_STATE, early_stopping=False, + max_iter=200, + max_depth=6, + learning_rate=0.05, + random_state=RevisionConfig.RANDOM_STATE, + early_stopping=False, ) bio_model.fit(X_train_bio, y_train) diff --git a/src/10_exploratory_data_analysis.py b/src/10_exploratory_data_analysis.py index 15a055b..02ea56e 100644 --- a/src/10_exploratory_data_analysis.py +++ b/src/10_exploratory_data_analysis.py @@ -1,8 +1,11 @@ +"""Exploratory data analysis: figures and summary statistics for the balanced dataset.""" + import pandas as pd from config import STAGE07_OUT, EDA_OUT import numpy as np import matplotlib -matplotlib.use('Agg') + +matplotlib.use("Agg") import matplotlib.pyplot as plt import seaborn as sns from scipy import stats @@ -36,13 +39,9 @@ # Load the dataset print("\n[1/10] Loading Dataset...") df = pd.read_csv(STAGE07_OUT / "Final_Dataset_Balanced.csv") -print( - f"✓ Dataset loaded successfully: {df.shape[0]:,} variants × {df.shape[1]} features" -) +print(f"✓ Dataset loaded successfully: {df.shape[0]:,} variants × {df.shape[1]} features") -# ============================================================================ # SECTION 1: DATASET OVERVIEW AND DESCRIPTIVE STATISTICS -# ============================================================================ print("\n" + "=" * 80) print("SECTION 1: DATASET OVERVIEW") print("=" * 80) @@ -100,9 +99,7 @@ print(f"✓ Dataset summary saved to {output_dir}/01_dataset_summary.txt") -# ============================================================================ # SECTION 2: TARGET VARIABLE ANALYSIS -# ============================================================================ print("\n" + "=" * 80) print("SECTION 2: TARGET VARIABLE DISTRIBUTION") print("=" * 80) @@ -115,12 +112,8 @@ f.write("=" * 80 + "\n\n") f.write("Class Distribution:\n") f.write("-" * 40 + "\n") - f.write( - f"Pathogenic variants (1): {target_counts.get(1, 0):,} ({target_pct.get(1, 0):.2f}%)\n" - ) - f.write( - f"Benign variants (0): {target_counts.get(0, 0):,} ({target_pct.get(0, 0):.2f}%)\n" - ) + f.write(f"Pathogenic variants (1): {target_counts.get(1, 0):,} ({target_pct.get(1, 0):.2f}%)\n") + f.write(f"Benign variants (0): {target_counts.get(0, 0):,} ({target_pct.get(0, 0):.2f}%)\n") f.write( f"\nClass ratio (Benign:Pathogenic): {target_counts.get(0, 0)/target_counts.get(1, 1):.2f}:1\n" ) @@ -142,9 +135,7 @@ ) axes[0].set_ylabel("Number of Variants", fontsize=11, fontweight="bold") axes[0].set_xlabel("Variant Class", fontsize=11, fontweight="bold") -axes[0].set_title( - "A) Distribution of Variant Classes", fontsize=12, fontweight="bold", pad=15 -) +axes[0].set_title("A) Distribution of Variant Classes", fontsize=12, fontweight="bold", pad=15) max_val = max(target_counts.get(0, 0), target_counts.get(1, 0)) axes[0].set_ylim(0, max_val * 1.2) axes[0].ticklabel_format(style="plain", axis="y") @@ -173,22 +164,16 @@ textprops={"fontsize": 11, "fontweight": "bold"}, wedgeprops={"edgecolor": "black", "linewidth": 1.5}, ) -axes[1].set_title( - "B) Proportion of Variant Classes", fontsize=12, fontweight="bold", pad=15 -) +axes[1].set_title("B) Proportion of Variant Classes", fontsize=12, fontweight="bold", pad=15) plt.tight_layout() -plt.savefig( - f"{EDA_OUT}/Figure1_Target_Distribution.png", dpi=300, bbox_inches="tight" -) +plt.savefig(f"{EDA_OUT}/Figure1_Target_Distribution.png", dpi=300, bbox_inches="tight") plt.savefig(f"{EDA_OUT}/Figure1_Target_Distribution.pdf", bbox_inches="tight") plt.close() print(f"✓ Figure 1 saved: Target distribution") -# ============================================================================ # SECTION 3: CHROMOSOMAL DISTRIBUTION -# ============================================================================ print("\n" + "=" * 80) print("SECTION 3: CHROMOSOMAL DISTRIBUTION ANALYSIS") print("=" * 80) @@ -272,17 +257,13 @@ axes[1].grid(False) plt.tight_layout() -plt.savefig( - f"{EDA_OUT}/Figure2_Chromosomal_Distribution.png", dpi=300, bbox_inches="tight" -) +plt.savefig(f"{EDA_OUT}/Figure2_Chromosomal_Distribution.png", dpi=300, bbox_inches="tight") plt.savefig(f"{EDA_OUT}/Figure2_Chromosomal_Distribution.pdf", bbox_inches="tight") plt.close() print(f"✓ Figure 2 saved: Chromosomal distribution") -# ============================================================================ # SECTION 4: NUCLEOTIDE SUBSTITUTION PATTERNS -# ============================================================================ print("\n" + "=" * 80) print("SECTION 4: NUCLEOTIDE SUBSTITUTION ANALYSIS") print("=" * 80) @@ -337,9 +318,7 @@ axes[0].set_yticks(range(len(top_subst))) axes[0].set_yticklabels(top_subst.index, fontsize=9) axes[0].set_xlabel("Number of Variants", fontsize=11, fontweight="bold") -axes[0].set_title( - "A) Top 12 Nucleotide Substitutions", fontsize=12, fontweight="bold", pad=15 -) +axes[0].set_title("A) Top 12 Nucleotide Substitutions", fontsize=12, fontweight="bold", pad=15) axes[0].grid(False) axes[0].invert_yaxis() @@ -353,22 +332,16 @@ linewidth=1.5, ) axes[1].set_ylabel("Number of Variants", fontsize=11, fontweight="bold") -axes[1].set_title( - "B) Transition vs Transversion", fontsize=12, fontweight="bold", pad=15 -) +axes[1].set_title("B) Transition vs Transversion", fontsize=12, fontweight="bold", pad=15) max_val2 = max(mutation_counts.get("Transition", 0), mutation_counts.get("Transversion", 0)) axes[1].set_ylim(0, max_val2 * 1.15) axes[1].grid(False) for i, k in enumerate(["Transition", "Transversion"]): v = mutation_counts.get(k, 0) - axes[1].text( - i, v + 1000, f"{v:,}", ha="center", va="bottom", fontweight="bold", fontsize=10 - ) + axes[1].text(i, v + 1000, f"{v:,}", ha="center", va="bottom", fontweight="bold", fontsize=10) # Panel C: Ti/Tv by pathogenicity -mutation_by_class_pct = ( - mutation_by_class.div(mutation_by_class.sum(axis=1), axis=0) * 100 -) +mutation_by_class_pct = mutation_by_class.div(mutation_by_class.sum(axis=1), axis=0) * 100 x_pos = np.arange(2) width = 0.35 axes[2].bar( @@ -392,26 +365,20 @@ linewidth=1, ) axes[2].set_ylabel("Percentage (%)", fontsize=11, fontweight="bold") -axes[2].set_title( - "C) Mutation Type by Pathogenicity", fontsize=12, fontweight="bold", pad=15 -) +axes[2].set_title("C) Mutation Type by Pathogenicity", fontsize=12, fontweight="bold", pad=15) axes[2].set_xticks(x_pos) axes[2].set_xticklabels(["Transition", "Transversion"]) axes[2].legend(frameon=True, fancybox=True, shadow=True) axes[2].grid(False) plt.tight_layout() -plt.savefig( - f"{EDA_OUT}/Figure3_Substitution_Patterns.png", dpi=300, bbox_inches="tight" -) +plt.savefig(f"{EDA_OUT}/Figure3_Substitution_Patterns.png", dpi=300, bbox_inches="tight") plt.savefig(f"{EDA_OUT}/Figure3_Substitution_Patterns.pdf", bbox_inches="tight") plt.close() print(f"✓ Figure 3 saved: Substitution patterns") -# ============================================================================ # SECTION 5: PATHOGENICITY SCORES DISTRIBUTION -# ============================================================================ print("\n" + "=" * 80) print("SECTION 5: PATHOGENICITY PREDICTION SCORES") print("=" * 80) @@ -453,9 +420,7 @@ ) score_stats_df = pd.DataFrame(score_stats) -score_stats_df.to_csv( - f"{EDA_OUT}/05_pathogenicity_scores_statistics.csv", index=False -) +score_stats_df.to_csv(f"{EDA_OUT}/05_pathogenicity_scores_statistics.csv", index=False) print(f"✓ Score statistics saved") # Figure 4: Score distributions @@ -496,9 +461,7 @@ axes[idx].set_xticks([0, 1]) axes[idx].set_xticklabels(["Benign", "Pathogenic"]) axes[idx].set_ylabel("Score Value", fontsize=10, fontweight="bold") - axes[idx].set_title( - f"{chr(65+idx)}) {score}", fontsize=11, fontweight="bold", pad=10 - ) + axes[idx].set_title(f"{chr(65+idx)}) {score}", fontsize=11, fontweight="bold", pad=10) axes[idx].grid(False) # Add p-value annotation @@ -523,16 +486,12 @@ dpi=300, bbox_inches="tight", ) -plt.savefig( - f"{EDA_OUT}/Figure4_Pathogenicity_Scores_Distribution.pdf", bbox_inches="tight" -) +plt.savefig(f"{EDA_OUT}/Figure4_Pathogenicity_Scores_Distribution.pdf", bbox_inches="tight") plt.close() print(f"✓ Figure 4 saved: Pathogenicity scores distribution") -# ============================================================================ # SECTION 6: CONSERVATION SCORES -# ============================================================================ print("\n" + "=" * 80) print("SECTION 6: EVOLUTIONARY CONSERVATION ANALYSIS") print("=" * 80) @@ -588,17 +547,13 @@ ) plt.tight_layout() -plt.savefig( - f"{EDA_OUT}/Figure5_Conservation_Scores.png", dpi=300, bbox_inches="tight" -) +plt.savefig(f"{EDA_OUT}/Figure5_Conservation_Scores.png", dpi=300, bbox_inches="tight") plt.savefig(f"{EDA_OUT}/Figure5_Conservation_Scores.pdf", bbox_inches="tight") plt.close() print(f"✓ Figure 5 saved: Conservation scores") -# ============================================================================ # SECTION 7: CLINICAL ANNOTATIONS -# ============================================================================ print("\n" + "=" * 80) print("SECTION 7: CLINICAL ANNOTATION ANALYSIS") print("=" * 80) @@ -638,9 +593,7 @@ ) clinical_stats_df = pd.DataFrame(clinical_stats) -clinical_stats_df.to_csv( - f"{EDA_OUT}/06_clinical_annotations_statistics.csv", index=False -) +clinical_stats_df.to_csv(f"{EDA_OUT}/06_clinical_annotations_statistics.csv", index=False) print(f"✓ Clinical annotation statistics saved") # Figure 6: Clinical annotations @@ -650,9 +603,7 @@ for idx, feature in enumerate(["IS_CANCER_GENE", "IS_ONCOGENE", "IS_TSG", "TIER"]): if feature == "TIER": # TIER distribution - tier_dist = ( - pd.crosstab(df["TIER"], df["LABEL_PATHOGENIC"], normalize="columns") * 100 - ) + tier_dist = pd.crosstab(df["TIER"], df["LABEL_PATHOGENIC"], normalize="columns") * 100 tier_dist.T.plot( kind="bar", ax=axes[idx], @@ -674,9 +625,7 @@ axes[idx].grid(False) else: # Binary features - contingency = ( - pd.crosstab(df["LABEL_PATHOGENIC"], df[feature], normalize="index") * 100 - ) + contingency = pd.crosstab(df["LABEL_PATHOGENIC"], df[feature], normalize="index") * 100 x_pos = np.arange(2) width = 0.35 @@ -703,18 +652,14 @@ axes[idx].set_ylabel("Percentage (%)", fontsize=11, fontweight="bold") axes[idx].set_xlabel("Variant Class", fontsize=11, fontweight="bold") - axes[idx].set_title( - f"{chr(65+idx)}) {feature}", fontsize=12, fontweight="bold", pad=15 - ) + axes[idx].set_title(f"{chr(65+idx)}) {feature}", fontsize=12, fontweight="bold", pad=15) axes[idx].set_xticks(x_pos) axes[idx].set_xticklabels(["Benign", "Pathogenic"]) axes[idx].legend(frameon=True, fancybox=True, shadow=True) axes[idx].grid(False) # Add p-value - pval = clinical_stats_df[clinical_stats_df["Feature"] == feature][ - "P_value" - ].values[0] + pval = clinical_stats_df[clinical_stats_df["Feature"] == feature]["P_value"].values[0] pval_text = f"p < 0.001" if pval < 0.001 else f"p = {pval:.3f}" axes[idx].text( 0.95, @@ -729,17 +674,13 @@ ) plt.tight_layout() -plt.savefig( - f"{EDA_OUT}/Figure6_Clinical_Annotations.png", dpi=300, bbox_inches="tight" -) +plt.savefig(f"{EDA_OUT}/Figure6_Clinical_Annotations.png", dpi=300, bbox_inches="tight") plt.savefig(f"{EDA_OUT}/Figure6_Clinical_Annotations.pdf", bbox_inches="tight") plt.close() print(f"✓ Figure 6 saved: Clinical annotations") -# ============================================================================ # SECTION 8: STRUCTURAL FEATURES -# ============================================================================ print("\n" + "=" * 80) print("SECTION 8: PROTEIN STRUCTURAL FEATURES") print("=" * 80) @@ -790,9 +731,7 @@ axes[idx].set_xlabel("Score Value", fontsize=11, fontweight="bold") axes[idx].set_ylabel("Density", fontsize=11, fontweight="bold") - axes[idx].set_title( - f"{chr(65+idx)}) {feature}", fontsize=12, fontweight="bold", pad=15 - ) + axes[idx].set_title(f"{chr(65+idx)}) {feature}", fontsize=12, fontweight="bold", pad=15) axes[idx].legend(frameon=True, fancybox=True, shadow=True) axes[idx].grid(False) @@ -817,9 +756,7 @@ dpi=300, bbox_inches="tight", ) -plt.savefig( - f"{EDA_OUT}/Figure7_Structural_Features_Continuous.pdf", bbox_inches="tight" -) +plt.savefig(f"{EDA_OUT}/Figure7_Structural_Features_Continuous.pdf", bbox_inches="tight") plt.close() print(f"✓ Figure 7 saved: Structural features (continuous)") @@ -829,9 +766,7 @@ axes = axes.flatten() for idx, feature in enumerate(structural_binary): - contingency = ( - pd.crosstab(df["LABEL_PATHOGENIC"], df[feature], normalize="index") * 100 - ) + contingency = pd.crosstab(df["LABEL_PATHOGENIC"], df[feature], normalize="index") * 100 x_pos = np.arange(2) width = 0.35 @@ -858,9 +793,7 @@ axes[idx].set_ylabel("Percentage (%)", fontsize=11, fontweight="bold") axes[idx].set_xlabel("Variant Class", fontsize=11, fontweight="bold") - axes[idx].set_title( - f"{chr(65+idx)}) {feature}", fontsize=12, fontweight="bold", pad=15 - ) + axes[idx].set_title(f"{chr(65+idx)}) {feature}", fontsize=12, fontweight="bold", pad=15) axes[idx].set_xticks(x_pos) axes[idx].set_xticklabels(["Benign", "Pathogenic"]) axes[idx].legend(frameon=True, fancybox=True, shadow=True) @@ -883,27 +816,20 @@ ) plt.tight_layout() -plt.savefig( - f"{EDA_OUT}/Figure8_Structural_Features_Binary.png", dpi=300, bbox_inches="tight" -) +plt.savefig(f"{EDA_OUT}/Figure8_Structural_Features_Binary.png", dpi=300, bbox_inches="tight") plt.savefig(f"{EDA_OUT}/Figure8_Structural_Features_Binary.pdf", bbox_inches="tight") plt.close() print(f"✓ Figure 8 saved: Structural features (binary)") -# ============================================================================ # SECTION 9: CORRELATION ANALYSIS -# ============================================================================ print("\n" + "=" * 80) print("SECTION 9: FEATURE CORRELATION ANALYSIS") print("=" * 80) # Select numerical features for correlation numerical_features = ( - pathogenicity_scores - + conservation_scores - + structural_continuous - + ["LABEL_PATHOGENIC"] + pathogenicity_scores + conservation_scores + structural_continuous + ["LABEL_PATHOGENIC"] ) correlation_matrix = df[numerical_features].corr() @@ -935,17 +861,13 @@ plt.xticks(rotation=45, ha="right", fontsize=14) plt.yticks(rotation=0, fontsize=14) plt.tight_layout() -plt.savefig( - f"{EDA_OUT}/Figure9_Correlation_Heatmap.png", dpi=300, bbox_inches="tight" -) +plt.savefig(f"{EDA_OUT}/Figure9_Correlation_Heatmap.png", dpi=300, bbox_inches="tight") plt.savefig(f"{EDA_OUT}/Figure9_Correlation_Heatmap.pdf", bbox_inches="tight") plt.close() print(f"✓ Figure 9 saved: Correlation heatmap") -# ============================================================================ # SECTION 10: COMPREHENSIVE SUMMARY STATISTICS -# ============================================================================ print("\n" + "=" * 80) print("SECTION 10: COMPREHENSIVE SUMMARY STATISTICS") print("=" * 80) @@ -975,15 +897,9 @@ f.write("-" * 80 + "\n") f.write(f"Total Variants: {len(df):,}\n") f.write(f"Features: {df.shape[1]}\n") - f.write( - f"Pathogenic Variants: {target_counts.get(1, 0):,} ({target_pct.get(1, 0):.2f}%)\n" - ) - f.write( - f"Benign Variants: {target_counts.get(0, 0):,} ({target_pct.get(0, 0):.2f}%)\n" - ) - f.write( - f"Class Balance Ratio: {target_counts.get(0, 0)/target_counts.get(1, 1):.2f}:1\n\n" - ) + f.write(f"Pathogenic Variants: {target_counts.get(1, 0):,} ({target_pct.get(1, 0):.2f}%)\n") + f.write(f"Benign Variants: {target_counts.get(0, 0):,} ({target_pct.get(0, 0):.2f}%)\n") + f.write(f"Class Balance Ratio: {target_counts.get(0, 0)/target_counts.get(1, 1):.2f}:1\n\n") f.write("GENOMIC CHARACTERISTICS\n") f.write("-" * 80 + "\n") @@ -1001,30 +917,20 @@ f.write( f"Variants in Oncogenes: {df['IS_ONCOGENE'].sum():,} ({df['IS_ONCOGENE'].sum()/len(df)*100:.2f}%)\n" ) - f.write( - f"Variants in TSGs: {df['IS_TSG'].sum():,} ({df['IS_TSG'].sum()/len(df)*100:.2f}%)\n" - ) + f.write(f"Variants in TSGs: {df['IS_TSG'].sum():,} ({df['IS_TSG'].sum()/len(df)*100:.2f}%)\n") f.write( f"Tier 1 Variants: {(df['TIER']==1).sum():,} ({(df['TIER']==1).sum()/len(df)*100:.2f}%)\n\n" ) f.write("KEY FINDINGS\n") f.write("-" * 80 + "\n") - f.write( - "1. All pathogenicity prediction scores show significant differences between\n" - ) + f.write("1. All pathogenicity prediction scores show significant differences between\n") f.write(" benign and pathogenic variants (p < 0.001)\n\n") - f.write( - "2. Conservation scores (GERP++, phyloP) demonstrate strong discriminative power\n" - ) + f.write("2. Conservation scores (GERP++, phyloP) demonstrate strong discriminative power\n") f.write(" with pathogenic variants showing higher conservation\n\n") - f.write( - "3. Clinical annotations show significant enrichment in pathogenic variants,\n" - ) + f.write("3. Clinical annotations show significant enrichment in pathogenic variants,\n") f.write(" particularly for cancer gene annotations\n\n") - f.write( - "4. Structural features reveal distinct patterns with pathogenic variants\n" - ) + f.write("4. Structural features reveal distinct patterns with pathogenic variants\n") f.write(" showing preferential localization in functional protein regions\n\n") f.write("OUTPUT FILES GENERATED\n") diff --git a/src/config.py b/src/config.py index 858ec38..dda7cdd 100644 --- a/src/config.py +++ b/src/config.py @@ -22,9 +22,19 @@ STAGE07_OUT = OUTPUT_DIR / "07_balancing" STAGE08_OUT = OUTPUT_DIR / "08_tda_fuzzy" STAGE09_OUT = OUTPUT_DIR / "09_experiments" -EDA_OUT = OUTPUT_DIR / "eda_figures" +EDA_OUT = OUTPUT_DIR / "eda_figures" # Ensure all output directories exist when config is imported -for d in [STAGE01_OUT, STAGE02_OUT, STAGE03_OUT, STAGE04_OUT, STAGE05_OUT, - STAGE06_OUT, STAGE07_OUT, STAGE08_OUT, STAGE09_OUT, EDA_OUT]: +for d in [ + STAGE01_OUT, + STAGE02_OUT, + STAGE03_OUT, + STAGE04_OUT, + STAGE05_OUT, + STAGE06_OUT, + STAGE07_OUT, + STAGE08_OUT, + STAGE09_OUT, + EDA_OUT, +]: d.mkdir(parents=True, exist_ok=True)