[Question about native contacts analysis] Can I only compute heavy atom pairs that their residue number is not same?

[narubi2]

Hi,

I'm working with MDAnalysis and have a question about 'native contact analysis'
I'm trying to get Q value (native contact value) of all heavy atoms.

I tried as follows:

# seletion
sel_ca_atom1 = "protein and not (name H*)"
sel_ca_atom2 = "protein and not (name H*)"
ref_ca_atom1 = u.select_atoms(sel_ca_atom1)
ref_ca_atom2 = u.select_atoms(sel_ca_atom2)

# calculate
ca = contacts.Contacts(u, select=(sel_ca_atom1, sel_ca_atom2),
                        refgroup=(ref_ca_atom1, ref_ca_atom2),
                        radius=4.5,
                        method='soft_cut',
                        kwargs={'beta': 5.0,
                                'lambda_constant': 1.5}).run()

It worked, but it took a long time to compute. Additionally it was hard to interpret the result.
I read the paper about Q value (https://doi.org/10.1073/pnas.1311599110) again, and found they compute only pairs of heavy atoms i and j belonging to residues θi and θj are in contact provided that |θi - θj| > 3

Is there a way to only count pairs of heavy atoms whose their residue number difference is larger than 3?

Any suggestion or advice would be appreciated.
Thanks

[orbeckst]

The current code for Contacts does not allow you to exclude neighbors and nearest neighbors. You could probably write your own version that does something like it. When you look at the source code's Contact._single_frame function

mdanalysis/package/MDAnalysis/analysis/contacts.py

Lines 486 to 499 in 0582265

	def _single_frame(self):
	self.results.timeseries[self._frame_index][0] = self._ts.frame
	# compute distance array for a frame
	d = distance_array(self.grA.positions, self.grB.positions,
	box=self._get_box(self._ts))
	for i, (initial_contacts, r0) in enumerate(zip(self.initial_contacts,
	self.r0), 1):
	# select only the contacts that were formed in the reference state
	r = d[initial_contacts]
	r0 = r0[initial_contacts]
	q = self.fraction_contacts(r, r0, **self.fraction_kwargs)
	self.results.timeseries[self._frame_index][i] = q

then you see that we calculate all distances with

     d = distance_array(self.grA.positions, self.grB.positions, 
                         box=self._get_box(self._ts)) 

but then filter these distances with initial_contacts.

initial_contacts is a boolean array of the same shape as d. So if you modify initial_contacts to exclude the contacts that you don't want to see then the Contacts code will produce what you want, e.g., in the special case that your selections A and B are a linear sequence of residues then you would set the diagonal and the two off-diagnals to False and they should be excluded from the calculations, perhaps something like this

import numpy as np
c = Contacts(...)
initial_contacts = c.initial_contacts[0]  # depending on how you run Contacts, can contain more
mask = np.ones(initial_contacts.shape, dtype=bool)
# keep second, first, and main diagonal at True, set upper and lower triangles to False
mask = np.triu(mask, k=-2) * np.tril(mask, k=2)
# set contacts on 2nd, 1st, and main diagonal to False
initial_contacts[mask] = False
# and update the class before running analysis
c.initial_contacts[0] = initial_contacts

# run analysis
c.run()

The above code will not just work — I only had a quick look at the source code and I think you can hack it using the ideas above. Perhaps the notes on native contact analysis are also useful.

You could also rewrite the Contacts code itself.

In any case, you'll have some coding to do in order to achieve what you want.

[narubi2]

Hi orbeckst.

Thank you for reply
I'll try your suggestions.