Welcome to Lahuta¶

Download a few test files¶

To get started, let's download a few files. We download a utility function that will handle the download for us.

In [2]:

from lahuta.api import download_structures

from lahuta.api import download_structures

In [4]:

download_structures(
    ['2RH1', '1KX2'], 
    pdb_or_cif='pdb', 
    dir_loc='data'
)

download_structures( ['2RH1', '1KX2'], pdb_or_cif='pdb', dir_loc='data' )

Out[4]:

{'2RH1': 'data/2rh1.pdb', '1KX2': 'data/1kx2.pdb'}

A simple, yet feature-rich, entry point¶

Luni serves as the entry point to Lahuta.

In [4]:

from lahuta import Luni

from lahuta import Luni

In [5]:

luni = Luni('data/2rh1.pdb')

luni = Luni('data/2rh1.pdb')

Acessing Atom, Residue, and Chain (ARC) objects¶

In [6]:

luni.arc.atoms, luni.arc.residues, luni.arc.chains

luni.arc.atoms, luni.arc.residues, luni.arc.chains

Out[6]:

(<lahuta.core.arc.Atoms at 0x7fd182a5d7b0>,
 <lahuta.core.arc.Residues at 0x7fd182a5df60>,
 <lahuta.core.arc.Chains at 0x7fd182a5e5c0>)

Accessing Atoms¶

In [7]:

atoms = luni.arc.atoms

atoms = luni.arc.atoms

In [8]:

names, elements, coords = atoms.names, atoms.elements, atoms.coordinates
names, elements, coords

names, elements, coords = atoms.names, atoms.elements, atoms.coordinates names, elements, coords

Out[8]:

(array(['N', 'CA', 'C', ..., 'O', 'O', 'O'], dtype='<U10'),
 array(['N', 'C', 'C', ..., 'O', 'O', 'O'], dtype='<U10'),
 array([[-52.822,  -1.611,  23.137],
        [-51.922,  -2.262,  22.148],
        [-52.178,  -1.713,  20.742],
        ...,
        [-37.397,  27.672,   6.12 ],
        [-19.199,  49.833,  44.41 ],
        [-35.618,  21.645,   8.827]], dtype=float32))

Accessing Residues¶

In [9]:

residues = luni.arc.residues

residues = luni.arc.residues

In [10]:

resids, resnames = residues.resids, residues.resnames
resids, resnames

resids, resnames = residues.resids, residues.resnames resids, resnames

Out[10]:

(array([ 29,  29,  29, ..., 546, 547, 548]),
 array(['ASP', 'ASP', 'ASP', ..., 'HOH', 'HOH', 'HOH'], dtype='<U10'))

Accessing Chains¶

In [11]:

chains = luni.arc.chains

chains = luni.arc.chains

In [12]:

auths, labels = chains.auths, chains.labels
auths, labels

auths, labels = chains.auths, chains.labels auths, labels

Out[12]:

(array(['A', 'A', 'A', ..., 'A', 'A', 'A'], dtype='<U10'),
 array(['Axp', 'Axp', 'Axp', ..., 'Axw', 'Axw', 'Axw'], dtype='<U10'))

Accessing the Underlying Objects¶

Lahuta relies on three libraries: MDAnalysis, Gemmi, and OpenBabel.

• MDAnalysis is used for loading data from MD simulations.
• Gemmi is used for parsing PDBx/mmCIF files.
• OpenBabel is used for SMARTS pattern matching and perception of chemical properties (bonds, aromaticity, etc.).

MDAnalysis AtomGroup and Universe¶

In [13]:

ag = luni.to("mda")
ag, ag.universe

ag = luni.to("mda") ag, ag.universe

Out[13]:

(<AtomGroup with 3804 atoms>, <Universe with 3804 atoms>)

OpenBabel OBMol object¶

Note that OpenBabel does not support vectorized creation of its objects, so this may be slow for large systems.

In [14]:

mol = luni.to("mol")
mol.NumAtoms()

mol = luni.to("mol") mol.NumAtoms()

Out[14]:

3804

Computing Neighbors¶

A "neighbor" refers to an atom placed within a defined proximity to another atom.

In Lahuta, a neighbor is determined by a user-defined distance criterion from another atom.
If one atom is within this specified distance, it's considered a neighbor. Lahuta captures this information in its NeighborPairs class

In [15]:

ns = luni.compute_neighbors(
    radius=5.0, # The cutoff distance for the neighbor search
    res_dif=2,  # The minimum residue difference between the atoms
)
ns

ns = luni.compute_neighbors( radius=5.0, # The cutoff distance for the neighbor search res_dif=2, # The minimum residue difference between the atoms ) ns

Out[15]:

<Lahuta NeighborPairs class containing 3648 atoms and 17616 pairs>

In [16]:

# atom indices that are neighbors given the cutoff distance and residue difference along with the distances
ns.pairs, ns.distances

# atom indices that are neighbors given the cutoff distance and residue difference along with the distances ns.pairs, ns.distances

Out[16]:

(array([[   0,   21],
        [   0,   25],
        [   1,   21],
        ...,
        [3787, 3803],
        [3789, 3803],
        [3797, 3800]]),
 array([4.61659716, 4.5690032 , 4.99562331, ..., 3.54749682, 2.78203339,
        3.3515611 ]))

In [17]:

# parner1, partner2 are the indices of the atoms in the pair (column 0 and 1)
ns.partner1, ns.partner2

# parner1, partner2 are the indices of the atoms in the pair (column 0 and 1) ns.partner1, ns.partner2

Out[17]:

(<AtomGroup with 17616 atoms>, <AtomGroup with 17616 atoms>)

Indexing¶

In [18]:

ns[0], ns[10:30]

ns[0], ns[10:30]

Out[18]:

(<Lahuta NeighborPairs class containing 2 atoms and 1 pairs>,
 <Lahuta NeighborPairs class containing 16 atoms and 20 pairs>)

Neighbors vs Contacts¶

In a sense, neighbor atoms can be considered to already be in contact with each other. Usually, however, we are interested in specific types of contacts or contacts that satisfy certain criteria. These criteria can be geometric, chemical, or both. For example, we might be interested in contacts between ionic or aromatic atoms or contacts between atoms that are within a certain distance of each other and satisfy a certain angle constraint. We may also want to apply restrictions on the difference between the Van der Waals radii of the atoms (e.g. avoid contacts between atoms that are too close to each other - this is useful for avoiding clashes).

In [19]:

from lahuta.contacts import F

from lahuta.contacts import F

In [20]:

# Compute aromatic interactions
F.aromatic_neighbors(ns)

# Compute aromatic interactions F.aromatic_neighbors(ns)

Out[20]:

<Lahuta NeighborPairs class containing 95 atoms and 91 pairs>

In [21]:

# Compute ionic interactions
F.ionic_neighbors(ns)

# Compute ionic interactions F.ionic_neighbors(ns)

Out[21]:

<Lahuta NeighborPairs class containing 65 atoms and 58 pairs>

In [23]:

# Compute hydrogen bonds
F.hbond_neighbors(ns)

# Compute hydrogen bonds F.hbond_neighbors(ns)

Out[23]:

<Lahuta NeighborPairs class containing 0 atoms and 0 pairs>

In [24]:

# Compute van der Waals interactions
F.vdw_neighbors(ns)

# Compute van der Waals interactions F.vdw_neighbors(ns)

Out[24]:

<Lahuta NeighborPairs class containing 492 atoms and 282 pairs>

In [25]:

hphobid_default_distance = F.hydrophobic_neighbors(ns)
hphobid_distance_6 = F.hydrophobic_neighbors(ns, distance=6.0)

hphobid_default_distance, hphobid_distance_6

hphobid_default_distance = F.hydrophobic_neighbors(ns) hphobid_distance_6 = F.hydrophobic_neighbors(ns, distance=6.0) hphobid_default_distance, hphobid_distance_6

Out[25]:

(<Lahuta NeighborPairs class containing 865 atoms and 1203 pairs>,
 <Lahuta NeighborPairs class containing 1049 atoms and 2297 pairs>)