diff --git a/backend/main/views_helper.py b/backend/main/views_helper.py index 281799b94..b63ed73ac 100644 --- a/backend/main/views_helper.py +++ b/backend/main/views_helper.py @@ -5,6 +5,10 @@ import pandas as pd from backend.protzilla.constants.paths import SETTINGS_PATH +from backend.protzilla.data_analysis.amino_acid_spheres import ( + calculate_amino_acid_spheres, + calculate_ptm_spheres, +) from backend.protzilla.disk_operator import YamlOperator from backend.protzilla.steps import Step from backend.protzilla.step_manager import StepManager @@ -211,6 +215,8 @@ def create_visualization( cif_df: pd.DataFrame, structure_entry_id: str, crosslinking_df: pd.DataFrame | None = None, + include_ptm_spheres: bool = False, + ignored_neighbors: int = 0, ) -> dict: """ Create visualization data, by packaging a mmCIF string (converted from a CIF DataFrame) with its structure entry ID. @@ -233,6 +239,14 @@ def create_visualization( result = {"structureEntryId": structure_entry_id, "cifString": cif_string} + if include_ptm_spheres: + result["trimeshMeshes"] = calculate_amino_acid_spheres( + cif_df, + only_intersecting_ptms=True, + ignored_neighbors=int(ignored_neighbors), + ) + result["trimeshMeshes"].extend(calculate_ptm_spheres(cif_df)) + if crosslinking_df is not None: result["crosslinks"] = extract_relevant_crosslink_information(crosslinking_df) diff --git a/backend/protzilla/all_steps.py b/backend/protzilla/all_steps.py index b5717dc1e..d36c3f12e 100644 --- a/backend/protzilla/all_steps.py +++ b/backend/protzilla/all_steps.py @@ -82,6 +82,7 @@ data_analysis.PTMDetailsVisualization, data_analysis.CrosslinkingValidationWithAngstromDeviation, data_analysis.CrosslinkingValidationWithAngstromDeviationForMultimer, + data_analysis.PtmValidation, data_preprocessing.ImputationByMinPerSample, data_integration.EnrichmentAnalysisGOAnalysisWithString, data_integration.EnrichmentAnalysisGOAnalysisWithEnrichr, diff --git a/backend/protzilla/constants/cif_columns.py b/backend/protzilla/constants/cif_columns.py index 5715f2f32..2d0f65afd 100644 --- a/backend/protzilla/constants/cif_columns.py +++ b/backend/protzilla/constants/cif_columns.py @@ -30,6 +30,7 @@ class ATOM_SITE_COLUMNS(StrEnum): AUTH_ASYM_ID = f"{ATOM_SITE_PREFIX}auth_asym_id" AUTH_ATOM_ID = f"{ATOM_SITE_PREFIX}auth_atom_id" PDBX_PDB_MODEL_NUM = f"{ATOM_SITE_PREFIX}pdbx_PDB_model_num" + PDBX_SIFTS_XREF_DB_ACC = f"{ATOM_SITE_PREFIX}pdbx_sifts_xref_db_acc" ATOM_SITE_LABEL_COMP_ID = ATOM_SITE_COLUMNS.LABEL_COMP_ID diff --git a/backend/protzilla/constants/van_der_waals.py b/backend/protzilla/constants/van_der_waals.py new file mode 100644 index 000000000..eb3519e6b --- /dev/null +++ b/backend/protzilla/constants/van_der_waals.py @@ -0,0 +1,123 @@ +vdw_radii = { + """ + Van der Waals radius in Ångström + """ + "H": 1.20, + "He": 1.43, + "Li": 2.12, + "Be": 1.98, + "B": 1.91, + "C": 1.77, + "N": 1.66, + "O": 1.50, + "F": 1.46, + "Ne": 1.58, + "Na": 2.50, + "Mg": 2.51, + "Al": 2.25, + "Si": 2.19, + "P": 1.90, + "S": 1.89, + "Cl": 1.82, + "Ar": 1.83, + "K": 2.73, + "Ca": 2.62, + "Sc": 2.58, + "Ti": 2.46, + "V": 2.42, + "Cr": 2.45, + "Mn": 2.45, + "Fe": 2.44, + "Co": 2.40, + "Ni": 2.40, + "Cu": 2.38, + "Zn": 2.39, + "Ga": 2.32, + "Ge": 2.29, + "As": 1.88, + "Se": 1.82, + "Br": 1.86, + "Kr": 2.25, + "Rb": 3.21, + "Sr": 2.84, + "Y": 2.75, + "Zr": 2.52, + "Nb": 2.56, + "Mo": 2.45, + "Tc": 2.44, + "Ru": 2.46, + "Rh": 2.44, + "Pd": 2.15, + "Ag": 2.53, + "Cd": 2.49, + "In": 2.43, + "Sn": 2.42, + "Sb": 2.47, + "Te": 1.99, + "I": 2.04, + "Xe": 2.06, + "Cs": 3.48, + "Ba": 3.03, + "La": 2.98, + "Ce": 2.88, + "Pr": 2.92, + "Nd": 2.95, + "Pm": None, + "Sm": 2.90, + "Eu": 2.87, + "Gd": 2.83, + "Tb": 2.79, + "Dy": 2.87, + "Ho": 2.81, + "Er": 2.83, + "Tm": 2.79, + "Yb": 2.80, + "Lu": 2.74, + "Hf": 2.63, + "Ta": 2.53, + "W": 2.57, + "Re": 2.49, + "Os": 2.48, + "Ir": 2.41, + "Pt": 2.29, + "Au": 2.32, + "Hg": 2.45, + "Tl": 2.47, + "Pb": 2.60, + "Bi": 2.54, + "Po": None, + "At": None, + "Rn": None, + "Fr": None, + "Ra": None, + "Ac": 2.80, + "Th": 2.93, + "Pa": 2.88, + "U": 2.71, + "Np": 2.82, + "Pu": 2.81, + "Am": 2.83, + "Cm": 3.05, + "Bk": 3.40, + "Cf": 3.05, + "Es": 2.70, + "Fm": None, + "Md": None, + "No": None, + "Lr": None, + "Rf": None, + "Db": None, + "Sg": None, + "Bh": None, + "Hs": None, + "Mt": None, + "Ds": None, + "Rg": None, + "Cn": None, + "Nh": None, + "Fl": None, + "Mc": None, + "Lv": None, + "Ts": None, + "Og": None, +} diff --git a/backend/protzilla/data_analysis/amino_acid_spheres.py b/backend/protzilla/data_analysis/amino_acid_spheres.py new file mode 100644 index 000000000..213d2bb43 --- /dev/null +++ b/backend/protzilla/data_analysis/amino_acid_spheres.py @@ -0,0 +1,194 @@ +from __future__ import annotations + +import numpy as np +import pandas as pd +import trimesh + +from backend.protzilla.constants.cif_columns import ( + ATOM_SITE_COLUMNS, + CHEM_COMP_COLUMNS, +) +from backend.protzilla.data_analysis.geometry_operations import ( + calculate_center_point, + extract_points_from_cif, + find_farthest_point_vdw, + find_intersecting_spheres, + mesh_to_polyhedron, + resolve_chain_column, +) +from backend.protzilla.utilities.ptm_helper import ( + get_all_ptm_atoms_with_coordinates, + get_center_points_and_radius_for_each_ptm, +) + + +def _calculate_residue_spheres( + cif_df: pd.DataFrame, chain_id: str | None = None +) -> list[dict]: + residue_df = cif_df + chain_column = resolve_chain_column(cif_df) + if chain_id is not None and chain_column is not None: + residue_df = cif_df[cif_df[chain_column] == chain_id] + + residue_positions = ( + pd.to_numeric(residue_df[ATOM_SITE_COLUMNS.LABEL_SEQ_ID], errors="coerce") + .dropna() + .astype(int) + .drop_duplicates() + .sort_values() + ) + + spheres = [] + for residue_position in residue_positions: + residue_points, residue_elements = extract_points_from_cif( + cif_df, + residue_range=(residue_position, residue_position), + chain_id=chain_id, + ) + center = calculate_center_point(residue_points) + _, radius = find_farthest_point_vdw(residue_points, residue_elements, center) + + spheres.append( + { + "chain": chain_id, + "position": residue_position, + "center": center.tolist(), + "radius": float(radius), + } + ) + + return spheres + + +def calculate_amino_acid_spheres( + cif_df: pd.DataFrame, + chain_id: str | None = None, + color: int = 0xFF8C00, + alpha: float = 0.25, + subdivisions: int = 1, + only_intersecting_ptms: bool = False, + ignored_neighbors: int = 0, +) -> list[dict]: + if only_intersecting_ptms: + residue_spheres = [] + seen_residues = set() + for collision in find_ptm_amino_acid_sphere_collisions( + cif_df, ignored_neighbors=ignored_neighbors, chain_id=chain_id + ): + for residue_sphere in collision["collisions"]: + residue_key = (residue_sphere["chain"], residue_sphere["position"]) + if residue_key not in seen_residues: + seen_residues.add(residue_key) + residue_spheres.append(residue_sphere) + else: + residue_spheres = _calculate_residue_spheres(cif_df, chain_id) + + spheres = [] + for residue_sphere in residue_spheres: + sphere = trimesh.creation.icosphere( + subdivisions=subdivisions, radius=residue_sphere["radius"] + ) + sphere.apply_translation(residue_sphere["center"]) + + spheres.append( + { + "label": f"Residue {residue_sphere['position']} sphere", + "mesh": mesh_to_polyhedron(sphere), + "color": color, + "alpha": alpha, + } + ) + + return spheres + + +def calculate_ptm_spheres( + cif_df: pd.DataFrame, + color: int = 0x00A6A6, + alpha: float = 0.35, + subdivisions: int = 1, +) -> list[dict]: + if CHEM_COMP_COLUMNS.MON_NSTD_FLAG not in cif_df.columns: + return [] + + ptms = get_center_points_and_radius_for_each_ptm( + get_all_ptm_atoms_with_coordinates(cif_df) + ) + + spheres = [] + for ptm in ptms: + center = np.array(ptm["center_point"], dtype=float) + radius = float(ptm["radius"]) + + sphere = trimesh.creation.icosphere(subdivisions=subdivisions, radius=radius) + sphere.apply_translation(center) + + spheres.append( + { + "label": f"PTM {ptm['ptm_name']} {ptm['chain']}:{ptm['position']} sphere", + "mesh": mesh_to_polyhedron(sphere), + "color": color, + "alpha": alpha, + } + ) + + return spheres + + +def find_ptm_amino_acid_sphere_collisions( + cif_df: pd.DataFrame, + ignored_neighbors: int = 0, + chain_id: str | None = None, +) -> list[dict]: + if CHEM_COMP_COLUMNS.MON_NSTD_FLAG not in cif_df.columns: + return [] + + amino_acid_df = cif_df[cif_df[CHEM_COMP_COLUMNS.MON_NSTD_FLAG] == True] + chain_column = resolve_chain_column(amino_acid_df) + if chain_id is not None: + chains = [chain_id] + elif chain_column is not None: + chains = amino_acid_df[chain_column].dropna().drop_duplicates() + else: + chains = [None] + amino_acid_spheres = [ + sphere + for chain in chains + for sphere in _calculate_residue_spheres(amino_acid_df, chain) + ] + + ptms = get_center_points_and_radius_for_each_ptm( + get_all_ptm_atoms_with_coordinates(cif_df) + ) + if chain_id is not None: + ptms = [ptm for ptm in ptms if ptm["chain"] == chain_id] + + collisions = [] + for ptm in ptms: + ptm_sphere = { + "ptm_name": ptm["ptm_name"], + "chain": ptm["chain"], + "position": int(ptm["position"]), + "center": ptm["center_point"], + "radius": float(ptm["radius"]), + } + checked_spheres = [ + sphere + for sphere in amino_acid_spheres + if not ( + sphere["chain"] == ptm_sphere["chain"] + and 0 + < abs(sphere["position"] - ptm_sphere["position"]) + <= ignored_neighbors + ) + ] + intersecting_spheres = find_intersecting_spheres(checked_spheres, ptm_sphere) + + collisions.append( + { + **ptm_sphere, + "collisions": intersecting_spheres, + } + ) + + return collisions diff --git a/backend/protzilla/data_analysis/crosslinking_validation.py b/backend/protzilla/data_analysis/crosslinking_validation.py index 8592e89f9..9b5af1401 100644 --- a/backend/protzilla/data_analysis/crosslinking_validation.py +++ b/backend/protzilla/data_analysis/crosslinking_validation.py @@ -36,6 +36,7 @@ PLOT_PRIMARY_COLOR, PLOT_SECONDARY_COLOR, ) +from backend.protzilla.constants.cif_columns import ATOM_SITE_COLUMNS def get_reactive_atom_of_amino_acid_residue(amino_acid_type: str) -> str: @@ -72,14 +73,14 @@ def get_coordinates_of_atom_crosslinker_bound_to( """ relevant_atom = get_reactive_atom_of_amino_acid_residue(amino_acid_type) - seq_ids = pd.to_numeric(cif_df["_atom_site.label_seq_id"], errors="coerce") + seq_ids = pd.to_numeric(cif_df[ATOM_SITE_COLUMNS.LABEL_SEQ_ID], errors="coerce") # Filter to the exact reactive atom of the amino acid residue # where the crosslinker is bound (e.g. CA at position 45) cif_df = cif_df[ - (cif_df["_atom_site.label_atom_id"] == relevant_atom) + (cif_df[ATOM_SITE_COLUMNS.LABEL_ATOM_ID] == relevant_atom) & (seq_ids == amino_acid_position_where_crosslinker_bound) - & (cif_df["_atom_site.auth_asym_id"] == chain_id) + & (cif_df[ATOM_SITE_COLUMNS.AUTH_ASYM_ID] == chain_id) ] if cif_df.empty: @@ -89,9 +90,9 @@ def get_coordinates_of_atom_crosslinker_bound_to( row = cif_df.iloc[0] - x = float(row["_atom_site.Cartn_x"]) - y = float(row["_atom_site.Cartn_y"]) - z = float(row["_atom_site.Cartn_z"]) + x = float(row[ATOM_SITE_COLUMNS.CARTN_X]) + y = float(row[ATOM_SITE_COLUMNS.CARTN_Y]) + z = float(row[ATOM_SITE_COLUMNS.CARTN_Z]) return x, y, z @@ -310,7 +311,7 @@ def get_chains( return [] target_ids_as_strings = [str(i) for i in target_ids] relevant_df = cif_df[cif_df[id_column_name].astype(str).isin(target_ids_as_strings)] - chain_ids = relevant_df["_atom_site.auth_asym_id"].dropna().unique().tolist() + chain_ids = relevant_df[ATOM_SITE_COLUMNS.AUTH_ASYM_ID].dropna().unique().tolist() return chain_ids @@ -402,7 +403,7 @@ def monomer_validation( pae_matrix=pae_matrix, plddt_df=plddt_df, valid_ids=valid_ids, - id_column_name="_atom_site.pdbx_sifts_xref_db_acc", + id_column_name=ATOM_SITE_COLUMNS.PDBX_SIFTS_XREF_DB_ACC, structures_to_validate=[protein_id], validation_criterion=validation_criterion, ) @@ -547,7 +548,7 @@ def multimer_validation( cif_df=cif_df, amino_acid_sequences_df=amino_acid_sequences_df, valid_ids=valid_ids, - id_column_name="_atom_site.label_entity_id", + id_column_name=ATOM_SITE_COLUMNS.LABEL_ENTITY_ID, structures_to_validate=structures_to_validate, pae_matrix=pae_matrix, plddt_df=plddt_df, diff --git a/backend/protzilla/data_analysis/geometry_operations.py b/backend/protzilla/data_analysis/geometry_operations.py new file mode 100644 index 000000000..1ef9bb5d8 --- /dev/null +++ b/backend/protzilla/data_analysis/geometry_operations.py @@ -0,0 +1,263 @@ +from __future__ import annotations + +import numpy as np +import pandas as pd +import trimesh + +from backend.protzilla.constants.cif_columns import ATOM_SITE_COLUMNS +from backend.protzilla.constants.van_der_waals import vdw_radii + +COORDINATE_COLUMNS = [ + ATOM_SITE_COLUMNS.CARTN_X, + ATOM_SITE_COLUMNS.CARTN_Y, + ATOM_SITE_COLUMNS.CARTN_Z, +] + + +def resolve_chain_column(cif_df: pd.DataFrame) -> str | None: + """ + Return the preferred chain identifier column if present in the CIF DataFrame. + + :param cif_df: DataFrame containing mmCIF atom_site data + :return: label chain column, author chain column as fallback, or None + """ + + if ATOM_SITE_COLUMNS.LABEL_ASYM_ID in cif_df.columns: + return ATOM_SITE_COLUMNS.LABEL_ASYM_ID + if ATOM_SITE_COLUMNS.AUTH_ASYM_ID in cif_df.columns: + return ATOM_SITE_COLUMNS.AUTH_ASYM_ID + return None + + +def extract_points_from_cif( + cif_df: pd.DataFrame, + residue_range: tuple[int, int], + chain_id: str | None = None, +) -> tuple[np.ndarray, np.ndarray]: + """ + Extract Cartesian atom coordinates in a residue range from a CIF DataFrame. + + :param cif_df: DataFrame containing mmCIF atom_site coordinates + :param residue_range: inclusive residue interval as [start, end] + :param chain_id: optional chain identifier. If given, the function will filter by + _atom_site.label_asym_id or _atom_site.auth_asym_id when available. + :return: array of n points und x, y, z coordinates. (n, 3) + :raises ValueError: if required columns are missing or too few points remain + """ + + required_columns = { + ATOM_SITE_COLUMNS.LABEL_SEQ_ID, + ATOM_SITE_COLUMNS.TYPE_SYMBOL, + *COORDINATE_COLUMNS, + } + missing_columns = sorted(required_columns - set(cif_df.columns)) + if missing_columns: + raise ValueError( + f"CIF DataFrame is missing required columns for 3D extraction: {missing_columns}" + ) + + start, end = residue_range + if start > end: + raise ValueError( + f"Invalid residue range {residue_range}. Start must be smaller than or equal to end." + ) + + filtered_df = cif_df.copy() + + chain_column = resolve_chain_column(filtered_df) + if chain_id is not None: + if chain_column is None: + raise ValueError( + "A chain_id was provided, but the CIF DataFrame has no chain identifier column." + ) + filtered_df = filtered_df[filtered_df[chain_column] == chain_id] + + residue_ids = filtered_df[ATOM_SITE_COLUMNS.LABEL_SEQ_ID].astype(int) + filtered_df = filtered_df[(residue_ids >= start) & (residue_ids <= end)] + + atom_data = filtered_df[ + [*COORDINATE_COLUMNS, ATOM_SITE_COLUMNS.TYPE_SYMBOL] + ].drop_duplicates() + + points = atom_data[COORDINATE_COLUMNS].astype(float).to_numpy() + elements = ( + atom_data[ATOM_SITE_COLUMNS.TYPE_SYMBOL] + .astype(str) + .str.strip() + .str.capitalize() + .to_numpy() + ) + + if len(points) == 0: + raise ValueError("No atom coordinates found.") + + return points, elements + + +def build_convex_hull(points: np.ndarray) -> trimesh.Trimesh: + """ + Build a 3D convex hull mesh from a point cloud. + + :param points: point coordinates with shape (n, 3) + :return: convex hull + :raises ValueError: if the param points have an invalid shape or contain fewer than + four points + """ + + if points.ndim != 2 or points.shape[1] != 3: + raise ValueError( + f"Expected points with shape (n, 3), got array with shape {points.shape}." + ) + + if len(points) < 4: + raise ValueError( + "At least four distinct points are required to build a 3D convex hull." + ) + + return trimesh.convex.convex_hull(points, qhull_options="QJ") + + +def mesh_to_polyhedron(mesh: trimesh.Trimesh) -> dict: + return { + "vertices": mesh.vertices.tolist(), + "faces": mesh.faces.tolist(), + } + + +def point_cloud_to_polyhedron(points: np.ndarray) -> dict: + return mesh_to_polyhedron(build_convex_hull(points)) + + +def calculate_center_point(points: np.ndarray) -> np.ndarray: + """ + Calculate the center point of a point cloud. + + :param points: point coordinates with shape (n, 3) + :return: arithmetic mean of all points + :raises ValueError: if the param points have an invalid shape or are empty + """ + if points.ndim != 2 or points.shape[1] != 3: + raise ValueError(f"Expected points with shape (n, 3), got {points.shape}.") + if len(points) == 0: + raise ValueError("At least one point is required to calculate a centroid.") + + return points.mean(axis=0) + + +def find_farthest_point(points: np.ndarray, reference_point: np.ndarray) -> np.ndarray: + """ + Calculate the farthest point inside a point cloud from a reference point. + + :param points: point coordinates with shape (n, 3) + :param reference_point: point from which distances are calculated + :return: point with the greatest distance from the reference point + :raises ValueError: if the reference point or the param points have an invalid shape or are empty + """ + if points.ndim != 2 or points.shape[1] != 3: + raise ValueError(f"Expected points with shape (n, 3), got {points.shape}.") + if len(points) == 0: + raise ValueError( + "At least one point is required to calculate a maximum distance." + ) + if reference_point.shape != (3,): + raise ValueError( + f"Expected reference_point with shape (3,), got {reference_point.shape}." + ) + + distances = np.linalg.norm(points - reference_point, axis=1) + return points[np.argmax(distances)] + + +def find_farthest_point_vdw( + points: np.ndarray, elements: np.ndarray, reference_point: np.ndarray +) -> tuple[np.ndarray, float]: + """ + Calculate the point, whose fdw-"Bubble" is most distant to a reference point. + + :param points: atom coordinates with shape (n, 3) + :param elements: element symbol mapping to each atom coordinate + :param reference_point: point from which distances are calculated + :return: farthest atom coordinate and its distance including the van der Waals radius + :raises ValueError: if the reference point or the param points have an invalid shape or are empty + """ + if points.ndim != 2 or points.shape[1] != 3: + raise ValueError(f"Expected points with shape (n, 3), got {points.shape}.") + if len(points) == 0: + raise ValueError( + "At least one point is required to calculate a maximum distance." + ) + if reference_point.shape != (3,): + raise ValueError( + f"Expected reference_point with shape (3,), got {reference_point.shape}." + ) + + radii = np.array([vdw_radii[element] for element in elements], dtype=float) + + distances = np.linalg.norm(points - reference_point, axis=1) + radii + max_index = np.argmax(distances) + return points[max_index], float(distances[max_index]) + + +def find_intersecting_spheres( + spheres: list[dict], + reference_sphere: dict, +) -> list[dict]: + """ + Find spheres that intersect or touch a reference sphere. + + :param spheres: spheres represented by center coordinates and radius + :param reference_sphere: sphere against which intersections are checked + :return: spheres intersecting or touching the reference sphere + """ + reference_center = np.array(reference_sphere["center"], dtype=float) + reference_radius = float(reference_sphere["radius"]) + + intersecting_spheres = [] + for sphere in spheres: + center = np.array(sphere["center"], dtype=float) + radius = float(sphere["radius"]) + distance = np.linalg.norm(center - reference_center) + + if distance <= reference_radius + radius: + intersecting_spheres.append(sphere) + + return intersecting_spheres + + +def meshes_intersect( + mesh_a: trimesh.Trimesh, + mesh_b: trimesh.Trimesh, + distance_tolerance: float = 1e-9, +) -> bool: + """ + Determine whether two meshes intersect or touch. + + :param mesh_a: first mesh + :param mesh_b: second mesh + :param distance_tolerance: maximum distance at which meshes count as touching + (standard is 1e-9 because of floating point error) + :return: whether the meshes intersect or are within the distance tolerance + """ + + from trimesh.collision import CollisionManager + + manager = CollisionManager() + manager.add_object("mesh_a", mesh_a) + return manager.min_distance_single(mesh_b) <= distance_tolerance + + +def meshes_distance(mesh_a: trimesh.Trimesh, mesh_b: trimesh.Trimesh) -> float: + """ + Calculate the minimum euclidean distance between two meshes. + + :param mesh_a: first mesh + :param mesh_b: second mesh + :return: minimum Euclidean distance between the meshes + """ + + from trimesh.collision import CollisionManager + + manager = CollisionManager() + manager.add_object("mesh_a", mesh_a) + distance = float(manager.min_distance_single(mesh_b)) + return max(distance, 0.0) diff --git a/backend/protzilla/data_analysis/ptm_analysis.py b/backend/protzilla/data_analysis/ptm_analysis.py index 7174e8a09..983332899 100644 --- a/backend/protzilla/data_analysis/ptm_analysis.py +++ b/backend/protzilla/data_analysis/ptm_analysis.py @@ -1,8 +1,13 @@ import numpy as np import pandas as pd import re +import logging from backend.protzilla.utilities.transform_dfs import long_to_wide +from backend.protzilla.data_analysis.amino_acid_spheres import ( + find_ptm_amino_acid_sphere_collisions, +) +from backend.protzilla.steps import OutputItem, OutputType def ptms_per_sample(psm_df: pd.DataFrame) -> dict: @@ -91,3 +96,53 @@ def from_string(mod_string: str) -> tuple[int, str]: name = name[1:] if name[0] == " " else name return amount, name + + +def ptm_validation( + cif_df: pd.DataFrame, + structure_metadata_df: pd.DataFrame, + ignored_neighbors: int = 0, +): + ignored_neighbors = int(ignored_neighbors) + collisions = find_ptm_amino_acid_sphere_collisions( + cif_df, ignored_neighbors=ignored_neighbors + ) + messages = [] + if not collisions: + messages.append( + dict( + level=logging.WARNING, + msg="No PTMs were found in the provided structure or PTM annotations are missing.", + ) + ) + collision_columns = [ + "ptm_name", + "ptm_chain", + "ptm_position", + "amino_acid_chain", + "amino_acid_position", + ] + collision_rows = [ + { + "ptm_name": collision["ptm_name"], + "ptm_chain": collision["chain"], + "ptm_position": collision["position"], + "amino_acid_chain": amino_acid["chain"], + "amino_acid_position": amino_acid["position"], + } + for collision in collisions + for amino_acid in collision["collisions"] + ] + data_for_visualization = { + "structure_entry_id": structure_metadata_df["entry_id"].iloc[0], + "cif_df": cif_df, + "include_ptm_spheres": True, + "ignored_neighbors": ignored_neighbors, + } + return { + "ptm_collisions_df": pd.DataFrame(collision_rows, columns=collision_columns), + "visualization": OutputItem( + output_type=OutputType.VISUALIZATION, value=data_for_visualization + ), + "messages": messages, + } diff --git a/backend/protzilla/methods/data_analysis.py b/backend/protzilla/methods/data_analysis.py index bab8043f1..bde74e40f 100644 --- a/backend/protzilla/methods/data_analysis.py +++ b/backend/protzilla/methods/data_analysis.py @@ -50,6 +50,7 @@ HEATMAP_HIGH_COLOR, ) from backend.protzilla.data_analysis.ptm_analysis import ( + ptm_validation, ptms_per_protein_and_sample, ptms_per_sample, ) @@ -2502,3 +2503,25 @@ def create_form(self): ), ], ) + + +class PtmValidation(PeptideAnalysisStep): + display_name = "PTM Validation" + method_description = "Validates PTMs in protein structure predictions." + output_keys = ["ptm_collisions_df"] + calc_method = staticmethod(ptm_validation) + + def create_form(self): + return Form( + label="PTM Validation", + input_fields=[ + NumberField( + name="ignored_neighbors", + label="Number of PTM residue neighbors to ignore", + value=0, + min=0, + step=1, + hasStepButtons=True, + ) + ], + ) diff --git a/backend/protzilla/utilities/ptm_helper.py b/backend/protzilla/utilities/ptm_helper.py new file mode 100644 index 000000000..be36a8738 --- /dev/null +++ b/backend/protzilla/utilities/ptm_helper.py @@ -0,0 +1,107 @@ +import pandas as pd +import numpy as np + +from backend.protzilla.constants.cif_columns import ( + ATOM_SITE_COLUMNS, + CHEM_COMP_COLUMNS, +) +from backend.protzilla.data_analysis.geometry_operations import find_farthest_point_vdw + + +def get_all_ptm_atoms_with_coordinates(cif_df: pd.DataFrame) -> list: + """ + Extracts all non-standard residues (PTMs) and their atomic coordinates + from a parsed mmCIF DataFrame. + """ + ptm_list = [] + + is_non_standard = cif_df[CHEM_COMP_COLUMNS.MON_NSTD_FLAG] == False + ptm_atoms_df = cif_df[is_non_standard] + + grouped_ptms = ptm_atoms_df.groupby( + [ATOM_SITE_COLUMNS.LABEL_ASYM_ID, ATOM_SITE_COLUMNS.LABEL_SEQ_ID] + ) + + for (chain_id, seq_id), atom_group in grouped_ptms: + comp_id = atom_group[ATOM_SITE_COLUMNS.LABEL_COMP_ID].iloc[0] + atoms_subset = atom_group[ + [ + ATOM_SITE_COLUMNS.LABEL_ATOM_ID, + ATOM_SITE_COLUMNS.TYPE_SYMBOL, + ATOM_SITE_COLUMNS.CARTN_X, + ATOM_SITE_COLUMNS.CARTN_Y, + ATOM_SITE_COLUMNS.CARTN_Z, + ] + ].rename( + columns={ + ATOM_SITE_COLUMNS.LABEL_ATOM_ID: "atom_name", + ATOM_SITE_COLUMNS.TYPE_SYMBOL: "element", + ATOM_SITE_COLUMNS.CARTN_X: "x", + ATOM_SITE_COLUMNS.CARTN_Y: "y", + ATOM_SITE_COLUMNS.CARTN_Z: "z", + } + ) + atoms_subset["element"] = ( + atoms_subset["element"].astype(str).str.strip().str.capitalize() + ) + + atoms_data = atoms_subset.to_dict(orient="records") + + ptm_list.append( + { + "ptm_name": comp_id, + "chain": chain_id, + "position": seq_id, + "atoms": atoms_data, + } + ) + + return ptm_list + + +def calculate_center_point(points: np.ndarray) -> np.ndarray: + """ + Calculate the center point of a point cloud. + """ + if points.ndim != 2 or points.shape[1] != 3: + raise ValueError(f"Expected points with shape (n, 3), got {points.shape}.") + if len(points) == 0: + raise ValueError("At least one point is required to calculate a centroid.") + + return points.mean(axis=0) + + +def find_farthest_point(points: np.ndarray, reference_point: np.ndarray) -> np.ndarray: + """ + Calculate the farthest point inside a point cloud from a reference point. + """ + if points.ndim != 2 or points.shape[1] != 3: + raise ValueError(f"Expected points with shape (n, 3), got {points.shape}.") + if len(points) == 0: + raise ValueError( + "At least one point is required to calculate a maximum distance." + ) + if reference_point.shape != (3,): + raise ValueError( + f"Expected reference_point with shape (3,), got {reference_point.shape}." + ) + + distances = np.linalg.norm(points - reference_point, axis=1) + return points[np.argmax(distances)] + + +def get_center_points_and_radius_for_each_ptm(ptm_list: list) -> list: + for ptm in ptm_list: + coords = [[atom["x"], atom["y"], atom["z"]] for atom in ptm["atoms"]] + elements = [atom["element"] for atom in ptm["atoms"]] + + coords_array = np.array(coords, dtype=np.float32) + elements_array = np.array(elements) + + center_point = calculate_center_point(coords_array) + _, radius = find_farthest_point_vdw(coords_array, elements_array, center_point) + + ptm["center_point"] = center_point.tolist() + ptm["radius"] = float(radius) + + return ptm_list diff --git a/backend/tests/main/test_views_helper.py b/backend/tests/main/test_views_helper.py index b8e38d9cd..9c03c8cf5 100644 --- a/backend/tests/main/test_views_helper.py +++ b/backend/tests/main/test_views_helper.py @@ -92,6 +92,7 @@ def test_get_all_possible_step_names(): "ArbitraryCSVImport", "CrosslinkingValidationWithAngstromDeviation", "CrosslinkingValidationWithAngstromDeviationForMultimer", + "PtmValidation", } steps = get_all_possible_steps() diff --git a/backend/tests/protzilla/utilities/test_ptm_helper.py b/backend/tests/protzilla/utilities/test_ptm_helper.py new file mode 100644 index 000000000..b25685614 --- /dev/null +++ b/backend/tests/protzilla/utilities/test_ptm_helper.py @@ -0,0 +1,108 @@ +import pandas as pd +from backend.protzilla.utilities.ptm_helper import ( + get_all_ptm_atoms_with_coordinates, + get_center_points_and_radius_for_each_ptm, +) +import pytest +import numpy as np + + +def test_get_all_ptm_atoms_with_coordinates(): + data = { + "_chem_comp.mon_nstd_flag": [True, True, False, False], + "_atom_site.label_asym_id": ["A", "A", "B", "B"], + "_atom_site.label_seq_id": [10, 10, 42, 42], + "_atom_site.label_comp_id": ["MET", "MET", "SEP", "SEP"], + "_atom_site.label_atom_id": ["N", "CA", "P", "O1P"], + "_atom_site.type_symbol": ["N", "C", "P", "O"], + "_atom_site.Cartn_x": [1.0, 2.0, 10.0, 11.0], + "_atom_site.Cartn_y": [1.1, 2.1, 10.1, 11.1], + "_atom_site.Cartn_z": [1.2, 2.2, 10.2, 11.2], + } + cif_df = pd.DataFrame(data) + + result = get_all_ptm_atoms_with_coordinates(cif_df) + + assert len(result) == 1 + + ptm = result[0] + assert ptm["ptm_name"] == "SEP" + assert ptm["chain"] == "B" + assert ptm["position"] == 42 + + atoms = ptm["atoms"] + assert len(atoms) == 2 + + assert atoms[0] == { + "atom_name": "P", + "element": "P", + "x": 10.0, + "y": 10.1, + "z": 10.2, + } + + assert atoms[1] == { + "atom_name": "O1P", + "element": "O", + "x": 11.0, + "y": 11.1, + "z": 11.2, + } + + +def test_get_center_points_and_radius_for_each_ptm(): + ptm_list = [ + { + "ptm_name": "SYM", + "chain": "A", + "position": 1, + "atoms": [ + {"atom_name": "A1", "x": 1.0, "y": 0.0, "z": 0.0}, + {"atom_name": "A2", "x": -1.0, "y": 0.0, "z": 0.0}, + {"atom_name": "A3", "x": 0.0, "y": 1.0, "z": 0.0}, + {"atom_name": "A4", "x": 0.0, "y": -1.0, "z": 0.0}, + ], + }, + { + "ptm_name": "SGL", + "chain": "B", + "position": 2, + "atoms": [ + {"atom_name": "A1", "x": 5.0, "y": 5.0, "z": 5.0}, + ], + }, + { + "ptm_name": "AB", + "chain": "C", + "position": 3, + "atoms": [ + {"atom_name": "A1", "x": 12.0, "y": 10.0, "z": 8.0}, + {"atom_name": "A2", "x": 8.0, "y": 10.0, "z": 7.0}, + {"atom_name": "A3", "x": 10.0, "y": 10.0, "z": 15.0}, + ], + }, + ] + for ptm in ptm_list: + for atom in ptm["atoms"]: + atom["element"] = "C" + + result_list = get_center_points_and_radius_for_each_ptm(ptm_list) + + ptm_1 = result_list[0] + + result_list = get_center_points_and_radius_for_each_ptm(ptm_list) + + ptm_1 = result_list[0] + np.testing.assert_allclose(ptm_1["center_point"], [0.0, 0.0, 0.0], atol=1e-6) + assert pytest.approx(ptm_1["radius"], 1e-6) == 2.77 + + ptm_2 = result_list[1] + np.testing.assert_allclose(ptm_2["center_point"], [5.0, 5.0, 5.0], atol=1e-6) + assert pytest.approx(ptm_2["radius"], 1e-6) == 1.77 + + ptm_3 = result_list[2] + + expected_center = [10.0, 10.0, 10.0] + np.testing.assert_allclose(ptm_3["center_point"], expected_center, atol=1e-6) + + assert pytest.approx(ptm_3["radius"], 1e-6) == 6.77 diff --git a/frontend/src/components/app/run-screen/run-screen.tsx b/frontend/src/components/app/run-screen/run-screen.tsx index 3cf3224b6..5e9774348 100644 --- a/frontend/src/components/app/run-screen/run-screen.tsx +++ b/frontend/src/components/app/run-screen/run-screen.tsx @@ -29,7 +29,6 @@ import { Col } from "react-grid-system"; import { useLocation, useNavigate } from "react-router-dom"; import { styled } from "styled-components"; -import { CrosslinkerInformation } from "../../core/shared/molstar-viewer/crosslinker-processing"; import { H3 } from "../../core/shared/text"; const StyledNavbar = styled(Navbar)` @@ -155,13 +154,14 @@ export const RunScreen: React.FC = () => { structureEntryId: response.data.structureEntryId, cifString: response.data.cifString, crosslinks: response.data.crosslinks, + trimeshMeshes: response.data.trimeshMeshes, }), [], ); const visualizations = useCertainStepOutputs< StepOutputInfo, ApiResponse, - { structureEntryId: string; cifString: string; crosslinks?: CrosslinkerInformation[] } + Visualization >({ available_outputs: availableVisualizations, endpoint: "get_step_visualizations/", @@ -380,7 +380,11 @@ export const RunScreen: React.FC = () => { {visualizations.length > 0 ? ( visualizations.map((viz) => ( - + )) ) : ( diff --git a/frontend/src/components/core/shared/molstar-viewer/index.ts b/frontend/src/components/core/shared/molstar-viewer/index.ts index 45fc7b6db..005c3ac68 100644 --- a/frontend/src/components/core/shared/molstar-viewer/index.ts +++ b/frontend/src/components/core/shared/molstar-viewer/index.ts @@ -1,2 +1,3 @@ export { default as MolstarViewer } from "./molstar-viewer"; export * from "./molstar-viewer.props"; +export * from "./molstar-trimesh-adapter"; diff --git a/frontend/src/components/core/shared/molstar-viewer/molstar-trimesh-adapter.ts b/frontend/src/components/core/shared/molstar-viewer/molstar-trimesh-adapter.ts new file mode 100644 index 000000000..859e060ca --- /dev/null +++ b/frontend/src/components/core/shared/molstar-viewer/molstar-trimesh-adapter.ts @@ -0,0 +1,103 @@ +import { Mesh } from "molstar/lib/mol-geo/geometry/mesh/mesh"; +import { MeshBuilder } from "molstar/lib/mol-geo/geometry/mesh/mesh-builder"; +import { Vec3 } from "molstar/lib/mol-math/linear-algebra"; +import { Shape } from "molstar/lib/mol-model/shape"; +import { PluginStateObject } from "molstar/lib/mol-plugin-state/objects"; +import { StateTransforms } from "molstar/lib/mol-plugin-state/transforms"; +import { PluginUIContext } from "molstar/lib/mol-plugin-ui/context"; +import { StateTransformer } from "molstar/lib/mol-state"; +import { Task } from "molstar/lib/mol-task"; +import { Color } from "molstar/lib/mol-util/color"; +import { ParamDefinition as PD } from "molstar/lib/mol-util/param-definition"; + +export interface TrimeshMesh { + vertices: number[][]; + faces: number[][]; +} + +export interface AddTrimeshMeshOptions { + color?: number; + alpha?: number; + label?: string; +} + +const ProtzillaTransforms = StateTransformer.builderFactory("protzilla"); + +const TrimeshMeshShape = ProtzillaTransforms({ + name: "trimesh-shape", + display: { name: "Trimesh Shape" }, + from: PluginStateObject.Root, + to: PluginStateObject.Shape.Provider, + params: { + mesh: PD.Value({ vertices: [], faces: [] }, { isHidden: true }), + color: PD.Color(Color(0xff8800)), + label: PD.Text("Trimesh Mesh"), + }, +})({ + canAutoUpdate: () => true, + apply({ params }) { + return Task.create("Create Trimesh Shape", () => { + return Promise.resolve( + new PluginStateObject.Shape.Provider( + { + label: params.label, + data: params.mesh, + params: Mesh.Params, + geometryUtils: Mesh.Utils, + getShape: (_ctx, data) => { + const builder = MeshBuilder.createState(data.vertices.length, data.faces.length); + builder.currentGroup = 0; + + for (const face of data.faces) { + if (face.length !== 3) continue; + const a = data.vertices[face[0]]; + const b = data.vertices[face[1]]; + const c = data.vertices[face[2]]; + if (a === undefined || b === undefined || c === undefined) continue; + + MeshBuilder.addTriangle( + builder, + Vec3.create(a[0], a[1], a[2]), + Vec3.create(b[0], b[1], b[2]), + Vec3.create(c[0], c[1], c[2]), + ); + } + + const mesh = MeshBuilder.getMesh(builder); + return Shape.create( + params.label, + data, + mesh, + () => params.color, + () => 1, + () => params.label, + ); + }, + }, + { label: params.label }, + ), + ); + }); + }, +}); + +export async function addTrimeshMesh( + plugin: PluginUIContext, + mesh: TrimeshMesh, + options: AddTrimeshMeshOptions = {}, +) { + if (mesh.vertices.length === 0 || mesh.faces.length === 0) return; + + const label = options.label ?? "Trimesh Mesh"; + const color = Color(options.color ?? 0xff8800); + const alpha = options.alpha ?? 0.35; + + await plugin + .build() + .toRoot() + .apply(TrimeshMeshShape, { mesh, color, label }) + .apply(StateTransforms.Representation.ShapeRepresentation3D, { alpha }) + .commit(); +} + +export const addTrimeshPolyhedron = addTrimeshMesh; diff --git a/frontend/src/components/core/shared/molstar-viewer/molstar-viewer.props.tsx b/frontend/src/components/core/shared/molstar-viewer/molstar-viewer.props.tsx index abf936f4d..c19f0391c 100644 --- a/frontend/src/components/core/shared/molstar-viewer/molstar-viewer.props.tsx +++ b/frontend/src/components/core/shared/molstar-viewer/molstar-viewer.props.tsx @@ -1,6 +1,15 @@ import { CrosslinkerInformation } from "./crosslinker-processing"; +import { TrimeshMesh } from "./molstar-trimesh-adapter"; + +export interface TrimeshShape { + label: string; + mesh: TrimeshMesh; + color?: number; + alpha?: number; +} export interface MolstarViewerProps { cifText: string; crosslinks?: CrosslinkerInformation[]; + trimeshMeshes?: TrimeshShape[]; } diff --git a/frontend/src/components/core/shared/molstar-viewer/molstar-viewer.tsx b/frontend/src/components/core/shared/molstar-viewer/molstar-viewer.tsx index 75f638bd5..acc08b3e6 100644 --- a/frontend/src/components/core/shared/molstar-viewer/molstar-viewer.tsx +++ b/frontend/src/components/core/shared/molstar-viewer/molstar-viewer.tsx @@ -5,6 +5,7 @@ import { PluginUIContext } from "molstar/lib/mol-plugin-ui/context"; import { renderReact18 } from "molstar/lib/mol-plugin-ui/react18"; import React, { useEffect, useRef, useState } from "react"; +import { addTrimeshMesh } from "./molstar-trimesh-adapter"; import { MolstarViewerProps } from "./molstar-viewer.props"; import { addCrosslinks, @@ -16,7 +17,7 @@ import { LegendOverlay } from "./molstar-viewer.ui"; import { CanvasWrapper, Container } from "./styles"; import "molstar/lib/mol-plugin-ui/skin/base/base.scss"; -const MolstarViewer: React.FC = ({ cifText, crosslinks }) => { +const MolstarViewer: React.FC = ({ cifText, crosslinks, trimeshMeshes }) => { const containerRef = useRef(null); const [isLoading, setIsLoading] = useState(true); const notify = useNotification(); @@ -56,6 +57,16 @@ const MolstarViewer: React.FC = ({ cifText, crosslinks }) => overrideLabels(plugin, crosslinkColors); } + if (trimeshMeshes !== undefined) { + for (const trimeshMesh of trimeshMeshes) { + await addTrimeshMesh(plugin, trimeshMesh.mesh, { + color: trimeshMesh.color, + alpha: trimeshMesh.alpha, + label: trimeshMesh.label, + }); + } + } + setIsLoading(false); } catch (error: unknown) { handleError(error, "MolstarViewer Error:", notify); @@ -74,7 +85,7 @@ const MolstarViewer: React.FC = ({ cifText, crosslinks }) => } } }; - }, [cifText, crosslinks, notify]); + }, [cifText, crosslinks, notify, trimeshMeshes]); return ( diff --git a/frontend/src/utils/protzilla-types.ts b/frontend/src/utils/protzilla-types.ts index 0b347d27a..c88f780eb 100644 --- a/frontend/src/utils/protzilla-types.ts +++ b/frontend/src/utils/protzilla-types.ts @@ -35,10 +35,23 @@ export interface Download { json_downloads: Record; } +export interface TrimeshMesh { + vertices: number[][]; + faces: number[][]; +} + +export interface TrimeshShape { + label: string; + mesh: TrimeshMesh; + color?: number; + alpha?: number; +} + export interface Visualization { structureEntryId: string; cifString: string; crosslinks?: CrosslinkerInformation[]; + trimeshMeshes?: TrimeshShape[]; } // We assume these are the only data types we receive for tables diff --git a/requirements.txt b/requirements.txt index 31426de7c..574dedc9c 100644 --- a/requirements.txt +++ b/requirements.txt @@ -32,6 +32,7 @@ pytest-django==4.11.1 pytest-order==1.3.0 PyYAML==6.0.3 py7zr!=1.1.2 +python-fcl==0.7.0.11 requests==2.32.5 restring==0.1.21 scikit-learn==1.7.2 @@ -40,6 +41,7 @@ seaborn==0.13.2 sphinx==8.2.3 sphinx-autoapi==3.6.1 statsmodels==0.14.5 +trimesh==4.11.5 tqdm==4.67.1 tzlocal==5.3.1 umap-learn==0.5.9.post2