Hbonds and assigning the segment of each atom

[Jason Held]

Hi,

I have a multiple protein universe, and each protein is defined as a different segment. I have worked through the excellent Hbond tutorial to generate a dataframe with the acceptor/donor atoms/residues, etc. But, for some reason assigning segments isn't working as I'd expect. Code below.

# build example universe, and add a second segment for testing multiple segments called CORE

import MDAnalysis
from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import HydrogenBondAnalysis as HBA
from MDAnalysis.tests.datafiles import CRD, PSF, DCD, DCD2

u = MDAnalysis.Universe(PSF, DCD)

# add an empty CORE segment
core_segment = u.add_Segment(segid='CORE')
core_segment.atoms

#assign atoms to the CORE segment
core_atoms = u.select_atoms('resid 1:29 or resid 60:121 or resid 160-214')
core_atoms.residues.segments = core_segment
core_segment.atoms

print(u.segments.segids)

['4AKE' 'CORE'] # there are two segments now for testing, as expected

# Analyze the hbonds

hbonds = HBA(universe=u)
hbonds.run()

# build a dataframe based on the tutorial

FRAME = 0
DONOR = 1
HYDROGEN = 2
ACCEPTOR = 3
DISTANCE = 4
ANGLE = 5

df = pd.DataFrame(hbonds.results.hbonds[:, :DISTANCE].astype(int),
                  columns=["frame", "donor_ix", "hydrogen_ix", "acceptor_ix",])

df["distances"] = hbonds.results.hbonds[:, DISTANCE]
df["angles"] = hbonds.results.hbonds[:, ANGLE]
df["donorResname"] = u.atoms[df.donor_ix].resnames
df["acceptorResname"] = u.atoms[df.acceptor_ix].resnames
df["donorResid"] = u.atoms[df.donor_ix].resids
df["acceptorResid"] = u.atoms[df.acceptor_ix].resids
df["donorName"] = u.atoms[df.donor_ix].names
df["acceptorName"] = u.atoms[df.acceptor_ix].names

# up till here works fine based on the tutorial, but the code below focused on defining the # segment the atoms is in is where the problem arises. Both of these assignments below # lead to '<Segment 4AKE>, <Segment CORE>' in the dataframe, not 1 segment or the

# other. In short, I can't figure out how to assign only the segment that the donor or acceptor atom # is in.
df['donorSegment'] = u.atoms[df.donor_ix].segments
df['acceptorSegment'] = u.atoms[df.acceptor_ix].segments

#And its clear that I'm using segments incorrectly since .names or .resnames or .resids

# is a list for each and all Hbonds

print(len(u.atoms[df.acceptor_ix].names))

10858

# but I only get 2 segments below when I examine the .segments, though

# I'd expect 10858 like .names or .resnames or .resids results in

print(len(u.atoms[df.acceptor_ix].residues.segments))

2

Any help would be appreciated.

Cheers,

Jason

[Jason Held]

For anyone that comes across this issue in the future, I figured out a solution. The issue is that there is an incompatibility between atoms and atom groups, and while donor_ix is an atom, you can't call an atom out of a trajectory using u.atoms[df.donor_ix].segment. So, you have to separately determine the segment for each row of the dataframe based on the donor_ix or acceptor_ix with a for loop once you buid the dataframe. Full code example below.

# Above is all from tutorial. New stuff beyond the tutorial to assign segments for each HBond below.

result = []
temp =[]

for value in df['donor_ix']:
    result = str(u.atoms[value].segid)
    temp.append(result)
df['donorSegment'] = temp

# do it again for acceptors
result = []
temp =[]

for value in df['acceptor_ix']:
    result = str(u.atoms[value].segid)
    temp.append(result)
df['acceptorSegment'] = temp

[Hugo Macdermott-Opeskin]

Glad you found a solution and sorry we weren't more of a help.

Cheers 

Hugo