Code quality, documentation, and repo hygiene cleanup

.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/

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.

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"] = (

src/02_remove_missing_values.py

@@ -1,88 +1,93 @@
-#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:
         print(f"\nError occurred: {e}")
 
 
 if __name__ == "__main__":
-    clean_dataset_selective()
+    clean_dataset_selective()