<jupyter_start><jupyter_text>Higher order Laplacians for protein atoms Authors: <NAME>1. [Protein 3-D point cloud of atoms](proteins)2. [Constructing simplex tree](tree)3. [Using ```HodgeLaplacinans()``` class](hodge) 4. [References](refs)  AbstractIn this example we will illustrate how to use ```HodgeLaplacians()``` class to get higher order Laplacians for a simplicial complex constructed from the protein data. To access Protein Data Bank we will use [BioPython Library](https://biopython.org/) and to efficiently compute the tree of simplices we will use [C++ GUDHI Library](http://gudhi.gforge.inria.fr/). A variety of alternative [libraries for Topological Data Analysis](https://en.wikipedia.org/wiki/Persistent_homologyComputation) are also available.   Protein point cloud data  [↑](toc)First, to install BioPython and nglview run```pythonpip3 install biopython nglview```To activate the nglview extension for notebook, run before starting Jupyter:```pythonjupyter-nbextension enable nglview --py --sys-prefixjupyter nbextension enable widgetsnbextension --user --py```There is a [nice video tutorial](https://www.youtube.com/watch?v=mL8NPpRxgJA) ([notebook for tutorial](https://nbviewer.jupyter.org/github/cgoliver/Notebooks/blob/master/COMP_364/L27/L27.ipynb)) on BioPython and Protein Data Bank format.<jupyter_code>from Bio.PDB import *
import numpy as np

import warnings
warnings.filterwarnings('ignore')

pdbl_1 = PDBList()
pdbl_1.retrieve_pdb_file('4XP1')

pdbl_2 = PDBList()
pdbl_2.retrieve_pdb_file('3I40')

parser_1 = MMCIFParser()
parser_2 = MMCIFParser()

structure_1 = parser_1.get_structure('4XP1', 'xp/4xp1.cif')
structure_2 = parser_2.get_structure('4XP1', 'i4/3i40.cif')<jupyter_output>WARNING: The default download format has changed from PDB to PDBx/mmCif
WARNING: The default download format has changed from PDB to PDBx/mmCif<jupyter_text>Let us visualize two proteins using the viewer form ```nglview``` module.<jupyter_code>import nglview as nv

view_1 = nv.show_biopython(structure_1)
view_1.clear_representations()
#view as ball and stick (atom and bond)
view_1.add_ball_and_stick()

view_1
view_2 = nv.show_biopython(structure_2)
view_2.clear_representations()
#view as ball and stick (atom and bond)
view_2.add_ball_and_stick()
view_2<jupyter_output><empty_output><jupyter_text>Next we extract the 3-D coordinates of individual atoms in each protein.<jupyter_code>points_1 = []
points_2 = []

for model in structure_1:
    for chain in model:
        for residue in chain:
            for atom in residue:
                points_1.append(tuple(atom.coord))
                
for model in structure_2:
    for chain in model:
        for residue in chain:
            for atom in residue:
                points_2.append(tuple(atom.coord))

points_1 = np.array(points_1)
points_2 = np.array(points_2)

print(f"There are {len(points_1)} atoms in the protein 1.")
print(f"There are {len(points_2)} atoms in the protein 2.")<jupyter_output>There are 7675 atoms in the protein 1.
There are 438 atoms in the protein 2.<jupyter_text>Constructing simplex tree  [↑](toc)In this section we will connect the neighboring points/atoms with edges thus obtaining a graph where each node/atom is connected to its neighbours in 3D. From the graph we derive the Vietoris-Rips simplicial complex simply by adding simplices associated to each clique of up to certain size in a graph.To get the connectivity of a graph we will use C++ GUDHI Library that will provides us simplex tree with filtration values as an output. We will use simplex tree data as input data to our ```HodgeLaplacians()``` class. To install GUDHI Library:1. [Download and compile C++ Sources](http://gudhi.gforge.inria.fr/doc/latest/installation.html)2. [Build Cython module for GUDHI](http://gudhi.gforge.inria.fr/python/latest/installation.html)P.S. Don't forget to install all C/C++ dependencies like Boost, CGAL, Eigen, GMP, CMake,...Let us assume that you've spent happy 30 minutes and compiled and installed the library, so we can proceed with the construction of connectivity and simplex tree to finally feed it into ```HodgeLaplacians()``` class.Construction of local neighbourhood connectivity is computationally intensive task that is why we rely on a separate compiled C++ library and do not incorporate this functionality inside ```HodgeLaplacians()``` class. In Python pairwise distance check for points could take a long time.<jupyter_code>import sys
sys.path.insert(0, "/home/mm/Download/library-gudhi_library_2.3.0/2018-09-04-14-25-00_GUDHI_2.3.0/build/cython")

import gudhi as gd
import time
start_time = time.time()


epsilon = 3 # Connect two points if distance is smaller than epsilon

skeleton_1  = gd.RipsComplex(points = points_1, max_edge_length=epsilon )
rips_simplex_tree_1 = skeleton_1.create_simplex_tree(max_dimension=3) 
rips_skeleton_gudhi_1 = rips_simplex_tree_1.get_skeleton(3)

skeleton_2  = gd.RipsComplex(points = points_2, max_edge_length=epsilon )
rips_simplex_tree_2 = skeleton_2.create_simplex_tree(max_dimension=3) 
rips_skeleton_gudhi_2 = rips_simplex_tree_2.get_skeleton(3)

end_time = time.time()

print(f"Done in {end_time-start_time} seconds")
print(len(rips_skeleton_gudhi_1))<jupyter_output>65155<jupyter_text>Using ```HodgeLaplacinans()``` to get diffusion operators, Ricci curvature and spectrum To install ```hodgelaplacians``` library ~~~pip3 install hodgelaplacians~~~ First let us import the main class ```HodgeLaplacians()``` and initialize it by simplex tree.<jupyter_code>from hodgelaplacians import HodgeLaplacians

start = time.time()

hl_1 = HodgeLaplacians(rips_skeleton_gudhi_1, maxdimension=3, mode='gudhi')
hl_2 = HodgeLaplacians(rips_skeleton_gudhi_2, maxdimension=3, mode='gudhi')


# # Boundary operators
#B0_1 = hl_1.getBoundaryOperator(0)
#B1_1 = hl_1.getBoundaryOperator(1)
#B2_1 = hl_1.getBoundaryOperator(2)
# B3_1 = hl_1.getBoundaryOperator(3)

# # Hodge Laplacians
# L0_1 = hl_1.getHodgeLaplacian(0)
# L1_1 = hl_1.getHodgeLaplacian(1)
# L2_1 = hl_1.getHodgeLaplacian(2)
# L3_1 = hl_1.getHodgeLaplacian(3)

# # Bochner Laplacian 
# # (row sum of absolute values equals to diagonal element)
# LB0_1 = hl_1.getBochnerLaplacian(0)
# LB1_1 = hl_1.getBochnerLaplacian(1)
# LB2_1 = hl_1.getBochnerLaplacian(2)
# LB3_1 = hl_1.getBochnerLaplacian(3)




end = time.time()

print(end-start)

# print(f"0-Boundary operator is \n {B0.toarray()} \n \n")
# print(f"1-Boundary operator is \n {B1.toarray()} \n \n")
# print(f"2-Boundary operator is \n {B2.toarray()} \n \n")

# print(f"0-Hodge Laplacian operator is \n {L0.toarray()} \n \n")
# print(f"1-Hodge Laplacian operator is \n {L1.toarray()} \n \n")
# print(f"2-Hodge Laplacian operator is \n {L2.toarray()} \n \n")

# print(hl.n_faces(0))
# print(hl.n_faces(1))
# print(hl.n_faces(2))<jupyter_output>6.887351989746094<jupyter_text>Hodge and Bochner Laplacian eigenvalues in an interval<jupyter_code>start = time.time()

evalsL1_1, evecsL1_1 = hl_1.getHodgeSpectrum(1, 40, around_point = 0.001)
evalsLB1_1, evecsLB1_1 = hl_1.getBochnerSpectrum(1, 40, around_point = 0.001)

evalsL1_2, evecsL1_2 = hl_2.getHodgeSpectrum(1, 40, around_point = 0.001)
evalsLB1_2, evecsLB1_2 = hl_2.getBochnerSpectrum(1, 40, around_point = 0.001)

end = time.time()
print(end - start)
%matplotlib inline
import pandas as pd
import matplotlib.pyplot as plt

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12,3))

df = pd.DataFrame({'Protein 1':evalsL1_1, 'Protein 2':evalsL1_2})
df[['Protein 1', 'Protein 2']].plot(ax=axes[0], title="Eigenvalues of 1-Hodge Laplacian", grid=True)

df = pd.DataFrame({'Protein 1':evalsLB1_1, 'Protein 2':evalsLB1_2})
df[['Protein 1', 'Protein 2']].plot(ax=axes[1], title="Eigenvalues of 1-Bochner Laplacian", grid=True)<jupyter_output><empty_output><jupyter_text>Distribution of combinatorial Ricci curvature of two proteins<jupyter_code># Combinatorial Ricci curvature
# (correction term for Laplacian to become Bochner Laplacian)
R0_1 = hl_1.getCombinatorialRicci(0)
R1_1 = hl_1.getCombinatorialRicci(1)
R2_1 = hl_1.getCombinatorialRicci(2)
R3_1 = hl_1.getCombinatorialRicci(3)

# # Combinatorial Ricci curvature
# # (correction term for Laplacian to become Bochner Laplacian)
R0_2 = hl_2.getCombinatorialRicci(0)
R1_2 = hl_2.getCombinatorialRicci(1)
R2_2 = hl_2.getCombinatorialRicci(2)
R3_2 = hl_2.getCombinatorialRicci(3)<jupyter_output><empty_output><jupyter_text>Now let us plot the historgrams side-by-side<jupyter_code>plt.figure(figsize=(12,3))
plt.subplot(1, 2, 1)
plt.suptitle('Distribution of combinatorial Ricci-1 curvature in two proteins \n')
plt.xlabel('Combinatorial Ricci Curvature')
plt.ylabel('Frequency')
plt.grid(axis='y', alpha=0.75)
plot_ricci1_1 = pd.Series(R1_1)
plot_ricci1_1.plot.hist(grid=True, bins=20, rwidth=0.9, color='#004c99')
plt.subplot(1, 2, 2)
plt.xlabel('Combinatorial Ricci Curvature')
plt.ylabel('Frequency')
plt.grid(axis='y', alpha=0.75)
plot_ricci1_2 = pd.Series(R1_2)
plot_ricci1_2.plot.hist(grid=True, bins=20, rwidth=0.9, color='#006633')


plt.figure(figsize=(12,3))
plt.suptitle('Distribution of combinatorial Ricci-2 curvature in two proteins \n')
plt.subplot(1, 2, 1)
plt.xlabel('Combinatorial Ricci Curvature')
plt.ylabel('Frequency')
plt.grid(axis='y', alpha=0.75)
plot_ricci1_1 = pd.Series(R2_1)
plot_ricci1_1.plot.hist(grid=True, bins=20, rwidth=0.9, color='#004c99')
plt.subplot(1, 2, 2)
plt.xlabel('Combinatorial Ricci Curvature')
plt.ylabel('Frequency')
plt.grid(axis='y', alpha=0.75)
plot_ricci1_2 = pd.Series(R2_2)
plot_ricci1_2.plot.hist(grid=True, bins=20, rwidth=0.9, color='#006633')


plt.figure(figsize=(12,3))
plt.suptitle('Distribution of combinatorial Ricci-3 curvature in two proteins \n')
plt.subplot(1, 2, 1)
plt.xlabel('Combinatorial Ricci Curvature')
plt.ylabel('Frequency')
plt.grid(axis='y', alpha=0.75)
plot_ricci1_1 = pd.Series(R3_1)
plot_ricci1_1.plot.hist(grid=True, bins=20, rwidth=0.9, color='#004c99')
plt.subplot(1, 2, 2)
plt.xlabel('Combinatorial Ricci Curvature')
plt.ylabel('Frequency')
plt.grid(axis='y', alpha=0.75)
plot_ricci1_2 = pd.Series(R3_2)
plot_ricci1_2.plot.hist(grid=True, bins=20, rwidth=0.9, color='#006633')<jupyter_output><empty_output><jupyter_text>kj  References  [↑](toc)1. kk* [Vector class](http://biopython.org/DIST/docs/api/Bio.PDB.vectors.Vector-class.html)* [Video introduction](https://www.youtube.com/watch?v=mL8NPpRxgJA)    * https://nbviewer.jupyter.org/github/cgoliver/Notebooks/blob/master/COMP_364/L27/L27.ipynb* https://stackoverflow.com/questions/10324674/parsing-a-pdb-file-in-python* https://github.com/maartenbreddels/talk-ipyvolume-scipy2018<jupyter_code>import numpy as np
import ipyvolume as ipv
V = np.zeros((128,128,128)) # our 3d array
# outer box
V[30:-30,30:-30,30:-30] = 0.75
V[35:-35,35:-35,35:-35] = 0.0
# inner box
V[50:-50,50:-50,50:-50] = 0.25
V[55:-55,55:-55,55:-55] = 0.0
ipv.quickvolshow(V, level=[0.25, 0.75], opacity=0.03, level_width=0.1, data_min=0, data_max=1)
import ipyvolume as ipv
import numpy as np
x, y, z = np.random.random((3, 10000))
ipv.quickscatter(x, y, z, size=1, marker="sphere")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12619433928249504, 'fraction_numerical': 0.05074815215431765, 'mean_word_length': 5.904169259489732, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 34, 'https://': 8, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '460038', 'license': 'MIT', 'n_tokens_mistral': 4618, 'n_tokens_neox': 4146, 'n_words': 972, 'path': 'examples/.ipynb_checkpoints/hodge_spectrum_protein_gudhi-checkpoint.ipynb', 'repo_name': 'tsitsvero/hodgelaplacians', 'size': '11105'}