diff --git a/contact_map/__init__.py b/contact_map/__init__.py index b46f5b6..d3f3f44 100644 --- a/contact_map/__init__.py +++ b/contact_map/__init__.py @@ -11,6 +11,8 @@ from .contact_count import ContactCount +from .contact_trajectory import ContactTrajectory, RollingContactFrequency + from .min_dist import NearestAtoms, MinimumDistanceCounter from .concurrence import ( diff --git a/contact_map/contact_map.py b/contact_map/contact_map.py index 428b87e..6af9b6f 100644 --- a/contact_map/contact_map.py +++ b/contact_map/contact_map.py @@ -174,6 +174,24 @@ def __init__(self, topology, query, haystack, cutoff, n_neighbors_ignored): } self._atom_idx_to_residue_idx = self._set_atom_idx_to_residue_idx() + @classmethod + def from_contacts(cls, atom_contacts, residue_contacts, topology, + query=None, haystack=None, cutoff=0.45, + n_neighbors_ignored=2): + obj = cls.__new__(cls) + super(cls, obj).__init__(topology, query, haystack, cutoff, + n_neighbors_ignored) + + def get_contact_counter(contact): + if isinstance(contact, ContactCount): + return contact.counter + else: + return contact + + obj._atom_contacts = get_contact_counter(atom_contacts) + obj._residue_contacts = get_contact_counter(residue_contacts) + return obj + def _set_atom_slice(self): """ Set atom slice logic """ if (self._class_use_atom_slice is None and @@ -774,6 +792,17 @@ def __init__(self, trajectory, query=None, haystack=None, cutoff=0.45, contacts = self._build_contact_map(trajectory) (self._atom_contacts, self._residue_contacts) = contacts + @classmethod + def from_contacts(cls, atom_contacts, residue_contacts, n_frames, + topology, query=None, haystack=None, cutoff=0.45, + n_neighbors_ignored=2): + obj = super(ContactFrequency, cls).from_contacts( + atom_contacts, residue_contacts, topology, query, haystack, + cutoff, n_neighbors_ignored + ) + obj._n_frames = n_frames + return obj + def __hash__(self): return hash((super(ContactFrequency, self).__hash__(), tuple(self._atom_contacts.items()), @@ -951,6 +980,10 @@ def __sub__(self, other): def contact_map(self, *args, **kwargs): #pylint: disable=W0221 raise NotImplementedError + @classmethod + def from_contacts(self, *args, **kwargs): #pylint: disable=W0221 + raise NotImplementedError + @property def atom_contacts(self): n_x = self.topology.n_atoms diff --git a/contact_map/contact_trajectory.py b/contact_map/contact_trajectory.py new file mode 100644 index 0000000..3bda7c8 --- /dev/null +++ b/contact_map/contact_trajectory.py @@ -0,0 +1,368 @@ +from collections import abc, Counter + +from .contact_map import ContactFrequency, ContactObject +import json + +class ContactTrajectory(ContactObject, abc.Sequence): + """Track all the contacts over a trajectory, frame-by-frame. + + Internally, this has a single-frame :class:`.ContactFrequency` for each + frame of the trajectory. + + Parameters + ---------- + trajectory : mdtraj.Trajectory + the trajectory to calculate contacts for + query : list of int + Indices of the atoms to be included as query. Default ``None`` + means all heavy, non-water atoms. + haystack : list of int + Indices of the atoms to be included as haystack. Default ``None`` + means all heavy, non-water atoms. + cutoff : float + Cutoff distance for contacts, in nanometers. Default 0.45. + n_neighbors_ignored : int + Number of neighboring residues (in the same chain) to ignore. + Default 2. + """ + _class_use_atom_slice = None + def __init__(self, trajectory, query=None, haystack=None, cutoff=0.45, + n_neighbors_ignored=2): + super(ContactTrajectory, self).__init__(trajectory.topology, query, + haystack, cutoff, + n_neighbors_ignored) + contacts = self._build_contacts(trajectory) + self._contact_maps = [ + ContactFrequency.from_contacts( + topology=self.topology, + query=self.query, + haystack=self.haystack, + cutoff=self.cutoff, + n_neighbors_ignored=self.n_neighbors_ignored, + atom_contacts=atom_contacts, + residue_contacts=residue_contacts, + n_frames=1 + ) + for atom_contacts, residue_contacts in zip(*contacts) + ] + + def __getitem__(self, num): + return self._contact_maps[num] + + def __len__(self): + return len(self._contact_maps) + + def __hash__(self): + return hash((super(ContactTrajectory, self).__hash__(), + tuple([frozenset(frame.counter.items()) + for frame in self.atom_contacts]), + tuple([frozenset(frame.counter.items()) + for frame in self.residue_contacts]))) + + def __eq__(self, other): + return hash(self) == hash(other) + + @classmethod + def from_contacts(cls, atom_contacts, residue_contacts, topology, + query=None, haystack=None, cutoff=0.45, + n_neighbors_ignored=2): + contact_maps = [ + ContactFrequency.from_contacts( + atom_cs, + res_cs, + n_frames=1, + topology=topology, + query=query, + haystack=haystack, + cutoff=cutoff, + n_neighbors_ignored=n_neighbors_ignored + ) + for atom_cs, res_cs in zip(atom_contacts, residue_contacts) + ] + return cls.from_contact_maps(contact_maps) + + def _build_contacts(self, trajectory): + # atom_contacts, residue_contacts = self._empty_contacts() + atom_contacts = [] + residue_contacts = [] + + residue_ignore_atom_idxs = self.residue_ignore_atom_idxs + residue_query_atom_idxs = self.residue_query_atom_idxs + used_trajectory = self.slice_trajectory(trajectory) + + # range(len(trajectory)) avoids recopying topology, as would occur + # in `for frame in trajectory` + for frame_num in range(len(trajectory)): + frame_contacts = self.contact_map(used_trajectory, frame_num, + residue_query_atom_idxs, + residue_ignore_atom_idxs) + frame_atom_contacts, frame_residue_contacts = frame_contacts + frame_atom_contacts = \ + self.convert_atom_contacts(frame_atom_contacts) + # TODO unify contact building with something like this? + # atom_contacts, residue_contact = self._update_contacts(...) + atom_contacts.append(frame_atom_contacts) + residue_contacts.append(frame_residue_contacts) + return atom_contacts, residue_contacts + + def contact_frequency(self): + """Create a :class:`.ContactFrequency` from this contact trajectory + """ + freq = ContactFrequency.from_contacts( + atom_contacts=Counter(), + residue_contacts=Counter(), + n_frames=0, + topology=self.topology, + query=self.query, + haystack=self.haystack, + cutoff=self.cutoff, + n_neighbors_ignored=self.n_neighbors_ignored + ) + for cmap in self._contact_maps: + # TODO: skipping compatibility checks would help performance; we + # know that everything in here *should* be compatible + freq.add_contact_frequency(cmap) + + return freq + + def to_dict(self): + return { + 'contact_maps': [cmap.to_dict() for cmap in self._contact_maps] + } + + @classmethod + def from_dict(cls, dct): + contact_maps = [ContactFrequency.from_dict(cmap) + for cmap in dct['contact_maps']] + obj = cls.from_contact_maps(contact_maps) + return obj + + @property + def atom_contacts(self): + return [cmap.atom_contacts for cmap in self._contact_maps] + + @property + def residue_contacts(self): + return [cmap.residue_contacts for cmap in self._contact_maps] + + @classmethod + def from_contact_maps(cls, maps): + obj = cls.__new__(cls) + super(cls, obj).__init__(maps[0].topology, maps[0].query, + maps[0].haystack, maps[0].cutoff, + maps[0].n_neighbors_ignored) + + for cmap in maps: + obj._check_compatibility(cmap) + + obj._contact_maps = maps + return obj + + @classmethod + def join(cls, others): + """Concatenate ContactTrajectory instances + + Parameters + ---------- + others : List[:class:.ContactTrajectory] + contact trajectories to concatenate + + Returns + ------- + :class:`.ContactTrajectory` : + concatenated contact trajectory + """ + contact_maps = sum([o._contact_maps for o in others], []) + return cls.from_contact_maps(contact_maps) + + def rolling_frequency(self, window_size=1, step=1): + """:class:`.RollingContactFrequency` iterator for this trajectory + + Parameters + ---------- + window_size : int + the number of frames in the window + step : int + the number of frames between successive starting points of the + window (like the ``step`` parameter in a Python slice object) + + Returns + ------- + :class:`.RollingContactFrequency` : + windowed iterator for this trajectory + """ + return RollingContactFrequency(self, width=window_size, step=step) + + +class MutableContactTrajectory(ContactTrajectory, abc.MutableSequence): + """Mutable version of :class:`.ContactTrajectory` + + Parameters + ---------- + trajectory : mdtraj.Trajectory + the trajectory to calculate contacts for + query : list of int + Indices of the atoms to be included as query. Default ``None`` + means all heavy, non-water atoms. + haystack : list of int + Indices of the atoms to be included as haystack. Default ``None`` + means all heavy, non-water atoms. + cutoff : float + Cutoff distance for contacts, in nanometers. Default 0.45. + n_neighbors_ignored : int + Number of neighboring residues (in the same chain) to ignore. + Default 2. + + """ + def __setitem__(self, key, value): + self._contact_maps[key] = value + + def __delitem__(self, key): + del self._contact_maps[key] + + def insert(self, key, value): + self._contact_maps.insert(key, value) + + def __hash__(self): + # mutable objects must have unique hashes + return id(self) + + +class WindowedIterator(abc.Iterator): + """ + Helper for windowed ("rolling average") iterators. + + The idea is that this is an easy and reusable code for getting windowed + quantitiies such as needed for rolling averages. This iterator itself + just returns sets of indices/slices to add/remove from whatever counter + is being tracked. The idea is that it will be used inside of another + iterator. + + + Parameters + ---------- + length : int + the length of the list windowed over + width : int + the number of items in the window + step : int + the number of items skipped between successive windows (as with the + ``step`` parameter in slices) + slow_build : bool + if True, the iterator builds up the window "step" objects at a time. + Otherwise, the first value is the full width of the window. + + Attributes + ---------- + min : int + the index of the first object in the cached window + max : int + the index of the last object in the cached window (note that this is + included in the window, unlike Python slices) + """ + def __init__(self, length, width, step, slow_build): + self.length = length + self.width = width + self.step = step + self.slow_build = slow_build + self.min = -1 + self.max = -1 + + def _startup(self): + to_sub = slice(0, 0) + self.min = max(self.min, 0) + if self.slow_build: + to_add = slice(self.max + 1, self.max + self.step + 1) + self.max += self.step + else: + self.max = self.width - 1 + to_add = slice(self.min, self.max + 1) + return to_add, to_sub + + def _normal(self): + self.min = max(0, self.min) + new_max = self.max + self.step + + if not self.slow_build: + new_max = max(new_max, self.width - 1) + + new_min = max(self.min, new_max - self.width + 1) + + to_sub = slice(self.min, new_min) + to_add = slice(self.max + 1, new_max + 1) + self.min = new_min + self.max = new_max + return to_add, to_sub + + def __next__(self): + # if self.max + self.step < self.width: + # to_add, to_sub = self._startup() + if self.max + self.step < self.length: + to_add, to_sub = self._normal() + else: + raise StopIteration + + return to_add, to_sub + + +class RollingContactFrequency(abc.Iterator): + """Iterator for "rolling-average" contact frequencies over a trajectory + + Parameters + ---------- + contact_trajectory : :class:`.ContactTrajectory` + input trajectory + width : int + the number of frames in the window + step : int + the number of frames between successive starting points of the + window (like the ``step`` parameter in a Python slice object) + """ + + _slow_build_iter = False + + def __init__(self, contact_trajectory, width=1, step=1): + self.trajectory = contact_trajectory + self.width = width + self.step = step + self.slow_build_iter = self._slow_build_iter + self._window_iter = None + self._contact_map = None + + def __iter__(self): + self._window_iter = WindowedIterator(length=len(self.trajectory), + width=self.width, + step=self.step, + slow_build=self.slow_build_iter) + self._contact_map = ContactFrequency.from_contacts( + Counter(), Counter(), + topology=self.trajectory.topology, + query=self.trajectory.query, + haystack=self.trajectory.haystack, + cutoff=self.trajectory.cutoff, + n_neighbors_ignored=self.trajectory.n_neighbors_ignored, + n_frames=0 + ) + return self + + def __next__(self): + to_add, to_sub = next(self._window_iter) + for frame in self.trajectory[to_add]: + self._contact_map.add_contact_frequency(frame) + for frame in self.trajectory[to_sub]: + self._contact_map.subtract_contact_frequency(frame) + + # need to make a copy in case the user does list(rolling_freq), + # otherwise they get copies of only the last version! + cmap = self._contact_map + map_copy = ContactFrequency.from_contacts( + cmap._atom_contacts.copy(), + cmap._residue_contacts.copy(), + topology=cmap.topology, + query=cmap.query, + haystack=cmap.haystack, + cutoff=cmap.cutoff, + n_neighbors_ignored=cmap.n_neighbors_ignored, + n_frames=cmap.n_frames + ) + return map_copy diff --git a/contact_map/tests/test_contact_map.py b/contact_map/tests/test_contact_map.py index 2d02445..11a2c58 100644 --- a/contact_map/tests/test_contact_map.py +++ b/contact_map/tests/test_contact_map.py @@ -12,7 +12,7 @@ # stuff to be testing in this file from contact_map.contact_map import * -from contact_map.contact_count import HAS_MATPLOTLIB +from contact_map.contact_count import HAS_MATPLOTLIB, ContactCount traj = md.load(find_testfile("trajectory.pdb")) @@ -158,6 +158,25 @@ def test_counters(self, idx): assert m._residue_contacts == expected assert m.residue_contacts.counter == expected + @pytest.mark.parametrize('contactcount', [True, False]) + def test_from_contacts(self, idx, contactcount): + expected = self.maps[idx] + atom_contact_list = self.expected_atom_contacts[expected] + residue_contact_list = self.expected_residue_contacts[expected] + atom_contacts = counter_of_inner_list(atom_contact_list) + residue_contacts = counter_of_inner_list(residue_contact_list) + if contactcount: + atom_contacts = ContactCount(atom_contacts, self.topology.atom, + 10, 10) + residue_contacts = ContactCount(residue_contacts, + self.topology.residue, 5, 5) + + cmap = ContactMap.from_contacts(atom_contacts, residue_contacts, + topology=self.topology, + cutoff=0.075, + n_neighbors_ignored=0) + _contact_object_compare(cmap, expected) + def test_to_dict(self, idx): m = self.maps[idx] dct = m.to_dict() @@ -374,6 +393,25 @@ def test_counters(self): def test_contacts_dict(self): _check_contacts_dict_names(self.map) + @pytest.mark.parametrize('contactcount', [True, False]) + def test_from_contacts(self, contactcount): + atom_contacts = self.expected_atom_contact_count + residue_contacts = self.expected_residue_contact_count + top = traj.topology + if contactcount: + atom_contacts = ContactCount(atom_contacts, top.atom, + 10, 10) + residue_contacts = ContactCount(residue_contacts, top.residue, + 5, 5) + + cmap = ContactFrequency.from_contacts(atom_contacts, + residue_contacts, + n_frames=5, + topology=top, + cutoff=0.075, + n_neighbors_ignored=0) + _contact_object_compare(cmap, self.map) + def test_check_compatibility_true(self): map2 = ContactFrequency(trajectory=traj[0:2], cutoff=0.075, diff --git a/contact_map/tests/test_contact_trajectory.py b/contact_map/tests/test_contact_trajectory.py new file mode 100644 index 0000000..c9575a4 --- /dev/null +++ b/contact_map/tests/test_contact_trajectory.py @@ -0,0 +1,348 @@ +# pylint: disable=wildcard-import, missing-docstring, protected-access +# pylint: disable=attribute-defined-outside-init, invalid-name, no-self-use +# pylint: disable=wrong-import-order, unused-wildcard-import + +from .utils import * +from .test_contact_map import ( + counter_of_inner_list, _contact_object_compare, traj_atom_contact_count, + traj_residue_contact_count +) + +import mdtraj as md + +from contact_map.contact_trajectory import * +from contact_map.contact_count import ContactCount + +TRAJ_ATOM_CONTACTS = [ + [[1, 4], [4, 6], [5, 6]], + [[1, 5], [4, 6], [5, 6]], + [[1, 4], [4, 6], [5, 6]], + [[1, 4], [4, 6], [4, 7], [5, 6], [5, 7]], + [[0, 9], [0, 8], [1, 8], [1, 9], [1, 4], [8, 4], [8, 5], [4, 6], [4, 7], + [5, 6], [5, 7]] +] + +TRAJ_RES_CONTACTS = [ + [[0, 2], [2, 3]], + [[0, 2], [2, 3]], + [[0, 2], [2, 3]], + [[0, 2], [2, 3]], + [[0, 2], [2, 3], [0, 4], [2, 4]] +] + +class TestContactTrajectory(object): + def setup(self): + self.traj = md.load(find_testfile("trajectory.pdb")) + self.map = ContactTrajectory(self.traj, cutoff=0.075, + n_neighbors_ignored=0) + self.expected_atom_contacts = TRAJ_ATOM_CONTACTS + self.expected_residue_contacts = TRAJ_RES_CONTACTS + + @pytest.mark.parametrize('contact_type', ['atom', 'residue']) + def test_contacts(self, contact_type): + assert len(self.map) == 5 + contacts = {'atom': self.map.atom_contacts, + 'residue': self.map.residue_contacts}[contact_type] + expected = {'atom': self.expected_atom_contacts, + 'residue': self.expected_residue_contacts}[contact_type] + + for contact, expect in zip(contacts, expected): + expected_counter = counter_of_inner_list(expect) + assert contact.counter == expected_counter + + @pytest.mark.parametrize('contact_type', ['atom', 'residue']) + def test_contacts_sliced(self, contact_type): + selected_atoms = [2, 3, 4, 5, 6, 7, 8, 9] + cmap = ContactTrajectory(self.traj, query=selected_atoms, + haystack=selected_atoms, cutoff=0.075, + n_neighbors_ignored=0) + contacts = {'atom': cmap.atom_contacts, + 'residue': cmap.residue_contacts}[contact_type] + expected = { + 'atom': [ + [[4, 6], [5, 6]], + [[4, 6], [5, 6]], + [[4, 6], [5, 6]], + [[4, 6], [4, 7], [5, 6], [5, 7]], + [[8, 4], [8, 5], [4, 6], [4, 7], [5, 6], [5, 7]] + ], + 'residue': [ + [[2, 3]], + [[2, 3]], + [[2, 3]], + [[2, 3]], + [[2, 3], [2, 4]] + ] + }[contact_type] + + for contact, expect in zip(contacts, expected): + expected_counter = counter_of_inner_list(expect) + assert contact.counter == expected_counter + + + @pytest.mark.parametrize('contactcount', [True, False]) + def test_from_contacts(self, contactcount): + atom_contacts = [ + counter_of_inner_list(frame_contacts) + for frame_contacts in self.expected_atom_contacts + ] + residue_contacts = [ + counter_of_inner_list(frame_contacts) + for frame_contacts in self.expected_residue_contacts + ] + top = self.traj.topology + if contactcount: + atom_contacts = [ContactCount(contact, top.atom, 10, 10) + for contact in atom_contacts] + residue_contacts = [ContactCount(contact, top.residue, 5, 5) + for contact in residue_contacts] + + cmap = ContactTrajectory.from_contacts(atom_contacts, + residue_contacts, + topology=top, + cutoff=0.075, + n_neighbors_ignored=0) + for truth, beauty in zip(self.map, cmap): + _contact_object_compare(truth, beauty) + assert truth == beauty + _contact_object_compare(cmap, self.map) + assert cmap == self.map + + + def test_contact_frequency(self): + freq = self.map.contact_frequency() + expected_atom_count = { + key: val / 5.0 for key, val in traj_atom_contact_count.items() + } + expected_res_count = { + key: val / 5.0 + for key, val in traj_residue_contact_count.items() + } + assert freq.atom_contacts.counter == expected_atom_count + assert freq.residue_contacts.counter == expected_res_count + + @pytest.mark.parametrize("intermediate", ["dict", "json"]) + def test_serialization_cycle(self, intermediate): + # NOTE: this is identical to TestContactFrequency; can probably + # abstract it out + serializer, deserializer = { + 'json': (self.map.to_json, ContactTrajectory.from_json), + 'dict': (self.map.to_dict, ContactTrajectory.from_dict) + }[intermediate] + + serialized = serializer() + reloaded = deserializer(serialized) + _contact_object_compare(self.map, reloaded) + assert self.map == reloaded + + def test_from_contact_maps(self): + maps = [ContactFrequency(frame, cutoff=0.075, n_neighbors_ignored=0) + for frame in self.traj] + cmap = ContactTrajectory.from_contact_maps(maps) + _contact_object_compare(self.map, cmap) + assert self.map == cmap + + def test_from_contact_maps_incompatible(self): + map0 = ContactFrequency(self.traj[0], cutoff=0.075, + n_neighbors_ignored=0) + maps = [map0] + [ContactFrequency(frame) for frame in self.traj[1:]] + with pytest.raises(AssertionError): + _ = ContactTrajectory.from_contact_maps(maps) + + def test_join(self): + segments = self.traj[0], self.traj[1:3], self.traj[3:] + assert [len(s) for s in segments] == [1, 2, 2] + assert md.join(segments) == self.traj + + cmaps = [ContactTrajectory(segment, cutoff=0.075, + n_neighbors_ignored=0) + for segment in segments] + + cmap = ContactTrajectory.join(cmaps) + + assert len(cmap) == len(self.map) + for i, (truth, beauty) in enumerate(zip(self.map, cmap)): + _contact_object_compare(truth, beauty) + assert truth == beauty + + _contact_object_compare(self.map, cmap) + assert self.map == cmap + + def test_rolling_frequency(self): + # smoke test; correctness is tested in tests for + # RollingContactFrequency + assert len(list(self.map.rolling_frequency(window_size=2))) == 4 + + +class TestMutableContactTrajectory(object): + def setup(self): + self.traj = md.load(find_testfile("trajectory.pdb")) + self.map = MutableContactTrajectory(self.traj, cutoff=0.075, + n_neighbors_ignored=0) + self.expected_atom_contacts = TRAJ_ATOM_CONTACTS.copy() + self.expected_residue_contacts = TRAJ_RES_CONTACTS.copy() + + def _test_expected_contacts(self, traj_map, exp_atoms, exp_res): + for cmap, exp_a, exp_r in zip(traj_map, exp_atoms, exp_res): + atom_counter = cmap.atom_contacts.counter + res_counter = cmap.residue_contacts.counter + assert atom_counter == counter_of_inner_list(exp_a) + assert res_counter == counter_of_inner_list(exp_r) + + def test_setitem(self): + cmap4 = ContactFrequency(self.traj[4], cutoff=0.075, + n_neighbors_ignored=0) + self.map[1] = cmap4 + expected_atoms = self.expected_atom_contacts + expected_atoms[1] = expected_atoms[4] + expected_res = self.expected_residue_contacts + expected_res[1] = expected_res[4] + self._test_expected_contacts(self.map, expected_atoms, expected_res) + + def test_delitem(self): + del self.map[1] + assert len(self.map) == 4 + expected_atoms = (self.expected_atom_contacts[:1] + + self.expected_atom_contacts[2:]) + expected_res = (self.expected_residue_contacts[:1] + + self.expected_residue_contacts[2:]) + self._test_expected_contacts(self.map, expected_atoms, expected_res) + + def test_insert(self): + cmap4 = self.map[4] + self.map.insert(0, cmap4) + expected_atoms = [TRAJ_ATOM_CONTACTS[4]] + TRAJ_ATOM_CONTACTS + expected_res = [TRAJ_RES_CONTACTS[4]] + TRAJ_RES_CONTACTS + self._test_expected_contacts(self.map, expected_atoms, expected_res) + + def test_hash_eq(self): + cmap = MutableContactTrajectory(self.traj, cutoff=0.075, + n_neighbors_ignored=0) + assert hash(cmap) != hash(self.map) + assert cmap != self.map + + +class TestWindowedIterator(object): + def setup(self): + self.iter = WindowedIterator(length=10, width=3, step=2, + slow_build=False) + + def test_startup_normal(self): + to_add, to_sub = self.iter._startup() + assert to_sub == slice(0, 0) + assert to_add == slice(0, 3) + assert self.iter.min == 0 + assert self.iter.max == 2 + + def test_startup_slow_build_step1(self): + itr = WindowedIterator(length=10, width=3, step=1, slow_build=True) + to_add, to_sub = itr._startup() + assert to_sub == slice(0, 0) + assert to_add == slice(0, 1) + assert itr.min == 0 + assert itr.max == 0 + + to_add, to_sub = itr._startup() + assert to_sub == slice(0, 0) + assert to_add == slice(1, 2) + assert itr.min == 0 + assert itr.max == 1 + + def test_normal(self): + self.iter.min = 0 + self.iter.max = 2 + to_add, to_sub = self.iter._normal() + assert to_sub == slice(0, 2) + assert to_add == slice(3, 5) + assert self.iter.min == 2 + assert self.iter.max == 4 + + @pytest.mark.parametrize('length,width,step,slow_build,expected', [ + (5, 3, 2, False, [(slice(0, 0), slice(0, 3), 0, 2), + (slice(0, 2), slice(3, 5), 2, 4)]), + (5, 3, 1, True, [(slice(0, 0), slice(0, 1), 0, 0), + (slice(0, 0), slice(1, 2), 0, 1), + (slice(0, 0), slice(2, 3), 0, 2), + (slice(0, 1), slice(3, 4), 1, 3), + (slice(1, 2), slice(4, 5), 2, 4)]), + (5, 3, 2, True, [(slice(0, 0), slice(0, 2), 0, 1), + (slice(0, 1), slice(2, 4), 1, 3)]), + (6, 3, 3, False, [(slice(0, 0), slice(0, 3), 0, 2), + (slice(0, 3), slice(3, 6), 3, 5)]), + (6, 3, 3, True, [(slice(0, 0), slice(0, 3), 0, 2), + (slice(0, 3), slice(3, 6), 3, 5)]), + ]) + def test_next(self, length, width, step, slow_build, expected): + itr = WindowedIterator(length, width, step, slow_build) + for expect in expected: + exp_sub, exp_add, exp_min, exp_max = expect + to_add, to_sub = next(itr) + assert to_add == exp_add + assert to_sub == exp_sub + assert itr.min == exp_min + assert itr.max == exp_max + with pytest.raises(StopIteration): + next(itr) + + +class TestRollingContactFrequency(object): + def setup(self): + self.traj = md.load(find_testfile("trajectory.pdb")) + self.map = ContactTrajectory(self.traj, cutoff=0.075, + n_neighbors_ignored=0) + self.rolling_freq = RollingContactFrequency(self.map, width=2, + step=1) + self.expected_atoms = [ + {frozenset([1, 4]): 0.5, frozenset([4, 6]): 1.0, + frozenset([5, 6]): 1.0, frozenset([1, 5]): 0.5}, + {frozenset([1, 5]): 0.5, frozenset([4, 6]): 1.0, + frozenset([5, 6]): 1.0, frozenset([1, 4]): 0.5}, + {frozenset([1, 4]): 1.0, frozenset([4, 6]): 1.0, + frozenset([5, 6]): 1.0, frozenset([4, 7]): 0.5, + frozenset([5, 7]): 0.5}, + {frozenset([0, 9]): 0.5, frozenset([0, 8]): 0.5, + frozenset([1, 8]): 0.5, frozenset([1, 9]): 0.5, + frozenset([1, 4]): 1.0, frozenset([8, 4]): 0.5, + frozenset([8, 5]): 0.5, frozenset([4, 6]): 1.0, + frozenset([4, 7]): 1.0, frozenset([5, 6]): 1.0, + frozenset([5, 7]): 1.0} + ] + self.expected_residues = [ + {frozenset([0, 2]): 1.0, frozenset([2, 3]): 1.0}, + {frozenset([0, 2]): 1.0, frozenset([2, 3]): 1.0}, + {frozenset([0, 2]): 1.0, frozenset([2, 3]): 1.0}, + {frozenset([0, 2]): 1.0, frozenset([2, 3]): 1.0, + frozenset([0, 4]): 0.5, frozenset([2, 4]): 0.5} + ] + + def test_normal_iteration(self): + results = list(freq for freq in self.rolling_freq) + assert len(results) == 4 + + atom_contacts = [r.atom_contacts.counter for r in results] + for beauty, truth in zip(atom_contacts, self.expected_atoms): + assert beauty == truth + + residue_contacts = [r.residue_contacts.counter for r in results] + for beauty, truth in zip(residue_contacts, self.expected_residues): + assert beauty == truth + + def test_slow_build_iteration(self): + self.rolling_freq.slow_build_iter = True + results = list(freq for freq in self.rolling_freq) + assert len(results) == 5 + + expected_atoms = [{frozenset([1, 4]): 1.0, frozenset([4, 6]): 1.0, + frozenset([5, 6]): 1.0}] + self.expected_atoms + expected_residues = ( + [{frozenset([0, 2]): 1.0, frozenset([2, 3]): 1.0}] + + self.expected_residues + ) + + atom_contacts = [r.atom_contacts.counter for r in results] + for beauty, truth in zip(atom_contacts, expected_atoms): + assert beauty == truth + + residue_contacts = [r.residue_contacts.counter for r in results] + for beauty, truth in zip(residue_contacts, expected_residues): + assert beauty == truth diff --git a/docs/api.rst b/docs/api.rst index 57e8502..7d0884a 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -17,6 +17,8 @@ Contact maps ContactMap ContactFrequency ContactDifference + ContactTrajectory + RollingContactFrequency Contact Concurrences -------------------- diff --git a/examples/contact_trajectory.ipynb b/examples/contact_trajectory.ipynb new file mode 100644 index 0000000..93592f4 --- /dev/null +++ b/examples/contact_trajectory.ipynb @@ -0,0 +1,342 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Contact Trajectories\n", + "\n", + "Sometimes you're interested in how contacts evolve in a trajectory, frame-by-frame. Contact Map Explorer provides the `ContactTrajectory` class for this purpose.\n", + "\n", + "We'll look at this using a trajectory of a specific inhibitor during its binding process to GSK3B. This system is also studied in the notebook on contact concurrences (with very similar initial discussion)." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n" + ] + } + ], + "source": [ + "from __future__ import print_function\n", + "%matplotlib inline\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "from contact_map import ContactTrajectory, RollingContactFrequency\n", + "import mdtraj as md\n", + "traj = md.load(\"data/gsk3b_example.h5\")\n", + "print(traj) # to see number of frames; size of system" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, we'll use MDTraj's [atom selection language](http://mdtraj.org/latest/atom_selection.html) to split out the protein and the ligand, which has residue name YYG in the input files. We're only interested in contacts between the protein and the ligand (not contacts within the protein). We'll also only look at heavy atom contacts." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "topology = traj.topology\n", + "yyg = topology.select('resname YYG and element != \"H\"')\n", + "protein = topology.select('protein and element != \"H\"')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Making an accessing a contact trajectory\n", + "\n", + "Contact trajectories have the same keyword arguments as other contact objects" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "contacts = ContactTrajectory(traj, query=yyg, haystack=protein)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Once the `ContactTrajectory` has been made, contacts for individual frames can be accessed either by taking the index of the `ContactTrajectory` itself, or by getting the list of contact (e.g., all the residue contacts frame-by-frame) and selecting the frame of interest." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[([YYG351, SER32], 1.0), ([YYG351, GLY31], 1.0), ([ASN30, YYG351], 1.0)]" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "contacts[0].residue_contacts.most_common()" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[([YYG351, SER32], 1.0), ([YYG351, GLY31], 1.0), ([ASN30, YYG351], 1.0)]" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "contacts.residue_contacts[0].most_common()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Advanced Python indexing is also allowed. In this example, note how the most common partners for YYG change! This is also what we see in the contact concurrences example." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[([VAL27, YYG351], 1.0), ([ARG107, YYG351], 1.0), ([ILE28, YYG351], 1.0)]\n", + "[([VAL27, YYG351], 1.0), ([ILE28, YYG351], 1.0), ([GLN151, YYG351], 1.0)]\n", + "[([VAL27, YYG351], 1.0), ([ASN30, YYG351], 1.0), ([GLY34, YYG351], 1.0)]\n", + "[([ASP166, YYG351], 1.0), ([PHE33, YYG351], 1.0), ([LYS149, YYG351], 1.0)]\n", + "[([YYG351, SER32], 1.0), ([VAL53, YYG351], 1.0), ([PHE33, YYG351], 1.0)]\n", + "[([GLU63, YYG351], 1.0), ([VAL53, YYG351], 1.0), ([PHE33, YYG351], 1.0)]\n", + "[([ASP166, YYG351], 1.0), ([VAL53, YYG351], 1.0), ([PHE33, YYG351], 1.0)]\n", + "[([YYG351, GLY168], 1.0), ([YYG351, SER32], 1.0), ([ASP166, YYG351], 1.0)]\n" + ] + } + ], + "source": [ + "for contact in contacts[50:80:4]:\n", + " print(contact.residue_contacts.most_common()[:3])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can easily turn the `ContactTrajectory` into `ContactFrequency`:" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEzCAYAAADHFzPDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAaY0lEQVR4nO3df7BdZX3v8feHgGDrLzQBY8BJ6g1VdGq0x8iMdor4g8B0mvqHNdpRaulN7UBbO/Zeoe0oHctU6697O6jciAj2tk1zW6zRiaXIlHqdyo/gRSRBaq5QOCSSRG0r/gATvvePvXLZHnbO3vtkr5yzDu/XzJ6z1trPevbzcDIfnvOsZ62dqkKS1E3HzHcDJElzZ4hLUocZ4pLUYYa4JHWYIS5JHWaIS1KHGeKSdJQkuTLJ3iR3HOb9JPmzJLuS3J7kxcPqNMQl6ei5Clg3y/vnAKub10bgo8MqNMQl6Sipqi8A356lyHrgk9VzI/C0JMtnq9MQl6SFYwVwX9/+dHPssI5ttTkjWrp0aa1cuXK+myFpkbr11lv3V9WyUcv/p6S+P4fP2QM7gB/2HdpUVZvGqCIDjs36bJQFEeIrV65k+803z3czJC1SWbLkX8cp/33gN+bwOZfAD6tqag6nHjINnNq3fwqwe7YTnE6RpIVjK/DmZpXKGcC/V9We2U5YECNxSVpIwuB5jSOuN/kr4ExgaZJp4F3AcQBVdTmwDTgX2EXvD4K3DKtzaIgnOQH4AnB8U/5vqupdSS4B/jOwryn6+1W1rTnnYuB84CDw21V17ci9lKQFoI0Qr6o3DHm/gAvGqXOUkfhDwFlV9WCS44AvJvlc896Hqur9/YWTnA5sAJ4PPAv4fJLTqurgOA2TpPnURoi3YeiceLNe8cFm97jmNdvV0vXA5qp6qKrupvdnwdojbqkkHUWZw2s+jHRhM8mSJLcBe4Hrquqm5q0Lm1tDr0xyYnNs7HWOkrTQLKoQr6qDVbWG3nKXtUleQO920OcAa4A9wAea4iOtc0yyMcn2JNv37ds34BRJmh9zCfAFHeKHVNW/ATcA66rqgSbcHwE+xqNTJiOtc6yqTVU1VVVTy5aNvAZfko6KRRPiSZYleVqz/UTgVcDXZtzP/1rg0FO5tgIbkhyfZBW9B7l4J4+kTulKiI+yOmU5cHWSJfRCf0tVfTbJnydZQ2+q5B6aG5yqakeSLcBO4ABwgStTJHVNV1anDA3xqrodeNGA42+a5ZxLgUuPrGmSNH8WTYhL0uPNfE6PjMsQl6QBDHFJ6jBDXJI6zBCXpA7rSoj7PHFJ6jBH4pI0g6tTJKnjDHFJ6jBDXJI6zBCXpA4zxCWpo7ywKUkdZ4hLUocZ4pLUYYa4JHWYIS5JHeWFTUnqOENckjrMEJekDjPEJanDDHFJ6qguXdj0SyEkqcMciUvSAF0ZiRvikjSAIS5JHWaIS1JHdenCpiEuSQN0JcSHrk5JckKSm5N8JcmOJH/UHH96kuuSfL35eWLfORcn2ZXkriRnt9kBSWpD5vCaD6MsMXwIOKuqXgisAdYlOQO4CLi+qlYD1zf7JDkd2AA8H1gHfCTJkjYaL0ltWTQhXj0PNrvHNa8C1gNXN8evBn6p2V4PbK6qh6rqbmAXsHairZakli2aEAdIsiTJbcBe4Lqqugk4uar2ADQ/T2qKrwDu6zt9ujkmSZ0wlwBf0CFeVQerag1wCrA2yQtmKT6oL/WYQsnGJNuTbN+3b99orZWko2RRhfghVfVvwA305rofSLIcoPm5tyk2DZzad9opwO4BdW2qqqmqmlq2bNkcmi5J7Vk0IZ5kWZKnNdtPBF4FfA3YCpzXFDsP+HSzvRXYkOT4JKuA1cDNk264JLWpKyE+yjrx5cDVzQqTY4AtVfXZJF8CtiQ5H7gXeB1AVe1IsgXYCRwALqiqg+00X5La0ZV14kNDvKpuB1404Pi3gFce5pxLgUuPuHWSNA+8Y1OSOq4rIe7zxCVpgDbmxJOsa+5k35XkogHvPzXJZ/rukH/LsDoNcUk6Cprrih8GzgFOB97Q3OHe7wJgZ3OH/JnAB5I8YbZ6DXFJGqCFkfhaYFdVfaOqHgY207vDvV8BT04S4EnAt+ktEDks58QlaYA5zokvTbK9b39TVW1qtgfdzf7SGedfRm+Z9m7gycDrq+qR2T7QEJekGY5gdcr+qpqapdqZZt7NfjZwG3AW8BzguiT/u6r+43Af6HSKJA3QwnTKKHezvwW4pnnw4C7gbuC5s1VqiEvSAC2E+C3A6iSrmouVG+hNnfS7l+b+myQnAz8NfGO2Sp1OkaQBJr1OvKoOJLkQuBZYAlzZ3OH+1ub9y4F3A1cl+WrThHdU1f7Z6jXEJWmANm72qaptwLYZxy7v294NvGacOg1xSZrB2+4lqeMMcUnqMENckjrMEJekDjPEJamjvLApSR1niEtShxniktRhXQlxn50iSR3mSFySBujKSNwQl6QZXJ0iSR1niEtShxniktRhhrgkdZghLkkd5YVNSeo4Q1ySOqwrIT70js0kpyb5xyR3JtmR5Hea45ckuT/Jbc3r3L5zLk6yK8ldSc5uswOS1IYWvu2+FaOMxA8Ab6+qLyd5MnBrkuua9z5UVe/vL5zkdGAD8HzgWcDnk5xWVQcn2XBJatOiGYlX1Z6q+nKz/V3gTmDFLKesBzZX1UNVdTewC1g7icZK0tEwl1H4fIX+WA/ASrISeBFwU3PowiS3J7kyyYnNsRXAfX2nTTN76EvSgrPoQjzJk4C/Bd5WVf8BfBR4DrAG2AN84FDRAafXgPo2JtmeZPu+ffvGbrgktWlRhXiS4+gF+F9U1TUAVfVAVR2sqkeAj/HolMk0cGrf6acAu2fWWVWbqmqqqqaWLVt2JH2QpIlbNCGeJMDHgTur6oN9x5f3FXstcEezvRXYkOT4JKuA1cDNk2uyJLWvKyE+yuqUlwFvAr6a5Lbm2O8Db0iyht5UyT3AbwBU1Y4kW4Cd9Fa2XODKFEldsqju2KyqLzK4P9tmOedS4NIjaJckaQTesSlJAyyakbgkPR4Z4pLUYYa4JHWYIS5JHbWoVqdI0uORIS5JHWaIS1KHGeKS1GGGuCR1lBc2JanjDHFJ6jBDXJI6zBCXpA4zxCWpo7ywKUkd15UQH+vb7iVJC4shLkkDtPEdm0nWJbkrya4kFx2mzJlJbkuyI8k/DavT6RRJGmDS0ylJlgAfBl4NTAO3JNlaVTv7yjwN+AiwrqruTXLSsHodiUvSAC2MxNcCu6rqG1X1MLAZWD+jzBuBa6rqXoCq2jusUkNckmaYS4CPEOIrgPv69qebY/1OA05MckOSW5O8eVilTqdI0gBznE5ZmmR73/6mqto0S5U1Y/9Y4GeBVwJPBL6U5Maq+pfDfaAhLkkDzDHE91fV1GHemwZO7ds/Bdg9oMz+qvoe8L0kXwBeCBw2xJ1OkaQBWphOuQVYnWRVkicAG4CtM8p8Gvi5JMcm+QngpcCds1XqSFySBpj06pSqOpDkQuBaYAlwZVXtSPLW5v3Lq+rOJH8P3A48AlxRVXfMVq8hLkkztHXbfVVtA7bNOHb5jP33Ae8btU5DXJIG6Mpt94a4JA1giEtSh3UlxIeuTklyapJ/THJncy//7zTHn57kuiRfb36e2HfOxc2zAe5KcnabHZCkNrTx7JQ2jLLE8ADw9qp6HnAGcEGS04GLgOurajVwfbNP894G4PnAOuAjzTMDJKkTWrpjsxVDQ7yq9lTVl5vt79Jbs7iC3j3/VzfFrgZ+qdleD2yuqoeq6m5gF71nBkhSZyyaEO+XZCXwIuAm4OSq2gO9oAcOPW1rlOcDSNKCtuhCPMmTgL8F3lZV/zFb0QHHZj4fgCQbk2xPsn3fvn2jNkOS1GekEE9yHL0A/4uquqY5/ECS5c37y4FDj0wc5fkAVNWmqpqqqqlly5bNtf2S1IpFMxJPEuDjwJ1V9cG+t7YC5zXb59G75//Q8Q1Jjk+yClgN3Dy5JktSu7p0YXOUdeIvA94EfDXJbc2x3wfeA2xJcj5wL/A6gOZZAFuAnfRWtlxQVQcn3nJJalFX1okPDfGq+iKH788rD3POpcClR9AuSZpXiybEJenxyBCXpA4zxCWpo+bzQuW4DHFJGsAQl6QOM8QlqcMMcUnqMENckjrKC5vj2r8frrpqvlshSf+fIT6OpUvhV391vlshabE6//yxTzHEJanDes/+G1M95qnbrTPEJWmQjoT4WN/sI0laWByJS9Igx8xhjPvII5NvxxCGuCTNlMwtxOeBIS5JgxjiktRRjsQlqeMM8TF4x6akhcYQH0MVHDgw362QpB6nU8b0pCfBGWfMdysk6VGG+BhOOAGe+9z5boUk9TgSH1MCxy6MpkgSYIhLUqcZ4pLUUU6nSFLHGeKS1FGOxCWp4wxxSeqwjoT40FYmuTLJ3iR39B27JMn9SW5rXuf2vXdxkl1J7kpydlsNlySNNhK/CrgM+OSM4x+qqvf3H0hyOrABeD7wLODzSU6rqoMTaKskHR2LaU68qr6QZOWI9a0HNlfVQ8DdSXYBa4EvzXrWj34E3/zmiB8hSUfBYgnxWVyY5M3AduDtVfUdYAVwY1+Z6ebY7O6/H/7wD4+gKZI0QYtpJH4YHwXeDVTz8wPArwGDvh564Nc/J9kIbAR49ooV8Md/PMemSNIQn/jE+Oe0EOJJ1gH/HVgCXFFV7zlMuZfQGxC/vqr+ZrY65xTiVfVA34d9DPhsszsNnNpX9BRg92Hq2ARsApiamiqe+cy5NEWS2jHhEE+yBPgw8Gp6WXlLkq1VtXNAufcC145S75xCPMnyqtrT7L4WOLRyZSvwl0k+SO/C5mrg5qEV/uAHcMcdQ4tJ0lHRznTKWmBXVX2j9xHZTO864s4Z5X4L+FvgJaNUOjTEk/wVcCawNMk08C7gzCRr6E2V3AP8BkBV7UiypWnUAeCCkVamPPGJ8IIXjNJeSTo6Jh/iK4D7+vangZf2F0iygt7A+CwmFeJV9YYBhz8+S/lLgUtH+XBJWpDmPhJfmmR73/6mZuoYRrtm+N+Ad1TVwWRQ8cfyjk1JGmRuIb6/qqYO894o1wyngM1NgC8Fzk1yoKr+7nAfaIhL0iCTn065BVidZBVwP70bI9/YX6CqVh3aTnIV8NnZAhwMcUl6rBYubFbVgSQX0lt1sgS4srmO+Nbm/cvnUq8hLkmDtLBOvKq2AdtmHBsY3lX1q6PUaYhL0kyPgzs2JWlxM8QlqcMMcUnqKKdTJKnjOhLi3WilJGkgR+KSNJPTKZLUcYa4JHWYIS5JHeV0iiR1nCEuSR3lSFySOs4Ql6SOciQuSR1niEtShxniktRRTqdIUscZ4pLUUY7EJanjDHFJ6jBDfAxVcODAfLdCknqcThlTAscujKZIUpcsjOTcuROmpua7FZL0qMUyEk9yJfALwN6qekFz7OnAXwMrgXuAX66q7zTvXQycDxwEfruqrh3aiuc9D770pbn1QJKGOf748covsumUq4DLgE/2HbsIuL6q3pPkomb/HUlOBzYAzweeBXw+yWlVdXDWT9i/H664Yg7Nl6SWLJYQr6ovJFk54/B64Mxm+2rgBuAdzfHNVfUQcHeSXcBaYPZh9sMPw/T0GM2WpJYtlhA/jJOrag9AVe1JclJzfAVwY1+56ebY7H7wA7jjjjk2RZImbJFNp4wjA47VwILJRmAjwLOf+lRYs2bCTZGkxmc+M/45izzEH0iyvBmFLwf2NsengVP7yp0C7B5UQVVtAjYBTE1NFZdcMsemSNIQ7373eOU7NBKfayu3Auc12+cBn+47viHJ8UlWAauBm4+siZI0D445ZvzXPBhlieFf0buIuTTJNPAu4D3AliTnA/cCrwOoqh1JtgA7gQPABUNXpkjSQtSRkfgoq1PecJi3XnmY8pcClx5JoyRpXnVoOmVh3LEpSQuNIS5JHeVIXJI6zhCXpA7LoNteFp6FEeKPPAI//OF8t0KSepxOGZPPE5e00BjiYzDEJWlOTE5JGsSRuCR1lHPiktRxhrgkdZQjcUnqOENckjqsIyHejVZK0tF0aDplws8TT7IuyV1JdjVfMj/z/V9Jcnvz+uckLxxWpyNxSRpkwiPxJEuADwOvpvctaLck2VpVO/uK3Q38fFV9J8k59L797KWz1WuIS9JM7VzYXAvsqqpv9D4im4H19L5EB4Cq+ue+8jfS+4rLWRnikjTI3EJ8aZLtffubmu8TBlgB3Nf33jSzj7LPBz437AMNcUkaZG4hvr+qpg7z3qDHItbAgskr6IX4y4d9oCEuSTO1M50yDZzat38KsPuxH52fAa4Azqmqbw2r1BCXpEEmH+K3AKuTrALuBzYAb+wvkOTZwDXAm6rqX0ap1BCXpJlaGIlX1YEkFwLXAkuAK6tqR5K3Nu9fDrwTeAbwkfS+lOLALNMzgCEuSYO1cLNPVW0Dts04dnnf9q8Dvz5OnYa4JM3ks1PG9KMfwTe/Od+tkKTOWRghftxx8MxnzncrJOlRjsTH4BclS1poDPExHHMMnHDCfLdCknqcE5ekjjPEJamjHi8j8ST3AN8FDtIsSk/ydOCvgZXAPcAvV9V3jqyZknSUPR5CvPGKqtrft38RcH1Vvad56PlFwDtmreHhh2F6egJNkaQJeRyF+EzrgTOb7auBGxgW4t//PmzfPmsRSTpqHi/TKfQeo/gPSQr4H81zc0+uqj0AVbUnyUlDa3nKU2DduiNsiiRN0OMkxF9WVbuboL4uyddGPTHJRmAjwLOf/WyXGEpaODo0Ej+iVlbV7ubnXuBT9L5+6IEkywGan3sPc+6mqpqqqqlly5YdSTMkafJa+KLkVpo51xOT/GSSJx/aBl4D3AFsBc5rip0HfPpIGylJR11HQvxIplNOBj7VPPP2WOAvq+rvk9wCbElyPnAv8Lojb6YkHUUdmk6Zc4g339j8wgHHvwW88kgaJUnzbrGHuCQtWo+HkbgkLWodCfFutFKSNJAjcUkapCMj8YUR4g8+CF/84ny3QpJ6nBMf0wknwHOfO9+tkKRHGeJjOHiwNxqXpIXAkfiYHnnEEJe0sBjiozuwcyf7XviY+4Ykaf4Y4qM79qd+imV/8ifz3QxJi9XrXz9eeadTxrR3L1x22Xy3QpIeZYiP4XvfgxtvnO9WSFKPI/ExHXssLF06362QtFjt2TP+OYb4GPx6Nklt+sQnxj/HEB/DypVwxRXz3QpJi9W4Ie50iiR1nCEuSR3lSFySOs4Ql6QOM8QlqaM6NJ3SjVZKkgZyJC5Jg3RkJG6IS9JMHZpOMcQlaRBDXJI6zBCXpI5yOkWSOs4Ql6SOciQuSR3XkRBvrZVJ1iW5K8muJBe19TmSNHGHRuLjvoZWO3supufPmvdvT/LiYXW2EuJJlgAfBs4BTgfekOT0Nj5Lklox4RAfMRfPAVY3r43AR4c2cy59G8FaYFdVfaOqHgY2A+tb+ixJmrzJj8RHycX1wCer50bgaUmWz9rMufZviBXAfX37080xSVr42plOGSUXx87Oti5sZsCx+rECyUZ6fy4APJQlS+5oqS0LwVJg/3w3omWLvY/2r9t+epzCt95667VZsmQu395+QpLtffubqmpTsz00F0cs82PaCvFp4NS+/VOA3f0Fmo5tAkiyvaqmWmrLvFvs/YPF30f7120zgnWoqmrjm9uH5uKIZX5MW9MptwCrk6xK8gRgA7C1pc+SpC4YJRe3Am9uVqmcAfx7Ve2ZrdJWRuJVdSDJhcC1wBLgyqra0cZnSVIXHC4Xk7y1ef9yYBtwLrAL+D7wlmH1tnazT1Vtaxo0ik3Di3TaYu8fLP4+2r9uWxD9G5SLTXgf2i7ggnHqTO8cSVIXdeO+UknSQBMP8SQnJLk5yVeS7EjyR83xS5Lcn+S25nVu3zkXN7eZ3pXk7MPU+/Qk1yX5evPzxEm3fRQt9u99Sb7W3Gr7qSRPO1p9mtGOVvrXV/b3klSSuSzfOmJt9i/JbzVldiT506PRnwFtaOvf55okNzbnbk+y9mj1aUY7xupfkmck+cckDya5bJZ6F0S+zElVTfRFb53jk5rt44CbgDOAS4DfG1D+dOArwPHAKuD/AksGlPtT4KJm+yLgvZNu+zz37zXAsc32exdb/5qyp9K7qPOvwNLF1D/gFcDngeOb/ZMWWf/+ATin2T4XuKEj/ftJ4OXAW4HLZql3QeTLXF4TH4lXz4PN7nHNa7aJ9/XA5qp6qKrupndVdtD/5dcDVzfbVwO/NKEmj6Wt/lXVP1TVgWb3RnrrQ4+6Fn9/AB8C/uuQ+lrVYv9+E3hPVT3UfM7eCTZ7ZC32r4CnNNtPZcja5baM27+q+l5VfRH44ZCqF0S+zEVrD8BKchuwF7iuqm5q3rqwmS64su/PlVFvMz25mvWSzc+T2mj7KFrqX79fAz430UaPoY3+JflF4P6q+kqbbR9FS7+/04CfS3JTkn9K8pLWOjBES/17G/C+JPcB7wcubqn5Q43Zv1EtmHwZVyshXlUHq2oNvdHk2iQvoPc0rucAa4A9wAea4mPfZjrf2uxfkj8ADgB/MdFGj2HS/UvyE8AfAO9srdFjaOn3dyxwIr0/7f8LsCXJoHNb11L/fhP43ao6Ffhd4OMTb/iIxuzfotfq6pSq+jfgBmBdVT3Q/Md/BPgYj/7JNuptpg+keZpX83Ne/lztN+H+keQ84BeAX6lmcm4+TbB/z6E33/qVJPc0Zb6c5JktNn+oCf/+poFrmj/3bwYeofc8knkz4f6dB1zTbP8vDj9ldtSM2L9RLbh8GVUbq1OWHVpZkeSJwKuAr+XHH6f4WuDQA6+2AhuSHJ9kFb3n6N48oOqt9P4h0fz89KTbPoq2+pdkHfAO4Ber6vtt9mE2bfSvqr5aVSdV1cqqWkkvOF5cVd9suTuP0eK/z78DzmrqPQ14AvPwQKkW+7cb+Plm+yzg6220f5g59G9UCyJf5mTSV0qBnwH+D3A7vf+Q72yO/znw1eb4VmB53zl/QO+q+F00V8Cb41cAU832M4Dr6f3juR54+qTbPs/920VvbvK25nX5YurfjM+4h/lbndLW7+8JwP9s6vwycNYi69/LgVvprWS5CfjZDvXvHuDbwIP0BhCnD+jfgsiXuby8Y1OSOsw7NiWpwwxxSeowQ1ySOswQl6QOM8QlqcMMcUnqMENckjrMEJekDvt/B6egDFZk3WQAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "freq = contacts.contact_frequency()\n", + "\n", + "fig, ax = plt.subplots(figsize=(5.5,5))\n", + "freq.residue_contacts.plot_axes(ax=ax)\n", + "ax.set_xlim(*contacts.query_residue_range);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Rolling Contact Frequencies\n", + "\n", + "A `ContactTrajectory` keeps all the time-dependent information about the contacts, whereas a `ContactFrequency`, as plotted above, loses all of it. What about something in between? For this, we have a `RollingContactFrequency`, which acts like a rolling average. It creates a contact frequency over a certain window of frames, with a certain step size between each window.\n", + "\n", + "This can be created either with the `RollingContactFrequency` object, or, more easily, with the `ContactTrajectory.rolling_frequency()` method." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "RollingContactFrequency(contacts, width=30, step=14)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "rolling_frequencies = contacts.rolling_frequency(window_size=30, step=14)\n", + "rolling_frequencies" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we'll plot each windowed frequency, and we will see the transition as some contacts fade out and others grow in." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtcAAAJDCAYAAADAeTgIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdf6xldX3v/+ebM4B6/cHgDHSE8c7EDm2pqSjT0aQ2tf6og+ntaFLraFOpF4M00Ns28V6gJtamIcHa9t5vI8odKQGvjXPHllunDS1FUmqaCsMMBWRAdAoUjjPlR/1VJcIdfN8/9jpfl+esc/Y+Z6119ufs/XwkO2fvtdfe+/2ZObx4z+ez1tqRmUiSJElq74RxFyBJkiRNCptrSZIkqSM215IkSVJHbK4lSZKkjthcS5IkSR2xuZYkSZI6YnMtSR2JiGsj4vGIuHeR5yMi/jgijkTEPRHxqtWuUZI00Fdm21xLUneuA3Yu8fx5wLbqdiHw8VWoSZLU7Dp6yGyba0nqSGZ+HvjaErvsAj6ZA7cBp0TEptWpTpJU11dm21xL0uo5A3i09ni22iZJKs+KMntdb+Usw4YNG3LLli3jLkPShDl06NCTmblx1P1/OCKfWuL5Y3AY+G5t057M3LOMkqJhWy7j9cUwtyV1bVIyu4jmesuWLRw8cGDcZUiaMDEz8y/L2f8p4H1LPP8h+G5mbm9R0iywufb4TOBoi/cbG3NbUtcmJbM9LESSKjHk1oH9wLurM9BfA3wzM49189aSNF1KzewiZq4lqRRtZhwi4tPA64ANETEL/A5wIkBmXg3cCLwFOMJg0uU9rYqVpClXYmbbXEtSTZugzsx3Dnk+gYtbfIQkqabEzLa5lqRKh0uJkqSelZrZNteSVOOJKJK0dpSY2TbXklRT4iyIJKlZiZltcy1JlaDMWRBJ0kKlZrbNtSTVlBjUkqRmJWa2zbUk1ZS4xChJalZiZttcS1Kl1CVGSdJCpWa2zbUk1ZQ4CyJJalZiZttcS1JNibMgkqRmJWa2zbUkVUpdYpQkLVRqZttcS1JNiUuMkqRmJWa2zbUk1ZQ4CyJJalZiZttcS1JNibMgkqRmJWa2zbUkVUo9fk+StFCpmW1zLUk1JQa1JKlZiZltcy1JNSUuMUqSmpWY2TbXklQpdYlRkrRQqZltcy1JNSXOgkiSmpWY2UMb/oh4TkQciIi7I+JwRPxutf3UiLg5Ir5S/Vxfe83lEXEkIh6IiDf3OQBJ6srcLMhit7XAzJY0LUrN7FE++2ng9Zn5CuAcYGdEvAa4DLglM7cBt1SPiYizgd3AjwM7gY9FxEwfxUtS10oM6mUysyVNjRIze+hn58C3q4cnVrcEdgHXV9uvB95a3d8F7M3MpzPzIeAIsKPTqiWpJ7HEbS0wsyVNkxIze6TGPiJmIuIu4HHg5sy8HTg9M48BVD9Pq3Y/A3i09vLZapskFa3UJcblMrMlTYNSM3ukz87MZzPzHOBMYEdEvHyJ3Zv+sZALdoq4MCIORsTBJ554YrRqJalnJc6CLFcfmQ3mtqTylJjZy2rsM/MbwK0Mjst7LCI2AVQ/H692mwU21152JnC04b32ZOb2zNy+cePGFZQuSd1rMwsSETurkwKPRMRlDc+/KCL+snay4Xs6LX6eLjO7ej9zW1JR2s5c95Hbo1wtZGNEnFLdfy7wRuBLwH7g/Gq384HPVvf3A7sj4uSI2ApsAw6MMkBJGqc2S4zVSYBXAecBZwPvrE4WrLsYuK862fB1wB9GxEmdDQAzW9L0aHtYSF+5Pcp1rjcB11cFnADsy8y/iogvAPsi4gLgEeDtAJl5OCL2AfcBx4GLM/PZET5HksauxVLiDuBIZj4IEBF7GZwseF9tnwReEBEBPB/4GoOc7JKZLWlqtDz8o5fcHtpcZ+Y9wCsbtv8b8IZFXnMFcMWw95ak0rQ4CabpxMBXz9vnowxmio8CLwDekZnfW/lHLmRmS5omLU9c7CW319IJ8JLUuyEnx2yYO6Gvul0476XzzT8x8M3AXcBLGFyD+qMR8cJOByBJU6RFZkNPue3Xn0tSZe74vSU8mZnbF3lulBMD3wNcmZkJHImIh4AfxWOcJWnZWmY29JTbzlxLUk2Lk2PuALZFxNbqZJfdDJYS6x6hOjQjIk4HfgR4sKPSJWnqtLxaSC+57cy1JNWs9OSYzDweEZcANwEzwLXVyYIXVc9fDfwecF1EfLH6qEsz88ku6pakadTmhMa+ctvmWpIqIywxLikzbwRunLft6tr9o8DPtfgISVKlbWZDP7ltcy1JNWvpmxgladqVmNk215JU44kokrR2lJjZNteSVOliiVGStDpKzWyba0mqKXGJUZLUrMTMtrmWpJoSZ0EkSc1KzGyba0mq1L7VS5JUuFIz2+ZakmpKnAWRJDUrMbNtriWppsSgliQ1KzGzba4lqVLqEqMkaaFSM9vmWpJqSpwFkSQ1KzGzba4lqabEWRBJUrMSM9vmWpIqpX4hgSRpoVIz2+ZakmpKDGpJUrMSM9vmWpJqSlxilCQ1KzGzba4lqVLqEqMkaaFSM9vmWpJqSpwFkSQ1KzGzba4lqabEWRBJUrMSM9vmWpIqpS4xSpIWKjWzba4lqabEJUZJUrMSM9vmWpJqSpwFkSQ1KzGzba4lqabEWRBJUrMSM9vmWpIqpR6/J0laqNTMtrmWpJoSg1qS1KzEzC6xJkkam1jiNvS1ETsj4oGIOBIRly2yz+si4q6IOBwRf99Z4ZI0hdpkNvST285cS1KlzRJjRMwAVwFvAmaBOyJif2beV9vnFOBjwM7MfCQiTmtbsyRNq7aHhfSV285cS1JNi1mQHcCRzHwwM58B9gK75u3zLuCGzHwEIDMf76xwSZpCLWeue8ltm2tJqjlhidsQZwCP1h7PVtvqzgLWR8StEXEoIt7dQcmSNLVaZDb0lNseFiJJlRGWGDdExMHa4z2Zuaf28vly3uN1wLnAG4DnAl+IiNsy88srKliSpljLzJ57i/la57bNtSTVDFlKfDIzty/y3Cywufb4TOBowz5PZuZ3gO9ExOeBVwA215K0Ai0yG3rKbQ8LkaSaFkuMdwDbImJrRJwE7Ab2z9vns8BPR8S6iHge8Grg/s6Kl6Qp0/KwkF5ye+hnR8TmiPi7iLi/ugTJb1TbT42ImyPiK9XP9bXXXF5d0uSBiHjzaOOTpPFa6sSYYSfHZOZx4BLgJgbBuy8zD0fERRFxUbXP/cDfAPcAB4BrMvPeTsdgZkuaEm0yG/rL7cicf2jJvB0iNgGbMvPOiHgBcAh4K/CrwNcy88rquoDrM/PSiDgb+DSDMzBfAnwOOCszn13sM7Zv354HDxwY+ocgScsRMzOHhiwJ/oBzIvJzSzy/EZb1fuOwGpkN5rak7k1KZg+duc7MY5l5Z3X/3xl09mcwuFTJ9dVu1zMIb6rtezPz6cx8CDjCILQlqXgtlxjHzsyWNE1KzOxlfXZEbAFeCdwOnJ6Zx2AQ5sDcRbVHuayJJBWn7RJjacxsSZOs1MweubmOiOcDfw78ZmZ+a6ldG7YtOPYkIi6MiIMRcfCJJ54YtQxJ6lWJsyAr0XVmV+9pbksqSomZPdJnR8SJDEL6TzPzhmrzY9WxfXPH+M19Y80olzUhM/dk5vbM3L5x48aV1i9JnSpxFmS5+shsMLcllafEzB7laiEB/Alwf2b+Ue2p/cD51f3zGVyqZG777og4OSK2AtsYnF0pSUWb+0KC0mZBlsPMljQtSs3sUb5E5qeAXwG+GBF3Vdt+G7gS2BcRFwCPAG8HqC5hsg+4DzgOXDzsrHNJKsVaaaKXYGZLmholZvbQ5joz/4HFZ9ffsMhrrgCuaFGXJI3FWjr8o4mZLWmalJjZfv25JFXmlhglSeUrNbNtriWppsRZEElSsxIz2+ZakmpKnAWRJDUrMbNtriWpUuoSoyRpoVIz2+ZakmpKXGKUJDUrMbNtriWpJk5YYh7ke99bvUIkSUOVmNk215JUV2BQS5IWUWBm21xL0pwIWLdELB4/vnq1SJKWVmhm21xL0pxhQS1JKkehmV1eRZI0TkstMUqSylJgZttcS9KcQmdBJEkNCs3s8iqSpHGJKHIWRJLUoNDMtrmWpLoCZ0EkSYsoMLPLq0iSxqXQJUZJUoNCM7u8iiRpXApdYpQkNSg0s22uJamuwFkQSdIiCszs8tp9SRqXuVmQxW5DXx47I+KBiDgSEZctsd9PRsSzEfGLndYvSdOkZWYP3qL73C6v3ZekcWlx/F5EzABXAW8CZoE7ImJ/Zt7XsN+HgZtaVitJ063lMdd95bYz15JUt27d4rel7QCOZOaDmfkMsBfY1bDfrwN/DjzebeGSNIVWntnQU27bXEvSnHZLjGcAj9Yez1bbam8fZwBvA67utG5JmkbtDwvpJbc9LESS5gxfYtwQEQdrj/dk5p65Vzfsn/Me/w/g0sx8NqJpd0nSyNplNvSU2zbXklS39GzHk5m5fZHnZoHNtcdnAkfn7bMd2FsF9AbgLRFxPDP/YoXVStJ0W3lmQ0+5bXMtSXPanRxzB7AtIrYCXwV2A++q75CZW7//UXEd8Fc21pK0Qu2/RKaX3La5lqQ5LYI6M49HxCUMziafAa7NzMMRcVH1vMdZS1KXWjbXfeW2zbUk1bX4tq/MvBG4cd62xnDOzF9d8QdJkgZafkNjH7ltcy1Jc9ovMUqSVkuhmV1eRZI0LoUGtSSpQaGZXV5FkjROLZcYJUmrqMDMtrmWpDmFzoJIkhoUmtnlVSRJ41TgLIgkaREFZrbNtSTNKXQWRJLUoNDMLq8iSRqXQoNaktSg0MwuryJJGqcClxglSYsoMLNtriVpTqGzIJKkBoVmdnkVSdK4RBQ5CyJJalBoZpfRXB8/Dt/4xrirkKQiZ0EkSYsoMLOHVhQR1wI/DzyemS+vtp0K/G9gC/Aw8EuZ+fXqucuBC4Bngf+SmTcNr2IdnHLKykYgSV0pdIlxuVYltyVp3ArN7FHm0q8Dds7bdhlwS2ZuA26pHhMRZwO7gR+vXvOxiJjprFpJ6tPcEuNit7XjOsxtSZOu0Mwe2u5n5ucjYsu8zbuA11X3rwduBS6ttu/NzKeBhyLiCLAD+MKSHzI7C5ddtoyyJaknBc6CLNdq5PZThw5xaMYeXNKYFZjZK63o9Mw8BpCZxyLitGr7GcBttf1mq21Le+YZePjhFZYiSR0p9OSYjnSa2ycCL+m8RElahkIzu+t2Pxq2ZeOOERcCFwK89KUvhb17Oy5F0tT7zGeWt3+hx+/1bMW5vemhh/qsS9K0We5qWKGZvdJ2/7GI2ARQ/Xy82j4LbK7tdyZwtOkNMnNPZm7PzO0bN25cYRmS1LF16xa/rW3mtqTJU2Bmr7S53g+cX90/H/hsbfvuiDg5IrYC24AD7UqUpFVS6MkxHTG3JU2WQjN7lEvxfZrBSTAbImIW+B3gSmBfRFwAPAK8HSAzD0fEPuA+4DhwcWY+21PtktStQpcYl8vcljQVCs3sUa4W8s5FnnrDIvtfAVzRpihJGpu1P0NtbkuaHgVmdhntvt/QKKkEEcs/oUaSNB6FZnYZzbXf0CipFAUuMUqSFlFgZpdXkSSNS6HXTJUkNSg0s8urSJLGZe7kmBVe1ikidkbEAxFxJCIWfO1sRPxyRNxT3f4xIl7RyzgkaRq0zOzBW3Sf285cS1LdCmdBImIGuAp4E4NrR98REfsz877abg8BP5OZX4+I84A9wKtbVixJ06vFzHVfuW1zLUlz2l3WaQdwJDMfHLxV7AV2MbjEHQCZ+Y+1/W9j8IUtkqSVaH8pvl5y2+Zakua0C+ozgEdrj2dZenbjAuCvV/phkjT12jfXveS2zbUk1S29xLghIg7WHu/JzD3V/WjYP5veJCJ+lkFIv3ZFNUqSBlae2dBTbttcS9Kc4bMgT2bm9kWemwU21x6fCRxd+BHxE8A1wHmZ+W8rLVWSpl67zIaecturhUhS3QknLH5b2h3AtojYGhEnAbuB/fUdIuKlwA3Ar2Tml3upX5KmycozG3rKbWeuJWlOi+P3MvN4RFwC3ATMANdm5uGIuKh6/mrgg8CLgY9FBMDxIbMqkqTFtDzmuq/ctrmWpDntg/pG4MZ5266u3X8v8N4Vf4Ak6fvan9DYS27bXEtSXYHf9iVJWkSBmW1zLUlzOpgFkSStkkIzu7yKJGlcIoqcBZEkNSg0s22uJamuwFkQSdIiCszs8iqSpHEpdIlRktSg0MwuryJJGpdClxglSQ0KzWyba0mqK3AWRJK0iAIzu7yKJGlcCp0FkSQ1KDSzba4laU6hx+9JkhoUmtnlVSRJ41RgUEuSFlFgZpdR0ewsXHbZuKuQNO0KXWIs0te/Dn/2Z+OuQtI0KzSzy2iun3kGHn543FVImnaFLjGW6OkHH+Sf3/GOcZchaZoVmtllVDQzA6ecMu4qJKnIWZASPQt8Y9xFSFKBmV1Gc/3CF8Ib3zjuKiRNmk98Ynn7FzoLUqLnnXsu5x44MO4yJE2SmZnl7V9oZhdR0bcefJCbXV6UNG6FBnWJ/vXQIT683P8RSlKXCs3sIio6Edg07iIkCYpcYizR84BXjbsISSows4torp977rm83OVFSV2bkCXGEr3w3HN5k7ktqUsTktnlVSRJ41TgLIgkaREFZrbNtSTNKXQWRJLUoNDMLq8iSRqXQoNaktSg0MwuryJJGqcClxglSYsoMLNtriVpTqGzIJKkBoVmdnkVSdI4FTgLIklaRIGZ3VtFEbEzIh6IiCMRcVlfnyNJnZmbBVnsNvTlS+deDPxx9fw9EVHMpaLNbElrTsvMHrxF97ndS3MdETPAVcB5wNnAOyPi7D4+S5I60yKoR8y984Bt1e1C4OPdD2L5zGxJa1L7CZFecruvmesdwJHMfDAznwH2Art6+ixJ6kbEYIlxsdvSRsm9XcAnc+A24JSIKOELas1sSWtPu8yGnnK7r+b6DODR2uPZapsklW3lsyCj5F6p2VhqXZK0tHaHhfSS232d0BgN2/IHdoi4kMH0OsDTMTNzb0+1lGAD8OS4i+jRJI9vkscGkz++H1nOzofuvPOmOPHEDUvs8pyIOFh7vCcz91T3h+beiPuMw0h1TVFuT/p/F5M8vkkeG0z++FYzs6Gn3O6ruZ4FNtcenwkcre9QDW4PQEQczMztPdUydo5v7ZrkscF0jG85+2fmzhYfNzT3RtxnHEaqa1pye5LHBpM9vkkeG0zH+Jazf8vMhp5yu6/DQu4AtkXE1og4CdgN7O/psySpBKPk3n7g3dXZ568BvpmZx1a70AZmtqRp1Etu9zJznZnHI+IS4CZgBrg2Mw/38VmSVILFci8iLqqevxq4EXgLcAR4CnjPuOqtM7MlTaO+cru3L5HJzBurgkaxZ/gua5rjW7smeWzg+DrVlHtVOM/dT+Di1axpVMvMbJjs351JHhtM9vgmeWzg+DrXR27H4DWSJEmS2irvOyMlSZKkNarz5joinhMRByLi7og4HBG/W23/UER8NSLuqm5vqb3m8uprJR+IiDcv8r6nRsTNEfGV6uf6rmsfpsexfSQivlR9reb/iYhTVmtM8+roZXy1fd8fERkRS102pzd9ji8ifr3a53BE/P5qjGfe5/f1u3lORNxWvfZgROxYrTHNq2NZ44uIF0fE30XEtyPio0u879hzZdwmObOrOiY2t81sM7vhfc3sEmRmpzcG1wN8fnX/ROB24DXAh4D3N+x/NnA3cDKwFfhnYKZhv98HLqvuXwZ8uOvaxzi2nwPWVfc/PI6x9Tm+at/NDE4Y+BdgwySND/hZ4HPAydXj0yZobH8LnFfdfwtw6xr5u/sPwGuBi4CPLvG+Y8+Vcd96/N0p4s+2x/GNPbf7Glu1r5m9NsdmZhdw63zmOge+XT08sbotdWD3LmBvZj6dmQ8xOBuz6V9au4Drq/vXA2/tqOSR9TW2zPzbzDxePbyNwTUUV12Pf3cA/x34b0Per1c9ju/XgCsz8+nqcx7vsOyR9Di2BF5Y3X8RY7om83LHl5nfycx/AL475K3HnivjNsmZDZOd22b2Ama2mV2EXo65joiZiLgLeBy4OTNvr566pFpCu7Y2lT/q10qentV1Baufp/VR+zA9ja3uPwN/3WnRy9DH+CLiF4CvZubdfdY+ip7+/s4Cfjoibo+Iv4+In+xtAEvoaWy/CXwkIh4F/gC4vKfyh1rm+EZVRK6M2yRnNkx2bpvZZvY8ZnYBemmuM/PZzDyHwb/kd0TEy4GPAy8DzgGOAX9Y7V7q1wE36nNsEfEB4Djwp50WvQxdjy8ingd8APhgb0UvQ09/f+uA9QyWvP4rsC8iml7bq57G9mvAb2XmZuC3gD/pvPARLXN8WoZJzmyY7Nw2s83seczsAvR6tZDM/AZwK7AzMx+r/qC/B3yC7y9njPq1ko9FxCaA6ueqL+PUdTw2IuJ84OeBX87Msf+PqsPxvYzB8WF3R8TD1T53RsQP9Vj+UB3//c0CN1TLYAeA7wFjOQEIOh/b+cAN1f3PsPgS8qoZcXyjKipXxm2SMxsmO7fNbMDMBjO7CH1cLWRjVGdNR8RzgTcCX5r7w6i8Dbi3ur8f2B0RJ0fEVmAbcKDhrfcz+KWh+vnZrmsfpq+xRcRO4FLgFzLzqT7HsJQ+xpeZX8zM0zJzS2ZuYRAQr8rMf+15OAv0+Lv5F8Drq/c9CzgJeLKfUTTrcWxHgZ+p7r8e+Eof9Q+zgvGNauy5Mm6TnNkw2bltZgNm9nxmdgmy+zNEfwL4J+AeBn9oH6y2/y/gi9X2/cCm2ms+wODM1weoznKttl8DbK/uvxi4hcEvyi3AqV3XPsaxHWFwLNVd1e3q1R5bn+Ob9xkPM74zz/v6+zsJ+FT1nncCr5+gsb0WOMTgLPXbgXPX0N/dw8DXgG8zaBDObhjf2HNl3Lcef3eK+LPtcXxjz+2+xjbvMx7GzF5LYzOzC7j5DY2SJElSR/yGRkmSJKkjNteSJElSR2yuJUmSpI7YXEuSJEkdsbmWJEmSOmJzLUmSJHXE5lqSJEnqiM21JEmS1BGba0mSJKkjNteSJElSR2yuJUmSpI7YXEuSJEkdsbmWpI5ExLUR8XhE3LvI8xERfxwRRyLinoh41WrXKEka6Cuzba4lqTvXATuXeP48YFt1uxD4+CrUJElqdh09ZLbNtSR1JDM/D3xtiV12AZ/MgduAUyJi0+pUJ0mq6yuzba4lafWcATxaezxbbZMklWdFmb2ut3KWYcOGDblly5ZxlyFpwhw6dOjJzNw46v4/HJFPLfH8MTgMfLe2aU9m7llGSdGwLZfx+mKY25K6NimZXURzvWXLFg4eODDuMiRNmJiZ+Zfl7P8U8L4lnv8QfDczt7coaRbYXHt8JnC0xfuNjbktqWuTktkeFiJJlWAQiovdOrAfeHd1BvprgG9m5rFu3lqSpkupmV3EzLUklaJpDXDk10Z8GngdsCEiZoHfAU4EyMyrgRuBtwBHGEy6vKdVsZI05UrMbJtrSappM9uRme8c8nwCF7f4CElSTYmZbXMtSZWg3SyIJGn1lJrZNteSVOOJKJK0dpSY2TbXklRTYlBLkpqVmNk215JUKXWJUZK0UKmZbXMtSTUlzoJIkpqVmNk215JUU+IsiCSpWYmZbXMtSZW5LySQJJWv1My2uZakmhKDWpLUrMTMtrmWpJoSlxglSc1KzGyba0mqlLrEKElaqNTMtrmWpJoSZ0EkSc1KzGyba0mqKXEWRJLUrMTMtrmWpEqpS4ySpIVKzWyba0mqKXGJUZLUrMTMtrmWpJoSZ0EkSc1KzGyba0mqKXEWRJLUrMTMtrmWpEqpx+9JkhYqNbNtriWppsSgliQ1KzGzh9YUEc+JiAMRcXdEHI6I3622nxoRN0fEV6qf62uvuTwijkTEAxHx5j4HIEldiSG3tcDMljQtSs3sURr+p4HXZ+YrgHOAnRHxGuAy4JbM3AbcUj0mIs4GdgM/DuwEPhYRM30UL0ldO2GJ2xphZkuaGiVm9tDPzoFvVw9PrG4J7AKur7ZfD7y1ur8L2JuZT2fmQ8ARYEenVUtST0qcBVkOM1vSNCkxs0dq7CNiJiLuAh4Hbs7M24HTM/MYQPXztGr3M4BHay+frbZJUtHmTo4pbRZkucxsSdOg1Mwe6bMz89nMPAc4E9gRES9fYvemfyzkgp0iLoyIgxFx8IknnhitWknqWYlBvVx9ZDaY25LKU2JmL+uzM/MbwK0Mjst7LCI2AVQ/H692mwU21152JnC04b32ZOb2zNy+cePGFZQuSd1rs8QYETurkwKPRMRlDc+/KCL+snay4Xs6LX6eLjO7ej9zW1JR2h4W0kduj3K1kI0RcUp1/7nAG4EvAfuB86vdzgc+W93fD+yOiJMjYiuwDTgwygAlaZzaLDFWJwFeBZwHnA28szpZsO5i4L7qZMPXAX8YESd1NgDMbEnTo+1hIX3l9ijXud4EXF8VcAKwLzP/KiK+AOyLiAuAR4C3A2Tm4YjYB9wHHAcuzsxnR/gcSRq7FifB7ACOZOaDABGxl8HJgvfV9kngBRERwPOBrzHIyS6Z2ZKmRssTF3vJ7aHNdWbeA7yyYfu/AW9Y5DVXAFcMe29JKk2L4/SaTgx89bx9Pspgpvgo8ALgHZn5vZV/5EJmtqRp0vLY6l5yey2doyNJvRphiXHD3Al91e3CeS+fb/6JgW8G7gJewuAa1B+NiBd2OQZJmhYtM3vuLeZrndt+/bkk1QxZYnwyM7cv8twoJwa+B7gyMxM4EhEPAT+KxzhL0oq0yGzoKbeduZakmhYnx9wBbIuIrdXJLrsZLCXWPUJ1aEZEnA78CPBgR6VL0tRpeSm+XnLbmWtJqlnpyTGZeTwiLgFuAmaAa6uTBS+qnr8a+D3guoj4YvVRl2bmk13ULUnTqM0JjX3lts21JFXmjt9bqcy8Ebhx3rara/ePAj/X4iMkSZW2mQ395LbNtSTVeKycJK0dJWa2zbUk1bS8ZqokaRWVmNk215JU6WKJUZK0OkrNbJtrSaopcRZEktSsxMy2uZakmhJnQSRJzUrMbJtrSaqUusQoSVqo1My2uZakmhKXGCVJzUrMbJtrSWCoMTAAACAASURBVKopcRZEktSsxMy2uZakSlDmLIgkaaFSM9vmWpJqSpwFkSQ1KzGzba4lqabEoJYkNSsxs22uJalS6hKjJGmhUjPb5lqSakqcBZEkNSsxs22uJammxFkQSVKzEjPb5lqSKqV+IYEkaaFSM9vmWpJqSgxqSVKzEjPb5lqSakpcYpQkNSsxs22uJalS6hKjJGmhUjPb5lqSakqcBZEkNSsxs22uJammxFkQSVKzEjPb5lqSKqUuMUqSFio1s22uJammxCVGSVKzEjPb5lqSakqcBZEkNSsxs0usSZLGJpa4DX1txM6IeCAijkTEZYvs87qIuCsiDkfE33dWuCRNoTaZDf3ktjPXklRpc/xeRMwAVwFvAmaBOyJif2beV9vnFOBjwM7MfCQiTmtbsyRNq7bHXPeV285cS1LNCUvchtgBHMnMBzPzGWAvsGvePu8CbsjMRwAy8/HOCpekKdQis6Gn3La5lqSaFkuMZwCP1h7PVtvqzgLWR8StEXEoIt7dQcmSNLVaHhbSS257WIgkVUZYYtwQEQdrj/dk5p7ay+fLeY/XAecCbwCeC3whIm7LzC+vqGBJmmItM3vuLeZrnds215JUM2S248nM3L7Ic7PA5trjM4GjDfs8mZnfAb4TEZ8HXgHYXEvSCrTIbOgptz0sRJJqWhy/dwewLSK2RsRJwG5g/7x9Pgv8dESsi4jnAa8G7u+seEmaMi2Pue4lt4d+dkRsjoi/i4j7q0uQ/Ea1/dSIuDkivlL9XF97zeXVJU0eiIg3jzY+SRqvuSXGlQR1Zh4HLgFuYhC8+zLzcERcFBEXVfvcD/wNcA9wALgmM+/tdAxmtqQp0Sazob/cjsz5h5bM2yFiE7ApM++MiBcAh4C3Ar8KfC0zr6yuC7g+My+NiLOBTzM4A/MlwOeAszLz2cU+Y/v27XnwwIGhfwiStBwxM3NoyJLgDzgnIj+3xPMbYVnvNw6rkdlgbkvq3qRk9tDGPjOPZead1f1/Z9DZn8HgUiXXV7tdzyC8qbbvzcynM/Mh4AiD0Jak4rVcYhw7M1vSNCkxs5f12RGxBXglcDtwemYeg0GYA3MX1R7lsiaSVJylLuk06rd9lcTMljTJSs3skZvriHg+8OfAb2bmt5batWHbgmNPIuLCiDgYEQefeOKJUcuQpF6VOAuyEl1ndvWe5rakopSY2SN9dkScyCCk/zQzb6g2P1Yd2zd3jN/cN9aMclkTMnNPZm7PzO0bN25caf2S1Jm2J8eUoo/MBnNbUllKzexRrhYSwJ8A92fmH9We2g+cX90/n8GlSua2746IkyNiK7CNwdmVklS8EpcYl8PMljRNSszsUb5E5qeAXwG+GBF3Vdt+G7gS2BcRFwCPAG8HqC5hsg+4DzgOXDzsrHNJKsVamqFehJktaWqUmNlDm+vM/AcW/wfAGxZ5zRXAFS3qkqSxWCsz1IsxsyVNkxIz268/l6TK3PF7kqTylZrZNteSVFNiUEuSmpWY2TbXklRT4hKjJKlZiZltcy1JlVKXGIv0r/8Kf/AH465C0hQrNbNtriWp5oRYYh4kG79bZTr90A/B+98/7iokTZJLL132S0rMbJtrSaorMKglSYsoMLNtriVpTgSsWyIWn3lm9WqRJC2t0My2uZakuhNKPIJPktSowMy2uZakOcNmQSRJ5Sg0s8urSJLGpdCgLtLXvw5/9mfjrkLSNCs0s8urSJLGqcAlxiIdPw7f+Ma4q5A07QrMbJtrSZpT6CxIkTZuhPe+d9xVSJok73vf8vYvNLPLq0iSxqnAWRBJ0iIKzGyba0maU+gsiCSpQaGZXV5FkjQuhQa1JKlBoZldXkWSNE4FLjFKkhZRYGbbXEvSnJazIBGxE/j/gBngmsy8cpH9fhK4DXhHZq7N69l5KT5J49bBzHUfuW1zLUlzIlY8CxIRM8BVwJuAWeCOiNifmfc17Pdh4KaW1Y7X+vXwi7847iokTbMWmT14eT+5Xd5cuiSN07p1i9+WtgM4kpkPZuYzwF5gV8N+vw78OfB4t4VL0hRaeWZDT7ntzLUkzWm3xHgG8Gjt8Szw6h98+zgDeBvweuAnV/pBkiS6OCykl9y2uZakOcOXGDdExMHa4z2ZuWfu1Q3757zH/wO4NDOfjWjaXZI0snaZDT3lts21JNUtPQvyZGZuX+S5WWBz7fGZwNF5+2wH9lYBvQF4S0Qcz8y/WGG1kjTdVp7Z0FNu21xL0px2J8fcAWyLiK3AV4HdwLvqO2Tm1u9/VFwH/JWNtSStUMsTGukpt22uJWlOi+P3MvN4RFzC4GzyGeDazDwcERdVz1/dXaGSpLbHXPeV2zbXklTXLqhvBG6ct60xnDPzV1f8QZKkgZbXue4jt22uJWlO+yVGSdJqKTSzba4laU4H3/YlSVolhWZ2eRVJ0jgVOAsiSVpEgZltcy1JcwqdBZEkNSg0s8urSJLGpdCgliQ1KDSzy6tIksapwCVGSdIiCsxsm2tJmlPoLIgkqUGhmV1eRZI0TgXOgkiSFlFgZttcS9KcQmdBJEkNCs3s8iqSpHEpNKglSQ0KzezyKpKkcSpwiVGStIgCM3tocx0R1wI/DzyemS+vtp0K/G9gC/Aw8EuZ+fXqucuBC4Bngf+SmTcN+4xvHjrEX87MrHAIktSRQmdBlms1cpsvfxl27uyjfEkaTaGZPUpF1wEfBT5Z23YZcEtmXhkRl1WPL42Is4HdwI8DLwE+FxFnZeazS33Ai849l/904MBK6pekxS33H+0RRc6CrMB19JzbnHUW/M3f9FG7pGk1IZk9tLnOzM9HxJZ5m3cBr6vuXw/cClxabd+bmU8DD0XEEWAH8IWlPuP/HjrEMWeuJZWgwFmQ5VqN3OaJJ+CaazqrWZJWpMDMXmlFp2fmMYDMPBYRp1XbzwBuq+03W21b0onPfS6bfvRHV1iKJC3in/5pefsXusTYkU5zmxe+EN74xs6LlKSRFZrZXVcUDduycceIC4ELAV760pfCwYMdlyJp6k3IEmPPVp7bW7b0WJYkDVFoZq+0osciYhNA9fPxavsssLm235nA0aY3yMw9mbk9M7dv3LhxhWVIUsfWrVv8traZ25ImT4GZvdLmej9wfnX/fOCzte27I+LkiNgKbAM8U1HS2jA3C7LYbW0ztyVNlkIze5RL8X2awUkwGyJiFvgd4EpgX0RcADwCvB0gMw9HxD7gPuA4cPHQM84lqRSFHr+3XOa2pKlQaGaPcrWQdy7y1BsW2f8K4Io2RUnS2EzAlYvMbUlTo8DMLq/dl6RxKfTkGElSg0Iz2+ZakuoKXGKUJC2iwMwuryJJGpdCZ0EkSQ0KzezyKpKkcZk7OWaFl3WKiJ0R8UBEHKm+Ynz+878cEfdUt3+MiFf0Mg5JmgYtM3vwFt3ntjPXklS3wiXGiJgBrgLexODa0XdExP7MvK+220PAz2Tm1yPiPGAP8OqWFUvS9GpxWEhfuW1zLUlz2i0x7gCOZOaDg7eKvcAuBpe4AyAz/7G2/20MvrBFkrQS7Q8L6SW3ba4laU67a6aeATxaezzL0rMbFwB/vdIPk6Sp1/46173kts21JNUtPQuyISIO1h7vycw91f1o2D+b3iQifpZBSL92RTVKkgZWntnQU27bXEvSnOGzIE9m5vZFnpsFNtcenwkcXfgR8RPANcB5mflvKy1VkqZeu8yGnnLb5lqS5rRbYrwD2BYRW4GvAruBd/3g28dLgRuAX8nML7cpVZKmXvvDQnrJbZtrSapb4ckxmXk8Ii4BbgJmgGsz83BEXFQ9fzXwQeDFwMciAuD4kFkVSdJSWpzQ2Fdu21xL0pyWsyCZeSNw47xtV9fuvxd474o/QJL0fe1nrnvJbZtrSaor8Nu+JEmLKDCzba4laU4HsyCSpFVSaGaXV5EkjUuhQS1JalBoZpdXkSSNU4FLjEV64gm45ppxVyFp2hWY2TbXkjSn0FmQIm3cCO/13ExJHXrf+5a3f6GZXV5FkjQuEUXOgkiSGhSa2TbXklRX4CyIJGkRBWZ2eRVJ0rgUusQoSWpQaGaXV5EkjUuhS4ySpAaFZrbNtSTVFTgLIklaRIGZXUZF3/oWfO5z465C0rQrdBakSE89BXfdNe4qJE2zQjO7jOb6ec+D7dvHXYWkaVfo8XtFmpmBU04ZdxWSplmhmV1ERf/37rs59uIXj7sMSSoyqEv09L338s8ve9m4y5A07QrM7CIqOvHHfoxNn/rUuMuQNGnOPXd5+xe6xFiik9ev52VvfOO4y5A0ST7zmeXtX2hmF9FcP3H//fzP5f5PUJK6VugSY5E8V0bSuBWa2UVUdCrwS+MuQtLEuWglLypwFqRI69fD29427iokTZJPfGL5rykws4tormfOPZf1Bw6MuwxJk2ZmZnn7FzoLUqT/+B/h6qvHXYWkSbLc5rrQzC6vIkkal0KDWpLUoNDMLq8iSRqnApcYJUmLKDCzba4laU6hsyCSpAaFZnZ5FUnSOBU4CyJJWkSBmW1zLUlzCp0FkSQ1KDSzy6tIksal0KCWJDUoNLN7m0uPiJ0R8UBEHImIy/r6HEnq1AknLH4bYljuxcAfV8/fExGv6mUMK2BmS1qTWmQ29JPbvTTXETEDXAWcB5wNvDMizu7jsySpM3OzIIvdlnzpSLl3HrCtul0IfLz7QSyfmS1pTWqR2YOX95Pbfc1c7wCOZOaDmfkMsBfY1dNnSVI3ItrMgoySe7uAT+bAbcApEbGp+4Esm5ktae1pl9nQU2731VyfATxaezxbbZOksq18FmSU3Cs1G0utS5KW1mLmmp5yu6+jwKNhW/7ADhEXMpheB3g6Zmbu7amWEmwAnhx3ET2a5PFN8thg8sf3I8vZ+dCdd94UJ564YYldnhMRB2uP92Tmnur+0NwbcZ9xGKmuKcrtSf/vYpLHN8ljg8kf32pmNvSU230117PA5trjM4Gj9R2qwe0BiIiDmbm9p1rGzvGtXZM8NpiO8S1n/8zc2eLjhubeiPuMw0h1TUtuT/LYYLLHN8ljg+kY33L2b5nZ0FNu93VYyB3AtojYGhEnAbuB/T19liSVYJTc2w+8uzr7/DXANzPz2GoX2sDMljSNesntXmauM/N4RFwC3ATMANdm5uE+PkuSSrBY7kXERdXzVwM3Am8BjgBPAe8ZV711ZrakadRXbvd25e3MvLEqaBR7hu+ypjm+tWuSxwaOr1NNuVeF89z9BC5ezZpGtczMhsn+3ZnkscFkj2+SxwaOr3N95HYMXiNJkiSprd6+oVGSJEmaNp031xHxnIg4EBF3R8ThiPjdavuHIuKrEXFXdXtL7TWXV18r+UBEvHmR9z01Im6OiK9UP9d3XfswPY7tIxHxpeprNf9PRJyyWmOaV0cv46vt+/6IyIhY6rI5velzfBHx69U+hyPi91djPPM+v6/fzXMi4rbqtQcjYsdqjWleHcsaX0S8OCL+LiK+HREfXeJ9x54r4zbJmV3VMbG5bWab2Q3va2aXIDM7vTG4HuDzq/snArcDrwE+BLy/Yf+zgbuBk4GtwD8DMw37/T5wWXX/MuDDXdc+xrH9HLCuuv/hcYytz/FV+25mcMLAvwAbJml8wM8CnwNOrh6fNkFj+1vgvOr+W4Bb18jf3X8AXgtcBHx0ifcde66M+9bj704Rf7Y9jm/sud3X2Kp9zey1OTYzu4Bb5zPXOfDt6uGJ1W2pA7t3AXsz8+nMfIjB2ZhN/9LaBVxf3b8eeGtHJY+sr7Fl5t9m5vHq4W0MrqG46nr8uwP478B/G/J+vepxfL8GXJmZT1ef83iHZY+kx7El8MLq/osY0zWZlzu+zPxOZv4D8N0hbz32XBm3Sc5smOzcNrMXMLPN7CL0csx1RMxExF3A48DNmXl79dQl1RLatbWp/FG/VvL0rK4rWP08rY/ah+lpbHX/GfjrTotehj7GFxG/AHw1M+/us/ZR9PT3dxbw0xFxe0T8fUT8ZG8DWEJPY/tN4CMR8SjwB8DlPZU/1DLHN6oicmXcJjmzYbJz28w2s+cxswvQS3Odmc9m5jkM/iW/IyJeDnwceBlwDnAM+MNq91K/DrhRn2OLiA8Ax4E/7bToZeh6fBHxPOADwAd7K3oZevr7WwesZ7Dk9V+BfRHR9Npe9TS2XwN+KzM3A78F/EnnhY9omePTMkxyZsNk57aZbWbPY2YXoNerhWTmN4BbgZ2Z+Vj1B/094BN8fzlj1K+VfCwiNgFUP1d9Gaeu47EREecDPw/8cmaO/X9UHY7vZQyOD7s7Ih6u9rkzIn6ox/KH6vjvbxa4oVoGOwB8DxjLCUDQ+djOB26o7n+GxZeQV82I4xtVUbkybpOc2TDZuW1mA2Y2mNlF6ONqIRujOms6Ip4LvBH40twfRuVtwL3V/f3A7og4OSK2AtuAAw1vvZ/BLw3Vz892XfswfY0tInYClwK/kJlP9TmGpfQxvsz8YmaelplbMnMLg4B4VWb+a8/DWaDH382/AF5fve9ZwEnAk/2MolmPYzsK/Ex1//XAV/qof5gVjG9UY8+VcZvkzIbJzm0zGzCz5zOzS5DdnyH6E8A/Afcw+EP7YLX9fwFfrLbvBzbVXvMBBme+PkB1lmu1/Rpge3X/xcAtDH5RbgFO7br2MY7tCINjqe6qblev9tj6HN+8z3iY8Z153tff30nAp6r3vBN4/QSN7bXAIQZnqd8OnLuG/u4eBr4GfJtBg3B2w/jGnivjvvX4u1PEn22P4xt7bvc1tnmf8TBm9loam5ldwM1vaJQkSZI64jc0SpIkSR2xuZYkSZI6YnMtSZIkdcTmWpIkSeqIzbUkSZLUEZtrSZIkqSM215IkSVJHbK4lSZKkjthcS5IkSR2xuZYkSZI6YnMtSZIkdcTmWpIkSeqIzbUkdSQiro2IxyPi3kWej4j444g4EhH3RMSrVrtGSdJAX5ltcy1J3bkO2LnE8+cB26rbhcDHV6EmSVKz6+ghs22uJakjmfl54GtL7LIL+GQO3AacEhGbVqc6SVJdX5ltcy1Jq+cM4NHa49lqmySpPCvK7HW9lbMMGzZsyC1btoy7DEkT5tChQ09m5sZR9//hiHxqieePwWHgu7VNezJzzzJKioZtuYzXF8PcltS1ScnsIprrLVu2cPDAgXGXIWnCxMzMvyxn/6eA9y3x/Ifgu5m5vUVJs8Dm2uMzgaMt3m9szG1JXZuUzPawEEmqBINQXOzWgf3Au6sz0F8DfDMzj3Xz1pI0XUrN7CJmriWpFE1rgCO/NuLTwOuADRExC/wOcCJAZl4N3Ai8BTjCYNLlPa2KlaQpV2Jm21xLUk2b2Y7MfOeQ5xO4uMVHSJJqSsxsm2tJqswtMUqSyldqZttcS1JNmyVGSdLqKjGzba4lqabEWRBJUrMSM9vmWpIqQZmzIJKkhUrNbJtrSaopcRZEktSsxMy2uZakmhKDWpLUrMTMtrmWpEqpS4ySpIVKzWyba0mqKXEWRJLUrMTMtrmWpJoSZ0EkSc1KzGyba0mqlPqFBJKkhUrNbJtrSaopMaglSc1KzGyba0mqKXGJUZLUrMTMtrmWpEqpS4ySpIVKzWyba0mqKXEWRJLUrMTMtrmWpJoSZ0EkSc1KzGyba0mqlLrEKElaqNTMtrmWpJoSlxglSc1KzGyba0mqKXEWRJLUrMTMHlpTRDwnIg5ExN0RcTgifrfafmpE3BwRX6l+rq+95vKIOBIRD0TEm/scgCR1JYbc1gIzW9K0KDWzR2n4nwZen5mvAM4BdkbEa4DLgFsycxtwS/WYiDgb2A38OLAT+FhEzPRRvCR17YQlbmuEmS1papSY2UM/Owe+XT08sbolsAu4vtp+PfDW6v4uYG9mPp2ZDwFHgB2dVi1JPSkxqJfDzJY0TUrM7JE+OyJmIuIu4HHg5sy8HTg9M48BVD9Pq3Y/A3i09vLZapskFa3UJcblMrMlTYNSM3uk5jozn83Mc4AzgR0R8fIldm8aTy7YKeLCiDgYEQefeOKJ0aqVpJ6VOAuyXH1kNpjbkspTYmYv67Mz8xvArQyOy3ssIjYBVD8fr3abBTbXXnYmcLThvfZk5vbM3L5x48YVlC5J3WszCxIRO6uTAo9ExGUNz78oIv6ydrLhezotfp4uM7t6P3NbUlHazlz3kdujXC1kY0ScUt1/LvBG4EvAfuD8arfzgc9W9/cDuyPi5IjYCmwDDowyQEkap7kvJFjJLEh1EuBVwHnA2cA7q5MF6y4G7qtONnwd8IcRcVJnA8DMljQ92mQ29Jfbo1znehNwfVXACcC+zPyriPgCsC8iLgAeAd4OkJmHI2IfcB9wHLg4M58d4XMkaexaLCXuAI5k5oMAEbGXwcmC99X2SeAFERHA84GvMcjJLpnZkqZGy8M/esntoc11Zt4DvLJh+78Bb1jkNVcAVwx7b0kqTYuTYJpODHz1vH0+ymCm+CjwAuAdmfm9lX/kQma2pGnS8sTFXnJ7LZ2jI0m9GmGJccPcCX3V7cJ5L59v/omBbwbuAl7C4BrUH42IF3Y5BkmaFi0ze+4t5mud2379uSTVDJkFeTIzty/y3CgnBr4HuDIzEzgSEQ8BP4rHOEvSirTIbOgpt525lqSaFifH3AFsi4it1ckuuxksJdY9QnVoRkScDvwI8GBHpUvS1Gl5Kb5ectuZa0mqzC0xrkRmHo+IS4CbgBng2upkwYuq568Gfg+4LiK+WH3cpZn5ZAelS9LUaZPZ0F9u21xLUk2bk2My80bgxnnbrq7dPwr8XIuPkCTVtP0mxj5y2+Zakmo8Vk6S1o4SM9vmWpJq2s6CSJJWT4mZbXMtSZW2x+9JklZPqZltcy1JNSUGtSSpWYmZbXMtSTUlLjFKkpqVmNk215JUKXWJUZK0UKmZbXMtSTUlzoJIkpqVmNk215JUU+IsiCSpWYmZbXMtSZVSlxglSQuVmtk215JUU+ISoySpWYmZbXMtSTUlzoJIkpqVmNk215JUCcqcBZEkLVRqZttcS1JNibMgkqRmJWa2zbUk1ZQY1JKkZiVmts21JFVKXWKUJC1UambbXEtSTYmzIJKkZiVmts21JNWUOAsiSWpWYmbbXEtSpdQvJJAkLVRqZttcS1JNiUEtSWpWYmbbXEtSTYlLjJKkZiVmts21JFVKXWKUJC1UambbXEtSTYlBLUlqVmJm21xLUk2JS4ySpGYlZnaJDb8kjcXcEuNit6Gvj9gZEQ9ExJGIuGyRfV4XEXdFxOGI+PuuapekadM2s6Gf3HbmWpJqVjoLEhEzwFXAm4BZ4I6I2J+Z99X2OQX4GLAzMx+JiNNaFyxJU6zNzHVfue3MtSTVtJgF2QEcycwHM/MZYC+wa94+7wJuyMxHADLz8c4Kl6Qp1HLmupfctrmWpErLJcYzgEdrj2erbXVnAesj4taIOBQR725ftSRNpw4OC+kltz0sRJJqhiwxboiIg7XHezJzzxIvzXmP1wHnAm8Angt8ISJuy8wvr6xaSZpuLTJ7sZe3zm2ba0mqGTLb8WRmbl/kuVlgc+3xmcDRhn2ezMzvAN+JiM8DrwBsriVpBVpkNvSU2x4WIkk1scRtiDuAbRGxNSJOAnYD++ft81ngpyNiXUQ8D3g1cH9nxUvSlGmR2dBTbg9triNic0T8XUTcX12C5Deq7adGxM0R8ZXq5/raay6vLmnyQES8ebTxSdJ4tTl+LzOPA5cANzEI3n2ZeTgiLoqIi6p97gf+BrgHOABck5n3djoGM1vSlGh7zHVfuR2Z8w8tmbdDxCZgU2beGREvAA4BbwV+FfhaZl5ZXRdwfWZeGhFnA59mcAbmS4DPAWdl5rOLfcb27dvz4IEDQ/8QJGk5Ymbm0JAlwR9wTkTessTzG2BZ7zcOq5HZYG5L6t6kZPbQxj4zj2XmndX9f2fQ2Z/B4FIl11e7Xc8gvKm2783MpzPzIeAIg9CWpOK1XGIcOzNb0jQpMbOXdcx1RGwBXgncDpyemcdgEObA3EW1R7msiSQVp4tv+yqJmS1pkpWa2SN/dkQ8H/hz4Dcz81tL7dqwbcGxJxFxYUQcjIiDTzzxxKhlSFKvSpwFWYmuM7t6T3NbUlFKzOyRmuuIOJFBSP9pZt5QbX6sOrZv7hi/uW+sGeWyJmTmnszcnpnbN27cuNL6Jakzpc6CLFcfmQ3mtqSylJrZo1wtJIA/Ae7PzD+qPbUfOL+6fz6DS5XMbd8dESdHxFZgG4OzKyWpeCUG9XKY2ZKmSYmZPcqXyPwU8CvAFyPirmrbbwNXAvsi4gLgEeDtANUlTPYB9wHHgYuHnXUuSaVYa4d/NDCzJU2NEjN7aHOdmf/A4rW/YZHXXAFc0aIuSVp1c0uMa5mZLWlalJrZfv25JNWUOAsiSWpWYmbbXEtSTYmzIJKkZiVmts21JFVKXWKUJC1UambbXEtSTYlLjJKkZiVmts21JNWcEEtEdTZ+t4okaUxKzGyba0mqKzCoJUmLKDCzba4laU4ErFsiFp95ZvVqkSQtrdDMtrmWpDmFBrUkqUGhmW1zLUl1J5R47rkkqVGBmV1Ec/3soUN8fWZm3GVImnbDZkH0/zO3JY1doZldREUzp5/O+ne/e9xlSJo0H/nI8l9T4CxIiWZ+7MdY/6lPjbsMSZPk3HOX/5oCM7uI5pozz4Qrrxx3FZImzXKb60JnQYr0vOfBOeeMuwpJ06zQzC6vIkkal0KDWpLUoNDMLq8i/b/27j9Wsvqs4/j72QvdWo3uwkIlyxI2m60JkkLaZds/MLVo7UJMtyQYtm0UtQ3SlCY2qRZCgvWPJtBq/JHS4oIE/BE3uGK7MdtQpFZjGn4sCJQlRVZY5QKCWDBB4jZLH/+YM2G4zNx7Z+Z87/neue9XMrlnzpyZ8zz3Dh+++z1z5kjqUoWHGCVJI1SY2Q6uJamv0lkQSdIQlWZ2fRVJUpcqnAWRJI1QYWbXMbh+6SXYv7/rKiStdVPOgkTELuCPyzx3aQAADgRJREFUgDng5swceqZ2RJwH3ANcmpmGnyRNooWZ6xK5XcfgeuNGuOSSrquQtNZNEdQRMQfcAHwAmAfuj4gDmfnYkO2uB+6cslpJWtumnxApktv1zaVLUpfWrRt9W9xO4EhmPpmZPwD2AbuHbPdp4G+AF9otXJLWoMkzGwrldh0z15JUg+lmQTYDTw/cnwfe88aXj83AxcAFwHmT7kiSRBsfCymS2w6uJakvYqnZjk0RcWjg/t7M3Nt/9pDtc8H9PwQ+l5mvRQzbXJK0bNNlNhTKbQfXkjRo8VmQFzNzx4jH5oEtA/dPB55dsM0OYF8T0JuAiyLieGZ+bcJqJWltmzyzoVBuO7iWpL7pDjHeD2yPiK3AM8Ae4KODG2Tm1td3FbcCf+fAWpImNP3HQorktoNrSepb+hDjSJl5PCKupHc2+RxwS2YejogrmsdvbK9QSdI0mQ3lctvBtSQNmmIWJDMPAgcXrBsazpn5qxPvSJLUM+X3XJfIbQfXktQ35SyIJGkFVZrZDq4lqa+Fq31JklZIpZldX0WS1KUKg1qSNEKFmV1fRZLUlUoPMUqShqg0sx1cS1JfpYcYJUlDVJrZ9VUkSV2qcBZEkjRChZnt4FqS+iqdBZEkDVFpZtdXkSR1pdKgliQNUWlm11eRJHWpwkOMkqQRKsxsB9eS1FfpLIgkaYhKM3vJ4X5E3BIRL0TEowPrToqIuyLiiebnxoHHro6IIxHxeER8sFThklTEunWjb6uEuS1pzagws5ez51uBXQvWXQXcnZnbgbub+0TEWcAe4Keb53wlIuZaq1aSSurPgoy6rR63Ym5LmnWVZvaSg+vM/Cfg+wtW7wZua5ZvAz48sH5fZh7LzKeAI8DOlmqVpLIqDepxmduS1oRKM3vSOfO3Z+ZzAM3PU5v1m4GnB7abb9ZJ0upQ4SHGlpjbkmZPhZnd9rA+hqzLoRtGXA5cDnDGGWe0XIYkTaDSk2MKM7clrU6VZvakw/rnI+I0gObnC836eWDLwHanA88Oe4HM3JuZOzJzxymnnDJhGZLUoogqZ0FaYm5Lmi2VZvakez4AXNYsXwZ8fWD9nohYHxFbge3AfdOVKEkrqMLP77XE3JY0eyrM7CX3HBF/BfwssCki5oHfAa4Dbo+IjwP/AfwSQGYejojbgceA48CnMvO1Jas4fhxefnnSHiSpHZUeYhzXiuT2Sy/B/v1lGpCk5ag0s5esKDM/MuKhnxux/ReAL4xXxQmwYcNYT5Gk1vUPMa5yK5LbGzfCJZeMWZkktajSzK5vuC9JXZrzK54ladWoMLMdXEtSX6WzIJKkISrNbAfXkjSows/vSZJGqDCz66tIkrpS6ckxkqQhKs3s+ubSJakrU35nakTsiojHI+JIRFw15PGPRcQjze07EXFOkT4kaS1o4XuuS+R2fcN9SerShLMgETEH3AB8gN6FWe6PiAOZ+djAZk8B78vMlyLiQmAv8J4pK5aktWuKmetSue3gWpL6pjs5ZidwJDOf7L1U7AN20/v+aAAy8zsD299D72qIkqRJTH9CY5HcdnAtSX3TfX5vM/D0wP15Fp/d+DjwjUl3Jklr3vSfuS6S2w6uJWnQ4kG9KSIODdzfm5l7m+UYsn0Oe5GIeD+9kD5/oholST2TZzYUym0H15LUt/Qhxhczc8eIx+aBLQP3TweeffMu4p3AzcCFmfnfk5YqSWvedJkNhXLbwbUk9U13iPF+YHtEbAWeAfYAH33jy8cZwB3AL2fmv05TqiStedN/LKRIbju4lqRBE54ck5nHI+JK4E5gDrglMw9HxBXN4zcC1wInA1+JCIDjS8yqSJIWM8UJjaVy28G1JPVNOQuSmQeBgwvW3Tiw/AngExPvQJL0uhYuIlMitx1cS1JfpVf7kiQNUWlm11eRJHVpuu9MlSStpAoz28G1JPVVOgsiSRqi0syuo6KHHoJNm7quQpKqnAWp0mOPwQ7PxZTUsQozu47B9dlnw7e+1XUVkmbNySePt32lsyBV2rYN9u/vugpJs2TbtvG2rzSz66johBNgw4auq5C01lUa1FVavx7OPLPrKiStZZVmdn0VSVKXKjzEKEkaocLMdnAtSX2VzoJIkoaoNLPrq0iSuhJR5SyIJGmISjPbwbUkDapwFkSSNEKFmV1HRceOwdGjXVchaa2r9BBjlcxtSV2rNLPrq0iSulLpIcYqvfoqHDrUdRWS1rJKM7uKwfV/PfoofzLudxtKUgkVzoLU6NUnn+SBSy/tugxJa12FmV1FRads2sRvXHxx12VImjFX3HTTeE+odBakRm875xze7cW/JLVpkgt/VZjZVQyu2bwZrruu6yokzZpJBtcVzoJUyYt/SepapZldX0WS1KUKg1qSNEKFmV1HRc6ASKpBpYcYJUlDVJrZdQyuJakGlR5ilCQNUWlm11eRJHWpwlkQSdIIFWa2g2tJ6qt0FkSSNESlmV1fRZLUlUqDWpI0RKWZXV9FktSlCg8xSpJGqDCzHVxLUl+lsyCSpCEqzexiw/2I2BURj0fEkYi4qtR+JKlV69aNvi1hqdyLnj9uHn8kIt5VpIcJmNmSVqUpMhvK5HaRwXVEzAE3ABcCZwEfiYizSuxLklrTnwUZdVv0qcvKvQuB7c3tcuCr7TcxPjNb0qo0RWb3nl4mt0vNXO8EjmTmk5n5A2AfsLvQviSpHdMF9XJybzfwZ9lzD7AhIk5rv5GxmdmSVp8pB9cUyu1Sg+vNwNMD9+ebdZJUt8kPMS4n92rNxlrrkqTFTfexkCK5XepT4DFkXb5hg4jL6U2vAxyLublHC9VSg03Ai10XUdAs9zfLvcHs9/dT42z8wIMP3hknnrhpkU3eGhGHBu7vzcy9zfKSubfMbbqwrLrWUG7P+n8Xs9zfLPcGs9/fSmY2FMrtUoPreWDLwP3TgWcHN2ia2wsQEYcyc0ehWjpnf6vXLPcGa6O/cbbPzF1T7G7J3FvmNl1YVl1rJbdnuTeY7f5muTdYG/2Ns/2UmQ2FcrvUx0LuB7ZHxNaIeAuwBzhQaF+SVIPl5N4B4Feas8/fC/xPZj630oUOYWZLWouK5HaRmevMPB4RVwJ3AnPALZl5uMS+JKkGo3IvIq5oHr8ROAhcBBwBXgV+rat6B5nZktaiUrld7Ju3M/NgU9By7F16k1XN/lavWe4N7K9Vw3KvCef+cgKfWsmalmvMzIbZfu/Mcm8w2/3Ncm9gf60rkdvRe44kSZKkadV3QXZJkiRplWp9cB0Rb42I+yLi4Yg4HBG/26z/fEQ8ExEPNbeLBp5zdXNZyccj4oMjXvekiLgrIp5ofm5su/alFOztSxHxveaymn8bERtWqqcFdRTpb2Dbz0ZERsRiX5tTTMn+IuLTzTaHI+KLK9HPgv2Xem+eGxH3NM89FBE7V6qnBXWM1V9EnBwR/xARr0TElxd53c5zpWuznNlNHTOb22a2mT3kdc3sGmRmqzd63wf4Y83yicC9wHuBzwOfHbL9WcDDwHpgK/BvwNyQ7b4IXNUsXwVc33btHfb2C8AJzfL1XfRWsr9m2y30Thj4d2DTLPUHvB/4e2B9c//UGertm8CFzfJFwLdXyd/uR4HzgSuALy/yup3nSte3gu+dKn63BfvrPLdL9dZsa2avzt7M7Apurc9cZ88rzd0Tm9tiH+zeDezLzGOZ+RS9szGH/UtrN3Bbs3wb8OGWSl62Ur1l5jcz83hz9x5636G44gr+7QD+APjtJV6vqIL9fRK4LjOPNft5ocWyl6Vgbwn8eLP8E3T0nczj9peZ/5uZ/wz83xIv3XmudG2WMxtmO7fN7Dcxs83sKhT5zHVEzEXEQ8ALwF2ZeW/z0JXNIbRbBqbyl3tZybdn872Czc9TS9S+lEK9Dfp14ButFj2GEv1FxIeAZzLz4ZK1L0ehv987gJ+JiHsj4h8j4rxiDSyiUG+/CXwpIp4Gfg+4ulD5Sxqzv+WqIle6NsuZDbOd22a2mb2AmV2BIoPrzHwtM8+l9y/5nRFxNvBVYBtwLvAc8PvN5rVeDniokr1FxDXAceAvWy16DG33FxFvA64Bri1W9BgK/f1OADbSO+T1W8DtETHsuUUV6u2TwGcycwvwGeBPWy98mcbsT2OY5cyG2c5tM9vMXsDMrkDRbwvJzJeBbwO7MvP55hf9Q+AmXj+csdzLSj4fEacBND9X/DDOoJZ7IyIuA34R+Fhmdv4/qhb720bv82EPR8TRZpsHI+InC5a/pJb/fvPAHc1hsPuAHwKdnAAErfd2GXBHs/zXjD6EvGKW2d9yVZUrXZvlzIbZzm0zGzCzwcyuQolvCzklmrOmI+JHgJ8Hvtf/ZTQuBh5tlg8AeyJifURsBbYD9w156QP03jQ0P7/edu1LKdVbROwCPgd8KDNfLdnDYkr0l5nfzcxTM/PMzDyTXkC8KzP/s3A7b1Lwvfk14ILmdd8BvAV4sUwXwxXs7Vngfc3yBcATJepfygT9LVfnudK1Wc5smO3cNrMBM3shM7sG2f4Zou8E/gV4hN4v7dpm/Z8D323WHwBOG3jONfTOfH2c5izXZv3NwI5m+WTgbnpvlLuBk9quvcPejtD7LNVDze3Gle6tZH8L9nGU7s48L/X3ewvwF81rPghcMEO9nQ88QO8s9XuBd6+iv91R4PvAK/QGCGcN6a/zXOn6VvC9U8XvtmB/ned2qd4W7OMoZvZq6s3MruDmFRolSZKklniFRkmSJKklDq4lSZKklji4liRJklri4FqSJElqiYNrSZIkqSUOriVJkqSWOLiWJEmSWuLgWpIkSWrJ/wMLSniuzsCMFAAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig, axs = plt.subplots(3, 2, figsize=(12, 10))\n", + "for ax, freq in zip(axs.flatten(), rolling_frequencies):\n", + " freq.residue_contacts.plot_axes(ax=ax)\n", + " ax.set_xlim(*contacts.query_residue_range);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.3" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": true, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": true + }, + "varInspector": { + "cols": { + "lenName": 16, + "lenType": 16, + "lenVar": 40 + }, + "kernels_config": { + "python": { + "delete_cmd_postfix": "", + "delete_cmd_prefix": "del ", + "library": "var_list.py", + "varRefreshCmd": "print(var_dic_list())" + }, + "r": { + "delete_cmd_postfix": ") ", + "delete_cmd_prefix": "rm(", + "library": "var_list.r", + "varRefreshCmd": "cat(var_dic_list()) " + } + }, + "types_to_exclude": [ + "module", + "function", + "builtin_function_or_method", + "instance", + "_Feature" + ], + "window_display": false + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}