Customizing the objective function¶
The genetic algorithm in cell2cell selects ligand-receptor pairs by maximizing something. By
default that something is how well the resulting cell-cell interaction scores correlate with a
reference matrix of phyisical distances. However, the GA search itself knows nothing about that, and any function of a candidate
set of pairs can be put in its place.
This notebook shows how we can customized our objective function and reference input.
Requirements. pip install cell2cell[ga]
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.spatial
import scipy.stats
from sklearn.metrics.pairwise import euclidean_distances
import cell2cell as c2c
%matplotlib inline
1. Writing your own objective function¶
The genetic algorithm does one thing: it proposes sets of ligand-receptor pairs and use them to evaluate an objective function (e.g. Spearman correlation between the aggregated CCI scores and phyisical distances). Then, it keeps a set of LR pairs that scores best.
The objective function drives that score, and therefore the selected candidates. In the analyses so far it has been the agreement between CCI scores and physical distances, but the search itself knows nothing about distances. Anything you can compute from a set of pairs can drive it.
A candidate set of LR pairs is written as a vector of 0s and 1s, where 1 means that the corresponding LR pair is in the set. Each generation the search proposes a whole population of them and expects one score back for each candidate set of LR pairs, with higher meaning better.
To costumize the objective function, cell2cell now implements an approach split in two parts:
A factory, which is called once to compute everything that does not depend on the selection of each LR pair: for every LR pair, it calculates the communication score between every pair of cells. With communication_score='expression_thresholding' those are presence/absence values — 1 where the pair's ligand is above the cutoff in one cell and its receptor is above it in the other. What the factory generates is then:
The scorer, which the search calls with the whole population. For each candidate set it aggregates the communication scores of the selected pairs into a cell-by-cell CCI matrix, and reduces that matrix to one number per pair of cells. With cci_score='bray_curtis', and counting only the selected pairs:
N(i,j) = pairs active from cell i to cell j (ligand on in i, receptor on in j) SA(i) = pairs whose ligand is on in i SB(j) = pairs whose receptor is on in j
CCI(i,j) = 2 N(i,j) / (SA(i) + SB(j)) the fraction of the two cells' active machinery that they share
fitness = correlation( CCI , reference distances ) – The objective function
THE PLAIN VERSION — one function, everything inside it (regular GA search performed by cell2cell)
candidate set (0/1 over the candidate pairs) │ ▼ ┌─────────────────────────────────────────────────┐ │ threshold the expression matrix │ identical for │ communication score of each pair, cell by cell │ every candidate │ ·············································· │ │ aggregate the selected pairs into a CCI matrix │ depends on the │ correlate the CCI matrix with the reference │ candidate └─────────────────────────────────────────────────┘ │ ▼ fitness repeated in full for every candidate: 200 individuals × 200 generations ≈ 40,000 times per runTHE SPLIT VERSION — the constant part taken out of the loop
pool of candidate LR pairs │ ▼ called ONCE per run ┌──────────────────────────────────────┐ │ factory │ │ threshold the expression matrix │ selection-independent, │ communication score of each pair, │ so computed a single time │ for every pair of cells │ └──────────────────┬───────────────────┘ │ returns a scorer holding those values ▼ ┌──────────────────────────────────────┐ │ scorer │◀──── masks │ add up the selected pairs → CCI │ │ correlate with the reference │────▶ fitness └──────────────────────────────────────┘ one value per candidate the whole population at once, one row per candidate LR pair The next run starts from the pairs this run kept, so the pool is smaller and the factory is called again to rebuild against it.
Why bother splitting it. A run scores tens of thousands of candidate sets, and the communication scores are the same in all of them — only which of them get added up changes. Computing them once turns scoring a candidate into "look up and add" rather than "rebuild the interaction space", which is the difference between seconds and hours. The reason it is a factory rather than a one-off precomputation is that the pool changes: each successive run starts from the pairs the previous one kept, so whatever was precomputed is indexed against a list that no longer applies and has to be rebuilt.
Three requirements, each of which fails quietly rather than raising: higher must mean better, since the search maximizes; the score must be deterministic for a given candidate set, or random_state guarantees nothing and the agreement between executions measures noise rather than reproducibility; and it must be finite for every candidate set, including the empty one.
CorrelationObjective is exactly this shape — built with the expression and reference matrices, called by the search with the pool, returning a scorer bound to it. A custom objective only has to follow the same division: whatever you can compute from the pool alone goes in the factory, whatever needs to know which pairs are on goes in the scorer.
2. Some data to work with¶
rng = np.random.default_rng(0)
n_celltypes, n_informative, n_noise = 12, 12, 48
positions = np.linspace(0, 100, n_celltypes)
celltypes = ['CT-{}'.format(i + 1) for i in range(n_celltypes)]
reference = pd.DataFrame(np.abs(positions[:, None] - positions[None, :]),
index=celltypes, columns=celltypes)
genes, profiles, pairs = [], [], []
for k in range(n_informative):
bump = 200 * np.exp(-((positions - rng.uniform(0, 100)) ** 2) / (2 * 22.0 ** 2))
for tag in ('L', 'R'):
genes.append('{}info{}'.format(tag, k))
profiles.append(bump * rng.uniform(0.8, 1.2, n_celltypes))
pairs.append(('Linfo{}'.format(k), 'Rinfo{}'.format(k)))
for k in range(n_noise):
for tag in ('L', 'R'):
genes.append('{}noise{}'.format(tag, k))
profiles.append(rng.uniform(0, 200, n_celltypes))
pairs.append(('Lnoise{}'.format(k), 'Rnoise{}'.format(k)))
rnaseq = pd.DataFrame(np.vstack(profiles), index=genes, columns=celltypes)
lr_pairs = pd.DataFrame(pairs, columns=['A', 'B'])
analysis_setup = {'communication_score': 'expression_thresholding',
'cci_score': 'bray_curtis', 'cci_type': 'undirected'}
cutoff_setup = {'type': 'constant_value', 'parameter': 50}
# Search settings that do not depend on the objective. `cutoff_setup` and
# `analysis_setup` are deliberately NOT in here: passing them alongside an
# `objective` is rejected, since the objective already carries them.
search = dict(ppi_data=lr_pairs, population_size=50, generations=30,
runs=1, random_state=888)
data = dict(rnaseq_data=rnaseq, reference_distances=reference,
cutoff_setup=cutoff_setup, analysis_setup=analysis_setup)
informative = lambda selected: selected['A'].str.startswith('Linfo').mean()
print('{} pairs, {} of them informative'.format(len(lr_pairs), n_informative))
60 pairs, 12 of them informative
3. The default, made explicit¶
CorrelationObjective is what the search builds when no objective is passed. Constructing it by
hand changes nothing — useful as the starting point for the examples below.
default = c2c.analysis.CorrelationObjective(rnaseq_data=rnaseq, reference_distances=reference,
cutoff_setup=cutoff_setup, analysis_setup=analysis_setup)
implicit = c2c.analysis.optimize_lr_pairs(**data, **search)
explicit = c2c.analysis.optimize_lr_pairs(objective=default, **search)
print('identical selections:', implicit['run1']['ppi_data'] == explicit['run1']['ppi_data'])
print('objective {:.4f} | {} pairs | {:.0%} informative'
.format(explicit['best_obj_fn'], len(explicit['best_ppi_data']),
informative(explicit['best_ppi_data'])))
identical selections: True objective 0.9603 | 30 pairs | 37% informative
Passing both an objective and the data it would be built from is rejected rather than silently resolved one way:
try:
c2c.analysis.optimize_lr_pairs(objective=default, **data, **search)
except ValueError as error:
print('ValueError:', error)
ValueError: Pass either `objective` or the data it would be built from (`rnaseq_data`, `reference_distances`, `cutoff_setup`, `analysis_setup`), not both -- otherwise it is ambiguous which one is in effect.
4. Example 1 — wrapping the default¶
The cheapest kind of customization: keep the default's scoring and add a term. Here, a penalty on the number of pairs, to prefer smaller explanations of the same quality.
Note the shape — the factory builds the inner objective once per run, and the returned closure does only arithmetic.
class SparsityPenalised:
"""Default objective, minus a penalty on the fraction of pairs kept."""
def __init__(self, inner, penalty=1.0):
self.inner, self.penalty = inner, penalty # Inner is another objective function to be added to the penalty
def __call__(self, pool):
bound = self.inner(pool) # built once per run
def objective(masks):
masks = np.atleast_2d(np.asarray(masks, dtype=float))
return bound(masks) - self.penalty * masks.sum(axis=1) / masks.shape[1]
return objective
for penalty in (0.0, 0.5, 1.0):
result = c2c.analysis.optimize_lr_pairs(
objective=SparsityPenalised(default, penalty), **search)
kept = result['best_ppi_data']
print('penalty {:>3}: {:>2} pairs | {:.0%} informative'
.format(penalty, len(kept), informative(kept)))
penalty 0.0: 30 pairs | 37% informative penalty 0.5: 15 pairs | 73% informative penalty 1.0: 11 pairs | 55% informative
Heavier penalties buy smaller sets, and on this data a purer one — the noise pairs are the first to go. That is data-dependent, not a general rule.
5. Example 2 — changing only the comparison¶
CorrelationObjective takes correlation and signed, so switching the statistic needs no new
class. correlation also accepts a callable, which is the escape hatch for anything scipy does not
provide.
def top_decile_agreement(distance_vector, reference_vector):
"""Fraction of the closest 10% of cell-type pairs that the candidate list also ranks closest (higher CCI score)."""
k = max(1, len(reference_vector) // 10)
closest_reference = set(np.argsort(reference_vector)[:k])
closest_candidate = set(np.argsort(distance_vector)[:k])
return len(closest_reference & closest_candidate) / k
variants = {
'spearman (default)': dict(correlation='spearman'),
'pearson': dict(correlation='pearson'),
'signed spearman': dict(correlation='spearman', signed=True),
'top-decile agreement': dict(correlation=top_decile_agreement),
}
for label, options in variants.items():
objective = c2c.analysis.CorrelationObjective(
rnaseq_data=rnaseq, reference_distances=reference,
cutoff_setup=cutoff_setup, analysis_setup=analysis_setup, **options) # Options here includes the correlation function, personalized in the top-decile agreement
result = c2c.analysis.optimize_lr_pairs(objective=objective, **search)
kept = result['best_ppi_data']
print('{:<22} fitness {:.4f} | {:>2} pairs | {:.0%} informative'
.format(label, result['best_obj_fn'], len(kept), informative(kept)))
spearman (default) fitness 0.9603 | 30 pairs | 37% informative pearson fitness 0.9656 | 31 pairs | 35% informative signed spearman fitness 0.9603 | 30 pairs | 37% informative top-decile agreement fitness 1.0000 | 35 pairs | 23% informative
signed=False — the absolute correlation — is the default deliberately. It is not known in advance
whether the search should favour pairs acting between cells that are close or pairs marking cells
that exclude each other, and both are real: the C. elegans results put semaphorin-plexin pairs,
classically repulsive guidance cues, among the most reproducible. signed=True restricts the search
to positive associations, and should be a decision rather than a default.
6. Example 3 — an objective that is not a correlation at all¶
Nothing requires the fitness to involve a reference matrix. The search only needs a number per candidate. Here the objective rewards pairs whose communication is specific to a few cell-type pairs rather than spread evenly — computed straight from the scorer, with no reference anywhere.
from cell2cell.core import PreparedCCIScorer
from cell2cell.core.interaction_space import InteractionSpace
from cell2cell.preprocessing import bidirectional_index, bidirectional_ppi_for_cci
class SpecificityObjective:
"""Prefers candidate sets whose CCI scores concentrate on few cell-type pairs."""
def __init__(self, rnaseq_data, cutoff_setup, analysis_setup):
self.rnaseq_data, self.cutoff_setup = rnaseq_data, cutoff_setup
self.analysis_setup = analysis_setup
def __call__(self, pool):
# Everything expensive happens here, once per run
space = InteractionSpace(rnaseq_data=self.rnaseq_data,
ppi_data=bidirectional_ppi_for_cci(pool, verbose=False),
gene_cutoffs=self.cutoff_setup, verbose=False,
**self.analysis_setup)
scorer = PreparedCCIScorer(space, cci_score=self.analysis_setup['cci_score'])
source = bidirectional_index(pool, verbose=False)
n = scorer.n_cells
off_diagonal = ~np.eye(n, dtype=bool)
def objective(masks):
masks = np.atleast_2d(np.asarray(masks, dtype=float))
scores = scorer.score_batch(masks[:, source])[:, off_diagonal]
# Gini-like concentration: high when a few cell-type pairs dominate
ordered = np.sort(scores, axis=1)
total = ordered.sum(axis=1)
weights = np.arange(1, ordered.shape[1] + 1)
gini = (2 * (ordered * weights).sum(axis=1) / (ordered.shape[1] * total)
- (ordered.shape[1] + 1) / ordered.shape[1])
return np.nan_to_num(gini)
return objective
result = c2c.analysis.optimize_lr_pairs(
objective=SpecificityObjective(rnaseq, cutoff_setup, analysis_setup), **search)
print('most specific set: fitness {:.4f} | {} pairs | {:.0%} informative'
.format(result['best_obj_fn'], len(result['best_ppi_data']),
informative(result['best_ppi_data'])))
most specific set: fitness 0.2746 | 9 pairs | 78% informative
It selects a different set from the distance-correlation objective, as it should — it is answering a different question. The point is that the search needed no modification to ask it.
7. Example 4 — combining datasets with a custom rule¶
CombinedObjective takes combine as a name or a callable over the (datasets, n) array of
per-dataset fitness, plus sd_penalty and per-dataset weights.
def make_donor(seed):
rng_d = np.random.default_rng(seed)
noisy = rnaseq * rng_d.uniform(0.7, 1.3, rnaseq.shape)
return c2c.analysis.CorrelationObjective(
rnaseq_data=noisy, reference_distances=reference,
cutoff_setup=cutoff_setup, analysis_setup=analysis_setup)
donors = [make_donor(s) for s in (11, 22, 33)]
rules = {
'mean (default)': dict(combine='mean'),
'mean - 1.0 x sd': dict(combine='mean', sd_penalty=1.0),
'worst donor': dict(combine='min'),
'weighted mean': dict(combine='mean', weights=[3.0, 1.0, 1.0]),
'harmonic mean': dict(combine=lambda values, weights=None:
values.shape[0] / np.sum(1.0 / np.clip(values, 1e-9, None), axis=0)),
}
for label, options in rules.items():
objective = c2c.analysis.CombinedObjective(donors, **options)
result = c2c.analysis.optimize_lr_pairs(objective=objective, **search)
mask = np.asarray(result[result['best_run']]['ppi_data'], dtype=float)
components = c2c.analysis.CombinedObjective(donors)(result['pool']) \
.evaluate_components(mask[None, :])[:, 0]
print('{:<16} per-donor {} mean {:.3f} sd {:.3f} | {:.0%} informative'
.format(label, np.round(components, 3), components.mean(), components.std(),
informative(result['best_ppi_data'])))
mean (default) per-donor [0.93 0.968 0.956] mean 0.951 sd 0.016 | 42% informative mean - 1.0 x sd per-donor [0.94 0.962 0.961] mean 0.954 sd 0.010 | 43% informative worst donor per-donor [0.94 0.951 0.945] mean 0.945 sd 0.005 | 37% informative weighted mean per-donor [0.952 0.948 0.958] mean 0.953 sd 0.004 | 44% informative harmonic mean per-donor [0.93 0.948 0.966] mean 0.948 sd 0.015 | 34% informative
8. The four ways to get this wrong¶
Each of these fails silently — the search still returns a plausible-looking answer.
Higher must be better¶
The search maximizes. An objective returning an error or a distance selects the worst set, and nothing will say so. Negate it.
It must be deterministic given the mask¶
If the objective samples, two identical candidates get different fitness, random_state stops
meaning anything, and the consensus across executions measures the objective's noise instead of
agreement between runs.
Precompute in the factory, not per call¶
The single most expensive mistake. Below, the same objective written both ways.
class RebuildsEveryTime:
"""The wrong shape: the interaction space is rebuilt for every candidate."""
def __init__(self, rnaseq_data, reference_distances, cutoff_setup, analysis_setup):
self.args = (rnaseq_data, reference_distances, cutoff_setup, analysis_setup)
def __call__(self, pool):
rnaseq_data, reference_distances, cutoff_setup, analysis_setup = self.args
cells = list(reference_distances.columns)
reference_vector = scipy.spatial.distance.squareform(
reference_distances.loc[cells, cells].values, checks=False)
def objective(masks):
out = []
for mask in np.atleast_2d(np.asarray(masks, dtype=float)):
kept = pool.loc[mask.astype(bool)]
space = InteractionSpace(
rnaseq_data=rnaseq_data[cells],
ppi_data=bidirectional_ppi_for_cci(kept, verbose=False),
gene_cutoffs=cutoff_setup, verbose=False, **analysis_setup)
space.compute_pairwise_cci_scores(verbose=False)
vector = scipy.spatial.distance.squareform(
space.distance_matrix.loc[cells, cells].values, checks=False)
out.append(c2c.analysis.correlation_fitness(vector, reference_vector))
return np.asarray(out)
return objective
pool = implicit['pool'] # the frame the masks are indexed against
population = np.random.default_rng(0).integers(0, 2, size=(20, len(pool))).astype(float)
good = default(pool)
slow = RebuildsEveryTime(rnaseq, reference, cutoff_setup, analysis_setup)(pool)
start = time.perf_counter(); good(population); t_good = time.perf_counter() - start
start = time.perf_counter(); slow(population); t_slow = time.perf_counter() - start
print('precomputed in the factory : {:.4f} s for 20 candidates'.format(t_good))
print('rebuilt per candidate : {:.4f} s -> {:.0f}x slower'.format(t_slow, t_slow / t_good))
print('\na 200 x 200 search is 40,000 candidates:')
print(' precomputed : {:>6.1f} s'.format(40000 * t_good / 20))
print(' rebuilt : {:>6.1f} h'.format(40000 * t_slow / 20 / 3600))
precomputed in the factory : 0.0027 s for 20 candidates rebuilt per candidate : 0.2350 s -> 88x slower a 200 x 200 search is 40,000 candidates: precomputed : 5.3 s rebuilt : 0.1 h
It must return a finite number for every mask¶
Including the all-zero mask, which the search will try. NaN fitness makes the genetic algorithm
behave unpredictably; the built-in objective floors non-finite correlations at 0.
Summary¶
| To change | Do this |
|---|---|
| The statistic | CorrelationObjective(correlation=..., signed=...) |
| Add a penalty | Wrap the default, as in Example 1 — wrapping the default |
| Optimize something else entirely | Write a factory, as in Example 3 — an objective that is not a correlation at all |
| Use several datasets | CombinedObjective(factories, combine=..., sd_penalty=...) |
| Tune the search rather than the objective | population_size, generations, runs, mutation_probability, keep_elitism |
The search itself never needs modifying for any of these.