In [1]:
%matplotlib inline
import numpy as np
import bomeba0 as bmb
import matplotlib.pyplot as plt
The first thing we do to use bomeba is to create a protein object. We can do this by passing a sequence (using one-letter code). At the moment only the 20 most common amino-acids are accepted.
In [2]:
prot = bmb.Protein('GAG')
A protein object has several attributes for example we can retrieve the sequence
In [3]:
prot.sequence
Out[3]:
or the cartesian coordinates as a (N, 3) numpy array
In [4]:
prot.coords
Out[4]:
We can also compute the energy of our protein. At the moment the energy is computed in vacuum, with a very simple force-field that only includes a Lenard-Jones term:
$$LJ_{ij} = \epsilon \left [ \left (\frac{\sigma_{ij}}{r_{ij}} \right)^{12} - 2 \left (\frac{\sigma_{ij}}{r_{ij}} \right)^{6} \right]$$$\sigma_{ij}$ is the distance at which the potential reaches its minimum
$\epsilon_{ij}$ is the depth of the potential well
$r_{ij}$ is the distance between the atoms
Only five atoms types are considered (C, H, O, N, S), the parameters are loosely based on the general purpose force field GAFF
In [5]:
prot.energy()
Out[5]:
We can also compute the radius of gyration, i.e. the root mean squared deviation of all atoms with respect to their common geometric center.
In [6]:
prot.rgyr()
Out[6]:
As with other Python object we can ask for the lenght of our protein. The answer we get is the number of amino acids.
In [7]:
len(prot)
Out[7]:
A protein object is an iterable object made of aminoacids:
In [8]:
for res in prot:
print(res)
Residues has several properties, including coordinates. One convenient way to access coordinates is just by indexing a residue with a number (the atom number inside the residue) or with a string (an atom name inside a residue). Thus, to the coordinates of the first two $\alpha$ carbons:
In [9]:
CA_0 = prot[0]['CA']
CA_1 = prot[1]['CA']
Passing 'CA' has the same effect as passing 1, because the 'CA' atom is also the atom 1 ('N' is 0)
In [10]:
CA_0 == prot[0][1]
Out[10]:
Aditionally we can use the strings "bb" and "sc", for the backbone and sidechains coordinates respectively.
In [11]:
prot[0]['bb']
Out[11]:
We can use the coordinates we extract to compute, for example, distances:
In [12]:
bmb.utils.dist(CA_0, CA_1)
Out[12]:
We can ask for all the torsional angles in our molecule, at the moment this only works for $\phi$ and $\psi$ torsional angles, sidechain is comming soon!
In [13]:
prot.get_torsionals(sidechain=False)
Out[13]:
and we can also change them!
In [14]:
prot.set_phi(1, -60)
prot.set_psi(1, -40)
prot.get_torsionals(sidechain=False)
Out[14]:
and we can plot a minimalistic (one-point) ramachadran!
In [15]:
bmb.plot_ramachandran(prot);
We can also initialize a protein object from a pdb file at the moment
In [16]:
pdb = bmb.Protein(pdb='../pdbs/1ubq.pdb')
pdb
Out[16]:
At the moment objects initilizased from pdbs have only a subset of features compared to those initialized from a sequence, this will be remeadiated in the near future. What features are not available? Basically computing the energy and setting torsional angles. All of the other features should work the same. For example retriving the coordinates or computing all $\phi, \psi$ torsionals.
In [17]:
pdb.coords
Out[17]:
In [18]:
bmb.plot_ramachandran(pdb)
Out[18]:
Bomeba allows us to retrive some of the information in the pdb files, like the b_factors. You can also compute the b_factors from a protein object initialized from a sequence but they will just a bunch os zeros.
In [19]:
plt.plot(pdb.bfactors);
As with distances we can also compute angles and torsionals angles from arbitrary points. For example
In [20]:
a = pdb[0]['N']
b = pdb[10]['CA']
c = pdb[20]['O']
d = pdb[30]['N']
bmb.utils.get_torsional(a, b, c, d) * bmb.constants.radians_to_degrees
Out[20]:
Notice that angles are returned in radians and we multiple the result by the constant radians_to_degrees to get the result in the more familiar degrees.
Bomeba uses nglview to embed in the notebook visualization of molecules. This feature is at very early stage of development
In [21]:
pdb.view3d()