Update app.py

app.py

@@ -1,19 +1,5 @@
-import joblib
-_scaler = joblib.load("Comb1_aac_ctd_RFE_selected_features_scaler.joblib")
-print("SCALER n_features_in_:", getattr(_scaler, "n_features_in_", "N/A"), flush=True)
-_names = getattr(_scaler, "feature_names_in_", None)
-if _names is not None:
-    print("SCALER FEATURE NAMES (%d):" % len(_names), flush=True)
-    print(list(_names), flush=True)
-else:
-    print("SCALER has NO feature_names_in_ (fit on numpy array)", flush=True)
-
-from tensorflow.keras.models import load_model
-_m = load_model("Comb1_aac_ctd_RFE_selected_features_model.keras")
-print("MODEL input_shape:", _m.input_shape, "output_shape:", _m.output_shape, flush=True)
-
 import os
-# Quiet TensorFlow logs and disable oneDNN nondeterminism notice
+# Quiet TensorFlow logs (must be set before importing tensorflow)
 os.environ.setdefault("TF_CPP_MIN_LOG_LEVEL", "3")
 os.environ.setdefault("TF_ENABLE_ONEDNN_OPTS", "0")
 
@@ -25,22 +11,18 @@ from propy import AAComposition, CTD
 from math import expm1
 
 # ---------------------------------------------------------------------------
-# LAZY LOADING
-# On the free 16GB Space, loading TensorFlow + PyTorch + ProtBert all at once
-# at import time causes an out-of-memory crash. We therefore load each heavy
-# component only when it is first needed, and cache it after that.
+# LAZY LOADING — keeps the free 16GB Space from OOM-ing at startup.
+# Heavy libs (TF, torch, ProtBert) load only when first needed.
 # ---------------------------------------------------------------------------
-
-_amp_model = None          # Keras AMP classifier
-_amp_scaler = None         # joblib scaler for AMP features
+_amp_model = None
+_amp_scaler = None
 _protbert_tokenizer = None
 _protbert_model = None
-_torch = None              # torch module, imported lazily
+_torch = None
 _device = None
 
 
 def get_amp_model():
-    """Load the Keras AMP classifier + scaler on first use."""
     global _amp_model, _amp_scaler
     if _amp_model is None:
         from tensorflow.keras.models import load_model
@@ -50,7 +32,6 @@ def get_amp_model():
 
 
 def get_protbert():
-    """Load ProtBert tokenizer + model on first use (only needed for MIC)."""
     global _protbert_tokenizer, _protbert_model, _torch, _device
     if _protbert_model is None:
         import torch
@@ -66,8 +47,9 @@ def get_protbert():
 
 
 # ---------------------------------------------------------------------------
-# Selected features (AAC + CTD, RFE-selected). 'Activity' is the target label
-# and is intentionally excluded from the input features.
+# The EXACT 343 features the scaler was fit on, IN THE EXACT TRAINING ORDER.
+# The scaler was fit on a numpy array (no stored names), so order is critical:
+# we must select these columns in this order BEFORE calling scaler.transform().
 # ---------------------------------------------------------------------------
 selected_features = [
     "_PolarizabilityC1", "_PolarizabilityC2", "_PolarizabilityC3",
@@ -77,7 +59,6 @@ selected_features = [
     "_PolarityC1", "_PolarityC2", "_PolarityC3",
     "_NormalizedVDWVC1", "_NormalizedVDWVC2", "_NormalizedVDWVC3",
     "_HydrophobicityC1", "_HydrophobicityC2", "_HydrophobicityC3",
-
     "_PolarizabilityT12", "_PolarizabilityT13", "_PolarizabilityT23",
     "_SolventAccessibilityT12", "_SolventAccessibilityT13", "_SolventAccessibilityT23",
     "_SecondaryStrT12", "_SecondaryStrT13", "_SecondaryStrT23",
@@ -85,14 +66,12 @@ selected_features = [
     "_PolarityT12", "_PolarityT13", "_PolarityT23",
     "_NormalizedVDWVT12", "_NormalizedVDWVT13", "_NormalizedVDWVT23",
     "_HydrophobicityT12", "_HydrophobicityT13", "_HydrophobicityT23",
-
     "_PolarizabilityD1001", "_PolarizabilityD1025", "_PolarizabilityD1050",
     "_PolarizabilityD1075", "_PolarizabilityD1100",
     "_PolarizabilityD2001", "_PolarizabilityD2025", "_PolarizabilityD2050",
     "_PolarizabilityD2075", "_PolarizabilityD2100",
     "_PolarizabilityD3001", "_PolarizabilityD3025", "_PolarizabilityD3050",
     "_PolarizabilityD3075", "_PolarizabilityD3100",
-
     "_SolventAccessibilityD1001", "_SolventAccessibilityD1025",
     "_SolventAccessibilityD1050", "_SolventAccessibilityD1075",
     "_SolventAccessibilityD1100",
@@ -102,28 +81,24 @@ selected_features = [
     "_SolventAccessibilityD3001", "_SolventAccessibilityD3025",
     "_SolventAccessibilityD3050", "_SolventAccessibilityD3075",
     "_SolventAccessibilityD3100",
-
     "_SecondaryStrD1001", "_SecondaryStrD1025", "_SecondaryStrD1050",
     "_SecondaryStrD1075", "_SecondaryStrD1100",
     "_SecondaryStrD2001", "_SecondaryStrD2025", "_SecondaryStrD2050",
     "_SecondaryStrD2075", "_SecondaryStrD2100",
     "_SecondaryStrD3001", "_SecondaryStrD3025", "_SecondaryStrD3050",
     "_SecondaryStrD3075", "_SecondaryStrD3100",
-
     "_ChargeD1001", "_ChargeD1025", "_ChargeD1050",
     "_ChargeD1075", "_ChargeD1100",
     "_ChargeD2001", "_ChargeD2025", "_ChargeD2050",
     "_ChargeD2075",
     "_ChargeD3001", "_ChargeD3025", "_ChargeD3050",
     "_ChargeD3075", "_ChargeD3100",
-
     "_PolarityD1001", "_PolarityD1025", "_PolarityD1050",
     "_PolarityD1075", "_PolarityD1100",
     "_PolarityD2001", "_PolarityD2025", "_PolarityD2050",
     "_PolarityD2075", "_PolarityD2100",
     "_PolarityD3001", "_PolarityD3025", "_PolarityD3050",
     "_PolarityD3075", "_PolarityD3100",
-
     "_NormalizedVDWVD1001", "_NormalizedVDWVD1025",
     "_NormalizedVDWVD1050", "_NormalizedVDWVD1075",
     "_NormalizedVDWVD1100",
@@ -133,7 +108,6 @@ selected_features = [
     "_NormalizedVDWVD3001", "_NormalizedVDWVD3025",
     "_NormalizedVDWVD3050", "_NormalizedVDWVD3075",
     "_NormalizedVDWVD3100",
-
     "_HydrophobicityD1001", "_HydrophobicityD1025",
     "_HydrophobicityD1050", "_HydrophobicityD1075",
     "_HydrophobicityD1100",
@@ -143,10 +117,8 @@ selected_features = [
     "_HydrophobicityD3001", "_HydrophobicityD3025",
     "_HydrophobicityD3050", "_HydrophobicityD3075",
     "_HydrophobicityD3100",
-
     "A", "R", "N", "D", "C", "E", "Q", "G", "H", "I",
     "L", "K", "M", "F", "P", "S", "T", "W", "Y", "V",
-
     "AR", "AD", "AQ", "AG", "AL", "AK", "AF", "AP", "AT", "AV",
     "RA", "RC", "RE", "RG", "RI", "RL", "RS", "RT", "RV",
     "NR", "NC", "NG", "NI", "NP", "NS", "NY", "NV",
@@ -171,21 +143,21 @@ selected_features = [
     "VA", "VR", "VD", "VC", "VE", "VG", "VI", "VL", "VK",
     "VS", "VT", "VY", "VV"
 ]
+assert len(selected_features) == 343, f"Expected 343 features, got {len(selected_features)}"
 
 
 def keras_predict_proba(X):
-    """Return probabilities as [P(Non-AMP), P(AMP)] for LIME."""
+    """Return probabilities as [P(Non-AMP), P(AMP)] for LIME (X already scaled)."""
     amp_model, _ = get_amp_model()
     preds = amp_model.predict(X, verbose=0)
     if preds.ndim == 1 or preds.shape[1] == 1:
         preds = preds.reshape(-1, 1)
-        # Assuming sigmoid output = P(AMP); adjust if your model is reversed.
-        return np.hstack([1 - preds, preds])
+        return np.hstack([1 - preds, preds])  # sigmoid output assumed = P(AMP)
     return preds
 
 
 def extract_features(sequence):
-    """Compute AAC (420) + CTD features, scale, and select RFE features."""
+    """Compute CTD + AAC, select the 343 training columns IN ORDER, then scale."""
     sequence = ''.join([aa for aa in sequence.upper() if aa in "ACDEFGHIKLMNPQRSTVWY"])
     if len(sequence) < 10:
         return "Error: Sequence too short."
@@ -193,33 +165,31 @@ def extract_features(sequence):
     try:
         _, amp_scaler = get_amp_model()
 
-        # AAC: 20 single AAs + 400 dipeptides = 420 features
-        dipeptide_features = AAComposition.CalculateAADipeptideComposition(sequence)
-        filtered_aac = {k: dipeptide_features[k] for k in list(dipeptide_features.keys())[:420]}
-
-        # CTD: Composition, Transition, Distribution
+        # Compute full feature pool
         ctd_features = CTD.CalculateCTD(sequence)
+        aac = AAComposition.CalculateAADipeptideComposition(sequence)
 
-        all_features_dict = {}
-        all_features_dict.update(ctd_features)
-        all_features_dict.update(filtered_aac)
+        # Merge everything into one lookup dict
+        pool = {}
+        pool.update(ctd_features)
+        pool.update(aac)
 
-        feature_df_all = pd.DataFrame([all_features_dict])
-        normalized_array = amp_scaler.transform(feature_df_all.values)
-        normalized_df = pd.DataFrame(normalized_array, columns=feature_df_all.columns)
+        # Verify all needed features are present
+        missing = [f for f in selected_features if f not in pool]
+        if missing:
+            return f"Error: Missing features from propy: {missing[:5]}..."
 
-        if not set(selected_features).issubset(normalized_df.columns):
-            missing = set(selected_features) - set(normalized_df.columns)
-            return f"Error: Missing features: {list(missing)[:5]}..."
+        # Build the 343-wide row IN THE EXACT TRAINING ORDER, THEN scale.
+        ordered_values = [pool[f] for f in selected_features]
+        feature_row = np.array(ordered_values, dtype=np.float64).reshape(1, -1)
 
-        selected_df = normalized_df[selected_features].fillna(0)
-        return selected_df.values.astype(np.float32)
+        scaled = amp_scaler.transform(feature_row)  # scaler expects exactly 343 cols
+        return scaled.astype(np.float32)
     except Exception as e:
         return f"Error in feature extraction: {str(e)}"
 
 
 def predictmic(sequence):
-    """Predict MIC values using ProtBert embeddings + per-bacterium models."""
     sequence = ''.join([aa for aa in sequence.upper() if aa in "ACDEFGHIKLMNPQRSTVWY"])
     if len(sequence) < 10:
         return {"Error": "Sequence too short or invalid."}
@@ -266,21 +236,20 @@ def full_prediction(sequence):
     amp_model, _ = get_amp_model()
     raw_pred = amp_model.predict(features, verbose=0)
 
-    # Handle sigmoid (1 output) vs softmax (>=2 outputs)
     if raw_pred.ndim == 1 or raw_pred.shape[1] == 1:
-        prob_amp = float(raw_pred.flatten()[0])  # assume output = P(AMP)
+        prob_amp = float(raw_pred.flatten()[0])  # sigmoid output assumed = P(AMP)
         if prob_amp >= 0.5:
-            prediction = 1  # AMP
+            prediction = 1
             confidence = round(prob_amp * 100, 2)
         else:
-            prediction = 0  # Non-AMP
+            prediction = 0
             confidence = round((1 - prob_amp) * 100, 2)
     else:
         class_idx = int(np.argmax(raw_pred[0]))
         prediction = class_idx
         confidence = round(float(raw_pred[0][class_idx]) * 100, 2)
 
-    # Label convention: 1 = AMP, 0 = Non-AMP (swap if your model is opposite)
+    # Label convention: 1 = AMP, 0 = Non-AMP (swap if your model is reversed)
     amp_result = "Antimicrobial Peptide (AMP)" if prediction == 1 else "Non-AMP"
     result = f"Prediction: {amp_result}\nConfidence: {confidence}%\n"
 
@@ -292,7 +261,6 @@ def full_prediction(sequence):
     else:
         result += "\nMIC prediction skipped for Non-AMP sequences.\n"
 
-    # LIME explanation (lazy import keeps startup light)
     try:
         from lime.lime_tabular import LimeTabularExplainer
         sample_data = np.random.rand(100, len(selected_features))
@@ -316,7 +284,6 @@ def full_prediction(sequence):
     return result
 
 
-# Gradio UI
 iface = gr.Interface(
     fn=full_prediction,
     inputs=gr.Textbox(label="Enter Protein Sequence"),