P5: Engineer features

This protocol engineers the CPP feature space: parts, splits, and numerical channels (embeddings or structure). Feature engineering is the stage between raw sequences and a CPP signature. Before Comparative Physicochemical Profiling (CPP) can contrast a test group (label=1) against a reference group (label=0) and read out what physicochemically distinguishes them, every protein has to be turned into a structured feature space. AAanalysis builds that space from one small grammar: a feature is always a Part-Split-Scale triple, a named region (part), a rule for which positions to average (split), and a value axis. Here we engineer that space at the domain level (dataset prefix DOM_), where the unit of comparison is the TMD-centric part set.

A feature = part × split × scale: where in the sequence (tmd, jmd_n_tmd_n, tmd_c_jmd_c), how the positions are picked (Segment / Pattern / PeriodicPattern), and what value is averaged. The value axis has two interchangeable sources: a physicochemical scale (sequence mode, run()) or a per-residue numerical tensor, a protein-language-model (PLM) embedding or a structure channel (numerical mode, run_num()). Same grammar, same df_feat output; only the value source changes. That is exactly what lets learned and structural representations ride the same interpretable CPP machinery.

When to use it. Reach for this protocol once you have labelled sequences (a df_seq with a binary label) and need to turn them into a CPP-ready feature space, the step that sits underneath every CPP run, whether the downstream task is determinant discovery or prediction. It answers the three engineering questions hiding under that grammar:

  • Where? Slice each protein into TMD-centric parts. By default get_df_parts returns the base part tmd plus the composites jmd_n_tmd_n and tmd_c_jmd_c (built from the base parts jmd_n / tmd / jmd_c).

  • Which positions? Pick a split type: Segment (a contiguous block), Pattern, or PeriodicPattern (spaced positions).

  • What value? Average a physicochemical scale (sequence mode, consumed by run()) or a per-residue numerical vector such as a PLM embedding or a structural channel (numerical mode, consumed by run_num()).

Both modes share the identical Part-Split grammar and emit the identical df_feat schema, so anything you learn in sequence mode transfers directly to embeddings and structure.

When not to use it. Skip the manual enumeration if you just want a signature with sensible defaults: run() builds parts and splits internally, so you only need this protocol when you want to control the feature space or feed in numerical channels. It is also not where you decide the compositional vs positional trade-off (that is P6: Compositional vs positional), and for residue-level tasks you first manufacture fixed-length windows (parts named by the Schechter-Berger convention) in Protocol 3 - Construct sets & sampling before engineering features.

Input. A df_seq with one row per protein: an entry identifier, a sequence, and a binary label (test group = 1 vs. reference group = 0). For a domain-level task it also carries tmd_start / tmd_stop, from which the part boundaries are derived. Upstream this df_seq typically arrives from Protocol 3 - Construct sets & sampling (residue level) or from a curated labelled set (domain level, as here).

The numerical mode adds one more input: a dict_num mapping each entry to a per-residue tensor of shape (L, D), where L equals the sequence length and D is the number of numerical channels (embedding or structure dimensions). Real PLM embeddings (ProtT5, ESM) are heavy to compute, so you generate them separately on a GPU (e.g. the companion Colab notebook), then save the (L, D) tensors and load them here; AAanalysis never runs the model itself. To keep this protocol offline and fast we instead build a small synthetic dict_num with NumPy. The real workflow would substitute a ProtT5 / ESM embedding or a DSSP / structure tensor of the same (L, D) shape. We use the bundled DOM_GSEC gamma-secretase dataset (substrates vs. non-substrates) throughout.

import aaanalysis as aa
import numpy as np

aa.options["verbose"] = False
aa.options["random_state"] = 42

# Domain-level gamma-secretase data: substrates (label=1) vs non-substrates (0).
# load_dataset(name=..., n=10) returns 2*n = 20 rows (10 per class).
df_seq = aa.load_dataset(name="DOM_GSEC", n=10,
                         random=False,                # take the first n per class, so the demo is stable
                         non_canonical_aa="remove",   # drop entries carrying non-canonical residues
                         min_len=None, max_len=None,  # no length window
                         aa_window_size=9,            # only used by residue-level (AA_*) datasets
                         verbose=False)
labels = df_seq["label"].to_list()          # use the label column, NOT len(df_seq)
aa.display_df(df=df_seq, n_rows=10, show_shape=True)

DataFrame shape: (20, 9)
  entry gene sequence label tmd_start tmd_stop jmd_n tmd jmd_c
1 Q14802 FXYD3 MQKVTLGLLVFLAGF...PGETPPLITPGSAQS 0 37 59 NSPFYYDWHS LQVGGLICAGVLCAMGIIIVMSA KCKCKFGQKS
2 Q86UE4 MTDH MAARSWQDELAQQAE...SPKQIKKKKKARRET 0 50 72 LGLEPKRYPG WVILVGTGALGLLLLFLLGYGWA AACAGARKKR
3 Q969W9 PMEPA1 MHRLMGVNSTAAAAA...AIWSKEKDKQKGHPL 0 41 63 FQSMEITELE FVQIIIIVVVMMVMVVVITCLLS HYKLSARSFI
4 P53801 PTTG1IP MAPGVARGPTPYWRL...GLFKEENPYARFENN 0 97 119 RWGVCWVNFE ALIITMSVVGGTLLLGIAICCCC CCRRKRSRKP
5 Q8IUW5 RELL1 MAPRALPGSAVLAAA...EVPATPVKRERSGTE 0 59 81 NDTGNGHPEY IAYALVPVFFIMGLFGVLICHLL KKKGYRCTTE
6 P01135 TGFA MVPSAGQLALFALGI...LLKGRTACCHSETVV 0 99 121 AVVAASQKKQ AITALVVVSIVALAVLIITCVLI HCCQVRKHCE
7 O43914 TYROBP MGGLEPCSRLLLLPL...SDVYSDLNTQRPYYK 0 42 64 DCSCSTVSPG VLAGIVMGDLVLTVLIALAVYFL GRLVPRGRGA
8 P05556 ITGB1 MNLQPIFWIGLISSV...KSAVTTVVNPKYEGK 0 729 751 ENPECPTGPD IIPIVAGVVAGIVLIGLALLLIW KLLMIIHDRR
9 P16234 PDGFRA MGTSHPAFLVLGCLL...DIGIDSSDLVEDSFL 0 527 549 VAPTLRSELT VAAAVLVLLVIVIISLIVLVVIW KQKPRYEIRW
10 P50895 BCAM MEPPDAPAQARGAPR...SGGARGGSGGFGDEC 0 549 571 TVSPQTSQAG VAVMAVAVSVGLLLLVVAVFYCV RRKGGPCCRQ

Run: sequence mode, scale-based features. Sequence mode turns df_seq into a CPP-ready feature space in three small steps: slice the parts, configure how to split, then enumerate the full grammar (see the SequenceFeature tutorial, tutorial3b_sequence_feature, for the function details):

  1. get_df_parts() slices each protein into the TMD-centric parts (by default the base part tmd plus the composites jmd_n_tmd_n and tmd_c_jmd_c), producing the df_parts that CPP consumes.

  2. get_split_kws() builds the split configuration, one entry per split type (Segment, Pattern, PeriodicPattern).

  3. get_features() enumerates the full Part-Split-Scale space as a list of feature ids.

sf = aa.SequenceFeature()

# 1) Parts: slice sequences into TMD-centric parts
df_parts = sf.get_df_parts(df_seq=df_seq,
                           list_parts=None,       # None -> tmd, jmd_n_tmd_n, tmd_c_jmd_c
                           all_parts=False,       # True -> all 9 parts
                           jmd_n_len=10, jmd_c_len=10,
                           tmd_len=None,          # only for sequence-based df_seq (no tmd_start/tmd_stop)
                           remove_entries_with_gaps=False,
                           replace_non_canonical_aa=False)

# 2) Split configuration: Segment / Pattern / PeriodicPattern
split_kws = sf.get_split_kws(split_types=["Segment", "Pattern", "PeriodicPattern"],
                             n_split_min=1, n_split_max=15,                     # Segment granularity
                             steps_pattern=None, n_min=2, n_max=4, len_max=15,  # Pattern grid
                             steps_periodicpattern=None,                        # PeriodicPattern steps
                             strategy=None)                 # "compositional" / "positional" shortcut

# Split types in the configuration: ['Segment', 'Pattern', 'PeriodicPattern']
list(split_kws)

['Segment', 'Pattern', 'PeriodicPattern']

# 3) Enumerate the PART-SPLIT-SCALE space. The full scale set explodes
#    combinatorially, so restrict to a 3-scale slice for the demo.
df_scales = aa.load_scales()
small_scales = list(df_scales.columns[:3])

features = sf.get_features(list_parts=None, all_parts=False,
                           split_kws=split_kws, list_scales=small_scales)
# Number enumerated (3 parts x splits x 3 scales): len(features)
len(features), features[:5]

(2970,
 ['TMD-Segment(1,1)-ANDN920101',
  'TMD-Segment(1,1)-ARGP820101',
  'TMD-Segment(1,1)-ARGP820102',
  'TMD-Segment(1,2)-ANDN920101',
  'TMD-Segment(1,2)-ARGP820101'])

# Optional: reduce redundancy. Build the feature matrix, then drop columns
# that correlate above max_cor with an already-kept feature.
# Note: filter_correlation keeps the FIRST-SEEN representative of each
# correlated group. Applied pre-CPP on the raw PART-SPLIT-SCALE enumeration
# order, it yields a redundancy-reduced but importance-agnostic set. To keep
# the most important representative, apply it to an importance-ordered
# CPP.run / df_feat output instead (see Protocol 1 - CPP signature).
X = sf.feature_matrix(features=features, df_parts=df_parts,
                      df_scales=None,        # None -> the default scale set
                      accept_gaps=False,
                      n_jobs=1,              # serial: avoids the spawn footgun on Py3.14 + macOS
                      batch=False,           # True -> chunk the columns to cap memory
                      df_seq=None, df_parts_kws=None)   # alternative: build df_parts on the fly
is_selected = aa.NumericalFeature.filter_correlation(X=X, max_cor=0.7)

# Features kept after redundancy filtering:
int(is_selected.sum())

368

Run: numerical mode, embedding or structure features. Numerical mode replaces the per-AA scale lookup with a per-residue tensor. The Part-Split grammar is unchanged; only the value source differs.

The numerical-mode tensor below is synthetic NumPy noise, so every ``abs_auc`` / ``mean_dif`` it produces is pure chance: treat it as a plumbing demo, not a signature. Swap in a real ProtT5 / ESM embedding or a DSSP / structure tensor of the same (L, D) shape for a meaningful readout.

  1. Build a dict_num of {entry: (L, D) ndarray}. We synthesise one with np.random.default_rng(42); rng.random already yields values in [0, 1), which is exactly the range run_num() expects (its max_std_test pre-filter is calibrated for [0, 1]), so no normalisation step is needed. With real PLM embeddings the values are unbounded, so you would first run them through encode() (embp.encode(df_seq=df_seq, embeddings=...), the sanctioned step from raw embeddings to dict_num that fits one [0, 1] normaliser per channel) before get_parts.

  2. :meth:`~aaanalysis.EmbeddingPreprocessor.build_scales` / ``build_cat`` name the D dimensions (df_scales of shape (20, D), df_cat of shape (D, ...)) so the CPP invariant D == len(df_scales.columns) == len(df_cat) holds. The per-AA scale values are not the per-residue value source for run_num (that is dict_num_parts), but the df_cat categories, derived from these pseudo-scales by clustering, do drive category-aware redundancy filtering and grouping. Deriving df_scales/df_cat from the same dict_num (rather than fabricating an arbitrary (20, D) / df_cat) is therefore what keeps both the invariant and the categories meaningful. build_scales(return_std=True) hands back the per-channel spread as well, which build_cat weights when it clusters the D channels.

  3. :meth:`~aaanalysis.NumericalFeature.get_parts` slices the sequence strings AND the tensor with the same boundaries, returning (df_parts, dict_num_parts).

  4. :meth:`~aaanalysis.CPP.run_num` runs the identical CPP algorithm on the sliced tensor.

D = 4                                   # tiny embedding dimensionality
rng = np.random.default_rng(42)

# Per-residue tensor: one (L, D) array per entry, L == sequence length, in [0, 1)
dict_num = {e: rng.random((len(s), D))
            for e, s in zip(df_seq["entry"], df_seq["sequence"])}

# Name the D dimensions. dict_num_parts (not these per-AA values) is the
# per-residue value source for run_num; but df_cat's categories - derived
# from these pseudo-scales - drive category-aware redundancy filtering, so we
# derive both from the same dict_num rather than fabricating them.
# Shapes: df_scales_num is (20, D); df_cat_num is (D, ...).
embp = aa.EmbeddingPreprocessor()
df_scales_num, df_stds_num = embp.build_scales(df_seq=df_seq, dict_num=dict_num,
                                               return_std=True)   # also the per-channel spread
df_cat_num = embp.build_cat(df_scales=df_scales_num,
                            df_stds=df_stds_num,     # spread-aware clustering of the D channels
                            cat_min_th=0.5,          # similarity threshold for a category
                            subcat_min_th=0.7,       # tighter threshold for a subcategory
                            metric="correlation",
                            random_state=42)
aa.display_df(df=df_scales_num, n_rows=10, show_shape=True)

DataFrame shape: (20, 4)
  dim_0 dim_1 dim_2 dim_3
A 0.492584 0.497712 0.485069 0.513540
C 0.502098 0.506278 0.518440 0.502096
D 0.487394 0.508945 0.499097 0.498795
E 0.482504 0.494474 0.511045 0.497436
F 0.508067 0.473032 0.475938 0.501151
G 0.500774 0.508791 0.498488 0.505346
H 0.513809 0.478418 0.495852 0.535164
I 0.501820 0.522848 0.511069 0.500499
K 0.507420 0.494591 0.520385 0.502633
L 0.482386 0.496649 0.504925 0.510880 
# Slice sequences AND tensor with shared part boundaries
df_parts_num, dict_num_parts = aa.NumericalFeature.get_parts(df_seq=df_seq, dict_num=dict_num,
                                                             list_parts=None, all_parts=False,
                                                             jmd_n_len=10, jmd_c_len=10,
                                                             tmd_len=None)

# Each part tensor has shape (n_samples, L_part_max, D), NaN-padded
{k: v.shape for k, v in dict_num_parts.items()}

{'tmd': (20, 23, 4), 'jmd_n_tmd_n': (20, 22, 4), 'tmd_c_jmd_c': (20, 21, 4)}

# Numerical-mode CPP: same algorithm and df_feat schema as CPP.run
cpp = aa.CPP(df_parts=df_parts_num,
             split_kws=split_kws,
             df_scales=df_scales_num,
             df_cat=df_cat_num,
             accept_gaps=False,
             verbose=False,
             random_state=42,
             bootstrap=False, bootstrap_kws=None)   # bootstrap=True adds selection_frequency

# run and run_num take the identical knobs, so set them once and reuse them for
# both modes (the sequence-mode run below the figure uses exactly this dict).
cpp_kws = dict(labels=labels, label_test=1, label_ref=0,
               n_filter=10,                            # size of the returned signature
               n_pre_filter=None, pct_pre_filter=5,    # pre-filter: keep the top 5% by effect size
               max_std_test=0.2,                       # drop features scattering within the test group
               max_overlap=0.5, max_cor=0.5,           # redundancy caps
               check_cat=True,                         # category-aware redundancy filtering
               parametric=False,                       # Mann-Whitney U instead of a t-test
               start=1, tmd_len=20, jmd_n_len=10, jmd_c_len=10,   # position annotation
               n_jobs=1,                               # serial: avoids spawn footgun on Py3.14 + macOS
               vectorized=True, n_batches=None, n_sample_batches=None,   # speed / memory knobs
               return_stats=False,                     # True -> (df_feat, dict_stats)
               redundancy="legacy")
df_feat = cpp.run_num(dict_num_parts=dict_num_parts, **cpp_kws)
aa.display_df(df=df_feat[["feature", "category", "subcategory", "mean_dif", "abs_auc"]],
              n_rows=10, show_shape=True)

DataFrame shape: (10, 5)
  feature category subcategory mean_dif abs_auc
1 TMD_C_JMD_C-Pat...C,4,7,10)-dim_1 Embeddings Embeddings_cat0_subcat1 -0.216000 0.360000
2 TMD_C_JMD_C-Pat...,9,12,15)-dim_1 Embeddings Embeddings_cat0_subcat1 -0.174000 0.360000
3 JMD_N_TMD_N-Pat...,6,10,14)-dim_0 Embeddings Embeddings_cat1_subcat0 0.174000 0.360000
4 JMD_N_TMD_N-Pat...,8,11,14)-dim_3 Embeddings Embeddings_cat0_subcat3 0.165000 0.350000
5 TMD-Pattern(N,2,5,8)-dim_3 Embeddings Embeddings_cat0_subcat3 0.152000 0.340000
6 TMD_C_JMD_C-Seg...ent(9,11)-dim_0 Embeddings Embeddings_cat1_subcat0 -0.207000 0.330000
7 TMD-Pattern(C,3,7,10,13)-dim_2 Embeddings Embeddings_cat2_subcat2 0.206000 0.320000
8 TMD_C_JMD_C-Seg...nt(10,10)-dim_1 Embeddings Embeddings_cat0_subcat1 0.188000 0.320000
9 TMD_C_JMD_C-Pat...2,5,8,11)-dim_2 Embeddings Embeddings_cat2_subcat2 0.143000 0.320000
10 TMD-Pattern(N,2,6,10)-dim_2 Embeddings Embeddings_cat2_subcat2 -0.174000 0.310000 
# Multiple per-residue sources combine along the D axis before get_parts.
# Here a second synthetic channel stands in for a structural tensor.
dict_struct = {e: rng.random((len(s), 2)) for e, s in zip(df_seq["entry"], df_seq["sequence"])}
dict_combined = aa.combine_dict_nums(dict_nums=[dict_num, dict_struct])
{e: dict_combined[e].shape for e in list(dict_combined)[:3]}   # D = 4 + 2 = 6

{'Q14802': (87, 6), 'Q86UE4': (582, 6), 'Q969W9': (287, 6)}

Output. Sequence mode produces three objects: - df_parts: one row per protein, one column per part (tmd, jmd_n_tmd_n, tmd_c_jmd_c); the primary input to CPP. - split_kws: the split configuration dict. - features: a list of Part-Split-Scale feature ids, optionally reduced by filter_correlation.

Numerical mode produces: - (df_parts_num, dict_num_parts) from get_parts; each dict_num_parts[part] is a NaN-padded array of shape (n_samples, L_part_max, D), aligned row-for-row with df_parts_num. - df_feat: the signature, with the same schema as run() (feature, category, subcategory, mean_dif, abs_auc, …). The value axis is now the embedding channels (dim_0 … dim_3) instead of physicochemical scales.

The bar plot below makes the engineered sequence-mode space concrete: it tallies the features enumeration by part (x-axis) and split type (colour). Notice how every part carries the same Segment / Pattern / PeriodicPattern stack: the grammar is uniform across parts, and what multiplies the count is simply how many scales (sequence mode) or numerical channels (numerical mode) you enumerate over. A second figure then puts the two modes side by side, so the difference between them is visible rather than asserted.

import matplotlib.pyplot as plt
import seaborn as sns

# Tally the sequence-mode enumeration: every feature id is a Part-Split-Scale triple,
# so we count features by part (x-axis) and split type (stacked colour).
parts_order = ["TMD", "JMD_N_TMD_N", "TMD_C_JMD_C"]
part_labels = ["TMD", "JMD_N\nTMD_N", "TMD_C\nJMD_C"]
split_types = list(split_kws)                 # ['Segment', 'Pattern', 'PeriodicPattern']
counts = {st: [0] * len(parts_order) for st in split_types}
for feat in features:
    tokens = feat.split("-")
    part, split = tokens[0], tokens[1]
    st = split.split("(")[0]
    if part in parts_order and st in counts:
        counts[st][parts_order.index(part)] += 1

# Publication styling (plotting_prelude): clean stacked bar with despine + tidy legend.
aa.plot_settings(font_scale=0.9, weight_bold=False)
colors = aa.plot_get_clist(n_colors=len(split_types))

fig, ax = plt.subplots(figsize=(6, 5))
x = range(len(parts_order))
bottom = [0] * len(parts_order)
dict_color = {}
for st, color in zip(split_types, colors):
    ax.bar(x, counts[st], bottom=bottom, label=st, color=color,
           width=0.62, edgecolor="white", linewidth=0.8)
    bottom = [b + c for b, c in zip(bottom, counts[st])]
    dict_color[st] = color

total = max(bottom)
ax.set_xticks(list(x))
ax.set_xticklabels(part_labels)
ax.set_ylim(0, total * 1.5)     # headroom for the legend above the bars
ax.set_xlabel("Part")
ax.set_ylabel("Number of features (over 3 scales)")
ax.set_title("Engineered feature space: Part x Split x Scale",
             size=aa.plot_gcfs() + 1)
sns.despine()

# Compact legend in the headroom above the bars (top right, clear of the y ticks and title).
aa.plot_legend(ax=ax, dict_color=dict_color, title="Split type",
               loc="upper right", x=1.0, y=1.0, n_cols=1, handletextpad=0.5)

plt.tight_layout()
plt.show()
../_images/protocol5_engineer_features_1_output_14_0.png

The two modes side by side. The bar plot above counts the engineered space; the figure below shows what actually changes when you switch modes. The two top panels show the same part of the same protein (the TMD of FXYD3), first as the three physicochemical scales that sequence mode averages, then as the four numerical channels that numerical mode averages. Both are per-residue values in [0, 1] and both are consumed by the identical Part-Split grammar, so the value axis is the only thing that differs.

The bottom panel runs each mode end to end at identical settings (the cpp_kws above) and shows the top five features of each signature. Read the feature ids: the first two tokens (part and split) come from the same grammar in both modes, and only the third token changes, a scale id (ARGP820101) in sequence mode, a channel id (dim_1) in numerical mode. Bar length is abs_auc, the separation strength. Sequence mode reaches the higher top effect size here because its value axis carries real physicochemistry, while these numerical channels are synthetic noise; with a real ProtT5 / ESM embedding or a structure tensor, the numerical bars are the ones to read.

# Sequence-mode signature, run with exactly the knobs numerical mode used above:
# same parts, same splits, same filters, so only the value source differs. The
# CPP tool is used sequentially here, so the canonical name is reassigned.
df_scales_seq = df_scales[small_scales]
cpp = aa.CPP(df_parts=df_parts, split_kws=split_kws, df_scales=df_scales_seq,
             df_cat=None, accept_gaps=False, verbose=False, random_state=42,
             bootstrap=False, bootstrap_kws=None)
df_feat_seq = cpp.run(**cpp_kws)

# The value source of each mode, on one part of one protein.
entry = df_parts.index[0]
gene = df_seq.set_index("entry").loc[entry, "gene"]
tmd = df_parts.loc[entry, "tmd"]
i_row = list(df_parts_num.index).index(entry)
M_seq = df_scales_seq.loc[list(tmd)].to_numpy().T          # (3 scales, L)
M_num = dict_num_parts["tmd"][i_row][:len(tmd)].T          # (4 channels, L)

aa.plot_settings(font_scale=0.8, weight_bold=False)
c_seq, c_num = aa.plot_get_clist(n_colors=2)

fig = plt.figure(figsize=(11, 7))
gs = fig.add_gridspec(2, 2, height_ratios=[1, 1.2])
ax_a, ax_b = fig.add_subplot(gs[0, 0]), fig.add_subplot(gs[0, 1])
for ax, M, names, title in [
        (ax_a, M_seq, small_scales, "Sequence mode\n3 physicochemical scales"),
        (ax_b, M_num, list(df_scales_num.columns), "Numerical mode\n4 embedding channels")]:
    sns.heatmap(M, ax=ax, cmap="viridis", vmin=0, vmax=1, xticklabels=list(tmd),
                yticklabels=names, linewidths=0.3, linecolor="white",
                cbar=ax is ax_b, cbar_kws={"label": "Value in [0, 1]", "pad": 0.02})
    ax.set_title(title, size=aa.plot_gcfs())
    ax.tick_params(axis="x", rotation=0, length=0, labelsize=aa.plot_gcfs() - 3)
    ax.tick_params(axis="y", rotation=0, length=0)
    ax.set_xlabel(f"TMD of {gene} ({entry})")

# The two signatures: identical Part-Split ids, different value token.
ax_c = fig.add_subplot(gs[1, :])
n_top = 5
feats = list(df_feat["feature"][:n_top])[::-1] + list(df_feat_seq["feature"][:n_top])[::-1]
vals = list(df_feat["abs_auc"][:n_top])[::-1] + list(df_feat_seq["abs_auc"][:n_top])[::-1]
ax_c.barh(range(len(feats)), vals, color=[c_num] * n_top + [c_seq] * n_top,
          edgecolor="white", linewidth=0.8, height=0.72)
ax_c.axhline(n_top - 0.5, color="0.7", lw=0.8)
ax_c.set_yticks(range(len(feats)))
ax_c.set_yticklabels(feats, size=aa.plot_gcfs() - 3)
ax_c.set_xlabel("abs_auc (separation strength)")
ax_c.set_xlim(0, max(vals) * 1.12)
ax_c.set_title("Top 5 features per mode: same Part-Split ids, different value token",
               size=aa.plot_gcfs())
sns.despine(ax=ax_c)
aa.plot_legend(ax=ax_c, dict_color={"Sequence mode (CPP.run)": c_seq,
                                    "Numerical mode (CPP.run_num)": c_num},
               n_cols=2, loc="upper left", loc_out=True, x=0, y=-0.42,
               fontsize=aa.plot_gcfs() - 2)

plt.tight_layout()
plt.show()
../_images/protocol5_engineer_features_2_output_16_0.png

How to interpret.

Output

Reading

a df_parts row

the default parts (the base part tmd plus the composites jmd_n_tmd_n and tmd_c_jmd_c, derived from the base parts jmd_n / tmd / jmd_c) of one protein

a split_kws entry

the rule for which positions are averaged (Segment = contiguous block, Pattern / PeriodicPattern = spaced positions)

a Part-Split-Scale id (sequence mode)

where + how + which physicochemical property

a Part-Split-dimension id (numerical mode)

where + how + which embedding / structure channel (dim_0 …)

mean_dif sign

direction of separation (positive = higher in the test group in that region)

abs_auc

effect size / separation strength of the feature

Both modes read the same way: a coherent block of high-abs_auc features from one property family or embedding channel, localized to one part, is the signature that distinguishes substrates from non-substrates. Numerical mode simply lets you profile learned / structural representations with the same interpretable grammar. One caveat for this demo: because the numerical-mode tensor is synthetic NumPy noise, the abs_auc / mean_dif values shown above are not biologically meaningful (the small effect sizes are pure chance). Swap in a real ProtT5 / ESM embedding or a DSSP / structure tensor of the same (L, D) shape and the identical readout surfaces the channels that actually separate the test group from the reference group.

The side-by-side figure above is the compact version of that reading: one part, two value axes, two signatures whose ids differ only in the final token.

Key takeaways - A feature is always part × split × scale: control where, how, and what value independently, and the engineered space follows mechanically. - Sequence mode (run(), physicochemical scales) and numerical mode (run_num(), per-residue tensors via dict_num) share one grammar and one df_feat schema, so embeddings and structure ride the same interpretable machinery. - The shape of the space is yours to set: enumerate broadly to explore, then reduce redundancy with filter_correlation before committing to a CPP run.

Common mistakes.

  • Passing ``df_seq`` into :class:`~aaanalysis.CPP`: CPP takes df_parts; build them with get_df_parts (sequence mode) or get_parts (numerical mode) first.

  • Feeding raw, unbounded embeddings into ``run_num``: values must be in [0, 1] (the max_std_test pre-filter is calibrated for that range); pass real PLM embeddings through encode() first. Our synthetic tensor is already in range, so it skips that step.

  • ``dict_num`` length mismatch: each dict_num[entry] must have L == len(sequence); otherwise get_parts cannot align the tensor to the parts.

  • ``D`` mismatch: D must equal len(df_scales.columns) and len(df_cat); deriving both with build_scales/build_cat from the same dict_num keeps the invariant.

  • Reading effect sizes across modes as a ranking of methods: in the figure above the numerical channels are synthetic noise, so their abs_auc is chance. Compare the two modes only when both value axes are real.

  • Using ``len(df_seq)`` for class sizes: use the label column; load_dataset(..., n=N) returns 2N rows.

  • Relying on the default ``n_jobs``: on Python 3.14 + macOS, multiprocessing spawn without a __main__ guard is fragile; pass n_jobs=1 in notebooks.

  • Enumerating the full scale set: get_features over hundreds of scales explodes combinatorially; slice list_scales for exploration.

Next step.

  • Compositional vs positional leads to P6: Compositional vs positional. You now know a feature is Part-Split-Scale; P6 shows how the split alone decides whether a feature is composition-like (whole-part Segment(1,1), position-agnostic) or position-resolved (Segment with n_split_max > 1, Pattern, PeriodicPattern), and how that choice maps onto the prediction level.

  • Run the signature by feeding the engineered df_parts straight into the CPP signature: see Protocol 1 - CPP signature.