how to select a binding site

[Nick-Mul]

Hello, this is a nice package. I think I am getting just getting my head round it, but fairly new to Julia. I am interested in selecting a binding site around the ligand.... at the moment my code is super ugly, I am selecting for all the atoms near the ligand, then selecting the residues, with lots of lists.

for i in ligand
lig, protein_atom, dist = (closest(i,protein))
push!(binding_site, (protein[protein_atom]))
end

I feel there should be a more elegant way of doing this and was wondering you had some pointers.

Many thanks, Nick

[lmiq]

I don´t think your approach is bad, with the tools of PDBTools. But maybe I would do something like:

julia> using PDBTools

julia> pdb = wget("1BSX", "chain A");

julia> protein = select(pdb, "protein");

julia> ligand = select(pdb, "resname T3");

julia> active_site = Residue[]
       for r in eachresidue(protein)
           if closest(r, ligand)[3] < 3.0
               push!(active_site, r)
           end
       end

julia> active_site
   Array{Residue,1} with 2 residues.

julia> active_site_atoms = [ at for r in active_site for at in r ]
   Array{Atoms,1} with 18 atoms with fields:
   index name resname chain   resnum  residue        x        y        z occup  beta model segname index_pdb
     866    N     ASN     A      331      113   19.724   33.389   26.541  1.00 69.68     1       -       866
     867   CA     ASN     A      331      113   18.925   32.240   26.919  1.00 69.68     1       -       867
     868    C     ASN     A      331      113   19.080   31.053   25.978  1.00 69.68     1       -       868
                                                       ⋮ 
    1685  CD2     HIS     A      435      217   21.566   44.624   33.553  1.00 44.44     1       -      1685
    1686  CE1     HIS     A      435      217   19.843   45.038   34.841  1.00 44.44     1       -      1686
    1687  NE2     HIS     A      435      217   20.213   44.449   33.723  1.00 44.44     1       -      1687

or if you want the list of atoms directly:

julia> active_site = Atom[]
       for r in eachresidue(protein)
           if closest(r, ligand)[3] < 3.0
               append!(active_site, at for at in r)
           end
       end

[Nick-Mul]

Oh that's much more elegant! Thanks you Leandro.

[lmiq]

A small note: you can use the distance function instead of the closest, because you don't need the indices of the atoms involved. So the same above is:

julia> active_site = Atom[]
       for r in eachresidue(protein)
           if distance(r, ligand) < 3.0
               append!(active_site, at for at in r)
           end
       end

julia> active_site
   Array{Atoms,1} with 18 atoms with fields:

[lmiq]

I've added this as an example in the docs (version 1.2.0 to be released at any moment). It is the only example so far, but a list of useful examples like this is probably a good idea to have. Thanks for the feedback.