Difference between revisions of "Team:DTU-Denmark/Model"

Background knowledge

To model the homology of promoter sequences between species, we needed a model that could take into account not only point mutations but also the possibility of insertions and deletions. For this purpose, we chose to use a hidden Markov model (HMM) which, as the name implies, is related to Markov processes, which can be thought of as a process which correlates the current state with previous states although it is not known what happened in the previous states. This “lack of memory” is called the Markov property. An HMM consists of a set of symbols on top of an underlying, but unobservable, Markov process. Each state in the process contains probabilities for emitting each symbol when the process enters that state. The adjective “hidden” refers to the fact that the sequence of states cannot be observed directly, only the sequence of emitted symbols.

HMMER structure

To tackle our original problem of modeling homology, while keeping it simple enough to handle computationally, we decided to use a HMM with the same structure as the one used in the HMMER software, illustrated. The HMMER structure can be thought of as an extension to a position-weight matrix, essentially considering our entire promoter sequence as a giant motif with structure on top.

Fig. 1: Hidden markov model

An overview of a HMM can be seen in figure 1. The structure contains 3 types of nodes: match nodes M_i, which represent bases in the motif, and have emission probabilities for each possible base; delete nodes D_i, which represent deletions, and don’t emit anything; and insert nodes I_i, that represent insertions, and have emission probabilities for each base.
From each match node M_i, there are only 3 non-zero transition probabilities: transitioning to the next match node M_i+1; to the insertion node I_i (which represents the beginning of an insertion); or to the delete state D_i+1 (representing the next base being deleted).
From each delete node, D_i, there are 2 non-zero transition probabilities: transitioning to M_i+1 (representing the next base being present), or to D_i+1 (representing the next base being deleted as well).
From each insert node I_i, there are 2 non-zero transition probabilities as well: transitioning to M_i+1 (representing the insertion being complete}; or to I_i (representing the insertion continuing).

In this way, we can model deletions of multiple bases always occurring concurrently, but insertions are modeled as being of geometrically distributed length, and each base in an insertion is drawn from the same probability distribution.

Our model

With this model structure in hand, we could determine the parameters needed to construct a specialized model for each set of promoter sequences stemming from the orthologs of a specific gene. For this, after constructing an initial multiple sequence alignment, we employed an expectation-maximization approach, alternating between building an HMM-based on alignment and aligning the sequences to the HMM.

In an attempt to minimize our promoter lengths, we decided to take advantage of the conceptualization previously mentioned of the model structure as being an extension to a position-weight matrix, by simply extracting the PWM, thus ridding the sequences of insertions.

Fig. 2: Sequence logo.

To determine the model parameters, we used a set of scripts we wrote to transform the genomic data of the entire set of accessible Aspergillus spp. genomes, into a format we could use as an input for the HMMER software suite, and then further transformed the results into the minimized form described above.
Using these parameters, we could then set up a mixed-integer linear program whose feasible region would represent the set of synthesizable sequences complying with the hard domestication constraints described on the design page. These would also comply with an objective function pointing towards the sequence of lowest surprisal (a measure of how unlikely a sequence is, given the parameters). Thus, the solution would be the minimally perturbed consensus sequence, that still lives up to our domestication requirements and is possible to de novo synthesize.

The process of representing our promoter sequence generation in the form of a linear program led to an intriguing insight. When representing the promoter sequences as position-weight matrices, we are considering each base as being drawn independently from their own categorical distribution over the bases A, C, G, and T, with the consensus sequence being the mode of the joint distribution. An alternative way of conceptualizing the consensus sequence arises from considering the individual distributions as being Gibbs distributions, with the probabilities p_x being proportional to exp(-e_x/(kT)) for the probability p_x of a given base, with e_x being a measure of some form of energy associated with that base.
We do not need to consider physical interpretation of these energies, but merely observe that as T goes to 0, the probability of sampling the consensus sequence approaches 1. We can then vary T to continuously move from the single consensus sequence to promoters with much higher diversity.

We hoped that this would allow us to sample promoters of varying strengths from a single model, thus expanding a model of the promoter for a single ortholog group into an entire library of promoters with similar expression characteristics.

For any model, validation is crucial. As Aspergillus is remarkably diverse, the possibility for validation was especially critical. [1] Thanks to the framework of sequence conservation and its relationship to information theory, we were able to procure a remarkably simple measure for “goodness of conservation”, namely the entropy of our models. Furthermore, we could use the score of the promoters of our specific target organisms against our models to give an indication of whether this model would be able to generate promoters that would work in A. niger.
However, despite this great promise, we soon realized that we had taken a LEAP into the unknown. The evolutionary conservation of promoter sequences has simply not been studied sufficiently to allow us to determine the biological meaning of a given level of entropy.
In fact, even from a purely mathematical point of view, the distribution of conservation of promoters is simply unknown
Not to let that stop us, we placed our trust in our models and reasoning, picked out genes with interesting expression levels, and using our software and IDT’s gBlocks, we translated our design intentions into physical DNA. We did not expect all computed promoters to work in vivo, but we had no easy way of knowing without trying. So in the spirit of innovation and experimentation, we decided to try anyway
This paid off handsomely: Not only did our homology modeling approach result in several novel synthetic promoters, which are compatible with multiple cloning standards and keep their dependencies on the growth phase. Our noise injection approach was validated, as we showed that we were able to generate “high-temperature” samples with promoter activity and to specifically weaken a strong promoter, gdpA_1, producing gdpA_2 with similar promoter dynamics.