Bio Structural Biology Geometric Analysis

from Bio.PDB import PDBParser
import numpy as np

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')

chain = structure[0]['A']
atom1, atom2 = chain[100]['CA'], chain[200]['CA']

distance = atom1 - atom2                              # Atom subtraction returns the distance directly
print(f'Distance: {distance:.2f} A')
print(np.linalg.norm(atom1.coord - atom2.coord))     # Equivalent via numpy on the .coord arrays

import numpy as np
from Bio.PDB import PDBParser

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')

ca_atoms = [r['CA'] for r in structure.get_residues() if r.has_id('CA') and r.id[0] == ' ']  # id[0]==' ' drops waters/hetero
n = len(ca_atoms)

dist = np.zeros((n, n))
for i in range(n):
    for j in range(i + 1, n):
        dist[i, j] = dist[j, i] = ca_atoms[i] - ca_atoms[j]
print(f'Distance matrix: {dist.shape}')

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

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')

res = structure[0]['A'][100]
angle = calc_angle(res['N'].get_vector(), res['CA'].get_vector(), res['C'].get_vector())  # calc_angle needs Vector, not .coord
print(f'N-CA-C angle: {np.degrees(angle):.1f} deg')

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

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')
chain = structure[0]['A']

prev, curr, nxt = chain[99], chain[100], chain[101]
phi = calc_dihedral(prev['C'].get_vector(), curr['N'].get_vector(), curr['CA'].get_vector(), curr['C'].get_vector())
psi = calc_dihedral(curr['N'].get_vector(), curr['CA'].get_vector(), curr['C'].get_vector(), nxt['N'].get_vector())
print(f'phi={np.degrees(phi):.1f}  psi={np.degrees(psi):.1f}')

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

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')

ppb = PPBuilder()                                    # PPBuilder builds peptides from connectivity, so chain breaks split them
rama = []
for pp in ppb.build_peptides(structure):
    for res, (phi, psi) in zip(pp, pp.get_phi_psi_list()):
        if phi is not None and psi is not None:      # None at termini and chain breaks by design - skip, do not fabricate
            rama.append((res.resname, np.degrees(phi), np.degrees(psi)))
print(f'{len(rama)} residues with phi/psi')

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

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')
res = structure[0]['A'][100]

if res.has_id('CB') and res.has_id('CG'):            # Gly/Ala lack CB/CG; chi atom quartets are residue-type specific
    chi1 = calc_dihedral(res['N'].get_vector(), res['CA'].get_vector(), res['CB'].get_vector(), res['CG'].get_vector())
    print(f'Chi1: {np.degrees(chi1):.1f} deg')

from Bio.PDB import PDBParser, Superimposer

parser = PDBParser(QUIET=True)
ref = parser.get_structure('ref', 'reference.pdb')
mob = parser.get_structure('mobile', 'mobile.pdb')

ref_ca = [r['CA'] for r in ref.get_residues() if r.has_id('CA') and r.id[0] == ' ']
mob_ca = [r['CA'] for r in mob.get_residues() if r.has_id('CA') and r.id[0] == ' ']
n = min(len(ref_ca), len(mob_ca))                    # Superimposer needs EQUAL-LENGTH ORDERED lists; it does NOT solve correspondence
ref_ca, mob_ca = ref_ca[:n], mob_ca[:n]              # Naive truncation is only valid when residues already correspond 1:1

sup = Superimposer()
sup.set_atoms(ref_ca, mob_ca)                        # Optimal rigid-body fit via SVD/Kabsch
print(f'RMSD (CA): {sup.rms:.2f} A')
rotation, translation = sup.rotran                   # The fitted transform, for reuse on other atoms
sup.apply(mob.get_atoms())                           # Mutates mob in place - copy first if the originals are still needed

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

parser = PDBParser(QUIET=True)
ref = parser.get_structure('ref', 'reference.pdb')
mob = parser.get_structure('mobile', 'mobile.pdb')

ref_ca = [r['CA'] for r in ref.get_residues() if r.has_id('CA') and r.id[0] == ' ']
mob_ca = [r['CA'] for r in mob.get_residues() if r.has_id('CA') and r.id[0] == ' ']
n = min(len(ref_ca), len(mob_ca))
ref_ca, mob_ca = ref_ca[:n], mob_ca[:n]

sup = Superimposer()
sup.set_atoms(ref_ca, mob_ca)
sup.apply([a for a in mob_ca])                        # Move only the paired CA set into the fitted frame
deviation = np.array([r - m for r, m in zip(ref_ca, mob_ca)])
print(f'core (<2A) residues: {(deviation < 2.0).sum()} / {n}')   # 2A is a common rigid-core cutoff, not a law
print(f'max deviation: {deviation.max():.2f} A at index {deviation.argmax()}')

import numpy as np
from Bio.PDB import PDBParser

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')

atoms = list(structure.get_atoms())
coords = np.array([a.coord for a in atoms])
masses = np.array([{'C': 12.0, 'N': 14.0, 'O': 16.0, 'S': 32.0, 'H': 1.0}.get(a.element, 12.0) for a in atoms])

com = (masses[:, None] * coords).sum(axis=0) / masses.sum()
print(f'Center of mass: {com}')

rg = np.sqrt(np.mean(np.sum((coords - coords.mean(axis=0)) ** 2, axis=1)))  # Unweighted radius of gyration
print(f'Radius of gyration: {rg:.2f} A')

from Bio.PDB import PDBParser

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')

v1 = structure[0]['A'][100]['CA'].get_vector()
v2 = structure[0]['A'][101]['CA'].get_vector()

diff = v2 - v1
print(f'length: {diff.norm():.2f}  unit: {diff.normalized()}')
cross = v1 ** v2                                      # ** is cross product on Vector objects
dot = v1 * v2                                         # * is dot product on Vector objects

from Bio.PDB import PDBParser
from Bio.PDB.SASA import ShrakeRupley

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')

sr = ShrakeRupley(probe_radius=1.40, n_points=100)   # 1.4A = water radius (convention); a SASA number is meaningless without its probe radius
sr.compute(structure, level='R')                     # level R attaches .sasa on each residue; children sum to parents
print(f'total SASA: {sum(r.sasa for r in structure.get_residues() if hasattr(r, "sasa")):.1f} A^2')

from Bio.PDB import PDBParser
from Bio.PDB.SASA import ShrakeRupley

MAX_ASA = {'ALA': 129.0, 'ARG': 274.0, 'ASN': 195.0, 'ASP': 193.0, 'CYS': 167.0, 'GLU': 223.0, 'GLN': 225.0, 'GLY': 104.0, 'HIS': 224.0, 'ILE': 197.0, 'LEU': 201.0, 'LYS': 236.0, 'MET': 224.0, 'PHE': 240.0, 'PRO': 159.0, 'SER': 155.0, 'THR': 172.0, 'TRP': 285.0, 'TYR': 263.0, 'VAL': 174.0}

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')
ShrakeRupley().compute(structure, level='R')

buried = 0
for res in structure.get_residues():
    if res.resname in MAX_ASA and hasattr(res, 'sasa'):
        rsa = res.sasa / MAX_ASA[res.resname]
        if rsa < 0.20:                               # RSA < 0.20 is the common buried heuristic (a rule of thumb, not a law)
            buried += 1
print(f'buried residues (RSA < 0.20): {buried}')

from Bio.PDB import PDBParser, DSSP

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')
model = structure[0]                                  # DSSP runs on ONE model only

dssp = DSSP(model, 'protein.pdb', dssp='mkdssp')      # pass the current binary name explicitly
codes = [dssp[k][2] for k in dssp.keys()]             # 8-state: H G I helix, E B strand, T S P - other
helix = sum(c in 'HGI' for c in codes)
strand = sum(c in 'EB' for c in codes)
print(f'helix {helix}, strand {strand}, other {len(codes) - helix - strand} of {len(codes)}')