<jupyter_start><jupyter_text>Clustering ExampleIn this example we cluster a short trajectory (1000 frames) of the disordered peptide[hiAPP](https://www.ncbi.nlm.nih.gov/pubmed/24021023)We create a normalized covariance matrix using four different metrics:    - Radius of gyration    - Exposed solvent surface    - Asphericity    - End-to-end distanceTable of ContentsLoad Environment Donwload Data  Loading the Trajectory  Quick Trajectory Visualization  Clustering  Quick View of the Clustering Tree  Extract PDB Files for Representative Structures  Calculation of X-Ray Profiles with CRYSOL  "Experimental" X-Ray profile  Fit the Tree Against the Experimental Profile  Analysis of the Tree Level with Best Fit to Experimental Profile  Weight of Each Cluster   (Top)Update EnvironmentPython packages specific to the example notebooks may be needed to install. Use either `pip` or `conda`.**Packages installed with conda**: [conda_requirements.yml](https://raw.githubusercontent.com/jmborr/idpflex/master/notebooks/conda_requirements.yml)    `conda env update -f conda_requirements.yml`**Packages installed with pip**: [pip_requirements.txt](https://raw.githubusercontent.com/jmborr/idpflex/master/notebooks/pip_requirements.txt)   `pip install -r pip_requirements.txt` (Top)Load Environment<jupyter_code>%matplotlib notebook
import matplotlib.pyplot as plt
plt.ion()

import os
import sys
import subprocess
import numpy as np
import MDAnalysis as mda
import nglview
from tqdm import tqdm
import pathos
import multiprocessing
import numpy as np
import scipy
from scipy.cluster.hierarchy import dendrogram
from distutils.spawn import find_executable

from idpflex.cnextend import load_tree
from idpflex.cluster import cluster_with_properties
from idpflex.properties import (RadiusOfGyration, EndToEnd, SaSa, Asphericity,
                               SaxsProperty, propagator_size_weighted_sum)
from idpflex.utils import write_frame
from idpflex.bayes import fit_to_depth<jupyter_output><empty_output><jupyter_text>(Top)Donwload DataIt's assumed git is installed in your system. Otherwise,[follow instructions](http://idpflex.readthedocs.io/en/latest/installation.htmltesting-tutorials-data)to download and unpack your data to /tmp/idpflex_data.<jupyter_code>%%bash
idpflex_data_dir="/tmp/idpflex_data"
if [ -d "${idpflex_data_dir}" ]; then
    cd ${idpflex_data_dir}
    git pull --rebase
else
    git clone https://github.com/jmborr/idpflex_data ${idpflex_data_dir}
fi
idpflex_data_dir = '/tmp/idpflex_data'
data_dir = os.path.join(idpflex_data_dir, 'data', 'simulation')
print(data_dir)<jupyter_output><empty_output><jupyter_text>(Top)Loading the Trajectory<jupyter_code>simulation = mda.Universe(os.path.join(data_dir, 'hiAPP.pdb'),
                          os.path.join(data_dir, 'hiAPP.xtc'))
print('Number of frames in trajectory is ', simulation.trajectory.n_frames)<jupyter_output><empty_output><jupyter_text>(Top)Quick Trajectory Visualization<jupyter_code>w_show = nglview.show_mdanalysis(simulation)
w_show<jupyter_output><empty_output><jupyter_text>(Top)ClusteringWe cluster usin four different scalar properties* Radius of gyration* End to end distance* Solvent accessible surface area* Asphericity<jupyter_code>properties = [RadiusOfGyration, EndToEnd, SaSa, Asphericity]
cl = cluster_with_properties(simulation, properties,
                             segment_length=100, n_representatives=100)<jupyter_output><empty_output><jupyter_text>(Top)Quick View of the Clustering Tree<jupyter_code>plt.title('Hierarchical Clustering Dendrogram')
plt.xlabel('sample index')
plt.ylabel('RMSD (Angstroms)')
dendrogram(cl.tree.z,
           truncate_mode='lastp',  # show only the last p merged clusters
           p=10,  # show this many cluster at the bottom of the tree
           show_leaf_counts=False,  # otherwise numbers in brackets are counts
           leaf_rotation=90.,
           leaf_font_size=12.,
           show_contracted=True,  # to get a distribution impression in truncated branches
          )<jupyter_output><empty_output><jupyter_text>(Top)Extract PDB Files for Representative StructuresWe extract PDB files for each of the 100 representatives and store under directory `/tmp/PDB`<jupyter_code>pdb_names = ['/tmp/PDB/conf_{}.pdb'.format(idx) for idx in cl.idx]

subprocess.call(['mkdir', '-p', '/tmp/PDB'])  # directory to store the PDB files
for idx, name in tqdm(list(zip(cl.idx, pdb_names))):
    write_frame(simulation, idx, name)<jupyter_output><empty_output><jupyter_text>(Top)Calculation of X-Ray Profiles with CRYSOLIt is assumed that `crysol` is installed in your computer, otherwise we fetch the output CRYSOL files from the `idpflex_data` repository. We store a profile for each representative in directory `/tmp/CRYSOL`.<jupyter_code>crysol_names = ['/tmp/CRYSOL/conf_{}.int'.format(idx) for idx in cl.idx]

if find_executable('crysol') is None:
    subprocess.call('cp /tmp/idpflex_data/data/simulation/CRYSOL.tar.gz /tmp'.split())
    subprocess.call('tar zxf /tmp/CRYSOL.tar.gz -C /tmp'.split())
    profiles = [SaxsProperty().from_crysol_int(name) for name in crysol_names]
else:
    pool = pathos.pools.ProcessPool(processes=multiprocessing.cpu_count())
    profiles = list(tqdm(pool.map(SaxsProperty().from_crysol_pdb, pdb_names), total=len(pdb_names)))
[profile.to_ascii(name) for profile, name in zip(profiles, crysol_names)]
propagator_size_weighted_sum(profiles, cl.tree)<jupyter_output><empty_output><jupyter_text>(Top)"Experimental" X-Ray profileWe do not have an experimental profile], so we are going to create a fake experimental profile using the profiles from some of the nodes. The task for the fit engine will be to identify which nodes did we use.Starting from the top of the tree (the root node), we will descend to `level=5`, which contains 6 nodes (the first level is the root node corresponding to `level=0`) We will assign different weights to each of the seven profiles and construct our profile with these weigths.The profile will be stored as a [SAXS property](http://idpflex.readthedocs.io/en/latest/idpflex/properties.htmlidpflex.properties.SaxsProperty)<jupyter_code>nodes = cl.tree.nodes_at_depth(5)
weights = np.asarray([0.00, 0.13, 0.00, 0.55, 0.32, 0.00])  # the weights add up to one
# x are the Q-values
x = nodes[0]['saxs'].x
# y are the intensities
y = np.sum(weights.reshape((6, 1)) * np.asarray([n['saxs'].y for n in nodes]), axis=0)
# Errors simple taken as 10% of the intensities
e = y * 0.1
# Now we create our X-Ray property
exp_saxs = SaxsProperty(qvalues=x, profile=y, errors=e)<jupyter_output><empty_output><jupyter_text>we can plot the property<jupyter_code>fig, ax = plt.subplots(1,1)
ax.plot(exp_saxs.x, exp_saxs.y)
ax.set_xlabel('Q', size=25)
ax.set_ylabel('Intensity', size=25)
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>(Top)Fit the Tree Against the Experimental ProfileStarting from the root node, we fit each tree level against the experimental profile, up to a maximum depth (in this case, `level=7`. Then we will inquire the goodnes of fit for each level<jupyter_code>fits = fit_to_depth(cl.tree, exp_saxs, exp_saxs.name, max_depth=7)<jupyter_output><empty_output><jupyter_text>`fits` is a list of [ModelResult](https://lmfit.github.io/lmfit-py/model.htmllmfit.model.ModelResult) instances, one result for every level. We extract the goodness of fit `\chi^2` and plot versus level<jupyter_code>chi2 = [fit.redchi for fit in fits]
fig, ax = plt.subplots(1,1)
ax.set_xlabel('level', size=25)
ax.set_ylabel('Chi-squared', size=25)
ax.set_yscale('log')
ax.plot(chi2)
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>the steep drop in orders of magnitude for $\chi^2$ at `level=5` indicates the fit engine successfully fitted the experimental profile. (Top)Analysis of the Tree Level with Best Fit to Experimental Profile<jupyter_code>best_fit = fits[5]<jupyter_output><empty_output><jupyter_text>(Top)Weight of Each ClusterWe inquire the weight that the fit engine assigned to each of the seven clusters of `level=6`<jupyter_code>for key in best_fit.best_values:
    if 'amplitude' in key:
        print(key, '{:4.2f}'.format(best_fit.best_values[key]))
print(['{:4.2f}'.format(x) for x in weights])  # weights used to construct the experimental profile<jupyter_output><empty_output><jupyter_text>The order in which the fitted weights are printed is different that the order of the experimental weight. Object `best_fit.best_values` is a python dictionary and order is not guaranteed for this type of object. However, we can use the node id in the amplitude name to sort the fitted weights from smaller to higher node id.The fit procedure correctly identified that only three out of the seven nodes are contributing to the experimental profile. (Top)Representative Structures of the NodesFind a representative structure for each of the three nodes contributing to the match of the experimental profile<jupyter_code>node_ids = [190, 192, 193]  # ID's for the clusters matching the experimental profile 
leafs = [cl.tree[id].representative(cl.rmsd) for id in node_ids]
repr_names = [pdb_names[l.id] for l in leafs]  # representative structures for each node
print(repr_names)
view = nglview.show_file(repr_names[0])
view.display()
view = nglview.show_file(repr_names[1])
view.display()
view = nglview.show_file(repr_names[2])
view.display()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09514687100893997, 'fraction_numerical': 0.010004257130693913, 'mean_word_length': 5.585143658023826, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 61, 'https://': 5, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '166954', 'license': 'MIT', 'n_tokens_mistral': 3127, 'n_tokens_neox': 2966, 'n_words': 890, 'path': 'notebooks/hiAPP_clustering_with_properties.ipynb', 'repo_name': 'jmborr/ipdflex', 'size': '9396'}