mirror of
https://github.com/biopython/biopython.git
synced 2025-10-20 21:53:47 +08:00
de0bb21fb3
$ ruff check --fix --select=I \ --config=lint.isort.force-single-line=true \ --config=lint.isort.order-by-type=false \ BioSQL/ Bio/ Tests/ Scripts/ Doc/ setup.py Using ruff version 0.4.10
665 lines
27 KiB
Python
665 lines
27 KiB
Python
# Copyright 2020-2022 by Robert T. Miller. All rights reserved.
|
|
# This file is part of the Biopython distribution and governed by your
|
|
# choice of the "Biopython License Agreement" or the "BSD 3-Clause License".
|
|
# Please see the LICENSE file that should have been included as part of this
|
|
# package.
|
|
|
|
"""Unit tests for internal_coords module of Bio.PDB."""
|
|
|
|
import copy
|
|
import re
|
|
import unittest
|
|
import warnings
|
|
|
|
try:
|
|
import numpy as np # noqa F401
|
|
except ImportError:
|
|
from Bio import MissingPythonDependencyError
|
|
|
|
raise MissingPythonDependencyError(
|
|
"Install NumPy if you want to use Bio.PDB."
|
|
) from None
|
|
|
|
from io import StringIO
|
|
|
|
from Bio.File import as_handle
|
|
from Bio.PDB.ic_rebuild import compare_residues
|
|
from Bio.PDB.ic_rebuild import IC_duplicate
|
|
from Bio.PDB.ic_rebuild import structure_rebuild_test
|
|
from Bio.PDB.internal_coords import AtomKey
|
|
from Bio.PDB.internal_coords import Dihedron
|
|
from Bio.PDB.internal_coords import IC_Chain
|
|
from Bio.PDB.internal_coords import IC_Residue
|
|
from Bio.PDB.MMCIFParser import MMCIFParser
|
|
from Bio.PDB.mmtf import MMTFParser
|
|
from Bio.PDB.PDBExceptions import PDBConstructionWarning
|
|
from Bio.PDB.PDBParser import PDBParser
|
|
from Bio.PDB.PICIO import read_PIC_seq
|
|
from Bio.PDB.SCADIO import write_SCAD
|
|
from Bio.Seq import Seq
|
|
from Bio.SeqRecord import SeqRecord
|
|
|
|
|
|
class Rebuild(unittest.TestCase):
|
|
"""Read PDB and mmCIF structures, convert to/from internal coordinates."""
|
|
|
|
PDB_parser = PDBParser(PERMISSIVE=True, QUIET=True)
|
|
CIF_parser = MMCIFParser(QUIET=True)
|
|
MMTF_parser = MMTFParser()
|
|
pdb_1LCD = PDB_parser.get_structure("1LCD", "PDB/1LCD.pdb")
|
|
# cif_1A7G = CIF_parser.get_structure("1A7G", "PDB/1A7G.cif")
|
|
# cif_1A7G2 = CIF_parser.get_structure("1A7G", "PDB/1A7G.cif")
|
|
pdb_2XHE = PDB_parser.get_structure("2XHE", "PDB/2XHE.pdb")
|
|
pdb_2XHE2 = PDB_parser.get_structure("2XHE", "PDB/2XHE.pdb")
|
|
pdb_1A8O = PDB_parser.get_structure("1A8O", "PDB/1A8O.pdb")
|
|
cif_3JQH = CIF_parser.get_structure("3JQH", "PDB/3JQH.cif")
|
|
cif_3JQH2 = CIF_parser.get_structure("3JQH", "PDB/3JQH.cif")
|
|
cif_4CUP = CIF_parser.get_structure("4CUP", "PDB/4CUP.cif")
|
|
cif_4CUP2 = CIF_parser.get_structure("4CUP", "PDB/4CUP.cif")
|
|
cif_4ZHL = CIF_parser.get_structure("4ZHL", "PDB/4ZHL.cif")
|
|
cif_4ZHL2 = CIF_parser.get_structure("4ZHL", "PDB/4ZHL.cif")
|
|
with warnings.catch_warnings(record=True) as w:
|
|
warnings.simplefilter("always", PDBConstructionWarning)
|
|
mmtf_1A8O = MMTF_parser.get_structure("PDB/1A8O.mmtf")
|
|
|
|
def test_mmtf(self):
|
|
chain = next(self.mmtf_1A8O.get_chains())
|
|
ic_chain = IC_Chain(chain)
|
|
self.assertEqual(len(ic_chain.ordered_aa_ic_list), 70)
|
|
|
|
def test_rebuild_1a8o(self):
|
|
"""Duplicate tutorial doctests which fail under linux."""
|
|
IC_Chain.MaxPeptideBond = 4.0
|
|
r = structure_rebuild_test(self.pdb_1A8O, False)
|
|
self.assertTrue(r["pass"])
|
|
|
|
cic = list(self.pdb_1A8O.get_chains())[0].internal_coord
|
|
distances = cic.distance_plot()
|
|
chirality = cic.dihedral_signs()
|
|
c2 = IC_duplicate(cic.chain)[0]["A"]
|
|
cic2 = c2.internal_coord
|
|
cic2.atomArray = np.zeros((cic2.AAsiz, 4), dtype=np.float64)
|
|
cic2.dihedraAngle[:] = 0.0
|
|
cic2.hedraAngle[:] = 0.0
|
|
cic2.hedraL12[:] = 0.0
|
|
cic2.hedraL23[:] = 0.0
|
|
cic2.copy_initNCaCs(cic)
|
|
cic2.distplot_to_dh_arrays(distances, chirality)
|
|
cic2.distance_to_internal_coordinates()
|
|
c2.internal_to_atom_coordinates()
|
|
np.allclose(cic2.atomArray, cic.atomArray)
|
|
|
|
def test_rebuild_multichain_missing(self):
|
|
"""Convert multichain missing atom struct to, from internal coords."""
|
|
# 2XHE has regions of missing chain, last residue has only N
|
|
r = structure_rebuild_test(self.pdb_2XHE, False)
|
|
self.assertEqual(r["residues"], 787)
|
|
self.assertEqual(r["rCount"], 835)
|
|
self.assertEqual(r["rMatchCount"], 835)
|
|
self.assertEqual(r["aCount"], 6267)
|
|
self.assertEqual(r["disAtmCount"], 0)
|
|
self.assertEqual(r["aCoordMatchCount"], 6267)
|
|
self.assertEqual(len(r["chains"]), 2)
|
|
self.assertTrue(r["pass"])
|
|
|
|
def test_rebuild_disordered_atoms_residues(self):
|
|
"""Convert disordered protein to internal coordinates and back."""
|
|
# 3jqh has both disordered residues
|
|
# and disordered atoms in ordered residues
|
|
|
|
with warnings.catch_warnings(record=True):
|
|
warnings.simplefilter("always", PDBConstructionWarning)
|
|
r = structure_rebuild_test(self.cif_3JQH, False)
|
|
# print(r)
|
|
self.assertEqual(r["residues"], 26)
|
|
self.assertEqual(r["rCount"], 47)
|
|
self.assertEqual(r["rMatchCount"], 47)
|
|
self.assertEqual(r["aCount"], 217)
|
|
self.assertEqual(r["disAtmCount"], 50)
|
|
self.assertEqual(r["aCoordMatchCount"], 217)
|
|
self.assertEqual(len(r["chains"]), 1)
|
|
self.assertTrue(r["pass"])
|
|
|
|
IC_Residue.no_altloc = True
|
|
with warnings.catch_warnings(record=True):
|
|
warnings.simplefilter("always", PDBConstructionWarning)
|
|
r = structure_rebuild_test(self.cif_3JQH2, verbose=False, quick=True)
|
|
self.assertEqual(r["aCoordMatchCount"], 167, msg="no_altloc fail")
|
|
self.assertTrue(r["pass"], msg="no_altloc fail")
|
|
IC_Residue.no_altloc = False
|
|
|
|
def test_no_crosstalk(self):
|
|
"""Deep copy, change few internal coords, test nothing else changes."""
|
|
# IC_Chain.ParallelAssembleResidues = False
|
|
self.cif_4CUP.atom_to_internal_coordinates()
|
|
cpy4cup = copy.deepcopy(self.cif_4CUP)
|
|
cic0 = self.cif_4CUP.child_list[0].child_list[0].internal_coord
|
|
cic1 = cpy4cup.child_list[0].child_list[0].internal_coord
|
|
alist = [
|
|
"omg",
|
|
"phi",
|
|
"psi",
|
|
"chi1",
|
|
"chi2",
|
|
"chi3",
|
|
"chi4",
|
|
"chi5",
|
|
"tau",
|
|
]
|
|
delta = 33 # degrees to change
|
|
tdelta = delta / 10.0 # more realistic for bond angle
|
|
targPos = 1
|
|
for ang in alist:
|
|
# skip by 2's with alist along original chain changing angle spec
|
|
ricTarg = cic0.chain.child_list[targPos].internal_coord
|
|
# print(targPos + 1, ricTarg.lc, ang)
|
|
targPos += 2
|
|
try:
|
|
edr = ricTarg.pick_angle(ang)
|
|
andx = edr.ndx
|
|
if ang == "tau":
|
|
cic0.hedraAngle[andx] += tdelta
|
|
cic0.hAtoms_needs_update[andx] = True
|
|
cic0.atomArrayValid[cic0.h2aa[andx]] = False
|
|
cic0.hAtoms_needs_update[:] = True
|
|
cic0.atomArrayValid[:] = False
|
|
cic0.dAtoms_needs_update[:] = True
|
|
else:
|
|
cic0.dihedraAngle[andx] += delta
|
|
if cic0.dihedraAngle[andx] > 180.0:
|
|
cic0.dihedraAngle[andx] -= 360.0
|
|
cic0.dihedraAngleRads[andx] = np.deg2rad(cic0.dihedraAngle[andx])
|
|
cic0.dAtoms_needs_update[andx] = True
|
|
cic0.atomArrayValid[cic0.d2aa[andx]] = False
|
|
# test Dihedron.bits()
|
|
pfd = IC_Residue.picFlagsDict
|
|
if ricTarg.rbase[2] == "P" and ang == "omg":
|
|
self.assertEqual(edr.bits(), (pfd["omg"] | pfd["pomg"]))
|
|
else:
|
|
self.assertEqual(edr.bits(), pfd[ang])
|
|
|
|
except AttributeError:
|
|
pass # skip if residue does not have e.g. chi5
|
|
cic0.internal_to_atom_coordinates() # move atoms
|
|
cic0.atom_to_internal_coordinates() # get new internal coords
|
|
# generate hdelta and ddelta difference arrays so can look for what
|
|
# changed
|
|
hdelta = cic0.hedraAngle - cic1.hedraAngle
|
|
hdelta[np.abs(hdelta) < 0.00001] = 0.0
|
|
ddelta = cic0.dihedraAngle - cic1.dihedraAngle
|
|
ddelta[np.abs(ddelta) < 0.00001] = 0.0
|
|
ddelta[ddelta < -180.0] += 360.0 # wrap around circle values
|
|
targPos = 1
|
|
for ang in alist:
|
|
# same skip along original chain looking at hdelta and ddelta
|
|
# if change is as specified, set difference to 0 then we can test
|
|
# for any remaining (spurious) changes
|
|
ricTarg = cic0.chain.child_list[targPos].internal_coord
|
|
# print(targPos + 1, ricTarg.lc, ang)
|
|
targPos += 2
|
|
try:
|
|
andx = ricTarg.pick_angle(ang).ndx
|
|
if ang == "tau":
|
|
self.assertAlmostEqual(hdelta[andx], tdelta, places=4)
|
|
hdelta[andx] = 0.0
|
|
# some other angle has to change to accommodate tau change
|
|
# N-Ca-Cb is artifact of choices in ic_data
|
|
# expected change so clear relevant hdelta here
|
|
adjAngNdx = ricTarg.pick_angle("N:CA:CB").ndx
|
|
self.assertNotAlmostEqual(hdelta[adjAngNdx], 0.0, places=1)
|
|
hdelta[adjAngNdx] = 0.0
|
|
else:
|
|
self.assertAlmostEqual(ddelta[andx], delta, places=4)
|
|
ddelta[andx] = 0.0
|
|
except AttributeError:
|
|
pass # if residue does not have e.g. chi5
|
|
|
|
hsum = hdelta.sum()
|
|
self.assertEqual(hsum, 0.0)
|
|
dsum = ddelta.sum()
|
|
self.assertEqual(dsum, 0.0)
|
|
|
|
# test hedron len12, angle, len23 setters and getters
|
|
hed = list(cic0.hedra.values())[10]
|
|
val = hed.len12 + 0.5
|
|
hed.len12 = val
|
|
self.assertEqual(hed.len12, val)
|
|
val = hed.len23 + 0.5
|
|
hed.len23 = val
|
|
self.assertEqual(hed.len23, val)
|
|
val = hed.angle + 1
|
|
hed.angle = val
|
|
self.assertEqual(hed.angle, val)
|
|
dihed = list(cic0.dihedra.values())[10]
|
|
val = dihed.angle + 196
|
|
dihed.angle = val
|
|
if val > 180.0:
|
|
val -= 360.0
|
|
if val < -180.0:
|
|
val += 360.0
|
|
self.assertEqual(dihed.angle, val)
|
|
|
|
def test_model_change_internal_coords(self):
|
|
"""Get model internal coords, modify psi and chi1 values and check."""
|
|
mdl = self.pdb_1LCD[1]
|
|
mdl.atom_to_internal_coordinates()
|
|
# other tests show can build with arbitrary internal coords
|
|
# build here so changes below trigger more complicated
|
|
# Atoms_needs_update mask arrays
|
|
mdl.internal_to_atom_coordinates()
|
|
nvt = {}
|
|
nvc1 = {}
|
|
nvpsi = {}
|
|
nvlen = {}
|
|
nvlen2 = {}
|
|
tcount = 0
|
|
l2count = 0
|
|
c1count = 0
|
|
psicount = 0
|
|
lcount = 0
|
|
for r in mdl.get_residues():
|
|
ric = r.internal_coord
|
|
if ric:
|
|
# hedra change
|
|
tau = ric.get_angle("tau")
|
|
if ric.rprev != [] and tau is not None:
|
|
tcount += 1
|
|
nv = tau + 0.5
|
|
ric.set_angle("tau", nv)
|
|
nvt[str(r)] = nv
|
|
l2count += 1
|
|
leng2 = ric.get_length("N:CA")
|
|
nv = leng2 + 0.05
|
|
ric.set_length("N:CA", nv)
|
|
nvlen2[str(r)] = nv
|
|
|
|
# sidechain dihedron change
|
|
chi1 = ric.get_angle("chi1")
|
|
if chi1 is not None:
|
|
c1count += 1
|
|
nv = chi1 + 90
|
|
if nv > 180.0:
|
|
nv -= 360.0
|
|
# ric.set_angle("chi1", nv)
|
|
ric.bond_set("chi1", nv)
|
|
nvc1[str(r)] = nv
|
|
# backbone dihedron change
|
|
psi = ric.get_angle("psi")
|
|
if psi is not None:
|
|
psicount += 1
|
|
nv = psi - 90
|
|
if nv < -180.0:
|
|
nv += 360.0
|
|
ric.set_angle("psi", nv)
|
|
nvpsi[str(r)] = nv
|
|
leng = ric.get_length("CA:CB")
|
|
if leng is not None:
|
|
lcount += 1
|
|
nv = leng + 0.05
|
|
ric.set_length("CA:CB", nv)
|
|
nvlen[str(r)] = nv
|
|
|
|
mdl.internal_to_atom_coordinates()
|
|
|
|
# prove not using stored results
|
|
for chn in mdl.get_chains():
|
|
if hasattr(chn, "hedraLen"):
|
|
delattr(chn.internal_coord, "hedraLen")
|
|
delattr(chn.internal_coord, "dihedraLen")
|
|
delattr(chn.internal_coord, "hedraAngle")
|
|
delattr(chn.internal_coord, "dihedraAngle")
|
|
for r in chn.get_residues():
|
|
r.internal_coord.hedra = {}
|
|
r.internal_coord.dihedra = {}
|
|
|
|
mdl.atom_to_internal_coordinates()
|
|
ttcount = 0
|
|
l2tcount = 0
|
|
c1tcount = 0
|
|
psitcount = 0
|
|
ltcount = 0
|
|
for r in mdl.get_residues():
|
|
ric = r.internal_coord
|
|
if ric:
|
|
tau = ric.get_angle("tau")
|
|
if ric.rprev != [] and tau is not None:
|
|
ttcount += 1
|
|
# print(str(r), "tau", tau, nvt[str(r)])
|
|
self.assertAlmostEqual(tau, nvt[str(r)], places=3)
|
|
l2tcount += 1
|
|
l2 = ric.get_length("N:CA")
|
|
self.assertAlmostEqual(l2, nvlen2[str(r)], places=3)
|
|
chi1 = ric.get_angle("chi1")
|
|
if chi1 is not None:
|
|
c1tcount += 1
|
|
# print(str(r), "chi1", chi1, nvc1[str(r)])
|
|
self.assertAlmostEqual(chi1, nvc1[str(r)], places=3)
|
|
psi = ric.get_angle("psi")
|
|
if psi is not None:
|
|
psitcount += 1
|
|
# print(str(r), "psi", psi, nvpsi[str(r)])
|
|
self.assertAlmostEqual(psi, nvpsi[str(r)], places=3)
|
|
leng = ric.get_length("CA:CB")
|
|
if leng is not None:
|
|
ltcount += 1
|
|
self.assertAlmostEqual(leng, nvlen[str(r)], places=3)
|
|
|
|
self.assertEqual(tcount, ttcount)
|
|
self.assertEqual(l2count, l2tcount)
|
|
self.assertEqual(c1count, c1tcount)
|
|
self.assertEqual(psicount, psitcount)
|
|
self.assertEqual(lcount, ltcount)
|
|
self.assertGreater(ttcount, 0)
|
|
self.assertGreater(c1count, 0)
|
|
self.assertGreater(psicount, 0)
|
|
self.assertGreater(lcount, 0)
|
|
|
|
def test_write_SCAD(self):
|
|
"""Check SCAD output plus MaxPeptideBond and Gly CB.
|
|
|
|
SCAD tests: scaling, transform mtx, extra bond created (allBonds)
|
|
"""
|
|
sf = StringIO()
|
|
write_SCAD(
|
|
self.cif_4CUP2,
|
|
sf,
|
|
10.0,
|
|
pdbid="4cup",
|
|
backboneOnly=True,
|
|
includeCode=False,
|
|
)
|
|
sf.seek(0)
|
|
next_one = False
|
|
with as_handle(sf, mode="r") as handle:
|
|
for aline in handle.readlines():
|
|
if "// (1856_S_CB, 1856_S_CA, 1856_S_C)" in aline:
|
|
m = re.search(r"\[\s+(\d+\.\d+)\,", aline)
|
|
if m:
|
|
# test correctly scaled atom bond length
|
|
self.assertAlmostEqual(float(m.group(1)), 15.30582, places=3)
|
|
else:
|
|
self.fail("scaled atom bond length not found")
|
|
elif '[ 1, "1857M",' in aline:
|
|
next_one = True
|
|
elif next_one:
|
|
next_one = False
|
|
# test last residue transform looks roughly correct
|
|
# some differences due to sorting issues on different
|
|
# python versions
|
|
target = [-12.413, -3.303, 35.771, 1.0]
|
|
ms = re.findall( # last column of each row
|
|
r"\s+(-?\d+\.\d+)\s+\]", aline
|
|
)
|
|
if ms:
|
|
for i in range(3):
|
|
self.assertAlmostEqual(float(ms[i]), target[i], places=0)
|
|
else:
|
|
self.fail("transform not found")
|
|
|
|
sf.seek(0)
|
|
IC_Residue.gly_Cbeta = True
|
|
IC_Chain.MaxPeptideBond = 100.0
|
|
chn = self.pdb_2XHE2[0]["A"]
|
|
chn.atom_to_internal_coordinates()
|
|
rt0 = chn.internal_coord.ordered_aa_ic_list[12]
|
|
rt1 = chn.internal_coord.ordered_aa_ic_list[16]
|
|
rt0.set_flexible()
|
|
rt1.set_hbond()
|
|
|
|
write_SCAD(
|
|
self.pdb_2XHE2[0]["A"],
|
|
sf,
|
|
10.0,
|
|
pdbid="2xhe",
|
|
# maxPeptideBond=100.0,
|
|
includeCode=False,
|
|
start=10,
|
|
fin=570,
|
|
)
|
|
sf.seek(0)
|
|
allBondsPass = False
|
|
maxPeptideBondPass = False
|
|
glyCbetaFound = False
|
|
startPass = True
|
|
finPass = True
|
|
flexPass = False
|
|
hbPass = False
|
|
with as_handle(sf, mode="r") as handle:
|
|
for aline in handle.readlines():
|
|
# test extra bond created in TRP (allBonds is True)
|
|
if '"Cres", 0, 0, 1, 0, StdBond, "W", 24, "CD2CE3CZ3"' in aline:
|
|
allBondsPass = True
|
|
# test 509_K-561_E long bond created
|
|
if "509_K" in aline and "561_E" in aline:
|
|
maxPeptideBondPass = True
|
|
if "(21_G_CB, 21_G_CA, 21_G_C)" in aline:
|
|
glyCbetaFound = True
|
|
target = [15.33630, 110.17513, 15.13861]
|
|
ms = re.findall(r"\s+(-?\d+\.\d+)", aline)
|
|
if ms:
|
|
for i in range(3):
|
|
self.assertAlmostEqual(float(ms[i]), target[i], places=0)
|
|
else:
|
|
self.fail("Cbeta internal coords not found")
|
|
if "8_K_CA" in aline:
|
|
startPass = False
|
|
if "572_N_CA" in aline:
|
|
finPass = False
|
|
if 'FemaleJoinBond, FemaleJoinBond, "N", 13, "NCAC"' in aline:
|
|
flexPass = True
|
|
if 'HBond, "R", 16, "CACO"' in aline:
|
|
hbPass = True
|
|
|
|
self.assertTrue(allBondsPass, msg="missing extra ring close bonds")
|
|
self.assertTrue(glyCbetaFound, msg="gly CB not created")
|
|
self.assertTrue(maxPeptideBondPass, msg="ignored maxPeptideBond setting")
|
|
self.assertTrue(startPass, msg="writeSCAD wrote residue before start")
|
|
self.assertTrue(finPass, msg="writeSCAD wrote residue past fin")
|
|
self.assertTrue(flexPass, msg="writeSCAD residue 12 not flexible")
|
|
self.assertTrue(hbPass, msg="writeSCAD residue 16 no hbond")
|
|
|
|
def test_i2a_start_fin(self):
|
|
"""Test assemble start/fin, default NCaC coordinates, IC_duplicate."""
|
|
chn = self.pdb_1LCD[2]["A"]
|
|
cpy = IC_duplicate(chn)[2]["A"] # generates internal coords as needed
|
|
cpy.internal_to_atom_coordinates(start=31, fin=45)
|
|
cdict = compare_residues(chn, cpy, quick=True)
|
|
self.assertFalse(cdict["pass"])
|
|
# transform source coordinates to put res 31 tau at origin like
|
|
# fragment
|
|
res = chn[31]
|
|
psi = res.internal_coord.pick_angle("psi")
|
|
cst = np.transpose(psi.cst)
|
|
chn.internal_coord.atomArray[:] = chn.internal_coord.atomArray.dot(cst)
|
|
cdict = compare_residues(chn, cpy, rtol=1e-03, atol=1e-05)
|
|
self.assertEqual(cdict["residues"], 51)
|
|
self.assertEqual(cdict["rMatchCount"], 77)
|
|
self.assertEqual(cdict["aCount"], 497)
|
|
self.assertEqual(cdict["disAtmCount"], 0)
|
|
self.assertEqual(cdict["aCoordMatchCount"], 140)
|
|
self.assertEqual(cdict["aFullIdMatchCount"], 140)
|
|
self.assertEqual(len(cdict["chains"]), 1)
|
|
self.assertEqual(cdict["rCount"], 77)
|
|
self.assertFalse(cdict["pass"])
|
|
|
|
def test_distplot_rebuild(self):
|
|
"""Build identical structure from distplot and chirality data."""
|
|
# load input chain
|
|
for _chn1 in self.cif_4ZHL.get_chains():
|
|
break
|
|
# create atomArray and compute distplot and dihedral signs array
|
|
_chn1.atom_to_internal_coordinates()
|
|
_c1ic = _chn1.internal_coord
|
|
atmNameNdx = AtomKey.fields.atm
|
|
CaSelect = [
|
|
_c1ic.atomArrayIndex.get(k)
|
|
for k in _c1ic.atomArrayIndex.keys()
|
|
if k.akl[atmNameNdx] == "CA"
|
|
]
|
|
dplot0 = _chn1.internal_coord.distance_plot(filter=CaSelect)
|
|
self.assertAlmostEqual(
|
|
dplot0[3, 9],
|
|
16.296,
|
|
places=3,
|
|
msg="fail generate distance plot with filter",
|
|
)
|
|
dplot1 = _chn1.internal_coord.distance_plot()
|
|
dsigns = _chn1.internal_coord.dihedral_signs()
|
|
|
|
# load second copy (same again) input chain
|
|
for _chn2 in self.cif_4ZHL2.get_chains():
|
|
break
|
|
# create internal coord structures but do not compute di/hedra
|
|
cic2 = _chn2.internal_coord = IC_Chain(_chn2)
|
|
cic2.init_edra()
|
|
# load relevant interatomic distances from chn1 distance plot
|
|
cic2.distplot_to_dh_arrays(dplot1, dsigns)
|
|
# compute di/hedra angles from dh_arrays
|
|
cic2.distance_to_internal_coordinates()
|
|
|
|
# clear chn2 atom coordinates
|
|
# cic2.atomArrayValid[:] = False # done in distance_to_internal_coordinates
|
|
|
|
# initialize values but this is redundant to Valid=False above
|
|
cic2.atomArray = np.zeros((cic2.AAsiz, 4), dtype=np.float64)
|
|
cic2.atomArray[:, 3] = 1.0
|
|
|
|
# 4zhl has chain breaks so copy initial coords of each segment
|
|
cic2.copy_initNCaCs(_chn1.internal_coord)
|
|
# compute chn2 atom coords from di/hedra data
|
|
cic2.internal_to_atom_coordinates()
|
|
|
|
# generate distance plot from second chain, confirm minimal distance
|
|
# from original
|
|
dp2 = cic2.distance_plot()
|
|
dpdiff = np.abs(dplot1 - dp2)
|
|
# print(np.amax(dpdiff))
|
|
self.assertTrue(np.amax(dpdiff) < 0.000001)
|
|
|
|
def test_seq_as_PIC(self):
|
|
"""Read seq as PIC data, extend chain, set each chi angle, check various."""
|
|
pdb_structure = read_PIC_seq(
|
|
SeqRecord(
|
|
Seq("GAVLIMFPSTCNQYWDEHKR"),
|
|
id="1RND",
|
|
description="my 20aa sequence",
|
|
)
|
|
)
|
|
chn = next(pdb_structure.get_chains())
|
|
cic = chn.internal_coord
|
|
cic.make_extended()
|
|
for chi in range(1, 6):
|
|
for ric in chn.internal_coord.ordered_aa_ic_list:
|
|
targ = "chi" + str(chi)
|
|
rchi = ric.get_angle(targ)
|
|
if rchi is not None:
|
|
nangl = rchi + 60.0
|
|
# nangl = nangl if (nangl <= 180.0) else nangl - 360.0
|
|
ric.set_angle(targ, nangl)
|
|
# ric.bond_set("chi1", nangl)
|
|
|
|
pdb_structure.internal_to_atom_coordinates()
|
|
|
|
self.assertEqual(
|
|
len(cic.atomArrayValid), 168, msg="wrong number atoms from sequence"
|
|
)
|
|
|
|
# test make_extended and each sequential chi rotation places selected
|
|
# atom from each sidechain where expected.
|
|
posnDict = {
|
|
"1_G_CA": [0.000, 0.000, 0.000, 1.000],
|
|
"2_A_CB": [-0.411, 2.505, 4.035, 1.000],
|
|
"3_V_CG2": [-2.246, -2.942, 4.865, 1.000],
|
|
"4_L_CD2": [-4.673, 2.124, 6.060, 1.000],
|
|
"5_I_CD1": [-4.260, -4.000, 10.546, 1.000],
|
|
"6_M_CE": [-10.035, 4.220, 12.322, 1.000],
|
|
"7_F_CE2": [-5.835, -1.076, 15.391, 1.000],
|
|
"8_P_CD": [-15.060, -2.368, 17.751, 1.000],
|
|
"9_S_OG": [-12.733, 0.388, 24.504, 1.000],
|
|
"10_T_CG2": [-19.084, -0.423, 22.650, 1.000],
|
|
"11_C_SG": [-15.455, 2.166, 27.172, 1.000],
|
|
"12_N_ND2": [-20.734, -3.990, 27.233, 1.000],
|
|
"13_Q_NE2": [-19.362, 4.465, 30.847, 1.000],
|
|
"14_Y_CE2": [-20.685, -3.571, 32.854, 1.000],
|
|
"15_W_CH2": [-23.295, 1.706, 32.883, 1.000],
|
|
"16_D_OD2": [-24.285, -3.683, 40.620, 1.000],
|
|
"17_E_OE2": [-31.719, 2.385, 38.887, 1.000],
|
|
"18_H_NE2": [-25.763, -0.071, 46.356, 1.000],
|
|
"19_K_NZ": [-36.073, -0.926, 42.383, 1.000],
|
|
"20_R_NH2": [-28.792, 3.687, 51.738, 1.000],
|
|
}
|
|
|
|
for k, v in posnDict.items():
|
|
atm = AtomKey(k)
|
|
ndx = cic.atomArrayIndex[atm]
|
|
coord = np.round(cic.atomArray[ndx], 3)
|
|
self.assertTrue(np.array_equal(coord, v), msg=f"position error on atom {k}")
|
|
|
|
cic.update_dCoordSpace()
|
|
rt = cic.ordered_aa_ic_list[5] # pick a residue
|
|
chi1 = rt.pick_angle("chi1") # chi1 coord space puts CA at origin
|
|
rt.applyMtx(chi1.cst)
|
|
coord = rt.residue.child_dict["CA"].coord # Biopython API Atom coords
|
|
self.assertTrue(
|
|
np.allclose(coord, [0.0, 0.0, 0.0]), msg="dCoordSpace transform error"
|
|
)
|
|
|
|
# test Dihedron repr and all that leads into it
|
|
psi = rt.pick_angle("psi")
|
|
|
|
self.assertEqual(
|
|
psi.__repr__(),
|
|
"4-6_M_N:6_M_CA:6_M_C:7_F_N MNMCAMCFN 123.0 ('1RND', 0, 'A', (' ', 6, ' '))",
|
|
msg="dihedron __repr__ error for M6 psi",
|
|
)
|
|
|
|
m = "Edron rich comparison failed"
|
|
self.assertTrue(chi1 != psi, msg=m)
|
|
self.assertFalse(chi1 == psi, msg=m)
|
|
self.assertTrue(psi < chi1, msg=m)
|
|
self.assertTrue(psi <= chi1, msg=m)
|
|
self.assertTrue(chi1 > psi, msg=m)
|
|
self.assertTrue(chi1 >= psi, msg=m)
|
|
|
|
self.assertTrue(
|
|
chi1.cre_class == "NsbCsbCsbCsb",
|
|
msg="covalent radii assignment error for chi1",
|
|
)
|
|
|
|
# dihedron atomkeys are all in residue atomkeys as expected, including
|
|
# i+1 N for psi.
|
|
self.assertTrue(all(ak in rt for ak in chi1.atomkeys))
|
|
self.assertFalse(all(ak in rt for ak in psi.atomkeys))
|
|
|
|
# test Hedron repr and all that leads into it
|
|
tau = rt.pick_angle("tau")
|
|
self.assertEqual(
|
|
tau.__repr__(),
|
|
"3-6_M_N:6_M_CA:6_M_C MNMCAMC 1.46091 110.97184 1.52499",
|
|
msg="hedron __repr__ error for M11 tau",
|
|
)
|
|
|
|
# some specific AtomKey comparisons missed in other tests
|
|
a0, a1 = tau.atomkeys[0], tau.atomkeys[1]
|
|
m = "AtomKey rich comparison failed"
|
|
self.assertTrue(a1 > a0, msg=m)
|
|
self.assertTrue(a1 >= a0, msg=m)
|
|
self.assertTrue(a0 <= a1, msg=m)
|
|
|
|
def test_angle_fns(self):
|
|
"""Test angle_dif and angle_avg across +/-180 boundaries."""
|
|
arr1 = np.array([179.0, 90.0, 88.0, 1.0])
|
|
arr2 = np.array([-179.0, -90.0, -91.0, -1.0])
|
|
assert (
|
|
Dihedron.angle_dif(arr1, arr2) == np.array([2.0, 180.0, -179.0, -2.0])
|
|
).all()
|
|
assert Dihedron.angle_avg(np.array([179.0, -179.0])) == 180.0
|
|
assert Dihedron.angle_avg(np.array([1.0, -1.0])) == 0.0
|
|
assert Dihedron.angle_avg(np.array([90.0, -90.0])) == 0.0
|
|
assert Dihedron.angle_avg(np.array([91.0, -91.0])) == 180.0
|
|
|
|
|
|
if __name__ == "__main__":
|
|
runner = unittest.TextTestRunner(verbosity=2)
|
|
unittest.main(testRunner=runner)
|