PyLipID

https://github.com/wlsong/pylipid
Admin
PyLipID is a python package for analyzing lipid interactions with membrane proteins from Molecular...
Tokens:17,413
Snippets:146
Trust Score:6.4
License:MIT
Update:4 months ago
Show doc for...
Tokens:
Raw
### Install PyLipID from source

Source: https://github.com/wlsong/pylipid/blob/master/docs/index.md

Install the PyLipID package by cloning the GitHub repository and running the setup script. This method is useful for developers or for obtaining the latest unreleased version.

```bash
git clone git://github.com/wlsong/PyLipID.git
python setup.py install

```

--------------------------------

### Install PyLipID from Source

Source: https://github.com/wlsong/pylipid/blob/master/README.rst

Clone the PyLipID repository from GitHub and install it from the source code. This method is an alternative to pip installation.

```bash
git clone https://github.com/wlsong/PyLipID
```

```bash
python setup.py install
```

--------------------------------

### Install PyLipID using pip

Source: https://github.com/wlsong/pylipid/blob/master/README.rst

Install the PyLipID package using pip. This is the recommended installation method.

```bash
pip install pylipid
```

--------------------------------

### Install PyLipID using pip

Source: https://github.com/wlsong/pylipid/blob/master/docs/INSTALL.md

Installs the PyLipID package and its dependencies using the pip package installer. Ensure all listed dependencies are available in their required versions for a successful installation.

```bash
$ pip install pylipid
```

--------------------------------

### PyLipID Lipid Interaction Analysis Setup

Source: https://github.com/wlsong/pylipid/blob/master/docs/demo.md

This script configures parameters for analyzing lipid interactions. It requires specifying trajectory and topology files, lipid type, interaction cutoffs, and output settings. Most parameters can be left at their default values after initial setup.

```python
import numpy as np
import matplotlib.pyplot as plt
from pylipid.api import LipidInteraction
from pylipid.util import check_dir

##################################################################
##### This part needs changes according to your setting ##########
##################################################################
trajfile_list = ["run1/md.xtc", "run2/md.xtc"]
topfile_list = ["run1/md.gro", "run2/md.gro"]  # topology file is needed when trajectory format does not
                                               # provide topology information. See mdtraj.load() for more
                                               # information.
dt_traj = None  # the timestep of trajectories. Need to use this param when trajectories are in a format
                # with no timestep information. Not necessary for trajectory formats of e.g. xtc, trr.
stride = 1   # tell pylipid to analyze every stride-th frame. Can be used to save computation memory
             # and speed up the calculation.

lipid = "CHOL"   # residue name in the topology.
lipid_atoms = None  # all lipid atoms will be considered for interaction calculation.
cutoffs = [0.5, 0.8]  # dual-cutoff scheme for coarse-grained simulations. Single-cutoff scheme can be
                      # achieved by using the same value for two cutoffs.

nprot = 1   # if the simulation system has N copies of receptors, "nprot=N" will report interactions
            # averaged from the N copies, but "nprot=1" will ask PyLipID to report interaction for
            # each copy.

binding_site_size = 4  # binding site should contain at least four residues.

n_top_poses = 3     # write out num. of representative bound poses for each binding site.
n_clusters = "auto"  # cluster the bound poses for a binding site into num. of clusters. PyLipID
                     # will write out a pose conformation for each of the cluster. By default, i.e. "auto", PyLipID
                     # will use a density based clusterer to find possible clusters.

save_dir = None  # save at current working directory if it is None.
save_pose_format = "gro"  # format that poses are written in
save_pose_traj = True  # save all the bound poses in a trajectory for each binding site. The generated
                       # trajectories can take some disk space (up to a couple GB depending on your system).
save_pose_traj_format = "xtc"  # The format for the saved pose trajectories. Can take any format that is supported
                               # by mdtraj.

timeunit = "us"  # micro-sec. "ns" is nanosecond. Time unit used for reporting the results.
resi_offset = 0  # shift the residue index, useful for MARTINI models.

radii = None  # Radii of protein atoms/beads. In the format of python dictionary {atom_name: radius}
              # Used for calculation of binding site surface area. The van der waals radii of common atoms were
              # defined by mdtraj (https://github.com/mdtraj/mdtraj/blob/master/mdtraj/geometry/sasa.py#L56).
              # The radii of MARTINI 2.2 beads were included in PyLipID.

pdb_file_to_map = None   # if a pdb coordinate of the receptor is provided, a python script
                         # "show_binding_site_info.py" will be generated which maps the binding
                         # site information to the structure in PyMol. As PyMol cannot recognize
                         # coarse-grained structures, an atomistic structure of the receptor is needed.

fig_format = "pdf"  # format for all pylipid produced figures. Allow for formats that are supported by
                    # matplotlib.pyplot.savefig().

num_cpus = None  # the number of cpu to use when functions are using multiprocessing. By default, i.e. None, the functions will use up all the cpus available. This can use up all the memory in some cases.

#####################################
###### no changes needed below ######
#####################################

```

--------------------------------

### Install PyLipID using pip

Source: https://github.com/wlsong/pylipid/blob/master/docs/index.md

Install the PyLipID package using pip, the Python package installer. This is the recommended method for most users.

```bash
pip install pylipid

```

--------------------------------

### Import Libraries and Print PyLipID Version

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/1-Distance_cutoff_determination.ipynb

Imports necessary libraries for PyLipID analysis and prints the installed version. Ensure you are using version 1.4 or later for compatibility.

```python
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import pylipid
from pylipid.api import LipidInteraction
print(pylipid.__version__)  # the following script works for versions later than 1.4
```

--------------------------------

### Get Binding Site Information

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Use the `binding_site` function to retrieve information about a specific binding site. This function returns a pandas DataFrame containing residue details for the specified binding site.

```python
df = li.binding_site(binding_site_id=3) # check the binding site No 3
df.head(10)
```

--------------------------------

### Initialize LipidInteraction Class

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Sets up the main LipidInteraction class with trajectory and topology files, lipid type, interaction cutoffs, and other simulation parameters. Data will be saved in the current directory by default.

```python
trajfile_list = ["./traj_data/A2a/run1/protein_lipids.xtc", "./traj_data/A2a/run2/protein_lipids.xtc"]
topfile_list = ["./traj_data/A2a/run1/protein_lipids.gro", "./traj_data/A2a/run2/protein_lipids.gro"]
lipid = "CHOL"
cutoffs = [0.475, 0.8] # use of dual-cutoff
nprot = 1  # num. of protein copies in the system. if the simulation system has N copies of receptors, 
           # "nprot=N" will report interactions averaged from the N copies, but "nprot=1"
           # will ask pylipid to report interaction data for each copy.
timeunit = 'us' # micro-second
save_dir = None  # if None, pylipid data will be saved at current working directory. 

# initialize 
li = LipidInteraction(trajfile_list, topfile_list=topfile_list, cutoffs=cutoffs, lipid=lipid, 
                      nprot=nprot, save_dir=save_dir)

```

--------------------------------

### Run GABAA Simulations and Save Test Cutoff Data

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/1-Distance_cutoff_determination.ipynb

This snippet sets up parameters for GABAA simulations, runs the `test_cutoffs` function to generate interaction data, and saves the results to a pickle file for later analysis. Ensure `check_dir`, `test_cutoffs`, `pickle`, `lipid_atoms`, `nprot`, `stride`, and `timeunit` are defined.

```python
save_dir = check_dir(suffix="test_cutoffs_GABAA", print_info=False)
lipid = "CHOL"
trajfile_list = ["./traj_data/GABAA/run1/protein_lipids.xtc", 
                 "./traj_data/GABAA/run2/protein_lipids.xtc"]
topfile_list = ["./traj_data/GABAA/run1/protein_lipids.gro", 
                 "./traj_data/GABAA/run2/protein_lipids.gro"]

num_of_binding_sites, duration_avgs, num_of_contacting_residues = test_cutoffs(
                                 cutoff_list, trajfile_list, topfile_list, lipid, lipid_atoms, 
                                 nprot=nprot, stride=stride, save_dir=save_dir, timeunit=timeunit)

test_data = {"num_of_binding_sites": num_of_binding_sites,
             "duration_avgs": duration_avgs,
             "num_of_contacting_residues": num_of_contacting_residues,
             "test_cutoff_list": cutoff_list}

with open(f"{save_dir}/test_cutoff_data_{lipid}.pickle", "wb") as f:
    pickle.dump(test_data, f, 2)
```

--------------------------------

### Detect and Analyze Binding Sites

Source: https://github.com/wlsong/pylipid/blob/master/docs/demo.md

Computes potential binding sites based on a threshold and analyzes their properties if any are detected. This includes site duration, occupancy, lipid count, and koff. Analysis is skipped if no binding sites are found.

```python
li.compute_binding_nodes(threshold=binding_site_size, print_data=False)
if len(li.node_list) == 0:
    print("*"*50)
    print("No binding site detected! Skip analysis for binding sites.")
    print("*"*50)
else:
    li.compute_site_duration(binding_site_id=None)
    li.compute_site_occupancy(binding_site_id=None)
    li.compute_site_lipidcount(binding_site_id=None)
    li.compute_site_koff(binding_site_id=None, plot_data=True, fig_close=True,
                         fig_format=fig_format, num_cpus=num_cpus)
```

--------------------------------

### Initialize LipidInteraction and Compute Residue Contacts

Source: https://github.com/wlsong/pylipid/blob/master/docs/demo.md

Initializes the LipidInteraction object with trajectory and topology files, then collects residue contact data. Ensure all necessary parameters like file paths, cutoffs, and lipid information are correctly provided.

```python
li = LipidInteraction(trajfile_list, topfile_list=topfile_list, cutoffs=cutoffs, lipid=lipid,
                      lipid_atoms=lipid_atoms, nprot=1, resi_offset=resi_offset,
                      timeunit=timeunit, save_dir=save_dir, stride=stride, dt_traj=dt_traj)
li.collect_residue_contacts()
```

--------------------------------

### Analyze Bound Poses and Save Trajectories

Source: https://github.com/wlsong/pylipid/blob/master/docs/demo.md

Analyzes the trajectories of bound poses within detected binding sites and saves them to disk if configured. This step requires binding sites to be detected.

```python
pose_traj, pose_rmsd_data = li.analyze_bound_poses(binding_site_id=None, pose_format=save_pose_format,
                                                       n_top_poses=n_top_poses, n_clusters=n_clusters,
                                                       fig_format=fig_format, num_cpus=num_cpus)
# save pose trajectories
if save_pose_traj:
    for bs_id in pose_traj.keys():
        pose_traj[bs_id].save("{}/Bound_Poses_{}/Pose_traj_BSid{}.{}".format(li.save_dir, li.lipid, bs_id,
                                                                          save_pose_traj_format))
# save memory space
del pose_traj
```

--------------------------------

### Display DataFrame Head with New Columns

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Displays the first 10 rows of the Pylipid dataset DataFrame. This shows how the 'Duration', 'Occupancy', and 'Lipid Count' columns are added after calculations.

```python
li.dataset.head(10) # show 10 lines
```

--------------------------------

### Calculate Binding Sites with Varying Cutoffs

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/1-Distance_cutoff_determination.ipynb

This script iterates through a list of dual-cutoff pairs to calculate and compare binding sites and their residence times. It utilizes the LipidInteraction API from PyLipID. Set plot_data=False to optimize for speed and disk space.

```python
from pylipid.api import LipidInteraction

lipid = "CHOL"
trajfile_list = ["./traj_data/GABAA/run1/protein_lipids.xtc", 
                 "./traj_data/GABAA/run2/protein_lipids.xtc"]
topfile_list = ["./traj_data/GABAA/run1/protein_lipids.gro", 
                 "./traj_data/GABAA/run2/protein_lipids.gro"]

atomistic_structure = "./traj_data/GABAA/GABAA_atomistic.pdb" # use `save_pymol_script` to map 
                                                              # binding sites to this atomistic structure. 

cutoffs_list = [[0.45, 0.7], [0.475, 0.7], [0.5, 0.7]]

for cutoffs in cutoffs_list:
    save_dir = "GABAA_cutoffs_{}".format("_"..join([str(cutoff) for cutoff in cutoffs]))
    li = LipidInteraction(trajfile_list, topfile_list=topfile_list, lipid=lipid, cutoffs=cutoffs, 
                         save_dir=save_dir, timeunit="us")
    li.collect_residue_contacts()
    li.compute_residue_koff(plot_data=False)  # plot_data=False to save time and disc space
    li.compute_binding_nodes(print_data=False) 
    li.compute_site_koff(plot_data=True)
    li.write_site_info()
    li.save_pymol_script(atomistic_structure)
    
```

--------------------------------

### Save Binding Site Data

Source: https://github.com/wlsong/pylipid/blob/master/docs/demo.md

Saves 'Duration BS' and 'Occupancy BS' data for analysis. Ensure 'li.node_list' is not empty before execution.

```python
if len(li.node_list) > 0:
    for item in ["Duration BS", "Occupancy BS"]:
        li.save_data(item=item)
```

--------------------------------

### Display Dataset Head

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Displays the first 10 rows of the dataset, which has been updated with koff calculation results. This includes columns like 'Residence Time', 'Koff', 'R Squared', and their bootstrap averages.

```python
li.dataset.head(10)
```

--------------------------------

### Displaying Binding Site Data

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

This snippet shows a sample of binding site data, indicating residue numbers and associated values. It's useful for understanding the structure of the output data.

```text
20                 0.998424                              0.998174  
300                0.998424                              0.998174  
296                0.998424                              0.998174  

[8 rows x 22 columns]
```

--------------------------------

### Save PyMOL Script for Binding Sites

Source: https://github.com/wlsong/pylipid/blob/master/docs/demo.md

Generates a PyMOL script to visualize binding sites, using a specified PDB file as a reference. This is conditional on `pdb_file_to_map` being provided.

```python
if pdb_file_to_map is not None:
    li.save_pymol_script(pdb_file_to_map)
```

--------------------------------

### Define Cutoff Ranges for Testing

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/1-Distance_cutoff_determination.ipynb

Sets up lists of lower and upper distance cutoff values to be used in the `test_cutoffs` function. The `product` function from `itertools` is then used to generate all possible pairs of these cutoffs for comprehensive testing.

```python
from itertools import product

lower_cutoff = [0.375, 0.40, 0.425, 0.45, 0.475, 0.5]
upper_cutoff = [0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85]
cutoff_list = list(product(lower_cutoff, upper_cutoff))
```

--------------------------------

### LipidInteraction Class Initialization

Source: https://github.com/wlsong/pylipid/blob/master/docs/api/index_LipidInteraction.md

Initializes the LipidInteraction class to set up trajectory analysis. This involves loading trajectory files, defining interaction cutoffs, specifying lipid types, and configuring various parameters for the analysis.

```APIDOC
## LipidInteraction Class

### Description
Initializes the LipidInteraction class to set up trajectory analysis. This involves loading trajectory files, defining interaction cutoffs, specifying lipid types, and configuring various parameters for the analysis.

### Parameters
* **trajfile_list** (str or list of str) – Trajectory filename(s). Read by mdtraj.load() to obtain trajectory information.
* **cutoffs** (list of two scalar or a scalar, default=[0.475, 0.7]) – Cutoff value(s) for defining contacts. When a list of two scalar are provided, the dual-cutoff scheme will be used. A contact in the dual-cutoff scheme starts when a lipid gets closer than the lower cutoff, and ends when this lipid moves farther than the upper cutoff. The duration between the two time points is the duration of this contact.
* **lipid** (str, default="CHOL") – Lipid name in topology.
* **topfile_list** (str or list of str, default=None) – Topology filename(s). Most trajectory formats do not contain topology information. Provide either the path to a RCSB PDB file, a trajectory, or a topology for each trajectory in trajfile_list for the topology information. See mdtraj.load(). for more information.
* **lipid_atoms** (list of str, default=None) – Lipid atom names. Only interactions of the provided atoms will be considered for the calculation of contacts. If None, all atoms of the lipid molecule will be used.
* **nprot** (int, default=1) – Number of protein copies in the system. If the system has N copies of the protein, ‘nprot=N’ will report averaged values from the N copies, but ‘nprot=1’ will report interaction values for each copy.
* **resi_offset** (int, default=0) – Shift residue index in the reported results from what is shown in the topology. Can be useful for MARTINI force field.
* **save_dir** (str, default=None) – The root directory to store the generated data. By default, a directory Interaction_{lipid} will be created in the current working directory, under which all the generated data are stored.
* **timeunit** ({"us", "ns"}, default="us") – The time unit used for reporting results. “us” is micro-second and “ns” is nanosecond.
* **stride** (int, default=1) – Only read every stride-th frame. The same stride in mdtraj.load().
* **dt_traj** (float, default=None) – Timestep of trajectories. It is required when trajectories do not have timestep information. Not needed for trajectory formats of e.g. xtc, trr etc. If None, timestep information will take from trajectories.
```

--------------------------------

### Saving Binding Site Information

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Writes the calculated binding site information to a text file. The order of data can be controlled using the 'sort_residue' parameter. Requires an initialized LipidInteraction object 'li'.

```python
li.write_site_info(sort_residue="Residence Time") # write binding site information in a txt file
```

--------------------------------

### Compare Simulation Data for A2aR and GABAA Receptors

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/1-Distance_cutoff_determination.ipynb

This function generates comparative graphs for simulation data (number of binding sites, average durations, and contacting residues) between A2aR and GABAA receptors. It requires `matplotlib.pyplot` and `numpy`. Load data from pickle files and ensure the cutoff lists match before plotting.

```python
import matplotlib.pyplot as plt

def graph_compare(cutoff_list, data1, data2, label1, label2, ylabel, title, fn):
    colors = ["#20639B", "#ED553B"]
    fig, ax = plt.subplots(1, 1, figsize=(len(cutoff_list)*0.42, 3.6))
    m1 = ax.scatter(np.arange(len(cutoff_list)), data1, s=50, color=colors[0])
    m2 = ax.scatter(np.arange(len(cutoff_list)), data2, s=50, color=colors[1])
    ax.set_xticks(np.arange(len(cutoff_list)))
    ax.set_xticklabels(cutoff_list, rotation=45, ha='right')
    ax.set_xlabel("Dual cutoff")
    ax.set_ylabel(ylabel)
    ax.set_title(title)
    ax.legend([m1, m2], [label1, label2], loc="upper left")
    plt.tight_layout()
    plt.savefig(fn, dpi=200)
    return

lipid = "CHOL"
with open(f"./test_cutoffs_A2a/test_cutoff_data_{lipid}.pickle", "rb") as f:
    data_A2a = pickle.load(f)
with open(f"./test_cutoffs_GABAA/test_cutoff_data_{lipid}.pickle", "rb") as f:
    data_GABAA = pickle.load(f)
    
if data_A2a["test_cutoff_list"] == data_GABAA["test_cutoff_list"]:
    cutoff_list = data_A2a["test_cutoff_list"]
    # plot 1 
    graph_compare(cutoff_list, 
                  [data_A2a["num_of_binding_sites"][cutoffs] for cutoffs in cutoff_list], 
                  [data_GABAA["num_of_binding_sites"][cutoffs] for cutoffs in cutoff_list],
                  "A2a", "GABAA", "Num. of binding sites",
                 lipid, f"cmp_num_of_bs_{lipid}.png")

    # plot 2
    graph_compare(cutoff_list, 
                  [data_A2a["duration_avgs"][cutoffs] for cutoffs in cutoff_list], 
                 [data_GABAA["duration_avgs"][cutoffs] for cutoffs in cutoff_list], 
                  "A2a", "GABAA", f"Average Durations ({timeunit})",
                 lipid, f"cmp_durations_{lipid}.png")

    # plot 3
    graph_compare(cutoff_list, 
                  [data_A2a["num_of_contacting_residues"][cutoffs] for cutoffs in cutoff_list],
                 [data_GABAA["num_of_contacting_residues"][cutoffs] for cutoffs in cutoff_list], 
                  "A2a", "GABAA", "num. of contacting residues",
                 lipid, f"cmp_num_of_contacting_residues_{lipid}.png")

```

--------------------------------

### Compute Binding Site Metrics

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Calculates the duration, occupancy, and lipid count for binding sites. This snippet demonstrates how to use the `compute_site_duration`, `compute_site_occupancy`, and `compute_site_lipidcount` methods.

```python
durations = li.compute_site_duration()
occupancies = li.compute_site_occupancy()
lipidcounts = li.compute_site_lipidcount()
print(len(node_list), len(li.trajfile_list))
print(len(durations), len(durations[0]))
print(len(occupancies), len(occupancies[0]))
print(len(lipidcounts), len(lipidcounts[0]))
```

--------------------------------

### Calculate and Save PIP2 Interaction Metrics

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/1-Distance_cutoff_determination.ipynb

Calculates interaction metrics for PIP2 with the A2a receptor using specified cutoffs and saves the results to a pickle file. This is useful for initial data generation before visualization.

```python
lipid = "POP2"
lipid_atoms = ["C1", "C2", "C3", "PO4", "P1", "P2"]
num_of_binding_sites, duration_avgs, num_of_contacting_residues = test_cutoffs(
                                 cutoff_list, trajfile_list, topfile_list, lipid, lipid_atoms, 
                                 nprot=nprot, stride=stride, save_dir=save_dir, timeunit=timeunit)
test_data = {"num_of_binding_sites": num_of_binding_sites,
             "duration_avgs": duration_avgs,
             "num_of_contacting_residues": num_of_contacting_residues,
             "test_cutoff_list": cutoff_list}

with open(f"{save_dir}/test_cutoff_data_{lipid}.pickle", "wb") as f:
    pickle.dump(test_data, f, 2)
```

--------------------------------

### Clone PyLipID from GitHub

Source: https://github.com/wlsong/pylipid/blob/master/docs/INSTALL.md

Clones the PyLipID repository from GitHub to build the package from source. This method is for users who wish to modify or inspect the source code.

```bash
$ git clone https://github.com/wlsong/PyLipID.git
```

--------------------------------

### Analyze Bound Poses

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Analyzes bound poses for binding sites, optionally specifying the pose format and the number of top poses or clusters. Set binding_site_id to None to analyze all residues.

```python
pose_pool, pose_rmsd = li.analyze_bound_poses(binding_site_id=None, pose_format="pdb",
                                             n_top_poses=3, n_clusters='auto') # None means analysis for all residues
```

--------------------------------

### Print Binding Site and Modularity Data

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Prints the calculated node list (binding site residue IDs) and the modularity score. Modularity indicates the quality of the detected binding sites.

```python
print(node_list)
print()
print(modularity)
```

--------------------------------

### Display Head of Surface Area DataFrame

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Displays the first few rows of the `surface_area` DataFrame, showing the time evolution of surface area for each binding site in nm^2.

```python
surface_area.head()
```

--------------------------------

### Compute Binding Site Koff and Residence Time

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Calculates the koff and residence time for each binding site. Set `plot_data=True` to generate a koff figure for each site.

```python
koff_BS, res_time_BS = li.compute_site_koff(plot_data=True)
```

--------------------------------

### Binding Site Analysis Results

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Displays the detailed results of binding site analysis, including residence time, koff, duration, occupancy, and lipid count for various residues. This table summarizes key interaction metrics.

```python
Result:
    Residue  Residue ID  Duration  Duration std  Occupancy  Occupancy std  \
282  283PRO         282  0.032000      0.024000   0.996016       0.597610   
293  294ARG         293  0.025000      0.021794   0.597610       0.597610   
286  287ALA         286  0.030000      0.030000   0.796813       0.796813   
16    17LEU          16  0.056923      0.047940  13.147410       0.796813   
297  298ARG         297  0.063750      0.060917   3.386454       0.199203   
20    21GLY          20  0.077143      0.090351   3.984064       0.796813   
300  301ILE         300  0.086809      0.071014  13.745020       2.589641   
296  297PHE         296  0.098095      0.099696  17.131474       6.772908   

     Lipid Count  Lipid Count std       Koff  Residence Time  ...  \
282     1.000000         0.000000  32.882928        0.030411  ...   
293     0.500000         0.500000  30.652487        0.032624  ...   
286     0.500000         0.500000  21.240404        0.047080  ...   
16      1.014286         0.014286  19.295543        0.051825  ...   
297     1.000000         0.000000  10.535137        0.094920  ...   
20      1.000000         0.000000   9.719108        0.102890  ...   
300     1.000000         0.000000   9.262659        0.107960  ...   
296     1.008333         0.008333   3.107937        0.321757  ...   

     R Squared Bootstrap avg  Binding Site ID  Binding Site Duration  \
282                 0.693587                3                0.08427   
293                 0.794795                3                0.08427   
286                 0.492717                3                0.08427   
16                  0.999266                3                0.08427   
297                 0.994440                3                0.08427   
20                  0.991340                3                0.08427   
300                 0.998140                3                0.08427   
296                 0.997092                3                0.08427   

     Binding Site Occupancy  Binding Site Lipid Count  Binding Site Koff  \
282               36.055777                  1.086764           2.489836   
293               36.055777                  1.086764           2.489836   
286               36.055777                  1.086764           2.489836   
16                36.055777                  1.086764           2.489836   
297               36.055777                  1.086764           2.489836   
20                36.055777                  1.086764           2.489836   
300               36.055777                  1.086764           2.489836   
296               36.055777                  1.086764           2.489836   

     Binding Site Koff Bootstrap avg  Binding Site Residence Time  \
282                          6.14596                     0.401633   
293                          6.14596                     0.401633   
286                          6.14596                     0.401633   
16                           6.14596                     0.401633   
297                          6.14596                     0.401633   
20                           6.14596                     0.401633   
300                          6.14596                     0.401633   
296                          6.14596                     0.401633   

     Binding Site R Squared  Binding Site R Squared Bootstrap avg  
282                0.998424                              0.998174  
293                0.998424                              0.998174  
286                0.998424                              0.998174  
16                 0.998424                              0.998174  
297                0.998424                              0.998174  
20                 0.998424                              0.998174  
300                0.998424                              0.998174  
296                0.998424                              0.998174  

```

--------------------------------

### Plot Binding Site Properties vs. ID

Source: https://github.com/wlsong/pylipid/blob/master/docs/demo.md

Generates scatter plots for Residence Time, Duration, Occupancy, and Lipid Count against Binding Site ID. Requires pre-defined 'li' object and 'plt' (matplotlib.pyplot).

```python
ylabel_timeunit = 'ns' if li.timeunit == "ns" else r"$mu$s"
ylabel_dict = {"Residence Time": "Residence Time ({})".format(ylabel_timeunit),
               "Duration": "Duration ({})".format(ylabel_timeunit),
               "Occupancy": "Occuoancy (100%)",
               "Lipid Count": "Lipid Count (num.)"}

# plot No. 1
binding_site_IDs = np.sort(
         [int(bs_id) for bs_id in li.dataset["Binding Site ID"].unique() if bs_id != -1])
for item in ["Residence Time", "Duration", "Occupancy", "Lipid Count"]:
    item_values = np.array(
              [li.dataset[li.dataset["Binding Site ID"]==bs_id]["Binding Site {}".format(item)].unique()[0]
               for bs_id in binding_site_IDs])
    fig, ax = plt.subplots(1, 1, figsize=(len(li.node_list)*0.5, 2.6))
    ax.scatter(np.arange(len(item_values)), np.sort(item_values)[::-1], s=50, color="red")
    ax.set_xticks(np.arange(len(item_values)))
    sorted_index = np.argsort(item_values)[::-1]
    ax.set_xticklabels(binding_site_IDs[sorted_index])
    ax.set_xlabel("Binding Site ID", fontsize=12)
    ax.set_ylabel(ylabel_dict[item], fontsize=12)
    for label in ax.xaxis.get_ticklabels()+ax.yaxis.get_ticklabels():
        plt.setp(label, fontsize=12, weight="normal")
    plt.tight_layout()
    plt.savefig("{}/{}_{}_v_binding_site.{}".format(li.save_dir, li.lipid, "_".join(item.split()), fig_format),
                dpi=200)
    plt.close()
```

--------------------------------

### Test Distance Cutoffs for Lipid-Protein Interactions

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/1-Distance_cutoff_determination.ipynb

This snippet sets up and runs simulations to test different distance cutoffs for lipid-protein interactions. It saves the simulation data and results for further analysis.

```python
import pickle

stride = 1   # use this to speed up if you have too many frames.
lipid = "CHOL"
lipid_atoms = None
nprot = 1
timeunit = "us"

# A2aR simulations
save_dir = check_dir(suffix="test_cutoffs_A2a", print_info=False)

trajfile_list = ["./traj_data/A2a/run1/protein_lipids.xtc", 
                 "./traj_data/A2a/run2/protein_lipids.xtc"]
topfile_list = ["./traj_data/A2a/run1/protein_lipids.gro", 
                 "./traj_data/A2a/run2/protein_lipids.gro"]
num_of_binding_sites, duration_avgs, num_of_contacting_residues = test_cutoffs(
                                 cutoff_list, trajfile_list, topfile_list, lipid, lipid_atoms, 
                                 nprot=nprot, stride=stride, save_dir=save_dir, timeunit=timeunit)


test_data = {"num_of_binding_sites": num_of_binding_sites,
             "duration_avgs": duration_avgs,
             "num_of_contacting_residues": num_of_contacting_residues,
             "test_cutoff_list": cutoff_list}

with open(f"{save_dir}/test_cutoff_data_{lipid}.pickle", "wb") as f:
    pickle.dump(test_data, f, 2)

```

--------------------------------

### Plot Surface Area vs. Binding Site ID

Source: https://github.com/wlsong/pylipid/blob/master/docs/demo.md

Generates a scatter plot of average surface area values against Binding Site ID. Requires 'surface_area_data' and 'li' objects. Assumes 'Binding Site ID' column exists and is numeric.

```python
# plot No. 3
surface_area_averages = np.array(
                       [surface_area_data["Binding Site {}".format(bs_id)].dropna(inplace=False).mean()
                        for bs_id in binding_site_IDs])
fig, ax = plt.subplots(1, 1, figsize=(len(li.node_list)*0.5, 2.6))
ax.scatter(np.arange(len(surface_area_averages)), np.sort(surface_area_averages)[::-1], s=50, color="red")
ax.set_xticks(np.arange(len(surface_area_averages)))
sorted_index = np.argsort(surface_area_averages)[::-1]
ax.set_xticklabels(binding_site_IDs[sorted_index])
ax.set_xlabel("Binding Site ID", fontsize=12)
ax.set_ylabel(r"Surface Area (nm$^2$)", fontsize=12)
for label in ax.xaxis.get_ticklabels()+ax.yaxis.get_ticklabels():
    plt.setp(label, fontsize=12, weight="normal")
plt.tight_layout()
plt.savefig("{}/{}_surface_area_v_binding_site.{}".format(li.save_dir, li.lipid, fig_format), dpi=200)
plt.close()
```

--------------------------------

### Plotting Interaction Metrics

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Iterates through a list of interaction metrics and plots each one as a function of protein residue index. Requires an initialized LipidInteraction object 'li'.

```python
for item in ["Residence Time", "Duration", "Occupancy", "Lipid Count"]:
    li.plot(item=item) # plot values as a function of time
```

--------------------------------

### Plotting Interaction Logo Figures

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Generates a logo figure for a specified interaction metric, highlighting amino acids that make contacts. Requires an initialized LipidInteraction object 'li'.

```python
li.plot_logo(item="Residence Time")
```

--------------------------------

### Compute Binding Sites

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Calculates binding sites by identifying community structures in the protein residue interaction network. The threshold parameter helps filter out spurious correlations by setting a minimum number of residues for a binding site.

```python
node_list, modularity = li.compute_binding_nodes(threshold=4, print_data=False)
```

--------------------------------

### PyLipID Calculation Progress Output

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/1-Distance_cutoff_determination.ipynb

This output shows the progress of PyLipID calculations for collecting interactions and computing binding sites across different trajectory files and cutoff values. It indicates the completion percentage and time taken for each step.

```text
COLLECT INTERACTIONS FROM TRAJECTORIES: 100%|██████████| 2/2 [02:58<00:00, 89.43s/it]
CALCULATE INTERACTIONS FOR BINDING SITES: 100%|██████████| 23/23 [00:13<00:00,  1.69it/s]
COLLECT INTERACTIONS FROM TRAJECTORIES: 100%|██████████| 2/2 [03:12<00:00, 96.04s/it]
CALCULATE INTERACTIONS FOR BINDING SITES: 100%|██████████| 16/16 [00:09<00:00,  1.62it/s]
COLLECT INTERACTIONS FROM TRAJECTORIES: 100%|██████████| 2/2 [03:06<00:00, 93.01s/it]
CALCULATE INTERACTIONS FOR BINDING SITES: 100%|██████████| 11/11 [00:07<00:00,  1.43it/s]

```

--------------------------------

### Compute Minimum Distances of PIP2 to A2aR

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/1-Distance_cutoff_determination.ipynb

Loads trajectory and topology files to compute and save the minimum distances between PIP2 headgroup beads and the A2aR protein. Requires PyTraj and PyLipID. The lipid_atoms parameter specifies which beads of the lipid to consider.

```python
import mdtraj as md
from pylipid import compute_minimum_distance
from pylipid.utils import check_dir

trajfile = "./traj_data/A2a/run1/protein_lipids.xtc"
topfile = "./traj_data/A2a/run1/protein_lipids.gro"
lipid = "POP2" # residue name of PIP2
lipid_atoms = ["C1", "C2", "C3", "PO4", "P1", "P2"] # PIP2 headgroup beads
nprot = 1
save_dir = "test_minimum_dist_A2a_{}".format(lipid)
fig_dir = check_dir(save_dir, "Figures", print_info=False)
contact_frames = 10  # will only plot data if the contact was formed over ${contact_frames} frames.
distance_threshold = 0.65

traj = md.load(trajfile, top=topfile)
minimum_distance_set = compute_minimum_distance(traj, lipid, fig_dir, lipid_atoms=lipid_atoms,
                                               contact_frames=10, distance_threshold=0.65)
```

--------------------------------

### Plot PIP2 Interaction Metrics vs. Cutoff

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/1-Distance_cutoff_determination.ipynb

Generates plots for the number of binding sites, average durations, and number of contacting residues as a function of distance cutoffs for PIP2 interactions. Use this to visualize how interaction metrics change with different cutoff values.

```python
graph(cutoff_list, [num_of_binding_sites[cutoffs] for cutoffs in cutoff_list], 
      "num. of binding sites", lipid, f"{save_dir}/test_cutoff_num_of_bs_{lipid}.png")

graph(cutoff_list, [duration_avgs[cutoffs] for cutoffs in cutoff_list], 
      f"Durations ({timeunit})", lipid, f"{save_dir}/test_cutoff_durations_{lipid}.png")

graph(cutoff_list, [num_of_contacting_residues[cutoffs] for cutoffs in cutoff_list], 
      "num. of contacting residues", lipid, 
      f"{save_dir}/test_cutoff_num_of_contacting_residues_{lipid}.png")
```

--------------------------------

### Plot RMSD vs. Binding Site ID

Source: https://github.com/wlsong/pylipid/blob/master/docs/demo.md

Generates a scatter plot of average RMSD values against Binding Site ID. Requires 'pose_rmsd_data' and 'li' objects. Assumes 'Binding Site ID' column exists and is numeric.

```python
# plot No. 2
binding_site_IDs_RMSD = np.sort([int(bs_id) for bs_id in binding_site_IDs
                                        if f"Binding Site {bs_id}" in pose_rmsd_data.columns])
RMSD_averages = np.array(
                     [pose_rmsd_data[f"Binding Site {bs_id}"].dropna(inplace=False).mean()
                      for bs_id in binding_site_IDs_RMSD])
fig, ax = plt.subplots(1, 1, figsize=(len(li.node_list)*0.5, 2.6))
ax.scatter(np.arange(len(RMSD_averages)), np.sort(RMSD_averages)[::-1], s=50, color="red")
ax.set_xticks(np.arange(len(RMSD_averages)))
sorted_index = np.argsort(RMSD_averages)[::-1]
ax.set_xticklabels(binding_site_IDs_RMSD[sorted_index])
ax.set_xlabel("Binding Site ID", fontsize=12)
ax.set_ylabel("RMSD (nm)", fontsize=12)
for label in ax.xaxis.get_ticklabels()+ax.yaxis.get_ticklabels():
    plt.setp(label, fontsize=12, weight="normal")
plt.tight_layout()
plt.savefig("{}/{}_RMSD_v_binding_site.{}".format(li.save_dir, li.lipid, fig_format), dpi=200)
plt.close()
```

--------------------------------

### Display Residue Dataset

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Prints the initial structure of the dataset attribute within the LipidInteraction object. This DataFrame contains residue names and IDs, and will be populated with interaction data after calculations.

```python
print(li.dataset)

```

--------------------------------

### Collect Residue Contact Data

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Prepares the data for residue interaction analysis by collecting contact information between lipids and protein residues. This step is necessary before calculating interaction metrics like duration or occupancy.

```python
li.collect_residue_contacts()

```

--------------------------------

### Display Top Poses RMSD Head

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Displays the first few rows of the pose_rmsd DataFrame, which contains the Root Mean Square Deviation (RMSD) of bound poses to the binding site average. This is useful for quickly inspecting the RMSD values.

```python
pose_rmsd.head()
```

--------------------------------

### Plot Residence Time vs. Surface Area

Source: https://github.com/wlsong/pylipid/blob/master/docs/demo.md

Generates a scatter plot comparing Residence Time and Surface Area for binding sites. Requires 'li' object and 'surface_area_data'.

```python
# plot No. 5
res_time_BS = np.array(
                  [li.dataset[li.dataset["Binding Site ID"]==bs_id]["Binding Site Residence Time"].unique()[0]
                   for bs_id in binding_site_IDs])
fig, ax = plt.subplots(1, 1, figsize=(len(li.node_list)*0.5, 2.6))
ax.scatter(res_time_BS, surface_area_averages, s=50, color="red")
```

--------------------------------

### Print Durations of First Residue

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Prints the detailed duration data for the first residue. The output shows a list of lists, where each inner list corresponds to the durations from a single trajectory.

```python
print(durations[0])
```

--------------------------------

### Compute Surface Area with Plotting

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Calculates the surface area for binding sites and generates time-series and violin plots if `plot_data` is True. The results are stored in a pandas DataFrame.

```python
surface_area = li.compute_surface_area(plot_data=True)
```

--------------------------------

### Plot Residence Time vs. RMSD

Source: https://github.com/wlsong/pylipid/blob/master/docs/demo.md

Generates a scatter plot comparing Residence Time and RMSD for binding sites. Requires 'li' object and 'pose_rmsd_data'.

```python
# plot No. 4
res_time_BS = np.array(
                  [li.dataset[li.dataset["Binding Site ID"]==bs_id]["Binding Site Residence Time"].unique()[0]
                   for bs_id in binding_site_IDs_RMSD])
fig, ax = plt.subplots(1, 1, figsize=(len(li.node_list)*0.5, 2.6))
ax.scatter(res_time_BS, RMSD_averages, s=50, color="red")
ax.set_xlabel(ylabel_dict["Residence Time"], fontsize=12)
ax.set_ylabel("RMSD (nm)", fontsize=12)
for label in ax.xaxis.get_ticklabels()+ax.yaxis.get_ticklabels():
    plt.setp(label, fontsize=12, weight="normal")
plt.tight_layout()
plt.savefig("{}/{}_Residence_Time_v_RMSD.{}".format(li.save_dir, li.lipid, fig_format), dpi=200)
plt.close()
```

--------------------------------

### Compute Surface Area and Save Data

Source: https://github.com/wlsong/pylipid/blob/master/docs/demo.md

Computes the surface area of binding sites and saves various analysis results, including pose RMSD data and surface area data, to CSV files. Ensure the save directory and lipid name are correctly set.

```python
surface_area_data = li.compute_surface_area(binding_site_id=None, radii=radii, fig_format=fig_format)
data_dir = check_dir(li.save_dir, "Dataset_{}".format(li.lipid))
pose_rmsd_data.to_csv("{}/Pose_RMSD_data.csv".format(data_dir), index=False, header=True)
surface_area_data.to_csv("{}/Surface_Area_data.csv".format(data_dir), index=True, header=True)
li.write_site_info(sort_residue="Residence Time")
```

--------------------------------

### Print Residue and Trajectory Counts

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/0-application-walk-through.ipynb

Prints the total number of residues and trajectory files. It also shows the dimensions of the computed durations, indicating the number of residues and the number of trajectories.

```python
print(len(li.residue_list), len(li.trajfile_list))
print(len(durations), len(durations[0]))
```

--------------------------------

### Plot PDF of Minimum Distances

Source: https://github.com/wlsong/pylipid/blob/master/docs/tutorials/1-Distance_cutoff_determination.ipynb

Generates and saves a Probability Density Function (PDF) plot for the computed minimum distances. This helps visualize the distribution and identify common interaction distances. Requires NumPy and PyLipID plotting functions.

```python
import numpy as np
from pylipid import plot_PDF

distance_set = np.concatenate(minimum_distance_set)
num_of_bins = 1000
fig_fn = "{}/dist_distribut_{}.png".format(save_dir, lipid)
plot_PDF(distance_set, num_of_bins, fig_fn)
```
Tokens:
Raw