UncleanCode commited on
Commit
2f40071
·
verified ·
1 Parent(s): 6726e6c

StaMP v2 release: 5-seed ensemble weights, ProteinMPNN backbone, inference, honest model card

Browse files
LICENSE ADDED
@@ -0,0 +1,21 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ MIT License
2
+
3
+ Copyright (c) 2026 Modular
4
+
5
+ Permission is hereby granted, free of charge, to any person obtaining a copy
6
+ of this software and associated documentation files (the "Software"), to deal
7
+ in the Software without restriction, including without limitation the rights
8
+ to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
9
+ copies of the Software, and to permit persons to whom the Software is
10
+ furnished to do so, subject to the following conditions:
11
+
12
+ The above copyright notice and this permission notice shall be included in all
13
+ copies or substantial portions of the Software.
14
+
15
+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16
+ IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17
+ FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
18
+ AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19
+ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
20
+ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
21
+ SOFTWARE.
README.md ADDED
@@ -0,0 +1,128 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ ---
2
+ license: mit
3
+ library_name: pytorch
4
+ tags:
5
+ - biology
6
+ - protein
7
+ - protein-stability
8
+ - ddg
9
+ - mutation-effect
10
+ - proteinmpnn
11
+ pipeline_tag: other
12
+ ---
13
+
14
+ # StaMP — Stability from Message-Passing
15
+
16
+ **StaMP** is a small, fast, structure-conditioned predictor of **single-point mutation stability**
17
+ (ΔΔG, kcal/mol) for proteins. Give it a backbone structure and it scores how much any single
18
+ amino-acid substitution stabilises or destabilises the fold, or scans **every** possible
19
+ substitution at every residue (a full in-silico deep-mutational scan) in a few seconds on CPU.
20
+
21
+ > **Convention:** positive ΔΔG = **destabilizing**, negative = stabilizing. ΔΔG(x→x) = 0 by construction.
22
+
23
+ It is built for the first filtering step of a protein-engineering campaign (enzymes, antibodies,
24
+ therapeutics, binders): rank thousands of candidate point mutations cheaply, then take only the
25
+ promising ones to the bench.
26
+
27
+ - **Trainable parameters:** 591,893 (~592k). Released as a **5-seed ensemble**.
28
+ - **Backbone:** frozen ProteinMPNN `v_48_020` (MIT, Dauparas et al. 2022), used only as a feature extractor — nothing in it is trained.
29
+ - **Speed:** a full single-mutation scan of a ~60-residue domain is ~3 s on 8 CPU threads (features dominate; the head is negligible).
30
+
31
+ ## How it works
32
+
33
+ All heavy components are **frozen**; only a tiny head is trained.
34
+
35
+ 1. **Frozen ProteinMPNN features**, computed once per chain → `R[L, 533]`:
36
+ - the 3 encoder-layer outputs (3×128, structure-only),
37
+ - the `W_s(wt)` sequence embedding (128),
38
+ - conditional log-probs averaged over 8 random decoding orders (21) — a zero-shot stability signal.
39
+ 2. **Structural neighbour context:** for each residue, the **48 nearest residues** on ProteinMPNN's own Cα k-NN graph are mean-pooled and concatenated with the centre features plus two burial counts (#neighbours within 8 Å and 12 Å).
40
+ 3. **Gated head (592k params):** a gated linear unit → LayerNorm → MLP emits one score per amino acid; the prediction is `o[mut] − o[wt]`, which forces ΔΔG(x→x)=0 exactly.
41
+
42
+ The neighbour context is what separates StaMP from a per-residue baseline: it adds **+0.026 Spearman**
43
+ on the Megascale test split and carries over to the external benchmarks.
44
+
45
+ ## Benchmarks
46
+
47
+ Trained and evaluated on the **Megascale** dataset (Tsuboyama et al., 2023), using the homology-clustered
48
+ train/val/test split released by the ThermoMPNN authors (mmseqs2 at 25% identity; 239/31/28 domains).
49
+ 5-seed ensemble, selected on validation Spearman (never on test).
50
+
51
+ ### Megascale held-out test (in-domain, n = 28,172)
52
+
53
+ | Metric | StaMP (v2, 5-seed) | ThermoMPNN (ref) |
54
+ |---|---|---|
55
+ | Spearman | **0.734** | ~0.73 |
56
+ | Pearson | **0.767** | ~0.75 |
57
+ | RMSE (kcal/mol) | 0.705 | ~0.71 |
58
+
59
+ Per-domain mean Spearman 0.786 / Pearson 0.823. StaMP slightly edges ThermoMPNN's headline
60
+ Megascale Spearman with **~5× fewer trainable parameters** (592k vs ~2.7M).
61
+
62
+ ### External generalisation (not seen in training)
63
+
64
+ Structures from RCSB; mutations mapped by PDB residue number. We also apply our **own homology filter**:
65
+ any benchmark protein with ≥25% identity to a Megascale *training* domain is removed (the unfiltered
66
+ numbers that are usually quoted are inflated by homology).
67
+
68
+ | Benchmark | n | Spearman | Pearson |
69
+ |---|---|---|---|
70
+ | S669 (all) | 667 | 0.453 | 0.415 |
71
+ | **S669 (homology-filtered)** | 488 | **0.384** | **0.368** |
72
+ | Ssym-direct | 342 | 0.655 | 0.666 |
73
+ | *ThermoMPNN, Ssym-direct (published)* | 342 | — | *0.72* |
74
+
75
+ ## ⚠️ Honest caveats — read before using
76
+
77
+ - **In-domain vs zero-shot.** The headline Megascale numbers are **in-domain**: train and test are different
78
+ proteins but the same assay and the same kind of small (≤72-residue) AlphaFold-modelled domains. This is
79
+ **not** a zero-shot number.
80
+ - **Performance drops out of domain.** On homology-filtered S669 (real crystal structures of larger, often
81
+ multi-domain proteins) Spearman falls to ~0.38. On Ssym StaMP (0.67 Pearson) is **below** ThermoMPNN's
82
+ published 0.72. StaMP trained only on AlphaFold models of tiny domains, so it is weaker on crystal
83
+ structures of full proteins. We tested coordinate-noise augmentation to close this gap and it did **not**
84
+ help — the gap is a training-data limitation, not a robustness trick.
85
+ - **Direction.** Only direct (forward) mutations are validated; inverse-mutation (Ssym-inverse) needs mutant
86
+ structures and was not evaluated.
87
+ - **Not a replacement for ΔΔG from physics or experiment.** Use StaMP as a cheap ranking / prescreen, not as a
88
+ calibrated free energy.
89
+
90
+ ## Usage
91
+
92
+ ```bash
93
+ pip install torch numpy
94
+ # predict one mutation (position is 1-based along the chain)
95
+ python stamp/predict.py --pdb your.pdb --chain A --mut A42G
96
+ # ...or use the PDB residue number
97
+ python stamp/predict.py --pdb your.pdb --chain A --mut A42G --pdb_numbering
98
+ # full single-mutation saturation scan -> TSV, ranked
99
+ python stamp/predict.py --pdb your.pdb --chain A --scan --out your_scan.tsv
100
+ ```
101
+
102
+ `predict.py` expects a **single-chain** backbone (N, CA, C, O); extract your chain of interest first.
103
+ The ProteinMPNN backbone weights (`weights/v_48_020.pt`) and the five head checkpoints
104
+ (`weights/stamp_s*.pt`) ship in this repo.
105
+
106
+ ## Files
107
+
108
+ ```
109
+ stamp/predict.py end-to-end inference (single mutation or full scan)
110
+ stamp/protein_mpnn_utils.py ProteinMPNN model + PDB featurizer (MIT, Dauparas et al.)
111
+ stamp/weights/stamp_s{0..4}.pt 5 trained head checkpoints (each stores its own cfg + feature norms)
112
+ stamp/weights/v_48_020.pt frozen ProteinMPNN backbone (MIT)
113
+ ```
114
+
115
+ ## License
116
+
117
+ MIT. StaMP's head and code are released under MIT. The bundled ProteinMPNN code and `v_48_020.pt`
118
+ weights are from the ProteinMPNN project (Dauparas et al., *Science* 2022), also MIT.
119
+
120
+ ## Citation
121
+
122
+ If you use StaMP, please also cite the work it builds on:
123
+
124
+ - Dauparas et al., "Robust deep learning–based protein sequence design using ProteinMPNN," *Science*, 2022.
125
+ - Tsuboyama et al., "Mega-scale experimental analysis of protein folding stability in biology and design," *Nature*, 2023.
126
+ - Dieckhaus et al. (ThermoMPNN), whose Megascale split and benchmarking protocol StaMP follows.
127
+
128
+ Built by [Modular](https://huggingface.co/Modularcomputing).
requirements.txt ADDED
@@ -0,0 +1,2 @@
 
 
 
1
+ torch>=2.0
2
+ numpy>=1.23
stamp/predict.py ADDED
@@ -0,0 +1,274 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ """StaMP — Stability from Message-Passing.
2
+
3
+ A small, fast structure-conditioned predictor of single-point mutation stability (ddG,
4
+ kcal/mol) built on frozen ProteinMPNN features. Positive ddG = destabilizing.
5
+
6
+ Given a PDB/mmCIF backbone, StaMP scores any single amino-acid substitution, or scans
7
+ every possible substitution at every residue (deep mutational scan) in a few seconds on CPU.
8
+
9
+ Pipeline (all frozen except the small head):
10
+ 1. ProteinMPNN v_48_020 (MIT, Dauparas et al. 2022) is run once per chain as a feature
11
+ extractor: 3 encoder-layer outputs (3x128), the W_s(wt) embedding (128), and
12
+ conditional log-probs averaged over 8 random decoding orders (21). -> R[L, 533]
13
+ 2. For each residue, the 48 nearest structural neighbours (ProteinMPNN's Ca k-NN graph)
14
+ are mean-pooled and concatenated with the centre features and two burial counts.
15
+ 3. A 592k-param gated head emits one score per amino acid; ddG = o[mut] - o[wt], so
16
+ ddG(x->x) = 0 exactly. The released model is a 5-seed ensemble (mean of the heads).
17
+
18
+ Usage:
19
+ python predict.py --pdb 1abc.pdb --chain A --mut A42G
20
+ python predict.py --pdb 1abc.pdb --chain A --scan # full single-mutation scan
21
+ See README.md for benchmarks and caveats.
22
+ """
23
+ import os, argparse, glob, numpy as np, torch
24
+ from protein_mpnn_utils import ProteinMPNN
25
+
26
+ ALPHABET = 'ACDEFGHIKLMNPQRSTVWYX'
27
+ HERE = os.path.dirname(os.path.abspath(__file__))
28
+ THREE = {'ALA': 'A', 'ARG': 'R', 'ASN': 'N', 'ASP': 'D', 'CYS': 'C', 'GLN': 'Q', 'GLU': 'E',
29
+ 'GLY': 'G', 'HIS': 'H', 'ILE': 'I', 'LEU': 'L', 'LYS': 'K', 'MET': 'M', 'PHE': 'F',
30
+ 'PRO': 'P', 'SER': 'S', 'THR': 'T', 'TRP': 'W', 'TYR': 'Y', 'VAL': 'V'}
31
+
32
+
33
+ def parse_chain(path, chain):
34
+ """Minimal backbone parser for one chain: returns (seq, resnums, X[L,4,3]).
35
+ Keeps only residues with all of N, CA, C, O. resnums are the PDB residue labels
36
+ (resseq+icode) in order, so --mut can be given in PDB numbering with --pdb_numbering."""
37
+ res = {}
38
+ for line in open(path):
39
+ if not line.startswith('ATOM'):
40
+ continue
41
+ if chain and line[21] != chain:
42
+ continue
43
+ if line[16] not in (' ', 'A'):
44
+ continue
45
+ name = line[12:16].strip()
46
+ if name not in ('N', 'CA', 'C', 'O'):
47
+ continue
48
+ rn = line[22:27].strip()
49
+ aa = THREE.get(line[17:20].strip())
50
+ if aa is None:
51
+ continue
52
+ xyz = [float(line[30:38]), float(line[38:46]), float(line[46:54])]
53
+ res.setdefault(rn, {'aa': aa, 'at': {}})['at'][name] = xyz
54
+ keys = [k for k in res if len(res[k]['at']) == 4]
55
+ keys.sort(key=lambda k: (int(''.join(c for c in k if c.isdigit() or c == '-')), k))
56
+ seq = ''.join(res[k]['aa'] for k in keys)
57
+ X = np.array([[res[k]['at'][a] for a in ('N', 'CA', 'C', 'O')] for k in keys], np.float32)
58
+ return seq, keys, X
59
+
60
+
61
+ # ----- frozen ProteinMPNN feature extractor -------------------------------------------
62
+ def load_mpnn(weights, device):
63
+ ckpt = torch.load(weights, map_location='cpu', weights_only=False)
64
+ m = ProteinMPNN(ca_only=False, num_letters=21, node_features=128, edge_features=128,
65
+ hidden_dim=128, num_encoder_layers=3, num_decoder_layers=3,
66
+ k_neighbors=ckpt['num_edges'], augment_eps=0.0)
67
+ m.load_state_dict(ckpt['model_state_dict'])
68
+ m.to(device).eval()
69
+ for p in m.parameters():
70
+ p.requires_grad = False
71
+ return m
72
+
73
+
74
+ @torch.no_grad()
75
+ def encode(model, X, S, mask, chain_M, residue_idx, chain_enc):
76
+ E, E_idx = model.features(X, mask, residue_idx, chain_enc)
77
+ h_V = torch.zeros((E.shape[0], E.shape[1], E.shape[-1]), device=E.device, dtype=E.dtype)
78
+ h_E = model.W_e(E)
79
+ ma = torch.gather(mask, 1, E_idx.view(E_idx.shape[0], -1)).view(E_idx.shape)
80
+ ma = mask.unsqueeze(-1) * ma
81
+ outs = []
82
+ for layer in model.encoder_layers:
83
+ h_V, h_E = layer(h_V, h_E, E_idx, mask, ma)
84
+ outs.append(h_V)
85
+ return torch.cat(outs, -1), model.W_s(S)
86
+
87
+
88
+ @torch.no_grad()
89
+ def cond_logprobs(model, X, S, mask, chain_M, residue_idx, chain_enc, n_orders, gen):
90
+ acc = None
91
+ for _ in range(n_orders):
92
+ randn = torch.randn(chain_M.shape, device=X.device, generator=gen)
93
+ lp = model.conditional_probs(X, S, mask, chain_M * mask, residue_idx, chain_enc, randn)
94
+ lp = lp[:, :, 0, :] if lp.dim() == 4 else lp
95
+ acc = lp if acc is None else torch.logaddexp(acc, lp)
96
+ return acc - np.log(n_orders)
97
+
98
+
99
+ def knn(X, k=48):
100
+ Ca = torch.as_tensor(X)[:, 1].float()
101
+ D = torch.sqrt(((Ca[:, None] - Ca[None]) ** 2).sum(-1) + 1e-6)
102
+ k = min(k, Ca.shape[0])
103
+ d, idx = torch.topk(D, k, dim=-1, largest=False)
104
+ return idx.numpy().astype(np.int64), d.numpy().astype(np.float32)
105
+
106
+
107
+ # ----- trained head -------------------------------------------------------------------
108
+ import torch.nn as nn
109
+
110
+
111
+ class MSHead(nn.Module):
112
+ def __init__(self, in_dim, gate_dim=256, hidden=(128, 64), vocab=21, dropout=0.0):
113
+ super().__init__()
114
+ self.W1 = nn.Linear(in_dim, gate_dim)
115
+ self.W2 = nn.Linear(in_dim, gate_dim)
116
+ self.norm = nn.LayerNorm(gate_dim)
117
+ layers, d = [], gate_dim
118
+ for h in hidden:
119
+ layers += [nn.Linear(d, h), nn.ReLU()]
120
+ if dropout:
121
+ layers += [nn.Dropout(dropout)]
122
+ d = h
123
+ layers += [nn.Linear(d, vocab)]
124
+ self.mlp = nn.Sequential(*layers)
125
+
126
+ def forward(self, x, wt_idx, mut_idx):
127
+ h = self.norm(self.W1(x) * torch.sigmoid(self.W2(x)))
128
+ o = self.mlp(h)
129
+ return o.gather(1, mut_idx[:, None]).squeeze(1) - o.gather(1, wt_idx[:, None]).squeeze(1)
130
+
131
+ def scores(self, x):
132
+ return self.mlp(self.norm(self.W1(x) * torch.sigmoid(self.W2(x))))
133
+
134
+
135
+ class MSNbrHead(nn.Module):
136
+ """Released config: mode='mean', k=48. (attn path kept for completeness.)"""
137
+ def __init__(self, xc_dim, r_dim, mode='mean', k=48, gate_dim=256, hidden=(128, 64),
138
+ dropout=0.0, heads=4, dh=32, ctx_dim=128):
139
+ super().__init__()
140
+ self.mode, self.k = mode, k
141
+ in_dim = xc_dim + r_dim + 2 if mode == 'mean' else xc_dim
142
+ self.head = MSHead(in_dim, gate_dim, hidden, dropout=dropout)
143
+
144
+ def context(self, xc, Rn, dist, valid):
145
+ burial = torch.stack([((dist < 8) & valid).float().sum(1),
146
+ ((dist < 12) & valid).float().sum(1)], -1) / 10.0
147
+ Rk, vk = Rn[:, 1:self.k + 1], valid[:, 1:self.k + 1]
148
+ w = vk.float()[..., None]
149
+ m = (Rk * w).sum(1) / w.sum(1).clamp(min=1)
150
+ return torch.cat([xc, m, burial], -1)
151
+
152
+ def forward(self, xc, Rn, dist, valid, wt_idx, mut_idx):
153
+ return self.head(self.context(xc, Rn, dist, valid), wt_idx, mut_idx)
154
+
155
+ def scores(self, xc, Rn, dist, valid):
156
+ return self.head.scores(self.context(xc, Rn, dist, valid))
157
+
158
+
159
+ def load_ensemble(ckpt_dir, device):
160
+ members = []
161
+ for cp in sorted(glob.glob(os.path.join(ckpt_dir, 'stamp_s*.pt'))):
162
+ ck = torch.load(cp, map_location='cpu', weights_only=False)
163
+ c = ck['cfg']
164
+ h = MSNbrHead(ck['mu'].shape[1], ck['r_mu'].shape[1], c['mode'], c['k'], c['gate_dim'],
165
+ tuple(int(x) for x in c['hidden'].split(',')), c['dropout'])
166
+ h.load_state_dict(ck['state_dict']); h.to(device).eval()
167
+ members.append((h, ck['mu'], ck['sd'], ck['r_mu'], ck['r_sd']))
168
+ if not members:
169
+ raise FileNotFoundError(f'no stamp_s*.pt checkpoints in {ckpt_dir}')
170
+ return members
171
+
172
+
173
+ # ----- end-to-end ---------------------------------------------------------------------
174
+ @torch.no_grad()
175
+ def featurize_chain(mpnn, pdb_path, chain, device):
176
+ seq, resnums, Xnp = parse_chain(pdb_path, chain)
177
+ if len(seq) < 1:
178
+ raise ValueError(f'no complete-backbone residues for chain {chain!r} in {pdb_path}')
179
+ L = len(seq)
180
+ X = torch.tensor(Xnp, device=device)[None]
181
+ S = torch.tensor([[ALPHABET.index(a) if a in ALPHABET else 20 for a in seq]], device=device)
182
+ mask = torch.ones((1, L), device=device)
183
+ chain_M = torch.ones((1, L), device=device)
184
+ residue_idx = torch.arange(L, device=device)[None]
185
+ chain_enc = torch.ones((1, L), dtype=torch.long, device=device)
186
+ enc, hseq = encode(mpnn, X, S, mask, chain_M, residue_idx, chain_enc)
187
+ gen = torch.Generator(device=device).manual_seed(0)
188
+ lp = cond_logprobs(mpnn, X, S, mask, chain_M, residue_idx, chain_enc, 8, gen)
189
+ enc, hseq, lp = enc[0].float().cpu().numpy(), hseq[0].float().cpu().numpy(), lp[0].float().cpu().numpy()
190
+ R = np.concatenate([enc, hseq, lp], -1).astype(np.float32) # [L, 533]
191
+ E, D = knn(Xnp, 48)
192
+ return seq, resnums, R, E, D
193
+
194
+
195
+ def predict_muts(members, seq, R, E, D, muts, device):
196
+ """muts: list of (wt, pos0, mut). Returns ensemble-mean ddG for each."""
197
+ xs, Rns, ds, wts, mus, valids = [], [], [], [], [], []
198
+ for wt, i, mut in muts:
199
+ wl, ml = R[i][-21:][ALPHABET.index(wt)], R[i][-21:][ALPHABET.index(mut)]
200
+ xs.append(np.concatenate([R[i], [wl, ml, ml - wl]]).astype(np.float32))
201
+ Rns.append(R[E[i]]); ds.append(D[i]); valids.append(D[i] < 98)
202
+ wts.append(ALPHABET.index(wt)); mus.append(ALPHABET.index(mut))
203
+ X = np.stack(xs); Rn = np.stack(Rns); Dn = np.stack(ds)
204
+ wt_i = torch.tensor(wts, device=device); mu_i = torch.tensor(mus, device=device)
205
+ valid = torch.tensor(np.stack(valids), device=device); Dt = torch.tensor(Dn, device=device)
206
+ preds = []
207
+ with torch.no_grad():
208
+ for h, m, s, rm, rs in members:
209
+ preds.append(h(torch.tensor((X - m) / s).float().to(device),
210
+ torch.tensor((Rn - rm[None]) / rs[None]).float().to(device),
211
+ Dt, valid, wt_i, mu_i).cpu().numpy())
212
+ return np.mean(preds, 0)
213
+
214
+
215
+ def main():
216
+ ap = argparse.ArgumentParser(description='StaMP single-mutation ddG predictor')
217
+ ap.add_argument('--pdb', required=True)
218
+ ap.add_argument('--chain', default='A')
219
+ ap.add_argument('--mut', help='e.g. A42G (wt, position, mut). Position is 1-based along the '
220
+ 'chain by default, or the PDB residue number with --pdb_numbering.')
221
+ ap.add_argument('--pdb_numbering', action='store_true',
222
+ help='interpret --mut position as the PDB residue number, not 1-based index')
223
+ ap.add_argument('--scan', action='store_true', help='score every single substitution')
224
+ ap.add_argument('--ckpt_dir', default=os.path.join(HERE, 'weights'))
225
+ ap.add_argument('--mpnn', default=os.path.join(HERE, 'weights', 'v_48_020.pt'))
226
+ ap.add_argument('--device', default='cuda' if torch.cuda.is_available() else 'cpu')
227
+ ap.add_argument('--out', default=None)
228
+ a = ap.parse_args()
229
+
230
+ dev = torch.device(a.device)
231
+ mpnn = load_mpnn(a.mpnn, dev)
232
+ members = load_ensemble(a.ckpt_dir, dev)
233
+ seq, resnums, R, E, D = featurize_chain(mpnn, a.pdb, a.chain, dev)
234
+
235
+ if a.mut:
236
+ wt, posstr, mut = a.mut[0], a.mut[1:-1], a.mut[-1]
237
+ if a.pdb_numbering:
238
+ assert posstr in resnums, f'residue {posstr} not found in chain {a.chain}'
239
+ i = resnums.index(posstr)
240
+ else:
241
+ i = int(posstr) - 1
242
+ assert 0 <= i < len(seq), f'position {posstr} out of range (len {len(seq)})'
243
+ assert seq[i] == wt, f'expected {wt} at position {posstr}, chain has {seq[i]}'
244
+ ddg = predict_muts(members, seq, R, E, D, [(wt, i, mut)], dev)[0]
245
+ print(f'{a.mut}\tddG = {ddg:+.3f} kcal/mol\t({"destabilizing" if ddg>0 else "stabilizing"})')
246
+ elif a.scan:
247
+ muts, labels = [], []
248
+ for i, wt in enumerate(seq):
249
+ if wt not in ALPHABET[:20]:
250
+ continue
251
+ for mut in ALPHABET[:20]:
252
+ if mut == wt:
253
+ continue
254
+ muts.append((wt, i, mut)); labels.append(f'{wt}{i+1}{mut}')
255
+ ddg = predict_muts(members, seq, R, E, D, muts, dev)
256
+ rows = sorted(zip(labels, ddg), key=lambda r: r[1])
257
+ out = a.out or (os.path.splitext(os.path.basename(a.pdb))[0] + f'_{a.chain}_scan.tsv')
258
+ with open(out, 'w') as f:
259
+ f.write('mutation\tddG_kcal_mol\n')
260
+ for lab, d in sorted(zip(labels, ddg), key=lambda r: r[0]):
261
+ f.write(f'{lab}\t{d:.4f}\n')
262
+ print(f'scanned {len(muts)} mutations over {len(seq)} residues -> {out}')
263
+ print('\nMost stabilizing:')
264
+ for lab, d in rows[:10]:
265
+ print(f' {lab}\t{d:+.3f}')
266
+ print('Most destabilizing:')
267
+ for lab, d in rows[-10:]:
268
+ print(f' {lab}\t{d:+.3f}')
269
+ else:
270
+ ap.error('give --mut A42G or --scan')
271
+
272
+
273
+ if __name__ == '__main__':
274
+ main()
stamp/protein_mpnn_utils.py ADDED
@@ -0,0 +1,1383 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ from __future__ import print_function
2
+ import json, time, os, sys, glob
3
+ import shutil
4
+ import numpy as np
5
+ import torch
6
+ from torch import optim
7
+ from torch.utils.data import DataLoader
8
+ from torch.utils.data.dataset import random_split, Subset
9
+
10
+ import copy
11
+ import torch.nn as nn
12
+ import torch.nn.functional as F
13
+ import random
14
+ import itertools
15
+
16
+ #A number of functions/classes are adopted from: https://github.com/jingraham/neurips19-graph-protein-design
17
+
18
+ def parse_fasta(filename,limit=-1, omit=[]):
19
+ header = []
20
+ sequence = []
21
+ lines = open(filename, "r")
22
+ for line in lines:
23
+ line = line.rstrip()
24
+ if line[0] == ">":
25
+ if len(header) == limit:
26
+ break
27
+ header.append(line[1:])
28
+ sequence.append([])
29
+ else:
30
+ if omit:
31
+ line = [item for item in line if item not in omit]
32
+ line = ''.join(line)
33
+ line = ''.join(line)
34
+ sequence[-1].append(line)
35
+ lines.close()
36
+ sequence = [''.join(seq) for seq in sequence]
37
+ return np.array(header), np.array(sequence)
38
+
39
+ def _scores(S, log_probs, mask):
40
+ """ Negative log probabilities """
41
+ criterion = torch.nn.NLLLoss(reduction='none')
42
+ loss = criterion(
43
+ log_probs.contiguous().view(-1,log_probs.size(-1)),
44
+ S.contiguous().view(-1)
45
+ ).view(S.size())
46
+ scores = torch.sum(loss * mask, dim=-1) / torch.sum(mask, dim=-1)
47
+ return scores
48
+
49
+ def _S_to_seq(S, mask):
50
+ alphabet = 'ACDEFGHIKLMNPQRSTVWYX'
51
+ seq = ''.join([alphabet[c] for c, m in zip(S.tolist(), mask.tolist()) if m > 0])
52
+ return seq
53
+
54
+ def parse_PDB_biounits(x, atoms=['N','CA','C'], chain=None):
55
+ '''
56
+ input: x = PDB filename
57
+ atoms = atoms to extract (optional)
58
+ output: (length, atoms, coords=(x,y,z)), sequence
59
+ '''
60
+
61
+ alpha_1 = list("ARNDCQEGHILKMFPSTWYV-")
62
+ states = len(alpha_1)
63
+ alpha_3 = ['ALA','ARG','ASN','ASP','CYS','GLN','GLU','GLY','HIS','ILE',
64
+ 'LEU','LYS','MET','PHE','PRO','SER','THR','TRP','TYR','VAL','GAP']
65
+
66
+ aa_1_N = {a:n for n,a in enumerate(alpha_1)}
67
+ aa_3_N = {a:n for n,a in enumerate(alpha_3)}
68
+ aa_N_1 = {n:a for n,a in enumerate(alpha_1)}
69
+ aa_1_3 = {a:b for a,b in zip(alpha_1,alpha_3)}
70
+ aa_3_1 = {b:a for a,b in zip(alpha_1,alpha_3)}
71
+
72
+ def AA_to_N(x):
73
+ # ["ARND"] -> [[0,1,2,3]]
74
+ x = np.array(x);
75
+ if x.ndim == 0: x = x[None]
76
+ return [[aa_1_N.get(a, states-1) for a in y] for y in x]
77
+
78
+ def N_to_AA(x):
79
+ # [[0,1,2,3]] -> ["ARND"]
80
+ x = np.array(x);
81
+ if x.ndim == 1: x = x[None]
82
+ return ["".join([aa_N_1.get(a,"-") for a in y]) for y in x]
83
+
84
+ xyz,seq,min_resn,max_resn = {},{},1e6,-1e6
85
+ for line in open(x,"rb"):
86
+ line = line.decode("utf-8","ignore").rstrip()
87
+
88
+ if line[:6] == "HETATM" and line[17:17+3] == "MSE":
89
+ line = line.replace("HETATM","ATOM ")
90
+ line = line.replace("MSE","MET")
91
+
92
+ if line[:4] == "ATOM":
93
+ ch = line[21:22]
94
+ if ch == chain or chain is None:
95
+ atom = line[12:12+4].strip()
96
+ resi = line[17:17+3]
97
+ resn = line[22:22+5].strip()
98
+ x,y,z = [float(line[i:(i+8)]) for i in [30,38,46]]
99
+
100
+ if resn[-1].isalpha():
101
+ resa,resn = resn[-1],int(resn[:-1])-1
102
+ else:
103
+ resa,resn = "",int(resn)-1
104
+ # resn = int(resn)
105
+ if resn < min_resn:
106
+ min_resn = resn
107
+ if resn > max_resn:
108
+ max_resn = resn
109
+ if resn not in xyz:
110
+ xyz[resn] = {}
111
+ if resa not in xyz[resn]:
112
+ xyz[resn][resa] = {}
113
+ if resn not in seq:
114
+ seq[resn] = {}
115
+ if resa not in seq[resn]:
116
+ seq[resn][resa] = resi
117
+
118
+ if atom not in xyz[resn][resa]:
119
+ xyz[resn][resa][atom] = np.array([x,y,z])
120
+
121
+ # convert to numpy arrays, fill in missing values
122
+ seq_,xyz_ = [],[]
123
+ try:
124
+ for resn in range(min_resn,max_resn+1):
125
+ if resn in seq:
126
+ for k in sorted(seq[resn]): seq_.append(aa_3_N.get(seq[resn][k],20))
127
+ else: seq_.append(20)
128
+ if resn in xyz:
129
+ for k in sorted(xyz[resn]):
130
+ for atom in atoms:
131
+ if atom in xyz[resn][k]: xyz_.append(xyz[resn][k][atom])
132
+ else: xyz_.append(np.full(3,np.nan))
133
+ else:
134
+ for atom in atoms: xyz_.append(np.full(3,np.nan))
135
+ return np.array(xyz_).reshape(-1,len(atoms),3), N_to_AA(np.array(seq_))
136
+ except TypeError:
137
+ return 'no_chain', 'no_chain'
138
+
139
+ def parse_PDB(path_to_pdb, input_chain_list=None, ca_only=False):
140
+ c=0
141
+ pdb_dict_list = []
142
+ init_alphabet = ['A', 'B', 'C', 'D', 'E', 'F', 'G','H', 'I', 'J','K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S', 'T','U', 'V','W','X', 'Y', 'Z', 'a', 'b', 'c', 'd', 'e', 'f', 'g','h', 'i', 'j','k', 'l', 'm', 'n', 'o', 'p', 'q', 'r', 's', 't','u', 'v','w','x', 'y', 'z']
143
+ extra_alphabet = [str(item) for item in list(np.arange(300))]
144
+ chain_alphabet = init_alphabet + extra_alphabet
145
+
146
+ if input_chain_list:
147
+ chain_alphabet = input_chain_list
148
+
149
+
150
+ biounit_names = [path_to_pdb]
151
+ for biounit in biounit_names:
152
+ my_dict = {}
153
+ s = 0
154
+ concat_seq = ''
155
+ concat_N = []
156
+ concat_CA = []
157
+ concat_C = []
158
+ concat_O = []
159
+ concat_mask = []
160
+ coords_dict = {}
161
+ for letter in chain_alphabet:
162
+ if ca_only:
163
+ sidechain_atoms = ['CA']
164
+ else:
165
+ sidechain_atoms = ['N', 'CA', 'C', 'O']
166
+ xyz, seq = parse_PDB_biounits(biounit, atoms=sidechain_atoms, chain=letter)
167
+ if type(xyz) != str:
168
+ concat_seq += seq[0]
169
+ my_dict['seq_chain_'+letter]=seq[0]
170
+ coords_dict_chain = {}
171
+ if ca_only:
172
+ coords_dict_chain['CA_chain_'+letter]=xyz.tolist()
173
+ else:
174
+ coords_dict_chain['N_chain_' + letter] = xyz[:, 0, :].tolist()
175
+ coords_dict_chain['CA_chain_' + letter] = xyz[:, 1, :].tolist()
176
+ coords_dict_chain['C_chain_' + letter] = xyz[:, 2, :].tolist()
177
+ coords_dict_chain['O_chain_' + letter] = xyz[:, 3, :].tolist()
178
+ my_dict['coords_chain_'+letter]=coords_dict_chain
179
+ s += 1
180
+ fi = biounit.rfind("/")
181
+ my_dict['name']=biounit[(fi+1):-4]
182
+ my_dict['num_of_chains'] = s
183
+ my_dict['seq'] = concat_seq
184
+ if s <= len(chain_alphabet):
185
+ pdb_dict_list.append(my_dict)
186
+ c+=1
187
+ return pdb_dict_list
188
+
189
+
190
+
191
+ def tied_featurize(batch, device, chain_dict, fixed_position_dict=None, omit_AA_dict=None, tied_positions_dict=None, pssm_dict=None, bias_by_res_dict=None, ca_only=False):
192
+ """ Pack and pad batch into torch tensors """
193
+ alphabet = 'ACDEFGHIKLMNPQRSTVWYX'
194
+ B = len(batch)
195
+ lengths = np.array([len(b['seq']) for b in batch], dtype=np.int32) #sum of chain seq lengths
196
+ L_max = max([len(b['seq']) for b in batch])
197
+ if ca_only:
198
+ X = np.zeros([B, L_max, 1, 3])
199
+ else:
200
+ X = np.zeros([B, L_max, 4, 3])
201
+ residue_idx = -100*np.ones([B, L_max], dtype=np.int32)
202
+ chain_M = np.zeros([B, L_max], dtype=np.int32) #1.0 for the bits that need to be predicted
203
+ pssm_coef_all = np.zeros([B, L_max], dtype=np.float32) #1.0 for the bits that need to be predicted
204
+ pssm_bias_all = np.zeros([B, L_max, 21], dtype=np.float32) #1.0 for the bits that need to be predicted
205
+ pssm_log_odds_all = 10000.0*np.ones([B, L_max, 21], dtype=np.float32) #1.0 for the bits that need to be predicted
206
+ chain_M_pos = np.zeros([B, L_max], dtype=np.int32) #1.0 for the bits that need to be predicted
207
+ bias_by_res_all = np.zeros([B, L_max, 21], dtype=np.float32)
208
+ chain_encoding_all = np.zeros([B, L_max], dtype=np.int32) #1.0 for the bits that need to be predicted
209
+ S = np.zeros([B, L_max], dtype=np.int32)
210
+ omit_AA_mask = np.zeros([B, L_max, len(alphabet)], dtype=np.int32)
211
+ # Build the batch
212
+ letter_list_list = []
213
+ visible_list_list = []
214
+ masked_list_list = []
215
+ masked_chain_length_list_list = []
216
+ tied_pos_list_of_lists_list = []
217
+ for i, b in enumerate(batch):
218
+ if chain_dict != None:
219
+ masked_chains, visible_chains = chain_dict[b['name']] #masked_chains a list of chain letters to predict [A, D, F]
220
+ else:
221
+ masked_chains = [item[-1:] for item in list(b) if item[:10]=='seq_chain_']
222
+ visible_chains = []
223
+ masked_chains.sort() #sort masked_chains
224
+ visible_chains.sort() #sort visible_chains
225
+ all_chains = masked_chains + visible_chains
226
+ for i, b in enumerate(batch):
227
+ mask_dict = {}
228
+ a = 0
229
+ x_chain_list = []
230
+ chain_mask_list = []
231
+ chain_seq_list = []
232
+ chain_encoding_list = []
233
+ c = 1
234
+ letter_list = []
235
+ global_idx_start_list = [0]
236
+ visible_list = []
237
+ masked_list = []
238
+ masked_chain_length_list = []
239
+ fixed_position_mask_list = []
240
+ omit_AA_mask_list = []
241
+ pssm_coef_list = []
242
+ pssm_bias_list = []
243
+ pssm_log_odds_list = []
244
+ bias_by_res_list = []
245
+ l0 = 0
246
+ l1 = 0
247
+ for step, letter in enumerate(all_chains):
248
+ if letter in visible_chains:
249
+ letter_list.append(letter)
250
+ visible_list.append(letter)
251
+ chain_seq = b[f'seq_chain_{letter}']
252
+ chain_seq = ''.join([a if a!='-' else 'X' for a in chain_seq])
253
+ chain_length = len(chain_seq)
254
+ global_idx_start_list.append(global_idx_start_list[-1]+chain_length)
255
+ chain_coords = b[f'coords_chain_{letter}'] #this is a dictionary
256
+ chain_mask = np.zeros(chain_length) #0.0 for visible chains
257
+ if ca_only:
258
+ x_chain = np.array(chain_coords[f'CA_chain_{letter}']) #[chain_lenght,1,3] #CA_diff
259
+ if len(x_chain.shape) == 2:
260
+ x_chain = x_chain[:,None,:]
261
+ else:
262
+ x_chain = np.stack([chain_coords[c] for c in [f'N_chain_{letter}', f'CA_chain_{letter}', f'C_chain_{letter}', f'O_chain_{letter}']], 1) #[chain_lenght,4,3]
263
+ x_chain_list.append(x_chain)
264
+ chain_mask_list.append(chain_mask)
265
+ chain_seq_list.append(chain_seq)
266
+ chain_encoding_list.append(c*np.ones(np.array(chain_mask).shape[0]))
267
+ l1 += chain_length
268
+ residue_idx[i, l0:l1] = 100*(c-1)+np.arange(l0, l1)
269
+ l0 += chain_length
270
+ c+=1
271
+ fixed_position_mask = np.ones(chain_length)
272
+ fixed_position_mask_list.append(fixed_position_mask)
273
+ omit_AA_mask_temp = np.zeros([chain_length, len(alphabet)], np.int32)
274
+ omit_AA_mask_list.append(omit_AA_mask_temp)
275
+ pssm_coef = np.zeros(chain_length)
276
+ pssm_bias = np.zeros([chain_length, 21])
277
+ pssm_log_odds = 10000.0*np.ones([chain_length, 21])
278
+ pssm_coef_list.append(pssm_coef)
279
+ pssm_bias_list.append(pssm_bias)
280
+ pssm_log_odds_list.append(pssm_log_odds)
281
+ bias_by_res_list.append(np.zeros([chain_length, 21]))
282
+ if letter in masked_chains:
283
+ masked_list.append(letter)
284
+ letter_list.append(letter)
285
+ chain_seq = b[f'seq_chain_{letter}']
286
+ chain_seq = ''.join([a if a!='-' else 'X' for a in chain_seq])
287
+ chain_length = len(chain_seq)
288
+ global_idx_start_list.append(global_idx_start_list[-1]+chain_length)
289
+ masked_chain_length_list.append(chain_length)
290
+ chain_coords = b[f'coords_chain_{letter}'] #this is a dictionary
291
+ chain_mask = np.ones(chain_length) #1.0 for masked
292
+ if ca_only:
293
+ x_chain = np.array(chain_coords[f'CA_chain_{letter}']) #[chain_lenght,1,3] #CA_diff
294
+ if len(x_chain.shape) == 2:
295
+ x_chain = x_chain[:,None,:]
296
+ else:
297
+ x_chain = np.stack([chain_coords[c] for c in [f'N_chain_{letter}', f'CA_chain_{letter}', f'C_chain_{letter}', f'O_chain_{letter}']], 1) #[chain_lenght,4,3]
298
+ x_chain_list.append(x_chain)
299
+ chain_mask_list.append(chain_mask)
300
+ chain_seq_list.append(chain_seq)
301
+ chain_encoding_list.append(c*np.ones(np.array(chain_mask).shape[0]))
302
+ l1 += chain_length
303
+ residue_idx[i, l0:l1] = 100*(c-1)+np.arange(l0, l1)
304
+ l0 += chain_length
305
+ c+=1
306
+ fixed_position_mask = np.ones(chain_length)
307
+ if fixed_position_dict!=None:
308
+ fixed_pos_list = fixed_position_dict[b['name']][letter]
309
+ if fixed_pos_list:
310
+ fixed_position_mask[np.array(fixed_pos_list)-1] = 0.0
311
+ fixed_position_mask_list.append(fixed_position_mask)
312
+ omit_AA_mask_temp = np.zeros([chain_length, len(alphabet)], np.int32)
313
+ if omit_AA_dict!=None:
314
+ for item in omit_AA_dict[b['name']][letter]:
315
+ idx_AA = np.array(item[0])-1
316
+ AA_idx = np.array([np.argwhere(np.array(list(alphabet))== AA)[0][0] for AA in item[1]]).repeat(idx_AA.shape[0])
317
+ idx_ = np.array([[a, b] for a in idx_AA for b in AA_idx])
318
+ omit_AA_mask_temp[idx_[:,0], idx_[:,1]] = 1
319
+ omit_AA_mask_list.append(omit_AA_mask_temp)
320
+ pssm_coef = np.zeros(chain_length)
321
+ pssm_bias = np.zeros([chain_length, 21])
322
+ pssm_log_odds = 10000.0*np.ones([chain_length, 21])
323
+ if pssm_dict:
324
+ if pssm_dict[b['name']][letter]:
325
+ pssm_coef = pssm_dict[b['name']][letter]['pssm_coef']
326
+ pssm_bias = pssm_dict[b['name']][letter]['pssm_bias']
327
+ pssm_log_odds = pssm_dict[b['name']][letter]['pssm_log_odds']
328
+ pssm_coef_list.append(pssm_coef)
329
+ pssm_bias_list.append(pssm_bias)
330
+ pssm_log_odds_list.append(pssm_log_odds)
331
+ if bias_by_res_dict:
332
+ bias_by_res_list.append(bias_by_res_dict[b['name']][letter])
333
+ else:
334
+ bias_by_res_list.append(np.zeros([chain_length, 21]))
335
+
336
+
337
+ letter_list_np = np.array(letter_list)
338
+ tied_pos_list_of_lists = []
339
+ tied_beta = np.ones(L_max)
340
+ if tied_positions_dict!=None:
341
+ tied_pos_list = tied_positions_dict[b['name']]
342
+ if tied_pos_list:
343
+ set_chains_tied = set(list(itertools.chain(*[list(item) for item in tied_pos_list])))
344
+ for tied_item in tied_pos_list:
345
+ one_list = []
346
+ for k, v in tied_item.items():
347
+ start_idx = global_idx_start_list[np.argwhere(letter_list_np == k)[0][0]]
348
+ if isinstance(v[0], list):
349
+ for v_count in range(len(v[0])):
350
+ one_list.append(start_idx+v[0][v_count]-1)#make 0 to be the first
351
+ tied_beta[start_idx+v[0][v_count]-1] = v[1][v_count]
352
+ else:
353
+ for v_ in v:
354
+ one_list.append(start_idx+v_-1)#make 0 to be the first
355
+ tied_pos_list_of_lists.append(one_list)
356
+ tied_pos_list_of_lists_list.append(tied_pos_list_of_lists)
357
+
358
+
359
+
360
+ x = np.concatenate(x_chain_list,0) #[L, 4, 3]
361
+ all_sequence = "".join(chain_seq_list)
362
+ m = np.concatenate(chain_mask_list,0) #[L,], 1.0 for places that need to be predicted
363
+ chain_encoding = np.concatenate(chain_encoding_list,0)
364
+ m_pos = np.concatenate(fixed_position_mask_list,0) #[L,], 1.0 for places that need to be predicted
365
+
366
+ pssm_coef_ = np.concatenate(pssm_coef_list,0) #[L,], 1.0 for places that need to be predicted
367
+ pssm_bias_ = np.concatenate(pssm_bias_list,0) #[L,], 1.0 for places that need to be predicted
368
+ pssm_log_odds_ = np.concatenate(pssm_log_odds_list,0) #[L,], 1.0 for places that need to be predicted
369
+
370
+ bias_by_res_ = np.concatenate(bias_by_res_list, 0) #[L,21], 0.0 for places where AA frequencies don't need to be tweaked
371
+
372
+ l = len(all_sequence)
373
+ x_pad = np.pad(x, [[0,L_max-l], [0,0], [0,0]], 'constant', constant_values=(np.nan, ))
374
+ X[i,:,:,:] = x_pad
375
+
376
+ m_pad = np.pad(m, [[0,L_max-l]], 'constant', constant_values=(0.0, ))
377
+ m_pos_pad = np.pad(m_pos, [[0,L_max-l]], 'constant', constant_values=(0.0, ))
378
+ omit_AA_mask_pad = np.pad(np.concatenate(omit_AA_mask_list,0), [[0,L_max-l]], 'constant', constant_values=(0.0, ))
379
+ chain_M[i,:] = m_pad
380
+ chain_M_pos[i,:] = m_pos_pad
381
+ omit_AA_mask[i,] = omit_AA_mask_pad
382
+
383
+ chain_encoding_pad = np.pad(chain_encoding, [[0,L_max-l]], 'constant', constant_values=(0.0, ))
384
+ chain_encoding_all[i,:] = chain_encoding_pad
385
+
386
+ pssm_coef_pad = np.pad(pssm_coef_, [[0,L_max-l]], 'constant', constant_values=(0.0, ))
387
+ pssm_bias_pad = np.pad(pssm_bias_, [[0,L_max-l], [0,0]], 'constant', constant_values=(0.0, ))
388
+ pssm_log_odds_pad = np.pad(pssm_log_odds_, [[0,L_max-l], [0,0]], 'constant', constant_values=(0.0, ))
389
+
390
+ pssm_coef_all[i,:] = pssm_coef_pad
391
+ pssm_bias_all[i,:] = pssm_bias_pad
392
+ pssm_log_odds_all[i,:] = pssm_log_odds_pad
393
+
394
+ bias_by_res_pad = np.pad(bias_by_res_, [[0,L_max-l], [0,0]], 'constant', constant_values=(0.0, ))
395
+ bias_by_res_all[i,:] = bias_by_res_pad
396
+
397
+ # Convert to labels
398
+ indices = np.asarray([alphabet.index(a) for a in all_sequence], dtype=np.int32)
399
+ S[i, :l] = indices
400
+ letter_list_list.append(letter_list)
401
+ visible_list_list.append(visible_list)
402
+ masked_list_list.append(masked_list)
403
+ masked_chain_length_list_list.append(masked_chain_length_list)
404
+
405
+
406
+ isnan = np.isnan(X)
407
+ mask = np.isfinite(np.sum(X,(2,3))).astype(np.float32)
408
+ X[isnan] = 0.
409
+
410
+ # Conversion
411
+ pssm_coef_all = torch.from_numpy(pssm_coef_all).to(dtype=torch.float32, device=device)
412
+ pssm_bias_all = torch.from_numpy(pssm_bias_all).to(dtype=torch.float32, device=device)
413
+ pssm_log_odds_all = torch.from_numpy(pssm_log_odds_all).to(dtype=torch.float32, device=device)
414
+
415
+ tied_beta = torch.from_numpy(tied_beta).to(dtype=torch.float32, device=device)
416
+
417
+ jumps = ((residue_idx[:,1:]-residue_idx[:,:-1])==1).astype(np.float32)
418
+ bias_by_res_all = torch.from_numpy(bias_by_res_all).to(dtype=torch.float32, device=device)
419
+ phi_mask = np.pad(jumps, [[0,0],[1,0]])
420
+ psi_mask = np.pad(jumps, [[0,0],[0,1]])
421
+ omega_mask = np.pad(jumps, [[0,0],[0,1]])
422
+ dihedral_mask = np.concatenate([phi_mask[:,:,None], psi_mask[:,:,None], omega_mask[:,:,None]], -1) #[B,L,3]
423
+ dihedral_mask = torch.from_numpy(dihedral_mask).to(dtype=torch.float32, device=device)
424
+ residue_idx = torch.from_numpy(residue_idx).to(dtype=torch.long,device=device)
425
+ S = torch.from_numpy(S).to(dtype=torch.long,device=device)
426
+ X = torch.from_numpy(X).to(dtype=torch.float32, device=device)
427
+ mask = torch.from_numpy(mask).to(dtype=torch.float32, device=device)
428
+ chain_M = torch.from_numpy(chain_M).to(dtype=torch.float32, device=device)
429
+ chain_M_pos = torch.from_numpy(chain_M_pos).to(dtype=torch.float32, device=device)
430
+ omit_AA_mask = torch.from_numpy(omit_AA_mask).to(dtype=torch.float32, device=device)
431
+ chain_encoding_all = torch.from_numpy(chain_encoding_all).to(dtype=torch.long, device=device)
432
+ if ca_only:
433
+ X_out = X[:,:,0]
434
+ else:
435
+ X_out = X
436
+ return X_out, S, mask, lengths, chain_M, chain_encoding_all, letter_list_list, visible_list_list, masked_list_list, masked_chain_length_list_list, chain_M_pos, omit_AA_mask, residue_idx, dihedral_mask, tied_pos_list_of_lists_list, pssm_coef_all, pssm_bias_all, pssm_log_odds_all, bias_by_res_all, tied_beta
437
+
438
+
439
+
440
+ def loss_nll(S, log_probs, mask):
441
+ """ Negative log probabilities """
442
+ criterion = torch.nn.NLLLoss(reduction='none')
443
+ loss = criterion(
444
+ log_probs.contiguous().view(-1, log_probs.size(-1)), S.contiguous().view(-1)
445
+ ).view(S.size())
446
+ loss_av = torch.sum(loss * mask) / torch.sum(mask)
447
+ return loss, loss_av
448
+
449
+
450
+ def loss_smoothed(S, log_probs, mask, weight=0.1):
451
+ """ Negative log probabilities """
452
+ S_onehot = torch.nn.functional.one_hot(S, 21).float()
453
+
454
+ # Label smoothing
455
+ S_onehot = S_onehot + weight / float(S_onehot.size(-1))
456
+ S_onehot = S_onehot / S_onehot.sum(-1, keepdim=True)
457
+
458
+ loss = -(S_onehot * log_probs).sum(-1)
459
+ loss_av = torch.sum(loss * mask) / torch.sum(mask)
460
+ return loss, loss_av
461
+
462
+ class StructureDataset():
463
+ def __init__(self, jsonl_file, verbose=True, truncate=None, max_length=100,
464
+ alphabet='ACDEFGHIKLMNPQRSTVWYX-'):
465
+ alphabet_set = set([a for a in alphabet])
466
+ discard_count = {
467
+ 'bad_chars': 0,
468
+ 'too_long': 0,
469
+ 'bad_seq_length': 0
470
+ }
471
+
472
+ with open(jsonl_file) as f:
473
+ self.data = []
474
+
475
+ lines = f.readlines()
476
+ start = time.time()
477
+ for i, line in enumerate(lines):
478
+ entry = json.loads(line)
479
+ seq = entry['seq']
480
+ name = entry['name']
481
+
482
+ # Convert raw coords to np arrays
483
+ #for key, val in entry['coords'].items():
484
+ # entry['coords'][key] = np.asarray(val)
485
+
486
+ # Check if in alphabet
487
+ bad_chars = set([s for s in seq]).difference(alphabet_set)
488
+ if len(bad_chars) == 0:
489
+ if len(entry['seq']) <= max_length:
490
+ if True:
491
+ self.data.append(entry)
492
+ else:
493
+ discard_count['bad_seq_length'] += 1
494
+ else:
495
+ discard_count['too_long'] += 1
496
+ else:
497
+ if verbose:
498
+ print(name, bad_chars, entry['seq'])
499
+ discard_count['bad_chars'] += 1
500
+
501
+ # Truncate early
502
+ if truncate is not None and len(self.data) == truncate:
503
+ return
504
+
505
+ if verbose and (i + 1) % 1000 == 0:
506
+ elapsed = time.time() - start
507
+ print('{} entries ({} loaded) in {:.1f} s'.format(len(self.data), i+1, elapsed))
508
+ if verbose:
509
+ print('discarded', discard_count)
510
+ def __len__(self):
511
+ return len(self.data)
512
+
513
+ def __getitem__(self, idx):
514
+ return self.data[idx]
515
+
516
+
517
+ class StructureDatasetPDB():
518
+ def __init__(self, pdb_dict_list, verbose=True, truncate=None, max_length=100,
519
+ alphabet='ACDEFGHIKLMNPQRSTVWYX-'):
520
+ alphabet_set = set([a for a in alphabet])
521
+ discard_count = {
522
+ 'bad_chars': 0,
523
+ 'too_long': 0,
524
+ 'bad_seq_length': 0
525
+ }
526
+
527
+ self.data = []
528
+
529
+ start = time.time()
530
+ for i, entry in enumerate(pdb_dict_list):
531
+ seq = entry['seq']
532
+ name = entry['name']
533
+
534
+ bad_chars = set([s for s in seq]).difference(alphabet_set)
535
+ if len(bad_chars) == 0:
536
+ if len(entry['seq']) <= max_length:
537
+ self.data.append(entry)
538
+ else:
539
+ discard_count['too_long'] += 1
540
+ else:
541
+ discard_count['bad_chars'] += 1
542
+
543
+ # Truncate early
544
+ if truncate is not None and len(self.data) == truncate:
545
+ return
546
+
547
+ if verbose and (i + 1) % 1000 == 0:
548
+ elapsed = time.time() - start
549
+
550
+ #print('Discarded', discard_count)
551
+ def __len__(self):
552
+ return len(self.data)
553
+
554
+ def __getitem__(self, idx):
555
+ return self.data[idx]
556
+
557
+
558
+
559
+ class StructureLoader():
560
+ def __init__(self, dataset, batch_size=100, shuffle=True,
561
+ collate_fn=lambda x:x, drop_last=False):
562
+ self.dataset = dataset
563
+ self.size = len(dataset)
564
+ self.lengths = [len(dataset[i]['seq']) for i in range(self.size)]
565
+ self.batch_size = batch_size
566
+ sorted_ix = np.argsort(self.lengths)
567
+
568
+ # Cluster into batches of similar sizes
569
+ clusters, batch = [], []
570
+ batch_max = 0
571
+ for ix in sorted_ix:
572
+ size = self.lengths[ix]
573
+ if size * (len(batch) + 1) <= self.batch_size:
574
+ batch.append(ix)
575
+ batch_max = size
576
+ else:
577
+ clusters.append(batch)
578
+ batch, batch_max = [], 0
579
+ if len(batch) > 0:
580
+ clusters.append(batch)
581
+ self.clusters = clusters
582
+
583
+ def __len__(self):
584
+ return len(self.clusters)
585
+
586
+ def __iter__(self):
587
+ np.random.shuffle(self.clusters)
588
+ for b_idx in self.clusters:
589
+ batch = [self.dataset[i] for i in b_idx]
590
+ yield batch
591
+
592
+
593
+
594
+ # The following gather functions
595
+ def gather_edges(edges, neighbor_idx):
596
+ # Features [B,N,N,C] at Neighbor indices [B,N,K] => Neighbor features [B,N,K,C]
597
+ neighbors = neighbor_idx.unsqueeze(-1).expand(-1, -1, -1, edges.size(-1))
598
+ edge_features = torch.gather(edges, 2, neighbors)
599
+ return edge_features
600
+
601
+ def gather_nodes(nodes, neighbor_idx):
602
+ # Features [B,N,C] at Neighbor indices [B,N,K] => [B,N,K,C]
603
+ # Flatten and expand indices per batch [B,N,K] => [B,NK] => [B,NK,C]
604
+ neighbors_flat = neighbor_idx.view((neighbor_idx.shape[0], -1))
605
+ neighbors_flat = neighbors_flat.unsqueeze(-1).expand(-1, -1, nodes.size(2))
606
+ # Gather and re-pack
607
+ neighbor_features = torch.gather(nodes, 1, neighbors_flat)
608
+ neighbor_features = neighbor_features.view(list(neighbor_idx.shape)[:3] + [-1])
609
+ return neighbor_features
610
+
611
+ def gather_nodes_t(nodes, neighbor_idx):
612
+ # Features [B,N,C] at Neighbor index [B,K] => Neighbor features[B,K,C]
613
+ idx_flat = neighbor_idx.unsqueeze(-1).expand(-1, -1, nodes.size(2))
614
+ neighbor_features = torch.gather(nodes, 1, idx_flat)
615
+ return neighbor_features
616
+
617
+ def cat_neighbors_nodes(h_nodes, h_neighbors, E_idx):
618
+ h_nodes = gather_nodes(h_nodes, E_idx)
619
+ h_nn = torch.cat([h_neighbors, h_nodes], -1)
620
+ return h_nn
621
+
622
+
623
+ class EncLayer(nn.Module):
624
+ def __init__(self, num_hidden, num_in, dropout=0.1, num_heads=None, scale=30):
625
+ super(EncLayer, self).__init__()
626
+ self.num_hidden = num_hidden
627
+ self.num_in = num_in
628
+ self.scale = scale
629
+ self.dropout1 = nn.Dropout(dropout)
630
+ self.dropout2 = nn.Dropout(dropout)
631
+ self.dropout3 = nn.Dropout(dropout)
632
+ self.norm1 = nn.LayerNorm(num_hidden)
633
+ self.norm2 = nn.LayerNorm(num_hidden)
634
+ self.norm3 = nn.LayerNorm(num_hidden)
635
+
636
+ self.W1 = nn.Linear(num_hidden + num_in, num_hidden, bias=True)
637
+ self.W2 = nn.Linear(num_hidden, num_hidden, bias=True)
638
+ self.W3 = nn.Linear(num_hidden, num_hidden, bias=True)
639
+ self.W11 = nn.Linear(num_hidden + num_in, num_hidden, bias=True)
640
+ self.W12 = nn.Linear(num_hidden, num_hidden, bias=True)
641
+ self.W13 = nn.Linear(num_hidden, num_hidden, bias=True)
642
+ self.act = torch.nn.GELU()
643
+ self.dense = PositionWiseFeedForward(num_hidden, num_hidden * 4)
644
+
645
+ def forward(self, h_V, h_E, E_idx, mask_V=None, mask_attend=None):
646
+ """ Parallel computation of full transformer layer """
647
+
648
+ h_EV = cat_neighbors_nodes(h_V, h_E, E_idx)
649
+ h_V_expand = h_V.unsqueeze(-2).expand(-1,-1,h_EV.size(-2),-1)
650
+ h_EV = torch.cat([h_V_expand, h_EV], -1)
651
+ h_message = self.W3(self.act(self.W2(self.act(self.W1(h_EV)))))
652
+ if mask_attend is not None:
653
+ h_message = mask_attend.unsqueeze(-1) * h_message
654
+ dh = torch.sum(h_message, -2) / self.scale
655
+ h_V = self.norm1(h_V + self.dropout1(dh))
656
+
657
+ dh = self.dense(h_V)
658
+ h_V = self.norm2(h_V + self.dropout2(dh))
659
+ if mask_V is not None:
660
+ mask_V = mask_V.unsqueeze(-1)
661
+ h_V = mask_V * h_V
662
+
663
+ h_EV = cat_neighbors_nodes(h_V, h_E, E_idx)
664
+ h_V_expand = h_V.unsqueeze(-2).expand(-1,-1,h_EV.size(-2),-1)
665
+ h_EV = torch.cat([h_V_expand, h_EV], -1)
666
+ h_message = self.W13(self.act(self.W12(self.act(self.W11(h_EV)))))
667
+ h_E = self.norm3(h_E + self.dropout3(h_message))
668
+ return h_V, h_E
669
+
670
+
671
+ class DecLayer(nn.Module):
672
+ def __init__(self, num_hidden, num_in, dropout=0.1, num_heads=None, scale=30):
673
+ super(DecLayer, self).__init__()
674
+ self.num_hidden = num_hidden
675
+ self.num_in = num_in
676
+ self.scale = scale
677
+ self.dropout1 = nn.Dropout(dropout)
678
+ self.dropout2 = nn.Dropout(dropout)
679
+ self.norm1 = nn.LayerNorm(num_hidden)
680
+ self.norm2 = nn.LayerNorm(num_hidden)
681
+
682
+ self.W1 = nn.Linear(num_hidden + num_in, num_hidden, bias=True)
683
+ self.W2 = nn.Linear(num_hidden, num_hidden, bias=True)
684
+ self.W3 = nn.Linear(num_hidden, num_hidden, bias=True)
685
+ self.act = torch.nn.GELU()
686
+ self.dense = PositionWiseFeedForward(num_hidden, num_hidden * 4)
687
+
688
+ def forward(self, h_V, h_E, mask_V=None, mask_attend=None):
689
+ """ Parallel computation of full transformer layer """
690
+
691
+ # Concatenate h_V_i to h_E_ij
692
+ h_V_expand = h_V.unsqueeze(-2).expand(-1,-1,h_E.size(-2),-1)
693
+ h_EV = torch.cat([h_V_expand, h_E], -1)
694
+
695
+ h_message = self.W3(self.act(self.W2(self.act(self.W1(h_EV)))))
696
+ if mask_attend is not None:
697
+ h_message = mask_attend.unsqueeze(-1) * h_message
698
+ dh = torch.sum(h_message, -2) / self.scale
699
+
700
+ h_V = self.norm1(h_V + self.dropout1(dh))
701
+
702
+ # Position-wise feedforward
703
+ dh = self.dense(h_V)
704
+ h_V = self.norm2(h_V + self.dropout2(dh))
705
+
706
+ if mask_V is not None:
707
+ mask_V = mask_V.unsqueeze(-1)
708
+ h_V = mask_V * h_V
709
+ return h_V
710
+
711
+
712
+
713
+ class PositionWiseFeedForward(nn.Module):
714
+ def __init__(self, num_hidden, num_ff):
715
+ super(PositionWiseFeedForward, self).__init__()
716
+ self.W_in = nn.Linear(num_hidden, num_ff, bias=True)
717
+ self.W_out = nn.Linear(num_ff, num_hidden, bias=True)
718
+ self.act = torch.nn.GELU()
719
+ def forward(self, h_V):
720
+ h = self.act(self.W_in(h_V))
721
+ h = self.W_out(h)
722
+ return h
723
+
724
+ class PositionalEncodings(nn.Module):
725
+ def __init__(self, num_embeddings, max_relative_feature=32):
726
+ super(PositionalEncodings, self).__init__()
727
+ self.num_embeddings = num_embeddings
728
+ self.max_relative_feature = max_relative_feature
729
+ self.linear = nn.Linear(2*max_relative_feature+1+1, num_embeddings)
730
+
731
+ def forward(self, offset, mask):
732
+ d = torch.clip(offset + self.max_relative_feature, 0, 2*self.max_relative_feature)*mask + (1-mask)*(2*self.max_relative_feature+1)
733
+ d_onehot = torch.nn.functional.one_hot(d, 2*self.max_relative_feature+1+1)
734
+ E = self.linear(d_onehot.float())
735
+ return E
736
+
737
+
738
+
739
+ class CA_ProteinFeatures(nn.Module):
740
+ def __init__(self, edge_features, node_features, num_positional_embeddings=16,
741
+ num_rbf=16, top_k=30, augment_eps=0., num_chain_embeddings=16):
742
+ """ Extract protein features """
743
+ super(CA_ProteinFeatures, self).__init__()
744
+ self.edge_features = edge_features
745
+ self.node_features = node_features
746
+ self.top_k = top_k
747
+ self.augment_eps = augment_eps
748
+ self.num_rbf = num_rbf
749
+ self.num_positional_embeddings = num_positional_embeddings
750
+
751
+ # Positional encoding
752
+ self.embeddings = PositionalEncodings(num_positional_embeddings)
753
+ # Normalization and embedding
754
+ node_in, edge_in = 3, num_positional_embeddings + num_rbf*9 + 7
755
+ self.node_embedding = nn.Linear(node_in, node_features, bias=False) #NOT USED
756
+ self.edge_embedding = nn.Linear(edge_in, edge_features, bias=False)
757
+ self.norm_nodes = nn.LayerNorm(node_features)
758
+ self.norm_edges = nn.LayerNorm(edge_features)
759
+
760
+
761
+ def _quaternions(self, R):
762
+ """ Convert a batch of 3D rotations [R] to quaternions [Q]
763
+ R [...,3,3]
764
+ Q [...,4]
765
+ """
766
+ # Simple Wikipedia version
767
+ # en.wikipedia.org/wiki/Rotation_matrix#Quaternion
768
+ # For other options see math.stackexchange.com/questions/2074316/calculating-rotation-axis-from-rotation-matrix
769
+ diag = torch.diagonal(R, dim1=-2, dim2=-1)
770
+ Rxx, Ryy, Rzz = diag.unbind(-1)
771
+ magnitudes = 0.5 * torch.sqrt(torch.abs(1 + torch.stack([
772
+ Rxx - Ryy - Rzz,
773
+ - Rxx + Ryy - Rzz,
774
+ - Rxx - Ryy + Rzz
775
+ ], -1)))
776
+ _R = lambda i,j: R[:,:,:,i,j]
777
+ signs = torch.sign(torch.stack([
778
+ _R(2,1) - _R(1,2),
779
+ _R(0,2) - _R(2,0),
780
+ _R(1,0) - _R(0,1)
781
+ ], -1))
782
+ xyz = signs * magnitudes
783
+ # The relu enforces a non-negative trace
784
+ w = torch.sqrt(F.relu(1 + diag.sum(-1, keepdim=True))) / 2.
785
+ Q = torch.cat((xyz, w), -1)
786
+ Q = F.normalize(Q, dim=-1)
787
+ return Q
788
+
789
+ def _orientations_coarse(self, X, E_idx, eps=1e-6):
790
+ dX = X[:,1:,:] - X[:,:-1,:]
791
+ dX_norm = torch.norm(dX,dim=-1)
792
+ dX_mask = (3.6<dX_norm) & (dX_norm<4.0) #exclude CA-CA jumps
793
+ dX = dX*dX_mask[:,:,None]
794
+ U = F.normalize(dX, dim=-1)
795
+ u_2 = U[:,:-2,:]
796
+ u_1 = U[:,1:-1,:]
797
+ u_0 = U[:,2:,:]
798
+ # Backbone normals
799
+ n_2 = F.normalize(torch.cross(u_2, u_1), dim=-1)
800
+ n_1 = F.normalize(torch.cross(u_1, u_0), dim=-1)
801
+
802
+ # Bond angle calculation
803
+ cosA = -(u_1 * u_0).sum(-1)
804
+ cosA = torch.clamp(cosA, -1+eps, 1-eps)
805
+ A = torch.acos(cosA)
806
+ # Angle between normals
807
+ cosD = (n_2 * n_1).sum(-1)
808
+ cosD = torch.clamp(cosD, -1+eps, 1-eps)
809
+ D = torch.sign((u_2 * n_1).sum(-1)) * torch.acos(cosD)
810
+ # Backbone features
811
+ AD_features = torch.stack((torch.cos(A), torch.sin(A) * torch.cos(D), torch.sin(A) * torch.sin(D)), 2)
812
+ AD_features = F.pad(AD_features, (0,0,1,2), 'constant', 0)
813
+
814
+ # Build relative orientations
815
+ o_1 = F.normalize(u_2 - u_1, dim=-1)
816
+ O = torch.stack((o_1, n_2, torch.cross(o_1, n_2)), 2)
817
+ O = O.view(list(O.shape[:2]) + [9])
818
+ O = F.pad(O, (0,0,1,2), 'constant', 0)
819
+ O_neighbors = gather_nodes(O, E_idx)
820
+ X_neighbors = gather_nodes(X, E_idx)
821
+
822
+ # Re-view as rotation matrices
823
+ O = O.view(list(O.shape[:2]) + [3,3])
824
+ O_neighbors = O_neighbors.view(list(O_neighbors.shape[:3]) + [3,3])
825
+
826
+ # Rotate into local reference frames
827
+ dX = X_neighbors - X.unsqueeze(-2)
828
+ dU = torch.matmul(O.unsqueeze(2), dX.unsqueeze(-1)).squeeze(-1)
829
+ dU = F.normalize(dU, dim=-1)
830
+ R = torch.matmul(O.unsqueeze(2).transpose(-1,-2), O_neighbors)
831
+ Q = self._quaternions(R)
832
+
833
+ # Orientation features
834
+ O_features = torch.cat((dU,Q), dim=-1)
835
+ return AD_features, O_features
836
+
837
+
838
+
839
+ def _dist(self, X, mask, eps=1E-6):
840
+ """ Pairwise euclidean distances """
841
+ # Convolutional network on NCHW
842
+ mask_2D = torch.unsqueeze(mask,1) * torch.unsqueeze(mask,2)
843
+ dX = torch.unsqueeze(X,1) - torch.unsqueeze(X,2)
844
+ D = mask_2D * torch.sqrt(torch.sum(dX**2, 3) + eps)
845
+
846
+ # Identify k nearest neighbors (including self)
847
+ D_max, _ = torch.max(D, -1, keepdim=True)
848
+ D_adjust = D + (1. - mask_2D) * D_max
849
+ D_neighbors, E_idx = torch.topk(D_adjust, np.minimum(self.top_k, X.shape[1]), dim=-1, largest=False)
850
+ mask_neighbors = gather_edges(mask_2D.unsqueeze(-1), E_idx)
851
+ return D_neighbors, E_idx, mask_neighbors
852
+
853
+ def _rbf(self, D):
854
+ # Distance radial basis function
855
+ device = D.device
856
+ D_min, D_max, D_count = 2., 22., self.num_rbf
857
+ D_mu = torch.linspace(D_min, D_max, D_count).to(device)
858
+ D_mu = D_mu.view([1,1,1,-1])
859
+ D_sigma = (D_max - D_min) / D_count
860
+ D_expand = torch.unsqueeze(D, -1)
861
+ RBF = torch.exp(-((D_expand - D_mu) / D_sigma)**2)
862
+ return RBF
863
+
864
+ def _get_rbf(self, A, B, E_idx):
865
+ D_A_B = torch.sqrt(torch.sum((A[:,:,None,:] - B[:,None,:,:])**2,-1) + 1e-6) #[B, L, L]
866
+ D_A_B_neighbors = gather_edges(D_A_B[:,:,:,None], E_idx)[:,:,:,0] #[B,L,K]
867
+ RBF_A_B = self._rbf(D_A_B_neighbors)
868
+ return RBF_A_B
869
+
870
+ def forward(self, Ca, mask, residue_idx, chain_labels):
871
+ """ Featurize coordinates as an attributed graph """
872
+ if self.augment_eps > 0:
873
+ Ca = Ca + self.augment_eps * torch.randn_like(Ca)
874
+
875
+ D_neighbors, E_idx, mask_neighbors = self._dist(Ca, mask)
876
+
877
+ Ca_0 = torch.zeros(Ca.shape, device=Ca.device)
878
+ Ca_2 = torch.zeros(Ca.shape, device=Ca.device)
879
+ Ca_0[:,1:,:] = Ca[:,:-1,:]
880
+ Ca_1 = Ca
881
+ Ca_2[:,:-1,:] = Ca[:,1:,:]
882
+
883
+ V, O_features = self._orientations_coarse(Ca, E_idx)
884
+
885
+ RBF_all = []
886
+ RBF_all.append(self._rbf(D_neighbors)) #Ca_1-Ca_1
887
+ RBF_all.append(self._get_rbf(Ca_0, Ca_0, E_idx))
888
+ RBF_all.append(self._get_rbf(Ca_2, Ca_2, E_idx))
889
+
890
+ RBF_all.append(self._get_rbf(Ca_0, Ca_1, E_idx))
891
+ RBF_all.append(self._get_rbf(Ca_0, Ca_2, E_idx))
892
+
893
+ RBF_all.append(self._get_rbf(Ca_1, Ca_0, E_idx))
894
+ RBF_all.append(self._get_rbf(Ca_1, Ca_2, E_idx))
895
+
896
+ RBF_all.append(self._get_rbf(Ca_2, Ca_0, E_idx))
897
+ RBF_all.append(self._get_rbf(Ca_2, Ca_1, E_idx))
898
+
899
+
900
+ RBF_all = torch.cat(tuple(RBF_all), dim=-1)
901
+
902
+
903
+ offset = residue_idx[:,:,None]-residue_idx[:,None,:]
904
+ offset = gather_edges(offset[:,:,:,None], E_idx)[:,:,:,0] #[B, L, K]
905
+
906
+ d_chains = ((chain_labels[:, :, None] - chain_labels[:,None,:])==0).long()
907
+ E_chains = gather_edges(d_chains[:,:,:,None], E_idx)[:,:,:,0]
908
+ E_positional = self.embeddings(offset.long(), E_chains)
909
+ E = torch.cat((E_positional, RBF_all, O_features), -1)
910
+
911
+
912
+ E = self.edge_embedding(E)
913
+ E = self.norm_edges(E)
914
+
915
+ return E, E_idx
916
+
917
+
918
+
919
+
920
+ class ProteinFeatures(nn.Module):
921
+ def __init__(self, edge_features, node_features, num_positional_embeddings=16,
922
+ num_rbf=16, top_k=30, augment_eps=0., num_chain_embeddings=16):
923
+ """ Extract protein features """
924
+ super(ProteinFeatures, self).__init__()
925
+ self.edge_features = edge_features
926
+ self.node_features = node_features
927
+ self.top_k = top_k
928
+ self.augment_eps = augment_eps
929
+ self.num_rbf = num_rbf
930
+ self.num_positional_embeddings = num_positional_embeddings
931
+
932
+ self.embeddings = PositionalEncodings(num_positional_embeddings)
933
+ node_in, edge_in = 6, num_positional_embeddings + num_rbf*25
934
+ self.edge_embedding = nn.Linear(edge_in, edge_features, bias=False)
935
+ self.norm_edges = nn.LayerNorm(edge_features)
936
+
937
+ def _dist(self, X, mask, eps=1E-6):
938
+ mask_2D = torch.unsqueeze(mask,1) * torch.unsqueeze(mask,2)
939
+ dX = torch.unsqueeze(X,1) - torch.unsqueeze(X,2)
940
+ D = mask_2D * torch.sqrt(torch.sum(dX**2, 3) + eps)
941
+ D_max, _ = torch.max(D, -1, keepdim=True)
942
+ D_adjust = D + (1. - mask_2D) * D_max
943
+ sampled_top_k = self.top_k
944
+ D_neighbors, E_idx = torch.topk(D_adjust, np.minimum(self.top_k, X.shape[1]), dim=-1, largest=False)
945
+ return D_neighbors, E_idx
946
+
947
+ def _rbf(self, D):
948
+ device = D.device
949
+ D_min, D_max, D_count = 2., 22., self.num_rbf
950
+ D_mu = torch.linspace(D_min, D_max, D_count, device=device)
951
+ D_mu = D_mu.view([1,1,1,-1])
952
+ D_sigma = (D_max - D_min) / D_count
953
+ D_expand = torch.unsqueeze(D, -1)
954
+ RBF = torch.exp(-((D_expand - D_mu) / D_sigma)**2)
955
+ return RBF
956
+
957
+ def _get_rbf(self, A, B, E_idx):
958
+ D_A_B = torch.sqrt(torch.sum((A[:,:,None,:] - B[:,None,:,:])**2,-1) + 1e-6) #[B, L, L]
959
+ D_A_B_neighbors = gather_edges(D_A_B[:,:,:,None], E_idx)[:,:,:,0] #[B,L,K]
960
+ RBF_A_B = self._rbf(D_A_B_neighbors)
961
+ return RBF_A_B
962
+
963
+ def forward(self, X, mask, residue_idx, chain_labels):
964
+ if self.augment_eps > 0:
965
+ X = X + self.augment_eps * torch.randn_like(X)
966
+
967
+ b = X[:,:,1,:] - X[:,:,0,:]
968
+ c = X[:,:,2,:] - X[:,:,1,:]
969
+ a = torch.cross(b, c, dim=-1)
970
+ Cb = -0.58273431*a + 0.56802827*b - 0.54067466*c + X[:,:,1,:]
971
+ Ca = X[:,:,1,:]
972
+ N = X[:,:,0,:]
973
+ C = X[:,:,2,:]
974
+ O = X[:,:,3,:]
975
+
976
+ D_neighbors, E_idx = self._dist(Ca, mask)
977
+
978
+ RBF_all = []
979
+ RBF_all.append(self._rbf(D_neighbors)) #Ca-Ca
980
+ RBF_all.append(self._get_rbf(N, N, E_idx)) #N-N
981
+ RBF_all.append(self._get_rbf(C, C, E_idx)) #C-C
982
+ RBF_all.append(self._get_rbf(O, O, E_idx)) #O-O
983
+ RBF_all.append(self._get_rbf(Cb, Cb, E_idx)) #Cb-Cb
984
+ RBF_all.append(self._get_rbf(Ca, N, E_idx)) #Ca-N
985
+ RBF_all.append(self._get_rbf(Ca, C, E_idx)) #Ca-C
986
+ RBF_all.append(self._get_rbf(Ca, O, E_idx)) #Ca-O
987
+ RBF_all.append(self._get_rbf(Ca, Cb, E_idx)) #Ca-Cb
988
+ RBF_all.append(self._get_rbf(N, C, E_idx)) #N-C
989
+ RBF_all.append(self._get_rbf(N, O, E_idx)) #N-O
990
+ RBF_all.append(self._get_rbf(N, Cb, E_idx)) #N-Cb
991
+ RBF_all.append(self._get_rbf(Cb, C, E_idx)) #Cb-C
992
+ RBF_all.append(self._get_rbf(Cb, O, E_idx)) #Cb-O
993
+ RBF_all.append(self._get_rbf(O, C, E_idx)) #O-C
994
+ RBF_all.append(self._get_rbf(N, Ca, E_idx)) #N-Ca
995
+ RBF_all.append(self._get_rbf(C, Ca, E_idx)) #C-Ca
996
+ RBF_all.append(self._get_rbf(O, Ca, E_idx)) #O-Ca
997
+ RBF_all.append(self._get_rbf(Cb, Ca, E_idx)) #Cb-Ca
998
+ RBF_all.append(self._get_rbf(C, N, E_idx)) #C-N
999
+ RBF_all.append(self._get_rbf(O, N, E_idx)) #O-N
1000
+ RBF_all.append(self._get_rbf(Cb, N, E_idx)) #Cb-N
1001
+ RBF_all.append(self._get_rbf(C, Cb, E_idx)) #C-Cb
1002
+ RBF_all.append(self._get_rbf(O, Cb, E_idx)) #O-Cb
1003
+ RBF_all.append(self._get_rbf(C, O, E_idx)) #C-O
1004
+ RBF_all = torch.cat(tuple(RBF_all), dim=-1)
1005
+
1006
+ offset = residue_idx[:,:,None]-residue_idx[:,None,:]
1007
+ offset = gather_edges(offset[:,:,:,None], E_idx)[:,:,:,0] #[B, L, K]
1008
+
1009
+ d_chains = ((chain_labels[:, :, None] - chain_labels[:,None,:])==0).long() #find self vs non-self interaction
1010
+ E_chains = gather_edges(d_chains[:,:,:,None], E_idx)[:,:,:,0]
1011
+ E_positional = self.embeddings(offset.long(), E_chains)
1012
+ E = torch.cat((E_positional, RBF_all), -1)
1013
+ E = self.edge_embedding(E)
1014
+ E = self.norm_edges(E)
1015
+ return E, E_idx
1016
+
1017
+
1018
+
1019
+ class ProteinMPNN(nn.Module):
1020
+ def __init__(self, num_letters, node_features, edge_features,
1021
+ hidden_dim, num_encoder_layers=3, num_decoder_layers=3,
1022
+ vocab=21, k_neighbors=64, augment_eps=0.05, dropout=0.1, ca_only=False):
1023
+ super(ProteinMPNN, self).__init__()
1024
+
1025
+ # Hyperparameters
1026
+ self.node_features = node_features
1027
+ self.edge_features = edge_features
1028
+ self.hidden_dim = hidden_dim
1029
+
1030
+ # Featurization layers
1031
+ if ca_only:
1032
+ self.features = CA_ProteinFeatures(node_features, edge_features, top_k=k_neighbors, augment_eps=augment_eps)
1033
+ self.W_v = nn.Linear(node_features, hidden_dim, bias=True)
1034
+ else:
1035
+ self.features = ProteinFeatures(node_features, edge_features, top_k=k_neighbors, augment_eps=augment_eps)
1036
+
1037
+ self.W_e = nn.Linear(edge_features, hidden_dim, bias=True)
1038
+ self.W_s = nn.Embedding(vocab, hidden_dim)
1039
+
1040
+ # Encoder layers
1041
+ self.encoder_layers = nn.ModuleList([
1042
+ EncLayer(hidden_dim, hidden_dim*2, dropout=dropout)
1043
+ for _ in range(num_encoder_layers)
1044
+ ])
1045
+
1046
+ # Decoder layers
1047
+ self.decoder_layers = nn.ModuleList([
1048
+ DecLayer(hidden_dim, hidden_dim*3, dropout=dropout)
1049
+ for _ in range(num_decoder_layers)
1050
+ ])
1051
+ self.W_out = nn.Linear(hidden_dim, num_letters, bias=True)
1052
+
1053
+ for p in self.parameters():
1054
+ if p.dim() > 1:
1055
+ nn.init.xavier_uniform_(p)
1056
+
1057
+ def forward(self, X, S, mask, chain_M, residue_idx, chain_encoding_all, randn, use_input_decoding_order=False, decoding_order=None):
1058
+ """ Graph-conditioned sequence model """
1059
+ device=X.device
1060
+ # Prepare node and edge embeddings
1061
+ E, E_idx = self.features(X, mask, residue_idx, chain_encoding_all)
1062
+ h_V = torch.zeros((E.shape[0], E.shape[1], E.shape[-1]), device=E.device)
1063
+ h_E = self.W_e(E)
1064
+
1065
+ # Encoder is unmasked self-attention
1066
+ mask_attend = gather_nodes(mask.unsqueeze(-1), E_idx).squeeze(-1)
1067
+ mask_attend = mask.unsqueeze(-1) * mask_attend
1068
+ for layer in self.encoder_layers:
1069
+ h_V, h_E = layer(h_V, h_E, E_idx, mask, mask_attend)
1070
+
1071
+ # Concatenate sequence embeddings for autoregressive decoder
1072
+ h_S = self.W_s(S)
1073
+ h_ES = cat_neighbors_nodes(h_S, h_E, E_idx)
1074
+
1075
+ # Build encoder embeddings
1076
+ h_EX_encoder = cat_neighbors_nodes(torch.zeros_like(h_S), h_E, E_idx)
1077
+ h_EXV_encoder = cat_neighbors_nodes(h_V, h_EX_encoder, E_idx)
1078
+
1079
+
1080
+ chain_M = chain_M*mask #update chain_M to include missing regions
1081
+ if not use_input_decoding_order:
1082
+ decoding_order = torch.argsort((chain_M+0.0001)*(torch.abs(randn))) #[numbers will be smaller for places where chain_M = 0.0 and higher for places where chain_M = 1.0]
1083
+ mask_size = E_idx.shape[1]
1084
+ permutation_matrix_reverse = torch.nn.functional.one_hot(decoding_order, num_classes=mask_size).float()
1085
+ order_mask_backward = torch.einsum('ij, biq, bjp->bqp',(1-torch.triu(torch.ones(mask_size,mask_size, device=device))), permutation_matrix_reverse, permutation_matrix_reverse)
1086
+ mask_attend = torch.gather(order_mask_backward, 2, E_idx).unsqueeze(-1)
1087
+ mask_1D = mask.view([mask.size(0), mask.size(1), 1, 1])
1088
+ mask_bw = mask_1D * mask_attend
1089
+ mask_fw = mask_1D * (1. - mask_attend)
1090
+
1091
+ h_EXV_encoder_fw = mask_fw * h_EXV_encoder
1092
+ for layer in self.decoder_layers:
1093
+ # Masked positions attend to encoder information, unmasked see.
1094
+ h_ESV = cat_neighbors_nodes(h_V, h_ES, E_idx)
1095
+ h_ESV = mask_bw * h_ESV + h_EXV_encoder_fw
1096
+ h_V = layer(h_V, h_ESV, mask)
1097
+
1098
+ logits = self.W_out(h_V)
1099
+ log_probs = F.log_softmax(logits, dim=-1)
1100
+ return log_probs
1101
+
1102
+
1103
+
1104
+ def sample(self, X, randn, S_true, chain_mask, chain_encoding_all, residue_idx, mask=None, temperature=1.0, omit_AAs_np=None, bias_AAs_np=None, chain_M_pos=None, omit_AA_mask=None, pssm_coef=None, pssm_bias=None, pssm_multi=None, pssm_log_odds_flag=None, pssm_log_odds_mask=None, pssm_bias_flag=None, bias_by_res=None):
1105
+ device = X.device
1106
+ # Prepare node and edge embeddings
1107
+ E, E_idx = self.features(X, mask, residue_idx, chain_encoding_all)
1108
+ h_V = torch.zeros((E.shape[0], E.shape[1], E.shape[-1]), device=device)
1109
+ h_E = self.W_e(E)
1110
+
1111
+ # Encoder is unmasked self-attention
1112
+ mask_attend = gather_nodes(mask.unsqueeze(-1), E_idx).squeeze(-1)
1113
+ mask_attend = mask.unsqueeze(-1) * mask_attend
1114
+ for layer in self.encoder_layers:
1115
+ h_V, h_E = layer(h_V, h_E, E_idx, mask, mask_attend)
1116
+
1117
+ # Decoder uses masked self-attention
1118
+ chain_mask = chain_mask*chain_M_pos*mask #update chain_M to include missing regions
1119
+ decoding_order = torch.argsort((chain_mask+0.0001)*(torch.abs(randn))) #[numbers will be smaller for places where chain_M = 0.0 and higher for places where chain_M = 1.0]
1120
+ mask_size = E_idx.shape[1]
1121
+ permutation_matrix_reverse = torch.nn.functional.one_hot(decoding_order, num_classes=mask_size).float()
1122
+ order_mask_backward = torch.einsum('ij, biq, bjp->bqp',(1-torch.triu(torch.ones(mask_size,mask_size, device=device))), permutation_matrix_reverse, permutation_matrix_reverse)
1123
+ mask_attend = torch.gather(order_mask_backward, 2, E_idx).unsqueeze(-1)
1124
+ mask_1D = mask.view([mask.size(0), mask.size(1), 1, 1])
1125
+ mask_bw = mask_1D * mask_attend
1126
+ mask_fw = mask_1D * (1. - mask_attend)
1127
+
1128
+ N_batch, N_nodes = X.size(0), X.size(1)
1129
+ log_probs = torch.zeros((N_batch, N_nodes, 21), device=device)
1130
+ all_probs = torch.zeros((N_batch, N_nodes, 21), device=device, dtype=torch.float32)
1131
+ h_S = torch.zeros_like(h_V, device=device)
1132
+ S = torch.zeros((N_batch, N_nodes), dtype=torch.int64, device=device)
1133
+ h_V_stack = [h_V] + [torch.zeros_like(h_V, device=device) for _ in range(len(self.decoder_layers))]
1134
+ constant = torch.tensor(omit_AAs_np, device=device)
1135
+ constant_bias = torch.tensor(bias_AAs_np, device=device)
1136
+ #chain_mask_combined = chain_mask*chain_M_pos
1137
+ omit_AA_mask_flag = omit_AA_mask != None
1138
+
1139
+
1140
+ h_EX_encoder = cat_neighbors_nodes(torch.zeros_like(h_S), h_E, E_idx)
1141
+ h_EXV_encoder = cat_neighbors_nodes(h_V, h_EX_encoder, E_idx)
1142
+ h_EXV_encoder_fw = mask_fw * h_EXV_encoder
1143
+ for t_ in range(N_nodes):
1144
+ t = decoding_order[:,t_] #[B]
1145
+ chain_mask_gathered = torch.gather(chain_mask, 1, t[:,None]) #[B]
1146
+ mask_gathered = torch.gather(mask, 1, t[:,None]) #[B]
1147
+ bias_by_res_gathered = torch.gather(bias_by_res, 1, t[:,None,None].repeat(1,1,21))[:,0,:] #[B, 21]
1148
+ if (mask_gathered==0).all(): #for padded or missing regions only
1149
+ S_t = torch.gather(S_true, 1, t[:,None])
1150
+ else:
1151
+ # Hidden layers
1152
+ E_idx_t = torch.gather(E_idx, 1, t[:,None,None].repeat(1,1,E_idx.shape[-1]))
1153
+ h_E_t = torch.gather(h_E, 1, t[:,None,None,None].repeat(1,1,h_E.shape[-2], h_E.shape[-1]))
1154
+ h_ES_t = cat_neighbors_nodes(h_S, h_E_t, E_idx_t)
1155
+ h_EXV_encoder_t = torch.gather(h_EXV_encoder_fw, 1, t[:,None,None,None].repeat(1,1,h_EXV_encoder_fw.shape[-2], h_EXV_encoder_fw.shape[-1]))
1156
+ mask_t = torch.gather(mask, 1, t[:,None])
1157
+ for l, layer in enumerate(self.decoder_layers):
1158
+ # Updated relational features for future states
1159
+ h_ESV_decoder_t = cat_neighbors_nodes(h_V_stack[l], h_ES_t, E_idx_t)
1160
+ h_V_t = torch.gather(h_V_stack[l], 1, t[:,None,None].repeat(1,1,h_V_stack[l].shape[-1]))
1161
+ h_ESV_t = torch.gather(mask_bw, 1, t[:,None,None,None].repeat(1,1,mask_bw.shape[-2], mask_bw.shape[-1])) * h_ESV_decoder_t + h_EXV_encoder_t
1162
+ h_V_stack[l+1].scatter_(1, t[:,None,None].repeat(1,1,h_V.shape[-1]), layer(h_V_t, h_ESV_t, mask_V=mask_t))
1163
+ # Sampling step
1164
+ h_V_t = torch.gather(h_V_stack[-1], 1, t[:,None,None].repeat(1,1,h_V_stack[-1].shape[-1]))[:,0]
1165
+ logits = self.W_out(h_V_t) / temperature
1166
+ probs = F.softmax(logits-constant[None,:]*1e8+constant_bias[None,:]/temperature+bias_by_res_gathered/temperature, dim=-1)
1167
+ if pssm_bias_flag:
1168
+ pssm_coef_gathered = torch.gather(pssm_coef, 1, t[:,None])[:,0]
1169
+ pssm_bias_gathered = torch.gather(pssm_bias, 1, t[:,None,None].repeat(1,1,pssm_bias.shape[-1]))[:,0]
1170
+ probs = (1-pssm_multi*pssm_coef_gathered[:,None])*probs + pssm_multi*pssm_coef_gathered[:,None]*pssm_bias_gathered
1171
+ if pssm_log_odds_flag:
1172
+ pssm_log_odds_mask_gathered = torch.gather(pssm_log_odds_mask, 1, t[:,None, None].repeat(1,1,pssm_log_odds_mask.shape[-1]))[:,0] #[B, 21]
1173
+ probs_masked = probs*pssm_log_odds_mask_gathered
1174
+ probs_masked += probs * 0.001
1175
+ probs = probs_masked/torch.sum(probs_masked, dim=-1, keepdim=True) #[B, 21]
1176
+ if omit_AA_mask_flag:
1177
+ omit_AA_mask_gathered = torch.gather(omit_AA_mask, 1, t[:,None, None].repeat(1,1,omit_AA_mask.shape[-1]))[:,0] #[B, 21]
1178
+ probs_masked = probs*(1.0-omit_AA_mask_gathered)
1179
+ probs = probs_masked/torch.sum(probs_masked, dim=-1, keepdim=True) #[B, 21]
1180
+ S_t = torch.multinomial(probs, 1)
1181
+ all_probs.scatter_(1, t[:,None,None].repeat(1,1,21), (chain_mask_gathered[:,:,None,]*probs[:,None,:]).float())
1182
+ S_true_gathered = torch.gather(S_true, 1, t[:,None])
1183
+ S_t = (S_t*chain_mask_gathered+S_true_gathered*(1.0-chain_mask_gathered)).long()
1184
+ temp1 = self.W_s(S_t)
1185
+ h_S.scatter_(1, t[:,None,None].repeat(1,1,temp1.shape[-1]), temp1)
1186
+ S.scatter_(1, t[:,None], S_t)
1187
+ output_dict = {"S": S, "probs": all_probs, "decoding_order": decoding_order}
1188
+ return output_dict
1189
+
1190
+
1191
+ def tied_sample(self, X, randn, S_true, chain_mask, chain_encoding_all, residue_idx, mask=None, temperature=1.0, omit_AAs_np=None, bias_AAs_np=None, chain_M_pos=None, omit_AA_mask=None, pssm_coef=None, pssm_bias=None, pssm_multi=None, pssm_log_odds_flag=None, pssm_log_odds_mask=None, pssm_bias_flag=None, tied_pos=None, tied_beta=None, bias_by_res=None):
1192
+ device = X.device
1193
+ # Prepare node and edge embeddings
1194
+ E, E_idx = self.features(X, mask, residue_idx, chain_encoding_all)
1195
+ h_V = torch.zeros((E.shape[0], E.shape[1], E.shape[-1]), device=device)
1196
+ h_E = self.W_e(E)
1197
+ # Encoder is unmasked self-attention
1198
+ mask_attend = gather_nodes(mask.unsqueeze(-1), E_idx).squeeze(-1)
1199
+ mask_attend = mask.unsqueeze(-1) * mask_attend
1200
+ for layer in self.encoder_layers:
1201
+ h_V, h_E = layer(h_V, h_E, E_idx, mask, mask_attend)
1202
+
1203
+ # Decoder uses masked self-attention
1204
+ chain_mask = chain_mask*chain_M_pos*mask #update chain_M to include missing regions
1205
+ decoding_order = torch.argsort((chain_mask+0.0001)*(torch.abs(randn))) #[numbers will be smaller for places where chain_M = 0.0 and higher for places where chain_M = 1.0]
1206
+
1207
+ new_decoding_order = []
1208
+ for t_dec in list(decoding_order[0,].cpu().data.numpy()):
1209
+ if t_dec not in list(itertools.chain(*new_decoding_order)):
1210
+ list_a = [item for item in tied_pos if t_dec in item]
1211
+ if list_a:
1212
+ new_decoding_order.append(list_a[0])
1213
+ else:
1214
+ new_decoding_order.append([t_dec])
1215
+ decoding_order = torch.tensor(list(itertools.chain(*new_decoding_order)), device=device)[None,].repeat(X.shape[0],1)
1216
+
1217
+ mask_size = E_idx.shape[1]
1218
+ permutation_matrix_reverse = torch.nn.functional.one_hot(decoding_order, num_classes=mask_size).float()
1219
+ order_mask_backward = torch.einsum('ij, biq, bjp->bqp',(1-torch.triu(torch.ones(mask_size,mask_size, device=device))), permutation_matrix_reverse, permutation_matrix_reverse)
1220
+ mask_attend = torch.gather(order_mask_backward, 2, E_idx).unsqueeze(-1)
1221
+ mask_1D = mask.view([mask.size(0), mask.size(1), 1, 1])
1222
+ mask_bw = mask_1D * mask_attend
1223
+ mask_fw = mask_1D * (1. - mask_attend)
1224
+
1225
+ N_batch, N_nodes = X.size(0), X.size(1)
1226
+ log_probs = torch.zeros((N_batch, N_nodes, 21), device=device)
1227
+ all_probs = torch.zeros((N_batch, N_nodes, 21), device=device, dtype=torch.float32)
1228
+ h_S = torch.zeros_like(h_V, device=device)
1229
+ S = torch.zeros((N_batch, N_nodes), dtype=torch.int64, device=device)
1230
+ h_V_stack = [h_V] + [torch.zeros_like(h_V, device=device) for _ in range(len(self.decoder_layers))]
1231
+ constant = torch.tensor(omit_AAs_np, device=device)
1232
+ constant_bias = torch.tensor(bias_AAs_np, device=device)
1233
+ omit_AA_mask_flag = omit_AA_mask != None
1234
+
1235
+ h_EX_encoder = cat_neighbors_nodes(torch.zeros_like(h_S), h_E, E_idx)
1236
+ h_EXV_encoder = cat_neighbors_nodes(h_V, h_EX_encoder, E_idx)
1237
+ h_EXV_encoder_fw = mask_fw * h_EXV_encoder
1238
+ for t_list in new_decoding_order:
1239
+ logits = 0.0
1240
+ logit_list = []
1241
+ done_flag = False
1242
+ for t in t_list:
1243
+ if (mask[:,t]==0).all():
1244
+ S_t = S_true[:,t]
1245
+ for t in t_list:
1246
+ h_S[:,t,:] = self.W_s(S_t)
1247
+ S[:,t] = S_t
1248
+ done_flag = True
1249
+ break
1250
+ else:
1251
+ E_idx_t = E_idx[:,t:t+1,:]
1252
+ h_E_t = h_E[:,t:t+1,:,:]
1253
+ h_ES_t = cat_neighbors_nodes(h_S, h_E_t, E_idx_t)
1254
+ h_EXV_encoder_t = h_EXV_encoder_fw[:,t:t+1,:,:]
1255
+ mask_t = mask[:,t:t+1]
1256
+ for l, layer in enumerate(self.decoder_layers):
1257
+ h_ESV_decoder_t = cat_neighbors_nodes(h_V_stack[l], h_ES_t, E_idx_t)
1258
+ h_V_t = h_V_stack[l][:,t:t+1,:]
1259
+ h_ESV_t = mask_bw[:,t:t+1,:,:] * h_ESV_decoder_t + h_EXV_encoder_t
1260
+ h_V_stack[l+1][:,t,:] = layer(h_V_t, h_ESV_t, mask_V=mask_t).squeeze(1)
1261
+ h_V_t = h_V_stack[-1][:,t,:]
1262
+ logit_list.append((self.W_out(h_V_t) / temperature)/len(t_list))
1263
+ logits += tied_beta[t]*(self.W_out(h_V_t) / temperature)/len(t_list)
1264
+ if done_flag:
1265
+ pass
1266
+ else:
1267
+ bias_by_res_gathered = bias_by_res[:,t,:] #[B, 21]
1268
+ probs = F.softmax(logits-constant[None,:]*1e8+constant_bias[None,:]/temperature+bias_by_res_gathered/temperature, dim=-1)
1269
+ if pssm_bias_flag:
1270
+ pssm_coef_gathered = pssm_coef[:,t]
1271
+ pssm_bias_gathered = pssm_bias[:,t]
1272
+ probs = (1-pssm_multi*pssm_coef_gathered[:,None])*probs + pssm_multi*pssm_coef_gathered[:,None]*pssm_bias_gathered
1273
+ if pssm_log_odds_flag:
1274
+ pssm_log_odds_mask_gathered = pssm_log_odds_mask[:,t]
1275
+ probs_masked = probs*pssm_log_odds_mask_gathered
1276
+ probs_masked += probs * 0.001
1277
+ probs = probs_masked/torch.sum(probs_masked, dim=-1, keepdim=True) #[B, 21]
1278
+ if omit_AA_mask_flag:
1279
+ omit_AA_mask_gathered = omit_AA_mask[:,t]
1280
+ probs_masked = probs*(1.0-omit_AA_mask_gathered)
1281
+ probs = probs_masked/torch.sum(probs_masked, dim=-1, keepdim=True) #[B, 21]
1282
+ S_t_repeat = torch.multinomial(probs, 1).squeeze(-1)
1283
+ S_t_repeat = (chain_mask[:,t]*S_t_repeat + (1-chain_mask[:,t])*S_true[:,t]).long() #hard pick fixed positions
1284
+ for t in t_list:
1285
+ h_S[:,t,:] = self.W_s(S_t_repeat)
1286
+ S[:,t] = S_t_repeat
1287
+ all_probs[:,t,:] = probs.float()
1288
+ output_dict = {"S": S, "probs": all_probs, "decoding_order": decoding_order}
1289
+ return output_dict
1290
+
1291
+
1292
+ def conditional_probs(self, X, S, mask, chain_M, residue_idx, chain_encoding_all, randn, backbone_only=False):
1293
+ """ Graph-conditioned sequence model """
1294
+ device=X.device
1295
+ # Prepare node and edge embeddings
1296
+ E, E_idx = self.features(X, mask, residue_idx, chain_encoding_all)
1297
+ h_V_enc = torch.zeros((E.shape[0], E.shape[1], E.shape[-1]), device=E.device)
1298
+ h_E = self.W_e(E)
1299
+
1300
+ # Encoder is unmasked self-attention
1301
+ mask_attend = gather_nodes(mask.unsqueeze(-1), E_idx).squeeze(-1)
1302
+ mask_attend = mask.unsqueeze(-1) * mask_attend
1303
+ for layer in self.encoder_layers:
1304
+ h_V_enc, h_E = layer(h_V_enc, h_E, E_idx, mask, mask_attend)
1305
+
1306
+ # Concatenate sequence embeddings for autoregressive decoder
1307
+ h_S = self.W_s(S)
1308
+ h_ES = cat_neighbors_nodes(h_S, h_E, E_idx)
1309
+
1310
+ # Build encoder embeddings
1311
+ h_EX_encoder = cat_neighbors_nodes(torch.zeros_like(h_S), h_E, E_idx)
1312
+ h_EXV_encoder = cat_neighbors_nodes(h_V_enc, h_EX_encoder, E_idx)
1313
+
1314
+
1315
+ chain_M = chain_M*mask #update chain_M to include missing regions
1316
+
1317
+ chain_M_np = chain_M.cpu().numpy()
1318
+ idx_to_loop = np.argwhere(chain_M_np[0,:]==1)[:,0]
1319
+ log_conditional_probs = torch.zeros([X.shape[0], chain_M.shape[1], 21], device=device).float()
1320
+
1321
+ for idx in idx_to_loop:
1322
+ h_V = torch.clone(h_V_enc)
1323
+ order_mask = torch.zeros(chain_M.shape[1], device=device).float()
1324
+ if backbone_only:
1325
+ order_mask = torch.ones(chain_M.shape[1], device=device).float()
1326
+ order_mask[idx] = 0.
1327
+ else:
1328
+ order_mask = torch.zeros(chain_M.shape[1], device=device).float()
1329
+ order_mask[idx] = 1.
1330
+ decoding_order = torch.argsort((order_mask[None,]+0.0001)*(torch.abs(randn))) #[numbers will be smaller for places where chain_M = 0.0 and higher for places where chain_M = 1.0]
1331
+ mask_size = E_idx.shape[1]
1332
+ permutation_matrix_reverse = torch.nn.functional.one_hot(decoding_order, num_classes=mask_size).float()
1333
+ order_mask_backward = torch.einsum('ij, biq, bjp->bqp',(1-torch.triu(torch.ones(mask_size,mask_size, device=device))), permutation_matrix_reverse, permutation_matrix_reverse)
1334
+ mask_attend = torch.gather(order_mask_backward, 2, E_idx).unsqueeze(-1)
1335
+ mask_1D = mask.view([mask.size(0), mask.size(1), 1, 1])
1336
+ mask_bw = mask_1D * mask_attend
1337
+ mask_fw = mask_1D * (1. - mask_attend)
1338
+
1339
+ h_EXV_encoder_fw = mask_fw * h_EXV_encoder
1340
+ for layer in self.decoder_layers:
1341
+ # Masked positions attend to encoder information, unmasked see.
1342
+ h_ESV = cat_neighbors_nodes(h_V, h_ES, E_idx)
1343
+ h_ESV = mask_bw * h_ESV + h_EXV_encoder_fw
1344
+ h_V = layer(h_V, h_ESV, mask)
1345
+
1346
+ logits = self.W_out(h_V)
1347
+ log_probs = F.log_softmax(logits, dim=-1)
1348
+ log_conditional_probs[:,idx,:] = log_probs[:,idx,:]
1349
+ return log_conditional_probs
1350
+
1351
+
1352
+ def unconditional_probs(self, X, mask, residue_idx, chain_encoding_all):
1353
+ """ Graph-conditioned sequence model """
1354
+ device=X.device
1355
+ # Prepare node and edge embeddings
1356
+ E, E_idx = self.features(X, mask, residue_idx, chain_encoding_all)
1357
+ h_V = torch.zeros((E.shape[0], E.shape[1], E.shape[-1]), device=E.device)
1358
+ h_E = self.W_e(E)
1359
+
1360
+ # Encoder is unmasked self-attention
1361
+ mask_attend = gather_nodes(mask.unsqueeze(-1), E_idx).squeeze(-1)
1362
+ mask_attend = mask.unsqueeze(-1) * mask_attend
1363
+ for layer in self.encoder_layers:
1364
+ h_V, h_E = layer(h_V, h_E, E_idx, mask, mask_attend)
1365
+
1366
+ # Build encoder embeddings
1367
+ h_EX_encoder = cat_neighbors_nodes(torch.zeros_like(h_V), h_E, E_idx)
1368
+ h_EXV_encoder = cat_neighbors_nodes(h_V, h_EX_encoder, E_idx)
1369
+
1370
+ order_mask_backward = torch.zeros([X.shape[0], X.shape[1], X.shape[1]], device=device)
1371
+ mask_attend = torch.gather(order_mask_backward, 2, E_idx).unsqueeze(-1)
1372
+ mask_1D = mask.view([mask.size(0), mask.size(1), 1, 1])
1373
+ mask_bw = mask_1D * mask_attend
1374
+ mask_fw = mask_1D * (1. - mask_attend)
1375
+
1376
+ h_EXV_encoder_fw = mask_fw * h_EXV_encoder
1377
+ for layer in self.decoder_layers:
1378
+ h_V = layer(h_V, h_EXV_encoder_fw, mask)
1379
+
1380
+ logits = self.W_out(h_V)
1381
+ log_probs = F.log_softmax(logits, dim=-1)
1382
+ return log_probs
1383
+
stamp/weights/stamp_s0.pt ADDED
@@ -0,0 +1,3 @@
 
 
 
 
1
+ version https://git-lfs.github.com/spec/v1
2
+ oid sha256:bb879d3a289e7c3d67a32ca7ad9f8ca0861ecb1427b6160e6d26c3622be51707
3
+ size 2384777
stamp/weights/stamp_s1.pt ADDED
@@ -0,0 +1,3 @@
 
 
 
 
1
+ version https://git-lfs.github.com/spec/v1
2
+ oid sha256:7623995807388936203413c1531bc2da539633799a7c8434eb0344ef7beefd46
3
+ size 2384777
stamp/weights/stamp_s2.pt ADDED
@@ -0,0 +1,3 @@
 
 
 
 
1
+ version https://git-lfs.github.com/spec/v1
2
+ oid sha256:ecb89d2c622c946fe86c04dcc187eb710e9d424ca6a7ecde0397b3485f9fc664
3
+ size 2384777
stamp/weights/stamp_s3.pt ADDED
@@ -0,0 +1,3 @@
 
 
 
 
1
+ version https://git-lfs.github.com/spec/v1
2
+ oid sha256:a27047c7a2b690df2a6738268208a9c6c9546c4106ca1f891403991d1776baa4
3
+ size 2384777
stamp/weights/stamp_s4.pt ADDED
@@ -0,0 +1,3 @@
 
 
 
 
1
+ version https://git-lfs.github.com/spec/v1
2
+ oid sha256:b6aa57361368c6b132d3ce33a1f88009a96bf0367eb702c7116aea02cf7aae18
3
+ size 2384777
stamp/weights/v_48_020.pt ADDED
@@ -0,0 +1,3 @@
 
 
 
 
1
+ version https://git-lfs.github.com/spec/v1
2
+ oid sha256:c9cb4a671d79604111231f8dbfc7c590e06f1197453b7a6854ac6661a642f5bd
3
+ size 6681301