Spaces:
Running
Running
| import os | |
| # Native-lib hygiene (prevents TF/PyTorch SIGSEGV when both load; harmless for RF) | |
| os.environ.setdefault("KMP_DUPLICATE_LIB_OK", "TRUE") | |
| os.environ.setdefault("OMP_NUM_THREADS", "1") | |
| os.environ.setdefault("MKL_NUM_THREADS", "1") | |
| os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") | |
| os.environ.setdefault("TOKENIZERS_PARALLELISM", "false") | |
| import sys | |
| import json | |
| import subprocess | |
| import joblib | |
| import numpy as np | |
| import pandas as pd | |
| from propy import AAComposition, Autocorrelation, CTD, PseudoAAC | |
| from lime.lime_tabular import LimeTabularExplainer | |
| import gradio as gr | |
| # --------------------------------------------------------------------------- | |
| # Load Random Forest AMP classifier + MinMax scaler (original files) | |
| # --------------------------------------------------------------------------- | |
| model = joblib.load("RF.joblib") | |
| scaler = joblib.load("norm (4).joblib") | |
| # --------------------------------------------------------------------------- | |
| # Original 138 RFE-selected features (CTD + AAC + Autocorrelation + APAAC) | |
| # --------------------------------------------------------------------------- | |
| selected_features = [ | |
| "_SolventAccessibilityC3", "_SecondaryStrC1", "_SecondaryStrC3", "_ChargeC1", "_PolarityC1", | |
| "_NormalizedVDWVC1", "_HydrophobicityC3", "_SecondaryStrT23", "_PolarizabilityD1001", "_PolarizabilityD2001", | |
| "_PolarizabilityD3001", "_SolventAccessibilityD1001", "_SolventAccessibilityD2001", "_SolventAccessibilityD3001", | |
| "_SecondaryStrD1001", "_SecondaryStrD1075", "_SecondaryStrD2001", "_SecondaryStrD3001", "_ChargeD1001", | |
| "_ChargeD1025", "_ChargeD2001", "_ChargeD3075", "_ChargeD3100", "_PolarityD1001", "_PolarityD1050", | |
| "_PolarityD2001", "_PolarityD3001", "_NormalizedVDWVD1001", "_NormalizedVDWVD2001", "_NormalizedVDWVD2025", | |
| "_NormalizedVDWVD2050", "_NormalizedVDWVD3001", "_HydrophobicityD1001", "_HydrophobicityD2001", | |
| "_HydrophobicityD3001", "_HydrophobicityD3025", "A", "R", "D", "C", "E", "Q", "H", "I", "M", "P", "Y", "V", | |
| "AR", "AV", "RC", "RL", "RV", "CR", "CC", "CL", "CK", "EE", "EI", "EL", "HC", "IA", "IL", "IV", "LA", "LC", "LE", | |
| "LI", "LT", "LV", "KC", "MA", "MS", "SC", "TC", "TV", "YC", "VC", "VE", "VL", "VK", "VV", | |
| "MoreauBrotoAuto_FreeEnergy30", "MoranAuto_Hydrophobicity2", "MoranAuto_Hydrophobicity4", | |
| "GearyAuto_Hydrophobicity20", "GearyAuto_Hydrophobicity24", "GearyAuto_Hydrophobicity26", | |
| "GearyAuto_Hydrophobicity27", "GearyAuto_Hydrophobicity28", "GearyAuto_Hydrophobicity29", | |
| "GearyAuto_Hydrophobicity30", "GearyAuto_AvFlexibility22", "GearyAuto_AvFlexibility26", | |
| "GearyAuto_AvFlexibility27", "GearyAuto_AvFlexibility28", "GearyAuto_AvFlexibility29", "GearyAuto_AvFlexibility30", | |
| "GearyAuto_Polarizability22", "GearyAuto_Polarizability24", "GearyAuto_Polarizability25", | |
| "GearyAuto_Polarizability27", "GearyAuto_Polarizability28", "GearyAuto_Polarizability29", | |
| "GearyAuto_Polarizability30", "GearyAuto_FreeEnergy24", "GearyAuto_FreeEnergy25", "GearyAuto_FreeEnergy30", | |
| "GearyAuto_ResidueASA21", "GearyAuto_ResidueASA22", "GearyAuto_ResidueASA23", "GearyAuto_ResidueASA24", | |
| "GearyAuto_ResidueASA30", "GearyAuto_ResidueVol21", "GearyAuto_ResidueVol24", "GearyAuto_ResidueVol25", | |
| "GearyAuto_ResidueVol26", "GearyAuto_ResidueVol28", "GearyAuto_ResidueVol29", "GearyAuto_ResidueVol30", | |
| "GearyAuto_Steric18", "GearyAuto_Steric21", "GearyAuto_Steric26", "GearyAuto_Steric27", "GearyAuto_Steric28", | |
| "GearyAuto_Steric29", "GearyAuto_Steric30", "GearyAuto_Mutability23", "GearyAuto_Mutability25", | |
| "GearyAuto_Mutability26", "GearyAuto_Mutability27", "GearyAuto_Mutability28", "GearyAuto_Mutability29", | |
| "GearyAuto_Mutability30", "APAAC1", "APAAC4", "APAAC5", "APAAC6", "APAAC8", "APAAC9", "APAAC12", "APAAC13", | |
| "APAAC15", "APAAC18", "APAAC19", "APAAC24" | |
| ] | |
| assert len(selected_features) == 138, f"Expected 138 features, got {len(selected_features)}" | |
| # --------------------------------------------------------------------------- | |
| # LIME explainer | |
| # Built ONCE at startup so explanations are reproducible across requests. | |
| # The training-data argument controls how LIME perturbs features around the | |
| # input. After MinMax scaling each feature lives in [0,1], so we use a small | |
| # uniform sample with a FIXED seed — that gives stable, repeatable weights. | |
| # (If you have a saved sample of real normalized training rows, swap it in | |
| # here and explanations will reflect the true feature distribution.) | |
| # --------------------------------------------------------------------------- | |
| _rng = np.random.default_rng(seed=42) | |
| _lime_background = _rng.uniform(low=0.0, high=1.0, size=(500, len(selected_features))) | |
| explainer = LimeTabularExplainer( | |
| training_data=_lime_background, | |
| feature_names=selected_features, | |
| class_names=["AMP", "Non-AMP"], | |
| mode="classification", | |
| discretize_continuous=True, | |
| random_state=42, # stable explanations | |
| ) | |
| # --------------------------------------------------------------------------- | |
| # Feature extraction — produces the full propy feature pool, scales it with | |
| # the saved MinMax scaler, then selects the 138 features the RF was trained on. | |
| # --------------------------------------------------------------------------- | |
| def extract_features(sequence): | |
| sequence = ''.join([aa for aa in sequence.upper() if aa in "ACDEFGHIKLMNPQRSTVWY"]) | |
| if len(sequence) < 10: | |
| return "Error: Sequence too short." | |
| try: | |
| # Original full pool: CTD + AAC(first 420) + Autocorrelation + PseudoAAC | |
| dipeptide_features = AAComposition.CalculateAADipeptideComposition(sequence) | |
| filtered_dipeptide_features = {k: dipeptide_features[k] for k in list(dipeptide_features.keys())[:420]} | |
| ctd_features = CTD.CalculateCTD(sequence) | |
| auto_features = Autocorrelation.CalculateAutoTotal(sequence) | |
| pseudo_features = PseudoAAC.GetAPseudoAAC(sequence, lamda=9) | |
| all_features_dict = {} | |
| all_features_dict.update(ctd_features) | |
| all_features_dict.update(filtered_dipeptide_features) | |
| all_features_dict.update(auto_features) | |
| all_features_dict.update(pseudo_features) | |
| feature_df_all = pd.DataFrame([all_features_dict]) | |
| normalized_array = scaler.transform(feature_df_all.values) | |
| normalized_df = pd.DataFrame(normalized_array, columns=feature_df_all.columns) | |
| 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]}..." | |
| selected_df = normalized_df[selected_features].fillna(0) | |
| return selected_df.values | |
| except Exception as e: | |
| return f"Error in feature extraction: {str(e)}" | |
| # --------------------------------------------------------------------------- | |
| # MIC prediction — runs in a SEPARATE process (mic_worker.py). | |
| # This isolates PyTorch/ProtBert from the main process and prevents the | |
| # native-library crash (exit 139) plus the OOM spike on the free tier. | |
| # --------------------------------------------------------------------------- | |
| def predictmic(sequence): | |
| sequence = ''.join([aa for aa in sequence.upper() if aa in "ACDEFGHIKLMNPQRSTVWY"]) | |
| if len(sequence) < 10: | |
| return {"Error": "Sequence too short or invalid."} | |
| try: | |
| proc = subprocess.run( | |
| [sys.executable, "mic_worker.py", sequence], | |
| capture_output=True, text=True, timeout=900 | |
| ) | |
| except subprocess.TimeoutExpired: | |
| return {"Error": "MIC prediction timed out (ProtBert may still be downloading; try again shortly)."} | |
| except Exception as e: | |
| return {"Error": f"Failed to start MIC worker: {str(e)}"} | |
| if proc.returncode != 0: | |
| tail = (proc.stderr or "").strip().splitlines()[-3:] | |
| return {"Error": f"MIC worker exited with code {proc.returncode}. {' '.join(tail)}"} | |
| out_lines = [ln for ln in (proc.stdout or "").splitlines() if ln.strip()] | |
| if not out_lines: | |
| return {"Error": "MIC worker produced no output."} | |
| try: | |
| return json.loads(out_lines[-1]) | |
| except Exception: | |
| return {"Error": f"Could not parse MIC worker output: {out_lines[-1][:200]}"} | |
| # --------------------------------------------------------------------------- | |
| # Main prediction pipeline | |
| # --------------------------------------------------------------------------- | |
| def full_prediction(sequence): | |
| features = extract_features(sequence) | |
| if isinstance(features, str): | |
| return features | |
| prediction = model.predict(features)[0] | |
| probabilities = model.predict_proba(features)[0] | |
| try: | |
| class_index = list(model.classes_).index(prediction) | |
| confidence = round(probabilities[class_index] * 100, 2) | |
| except Exception: | |
| confidence = "Unknown" | |
| amp_result = "Antimicrobial Peptide (AMP)" if prediction == 0 else "Non-AMP" | |
| result = f"Prediction: {amp_result}\nConfidence: {confidence}%\n" | |
| # ---- LIME first (per your spec: LIME before SHAP in the report) ---- | |
| try: | |
| explanation = explainer.explain_instance( | |
| data_row=features[0], # <-- explicitly the single input sequence | |
| predict_fn=model.predict_proba, | |
| num_features=10, | |
| num_samples=2000, # perturbations around this single input | |
| ) | |
| result += "\nTop Features Influencing Prediction (LIME):\n" | |
| for feat, weight in explanation.as_list(): | |
| result += f"- {feat}: {round(weight, 4)}\n" | |
| except Exception as e: | |
| result += f"\nLIME explanation failed: {str(e)}\n" | |
| # ---- MIC (only for AMPs) ---- | |
| if prediction == 0: | |
| mic_values = predictmic(sequence) | |
| result += "\nPredicted MIC Values (μM):\n" | |
| for org, mic in mic_values.items(): | |
| result += f"- {org}: {mic}\n" | |
| else: | |
| result += "\nMIC prediction skipped for Non-AMP sequences.\n" | |
| return result | |
| # Gradio UI | |
| iface = gr.Interface( | |
| fn=full_prediction, | |
| inputs=gr.Textbox(label="Enter Protein Sequence"), | |
| outputs=gr.Textbox(label="Results"), | |
| title="AMP & MIC Predictor + LIME Explanation", | |
| description="Paste an amino acid sequence (≥10 characters). Get AMP classification, MIC predictions, and LIME interpretability insights." | |
| ) | |
| iface.launch() |