Tutorial

GaugeFixer provides a simple interface for fixing the gauge of parameters of linear models that describe sequence-function relationships. These models are implemented as subclasses of HierarchicalModel.

Defining a linear sequence-function model

import numpy as np
import pandas as pd

from gaugefixer import AllOrderModel, PairwiseModel
from gaugefixer.utils import get_all_seqs

To create the most general HierarchicalModel object you must specify the configuration of the sequence space describing the sequence→function relationship:

• L: sequence length (integer).

• alphabet: sequence alphabet. Options:

  • a named biological alphabet: 'dna', 'rna', or 'protein';

  • an explicit list of alleles: ['A','C','G','T'];

  • position-specific alphabets: e.g. [['A','C'], ['A','C','G']].

• generating_orbits: set of generating orbits (tuples of position indices). Orbits determine which subsequences share parameters. More specific model classes (e.g., AllOrderModel, KorderModel, KadjacentModel, PairwiseModel) provide pre-defined generating orbits for common choices (all orders, up to order k, adjacent-site interactions, pairwise, etc.).

Example: define a 3-nucleotide sequence space with interactions across all possible position subsets (use the AllOrder model):

model = AllOrderModel(L=3, alphabet_name='dna')
model

AllOrderModel(L=3,alphabet_name=dna,n_features=125,n_orbits=8)

See some model properties

Printing the model already gives a short summary of the model configuration, but we can access other model properties. For instance, we can see the set of generating orbits defining a given hierarchical model

model.get_generating_orbits()

[(0, 1, 2)]

We can also access the complete set of orbits and see what is the pattern of interactions across positions in this linear model. In this case, the orbits show that there are interactions of all orders across all possible sets of positions.

model.get_orbits()

[(), (0,), (1,), (2,), (0, 1), (0, 2), (1, 2), (0, 1, 2)]

Then, we can extract the sequence features that are associated to each of the orbits. Features are defined by a tuple containing a tuple of the subset of positions and the subsequence at that position, e.g., ((0, 2), 'AG'). As there are many features, lets just show a few of them.

model.get_features()[:15]

[((), ''),
 ((0,), 'A'),
 ((0,), 'C'),
 ((0,), 'G'),
 ((0,), 'T'),
 ((1,), 'A'),
 ((1,), 'C'),
 ((1,), 'G'),
 ((1,), 'T'),
 ((2,), 'A'),
 ((2,), 'C'),
 ((2,), 'G'),
 ((2,), 'T'),
 ((0, 1), 'AA'),
 ((0, 1), 'AC')]

Defining model parameters

The next step in model definition is to define the values of the model parameters \(\theta\). For simplicity, lets just create some random parameter values which have to be provided as a pd.Series with the parameter values indexed by the model features.

model.set_random_params()

Once the model parameters are defined, we can evaluate it at any sequence of interest by calling the model directly

model(['ACG'])

array([-2.00383214])

In the case of an AllOrderModel, we can also define the model parameters from the function values across all possible sequences. Lets initialize the model with some random function values f and evaluate the function at some sequences:

f_values = np.random.normal(size=model.size)
seqs = get_all_seqs(model.alphabet_list)
f = pd.Series(f_values, index=seqs)
model.set_landscape(f)
model(['ACG']), f.loc['ACG']

(array([0.32196665]), np.float64(0.32196664793539864))

Fixing the Gauge

Now that the model is completely specified, we may want to interpret the parameter values. However, there are subspaces of parameter values that encode the same sequence-to-function model. Fixing the gauge means choosing one among the many set of parameter values that encode the same model by specifying certain properties of the parameter values e.g. in the hierarchical gauge the parameters associated to low order subsequences explain as much variance as possible.

We can do this operation easily by using the method fix_gauge. The gauge defined either by their specific names, e.g., wild-type, zero-sum, or hierarchical, or by directly defining the \(\lambda\) value as lda and the site- and allele-specific probabilities pi_lc, which define a specific gauge in the \(\lambda-\pi\) family of linear gauges.

For instance, lets express our parameters in the commonly used zero-sum gauge.

theta_fixed = model.get_fixed_params(gauge='zero-sum')
theta_fixed

((), )             -0.039335
((0,), A)          -0.106147
((0,), C)          -0.469653
((0,), G)           0.174317
((0,), T)           0.401483
                      ...
((0, 1, 2), TGT)    0.524674
((0, 1, 2), TTA)   -0.636448
((0, 1, 2), TTC)    0.091942
((0, 1, 2), TTG)    1.196552
((0, 1, 2), TTT)   -0.652047
Length: 125, dtype: float64

We can verify that, while the new parameter values are different, the two models encode the same function

f1 = model(seqs)
theta = model.theta.copy()
model.set_params(theta_fixed)
f2 = model(seqs)
assert not np.allclose(theta, theta_fixed)
assert np.allclose(f1, f2)

Now the parameters can be more easily interpreted. For example, the parameter for the empty subsequence represents the average across all posible sequences

theta_fixed.loc[[((), '')]]

((), )   -0.039335
dtype: float64

whereas the parameter associated to nucleotide A at position 3 represents the average effect of placing an A at that position across all possible sequence.

theta_fixed.loc[[((2,), 'A')]]

((2,), A)    0.055764
dtype: float64

But maybe we are interested in the effect of placing an A at position 2 only in sequences with A or G at position one. We can do this easily using the hierarchical Gauge with a under a site-independent probability distribution pi_lc

pi_lc = [np.array([0.5, 0, 0.5, 0]), np.array([0.25, 0.25, 0.25, 0.25]), np.array([0.25, 0.25, 0.25, 0.25])]
theta_fixed = model.get_fixed_params(gauge='hierarchical', pi_lc=pi_lc)
theta_fixed.loc[[((2,), "A")]]

((2,), A)    0.026694
dtype: float64

Hierarchical models

In the previous case, we have fixed the Gauge in a model in which we can enumerate all possible sequences in the space and have parameter values of all possible orders, which is possible only for sequences of limited length. When having longer sequences, we often fit models that only contain additive, pairwise or, more generally, interactions up to order k, or restrict interactions within a window of adjacent sites in the sequence. These models are specific cases of HierarchicalModel with predefined sets of orbits and are implemented as specific classes in GaugeFixer, such as KorderModel or KadjacentModel, which are defined by the highest order interaction k.

Let’s try with a pairwise model for a 20 aminoacid sequence

model = PairwiseModel(L=20, alphabet_name='protein')

As before, we can generate random parameter values

model.set_random_params()

and express this model in the zero-sum Gauge

theta_fixed = model.get_fixed_params(gauge='zero-sum')
theta_fixed

((), )           -0.940365
((0,), A)        -1.810552
((0,), C)         0.478341
((0,), D)         1.185330
((0,), E)         3.573782
                    ...
((18, 19), YS)   -0.610906
((18, 19), YT)    0.617389
((18, 19), YV)   -0.678849
((18, 19), YW)    0.504154
((18, 19), YY)    0.920631
Length: 76401, dtype: float64

GaugeFixer’s new algorithm enables efficient Gauge fixing even in models with hundreds of thousands to few millions parameters. Let’s see how long it takes to fix the Gauge of a pairwise model on 100-aminoacid sequences with nearly 2 million parameters.

model = PairwiseModel(L=100, alphabet_name="protein")
model.set_random_params()

%time theta_fixed = model.get_fixed_params(gauge='zero-sum')

CPU times: user 2.11 s, sys: 23.3 ms, total: 2.13 s
Wall time: 2.13 s