Lesson 4 of 13 · 12 min

Structural, Physicochemical, and Expression Features

Sequence tokens are not the only signal a biological model can learn from. Structural geometry, per-residue physicochemical properties, and transcriptome-wide expression carry orthogonal information that often lifts performance well beyond one-hot or k-mer encodings. This lesson builds fixed-length feature matrices from each of these sources and handles the scaling choices that make or break an omics model.

Amino-Acid Physicochemical Scales

A physicochemical scale maps each of the 20 amino acids to a scalar, such as Kyte-Doolittle hydropathy, Eisenberg mean hydrophobicity, or a volume, charge, or flexibility index from the AAindex database. To turn a variable-length protein into a fixed vector you aggregate these per-residue values, typically the mean, plus composition and a few global descriptors. Biopython's ProteinAnalysis exposes GRAVY (mean Kyte-Doolittle), isoelectric point, instability index, aromaticity, and helix/turn/sheet propensity fractions directly.

python

import numpy as np
from Bio.SeqUtils.ProtParam import ProteinAnalysis

np.set_printoptions(suppress=True)

seqs = {
    "P0DTC2": "MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHS",
    "P69905": "MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLS",
}
AA = "ACDEFGHIKLMNPQRSTVWY"

def physchem_features(seq):
    pa = ProteinAnalysis(seq)
    comp = pa.amino_acids_percent                # 20-dim composition (percent)
    helix, turn, sheet = pa.secondary_structure_fraction()
    scalars = [pa.gravy(), pa.isoelectric_point(),
               pa.instability_index(), pa.aromaticity(),
               helix, turn, sheet, pa.molecular_weight()]
    return np.array([comp[a] for a in AA] + scalars)

X = np.vstack([physchem_features(s) for s in seqs.values()])
print(X.shape)
print(np.round(X[0, 20:], 3))               # the 8 global descriptors

(2, 28)
[   0.326    9.623   33.15     0.14     0.2      0.28     0.52  5694.582]

Structure And Accessibility

python

from Bio.PDB import PDBParser, DSSP
import numpy as np

parser = PDBParser(QUIET=True)
model = parser.get_structure("1UBQ", "1ubq.pdb")[0]
dssp = DSSP(model, "1ubq.pdb", dssp="mkdssp")   # needs the mkdssp binary on PATH

# collapse DSSP 8-state code (index 2) to 3-state; relative ASA is index 3
three = {"H": "H", "G": "H", "I": "H", "E": "E", "B": "E",
         "T": "C", "S": "C", "P": "C", "-": "C"}
ss  = [three[dssp[k][2]] for k in dssp.keys()]
rsa = np.array([dssp[k][3] for k in dssp.keys()], dtype=float)
n = len(ss)

feat = {
    "frac_H": ss.count("H") / n,
    "frac_E": ss.count("E") / n,
    "frac_C": ss.count("C") / n,
    "mean_rsa": float(np.nanmean(rsa)),
    "frac_buried": float(np.mean(rsa < 0.25)),   # RSA<0.25 = core residue
}
print({k: round(v, 3) for k, v in feat.items()})   # frac_H + frac_E + frac_C == 1.0

Expression And Scaling

python

import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.feature_selection import VarianceThreshold

rng = np.random.default_rng(0)
counts = rng.negative_binomial(5, 0.3, size=(40, 20000)).astype(float)  # n x p, p >> n
counts[:, :4000] = 0.0            # 4000 genes are unexpressed in every sample
y = rng.integers(0, 2, size=40)

# counts -> CPM (per-sample library-size scaling) -> log1p(CPM)
lib = counts.sum(axis=1, keepdims=True)
logcpm = np.log1p(counts / lib * 1e6)

X_tr, X_te, y_tr, y_te = train_test_split(
    logcpm, y, test_size=0.25, stratify=y, random_state=0)

# fit EVERY transform on train only, then apply to test (no leakage)
vt = VarianceThreshold(1e-3).fit(X_tr)          # drop near-constant genes
X_tr, X_te = vt.transform(X_tr), vt.transform(X_te)
sc = StandardScaler().fit(X_tr)                  # per-gene z-score, train stats
X_tr, X_te = sc.transform(X_tr), sc.transform(X_te)

print("train", X_tr.shape, "test", X_te.shape)
print("train mean", round(X_tr.mean(), 6), "std", round(X_tr.std(), 3))
print("test mean", round(X_te.mean(), 3))

train (30, 16000) test (10, 16000)
train mean 0.0 std 1.0
test mean -0.02

Warning: with p >> n (thousands of genes, dozens of samples) any transform that learns cross-sample statistics or the labels leaks test information and inflates accuracy. Per-sample CPM/log can be applied up front because it never crosses samples, but variance filtering, univariate selection, and StandardScaler must be fit inside each cross-validation fold via a Pipeline, and you should lean on regularization (L1/L2 or elastic net) or supervised dimensionality reduction rather than adding raw features. Report nested-CV scores, not a single split.

Try it yourself: assemble a labeled set of proteins from two families (say ~100 enzymes vs ~100 non-enzymes with solved structures), run mkdssp on each PDB, and concatenate every protein's physchem vector with its DSSP structure fractions and mean_rsa into one feature matrix. Fit an L2-regularized LogisticRegression inside a StratifiedKFold and record ROC-AUC. Then measure how the score changes when you (a) drop StandardScaler entirely and (b) fit the scaler on the whole matrix before the split instead of inside each fold, and explain why only (b) is data leakage.

Encoding Sequences: k-mers and One-Hot Train/Val/Test Splits and Sequence-Similarity Leakage