HSExposure list of residues?

Question

[johnbmcarthur]

from Bio import PDB

pdbfile=open('2wqq.pdb')
p=PDB.PDBParser()
s=p.get_structure('X', pdbfile)
m=s[0]
RADIUS=10.0

hse=PDB.HSExposureCA(m, RADIUS) 
residue_list=PDB.Selection.unfold_entities(m,'R')
for r in residue_list[:]:
    print(r.get_resname(), r.xtra)

This returns a dictionary for each residue where the keys are 'EXP_HSE_A_U', 'EXP_HSE_A_D', and 'EXP_CB_PCB_ANGLE'. I would instead like a dictionary where the key is the residue id number of r and the values are the residue id numbers of residues within the 'up' half shell. Does anybody know how I could do this?

[johnbmcarthur]

Update:

Got it almost sort of working, but the output doesn't seem to be correct. The following code is supposed to make a dictionary of residues that are nearby (within defined radius), the angle between res1 Ca-Cb and res1 Ca-res2 Ca is less than 90 degrees, and the angle between res 2 Ca-Cb and res2 Ca - res1 Ca is less than 90 degrees. I'm sorry if it's a mess. I don't know what I am doing. I think there's something wrong with the angle calculation...

from Bio.PDB import *
import warnings
import math
from Bio import PDB
import numpy 
from numpy.linalg import norm
from collections import defaultdict
from math import pi

residues=[]
pdbfile=open('2wqq.pdb')
p=PDB.PDBParser()
s=p.get_structure('X', pdbfile)
m=s[0]
RADIUS=10.0
frontcontacts= defaultdict(list)
hse=PDB.HSExposureCA(m, RADIUS)
residue_list=PDB.Selection.unfold_entities(m,'R')
for model in s.get_list():
    for chain in model.get_list():
         for residue in chain.get_list():
            residues.append(residue)
for c1 in residue_list[:]:
    if not is_aa(c1) or not c1.has_id('CA'):
        continue
    one = c1['CA'].get_coord()
    for c2 in residues:
        if not is_aa(c2) or not c2.has_id('CA') or c2==c1:
            continue
        two = c2['CA'].get_coord()
        if numpy.linalg.norm(one-two) < RADIUS:
            angle1 = calc_angle(c2['CA'].get_vector(), c1['CA'].get_vector(), hse._get_gly_cb_vector(c1))
            angle2 = calc_angle(c1['CA'].get_vector(), c2['CA'].get_vector(), hse._get_gly_cb_vector(c2))
            if angle1*180/pi < 90 and angle2*180/pi < 90:
                frontcontacts[str(c1.id[1])].append(str(c2.id[1]))
                #print(c1.id[1],c2.id[1],angle1*180/pi)
print(frontcontacts)

Answer 1

[Mensur Dlakic]

Pretty sure you can't do this without modifying the original code - and even that may be questionable. I have an old version of this program called hsexpo before it was part of BioPython, but it was >15 years ago when the original paper was published so I don't remember how I got it. You may want to contact Thomas Hamelryck directly ( https://www1.bio.ku.dk/english/staff/?pure=en/persons/271052 ). Beware that the original program also doesn't have the exact functionality you seek. Its output is below, and it is residues and their numbers followed by HSEb-up and HSEb-down values.

MET 1  14  15
LYS 2  12  18
PHE 3  24  22
VAL 4  19  26
SER 5  24  26
PHE 6  25  25
ASN 7  22  24
ILE 8  23  21
ASN 9  10  23
GLY 10   0   0