avGFP Fluorescence Optimisation with ALSEBO

This tutorial walks through a complete ALSEBO run end-to-end: from a FASTA file of VAE-generated sequences to a recommended set of high-fluorescence avGFP variants, using only ~65 simulated experiments.

---

Biological context

avGFP (Aequorea victoria Green Fluorescent Protein) is one of the most widely used reporter proteins in cell biology. Engineering brighter variants is a classic protein engineering benchmark.

We use the deep mutational scanning dataset from:

Sarkisyan et al. (2016). Local fitness landscape of the green fluorescent protein. Nature, 533, 397–401. https://doi.org/10.1038/nature17995

A pre-trained Support Vector Regression (SVR) model (best_svr_model.pkl) trained on this dataset acts as our fitness oracle — standing in for wet-lab fluorescence measurements.

---

Files in this tutorial

File	Description
generated_seqs.fasta	VAE-generated avGFP variant sequences (candidate pool)
avgfp_jhmmer.fasta	Multiple sequence alignment used for DCA featurisation
best_svr_model.pkl	Pre-trained SVR fitness oracle

Pipeline

generated_seqs.fasta  +  avgfp_jhmmer.fasta
           │
           ▼
  Step 1: Featurise sequence space  (DCA → seq_space.csv)
           │
           ▼
  Step 2: Explore fitness landscape  (t-SNE + SVR oracle)
           │
           ▼
  Step 3: Sample initial training set  (k-means diversity, n=15)
           │
           ▼
  Step 4: Simulate experiments  (SVR oracle)
           │
           ▼
  Step 5: Bayesian Optimisation loop  (GPR + UCB, 10 rounds × 5 sequences)
           │
           ▼
  Step 6: Analyse results  (convergence + trajectory)

Setup

import warnings
warnings.filterwarnings("ignore")

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import joblib
from sklearn.manifold import TSNE

from alsebo.seq_space import generate_seq_space
from alsebo.optimizer import (
    read_seq_files,
    gpr,
    seq_space_prediction,
    get_next_seq_bo,
    save_next_batch_results,
)
from alsebo.training_space import (
    sample_initial_training_sequnces,
    generate_sequence_training_file,
)

plt.rcParams.update({
    "font.size": 11,
    "axes.linewidth": 1.2,
    "axes.labelweight": "bold",
    "figure.dpi": 120,
})

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 11
      7 import matplotlib.pyplot as plt
      8 import joblib
      9 from sklearn.manifold import TSNE
     10 
---> 11 from alsebo.seq_space import generate_seq_space
     12 from alsebo.optimizer import (
     13     read_seq_files,
     14     gpr,

ModuleNotFoundError: No module named 'alsebo'

# ── Configuration ──────────────────────────────────────────────────────────
EXP_DIR     = "./"                   # directory containing all data files
MSA_FILE    = "avgfp_jhmmer.fasta"   # MSA for DCA featurisation
GEN_SEQS    = "generated_seqs.fasta" # VAE-generated candidate sequences
ORACLE_PATH = "best_svr_model.pkl"

INIT_BATCH  = 15    # sequences in the initial training set
BO_ROUNDS   = 10    # number of Bayesian optimisation rounds
BATCH_SIZE  = 5     # sequences recommended per round
BETA        = 2.0   # UCB exploration-exploitation trade-off

obj_config = {
    "names":      ["fitness"],
    "directions": ["max"],
}
# ───────────────────────────────────────────────────────────────────────────

# Remove leftovers from a previous run so we start clean
for fname in ["seq_exp_data.csv", "training_seqs.csv", "seq_space.csv"]:
    fpath = os.path.join(EXP_DIR, fname)
    if os.path.exists(fpath):
        os.remove(fpath)

oracle = joblib.load(ORACLE_PATH)
print("Oracle loaded. Ready to start.")

Step 1 — Featurise the Sequence Space

We first convert each candidate sequence in generated_seqs.fasta into a numerical feature vector using Direct Coupling Analysis (DCA).

DCA fits a maximum-entropy statistical model to the MSA and computes per-position features that capture both single-site preferences and residue–residue co-evolutionary couplings.

The result is written to seq_space.csv — one row per sequence.

Note: This step fits a DCA model on the MSA and is the slowest part of the pipeline (~1–3 min depending on MSA depth). It only needs to run once.

generate_seq_space(
    exp_dir=EXP_DIR,
    msa_fname=MSA_FILE,
    gen_seq_fasta_fname=GEN_SEQS,
    featuarization_method="DCA",
)

seq_df = pd.read_csv(f"{EXP_DIR}seq_space.csv")
print(f"Sequence space: {len(seq_df):,} sequences × {seq_df.shape[1]-1} DCA features")

Step 2 — Explore the Fitness Landscape

We use the SVR oracle to predict the fluorescence fitness of every sequence in the candidate pool, then project the high-dimensional DCA feature space down to 2-D with t-SNE to visualise the landscape.

features    = seq_df.drop("seq_id", axis=1)
all_seqs    = seq_df["seq_id"].tolist()
fitness_all = oracle.predict(features.values)

print("Running t-SNE projection...")
tsne_2d = TSNE(n_components=2, random_state=42).fit_transform(features)

fig, ax = plt.subplots(figsize=(7, 5))
sc = ax.scatter(tsne_2d[:, 0], tsne_2d[:, 1],
                c=fitness_all, cmap="magma", alpha=0.75, s=8)
plt.colorbar(sc, ax=ax, label="Predicted fluorescence fitness")
ax.set_xlabel("t-SNE 1")
ax.set_ylabel("t-SNE 2")
ax.set_title("avGFP Sequence Space — Predicted Fitness Landscape")
plt.tight_layout()
plt.show()

print(f"Fitness range : {fitness_all.min():.2f}  →  {fitness_all.max():.2f}")

Step 3 — Sample a Diverse Initial Training Set

Before running BO we need a small starting set of measured sequences. A random selection risks clustering in one region of the landscape.

ALSEBO avoids this with a t-SNE + k-means strategy:

1. Project feature space to 2-D

2. Partition into k = INIT_BATCH clusters

3. Pick the sequence closest to each centroid

This ensures the 15 initial sequences span the full landscape.

sample_initial_training_sequnces(EXP_DIR, training_seq_size=INIT_BATCH, manipold="TSNE")

training_df = pd.read_csv(f"{EXP_DIR}training_seqs.csv")
train_idx   = [all_seqs.index(s) for s in training_df["seq_id"] if s in all_seqs]

fig, ax = plt.subplots(figsize=(7, 5))
sc = ax.scatter(tsne_2d[:, 0], tsne_2d[:, 1],
                c=fitness_all, cmap="magma", alpha=0.45, s=8, label="All sequences")
ax.scatter(tsne_2d[train_idx, 0], tsne_2d[train_idx, 1],
           c="cyan", edgecolors="black", s=70, zorder=5,
           label=f"Initial batch (n={INIT_BATCH})")
plt.colorbar(sc, ax=ax, label="Predicted fluorescence fitness")
ax.set_xlabel("t-SNE 1")
ax.set_ylabel("t-SNE 2")
ax.set_title("Diverse Initial Training Set")
ax.legend(frameon=False)
plt.tight_layout()
plt.show()

Step 4 — Simulate Initial Experiments

In a real campaign you would now synthesise and measure the 15 selected sequences in the lab. Here the SVR oracle simulates those measurements instantly. The results are written to seq_exp_data.csv — the training log that ALSEBO appends to after every BO round.

train_features = training_df.drop("seq_id", axis=1)
init_fitness   = oracle.predict(train_features.values)

# generate_sequence_training_file expects list[list[float]]
obj_values = [[float(v)] for v in init_fitness]
generate_sequence_training_file(EXP_DIR, obj_config, obj_values)

print("Initial training data written to seq_exp_data.csv")
print(f"  Sequences measured : {INIT_BATCH}")
print(f"  Mean fitness       : {np.mean(init_fitness):.3f}")
print(f"  Best fitness       : {np.max(init_fitness):.3f}")

Step 5 — Bayesian Optimisation Loop

Each round ALSEBO:

1. Fits a Gaussian Process surrogate on all sequences measured so far

2. Predicts posterior mean and uncertainty across every untested candidate

3. Scores each candidate with UCB: score = mean + β × std

4. Selects the top-5 and evaluates them with the oracle

5. Appends results and repeats

After 10 rounds we will have evaluated 15 + 10 × 5 = 65 sequences — just ~3% of the candidate pool.

history = []  # stores (round, mean_fitness, max_fitness)

for round_idx in range(BO_ROUNDS):

    # Load current training data + remaining candidates
    x_train, y_train, x_space, seq_ids_remaining = read_seq_files(EXP_DIR, obj_config)

    # Fit GPR surrogate — one model per objective
    models = gpr(x_train, y_train)

    # Predict mean + uncertainty over the full candidate space
    preds = seq_space_prediction(models, x_space)

    # UCB acquisition → select top-k candidates
    next_seqs, scores, idx = get_next_seq_bo(
        seq_ids_remaining, obj_config, preds,
        top_k=BATCH_SIZE, strategy="UCB", beta=BETA,
    )

    # Oracle evaluation (simulates wet-lab measurement)
    new_fitness = oracle.predict(x_space.iloc[idx].values)

    # Append results to seq_exp_data.csv
    save_next_batch_results(
        EXP_DIR, next_seqs, idx, x_space, obj_config,
        obj_values=new_fitness.tolist(),
    )

    history.append({"round": round_idx + 1,
                    "mean":  float(np.mean(new_fitness)),
                    "max":   float(np.max(new_fitness))})

    best_so_far = max(h["max"] for h in history)
    print(f"Round {round_idx+1:2d}  |  batch max = {np.max(new_fitness):.3f}  "
          f"|  best so far = {best_so_far:.3f}")

print("\nOptimisation complete.")

Step 6 — Results

Convergence

Mean and best fitness per round, compared against the oracle maximum.

exp_df = pd.read_csv(f"{EXP_DIR}seq_exp_data.csv")

rounds_labels = ["Init"] + [f"R{h['round']}" for h in history]
mean_by_round = [exp_df["fitness"][:INIT_BATCH].mean()] + [h["mean"] for h in history]
max_by_round  = [exp_df["fitness"][:INIT_BATCH].max()]  + [h["max"]  for h in history]

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(rounds_labels, mean_by_round, "-o", label="Mean fitness", lw=2)
ax.plot(rounds_labels, max_by_round,  "-o", label="Best fitness",  lw=2)
ax.axhline(fitness_all.max(), color="grey", linestyle="--",
           label=f"Oracle max ({fitness_all.max():.2f})")
ax.set_xlabel("Optimisation round")
ax.set_ylabel("Predicted fluorescence fitness")
ax.set_title("ALSEBO Convergence on avGFP")
ax.legend(frameon=False)
plt.tight_layout()
plt.show()

Optimisation Trajectory

Best sequence per round overlaid on the fitness landscape.

exp_seqs = exp_df["seq_id"].tolist()

opt_tsne = []
for i, h in enumerate(history):
    start     = INIT_BATCH + i * BATCH_SIZE
    batch_fit = exp_df["fitness"].iloc[start : start + BATCH_SIZE].values
    best_seq  = exp_seqs[start + int(np.argmax(batch_fit))]
    if best_seq in all_seqs:
        opt_tsne.append(tsne_2d[all_seqs.index(best_seq)])

opt_tsne = np.array(opt_tsne)

fig, ax = plt.subplots(figsize=(7, 5))
sc = ax.scatter(tsne_2d[:, 0], tsne_2d[:, 1],
                c=fitness_all, cmap="magma", alpha=0.45, s=8)
ax.scatter(tsne_2d[train_idx, 0], tsne_2d[train_idx, 1],
           c="cyan", edgecolors="black", s=50, zorder=4, label="Initial training")

if len(opt_tsne) > 1:
    ax.plot(opt_tsne[:, 0], opt_tsne[:, 1], "w-", lw=2, zorder=5)
if len(opt_tsne) > 0:
    ax.scatter(opt_tsne[:, 0], opt_tsne[:, 1],
               c="white", edgecolors="black", s=50, zorder=6, label="Best per round")
    ax.scatter(opt_tsne[-1, 0], opt_tsne[-1, 1],
               c="red", edgecolors="black", s=90, zorder=7, label="Final best")

plt.colorbar(sc, ax=ax, label="Predicted fluorescence fitness")
ax.set_xlabel("t-SNE 1")
ax.set_ylabel("t-SNE 2")
ax.set_title("Optimisation Trajectory on the Fitness Landscape")
ax.legend(frameon=False, fontsize=9)
plt.tight_layout()
plt.show()

Summary

Initial training set	15 sequences (k-means diversity sampling)
BO rounds	10
Batch size per round	5
Total sequences evaluated	65 out of the full candidate pool (~3%)

ALSEBO navigated the avGFP fitness landscape using only a small fraction of the candidate pool, guided by the GPR surrogate and UCB acquisition function.

---

Things to try

• manipold="PCA" — swap t-SNE for PCA in the initial sampling step (faster, linear)

• beta=0.5 — more exploitative; converges faster but may miss the global optimum

• beta=5.0 — more exploratory; better for rugged or multi-modal landscapes

• ESM features — set featuarization_method="ESM" in generate_seq_space() for protein language model embeddings

• Multi-objective — add a second objective (e.g. thermostability) to obj_config and provide two values per sequence in obj_values