diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
new file mode 100644
index 0000000..5c862e6
--- /dev/null
+++ b/.github/workflows/test.yml
@@ -0,0 +1,67 @@
+name: Tests
+
+on:
+ push:
+ branches: [master, main, new_version]
+ pull_request:
+ # Runs the unpinned job below, so a new release of a dependency is found on a
+ # schedule instead of on an unrelated push.
+ schedule:
+ - cron: '0 6 * * 1'
+ workflow_dispatch:
+
+permissions:
+ contents: read
+
+jobs:
+ test:
+
+ runs-on: ubuntu-latest
+ strategy:
+ fail-fast: false
+ matrix:
+ python-version: ['3.10', '3.12']
+
+ steps:
+ - uses: actions/checkout@v4
+ - name: Set up Python ${{ matrix.python-version }}
+ uses: actions/setup-python@v5
+ with:
+ python-version: ${{ matrix.python-version }}
+ cache: pip
+ - name: Install dependencies
+ # Constrained so that a dependency release cannot turn this red without a commit
+ # here. The pins only apply to CI; setup.py stays open for anyone installing the
+ # package. See ci/constraints-py*.txt for how to move them forward.
+ run: |
+ python -m pip install --upgrade pip
+ pip install -e .[test] -c ci/constraints-py${{ matrix.python-version }}.txt
+ - name: Run tests
+ run: pytest -v --cov=cell2cell --cov-report=term-missing
+
+ latest-dependencies:
+ # The same suite with nothing pinned, to find out when the ecosystem moves ahead of
+ # the constraints. A failure here means a dependency needs attention, not that the
+ # commit is broken, which is why it does not run on pushes or pull requests.
+ if: github.event_name == 'schedule' || github.event_name == 'workflow_dispatch'
+
+ runs-on: ubuntu-latest
+ strategy:
+ fail-fast: false
+ matrix:
+ python-version: ['3.10', '3.12']
+
+ steps:
+ - uses: actions/checkout@v4
+ - name: Set up Python ${{ matrix.python-version }}
+ uses: actions/setup-python@v5
+ with:
+ python-version: ${{ matrix.python-version }}
+ cache: pip
+ - name: Install the newest dependencies
+ run: |
+ python -m pip install --upgrade pip
+ pip install -e .[test]
+ pip list
+ - name: Run tests
+ run: pytest -v
diff --git a/README.md b/README.md
index 36b4873..15be95a 100644
--- a/README.md
+++ b/README.md
@@ -44,7 +44,7 @@ pip install cell2cell
| cell2cell Examples | Tensor-cell2cell Examples |
| --- |-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
|  |  |
-| - [Step-by-step Pipeline](https://github.com/earmingol/cell2cell/blob/master/examples/cell2cell/Toy-Example.ipynb) - [Interaction Pipeline for Bulk Data](https://earmingol.github.io/cell2cell/tutorials/Toy-Example-BulkPipeline) - [Interaction Pipeline for Single-Cell Data](https://earmingol.github.io/cell2cell/tutorials/Toy-Example-SingleCellPipeline) - [Whole Body of *C. elegans*](https://github.com/LewisLabUCSD/Celegans-cell2cell) | - [Obtaining patterns of cell-cell communication](https://earmingol.github.io/cell2cell/tutorials/ASD/01-Tensor-Factorization-ASD/) - [Downstream 1: Factor-specific analyses](https://earmingol.github.io/cell2cell/tutorials/ASD/02-Factor-Specific-ASD/) - [Downstream 2: Patterns to functions (GSEA)](https://earmingol.github.io/cell2cell/tutorials/ASD/03-GSEA-ASD/) - [Tensor-cell2cell in Google Colab (**GPU**)](https://colab.research.google.com/drive/1T6MUoxafTHYhjvenDbEtQoveIlHT2U6_?usp=sharing) - [Communication patterns in **Spatial Transcriptomics**](https://earmingol.github.io/cell2cell/tutorials/Tensor-cell2cell-Spatial/) - [Multi-modal communication patterns with **Coupled Tensor Component Analysis**](https://earmingol.github.io/cell2cell/tutorials/Version2/Tensor-cell2cell-CTCA/) |
+| - [Step-by-step Pipeline](https://github.com/earmingol/cell2cell/blob/master/examples/cell2cell/Toy-Example.ipynb) - [Interaction Pipeline for Bulk Data](https://earmingol.github.io/cell2cell/tutorials/Toy-Example-BulkPipeline) - [Interaction Pipeline for Single-Cell Data](https://earmingol.github.io/cell2cell/tutorials/Toy-Example-SingleCellPipeline) - [Whole Body of *C. elegans*](https://github.com/LewisLabUCSD/Celegans-cell2cell) - [Identifying a **spatial code** of ligand-receptor pairs (Genetic Algorithm)](https://earmingol.github.io/cell2cell/tutorials/Genetic-Algorithm-LR-Selection) | - [Obtaining patterns of cell-cell communication](https://earmingol.github.io/cell2cell/tutorials/ASD/01-Tensor-Factorization-ASD/) - [Downstream 1: Factor-specific analyses](https://earmingol.github.io/cell2cell/tutorials/ASD/02-Factor-Specific-ASD/) - [Downstream 2: Patterns to functions (GSEA)](https://earmingol.github.io/cell2cell/tutorials/ASD/03-GSEA-ASD/) - [Tensor-cell2cell in Google Colab (**GPU**)](https://colab.research.google.com/drive/1T6MUoxafTHYhjvenDbEtQoveIlHT2U6_?usp=sharing) - [Communication patterns in **Spatial Transcriptomics**](https://earmingol.github.io/cell2cell/tutorials/Tensor-cell2cell-Spatial/) - [Multi-modal communication patterns with **Coupled Tensor Component Analysis**](https://earmingol.github.io/cell2cell/tutorials/Version2/Tensor-cell2cell-CTCA/) |
Reproducible runs of the analyses in the [Tensor-cell2cell paper](https://doi.org/10.1038/s41467-022-31369-2) are available at [CodeOcean.com](https://doi.org/10.24433/CO.0051950.v2)
diff --git a/cell2cell/__init__.py b/cell2cell/__init__.py
index 3cab73b..e534ceb 100644
--- a/cell2cell/__init__.py
+++ b/cell2cell/__init__.py
@@ -13,4 +13,4 @@
from cell2cell import tensor
from cell2cell import utils
-__version__ = "0.8.4"
+__version__ = "0.9.0"
diff --git a/cell2cell/analysis/__init__.py b/cell2cell/analysis/__init__.py
index d46a775..1b77e62 100644
--- a/cell2cell/analysis/__init__.py
+++ b/cell2cell/analysis/__init__.py
@@ -1,4 +1,8 @@
from cell2cell.analysis.cell2cell_pipelines import (initialize_interaction_space, BulkInteractions, SingleCellInteractions)
+from cell2cell.analysis.genetic_algorithm import (optimize_lr_pairs, PreparedCCIScorer,
+ lr_selection_frequency, lr_cooccurrence,
+ consensus_from_cooccurrence,
+ consensus_from_frequency)
from cell2cell.analysis.tensor_pipelines import (run_tensor_cell2cell_pipeline)
import cell2cell.analysis.tensor_downstream as tensor_downstream
diff --git a/cell2cell/analysis/cell2cell_pipelines.py b/cell2cell/analysis/cell2cell_pipelines.py
index ac36fed..1c54600 100644
--- a/cell2cell/analysis/cell2cell_pipelines.py
+++ b/cell2cell/analysis/cell2cell_pipelines.py
@@ -5,6 +5,7 @@
import pandas as pd
import scanpy
import numpy as np
+from natsort import natsorted
from tqdm import tqdm
from cell2cell.core import interaction_space as ispace
@@ -462,6 +463,11 @@ class SingleCellInteractions:
of a given gene.
- 'average' : Computes the average gene expression among the single cells
composing a cell type for a given gene.
+ - 'trimean' : Computes the Tukey's trimean of the gene expression among the
+ single cells composing a cell type for a given gene. It is a weighted
+ average of the median and the first and third quartiles
+ (0.5 * Q2 + 0.25 * (Q1 + Q3)), so it is more robust to outliers than
+ the average while still accounting for the spread of the distribution.
barcode_col : str, default='barcodes'
Column-name for the single cells in the metadata.
@@ -609,6 +615,11 @@ class SingleCellInteractions:
of a given gene.
- 'average' : Computes the average gene expression among the single cells
composing a cell type for a given gene.
+ - 'trimean' : Computes the Tukey's trimean of the gene expression among the
+ single cells composing a cell type for a given gene. It is a weighted
+ average of the median and the first and third quartiles
+ (0.5 * Q2 + 0.25 * (Q1 + Q3)), so it is more robust to outliers than
+ the average while still accounting for the spread of the distribution.
ccc_permutation_pvalues : pandas.DataFrame
Contains the P-values of the permutation analysis on the
@@ -935,7 +946,7 @@ def initialize_interaction_space(rnaseq_data, ppi_data, cutoff_setup, analysis_s
if excluded_cells is None:
excluded_cells = []
- included_cells = sorted(list((set(rnaseq_data.columns) - set(excluded_cells))))
+ included_cells = natsorted(set(rnaseq_data.columns) - set(excluded_cells))
interaction_space = ispace.InteractionSpace(rnaseq_data=rnaseq_data[included_cells],
ppi_data=ppi_data,
diff --git a/cell2cell/analysis/genetic_algorithm.py b/cell2cell/analysis/genetic_algorithm.py
new file mode 100644
index 0000000..c7cfbae
--- /dev/null
+++ b/cell2cell/analysis/genetic_algorithm.py
@@ -0,0 +1,1091 @@
+# -*- coding: utf-8 -*-
+
+'''Selection of ligand-receptor pairs with a genetic algorithm.
+
+This module reimplements, as part of the package, the analysis in
+https://github.com/LewisLabUCSD/Celegans-cell2cell (`code/genetic_algorithm.py`),
+which searches for the subset of ligand-receptor pairs whose cell-cell interaction
+scores best reproduce a reference distance between cells, for example the physical
+distances measured in a 3D map.
+
+Reference
+---------
+Armingol E, Ghaddar A, Joshi CJ, Baghdassarian H, Shamie I, Chan J, et al. (2022)
+Inferring a spatial code of cell-cell interactions across a whole animal body.
+PLOS Computational Biology 18(11): e1010715.
+https://doi.org/10.1371/journal.pcbi.1010715
+
+The objective function is the one used there: the absolute Spearman correlation
+between the CCI distance matrix and the reference distance matrix. The original
+implementation used `pyevolve`, which only supports Python 2; this one uses
+`pygad`, and evaluates the objective in a vectorized way that is orders of
+magnitude faster (see `optimize_lr_pairs` for the details).
+'''
+
+from __future__ import absolute_import
+
+import warnings
+
+import numpy as np
+import pandas as pd
+import scipy.spatial
+import scipy.stats
+
+from types import ModuleType
+
+from cell2cell.clustering.cluster_interactions import get_clusters_from_linkage
+from cell2cell.core.interaction_space import InteractionSpace
+from cell2cell.preprocessing.manipulate_dataframes import check_symmetry
+from cell2cell.preprocessing.ppi import bidirectional_ppi_for_cci, remove_ppi_bidirectionality
+
+
+# CCI scores whose value is a function of three quantities that are linear in the
+# PPI weights, which is what makes the vectorized objective possible. See
+# `PreparedCCIScorer` for the derivation.
+LINEAR_CCI_SCORES = ('bray_curtis', 'jaccard', 'count', 'icellnet')
+
+# Scores that are not bounded between 0 and 1, and whose distance matrix is
+# therefore computed with the regularized formula in `InteractionSpace`.
+UNBOUNDED_CCI_SCORES = ('count', 'icellnet')
+
+
+def _check_if_pygad() -> ModuleType:
+ try:
+ import pygad
+
+ except Exception:
+ raise ImportError('pygad is not installed. Please install it with: '
+ 'pip install pygad'
+ )
+ return pygad
+
+
+class PreparedCCIScorer:
+ '''
+ Precomputes the expression-dependent part of the CCI scores, so that scoring
+ the same cells again under a different set of PPI weights becomes a matrix
+ multiplication instead of a rebuild of the interaction space.
+
+ Every CCI score in `cell2cell.core.cci_scores` is a function of three
+ quantities, each of which is a weighted sum over the PPIs and therefore
+ **linear in the PPI weights** `w`:
+
+ .. code-block:: text
+
+ N(i, j) = sum_k w_k * A_ki * B_kj (the ligand-receptor product)
+ SA(i) = sum_k w_k * A_ki^2
+ SB(j) = sum_k w_k * B_kj^2
+
+ bray_curtis = 2 * N / (SA + SB)
+ jaccard = N / (SA + SB - N)
+ icellnet = N
+ count = sum_k [w_k * A_ki * B_kj != 0]
+
+ Since `A` and `B` depend only on the expression data, they can be computed
+ once and reused for any number of weight vectors. Scoring a whole population
+ of weight vectors is then a single matrix product.
+
+ Parameters
+ ----------
+ interaction_space : cell2cell.core.interaction_space.InteractionSpace
+ A built interaction space. Only its expression-derived matrices are used;
+ it is not modified.
+
+ cci_score : str, default=None
+ CCI score to compute. If None, the one the interaction space was built
+ with is used. Must be one of 'bray_curtis', 'jaccard', 'count' or
+ 'icellnet'.
+
+ max_memory_mb : float, default=512
+ Budget for the precomputed ligand-receptor outer products, which are what
+ make the batched path possible. If the array would exceed this, only the
+ per-vector path is prepared, which is still much faster than rebuilding
+ the interaction space but does not gain from batching.
+
+ Attributes
+ ----------
+ A, B : numpy.ndarray
+ Weighted ligand and receptor expression, of shape (PPIs, cells).
+
+ cell_names : list
+ Cell names, in the order they have in the interaction space. Rows and
+ columns of every returned matrix follow this order.
+
+ batched : boolean
+ Whether the precomputed outer products fitted in `max_memory_mb`.
+ '''
+
+ def __init__(self, interaction_space, cci_score=None, max_memory_mb=512):
+ if cci_score is None:
+ cci_score = interaction_space.cci_score
+ if cci_score not in LINEAR_CCI_SCORES:
+ raise NotImplementedError(
+ "'{}' is not supported by the vectorized scorer. Use one of {}, or "
+ "pass fast=False to fall back to the reference implementation."
+ .format(cci_score, list(LINEAR_CCI_SCORES)))
+
+ self.cci_score = cci_score
+ self.cci_type = interaction_space.cci_type
+ self.cell_names = list(interaction_space.interaction_elements['cell_names'])
+
+ cells = interaction_space.interaction_elements['cells']
+ self.A = np.column_stack([cells[c].weighted_ppi['A'].values for c in self.cell_names])
+ self.B = np.column_stack([cells[c].weighted_ppi['B'].values for c in self.cell_names])
+
+ # `nansum` in the scalar scores treats missing values as zero. Doing the
+ # substitution once here reproduces that without a NaN-aware reduction in
+ # the inner loop.
+ self._has_nans = bool(np.isnan(self.A).any() or np.isnan(self.B).any())
+ self.A = np.nan_to_num(self.A)
+ self.B = np.nan_to_num(self.B)
+
+ self.n_ppi, self.n_cells = self.A.shape
+ self._A2 = self.A * self.A
+ self._B2 = self.B * self.B
+
+ # Outer product of every ligand-receptor pair, flattened over the two cell
+ # axes so the contraction over PPIs is a plain matrix product.
+ outer_mb = self.n_ppi * self.n_cells * self.n_cells * 8 / 1e6
+ self.batched = outer_mb <= max_memory_mb
+ if self.batched:
+ outer = self.A[:, :, None] * self.B[:, None, :]
+ if self.cci_score == 'count':
+ # 'count' counts the non-zero products rather than adding them up,
+ # which is linear in a binary weight vector but not in a general one.
+ self._P = (outer != 0).astype(float).reshape(self.n_ppi, -1)
+ else:
+ self._P = outer.reshape(self.n_ppi, -1)
+ else:
+ self._P = None
+
+ def _terms(self, W):
+ '''Computes N, SA and SB for a stack of weight vectors W of shape (n, PPIs).'''
+ SA = W @ self._A2
+ SB = W @ self._B2
+ if self.batched:
+ N = (W @ self._P).reshape(-1, self.n_cells, self.n_cells)
+ else:
+ N = np.stack([(self.A * w[:, None]).T @ self.B for w in W])
+ return N, SA, SB
+
+ def _combine(self, N, SA, SB):
+ '''Applies the score-specific formula. Shapes: N (n, C, C), SA/SB (n, C).'''
+ if self.cci_score == 'icellnet' or self.cci_score == 'count':
+ return N
+
+ denominator = SA[:, :, None] + SB[:, None, :]
+ if self.cci_score == 'jaccard':
+ denominator = denominator - N
+
+ with np.errstate(divide='ignore', invalid='ignore'):
+ if self.cci_score == 'bray_curtis':
+ scores = np.divide(2.0 * N, denominator)
+ else:
+ scores = np.divide(N, denominator)
+ # The scalar implementations return 0.0 when the denominator is zero
+ scores[denominator == 0.0] = 0.0
+ return scores
+
+ def score_batch(self, W):
+ '''
+ Computes the CCI matrix for each of several PPI weight vectors.
+
+ Parameters
+ ----------
+ W : numpy.ndarray
+ Weights, of shape (n, PPIs). One row per weight vector.
+
+ Returns
+ -------
+ scores : numpy.ndarray
+ CCI matrices, of shape (n, cells, cells). Rows and columns follow
+ `self.cell_names`. For an undirected interaction space the matrices
+ are symmetrized the same way `compute_pairwise_cci_scores` does, by
+ mirroring the upper triangle.
+ '''
+ W = np.atleast_2d(np.asarray(W, dtype=float))
+ if W.shape[1] != self.n_ppi:
+ raise ValueError('Expected weight vectors of length {}, got {}'
+ .format(self.n_ppi, W.shape[1]))
+
+ if self.batched and self.cci_score == 'count' and not np.isin(W, (0.0, 1.0)).all():
+ raise ValueError("The 'count' score is only vectorized for binary weights, "
+ "because it counts non-zero products rather than adding them.")
+
+ scores = self._combine(*self._terms(W))
+
+ if self.cci_type == 'undirected':
+ # `generate_pairs` yields the upper triangle plus the diagonal, and the
+ # scoring loop mirrors each value. The lower triangle of the directed
+ # result is therefore never used.
+ upper = np.triu(scores)
+ scores = upper + np.triu(scores, k=1).transpose(0, 2, 1)
+ return scores
+
+ def score(self, ppi_score):
+ '''
+ Computes the CCI matrix for a single PPI weight vector.
+
+ Parameters
+ ----------
+ ppi_score : array-like
+ Weights, one per PPI.
+
+ Returns
+ -------
+ cci_matrix : pandas.DataFrame
+ CCI scores, with cells as rows and columns.
+ '''
+ scores = self.score_batch(np.asarray(ppi_score, dtype=float)[None, :])[0]
+ return pd.DataFrame(scores, index=self.cell_names, columns=self.cell_names)
+
+ def distance_batch(self, W):
+ '''
+ Computes the distance matrix for each of several PPI weight vectors,
+ reproducing what `InteractionSpace.compute_pairwise_cci_scores` derives.
+
+ Bounded scores use `1 - score`; the unbounded ones ('count', 'icellnet')
+ use the regularized `1 - score / (score + mean)`, where the mean is taken
+ over the whole CCI matrix of that weight vector. The diagonal is zeroed.
+
+ Parameters
+ ----------
+ W : numpy.ndarray
+ Weights, of shape (n, PPIs).
+
+ Returns
+ -------
+ distances : numpy.ndarray
+ Distance matrices, of shape (n, cells, cells).
+ '''
+ scores = self.score_batch(W)
+ if self.cci_score in UNBOUNDED_CCI_SCORES:
+ means = np.nanmean(scores, axis=(1, 2))[:, None, None]
+ with np.errstate(divide='ignore', invalid='ignore'):
+ distances = 1.0 - np.divide(scores, scores + means)
+ else:
+ distances = 1.0 - scores
+
+ idx = np.arange(self.n_cells)
+ distances[:, idx, idx] = 0.0
+ return distances
+
+
+def _reference_distance_matrix(interaction_space, ppi_score, cells):
+ '''Distance matrix through the unmodified `InteractionSpace` code path.'''
+ interaction_space.ppi_data['score'] = np.asarray(ppi_score, dtype=float)
+ interaction_space.interaction_elements['ppi_score'] = interaction_space.ppi_data['score'].values
+ interaction_space.compute_pairwise_cci_scores(use_ppi_score=True, verbose=False)
+ return interaction_space.distance_matrix.loc[cells, cells]
+
+
+def _bidirectional_index(ppi_data, interaction_columns=('A', 'B'), verbose=False):
+ '''
+ Maps every row of the bidirectional PPI table back to the row of `ppi_data`
+ it came from.
+
+ `bidirectional_ppi_for_cci` duplicates every interaction with its partners
+ swapped and then drops duplicates, which collapses self-interactions back to a
+ single copy. Rather than reimplementing that, this runs it once on a table
+ whose score column holds each row's position, and reads the positions back.
+
+ Returns
+ -------
+ source : numpy.ndarray
+ For each row of the bidirectional table, the index of the row of
+ `ppi_data` it originates from.
+ '''
+ probe = ppi_data.copy()
+ probe['score'] = np.arange(len(probe), dtype=float)
+ bi_probe = bidirectional_ppi_for_cci(ppi_data=probe,
+ interaction_columns=interaction_columns,
+ verbose=verbose)
+ source = bi_probe['score'].values.astype(int)
+
+ # `drop_duplicates` acts on (A, B, score), so if the table still contains a pair
+ # and its reciprocal, two rows that differ only by score stop being duplicates and
+ # the bidirectional table changes length with the weights. The mapping is then not
+ # well defined -- and neither is assigning that column to a fixed interaction space.
+ # The interaction space is built from the all-ones table, so that one has to
+ # match too, not just an arbitrary binary vector.
+ rng = np.random.default_rng(0)
+ probes = [np.ones(len(ppi_data)),
+ (np.arange(len(ppi_data)) % 2).astype(float),
+ rng.integers(0, 2, size=len(ppi_data)).astype(float)]
+ lengths = {len(bidirectional_ppi_for_cci(ppi_data=ppi_data.assign(score=p),
+ interaction_columns=interaction_columns,
+ verbose=verbose))
+ for p in probes}
+ if lengths != {len(source)}:
+ raise ValueError(
+ 'The number of bidirectional interactions depends on the weights, so a '
+ 'ligand-receptor pair cannot be mapped onto a fixed set of rows. This '
+ 'happens when `ppi_data` holds a pair and its reciprocal as separate rows, '
+ 'or the exact same pair more than once. Deduplicate it first -- '
+ 'cell2cell.preprocessing.remove_ppi_bidirectionality() followed by '
+ 'drop_duplicates() on the interaction columns.')
+ return source
+
+
+def _correlation(distance_vector, reference_vector, method='spearman'):
+ if method == 'spearman':
+ corr = scipy.stats.spearmanr(distance_vector, reference_vector)[0]
+ elif method == 'pearson':
+ corr = scipy.stats.pearsonr(distance_vector, reference_vector)[0]
+ else:
+ raise ValueError("`method` must be either 'spearman' or 'pearson'")
+ return abs(np.nan_to_num(corr))
+
+
+def lr_selection_frequency(selection_masks):
+ '''
+ Fraction of independent genetic-algorithm executions that selected each pair.
+
+ Parameters
+ ----------
+ selection_masks : array-like
+ Binary matrix of shape (executions, LR pairs). One row per independent run
+ of the genetic algorithm, holding the 0/1 mask it converged to.
+
+ Returns
+ -------
+ frequency : numpy.ndarray
+ Value in [0, 1] per ligand-receptor pair.
+ '''
+ masks = np.atleast_2d(np.asarray(selection_masks, dtype=float))
+ return np.nansum(masks, axis=0) / masks.shape[0]
+
+
+def lr_cooccurrence(selection_masks, labels=None):
+ '''
+ Co-occurrence of ligand-receptor pairs across independent genetic-algorithm runs.
+
+ Two pairs co-occur when the same run selected both. The value reported is the
+ Jaccard index of their selection patterns -- the number of runs that selected
+ both, divided by the number that selected either -- so a pair that is chosen
+ rarely can still co-occur strongly with another it is always chosen alongside.
+
+ Parameters
+ ----------
+ selection_masks : array-like
+ Binary matrix of shape (executions, LR pairs), one row per independent run.
+
+ labels : list, default=None
+ Names for the ligand-receptor pairs, used as the index and columns of the
+ result. If None, positional integers are used.
+
+ Returns
+ -------
+ cooccurrence : pandas.DataFrame
+ Symmetric matrix of Jaccard indexes, with ones on the diagonal for pairs
+ selected at least once and zeros for pairs never selected.
+ '''
+ masks = np.atleast_2d(np.asarray(selection_masks)).astype(bool)
+ intersection = (masks.astype(float).T @ masks.astype(float))
+ counts = masks.sum(axis=0)
+ union = counts[:, None] + counts[None, :] - intersection
+
+ with np.errstate(divide='ignore', invalid='ignore'):
+ cooccurrence = np.divide(intersection, union)
+ # A pair of LRs that no run ever selected has an empty union
+ cooccurrence[union == 0] = 0.0
+
+ if labels is None:
+ labels = list(range(masks.shape[1]))
+ return pd.DataFrame(cooccurrence, index=labels, columns=labels)
+
+
+def consensus_from_cooccurrence(cooccurrence, n_clusters=2, method='ward',
+ select='cooccurrence', frequency=None, min_frequency=0.0):
+ '''
+ Picks the group of ligand-receptor pairs that keep being selected together.
+
+ Clusters the co-occurrence matrix and returns one cluster: the tight group of
+ pairs that are chosen alongside each other across independent runs of the
+ genetic algorithm. This is how the published selection for *C. elegans* was
+ produced.
+
+ Parameters
+ ----------
+ cooccurrence : pandas.DataFrame
+ Square co-occurrence matrix, as returned by `lr_cooccurrence`. Pairs that no
+ run ever selected (an all-zero row and column) are dropped first.
+
+ n_clusters : int, default=2
+ Number of clusters to cut the dendrogram into.
+
+ method : str, default='ward'
+ Linkage method.
+
+ frequency : array-like, default=None
+ Fraction of runs that selected each pair, aligned with `cooccurrence`. Only
+ needed for `min_frequency`.
+
+ min_frequency : float, default=0.0
+ Drop pairs selected in fewer than this fraction of runs before clustering.
+ Two pairs selected once, in the same run, have a Jaccard index of 1.0 even
+ though a single run is no evidence that they belong together, and with few
+ runs a cluster of such pairs can outscore the genuinely reproducible one.
+ This removes them, but it is not a substitute for running the search enough
+ times: **the clustering needs on the order of 30 or more runs to be stable**,
+ and below that the frequency route is the more reliable of the two. Requires
+ `frequency`.
+
+ select : str, default='cooccurrence'
+ Which cluster to return.
+
+ - 'cooccurrence' : the one with the highest mean co-occurrence among its own
+ members, i.e. the most consistently co-selected group. This is the
+ intent of the analysis and does not depend on cluster sizes.
+ - 'smallest' : the one with fewest members. This is literally what the
+ reference notebook did, and on that data it is also the highest-
+ co-occurrence one; it is kept for exact reproducibility.
+
+ Returns
+ -------
+ selected : list
+ Labels of the pairs in the chosen cluster.
+
+ clusters : dict
+ Every cluster, keyed by its scipy cluster id, so the others can be inspected.
+
+ scores : dict
+ Mean intra-cluster co-occurrence per cluster id, which is what `select`
+ ranks on.
+ '''
+ import scipy.cluster.hierarchy as hc
+ from sklearn.metrics import pairwise_distances
+
+ keep = (cooccurrence != 0).any(axis=0)
+ if min_frequency > 0.0:
+ if frequency is None:
+ raise ValueError('`frequency` is required when `min_frequency` is set')
+ frequency = np.asarray(frequency, dtype=float)
+ if len(frequency) != cooccurrence.shape[0]:
+ raise ValueError('`frequency` must have one value per pair in `cooccurrence`')
+ keep = keep & (frequency >= min_frequency)
+ data = cooccurrence.loc[keep, keep]
+ if data.shape[0] < n_clusters:
+ raise ValueError('Only {} pairs passed the filters, which is fewer than the {} '
+ 'clusters requested. Lower `min_frequency`, or run more '
+ 'executions.'.format(data.shape[0], n_clusters))
+
+ # Distances between co-occurrence profiles, then Ward on those. `hc.linkage`
+ # reads a square array as observations x features, so it clusters the rows of
+ # the distance matrix -- which is what the reference analysis did.
+ distances = pairwise_distances(data.values)
+ with warnings.catch_warnings():
+ # scipy notices that a square hollow matrix was passed where it usually takes
+ # a condensed one, and warns. That is what is meant here: the rows of the
+ # distance matrix are the observations, matching the reference analysis.
+ warnings.simplefilter('ignore', hc.ClusterWarning)
+ linkage = hc.linkage(distances, method=method, optimal_ordering=True)
+
+ clusters = get_clusters_from_linkage(linkage, n_clusters, criterion='maxclust',
+ labels=list(data.index))
+
+ # Mean co-occurrence between distinct members of each cluster
+ scores = {}
+ for key, members in clusters.items():
+ if len(members) < 2:
+ scores[key] = 0.0
+ continue
+ block = data.loc[members, members].values
+ off_diagonal = block[~np.eye(len(members), dtype=bool)]
+ scores[key] = float(np.nanmean(off_diagonal))
+
+ if select == 'cooccurrence':
+ chosen = max(scores, key=lambda k: scores[k])
+ elif select == 'smallest':
+ chosen = min(clusters, key=lambda k: len(clusters[k]))
+ else:
+ raise ValueError("`select` must be either 'cooccurrence' or 'smallest'")
+ return clusters[chosen], clusters, scores
+
+
+def consensus_from_frequency(frequency, percentile=90):
+ '''
+ Keeps the ligand-receptor pairs selected most often across independent runs.
+
+ The simpler alternative to `consensus_from_cooccurrence`: rather than asking
+ which pairs are chosen *together*, it asks which are chosen *often*. Cheaper to
+ reason about, but it cannot separate two groups of pairs that are each
+ self-consistent yet rarely co-selected.
+
+ Parameters
+ ----------
+ frequency : array-like
+ Fraction of runs that selected each pair, from `lr_selection_frequency`.
+
+ percentile : float, default=90
+ Percentile of the frequency distribution used as the cutoff. Pairs strictly
+ above it are kept. The reference analysis used the 90th percentile.
+
+ Returns
+ -------
+ mask : numpy.ndarray
+ Boolean, True for the pairs that are kept.
+
+ threshold : float
+ The cutoff value.
+ '''
+ values = np.asarray(frequency, dtype=float)
+ threshold = float(np.percentile(values, percentile))
+ return values > threshold, threshold
+
+
+
+def _optimize_once(rnaseq_data, ppi_data, reference_distances, cutoff_setup, analysis_setup,
+ included_cells=None, population_size=200, generations=200, runs=None,
+ inc_percentage=0.025, max_runs=100, correlation='spearman',
+ mutation_probability=0.05, keep_elitism=1, random_state=None,
+ interaction_columns=('A', 'B'), complex_sep=None, complex_agg_method='min',
+ fast=True, validate_fast=True, max_memory_mb=512, deduplicate=True,
+ verbose=False):
+ '''
+ Selects the subset of ligand-receptor pairs whose cell-cell interaction scores
+ best reproduce a reference distance between cells, using a genetic algorithm.
+
+ Each individual is a binary vector with one entry per ligand-receptor pair,
+ indicating whether it is included. The objective function is the absolute
+ Spearman correlation between the resulting CCI distance matrix and
+ `reference_distances`, as in Armingol et al. (2022) on the whole body of *C. elegans*.
+
+ The search is repeated in successive runs: each run keeps only the pairs that
+ the previous one selected, so the set shrinks until the objective stops
+ improving by at least `inc_percentage`.
+
+ Parameters
+ ----------
+ rnaseq_data : pandas.DataFrame
+ Gene expression matrix, with genes as rows and cells as columns.
+
+ ppi_data : pandas.DataFrame
+ List of ligand-receptor pairs. A 'score' column is added if missing.
+
+ reference_distances : pandas.DataFrame
+ Square, symmetric matrix of reference distances between cells, for example
+ physical distances. Rows and columns are cell names.
+
+ cutoff_setup : dict
+ Cutoff setup, as in `cell2cell.analysis.initialize_interaction_space`.
+
+ analysis_setup : dict
+ Analysis setup with the keys 'communication_score', 'cci_score' and
+ 'cci_type', as in `cell2cell.analysis.initialize_interaction_space`.
+ `cci_type` must be 'undirected', since the objective compares the
+ condensed form of a symmetric distance matrix.
+
+ included_cells : list, default=None
+ Cells to consider. If None, the cells present in both `rnaseq_data` and
+ `reference_distances` are used.
+
+ population_size : int, default=200
+ Number of individuals per generation.
+
+ generations : int, default=200
+ Number of generations per run.
+
+ runs : int, default=None
+ Number of runs. If None, runs continue until the objective improves by
+ less than `inc_percentage` with respect to the previous run, or until
+ `max_runs` is reached.
+
+ inc_percentage : float, default=0.025
+ Minimum relative improvement of the objective for another run to start.
+ Only used when `runs` is None.
+
+ max_runs : int, default=100
+ Upper bound on the number of runs when `runs` is None.
+
+ correlation : str, default='spearman'
+ Correlation between the CCI distances and the reference distances, either
+ 'spearman' or 'pearson'.
+
+ mutation_probability : float, default=0.05
+ Probability of flipping each gene, equivalent to the flip mutator of the
+ original implementation.
+
+ keep_elitism : int, default=1
+ Number of best individuals carried over to the next generation.
+
+ random_state : int, default=None
+ Seed for reproducibility.
+
+ interaction_columns : tuple, default=('A', 'B')
+ Columns of `ppi_data` holding the ligands and the receptors.
+
+ complex_sep : str, default=None
+ Separator of the subunits of a protein complex, if any.
+
+ complex_agg_method : str, default='min'
+ Method to aggregate the expression of the subunits of a complex.
+
+ fast : boolean, default=True
+ Whether to evaluate the objective with `PreparedCCIScorer`, which is
+ equivalent but vectorized. If False, every individual is evaluated by
+ rebuilding the CCI scores through `InteractionSpace`, exactly as the
+ original implementation did.
+
+ validate_fast : boolean, default=True
+ Whether to check the vectorized objective against the reference one on a
+ few random individuals before starting. Cheap, and it catches the case
+ where the two paths would disagree.
+
+ max_memory_mb : float, default=512
+ Memory budget for the precomputed ligand-receptor outer products.
+
+ deduplicate : boolean, default=True
+ Whether to collapse each interaction and its reciprocal into a single row
+ with `remove_ppi_bidirectionality` before the search. This is required for
+ a pair to map onto a fixed set of bidirectional rows: `bidirectional_ppi_for_cci`
+ drops duplicates on (A, B, score), so if both directions are present as
+ separate rows, the number of bidirectional rows depends on the candidate
+ solution. Pairs loaded with `cell2cell.io.load_ppi` are already deduplicated
+ by `preprocess_ppi_data`, so this is a no-op for them; it matters when the
+ table comes from somewhere else. Turn it off only if the input is already
+ deduplicated -- the ambiguous case is rejected either way.
+ Note that the returned masks are then indexed against the deduplicated
+ table, which is also what 'best_ppi_data' contains.
+
+ verbose : boolean, default=False
+ Whether to print the progress of each run.
+
+ Returns
+ -------
+ results : dict
+ Dictionary with one entry per run, keyed 'run1', 'run2', ..., each with:
+
+ - 'obj_fn' : the objective function of the best individual.
+ - 'ppi_data' : list of 0/1, one per row of the original `ppi_data`,
+ indicating the pairs selected in that run.
+ - 'drop_fraction' : fraction of the pairs available to that run that were
+ dropped.
+ - 'n_selected' : number of pairs selected.
+
+ The dictionary also holds 'best_run', 'best_obj_fn' and 'best_ppi_data',
+ the last being a copy of `ppi_data` restricted to the selected pairs.
+
+ Examples
+ --------
+ >>> import cell2cell as c2c
+ >>> results = c2c.analysis.optimize_lr_pairs(rnaseq_data=rnaseq,
+ ... ppi_data=lr_pairs,
+ ... reference_distances=physical_distances,
+ ... cutoff_setup={'type': 'constant_value',
+ ... 'parameter': 10},
+ ... analysis_setup={'communication_score': 'expression_thresholding',
+ ... 'cci_score': 'bray_curtis',
+ ... 'cci_type': 'undirected'},
+ ... random_state=888)
+ >>> selected = results['best_ppi_data']
+ '''
+ pygad = _check_if_pygad()
+
+ if analysis_setup['cci_type'] != 'undirected':
+ raise NotImplementedError("Only 'undirected' interactions are supported, because the "
+ "objective compares condensed symmetric distance matrices.")
+
+ reference_distances = _as_symmetric(reference_distances)
+
+ # Cells shared by the expression data and the reference distances
+ if included_cells is None:
+ included_cells = sorted(set(rnaseq_data.columns) & set(reference_distances.columns))
+ included_cells = list(included_cells)
+ if len(included_cells) < 3:
+ raise ValueError('At least three cells are needed to correlate distances')
+
+ reference_vector = scipy.spatial.distance.squareform(
+ np.asarray(reference_distances.loc[included_cells, included_cells].values, dtype=float),
+ checks=False)
+
+ if deduplicate:
+ # Required for the pair-to-bidirectional-row mapping to be well defined; see
+ # `_bidirectional_index`. Also what the reference analysis effectively had,
+ # since its LR list held each interaction once.
+ ppi_data = remove_ppi_bidirectionality(ppi_data=ppi_data,
+ interaction_columns=interaction_columns,
+ verbose=verbose)
+ ppi_data = ppi_data.drop_duplicates(subset=list(interaction_columns))
+ ppi_data = ppi_data.reset_index(drop=True)
+
+ theta_ppi_data = ppi_data.copy()
+ if 'score' not in theta_ppi_data.columns:
+ theta_ppi_data = theta_ppi_data.assign(score=1.0)
+
+ prot_a, prot_b = interaction_columns
+ results = dict()
+ run = 1
+ previous_obj = None
+
+ while True:
+ if runs is None:
+ if run > max_runs:
+ break
+ elif run > runs:
+ break
+
+ # Each run searches only among the pairs the previous run kept
+ theta_ppi_data = theta_ppi_data.loc[theta_ppi_data['score'] == 1].reset_index(drop=True)
+ n_ppi = len(theta_ppi_data)
+ if n_ppi == 0:
+ break
+
+ bi_ppi_data = bidirectional_ppi_for_cci(ppi_data=theta_ppi_data,
+ interaction_columns=interaction_columns,
+ verbose=verbose)
+ interaction_space = InteractionSpace(rnaseq_data=rnaseq_data[included_cells],
+ ppi_data=bi_ppi_data,
+ gene_cutoffs=cutoff_setup,
+ communication_score=analysis_setup['communication_score'],
+ cci_score=analysis_setup['cci_score'],
+ cci_type=analysis_setup['cci_type'],
+ complex_sep=complex_sep,
+ complex_agg_method=complex_agg_method,
+ interaction_columns=interaction_columns,
+ verbose=verbose)
+
+ # Position in `theta_ppi_data` that each bidirectional row comes from, so a
+ # candidate solution can be expanded to the weights the scorer expects
+ source = _bidirectional_index(theta_ppi_data,
+ interaction_columns=interaction_columns,
+ verbose=verbose)
+
+ space_cells = list(interaction_space.interaction_elements['cell_names'])
+ take = [space_cells.index(c) for c in included_cells]
+
+ def reference_objective(theta):
+ weights = np.asarray(theta, dtype=float)[source]
+ distances = _reference_distance_matrix(interaction_space, weights, included_cells)
+ vector = scipy.spatial.distance.squareform(np.asarray(distances.values, dtype=float),
+ checks=False)
+ return _correlation(vector, reference_vector, method=correlation)
+
+ use_fast = fast
+ if use_fast:
+ try:
+ scorer = PreparedCCIScorer(interaction_space,
+ cci_score=analysis_setup['cci_score'],
+ max_memory_mb=max_memory_mb)
+ except NotImplementedError:
+ if verbose:
+ print('Falling back to the reference objective for this CCI score')
+ use_fast = False
+
+ if use_fast and scorer._has_nans and analysis_setup['cci_score'] == 'count':
+ # `count` treats a NaN product as active, which the substitution above
+ # does not reproduce. Only this score is affected.
+ use_fast = False
+
+ if use_fast:
+ def batch_objective(THETA):
+ W = np.asarray(THETA, dtype=float)[:, source]
+ distances = scorer.distance_batch(W)[:, take][:, :, take]
+ out = np.empty(len(W))
+ for n, d in enumerate(distances):
+ vector = scipy.spatial.distance.squareform(d, checks=False)
+ out[n] = _correlation(vector, reference_vector, method=correlation)
+ return out
+
+ if validate_fast:
+ rng = np.random.default_rng(random_state)
+ probes = rng.integers(0, 2, size=(2, n_ppi)).astype(float)
+ fast_values = batch_objective(probes)
+ ref_values = np.array([reference_objective(p) for p in probes])
+ if not np.allclose(fast_values, ref_values, rtol=1e-9, atol=1e-9):
+ raise RuntimeError(
+ 'The vectorized objective disagrees with the reference one '
+ '({} vs {}). Please report this, and use fast=False meanwhile.'
+ .format(fast_values, ref_values))
+
+ def fitness_func(ga_instance, solution, solution_idx):
+ return float(batch_objective(np.atleast_2d(solution))[0])
+ else:
+ def fitness_func(ga_instance, solution, solution_idx):
+ return float(reference_objective(solution))
+
+ ga = pygad.GA(num_generations=generations,
+ num_parents_mating=max(2, population_size // 2),
+ fitness_func=fitness_func,
+ sol_per_pop=population_size,
+ num_genes=n_ppi,
+ gene_type=int,
+ init_range_low=0,
+ init_range_high=2,
+ gene_space=[0, 1],
+ parent_selection_type='tournament',
+ keep_elitism=keep_elitism,
+ mutation_type='random',
+ mutation_probability=mutation_probability,
+ random_seed=random_state if random_state is None else random_state + run,
+ suppress_warnings=True,
+ )
+ if use_fast:
+ # Evaluate the whole generation with one matrix product
+ ga.fitness_batch_size = population_size
+
+ def batch_fitness(ga_instance, solutions, solutions_indices):
+ return list(batch_objective(np.atleast_2d(solutions)))
+
+ ga.fitness_func = batch_fitness
+
+ ga.run()
+
+ best_solution, best_fitness, _ = ga.best_solution()
+ best = np.asarray(best_solution, dtype=int)
+
+ theta_ppi_data['score'] = best.astype(float)
+ drop_fraction = 1.0 - best.sum() / len(best)
+
+ # Map the selection back onto the rows of the original ppi_data
+ selected = theta_ppi_data.loc[theta_ppi_data['score'] == 1, [prot_a, prot_b]]
+ selected_pairs = set(map(tuple, selected.values))
+ mask = [1 if tuple(row) in selected_pairs else 0
+ for row in ppi_data[[prot_a, prot_b]].values]
+
+ results['run{}'.format(run)] = {'obj_fn': float(best_fitness),
+ 'ppi_data': mask,
+ 'drop_fraction': float(drop_fraction),
+ 'n_selected': int(best.sum()),
+ }
+ if verbose:
+ print('Run {}: objective {:.4f}, {} of {} pairs kept'
+ .format(run, best_fitness, int(best.sum()), len(best)))
+
+ if runs is None and previous_obj is not None:
+ if (best_fitness - previous_obj) / previous_obj < inc_percentage:
+ run += 1
+ break
+ previous_obj = best_fitness
+ run += 1
+
+ if not results:
+ raise RuntimeError('The genetic algorithm produced no results')
+
+ best_key = max(results, key=lambda k: results[k]['obj_fn'])
+ best_mask = np.asarray(results[best_key]['ppi_data'], dtype=bool)
+ results['best_run'] = best_key
+ results['best_obj_fn'] = results[best_key]['obj_fn']
+ results['best_ppi_data'] = ppi_data.loc[best_mask].reset_index(drop=True)
+ return results
+
+def optimize_lr_pairs(rnaseq_data, ppi_data, reference_distances, cutoff_setup, analysis_setup,
+ executions=1, random_state=None, consensus_method='cooccurrence',
+ n_clusters=2, cluster_selection='cooccurrence', min_frequency=0.0,
+ frequency_percentile=90, verbose=False, **kwargs):
+ '''
+ Selects ligand-receptor pairs whose cell-cell interaction scores best reproduce a
+ reference distance between cells, using a genetic algorithm.
+
+ A genetic algorithm converges to a *local* optimum, so a single execution is not
+ conclusive: different seeds settle on different, largely overlapping sets of
+ pairs. With `executions > 1` the search is repeated independently and the results
+ are integrated the way the reference analysis did -- by how often each pair is
+ selected, and by which pairs are selected *together* -- which separates the
+ reproducible core from the noise of any one execution.
+
+ Parameters
+ ----------
+ rnaseq_data : pandas.DataFrame
+ Gene expression matrix, with genes as rows and cells as columns.
+
+ ppi_data : pandas.DataFrame
+ List of ligand-receptor pairs. A 'score' column is added if missing.
+
+ reference_distances : pandas.DataFrame
+ Square, symmetric matrix of reference distances between cells, for example
+ physical distances. Rows and columns are cell names.
+
+ cutoff_setup : dict
+ Cutoff setup, as in `cell2cell.analysis.initialize_interaction_space`.
+
+ analysis_setup : dict
+ Analysis setup with the keys 'communication_score', 'cci_score' and
+ 'cci_type'. `cci_type` must be 'undirected'.
+
+ executions : int, default=1
+ Number of independent runs of the genetic algorithm. Each uses a different
+ seed derived from `random_state`. With more than one, the consensus outputs
+ described below are added to the result. **Around 30 or more is needed for
+ the co-occurrence clustering to be stable**; the reference analysis used
+ about a hundred. The frequency route tolerates fewer.
+
+ random_state : int, default=None
+ Seed. Execution *i* uses `random_state + i`, so the whole set is reproducible.
+
+ consensus_method : str, default='cooccurrence'
+ How to integrate the executions when there is more than one.
+
+ - 'cooccurrence' : cluster the pairs by how often they are selected
+ *together* and keep one cluster. This is what the reference analysis
+ did, and it is able to tell apart groups of pairs that are each
+ self-consistent.
+ - 'frequency' : keep the pairs selected most often, above
+ `frequency_percentile`. Simpler, and blind to which pairs go together.
+
+ n_clusters : int, default=2
+ Number of clusters to cut the co-occurrence dendrogram into.
+
+ cluster_selection : str, default='cooccurrence'
+ Which cluster to keep: 'cooccurrence' for the one whose members co-occur
+ most with each other, or 'smallest' for the fewest members, which is
+ literally what the reference notebook did. See `consensus_from_cooccurrence`.
+
+ min_frequency : float, default=0.0
+ Drop pairs selected in fewer than this fraction of executions before building
+ the co-occurrence clusters, removing pairs that co-occur perfectly only
+ because they were each chosen once, in the same execution. Note the
+ co-occurrence route needs roughly 30 or more executions to be stable; with
+ fewer, prefer `consensus_method='frequency'`. See
+ `consensus_from_cooccurrence`.
+
+ frequency_percentile : float, default=90
+ Cutoff percentile when `consensus_method='frequency'`.
+
+ verbose : boolean, default=False
+ Whether to print the progress of each execution.
+
+ **kwargs
+ Passed to each individual search: `population_size`, `generations`, `runs`,
+ `inc_percentage`, `max_runs`, `correlation`, `mutation_probability`,
+ `keep_elitism`, `included_cells`, `interaction_columns`, `complex_sep`,
+ `complex_agg_method`, `fast`, `validate_fast`, `max_memory_mb` and
+ `deduplicate`. See `_optimize_once` for their meaning.
+
+ Returns
+ -------
+ results : dict
+ With a single execution, the result of that search: one entry per run
+ ('run1', 'run2', ...) with 'obj_fn', 'ppi_data' (a 0/1 mask over the rows of
+ `ppi_data`), 'drop_fraction' and 'n_selected', plus 'best_run',
+ 'best_obj_fn' and 'best_ppi_data'.
+
+ With several executions, the same keys report the single best execution, and
+ these are added:
+
+ - 'executions' : the full result of each execution, keyed 'execution1', ...
+ - 'selection_masks' : binary array of shape (executions, LR pairs), the mask
+ each execution converged to.
+ - 'selection_frequency' : dataframe of the pairs with the fraction of
+ executions that selected each one.
+ - 'cooccurrence' : Jaccard co-occurrence between pairs across executions.
+ - 'consensus_ppi_data' : the consensus selection -- with the default method,
+ the pairs of the co-occurrence cluster whose members are most consistently
+ chosen together. **This is the recommended output.**
+ - 'consensus_clusters' : every cluster, so the others can be inspected.
+ - 'consensus_cluster_scores' : mean intra-cluster co-occurrence per cluster.
+
+ With `consensus_method='frequency'`, 'consensus_threshold' holds the
+ frequency cutoff instead of the cluster keys.
+
+ Examples
+ --------
+ >>> import cell2cell as c2c
+ >>> results = c2c.analysis.optimize_lr_pairs(rnaseq_data=rnaseq,
+ ... ppi_data=lr_pairs,
+ ... reference_distances=physical_distances,
+ ... cutoff_setup={'type': 'constant_value',
+ ... 'parameter': 10},
+ ... analysis_setup={'communication_score': 'expression_thresholding',
+ ... 'cci_score': 'bray_curtis',
+ ... 'cci_type': 'undirected'},
+ ... executions=20, random_state=888)
+ >>> results['consensus_ppi_data']
+ '''
+ if executions < 1:
+ raise ValueError('`executions` must be at least 1')
+
+ common = dict(rnaseq_data=rnaseq_data, ppi_data=ppi_data,
+ reference_distances=reference_distances, cutoff_setup=cutoff_setup,
+ analysis_setup=analysis_setup, verbose=verbose, **kwargs)
+
+ if executions == 1:
+ return _optimize_once(random_state=random_state, **common)
+
+ interaction_columns = kwargs.get('interaction_columns', ('A', 'B'))
+ prot_a, prot_b = interaction_columns
+
+ # The pool the masks are indexed against, matching what each execution searches
+ pool = ppi_data
+ if kwargs.get('deduplicate', True):
+ pool = remove_ppi_bidirectionality(ppi_data=ppi_data,
+ interaction_columns=interaction_columns,
+ verbose=False)
+ pool = pool.drop_duplicates(subset=list(interaction_columns)).reset_index(drop=True)
+
+ all_executions, masks = {}, []
+ for i in range(executions):
+ seed = None if random_state is None else random_state + i
+ result = _optimize_once(random_state=seed, **common)
+ all_executions['execution{}'.format(i + 1)] = result
+ # The pairs that execution converged to, i.e. its last run
+ last = max((k for k in result if k.startswith('run')),
+ key=lambda k: int(k[3:]))
+ masks.append(result[last]['ppi_data'])
+ if verbose:
+ print('Execution {}: best objective {:.4f}, {} pairs'
+ .format(i + 1, result['best_obj_fn'], sum(result[last]['ppi_data'])))
+
+ masks = np.asarray(masks, dtype=int)
+ labels = ['{}^{}'.format(a, b) for a, b in pool[[prot_a, prot_b]].values]
+
+ frequency = pd.DataFrame({prot_a: pool[prot_a].values, prot_b: pool[prot_b].values,
+ 'frequency': lr_selection_frequency(masks)})
+ cooccurrence = lr_cooccurrence(masks, labels=labels)
+
+ best_key = max(all_executions, key=lambda k: all_executions[k]['best_obj_fn'])
+ results = dict(all_executions[best_key])
+ results['executions'] = all_executions
+ results['best_execution'] = best_key
+ results['selection_masks'] = masks
+ results['selection_frequency'] = frequency.sort_values('frequency', ascending=False)
+ results['cooccurrence'] = cooccurrence
+
+ if consensus_method == 'frequency':
+ mask, threshold = consensus_from_frequency(frequency['frequency'].values,
+ percentile=frequency_percentile)
+ results['consensus_ppi_data'] = pool.loc[mask].reset_index(drop=True)
+ results['consensus_threshold'] = threshold
+ if verbose:
+ print('Frequency consensus: {} pairs above {:.3f}'
+ .format(int(mask.sum()), threshold))
+ elif consensus_method == 'cooccurrence':
+ try:
+ selected, clusters, scores = consensus_from_cooccurrence(
+ cooccurrence, n_clusters=n_clusters, select=cluster_selection,
+ frequency=frequency['frequency'].values, min_frequency=min_frequency)
+ chosen = set(selected)
+ results['consensus_ppi_data'] = pool.loc[[l in chosen for l in labels]].reset_index(drop=True)
+ results['consensus_clusters'] = clusters
+ results['consensus_cluster_scores'] = scores
+ if verbose:
+ print('Co-occurrence consensus: {} pairs, cluster sizes {}, mean co-occurrence {}'
+ .format(len(chosen), {k: len(v) for k, v in clusters.items()},
+ {k: round(v, 3) for k, v in scores.items()}))
+ except ValueError as error:
+ if verbose:
+ print('No consensus could be built: {}'.format(error))
+ results['consensus_ppi_data'] = None
+ results['consensus_clusters'] = None
+ results['consensus_cluster_scores'] = None
+ else:
+ raise ValueError("`consensus_method` must be either 'cooccurrence' or 'frequency'")
+ return results
+
+
+def _as_symmetric(matrix):
+ '''
+ Validates a reference distance matrix and makes it exactly symmetric.
+
+ `check_symmetry` compares with exact equality, which a matrix produced by a
+ distance function often fails by a few ULP. Rather than rejecting those, they
+ are averaged with their transpose; genuinely asymmetric input still raises.
+ '''
+ values = np.asarray(matrix.values, dtype=float)
+ if values.shape[0] != values.shape[1]:
+ raise ValueError('`reference_distances` must be a square matrix')
+ if list(matrix.index) != list(matrix.columns):
+ raise ValueError('`reference_distances` must have the same cells as rows and columns')
+ if not np.allclose(values, values.T, rtol=1e-8, atol=1e-8, equal_nan=True):
+ raise ValueError('`reference_distances` must be a symmetric matrix')
+ return pd.DataFrame((values + values.T) / 2.0, index=matrix.index, columns=matrix.columns)
diff --git a/cell2cell/analysis/tensor_downstream.py b/cell2cell/analysis/tensor_downstream.py
index 34abd66..22dcac3 100644
--- a/cell2cell/analysis/tensor_downstream.py
+++ b/cell2cell/analysis/tensor_downstream.py
@@ -3,6 +3,8 @@
import numpy as np
import pandas as pd
+from natsort import natsorted
+
from cell2cell.stats import gini_coefficient
@@ -97,7 +99,7 @@ def get_factor_specific_ccc_networks(result, sender_label='Sender Cells', receiv
else:
raise ValueError('result is not of a valid type. It must be an InteractionTensor or a dict.')
- factors = sorted(list(set(result[sender_label].columns) & set(result[receiver_label].columns)))
+ factors = natsorted(set(result[sender_label].columns) & set(result[receiver_label].columns))
networks = dict()
for f in factors:
@@ -113,7 +115,7 @@ def flatten_factor_ccc_networks(networks, orderby='senders'):
'''
Flattens all adjacency matrices in the factor-specific
cell-cell communication networks. It generates a matrix
- where rows are factors and columns are cell-cell pairs.
+ where rows are cell-cell pairs and columns are factors.
Parameters
----------
@@ -131,11 +133,19 @@ def flatten_factor_ccc_networks(networks, orderby='senders'):
Returns
-------
flatten_networks : pandas.DataFrame
- A dataframe wherein rows contains a factor-specific network. Columns are
- the directed cell-cell pairs.
+ A dataframe wherein each column contains a factor-specific network. Rows are
+ the directed cell-cell pairs, named as ' --> '. Cells keep
+ the order they have in the tensor dimensions.
'''
- senders = sorted(set.intersection(*[set(v.index) for v in networks.values()]))
- receivers = sorted(set.intersection(*[set(v.columns) for v in networks.values()]))
+ net_list = list(networks.values())
+ common_senders = set.intersection(*[set(v.index) for v in net_list])
+ common_receivers = set.intersection(*[set(v.columns) for v in net_list])
+
+ # Keep the order of the elements in the tensor dimensions instead of sorting them.
+ # Sorting the names here without reordering the data assigns loadings to the wrong
+ # cell-cell pair whenever the tensor elements are not alphabetically sorted.
+ senders = [s for s in net_list[0].index if s in common_senders]
+ receivers = [r for r in net_list[0].columns if r in common_receivers]
if orderby == 'senders':
cell_pairs = [s + ' --> ' + r for s in senders for r in receivers]
@@ -146,7 +156,10 @@ def flatten_factor_ccc_networks(networks, orderby='senders'):
else:
raise ValueError("`orderby` must be either 'senders' or 'receivers'.")
- data = np.asarray([v.values.flatten(flatten_order) for v in networks.values()]).T
+ # Reindexing guarantees that every network is flattened in the same order the
+ # cell-cell pair names were built from.
+ data = np.asarray([v.reindex(index=senders, columns=receivers).values.flatten(flatten_order)
+ for v in net_list]).T
flatten_networks = pd.DataFrame(data=data,
index=cell_pairs,
columns=list(networks.keys())
@@ -191,7 +204,7 @@ def compute_gini_coefficients(result, sender_label='Sender Cells', receiver_labe
else:
raise ValueError('result is not of a valid type. It must be an InteractionTensor or a dict.')
- factors = sorted(list(set(result[sender_label].columns) & set(result[receiver_label].columns)))
+ factors = natsorted(set(result[sender_label].columns) & set(result[receiver_label].columns))
ginis = []
for f in factors:
@@ -276,7 +289,7 @@ def get_lr_by_cell_pairs(result, lr_label, sender_label, receiver_label, order_c
assert receiver_label in result.keys(), 'The specified dimension ' + receiver_label + ' is not present in the `result` input'
# Sort factors
- sorted_factors = sorted(result[lr_label].columns, key=lambda x: int(x.split(' ')[1]))
+ sorted_factors = natsorted(result[lr_label].columns)
# Get CCI network per factor
networks = get_factor_specific_ccc_networks(result=result,
diff --git a/cell2cell/clustering/cluster_interactions.py b/cell2cell/clustering/cluster_interactions.py
index 0c14a13..62a216a 100644
--- a/cell2cell/clustering/cluster_interactions.py
+++ b/cell2cell/clustering/cluster_interactions.py
@@ -86,10 +86,10 @@ def compute_linkage(distance_matrix, method='ward', optimal_ordering=True):
Z : numpy.ndarray
The hierarchical clustering encoded as a linkage matrix.
'''
- if (type(distance_matrix) is pd.core.frame.DataFrame):
- data = distance_matrix.values
- else:
- data = distance_matrix.copy()
+ # `np.array` accepts both dataframes and arrays, and always copies, so the diagonal
+ # can be zeroed below. The array behind `DataFrame.values` is read-only under the
+ # copy-on-write of pandas >= 3.0.
+ data = np.array(distance_matrix, dtype=float)
if ~(data.transpose() == data).all():
raise ValueError('The matrix is not symmetric')
diff --git a/cell2cell/core/interaction_space.py b/cell2cell/core/interaction_space.py
index d9c946b..ec377df 100644
--- a/cell2cell/core/interaction_space.py
+++ b/cell2cell/core/interaction_space.py
@@ -2,7 +2,8 @@
from __future__ import absolute_import
-from cell2cell.preprocessing import integrate_data, cutoffs, get_genes_from_complexes, add_complexes_to_expression
+from cell2cell.preprocessing import (integrate_data, cutoffs, get_genes_from_complexes, add_complexes_to_expression,
+ zero_diagonal)
from cell2cell.core import cell, cci_scores, communication_scores
import itertools
@@ -64,7 +65,9 @@ def generate_pairs(cells, cci_type, self_interaction=True, remove_duplicates=Tru
else:
raise NotImplementedError("CCI type has to be directed or undirected")
if remove_duplicates:
- pairs = list(set(pairs)) # Remove duplicates
+ # `dict.fromkeys` removes duplicates while keeping the order given by the
+ # list of cells, so the resulting pairs are reproducible across runs.
+ pairs = list(dict.fromkeys(pairs))
return pairs
@@ -510,14 +513,15 @@ def compute_pairwise_cci_scores(self, cci_score=None, use_ppi_score=False, verbo
# )
# Generate distance matrix
- if ~(cci_score in ['count', 'icellnet']):
+ if cci_score not in ['count', 'icellnet']:
self.distance_matrix = self.interaction_elements['cci_matrix'].apply(lambda x: 1 - x)
else:
#self.distance_matrix = self.interaction_elements['cci_matrix'].div(self.interaction_elements['cci_matrix'].max().max()).apply(lambda x: 1 - x)
# Regularized distance
mean = np.nanmean(self.interaction_elements['cci_matrix'])
self.distance_matrix = self.interaction_elements['cci_matrix'].div(self.interaction_elements['cci_matrix'] + mean).apply(lambda x: 1 - x)
- np.fill_diagonal(self.distance_matrix.values, 0.0) # Make diagonal zero (delete autocrine-interactions)
+ # Make diagonal zero (delete autocrine-interactions)
+ self.distance_matrix = zero_diagonal(self.distance_matrix)
def pair_communication_score(self, cell1, cell2, communication_score='expression_thresholding',
use_ppi_score=False, verbose=True):
diff --git a/cell2cell/datasets/__init__.py b/cell2cell/datasets/__init__.py
index bc7ec62..cb377c1 100644
--- a/cell2cell/datasets/__init__.py
+++ b/cell2cell/datasets/__init__.py
@@ -3,4 +3,6 @@
from cell2cell.datasets.heuristic_data import (HeuristicGOTerms)
from cell2cell.datasets.random_data import (generate_random_rnaseq, generate_random_ppi, generate_random_cci_scores,
generate_random_metadata)
-from cell2cell.datasets.toy_data import (generate_toy_distance, generate_toy_rnaseq, generate_toy_ppi, generate_toy_metadata)
\ No newline at end of file
+from cell2cell.datasets.toy_data import (generate_toy_distance, generate_toy_rnaseq, generate_toy_ppi, generate_toy_metadata,
+ generate_toy_contexts, generate_toy_single_cells, generate_toy_coordinates,
+ generate_toy_spatial_adata, generate_toy_liana_output)
\ No newline at end of file
diff --git a/cell2cell/datasets/toy_data.py b/cell2cell/datasets/toy_data.py
index ca2b803..dc916e1 100644
--- a/cell2cell/datasets/toy_data.py
+++ b/cell2cell/datasets/toy_data.py
@@ -111,4 +111,235 @@ def generate_toy_distance():
index=['C1', 'C2', 'C3', 'C4', 'C5'],
columns=['C1', 'C2', 'C3', 'C4', 'C5']
)
- return distance
\ No newline at end of file
+ return distance
+
+
+def generate_toy_contexts(n_contexts=4, context_names=None):
+ '''Generates a toy RNA-seq dataset for multiple contexts.
+
+ Each context contains the same genes and cells as the dataset generated by
+ `generate_toy_rnaseq()`, but with different expression values, so the contexts
+ can be used to build a 4D communication tensor.
+
+ Parameters
+ ----------
+ n_contexts : int, default=4
+ Number of contexts to generate. Contexts are named 'Context-1' to
+ 'Context-N'. Using a number equal to or greater than 10 results in
+ context names whose alphabetical order differs from their natural
+ order (e.g. 'Context-2' and 'Context-10').
+
+ context_names : list, default=None
+ Names to use for the contexts. If None, contexts are named 'Context-1'
+ to 'Context-N'. Its length must be equal to `n_contexts`.
+
+ Returns
+ -------
+ contexts : dict
+ Dictionary where keys are the context names and values are
+ pandas.DataFrame objects containing the gene expression of that context.
+ Columns are cells and rows are genes.
+ '''
+ if context_names is None:
+ context_names = ['Context-{}'.format(i) for i in range(1, n_contexts + 1)]
+ else:
+ assert len(context_names) == n_contexts, \
+ "The length of `context_names` must be equal to `n_contexts`"
+
+ base = generate_toy_rnaseq()
+
+ contexts = dict()
+ for i, name in enumerate(context_names):
+ # Shifting the genes and scaling the values generates a different
+ # communication pattern per context, while keeping it deterministic.
+ data = np.roll(base.values, i, axis=0) * (1.0 + 0.25 * i)
+ df = pd.DataFrame(data, index=base.index, columns=base.columns)
+ df.index.name = base.index.name
+ contexts[name] = df
+ return contexts
+
+
+def generate_toy_single_cells(n_cell_types=3, n_cells_per_type=4):
+ '''Generates a toy single-cell RNA-seq dataset with its metadata.
+
+ Parameters
+ ----------
+ n_cell_types : int, default=3
+ Number of cell types to generate. Cell types are named 'CT-1' to 'CT-N'.
+ Using a number equal to or greater than 10 results in cell-type names
+ whose alphabetical order differs from their natural order (e.g. 'CT-2'
+ and 'CT-10').
+
+ n_cells_per_type : int, default=4
+ Number of single cells to generate for each cell type.
+
+ Returns
+ -------
+ rnaseq : pandas.DataFrame
+ Gene expression of the single cells. Columns are single cells and rows
+ are genes, as in `generate_toy_rnaseq()`. To aggregate it with
+ `cell2cell.preprocessing.aggregate_single_cells`, either pass
+ `rnaseq.T` with `transposed=True` or `rnaseq` with `transposed=False`.
+
+ metadata : pandas.DataFrame
+ Metadata of the single cells. Contains the columns 'barcodes' and
+ 'cell_types', matching the default parameters of
+ `cell2cell.preprocessing.aggregate_single_cells`.
+ '''
+ base = generate_toy_rnaseq()
+
+ barcodes = []
+ cell_types = []
+ columns = []
+ for t in range(1, n_cell_types + 1):
+ cell_type = 'CT-{}'.format(t)
+ # Each cell type is based on one of the cells in the toy RNA-seq dataset
+ profile = base.iloc[:, (t - 1) % base.shape[1]].values.astype(float)
+ for c in range(1, n_cells_per_type + 1):
+ barcodes.append('{}-cell-{}'.format(cell_type, c))
+ cell_types.append(cell_type)
+ # Deterministic variation across the single cells of a cell type
+ columns.append(profile * (1.0 + 0.1 * (c - 1)) + (c - 1))
+
+ rnaseq = pd.DataFrame(np.asarray(columns).T, index=base.index, columns=barcodes)
+ rnaseq.index.name = base.index.name
+
+ metadata = pd.DataFrame({'barcodes': barcodes, 'cell_types': cell_types})
+ return rnaseq, metadata
+
+
+def generate_toy_coordinates(n_cell_types=3, n_cells_per_type=5):
+ '''Generates toy spatial coordinates for cells of different cell types.
+
+ Cells of a same cell type are placed close to each other, so the distances
+ between cell types are meaningful.
+
+ Parameters
+ ----------
+ n_cell_types : int, default=3
+ Number of cell types to generate. Cell types are named 'CT-1' to 'CT-N'.
+
+ n_cells_per_type : int, default=5
+ Number of cells to generate for each cell type.
+
+ Returns
+ -------
+ coordinates : pandas.DataFrame
+ DataFrame containing the columns 'X' and 'Y' with the coordinates of
+ each cell, and the column 'celltype' with its cell type. Rows are
+ indexed by the cell barcodes.
+ '''
+ records = []
+ for t in range(1, n_cell_types + 1):
+ cell_type = 'CT-{}'.format(t)
+ # Cell types are centered along a diagonal, separated by 50 units
+ center_x = 50.0 * ((t - 1) % 3)
+ center_y = 50.0 * ((t - 1) // 3)
+ for c in range(1, n_cells_per_type + 1):
+ records.append({'barcode': '{}-cell-{}'.format(cell_type, c),
+ 'X': center_x + (c - 1) * 2.0,
+ 'Y': center_y + ((c - 1) % 3) * 2.0,
+ 'celltype': cell_type})
+
+ coordinates = pd.DataFrame.from_records(records).set_index('barcode')
+ coordinates.index.name = 'barcode'
+ return coordinates
+
+
+def generate_toy_spatial_adata(num_cells=225, n_cell_types=3):
+ '''Generates a toy AnnData object containing spatial coordinates.
+
+ The cells are placed on a regular square lattice covering coordinates from
+ 0 to 100 in both dimensions, which makes the object usable with the
+ functions in `cell2cell.spatial`.
+
+ Parameters
+ ----------
+ num_cells : int, default=225
+ Number of cells/spots to generate.
+
+ n_cell_types : int, default=3
+ Number of cell types to assign to the cells. Cell types are named
+ 'CT-1' to 'CT-N'.
+
+ Returns
+ -------
+ adata : AnnData
+ Annotated data matrix with the toy genes as variables and the cells as
+ observations. Spatial coordinates are stored in `adata.obsm['spatial']`
+ and cell types in the 'celltype' column of `adata.obs`.
+ '''
+ import anndata
+
+ base = generate_toy_rnaseq()
+
+ # Regular lattice of coordinates, so the windows and grids are reproducible
+ side = int(np.ceil(np.sqrt(num_cells)))
+ x_coords, y_coords = np.meshgrid(np.linspace(0., 100., side),
+ np.linspace(0., 100., side))
+ coordinates = np.column_stack([x_coords.ravel(), y_coords.ravel()])[:num_cells]
+
+ n_genes = base.shape[0]
+ n_profiles = base.shape[1]
+ expression = np.zeros((num_cells, n_genes))
+ barcodes = []
+ cell_types = []
+ for i in range(num_cells):
+ profile = base.iloc[:, i % n_profiles].values.astype(float)
+ # Expression depends on the position, generating a spatial pattern
+ expression[i, :] = profile * (1.0 + coordinates[i, 0] / 100.)
+ barcodes.append('spot-{}'.format(i + 1))
+ cell_types.append('CT-{}'.format(i % n_cell_types + 1))
+
+ obs = pd.DataFrame({'celltype': cell_types}, index=barcodes)
+ obs.index.name = 'barcode'
+ var = pd.DataFrame(index=list(base.index))
+ var.index.name = 'gene_id'
+
+ adata = anndata.AnnData(X=expression, obs=obs, var=var)
+ adata.obsm['spatial'] = coordinates
+ return adata
+
+
+def generate_toy_liana_output(n_contexts=3, n_cell_types=3):
+ '''Generates a toy output resembling the one obtained from LIANA.
+
+ It contains the communication scores of each ligand-receptor pair for each
+ pair of sender-receiver cells, across multiple contexts, in a long format.
+
+ Parameters
+ ----------
+ n_contexts : int, default=3
+ Number of contexts to generate. Contexts are named 'Context-1' to
+ 'Context-N'.
+
+ n_cell_types : int, default=3
+ Number of cell types to generate. Cell types are named 'CT-1' to 'CT-N'.
+
+ Returns
+ -------
+ liana_outputs : pandas.DataFrame
+ Dataframe in a long format, containing the columns 'context', 'source',
+ 'target', 'ligand', 'receptor' and 'score'. Grouping it by the 'context'
+ column generates the dictionary that
+ `cell2cell.tensor.dataframes_to_tensor` expects.
+ '''
+ ppi = generate_toy_ppi(prot_complex=False)
+ lr_pairs = list(zip(ppi['A'], ppi['B']))
+
+ records = []
+ for i in range(1, n_contexts + 1):
+ context = 'Context-{}'.format(i)
+ for s in range(1, n_cell_types + 1):
+ for t in range(1, n_cell_types + 1):
+ for k, (ligand, receptor) in enumerate(lr_pairs):
+ # Deterministic score, different for each combination
+ score = ((i * 7 + s * 13 + t * 17 + k * 3) % 100) / 100.
+ records.append({'context': context,
+ 'source': 'CT-{}'.format(s),
+ 'target': 'CT-{}'.format(t),
+ 'ligand': ligand,
+ 'receptor': receptor,
+ 'score': score})
+ liana_outputs = pd.DataFrame.from_records(records)
+ return liana_outputs
\ No newline at end of file
diff --git a/cell2cell/external/pcoa.py b/cell2cell/external/pcoa.py
index fb0ac83..18a3c09 100644
--- a/cell2cell/external/pcoa.py
+++ b/cell2cell/external/pcoa.py
@@ -87,8 +87,10 @@ def pcoa(distance_matrix, method="eigh", number_of_dimensions=0,
"""
distance_matrix = convert_to_distance_matrix(distance_matrix)
- # Center distance matrix, a requirement for PCoA here
- matrix_data = center_distance_matrix(distance_matrix.values, inplace=inplace)
+ # Center distance matrix, a requirement for PCoA here. `np.array` always copies, so
+ # the `inplace` option gets a writable array, unlike `DataFrame.values` under the
+ # copy-on-write of pandas >= 3.0.
+ matrix_data = center_distance_matrix(np.array(distance_matrix, dtype=float), inplace=inplace)
# If no dimension specified, by default will compute all eigenvectors
# and eigenvalues
@@ -357,7 +359,9 @@ def pcoa_biplot(ordination, y):
raise ValueError('The eigenvectors and the descriptors must describe '
'the same samples.')
- eigvals = ordination['eigvals']
+ # Converted to a numpy array because using np.power() with the `where` argument
+ # on a pandas Series recurses through pandas' ufunc handling.
+ eigvals = np.asarray(ordination['eigvals'])
coordinates = ordination['samples']
N = coordinates.shape[0]
@@ -373,8 +377,11 @@ def pcoa_biplot(ordination, y):
#
# Only get the power of non-zero values, otherwise this will raise a
# divide by zero warning. There shouldn't be negative eigenvalues(?)
- Uproj = np.sqrt(N - 1) * spc.dot(np.diag(np.power(eigvals, -0.5,
- where=eigvals > 0)))
+ # `out` is needed so the entries excluded by `where` are zero instead of
+ # whatever was left in the uninitialized output array.
+ inverse_sqrt = np.power(eigvals, -0.5, where=eigvals > 0,
+ out=np.zeros_like(eigvals, dtype=float))
+ Uproj = np.sqrt(N - 1) * spc.dot(np.diag(inverse_sqrt))
ordination['features'] = pd.DataFrame(data=Uproj,
index=y.columns.copy(),
diff --git a/cell2cell/external/pcoa_utils.py b/cell2cell/external/pcoa_utils.py
index 39d5023..9703f3a 100644
--- a/cell2cell/external/pcoa_utils.py
+++ b/cell2cell/external/pcoa_utils.py
@@ -107,9 +107,12 @@ def scale(a, weights=None, with_mean=True, with_std=True, ddof=0, copy=True):
Wherever std equals 0, it is replaced by 1 in order to avoid
division by zero.
"""
- if copy:
- a = a.copy()
a = np.asarray(a, dtype=np.float64)
+ # `a` is standardized in place below. Copying it beforehand is not only what `copy`
+ # asks for, but also a requirement when the array is not writable, which is the case
+ # for the one behind `DataFrame.values` under the copy-on-write of pandas >= 3.0.
+ if copy or not a.flags.writeable:
+ a = a.copy()
avg, std = mean_and_std(a, axis=0, weights=weights, with_mean=with_mean,
with_std=with_std, ddof=ddof)
if with_mean:
@@ -218,7 +221,7 @@ def _e_matrix_inplace(distance_matrix):
distance_matrix : 2D array_like
Distance matrix.
"""
- distance_matrix = distance_matrix.astype(np.float)
+ distance_matrix = distance_matrix.astype(float)
for i in np.arange(len(distance_matrix)):
distance_matrix[i] = (distance_matrix[i] * distance_matrix[i]) / -2
@@ -238,7 +241,7 @@ def _f_matrix_inplace(e_matrix):
e_matrix : 2D array_like
A matrix representing the "E matrix" as described above.
"""
- e_matrix = e_matrix.astype(np.float)
+ e_matrix = e_matrix.astype(float)
row_means = np.zeros(len(e_matrix), dtype=float)
col_means = np.zeros(len(e_matrix), dtype=float)
diff --git a/cell2cell/io/directories.py b/cell2cell/io/directories.py
index 0231b87..ffbecb1 100644
--- a/cell2cell/io/directories.py
+++ b/cell2cell/io/directories.py
@@ -2,6 +2,8 @@
import os
+from natsort import natsorted
+
def create_directory(pathname):
'''Creates a directory.
@@ -37,8 +39,10 @@ def get_files_from_directory(pathname, dir_in_filepath=False):
-------
filenames : list
A list containing the names (strings) of the files
- in the folder.
+ in the folder, naturally sorted by filename.
'''
directory = os.fsencode(pathname)
- filenames = [pathname + '/' + os.fsdecode(file) if dir_in_filepath else os.fsdecode(file) for file in os.listdir(directory)]
+ # Naturally sorted to avoid a filesystem-dependent order of the files
+ files = natsorted([os.fsdecode(file) for file in os.listdir(directory)])
+ filenames = [pathname + '/' + file if dir_in_filepath else file for file in files]
return filenames
diff --git a/cell2cell/plotting/aesthetics.py b/cell2cell/plotting/aesthetics.py
index 5c352f4..fe1bd1b 100644
--- a/cell2cell/plotting/aesthetics.py
+++ b/cell2cell/plotting/aesthetics.py
@@ -6,6 +6,8 @@
import matplotlib.patches as patches
import numpy as np
+from natsort import natsorted
+
def get_colors_from_labels(labels, cmap='gist_rainbow', factor=1):
'''Generates colors for each label in a list given a colormap
@@ -158,7 +160,7 @@ def generate_legend(color_dict, loc='center left', bbox_to_anchor=(1.01, 0.5), n
'''
color_patches = []
if sorted_labels:
- iteritems = sorted(color_dict.items())
+ iteritems = natsorted(color_dict.items(), key=lambda x: x[0])
else:
iteritems = color_dict.items()
for k, v in iteritems:
diff --git a/cell2cell/plotting/circular_plot.py b/cell2cell/plotting/circular_plot.py
index 5cdbbab..978b270 100644
--- a/cell2cell/plotting/circular_plot.py
+++ b/cell2cell/plotting/circular_plot.py
@@ -9,6 +9,8 @@
import numpy as np
import pandas as pd
+from natsort import natsorted
+
from cell2cell.plotting.aesthetics import get_colors_from_labels, generate_legend
@@ -165,7 +167,7 @@ def circos_plot(interaction_space, sender_cells, receiver_cells, ligands, recept
small_R = determine_small_radius(edges_dict)
# Colors
- cells = list(set(sender_cells+receiver_cells))
+ cells = natsorted(set(sender_cells + receiver_cells))
if metadata is not None:
meta = metadata.set_index(sample_col).reindex(cells)
meta = meta[[group_col]].fillna('NA')
@@ -303,7 +305,7 @@ def get_arc_angles(G, sorting_feature=None):
values are tuples with angles for the start and end of
the arc that represents a node.
'''
- elements = list(set(G.nodes()))
+ elements = natsorted(G.nodes())
n_elements = len(elements)
if sorting_feature is not None:
diff --git a/cell2cell/plotting/factor_plot.py b/cell2cell/plotting/factor_plot.py
index e0742b2..c3c1829 100644
--- a/cell2cell/plotting/factor_plot.py
+++ b/cell2cell/plotting/factor_plot.py
@@ -6,6 +6,7 @@
from matplotlib import pyplot as plt
from statannotations.Annotator import Annotator
from scipy.stats import zscore
+from natsort import natsorted
from cell2cell.clustering.cluster_interactions import compute_distance, compute_linkage
from cell2cell.analysis.tensor_downstream import get_factor_specific_ccc_networks
@@ -116,7 +117,7 @@ def context_boxplot(context_loadings, metadict, included_factors=None, group_ord
if group_order is not None:
assert len(set(group_order) & set(metadict.values())) == len(set(metadict.values())), "All groups in `metadict` must be contained in `group_order`"
else:
- group_order = list(set(metadict.values()))
+ group_order = natsorted(set(metadict.values()))
df = context_loadings.copy()
if included_factors is None:
@@ -143,16 +144,25 @@ def context_boxplot(context_loadings, metadict, included_factors=None, group_ord
order = group_order
# Plot the boxes
+ # `hue` repeats `x` because seaborn 0.14 removes the option of passing a palette
+ # without one. `dodge=False` keeps one box per group instead of splitting them,
+ # which is what makes this identical to passing the palette on its own. It is
+ # preferred over the `legend` argument, which only exists from seaborn 0.13.
ax = sns.boxplot(x=x,
y=y,
data=df,
order=order,
whis=[0, 100],
width=.6,
+ hue=x,
palette=cmap,
+ dodge=False,
boxprops=dict(alpha=.5),
ax=ax
)
+ # The x axis already labels the groups, so the legend the hue brings is redundant
+ if ax.get_legend() is not None:
+ ax.get_legend().remove()
# Plot the dots
sns.stripplot(x=x,
diff --git a/cell2cell/plotting/pcoa_plot.py b/cell2cell/plotting/pcoa_plot.py
index fb14ca4..c31f23a 100644
--- a/cell2cell/plotting/pcoa_plot.py
+++ b/cell2cell/plotting/pcoa_plot.py
@@ -5,6 +5,8 @@
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
+from natsort import natsorted
+
from cell2cell.external import pcoa, _check_ordination
from cell2cell.plotting.aesthetics import get_colors_from_labels
@@ -138,7 +140,7 @@ def pcoa_3dplot(interaction_space, metadata=None, sample_col='#SampleID', group_
assert all(elem in colors.keys() for elem in set(labels))
# Plot each data point with respective color
- for i, cell_type in enumerate(sorted(meta_[group_col].unique())):
+ for i, cell_type in enumerate(natsorted(meta_[group_col].unique())):
cells = list(meta_.loc[meta_[group_col] == cell_type].index)
if colors is not None:
ax.scatter(ordination['samples'].loc[cells, 'PC1'],
diff --git a/cell2cell/plotting/pval_plot.py b/cell2cell/plotting/pval_plot.py
index c204930..cf92869 100644
--- a/cell2cell/plotting/pval_plot.py
+++ b/cell2cell/plotting/pval_plot.py
@@ -189,7 +189,8 @@ def generate_dot_plot(pval_df, score_df, significance=0.05, xlabel='', ylabel=''
# Drop all zeros
df = df.loc[(df != 0).any(axis=1)]
df = df.T.loc[(df != 0).any(axis=0)].T
- pval_df = pval_df[df.columns].loc[df.index].applymap(lambda x: -1. * np.log10(x + 1e-9))
+ # Vectorized instead of elementwise, since `DataFrame.applymap` was removed in pandas 3.0
+ pval_df = -1. * np.log10(pval_df[df.columns].loc[df.index] + 1e-9)
n_rows = len(pval_df.index)
n_cols = len(pval_df.columns)
@@ -211,8 +212,9 @@ def generate_dot_plot(pval_df, score_df, significance=0.05, xlabel='', ylabel=''
norm = mpl.colors.Normalize(vmin=-1. * max_abs, vmax=max_abs)
max_size = mpl.colors.Normalize(vmin=0., vmax=3)
- # Colormap
- cmap = mpl.cm.get_cmap(cmap)
+ # Colormap. `matplotlib.cm.get_cmap` was removed in matplotlib 3.11, and the pyplot
+ # function is the one used elsewhere in this subpackage. It also takes a Colormap.
+ cmap = plt.get_cmap(cmap)
# Create figure with proper height ratios
# Use height_ratios based on actual inches rather than arbitrary numbers
diff --git a/cell2cell/plotting/tensor_plot.py b/cell2cell/plotting/tensor_plot.py
index fe6cf96..66f1be1 100644
--- a/cell2cell/plotting/tensor_plot.py
+++ b/cell2cell/plotting/tensor_plot.py
@@ -469,7 +469,8 @@ def reorder_dimension_elements(factors, reorder_elements, metadata=None):
assert all((len(set(factors[key].index).difference(set(reorder_elements[key]))) == 0) for key in reorder_elements.keys()), "All elements of each dimension included should be present"
reordered_factors = factors.copy()
- new_metadata = metadata.copy()
+ # `metadata` is optional, so it is only copied when it was actually provided
+ new_metadata = metadata.copy() if metadata is not None else None
i = 0
for k, df in reordered_factors.items():
@@ -596,7 +597,7 @@ def plot_multiple_run_elbow(all_loss, elbow=None, ci='95%', figsize=(4, 2.25), y
raise ValueError("Specify a correct ci. Either '95%' or 'std'")
plt.fill_between(x, mean - coeff * std, mean + coeff * std, color='steelblue', alpha=.2,
- label='$\pm$ 1 std')
+ label=r'$\pm$ 1 std')
plt.tick_params(axis='both', labelsize=fontsize)
plt.xlabel('Rank', fontsize=int(1.2 * fontsize))
diff --git a/cell2cell/plotting/umap_plot.py b/cell2cell/plotting/umap_plot.py
index adf0bfb..45c1365 100644
--- a/cell2cell/plotting/umap_plot.py
+++ b/cell2cell/plotting/umap_plot.py
@@ -54,11 +54,15 @@ def umap_biplot(umap_df, figsize=(8 ,8), ax=None, show_axes=True, show_legend=Tr
if ax is None:
fig = plt.figure(figsize=figsize)
+ # A palette is only meaningful together with a hue. Passing it without one was
+ # ignored with a warning by seaborn, and is removed in seaborn 0.14.
+ palette = cmap if hue is not None else None
+
ax = sns.scatterplot(x='umap1',
y='umap2',
data=umap_df,
hue=hue,
- palette=cmap,
+ palette=palette,
ax=ax
)
diff --git a/cell2cell/preprocessing/__init__.py b/cell2cell/preprocessing/__init__.py
index 0e1d7e6..3067e65 100644
--- a/cell2cell/preprocessing/__init__.py
+++ b/cell2cell/preprocessing/__init__.py
@@ -6,7 +6,8 @@
from cell2cell.preprocessing.integrate_data import (get_thresholded_rnaseq, get_modified_rnaseq, get_ppi_dict_from_go_terms,
get_ppi_dict_from_proteins, get_weighted_ppi)
from cell2cell.preprocessing.manipulate_dataframes import (check_presence_in_dataframe, shuffle_cols_in_df, shuffle_rows_in_df,
- shuffle_dataframe, subsample_dataframe)
+ shuffle_dataframe, subsample_dataframe, check_symmetry,
+ convert_to_distance_matrix, zero_diagonal)
from cell2cell.preprocessing.ppi import (bidirectional_ppi_for_cci, filter_ppi_by_proteins, filter_ppi_network,
get_all_to_all_ppi, get_filtered_ppi_network, get_one_group_to_other_ppi,
remove_ppi_bidirectionality, simplify_ppi, filter_complex_ppi_by_proteins,
diff --git a/cell2cell/preprocessing/find_elements.py b/cell2cell/preprocessing/find_elements.py
index d28d246..2e6f7a2 100644
--- a/cell2cell/preprocessing/find_elements.py
+++ b/cell2cell/preprocessing/find_elements.py
@@ -44,9 +44,12 @@ def get_element_abundances(element_lists):
abundance_dict : dict
Dictionary containing the number of times that an
element was present, divided by the total number of
- lists in `element_lists`.
+ lists in `element_lists`. Keys keep the order in which
+ the elements were first found across `element_lists`.
'''
- abundance_dict = Counter(itertools.chain(*map(set, element_lists)))
+ # `dict.fromkeys` removes duplicates within each list while keeping their order,
+ # so that the resulting keys are reproducible across runs (unlike using sets).
+ abundance_dict = Counter(itertools.chain(*[dict.fromkeys(l) for l in element_lists]))
total = len(element_lists)
abundance_dict = {k : v/total for k, v in abundance_dict.items()}
return abundance_dict
diff --git a/cell2cell/preprocessing/gene_ontology.py b/cell2cell/preprocessing/gene_ontology.py
index f0c693f..a2eff5e 100644
--- a/cell2cell/preprocessing/gene_ontology.py
+++ b/cell2cell/preprocessing/gene_ontology.py
@@ -5,6 +5,8 @@
import numpy as np
import networkx
+from natsort import natsorted
+
def get_genes_from_go_terms(go_annotations, go_filter, go_header='GO', gene_header='Gene', verbose=True):
'''
@@ -80,7 +82,7 @@ def get_genes_from_go_hierarchy(go_annotations, go_terms, go_filter, go_header='
iter = len(go_hierarchy)
for i in range(iter):
find_all_children_of_go_term(go_terms, go_hierarchy[i], go_hierarchy, verbose=verbose)
- go_hierarchy = list(set(go_hierarchy))
+ go_hierarchy = natsorted(set(go_hierarchy))
genes = get_genes_from_go_terms(go_annotations=go_annotations,
go_filter=go_hierarchy,
go_header=go_header,
diff --git a/cell2cell/preprocessing/manipulate_dataframes.py b/cell2cell/preprocessing/manipulate_dataframes.py
index d7d5f81..b6d16a3 100644
--- a/cell2cell/preprocessing/manipulate_dataframes.py
+++ b/cell2cell/preprocessing/manipulate_dataframes.py
@@ -3,6 +3,8 @@
from __future__ import absolute_import
import random
+import warnings
+
import numpy as np
import pandas as pd
@@ -33,9 +35,14 @@ def check_presence_in_dataframe(df, elements, columns=None):
'''
if columns is None:
columns = list(df.columns)
- df_elements = pd.Series(np.unique(df[columns].values.flatten()))
- df_elements = df_elements.loc[df_elements.isin(elements)].values
- found_elements = list(df_elements)
+ elif isinstance(columns, str):
+ # A string would be read by pandas as a single column name, returning a Series
+ # instead of the dataframe the code below expects, as in `shuffle_cols_in_df`.
+ columns = [columns]
+ # `pd.unique` does not sort the values, so it also works when the considered
+ # columns contain a mix of data types (e.g. gene names and scores).
+ df_elements = pd.Series(pd.unique(df[columns].to_numpy(dtype=object).ravel()))
+ found_elements = df_elements.loc[df_elements.isin(elements)].tolist()
return found_elements
@@ -136,7 +143,9 @@ def shuffle_dataframe(df, shuffling_number=1, axis=0, random_state=None):
'''
df_ = df.copy()
axis = int(not axis) # pandas.DataFrame is always 2D
- to_shuffle = np.rollaxis(df_.values, axis)
+ # `to_numpy(copy=True)` is shuffled in place below, so it must be writable. The array
+ # behind `DataFrame.values` is read-only under the copy-on-write of pandas >= 3.0.
+ to_shuffle = np.rollaxis(df_.to_numpy(copy=True), axis)
for _ in range(shuffling_number):
for i, view in enumerate(to_shuffle):
if random_state is not None:
@@ -202,6 +211,29 @@ def check_symmetry(df):
return symmetric
+def zero_diagonal(df):
+ '''
+ Sets all diagonal elements of a square dataframe to zero.
+
+ Parameters
+ ----------
+ df : pandas.DataFrame
+ A square dataframe.
+
+ Returns
+ -------
+ df_ : pandas.DataFrame
+ A copy of df, but with all diagonal elements with a
+ value of zero.
+ '''
+ # `np.array` always copies, so the result is writable. The array behind
+ # `DataFrame.values` is read-only under the copy-on-write of pandas >= 3.0,
+ # which makes it unusable with in-place functions such as `np.fill_diagonal`.
+ values = np.array(df, dtype=float)
+ np.fill_diagonal(values, 0.0)
+ return pd.DataFrame(values, index=df.index, columns=df.columns)
+
+
def convert_to_distance_matrix(df):
'''
Converts a symmetric dataframe into a distance dataframe.
@@ -219,10 +251,10 @@ def convert_to_distance_matrix(df):
value of zero.
'''
if check_symmetry(df):
- df_ = df.copy()
- if np.trace(df_.values,) != 0.0:
- raise Warning("Diagonal elements are not zero. Automatically replaced by zeros")
- np.fill_diagonal(df_.values, 0.0)
+ if np.trace(df.values,) != 0.0:
+ # Warned instead of raised, so the diagonal is actually replaced below
+ warnings.warn("Diagonal elements are not zero. Automatically replaced by zeros")
+ df_ = zero_diagonal(df)
else:
raise ValueError('The DataFrame is not symmetric')
return df_
diff --git a/cell2cell/preprocessing/ppi.py b/cell2cell/preprocessing/ppi.py
index beae77b..3158b36 100644
--- a/cell2cell/preprocessing/ppi.py
+++ b/cell2cell/preprocessing/ppi.py
@@ -7,6 +7,8 @@
from itertools import combinations
+from natsort import natsorted
+
### Preprocess a PPI table from a known list
def preprocess_ppi_data(ppi_data, interaction_columns, sort_values=None, score=None, rnaseq_genes=None, complex_sep=None,
@@ -596,7 +598,7 @@ def get_filtered_ppi_network(ppi_data, contact_proteins, mediator_proteins=None,
interaction_columns=interaction_columns)
elif interaction_type == 'complete':
- total_proteins = list(set(contact_proteins + mediator_proteins))
+ total_proteins = natsorted(set(contact_proteins + mediator_proteins))
new_ppi_data = get_all_to_all_ppi(ppi_data=ppi_data,
proteins=total_proteins,
diff --git a/cell2cell/preprocessing/rnaseq.py b/cell2cell/preprocessing/rnaseq.py
index 4daef0d..763e47d 100644
--- a/cell2cell/preprocessing/rnaseq.py
+++ b/cell2cell/preprocessing/rnaseq.py
@@ -5,6 +5,8 @@
import numpy as np
import pandas as pd
+from natsort import natsorted
+
### Pre-process RNAseq datasets
def drop_empty_genes(rnaseq_data):
@@ -82,7 +84,9 @@ def scale_expression_by_sum(rnaseq_data, axis=0, sum_value=1e6):
cell-types/tissues/samples and rows are genes.
'''
data = rnaseq_data.values
- data = sum_value * np.divide(data, np.nansum(data, axis=axis))
+ # `keepdims` is needed so the sums broadcast back along the specified axis.
+ # Without it, normalizing across columns (axis=1) raises a broadcasting error.
+ data = sum_value * np.divide(data, np.nansum(data, axis=axis, keepdims=True))
scaled_data = pd.DataFrame(data, index=rnaseq_data.index, columns=rnaseq_data.columns)
return scaled_data
@@ -196,9 +200,33 @@ def add_complexes_to_expression(rnaseq_data, complexes, agg_method='min'):
return tmp_rna
+def _trimean(x, axis):
+ '''
+ Computes the trimean of the data along the specified axis.
+
+ Parameters
+ ----------
+ x : numpy.ndarray
+ The input data for which the trimean is to be computed.
+
+ axis : int
+ The axis along which to compute the trimean. Use 0 for columns, 1 for rows.
+
+ Returns
+ -------
+ trimean : numpy.ndarray
+ An array containing the trimean values for each row or column, depending on
+ the specified axis.
+ '''
+ q1, q2, q3 = np.nanpercentile(x, [25, 50, 75], axis=axis)
+ trimean = 0.5 * q2 + 0.25 * (q1 + q3)
+ return trimean
+
+
def aggregate_single_cells(rnaseq_data, metadata, barcode_col='barcodes', celltype_col='cell_types', method='average',
transposed=True):
- '''Aggregates gene expression of single cells into cell types for each gene.
+ '''
+ Aggregates gene expression of single cells into cell types for each gene.
Parameters
----------
@@ -229,6 +257,11 @@ def aggregate_single_cells(rnaseq_data, metadata, barcode_col='barcodes', cellty
of a given gene.
- 'average' : Computes the average gene expression among the single cells
composing a cell type for a given gene.
+ - 'trimean' : Computes the Tukey's trimean of the gene expression among the
+ single cells composing a cell type for a given gene. It is a weighted
+ average of the median and the first and third quartiles
+ (0.5 * Q2 + 0.25 * (Q1 + Q3)), so it is more robust to outliers than
+ the average while still accounting for the spread of the distribution.
transposed : boolean, default=True
Whether the rnaseq_data is organized with columns as
@@ -241,7 +274,7 @@ def aggregate_single_cells(rnaseq_data, metadata, barcode_col='barcodes', cellty
by cell types. Columns are cell types and rows are genes.
'''
assert metadata is not None, "Please provide metadata containing the barcodes and cell-type annotation."
- assert method in ['average', 'nn_cell_fraction'], "{} is not a valid option for method".format(method)
+ assert method in ['average', 'nn_cell_fraction', 'trimean'], "{} is not a valid option for method".format(method)
meta = metadata.reset_index()
meta = meta[[barcode_col, celltype_col]].set_index(barcode_col)
@@ -251,18 +284,24 @@ def aggregate_single_cells(rnaseq_data, metadata, barcode_col='barcodes', cellty
df = rnaseq_data
else:
df = rnaseq_data.T
- df.index = [mapper[c] for c in df.index]
- df.index.name = 'celltype'
- df.reset_index(inplace=True)
- agg_df = pd.DataFrame(index=df.columns).drop('celltype')
+ # Grouping by an external list of cell types, instead of replacing the index of
+ # `df` and adding a column to it, avoids modifying the dataframe passed by the user.
+ celltypes = [mapper[c] for c in df.index]
+
+ agg_df = pd.DataFrame(index=df.columns)
- for celltype, ct_df in df.groupby('celltype'):
- ct_df = ct_df.drop('celltype', axis=1)
+ for celltype, ct_df in df.groupby(celltypes):
if method == 'average':
agg = ct_df.mean()
elif method == 'nn_cell_fraction':
agg = ((ct_df > 0).sum() / ct_df.shape[0])
+ elif method == 'trimean':
+ agg = pd.Series(_trimean(ct_df.values, axis=0), index=ct_df.columns)
agg_df[celltype] = agg
+
+ # Naturally sorted to avoid a lexicographic order of the cell types (e.g. to obtain
+ # 'CT-1', 'CT-2', 'CT-10' instead of 'CT-1', 'CT-10', 'CT-2')
+ agg_df = agg_df[natsorted(agg_df.columns)]
return agg_df
diff --git a/cell2cell/spatial/__init__.py b/cell2cell/spatial/__init__.py
index 45f6da4..74b9088 100644
--- a/cell2cell/spatial/__init__.py
+++ b/cell2cell/spatial/__init__.py
@@ -1,3 +1,5 @@
-from cell2cell.spatial.distances import (celltype_pair_distance, pairwise_celltype_distances)
+from cell2cell.spatial.distances import (celltype_pair_distance, pairwise_celltype_distances,
+ get_spatial_coordinates, celltype_centroids,
+ celltype_centroid_distances, celltype_distances)
from cell2cell.spatial.filtering import (dist_filter_liana, dist_filter_tensor)
-from cell2cell.spatial.neighborhoods import (create_spatial_grid, create_sliding_windows, calculate_window_size, add_sliding_window_info_to_adata)
\ No newline at end of file
+from cell2cell.spatial.neighborhoods import (create_spatial_grid, create_sliding_windows, calculate_window_size, add_sliding_window_info_to_adata)
diff --git a/cell2cell/spatial/distances.py b/cell2cell/spatial/distances.py
index 5e20043..0681cb9 100644
--- a/cell2cell/spatial/distances.py
+++ b/cell2cell/spatial/distances.py
@@ -2,6 +2,7 @@
import itertools
import numpy as np
import pandas as pd
+from natsort import natsorted
from sklearn.metrics.pairwise import euclidean_distances, manhattan_distances
@@ -22,7 +23,7 @@ def celltype_pair_distance(df1, df2, method='min', distance='euclidean'):
method : str, default='min'
The aggregation method for the calculated distances. It can be one of 'min',
- 'max', or 'mean'.
+ 'max', 'mean', or 'median'.
distance : str, default='euclidean'
The distance metric to use. It can be 'euclidean' or 'manhattan'.
@@ -46,6 +47,8 @@ def celltype_pair_distance(df1, df2, method='min', distance='euclidean'):
agg_dist = np.nanmax(distances)
elif method == 'mean':
agg_dist = np.nanmean(distances)
+ elif method == 'median':
+ agg_dist = np.nanmedian(distances)
else:
raise NotImplementedError('Method {} is not implemented.'.format(method))
return agg_dist
@@ -98,4 +101,283 @@ def pairwise_celltype_distances(df, group_col, coord_cols=['X', 'Y'],
)
distances.loc[pair[0], pair[1]] = dist
distances.loc[pair[1], pair[0]] = dist
- return distances
\ No newline at end of file
+ return distances
+
+def get_spatial_coordinates(adata, spatial_key='spatial', coord_names=None):
+ '''
+ Extracts the spatial coordinates of an AnnData object as a dataframe.
+
+ Parameters
+ ----------
+ adata : anndata.AnnData
+ Object containing the spatial coordinates of each single cell.
+
+ spatial_key : str, default='spatial'
+ Key in `adata.obsm` where the coordinates are stored. Objects written by
+ different tools use different keys (e.g. 'spatial', 'X_spatial',
+ 'X_umap'), so it can be changed here.
+
+ coord_names : list, default=None
+ Names to give to the coordinate columns. If None, they are named 'X', 'Y'
+ and 'Z' for the first three dimensions, and 'Dim4', 'Dim5', ... beyond
+ that, so the result works with the `coord_cols` parameter of the other
+ functions in this module.
+
+ Returns
+ -------
+ coordinates : pandas.DataFrame
+ Coordinates of each single cell. Rows are the observation names of
+ `adata`, in the same order, and columns are the dimensions.
+ '''
+ if spatial_key not in adata.obsm.keys():
+ raise KeyError("'{}' is not in adata.obsm. Available keys are: {}"
+ .format(spatial_key, list(adata.obsm.keys())))
+
+ coords = np.asarray(adata.obsm[spatial_key])
+ if coords.ndim != 2:
+ raise ValueError('The coordinates in adata.obsm[\'{}\'] must be two-dimensional'
+ .format(spatial_key))
+
+ if coord_names is None:
+ default = ['X', 'Y', 'Z']
+ coord_names = [default[i] if i < len(default) else 'Dim{}'.format(i + 1)
+ for i in range(coords.shape[1])]
+ elif len(coord_names) != coords.shape[1]:
+ raise ValueError('`coord_names` must have one name per dimension ({})'
+ .format(coords.shape[1]))
+
+ return pd.DataFrame(coords, index=adata.obs_names, columns=coord_names)
+
+
+def celltype_centroids(adata, group_col, spatial_key='spatial', coord_names=None,
+ method='mean'):
+ '''
+ Computes the centroid of each cell type from the coordinates of its single cells.
+
+ Parameters
+ ----------
+ adata : anndata.AnnData or pandas.DataFrame
+ Either an AnnData object with coordinates in `adata.obsm[spatial_key]` and
+ the cell-type annotation in `adata.obs[group_col]`, or a dataframe with one
+ row per single cell containing both the coordinates and the annotation.
+
+ group_col : str
+ Column with the cell-type annotation. Taken from `adata.obs` for an AnnData
+ object, and from the dataframe itself otherwise.
+
+ spatial_key : str, default='spatial'
+ Key in `adata.obsm` where the coordinates are stored. Ignored when a
+ dataframe is passed.
+
+ coord_names : list, default=None
+ Names of the coordinate columns. For an AnnData object they name the
+ extracted dimensions; for a dataframe they select which columns to use. If
+ None, an AnnData is named 'X', 'Y', 'Z', ... and for a dataframe every
+ column other than `group_col` is used.
+
+ method : str, default='mean'
+ How to summarize the coordinates of the single cells of a cell type. It can
+ be 'mean' (the centroid proper) or 'median' (the component-wise median,
+ which is robust to cells scattered far from the rest of their type).
+
+ Returns
+ -------
+ centroids : pandas.DataFrame
+ One row per cell type and one column per dimension. Cell types are
+ naturally sorted, so 'CT-2' comes before 'CT-10'.
+
+ Examples
+ --------
+ >>> import cell2cell as c2c
+ >>> adata = c2c.datasets.generate_toy_spatial_adata()
+ >>> centroids = c2c.spatial.celltype_centroids(adata, group_col='cell_type')
+ '''
+ coords, groups = _coordinates_and_groups(adata, group_col, spatial_key, coord_names)
+
+ if method == 'mean':
+ centroids = coords.groupby(groups, observed=True).mean()
+ elif method == 'median':
+ centroids = coords.groupby(groups, observed=True).median()
+ else:
+ raise NotImplementedError("Method {} is not implemented. Use 'mean' or 'median'."
+ .format(method))
+ centroids.index.name = group_col
+ return centroids.loc[natsorted(centroids.index)]
+
+
+def celltype_centroid_distances(adata, group_col, spatial_key='spatial', coord_names=None,
+ centroid_method='mean', distance='euclidean'):
+ '''
+ Computes the distances between the centroids of every pair of cell types.
+
+ This summarizes each cell type by one point before measuring distances, so its
+ cost does not depend on how many single cells each type contains. That makes it
+ the option to use on large datasets, where the all-versus-all single-cell
+ distances of `pairwise_celltype_distances` become prohibitive.
+
+ Parameters
+ ----------
+ adata : anndata.AnnData or pandas.DataFrame
+ Object or dataframe containing the coordinates and the cell-type annotation.
+
+ group_col : str
+ Column with the cell-type annotation.
+
+ spatial_key : str, default='spatial'
+ Key in `adata.obsm` where the coordinates are stored.
+
+ coord_names : list, default=None
+ Names of the coordinate columns.
+
+ centroid_method : str, default='mean'
+ How to summarize the coordinates of each cell type, 'mean' or 'median'.
+
+ distance : str, default='euclidean'
+ The distance metric to use. It can be 'euclidean' or 'manhattan'.
+
+ Returns
+ -------
+ distances : pandas.DataFrame
+ Symmetric matrix with a zero diagonal, where rows and columns are the cell
+ types, naturally sorted.
+
+ Examples
+ --------
+ >>> import cell2cell as c2c
+ >>> adata = c2c.datasets.generate_toy_spatial_adata()
+ >>> distances = c2c.spatial.celltype_centroid_distances(adata, group_col='cell_type')
+ '''
+ centroids = celltype_centroids(adata, group_col, spatial_key=spatial_key,
+ coord_names=coord_names, method=centroid_method)
+
+ if distance == 'euclidean':
+ matrix = euclidean_distances(centroids.values, centroids.values)
+ elif distance == 'manhattan':
+ matrix = manhattan_distances(centroids.values, centroids.values)
+ else:
+ raise NotImplementedError("{} distance is not implemented.".format(distance.capitalize()))
+
+ # Forced rather than assumed, so the result always satisfies `check_symmetry`
+ # and `squareform`, whatever rounding the metric introduced
+ matrix = (matrix + matrix.T) / 2.0
+ np.fill_diagonal(matrix, 0.0)
+ return pd.DataFrame(matrix, index=centroids.index, columns=centroids.index)
+
+
+def celltype_distances(adata, group_col, spatial_key='spatial', coord_names=None,
+ method='centroid', distance='euclidean', centroid_method='mean',
+ pairs=None, verbose=False):
+ '''
+ Computes a distance between every pair of cell types from the coordinates of
+ their single cells.
+
+ Single entry point for the two ways of summarizing the distance between two
+ cell types: aggregating the distances between all of their single cells
+ ('min', 'max', 'mean', 'median'), or measuring between their centroids
+ ('centroid').
+
+ Parameters
+ ----------
+ adata : anndata.AnnData or pandas.DataFrame
+ Object or dataframe containing the coordinates and the cell-type annotation.
+
+ group_col : str
+ Column with the cell-type annotation.
+
+ spatial_key : str, default='spatial'
+ Key in `adata.obsm` where the coordinates are stored. Ignored when a
+ dataframe is passed.
+
+ coord_names : list, default=None
+ Names of the coordinate columns.
+
+ method : str, default='centroid'
+ How to summarize the distance between two cell types:
+
+ - 'centroid' : distance between the centroids of the two cell types. Cost
+ is independent of the number of single cells, so this is the one to use
+ on large datasets.
+ - 'min' : smallest distance between any two of their single cells, i.e.
+ how close the two types get to each other.
+ - 'max' : largest distance between any two of their single cells.
+ - 'mean' : average over all pairs of their single cells.
+ - 'median' : median over all pairs of their single cells, less sensitive
+ to a few distant cells than 'mean'.
+
+ Every option other than 'centroid' evaluates all pairs of single cells of
+ the two types, so its cost grows with the product of their sizes.
+
+ distance : str, default='euclidean'
+ The distance metric to use. It can be 'euclidean' or 'manhattan'.
+
+ centroid_method : str, default='mean'
+ How to summarize the coordinates of each cell type when
+ `method='centroid'`, either 'mean' or 'median'.
+
+ pairs : list, default=None
+ Specific pairs of cell types to compute. If None, all combinations are
+ used. Ignored when `method='centroid'`, which computes all of them at once.
+
+ verbose : boolean, default=False
+ Whether to warn when the all-versus-all computation is going to be large.
+
+ Returns
+ -------
+ distances : pandas.DataFrame
+ Symmetric matrix with a zero diagonal, where rows and columns are the cell
+ types, naturally sorted.
+
+ Examples
+ --------
+ >>> import cell2cell as c2c
+ >>> adata = c2c.datasets.generate_toy_spatial_adata()
+ >>> # Fast, and the sensible default on large data
+ >>> distances = c2c.spatial.celltype_distances(adata, group_col='cell_type')
+ >>> # How close the two cell types get to each other
+ >>> distances = c2c.spatial.celltype_distances(adata, group_col='cell_type',
+ ... method='min')
+ '''
+ if method == 'centroid':
+ return celltype_centroid_distances(adata, group_col, spatial_key=spatial_key,
+ coord_names=coord_names,
+ centroid_method=centroid_method,
+ distance=distance)
+
+ coords, groups = _coordinates_and_groups(adata, group_col, spatial_key, coord_names)
+
+ counts = pd.Series(groups).value_counts()
+ if verbose:
+ worst = int(counts.max()) ** 2
+ if worst > 1e8:
+ print('Computing all-versus-all distances for up to {:.1e} pairs of single '
+ "cells per cell-type pair. Consider method='centroid'.".format(worst))
+
+ df = coords.copy()
+ df[group_col] = groups
+ return pairwise_celltype_distances(df, group_col=group_col,
+ coord_cols=list(coords.columns),
+ method=method, distance=distance, pairs=pairs)
+
+
+def _coordinates_and_groups(adata, group_col, spatial_key, coord_names):
+ '''
+ Normalizes the two accepted inputs into a coordinates dataframe and a list of
+ cell-type labels aligned with it.
+ '''
+ if hasattr(adata, 'obsm'):
+ if group_col not in adata.obs.columns:
+ raise KeyError("'{}' is not a column of adata.obs".format(group_col))
+ coords = get_spatial_coordinates(adata, spatial_key=spatial_key,
+ coord_names=coord_names)
+ groups = np.asarray(adata.obs[group_col].values)
+ elif isinstance(adata, pd.DataFrame):
+ if group_col not in adata.columns:
+ raise KeyError("'{}' is not a column of the dataframe".format(group_col))
+ if coord_names is None:
+ coord_names = [c for c in adata.columns if c != group_col]
+ coords = adata[list(coord_names)]
+ groups = np.asarray(adata[group_col].values)
+ else:
+ raise TypeError('`adata` must be an AnnData object or a pandas DataFrame, got {}'
+ .format(type(adata).__name__))
+ return coords, groups
diff --git a/cell2cell/spatial/neighborhoods.py b/cell2cell/spatial/neighborhoods.py
index 26d3a85..d2cc62d 100644
--- a/cell2cell/spatial/neighborhoods.py
+++ b/cell2cell/spatial/neighborhoods.py
@@ -2,6 +2,8 @@
import numpy as np
import pandas as pd
+from natsort import natsorted
+
def create_spatial_grid(adata, num_bins, copy=False):
"""
@@ -162,9 +164,10 @@ def add_sliding_window_info_to_adata(adata, window_mapping):
"""
# Initialize all window columns to 0.0
- for window in sorted(window_mapping.keys()):
+ for window in natsorted(window_mapping.keys()):
adata.obs[window] = 0.0
# Mark cells that belong to each window
for window, barcode_indeces in window_mapping.items():
- adata.obs.loc[barcode_indeces, window] = 1.0
\ No newline at end of file
+ # Converted to a list because pandas does not accept a set as an indexer
+ adata.obs.loc[list(barcode_indeces), window] = 1.0
\ No newline at end of file
diff --git a/cell2cell/stats/permutation.py b/cell2cell/stats/permutation.py
index 656eb40..0488162 100644
--- a/cell2cell/stats/permutation.py
+++ b/cell2cell/stats/permutation.py
@@ -10,8 +10,9 @@
import seaborn as sns
import cell2cell.core.interaction_space as ispace
-from cell2cell.preprocessing import shuffle_rows_in_df
+from cell2cell.preprocessing import shuffle_rows_in_df, zero_diagonal
+from natsort import natsorted
from sklearn.utils import shuffle
from tqdm import tqdm
@@ -142,7 +143,10 @@ def pvalue_from_dist(obs_value, dist, label='', consider_size=False, comparison=
label_ = label + ' - p-val: <{:g}'.format(float('{:.1g}'.format(1. / len(dist))))
else:
label_ = label + ' - p-val: {0:.2E}'.format(pval)
- fig = sns.distplot(dist, hist=True, kde=True, norm_hist=False, rug=False, label=label_)
+ # `sns.distplot` is deprecated and removed in seaborn 0.14. It normalized the
+ # histogram to a density whenever a KDE was drawn, regardless of `norm_hist`,
+ # so `stat='density'` reproduces what it did here.
+ fig = sns.histplot(dist, kde=True, stat='density', label=label_)
fig.axvline(x=obs_value, color=fig.get_lines()[-1].get_c(), ls='--')
fig.tick_params(axis='both', which='major', labelsize=16)
@@ -196,24 +200,31 @@ def random_switching_ppi_labels(ppi_data, genes=None, random_state=None, interac
prot_b = interaction_columns[1]
if permuted_column == 'both':
if genes is None:
- genes = list(np.unique(ppi_data_[interaction_columns].values.flatten()))
+ # `interaction_columns` is a tuple, which pandas would treat as a single
+ # column name, so it is converted into a list before selecting them.
+ # An object dtype is requested because the values are protein names, which
+ # pandas >= 3.0 returns as an extension array that has no `.ravel()`.
+ genes = list(np.unique(ppi_data_[list(interaction_columns)].to_numpy(dtype=object).ravel()))
else:
- genes = list(set(genes))
+ # Sorted to make the permutation reproducible for a given random_state
+ genes = natsorted(set(genes))
mapper = dict(zip(genes, shuffle(genes, random_state=random_state)))
ppi_data_[prot_a] = ppi_data_[prot_a].apply(lambda x: mapper[x])
ppi_data_[prot_b] = ppi_data_[prot_b].apply(lambda x: mapper[x])
elif permuted_column == 'first':
if genes is None:
- genes = list(np.unique(ppi_data_[prot_a].values.flatten()))
+ genes = list(np.unique(ppi_data_[prot_a].to_numpy(dtype=object)))
else:
- genes = list(set(genes))
+ # Sorted to make the permutation reproducible for a given random_state
+ genes = natsorted(set(genes))
mapper = dict(zip(genes, shuffle(genes, random_state=random_state)))
ppi_data_[prot_a] = ppi_data_[prot_a].apply(lambda x: mapper[x])
elif permuted_column == 'second':
if genes is None:
- genes = list(np.unique(ppi_data_[prot_b].values.flatten()))
+ genes = list(np.unique(ppi_data_[prot_b].to_numpy(dtype=object)))
else:
- genes = list(set(genes))
+ # Sorted to make the permutation reproducible for a given random_state
+ genes = natsorted(set(genes))
mapper = dict(zip(genes, shuffle(genes, random_state=random_state)))
ppi_data_[prot_b] = ppi_data_[prot_b].apply(lambda x: mapper[x])
else: raise ValueError('Not valid option')
@@ -331,9 +342,9 @@ def run_label_permutation(rnaseq_data, ppi_data, genes, analysis_setup, cutoff_s
genes = list(rnaseq_data.index)
if excluded_cells is not None:
- included_cells = sorted(list(set(rnaseq_data.columns) - set(excluded_cells)))
+ included_cells = natsorted(set(rnaseq_data.columns) - set(excluded_cells))
else:
- included_cells = sorted(list(set(rnaseq_data.columns)))
+ included_cells = natsorted(set(rnaseq_data.columns))
rnaseq_data_ = rnaseq_data.loc[genes, included_cells]
@@ -367,12 +378,15 @@ def run_label_permutation(rnaseq_data, ppi_data, genes, analysis_setup, cutoff_s
cci_type=analysis_setup['cci_type'],
verbose=verbose)
+ # The CCI matrix is only filled by this method. Without it, the scores below are the
+ # zeros the interaction space is initialized with.
+ interaction_space.compute_pairwise_cci_scores(verbose=verbose)
+
# Keep scores
cci = interaction_space.interaction_elements['cci_matrix'].loc[included_cells, included_cells]
cci_diag = np.diag(cci).copy()
- np.fill_diagonal(cci.values, 0.0)
- iter_scores = scipy.spatial.distance.squareform(cci)
+ iter_scores = scipy.spatial.distance.squareform(zero_diagonal(cci))
iter_scores = np.reshape(iter_scores, (len(iter_scores), 1)).T
iter_diag = np.reshape(cci_diag, (len(cci_diag), 1)).T
@@ -396,12 +410,13 @@ def run_label_permutation(rnaseq_data, ppi_data, genes, analysis_setup, cutoff_s
cci_type=analysis_setup['cci_type'],
verbose=verbose)
+ base_interaction_space.compute_pairwise_cci_scores(verbose=verbose)
+
# Keep scores
base_cci = base_interaction_space.interaction_elements['cci_matrix'].loc[included_cells, included_cells]
base_cci_diag = np.diag(base_cci).copy()
- np.fill_diagonal(base_cci.values, 0.0)
- base_scores = scipy.spatial.distance.squareform(base_cci)
+ base_scores = scipy.spatial.distance.squareform(zero_diagonal(base_cci))
# P-values
pvals = np.zeros((scores.shape[1], 1))
diff --git a/cell2cell/tensor/external_scores.py b/cell2cell/tensor/external_scores.py
index 18e7995..c075430 100644
--- a/cell2cell/tensor/external_scores.py
+++ b/cell2cell/tensor/external_scores.py
@@ -4,11 +4,30 @@
import pandas as pd
from collections import defaultdict
+from natsort import natsorted
from tqdm import tqdm
from cell2cell.preprocessing.find_elements import get_element_abundances, get_elements_over_fraction
from cell2cell.tensor.tensor import PreBuiltTensor
+def _ordered_intersection(element_lists):
+ '''Intersects multiple lists of elements, keeping the order in which the
+ elements appear in the first list.
+
+ Parameters
+ ----------
+ element_lists : list
+ A list containing lists of elements, one per context.
+
+ Returns
+ -------
+ elements : list
+ Elements present in all lists, ordered as in `element_lists[0]`.
+ '''
+ common = set.intersection(*map(set, element_lists))
+ return [e for e in element_lists[0] if e in common]
+
+
def dataframes_to_tensor(context_df_dict, sender_col, receiver_col, ligand_col, receptor_col, score_col, how='inner',
outer_fraction=0.0, lr_fill=np.nan, cell_fill=np.nan, lr_sep='^', dup_aggregation='max',
context_order=None, order_labels=None, sort_elements=True, device=None):
@@ -123,10 +142,12 @@ def dataframes_to_tensor(context_df_dict, sender_col, receiver_col, ligand_col,
if order_labels is None:
order_labels = ['Contexts', 'Ligand-Receptor Pairs', 'Sender Cells', 'Receiver Cells']
- # Find all existing LR pairs, sender and receiver cells across contexts
- lr_dict = defaultdict(set)
- sender_dict = defaultdict(set)
- receiver_dict = defaultdict(set)
+ # Find all existing LR pairs, sender and receiver cells across contexts.
+ # Lists (instead of sets) are used to keep the order in which elements are found,
+ # making the tensor reproducible across runs when `sort_elements=False`.
+ lr_dict = defaultdict(list)
+ sender_dict = defaultdict(list)
+ receiver_dict = defaultdict(list)
for k, df in cont_dict.items():
df['LRs'] = df.apply(lambda row: row[ligand_col] + lr_sep + row[receptor_col], axis=1)
@@ -136,9 +157,9 @@ def dataframes_to_tensor(context_df_dict, sender_col, receiver_col, ligand_col,
# ccc_df = ccc_df.dropna(how='any')
# lr_dict[k].update(list(ccc_df.index))
# else:
- lr_dict[k].update(df['LRs'].unique().tolist())
- sender_dict[k].update(df[sender_col].unique().tolist())
- receiver_dict[k].update(df[receiver_col].unique().tolist())
+ lr_dict[k] = list(dict.fromkeys(lr_dict[k] + df['LRs'].unique().tolist()))
+ sender_dict[k] = list(dict.fromkeys(sender_dict[k] + df[sender_col].unique().tolist()))
+ receiver_dict[k] = list(dict.fromkeys(receiver_dict[k] + df[receiver_col].unique().tolist()))
# Subset LR pairs, sender and receiver cells given parameter 'how'
df_lrs = [list(lr_dict[k]) for k in context_order]
@@ -146,9 +167,9 @@ def dataframes_to_tensor(context_df_dict, sender_col, receiver_col, ligand_col,
df_receivers = [list(receiver_dict[k]) for k in context_order]
if how == 'inner':
- lr_pairs = list(set.intersection(*map(set, df_lrs)))
- sender_cells = list(set.intersection(*map(set, df_senders)))
- receiver_cells = list(set.intersection(*map(set, df_receivers)))
+ lr_pairs = _ordered_intersection(df_lrs)
+ sender_cells = _ordered_intersection(df_senders)
+ receiver_cells = _ordered_intersection(df_receivers)
elif how == 'outer':
lr_pairs = get_elements_over_fraction(abundance_dict=get_element_abundances(element_lists=df_lrs),
fraction=outer_fraction)
@@ -159,10 +180,10 @@ def dataframes_to_tensor(context_df_dict, sender_col, receiver_col, ligand_col,
elif how == 'outer_lrs':
lr_pairs = get_elements_over_fraction(abundance_dict=get_element_abundances(element_lists=df_lrs),
fraction=outer_fraction)
- sender_cells = list(set.intersection(*map(set, df_senders)))
- receiver_cells = list(set.intersection(*map(set, df_receivers)))
+ sender_cells = _ordered_intersection(df_senders)
+ receiver_cells = _ordered_intersection(df_receivers)
elif how == 'outer_cells':
- lr_pairs = list(set.intersection(*map(set, df_lrs)))
+ lr_pairs = _ordered_intersection(df_lrs)
sender_cells = get_elements_over_fraction(abundance_dict=get_element_abundances(element_lists=df_senders),
fraction=outer_fraction)
receiver_cells = get_elements_over_fraction(abundance_dict=get_element_abundances(element_lists=df_receivers),
@@ -172,10 +193,10 @@ def dataframes_to_tensor(context_df_dict, sender_col, receiver_col, ligand_col,
if sort_elements:
if sort_context:
- context_order = sorted(context_order)
- lr_pairs = sorted(lr_pairs)
- sender_cells = sorted(sender_cells)
- receiver_cells = sorted(receiver_cells)
+ context_order = natsorted(context_order)
+ lr_pairs = natsorted(lr_pairs)
+ sender_cells = natsorted(sender_cells)
+ receiver_cells = natsorted(receiver_cells)
# Build temporal tensor to pass to PreBuiltTensor
tmp_tensor = []
diff --git a/cell2cell/tensor/tensor.py b/cell2cell/tensor/tensor.py
index f4b2890..52d4bac 100644
--- a/cell2cell/tensor/tensor.py
+++ b/cell2cell/tensor/tensor.py
@@ -5,6 +5,7 @@
import tensorly as tl
from collections import OrderedDict
+from natsort import natsorted
from tqdm import tqdm
from cell2cell.core.communication_scores import compute_ccc_matrix, aggregate_ccc_matrices
@@ -1142,12 +1143,12 @@ def build_context_ccc_tensor(rnaseq_matrices, ppi_data, how='inner', outer_fract
if set(df_idxs[0]) == genes:
genes = df_idxs[0]
else:
- genes = sorted(list(genes))
+ genes = natsorted(genes)
if set(df_cols[0]) == cells:
cells = df_cols[0]
else:
- cells = sorted(list(cells))
+ cells = natsorted(cells)
# Filter PPI data for
ppi_data_ = filter_ppi_by_proteins(ppi_data=ppi_data,
diff --git a/cell2cell/tensor/tensor_manipulation.py b/cell2cell/tensor/tensor_manipulation.py
index e62ff77..2d960af 100644
--- a/cell2cell/tensor/tensor_manipulation.py
+++ b/cell2cell/tensor/tensor_manipulation.py
@@ -65,19 +65,24 @@ def concatenate_interaction_tensors(interaction_tensors, axis, order_labels, rem
except:
context = {'dtype': interaction_tensors[0].tensor.dtype, 'device' : None}
- # Concatenate tensors
- concat_tensor = tl.concatenate([tensor.tensor.to('cpu') for tensor in interaction_tensors], axis=axis)
+ # Concatenate tensors. `.to('cpu')` only exists in backends such as pytorch, so
+ # it is skipped for backends whose tensors are numpy arrays (the default one).
+ def to_cpu(data):
+ return data.to('cpu') if hasattr(data, 'to') else data
+
+ concat_tensor = tl.concatenate([to_cpu(tensor.tensor) for tensor in interaction_tensors], axis=axis)
if mask is not None:
assert mask.shape == concat_tensor.shape, "Mask must have the same shape of the concatenated tensor. Here: {}".format(concat_tensor.shape)
else: # Generate a new mask from all previous masks if all are not None
if all([tensor.mask is not None for tensor in interaction_tensors]):
- mask = tl.concatenate([tensor.mask.to('cpu') for tensor in interaction_tensors], axis=axis)
+ mask = tl.concatenate([to_cpu(tensor.mask) for tensor in interaction_tensors], axis=axis)
else:
mask = None
- concat_tensor = tl.tensor(concat_tensor, device=context['device'])
+ # The context of a numpy-backed tensor does not include a 'device' key
+ concat_tensor = tl.tensor(concat_tensor, device=context.get('device', None))
if mask is not None:
- mask = tl.tensor(mask, device=context['device'])
+ mask = tl.tensor(mask, device=context.get('device', None))
# Concatenate names of elements for the given axis but keep the others as in one tensor
order_names = []
diff --git a/ci/constraints-py3.10.txt b/ci/constraints-py3.10.txt
new file mode 100644
index 0000000..c8d058f
--- /dev/null
+++ b/ci/constraints-py3.10.txt
@@ -0,0 +1,65 @@
+# Versions the CI job resolved on Python 3.10, pinned so a new release of a
+# dependency cannot turn a build red without a commit to this repository.
+#
+# This constrains CI only. setup.py stays open, so nothing here limits what a
+# user installing cell2cell gets.
+#
+# To move to newer dependencies, regenerate on Python 3.10 and commit the diff:
+# python -m venv /tmp/c && /tmp/c/bin/pip install -e '.[test]'
+# /tmp/c/bin/pip freeze --exclude-editable > ci/constraints-py3.10.txt
+# then re-add this header. Review the diff and run the suite before pushing.
+anndata==0.11.4
+array-api-compat==1.15.0
+certifi==2026.7.22
+charset-normalizer==3.4.9
+contourpy==1.3.2
+coverage==7.15.3
+cycler==0.12.1
+et_xmlfile==2.0.0
+exceptiongroup==1.3.1
+fonttools==4.63.0
+gseapy==1.3.1
+h5py==3.16.0
+idna==3.18
+iniconfig==2.3.0
+joblib==1.5.3
+kiwisolver==1.5.0
+kneed==0.8.6
+legacy-api-wrap==1.5
+llvmlite==0.48.0
+matplotlib==3.10.9
+natsort==8.4.0
+networkx==3.4.2
+numba==0.66.0
+numpy==2.2.6
+openpyxl==3.1.5
+packaging==26.3
+pandas==2.3.3
+patsy==1.0.2
+pillow==12.3.0
+pluggy==1.6.0
+Pygments==2.20.0
+pynndescent==0.6.0
+pyparsing==3.3.2
+pytest-cov==7.1.0
+pytest==9.1.1
+python-dateutil==2.9.0.post0
+pytz==2026.3.post1
+requests==2.34.2
+scanpy==1.11.5
+scikit-learn==1.7.2
+scipy==1.15.3
+seaborn==0.13.2
+session-info2==0.4.1
+six==1.17.0
+statannotations==0.7.2
+statsmodels==0.14.6
+tensorly==0.9.0
+threadpoolctl==3.6.0
+tomli==2.4.1
+tqdm==4.70.0
+typing_extensions==4.16.0
+tzdata==2026.3
+umap-learn==0.5.12
+urllib3==2.7.0
+xlrd==2.0.2
diff --git a/ci/constraints-py3.12.txt b/ci/constraints-py3.12.txt
new file mode 100644
index 0000000..6e9fc63
--- /dev/null
+++ b/ci/constraints-py3.12.txt
@@ -0,0 +1,75 @@
+# Versions the CI job resolved on Python 3.12, pinned so a new release of a
+# dependency cannot turn a build red without a commit to this repository.
+#
+# This constrains CI only. setup.py stays open, so nothing here limits what a
+# user installing cell2cell gets.
+#
+# To move to newer dependencies, regenerate on Python 3.12 and commit the diff:
+# python -m venv /tmp/c && /tmp/c/bin/pip install -e '.[test]'
+# /tmp/c/bin/pip freeze --exclude-editable > ci/constraints-py3.12.txt
+# then re-add this header. Review the diff and run the suite before pushing.
+anndata==0.13.2
+annotated-types==0.8.0
+array-api-compat==1.15.0
+certifi==2026.7.22
+charset-normalizer==3.4.9
+contourpy==1.3.3
+coverage==7.15.3
+cycler==0.12.1
+donfig==0.8.1.post1
+et_xmlfile==2.0.0
+fast-array-utils==1.5
+fonttools==4.63.0
+google-crc32c==1.8.0
+gseapy==1.3.1
+h5py==3.16.0
+idna==3.18
+iniconfig==2.3.0
+joblib==1.5.3
+kiwisolver==1.5.0
+kneed==0.8.6
+legacy-api-wrap==1.5
+llvmlite==0.48.0
+matplotlib==3.11.1
+narwhals==2.24.0
+natsort==8.4.0
+networkx==3.6.1
+numba==0.66.0
+numcodecs==0.16.5
+numpy==2.4.6
+openpyxl==3.1.5
+packaging==26.3
+pandas==3.0.5
+patsy==1.0.2
+pillow==12.3.0
+pluggy==1.6.0
+pydantic_core==2.46.4
+pydantic-settings==2.14.2
+pydantic==2.13.4
+Pygments==2.20.0
+pynndescent==0.6.0
+pyparsing==3.3.2
+pytest-cov==7.1.0
+pytest==9.1.1
+python-dateutil==2.9.0.post0
+python-dotenv==1.2.2
+PyYAML==6.0.3
+requests==2.34.2
+scanpy==1.12.3
+scikit-learn==1.9.0
+scipy==1.18.0
+scverse-misc==0.1.3
+seaborn==0.13.2
+session-info2==0.4.2
+six==1.17.0
+statannotations==0.7.2
+statsmodels==0.14.6
+tensorly==0.9.0
+threadpoolctl==3.6.0
+tqdm==4.70.0
+typing_extensions==4.16.0
+typing-inspection==0.4.2
+umap-learn==0.5.12
+urllib3==2.7.0
+xlrd==2.0.2
+zarr==3.3.0
diff --git a/docs/index.md b/docs/index.md
index 61aa7bd..8e370a6 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -42,7 +42,7 @@ pip install cell2cell
| cell2cell Examples | Tensor-cell2cell Examples |
|---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
|  |  |
-| - [Step-by-step Pipeline](https://github.com/earmingol/cell2cell/blob/master/examples/cell2cell/Toy-Example.ipynb) - [Interaction Pipeline for Bulk Data](./tutorials/Toy-Example-BulkPipeline) - [Interaction Pipeline for Single-Cell Data](./tutorials/Toy-Example-SingleCellPipeline) - [Whole Body of *C. elegans*](https://github.com/LewisLabUCSD/Celegans-cell2cell) | - [Obtaining patterns of cell-cell communication](./tutorials/ASD/01-Tensor-Factorization-ASD/) - [Downstream 1: Factor-specific analyses](./tutorials/ASD/02-Factor-Specific-ASD/) - [Downstream 2: Patterns to functions (GSEA)](./tutorials/ASD/03-GSEA-ASD/) - [Tensor-cell2cell in Google Colab (**GPU**)](https://colab.research.google.com/drive/1T6MUoxafTHYhjvenDbEtQoveIlHT2U6_?usp=sharing) - [Communication patterns in **Spatial Transcriptomics**](./tutorials/Tensor-cell2cell-Spatial/) - [Multi-modal communication patterns with **Coupled Tensor Component Analysis**](./tutorials/Version2/Tensor-cell2cell-CTCA/) |
+| - [Step-by-step Pipeline](https://github.com/earmingol/cell2cell/blob/master/examples/cell2cell/Toy-Example.ipynb) - [Interaction Pipeline for Bulk Data](./tutorials/Toy-Example-BulkPipeline) - [Interaction Pipeline for Single-Cell Data](./tutorials/Toy-Example-SingleCellPipeline) - [Whole Body of *C. elegans*](https://github.com/LewisLabUCSD/Celegans-cell2cell) - [Identifying a **spatial code** of ligand-receptor pairs (Genetic Algorithm)](./tutorials/Genetic-Algorithm-LR-Selection/) | - [Obtaining patterns of cell-cell communication](./tutorials/ASD/01-Tensor-Factorization-ASD/) - [Downstream 1: Factor-specific analyses](./tutorials/ASD/02-Factor-Specific-ASD/) - [Downstream 2: Patterns to functions (GSEA)](./tutorials/ASD/03-GSEA-ASD/) - [Tensor-cell2cell in Google Colab (**GPU**)](https://colab.research.google.com/drive/1T6MUoxafTHYhjvenDbEtQoveIlHT2U6_?usp=sharing) - [Communication patterns in **Spatial Transcriptomics**](./tutorials/Tensor-cell2cell-Spatial/) - [Multi-modal communication patterns with **Coupled Tensor Component Analysis**](./tutorials/Version2/Tensor-cell2cell-CTCA/) |
Reproducible runs of the analyses in the [Tensor-cell2cell paper](https://doi.org/10.1038/s41467-022-31369-2) are available at [CodeOcean.com](https://doi.org/10.24433/CO.0051950.v2)
diff --git a/docs/tutorials/Genetic-Algorithm-LR-Selection.ipynb b/docs/tutorials/Genetic-Algorithm-LR-Selection.ipynb
new file mode 100644
index 0000000..3372e65
--- /dev/null
+++ b/docs/tutorials/Genetic-Algorithm-LR-Selection.ipynb
@@ -0,0 +1,1714 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "3984da0a",
+ "metadata": {},
+ "source": [
+ "# Identifying a spatial code of ligand-receptor pairs with a genetic algorithm\n",
+ "\n",
+ "Cells that sit close together in a tissue should be talking to each other — but only some of the\n",
+ "ligand-receptor pairs in a resource carry that spatial signal, and the rest dilute it. This tutorial\n",
+ "searches for the subset whose cell-cell interaction scores best reproduce **how far apart the cells\n",
+ "actually are**, using a genetic algorithm.\n",
+ "\n",
+ "The pairs it recovers are a *spatial code*: the interactions whose activity tracks physical\n",
+ "proximity, and which can therefore be read as a signature of where cells are relative to one another.\n",
+ "\n",
+ "This is the analysis of\n",
+ "[Armingol et al. (2022)](https://doi.org/10.1371/journal.pcbi.1010715) on the whole body of\n",
+ "*C. elegans*, where the reference was a digital 3D map of the animal\n",
+ "([Celegans-cell2cell](https://github.com/LewisLabUCSD/Celegans-cell2cell)). It is now part of the\n",
+ "package as `cell2cell.analysis.optimize_lr_pairs`, and §6 reproduces it end to end.\n",
+ "\n",
+ "The reference does not have to be a physical distance — any square matrix of distances between the\n",
+ "same cells works, so the same search applies to developmental, functional or phenotypic similarity.\n",
+ "Here we start from spatial transcriptomics, using the helpers that turn single-cell coordinates into\n",
+ "cell-type distances, and finish on the *C. elegans* 3D map.\n",
+ "\n",
+ "**Requirements.** The genetic algorithm needs the optional dependency `pygad`:\n",
+ "\n",
+ "```\n",
+ "pip install cell2cell[ga]\n",
+ "```"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "44fd2f26",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-08-23T22:44:19.971590Z",
+ "iopub.status.busy": "2026-08-23T22:44:19.971388Z",
+ "iopub.status.idle": "2026-08-23T22:44:26.339568Z",
+ "shell.execute_reply": "2026-08-23T22:44:26.339079Z"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/private/tmp/claude-503/-Users-eg22-Repos-cell2cell/fa1914c8-fe56-40d2-b837-7db6feb472f8/scratchpad/venv-ci312/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n",
+ " from .autonotebook import tqdm as notebook_tqdm\n"
+ ]
+ }
+ ],
+ "source": [
+ "import numpy as np\n",
+ "import pandas as pd\n",
+ "import matplotlib.pyplot as plt\n",
+ "\n",
+ "import cell2cell as c2c\n",
+ "\n",
+ "%matplotlib inline"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "21a766fe",
+ "metadata": {},
+ "source": [
+ "## 1. Data\n",
+ "\n",
+ "We build a small synthetic system where we know the answer in advance, so the search can be judged.\n",
+ "\n",
+ "Twelve cell types sit at different positions along an axis. Around each one we scatter single cells,\n",
+ "which is what a spatial transcriptomics experiment would give us."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "927a32e7",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-08-23T22:44:26.341694Z",
+ "iopub.status.busy": "2026-08-23T22:44:26.341470Z",
+ "iopub.status.idle": "2026-08-23T22:44:26.346324Z",
+ "shell.execute_reply": "2026-08-23T22:44:26.345922Z"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "AnnData object with n_obs × n_vars = 480 × 1\n",
+ " obs: 'celltype'\n",
+ " obsm: 'spatial'\n",
+ " layers: None (.X)"
+ ]
+ },
+ "execution_count": 2,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "import anndata\n",
+ "\n",
+ "rng = np.random.default_rng(0)\n",
+ "\n",
+ "n_celltypes = 12\n",
+ "cells_per_type = 40\n",
+ "positions = np.linspace(0, 100, n_celltypes)\n",
+ "celltypes = ['CT-{}'.format(i + 1) for i in range(n_celltypes)]\n",
+ "\n",
+ "coords, labels = [], []\n",
+ "for name, x in zip(celltypes, positions):\n",
+ " coords.append(np.column_stack([rng.normal(x, 2.5, cells_per_type),\n",
+ " rng.normal(0, 2.5, cells_per_type)]))\n",
+ " labels += [name] * cells_per_type\n",
+ "\n",
+ "coords = np.vstack(coords)\n",
+ "adata = anndata.AnnData(X=np.zeros((len(coords), 1), dtype=float),\n",
+ " obs=pd.DataFrame({'celltype': labels},\n",
+ " index=['cell-{}'.format(i) for i in range(len(coords))]))\n",
+ "adata.obsm['spatial'] = coords\n",
+ "adata"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "bf559184",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-08-23T22:44:26.347622Z",
+ "iopub.status.busy": "2026-08-23T22:44:26.347532Z",
+ "iopub.status.idle": "2026-08-23T22:44:26.444169Z",
+ "shell.execute_reply": "2026-08-23T22:44:26.443809Z"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwgAAAFICAYAAAD0wtlSAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAjvFJREFUeJztnQeYFFXWhs/kGWbISBJEARURMIGCCIqKrJjTmtc1J1RQjL+65oiKCRUjZnSNa9gVFQQUVIwgZhBBkQGGNDnV/3y3u5rq6qruqurqON/7PENxb1dXV1enc+453zk5mqZpQgghhBBCCCEiksurQAghhBBCCNGhg0AIIYQQQggJQQeBEEIIIYQQEoIOAiGEEEIIISQEHQRCCCGEEEJICDoIhBBCCCGEkBB0EAghhBBCCCEh6CAQQgghhBBCQtBBIISkDUOHDpWDDjoo6x7Lj3NL1/NN1/NySlVVleTn58stt9wSdb9PP/1U7ffmm29KpvDjjz+qc37++eejzhFCiBk6CISQhNLQ0CCPPvqoDBs2TDp37iydOnVS/7/99ttl9erVYfs2NjZKU1NTUl6RZD6WH+eWruebruflFE3T1Pk3Nzf7sl86YXXOTp/HN998oxyJl156KQlnSghJN/JTfQKEkOzm+OOPV6uukyZNkqOOOkpKS0vlf//7n1x66aUyf/58ee2118JWaXNyclJ6voSQzHSICCH+wQgCISRhfP/99/LKK6/IueeeKxdeeKFsueWW0q5dOzn22GPl66+/lt133z1s/7y8PMnN5dcSIYQQkkr4S0wISRhr165V2x49ekTc1qFDB7nyyitj5rP369dPORTInT7ggAOkTZs2svXWW8sdd9wRcUyseN5www2yzTbbqP323XdfWbRokYpcDBgwwNE5v/322+p+7du3V8fYZ599ZNasWY7vO3r0aOnYsaNKpRo7dqx89tlnvh3fDO6Hx+vSpYt6zJEjRyqHzK9zffjhh2WXXXaRsrIylR6G6/jdd985Or6T++K2Sy65JOK+hx12mOy8886W7wMcA++Dtm3bykUXXaRuq62tleuvv1769++vIlTdu3eXM844Q8rLy8OOsXz5cvn73/+urj2u2QUXXCB1dXXilkceeUT69u2rzn+vvfaSuXPnhm5D2hzOwep5rV+/Xr3m48aNi3p8PJ8bb7xRBg4cqB6jT58+Mn78eKmoqAjbx8lz9gIifoMHD1b/P/HEE1WqEf7w2XLz/Iy6jWjXLBnPiRDiDjoIhJCEAaO8devW8uSTT8qSJUs8597DQLjiiivktttuk6VLl8rZZ58tl19+ubz66qth+yJKAbHp//3f/6nHQ1rTZZddpowaHCcWDzzwgBxyyCEyfPhwFeH46aefZMSIEbL//vvLjBkzot538uTJ6r4w6ubNm6fui/OB1sKP45v5+eef5cADD1TG1JdffqmeL67Piy++KMuWLYv7XJECBkPv5JNPVsf+8MMPlYEKJw5OVzSc3tdOv4A58+uF8apVq2TixInqNV68eLG6jvX19TJmzBiZOnWq3HzzzcoJ+O9//6tuh8NUWVmp7r9hwwbZe++91eO/88476vYhQ4Yow9sNzzzzjPz666/y0UcfqdcQjgZev88//1zdvsUWWygHB+/5mpqasPtibtOmTfLPf/7T9vhwWHC8+++/XznQuH7vv/++9O7dWx5//HG1j9Pn7BW8N2Dc688Xhjv+rr76alfPT09TmjZtWtRrloznRAhxiUYIIQnkzTff1LbYYgstNzdXGzx4sHbWWWdpzzzzjLZ69eqIfXfbbTdtzJgxYXN9+vTRioqKtBUrVoTN9+vXTzvggANC46VLl6rHuPzyy8P2w3x+fr62/fbbR32s8vJyrbi4WDv11FMjzguPs8suu9g+x5UrV6pzPPHEE233cXN8q+tgnnvqqac0fIX//vvvto/p9Vx/++03dS3POeecsPl169Zpbdq00Q455BDb83JzX5zHRRddFPH4Bx10kLbjjjtGvA/y8vLU62lkypQp6jp88MEHYfPLly9Xx7/rrrvU+JZbblH7ff3112H73XDDDWr+xhtv1KIxb948td++++4bNl9bW6t1795d23///UNzc+fOVfs+/vjjobnm5mZtu+220wYNGhT1ce677z5133fffdd2H6fP+fvvv1f74fOmYzVnxVdffaX2e+GFFyJuc/r83Fwzp8+JEJIcGEEghCQUrEZiNfA///mPSh9auXKligAgDeiJJ55wdIxdd91V6ReMYPXbGJVAug0ElYceemjYfkhHwr6xeO+999QqKVJQzCAV56uvvlKr0FZg9R8rv1gxT8Tx7aIzEHSfeeaZ8sEHHzhOlXFyrjNnzlTXEmlBRqAfwbni8RJx31gg7QivpxGkryBNCmlbRpDWhrQkPX0Lj4v33E477RS235FHHunqHLB6bqSoqEilZ2F1HBW7ACIbeJyHHnootB8eH5Ga008/Perx33rrLZXy9be//c12H6fPOVG4fX5OrlmqnxMhJBxWMSKEJBzdIMAfQKoInAU4CnvuuacyAKLRrVu3iDnkOiPnWWfNmjVqi/QFM5hDalI0/vrrr5AxA8Mb6RH4A3oaDDQVyH03o+dIm50Yv45vxW677abSiW699VaVrlFYWKgMNwjCjznmGNv7OTlXXTtidd27du0q1dXVqn8A8sT9vK+Ofl3MWJ0zrivSl4qLi0P31f9wXVu1ahV6f1i9N3BObrA7BgzdjRs3KiMXnHfeeer9jTQapDJNmTJFfQ6Q0x8NvD7RXhs3zzmRuHl+Tq5ZOjwnQshmGEEghCQdGAzIUUde+SeffBJzf7vSp0ZDEqJnYCVodCJy1A07rLAj3xlGLIxZ/GHFHcYM8sCtQF42+PPPPxNyfDsQjUDkARqL6dOnKwMNc9Fq1zs5V/1awpEzgzkYcXYGm5v7whmyyi23O7eCggLL69qrVy91HON1RX488trnzJkTOi+r94HVeUbD7hgQ48Jp1YGhjOeHVfY//vhDRdAOP/zw0Psg2usT7bVx85wTiZvn5+SapcNzIoRshg4CISRhLFiwwLZjK1KNgNMV81hAgApH4t133w2bX7FihSxcuDDm/ZH+glV49GXQq7aY/+zQV/BfeOGFhBw/FkhJgXH2xhtvqDKxSN2I51z1a4njGYEAFakko0aNsnXa3NwXDpFZ8AyBtdNKSeDggw+W3377TYlfra4pSucCPC6iSBC9GnHbGRkpQEbg2EFMi8o8RgcGEZJTTjlFRXkg/oYzHCu9CCDKBocvWiqW0+ccD3qEB8a53e1On5+Ta5aM50QIcUGStA6EkBbInDlzlPDw2GOP1T777DOtvr5eW7t2rTZt2jStrKxM22GHHbTq6uqYIuWjjjoq4tinn3661rFjx7C5M844Q2vVqpX23HPPaZWVlUqQeeihh2p77rlnTJEygBASAttrr71WW7ZsmTrfX375RXv00Ue1k046Kepzve2229R9r7zySiUcrqqq0mbNmqWdcMIJro/vRKR8//33a5dddpkS3eIabty4UbvnnntshaVuz/X888/XCgsLtYcffljbsGGDOs+//e1vWklJifbll19GPVen933sscfU+U6dOlWdwzfffKMdfvjh6vWyEilbvQ9qamq0YcOGab169dJeffVVraKiQtu0aZP2xRdfaBdffHFISIv33ZZbbqntvPPO6prh8V5++WX1nN2IlCGyvu6667Q1a9Yo4TzE3hDB471u5ocfflD3wR/Or6mpKepj6M9nyJAhWrdu3dTzwfWDsBzXCK+vm+ccj0gZQuLWrVtrJ598sno8K2I9PzfXzOlzIoQkBzoIhJCE0djYqL3zzjvKCIOBV1BQoAxHVCCaOHGiMtqMxOsgNDQ0KKO3c+fO6rGGDh2qDNIDDzwwonqM1WOBt956S1VdgXGEY/Tt21dV5IExFIuXXnpJPSaqrqBiD44ze/Zs18d34iCsX79emzRpkqp+BMMbj4fHfv7552Oep5NzhbEHhwaOFaoHlZaWKiMfBlus6+j0vtgPr1enTp3UtRg5cqS6DnZVjKzeB7pxiWpEcDhxHDwfVMyaPHmyMjJ14KgcfPDBqpoUrj8qSq1atcqVg/Daa6+p644qPHhuuP7vvfee7f32228/dT8YyE6Bc4tqXFtvvbV6jB49eqj3CM7VzXOOx0EAcDRx3eFM4jyuv/56V8/P7TVz+joSQhJPDv5xE3EghBCv4OvGLjUFQIyI243dlK3mACrl4HhOUg8GDRqkRLP/+9//Yh43HXBzHVJNIs7L6rVN5PNHegweK9p7081+RiAYR78OVNxCjn2ysTpnt89DFwrj2puvf7TnN3/+fBk2bJhKq0MKHCEkc2AVI0JI0ohlkFgZ+3YOgFNDER2YkdOOLrxOjpsOuLkOqSYR52X12iby+TvVf7jViUBsi6Zs0J+kwjmwO2e3zwOfW6v7pMPzI4QkhvRaiiKEkDiAGPLOO+9UYkf0HECFJKxwdu7cWZVkJCRZYJUeHZ9Rieeqq67Kuguf7c+PkJYOHQRCSNaASjWopX7AAQdI69atVWUU9FiYO3euqvRDSDLQewI89thjctddd8nIkSOz6sJn+/MjhIhQg0AIIYT4CHL2oaNI17SwZD4/L7oNQkjqoYNACCGEEEIICcEUI0IIIYQQQkgIVjEygbAp2twjf5khUUIIIYQQkk0pguhs371796jVAOkgmIBz0LNnz0S/PoQQQgghhKSE5cuXS48ePWxvp4NgApED/cK1adMmsa8OIYQQQgghSWLjxo1qIVy3d+2gg2BCTyuCc0AHgRBCCCGEZBux0ugpUiaEEEIIIYSEoINACCGEEEIICUEHgRBCCCGEEBKCDgIhGURdU50s37RcbQkhhBBCEgFFyoRkCPNXzpfxM8dLVUOVlBWUyT2j7pGh3Yam+rQIIYQQkmUwgkBIBoCIAZyD6oZqNYaTMGHmBEYSCCGEEOI7dBAIyQDKq8uVU6CJpsbYVjZUqnlCCCGEED+hg0BIBtC5VWeVVpQjgbrF2GKMeUIIIYQQP6GDQEgGUJRXpDQHpQWlaowtxpgnhBBCCPETipQJyRAgSJ517CyVVoTIAZ0DQgghhCQCOgiEZBBwCnq27pnq0yCEEEJIFsMUI0IIIYQQQkgIOgiEEEIIIYSQEHQQCCGEEEIIISHoIBBCCCGEEEJC0EEghBBCCCGEhKCDQAghhBBCCAlBB4EQQgghxCO1Tc2yrKZObUma01ArUrE0sCVRYR8EQgghhBAPzKnYJKctWiqbmpqldV6uPDFgGxnRoTWvZTqyZJbI9JNE6jaJFLUWOfZZkd77pPqs0hZGEAghhBBCXIKIAZyDymDkoDI4ZiQhDUHEQDkHlYExthgzkmALHQRCCCGEEJesqm9QkQMtOMYWY8yTNGPTykDkwPhqYYx5YgkdBEIymLqmOlm+abnaEkIISR5dCgtUWlFOcIwtxpgnaUbrboG0IuOrhTHmiSV0EAjJUOavnC97T99bxr46VvaZvo8aE0JIImloaJCKigq1bekUBzUHZXkBU6osOMY8STMKigOag6KywBhbjDFPLMnRNE2PtxAR2bhxo7Rt21Y2bNggbdq04TUhaQkiBnAOqhuqRRNNciRHSgtKZdaxs6QoryjVp0cIyUKWLFki06dPl7q6OikqKpJjjz1WevfuLS2d2mBaESIH2eoc1DY0SfnGOuncpkiKC/IkY4HmAGlFiBy0UOdgo0M7NzvfyYRkOeXV5VLVUKWcA4BtZUOlmieEEL9BxEB3DgC2GDOSEIgk9Copylrn4ONf1sjgm96XkXfOVFuMMxY4BR22abHOgRuy891MSJbTuVVnKSsoU5EDgC3GmCeEEL/ZtGlTyDnQwRjzJHtB5ODsZ76QqvpGNcYWY8yT7Car+iA0NTXJzJkzI+Z33HFH6daNQhSSPSCN6J5R98iEmRNU5ADpRRgzvYgQkghat26t0oqMTgLGmCfZm+KEtKLKuoBzAJCUjjHmt+rYytfHIulFVjkINTU1Mnr0aBkyZEhYXtVll11GB4FkHUO7DVWaA6QVIXJA54AQkigKCgqU5sCsQcA8yd5GbdAclBXlq8gBnIOcHJHSwnw1T7KbrBIpV1ZWqtWMefPmydChQz0dgyLlxJE1IidCDDQ2NEnV+nopbVco+XxfJ5TG+nqpXFchZe07SH5hId+HKQCaA6QV4beWzkHqIwcDP16kGrTBkMsJVlJaOHyAr5EEaA6QVoTIAZyFR07eTYb37eTb8UlycWrnZlUEwVhpAelGffr0ka5du6b6dFoM0RwAfsGQbGT5DxXy34cXSn1tkxQW58nfzhkoPft1SPVpZSXLFn4tb951i9TXVEthSSs59JKrpNfAnVN9Wi0OOAUdOvA9nk6N2sSiURtE034BZ2DB1ftzga+FkZUi5csvv1zGjx8v22yzjRx22GGydu1a230RKoU3Zfwj/lY5oMiJZGvkQDkHdQGxHrYYY574fK3r6wPOQW1N4FrX1qgx5glpqSSzURsW/aA5YPS/5ZBVDkJ+fr689NJLsnz5cvn888/lp59+ksWLF8u5555re59bb71VhVr0v549eyb1nLOBWA6ALnLSk9mMIidCMhWkFSFyEKw0q7YYY574C9KKEDkwfolgjHlCWips1EYSSVY5CMXFxXLMMceExjD2IVB+/fXXbWs1X3nllSoPS/+Dc0HcEcsB0EVOEDcBbDGmyIlkMtAcIK3IuHyHMeaJv0BzgLQi45cIxpgnpCUDQTI0B58O3UFt/RQok5ZNVjkIVnTq1Ek5B3ZpRqjEAJGG8S9bwYr+72urXdcvjnW/WA4AQpIQNZUWBiQv2GLMUCXJZCBIhuagsCigt8EWYwqVE3CtCwuV5qCwuCRwrYtL1DiWUFlraJaGVVXSsKpa/Z94B7+jFRUVbIyWhmR7ozaSGrKqihEiAEgTMnLmmWfKW2+9JX/++afk6BZsC6xi5FUk7PR+TvZjFSOSjbCKUXpWMar9ZZ2smbZYJOgY5BTkSsdT+ktx3/ZJOtvsKvxhLm/au3fvVJ8WIZlBQ63IppUirbulRQdnp3ZuVjkIDz74oMyYMUMOPfRQadeunbzzzjvy9NNPy5NPPiknnniio2Nko4MAwxzCYb2OMSgtzJMvrhkddRXffD+9/jGqGVjdjw4AISQdQLTgjxvnidSHRw1yivKk+9VDlbOQCaRDSVGcw6RJkyIapE2cOJFlTklUaBOIyJJZItNPEqnbJFLUWuTYZ0V675PSd06LLHN6/vnny3bbbScvv/yyrF69Wq1wfPPNN7LDDjtIS8bcCRFU1TfJuwtXyhG79vCtg6Je5YAQQlJJE/RPJucAaHVN6rb8joFUpXQmXVbt4aAYnQOAMeZZ7pTYwdLmEogcKOegMvjBqQyMJ/6SFpGEWGTGMooL0El56tSp8tprr8ldd93V4p0DAC0AIgZmrnnju6h6BIqLCSGZSB70T4WRP2+IIKjb0hys2uvOAcAWY7tiG4kE0Qs4KEYwxny60dRUJzU1v6stSR0sbR4EaUWIHBhL3WGM+Qwg6xwEIpYr+zcdPiBiPlapUYqL05u6pjpZvmm52pL00SNsWF3DXggpBilEnf7RX8SQSqQ0CCfvkBHpRdFW7ZMNUpsQvdCdBD2akW5dlCsqPpY5c3eXT+aNUluMSWpgafMg0BwgrchY6g5jzGcAWZViROw5cGA3ufr1RSq1COh6glilRtlBMT2Zv3K+jJ85XqoaqqSsoEzuGXWPDO02NNWn1aJhV+X0AmLkLa8dJo0VaK6WI/kdijPCOTCu2pvz/jGfCl0CUpugOUi1HsIORAy+XXieNDVVBcdVajxir88kLy/9I0bZhp59YNYvtrjS5gXFAc1BSINQFhhnQHoRyIxvSxI3XyxbJ80GOXpRfq5lqVGrkqbsoJheIGIA56C6oVqN4SRMmDmBkYQUwq7KputRXy/rV/2V8k7HcAgKupRKQZdWGeMcRFu1R58eCIbvu+8+tYVOIZnnBM1BujkHoL5+lTQ1VYalcmCMeZJ8mH1gAIJkaA4u/DqwTbFA2Q2MILSgfMDaxmD0AC98bq7s1iu83B9FRZlBeXW5cgp0NNGksqFSzfdszU7gKe2qLJFdldtukf6CWD9ZtvBrefOuW1SnYzQzQ7+CXgN3TvVpZRzmVXtgrCak6xJYTUiksLCL5OWVBSMIcBJyJC+vVM2T1MDsAwOIGHTYRjKNzFlSIf7lA1roDxIpKvLaoI1Y07lVZ5VWlBPMa8QWY8yT1MCuygEQMVDOQS3SeuAk1ahxqiMJTkujNq6tSauGasZV+3TSJaQbSCMaNHCKcgoC49Lg2H1KS21TsyyrqVNbEh/MPshsGEFoATjJB3Rb0tQpaRmVSLOmJW4pyitSmgOkFSFyUFpQqsaYd5qihGgDHAqn9yHOuir/9+GFKnLQUrsqo4kZIgchNE2NMd+uS1fXzc6SBZqqrX3me1UGFZWOIGZOt4Zq0XQJLQ1oDpA+FIgcBL7DOnQYrjQH5nk3zKnYJKctWiqbmpqldV6uPDFgGxnRoeVdX0KyrlGaH2RjozQnhrrbpmj6feBAwNGwa5zm9pgtsWmJV2IZ+la3U9ycnK7KRaX5UlfVqCILLclJgPH/0FknBSIIwQ99YXGJnDv1WeUMpGP6ESIGf940XzQUcAhkp0hOYXo2VEuX3gipBNWJAoLkSpVWhEgBnIN4QcRg4MeLpLKpWX8bSFleriwcPkCK89LrfUBIPLTITsot3UGwMtiNcyCaQe9mtd/JvkgrGnnnzIj7zr50lH8N1dxEA7DvpL7BpiXBnwBUFciQpiVusHIEdum8i+w9fW8lboZuAalJiD7MOnYWIwk+0tKrGdk5AbGch1SBtKK/7lwQMd/10sFp2VAtHborJysqYLUPSpiatQaxqhU5OTbSivaY/33E/KdDd5BeJYy0kuyhRXZSbslYGezATXqPU1GRnV7BHBlIeKkzt9GAUNMSiWxakoECIrdVjp4d+yzFzSmqZnTqnXu1mEgCnAEY/eY0og3lf8VMP4q2yo8OyGhy5veqvjpmUV5EBCFdG6rpuoRsiwp88+250txcJbm5pbLToIcsowKbqxVJWLWijRu/kjZtdrE0/p1GHLoUFqi0InMEAfOEtEQYN8sCrAz2s55e4El07ERU5LQJSlylzrDaX7E0sHXTwtxu/yxoWuK2yhGiBMYqRyAV4ubmujqpX75cbVtMNSMtsppRSwJOAYx+3TlAVOG5/7skfCdEEEpaKScilj4AKUBY5ccWYz8JNVALdpvHNpkN1RARqKioSEiX5EQe2y+wuv/Nt2cp5wBgi7FVN2S9WtHm7/AAX351omVzNLv+CFbHLg5qDuAUgLLgmOlFpKXCCEIWYCUw1hui+S06dhsZ8FTqzElkwEs0IMOblritcqQ7CXoqUY/WPcLEza3yW8lVe1yV0HOpmjdPVlxwoTRXVkpOq1ay5eR7pPXIkZLt1YxUBCG4DAnBMuZbKrWVlfLGnTdJQ124815YXKzSj6KlFyFyoMTDwe8zbDH2Wx8AQTKOmagoRSo0BZmiV6itXS7NzeHvDYwxX1ra17JakR4RCJATjCRENkezizhgvqRkq4hzgSAZmoNV9Q0qckDnwOVrGUOX6Nd9SHJgBCEL0A12GOogYLDnRcxh7Ed6j9vIgKtSZ04jAxHRAFgcpSIl7bO2aYnbKkdwCoCxyhG6LUNzcMtet6hLd9Xcq2Sf6fsozYLfIGKgnIOqwOqdVl0tK846Wypnz5Z0SgnasLpGbf2sZgSnAGRDNaN4mp4hcjD1/H9GOAfg+JvujilQhsGOykLGiAzGmPcbOAXQHCQzcqAb8Ma+Bn6s9ify2H5jp4K0m9erFe26y3P6nrbN0SIjDtAslEXtjwCnAJqDVDoHmVhqFWnOKEoC3SG2GCfiPq6zDYhnGEHIAnSDXdcblAYNdmCe88tDT1gTFKeRAXM0ANRXiUweEFuLkKFNS9ygOwJ2VY5u/vTmCI2C32LlxvJyFTkws2L8BNlu3ieSG+wSm21iYhwDmgOkFWV6FaN4qg7pPREaLFLLcKx2nbtknT7ADdH6GsSrMUjksf2mpKSn5OaWSHNzoHcGwBjzdiBCENAcRG+OZo44xNMfIVlkYqlVp7rEeO+TzVUJ0xFGELIE3WBHhSBsMbaaS/smKG50AvgiGL8oEDnQ93eiRWghqw8w9tFZ2Wz022kUMO8n+Z07q7QiM4gkwHlIRzGxn5EEdFBOF+fASxQgnqZn2OfPn38MiJJNS8EFRbFTi6LpA4oO7ipN2uaUykxF72tgxK++Bok8tt/AWN9p0CPBlX6My4LjIl+ao+kRhz2HzVRbP0qiJgpEDOAcQCgNKoPjdI8kONUlxnufuHWIxBWMIGQRusEeay6tcwHd6gRq1gUiB14qEzldfcjwxmpONQp+ipWRXgQnoPudd8gf548Luy23rEw5D2khJpZIMTEM+2zCaxTASdOzWI9nzBFHnmNBUZGc9eBTUlwWMAbd6AN+//xreWvqnVI7aVPa9FCItxoRdAFmnYAfpUsTeexE0LbtYBky+BVlJCJy4HSF32lzNMxbaQ7SBTgA0D3UNTeryIEObGeMv9xYJbu2KU1bTYSXioVxVzlsIVUJUwn7IGRRH4R4SLuOx06Ncq+9DZzeL8UhzER1PYbmQBcr630SkJbktzAZzkDHM8+QNY9MVZEDjHvcf5+UDhsmqQSRgicvnRshJs62cqTx9B7wct+I+xiIx6hP1x4K6d7XIBN6JiSq8VmmYEwpQuWkJk2T2mY9vruZdE83MtoQrQrz5METdpVR/Tonzu5oQX2N/IaN0hJ84bIJc8djAJHzF9eMzoyqAl6MeKQV3WdhqEC8rK8+pPgLKNFdjxPhfCBy8PPwvQLC5KAhl1taKn0+/ECaN2xQkYNUaw/cahD07sh+6gkScUwrkFb0+IVnRMyfft9jMXsPeIk+2D3eMdfeKt233d6zMb/ujz9l+mWXSk1TpTRrTa6fB0lPvDY+yxasujcX5+YIvhEqmzVj/C0jOjvP/KFczn/+S6mub3Js8MeVuUANgifYKI14LpMKUCb13YUr5Yhde6T/ldQrE7lJA9K1Dmbj36h1SGEI067ZmZ9CYl2j4CcRwmRNU2M4B4U9/X2sZIiJEyFkTmanZfQYgGFvXnmP1XsgVtMzt48Xj3OAvgfVz/wmB/c8Rxqa62TuqtekvO53V8+jJaFHDUpKSqSmpiYsepBuEQW3ZUi9dklOV5BWZE4pqmnWZPaQ7WVNQ6Mc+fWvxgJeal/cJx07O8PQv+CFr6QmqONyKjp2kwbty28/cQw1CER57ogYmHsnXPPGd3LgwG7pH0XwohFwonVw4kQkCF1IrGMUEvtt1PuJihCUlUVEEFKtOYglJk5WV+Rkd1qGUY5V/1AUoLjEsUA4HR7P3AchP6dA9upyhLxT/rgcfMnlKUkvSmRXZz97H+jo+gOQbn0R9DKk0SoRZXN6kl335q1KitRfJnV2turH5Ffvpai0gKqEqYIOAlEOwE2HD5AJL30TdjWS8uGOl3hCjLFWH1LYWC0ZQuJEpCMhfQgag5AGobRUjdMlrchNqtT6H//wJGSOlj6UCnG02yhAvALneB7Ptg9CkJycXCnIKZLTb3tEiru2SbqRj2iGcljqmlT5VVRYgohaP2ZjBSo+5Uh+h+KkOw/m3gc6eg8EDSLzYAUqfW7ixIkpjSTolYi++fbsYJlTTTStUTZsWBDV0EfUAI3UrLokZ1J6kt692ahBMHZvjnZbKoiWDhQSHaMykWH+19WV6W1DEFvoIBAFIgWIGOgfbr2iQF1jk/pSSMsogl2ZMycCZaNTEG31IUUhTL3ZmS4kNjY7SzZutRAQIG/78VyVbpROmoNY6Ia9/PiV/DX+ImmorpW8PW+Xpvzga+6gK3Ks9KFUdFqGwNeLsW5V5vSNSTfLiTffJW07d416LNzmhzbArg9CUUfnVZDcGPnRiNbVuW7ZBlkzbbFIQ/Pm8qyn9Hd0XL+w6n2gYzWPud9++0223nrrlDoJqGBkrLje3FwX1dA3Rg3CcZ+elA5E697sV2dnvUpSPMeIJSiGjXD/8bvIqU99HnY/pB256m2QLjRkV/VCL6RXfJSkjFB35KKAz1icnyeNzc0y+p7Z3jscJpqQRkCL1AhEizhAeAyBMrYYx0J3IpL8JaE3O3vnyHfUFmOs5C/ftFxtk4GdFsL8+Fhxr1++XG0BnAI4B0qTYGO0pBMw7FHV6Nlr5skLT62XNQVbSl5zowz8bqrkNgbqahcU5kbtiuykt4Jdp2XgZ0dnYwQA1X8gHMYWY6eEypwaCpU31NbIU5ec5/pYXrHqg6DGDlfn1ar+2hq1jWbk67d76ercWFEra57e7ByEPY6D4/qFVe8DHcwXWjh0zz33nNxxxx0qNSlVwKBvbg4vVW3uiGyMHBijBuHE7pKcThi7JUfr3hxvZ2dUSYIQeo/536stxq7P1aapGeaN9Nki0nF31dsgXfBiJ2QhdBCI+pD/vrZaBnRvI6+dt6e8fcFwyc0RqWtsjvplkHLcNFXL0MYqxmZnWMkf+eJIGfvqWNn7xb3VONE4aaqGsqaoXPTr6APUFmOruVRhdl5iGfZNeYWycMczpSk3P1C+X7cIEVaLQih9SItMH7ISR5904zC1Bbpzgi2cFT+Ip9EZqPhzue1tbo8VD3ofhK6XDlZbp6vyiBT8edN8+evOBWqLsZ2Rj3mn0Qzj140aqxc60hFwely/0HsfWDVIw/xxxx1n6UDoqUnYplKHYLywdob+ZlGzuQiofaO0dMSp0W50IpLdeE23C/S0IidNzdqW5Cs9o/5ViS2iDY57G6QDGWgnJAqmGLVg8MFHpaKrX18UJlA2C5aTJjZyi1uNQAY3VsGK/QUfXCC1TYEvqarGKjWee/zchKYdxdJCwOhWegOIkjGuqpLl4y5QP/XN1dWhOeyDtCNEFvQmaolIPzIf29yTwaoHQ4QuICdXmvJLpKa4oywccKY0Ba9vQwxBsZv0IV0cnUjRstdGZwCG/9v33mm/g4tj+QEiBvkdnWs07CIF3S7f3TJlCca/k3NA9CKUnhSMZuR3KBEpzI1wEvA4To7rJxAdQ1dgV8Xo9NNPlylTplimG+E+HTp0SJkOYbPY2N7QtxI15+a2kiGDX3XVYC2V2Bnt5vKlxv4IXnsgWFVJclIJyZxOhNQhq6ZmcAjgRMAB+GLZumCUwWhLBFKRUppe5DZVKIPtBL+hg9BCMX4BmEENYyOuOxwmEzcagRRWJYqXFZtWhJwDHYwx36ddn5RpIazKmmpVcCYkotQp9m1YsSKmwe4VszOw5d13yR8XXxLmvBgdFVvDHkZLY63kaKIcBaeCYj19KKRBCKYPRTP0EylajqfEaYRzYfwy8FAuNdmYxc16pKC5psHSyHeasqRHM8wC507/6B+pQXBxXD+BM6Ab+nASjLRv316lGuliZR3MmfdNJm46Ils5E2VlfZNynn7k8jsx2p06EV6rJEWrhGSVTgQdAZwEbGEzlBbmywX79pXht88MORFISdazDgCapX18+Shp26owfYuYWDkPGWwn+A0dhBZI6AvAwjkAmuEDDmehNB1WAfwoc5bCqkTJxO8GaLoWwuqYVmVNpaQELdpFqwmktuilTnPbtpUVRxwZ02D3EjFQxzZFMlaMn6C6N1s5KsaeDBGGfXGejDltN+nUepp8ee9vKnLgVFDspLeCkUSKluMpORpyLgxOQn5hkeTm5Up9TU3Cy6XGS8Pq4HvPYkUfkQgrIz+eaAYchy2vHZbSKkZOnQekGr3wwgthKUVwGCZPnpzS0qcw/qOJi/V+BxA1O3EmnOK0j4LXFX2zU+HEaPe68u+2SpKbcqXQF0BsjNsROYBzEHIiTJWLAGyHDTWNqXMQYhUxsXMeYtkJDS1HvJyjofYZaVGdlBESHHnnzKj7IM3okyv2VR9wTx0O05kM/IDD6N/rhb0ioghT9psiI3qMSFr35Vgr9znFxeoXT6vZfJ56pKCgRw+lSTDTZ8Z7npqohT1uq1bhzkAQNQ9HxdCTwc4hsSpPmoymZl4ew00nZq9VjD578xWZ89yTofGIE0+VXf92iC/lSxMJ0ougOQiLIIhIx3/2l5J+HVN2XukEnIPVq1fLU089FRZNgEYh1aVP3fQ7iLdJmtM+ClYdj81dja2iC3ZORSxnw8njJSrygQVEFCYxpxMZKxHFsiGs7uOmZGpc3ZV1KpYGRMZmLvw68NsP8bE5SmCsgNhgsBMA/r/mZ5FXTvNWVj0D7Vw6CB4vXCZj/gKw4sl/DpFR/dKzuVVLZc6KOXLeB+eFzcEJ0Lsrw4nYe/requKQUS/gZ/flqELgFStk2d+PDWgPgr8sOSUl0nfWTMlv00btA8GyuYmalwiC1bHC3szBY3e/+y75E2lGcaQ0uTHGveLmMZLhtMCpQKUic3oS+hukq2Ogg6pFECabgcjZjY4hG4jWObmiokLuu+++iPtceOGFKdEi2AEnYM7c3SOaqQ3YcbIs+m68rXFv5TwY54DVca3Kq0IoDEGxmU+H7qBW9K0M/iFtS6Ma+bGMdj80CF6JVdLUyolA5UMUN4EGweo+To8f67EdAwPfzgmAsW/nPJizEZYYIg2KHHunIsvsXKYYtUD0esXjnv9SfZhLCvKkWdNU/iDSih48YVc6B2kYcdi67dYRc8buyqnsvgwDP7ewMJTiEzgBTa3qN2/YINKmja9N1Ky0D0CPJOjOQdE220ifDz9Q5+BVFB2t27JfWD2GldOQrE7M68tXeRY4O8FrVCOe3gnJFgynU1dlq87JdpqDVGoRrNhcuSi8DOrCRRcGm6tFNkmzigwA49z2219neVyrPgrR0oLs9AJv77pt1DQhvXypHX71QPACDHI9nchqFV8vi64b8qXBNOTderWPufJvVzIVjwfsbos4Xqzf5WipQk51Bg2mNCWF1mLEy3QQWhB62A6dDSE2gnOgOwTD+nSMP6SXjYZ6PJ2ak1xRKNXdl630CDDYm+vr1Yo/jHO/mqhZPRacAt0ZqF+6NCBQNkQOvKQxpQq7KEEyOjEHuiffHD7poyjZS3dmN9hVG7LTBHjtqpzOmLsqW3VORoUjKzCfTk5CXn47y/nm5mpL4x7RAXOH5W++PVdycnLC5n788TrJzS0NHmdzBMGqvGq0XH5EF6wcAdicbgXCVo/rRnPgJ7ADolUttHMiYlU6tNM46CVT7W4LO+6SWaK9eJLk1G8Srai15Nj9LtsVMXGqR9xkrmhkJPvFy9nxbUhigrAdQoLIG0SnQ12gXNPQpJwFgA9g2joHThuXwIlA7qFes9g8zuB6yHpFIRj9wFxRKNbtiUaPEMBQBzkw/pubZOlBB4f1QsB+MNbjKXFqfiw9GoFUJjgPVtWL0qVhG6IA5qZoxrna6gZ59yHrhmu6qNlYix9jvzoxb+6dEP4eLygq9kWUHG9vhnh6J5ibptn1SrDaT59rrm6MuC0Rxj3Sf+LpS2DVVVkvZRqtsRrGmPcbpPbU1Pyutq7v27jecj4vr5Vl74TIXgmaasRmnsO43/bXK6fASR8FfUUfaUXY6uk+enTB8JFU486FBXLrdj2kNLjy70QgnAkY+yPoToQbuwHOBFKHrHolRLstREOtNL5womjB32WttlKNbX+X7Rqd6s4D0oqwtXIwWpt7LRnI0iInRhhBaAGYQ3pi/Jr0sceBL8IiL9UI7Fb7R14qMvtOZ6v/VtGJdb9Fr4ecgohGtIpCTm5PNHqEIEyP4KFikZNeCXbRCKv0I6vqRekSGQD6XH5Bror+NDVotlECt6VU3WBX3vTQi6+ULbfvH3eKUDy9GdxirDYEwz8UUSgKRBSKerWN7JUwbbFosAXqm0P7Af2+oWMHb3PatM2vtCCn6Ma/0UkwG/96YzXz4/ktUHYqBLZjc++DzZ9pjAfseK8s+u6iiN4Jdr0SNkcQNkcLOnceq/6cCp2tVvStogsTenWRIfMWB8a5OfJAv63k4M7tMt458EMfYJeepNsM0W4Ddev+kKKGypDNnpujSW5DZWC+cx9/KyAWmCMNrUWOfkKk47aRv/tpkorsJxQptwCRcqyKA/igR6s24ATfhEVuqxHoH+4IQZI4FxRZpRGBF08SqTc6CMFjjF8k8tN/Rd6eGLgd9xk7SaT/4VnzxRAv6FzspWIRHINN//2vrLzhRtVPwYuw2E8xtJ8gAoBOycaSpgVosJWTs7mUqhXB0qdGnUGihNMR4uTQ50dC6UDAa4pQssXPauW/okbKp3wTrkkoyJUtzt5Jyu8PRE8tCWoXtOZmEaPDZrgN0Qm/0pIQMZg0aVKEUe+2qpAuTF6zZo288sorMZ0N7L9u3bpQnwQ/HQQ7gbGVENhrFaOamuUBkWzx5kZpTjQIbh2VWOii47b5eco58KsCUbrgpLqR0+PopVLtqiRGW2z8fVWFtJ/SX0qlVjkHzVqOVEmxrDtvsWzVJUHi+gaT8W8ep1EqshMoUiYh9LCdVa1icPLQreJyDqKJjnyJJDgRFNnmCsYQFFlFJ148MfAY9UZhEiyk0kBU4p4dw2/DcV87W+TtS0SOey6tvxiShZ1GAPNRy5aOuyBM6OylV4KfYmg/sdIPNNTFTlUpKMyT0WfsGOYQJEo4be6dYPRaAulAN6uXs6GuNixFyKmBH09vBrcYowZhoGlafbOUP/y1Mu41NHey+mIMNlezJHgbtAt+VUaKlhbktKqQOQJx9NFHS8eOHS2rGOksX77cl6iFG4GxlRDYSyO1DRsWWBr9dvu76Z/gtimaHl2w0yS47WGQbkTTDjjNPrBaSLS6bzT9Q+cObeV8mSiTZZK0lhrlHIyXifJgh7aSMAoMkQazM4CIwr9Pi53hkIFkrjtLXIf0SousM8qenrdMGfnxfnHoVSbNoqO40cN8cArscv9scwUNmaG43SwoCjkWBkcCxr+KHJishlPfDaQsmR0HnfrUahTSCTuNgJ2RjlV/Y6Mzq/QgN+jpR4hYYOtXt+Z4sNIPFBTlSoFxzgQiDPuduoPMeOw7efaaeSoCgTSlRIJoAAz+Y669NfwGlQ5UIw2h6EJ4ipDb459+32Nq66dA2agVMKYPWdKAZK5AJEBPGxJEA2xei4gIQrDxml/EqwmwEib/+9//juoc2ImZrfQPXrQRerqPlVbAayM1Y8lSsxg5MK6z3N9uzgqUGEWJUpQ2xRZjp9hpEtyIk9MRsz5A75cUpg/wsJDo1vaAPXPayf+UveVRGVF3j9pinBT9ZIPFguLL/4y0IfTFyAyHDkILQa848H9j+0XchmpG8RjzjoRF8RJLUBThRLQWGX1DdKfC0rFAPkdZ5BzGeYVRKhoEScEXA/ofLN+0XG2Tcb9EGOkRugEDiEREizzY4YcY2k/0js1IFwLYHnjuIDnQMAcNQj7SjnB7cZ4ccNYA+fCp7y0Fywk918JC6b7t9iqFyGgRFBQVSX4RmuFt/rBjH7fVjXB8aA78jBwYBccrb/s0sPofqw1oQ7N0Pm+nkJi50yn9Qw6DBF8HK2JVRnKDbngDrN7rToJbTYATYbLX+yAygfQn9E3AFmMnwBAPaAPshcBGAbOdmNlq3kqMrEcn4sGubCnmnaBrEpBWFK84GY+JiITTx07GQmNR/ubn0ayJfLEskJ6WzIVE2DOfXD1Wnpt4nNr6ls4ci01WC4pVgeyCWIuRGUjWipTRIBqiJBL+AT952NYy+f2flVMA9DzCeIz5WKIj34glKLIqabb72d5qJQPzXPutTalOqS975rVzcrI6LutGesz92tqEh4uLVT+DdDHy4wWlSk++ZU9Zs7xS2nVtJU31zSqyAH2BnkIE9P8no6yp83QjpBbVST6chIJCaayv85wi5HcfBEQOwgTHVlWGYPCr/CgtTEeQ36EkZOjr1Y+QOpRbUiArb/8sQpwMh8J4H79FyePHj5e//vpLunbt6qrcqBNhspf7OCmZGg27dB+zViA3tyRYcag2LF3ITn9gJUa2K1PqBqQCxZsi5EcPA7tGaW5Tn/xkQPc2kmuwq2obmxynE4dSnU0aBq+2R6wyrAmhtU2681FPGDosZ091o6wTKS9dulTOPvtsmTVrlhQWFsrf//53uf/++6U0mOrQkkTKdkKfRAmKE1bFyC3RqgnY3WY1bzVnzD9EpGHIGSKfPxZIL0qyOMlr52S/Oi47qTQUr6gZeO2AnA6YxcTGKkY60bohWwmbzYLlRFNbWSlTz/+ncg42C4uL5fib7pZ2nbu4NvAT0QfBrnuykTYHbi2FW5ZFVDOKVolo40fLZeO7v4Udo83ePRMmSobBjYWt+vp6T3oAL1WQYt0nUR2XIwXMRgLG/p57zpFPPhlhK3KOt0KSFTDAo3VATgZ25/Bw/15yzuJlSe2ubOyfdP7zX0q1Rdre7EtHOTLWE1rMJFkssREkw15A5UOAxcQ0dhCc2rlZ5SDgC3fQoEGy3XbbyeOPP64qMxx00EGy++67y7PPBleFW4iD4KRVul/GfNo4BiBaNQG/Kg3EqmiQJJAeNPbVsRHz7xz5TtTOyV7vFyEo1kXAPhjwEZWHjKRJFaJ4S5pCaAwtQcjYN4Db7Yx+u6ZpyWL9qr/k8QvPiJiHhsBtadJEVTFCxADpRcZKReZiZnrVIeCkMZrVMf2sXGRneBsxVjHSqxPpmgLzWMduPhrR7uNXdSUzSBn6ZN6oqPvsustz8uVXKBoRzp7DZoZEznA0nAqPM2X1HmlF0D+YQU+F6iQ6LkY7Qiw+Vm6rIMJWWF4RiEb27JDGfZdi0RBj8TDNKxm1yCpG7777rvzwww/y3nvvSadOndTfDTfcICeeeKL6gkPINt3xw9h2UlXIr/BcWq0IROuXACKqFZ0gMmGxSEk79y3cjalOsVKfgqv2fvcm8No5Od6Oy2ZBsZdKQ7EqD4Vh08fAzwiG32DlXxn1Bu3Ae49+FyhnakG0tCE4A8YUpGRFDnSQBoSVfrNRj3m3qUKJ6oMQ0T0Z1Ynqm22rDpkrD1l1U8Y4rIqRz5WLrNJ7zOh6gPXr14et8o8cOVJmz55tueoPox3HduMkYB+7aECi+iVEpggZCUQKSsv6x0wj0oXHfmKVImTnNCQCXehsjCC0ys2RqiRWRwrZEcbKRaZ9WhUG0oud2CppZSvES4HhN1+PHKiy6NlVySirRMrz5s1TX5I9DUbEqFGjpLm5WT777DPJpG7H2GKcllWFfK5KkFABkS4athMX3dM/vCuz047NLkC+P1J6sGq/z/R91NgPvHZOjrfjckhQbHiDeak0ZCdq3ubttyQHKYEGIaxZqIwIBiIOSEsydmlOF0LaAcPbDc5BoO9B5P6xuiHrZU2T7RwYtQhwCoCuO/jjx8UqGoDoArZIHXLqbMQrco7VPbnbFXsEqhLlxK46ZNVNGShnweExvKAb3kZRstnoxhx0CGYNwIwZM2wrDxkFxXfeead88803cXVlBvhdRcQAaUXY+lEG1SxghgYhN7c4TMxcWNAmpsg5UehlS7H1IlyOR2BsJXSeuuPWSa2OFLIjLG7LCVYwmnfFvhFGvrHTsnHOGInANqW2gl8sCdoLU/YwVT7MjkpGWRVBWLVqlWyxxRZhc4giIKcTt1mBL1fjCg5CL6nAz14CfouBElkX2Vdi9UuwEhjDSdA9/cbaQFQBcwAf8DhXARA5gBgY+f4Aq/YTZk5wne9vF4nw2jk5no7LXnocOAWRgOI+faTnA/fb9jFIRAQjUSVNzdqBUJqRhQZBtCZZv2q1b8JdP9FLk+rRAhBKFXLZD2Hf086RDx5/SJVL9bsPgrF7clhEwabqUIS4Gd2Un/k+lEbk5BhWWEUkohne+mq/VU+CmpqaqFEGY6QBxzA6E9AyvPbaa/LOO+/E3d8gWpTBDcaUILOAOXDO4elC0UTOycKtcNmPaINVFMPcsdlrdSQnwF6AE6AXNDFHDqb+Y7C0bVXoKEpgthVASm2FRGQsWFHib6f1ZJNVDgJAtMCILrGwq2h06623yvXXXy+pxk9jO1lVhZLiiLjJ77erSKTfD/83OgA62Hfx6yJvX2x9m7m5mgtggMMp0EFKT2VDpZp3mu8fq/IQjHu3xwLm+zlNg/LaiMxNSpAeTbDaP6Ikqk0KUjqUNA1pB4o2awf0dKGi0nypq2pUzsQfPyyUh87yV7ibSLykCoWLk0vkwPMvke2GDk+YM2SsSGRnqMdKI3JyjGgN2pwIoY2Gt9lh0LUGsVKR9MpDVmVLvVQeShR2omJjipBVulAi0ojiTfkps1m9t4s2eNEK6FEMP6sjOX7sgjx54IRd5dSnPo9wDj65Yt8I5wCLnGdMWyA1wagAbA+M3xw3XHVNtsJuPiPYZNec1UDNusgU5gwiq1KMunXrJuWmNIfVq1crJ8FOf3DllVcqoYb+hxWcVOB3LwG97wGqC2CbiFy/UAO2wsCH3HdHxEu6T7R+Cfg/NAeqZrGE9z14e2Kkc+DDKoCe7488/8Cj5aix03z/WJGIWD0MnPY6cJsG5bYRmZeUILs+BnoEI1oKUjqgOwMn3ThMbXVhsZ4uVNyqIKA50JoChrNpNR75/ekCjHtjOlHFnytcpQrhuRjLpaLh2gdPPJTw89YjCnaGvZM0oljHcBKRsCy7GsNh0A15q1Sk0aNHW/ZLsGq25rQngh+N0aIRrbGZXf+DdMFNbwM92qBZRBsyjVH9OsuT/xyinAIAm+RRi8gBgPhYdw50MB59z2wZNekjy+NvqAmPKqQULEZWLA1vdNpgMRezOWuQLOiFkMHuWyR77bWX3HLLLSoHUw+lfvDBB5KXlydDh1rXeMeXqd0XajJJxKp/MuoE646I71WMogmOnUQS7Fb84c0f93x4lGHsJJHXzk7IKoCe7w9jHpEDt/n+8UQinPY68JoG5bTHgd8pQeYIRk5JSVr0SjCXNDU6A6kQ7vpFyLg3ODBv33uHHHTRpfL2vXcGIgIxUoUinqNyEqI/R6dpOvFglUbU7rA+no+XKGGzVWQB1fnMQmTdmXjxxRdVepHXrsxeyqXGYnNjMwlrbFZe/o78+NN1vpYqTQROV+/N0QadX6pq4xYTJ1MorTOsT0d54/zhcVUeqjY5DolKe/a1AhGIVpXInLFQ0CooNqvZvH8GC5SzrsxpU1OTDB48WDp27CiPPPKIWv04+uijZcyYMTJ16tSMKHOaTiVDU3ou8NoROTCDyIDHdB/bLwmACEVEuDCoYfChEkG8VYxwf6zsmysP2Rnxbnod+FH21EuPA0Qf4kkJqpw9W1aMnyBadXXSeiXYpUnFU4o0UaU/k1HmFBEDJ1WM0EvhwdOPi5g///EX1f3Mx3CbphMvcEaqF66Wda//IlLf7PkxcZw/bpynjqHwuTSq0xKm2G/x4sXy1ltvhSIA2P/444+PaegnqqxpZN+DQEUimCDNzdWWfQ6M902lBsEtH6zZICcuXBo21zrOkqSp6M/gpvLQhup62emGGVGPh0gE+iikrIqRlSPQc2jw99+kXdSChUyMc0ZbQLchkGGARUTdljDPWdkOKSqN7tbOzaoUI0QK3n77bWnfvr1aXTnssMPkyCOPVI3SMgV91T/VzoFVRSWr6gQJIyJ8F6N9ebRQYLQoA7b6SoB6PAM+dkTU8/29ljh1W3lIjzjAOTBHHBKVBhWtSzJW+UP4kBIEQ/2Piy8RraYmLCqB+URhlyZlVdIUY8zHUyUoHZyDWJWHcI6IAMQ619oqayHfki8/j6iE5EeajhfWv/GrSPAxvD5m3bINkmNccsvPtRVHo7lbrOMb03yM1YmwxdgOGPL9+/cP097pHZFjpQxZ6RicpidFSxUyVy3Cdvvtr5PmZmOZ00BUAc6AUbcAxwI9E7DF2AmpTFvqWxr5m4FV/99r6jxXNkp26pLbKoWx0oX0qkeJTHv2lJWAEqVW1Q/rK+2rEsHRuLNPYBHznh1F1i/bbEvg/5MH2KdGJ6BSYtqkGOHLJZUip1h0795dXn755VSfRkZj9cUAsVFuDv6fJO9fN9pfPDHwQYVuwM5Y96NBia5dcOL9pwg3lYfc9DrwKw3KChjRy887T7SazY5bTnEgRSielKBkC5WjpUlVbWwOq0qkKuhG6WvgpEpQspwDYx8DYPX4ugMTEhh7cGBCvRRMaUYzHn1AGhvqw7QXZ972WEL7D1ilMWmNzZaP2VhRIwVdSh0fSzkVjZuNv5zcHCnq1TZsP6fREXOaDwpw6Ma9E9ExDHpzmpFu6EerRmTVn8FJepKTrsYYoztyVeXiQI+D3CL5Me86yz4HMOxra5db6hbMEQYv5xIv0ZqmWYmai3NzZOwXP0lls+YpPciNUDoVhVNiCY6bNZFFf260tBuSkq0QISoOGv3AqvqhZhFBgE0AR+OF4wKpRAD2CcaXBbsoR0uNjid1OgW4jiBg9QJNWp566impCv5YkuzCqo8CxEZ6ubPk9jswqLatsPvAOY0kWEUVoDfQowtphtNIhNuIg+58IK0IWyutgiejetwFYc4B0HJypWTXXdXtSD/ysuqfCKEyVvw3rK6xXPmP1vtBL2lqDHZF62tg9zhOV+NdPaf6epUiZCV4NgqPHzzjeHnw9ONtexroDgzSirB1W2EJz+mgiy6zOL+6sGsKB6KmqTKh/Qes+h+UT/lGBKv8pq8ZzOt9ERzrDwwRBN2xCY0dRkf01X5jrwPzyn+sVX0rwbITQ99KFB2rMVo0AbLZcP/kkxGqOzK2GzYssOxzgHlEC+Z/OiaoWwiPMMBxsIsOOD2XeLUASPdBt2NsMTYCh+Hh/r1UczO9AzKogpXssI9CPELpVBROiRVBgA1hZTf41f/Jc1ZC+62DGQRl4ZkDxz0XOQebABEH3TnQwRjz0XoxgVi3Z3oE4fDDD5dff/1VJkyYoJqm4IvjtNNOk2EJzv0lyQEf3rrGpvDypaYOiknpd6Ab/rE6E9qtCsRRmjRbcNvrwGu5VDuUUW21iFBdLZv++1/568abAmVSPegHvJZatSOWhiBa74dcm5KmVk3N4tEquCW8rGh46VSz8LjR4KTZ9TTQHRivdOjew/oGWB4G7UXrLTpJycklnvoPOMXKUM/JzxUtP0ekYfO3nbkvQjT0ikjqmMFFR5y70bFxKmK2K1dq1UQNKUhWmoR4OiBbiaK9CJAxr5cntTPcEQ0w90LYrFUwkyO5uUXy+YKjbKMDTs4lHpyUMYXDcM7iZcohgHNwY9/uMuHHFXF3QU52mVM3hVPMZc+twHFQ7agoPy/kaPjV/ykm0cqgGzMIjJkDEy3moqH2ayUSLPgReNxWm1OjY/VqynQHoU+fPjJlyhS5++675dVXX5UnnnhChg8fLv369VOOwsknnyxdumxug04yB6MgqbggV4rz85TXX1qUL43NzVKHMHwCG6+F4dTwz7APnJ9YCZ/Nc/EY/fEKq5XxXFoa4STktGolK2+4UYmL46lqZO6VABCRcNJnwUiEhqC2Sd596Fs5bdKIkJEfyyEx9jcwVjGK+jhBrQLu53eH5NqqGnlj0s3SUFdrafRbVRVKdBWlUJqRQYydX1AouXm5quypMXUpv2+h6/4DbrA01BuapeOpA2Ttk4s8pTc5aazmxImwS/OBkQ5NAdKGcBsi+ZMnT45aacitoe+1MVpAPFxmmSrk1HDXjfdAZMBas5KbiwWp5qgpR07OJZFN08wORHVTs1zzy59SlpujHAY9PajUY3qQuTdCKqoUYiERhr6xspHZobCiKD9XjpjySUj0fONhO/rbbDWW+NfOEbCrflhgMYeIg5UTgPnAszA9qOa8V1Oa4flbt7i4WE444QR5//33VUQBHYwvvfRS6dGjhxxzzDHy9dfhIWqSWboDOAOIjs6YMFJ9QTx+yhBpFfxywDYRjdc8iZRDAmOLUKAdbgXNKcaql4FV3wK3vQyi4cexlFH9wP2SU7L5tcgpLpYu/3eVaPpKvCldx8tjQHNQ8+WXrvss6MCoVxoCw/d4Q12z/PLFale9H/SSpnbGfsTjGLQKfoIoxbTL31Odis3pOzD6LYXHRmL0NPCKlRj78MuukXOnPmeZuuS0/4BTwa+T/gdFPVuLFBoez2V6k95Yreulg9XWrC3QnQg4BWpsEx2xSvNBBSL8xiJyP378eJk9e3ZYCpKdANncVyERWAmQA+PN10033I0XHWOz4W61X25uqeyx+/9kyOBXpbm5Jqqo2cm5xIOuBTC8ddRYN/atxMRwFm7brqfSIehzTZomn2+InaYNh8OrsDkRhVOwkLjLDTNUfwP87XzDe2FpQdEKY+bm5IRFC6554zslXval/5NT8a+xQIkXCopFjn9hc1ETbDHGPBwPq/QjYwpRtF5N2VTmFCKqJ598UqZNmybr16+X4447TjVwgUj4zTfflJkzZ2Zc6lGqy5ymClQnQg6gGVQcwJeDm3JnvhFLfGwuVeokFOiHoDmJWPUy2KXzLhElTFvlw9gTR2VNY0UH3JRIdYLSGqwIhNcb//hD/phwcXhUIZiu47UvAo4Pp8Cc/uP0eFjZf2LiHOUUGCkozpPTfFzdx+M8eencQAQhuIyIdKRTE/AYdbV1UrfuYazbBm4wlE4FcBTQ7Az9DOA45OM6aQFNQKI7ORuF0fHqLeIph2p1X7B22uLNzkZBrnQ6pb/vJVad9niwK2uKtCLoAc3AeXC66p8IYpUjdSoettvPrlSqlWg5kaVRo/UjsCtH+vmw/jL4k++USFkclilNRd+DWAuJ0AmYIwSlwe7Kw2+faZlihNuvPaS/XP7Kwohj3vP3nZSjEJdtAVvAqlRpIsW/DRbRimjnAVJY2tSLnes6xaimpkb+/e9/q9Sijz76SJUT/de//qWcg9LSgMeOCAJWOFByNNMchJaKOX/QmEZkV+4sIXmCTsOBXgz9DKsgYNXA7KIPL5Kn/vZURNO0qsbwlahojdRiNVDz0pQtGjDSi/v0UYb8smOPk+ZgapEx5Sge/UC8FY1gnI88YXv54Mnvw+YbLCoRWTVDCzuXKLfnu9AqeEWPUuRIvhSUHSoNlW8qJ6GgqFgZ/X/8uDhMl4BmZx2694xaxchv4tUyxBL8Ou05oK/264Y6gGg5VhUiP9CjI17TfLxWGooHJwY35qPl+cPIN+oN7I5jt58eHdjsPNhHB2KdSzyVi4a0LbXVAuhiYt2w18XEGxqbQs6BEx2CE61DsjFXNtJBAZPvV26yvO3/xu4gJw/rpf5/41vfR9gYBw7spv7iqmLklxbRTX+CApuUJKsUouXzM2ph0rOD8OCDD8qtt96qtAYPPPCA7Ljjjpb7HXDAAco7IZlBNEESogtWeYLIQdy2S4JXM6w+hF4N/QwTNFsZ6tWN1XLKu6dIcV6xciBiRRDaFrVV6Ul6pMBJ12Q3JVLjMuSD9HzicSno0EE5EF6chGgCYqf03bWzzH7hJ2kwre4bKxHFEhj/tmiNvPfod+oYxtuNToMTrUI86BWVEKXIK9hKctufIwX5NXLK7Qeoxwo1ZAt1RL4zTIycDp2bk9m12GioqzQl4/EkceVV4yUeAbIX/Cwb6tRwt9vPqZNhh9fIgpvVfCsx8fqGRikN6hDEQZnSWFqHVC4kWkUQdujWWm31aoc6937ws3IQYome4yp44ocW0a/sgt6mhU0Qiiqk/8KkEdduKPIg//zzTyWOsnMOwNixY9W+JL2I1uxMFySZG5mYy53pQGxkVZIs4Q3VvJYKc9t8LcWYG5jp1DQFDLxWqrV7oITp5H0ny+RRk8PKmp424DQZ/fLoMB2BkwZqbkukei5NGuT3Y4/zpB3Q0QXEcArU2ENFIxjPB54zUDkFwLy6H6sZGpyDtx/4NuBgBIXOuB3zSPl59pp5agsnI5ZWIR70KIX+PIqKi2TsuBFSXFqyWZhso0vINOx0BF7Lofp9PL8wNkuzEiAjrQhbo0DZ7j5eSEbZULfozoNb58Br0zW71fxougBdTIwtnIsh8xaHnAMnZUqttA64T11zc0r1CDDqSwzfXShoMvUfg6Vtq0K56fABEffRRcfRbIy48aJFTFS5dLPOIcNKm/qmQchGslmDEI+OAPc96+kFYasDeojQmGqUFK1CPPmGCdQgxFv1xwoY9UgrQuTAzOuHvS6FeYWWVYx+rvhZLpp1UWhfPQrwv6P/J2P+PSYiOmClL0jE81GN09Abwar8qQ9aBL2iked0JZsUIfQtgJFv5qQbh6l9n7hkjjTUR/5oFxTlSUN94jQHbp4Hcv9DEQRDWVFzOdNMIh4NQjKOFy/mZmlWlYr8uE80UFUIBrWZPYfN9D2FJ5G40S+YgUgYPQ/MfDp0h5ir+VaahFZ5ufLFsP7SriDfcdSiJChwrvHYaM1PrKoYGTUKVeijZGMjJJRoKULRbkPREoibzVz4dfzZBanQR/hk56YmkY0kHbdt083AyH/9/PCQsrEkmR+PkZTVglgVBDxWOPKzgpARaANmHDMjkEIUBEY9Igs9WveIaJqG/8Ogv2LuFWHH0SMFG+o2OI4OOG3K5gZU/tn6penWN8ZRzchY0SieDs12q/vRmqHBELdyDvILczenLCWwapHT52FVRchtR+R0I1bVILfVjtwczwtuKi5ZNUuzq1QUz31i4bT6ULqzucyqfQUkr5WLomFV1aiqqVlpEmI5FluVFCpx8+wh2ws+ybVxNFrzsyISjH2kF+PPaPjrEQaURgelMXon+I5dhaJYFY5K2osUGjum+5hdEG90I5M0CCQzcds23QqsFNgJmWM+RptcfxX80QTMXnQNcUQXnOT1x0ObwjZy7773qmPCyI+V8oNV/5pGU6k1rFrlt1LOA/6eHRuoZgMnw08HwAmFPXoENANmPUIwgpDbtq2nfgaJJJrAOJT3jxKmBkb9o5989OyPEfPrVlWFCZ/tiCWI9gIqE535wBNSvmypdO61jRQj5SvDcSL4dRMZcCogTnR0wqpZmt492a5SkZf7xGLZsj/ku0UjpO+270t+foPk5JT4WjbUTKKqD8XTH8FOeOxELKw7F+aqRtGcC7Pe4dZte7gSODs5ZqIiEHa9E1JGLM2i/rtfb4hq+23E947DXkkhdBBaCNGqFDkllsjI7jG6rp0vMvUf/qf12Bn6XioTxFHhyO+qP/F2Rda1CzgHI5P2niRflX8VtYJRMjA3HQvNl5ZKp3PPkV/33c9zh+VEYicwNjsPOnAOdhmzlXz6xtKw48x47LuYaUaJ6rhs1V15y+37e65e5GfJ0kQRb7WjVJ2Dl0pFflc32hyR6CRr1hwjRUVYBGknI/baXdJdDG3GTQUkP7sYu3UurPQOV/y0XN2vyoWTkcqKSHrvhHhA5oEvTka04iSwCYy/+yo03Epk/CKRknbiK07tlTSCGgQT1CDE98E1axCmnjBA9nx1j9Tl3zmNCsSRg4gIAtKKnOT1JwukOOkRB0QO4Bzs3m13X/sbxIuuGUDEoHnDBrVVzoHHfgappra6QZ6+8pMwzUFBYV5IuGzWLthFEZz0S/ASXbDSIIR1Mjb1P4hl/Fs5G4nqnRAPSOn5684FEfNII0pWdSKv55BqDUIy+y3EoxFw+ziJ6o/gpERqLOfCTu/wwA5byZU/rfAUAYhHQ5E0Q96Ar1rGaBoAOAnRfvcbXJQ9zSAS1geBZC5+hf6irQ5EPMam31NXWtRNVCCOMmkwrs8adJbc/cXdagzjG+NUOQd2EQeUO01EpMOrOFjXDCjatFFpRfH0M0g1dVWN4c6Aho7MTVIALQJyzm3Kp9p2XJZI7QKcCq/RhVAVo9BxNdUYTa8qBccBBj9Ey+Z+CWbjH86Dut1QMlW/b7pFEvTqRGr1PvgaoINxMqsTeT0HvVKRVbM0P+9jRzL7LWzWCEiERsBPMXQi+iM4Qa9q5DUl6eAt2qk/txGMaMd0GoGwcwQSUZTE975Ldr0JMB/td3+JzeKi04yELHAsKFJuYZjbpif8MVJZWtRNebE4hESIIDzy7SNhc1O/narm4wH3h1Hv9ThmkbG5bKoudo6nvwGqEqE8aTxlSm3LoCKCUFYWtZ+B6tIMx8KUd50KIoTMQaHyfqf1ty2f6ug4BkF0rHKr0UAkAMZ+RL1iU9nTDeV/WRr/cAp0MqlkKlJ4kO8Pg1yNCwP5/8lKL4r3HPRmaW4MfS/3sTsOIhBwCoDXfgtYtUc1pGilUbNFDB0vekoSDHhgTEkylk7165hOgCOA6kQj75ypthgnqiiJrmU0fLWEFUPxhF1xErvffTDdYnHx5/eii52dCKIzCKYYtaAUo5SRwNKivpcX8+D5w4hH9SIz7xz5jueVeWO3Yz1FaESPEZ6OZZd6FK8GAUY5nAI/U4LgYOjahFgaBDf7xoObCAlW99958FtpNFSpyS/IlTFnD5D2XUodpwTZRQnsyq0eNmEX6dq7TcxjG9OCCopLRNOaA4a/oezp8TfdJdMuOTfivqff91iokVomlkyFDkDvmpxM5yDdzsGrFsFrRMKNrsCLBiFVKUPpkpKU6GOGSpeatIWvnbenjL5ndsT+6HEQj/7A7vESWirV/Ltvl3JcWCpSX21vT6RhSdN47Fw6CB4vHHFJskNu+uOt+VnkldMS6pz4rUHA8YxaAZ0p+03xxUnwq78BVu4ROTDTZ8Z7caUEOTHIE+Gc+OGEYCX/iUvnSoOpclFBUa6cNmmEq2pElr0MzPoEA07TjYzaAqtUIoiWnRj/XjUImSBsJqnTFbgx+BMpak6VEZ/Kc7K6L5qeInJgZsaEkapZaiIM+aT0U4pmm1gZ+oWtwisdWekUy38QmbJH9H3sHjOJ0EFI8IUjGRSxOPoJkY7bJvTDOP+LR2TCt/dLZW6OlDVrcs+gC2Tobmd7OpZdRAKRhNnHzU6ptsFPIz2exmaJck7snl9TTp7UF7WT4vxG6Td3lu352q3wG5uqwegvKs1XmgUvJUyN0YUwPDZiszLYnRr/bo39TBE2E39IZJM1P0TN8UQfklU21I1x/1b5elXxCCVR3Z6T3fOJtqL/xbJ1CTPkEyF+dpXdsMR021FPBBcbbaID2P/FE0XqjToaiwhCqjIqDNBBSPCFI2lKKkJ8wcesq6+U8rxc6dzULEWF3h8TK/wjXxxp2T3ZbdpSIrohG1n7+ONSfuek0LjzpROl4+mnx7xfvOlByYgg6E5IRbvtZeGAM6Upv0TyGmtkzD/6yjYjt3cdQTjgzAGqxKnRsPdawhSP89eSjfLGPV+5qpDk6jF8Xun3Oy0pU9N1WhKJrEwUr/MRT/TBqjsy8voTVTbUiXF/qqGMqbg8p1jPJ9qKflIM+VTYCg21Iut+C+zbfuvwngkw7pFudPRTItsdYHGsIGYHIE1SkNhJmSQEfBkg5OhGiOTlPkkRJsfCaVfl4GMWaZr0bGxSW8+PGRQXQ3Ngxq2gON7uzrFE0jDS1zz0cNgcxrEEw7hdOQcw7jGuqgqMTfeLJkDWeynAKVDj0tLA2Mf0IkQ2tNbtAs5B0JjBdsarf9mKgrFyfyBEyIW5YRqEkHNgKnkKZ8GpyNj8ONAc2ImZ/QBGOzQHfqUB+SlsRtOxP2+ar0qHYotxtgMdAMqOxtMVOdnovQfgFATG0XsPOBEz+yFqxvEDzkHgOwjbwNiZENaqO7LeuMxPnHQ+1nscoEeCETfnFOv56NUJoS/A1hglSEbhk6TbCrrQeMoeIo/vL7I8+NsJQx8ZCUqLUBWIKGDfiGMFOf398OiAn/ZJEmCZU5LQkmZJzSWMs1xpGG7CgH49pgFoDaA5mPjRRBVJ0AXFTqMA8XZ3Noqk7cTMSA8yd0N2UpY04n4W5UydRBgwRsTAa5pSLHC8DjdPkqZXDUZDTq6KDuglR61ANOD0u0bIxjVYKc+RNlsUR5YuNWAsYepXd+d0RK+iZI4gYD4dG5+lU4QiVm+DeETEiQar8ogYxErlcbuiH0/jM6clVe3y+eMpG+pUI+A0hUk37q1o7fCcnDwfP5qfpRV2v9sl7UUeGWFdHh38+7SgUNlw27gF1gJmRB6cPGYyqjp6gHFZErWk2W43zpCZP5R7KmmWqDJoUYmjXGnM/gl2kQQ/HtPGSYDmAGlFMOx36byL47KnendnXeRs7Hng1bkwP66TsqRWUYBY93MaYTD2UkhUI7VOo4ZKAVbpxd0qPYz0Dt3KpEP3UvV/qxKoOvGs+uvdnZFWhK0f3ZYTBSIR0BzAKQDYYuw2QgGjXTOKtDVRY8z7RTwRCjgWaIyGrR9s7mYceH7YYqxHEuA8TJo0STU0wxbjdEPvPRAtcuBlRV93PpBWhK3TFCEn0QcY6Ei7QYMxbDGOt2xotGM66XxsFUnQjXvzV0tZbo7jUqbxlkHNSOx+t2vW2a/yb7KJANy/a1DArEW3ARJkKyQKRhCyjHjzAUNGfZ1u1DfJqU99LrcfNVBFAXSMtYntVhX0esZu7uMLes1jr1UCorVmt2vuFuMxvWoB9F4GTlb0jeg9D8yVlZykKOnORayGanqaT2il35TmYxcFiHU/JxGGZKGnDMW7Sm9e7TdrEOJZ9cd9bbsyp1nFIAiSoTmI55wS3fgsnggFHAl137omdY7odVDct31c54PIgLFRGcBYjxhYOQ9olJZukYRENUkzNz5zIjyOFX2wM9CN+fxYzcfYacUgJ8fUowt1zc1hUQFjyo+52Zpu3OvRBhj3t23XQzVVc2Pgu30+WYHV7zYWAqOt8heZblOdLA3VjRBJGL9IpKSd88dMU+ggZBF+pPOYjXqd6978TkoL85TDAPRKBnBE7MBtOA9z9YNo9/ENfOi8dmr2Gga0eUy3xr0f6UKYx+PoPQ/gHDhNUXLjXNil+dhFAXQRcbT0ID3CYBYgR2uYlkj0VXpzydF4jhNPFSOnpGvFIF3bEG/TsZAh7nPjs1CEQiIjFPkdS5Ke+hStm3E05wGN0uIhmWlL+oq+WczstkmamzSlaKlP5rQdOwPdaXdkJ8c0phTByC/JzZHaZs1RCpNfxr2b55M1mH+3o3VeBsbbrEqfYowohJ2DEK99kkRagIvYMvArnQfGOxwBMzXIx23eLMApys9VDki0KAVuwz6lhQE/FNtY90kLfAwDOk3XSUS6EJwQOBF6ipJTp0R3LuAUgFjOhVWaTygKYBCk6lGAaPdLlgDZ6yp9vMb85uMk9qsXkYNY3ZAzGazKw/DueulgtY13ld4qQmEUf2McK0KRqNSnaN2MdefBiO48xEOy05bcipn9SlOyS30yp+3kuMjntyPaMc3RBV1wXJqb4zjlx2uXZRJllf+8TwNC455DrbsyT1gcWEw0vqoYp6mmwC2MIGQJfqXzwHh/4IRdVVqRmbrglxY+Cvm5ubJbr9g/ynr1A1/KoCW6uYjx+D6FAZ2m60QjnnQhPUXJLbpz4bVEqpcogLEvQqIFyKnErkuyn4QqBukYKgbFs3qfTmBVPtqKfrIjFF5Tn5ys1EOQjLQh837Yjhw5UmbMmBHaF+N4VvztNA+JTltyKmZORJqSGau0nXhz8qMdE5WKzNGFmmZNZg/ZXoqCTgQN/ySzfL59oRJjBCBatCHDoYOQJfiZzjOqX2d58p9D5Pznv5Tq+iZpVZintqHFYHHnfPhS/SDRzUXsjh9nGNBs3OsNz9oWtXV8jHjSheLBq3MBYukMzNjpFZKtOUg0tdUN8u6Ub6WhPmAMoOwpnAW3zc2SVTGopaJHKNxUMfLiWMSqTmQExrk5bQjG/OzZs8PmMN599909G/OJTFuKhVlPkIo0pUTm5Nsd066K0FaMCKQGu0IlEy36FWSQpsAtjEVlCX6n88BJ+PKa0aru8bwr9lXOh6HojBonRUsQ7cNas95ZnwKvx4/3uBbpOgBlS8f8e4yrvgRe04VSiR4FQHdjbO0aobmpWpTpkYOnr/wk5BwotM1lTtOxYlBLRo9QuNEPuEl9ilWdKF5j3iuJSltKNH6kKSUjbcfqmC2yilA647ZfQUEwopBFzgFgBCGL8DWdx7TyD2dDF0CXJltLYFdV6J7+AUFQvBEFL1WLXABj/n9H/09GvzxaahprVCTBbV+CeFf0U4WuM4hGOlUtShRohIZIQUNQwGrEz+ZmflcMIolLffJjpT6agNkruubBHNnIhKpI8aYppZIWWUUoXTCnL2dYv4JEQQchy0hUMxO/nQ9XRHxYg5iblXhtV56EL4MNdRtU5MCtFsFredRMIt2qFiUCu0ZpBYW5CW1uFm/FoGwm1Y3Q/DDuE2XM22keMoF40pRSTbZUEYq33HpSMaYXF5aJHHSXSP/Ds1pb4JQcTdMzywnYuHGjtG3bVjZs2CBt2rThRUnLD3GwzbkZVBXwuuKfYI0DDP19pu8TITSOFkGItzxqJuGkc3KmRxCevHSu0hzoPmhBYZ7849Y9pbhV5hhf2UIi+hV4wY0GIRrp3EmZtCz8KLee1MjBpL6Ri4+6DdBzaFZqC5zauXQQPF44ksIwIFqhTx4QueLvNYJgPn6Cvgxg8OtC41gGPxyKvafvrcqjOnUoMh1jFaNsqlqUzOpFxFnkAN2RzdWG4u1X4BUa9yRbQORg8E3vRxRLQfZBWkYSoGG8z6o3jE82RYbbuUwxIplDokuLJbh5iZvSoX6UR800nOgVMhm/Gq6R5DVCS0YaklV1IjvoTBA79C7MqdQv+FVuPWmE0ovNon5/dYiZCh0EkplkaGkxOAVwDmI5CfH0PiDpi94ojaQOp/0K0iUNSYfpSMQOYxfm1sEKSBA9Z3K59aQ2RX3xRJF6Yw+NlilKNkOZPMlcMrC0GNKMkDo09tWxSpNgV+rUbTdjQogz9H4FcArU2KJfASIHyjkIVp3CVo0bDCVqk4gfJVFT0SGZJB5zF+bK4Bjzro7T0CS/r61W23Qpt560xcZLfxU54pFgV2RpsaJkM9QgmKAGgSQKL7qCllDFiJBUEC19qHFtjfx154KI+6C/QSK6N8eioqJCGfVmLrzwQlfpSXAKzFWTEt0hmSQWdGHeY/73EfOfDt3BcUUkv4XFGVXFKIk6xEyzcxlBICRJ6LoCvaOyUVcQq/cBnQNCktcITU9DUi1t1c6ixuY0pGThR/OyRDRV85umpjqpqfldbYkz9C7MhreqGmPeqTEP5wBpQQBbjOONJEBzkFHOQYZmJSSSrHIQampqZOjQoRF/b7/9dqpPjSTS4/ejm3IS0HUFiBwAbDGmroC0dBrr62X9qr/UNlPSkOKNXiBK4TRlSe93oDsJXvodpHuH5IqKj2XO3N3lk3mj1BZjEpt4uzDrwmK94L1RWExaNlmVYlRZWam+7B577DHZcccdQ/N9+vSRLbbYwtExmGKUQSS4d0EicFPqlJCWwLKFX8ubd90i9TXVUljSSg695CrVBTodSEQVo3jEz/FWMfJL6Ow3iBjAKWhqQuW2gHI8L69UdUXOpG7ImVjFKONKk5K4aZF9EHQHYd68eSpy4AU6CBlCRIOTzKlbbNQVAGoMSEsFEYOHzjpJ6mtrQl20C4tL5Nypz6ou0NlGOvRgSMdSqUgrQuTAzJ7DZmZsV+RMIqOam5G4adF9EC6//HLJzc1VkYMzzzxT9thjj1SfEvEbCInCahdnTt1iXVfQkjolE2JF5boKFTkIoWlqjPl2Xbq26B4MicJN34VkUVjYRfLyyiIiCJhPdOSivn5V8PFbbqQCzgAiBhkpLCYJI6s0CGDw4MFy6qmnyhVXXCFlZWWy1157yUsvvWS7P0Kt8KaMfyQD0BucGKVZGGdI3WJEEeAcoKIRgJOA1CPME9JSKGvfQaUVqbwGgAhCSSs1n42km/g5XYBxPmjgFOUUBMalwXHirgs1D1kiLCYJI61TjBD+GDNmTNR99t57b7n99tvV/5ubm6WpqSksbDpu3Dh5/fXXZcWKFZb3v+666+T666+3fOxooReSBmSgBkFn+ablqheCmXeOfCdrOyUTkmkahESQbg3Y0olkrehT80BaMhuzQYPQ2NgoCxZE1qI2glDpdtttZ3s7nIMjjjhC1qxZIx07drSMIBhLv+HC9ezZkw5CppChdYsRKUCjNHOn5Gg9EQjJZi0C0ooQOchG7UEyxM/EOdQ8kJbMxmzQIOTn53sWG+v89ddfSo9QXGxtPKKSg7n0G8mgpid63eIMQ++UrFc0Yqdk0pKBU5CNmoNYPRhIy9I8EJJJpHUEwS1vvvmmdOvWTYYMGaLGP//8s4wePVp23nlnFUlwAqsYhcPqBomFnZIJIST5QIPw7cLzpKmpUjkL0Dx06DCcLwXJejZmQ4qRW7777js577zz5JdfflFP/tdff1V1nu+9915p395ZjicdhM2wPjIhhJBshVWMSEtkYzakGLkFzdE++ugjKS8vl9WrV8s222wjrVq1SvVpZSx6h0UdY4dFVDsghBBCMhUIodlngZAW4CDodO7cWf2ROK9jmyLVNMXcYRHzSSFDBciEEEKIH52OCUkVfJcSWyBIRkfF0sKAH4ktxkkRKqOEKTol37dzYItxqoHDUrE0sCWEkDQDXZIrKirUlqQPcyo2ycCPF8ke879XW4wJSXeySoPgB9QgpEEVIxjgcArqKkMVJqSoTGTiL6mLJGRwzwVCSPazZMkSmT59uirbjcp80N/17t071afV4kHkAE5BZVOz/msmZXm5snD4AEYSSFrbuYwgkPTrsIi0Ihji6usUaIEx5lMBHBblHFQGSw9VBsaMJBCStD4J61f9pbbE4iuqoSHkHKivqLo6NWYkIbECZ/RTwDYaSCvaFHQOALYYY56QdCYrNQgkw4HmAKv05ggC5lNByGGRSIclA3swEJJJtLROy17YtGlTWMNPgDHm0UyUpK5EKjQHrfNyIyIImCcknWEEgaQfSCNCCg+cAoAtxqlKL9IdFvXVLkGHpXXqHBZCWgiIGCjnoLZGjbHFmJGEcFq3bh3R8BNjzBN/QcQg4BxUBcdVwbF1JAGC5CcGbKOcAlAWHFOoTNIdRhBIeoL8fmgO0qGKke6whDQIKXZYCGkhVK6rUJGDEJqmxphvSZ2XY1FQUKA0B2YNAuaJv9TXr1KRg81oaox5u5KpIzq0VpoDVjEimQQdBJK+wABPlxSedHJYCGkhlLXvoNKKVAQhWGu5sLhEzZNwIEieOHGiSitC5IDOQWIoLOyi0ooCEYRA0lBeXqmajwYiBr1KklQinBAfYIoRIW4dFjoHhCSF/MJCpTmAUwCwxRjzxOIrqqBAaQ7oHCS2uRo0B3AKAuPS4JjGP8kuWObUBMucEkJIegHNAdKKEDmgc0DSAWgOkFYUiCjQOSDZZ+cyxYgQQkhaA6eAmgOSTsApsNMcEJINMMWIEEIIIYQQEoIOAiGEEEIIISQEHQRCCCGEEEJAQ61IxdLAtgVDDQIhpMXRXFcnjeXlkt+5s+SaGkwRQghpoSyZZeh51DrQ8whlzlsgjCAQQloUVfPmyc/D95JfRx+gthgTQghp4SBioJyDYCM8bDFuoZEEOgiEkBYVOVhxwYXSXIUmR6K2alxXl+pTI4QQkkrQCBWRA9UAD2iBMeZbIHQQCCEtBqQVNVdWBrryAk1TY8wTQghpwbTuFkgrkpzgRE5gjPkWCB0EQtKYuqY6Wb5pudqS+FGag7IykZzgD0BOjhpjnhBCSAumoDigOSgqC4yxxRjzLRCKlAlJU+avnC/jZ46XqoYqKSsok3tG3SNDuw1N9WllNBAk97j/vkBaUWWl5JaWqjGFyoQQQpQgeeIvgbQiRA5aqHMAcjRNj7UTNy2oCUkkiBjsPX1vqW6oFk00yZEcKS0olVnHzpKiPFbdiRdWMSKEENIS2ejQzmWKESFpSHl1uYocwDkA2FY2VKp5Ej+IGBT27MnIASGEEGIBHQRC0pDOrTqrtCJEDgC2GGOeEEIIISSR0EEgJA1BGhE0B0grAthizPQiQgghhCQaipQJSVMgSIbmAGlFiBzQOSCEEEJIMqCDQEgaA6egZ+ueqT4NQgghhLQgmGJECCGEEEIICUEHgRBCCCGEEBKCDgIhhBBCCCEkBB0EQgghhBBCSAg6CIQQQgghhJAQdBAIIYQQQgghIeggEEIIIYQQQkLQQSCEEEIIIYSEoINACEkJzXV1Ur98udoSQgghJH2gg0AISTpV8+bJz8P3kl9HH6C2GBNCCCEkPaCDQAhJKogYrLjgQmmuqgqMq6oCY0YSCCGEkLQgXzKMhQsXyksvvSSdOnWSiy66yHKf2bNnywcffCDFxcVyxBFHSL9+/ZJ+noQQaxrLy6W5snLzhKapMeYLe/bkZUtzGuvrpXJdhZS17yD5hYWpPh1CCCEtOYLQ2NgoI0aMkOOPP17effddeeaZZyz3u+666+Tggw+WDRs2yA8//CA77bST/Oc//0n6+RJCrMnv3Flyy8pEcnICEzk5aox5kt4sW/i1PHTWSfL4hWeoLcaEEEKyj4xxEHJycuSWW26RRYsWyV577WW5z6+//io33XSTPPnkkzJ58mSZNm2anHfeeXLuuedKU1NT0s+ZEBJJblGR9Lj/PsktLQ2MS0sD46IiXq40jxy8edctUl9bo8bYYox5Qggh2UXGOAh5eXkqghANRArKysrksMMOC83985//lD/++EMWLFiQhLMkhDihdNgw2fbjudJnxntqizFJb5BWVF9TrVLCFJqmxpgnhBCSXWSMg+CEH3/8UbbaaivJz98srejTp4/a/vTTT5b3qaurk40bN4b9EUISDyIG0BwwcpAZQHNQWNIqLDUMY8wTQgjJLlImUq6urpZrr7026j4DBgxQEQA3x2zTpk3YHCIKiD5UBSummLn11lvl+uuvd/wYhBDSEoEg+dBLrgqkGdVUS2FxiRpTqEwIIdlHfio1BV27do26T/v27V0dE87A+vXrw+Y2bdqk9AetW7e2vM+VV14pF198cWiMCEJPVlIhhJAIeg3cWc6d+iyrGBFCSJaTMgehpKREJk6c6Osx+/fvL08//bRKGyoKCh5RyQjssMMOlvfBfvq+hBBCooOIQbsu0Rd3CCGEZDZZpUGAOLm+vl6ee+650NzDDz8sffv2lV122SWl50YIIYQQQkgmkFGN0u666y5ZuXKlzJkzR1Um0iMQKH9aWFgoPXr0UPuMGzdOZsyYIRUVFTJv3jx56623VEoTIYQQQgghJIschI4dO4qmaapZmhGj8Q/nYL/99pNZs2ap1CH0QoildTCC4wNWMyKEEEIIIdmEbt/q9q4dOVqsPVoYK1asoEiZEEIIIYRkLcuXL1eZN3bQQTDR3Nwsf/75p6p61JLTkvRqTngDmUvHksyDr2d2wdczu+DrmV3w9cw+NmaRTYS4ACp8du/eXXJzc7MjxSgZ4GJF86haGvggZPqHgWyGr2d2wdczu+DrmV3w9cw+2mSJTdS2bduWVcWIEEIIIYQQEh90EAghhBBCCCEh6CAQS1AB6l//+hebyGUJfD2zC76e2QVfz+yCr2f2UdQCbSKKlAkhhBBCCCEhGEEghBBCCCGEhKCDQAghhBBCCAlBB4EQQgghhBASgn0QSBi1tbXy7LPPyvfff6/6QZx88snSqVMnXqUMoKmpSd5++2357LPPpFWrVrLvvvvK0KFDI/b77bff5IUXXpCKigoZPHiwHHPMMVGbpZDU88MPP8gDDzwgu+yyi5x++ulht61bt06eeeYZWbZsmWy77bbyj3/8Q73+JD358ccf5bXXXlONiv72t7/JiBEjwm6vr69X38GLFy+Wbt26qe/gzp07p+x8iT2VlZUyffp0+emnn6S0tFRGjhwp++yzT8R+eL0/+eQTVXv+2GOPVZ9Tkh72zksvvaR+M0888UQZNmyYZYO0p59+WpYuXSp9+vRR369lZWWu98lEaBWQsA8LvuDuvfdeadeunfznP/+RnXbaSXUOJOkNvqAGDBggTzzxhBQXF8uaNWtk9OjRctVVV4Xt98UXX8jAgQPVFkbkpZdeKkcffXTKzps4+1zCqMAP2bvvvht226pVq2TXXXeVl19+WX1mp06dqpxCGC4k/XjkkUeUk7dkyRLp0KGD3HDDDeo1MzoHo0aNkkmTJqlmTP/73//U5xWGB0kvVq9erb5z8ZritYTDd+ihh8oll1wStt+pp54q5513nvpehqM/aNAg+eijj1J23iTAO++8o4z5GTNmqNfwu+++i7g0a9euld12202ee+459f365JNPypAhQ2T9+vWu9slYNEKCTJ48WWvXrp22du1aNa6vr9cGDRqknXrqqbxGaU51dbX2008/hc298MILGj7iK1asCM3ttdde2hFHHBEaL1q0SO3z9ttvJ/V8iXPOPfdc9XfQQQdpRx11VNht48aN0/r166fV1taq8fr167UttthCu/nmm3mJ04xvv/1Wy8/P155++umw+aVLl4b+//DDD2ulpaXaqlWr1LixsVEbMmSIdvzxxyf9fEl0pk6dqhUVFWmVlZWhufvvv1/NNTQ0qPHcuXPV9+unn34a2ufkk0/WBg4cyMubYvB7qds6eM0effTRiH0mTpyo9e7dW/2+gk2bNmndu3fXrrnmGlf7ZCqMIJAQb7zxhhx44IFqNQQUFBTIcccdp+ZJelNSUhIRtt5xxx3V9q+//gqtdHz88cdy0kknhe2DFWi+xunJq6++KrNmzZK77rrL8na8bkgR02tzI4XhkEMO4euZhiBSoKdtGtl6663DXs8DDjgglFKUl5cnxx9/vLz55ptYzEv6ORN7unbtqtI6q6qqQnNYNcZrl5+fH3o9+/btK7vvvntoH7z+CxcuZFQoxeD3Urd17HjjjTfkqKOOUr+vAGlDhx9+eNj3q5N9MhU6CCTEzz//LNtss03YFendu7fKVccfySyQbtSxY8eQo/DLL78oI8PqNcZrT9KL33//XaUmIHSt//gYqaurU+l/fD0zA6T1DR8+XOU7X3HFFXLTTTfJvHnzHH0HwwhduXJlks+YRAOO+I033qi0XqeddpocdthhSgP2+uuvh72eeP2M6GN+56Y3mqbJr7/+GvX71ck+mQwdBBKipqZGWrduHXZFkAcLqqureaUyCOSk33ffffLQQw+p3Ff99QVWrzFf3/SisbFRrRxffPHFKmfdCr6emacT+vzzz+Xss89WnznksO+3335y6623hvbhd3BmvZ6zZ88ORW8RKUChALzGOnw9M5e6ujppbm62/L3E6wrnwMk+mQyrGJEQCI2ZhTWokALMHwCSvrz11lsqjQgOAtJPdPSqClavse4IkvQA4tQFCxaoIgHjxo1TcxDRIeUEY4hb8ZnMycnh65kh4PXCawgjsn379moOhuWECRPkwgsvVFVwon0H8zOaXtx5553y7bffqpViPcK3xx57qGo4Y8eOlZ49e/L1zGCKiopUqpjV51H/7nWyTybDCAIJgVQUlDc1gvGWW26pcptJZlRmQFUiVEE5//zzw27bfvvtlYFp9Rr3798/yWdKooHXA7oDbPv166f+YEDiRwf/hz4IfzAw+XpmBqh4A72B7hwARIcQLVqxYkXU72CkCrLUaXqB0qb4fBrT/1DNBq8n0k701xOVi4wryfrru8MOO6TgrIlTcnJy1Osb7fvVyT4ZTapV0iR9mDZtmlZSUqL9/PPPoYooW2+9tXbxxRen+tSIA959912tuLhYu++++2z3OeSQQ7QRI0aEqmx88MEHqsrGJ598wmuc5lhVMUKljC233DJUjeO3337TysrKtClTpqToLIkdH374ofp8LlmyJDR37bXXam3btlUV48D06dO1goIC7fvvvw9VRNluu+208847jxc2zbjxxhu1Dh06aOXl5aG5Bx54QFWqWrlypRovXLhQy8nJ0V5//XU1bmpq0saMGaONHDkyZedNIrGrYnTzzTdrXbp0CVUVQ0VAVHq8++67Xe2TqeTgn1Q7KSQ9QC4d6q0jr3L//feXTz/9VEUOPvjgA1Xfl6Qvf/75p6rp3KVLFzn44IPDbjvrrLNU7W29SRoa+WAVExEFRBzOOeccueOOO1J05sQpeF2hJ/n3v/8dmoN4Ff0uIGDdc889ZebMmaoqFRozIcJA0gv0HXnqqadUCgp6WKB51rRp0+SII45Qt+PnGFVukGKGakZIM0MaA15XRBFI+oDPHn4n8Z06ZswY1Xvm/fffl9tuu03Gjx8f2u/mm29WOhO85tgXhQXwm5oVK8wZDHqR3H333er/Dz/8sOoBhdcEFafQ6AxAR4DKjuhDstdee6n+FdgHPaL0ynFO9slU6CCQCObMmaPCokgtgvFBQyP9QQ4kuq9acdBBB4VVWcAPGwwQvZPyzjvvnMQzJV5BhRSkiKH7rhGUWoRhgqpHSDnae++9Mz73NZtB3jqErHDSYVBYpQ7BcYBeAZ2U4SgUFham5FxJdODQoXQ0dAjQj0CDsNVWW0Xsh9dy/vz5SkcCZ4J6ktSDRZVXXnklYh4LZ7B7jAuncOj0LsmoWmX+fnWyTyZCB4EQQgghhBASgiJlQgghhBBCSAg6CIQQQgghhJAQdBAIIYQQQgghIeggEEIIIYQQQkLQQSCEEEIIIYSEoINACCGEEEIICUEHgRBCCCGEEBKCDgIhhBBCCCEkBB0EQgghhBBCSAg6CIQQQpJKZWWl/Pbbb9LU1BQ2v3z5cqmoqOCrQQghKYYOAiGEkKRSW1sre+65p9xwww2huVdffVX69u0rv//+O18NQghJMTmapmmpPglCCCEtixkzZsjYsWNl1qxZsvXWW8ugQYPkyiuvlIkTJ6b61AghpMVDB4EQQkhKuOSSS+SVV16RXr16SUFBgXIacnJy+GoQQkiKoYNACCEkJdTX18u2224r5eXl8ssvv8iWW27JV4IQQtIAahAIIYSkhA8//FBWrlyp/j9z5ky+CoQQkiYwgkAIISTprF69WukOzjvvPGnXrp1cffXV8s033yg9AiGEkNRCB4EQQkjSOeyww2TNmjUye/ZsycvLU4LlTZs2KdEyxoQQQlIHU4wIIYQklddee01++OEHefbZZ0POwJNPPinr16+XadOm8dUghJAUwwgCIYQQQgghJAQjCIQQQgghhJAQdBAIIYQQQgghIeggEEIIIYQQQkLQQSCEEEIIIYSEoINACCGEEEIICUEHgRBCCCGEEBKCDgIhhBBCCCEkBB0EQgghhBBCSAg6CIQQQgghhJAQdBAIIYQQQgghIeggEEIIIYQQQkLQQSCEEEIIIYSEoINACCGEEEIICUEHgRBCCCGEEBKCDgIhhBBCCCEkBB0EQgghhBBCSAg6CIQQQgghhJAQdBAIIYQQQgghIeggEEIIIYQQQkLQQSCEEEIIIYSEoINACCGEEEIICUEHgRBCCCGEEBKCDgIhhBBCCCEkBB0EQgghhBBCSAg6CIQQQgghhJAQdBAIIYQQQgghIeggEEIIIYQQQkLQQSCEEEIIIYSEyN/8X0JIvDQ1NUlDQwMvJCGEkIRRUFAgeXl5vMIkYdBBIMQnKisrZcWKFaJpGq8pIYSQhJGTkyM9evSQsrIyXmWSEHI0WjOE+BI5+Pnnn6VVq1ayxRZbqC9vQgghxG9gtq1evVqqq6tl2223ZSSBJARGEAjxAaQV4UsbzkFJSQmvKSGEkISB35rffvtN/fYw1YgkAoqUCfERRg4IIYQkGv7WkETDCAIhWQSiGJ988olKd9pmm21k5MiR8vnnn8tnn30WsW/fvn3lb3/7W2gM/cTrr7+u/r/nnnvKrrvuKtlATU2NzJ07V/7880/ZZZddZNCgQfKf//xHli1bFrHvXnvtJTvvvHPY3I8//ijz5s2TAQMGyODBgyUbKC8vV8+pqqpKRowYIT179pSpU6dKfX19xL7HHHOMdOnSJTR+4403ZPny5aHxySefLG3btpVM55dffpEFCxZIUVGRjBo1Spqbm+X555+33HfcuHGh/yPVY/r06WG3Dx06NCveK1988YV899130q1bN9lnn33U98qHH34YsV/Xrl3l6KOPDptbsmSJzJkzR902evRoyc3N3vXIxsZG9VzxndK/f3/Zfffd5YMPPpDvv/8+Yl98Bw0fPjw0Xrx4ceiajhkzRqUMRZsnJFnQQSAkS9i0aZMcdNBBUlFRIUOGDJFp06YpY+fCCy+UH374Qe3z5JNPyvHHHy/FxcURqVAwFrEfnAkIrrPBQYBxc+ihh8qWW24pffr0kcmTJ8shhxyijBY8V1wzGLwnnXSS2n+HHXYIu/9TTz2l7rPbbrvJtddeKxdccIFceumlksm8+uqrcvbZZysnsHXr1nLDDTfI7bffrhyhuro6+emnn5QDsN9++4UcLCP333+/dO7cWTp06KDG2VC16/rrr5cpU6YoxwDG3tVXXy0PP/xw6HMze/Zs6dixo+y4444R98X10fcDL7/8sjzzzDOSyeAawOBfuHChcpr/+usv9f6/9dZbQ8/1pZdekn333Vc6deqk3jdGZs2apRxLfB/hvvgcvfDCC5KN4LMyduxY9X2K98cjjzwiAwcOVIszuFbQpz3++ONy1llnqf0hLDayYcMGtR8cClxL3RGwmyckaUCkTAiJj5qaGm3x4sVq6+p+9Y3asjVVahsvF198sXbSSSdpTU1Nobk5c+aE7dOlSxdt9erVUY/zr3/9S7v11lu1VFHbWKv9vvF3tY2XPffcU7vvvvvC5ozXZOnSpVqfPn1s7//FF19ozc3N6v8fffSRNnz4cC0VNNQ3auvLq9U2HioqKrROnTppn332WWgO79nPP/88NH7mmWe0U045xfYY++23n3b77bdrzz//vFZeXq6lgub6Jq1hTbXaxsvHH3+sbbnllmGfi5UrV2o///xzaHz++edrDz30UMxj/fHHH9p2220Xes8kk8bGWq26epnaxsvdd9+t7b///lp9fX1o7ssvv9QqKytD491220376quvLO8/btw4dQydQYMGaYsWLdLSBT+/d4888kjtqquusv2OweertLQ05nHwmXvhhRccz3v9zSHEKYwgEJIiPv5ljZz9zBdSWdcoZUX58sjJu8nwvp08Hw8rTVgFNYbysfqXScxfOV/GzxwvVQ1VUlZQJveMukeGdhvq6VhY2UV61cyZM8Pm3VwTRFGQnoTjvPvuu6FVwGSy/IcK+e/DC6W+tkkKi/Pkb+cMlJ79Aqv3XlJGsBKJCJMOoklu0mEOP/xwtbL58ccfq4gKrk+/fv0kWdT+sk7WPvO9aHVNklOUJx1P3kGK+7b3fDykcfz9739Xq7Q6iDB5AavHp556atLzwysqPpZvF54nTU2VkpdXJoMGTpEOHTansXj5LjnzzDNVrX1jaoxThg0bpqIwSFlbuXKlSsNatGiRZQQm07938f5BVM1Ipn3vEmIFHQRCUkBtQ5P6kaqqb1RjbDFecPX+UlzgrfkNfowLCwsd7//000/Lxo0bpXv37nLkkUdKqqlrqlPOQXVDtRrDSZgwc4LMOnaWFOUVuT4eUl/gLDmt8IGSgU888UREnvCqVauUQYyQP65XMmlsaAo4B3VNaowtxqfeuZfke3ifuH2P/Prrr8oxAkgXga7FmH+P9KSHHnpI7r33XkkGWkNzwDmoD1wPbDHufvVQySnwluPu9prAMfrqq6/U/0877TRV2lh/vyGVxkrvk0iamuqCzkFVcFylxiP2+kzyPHxuvFwTpK1B49OmTRv5xz/+ISeccIK0a9dOPvroI+V8In0tG793oflye60eeOABtd1+++2VNoOQdCV7VUOEpDHlG+vUCpbeUw1bjDHvlT322ENeeeWVsDmUwYtm/MHwtRLrpoLy6nLlFGgSuCjYVjZUqnkvwFiBQfvWW285uibIFcb1wB9WPMEff/whRx11lFoZfvvtt1UedjKpWl+vIgfBS6K2GGPeC3B8vv32W/n999/D5u3eA9Bo6NcE/zfTvn37iPzzRNK0sU5FDozXA2PMx/O5gWjdKNCura1VufNW4L2hXxO8Z4xGMsTJRkF3MqivX6UiB8aLgjHm/fwuwfWwe63x/sH1wHeKDvLyoW1BIQRoOMzi/2z43kWkCIJkN9+7+nsHkRVC0hlGEAhJAZ3bFKnwNlaw8COFjITSwnw175Wbb75ZCU+xWocfePxooyqL3YomhJlGYABC2Iz9kVqAla4zzjhDpaAkg86tOqu0It1JyJEcKS0oVfNeweo2xJLvvPOO9O7dW77++mspLS2Vxx57LGJfCHb11T2diRMnKkEuqvy89957SsSaTErbFaq0IhVBgFGTI1JYlKfmvYBqNNdcc40yZFF9CE7UjBkz1P9PP/30iP1h1BmvCVbJ4Szp4kxUPoJxnSzy2hSptCIVQQhej5zCPDXvFURGnn32WSVEhzAXzxHvl0cffdQy1QgpVvgzg/S+f/3rX5JsCgu7qLSiQAQhcFHy8krVvFcuu+wyFUE74IAD1HseUTSk6iF6gsIHZiZMmBA2Xr9+vbqmKHyASCWiClgxz8bvXRQxQJUhVI9DChUqF6FQBIofWGH+jsEixGuvvaacBkQ8EaVEKqPdPCHJghEEQlIAwtnIfS0tDPjo2GLsNb0IIFUIVXtg8Kxbt04Zd+b8e6RE2DVyg2GEHyMY0jCIzSukiQZpRNAcwCkA2GLsJb1IB6UZkQ6CH25cE6RSwfDTgYEM49gOOEwoW4hVY+xnV/YyUSCNCJoDOAUAW4y9pBcZjTlEVfTV/5tuuinMOYAht//++9umVOB9AScUjhZWhpOZb400ImgO4BSocWFAg+A1vUgHFXbuuOMO9fzwvF588UXlMOjsvffeqsytHaj6Bac82Q4kQBoRNAdwCgLj0uDY++cGzjL0KqeccopaOEAFMHyXlJWVhfY59thjVbMuK/C+wvtk7dq1oTS0bP3e3WmnnZS+Aq8/HANUdjJGFPLz81XVsGipjbhWSMVCupqejmU3T0iyyIFSOWmPRkiWgpSEpUuXqpQWNyvuyIlFeBsrWPH8SGUT0CIgrQiRg3icg2wCWgSkFSFyEI9zkC1Ai4C0IhVRiNM5yBagRUBaUSCiwM9Ntn/vev3NIcQpTDEiJIXgx2mrjgGRIwkAp6Bn6568HAbgFLTdwjry0xKBU5DfkdfDCJyCkpKtUvaaZBL83iUkNlx6IYQQQgghhISgg0AIIYQQQggJQQeBEEIIIYQQEoIOAiGEEEIIISQEHQRCCCGEEEJICFYxIiTLQP3sJUuWSK9evVQ98xUrVqg/Mx07dpRtt902bK6xsVHV20ZPBdTJzxbQA+HPP/9Uz7ewsFDVF0czJzO4ZmgmZlXXHf0UBg0apGqSZwNoxIRGVrgm6Aj7+eefW/a9QP1/Y/17gOrYeI+hP4TdNcs08NxRNhKNwNAHpKamRr755hvLfdFozgzui8/edtttpxoNZgN436M7MhrGdejQQb3exm7JOvieQa8RYwlONCUEuO/WW28d0WkY1zkb3jc6aGSGDuXoGYFeM7hOekd2I/hu3WqrzdWmjNcUn0V8Lxuv/88//xz6LickqaAPAiEkPmpqarTFixerbapobGzUxo8fr5WVlWkDBgzQ2rdvr1122WXao48+qu2xxx7qr6CgQNttt93U/6+88sqw+8+ePVvbZptttB133FFr06aNNnnyZC3TKS8v1w455BCtXbt26pp06tRJe/bZZ9V1wjXYeeedtaKiotD1ee655yyPc+GFF2odOnTQvvrqKy3TwXPAtejevbu2/fbba3369NE+//xzbdSoUeoaYIzrpF+TRYsWhd1/yZIl2qBBg7SePXtGvWaZxCuvvKJ17dpVPfdevXppI0aM0ObPnx+6Bl26dFGfDX1spLm5WTviiCO0bt26af369VP7/f7771qmc9NNN6nvgf79+6v3w6mnnqq9+uqroWtQWlqqDRw4UP3/9NNPD7vvb7/9puZ79+6tnX/++aH5pqYm7fDDD9d69OihbbHFFtrEiRO1TGfTpk3aSSedpLVu3Vp9rvA9ce+992o33HCDuga77767lpubG7pu9913X9j933jjDTWPa/zCCy+E5nGt8Z7CMdu2bau9/PLLafebQ7IbOgiE+IDnL+v6Gk1buySwjZNJkyZpw4cP19asWRP6Mb7//vvD9oGhs3r1asv7T58+XVu+fLn6/4IFC9S+qaCptlar+/13tY2Xo48+WjvrrLO0urq60I/5Y489Frp96dKlyiiMZTzefffdyrFKlYPQUFenrftrpdrGA96fMGDhNBqvweuvvx4aP/PMM9opp5xie4wDDzxQO/fcc7VUUl9fr61du1Zt4+XHH39URh0cAp25c+eqz4AOjNyHHnrI8v6//PKLcpb0cznttNOUgZhsahqbtN+qa9U2Xl588UVthx12CHN0nnrqKW39+vWhsZPPA95nRgfhzTffVMZwQ0ODtmHDBm2rrbbSvv/+ey3p+Pi9O27cOO2oo47SKisr1RjfNcb3Cj5zcKZigc+c0UF44okntIqKCvX/f//73xGOKR0EkmiYYkRIqlgyS2T6SSJ1m0SKWosc+6xI7308H+7ll1+Wm266KRSizs3NlXHjxjm+/9///vfQ/5F6stNOO0myqZo3T1ZccKE0V1ZKblmZ9Lj/PikdNszTsRoaGuTNN99UIXykFQGkypx++umOj4GUgbfeekueeOIJee655yQVLFv4tbx51y1SX1MthSWt5NBLrpJeA3f2dKwFCxZImzZt5IwzzgjNIf3DnAJiB1LQZsyYodKLvv32W5WKk+xUNDz29OnTVfoF0lSOPfZY6d27t+fj4T1yzDHHyB577BGaGz58uOP7b7nlluoavPTSS7LFFlvIokWL5Nxzz5VkMqdik5y2aKlsamqW1nm58sSAbWREh9ZxfZdcdtll6vXVOeWUU+I+z88++0wOPfRQyc/PV+/DAw44QM3169dPMvl7d+7cuVJaWqrG+K4555xz4j7NU089NeXfx6RlQ5EyIamgoTb4I1UZGGOLMeY9snbtWmWgOAU59fPnz5fFixeHzWN81VVXyWOPPSbJpLmuLuAcVFUFxlVVgXFdneec4Ly8PGWIOHUocD3wB8cAOelnn322cpwwhx/phQsXSmVl8DVLAo319QHnoLZGjbHFGPPJeI9UVFSErgl0HBs2bJDi4mLlYMBghGH+5JNPSrLAa6Q7BwBbjDHvFbfXBO8N/ZrgceGkjB07Vv7v//5PLr74YmnXrp1ss802kixqm5qVc1DZ1KzGlcEx5pN1TfCdgeuB75RobNq0KUzPgv9jLlO/d5GFgc+Im2ulv3estBxWzJ49W5566im54447PJ0jIV6hg0BIKti0MrCCJVpwQguMMe8RiATnzJkTNgexoB033HCDjB8/Xu69997QHFbCYBRjNdS4epgMGsvLVeRAtOA10TQ1xrwXEEmBc4BVcyfXBIYKrgf+3njjDWUMwyi+7rrr1Nzy5cvlzjvvlGXLlkmyqFxXoSIHxmuCMea90L9/fyUeNRtldtcEIl39msA5glAVjtPdd9+tjMH3339fbrnlFkkWOG/dOdDBOB4j0+pzE+2a4L2hXxM87nvvvScff/yxEikjerDPPvvIzTffLMliVX2DihwYvknUGPN+XhNEj6xE7ADfIbge+E6JFW1BBEgHRnKPHj0kU793Ie7HZ8rN967+3nnmmWdiHv+VV15R76XXX39d2rZt6+kcCfEKU4wISQWtuwXC22olCz9WOSJFZYF5j/zrX/+S0aNHq5VzpEvAkH300Udl3rx5lvu/9tprEYbPhRdeKFOmTFGrpPizqtaSKPI7d1ZpRSqCAIM4J0dyS0vVvFdgvB511FFy9dVXq9VuGMdffvmlZboQjF+s7BkxjgcPHqyiKsZqLYmmrH0HlVakIgjBa1JYXKLmvYAqKYcccoiMGTNGJkyYoBwoGLidO3eWyy+/PGL/UaNGRVyTE044QRk3WDWH4RJPeo9bUMkFK/ZGJwHjeCq8IEVp0qRJcuKJJ6rnhqgA3h9IE9lvv/0i9r/gggvUn9ERRaWZV199VRlxSMHaf//9JVl0KSxQaUWIHAS/SaQsL1fNewXpRXi/I00R74FVq1bJQw89JG+//baKkJh55JFHIubwvoEzgPvi/6gAhs/i7rvvrr5XUEUMlbNefPFFyeTvXaR1IqKGCBK+G77//nt555131LWywvx5wvc1IjBIhfzll1/UNRkyZIi63nC8sMXteJ/vsssuns+TELcwgkBIKigoDuS+4scJYIsx5j2CHw+sZGKl+9Zbb1U/RE8//XTYPrvttpttCUbklKPs4I033hha5UomuUVFSnMAp0CNS0sD46Iiz8c87bTTlH4AK3y4JuXl5fLggw+Gbnfzoztw4MBQnnGyyC8sVJoDOAUAW4wx75XHH39cXRek5sAAQcnFiRMnhm5HugRKNdpx//33S3Nzs1x//fVqpXTatGmSLPDehUGP1w3oGoR4yorivvjcIA8exhicHzgKRucAKUNdunSxvD8MaTgYSLW67bbbZMSIEXLJJZdIsigOag7gFICy4BjzXsF7Ao40Sr0iaganB6+70TmI9nlApAHfHx9++KH6PsL/UVYXzuTzzz+vnALc9u677ya3bHACvncPPvhgpWNBRA0LEj/99JP6jOnAyYJTZAecAlwfpHVB76Q76ojY4XpfeeWVjqIzhPhNDpTKvh+VkBYGDCWkGMCQQI62Y5D7ivA2VrDi+JHKJqA5QFqRiijE4RxkE9AcIK0IkYN4nINsAav8SO9B5CBbeg7ES20wrQiRg3icgxZBFnzvev7NIcQhTDEiJJXgx6lD8gSNmQCcgsIk6x/SHTgF7bp0TfVppA1wCpASRjYDp6BXCR1qZ28gfu8SEgsuMxBCCCGEEEJC0EEghBBCCCGEhKCDQIiPUNJDCCEk0fC3hiQaahAI8SknGjWxUaoOVWDwf0IIISQRzgF+a/A7Q5E+SRSsYkSIT6DD7ooVK7iyQwghJKHAOUCTOWNnakL8hA4CIT6C+t8owUgIIYQkCkQO8vLyeIFJwqCDQAghhBBCCAlBkTIhhBBCCCEkBB0EQgghhBBCSAg6CIQQQgghhJAQdBAIIYQQQgghIeggEEIIIYQQQkLQQSCEEEIIIYSEoINACCGEEEIIEZ3/B8SUHcKNduWYAAAAAElFTkSuQmCC",
+ "text/plain": [
+ "
"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "fig, ax = plt.subplots(figsize=(9, 2.4))\n",
+ "for name in celltypes:\n",
+ " mask = adata.obs['celltype'] == name\n",
+ " ax.scatter(adata.obsm['spatial'][mask, 0], adata.obsm['spatial'][mask, 1], s=8, label=name)\n",
+ "ax.set_xlabel('x'); ax.set_ylabel('y'); ax.set_title('Single cells coloured by cell type')\n",
+ "ax.legend(ncol=6, fontsize=7, loc='upper center', bbox_to_anchor=(0.5, -0.35))\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e99e39e0",
+ "metadata": {},
+ "source": [
+ "### Gene expression\n",
+ "\n",
+ "Half of the ligand-receptor pairs are **informative**: the expression of their ligand and receptor\n",
+ "follows a bump centred somewhere along the axis, so cell types close to each other co-express them.\n",
+ "The other pairs are **noise**, with expression unrelated to position.\n",
+ "\n",
+ "A perfect search would keep the informative pairs and drop the rest."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "a67ee92f",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-08-23T22:44:26.445368Z",
+ "iopub.status.busy": "2026-08-23T22:44:26.445283Z",
+ "iopub.status.idle": "2026-08-23T22:44:26.453451Z",
+ "shell.execute_reply": "2026-08-23T22:44:26.453028Z"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "expression : (120, 12)\n",
+ "LR pairs : 60 (12 informative, 48 noise)\n"
+ ]
+ },
+ {
+ "data": {
+ "text/html": [
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "
CT-1
\n",
+ "
CT-2
\n",
+ "
CT-3
\n",
+ "
CT-4
\n",
+ "
CT-5
\n",
+ "
CT-6
\n",
+ "
\n",
+ " \n",
+ " \n",
+ "
\n",
+ "
Linfo0
\n",
+ "
162.9
\n",
+ "
222.6
\n",
+ "
171.7
\n",
+ "
124.8
\n",
+ "
108.6
\n",
+ "
64.2
\n",
+ "
\n",
+ "
\n",
+ "
Rinfo0
\n",
+ "
147.8
\n",
+ "
190.2
\n",
+ "
163.4
\n",
+ "
143.3
\n",
+ "
81.5
\n",
+ "
64.9
\n",
+ "
\n",
+ "
\n",
+ "
Linfo1
\n",
+ "
176.8
\n",
+ "
162.4
\n",
+ "
136.3
\n",
+ "
122.0
\n",
+ "
63.0
\n",
+ "
30.5
\n",
+ "
\n",
+ "
\n",
+ "
Rinfo1
\n",
+ "
199.9
\n",
+ "
162.7
\n",
+ "
189.7
\n",
+ "
115.6
\n",
+ "
69.1
\n",
+ "
31.4
\n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
"
+ ],
+ "text/plain": [
+ " CT-1 CT-2 CT-3 CT-4 CT-5 CT-6\n",
+ "Linfo0 162.9 222.6 171.7 124.8 108.6 64.2\n",
+ "Rinfo0 147.8 190.2 163.4 143.3 81.5 64.9\n",
+ "Linfo1 176.8 162.4 136.3 122.0 63.0 30.5\n",
+ "Rinfo1 199.9 162.7 189.7 115.6 69.1 31.4"
+ ]
+ },
+ "execution_count": 4,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "n_informative, n_noise = 12, 48\n",
+ "\n",
+ "genes, profiles, pairs = [], [], []\n",
+ "\n",
+ "for k in range(n_informative):\n",
+ " centre = rng.uniform(0, 100)\n",
+ " bump = 200 * np.exp(-((positions - centre) ** 2) / (2 * 22.0 ** 2))\n",
+ " for tag in ('L', 'R'):\n",
+ " genes.append('{}info{}'.format(tag, k))\n",
+ " profiles.append(bump * rng.uniform(0.8, 1.2, n_celltypes))\n",
+ " pairs.append(('Linfo{}'.format(k), 'Rinfo{}'.format(k)))\n",
+ "\n",
+ "for k in range(n_noise):\n",
+ " for tag in ('L', 'R'):\n",
+ " genes.append('{}noise{}'.format(tag, k))\n",
+ " profiles.append(rng.uniform(0, 200, n_celltypes))\n",
+ " pairs.append(('Lnoise{}'.format(k), 'Rnoise{}'.format(k)))\n",
+ "\n",
+ "rnaseq = pd.DataFrame(np.vstack(profiles), index=genes, columns=celltypes)\n",
+ "lr_pairs = pd.DataFrame(pairs, columns=['A', 'B'])\n",
+ "\n",
+ "print('expression :', rnaseq.shape)\n",
+ "print('LR pairs :', lr_pairs.shape[0], '({} informative, {} noise)'.format(n_informative, n_noise))\n",
+ "rnaseq.iloc[:4, :6].round(1)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "0b7126df",
+ "metadata": {},
+ "source": [
+ "## 2. Reference distances from the single-cell coordinates\n",
+ "\n",
+ "`cell2cell.spatial.celltype_distances` summarizes the distance between two cell types from the\n",
+ "coordinates of their single cells. It takes an `AnnData` (reading `obsm`) or a plain dataframe, and\n",
+ "offers several ways to summarize:\n",
+ "\n",
+ "| `method` | what it measures |\n",
+ "|---|---|\n",
+ "| `'centroid'` | distance between the centroids of the two cell types |\n",
+ "| `'min'` | how close the two types get to each other |\n",
+ "| `'max'` | how far apart their most distant cells are |\n",
+ "| `'mean'` / `'median'` | average / median over all pairs of their single cells |\n",
+ "\n",
+ "Only `'centroid'` has a cost independent of how many single cells there are, so it is the one to\n",
+ "reach for on large datasets — the others evaluate all pairs of cells.\n",
+ "\n",
+ "The key holding the coordinates is configurable with `spatial_key`, since objects written by\n",
+ "different tools use different ones (`'spatial'`, `'X_spatial'`, ...)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "7a438459",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-08-23T22:44:26.454725Z",
+ "iopub.status.busy": "2026-08-23T22:44:26.454633Z",
+ "iopub.status.idle": "2026-08-23T22:44:26.459430Z",
+ "shell.execute_reply": "2026-08-23T22:44:26.458991Z"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "