From Statistical Physics to Data-Driven modelling in Biology - session 9

Um so this will be the last uh session and uh I wanted to also take the opportunity to advertise our sponsor Jstat. I guess you know Jstat. Uh but yeah it's it's a very important journal for our community. I am an editor so I might be in a conflict of interest but I think it's an important journal for our community and it's owned by the community somehow. So yeah you are all welcome to submit your nice papers to these journals and uh I think we are all committed to do the best to provide a good uh you know editorial service to the to the community. This said, let's go back to what we were discussing uh yesterday. So I was I would like to conclude this series of lectures by giving you an application to proteins. uh and what I was saying yesterday is that the goal in in uh in uh protein uh uh science let's say is to understand given a sequence of a protein what is the function that the protein performs and uh structure is an intermediate step that is also very important. So for going to se from sequence to structure we we have now alpha fold and similar machines but going from sequence to function or from structure to function is complicated because function is not a very well- definfined concept. It depends a lot on the environment on the context and so on and also because we have few data. Okay. So I I made yesterday an overview of the kind of data we have. I will not repeat everything but just very quickly. We have data from natural evolution from which we can um obtain information at large scale on the on the space of possible sequences that are associated to a given function. We have very local scale information coming from deep mutational scans where you take one sequence and you do all the possible single mutations. We have experiments where people do evolution in vitro. So they start from one sequence and they let it evolve and diversify. And uh these are typically confined to short time scales. So they can do 20 maybe 100 generations hopefully in the future maybe a bit more. Um and we have experiments where you can take two sequences that are um that two natural sequences for example or maybe one sequence and one ancestor that you reconstruct in some way. And then you can study parts of single mutations that go from one sequence to the other. Or if the two sequences are not too far away, you can study all the combinatorial space of intermediates. So in this experiment, they have two proteins that differ by 13 mutations. And so you can do all the two to the 13 intermediates. In this experiment, they could only reconstruct one path, but they studied the evolution along that path. And there are many data of this type that are appearing now. so that you would like to study. Okay. So now what can we get from this data? The first example I want to show is how you can infer some information on the structure from the sequence. Now today there is alpha fold that is doing that. But what kind of what is the idea behind alpha fold? This is something that comes from uh previous studies and in particular um so the idea is okay I have I have an alignment of sequences that are all doing the same uh thing and so they have also in particular they have the same structure. So this is an example. This is the okay. So these these um proteins are organized in families. Families are typically constructed by aligning sequences with standard alignment techniques that I'm not going to discuss but so you can find them. And so this is an example. This is this domain. This is protein called trips inhibitor. And it's something that you find in many different species from insects to mammals to bacteria and so on. So you have all these proteins coming from different species and you can align them. And what this is this is what you see. This is this is an example. It's a 53 amino acid long alignment. You see that you have gaps but you see that you have some specific patterns that you can recognize that are also used for aligning these sequences. For example, you see that the this one it's very often a C almost always a C a C stain and also the second one is very often a C. Okay. So you have for example in this family you have these two specific C and you have other patterns that you can exploit to align the sequences. So this is typically called conservation which means on a certain site along the sequence there is very often one amino acid or two amino acids or a few amino acids that are present. So to to to study conservation what you can do is to compute the frequencies of amino acids at each site. A simple way to do so is shown here and it's what we will use in the code. I will show you a code that does this. So you can do what is called the one or 10 coding of the sequences. So you represent your symbols as strings of bits and as you see here I'm sorry that maybe people in the back cannot really see very well but the screen is not very big but okay I can describe what is there so for example the first sequence has a gap at the beginning so we will encode the gap by one 0 0 0 okay so this will be and then if you have an a for example the Third sequence starts with an A then it will be 0 1 0 0 0. So if you put the one in the first position it will be a gap. If you put the one in the second position it will be an A. If you put the one in the third position like here it will be a C. There is no B. B is not an amino acid. And so on and so forth and so you you have 21 in principal strings that have all zeros and a one in the position that corresponds to the letter. In reality here, because of what I'm going to show you next, we are going to use 20 bits instead of 21. Why? Because, as I told you yesterday, there is a redundancy, of course, because there must be one of the 21 symbols. So, if if the first 20 are not are not there, then the last one must be there. So, you can use 20 zeros and no one to encode the last uh the last symbol. Okay. And in if you don't do that in your data matrix, if you use 21 symbols, then the last column is going to be one minus the sum of the first 20. And so you have a linear dependency between the columns of your matrix that will give you problems when you will invert the matrix because we are going to invert the the correlation matrix to get information. Okay, so I I'll come back to that. But anyways, so we do a one not encoding but with the specification that we only use 20 bits to for the one not encoding and the 21 is basically dropped. And if you have all zeros then it means that there would be a one in the 21 position. So you have you have 53 amino acids. So you will get 53 * 20 bits and this will be the lines of your data in one not encoding. So if you want to count the frequency of for example a C in the second position, what you have to do is to count the number of ones that you have in the third um column of the second position. This is the first position, second position, third position. So you count how many ones you have and you get you divide by the number of lines and you get the frequency of a C which is basically here you see it's basically one. Okay. So what you have to do is the operation that is done here. You just sum the the the the the columns of the matrix. You divide by the number of sequences and you get a vector of probabilities that gives you the conservation profile. So this is um well this is what I called actually FIA in my lecture yesterday. Here it's called PIA but it's the frequency that you get from data. Okay, so this will tell you information on which amino acids are present on which sites. But you also have uh another kind of information which is co-evolution. So you can have pairs of sites that are correlated. So for example, I don't know if there is a P here, maybe it's more frequent that there is a Y here or something like this. So what you can do you can count you can compute the joint probability of having amino acid A on site I and amino acid B on site J and uh and this is what we call P I A JB. So what you have to do is to count how many times you get the one here and the one here and you average over the lines and you get this and what we are interested in is the correlation. So the correlation is going to be the difference between P I A JB minus the product of the P because if the columns were independent then the joint probability would be the product of probabilities and the correlation would be zero. So if the correlation is non zero it means they are not independent. Okay. Now from from this kind of so this is the basic uh information that you can get. Then of course from this data you can get other quantities other statistics but we we will stick to these two for the moment. >> Excuse me. >> Yes. >> In the first step how did you group the proteins? You pick the proteins for example in the different regions that have similar function. >> So yeah this is um this is the problem of aligning the sequences to construct these families which I didn't discuss in this lecture. So the idea is that you start from a seed. So you start from a few sequences that you know to be trips in inhibitors for example and this you can know because you did some experiment or because you look at the sequence and you are an expert biologist and you recognize some pattern. So this is done this is an art of making the seed. Okay. So you start from a seed. You can do this also by aligning pairwise. you take pairs of sequences and you align them. So at some point by some process that depends on the protein you have I don't know 100 sequences or 50 sequences or even 10 you have a few sequences that you can align then you learn a model from this seed and the model is typically what is called a hidden mark of model which is a model that basically uses essentially the conservation So it's a model that looks at the conservation profile of the sites, but it also allows to insert gaps or to insert amino acids. So it's a model where you scan your sequence and you look for conservation patterns allowing for a bit of flexibility because you know that some amino acids can be removed or inserted by evolution. So basically from the seed you learn this model and then once you have the model you scan all the proteins that have been sequenced by by humanity and you look for hits. What are called hits are strings I mean sequences that the model says these are good they have same similar properties to the seed. Then the model will give you a score and you put a threshold and you say if the score is above my threshold I say okay this is a heat and I put it in my alignment and then you get your alignment. This is how it is done and it's a non-trivial I mean there is a whole uh domain of bioinformatics that is devoted to how to do this in a proper way and uh there is research also on this but this is another this would be another series of lectures so I'm not so let's say this is given basically today if you go to interpro which is a database you have this families and someone did the job and you can download this already this is already done. >> So we consider this as a given uh piece of information. >> Another question. Sure. You said that we have we have 20 different amino acids. >> Yes. >> And we can have also u we can not have amino acids. So it's in total with 21 different codes. Yes. >> But you said we just need 20 because the the last one is would be dependent to >> Yes. When I do the one not encoding, when I do the one not encoding, for example, if you have two symbols and you want to do a one not 10 encoding, it's enough to have one bit because you put one if it is the first and zero it's not the first then it's the second. So if you have four symbols you do you do three you put three bits and you say one for the first symbol in the I mean you can say okay 1 0 0 is the first symbol 0 1 0 is the second and this is the third and the fourth one if I if I have 0 0 0 then it means there is a one the the reason I mean you could do this this would be the full one of the encoding but then this line this column is a linear combin combination of of these columns because this column here is one minus the sum of the first three columns, right? Because well in in uh uh yes it's one minus the sum of the first three columns because if if there is a one it's zero and if there is no one then it's one. So this gives you a matrix of data that has some linear dependency between the columns and so it's typically not invertible and since we are going to invert the matrix this is going to create problems. Okay. Okay. So what what was the first idea that people had? This kind dates back to the to the '9s. The idea was okay we can look once we have these probabilities P I A and P I A JB we can look at correlations between sides and the idea is that maybe if two sides are strongly co-evolving. So if they they have strong correlations between the amino acid on one side and the amino acid on the other side this could be a hint that this pair of sides is in contact. So if you have your protein that is at the beginning the protein is a linear chain of amino acids. So you will have all the amino acids that are organized like this. Then the protein is going to fold. So maybe it will fold and maybe these two sides I don't know this one and this one will be in contact and then okay my drawing is not very good but whatever the one two three four. So this one and this one become a contact. And so if there is a contact, it means that these amino acids are somehow interacting in some way. And so if you if you change one of the two amino acids, it's possible that you will have to change also the other one to keep the interaction uh strong enough to maintain the contact. Okay. So the idea was okay let's look at these patterns of correlation and let's try to see if this gives information on the contacts. So okay so what you want to do is to fix the the I and J indices that are the the ones that give you the site but then you want to average over the amino acids that you can have on that pair. So a way of doing that is using mutual information that we defined at the beginning. So this is in this case the mutual information between site I and site J. So you fix u i and j and then you sum over all possible pairs of amino acids that you can put on these two sides and you compare the logarithm of p i a jb with the product of the probabilities. If the sides are independent this quantity is zero and otherwise it is a positive quantity that uh that will be um the more you have correlation and the bigger it will be the mutual information. Okay. So now what you do is you take all the pairs of sides. For each pair you have one number which is the mutual information of that pair and then you sort the pairs from the one that has the biggest mutual information to the to the smallest. So what is shown here in this case for this protein for the trips inhibitor uh the structure is known. So this is the structure it is shown here and what is shown um in this plot is the following thing. So you take your um pairs. So for each pair for example one two will have some mutual information. So there will be m12 then you have one three and you have one m13 blah blah blah and then you will have 52 53 that have some m. Now you you rank according to M. So you take the pair that has the biggest M, the second biggest and so on and you take the first K pairs, the one that has the the top motor information. So the Xaxis here is the number of pairs okay that that have the biggest uh m. So for example when it is 10 here you are taking the 10 biggest mutual informations pairs and you what is in the yaxis is the fraction of those pairs that correspond to contacts. Okay, what is a contact? There is a convention because what is exactly a contact is it's not well defined. But typically one says that two amino acids are in contact if the they the they the distance is smaller than eight angstrom the distance between the carbon atoms in the center. Anyways, what is it? I mean if they are close enough you say that they are a contact. So from the structure you you can get the structure from experiments and from the structure you can get uh the position of the atoms and then you can decide that some uh sides are in contact. So among the 10 the top 10 mutual information pairs all of them are contacts. Then you increase the number of pairs. So you take the top 20 here. This is log scale. So you take the top 20 and you find that 80% are contacts. Then you take the top 100 and you find that 50% are contacts and and you can keep going. So you see at the beginning you all the top mutual information pairs are are indeed contacts. This is what you expect and it's it's so this is an indication that the hypothesis we made is correct. So if two sides are in contact they are likely to be strongly correlated. I mean the amino acids that you put on those sides are likely to be strongly correlated. So the two curves here are the the full line is I mean the dash line is all contacts. The full line is you remove you take only contacts that are a distance bigger than four along the chain. Why? Because of course I mean amino acids that are close that are neighbors on the chain are are in contact. So it's a bit trivial right these are in contact. So when you try to predict contacts you want to predict non-trivial contacts that are contacts between sides that are distant along the chain. So you typically these are also again a bit arbitrary but you typically restrict to things that are bigger than that are a distance bigger than four along the chain. You see it doesn't change much in this case but okay. Um now the um here it's the same thing but with a correction that in this case it's not really useful but uh that that sometimes is used. It's called average product correction where you construct a score that is given by the mutual information minus here this notation means m i dot means that you sum over j mi j and here you sum over i mig j and here you sum over both and the idea of this score is of this correction is that basically one you know that in three dimensions one sphere can have 13 neighbors I think and uh so you cannot have more than a certain number of contacts just by geometry. So uh basically if the if a site I has too many candidate partners so if a site I has strong MIJ with many J there is something wrong. So you want to lower the score. So this is a way you can think a little bit and this is a way to remove to downgrade a little bit the the sides that have a lot of candidate connections. Okay, just to take into account of this constraint but it's a an empirical recipe that is uh well motivated in this way but okay anyways in this particular example it doesn't change much. So you can see that the two curves are very are very similar. What you see here to get an idea, a more visual idea is the structure. And here you have the top 30 mutual information pairs. So you stay here at 30. When you are at 30, you know that 60% of the top are contacts. So here the color code is the non-intuitive one. So red is good. It's a contact and green is not a contact. So you see that among the top 30 there are some contacts but there are also some a lot of false positives. So can you do better than that? Yes. How do you do better? By doing what we discussed in uh the the previous uh lectures you should take into account the fact that correlation is not causation. So two sides can be strong. There could be I don't know maybe this site can be in contact with some other site here and so this one and this one might be correlated but just because they are both in contact with the third one. So you want to go from correlations to couplings. So this is where inferring a graphical model is useful. And so we have seen that one way to do that is to assume that the variables are gausian and then invert the matrix or assuming that the variables are categorical which is correct in this case and then do field and in both cases what you get is that the matrix of couplings uh so the graphical model is given by the inverse of the matrix of correlations. So in this specific case what you're what you have to do the correlation matrix is a matrix that is of size 53 * 20 * 53 * 20. So you consider your one not encoding here as one data point. So it's 53 * 20 and so this this is considered as a single uh index right and then you do uh so then your if you call y I a m sorry y i a m is going to be m is the index of the row and y a is the index of the of the columns. Then this defines this matrix Y. This is a matrix of size M times L u I mean times 53 * 20 where this is the size this is L the length of our uh strings of bits and then the correlation matrix is just uh 1 / m y transpose y as we did for in in the lectures. I mean you have to subtract the the mean first and divide by the you have to standardize the the the y. So you have to divide by the the mean and uh sorry you have to remove the mean and divide by the variance. Okay but the variance is just equal to the mean because these are binary 01 variables. Okay. So uh once you do that you can invert the matrix. So you uh you can invert the matrix and you get J is minus C2 minus one and um so J is a is a matrix of size L by L. So 53 * 20 * 53 * 20. So it's going to be a matrix J I A JB. And again we are not interested in the in the coupling for each amino acid pair but we're interested in a coupling between pairs of sides. So in order to average over the amino acids what you do you define this this um uh object that's called the forbinous norm. So you just sum the squares of the J because the J can have positive or negative signs. So you can have if the J is positive it means that these two amino acids are likely to co co co co co co co co co co co co co co co co co co co co co co co co co co co co co co co co co co co app appear together to cocore together if if it's negative they are likely to not be present together but both senials are good so you want to average the square >> uh it's just a matter of convention here I I mean the minus can be or not be there I mean when I put the minus sign is that I'm considering that the model is a Bzman distribution and J is so if my probability is the exponential of plus J sig J then I have um sigma I a okay if I call if I call sigma i a variable that is one when amino acid a is present on site i. If I put the plus here then my Hamiltonian is minus this and this means that positive coupling means ferommonetic interaction so they are likely to be together negative coupling is so this is my usual convention but sometimes you can put the minus and then you have a plus anyways since for this particular application since I'm going to average the square of the couplings it doesn't really matter but usually we use the convention of physics So positive coupling is ferommenatic. You want to be together and negative is antifromagnetic. Okay. So now for each pair of sites you have another score that is not the mutual information but the um average of the square the sum of the squares of the couplings uh of that uh that pair of sides and you can do the same thing. So you rank the pairs by um this score and you fix a number here uh for example 10. So you take the top 10 sides that have the biggest couplings and you ask is this a contact and yes it is. And you see that okay this is for all contacts and this is for contacts that are far away sites that are far away. And well if you just use f you don't do much much better than mutual information in this case but if you correct uh with this APC so you do fig j minus the sum over j of fi j sum over i of fi j divid you see that you get a big improvement. So now here out of the first top 40 contacts sorry the top 40 side pairs according to this score now are all contacts and then it decays uh but you see then here the top 30 from norm pairs are all contacts. Okay so you can have a big improvement. >> Yes. >> Can you repeat what? >> Yes. So I look I fix a number that is the x-axis for example 20. I take the top the I take the 20 pairs I J that have the biggest uh f and this is my 20 predictions for contacts. Then I take the structure and I look are they contacts. So are they at distance smaller than 8 anstrom? And if yes, then the y-axis is the fraction of sides that are true contacts. So it's the fraction of true positives. Okay. So when this until up until the point where this is one, all your predictions are correct. So here the first 40 something predictions are correct. And then like here if you do 100 predictions you are tw you have 20% of uh false uh positives. Okay. Uh okay. So you see that you have um a big improvement. There are some technicalities that we can see later in the in the code. You you you have to use episode count. It's sort of regularization. So episod count means that you okay I will show you and you have to because you have this invariance of the model under uh this u uh over parameterization thing this gauge transformation as we call it you have to choose a it's better to choose a gauge that minimizes this f so you do a transformation to what is called zero sum gauge that is a a a choice of these js such that sum over a of J I A J B is equal to zero. Okay, so you fix you I told you that you can add a constant to each line and column of the matrix. So you add the constant such that you set to zero the average of J over A or B in such a way that your matrix is as centered as possible and the the part that you remove goes in the fields because it's not a true um co-evolution senior. So these are technical things. Now before I show you the code so that you can do it yourself I just want to mention that okay this was done in the 2009 the first paper I think is from Martin B and co-workers Tiwa and uh others um and uh there was no alpha fold at the time but so it was it was a promising uh result uh And alpha fold is is doing many things but one of the first thing that it does when you give to alpha fold an input sequence it will look for an MSA. So it will look it will do a nomology search to construct an alignment and uh from the MSA it will compute C the matrix C and it will invert it to get J and then this J is fed into the rest of the of the layers. Okay. So alpha fold is using this information and indeed when you give to alpha fold the sequence that has no homologous ones. So if you give a sequence where there is no family the performance is um reduced. Okay. So this this information is really important this co-evolution information is really important for alpha fold to predict the structure because basically I mean if you know like in this case if you know that many sites are in contact that then this gives you a lot of constraints on the structure. So you it helps a lot. So okay this is the message. So the the statistical properties of these um natural sequences contain information about the structure and I want to show you how you can do it yourself. So if you're interested, you can go to this notebook that you have in the folder. And so here you have so you can input a fast file that is for example here this PF 14. So let's see if it works now. Uh anyways while it uh so this is the kind of uh yes okay so this is how it looks like when you read the faster so it's the file I was showing before here you see that you have um uh how many sequences you okay I will not run all of it But um you have 800 8,871 sequences of length 53. So you you can uh okay in in the fast format you have a header that that that is this one and it typically tells you something like where does it come from. So this is a domain. Uh so here you see it's it's from 9 990 to 10 uh to to 142. So it's a domain in a longer uh protein. Um so you can drop the header. This is how your alignment will look like. You you convert into one not and this is what you get. And then the only thing you have to do is to compute the frequencies and the so the P is the sum you count the ones in each column and the the P I A JB is this sorry this is actually this is P because you have to remove the mean. So then you you to get C you do just this operation you take P you take this and you subtract the outer product of P. So okay you can think a little bit about the indices and you will uh find it's correct. Now I was saying you have to put episode count. Uh this is a regularization. Uh so what you do is basically you you have all the lines of your matrix and what what you do is to add to to say that PI A is replaced by some alpha P I A plus 1 minus alpha times a uniform distribution 1 / Q where Q is 21 and then you you can tune alpha, okay, to to so alpha is a kind of regularizer. You add a little bit of uh uniformity because the idea is that maybe because you have few data, maybe some probabilities can be zero just because you never see an amino acid, but it's because you don't have enough data. So you add a little bit of uh uniformity to your distribution. Um so it's a sort of L2 regularization where you try to avoid big um big u parameters in your model. Here actually this has been found empirically that uh you actually need quite big regularization. So this alpha is chosen to be one alpha which is a bit strange and this is because the mean field approximation is not a very good one. So you need a strong regularization. Actually if you do the likelihood maximization you don't need this big regularization and you can get uh similar results actually better results. Uh anyway, so this this does the psel count. You have to do this for the P and for the P2 for the the the joint. Um you have to do a regularization a bit more carefully, but you can look here and you can it's it's also described in the in the book. So I will not go into the details. And then here you have all the routines that do the mutual information calculation. And then you just invert the matrix here. uh you just invert and you compute this probinous norm and you and you can get these other files these other figures. So if you want you can play with this you can say you can look at what changes if you change for example the regularization you can take other families you can do what you want. Um so this is one uh one notebook that you have in the in the folder here now. Okay. So is is it uh okay up to this point? Do you have any questions on this part? So this is a first hint that you can do something with this data. You can use this data to get information on structure. Okay. Yes. two proteins are touching and then they >> say that two sides because they're touching because >> yeah the idea is that basically if the two sides are touching then it means that there is a physical or chemical interaction between the two amino acids. So if you mutate one of the two for example you you can change the the charge. So if you have a positive amino acid and you mutate into a negative one negatively charged from positive then on the other side maybe you want to compensate because if you have two positively charged amino acids they will repel. So if you >> But if not >> exactly but the point is that we are using our data are natural proteins that are functional. >> But we can also have natural proteins that are nonunctional. >> No no no no. The idea is that the the idea is that the proteins that I'm putting in the alignment are all proteins that are known to be well not functional but they are known to fold into this structure. So at least they fold and then they are also functional. But now we are just looking at the structure. So all these sequences fold into this structure the same structure. So if they fold it means that the contact is there. We have to check this like we have to check. >> Uh no in because we are we are assuming here that these are good sequences that nature already selected for us. So we don't need to check. I mean okay you could suspect that some of these are uh like false hits but in in general in these databases okay in these databases you have many sequences uh you can have or not a structure in this case you have a structure you you don't have this here you have 80 8,000 9,000 sequences you don't have the structure experimentally for all of them but you have a few that that are measured and they have the same structure and the others you assume that they have the same structure. So the assumption is all these sequences if you could test them they would have the same structure. So the if if these two sides are in contact in one of them they are in contact in all of them. And so if there is a mutation on one of the two you should expect a compensatory mutation. Did I get I'm not sure I got your question right but >> uh no it was because last year we had a project on this one function and function. >> Ah well then it depends on what you want to do. Yes you can have okay I know what you're talking about but this is another problem. So okay this is function now we are doing structure >> we did the same but >> but the yes but it was a different uh okay I I will maybe I will come to that later for the moment the only thing we are doing we are not talking about function we are just talking about structure and we are assuming that all the natural sequences are working well are working well in the sense they are folding they are folding in the species that uses them. So if if for example if this one comes from an insect I'm assuming that in that insect it is folding. Now we will come to what you are saying. So yeah that's the that's the second point. The second point is now we we have our couplings and so we can give a probability to a a sequence. Okay. So we as we said um yesterday we say I can say okay the probability of of a string of uh amino acids I can can be written as a exponential of sum over i a of delta sorry h i a delta i sorry delta sigma i a Plus times this divided by partition function and I can so if I learn from my data the couplings and the fields then I can assign a probability to sequence of the same length. So what what can I do with this probability? I can try to check if this model knows something about whether these sequences are functional or not. And so this is generative modeling. The idea is in generative modeling, I mean you know it I I start from a training set of examples that I think are good. For example, I have many images of human faces. So I I'd say okay this is my this is what I want. So these are good images. I learn a model typically by maximizing the likelihood of my training data. So now I have a probability distribution and then I can generate one new sample that should be an image of a human face. Um how do I evaluate the quality of the generation? This is very tricky because um I mean in images you can just look at the image and do you think this is a human face? Well, yes. Um but in proteins is not totally trivial to check what you're doing. And also you would like to be if you have a good generative model you would like the samples that you generate to be different from the training sample because of course a trivial way of generating new data is just to generate data that are identical to the training set but they're not new. They are you generate human faces but they are not new. So there is you have to check whether your data are um look like human faces but are also different enough from the ones that you have in the training set and uh this is tricky. So in the case of proteins uh how do we do that? I I already said that we we start from our MSA in for example the PF14 the trips inhibitor. We learn the model. Uh the model has this form. We can learn the model in many ways as we discussed yesterday. But and so for the contact prediction we can do mean field or we can do the likelihood maximization. If you want to generate things mean field is not good. Uh and likelihood maximization is not good. So you have to do Bzman learning but okay you do Bzman learning and you can check that if you do Bman learning you can generate sequences that look statistically very similar to the ones in the training set. I will show you measures for this and also that are distinct enough from the ones in the training set. We we are going to use these um Boltsman machines because they are simple. So I mean now this was just to you know typically people will say okay but why don't you use language models or more deep nets or they they are better no I mean they are not necessarily better I mean they can be better for some things but they're not always better so in this case they are not for what at least for what I've seen up to this point they are not really better than a simple bulsa machine. But the Balsson machine has the advantage that is um more interpretable, easier to train, so you save energy and it's very simple to implement. So you can have a code that everyone can can can use. So I don't see for for the moment the need to invest uh uh much more energy or money to to to have different models, but maybe I'm wrong and uh it will be the case in the future. So for the moment we will stick to BSON machines. And so this is um these are data taken from this paper uh from a few years ago from uh Martin B and also Remy Simon Ramaran Ganatan and many other collaborators Russ was the first author that carried on the the study. Um these are maybe the data that you have seen last year. So the point here is different from structure. So the point is this is another protein. This is the corismate mutase. So the corismate mutase is an enzyme. It's 96 u amino acids long. Um and it is used by E.coli in a pathway that is needed to produce some some amino acid. So if you grow the bacteria in a medium where this amino acid is not present, the bacteria have to produce it and they need a functioning chismate mutase. Otherwise if you give to the bacteria the amino acid in the medium they can eat it and they don't need the corismate mutase. So you can select your bacteria by starving them of that amino acid in such a way that they have to produce it. But now this is really function. this is not structure. So you you can have a protein that is able to fold in the correct uh structure but maybe for some reason it doesn't work in E.coli. So this is shown here. So what they did here as a first uh test of the platform. So here you have again you have an alignment of sequences. you have uh I think in this experiment at the time they had uh of the order of 1,000 sequences. So here you have 1,000 sequences that come from different species. You take you take the the DNA I mean you you synthesize the the gene you put it in a plasmid the plasmid is put into the bacteria. So now the bacteria have this gene instead of the the original one of E.coli. They have another one that you take from another species and you let them grow in this medium where they don't have the the amino acid that is produced by the githic mutase. So you they have to produce it. So if the chismetic mutase works but now works is not just folding into a given structure. It means it it it falls it can interact with the other um partners that that are needed in this um in this uh pathway and so on and so forth. Uh then the bacteria will survive. So what you have here on the x-axis is the enrichment. So they put all the they put all the mutants together. So you have a library of here 1,000 natural corism mutasis taken from all possible species. You put them you make a culture with bacteria that have different genes different corism mutasis and you let them grow and you measure after some time the uh enrichment of each of the mutants in the population. So if the enrichment is uh let's say this is normalized to the wild type the wild type means the ecoli will be its own corismate mutase. So if it is close to zero it means this mutant is as good as the wild type. If it is minus three here it means that basically that mutant disappeared from the population after some time. So it's not working. Okay. If you look to if you make a histogram of this data, you see that it's quite bimodal. Okay? So you have a I mean you see it from the scatter plot but you see it better in a histogram. So they put a threshold here and they say okay below 10 the minus1 it's uh okay actually I'm not sure if it is 10 - 1 or e to the minus one but anyway you put the threshold here should be 10 to the minus one I guess but uh if you put the threshold here and below the threshold you say okay this mutant is not working so in this sense now this is nonfunctional in the sense that you are taking a protein that comes from another species that is able to fold and function in that species because it was selected by nature. But once you take it and you put it in ecoli in that specific experimental setting, the ecoli dies and you don't know really why it's uh it's not working. It could be because the machinery of ecoli is not able to express that protein or maybe it it is expressed it falls but then it's not binding to the other protein of ecoli that is in the pathway or any sort of reason. Okay. So there there can be many reasons. Anyways, 30% roughly of the natural proteins in the alignment work in E.coli. Okay. Uh and this is the x-axis. The the y-axis is the energy that the model gives to these sequences. Now the model is trained on these sequences. So because you train the model by maximizing the likelihood so minimizing the energy of the natural sequences as you expect all the natural sequences have low energy. Okay. So here you have low energy 30% of the sequences work. Now what they did here was to now you to you sample new sequences from this probability. when you sample you can change the temperature you can put a temperature here in front and in in this way as usual in statistical physics you can change the energy okay so let's say let's say I I I say this is my energy so this is e to the minus energy so now I can sample from e to the minus beta times the energy And by changing beta I can change the energy of the sequences I'm sampling. So so in this way they could sample sequences at different energies and then they tested them in the same experiment and what you see is that among the sequences that that that have low energy where low here there is a threshold that is more or less the biggest energy of the natural sequences. So among these again 30% are functioning. Um is this trivial? It is not because if you think to the number of sequences that you have, you have 21 to the 96 sequences. So it's a huge space. If you take a random sequence of of 96 amino acids, it will never work. Um so the fact that the model is able to generate random sequences such that 30% of them work is uh is really non-trivial and in particular if what they show in this experiment is also that if you don't put the J and you just use the conservation profile the success rate is zero. how many sequences you can generate this to to have an idea what you can do is to compute the the entropy of the model. Okay. So from the okay let's say you fix the temperature to one one is more or less the temperature at which you sample in this region here because it's the training temperature and you can compute the entropy. If you compute the entropy, what you find is that this entropy is of the order of one the entropy divided by the number of uh state of um amino acids is about 1.5. Now the uniform distribution has an entropy equal to logarithm of Q which is 21 and 1.5 is about half half that is about logarithm of 21 / two. So this means that the number of sequences that you can generate from this model is about so the the number of sequences that that you can generate that have the right energy is of the order of uh Q to the L / 2 because the entropy is divided by two. So the number which is the exponential of the entropy is the square root. This is still huge. Okay, if you put the numbers, you get 10 to the 65. Okay, so you have of the order the model can generate of the order of 10 to the 65 sequences out of which the experiment tells you 30% are functional. So the number this is the number of generated of the sequence this is the number of sequences you can generate the number of functional sequences is going to be 0.3 * 10 65 which means still huge. Okay. So the point is if the success rate is of order one then compared to the entropy it's nothing. So this means that >> sorry >> yes >> how did they test they are functional >> this is the test so you generate 1,000 >> and you and you see how many of them are functional and for >> they actually synthesize them >> yes so you you gen you sample from the model now you have sequences you synthesize them as you as you did for the natural ones the natural ones you have to you know the sequence because someone sequenced the corismate mutase in some other organism that is not ecoli. So what they did here is you synthesize the the this natural sequence and you put it in a plasmid and the plasmid is put in ecoli and then you run the experiment. So you do exactly the same thing here and here except that here the sequences are the training ones that you get from nature and here you're generating new ones and you do it for 1,000 and because the success rate is of order one you it converges to a value that you can estimate it from 1,000 but then you can say okay if I did 10,000 it should be the same. Okay. So you conclude that there is a huge space of sequences that work as well as the natural ones. Okay, this is interesting because it tells you that nature is not really you know finding is not solving a very difficult problem. I mean it's there are a lot of solutions to the problem of of making a functional charismatic mutase. So if you sample if you have an evolutionary process that is going to sample the space it's going to converge pretty fast to to some solution. Um but now the question is how different are the artificial ones from the natural ones because of course one could say okay but if the model is just proposing one or two mutations. So this is this plot. In this plot you have in the x-axis the identity to the closest natural sequence and what you see is that you have a range I told you that if you take two typical natural sequences so if you take two of these sequences at random the identity is close to 20%. They have 80% of the sides that are different. So you are more or less somewhere here. Now when you generate from the model and you look to the closest one uh you see something like this which means that the if you make a sketch of this it means really that the energy landscape if you make a plot of the energy as a function of the sequence this is of course the sequence is a high dimensional thing but this the sketch, but you really expect to have something like this where the natural sequences are staying like in the low energy part and when you generate artificial sequences, the more you go far away from the natural sequences, this is the energy and the more the energy goes up. So if you move away from the natural sequences, you go up in energy, which is not surprising because the model is constructed to give low energy to the natural sequences. And the more you go up in energy, the more it's unlikely that your sequences work. So the red dots are the one that work and the gray dots are the one that don't work. But you see that you can go as I mean that this these wells are flat enough here. I I my drawing is not very good. They're flat. So you can go up to 40% mutations away from the from the closest natural. Here you have a few that are even at 60%. So you are generating things that are very different from the natural. Typically here this is a 100 sequence long. So if you have here you have 40 mutations and you can go up to 60 mutations to the closest natural and still um it works. Okay. So this is this is something that I think you should always do when you try to make a generative model. Check the distance of what you generate from the closest training point. And people not always do it. So it's important to have an idea of how far you are from the train set and here you can go quite far. Okay. Is uh is this clear? So in this sense you then what I think you did in your project last year was that once you have this data you can use it to improve the model. So you can say okay now I tested these sequences. I know that these ones are functional and these ones are not. So I can improve my model by integrating this information. So just to give you an idea of the kind of things you can do with the ah sorry I forgot maybe before going to this let's I I wanted to show you if you want to use this you you have this other notebook in the drive. So, okay, maybe I will not uh well, let's try to run it. Let's see if we can run it um for PF14 just to give you an idea. So we have this we have this package um to do the bolts learning that is uh a Python uh package that you can uh use. So ah something went wrong. What's the issue? Yes. So I I read the same data. I read the the PF14 alignment. And then you you can install uh this package is called ADA BMDCA PI where DCA is direct coupling analysis. BM is BA machine. ADA stands for adaptive because you can also in this package you can also use this procedure of adding or removing couplings to make your model sparse which we are not going to do here but okay so you can you can install your uh your u model and then okay it's done. I I just want to run it to show you what happens when you run it. But so, okay, hopefully you can get the GPU from from Google. Um, and then, okay, you can import the the data. You see here it tells you that you have 8,871 sequences of length 53 with Q equal 21. And now um you prepare the system. So uh something that I forgot to say is that when we do this Bzman learning we give a weight to the training set and the weight is typically calculated by checking how many sequences you have around you at distance smaller than 20% which is a convention but it's found to work well. So you count how many sequences you have that are closer to 20% identity to you and the weight is one over this number. And the idea is that you want to give less weight to sequences that are close to many other sequences because like if you have many sequences here, you want to downweight them a bit because the idea is that these are all maybe related by fogyny, you can have bias in the in the alignment. You can have a fellow genetic bias that maybe there is a big cluster of sequences that are close or maybe for some reason people sequenced more of this because it's interesting. I don't know it comes from a species that is well studied. So you have more sequences. So you have this kind of biases. So you want to downweight a bit the regions that are too populated. So this will this is done here and then you get an effective number of sequences that is the the average weight that you have one over the average weight. So you get the an effective number that is here it's about half the original one. Okay, we put a small P count that is in this case is very small is one over 4,000 compared to mean field where you put the P count of order 1/2 and then with this P count you compute the FIA and the you compute the empirical frequencies. Now I will call them fi a and fi a jb. So these are the expectation over the data of delta sigma i a. So you do this once at the beginning of the learning and then you you you run. So here you can choose how many Monte Carlo chains you you want to use to sample from the to compute because you you remember in order to do the BSA learning you have to compute these quantities from the model and to do that you need Monte Carlo. So you you can choose the number of chains. You can choose how many Monte Carlo sweeps you do at each gradient update which is an important parameter because when you do Monte Carlo you want to be as uh close as possible to a good sampling. So you want to be in equilibrium and in order to equilibrate you need to do a certain number of sweeps. So this number depends a lot on the family and if you're not careful and if if this number is too small you can be uh you can have a problem of convergence of the BS and learning but okay this would require a long discussion and again I you can check I can give you references if you want but anyways for the for this PF14 the Monte Carlo converges fast enough so 10 sweeps is good you can fix a target Pearson which is the stopping criterion. So the the learning will stop when the CI of the data and the CI of the model have a Pearson coefficient of so as soon as it becomes bigger than here.9 you stop. Usually we do 0 95 but if I do 0 95 it will take half an hour to run. So uh well and then you have to fix a learning rate that you can choose here is 5% is good enough and finally you can choose a sampler you can use Gibbs or Metropolis here and well it yes >> can you repeat why we're doing >> yes because when we do Bzman learning we want to update the parameters of the model so we start from in this code we typically start from a model that has zero couplings and the fields are fixed by the frequencies of the amino acids. And then the update equation is something like h i a t + one is equal to h i a of t plus a learning rate times the expectation of this let's say f i a minus p i a of t where fia a is the expectation over the data and pia is the expectation over the model. So you need to compute these. So you need to compute these are essentially the magnetization of the model. So you need to compute an average from a Bzman distribution that is complicated because you have uh parameters that depend on the sites and so on and so forth. So you cannot compute it exactly. So you need to sample. So the Monte Carlo is is done in order to compute these quantities. Basically you generate here you have um 10,000 chains. So you generate you run the Monte Carlo with 10,000 chains. You generate 10,000 sequences and then you compute the probability that over this 10,000 sequences that you have amino acid A on site I and this is your estimate for P. Okay. Okay. Now the training works like this. Okay. I am perfect. It's not working. Missing two required positional arguments, but they are there. H maybe there was a let's see sometimes the we update the package. Okay, there was an update of the package. Okay, so you see here how it works. You the the training runs with epochs. Each epoch is one of this t. So now we are t is 35 3 5 3 5 3 5 3 5 3 5 3 5 3 5 3 5 3 5 3 5 3 5 3 5 3 5 3 5 3 5 3 5 3 5 3 5 3 5 3 5 3 Every epo is an update of the parameters. At each epoch you do the Monte Carlo and what you have here is the log likelihood that is given by the log of this averaged over the data and you have the Pearson here. So you want to reach 0.9. So here we are at 0.5 and it's increasing. Um to the log likelihood is actually the true log likelihood because in this code we are also able to track zed and the way we do zed is that because you change a little bit the parameters at each step. You can write you have an expression for log zed at time t + 1 which is log zed at time t plus a correction and you can compute the correction. So while you learn you can also track zed. So you can have the the full likelihood and the entropy of your model as a function of time. So okay something that I can tell you is that if you plot here the plot will come later but if you plot the person if you plot one minus the person versus time you typically get something like this. So the you want the person to be as big as possible. You want it to go to one ideally. So at the beginning it is zero because you don't have any coupling. So the correlations are random. So it typically does something like this is in log scale and this is in log scale. It will typically do something like this. So it starts the person will start growing and so this goes down and then it will go to one. So one minus row will go to zero as a power law of time. So at the beginning it's fast and you see okay it's converging very fast. But then you see that now it's saturating a bit and it takes quite a lot of time to go down. Okay, you can play with the learning rate, you can do things. We try to optimize uh the code as much as possible but we didn't find the way to go faster. So you have to wait and so well you see now we are about 08 and uh in in two or three more minutes it will converge. So uh once you convert okay this is a plot okay this this is the plot I just uh made. So this is these are the epochs this is for another run but okay and this is one minus the person. So you see you get this power law behavior here. Um we can maybe well this will take two or three more minutes to converge. So in the meanwhile I can show you what what you can do once you have the the data. So first of all you can check once you learn your model you can check that you did right by uh running the Monte Carlo again and computing the time it takes to the correlate from initial state. So this is shown here. uh and you should check that the number of sweeps that you used in the train. So I forgot to say one Monte Carlo sweep means that you try to change each amino acid once. So if you have 100 amino acid, one Monte Carlo sweep is 100 Monte Carlo steps one per site. So you should check that the number of sweeps that you used in the training is this is what is shown here is the sequence identity between the sequence at time zero and the sequence at time t. Uh and this is compared with the sequence identity between the sequence at time t /2 and the sequence at time t. So if you the correlate then zero t /2 and t should be independent and so these two quantities should be similar. You see at the beginning they are not similar and t when time increases it goes uh it goes down. Okay. So uh sorry I I no I I said it in a wrong way. So sequence identity at time t here is the sequence identity between two different chains at time t. And this is the same chain at time t and t /2. So if you have if you have you're doing your monte carlo you start from the chains typically start from either from the natural data or from some sample that you give it depends I mean but say you start your chains then your chains will evolve in time. and you look at the sequence identity between two different chains at time t and you look at the sequence identity between the same chain at time t and at time t /2. So if your time is bigger than the mixing time this should be an equilibrium sample this should be an equilibrium sample and this should be an equilibrium sample and they should be all all independent. So if this then these two quantities should be the same. So this is if the time is long enough you see that they are the same but this is for another protein. So let's see if I can update it's still it will take a few uh maybe one more minute and then we get our result. So once we have that you you can measure a few things to check that you're generating right. So this is a PCA plot where we have the in black we have the natural sequences. So we we take the natural sequences we compute the C i J A B matrix we do PCA and we project on the first two principal components and then we do the same with the uh generated sequences that are projected on the principal components of the natural sequences. So the principal components are the same here and I'm projecting the two data sets and you see that they are very similar. Also you can make a histogram of the principal components and you see that they are similar and then you can compute the one point frequencies and it's not shown because they are basically perfectly correlated because the fields are always very easy to fit. So the conservation is always very easy. What is more difficult is to fit the CIG. So this is what you get typically at the end of the training. You get this is the correlation between the CI of your natural data, the training set and the CJ of the generated data. And then you can compute and these are this this is um this is fit by the model, right? So this should be close to the diagonal by construction. But then this is not a fit because these are free sight correlations. So you take I J K and you ask what is the probability that I have three amino acids in three different sites. You do it for many triplets. Here we choose because the there are too many triplets. We just choose a few at random. And you can see that you have a good fit also of the three-point correlations that you didn't put in the model. So this is just the model is able to catch. >> Sorry. >> Yes. >> Yes. The package is there is a routine in the package that does it. And then here you get a histogram of the energies. I want to show this because this is important. This is the histogram of the energies of the in blue you have the natural ones and in yellow you have the ones that you generate and you see this is not perfectly matched because by construction the model is trained to maximize the likelihood of the natural sequences. So it always gives a slightly lower energy to the natural sequences which means the one the natural ones will be always a bit on the bottom of these wells and the ones that you generate will always be a little bit higher in energy. You can regularize this a little bit. So if you make the model sparse, this will improve. You get a better match between the two. But it's never perfect because the model is always a bit aware of the fact that the blue are the training set. So by construction of the maximum likely uh oh this is done. So we can see quickly and then we conclude uh so here it uh it finished. So we can do this plot. Okay. This is for the PF14. This is the Pearson versus the number of epochs. And you see the power law. It's even more visible than for the other family. And now you can sample. So when you sample, you have to tell how many sweeps you want to do to now we are fixing the model. The model has been trained and we are generating new artificial sequences from from the model. So this will uh run the sampling for a certain number of sweeps and so at the end you get your artificial sequences that you can compare with the natural ones. Now this will take a little bit of time. So while it runs we can uh see the last things. So this this I was saying is the histogram of the energies and here you can check the distances between sequences that is important to know if you're generating well. So here you have MSA1 is the training one the natural ones and MSA 2 is the generated ones. So this is the histogram of pairwise distances this is not the PF14 so it's longer than 53. uh but when we will redo it with the 50 the PF14 it will be different but anyways you see that if you take two training sequences at random they will have a certain distance and if you take two generated sequence at random they have the same distribution of distances but what is more important is that you can compare you take one of the sequences you generate and you look to the distance between this sequence and the ones in the training set and you get a histogram like this. And then you do the same but you take one sequence of the generated set and you look to the closest sequence in the generated set and if you're generating well these two should be the same and they are very similar. Okay. So it means that the things that you are generating are between themselves at the same distance than they have from the training set which means you are not your that are not closer to the training set than they are uh between themselves. So this is a good measure of um generation quality. And here the last plot that you have in this notebook is the same as the one I was showing before. This is the distance to the closest natural and the energy. So you see that the more you move away from the closest natural, the more the energy goes up. This is for the generated ones. And of course for the natural ones, the distance to the closest natural is always zero because it's the is each sequence is a distance zero to itself. But this gives you the range of energies that is comparable. So the generation should be done. I can show you just the results for the PF. I know it's not it's not done. No, it's done. Well, okay. You see this is the this is the autocorrelation curve. This one for the PF14. And you see that after 50 time steps you are very close to the to the asytoic limit. In the training we use 10 which is a bit short but because you your learning rate is small it has enough time to decorrelate so it's a good u and well once you ah no I should not do this but okay this will run again. Okay anyways you can run the notebook by yourself and you can check the the results. So with this you can choose your favorite protein family and you can run BMDCA you can get the parameters. So you can do contact prediction if you don't if you want to check uh the history of these things because now you do contact prediction with alpha fold but okay uh but more importantly you can generate artificial sequences and you can do what you want with them. Um the last thing that you can find in in the here that I didn't have time to discuss but uh can be done is to use these models to do Monte Carlo short times and compare with experiments in in in uh in vitro. So we have a few things about it but uh I I didn't have time to discuss but you find the notebook in the in the folder that does this so you you can run the simulation for short times and you can compare with experimental data and well I I conclude here but if you have any questions on these things uh I'm always happy to chat so send me an email and uh we can discuss very very good. >> Ah uh and Sergio said he will correct the exam at 4. >> Okay. Ah yes. So if you if you want to know >> Sergio will be here if you want to correct a correction discussion of his exam uh you can come and discuss with him. Okay. Very good. So uh so this is the end of the lectures in the spring college. Tomorrow we will have the exam of um uh Francesco and uh so I I would like to uh um propose uh one thing because uh while we were here in this four weeks uh looks like the world outside went crazy and uh and So um maybe I mean so these are things uh there are people here coming from different parts of the world. So and uh so what happens typically is that uh uh we sort of tend to maybe uh discuss these things with people who uh who think like us. So I wanted to see maybe this is a crazy idea but to have a um a session tomorrow at five who wants to come just as a uh like uh re-equilibibration chamber just to express uh different views how we see what is going on how we live these things because uh maybe this is a opportunity we don't have normally I mean to to to have such a broad range of people of course this is completely unrelated to the spring college this is completely informal uh don't yeah so um has nothing to do with physics but uh but I think it can be useful okay especially because uh well when these things these crazy things happen we tend to be very much focused on the short term but uh especially for you who are younger I think what is important is the long term okay and and then not to be too much uh uh I mean your trajectory not to be too much affected by I mean I mean uh as little as possible by by this thing. Okay. So tomorrow at 5 um for those of who you want we will be here and uh it will be just uh I mean the the only thing is uh of course uh we want to avoid any type of confrontation. So uh no um uh reactive uh things. I mean to uh we just want to I mean uh have a have a safe space where everybody can say uh what he thinks. Okay. Okay. That's it. So thank you very much. See you tomorrow at 9. Okay.