nonzeroexit commited on
Commit
d750a2d
·
verified ·
1 Parent(s): b8e2131

Update app.py

Browse files
Files changed (1) hide show
  1. app.py +405 -552
app.py CHANGED
@@ -1,157 +1,71 @@
1
- import gradio as gr
2
  import joblib
3
  import numpy as np
4
  import pandas as pd
 
5
  import torch
6
 
7
- from propy import AAComposition, Autocorrelation, CTD, PseudoAAC
8
- from transformers import BertTokenizer, BertModel
9
- from lime.lime_tabular import LimeTabularExplainer
10
  from math import expm1
11
 
 
 
 
 
12
 
13
- # ============================================================
14
- # CONFIG
15
- # ============================================================
16
-
17
- MODEL_PATH = "RF.joblib"
18
- SCALER_PATH = "norm (4).joblib"
19
 
20
- # REAL training data used by the 146-feature classifier
21
- LIME_DATA_PATH = "lime_training_data.csv"
22
 
23
 
24
  # ============================================================
25
- # LOAD CLASSIFIER
26
  # ============================================================
27
 
28
- print("Loading Random Forest...")
29
-
30
- model = joblib.load(MODEL_PATH)
31
-
32
- print(
33
- "Random Forest loaded."
34
- )
35
-
36
- print(
37
- "RF expects:",
38
- getattr(model, "n_features_in_", "unknown"),
39
- "features"
40
- )
41
 
42
 
43
  # ============================================================
44
- # LOAD SCALER
45
  # ============================================================
46
 
47
- print("Loading scaler...")
48
-
49
- amp_scaler = joblib.load(SCALER_PATH)
50
-
51
- print(
52
- "Scaler loaded."
53
- )
54
-
55
- print(
56
- "Scaler expects:",
57
- getattr(
58
- amp_scaler,
59
- "n_features_in_",
60
- "unknown"
61
- ),
62
- "features"
63
- )
64
 
65
 
66
  # ============================================================
67
- # YOUR 146 FEATURES
68
  # ============================================================
69
 
70
  selected_features = [
71
- "_SolventAccessibilityC3",
72
- "_SecondaryStrC1",
73
- "_SecondaryStrC3",
74
- "_ChargeC1",
75
- "_PolarityC1",
76
- "_NormalizedVDWVC1",
77
- "_HydrophobicityC3",
78
- "_SecondaryStrT23",
79
- "_PolarizabilityD1001",
80
- "_PolarizabilityD2001",
81
- "_PolarizabilityD3001",
82
  "_SolventAccessibilityD1001",
83
  "_SolventAccessibilityD2001",
84
- "_SolventAccessibilityD3001",
85
- "_SecondaryStrD1001",
86
- "_SecondaryStrD1075",
87
- "_SecondaryStrD2001",
88
- "_SecondaryStrD3001",
89
- "_ChargeD1001",
90
- "_ChargeD1025",
91
- "_ChargeD2001",
92
- "_ChargeD3075",
93
- "_ChargeD3100",
94
- "_PolarityD1001",
95
- "_PolarityD1050",
96
- "_PolarityD2001",
97
- "_PolarityD3001",
98
- "_NormalizedVDWVD1001",
99
- "_NormalizedVDWVD2001",
100
- "_NormalizedVDWVD2025",
101
- "_NormalizedVDWVD2050",
102
- "_NormalizedVDWVD3001",
103
- "_HydrophobicityD1001",
104
- "_HydrophobicityD2001",
105
- "_HydrophobicityD3001",
106
  "_HydrophobicityD3025",
107
 
108
- "A",
109
- "R",
110
- "D",
111
- "C",
112
- "E",
113
- "Q",
114
- "H",
115
- "I",
116
- "M",
117
- "P",
118
- "Y",
119
- "V",
120
-
121
- "AR",
122
- "AV",
123
- "RC",
124
- "RL",
125
- "RV",
126
- "CR",
127
- "CC",
128
- "CL",
129
- "CK",
130
- "EE",
131
- "EI",
132
- "EL",
133
- "HC",
134
- "IA",
135
- "IL",
136
- "IV",
137
- "LA",
138
- "LC",
139
- "LE",
140
- "LI",
141
- "LT",
142
- "LV",
143
- "KC",
144
- "MA",
145
- "MS",
146
- "SC",
147
- "TC",
148
- "TV",
149
- "YC",
150
- "VC",
151
- "VE",
152
- "VL",
153
- "VK",
154
- "VV",
155
 
156
  "MoreauBrotoAuto_FreeEnergy30",
157
 
@@ -210,176 +124,292 @@ selected_features = [
210
  "GearyAuto_Mutability23",
211
  "GearyAuto_Mutability25",
212
  "GearyAuto_Mutability26",
213
- "GearyAuto_Mutability27",
214
  "GearyAuto_Mutability28",
215
  "GearyAuto_Mutability29",
216
  "GearyAuto_Mutability30",
217
 
218
- "APAAC1",
219
- "APAAC4",
220
- "APAAC5",
221
- "APAAC6",
222
- "APAAC8",
223
- "APAAC9",
224
- "APAAC12",
225
- "APAAC13",
226
- "APAAC15",
227
- "APAAC18",
228
- "APAAC19",
229
- "APAAC24"
230
  ]
231
 
232
 
 
 
 
 
 
 
 
 
233
  print(
234
- "Selected feature count:",
235
  len(selected_features)
236
  )
237
 
238
 
239
  # ============================================================
240
- # VERIFY CLASSIFIER
241
  # ============================================================
242
 
243
- if hasattr(model, "n_features_in_"):
244
-
245
- if model.n_features_in_ != len(selected_features):
246
-
247
- raise ValueError(
248
- "\n\nMODEL / FEATURE MISMATCH\n"
249
- f"RF expects {model.n_features_in_} features.\n"
250
- f"Your selected feature list contains "
251
- f"{len(selected_features)} features.\n\n"
252
- "Do not continue until these match."
253
- )
254
 
 
 
 
255
 
256
- # ============================================================
257
- # LOAD REAL LIME TRAINING DATA
258
- # ============================================================
259
 
260
  print(
261
- "\nLoading REAL training data for LIME..."
262
- )
263
-
264
- training_df = pd.read_csv(
265
- LIME_DATA_PATH
266
  )
267
 
268
 
269
- # Remove accidental pandas index
270
- if "Unnamed: 0" in training_df.columns:
 
271
 
272
- training_df = training_df.drop(
273
- "Unnamed: 0",
274
- axis=1
275
  )
276
 
277
 
 
 
 
 
 
 
 
 
 
 
 
 
278
  print(
279
- "Training CSV shape:",
280
- training_df.shape
 
281
  )
282
 
283
 
284
  # ============================================================
285
- # CHECK ALL 146 FEATURES EXIST
286
  # ============================================================
287
 
288
- missing_features = [
289
- feature
290
- for feature in selected_features
291
- if feature not in training_df.columns
292
- ]
293
 
 
 
 
 
294
 
295
- if missing_features:
296
 
297
- raise ValueError(
298
- "\n\nLIME TRAINING DATA IS MISSING FEATURES:\n"
299
- + "\n".join(
300
- missing_features[:20]
 
 
301
  )
 
 
 
 
302
  )
303
 
304
 
305
- # ============================================================
306
- # PREPARE REAL LIME DATA
307
- # ============================================================
308
 
309
- lime_training_df = (
310
- training_df[
311
- selected_features
312
- ]
313
- .copy()
314
- )
315
 
316
 
317
- # Convert to numbers
318
- lime_training_df = (
319
- lime_training_df.apply(
320
- pd.to_numeric,
321
- errors="coerce"
322
  )
323
- )
324
 
325
 
326
- # Missing numerical values
327
- lime_training_df = (
328
- lime_training_df.fillna(0)
329
- )
330
 
 
331
 
332
- # Convert to numpy
333
- lime_training_data = (
334
- lime_training_df
335
- .values
336
- .astype(np.float64)
337
- )
338
 
 
 
 
339
 
340
- print(
341
- "Real LIME data:",
342
- lime_training_data.shape
343
- )
344
 
 
 
 
345
 
346
- # ============================================================
347
- # LIMIT LIME BACKGROUND
348
- # ============================================================
349
 
350
- # This is NOT generating artificial data.
351
- # It only selects real rows from your training dataset.
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
352
 
353
- MAX_LIME_BACKGROUND = 1000
354
 
 
 
 
355
 
356
- if len(lime_training_data) > MAX_LIME_BACKGROUND:
357
 
358
- rng = np.random.RandomState(42)
 
 
 
 
359
 
360
- indices = rng.choice(
361
- len(lime_training_data),
362
- MAX_LIME_BACKGROUND,
363
- replace=False
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
364
  )
365
 
366
- lime_training_data = (
367
- lime_training_data[
368
- indices
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
369
  ]
 
 
 
370
  )
371
 
372
 
373
- print(
374
- "LIME background:",
375
- lime_training_data.shape
376
- )
 
 
 
 
 
 
 
 
 
 
 
377
 
378
 
379
  # ============================================================
380
  # CREATE LIME EXPLAINER
381
  # ============================================================
382
 
 
 
 
 
 
383
  explainer = LimeTabularExplainer(
384
 
385
  training_data=lime_training_data,
@@ -399,53 +429,44 @@ explainer = LimeTabularExplainer(
399
  )
400
 
401
 
402
- print(
403
- "REAL LIME EXPLAINER READY."
404
- )
405
 
406
 
407
  # ============================================================
408
- # LOAD PROTBERT
409
  # ============================================================
410
 
411
- print(
412
- "\nLoading ProtBert..."
413
- )
414
-
415
 
416
  tokenizer = BertTokenizer.from_pretrained(
417
  "Rostlab/prot_bert",
418
  do_lower_case=False
419
  )
420
 
421
-
422
  protbert_model = BertModel.from_pretrained(
423
  "Rostlab/prot_bert"
424
  )
425
 
426
-
427
  device = torch.device(
428
  "cuda"
429
  if torch.cuda.is_available()
430
  else "cpu"
431
  )
432
 
433
-
434
  protbert_model = (
435
  protbert_model
436
  .to(device)
437
  .eval()
438
  )
439
 
440
-
441
  print(
442
- "ProtBert loaded on:",
443
  device
444
  )
445
 
446
 
447
  # ============================================================
448
- # VALID AMINO ACIDS
449
  # ============================================================
450
 
451
  VALID_AA = set(
@@ -453,29 +474,21 @@ VALID_AA = set(
453
  )
454
 
455
 
456
- # ============================================================
457
- # CLEAN SEQUENCE
458
- # ============================================================
459
-
460
  def clean_sequence(sequence):
461
 
462
  sequence = str(
463
  sequence
464
  ).upper().strip()
465
 
466
-
467
- sequence = "".join(
468
  aa
469
  for aa in sequence
470
  if aa in VALID_AA
471
  )
472
 
473
 
474
- return sequence
475
-
476
-
477
  # ============================================================
478
- # FEATURE EXTRACTION
479
  # ============================================================
480
 
481
  def extract_features(sequence):
@@ -484,65 +497,43 @@ def extract_features(sequence):
484
  sequence
485
  )
486
 
487
-
488
  if len(sequence) < 10:
489
 
490
- return (
491
- None,
492
- "Sequence must contain at least 10 valid amino acids."
493
- )
494
 
495
 
496
  try:
497
 
498
- # ----------------------------------------------------
499
- # DIPEPTIDE
500
- # ----------------------------------------------------
501
-
502
- dipeptide_features = (
503
- AAComposition.CalculateAADipeptideComposition(
504
  sequence
505
  )
506
  )
507
 
508
 
509
- filtered_dipeptide_features = {
510
-
511
- k: dipeptide_features[k]
512
-
513
  for k in list(
514
- dipeptide_features.keys()
515
  )[:420]
516
  }
517
 
518
 
519
- # ----------------------------------------------------
520
- # CTD
521
- # ----------------------------------------------------
522
-
523
- ctd_features = (
524
- CTD.CalculateCTD(
525
- sequence
526
- )
527
  )
528
 
529
 
530
- # ----------------------------------------------------
531
- # AUTOCORRELATION
532
- # ----------------------------------------------------
533
-
534
- auto_features = (
535
- Autocorrelation.CalculateAutoTotal(
536
  sequence
537
  )
538
  )
539
 
540
 
541
- # ----------------------------------------------------
542
- # PSEUDO AAC
543
- # ----------------------------------------------------
544
-
545
- pseudo_features = (
546
  PseudoAAC.GetAPseudoAAC(
547
  sequence,
548
  lamda=9
@@ -550,27 +541,12 @@ def extract_features(sequence):
550
  )
551
 
552
 
553
- # ----------------------------------------------------
554
- # COMBINE
555
- # ----------------------------------------------------
556
-
557
  all_features = {}
558
 
559
- all_features.update(
560
- ctd_features
561
- )
562
-
563
- all_features.update(
564
- filtered_dipeptide_features
565
- )
566
-
567
- all_features.update(
568
- auto_features
569
- )
570
-
571
- all_features.update(
572
- pseudo_features
573
- )
574
 
575
 
576
  feature_df = pd.DataFrame(
@@ -578,34 +554,32 @@ def extract_features(sequence):
578
  )
579
 
580
 
581
- # ----------------------------------------------------
582
- # ALIGN WITH SCALER
583
- # ----------------------------------------------------
584
 
585
  if hasattr(
586
- amp_scaler,
587
  "feature_names_in_"
588
  ):
589
 
590
  scaler_features = list(
591
- amp_scaler.feature_names_in_
592
  )
593
 
594
 
595
- missing_scaler = [
596
  f
597
  for f in scaler_features
598
  if f not in feature_df.columns
599
  ]
600
 
601
 
602
- if missing_scaler:
603
 
604
  return (
605
  None,
606
- "Missing features required by scaler: "
607
  + str(
608
- missing_scaler[:10]
609
  )
610
  )
611
 
@@ -617,69 +591,52 @@ def extract_features(sequence):
617
  )
618
 
619
 
620
- # ----------------------------------------------------
621
- # NORMALIZATION
622
- # ----------------------------------------------------
623
 
624
- normalized_array = (
625
- amp_scaler.transform(
626
- feature_df
627
- )
628
  )
629
 
630
 
631
  normalized_df = pd.DataFrame(
632
- normalized_array,
633
  columns=feature_df.columns
634
  )
635
 
636
 
637
- # ----------------------------------------------------
638
- # SELECT THE SAME 146 FEATURES
639
- # ----------------------------------------------------
640
 
641
- missing_selected = [
642
  f
643
  for f in selected_features
644
  if f not in normalized_df.columns
645
  ]
646
 
647
 
648
- if missing_selected:
649
 
650
  return (
651
  None,
652
  "Missing selected features: "
653
- + str(
654
- missing_selected[:10]
655
- )
656
  )
657
 
658
 
659
- selected_df = (
660
  normalized_df[
661
  selected_features
662
  ]
663
  .fillna(0)
664
- )
665
-
666
-
667
- X = (
668
- selected_df
669
  .values
670
- .astype(np.float64)
671
  )
672
 
673
 
674
- # Final check
675
- if X.shape[1] != len(
676
- selected_features
677
- ):
678
 
679
  return (
680
  None,
681
- f"Feature count error: "
682
- f"{X.shape[1]} generated."
683
  )
684
 
685
 
@@ -695,7 +652,7 @@ def extract_features(sequence):
695
 
696
 
697
  # ============================================================
698
- # MIC PREDICTION
699
  # ============================================================
700
 
701
  def predict_mic(sequence):
@@ -705,17 +662,6 @@ def predict_mic(sequence):
705
  )
706
 
707
 
708
- if len(sequence) < 10:
709
-
710
- return {
711
- "Error": "Sequence too short."
712
- }
713
-
714
-
715
- # --------------------------------------------------------
716
- # PROTBERT
717
- # --------------------------------------------------------
718
-
719
  seq_spaced = " ".join(
720
  list(sequence)
721
  )
@@ -746,8 +692,7 @@ def predict_mic(sequence):
746
 
747
 
748
  embedding = (
749
- outputs
750
- .last_hidden_state
751
  .mean(dim=1)
752
  .squeeze()
753
  .cpu()
@@ -756,77 +701,55 @@ def predict_mic(sequence):
756
  )
757
 
758
 
759
- # --------------------------------------------------------
760
- # BACTERIAL MODELS
761
- # --------------------------------------------------------
762
-
763
  bacteria_config = {
764
 
765
- "E. coli": {
766
- "model":
767
- "coli_xgboost_model.pkl",
768
-
769
- "scaler":
770
- "coli_scaler.pkl",
771
-
772
- "pca":
773
- None
774
- },
775
-
776
-
777
- "S. aureus": {
778
- "model":
779
- "aur_xgboost_model.pkl",
780
-
781
- "scaler":
782
- "aur_scaler.pkl",
783
-
784
- "pca":
785
- None
786
- },
787
-
788
-
789
- "P. aeruginosa": {
790
- "model":
791
- "arg_xgboost_model.pkl",
792
-
793
- "scaler":
794
- "arg_scaler.pkl",
795
-
796
- "pca":
797
- None
798
- },
799
-
800
 
801
- "K. pneumoniae": {
802
- "model":
803
- "pne_mlp_model.pkl",
 
 
804
 
805
- "scaler":
806
- "pne_scaler.pkl",
 
 
 
807
 
808
- "pca":
809
- "pne_pca.pkl"
810
- }
 
 
811
  }
812
 
813
 
814
- mic_results = {}
815
 
816
 
817
- for bacterium, config in (
818
  bacteria_config.items()
819
  ):
820
 
821
  try:
822
 
823
- mic_scaler = joblib.load(
824
- config["scaler"]
825
  )
826
 
827
 
828
  mic_model = joblib.load(
829
- config["model"]
 
 
 
 
 
830
  )
831
 
832
 
@@ -837,15 +760,14 @@ def predict_mic(sequence):
837
  )
838
 
839
 
840
- if config["pca"]:
841
 
842
- mic_pca = joblib.load(
843
- config["pca"]
844
  )
845
 
846
-
847
  transformed = (
848
- mic_pca.transform(
849
  scaled
850
  )
851
  )
@@ -855,60 +777,37 @@ def predict_mic(sequence):
855
  transformed = scaled
856
 
857
 
858
- mic_log = (
859
  mic_model.predict(
860
  transformed
861
  )[0]
862
  )
863
 
864
 
865
- mic = round(
866
- expm1(mic_log),
 
 
867
  3
868
  )
869
 
870
 
871
- mic_results[
872
- bacterium
873
- ] = mic
874
-
875
-
876
  except Exception as e:
877
 
878
- mic_results[
879
- bacterium
880
- ] = (
881
  f"Error: {str(e)}"
882
  )
883
 
884
 
885
- return mic_results
886
 
887
 
888
  # ============================================================
889
- # FULL PREDICTION
890
  # ============================================================
891
 
892
  def full_prediction(sequence):
893
 
894
- sequence = clean_sequence(
895
- sequence
896
- )
897
-
898
-
899
- if len(sequence) < 10:
900
-
901
- return (
902
- "Invalid sequence.\n\n"
903
- "Please enter at least 10 valid "
904
- "amino acids."
905
- )
906
-
907
-
908
- # ========================================================
909
- # EXTRACT FEATURES
910
- # ========================================================
911
-
912
  features, error = (
913
  extract_features(
914
  sequence
@@ -924,109 +823,69 @@ def full_prediction(sequence):
924
  )
925
 
926
 
927
- # ========================================================
928
- # RANDOM FOREST PREDICTION
929
- # ========================================================
930
-
931
- try:
932
-
933
- prediction = (
934
- model.predict(
935
- features
936
- )[0]
937
- )
938
-
939
-
940
- probabilities = (
941
- model.predict_proba(
942
- features
943
- )[0]
944
- )
945
-
946
-
947
- except Exception as e:
948
-
949
- return (
950
- "AMP prediction failed:\n\n"
951
- f"{type(e).__name__}: {e}"
952
- )
953
-
954
-
955
- # ========================================================
956
- # CONFIDENCE
957
- # ========================================================
958
-
959
- try:
960
-
961
- class_index = list(
962
- model.classes_
963
- ).index(
964
- prediction
965
- )
966
 
 
 
 
967
 
968
- confidence = round(
969
- probabilities[
970
- class_index
971
- ] * 100,
972
- 2
973
- )
974
 
 
 
 
 
 
975
 
976
- except Exception:
977
-
978
- confidence = None
979
 
 
 
 
 
 
980
 
981
- # ========================================================
982
- # AMP LABEL
983
- # ========================================================
984
 
985
- if prediction == 0:
 
 
 
986
 
987
- amp_result = (
988
- "Antimicrobial Peptide (AMP)"
989
- )
990
 
991
- else:
992
-
993
- amp_result = (
994
- "Non-AMP"
995
- )
 
996
 
997
 
998
  result = (
999
- "====================================\n"
1000
  "AMP CLASSIFICATION\n"
1001
- "====================================\n\n"
1002
- f"Prediction: {amp_result}\n"
1003
  f"Confidence: {confidence}%\n"
1004
  )
1005
 
1006
 
1007
- # ========================================================
1008
  # LIME
1009
- # ========================================================
1010
 
1011
  result += (
1012
- "\n====================================\n"
1013
- "LIME INTERPRETABILITY\n"
1014
- "====================================\n"
1015
  )
1016
 
1017
 
1018
  try:
1019
 
1020
- # Reproducibility
1021
- np.random.seed(42)
1022
-
1023
-
1024
  explanation = (
1025
  explainer.explain_instance(
1026
 
1027
- data_row=features[0],
1028
 
1029
- predict_fn=model.predict_proba,
1030
 
1031
  labels=[
1032
  int(prediction)
@@ -1034,35 +893,28 @@ def full_prediction(sequence):
1034
 
1035
  num_features=10,
1036
 
1037
- num_samples=3000
1038
- )
1039
- )
1040
-
1041
-
1042
- lime_results = (
1043
- explanation.as_list(
1044
- label=int(prediction)
1045
  )
1046
  )
1047
 
1048
 
1049
  result += (
1050
- "\nTop features contributing "
1051
- "to this prediction:\n"
1052
  )
1053
 
1054
 
1055
  for feature, weight in (
1056
- lime_results
 
 
1057
  ):
1058
 
1059
- if weight > 0:
1060
-
1061
- direction = "↑"
1062
-
1063
- else:
1064
-
1065
- direction = "↓"
1066
 
1067
 
1068
  result += (
@@ -1075,21 +927,21 @@ def full_prediction(sequence):
1075
  except Exception as e:
1076
 
1077
  result += (
1078
- "\n\nLIME failed:\n"
1079
- f"{type(e).__name__}: {e}"
 
1080
  )
1081
 
1082
 
1083
- # ========================================================
1084
  # MIC
1085
- # ========================================================
1086
 
1087
  if prediction == 0:
1088
 
1089
  result += (
1090
- "\n\n====================================\n"
1091
- "PREDICTED MIC VALUES\n"
1092
- "====================================\n"
1093
  )
1094
 
1095
 
@@ -1112,7 +964,7 @@ def full_prediction(sequence):
1112
 
1113
  result += (
1114
  "\n\nMIC prediction skipped "
1115
- "for Non-AMP sequences."
1116
  )
1117
 
1118
 
@@ -1129,10 +981,10 @@ iface = gr.Interface(
1129
 
1130
  inputs=gr.Textbox(
1131
  label="Protein Sequence",
 
1132
  placeholder=(
1133
  "Enter an amino acid sequence..."
1134
- ),
1135
- lines=5
1136
  ),
1137
 
1138
  outputs=gr.Textbox(
@@ -1141,15 +993,14 @@ iface = gr.Interface(
1141
  ),
1142
 
1143
  title=(
1144
- "AMP & MIC Predictor + "
1145
- "LIME Interpretability"
1146
  ),
1147
 
1148
  description=(
1149
- "Enter an amino acid sequence "
1150
- "(≥10 amino acids) to obtain AMP "
1151
- "classification, confidence, LIME "
1152
- "interpretability, and MIC predictions."
1153
  ),
1154
 
1155
  examples=[
@@ -1163,9 +1014,11 @@ iface = gr.Interface(
1163
 
1164
 
1165
  # ============================================================
1166
- # LAUNCH
1167
  # ============================================================
1168
 
1169
  if __name__ == "__main__":
1170
 
 
 
1171
  iface.launch()
 
1
+ import os
2
  import joblib
3
  import numpy as np
4
  import pandas as pd
5
+ import gradio as gr
6
  import torch
7
 
 
 
 
8
  from math import expm1
9
 
10
+ from propy import AAComposition
11
+ from propy import Autocorrelation
12
+ from propy import CTD
13
+ from propy import PseudoAAC
14
 
15
+ from transformers import BertTokenizer, BertModel
 
 
 
 
 
16
 
17
+ from lime.lime_tabular import LimeTabularExplainer
 
18
 
19
 
20
  # ============================================================
21
+ # FILES
22
  # ============================================================
23
 
24
+ RF_MODEL_PATH = "RF.joblib"
25
+ SCALER_PATH = "norm (4).joblib"
26
+ TRAINING_DATA_PATH = "training_data.csv"
 
 
 
 
 
 
 
 
 
 
27
 
28
 
29
  # ============================================================
30
+ # MODEL SETTINGS
31
  # ============================================================
32
 
33
+ LIME_BACKGROUND_SIZE = 1000
34
+ LIME_NUM_SAMPLES = 3000
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
35
 
36
 
37
  # ============================================================
38
+ # EXACT 146 RF FEATURES
39
  # ============================================================
40
 
41
  selected_features = [
42
+ "_SolventAccessibilityC3", "_SecondaryStrC1",
43
+ "_SecondaryStrC3", "_ChargeC1", "_PolarityC1",
44
+ "_NormalizedVDWVC1", "_HydrophobicityC3",
45
+ "_SecondaryStrT23", "_PolarizabilityD1001",
46
+ "_PolarizabilityD2001", "_PolarizabilityD3001",
 
 
 
 
 
 
47
  "_SolventAccessibilityD1001",
48
  "_SolventAccessibilityD2001",
49
+ "_SolventAccessibilityD3001", "_SecondaryStrD1001",
50
+ "_SecondaryStrD1075", "_SecondaryStrD2001",
51
+ "_SecondaryStrD3001", "_ChargeD1001",
52
+ "_ChargeD1025", "_ChargeD2001", "_ChargeD3075",
53
+ "_ChargeD3100", "_PolarityD1001", "_PolarityD1050",
54
+ "_PolarityD2001", "_PolarityD3001",
55
+ "_NormalizedVDWVD1001", "_NormalizedVDWVD2001",
56
+ "_NormalizedVDWVD2025", "_NormalizedVDWVD2050",
57
+ "_NormalizedVDWVD3001", "_HydrophobicityD1001",
58
+ "_HydrophobicityD2001", "_HydrophobicityD3001",
 
 
 
 
 
 
 
 
 
 
 
 
59
  "_HydrophobicityD3025",
60
 
61
+ "A", "R", "D", "C", "E", "Q", "H", "I",
62
+ "M", "P", "Y", "V",
63
+
64
+ "AR", "AV", "RC", "RL", "RV", "CR", "CC",
65
+ "CL", "CK", "EE", "EI", "EL", "HC", "IA",
66
+ "IL", "IV", "LA", "LC", "LE", "LI", "LT",
67
+ "LV", "KC", "MA", "MS", "SC", "TC", "TV",
68
+ "YC", "VC", "VE", "VL", "VK", "VV",
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
69
 
70
  "MoreauBrotoAuto_FreeEnergy30",
71
 
 
124
  "GearyAuto_Mutability23",
125
  "GearyAuto_Mutability25",
126
  "GearyAuto_Mutability26",
 
127
  "GearyAuto_Mutability28",
128
  "GearyAuto_Mutability29",
129
  "GearyAuto_Mutability30",
130
 
131
+ "APAAC1", "APAAC4", "APAAC5", "APAAC6",
132
+ "APAAC8", "APAAC9", "APAAC12", "APAAC13",
133
+ "APAAC15", "APAAC18", "APAAC19", "APAAC24"
 
 
 
 
 
 
 
 
 
134
  ]
135
 
136
 
137
+ # ============================================================
138
+ # CHECK FEATURE COUNT
139
+ # ============================================================
140
+
141
+ print("=" * 70)
142
+ print("AMP + LIME APPLICATION")
143
+ print("=" * 70)
144
+
145
  print(
146
+ "Selected RF feature count:",
147
  len(selected_features)
148
  )
149
 
150
 
151
  # ============================================================
152
+ # LOAD RF
153
  # ============================================================
154
 
155
+ print("\nLoading Random Forest...")
 
 
 
 
 
 
 
 
 
 
156
 
157
+ model = joblib.load(
158
+ RF_MODEL_PATH
159
+ )
160
 
161
+ print("Random Forest loaded.")
 
 
162
 
163
  print(
164
+ "RF expects:",
165
+ model.n_features_in_,
166
+ "features"
 
 
167
  )
168
 
169
 
170
+ if model.n_features_in_ != len(
171
+ selected_features
172
+ ):
173
 
174
+ raise ValueError(
175
+ "RF feature count does not match "
176
+ "selected feature list."
177
  )
178
 
179
 
180
+ # ============================================================
181
+ # LOAD SCALER
182
+ # ============================================================
183
+
184
+ print("\nLoading scaler...")
185
+
186
+ scaler = joblib.load(
187
+ SCALER_PATH
188
+ )
189
+
190
+ print("Scaler loaded.")
191
+
192
  print(
193
+ "Scaler expects:",
194
+ scaler.n_features_in_,
195
+ "features"
196
  )
197
 
198
 
199
  # ============================================================
200
+ # BUILD REAL LIME DATA
201
  # ============================================================
202
 
203
+ def build_lime_background():
 
 
 
 
204
 
205
+ print("\n")
206
+ print("=" * 70)
207
+ print("LOADING ORIGINAL TRAINING DATA")
208
+ print("=" * 70)
209
 
 
210
 
211
+ if not os.path.exists(
212
+ TRAINING_DATA_PATH
213
+ ):
214
+
215
+ raise FileNotFoundError(
216
+ f"Missing {TRAINING_DATA_PATH}"
217
  )
218
+
219
+
220
+ df = pd.read_csv(
221
+ TRAINING_DATA_PATH
222
  )
223
 
224
 
225
+ if "Unnamed: 0" in df.columns:
 
 
226
 
227
+ df = df.drop(
228
+ "Unnamed: 0",
229
+ axis=1
230
+ )
 
 
231
 
232
 
233
+ print(
234
+ "Training data shape:",
235
+ df.shape
 
 
236
  )
 
237
 
238
 
239
+ # --------------------------------------------------------
240
+ # Remove non-feature columns
241
+ # --------------------------------------------------------
 
242
 
243
+ non_feature_columns = []
244
 
245
+ if "Sequence" in df.columns:
 
 
 
 
 
246
 
247
+ non_feature_columns.append(
248
+ "Sequence"
249
+ )
250
 
251
+ if "Activity" in df.columns:
 
 
 
252
 
253
+ non_feature_columns.append(
254
+ "Activity"
255
+ )
256
 
 
 
 
257
 
258
+ X_train = df.drop(
259
+ columns=non_feature_columns,
260
+ errors="ignore"
261
+ )
262
+
263
+
264
+ print(
265
+ "Feature matrix shape:",
266
+ X_train.shape
267
+ )
268
+
269
+
270
+ # --------------------------------------------------------
271
+ # Match scaler columns
272
+ # --------------------------------------------------------
273
+
274
+ if hasattr(
275
+ scaler,
276
+ "feature_names_in_"
277
+ ):
278
+
279
+ scaler_features = list(
280
+ scaler.feature_names_in_
281
+ )
282
+
283
+ missing = [
284
+ col
285
+ for col in scaler_features
286
+ if col not in X_train.columns
287
+ ]
288
+
289
+
290
+ if missing:
291
+
292
+ raise ValueError(
293
+ "Training data is missing "
294
+ "features required by scaler:\n\n"
295
+ + "\n".join(
296
+ missing[:30]
297
+ )
298
+ )
299
+
300
 
301
+ # EXACT scaler order
302
 
303
+ X_train = X_train[
304
+ scaler_features
305
+ ]
306
 
 
307
 
308
+ else:
309
+
310
+ if X_train.shape[1] != (
311
+ scaler.n_features_in_
312
+ ):
313
 
314
+ raise ValueError(
315
+ "Training feature count does "
316
+ "not match scaler."
317
+ )
318
+
319
+
320
+ # --------------------------------------------------------
321
+ # Select reproducible subset
322
+ # --------------------------------------------------------
323
+
324
+ if len(X_train) > LIME_BACKGROUND_SIZE:
325
+
326
+ X_train = (
327
+ X_train
328
+ .sample(
329
+ n=LIME_BACKGROUND_SIZE,
330
+ random_state=42
331
+ )
332
+ )
333
+
334
+
335
+ # --------------------------------------------------------
336
+ # SAME SCALING USED BY MODEL
337
+ # --------------------------------------------------------
338
+
339
+ print(
340
+ "Applying saved scaler..."
341
  )
342
 
343
+
344
+ X_scaled = scaler.transform(
345
+ X_train
346
+ )
347
+
348
+
349
+ X_scaled = pd.DataFrame(
350
+ X_scaled,
351
+ columns=X_train.columns
352
+ )
353
+
354
+
355
+ # --------------------------------------------------------
356
+ # EXACT 146 RF FEATURES
357
+ # --------------------------------------------------------
358
+
359
+ missing_selected = [
360
+ f
361
+ for f in selected_features
362
+ if f not in X_scaled.columns
363
+ ]
364
+
365
+
366
+ if missing_selected:
367
+
368
+ raise ValueError(
369
+ "The scaled training data is "
370
+ "missing RF features:\n\n"
371
+ + "\n".join(
372
+ missing_selected
373
+ )
374
+ )
375
+
376
+
377
+ lime_background = (
378
+ X_scaled[
379
+ selected_features
380
  ]
381
+ .fillna(0)
382
+ .astype(np.float64)
383
+ .values
384
  )
385
 
386
 
387
+ print(
388
+ "REAL LIME background shape:",
389
+ lime_background.shape
390
+ )
391
+
392
+
393
+ if lime_background.shape[1] != 146:
394
+
395
+ raise ValueError(
396
+ "LIME background does not contain "
397
+ "146 features."
398
+ )
399
+
400
+
401
+ return lime_background
402
 
403
 
404
  # ============================================================
405
  # CREATE LIME EXPLAINER
406
  # ============================================================
407
 
408
+ lime_training_data = (
409
+ build_lime_background()
410
+ )
411
+
412
+
413
  explainer = LimeTabularExplainer(
414
 
415
  training_data=lime_training_data,
 
429
  )
430
 
431
 
432
+ print("\nREAL LIME EXPLAINER READY.")
 
 
433
 
434
 
435
  # ============================================================
436
+ # PROTBERT
437
  # ============================================================
438
 
439
+ print("\nLoading ProtBert...")
 
 
 
440
 
441
  tokenizer = BertTokenizer.from_pretrained(
442
  "Rostlab/prot_bert",
443
  do_lower_case=False
444
  )
445
 
 
446
  protbert_model = BertModel.from_pretrained(
447
  "Rostlab/prot_bert"
448
  )
449
 
 
450
  device = torch.device(
451
  "cuda"
452
  if torch.cuda.is_available()
453
  else "cpu"
454
  )
455
 
 
456
  protbert_model = (
457
  protbert_model
458
  .to(device)
459
  .eval()
460
  )
461
 
 
462
  print(
463
+ "ProtBert device:",
464
  device
465
  )
466
 
467
 
468
  # ============================================================
469
+ # SEQUENCE CLEANING
470
  # ============================================================
471
 
472
  VALID_AA = set(
 
474
  )
475
 
476
 
 
 
 
 
477
  def clean_sequence(sequence):
478
 
479
  sequence = str(
480
  sequence
481
  ).upper().strip()
482
 
483
+ return "".join(
 
484
  aa
485
  for aa in sequence
486
  if aa in VALID_AA
487
  )
488
 
489
 
 
 
 
490
  # ============================================================
491
+ # FEATURE EXTRACTION FOR NEW INPUT
492
  # ============================================================
493
 
494
  def extract_features(sequence):
 
497
  sequence
498
  )
499
 
 
500
  if len(sequence) < 10:
501
 
502
+ return None, "Sequence too short."
 
 
 
503
 
504
 
505
  try:
506
 
507
+ dipeptide = (
508
+ AAComposition
509
+ .CalculateAADipeptideComposition(
 
 
 
510
  sequence
511
  )
512
  )
513
 
514
 
515
+ filtered_dipeptide = {
516
+ k: dipeptide[k]
 
 
517
  for k in list(
518
+ dipeptide.keys()
519
  )[:420]
520
  }
521
 
522
 
523
+ ctd = CTD.CalculateCTD(
524
+ sequence
 
 
 
 
 
 
525
  )
526
 
527
 
528
+ autocorr = (
529
+ Autocorrelation
530
+ .CalculateAutoTotal(
 
 
 
531
  sequence
532
  )
533
  )
534
 
535
 
536
+ pseaac = (
 
 
 
 
537
  PseudoAAC.GetAPseudoAAC(
538
  sequence,
539
  lamda=9
 
541
  )
542
 
543
 
 
 
 
 
544
  all_features = {}
545
 
546
+ all_features.update(ctd)
547
+ all_features.update(filtered_dipeptide)
548
+ all_features.update(autocorr)
549
+ all_features.update(pseaac)
 
 
 
 
 
 
 
 
 
 
 
550
 
551
 
552
  feature_df = pd.DataFrame(
 
554
  )
555
 
556
 
557
+ # Match scaler order
 
 
558
 
559
  if hasattr(
560
+ scaler,
561
  "feature_names_in_"
562
  ):
563
 
564
  scaler_features = list(
565
+ scaler.feature_names_in_
566
  )
567
 
568
 
569
+ missing = [
570
  f
571
  for f in scaler_features
572
  if f not in feature_df.columns
573
  ]
574
 
575
 
576
+ if missing:
577
 
578
  return (
579
  None,
580
+ "Missing scaler features: "
581
  + str(
582
+ missing[:20]
583
  )
584
  )
585
 
 
591
  )
592
 
593
 
594
+ # Same normalization
 
 
595
 
596
+ normalized = scaler.transform(
597
+ feature_df
 
 
598
  )
599
 
600
 
601
  normalized_df = pd.DataFrame(
602
+ normalized,
603
  columns=feature_df.columns
604
  )
605
 
606
 
607
+ # Same 146 selected features
 
 
608
 
609
+ missing = [
610
  f
611
  for f in selected_features
612
  if f not in normalized_df.columns
613
  ]
614
 
615
 
616
+ if missing:
617
 
618
  return (
619
  None,
620
  "Missing selected features: "
621
+ + str(missing[:20])
 
 
622
  )
623
 
624
 
625
+ X = (
626
  normalized_df[
627
  selected_features
628
  ]
629
  .fillna(0)
630
+ .astype(float)
 
 
 
 
631
  .values
 
632
  )
633
 
634
 
635
+ if X.shape != (1, 146):
 
 
 
636
 
637
  return (
638
  None,
639
+ f"Wrong final shape: {X.shape}"
 
640
  )
641
 
642
 
 
652
 
653
 
654
  # ============================================================
655
+ # MIC
656
  # ============================================================
657
 
658
  def predict_mic(sequence):
 
662
  )
663
 
664
 
 
 
 
 
 
 
 
 
 
 
 
665
  seq_spaced = " ".join(
666
  list(sequence)
667
  )
 
692
 
693
 
694
  embedding = (
695
+ outputs.last_hidden_state
 
696
  .mean(dim=1)
697
  .squeeze()
698
  .cpu()
 
701
  )
702
 
703
 
 
 
 
 
704
  bacteria_config = {
705
 
706
+ "E. coli": (
707
+ "coli_xgboost_model.pkl",
708
+ "coli_scaler.pkl",
709
+ None
710
+ ),
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
711
 
712
+ "S. aureus": (
713
+ "aur_xgboost_model.pkl",
714
+ "aur_scaler.pkl",
715
+ None
716
+ ),
717
 
718
+ "P. aeruginosa": (
719
+ "arg_xgboost_model.pkl",
720
+ "arg_scaler.pkl",
721
+ None
722
+ ),
723
 
724
+ "K. pneumoniae": (
725
+ "pne_mlp_model.pkl",
726
+ "pne_scaler.pkl",
727
+ "pne_pca.pkl"
728
+ )
729
  }
730
 
731
 
732
+ results = {}
733
 
734
 
735
+ for organism, config in (
736
  bacteria_config.items()
737
  ):
738
 
739
  try:
740
 
741
+ model_path, scaler_path, pca_path = (
742
+ config
743
  )
744
 
745
 
746
  mic_model = joblib.load(
747
+ model_path
748
+ )
749
+
750
+
751
+ mic_scaler = joblib.load(
752
+ scaler_path
753
  )
754
 
755
 
 
760
  )
761
 
762
 
763
+ if pca_path:
764
 
765
+ pca = joblib.load(
766
+ pca_path
767
  )
768
 
 
769
  transformed = (
770
+ pca.transform(
771
  scaled
772
  )
773
  )
 
777
  transformed = scaled
778
 
779
 
780
+ log_value = (
781
  mic_model.predict(
782
  transformed
783
  )[0]
784
  )
785
 
786
 
787
+ results[organism] = round(
788
+ expm1(
789
+ log_value
790
+ ),
791
  3
792
  )
793
 
794
 
 
 
 
 
 
795
  except Exception as e:
796
 
797
+ results[organism] = (
 
 
798
  f"Error: {str(e)}"
799
  )
800
 
801
 
802
+ return results
803
 
804
 
805
  # ============================================================
806
+ # MAIN PREDICTION
807
  # ============================================================
808
 
809
  def full_prediction(sequence):
810
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
811
  features, error = (
812
  extract_features(
813
  sequence
 
823
  )
824
 
825
 
826
+ # --------------------------------------------------------
827
+ # RF
828
+ # --------------------------------------------------------
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
829
 
830
+ prediction = model.predict(
831
+ features
832
+ )[0]
833
 
 
 
 
 
 
 
834
 
835
+ probabilities = (
836
+ model.predict_proba(
837
+ features
838
+ )[0]
839
+ )
840
 
 
 
 
841
 
842
+ class_index = list(
843
+ model.classes_
844
+ ).index(
845
+ prediction
846
+ )
847
 
 
 
 
848
 
849
+ confidence = round(
850
+ probabilities[class_index] * 100,
851
+ 2
852
+ )
853
 
 
 
 
854
 
855
+ label = (
856
+ "Antimicrobial Peptide (AMP)"
857
+ if prediction == 0
858
+ else
859
+ "Non-AMP"
860
+ )
861
 
862
 
863
  result = (
 
864
  "AMP CLASSIFICATION\n"
865
+ "===================\n\n"
866
+ f"Prediction: {label}\n"
867
  f"Confidence: {confidence}%\n"
868
  )
869
 
870
 
871
+ # --------------------------------------------------------
872
  # LIME
873
+ # --------------------------------------------------------
874
 
875
  result += (
876
+ "\n\nLIME EXPLANATION\n"
877
+ "===================\n"
 
878
  )
879
 
880
 
881
  try:
882
 
 
 
 
 
883
  explanation = (
884
  explainer.explain_instance(
885
 
886
+ features[0],
887
 
888
+ model.predict_proba,
889
 
890
  labels=[
891
  int(prediction)
 
893
 
894
  num_features=10,
895
 
896
+ num_samples=LIME_NUM_SAMPLES
 
 
 
 
 
 
 
897
  )
898
  )
899
 
900
 
901
  result += (
902
+ "\nTop features influencing "
903
+ "this prediction:\n"
904
  )
905
 
906
 
907
  for feature, weight in (
908
+ explanation.as_list(
909
+ label=int(prediction)
910
+ )
911
  ):
912
 
913
+ direction = (
914
+ "↑"
915
+ if weight > 0
916
+ else "↓"
917
+ )
 
 
918
 
919
 
920
  result += (
 
927
  except Exception as e:
928
 
929
  result += (
930
+ "\nLIME failed: "
931
+ f"{type(e).__name__}: "
932
+ f"{str(e)}"
933
  )
934
 
935
 
936
+ # --------------------------------------------------------
937
  # MIC
938
+ # --------------------------------------------------------
939
 
940
  if prediction == 0:
941
 
942
  result += (
943
+ "\n\n\nPREDICTED MIC\n"
944
+ "===================\n"
 
945
  )
946
 
947
 
 
964
 
965
  result += (
966
  "\n\nMIC prediction skipped "
967
+ "for Non-AMP."
968
  )
969
 
970
 
 
981
 
982
  inputs=gr.Textbox(
983
  label="Protein Sequence",
984
+ lines=5,
985
  placeholder=(
986
  "Enter an amino acid sequence..."
987
+ )
 
988
  ),
989
 
990
  outputs=gr.Textbox(
 
993
  ),
994
 
995
  title=(
996
+ "AMP & MIC Predictor "
997
+ "+ LIME Interpretability"
998
  ),
999
 
1000
  description=(
1001
+ "AMP classification, confidence, "
1002
+ "real LIME interpretation, and "
1003
+ "MIC prediction."
 
1004
  ),
1005
 
1006
  examples=[
 
1014
 
1015
 
1016
  # ============================================================
1017
+ # START
1018
  # ============================================================
1019
 
1020
  if __name__ == "__main__":
1021
 
1022
+ print("\nApplication starting...")
1023
+
1024
  iface.launch()