SeqOpt.run

SeqOpt.run(df_seq, df_feat, objectives, algorithm='nsga2', pop_size=50, n_gen=20, crossover='uniform', mutation='substitution', cx_prob=0.5, mut_prob=0.2, survival='mu_plus_lambda', variation='and', constraints=None, penalty='delta', hof_size=10, n_mut_max=5, region=None, to_aa=None, init='random', seed=None, jmd_n_len=10, jmd_c_len=10)[source]

Run multi-objective directed evolution and return the Pareto front of variants.

A population of multi-mutation variants of the single wild-type df_seq is evolved by NSGA-II (algorithm='nsga2') or an importance-ordered greedy walk (algorithm='greedy'), scored on every objective, and reduced to the non-dominated trade-off front.

Parameters:

  • df_seq (pd.DataFrame, shape (1, n_seq_info)) – DataFrame containing an entry column with unique protein identifiers, in the position-based format (sequence, tmd_start, tmd_stop). See SequenceFeature.get_df_parts() for the full df_seq format specification. Must hold exactly one wild-type sequence (SeqOpt optimizes one sequence per run).

  • df_feat (pd.DataFrame) – CPP feature set (output of CPP.run()) defining the features and the residue attribution (feat_importance / feat_impact, positions) the search reads.

  • objectives (list of (str, str, object)) – (name, goal, source) per objective; goal in {'max','min'} and source in {'delta_pred','delta_cpp','shift_score','n_mut'} or a callable(sequence) -> float. The callable receives the variant sequence and returns a scalar, so any external predictor (scikit / torch model, or a sequence-level tool / web API such as a topology or signal-peptide predictor) can be optimized as an objective; its result is cached per distinct variant. At least two objectives.

  • algorithm (str, default='nsga2') – 'nsga2' (population) or 'greedy' (importance-ordered single path).

  • pop_size (int, default=50) – Population size (NSGA-II only).

  • n_gen (int, default=20) – Number of generations (NSGA-II only).

  • crossover (str, default='uniform') – Crossover operator: 'uniform' / 'one_point' / 'two_point'.

  • mutation (str, default='substitution') – Mutation operator: 'substitution' or 'shift'.

  • cx_prob (float, default=0.5) – Per-pair crossover probability.

  • mut_prob (float, default=0.2) – Per-individual mutation probability.

  • survival (str, default='mu_plus_lambda') – Survival scheme: 'mu_plus_lambda' (elitist), 'mu_comma_lambda' or 'ea_simple' (generational replacement).

  • variation (str, default='and') – Variation scheme: 'and' (varAnd — crossover and mutation) or 'or' (varOr — each offspring is crossover or mutation or reproduction; needs cx_prob+mut_prob<=1).

  • constraints (list of callable, optional) – Feasibility predicates genome -> bool (True = feasible). Infeasible variants are penalized so the search avoids them.

  • penalty (str, default='delta') – Penalty applied to infeasible variants: 'delta' (fixed worst objective) or 'closest_valid' (penalty scaled by the number of violated constraints).

  • hof_size (int, default=10) – Size of the single-objective Hall of Fame (SeqOpt.hall_of_fame_) accumulated across generations.

  • n_mut_max (int, default=5) – Maximum number of point mutations per variant.

  • region (str or list of int, optional) – Restrict the mutable span (see SeqMut.scan()).

  • to_aa (list of str, optional) – Substitution alphabet (see SeqMut.scan()).

  • init (str, default='random') – Population initialization: 'random' or 'suggest' (warm start from the top single mutations).

  • seed (int, optional) – Per-call seed; overrides the constructor random_state.

  • jmd_n_len (int, default=10) – Length of JMD-N in number of amino acids.

  • jmd_c_len (int, default=10) – Length of JMD-C in number of amino acids.

Returns:

df_pareto – One row per final-population variant with entry, variant, n_mut, sequence_mut, one column per objective (named by objectives), the non-dominated rank (0 = best front) and the crowding distance, sorted by rank then descending crowding.

Return type:

pd.DataFrame

Examples

import numpy as np, pandas as pd
import matplotlib.pyplot as plt
from sklearn.ensemble import RandomForestClassifier
import aaanalysis as aa
aa.options["verbose"] = False

# Gamma-secretase (GSEC) substrate data + the bundled interpretable CPP feature set.
df_feat = aa.load_features(name="DOM_GSEC")           # 150 CPP features (with positions, feat_importance)
df_seq  = aa.load_dataset(name="DOM_GSEC", n=50)      # 100 TMD sequences, label 1 = GSEC substrate
labels  = df_seq["label"].to_list()

# A simple RandomForest substrate classifier on the CPP feature matrix.
sf = aa.SequenceFeature()
X = np.asarray(sf.feature_matrix(features=df_feat["feature"],
                                 df_parts=sf.get_df_parts(df_seq=df_seq),
                                 df_scales=aa.load_scales()), dtype=float)
model = RandomForestClassifier(n_estimators=100, random_state=0).fit(X, labels)

# Pick a NON-substrate as the wild-type and design a "super substrate": mutate its TMD to
# maximize the predicted substrate probability with as few mutations as possible.
wt = df_seq[df_seq["label"] == 0].iloc[[0]].reset_index(drop=True)
objectives = [("substrate", "max", "delta_pred"),     # raise P(GSEC substrate) (RF prediction shift)
              ("parsimony", "min", "n_mut")]          # with as few mutations as possible

Run NSGA-II and read off the Pareto front (rank=0 = the non-dominated trade-offs):

seqopt = aa.SeqOpt(mode="importance", model=model, target_class=1, random_state=42)
df_pareto = seqopt.run(df_seq=wt, df_feat=df_feat, objectives=objectives,
                       algorithm="nsga2", pop_size=40, n_gen=20, n_mut_max=5, region="tmd")
aa.display_df(df_pareto, n_rows=10, show_shape=True)

DataFrame shape: (7, 8)
  entry variant n_mut sequence_mut substrate parsimony rank crowding
1 Q14802 0 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 0.000000 0.000000 0 inf
2 Q14802 L37T+I55T+V56L+S58Q+A59R 5 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 37.000000 5.000000 0 inf
3 Q14802 L37G+A59R 2 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 29.000000 2.000000 0 0.429730
4 Q14802 A59R 1 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 16.000000 1.000000 0 0.316216
5 Q14802 L37G+V56L+A59R 3 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 33.000000 3.000000 0 0.281081
6 Q14802 A59K 1 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 16.000000 1.000000 0 0.275676
7 Q14802 L37G+V56L+S58Q+A59R 4 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 35.000000 4.000000 0 0.254054

Guidance mode. Above used mode=\"importance\" (static feat_importance). The headline mode=\"impact\" refits a ShapModel every generation (fuzzy labeling) to mutate the residues SHAP deems most important — pass the labeled reference set:

df_shap = aa.SeqOpt(mode="impact", model=model, target_class=1, df_seq_ref=df_seq,
                    labels=labels, random_state=0).run(
    df_seq=wt, df_feat=df_feat, objectives=objectives, pop_size=12, n_gen=4,
    n_mut_max=4, region="tmd")
aa.display_df(df_shap, n_rows=10, show_shape=True)

DataFrame shape: (9, 8)
  entry variant n_mut sequence_mut substrate parsimony rank crowding
1 Q14802 0 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 0.000000 0.000000 0 inf
2 Q14802 C49E+S58R 2 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 21.000000 2.000000 0 inf
3 Q14802 S58F 1 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 8.000000 1.000000 0 0.559524
4 Q14802 V56M 1 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 8.000000 1.000000 0 0.440476
5 Q14802 C49E 1 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 8.000000 1.000000 0 0.000000
6 Q14802 C49E+S58R+A59W 3 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 20.000000 3.000000 1 inf
7 Q14802 A50G 1 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 0.000000 1.000000 1 inf
8 Q14802 A50G+S58R 2 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 14.000000 2.000000 1 1.000000
9 Q14802 C49Y+A50G+S58R+A59W 4 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 20.000000 4.000000 2 inf

Algorithm + operators. Swap the search and the evolutionary operators: a deterministic importance-ordered greedy baseline, or NSGA-II with varOr variation, ea_simple survival, two-point crossover, shift mutation.

df_greedy = seqopt.run(df_seq=wt, df_feat=df_feat, objectives=objectives,
                       algorithm="greedy", n_mut_max=5, region="tmd")
df_or = seqopt.run(df_seq=wt, df_feat=df_feat, objectives=objectives, pop_size=40,
                   n_gen=20, n_mut_max=5, region="tmd", variation="or",
                   survival="ea_simple", crossover="two_point", mutation="shift",
                   cx_prob=0.5, mut_prob=0.3)
aa.display_df(df_greedy, n_rows=10, show_shape=True)

DataFrame shape: (5, 8)
  entry variant n_mut sequence_mut substrate parsimony rank crowding
1 Q14802 A59K 1 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 16.000000 1.000000 0 inf
2 Q14802 I55A+V56L+M57A+S58K+A59K 5 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 29.000000 5.000000 0 inf
3 Q14802 S58K+A59K 2 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 23.000000 2.000000 0 0.596154
4 Q14802 M57A+S58K+A59K 3 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 25.000000 3.000000 0 0.442308
5 Q14802 I55A+M57A+S58K+A59K 4 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 28.000000 4.000000 0 0.403846

Constraints + Hall of Fame. Forbid mutating the first TMD residue (a feasibility callable, penalized), and read the single-objective best-k variants:

tmd_start = int(wt["tmd_start"].iloc[0])
df_con = seqopt.run(df_seq=wt, df_feat=df_feat, objectives=objectives, pop_size=40,
                    n_gen=20, n_mut_max=5, region="tmd",
                    constraints=[lambda g: tmd_start not in g], penalty="delta",
                    hof_size=5)
print("Hall of Fame (best by substrate gain):", seqopt.hall_of_fame_)
aa.display_df(df_con, n_rows=10, show_shape=True)

Hall of Fame (best by substrate gain): ['C49L+S58R+A59R', 'C49L+I55C+S58R+A59R', 'I43V+C49L+I55W+A59R', 'I43V+C49L+A50G+I55W+A59R', 'I43V+C49L+A50S+I55W+A59R']
DataFrame shape: (6, 8)
  entry variant n_mut sequence_mut substrate parsimony rank crowding
1 Q14802 0 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 0.000000 0.000000 0 inf
2 Q14802 C49L+S58R+A59R 3 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 36.000000 3.000000 0 inf
3 Q14802 A59R 1 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 16.000000 1.000000 0 0.388889
4 Q14802 C49M+A59R 2 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 27.000000 2.000000 0 0.319444
5 Q14802 A59K 1 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 16.000000 1.000000 0 0.319444
6 Q14802 S58R+A59R 2 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 27.000000 2.000000 0 0.291667

Any predictor as an objective. Besides the features-model (delta_pred), an objective can be any callable(sequence) -> float — a scikit/torch model or a sequence-level tool / web API (e.g. a topology or signal-peptide predictor such as Phobius / TMHMM). Results are cached per variant. Here a simple aromatic-stretch proxy:

import re
def aromatic_stretch(sequence):
    return float(max((len(r) for r in re.findall(r"[FWY]+", sequence)), default=0))
obj_ext = [("substrate", "max", "delta_pred"),
           ("aromatic", "min", aromatic_stretch),  # external predictor / API goes here
           ("parsimony", "min", "n_mut")]
df_ext = seqopt.run(df_seq=wt, df_feat=df_feat, objectives=obj_ext, pop_size=30,
                    n_gen=15, n_mut_max=4, region="tmd")
aa.display_df(df_ext, n_rows=10, show_shape=True)

DataFrame shape: (6, 9)
  entry variant n_mut sequence_mut substrate aromatic parsimony rank crowding
1 Q14802 0 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 0.000000 3.000000 0.000000 0 inf
2 Q14802 C49F+G52L+S58R+A59K 4 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 36.000000 3.000000 4.000000 0 inf
3 Q14802 A59R 1 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 16.000000 3.000000 1.000000 0 inf
4 Q14802 C49F+A59K 2 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 26.000000 3.000000 2.000000 0 0.305556
5 Q14802 C49F+S58R+A59K 3 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 31.000000 3.000000 3.000000 0 0.259259
6 Q14802 A59K 1 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 16.000000 3.000000 1.000000 0 0.231481

Further parameters. init chooses the population start ('random' or 'suggest', a warm start from the top single mutations), to_aa restricts the substitution alphabet, seed makes a single run reproducible (overriding the constructor random_state), and jmd_n_len / jmd_c_len set the split geometry used to place mutations and build the CPP features. Kept small (few generations) so it finishes fast:

df_further = seqopt.run(df_seq=wt, df_feat=df_feat, objectives=objectives,
                        algorithm="nsga2", pop_size=12, n_gen=4, n_mut_max=4,
                        region="tmd",
                        init="suggest",                   # warm start from top single mutations
                        to_aa=["A", "L", "V", "K", "R"],  # restrict substitution alphabet
                        seed=1,                           # reproducible run (overrides random_state)
                        jmd_n_len=10, jmd_c_len=10)        # split geometry
aa.display_df(df_further, n_rows=10, show_shape=True)

DataFrame shape: (4, 8)
  entry variant n_mut sequence_mut substrate parsimony rank crowding
1 Q14802 S58R+A59K 2 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 23.000000 2.000000 0 inf
2 Q14802 0 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 0.000000 0.000000 0 inf
3 Q14802 A59R 1 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 16.000000 1.000000 0 0.597826
4 Q14802 A59K 1 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 16.000000 1.000000 0 0.402174