CaviTracer: Pores reconstruction

[karolamik13]

@briza81 @MatthewLicht
I prepared the reconstruction of pores from identified channels via calcChannels().

Additions:

• return_details is added to calcChannels() because additional information is required to reconstruct pores.
• calcPoresFromChannels
• showPores()

As an example, I used the 2FGQ used by Jan. In the structure, we have a very clear and nice pore (see the comparison below between Jan's channel prediction and the further reconstruction of two pores).

The code should be adapted to multi-model PDBs and trajectories
(I will do that in the near future; probably adding new commits here).

PORE 1:

PORE 2:

How to use it:
from prody import *
p = parsePDB('2fgq')
protein = p.select('protein')
channels, surface, details = calcChannels(protein, r1=20, r2=1, start_point = [39.277, 43.995, -0.961],
max_deviation=0.1, return_details=True)

vmd_path = '/usr/local/bin/vmd'
model = getVmdModel(vmd_path, protein)

pores = calcPoresFromChannels(channels, details)
showPores(pores[0], model=model)
showPores(pores[1], model=model)

[jamesmkrieger]

Looks good to me. The only thing is it seems redundant to have an extra show function that doesn’t seem to do anything except call another show function with all the parameters the same. You could just alias it instead

[karolamik13]

Channels and pores use the same function to obtain residues that form those channels/pores, but there are some differences, such as LOGGER info/names and output files having different names. Also, the docs would be confusing if for pores we had channel information. That is why I didn't use an alias and instead created a supporting function that is used by both channels and pores.

[karolamik13]

After those improvements, those functions are available:

getPoreResidueNames(protein, pores)
getPoreParameters(pores)

[jamesmkrieger]

Channels and pores use the same function to obtain residues that form those channels/pores, but there are some differences, such as LOGGER info/names and output files having different names. Also, the docs would be confusing if for pores we had channel information. That is why I didn't use an alias and instead created a supporting function that is used by both channels and pores.

yes, that makes sense. It's just the show function that was the same

[karolamik13]

After testing the GPCR protein and obtaining ~400 pores, which were crashing my Open3D window, I decided that the filters inside calcPoresFromChannels() are needed. I added several filters by end-to-end distance of openings of pore, length, volume, and bottleneck. With those filters, I was able to exclude useless pores, such as those with openings close together. From ~400 (without filters), I had ~70 for pores2 and 17 pores for pores1 (criteria below). Protein is taken from caseStudy2.pdb from the WatFinder tutorial (r2=0.8 for channels prediction).

Example of usage:
pores1 = calcPoresFromChannels(channels, details, min_end_to_end=45, min_bottleneck=0.6)

pores2 = calcPoresFromChannels(channels, details, min_end_to_end=15, min_volume=600, min_length=60)

The orientation is probably different, but this is what we showed in the WatFinder paper (water interaction)

Single pores:
PORE 1:

PORE 2:

[karolamik13]

I added functions that can handle multi-model PDB and trajectory for pores:

• calcPoresFromChannelsMultipleFrames()
• getPoreParametersMultipleFrames()
• getPoreResidueNamesMultipleFrames()

I didn't want to multiply the code, so I wrote:

• getObjectResidueNames()
• getObjectResidueNamesMultipleFrames()
  which are used by both channels and pores, but we need to define which object type we want to analyze (object_type='channel' or 'pore').

Additionally, I modified calcPoresFromChannels() and calcPoresFromChannelsMultipleFrames() to add output_path to save files and PQR/PDB files.

I tested both channels, pore, multi-model PDB, and trajectories. Everything looks good.

[karolamik13]

The only thing still missing here is parallel computing for multi-model PDBs and trajectories. I will take a look at that later, and hopefully that will be the only thing I will add/change.

[karolamik13]

Now, we have multiprocessing included for:

• calcPoresFromChannelsMultipleFrames()
• calcChannelsMultipleFrames()
• calcChannelSurfaceOverlaps() - that was done before, but this time with protection when someone is using Windows or macOS (mp_context can be used to use 'spawn' instead of "fork"), and handling empty files when channels are not detected.

Tests were done with the files below (results are also included):

Multi-model PDB:
multimodelPDB_script.py
multimodelPDB_test.txt

Trajectory (WatFinder case study 2; first 10 frames or so):
dcd_test_script.py
dcd_test.txt

I am done with code improvements. @jamesmkrieger, let me know if you have any suggestions.