Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def test_connect_via_residue_names(single_model):
"""
Test whether the created bond list is equal to the bonds deposited
in the MMTF file.
"""
# Structure with peptide, nucleotide, small molecules and water
file = mmtf.MMTFFile.read(join(data_dir("structure"), "5ugo.mmtf"))
if single_model:
atoms = mmtf.get_structure(file, include_bonds=True, model=1)
else:
atoms = mmtf.get_structure(file, include_bonds=True)
ref_bonds = atoms.bonds
test_bonds = struc.connect_via_residue_names(atoms)
assert test_bonds == ref_bonds
def test_array_conversion(path, single_model):
model = 1 if single_model else None
mmtf_file = mmtf.MMTFFile.read(path)
a1 = mmtf.get_structure(mmtf_file, model=model, include_bonds=True)
mmtf_file = mmtf.MMTFFile()
mmtf.set_structure(mmtf_file, a1)
temp = TemporaryFile("w+b")
mmtf_file.write(temp)
temp.seek(0)
mmtf_file = mmtf.MMTFFile.read(temp)
temp.close()
a2 = mmtf.get_structure(mmtf_file, model=model, include_bonds=True)
for category in a1.get_annotation_categories():
assert a1.get_annotation(category).tolist() == \
a2.get_annotation(category).tolist()
assert a1.coord.flatten().tolist() == \
approx(a2.coord.flatten().tolist(), abs=1e-3)
def test_dihedral_backbone_result(file_name):
import mdtraj
mmtf_file = mmtf.MMTFFile.read(file_name)
array = mmtf.get_structure(mmtf_file, model=1)
array = array[struc.filter_amino_acids(array)]
if array.array_length() == 0:
# Structure contains no protein
# -> determination of backbone angles makes no sense
return
for chain in struc.chain_iter(array):
print("Chain: ", chain.chain_id[0])
if len(struc.check_res_id_continuity(chain)) != 0:
# Do not test discontinuous chains
return
test_phi, test_psi, test_ome = struc.dihedral_backbone(chain)
temp = NamedTemporaryFile("w+", suffix=".pdb")
strucio.save_structure(temp.name, chain)
traj = mdtraj.load(temp.name)
def test_fetch(format, as_file_like):
path = None if as_file_like else tempfile.gettempdir()
file_path_or_obj = rcsb.fetch("1l2y", format, path, overwrite=True)
if format == "pdb":
file = pdb.PDBFile.read(file_path_or_obj)
pdb.get_structure(file)
elif format == "pdbx":
file = pdbx.PDBxFile.read(file_path_or_obj)
pdbx.get_structure(file)
elif format == "mmtf":
file = mmtf.MMTFFile.read(file_path_or_obj)
mmtf.get_structure(file)
elif format == "fasta":
file = fasta.FastaFile.read(file_path_or_obj)
# Test if the file contains any sequences
assert len(fasta.get_sequences(file)) > 0
# Converter for the DSSP secondary structure elements
# to the classical ones
dssp_to_abc = {"I" : "c",
"S" : "c",
"H" : "a",
"E" : "b",
"G" : "c",
"B" : "b",
"T" : "c",
"C" : "c"}
# Fetch and load structure
file_name = rcsb.fetch("1QGD", "mmtf", gettempdir())
mmtf_file = mmtf.MMTFFile.read(file_name)
array = mmtf.get_structure(mmtf_file, model=1)
# Transketolase homodimer
tk_dimer = array[struc.filter_amino_acids(array)]
# Transketolase monomer
tk_mono = tk_dimer[tk_dimer.chain_id == "A"]
# The chain ID corresponding to each residue
chain_id_per_res = array.chain_id[struc.get_residue_starts(tk_dimer)]
sse = mmtf_file["secStructList"]
sse = sse[sse != -1]
sse = sse[chain_id_per_res == "A"]
sse = np.array([sec_struct_codes[code] for code in sse if code != -1],
dtype="U1")
sse = np.array([dssp_to_abc[e] for e in sse], dtype="U1")
# Helper function to convert secondary structure array to annotation
# and visualize it
"""
# Code source: Patrick Kunzmann
# License: BSD 3 clause
import numpy as np
import matplotlib.pyplot as plt
import biotite.structure as struc
import biotite.structure.io.mmtf as mmtf
import biotite.structure.graphics as graphics
import biotite.database.rcsb as rcsb
# Structure of a DNA double helix
mmtf_file = mmtf.MMTFFile.read(rcsb.fetch("1qxb", "mmtf"))
structure = mmtf.get_structure(mmtf_file, model=1, include_bonds=True)
nucleotides = structure[struc.filter_nucleotides(structure)]
# Choose one adenine-thymine and one guanine-cytosine base pair
base_pairs = struc.base_pairs(nucleotides)
for i, j in base_pairs:
if (nucleotides.res_name[i], nucleotides.res_name[j]) == ("DG", "DC"):
guanine, cytosine = [nucleotides[mask] for mask
in struc.get_residue_masks(nucleotides, [i, j])]
break
for i, j in base_pairs:
if (nucleotides.res_name[i], nucleotides.res_name[j]) == ("DA", "DT"):
adenine, thymine = [nucleotides[mask] for mask
in struc.get_residue_masks(nucleotides, [i, j])]
break
pairs = [(guanine, cytosine), (adenine, thymine)]
# The corresponding structure
PDB_ID = "1MUG"
# The normal mode to be visualized
# '-1' is the last (and most significant) one
MODE = -1
# The amount of frames (models) per oscillation
FRAMES = 60
# The maximum oscillation amplitude for an atom
# (The length of the ANM's eigenvectors make only sense when compared
# relative to each other, the absolute values have no significance)
MAX_AMPLITUDE = 5
# Load structure
mmtf_file = mmtf.MMTFFile.read(rcsb.fetch(PDB_ID, "mmtf"))
structure = mmtf.get_structure(mmtf_file, model=1)
# Filter first peptide chain
protein_chain = structure[
struc.filter_amino_acids(structure)
& (structure.chain_id == structure.chain_id[0])
]
# Filter CA atoms
ca = protein_chain[protein_chain.atom_name == "CA"]
# Load eigenvectors for CA atoms
# The first axis indicates the mode,
# the second axis indicates the vector component
vectors = np.loadtxt(VECTOR_FILE, delimiter=",").transpose()
# Discard the last 6 modes, as these are movements of the entire system:
########################################################################
# Structure manipulation
# ----------------------
#
# The most basic way to manipulate a structure is to edit the
# annotation arrays or coordinates directly.
from tempfile import gettempdir
import biotite.database.rcsb as rcsb
import biotite.structure as struc
import biotite.structure.io.mmtf as mmtf
file_path = rcsb.fetch("1l2y", "mmtf", gettempdir())
mmtf_file = mmtf.MMTFFile.read(file_path)
structure = mmtf.get_structure(mmtf_file, model=1)
print("Before:")
print(structure[structure.res_id == 1])
print()
structure.coord += 100
print("After:")
print(structure[structure.res_id == 1])
########################################################################
# *Biotite* provides also some transformation functions, for example
# :func:`rotate()` for rotations about the *x*-, *y*- or *z*-axis.
structure = mmtf.get_structure(mmtf_file, model=1)
print("Before:")
print(structure[structure.res_id == 1])
print()
# Rotation about z-axis by 90 degrees
])
# Fetch animal lysoyzme structures
lyso_files = rcsb.fetch(
["1REX", "1AKI", "1DKJ", "1GD6"],
format="mmtf", target_path=biotite.temp_dir()
)
organisms = ["H. sapiens", "G. gallus", "C. viginianus", "B. mori"]
# Create a PB sequence from each structure
pb_seqs = []
for file_name in lyso_files:
file = mmtf.MMTFFile.read(file_name)
# Take only the first model into account
array = mmtf.get_structure(file, model=1)
# Remove everything but the first protein chain
array = array[struc.filter_amino_acids(array)]
array = array[array.chain_id == array.chain_id[0]]
# Calculate backbone dihedral angles,
# as the PBs are determined from them
phi, psi, omega = struc.dihedral_backbone(array)
# A PB requires the 8 phi/psi angles of 5 amino acids,
# centered on the amino acid to calculate the PB for
# Hence, the PBs are not defined for the two amino acids
# at each terminus
pb_angles = np.full((len(phi)-4, 8), np.nan)
pb_angles[:, 0] = psi[ : -4]
pb_angles[:, 1] = phi[1 : -3]
pb_angles[:, 2] = psi[1 : -3]
pb_angles[:, 3] = phi[2 : -2]