Import all the Python modules that we are going to use throughout the script

In [1]:
from cresset import flare
from rdkit.Geometry.rdGeometry import Point3D
import urllib.request

Download a PDB structure with ID code from the Protein Data Bank into the Flare project project

In [2]:
def download_from_pdb(project, code):
    url = 'http://files.rcsb.org/view/{0:s}.pdb'.format(code)
    response = urllib.request.urlopen(url)
    text = response.read().decode('utf-8')
    project.proteins.extend(flare.read_string(text, 'pdb'))
    return project.proteins[-1]

Create a new project
If we are running in Flare GUI, swap the Project in main_window with the one we have just created

In [3]:
p = flare.Project()
if (flare.main_window()):
    flare.main_window().project = p

Download 1udt and 1udu from the Protein Data Bank

In [4]:
protein_dict = {'1udt': download_from_pdb(p, '1udt'),
               '1udu': download_from_pdb(p, '1udu')}

Delete chain B

In [5]:
for seq_list in [prot.sequences for prot in p.proteins]:
    for i, seq in reversed(tuple(enumerate(seq_list))):
        if seq.chain == 'B':
            del seq_list[i]

Carry out protein preparation on all proteins and wait for it to be finished

In [6]:
prep = flare.ProteinPrep()
prep.proteins = p.proteins
prep.start()
prep.wait()

Find natural amino acid sequences

In [7]:
nat_seq_dict = {}
for seq_list in [prot.sequences for prot in p.proteins]:
    for seq in seq_list:
        if seq.type == flare.Sequence.Type.Natural:
            nat_seq_dict[seq.molecule.title.lower()] = seq

Create a dictionary of Calpha atom lists

In [8]:
atoms = {'1udt': [], '1udu': []}
for i in range(len(nat_seq_dict['1udt'].alignment)):
    if (nat_seq_dict['1udt'].alignment[i] == '-'
        or nat_seq_dict['1udu'].alignment[i] == '-'):
        continue
    for name in atoms.keys():
        atoms[name].extend([a for a in nat_seq_dict[
            name][i].atoms if a.name == 'CA'])

RMS fit Calphas of 1udu onto 1udt and apply the transformations to 1udu

In [9]:
trans, rms = flare.fit_atoms(atoms['1udu'], atoms['1udt'])
In [10]:
protein_dict['1udu'].apply_transform(trans)

Extract all sequences which are of Ligand type and move them to the Ligand table

In [11]:
for seq_list in [prot.sequences for prot in p.proteins]:
    for i, seq in reversed(tuple(enumerate(seq_list))):
        if (seq.type == flare.Sequence.Type.Ligand):
            p.ligands.extend(seq)
            del seq_list[i]

Zoom and center on ligands

In [12]:
if (flare.main_window()):
    flare.main_window().selected_ligands = p.ligands
    flare.main_window().selected_proteins = p.proteins
    flare.main_window().camera.center_on(p.ligands)
    flare.main_window().camera.fit_to_window(p.ligands)
In [13]:
def residues_near(what, nat_seq_dict, radius):
    residues = set()
    for al in what.atoms:
        if (al.atomic_number == 1):
            continue
        for ap in nat_seq_dict[what.protein.title.lower()].atoms:
            if (ap.atomic_number == 1):
                continue
            d = Point3D(*al.pos).Distance(Point3D(*ap.pos))
            if (d < radius):
                residues.add(str(ap.residue))
    return residues
In [14]:
contacts = [residues_near(l, nat_seq_dict, 3.5) for l in p.ligands]

Common contacts

In [15]:
common_contacts = sorted(list(contacts[0].intersection(contacts[1])),
       key=lambda i: int(i.split()[-1]))

Color common contacts in orange, and display them as thin sticks

In [16]:
for seq in nat_seq_dict.values():
    for r in seq:
        if (str(r) in common_contacts):
            for a in r.atoms:
                if (a.atomic_number == 6):
                    a.style.color = (1, 0.5, 0.1)
                a.style.display_style = flare.DisplayStyle.CPK

Viagra-only contacts

In [17]:
viagra_contacts = sorted(list(contacts[0].difference(contacts[1])),
       key=lambda i: int(i.split()[-1]))

Color Viagra-only contacts with the same color as the Viagra ligand, and display them as CPK

In [18]:
for r in nat_seq_dict['1udt']:
    if (str(r) in viagra_contacts):
        for a in r.atoms:
            if (a.atomic_number == 6):
                a.style.color = p.ligands[0].style.color
            a.style.display_style = flare.DisplayStyle.CPK

Cialis-only contacts

In [19]:
cialis_contacts = sorted(list(contacts[1].difference(contacts[0])),
       key=lambda i: int(i.split()[-1]))

Color Cialis-only contacts with the same color as the Cialis ligand, and display them as CPK

In [20]:
for r in nat_seq_dict['1udu']:
    if (str(r) in cialis_contacts):
        for a in r.atoms:
            if (a.atomic_number == 6):
                a.style.color = p.ligands[1].style.color
            a.style.display_style = flare.DisplayStyle.CPK

View each ligand-protien complex in a grid

In [21]:
if (flare.main_window()):
    flare.main_window().camera.grid_enabled=True

Make the view less cluttered by clipping the front and back of the structures

In [22]:
if (flare.main_window()):
    flare.main_window().camera.front_clip=-9
    flare.main_window().camera.back_clip=9
In [23]:
if (flare.main_window()):
    flare.main_window().camera.fit_to_window(flare.atom_pick.active_site(p.ligands,p.proteins,5.0))

Print out the residues

In [24]:
print(f"Viagra Only =",viagra_contacts)
print(f"Cialis Only =",cialis_contacts)
print(f"Common =",common_contacts)
Viagra Only = ['A TYR 664', 'A ALA 779']
Cialis Only = ['A SER 663', 'A ALA 783', 'A PHE 786']
Common = ['A TYR 612', 'A LEU 804', 'A MET 816', 'A GLN 817', 'A PHE 820']