Chainsaw / model /domain_chop.py
wuxing0105's picture
Upload folder using huggingface_hub
80a72c3 verified
Raw History Blame Contribute Delete
13.9 kB
"""Domain predictor classes.
"""
import os
import glob
import numpy as np
import pandas as pd
import torch
from torch import nn
from model.domain_assignment.util import make_pair_labels, make_domain_mapping_dict
import logging
LOG = logging.getLogger(__name__)
def get_checkpoint_epoch(checkpoint_file):
return int(os.path.splitext(checkpoint_file)[0].split(".")[-1])
class PairwiseDomainPredictor(nn.Module):
"""Wrapper for a pairwise domain co-membership predictor, adding in domain prediction post-processing."""
def __init__(
self,
model,
domain_caller,
device,
loss="bce",
x_has_padding_mask=True,
mask_padding=True,
n_checkpoints_to_average=1,
checkpoint_dir=None,
load_checkpoint_if_exists=False,
save_val_best=True,
max_recycles=0,
post_process_domains=True,
min_ss_components=2,
min_domain_length=30,
remove_disordered_domain_threshold=0,
trim_each_domain=True,
dist_transform_type="min_replace_inverse",
distance_denominator=10,
):
super().__init__()
self._train_model = model # we want to keep this hanging around so that optimizer references dont break
self.model = self._train_model
self.domain_caller = domain_caller
self.device = device
self.x_has_padding_mask = x_has_padding_mask
self.mask_padding = mask_padding # if True use padding mask to mask loss
self.n_checkpoints_to_average = n_checkpoints_to_average
self.checkpoint_dir = checkpoint_dir
self._epoch = 0
self.save_val_best = save_val_best
self.best_val_metrics = {}
self.max_recycles = max_recycles
self.post_process_domains = post_process_domains
self.remove_disordered_domain_threshold = remove_disordered_domain_threshold
self.trim_each_domain = trim_each_domain
self.min_domain_length = min_domain_length
self.min_ss_components = min_ss_components
self.dist_transform_type = dist_transform_type
self.distance_denominator = distance_denominator
if load_checkpoint_if_exists:
checkpoint_files = sorted(
glob.glob(os.path.join(self.checkpoint_dir, "weights*")),
key=get_checkpoint_epoch,
reverse=True,
)
if len(checkpoint_files) > 0:
self._epoch = get_checkpoint_epoch(checkpoint_files[0])
LOG.info(f"Loading saved checkpoint(s) ending at epoch {self._epoch}")
self.load_checkpoints(average=True)
self.load_checkpoints()
else:
LOG.info("No checkpoints found to load")
if loss == "bce":
self.loss_function = nn.BCELoss(reduction="none")
elif loss == "mse":
self.loss_function = nn.MSELoss(reduction="none")
def load_checkpoints(self, average=False, old_style=False):
if self.n_checkpoints_to_average == 1:
weights_file = os.path.join(self.checkpoint_dir, "weights.pt")
else:
# for e.g. resetting training weights for next training epoch after testing with avg
LOG.info(f"Loading last checkpoint (epoch {self._epoch})")
weights_file = os.path.join(self.checkpoint_dir, f"weights.{self._epoch}.pt")
LOG.info(f"Loading weights from: {weights_file}")
state_dict = torch.load(weights_file, map_location=self.device)
if old_style:
self.load_state_dict(state_dict, strict=False)
else:
self.model.load_state_dict(state_dict)
def predict_pairwise(self, x):
x = x.to(self.device)
if np.isnan(x.cpu().numpy()).any():
raise Exception('NAN values in data')
y_pred = self.model(x).squeeze(1) # b, L, L
assert y_pred.ndim == 3
return y_pred
def get_mask(self, x):
"""Binary mask 1 for observed, 0 for padding."""
x = x.to(self.device)
if self.x_has_padding_mask:
mask = 1 - x[:, -1] # b, L, L
else:
mask = None
return mask
def epoch_start(self):
self.model = self._train_model
self.model.train()
self._epoch += 1
def test_begin(self):
if self.n_checkpoints_to_average > 1:
torch.save(self.model.state_dict(), os.path.join(self.checkpoint_dir, f'weights.{self._epoch}.pt'))
start_idx = self._epoch - self.n_checkpoints_to_average
if start_idx >= 2:
os.remove(os.path.join(self.checkpoint_dir, f"weights.{start_idx-1}.pt"))
else:
torch.save(self.model.state_dict(), os.path.join(self.checkpoint_dir, 'weights.pt'))
if self.n_checkpoints_to_average > 1:
# self.model.to("cpu") # free up gpu memory for average model
self.load_checkpoints(average=True)
self.model.eval()
def forward(self, x, y, batch_average=True):
"""A training step."""
x, y = x.to(self.device), y.to(self.device)
y_pred = self.predict_pairwise(x)
mask = self.get_mask(x)
return self.compute_loss(y_pred, y, mask=mask)
def compute_loss(self, y_pred, y, mask=None, batch_average=True):
y_pred, y = y_pred.to(self.device), y.to(self.device)
if mask is None or not self.mask_padding:
mask = torch.ones_like(y)
# mask is b, L, L. To normalise correctly, we need to divide by number of observations
loss = (self.loss_function(y_pred, y)*mask).sum((-1,-2)) / mask.sum((-1,-2))
# metrics characterising inputs: how many residues, how many with domain assignments.
labelled_residues = ((y*mask).sum(-1) > 0).sum(-1) # b
non_padding_residues = (mask.sum(-1) > 0).sum(-1) # b
labelled_frac = labelled_residues / non_padding_residues
metrics = {
"labelled_residues": labelled_residues.detach().cpu().numpy(),
"residues": non_padding_residues.detach().cpu().numpy(),
"labelled_frac": labelled_frac.detach().cpu().numpy(),
"loss": loss.detach().cpu().numpy(),
}
if batch_average:
loss = loss.mean(0)
metrics = {k: np.mean(v) for k, v in metrics.items()}
return loss, metrics
def domains_from_pairwise(self, y_pred):
assert y_pred.ndim == 3
domain_preds = []
confidence_list = []
for pred_single in y_pred.cpu().numpy():
single_domains, confidence = self.domain_caller(pred_single)
domain_preds.append(single_domains)
confidence_list.append(confidence)
return domain_preds, confidence_list
def distance_transform(self, x):
dist_chan = x[0, 0]
# Find the minimum non-zero value in the channel
min_nonzero = dist_chan[dist_chan > 0].min()
# Replace zero values in the channel with the minimum non-zero value
dist_chan[dist_chan == 0] = min_nonzero
if self.dist_transform_type == "min_replace_inverse":
# replace zero values and then invert.
dist_chan = dist_chan ** (-1)
x[0, 0] = dist_chan
return x
elif self.dist_transform_type == "unidoc_exponent": # replace zero values in pae / distance
spread = self.distance_denominator
dist_chan = (1 + np.exp((dist_chan - 8) / spread)) ** -1
x[0,0] = dist_chan
return x
@torch.no_grad()
def predict(self, x, return_pairwise=True):
x = self.distance_transform(x)
if self.max_recycles > 0:
for i in range(self.max_recycles):
# add recycling channels
n_res = x.shape[-1]
recycle_channels = torch.zeros(1, 2, n_res, n_res)
# Concatenate the original tensor and the zeros tensor along the second dimension
x = torch.cat((x, recycle_channels), dim=1)
x = self.recycle_predict(x)
y_pred = self.predict_pairwise(x)
domain_dicts, confidence = self.domains_from_pairwise(y_pred)
if self.post_process_domains:
domain_dicts = self.post_process(domain_dicts, x) # todo move this to domains from pairwise function
if return_pairwise:
return y_pred, domain_dicts, confidence
else:
return domain_dicts, confidence
@torch.no_grad()
def recycle_predict(self, x):
x = x.to(self.device)
y_pred = self.predict_pairwise(x)
domain_dicts, confidence = self.domains_from_pairwise(y_pred)
y_pred_from_domains = np.array(
[make_pair_labels(n_res=x.shape[-1], domain_dict=d_dict) for d_dict in domain_dicts])
y_pred_from_domains = torch.tensor(y_pred_from_domains).to(self.device)
if self.x_has_padding_mask:
x[:, -2, :, :] = y_pred # assumes that last dimension is padding mask
x[:, -3, :, :] = y_pred_from_domains
else:
x[:, -1, :, :] = y_pred
x[:, -2, :, :] = y_pred_from_domains
return x
def post_process(self, domain_dicts, x_batch):
new_domain_dicts = []
for domain_dict, x in zip(domain_dicts, x_batch):
x = x.cpu().numpy()
domain_dict = {k: list(v) for k, v in domain_dict.items()}
helix, sheet = x[1], x[2]
diag_helix = np.diagonal(helix)
diag_sheet = np.diagonal(sheet)
ss_residues = list(np.where(diag_helix > 0)[0]) + list(np.where(diag_sheet > 0)[0])
domain_dict = self.trim_disordered_boundaries(domain_dict, ss_residues)
if self.remove_disordered_domain_threshold > 0:
domain_dict = self.remove_disordered_domains(domain_dict, ss_residues)
if self.min_ss_components > 0:
domain_dict = self.remove_domains_with_few_ss_components(domain_dict, x)
if self.min_domain_length > 0:
domain_dict = self.remove_domains_with_short_length(domain_dict)
new_domain_dicts.append(domain_dict)
return new_domain_dicts
def trim_disordered_boundaries(self, domain_dict, ss_residues):
if not self.trim_each_domain:
start = min(ss_residues)
end = max(ss_residues)
for dname, res in domain_dict.items():
if dname == "linker":
continue
if self.trim_each_domain:
domain_specific_ss = set(ss_residues).intersection(set(res))
if len(domain_specific_ss) == 0:
continue
start = min(domain_specific_ss)
end = max(domain_specific_ss)
domain_dict["linker"] += [r for r in res if r < start or r > end]
domain_dict[dname] = [r for r in res if r >= start and r <= end]
return domain_dict
def remove_disordered_domains(self, domain_dict, ss_residues):
new_domain_dict = {}
for dname, res in domain_dict.items():
if dname == "linker":
continue
if len(res) == 0:
continue
if len(set(res).intersection(set(ss_residues))) / len(res) < self.remove_disordered_domain_threshold:
domain_dict["linker"] += res
else:
new_domain_dict[dname] = res
new_domain_dict["linker"] = domain_dict["linker"]
return new_domain_dict
def remove_domains_with_few_ss_components(self, domain_dict, x):
"""
Remove domains where number of ss components is less than minimum
eg if self.min_ss_components=2 domains made of only a single helix or sheet are removed
achieve this by counting the number of unique string hashes in domain rows of x
"""
new_domain_dict = {}
for dname, res in domain_dict.items():
if dname == "linker":
continue
res = sorted(res)
helix = x[1][res, :][:, res]
strand = x[2][res, :][:, res]
helix = helix[np.any(helix, axis=1)]
strand = strand[np.any(strand, axis=1)]
# residues in the same secondary structure component have the same representation in the helix or strand matrix
n_helix = len(set(["".join([str(int(i)) for i in row]) for row in helix]))
n_sheet = len(set(["".join([str(int(i)) for i in row]) for row in strand]))
if len(res) == 0:
continue
if n_helix + n_sheet < self.min_ss_components:
domain_dict["linker"] += res
else:
new_domain_dict[dname] = res
new_domain_dict["linker"] = domain_dict["linker"]
return new_domain_dict
def remove_domains_with_short_length(self, domain_dict):
"""
Remove domains where length is less than minimum
"""
new_domain_dict = {}
for dname, res in domain_dict.items():
if dname == "linker":
continue
if len(res) < self.min_domain_length:
domain_dict["linker"] += res
else:
new_domain_dict[dname] = res
new_domain_dict["linker"] = domain_dict["linker"]
return new_domain_dict
class CSVDomainPredictor:
def __init__(self, csv_predictions):
self.csv_filepath = csv_predictions
self.predictions = pd.read_csv(csv_predictions)
def predict(self, x):
# x should be a chain_id
one_pred = self.predictions[self.predictions.chain_id == x[0][:5]]
if len(one_pred) == 0:
return None
domain_dicts = [make_domain_mapping_dict(one_pred.iloc[0])]
return None, domain_dicts