SpheronizaTor: Spherical Voxelization for Interpretable Protein Microenvironment Modeling.

Code ↔ Paper

5 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 5 matches

[1] § Methods › α-Helix classification use case ↔ examples/3DCNN-SpheronizaTor.py, lines 791–913 · score 0.90 · weight decay, AdamW, Dice loss, BCE, optimizer, IoU
[2] § Methods › α-Helix classification use case ↔ examples/3DCNN-SpheronizaTor.py, lines 188–237 · score 0.81 · linear layers, LeakyReLU, dilation, head, stem, padding
[3] § Methods › Construction of residue-centered spatial representations ↔ src/spheronizator/voxelBuilder.py, lines 143–154 · score 0.57 · buildSphere, get_boxProjection, position, spheres, voxel, atoms
[4] § Methods › Protein structure acquisition and preprocessing ↔ src/spheronizator/mol2parser.py, lines 13–139 · score 0.56 · mol2parser, Tripos, Biopython, parsed, matching, atomic
[5] § Methods › Construction of residue-centered spatial representations ↔ src/spheronizator/functions.py, lines 168–221 · score 0.50 · get_boxProjection, vector, distance, position, atoms, residue

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Python · 1,095 lines · 35 KB · BSD-3-Clause · 2 matches

#!/usr/bin/env python
# coding: utf-8
# In[1]:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
from __future__ import annotations
import os
import json
import warnings
import argparse
from dataclasses import dataclass, asdict
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import torch
import torch.nn as nn
import torch.optim as optim
import torch.nn.functional as F
from torch.utils.data import Dataset, DataLoader
import MDAnalysis as mda
from MDAnalysis.analysis.dssp import DSSP
from pyuul import utils, VolumeMaker
from sklearn.metrics import (
roc_curve,
auc,
confusion_matrix,
ConfusionMatrixDisplay,
roc_auc_score,
precision_score,
recall_score,
f1_score,
matthews_corrcoef,
accuracy_score,
)
from sklearn.model_selection import train_test_split
# =========================================================
# CONFIG
# =========================================================
@dataclass
class Config:
backend: str = "spheronizator" # pyuul | spheronizator | compare
pdb_dir: Path = Path("data/classification_data")
pdb_pattern: str = "*fixed.pdb"
sphero_dataset_dir: Path = Path(".")
sphero_atoms_dir: Path = Path("output_vox_atoms")
sphero_bonds_dir: Path = Path("output_vox_bonds")
resolution: float = 1.0
epochs: int = 10
learning_rate: float = 1e-3
batch_size: int = 1
num_workers: int = 0
pin_memory: bool = False
use_amp: bool = True
plot_first_sample: bool = False
save_checkpoint_every: int = 1
sphero_aggregate: str = "sum" # sum | mean | max
device: str = "cuda" if torch.cuda.is_available() else "cpu"
target_build_device: str = "cuda" if torch.cuda.is_available() else "cpu"
debug: bool = False
test_size: float = 0.2
random_state: int = 42
cls_threshold: float = 0.5
output_dir: Path = Path("outputs_classifier")
def parse_args():
parser = argparse.ArgumentParser(description="3D CNN protein classifier with pyuul vs spheronizator comparison")
parser.add_argument("--backend", default="spheronizator", choices=["spheronizator", "pyuul", "compare"])
parser.add_argument("--epochs", type=int, default=10)
parser.add_argument("--lr", type=float, default=1e-2)
parser.add_argument("--resolution", type=float, default=1)
parser.add_argument("--aggregate", default="sum", choices=["sum", "mean", "max"])
parser.add_argument("--batch_size", type=int, default=1)
parser.add_argument("--num_workers", type=int, default=0)
parser.add_argument("--plot_sample", action="store_true")
parser.add_argument("--debug", action="store_true")
parser.add_argument("--test_size", type=float, default=0.2)
parser.add_argument("--cls_threshold", type=float, default=0.25)
parser.add_argument("--output_dir", type=str, default="outputs_classifier")
return parser.parse_args()
def build_config():
args = parse_args()
cfg = Config()
cfg.backend = args.backend
cfg.epochs = args.epochs
cfg.learning_rate = args.lr
cfg.resolution = args.resolution
cfg.sphero_aggregate = args.aggregate
cfg.batch_size = args.batch_size
cfg.num_workers = args.num_workers
cfg.plot_first_sample = args.plot_sample
cfg.debug = args.debug
cfg.test_size = args.test_size
cfg.cls_threshold = args.cls_threshold
cfg.output_dir = Path(args.output_dir)
cfg.sphero_atoms_dir = cfg.sphero_dataset_dir / "output_vox_atoms"
cfg.sphero_bonds_dir = cfg.sphero_dataset_dir / "output_vox_bonds"
os.environ["PYTORCH_CUDA_ALLOC_CONF"] = "expandable_segments:True"
cfg.output_dir.mkdir(parents=True, exist_ok=True)
return cfg
def backend_output_dir(cfg: Config, backend: str) -> Path:
out = cfg.output_dir / backend
out.mkdir(parents=True, exist_ok=True)
return out
def debug_print(cfg: Config, *args):
if cfg.debug:
print(*args)
# =========================================================
# MODEL
# =========================================================
class ResidualBlock3D(nn.Module):
def __init__(self, in_channels: int, out_channels: int, groups: int = 4, dropout: float = 0.0):
super().__init__()
g = min(groups, out_channels)
while out_channels % g != 0:
g -= 1
self.conv1 = nn.Conv3d(in_channels, out_channels, kernel_size=3, padding=1, bias=False)
self.norm1 = nn.GroupNorm(g, out_channels)
self.act1 = nn.LeakyReLU(0.1, inplace=True)
self.conv2 = nn.Conv3d(out_channels, out_channels, kernel_size=3, padding=1, bias=False)
self.norm2 = nn.GroupNorm(g, out_channels)
self.dropout = nn.Dropout3d(dropout) if dropout > 0 else nn.Identity()
if in_channels != out_channels:
self.skip = nn.Conv3d(in_channels, out_channels, kernel_size=1, bias=False)
else:
self.skip = nn.Identity()
self.act2 = nn.LeakyReLU(0.1, inplace=True)
def forward(self, x: torch.Tensor) -> torch.Tensor:
identity = self.skip(x)
out = self.conv1(x)
out = self.norm1(out)
out = self.act1(out)
out = self.dropout(out)
out = self.conv2(out)
out = self.norm2(out)
out = out + identity
out = self.act2(out)
return out
class Conv3dVoxelClassifier(nn.Module):
def __init__(self, ch_in: int):
super().__init__()
self.stem = nn.Sequential(
nn.Conv3d(ch_in, 16, kernel_size=3, padding=1, bias=False),
nn.GroupNorm(4, 16),
nn.LeakyReLU(0.1, inplace=True),
)
self.features = nn.Sequential(
ResidualBlock3D(16, 32, groups=4, dropout=0.05),
ResidualBlock3D(32, 32, groups=4, dropout=0.10),
ResidualBlock3D(32, 32, groups=4, dropout=0.10),
ResidualBlock3D(32, 16, groups=4, dropout=0.05),
)
self.context = nn.Sequential(
nn.Conv3d(16, 16, kernel_size=3, padding=2, dilation=2, bias=False),
nn.GroupNorm(4, 16),
nn.LeakyReLU(0.1, inplace=True),
)
self.head = nn.Sequential(
nn.Conv3d(32, 16, kernel_size=3, padding=1, bias=False),
nn.GroupNorm(4, 16),
nn.LeakyReLU(0.1, inplace=True),
nn.Conv3d(16, 1, kernel_size=1),
)
#linear layer
self.classifier = nn.Sequential(
nn.Linear(1, 8),
nn.LeakyReLU(0.1, inplace=True),
nn.Dropout(0.1),
nn.Linear(8, 1)
)
def forward(self, x: torch.Tensor) -> torch.Tensor:
x = self.stem(x)
feat = self.features(x)
ctx = self.context(feat)
x = torch.cat([feat, ctx], dim=1)
x = self.head(x)
# NOVO: global pooling
x = x.mean(dim=(2, 3, 4)) # (N, 1)
x = self.classifier(x) # (N, 1)
return x
# =========================================================
# HELPERS
# =========================================================
def ensure_coords_shape(coords: torch.Tensor) -> torch.Tensor:
if coords.ndim == 2:
if coords.shape[1] != 3:
raise ValueError(f"coords expected [N,3], got {tuple(coords.shape)}")
coords = coords.unsqueeze(0)
elif coords.ndim == 3:
if coords.shape[-1] != 3:
raise ValueError(f"coords expected [B,N,3], got {tuple(coords.shape)}")
else:
raise ValueError(f"coords must have 2 or 3 dims, got {coords.ndim}")
return coords.float()
def ensure_radii_shape(radii: torch.Tensor) -> torch.Tensor:
if radii.ndim == 1:
radii = radii.unsqueeze(0)
elif radii.ndim == 2:
if radii.shape[0] == 1:
pass
elif radii.shape[1] == 1:
radii = radii.squeeze(1).unsqueeze(0)
else:
raise ValueError(f"Unexpected radii shape: {tuple(radii.shape)}")
elif radii.ndim == 3:
if radii.shape[0] == 1 and radii.shape[2] == 1:
radii = radii.squeeze(-1)
else:
raise ValueError(f"Unexpected radii shape: {tuple(radii.shape)}")
else:
raise ValueError(f"radii must have 1, 2 or 3 dims, got {radii.ndim}")
return radii.float()
def ensure_atomwise_2d(x: torch.Tensor, n_atoms_hint: int | None = None, name: str = "tensor") -> torch.Tensor:
if x.ndim == 1:
x = x.unsqueeze(0)
elif x.ndim == 2:
if x.shape[0] == 1:
pass
elif x.shape[1] == 1:
x = x.squeeze(1).unsqueeze(0)
else:
if n_atoms_hint is not None:
if x.shape[0] == n_atoms_hint:
x = x.unsqueeze(0)
elif x.shape[1] == n_atoms_hint:
pass
else:
raise ValueError(f"Unexpected {name} shape: {tuple(x.shape)}")
else:
raise ValueError(f"Unexpected {name} shape: {tuple(x.shape)}")
elif x.ndim == 3:
if x.shape[0] != 1:
raise ValueError(f"Unexpected {name} batch shape: {tuple(x.shape)}")
if x.shape[2] == 1:
x = x.squeeze(-1)
elif x.shape[1] == 1:
x = x.squeeze(1)
else:
if n_atoms_hint is not None:
if x.shape[1] == n_atoms_hint:
x = x[:, :, 0]
elif x.shape[2] == n_atoms_hint:
x = x[:, 0, :]
else:
raise ValueError(f"Unexpected {name} shape: {tuple(x.shape)}")
else:
raise ValueError(f"Unexpected {name} shape: {tuple(x.shape)}")
else:
raise ValueError(f"{name} must have 1, 2, or 3 dims, got {x.ndim}")
if x.ndim != 2:
raise ValueError(f"{name} could not be normalized to [B,N], got {tuple(x.shape)}")
return x.float()
def align_atomwise_tensors(coords, atom_channels, radii, atom_labels):
n = min(coords.shape[1], atom_channels.shape[1], radii.shape[1], atom_labels.shape[1])
coords = coords[:, :n, :]
atom_channels = atom_channels[:, :n]
radii = radii[:, :n]
atom_labels = atom_labels[:, :n]
return coords, atom_channels, radii, atom_labels
def ensure_volume_shape_for_cnn(volume: torch.Tensor) -> torch.Tensor:
if volume.ndim == 4:
volume = volume.unsqueeze(1)
elif volume.ndim == 5:
if volume.shape[-1] <= 64 and volume.shape[1] > 64:
volume = volume.permute(0, 4, 1, 2, 3).contiguous()
else:
raise ValueError(f"Unexpected volume shape: {tuple(volume.shape)}")
return volume.float()
# =========================================================
# LABELS
# =========================================================
"""
def extract_atom_labels_from_dssp(pdb_file: Path) -> np.ndarray:
with warnings.catch_warnings():
warnings.simplefilter("ignore")
u = mda.Universe(str(pdb_file))
protein = u.select_atoms("protein")
residues = protein.residues
dssp = DSSP(u).run()
res_labels = dssp.results.dssp[0]
n = min(len(residues), len(res_labels))
atom_labels = []
for i in range(n):
residue = residues[i]
ss = res_labels[i]
label = 1.0 if ss == "H" else 0.0
for _atom in residue.atoms:
atom_labels.append(label)
return np.asarray(atom_labels, dtype=np.float32)
"""
def extract_atom_labels_from_dssp(pdb_file: Path) -> np.ndarray:
with warnings.catch_warnings():
warnings.simplefilter("ignore")
u = mda.Universe(str(pdb_file))
protein = u.select_atoms("protein")
residues = protein.residues
dssp = DSSP(u).run()
res_labels = dssp.results.dssp[0]
n = min(len(residues), len(res_labels))
atom_labels = []
for i in range(n):
residue = residues[i]
ss = res_labels[i]
label = 1.0 if ss == "H" else 0.0
#for _atom in residue.atoms:
atom_labels.append(label)
return np.asarray(atom_labels, dtype=np.float32)
def get_pdb_files(pdb_dir: Path, pattern: str) -> list[Path]:
files = sorted(pdb_dir.glob(pattern))
if not files:
raise FileNotFoundError(f"No PDB files found in: {pdb_dir.resolve()}")
return files
# =========================================================
# BUILDERS
# =========================================================
"""
def build_pyuul_input_and_target(pdb_file: Path, device: str = "cpu", resolution: float = 0.5):
coords, atname = utils.parsePDB(str(pdb_file))
atom_channels = utils.atomlistToChannels(atname)
radii = utils.atomlistToRadius(atname)
coords = ensure_coords_shape(coords.to(device))
radii = ensure_radii_shape(radii.to(device))
atom_labels_np = extract_atom_labels_from_dssp(pdb_file)
atom_labels = torch.tensor(atom_labels_np, dtype=torch.float32, device=device)
n_atoms_hint = coords.shape[1]
atom_channels = ensure_atomwise_2d(atom_channels.to(device), n_atoms_hint=n_atoms_hint, name="atom_channels")
#atom_labels = ensure_atomwise_2d(atom_labels, n_atoms_hint=n_atoms_hint, name="atom_labels")
residue_labels_np = extract_atom_labels_from_dssp(pdb_file) # [N]
target_volume = torch.tensor(residue_labels_np, dtype=torch.float32).unsqueeze(1) # [N, 1]
coords, atom_channels, radii, atom_labels = align_atomwise_tensors(coords, atom_channels, radii, atom_labels)
volmaker = VolumeMaker.Voxels(device=device)
input_volume = volmaker(coords, radii, atom_channels, resolution=resolution).to_dense().float()
input_volume = ensure_volume_shape_for_cnn(input_volume)
if input_volume.shape[0] != target_volume.shape[0]:
n = min(input_volume.shape[0], target_volume.shape[0])
input_volume = input_volume[:n]
target_volume = target_volume[:n]
target_volume = (target_volume > 0).float()
# NOVO: reduzir voxel → label global
target_volume = target_volume.amax(dim=(2, 3, 4)) # (N, 1)
return input_volume, target_volume
"""
def build_pyuul_input_and_target(
pdb_file: Path,
device: str = "cpu",
resolution: float = 1.0,
box_size: int = 41,
):
"""
Build per-residue local PyUUL volumes for residue-level classification.
Returns
-------
input_volume : torch.Tensor
Shape [N_res, C, box_size, box_size, box_size]
target : torch.Tensor
Shape [N_res, 1]
"""
# --------------------------------------------------
# Parse structure with PyUUL
# --------------------------------------------------
coords, atname = utils.parsePDB(str(pdb_file))
atom_channels = utils.atomlistToChannels(atname)
radii = utils.atomlistToRadius(atname)
coords = ensure_coords_shape(coords.to(device)) # [1, N_atoms, 3]
radii = ensure_radii_shape(radii.to(device)) # [1, N_atoms]
n_atoms_hint = coords.shape[1]
atom_channels = ensure_atomwise_2d(
atom_channels.to(device),
n_atoms_hint=n_atoms_hint,
name="atom_channels"
) # [1, N_atoms]
# number of atom-type channels expected by the model
n_input_channels = int(atom_channels.max().item()) + 1
# --------------------------------------------------
# Residue-level labels from DSSP
# IMPORTANT: extract_atom_labels_from_dssp must return
# one label per residue for this pipeline
# --------------------------------------------------
residue_labels_np = extract_atom_labels_from_dssp(pdb_file) # [N_res]
target = torch.tensor(
residue_labels_np,
dtype=torch.float32,
device=device
).unsqueeze(1) # [N_res, 1]
# --------------------------------------------------
# Load residues with MDAnalysis
# --------------------------------------------------
with warnings.catch_warnings():
warnings.simplefilter("ignore")
u = mda.Universe(str(pdb_file))
protein = u.select_atoms("protein")
residues = protein.residues
n_res = min(len(residues), target.shape[0])
if n_res == 0:
raise ValueError(f"No valid residues found for {pdb_file}")
# trim target if needed
target = target[:n_res]
volmaker = VolumeMaker.Voxels(device=device)
half_box_ang = (box_size * resolution) / 2.0
local_volumes = []
# global coords tensor without batch dim -> [N_atoms, 3]
global_coords = coords[0]
# --------------------------------------------------
# Build one local volume per residue
# --------------------------------------------------
for i in range(n_res):
res = residues[i]
# Prefer backbone centroid (N, CA, C); fallback to residue center
backbone = res.atoms.select_atoms("name N CA C")
if len(backbone) > 0:
center_np = backbone.positions.mean(axis=0)
else:
center_np = res.atoms.center_of_geometry()
center = torch.tensor(center_np, dtype=torch.float32, device=device)
# Local coordinates centered on the residue
local_coords = global_coords - center # [N_atoms, 3]
# Cubic local crop around the residue center
mask = (
(local_coords[:, 0].abs() <= half_box_ang) &
(local_coords[:, 1].abs() <= half_box_ang) &
(local_coords[:, 2].abs() <= half_box_ang)
)
# Fallback: empty local box with correct channel count
if mask.sum().item() == 0:
local_volume = torch.zeros(
(1, n_input_channels, box_size, box_size, box_size),
dtype=torch.float32,
device=device
)
local_volumes.append(local_volume)
continue
sel_coords = local_coords[mask].unsqueeze(0) # [1, n_sel, 3]
sel_channels = atom_channels[:, mask] # [1, n_sel]
sel_radii = radii[:, mask] # [1, n_sel]
# Voxelize local environment
local_volume = volmaker(
sel_coords,
sel_radii,
sel_channels,
resolution=resolution
).to_dense().float()
# Normalize tensor layout to [1, C, D, H, W]
if local_volume.ndim == 5:
# Sometimes dense output may come as [1, D, H, W, C]
if local_volume.shape[1] != n_input_channels and local_volume.shape[-1] == n_input_channels:
local_volume = local_volume.permute(0, 4, 1, 2, 3).contiguous()
elif local_volume.ndim == 4:
local_volume = local_volume.unsqueeze(1)
else:
raise ValueError(
f"Unexpected local_volume shape for residue {i} in {pdb_file}: "
f"{tuple(local_volume.shape)}"
)
# If channel dimension is still wrong, force a safe fix when possible
if local_volume.shape[1] != n_input_channels:
if local_volume.shape[1] == 1 and local_volume.shape[-1] == n_input_channels:
local_volume = local_volume.permute(0, 4, 1, 2, 3).contiguous()
elif local_volume.shape[1] != n_input_channels:
raise ValueError(
f"Channel mismatch for residue {i} in {pdb_file}: "
f"expected {n_input_channels}, got {local_volume.shape[1]}"
)
# Force fixed spatial size
if local_volume.shape[-3:] != (box_size, box_size, box_size):
local_volume = F.interpolate(
local_volume,
size=(box_size, box_size, box_size),
mode="nearest"
)
local_volumes.append(local_volume)
input_volume = torch.cat(local_volumes, dim=0) # [N_res, C, box, box, box]
# Final safety check
if input_volume.shape[0] != target.shape[0]:
n = min(input_volume.shape[0], target.shape[0])
input_volume = input_volume[:n]
target = target[:n]
return input_volume.float(), target.float()
def build_spheronizator_input_and_target(
pdb_file: Path,
atoms_dir: Path,
bonds_dir: Path,
target_device: str = "cpu",
resolution: float = 0.5,
aggregate: str = "sum",
):
input_volume = load_spheronizator_volume(pdb_file, atoms_dir, bonds_dir) # [N, C, D, H, W]
residue_labels_np = extract_atom_labels_from_dssp(pdb_file) # [N]
target_volume = torch.tensor(residue_labels_np, dtype=torch.float32, device=target_device).unsqueeze(1) # [N, 1]
if input_volume.shape[0] != target_volume.shape[0]:
n = min(input_volume.shape[0], target_volume.shape[0])
input_volume = input_volume[:n]
target_volume = target_volume[:n]
return input_volume, target_volume
# =========================================================
# DATASETS
# =========================================================
class PyuulDataset(Dataset):
def __init__(self, pdb_files: list[Path], resolution: float = 0.5, build_device: str = "cpu"):
self.pdb_files = list(pdb_files)
self.resolution = resolution
self.build_device = build_device
def __len__(self):
return len(self.pdb_files)
def __getitem__(self, idx):
pdb_file = self.pdb_files[idx]
x, y = build_pyuul_input_and_target(pdb_file, device=self.build_device, resolution=self.resolution)
return x, y, pdb_file.name
class SpheronizatorDataset(Dataset):
def __init__(
self,
pdb_files: list[Path],
atoms_dir: Path,
bonds_dir: Path,
resolution: float = 0.5,
target_build_device: str = "cuda",
aggregate: str = "sum",
):
self.pdb_files = list(pdb_files)
self.atoms_dir = Path(atoms_dir)
self.bonds_dir = Path(bonds_dir)
self.resolution = resolution
self.target_build_device = target_build_device
self.aggregate = aggregate
def __len__(self):
return len(self.pdb_files)
def __getitem__(self, idx):
pdb_file = self.pdb_files[idx]
x, y = build_spheronizator_input_and_target(
pdb_file=pdb_file,
atoms_dir=self.atoms_dir,
bonds_dir=self.bonds_dir,
target_device=self.target_build_device,
resolution=self.resolution,
aggregate=self.aggregate,
)
return x, y, pdb_file.name
def single_item_collate(batch):
return batch[0]
def build_loader(backend: str, pdb_files: list[Path], cfg: Config) -> DataLoader:
if backend == "pyuul":
dataset = PyuulDataset(
pdb_files=pdb_files,
resolution=cfg.resolution,
build_device=cfg.target_build_device,
)
elif backend == "spheronizator":
dataset = SpheronizatorDataset(
pdb_files=pdb_files,
atoms_dir=cfg.sphero_atoms_dir,
bonds_dir=cfg.sphero_bonds_dir,
resolution=cfg.resolution,
target_build_device=cfg.target_build_device,
aggregate=cfg.sphero_aggregate,
)
else:
raise ValueError(f"Unknown backend: {backend}")
return DataLoader(
dataset,
batch_size=cfg.batch_size,
shuffle=True,
num_workers=cfg.num_workers,
pin_memory=cfg.pin_memory,
collate_fn=single_item_collate,
)
# =========================================================
# METRICS
# =========================================================
def compute_voxel_metrics(logits, y, threshold=0.5, eps=1e-8):
probs = torch.sigmoid(logits)
preds = (probs > threshold).float()
y = y.float()
tp = ((preds == 1) & (y == 1)).sum().float()
tn = ((preds == 0) & (y == 0)).sum().float()
fp = ((preds == 1) & (y == 0)).sum().float()
fn = ((preds == 0) & (y == 1)).sum().float()
acc = (tp + tn) / (tp + tn + fp + fn + eps)
precision = tp / (tp + fp + eps)
recall = tp / (tp + fn + eps)
f1 = 2 * precision * recall / (precision + recall + eps)
dice = (2 * tp) / (2 * tp + fp + fn + eps)
iou = tp / (tp + fp + fn + eps)
return {
"acc": acc.item(),
"precision": precision.item(),
"recall": recall.item(),
"f1": f1.item(),
"dice": dice.item(),
"iou": iou.item(),
}
def dice_loss_from_logits(logits, targets, eps=1e-6):
probs = torch.sigmoid(logits)
probs = probs.reshape(probs.shape[0], -1)
targets = targets.reshape(targets.shape[0], -1).float()
intersection = (probs * targets).sum(dim=1)
union = probs.sum(dim=1) + targets.sum(dim=1)
dice = (2 * intersection + eps) / (union + eps)
return 1 - dice.mean()
def plot_roc_curve(y_true, y_score, output_path="roc_curve.pdf", title="ROC Curve"):
y_true = np.asarray(y_true).astype(int)
y_score = np.asarray(y_score).astype(float)
if len(np.unique(y_true)) < 2:
print(f"[WARNING] ROC curve skipped for {title}: only one class present.")
return float("nan")
fpr, tpr, _ = roc_curve(y_true, y_score)
roc_auc = auc(fpr, tpr)
plt.figure(figsize=(6, 6))
plt.plot(fpr, tpr, label=f"AUC = {roc_auc:.4f}")
plt.plot([0, 1], [0, 1], linestyle="--")
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
plt.title(title)
plt.legend(loc="lower right")
plt.tight_layout()
plt.savefig(output_path)
plt.close()
return roc_auc
def plot_confusion_matrix(y_true, y_pred, output_path="confusion_matrix.pdf", title="Confusion Matrix"):
y_true = np.asarray(y_true).astype(int)
y_pred = np.asarray(y_pred).astype(int)
#cm = confusion_matrix(y_true, y_pred)
cm = confusion_matrix(y_true, y_pred, labels=[0, 1])
fig, ax = plt.subplots(figsize=(5, 5))
disp = ConfusionMatrixDisplay(confusion_matrix=cm, display_labels=[0, 1])
disp.plot(ax=ax, colorbar=False)
ax.set_title(title)
plt.tight_layout()
plt.savefig(output_path)
plt.close()
return cm
# =========================================================
# TRAIN / EVAL
# =========================================================
def train_voxel_classification(
model,
train_loader,
test_loader,
epochs,
lr,
device,
use_amp,
save_checkpoint_every,
backend,
output_dir,
):
bce = nn.BCEWithLogitsLoss()
optimizer = torch.optim.AdamW(model.parameters(), lr=1e-3, weight_decay=1e-4)
use_amp = use_amp and device.startswith("cuda")
scaler = torch.amp.GradScaler("cuda", enabled=use_amp)
history = {
"train_loss": [],
"test_loss": [],
"accuracy": [],
"precision": [],
"recall": [],
"f1": [],
"dice": [],
"iou": [],
}
for epoch in range(1, epochs + 1):
model.train()
train_loss = 0.0
for x, y, _ in train_loader:
x = x.to(device, non_blocking=True)
y = y.to(device, non_blocking=True).float()
optimizer.zero_grad(set_to_none=True)
with torch.amp.autocast(device_type="cuda", enabled=use_amp):
logits = model(x)
if logits.shape != y.shape:
raise RuntimeError(
f"Shape mismatch: output={tuple(logits.shape)} target={tuple(y.shape)}"
)
loss = bce(logits, y) + dice_loss_from_logits(logits, y)
scaler.scale(loss).backward()
scaler.step(optimizer)
scaler.update()
train_loss += loss.item()
train_loss /= max(len(train_loader), 1)
model.eval()
test_loss = 0.0
metric_sums = {
"acc": 0.0,
"precision": 0.0,
"recall": 0.0,
"f1": 0.0,
"dice": 0.0,
"iou": 0.0,
}
n_test = 0
with torch.no_grad():
for x, y, _ in test_loader:
x = x.to(device, non_blocking=True)
y = y.to(device, non_blocking=True).float()
with torch.amp.autocast(device_type="cuda", enabled=use_amp):
logits = model(x)
loss = bce(logits, y) + dice_loss_from_logits(logits, y)
test_loss += loss.item()
metrics = compute_voxel_metrics(logits, y)
for k, v in metrics.items():
metric_sums[k] += v
n_test += 1
test_loss /= max(n_test, 1)
epoch_metrics = {k: v / max(n_test, 1) for k, v in metric_sums.items()}
history["train_loss"].append(train_loss)
history["test_loss"].append(test_loss)
history["accuracy"].append(epoch_metrics["acc"])
history["precision"].append(epoch_metrics["precision"])
history["recall"].append(epoch_metrics["recall"])
history["f1"].append(epoch_metrics["f1"])
history["dice"].append(epoch_metrics["dice"])
history["iou"].append(epoch_metrics["iou"])
print(
f"[{backend}] Epoch {epoch:02d}/{epochs:02d} | "
f"train_loss={train_loss:.4f} | "
f"test_loss={test_loss:.4f} | "
f"acc={epoch_metrics['acc']:.4f} | "
f"prec={epoch_metrics['precision']:.4f} | "
f"rec={epoch_metrics['recall']:.4f} | "
f"f1={epoch_metrics['f1']:.4f} | "
f"dice={epoch_metrics['dice']:.4f} | "
f"iou={epoch_metrics['iou']:.4f}"
)
if save_checkpoint_every > 0 and epoch % save_checkpoint_every == 0:
torch.save(
{
"epoch": epoch,
"model_state_dict": model.state_dict(),
"optimizer_state_dict": optimizer.state_dict(),
"train_loss": train_loss,
"test_loss": test_loss,
"metrics": epoch_metrics,
"backend": backend,
},
output_dir / f"checkpoint_voxel_{backend}_epoch_{epoch:02d}.pt",
)
return history
# =========================================================
# PLOTS
# =========================================================
def plot_sample_voxels(input_volume: torch.Tensor, target_volume: torch.Tensor, title: str = "", output_dir: Path | None = None):
x = input_volume.detach().cpu()
y = target_volume.detach().cpu()
input_mask = (x[0].sum(dim=0) > 0)
target_mask = (y[0, 0] > 0)
fig = plt.figure(figsize=(12, 5))
ax1 = fig.add_subplot(1, 2, 1, projection="3d")
ax1.voxels(input_mask.numpy())
ax1.set_title(f"{title} - Input")
ax2 = fig.add_subplot(1, 2, 2, projection="3d")
ax2.voxels(target_mask.numpy())
ax2.set_title(f"{title} - Target")
plt.tight_layout()
out = (output_dir / f"{title}_sample_voxels.png") if output_dir else Path(f"{title}_sample_voxels.png")
plt.savefig(out, dpi=200)
plt.close(fig)
def plot_classification_history(history: dict, prefix: str, output_dir: Path):
plt.figure(figsize=(7, 4))
plt.plot(history["train_loss"], marker="o", label="train_loss")
plt.plot(history["test_loss"], marker="o", label="test_loss")
plt.xlabel("Epoch")
plt.ylabel("Loss")
plt.title(f"Classification Loss - {prefix}")
plt.legend()
plt.tight_layout()
plt.savefig(output_dir / f"{prefix}_cls_loss.pdf")
plt.close()
metrics_to_plot = ["accuracy", "precision", "recall", "f1"]
plt.figure(figsize=(10, 5))
for metric in metrics_to_plot:
plt.plot(history[metric], marker="o", label=metric)
plt.xlabel("Epoch")
plt.ylabel("Score")
plt.title(f"Classification Metrics - {prefix}")
plt.legend()
plt.tight_layout()
plt.savefig(output_dir / f"{prefix}_cls_metrics.pdf")
plt.close()
def plot_backend_comparison(results_df: pd.DataFrame, output_path: Path):
metrics = ["accuracy", "precision", "recall", "f1"]
x = np.arange(len(matrics))
width = 0.12
offsets = np.linspace(-2.5 * width, 2.5 * width, len(metrics))
plt.figure(figsize=(12, 6))
for metric, offset in zip(metrics, offsets):
plt.bar(x + offset, results_df[metric].values, width=width, label=metric)
plt.xticks(x, backends)
plt.ylim(0, 1)
plt.ylabel("Score")
plt.title("Backend Comparison - Final Test Metrics")
plt.legend()
plt.tight_layout()
plt.savefig(output_path)
plt.close()
# =========================================================
# EXPERIMENT RUNNER
# =========================================================
def run_backend_classification(cfg: Config, backend: str, train_files: list[Path], test_files: list[Path]):
out_dir = backend_output_dir(cfg, backend)
train_loader = build_loader(backend, train_files, cfg)
test_loader = build_loader(backend, test_files, cfg)
x0, y0, name0 = next(iter(train_loader))
print(f"\n=== Backend: {backend} ===")
print("First sample:")
print(" name :", name0)
print(" x :", tuple(x0.shape))
print(" y :", tuple(y0.shape))
if cfg.plot_first_sample:
plot_sample_voxels(x0, y0, title=name0, output_dir=out_dir)
ch_in = x0.shape[1]
model = Conv3dVoxelClassifier(ch_in=ch_in).to(cfg.device)
history = train_voxel_classification(
model=model,
train_loader=train_loader,
test_loader=test_loader,
epochs=cfg.epochs,
lr=cfg.learning_rate,
device=cfg.device,
use_amp=cfg.use_amp,
save_checkpoint_every=cfg.save_checkpoint_every,
backend=backend,
output_dir=out_dir)
plot_classification_history(history, prefix=backend, output_dir=out_dir)
model_path = out_dir / f"conv3d_classifier_{backend}.pt"
torch.save(model.state_dict(), model_path)
print(f"Model saved to: {model_path}")
with open(out_dir / "config.json", "w") as f:
json.dump(
{k: str(v) if isinstance(v, Path) else v for k, v in asdict(cfg).items()},
f,
indent=2,
)
return history
# =========================================================
# MAIN
# =========================================================
def main():
cfg = build_config()
print("Device:", cfg.device)
print("Backend mode:", cfg.backend)
print("Resolution:", cfg.resolution)
print("Output dir:", cfg.output_dir)
pdb_files = get_pdb_files(cfg.pdb_dir, cfg.pdb_pattern)
print(f"Found {len(pdb_files)} PDB files")
train_files, test_files = train_test_split(
pdb_files,
test_size=cfg.test_size,
random_state=cfg.random_state,
shuffle=True,
)
if cfg.backend == "compare":
summaries = []
for backend in ["pyuul","spheronizator"]:
summary = run_backend_classification(cfg, backend, train_files, test_files)
summaries.append(summary)
results_df = pd.DataFrame(summaries)
results_df.to_csv(cfg.output_dir / "backend_comparison.csv", index=False)
plot_backend_comparison(
results_df,
output_path=cfg.output_dir / "backend_comparison.pdf",
)
print("\n=== Final Backend Comparison ===")
print(results_df.to_string(index=False))
else:
summary = run_backend_classification(cfg, cfg.backend, train_files, test_files)
summary_df = pd.DataFrame([summary])
summary_df.to_csv(cfg.output_dir / f"{cfg.backend}_summary.csv", index=False)
print("\n=== Final Summary ===")
print(summary_df.to_string(index=False))
if __name__ == "__main__":
main()
# In[ ]:
# python train_classifier.py --backend compare
# python train_classifier.py --backend spheronizator
# python train_classifier.py --backend compare

3DCNN-SpheronizaTor.py at commit 36b80d6, under BSD-3-Clause · at the source

Overview

Authors: Jose Cleydson Ferreira Silva¹, Matthew Richardson², José D. D. Cediel-Becerra¹, Layla Schuster¹, Raquel Dias¹

ORCID iDs: Jose Cleydson Ferreira Silva, José D. D. Cediel-Becerra, Layla Schuster, Raquel Dias

Institutions: University of Florida (United States); Institute of Food and Agricultural Sciences (United States)

Journal: Computational and structural biotechnology journal, volume 35, issue 1, article 0076

Dates: received 17 December 2025; accepted 9 April 2026; published online 7 May 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.34133/csbj.0076 · PMID 42110217 · PMCID PMC13150068 · OpenAlex W7153140616

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: cellular / molecular (subfield)

Methods: Machine learning

Journal subjects: Software/Web Server Article, General

Topic: Cell Image Analysis Techniques (Biophysics, Biochemistry, Genetics and Molecular Biology), according to OpenAlex

Funding: University of Florida

Citations: not cited yet (Europe PMC); 32 references in the paper

Abstract

Artificial intelligence (AI) has expanded the reach of structural biology by enabling models to extract biochemical and geometric features directly from 3-dimensional (3D) protein structures. Yet, the effectiveness of these models depends critically on how protein environments are encoded. Most existing volumetric representations rely on Cartesian voxel grids derived from smoothed atomic densities, an approach that offers broad applicability but struggles to reconcile rotational invariance, residue-level specificity, and explicit biochemical detail. We present SpheronizaTor, a residue-centered voxelization framework for protein structures that builds local spherical voxel maps centered on each residue. Each spherical map encodes atom types, covalent bonding information, and whether atoms belong to the central residue or neighboring residues. By producing one voxel representation per residue, SpheronizaTor emphasizes the structural and functional granularity through which proteins organize catalysis, recognition, and stability. The combination of spherical alignment with chemically explicit feature channels enables richer interpretability and enhances compatibility with 3D convolutional and hybrid neural architectures. Designed specifically for proteins and engineered for extensibility, SpheronizaTor provides a voxelization strategy that is both chemically realistic and computationally efficient. The residue-centric approach bridges the gap between global volumetric encoders and graph-based models, offering a versatile foundation for downstream tasks such as mutation effect prediction, binding site analysis, and structural comparison across protein families.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 5 matches between paragraphs and lines of code.

Dias-Lab/spheronizator

License: BSD-3-Clause

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: 36b80d66e065474d30036a931dff43bc4ad83291, 8 April 2026

Languages: Python (9), Jupyter (1)

Size: 45 files, 10 scripts

Software Heritage: not archived

Found in: “Data Availability”

Holds: README, license file, environment (pyproject.toml, setup.py), tests, continuous integration, documentation, 1 notebook

Not found: CITATION.cff

Tools: NumPy (7 files), Biopython (5 files), pandas (2 files), Matplotlib (1 file), MDAnalysis (1 file), PyTorch (1 file), scikit-learn (1 file)

Availability: 1 check, the latest on 28 September 2026: the link answers

28 September 2026: the link answers

12 files

examples/3DCNN-SpheronizaTor.py, Python, 1,095 lines, 2 matches
examples/demo.ipynb, Jupyter, 227 lines
setup.py, Python, 17 lines
src/spheronizator/__init__.py, Python, 2 lines
src/spheronizator/functions.py, Python, 253 lines, 1 match
src/spheronizator/mol2parser.py, Python, 142 lines, 1 match
src/spheronizator/voxelBuilder.py, Python, 231 lines, 1 match
src/spheronizator/voxelize.py, Python, 152 lines
tests/dataWrapper.py, Python, 46 lines
tests/debug.exportBoxes.py, Python, 58 lines
LICENSE, License, 28 lines
README.md, Text, 77 lines

At the source: github.com/Dias-Lab/spheronizator

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
10 scripts, each with its path and the digest of its content;
5 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability

The SpheronizaTor software is freely available on GitHub at https://github.com/Dias-Lab/spheronizator and can also be installed via pip (https://pypi.org/project/spheronizator/). Detailed documentation is available at https://dias-lab.github.io/spheronizator/.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 1 funder, 26 references.

Cite

This paper

Silva, J. C. F., Richardson, M., Cediel-Becerra, J. D. D., Schuster, L., & Dias, R. (2026). SpheronizaTor: Spherical Voxelization for Interpretable Protein Microenvironment Modeling. Computational and structural biotechnology journal, 35(1), 0076. https://doi.org/10.34133/csbj.0076

BibTeX

@article{silva2026spheronizator, author = {Silva, Jose Cleydson Ferreira and Richardson, Matthew and Cediel-Becerra, José D. D. and Schuster, Layla and Dias, Raquel}, title = {{SpheronizaTor: Spherical Voxelization for Interpretable Protein Microenvironment Modeling}}, journal = {Computational and structural biotechnology journal}, year = {2026}, month = may, volume = {35}, number = {1}, pages = {0076}, publisher = {AAAS Science Partner Journal Program}, issn = {2001-0370}, doi = {10.34133/csbj.0076}, url = {https://doi.org/10.34133/csbj.0076}, pmid = {42110217}, pmcid = {PMC13150068} }

RIS

TY - JOUR AU - Silva, Jose Cleydson Ferreira AU - Richardson, Matthew AU - Cediel-Becerra, José D. D. AU - Schuster, Layla AU - Dias, Raquel TI - SpheronizaTor: Spherical Voxelization for Interpretable Protein Microenvironment Modeling T2 - Computational and structural biotechnology journal J2 - Comput Struct Biotechnol J PY - 2026 DA - 2026/05/07 VL - 35 IS - 1 SP - 0076 SN - 2001-0370 PB - AAAS Science Partner Journal Program DO - 10.34133/csbj.0076 UR - https://doi.org/10.34133/csbj.0076 LA - en ER -

CSL-JSON

{ "id": "10.34133/csbj.0076", "type": "article-journal", "title": "SpheronizaTor: Spherical Voxelization for Interpretable Protein Microenvironment Modeling", "container-title": "Computational and structural biotechnology journal", "author": [ { "family": "Silva", "given": "Jose Cleydson Ferreira" }, { "family": "Richardson", "given": "Matthew" }, { "family": "Cediel-Becerra", "given": "José D. D." }, { "family": "Schuster", "given": "Layla" }, { "family": "Dias", "given": "Raquel" } ], "container-title-short": "Comput Struct Biotechnol J", "volume": "35", "issue": "1", "page": "0076", "DOI": "10.34133/csbj.0076", "PMID": "42110217", "PMCID": "PMC13150068", "ISSN": "2001-0370", "publisher": "AAAS Science Partner Journal Program", "URL": "https://doi.org/10.34133/csbj.0076", "language": "en", "issued": { "date-parts": [ [ 2026, 5, 7 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1016/j.molcel.2026.07.006 [code]: DeorphaNN: Virtual screening of GPCR peptide agonists using AlphaFold-predicted active-state complexes and deep learning embeddings.
Journal: Molecular cell
In common: Biopython, PyTorch, scikit-learn, 3 other tools, cellular / molecular, 2 references
[2] doi:10.1038/s41467-026-75444-4 [code]: Structural insights enable drug discovery for the neuronal NBCn2 carbonate transporter.
Journal: Nature communications
In common: MDAnalysis, Biopython, pandas, 2 other tools, cellular / molecular
[3] doi:10.1002/advs.202523984 [code]: INB3P: A Multi-Modal and Interpretable Co-Attention Framework Integrating Property-Aware Explanations and Memory-Bank Contrastive Fusion for Blood-Brain Barrier Penetrating Peptide Discovery.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: Biopython, PyTorch, scikit-learn, 3 other tools, 1 reference
[4] doi:10.1038/s41586-026-10658-6 [code]: An AI system to help scientists write expert-level empirical software.
Journal: Nature
In common: PyTorch, scikit-learn, pandas, 2 other tools, 2 references
[5] doi:10.1038/s42003-026-10957-8 [code]: Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: Biopython, PyTorch, scikit-learn, 3 other tools, cellular / molecular
[6] doi:10.3390/ijms27156614 [code]: Candidalysin Inhibits Porphyromonas gingivalis Lipoprotein-Induced IL-1β Production in BV-2 Microglia via Hydrophobic Microbial Interactions.
Journal: International journal of molecular sciences
In common: Biopython, PyTorch, scikit-learn, 3 other tools, cellular / molecular
[7] doi:10.1038/s41598-026-53415-5 [code]: Computational design and immunoinformatics validation of a T cell multi-epitope vaccine targeting glioblastoma stem cells.
Journal: Scientific reports
In common: Biopython, PyTorch, scikit-learn, 3 other tools, cellular / molecular
[8] doi:10.1038/s41586-026-10391-0 [code]: Cell-type-targeted mitochondrial transplantation rescues cell degeneration.
Journal: Nature
In common: Biopython, PyTorch, scikit-learn, 3 other tools, cellular / molecular
[9] doi:10.1016/j.xops.2026.101249 [code]: A Disorder-Aware Computational Framework to Identify Structurally Tractable Targets in Proliferative Vitreoretinopathy.
Journal: Ophthalmology science
In common: Biopython, PyTorch, pandas, 2 other tools, 1 reference
[10] doi:10.1523/eneuro.0362-25.2026 [code]: Similarities between Ciona Dorsal Motor Ganglion and Vertebrate Cerebellum: Did a Chordate Ancestor Already Show D/V Subdivision within a Hindbrain Precursor?
Journal: eNeuro
In common: Biopython, PyTorch, scikit-learn, 3 other tools

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.