StackAPP Recreated Code
Color-coded + heavily commented Python implementation
Legend: Blue = Python keywords, Green = strings, Gray = comments. This PDF contains a full
runnable script that recreates the paper's pipeline (AAC + APAAC feature fusion + stacking
ensemble).
Run it like this: python [Link] --data_csv your_data.csv --lam 10 --w 0.5
Generated PDF of the StackAPP recreation code
StackAPP Recreated Code (continued) Page 1
1 #!/usr/bin/env python3
2 """
3 StackAPP recreation (from the paper "StackAPP: Advancing autophagy protein identificati
on with ensemble learning")
4
5 What this program does:
6 1) Reads protein sequences and labels (APP = 1, non-APP = 0)
7 2) Builds features:
8 - AAC : counts of amino acids (20 numbers)
9 - APAAC: AAC + extra "order" features (how amino acids relate across the sequence)
10 3) Fuses features: [AAC + APAAC]
11 4) Trains a stacking model:
12 Base learners : RandomForest, XGBoost, LightGBM, GradientBoosting
13 Meta learner : LogisticRegression
14 5) Evaluates with:
15 - 10-fold cross validation (on training set)
16 - independent test (held-out set)
17
18 Tip:
19 If you have the official train/test files from the paper, add a 'split' column in the C
SV
20 with values 'train' or 'test' and run with --use_official_split.
21 """
22
23 from __future__ import annotations
24
25 import argparse
26 from typing import Dict, List, Tuple
27
28 import numpy as np
29 import pandas as pd
30
31 from [Link] import RandomForestClassifier, GradientBoostingClassifier, Stacki
ngClassifier
32 from sklearn.linear_model import LogisticRegression
33 from sklearn.model_selection import StratifiedKFold, train_test_split
34 from [Link] import (
35 accuracy_score, precision_score, f1_score, matthews_corrcoef,
36 confusion_matrix, cohen_kappa_score
37 )
38
39 # Optional dependencies (install if missing):
40 # pip install xgboost lightgbm
41 try:
42 from xgboost import XGBClassifier
43 except Exception:
44 XGBClassifier = None
45
46 try:
47 from lightgbm import LGBMClassifier
48 except Exception:
49 LGBMClassifier = None
50
51
52 # ---------------------------------------------------------------------------
53 # 1) Amino acid basics
54 # ---------------------------------------------------------------------------
55
StackAPP Recreated Code (continued) Page 2
56 # The 20 standard amino acids (we will ignore any others like X, B, Z, U, etc.)
57 AA_ORDER = list("ACDEFGHIKLMNPQRSTVWY")
58 AA_TO_IDX = {aa: i for i, aa in enumerate(AA_ORDER)}
59
60
61 # ---------------------------------------------------------------------------
62 # 2) Two properties needed for APAAC:
63 # - hydrophobicity (likes oil)
64 # - hydrophilicity (likes water)
65 #
66 # We will normalize these values so they are centered around 0.
67 # ---------------------------------------------------------------------------
68
69 RAW_HYDROPHOBICITY = {
70 "A": 0.62, "C": 0.29, "D": -0.90, "E": -0.74, "F": 1.19,
71 "G": 0.48, "H": -0.40, "I": 1.38, "K": -1.50, "L": 1.06,
72 "M": 0.64, "N": -0.78, "P": 0.12, "Q": -0.85, "R": -2.53,
73 "S": -0.18, "T": -0.05, "V": 1.08, "W": 0.81, "Y": 0.26
74 }
75
76 RAW_HYDROPHILICITY = {
77 "A": -0.50, "C": -1.00, "D": 3.00, "E": 3.00, "F": -2.50,
78 "G": 0.00, "H": -0.50, "I": -1.80, "K": 3.00, "L": -1.80,
79 "M": -1.30, "N": 0.20, "P": 0.00, "Q": 0.20, "R": 3.00,
80 "S": 0.30, "T": -0.40, "V": -1.50, "W": -3.40, "Y": -2.30
81 }
82
83
84 def normalize_property(prop: Dict[str, float]) -> Dict[str, float]:
85 """
86 Make the property values have:
87 - mean (average) = 0
88 - std (spread) = 1
89 This is a common trick so numbers are comparable.
90 """
91 vals = [Link]([prop[aa] for aa in AA_ORDER], dtype=float)
92 mean = float([Link]())
93 std = float([Link](ddof=0))
94 if std == 0:
95 raise ValueError("Property has zero standard deviation; cannot normalize.")
96 return {aa: (prop[aa] - mean) / std for aa in AA_ORDER}
97
98
99 HYDROPHOBICITY = normalize_property(RAW_HYDROPHOBICITY)
100 HYDROPHILICITY = normalize_property(RAW_HYDROPHILICITY)
101
102
103 # ---------------------------------------------------------------------------
104 # 3) Feature functions (AAC and APAAC)
105 # ---------------------------------------------------------------------------
106
107 def clean_sequence(seq: str) -> str:
108 """
109 Keep only the 20 standard amino acids.
110 Anything else is removed.
111 """
112 seq = [Link]().upper()
113 return "".join([c for c in seq if c in AA_TO_IDX])
StackAPP Recreated Code (continued) Page 3
114
115
116 def aac(seq: str) -> [Link]:
117 """
118 AAC = Amino Acid Composition
119 Output: 20 numbers (frequencies that sum to 1).
120 """
121 seq = clean_sequence(seq)
122 L = len(seq)
123 if L == 0:
124 return [Link](20, dtype=float)
125
126 counts = [Link](20, dtype=float)
127 for c in seq:
128 counts[AA_TO_IDX[c]] += 1.0
129 return counts / L
130
131
132 def apaac(seq: str, lam: int = 10, w: float = 0.5) -> [Link]:
133 """
134 APAAC = Amphiphilic Pseudo-Amino Acid Composition
135
136 It has:
137 - 20 AAC-like parts
138 - 2*lam extra "order" parts (lam for hydrophobicity + lam for hydrophilicity)
139
140 Big idea:
141 Look at pairs of amino acids that are k steps apart and compute correlations.
142 """
143 seq = clean_sequence(seq)
144 L = len(seq)
145 if L == 0:
146 return [Link](20 + 2 * lam, dtype=float)
147
148 # Start from AAC (these sum to 1)
149 freq = aac(seq)
150
151 # Compute correlation factors (taus)
152 taus: List[float] = []
153 lam_eff = min(lam, max(1, L - 1)) # cannot look farther than sequence length
154
155 for k in range(1, lam_eff + 1):
156 s_hydro = 0.0
157 s_hphil = 0.0
158 for i in range(L - k):
159 a1 = seq[i]
160 a2 = seq[i + k]
161 s_hydro += HYDROPHOBICITY[a1] * HYDROPHOBICITY[a2]
162 s_hphil += HYDROPHILICITY[a1] * HYDROPHILICITY[a2]
163
164 [Link](s_hydro / (L - k))
165 [Link](s_hphil / (L - k))
166
167 # If sequence is short, we pad remaining taus with zeros
168 if lam_eff < lam:
169 [Link]([0.0] * (2 * (lam - lam_eff)))
170
171 taus_arr = [Link](taus, dtype=float)
StackAPP Recreated Code (continued) Page 4
172
173 # Combine AAC and taus using a weight w
174 denom = 1.0 + w * float(taus_arr.sum())
175 if denom <= 1e-12:
176 denom = 1e-12
177
178 p_aa = freq / denom
179 p_tau = (w * taus_arr) / denom
180 return [Link]([p_aa, p_tau], axis=0)
181
182
183 def fuse_features(seqs: List[str], lam: int = 10, w: float = 0.5) -> [Link]:
184 """
185 Feature fusion: AAC + APAAC (concatenate).
186 """
187 feats = []
188 for s in seqs:
189 [Link]([Link]([aac(s), apaac(s, lam=lam, w=w)], axis=0))
190 return [Link](feats)
191
192
193 # ---------------------------------------------------------------------------
194 # 4) Loading data
195 # ---------------------------------------------------------------------------
196
197 def load_csv(path: str) -> Tuple[List[str], [Link], [Link]]:
198 """
199 CSV must have:
200 - sequence : protein sequence (string)
201 - label : 1 for APP, 0 for non-APP
202 Optional:
203 - split : 'train' or 'test'
204 """
205 df = pd.read_csv(path)
206 if "sequence" not in [Link] or "label" not in [Link]:
207 raise ValueError("CSV must contain columns: sequence,label")
208 seqs = df["sequence"].astype(str).tolist()
209 y = df["label"].astype(int).to_numpy()
210 return seqs, y, df
211
212
213 # ---------------------------------------------------------------------------
214 # 5) Metrics helper (same set used in the paper)
215 # ---------------------------------------------------------------------------
216
217 def compute_metrics(y_true: [Link], y_pred: [Link]) -> Dict[str, float]:
218 """
219 Return a dictionary of metrics:
220 accuracy, sensitivity, specificity, precision, f1, mcc, kappa
221 """
222 tn, fp, fn, tp = confusion_matrix(y_true, y_pred, labels=[0, 1]).ravel()
223
224 acc = accuracy_score(y_true, y_pred)
225 sens = tp / (tp + fn) if (tp + fn) else 0.0 # also called recall for posit
ive class
226 spec = tn / (tn + fp) if (tn + fp) else 0.0
227 prec = precision_score(y_true, y_pred, zero_division=0)
228 f1 = f1_score(y_true, y_pred, zero_division=0)
StackAPP Recreated Code (continued) Page 5
229 mcc = matthews_corrcoef(y_true, y_pred)
230 kappa = cohen_kappa_score(y_true, y_pred)
231
232 return {
233 "accuracy": acc,
234 "sensitivity": sens,
235 "specificity": spec,
236 "precision": prec,
237 "f1": f1,
238 "mcc": mcc,
239 "kappa": kappa,
240 "tn": float(tn),
241 "fp": float(fp),
242 "fn": float(fn),
243 "tp": float(tp),
244 }
245
246
247 # ---------------------------------------------------------------------------
248 # 6) Main training + evaluation
249 # ---------------------------------------------------------------------------
250
251 def build_model(seed: int = 42):
252 """
253 Build the stacking model:
254 Base learners: RF, XGB, LGBM, GBC
255 Meta learner : Logistic Regression
256
257 Stacking uses predict_proba, so the meta learner sees probability scores.
258 """
259 rf = RandomForestClassifier(n_estimators=500, random_state=seed, n_jobs=-1)
260
261 if XGBClassifier is None:
262 raise ImportError("xgboost is not installed. Run: pip install xgboost")
263 xgb = XGBClassifier(
264 n_estimators=500,
265 max_depth=5,
266 learning_rate=0.05,
267 subsample=0.9,
268 colsample_bytree=0.9,
269 reg_lambda=1.0,
270 random_state=seed,
271 eval_metric="logloss",
272 n_jobs=-1
273 )
274
275 if LGBMClassifier is None:
276 raise ImportError("lightgbm is not installed. Run: pip install lightgbm")
277 lgbm = LGBMClassifier(
278 n_estimators=800,
279 learning_rate=0.03,
280 num_leaves=31,
281 subsample=0.9,
282 colsample_bytree=0.9,
283 random_state=seed
284 )
285
286 gbc = GradientBoostingClassifier(random_state=seed)
StackAPP Recreated Code (continued) Page 6
287
288 meta = LogisticRegression(max_iter=5000, solver="lbfgs")
289
290 model = StackingClassifier(
291 estimators=[("rf", rf), ("xgb", xgb), ("lgbm", lgbm), ("gbc", gbc)],
292 final_estimator=meta,
293 stack_method="predict_proba",
294 passthrough=False,
295 cv=5,
296 n_jobs=-1
297 )
298 return model
299
300
301 def main():
302 ap = [Link]()
303 ap.add_argument("--data_csv", type=str, required=True,
304 help="CSV with columns: sequence,label and optional split=train/tes
t")
305 ap.add_argument("--use_official_split", action="store_true",
306 help="Use split column for train/test instead of random split")
307 ap.add_argument("--split_col", type=str, default="split")
308 ap.add_argument("--lam", type=int, default=10, help="APAAC lambda (max lag).")
309 ap.add_argument("--w", type=float, default=0.5, help="APAAC weight w (paper uses 0.
5).")
310 ap.add_argument("--seed", type=int, default=42)
311 args = ap.parse_args()
312
313 _, _, df = load_csv(args.data_csv)
314
315 # Choose train/test
316 if args.use_official_split:
317 if args.split_col not in [Link]:
318 raise ValueError(f"Missing split column '{args.split_col}' in CSV.")
319 train_df = df[df[args.split_col].astype(str).[Link]() == "train"].copy()
320 test_df = df[df[args.split_col].astype(str).[Link]() == "test"].copy()
321 if len(train_df) == 0 or len(test_df) == 0:
322 raise ValueError("Split column must contain 'train' and 'test'.")
323 else:
324 # A practical stratified split (keeps label balance).
325 train_df, test_df = train_test_split(
326 df,
327 test_size=0.284,
328 random_state=[Link],
329 stratify=df["label"].astype(int)
330 )
331
332 train_seqs = train_df["sequence"].astype(str).tolist()
333 y_train = train_df["label"].astype(int).to_numpy()
334
335 test_seqs = test_df["sequence"].astype(str).tolist()
336 y_test = test_df["label"].astype(int).to_numpy()
337
338 # Build features
339 X_train = fuse_features(train_seqs, lam=[Link], w=args.w)
340 X_test = fuse_features(test_seqs, lam=[Link], w=args.w)
341
342 model = build_model(seed=[Link])
StackAPP Recreated Code (continued) Page 7
343
344 # --------------------- 10-fold CV on training set ---------------------
345 skf = StratifiedKFold(n_splits=10, shuffle=True, random_state=[Link])
346 fold_metrics = []
347
348 for fold, (tr_idx, va_idx) in enumerate([Link](X_train, y_train), start=1):
349 X_tr, X_va = X_train[tr_idx], X_train[va_idx]
350 y_tr, y_va = y_train[tr_idx], y_train[va_idx]
351
352 [Link](X_tr, y_tr)
353 pred = [Link](X_va)
354
355 m = compute_metrics(y_va, pred)
356 m["fold"] = fold
357 fold_metrics.append(m)
358
359 cv_df = [Link](fold_metrics)
360 print("\n=== 10-fold cross validation (training set) ===")
361 print(cv_df[["accuracy","sensitivity","specificity","precision","f1","mcc","kappa"]].d
escribe().loc[["mean","std"]].round(4))
362
363 # --------------------- Independent test ---------------------
364 [Link](X_train, y_train)
365 test_pred = [Link](X_test)
366 test_m = compute_metrics(y_test, test_pred)
367
368 print("\n=== Independent test ===")
369 print([Link]([test_m])[["accuracy","sensitivity","specificity","precision","f
1","mcc","kappa"]].round(4).to_string(index=False))
370
371 print("\nConfusion matrix (rows=true, cols=pred), labels=[0,1]:")
372 print([Link]([[int(test_m["tn"]), int(test_m["fp"])],
373 [int(test_m["fn"]), int(test_m["tp"])]]))
374
375
376 if __name__ == "__main__":
377 main()