If they are independent the residues in one of the networks mutate in sync with each other, but the residues from the other network mutate in a different rate. For example, If you compare two subfamilies the positions of one network all have a different residue in each of the two subfamilies (so they mutate between the two subfamilies), and the positions of the other network might be 100% conserved over the two subfamilies.

You can explore this idea by making an alignment with two groups of proteins. Each group should contain functionally similar proteins but this function differs between the groups. For instance, a group can be made with enzymes all converting the same substrate, and a second group that convert another substrate. The correlated mutations then often reveal hotspots that determine the substrate specificity change. A very important note is that correlated mutations reflect functional changes and not environmental changes. This can easily be explained: You cannot use correlated mutations to find hotspots for thermostability. You could make an alignment of two groups: one composed of proteins from thermophilic organisms and another made from non-thermophilic organisms. Then the differences between these groups should show up as correlated mutations, right? This is not the case, because in fact these two groups cover all sequences of the complete superfamily alignment. You can investigate the differences between these groups just by comparing the amino acid contents between the groups. This concept will be discussed later. One more important note: The problem with functional grouping is that the functional annotations of proteins is often based only on sequence similarity and can therefore be wrong. You need to make sure that at least a high percentage of the sequences share the same function to make functional based hotspot finding work. There are examples though were the alignment was successfully grouped simply by using keyword searches in the protein descriptions and used to find novel enzymes. The group of Prof Uwe Bornscheuer has found R-selective amine transferases (Only S-selective enzymes where known) by deleting from the protein family groups of proteins that had other activities based on motifs that were associated with the different functions. This exercise was repeated with 3DM. The protein family was first grouped by giving keywords that were related to enzymes with unwanted reaction mechanism. The keywords lyase, for instance, was used to select of enzymes with likely the unwanted lyase activity. Then, this set of sequences was used to find a sequence motif specific for lyase activity. All sequences that had this lyase specific sequence motif were subsequently deleted from the alignment. This step ensures that all lyases that were not annotated as such still get deleted from the alignment. This procedure was repeated for all different unwanted enzyme activities. This exercise resulted in 42 sequences that most likely are all R-selective amine transferases.

There are multiple ways of doing this. One easy way would be to look at the amino acids distributions in each network. If the number of different motifs/different groups in each network differ (e.g. the positions of the first network all have only two different residue types and the positions of the other network all have five different residues) you know the two networks independently mutate. See question 6. Structural separation of the two networks might give a clue too.

Yes, one network clearly contains hotspots for specificity. This means that changes in specificity where the driving evolutionary force behind the correlated mutation behavior of this sub-network. For the other network there is no literature in the 3DM database describing mutations effecting specificity. As this database contains huge amounts of mutation data it is unlikely that these positions are hotspots for specificity.

Yes, it is very likely they are also important for specificity. You could have a look in the 3D structure to see if they are close to the other members of this network (and/or close to the ligand binding pocket). But the only way to really make sure is to make mutations in the lab and test the binding affinity of the ligand. In enzymes such mutations would have an effect on the Km but not so much on the Kcat. Such mutations typically indicate a hotspot for specificity. The Km on other substrates might improve at such positions.

Usually 3DM alignments are composed of sequences from organisms from the entire tree of live. If you would make a phylogenic tree of the whole alignment and you see that it nicely overlays the tree of live, then correlated mutations can reflect the different domains of the tree of live. But most of the time multiple sequences from the same organism are in a 3DM alignment and these sequences can sequentially be very different. They almost always have different functions. Because there are proteins with different functions in an alignment the different groups of the phylogenetic tree are function based. This means that proteins with the same function form the large groups/clusters in the tree instead of the different domains or kingdoms of the tree of live. These different groups in such alignments result in correlated mutations that reflect the functional changes during evolution. Note that you can have both types of differentiation in one alignment.

The pie chart of position 74 shows that there are two residues (W and F) very abundant at this position. Note that the * in the pie chart indicates all residues that are less than 5% present in the alignment. The variability at position 88 is much higher. 8 different residue types are available in more than 5% of the sequences. The two networks clearly show different mutation rates. To get the exact distributions you can click on "open position info page" below the pie chart.

Yes, most positions in the network of 74 have only two types of residues. The positions of the other network all have many more. This clearly shows that these two networks are independently co-evolving.

The amino acid distribution of these two positions indicates that it is possible that they form a salt-bridge. In a large portion of the sequences they have opposite charge indicating they might structurally be close to each other and being able to form a salt-bridge. On the other hand there are quite some sequences that have E-E. This might indicate that they don't form a salt-bridge. You can simply check this if you select both positions and then click on "selected nodes" below the pie chart.

The E-score is quite low. An E-score of 1 means that there is no enrichment for the specificity related mutations in the network (just as likely to find specificity mutations at a position in the network as a position outside the network). There is nothing wrong, however, because not all positions of the complete network are hotspots for specificity. Only one small sub-network seems specificity related and thus the E-score for the complete network is low.

The E-score for the sub-network is 5.0 (make sure the correlation cut-off is still set on the default of 1). An E-score of 5.0 indicates that the chance that the positions in this sub-network are specificity hotspots is 5 times higher than random. Or better, the chance they are specificity hotspots is 5 times higher than other positions. In other words: because these positions are important for specificity there has been an evolutionary pressure at these positions that has resulted in the alignment in this correlated mutational behavior. This is why correlated mutations are often referred to as "co-evolution".

Fig 1. This figure shows the E-scores of different keywords. These figures are always availabe below the CorNet network. At high correlated mutation cut-offs specificity mutations are clearly over-represented in the GPCR CorNet network.

If set on 2 then 3DM will only show positions in the network for which at least two different mutations have been published that had an effect on specificity. This is to ensure no false positives are included. In protein families for which a lot of mutations are available it is usually smart to leave this at at least two. If there is not so much data available it might be better to set it on 1 to get reliable E-scores.

As position 93 is part of the sub-network that has a high E-score for specificity it is very likely that 93 also is a hotspot for specificity. For position 248 on the other hand this is less clear. Note that 3DM "just" gives you the data. To make sure, always read the article. Especially if there is only one article found for a position.

The specificity related network does cluster in the structure. The other big network doesn't cluster structurally and the function remains unclear based on structural data.

Considering the 3D location (near the ligand binding pocket and close to the specificity related network) position 60 probably is also a hotspot for specificity.

Just like 60, 247 and 248 are located close to the light blue positions of the specificity related network and are therefore likely specificity related positions. The role of the other residues are not clear based on this exercise.

There are several compounds crystalized in 1F88A. One of them is located inside the GPCR. This ligand is located exactly in between all the specificity related positions. This is a good example of how to use 3DM. With 3DM you are looking for correlations between different data types. In this case there is a good correlation between the CorNet network and structural location. If you find something like this your alarm bells should start ringing. In this case this screams specificity. If correlated mutations are located around the ligand (or substrate in case of enzymes) they almost always are specificity hotspots. Other examples of correlation between different data are: 

1. High B-factor & High RMSD → hotspots for thermostability

2. Motifs & sequence groups. For example: A sequence motif that is highly conserved in a subset composed of sequences with the same function and a different, again conserved, motif is found in a subgroup of sequences with another function. Often these positions determine the functional swap between the two subgroups (e.g. specificity hotspots)

3. Ligand contacts & activity data of ligands. Some positions are bound only by inhibiting- or activating ligands → hotspots for inhibition and/or activation. A nice example was shown in the NR practical.

4. Can you think of another example?

38,324

45,796 for the full alignment. 2,940 for class B, 3,978 for class C, and 4,008 for class F. Note that these numbers might be outdated when the database has been updated.

There can be different reasons. 3DM sometimes deletes sequences during the alignment generation procedure if they can't reliably be aligned. It might be that some sequences can reliably be aligned in the smaller family alignments because they have a bigger core. Meaning there are more residues that can be used to make the alignment resulting in more reliably aligned sequences. In the complete superfamily 3DM alignment only the trans-membrane helices are aligned (as those are the only parts that are structurally conserved between the different families). This means that there can be more uncertainty on the quality of the alignments. Another reason might be that some sequences are aligned in more then one class.

The first conserved residue (>95%) is Cys at position 68 in the full 3DM alignment. The same conserved Cys has number 81 in the class A family. There are less structural conserved positions if a structural conserved core is made for all the GPCRs. There are more positions structural conserved if only class A GPCRs are used to determine the structural conserved core.

They have the same number using the class A numbering scheme: 3.25a. In the class B alignment this residue has number 3.29b. Obviously there are more structural conserved positions before this residue in the class B alignment compared to the class A alignment. When this common numbering scheme was designed they didn't realize that they could have made a numbering scheme that could be applied to all classes. As everyone uses this numbering scheme now, it is too late to synchronize the numbering schemes of the different families.

If you have searched in the titles of the articles you will find that there are two human proteins known to be involved in Hirschsprung disease: J7HDH6 and P24530.

Alignment numbering: This uses the 3D numbers. Note that the active numbering scheme is used. This can also be a user defined numbering scheme. Later you will play with the different numbering schemes already available in the GPCR protein family. 


Residue numbering: This option uses the numbering of the sequence that is being modeled. 


Template numbering: This option uses the numbers of the template pdb structure.

If you search with the keyword "Hirschsprung" with the literature cut-off set on 2 you will find that at 7 positions the literate reports Hirschsprung related mutations. In Pymol you can use the "Show scene content details" to get a list of the search results. Note that using "hirschsprung disease" as keyword doesn't change anything. But, sometimes you have to think about which keyword exactly to use to find what you need. For instance, if you search for mutations that resulted in a more active enzyme then you can use "increase activity" but also "high activity". The search works such that the words don't necessary have to be adjacent. So, a sentence like: "the A345Q mutant resulted in a enzyme of which the activity was increased", will also be found with the keyword "increase activity". If you put the literature cut-off on 1 you will see that there are many more positions mentioned in the literature in relation to Hirschsprung disease. 





Take home message: Always think about which keyword (combinations) and cut-off best to use.

In this picture the first three templates were selected, the top five correlating positions, the first three ligands and all positions >90% conserved. You can see that the correlated mutation nicely surround the ligands located inside the GPCRs.

SNPs are less common at conserved positions. With these cut-offs you see that the Enrichment is 0,55. These means SNPs are found less at conserved positions than at other positions. Note that it is important to set the cut-offs in this exercise because otherwise too many positions fall within the cut-offs to calculate a meaningful Enrichment-score.

Not really. Maybe position 84. As SNPs are from genome sequencing projects of healthy humans you don't expect many SNPs at important positions as they might be damaging to the functioning of the GPCR, which, in turn, can result in a disease. This information can also be used the other way around. At positions where you find many SNPs it is likely that mutations are allowed. This means that they might be good candidates if you want to change the thermostability, because you can be quite sure that mutations won't harm the functioning of the protein.

There are 7 positions: 29, 85, 130, 134, 145, 195, and 205. Position 85 is a common SNP closest to the ligand binding pocket.

The result should be the same.

If all OK, you should have these positions in your basket: 84 88 92 93 98 109 218 245 248

The B-factor method doesn't seem to be a good method for the GPCRs. Likely this is because the transmembrane helices are so tightly bound that automatically the rest of the protein is much more flexible. It is known that making mutations in the transmembrane helices can make GPCRs more stable. This is nicely demonstrated by the "thermostability" literature search. Clearly the two methods do not overlap indicating that the B-factor method does not apply to GPCRs. It is in this protein family probably much better to first try the positions that are already published to have an effect on the stability. The best way to do this is to find amino acids described in literature that are known to stabilize GPCRs. If the stabilizing residue is not in your target GPCR then it might be a good idea to try that residue. If your target sequence has a different residue than the consensus it might be smart to try the consensus residue as there are several papers demonstrating that making the consensus often has a beneficial effect on the stability. (Note that if it is easy to screen your protein for thermostability than randomizing each hotspot might be smarter). There are many other tricks too. Here are some examples: 

1. Introducing prolines at positions where a proline is common and your target doesn't have a proline

2. Inserting negative charges at the N side of a helix (if there isn't any). Or a positive residue at the C-side. This is known as helix capping because the N terminus of a helix is slightly positive and the C-terminus is slightly negative.

3. Creating salt-bridges by inserting positive or negative residues at positions on the outside of the protein. Always check if the residue you want to use is actually common in the alignment at that position.

4. Replacing glycines where, according to the alignment, glycines can be replaced by something else.

1267 groups.

Motifs for which only one sequence is available are not put in a group. The hypothesis behind this is that those might be sequence errors. Also sequences that have a gap at a selected position are not grouped. If you put "minimum group size" on 1 and deselect "hide motifs with gaps" the number of motifs is equal to the number of groups.

There are now 1436 total groups.

In the abstract you can find R124. This R is at position 30. Remember that articles, of course, do not use the 3D numbers. In 3DM you can see which number is used in the article. If you go back to the 3DM page that gave the link to the paper you see that at position 30 mutation R124Y has been found. Remember the residue number when you go reading a paper. This makes it easy to find back the residue in the paper. Clearly this is a highly variable position and likely you can mutate this residue without destroying the function too much. Furthermore, R is not a very common residue. Even more "evidence" that the R can be mutated to something else. There is literature available that making the consensus often has a positive effect on the stability of proteins. Although we don't see this as a very common rule, in this case we do see that mutating to the most common residue has a positive effect on the stability.

• Go to the protein detail page of P24530, make a model of this sequence and open it in Yasara or Pymol. Remember you can find a generated model at the "Visualize" option of 3DM.

• In Yasara: look at R124 → use F3 to see all the atoms. Make the R124Y mutation (right click residue → swap → residue).

• In Pymol follow the following steps:

  • Click on the ‘H' next to object in the object list. Select "Everything". The molecule will disappear.

  • Click on the ‘S' next to the object. Select "Lines". The structure will reappear with individual atoms visible.

  • In the top menu, click one Wizard > Mutagenesis. Click on the residue you wish to mutate.

  • In the Mutagenesis menu in the bottom right, click on the first button, which will normally say "No change". Here you can select the residue you wish to mutate to and click "Done".

• Do you understand why this mutation might be stabilizing? Investigate all other mutations that are described in the abstract of the paper. Look at conservation and structural location to see if you would have picked these positions/mutations in order to change thermostability.

It is almost always very difficult by simply looking at a structure to explain the effect of a mutation. It is even more difficult to predict the effect by just looking in a structure. This is why it is so handy to have the 3DM data behind the structure. But let's just give it a try:

R124 (3D number 30) might be a bit more unstable because there is another positive residue (K34) very close to R30. K34 is actually really been pushed away from R30. The mutation to the Y makes it more stable. Also because the hydrophobic part of K34 nicely packs to the hydrophobic ring of Y30 (see figure).

D154 is not part of the core and doesn't have a 3D number. You can find this residue by looking in 3DM what 3D number residues have before or after this residue. 3D number Ile61 comes after the D154. Looking in the model you will find that D154 could not be modeled because the template doesn't have this residue. This could be solved by finding a template that does have a residue there, but 3DM cannot do this automatically and is therefore beyond the scope of this course. If you are interested in how to do this, find a good homology modeling course.

K270 (in 3DM: K163). This position is also highly variable. Although Ala is not the most common residue clearly the hydrophobic residues are by far the most common type of residue. So based on the amino acid distribution data mutating this Lysine seems a smart choice. Looking in the structure K163 doesn't really fit there and it is close to other positively charged amino acids. It's looks like mutating it to something hydrophobic will indeed have a positive effect on the stability and a smaller residue might fit there better as the side chain of K163 bumps into the side chain of L225.

S342A (3D number S220): This position is again highly variable. S in not very common at this position and just like K270 the larger hydrophobic residues are the most common residues. So mutating this residue seems a logical choice if you want to target thermostability. Looking in the structure there is not clear reason why Ala would be more stable than Ser.

I381A (3D: I250): This position contains mainly large hydrophobic residues. So mutating to Ala doesn't seem the best choice unless in our model the Ile doesn't really fit in the structure. Looking in the structure this is not the case. So why the I250A mutation is beneficial to the stability is unclear based on just these analyses.

The 'measured data' is randomly generated, so the best and worst performing mutations will be slightly different for everybody.

Your base sequence will contain some mutations from previous rounds and in this case, there's a 'V' at position 279.