Difference between revisions of "BIN-SX-Homology modelling"

From "A B C"
Jump to navigation Jump to search
m
m
Line 19: Line 19:
  
  
{{DEV}}
+
{{LIVE}}
  
 
{{Vspace}}
 
{{Vspace}}
Line 134: Line 134:
 
# To emphasize this better, select the protein atoms and display them as a '''sphere''' model to better appreciate the packing, i.e. the Van der Waals contacts. Use the '''Favorites → Side view''' panel to move the clipping plane and see a section through the protein. Study the packing, in particular, note that the additional methyl groups of a valine or isoleucine would not have enough space in the structure. Then restore the clipping planes so you can see the whole molecule.
 
# To emphasize this better, select the protein atoms and display them as a '''sphere''' model to better appreciate the packing, i.e. the Van der Waals contacts. Use the '''Favorites → Side view''' panel to move the clipping plane and see a section through the protein. Study the packing, in particular, note that the additional methyl groups of a valine or isoleucine would not have enough space in the structure. Then restore the clipping planes so you can see the whole molecule.
 
# Lets simplify the view: choose '''Actions &rarr; Atoms/Bonds &rarr; show''' and '''Actions &rarr; Atoms/Bonds &rarr; backbone&nbsp;only &rarr; chain&nbsp;trace'''. Then select <code>A&nbsp;42</code> again in the sequence window and choose '''Actions &rarr; Atoms/Bonds &rarr; show'''.
 
# Lets simplify the view: choose '''Actions &rarr; Atoms/Bonds &rarr; show''' and '''Actions &rarr; Atoms/Bonds &rarr; backbone&nbsp;only &rarr; chain&nbsp;trace'''. Then select <code>A&nbsp;42</code> again in the sequence window and choose '''Actions &rarr; Atoms/Bonds &rarr; show'''.
# Add the surrounding residues: choose '''Select &rarr; Zone...'''. In the window, see that the box is checked that selects all atoms at a distance of less then 5&Aring; to the current selection, and check the lower box to select the whole residue of any atom that matches the distance cutoff criterion. Click '''OK''' and choose '''Actions &rarr; Atoms/Bonds &rarr; show'''. You now have a very clear scene of the alanine residue in red, the surrounding side chains, and the rest of the structure as a C-alpha trace. You also see three water molecules. Spend a bit of time again, to get a sense for the spatial context<ref>Chimera uses a default '''distance to screen''' that is too close and that exaggerates the depth of the scene to a degree that it is difficult to fuse the stereo pairs. Choose '''Tools &rarr; Viewing controls &rarr; Camera''' and set the distance to screen to 50 cm. This will make stereo viewing easier and will also give a better match between distance estimates in all three dimensions.</ref>.
+
# Add the surrounding residues: choose '''Select &rarr; Zone...'''. In the window, see that the box is checked that selects all atoms at a distance of less then 5&Aring; to the current selection, and check the lower box to select the whole residue of any atom that matches the distance cutoff criterion. Click '''OK''' and choose '''Actions &rarr; Atoms/Bonds &rarr; show'''. You now have a very clear scene of the alanine residue in red, the surrounding side chains, and the rest of the structure as a C-alpha trace. You also see three water molecules. Spend a bit of time again, to get a sense for the spatial context<ref>Chimera uses a default '''distance to screen''' that is too close and that exaggerates the depth of the scene to a degree that it is difficult to fuse the stereo pairs. Choose '''Tools &rarr; Viewing controls &rarr; Camera''' and set the distance to screen to 50 cm. This will make stereo viewing easier and will also give a better sense of distance estimates in all three dimensions.</ref>.
 
#Select <code>A&nbsp;42</code> again: '''left-click''' (control click) on any atom of the alanine to select the atom, then '''up-arrow''' to select the entire residue. Now let's mutate this residue to isoleucine.
 
#Select <code>A&nbsp;42</code> again: '''left-click''' (control click) on any atom of the alanine to select the atom, then '''up-arrow''' to select the entire residue. Now let's mutate this residue to isoleucine.
 
#Choose '''Tools &rarr; Structure&nbsp;Editing &rarr; Rotamers''' and select <code>ILE</code> as the rotamer type. Click '''OK''', a window will pop up that shows you the possible rotamers for isoleucine together with their database-derived probabilities; you can select them in the window and cycle through them with your arrow keys. But note that the probabilities are '''very''' different - and thus show you high-energy and low-energy rotamers to choose from. Therefore, unless you have compelling reasons to do otherwise, try to find the highest-probability rotamer that may fit. This is where your stereo viewing practice becomes important, if not essential. It is really, really hard to do this reasonably in a 2D image! It becomes quite obvious in 3D.
 
#Choose '''Tools &rarr; Structure&nbsp;Editing &rarr; Rotamers''' and select <code>ILE</code> as the rotamer type. Click '''OK''', a window will pop up that shows you the possible rotamers for isoleucine together with their database-derived probabilities; you can select them in the window and cycle through them with your arrow keys. But note that the probabilities are '''very''' different - and thus show you high-energy and low-energy rotamers to choose from. Therefore, unless you have compelling reasons to do otherwise, try to find the highest-probability rotamer that may fit. This is where your stereo viewing practice becomes important, if not essential. It is really, really hard to do this reasonably in a 2D image! It becomes quite obvious in 3D.
Line 419: Line 419:
  
  
We have spent a significant amount of time to prepare data for the analysis and in practice it usually seems to turn out that way, that the preparation of data occupies the greatest part of our efforts. The actual computational analysis is generally quite fast. And, unfortunately, the '''interpretation of results''' is often somewhat neglected. Don't be that way. Data does not explain itself. The interpretation of your computational results is the most important part. The integrator
+
We have spent a significant amount of time to prepare data for the analysis and in practice it usually seems to turn out that way, that the preparation of data occupies the greatest part of our efforts. The actual computational analysis is generally quite fast. And, unfortunately, the '''interpretation of results''' is often somewhat neglected. Don't be that way. Data does not explain itself. The interpretation of your computational results in a bilological context is the most important part.
  
 
{{Vspace}}
 
{{Vspace}}
Line 503: Line 503:
 
# Start Chimera and load the '''model''' coordinates that you have just renumbered.
 
# Start Chimera and load the '''model''' coordinates that you have just renumbered.
 
# Select all, hide Ribbons and show Atoms, bonds to view the entire model structure.
 
# Select all, hide Ribbons and show Atoms, bonds to view the entire model structure.
# Choose '''Tools &rarr; Depiction &rarr; Render by attribute''' and select '''attributes of atoms''', '''Attribute: bfactor''', check '''color atoms''' and click '''OK'''.
+
# Choose '''Tools &rarr; Depiction &rarr; Render by attribute''' and select '''attributes of atoms''', '''Attribute: bfactor''', check '''color atoms''' and click '''Apply'''. Note that you can change the way the spectrum is mapped to the values by moving the blue, white and red bars over the histogram with your mouse.
# Study the result: It seems that residues in the core of the protein have better energies (higher values) than residues at the surface. Why could that be the case?
+
# Study the result: It seems that residues in the core of the protein have better energies (higher values) than residues at the surface, i.e. Swiss-Model was more confident in the predicted conformationstes. Why could that be the case?
 
}}
 
}}
  

Revision as of 14:19, 31 October 2017

Abstract

This unit introduces the principles of modelling structures based on the known coordinates of a homologue. The key to sucessful modelling is a carfully done multiple sequence alignment.


 


This unit ...

Prerequisites

You need to complete the following units before beginning this one:


 


Objectives

This unit will ...

  • ... introduce the principles behind homology modeling of structurs;
  • ... teach how to produce a structural model of the MBP1_MYSPE APSES domain;
  • ... demonstrate how to analyze the model;


 


Outcomes

After working through this unit you ...

  • ... can produce a homology model using the Swiss-Model server;
  • ... can work with Chimera to analyze its structural details.


 


Deliverables

  • Time management: Before you begin, estimate how long it will take you to complete this unit. Then, record in your course journal: the number of hours you estimated, the number of hours you worked on the unit, and the amount of time that passed between start and completion of this unit.
  • Journal: Document your progress in your Course Journal. Some tasks may ask you to include specific items in your journal. Don't overlook these.
  • Insights: If you find something particularly noteworthy about this unit, make a note in your insights! page.


 


Evaluation

Evaluation: NA

This unit is not evaluated for course marks.


 


Contents

Task:

  • Read:
Biasini et al. (2014) SWISS-MODEL: modelling protein tertiary and quaternary structure using evolutionary information. Nucleic Acids Res 42:W252-8. (pmid: 24782522)

PubMed ] [ DOI ] Protein structure homology modelling has become a routine technique to generate 3D models for proteins when experimental structures are not available. Fully automated servers such as SWISS-MODEL with user-friendly web interfaces generate reliable models without the need for complex software packages or downloading large databases. Here, we describe the latest version of the SWISS-MODEL expert system for protein structure modelling. The SWISS-MODEL template library provides annotation of quaternary structure and essential ligands and co-factors to allow for building of complete structural models, including their oligomeric structure. The improved SWISS-MODEL pipeline makes extensive use of model quality estimation for selection of the most suitable templates and provides estimates of the expected accuracy of the resulting models. The accuracy of the models generated by SWISS-MODEL is continuously evaluated by the CAMEO system. The new web site allows users to interactively search for templates, cluster them by sequence similarity, structurally compare alternative templates and select the ones to be used for model building. In cases where multiple alternative template structures are available for a protein of interest, a user-guided template selection step allows building models in different functional states. SWISS-MODEL is available at http://swissmodel.expasy.org/.


 


Introduction

In order to understand how specific residues in the sequence contribute to the putative function of the protein, and why and how they are conserved throughout evolution, we would need to study an explicit molecular model of an APSES domain protein, bound to its cognate DNA sequence. Explanations of a protein's observed properties and functions can't rely on the general fact that it binds DNA, we need to consider details in terms of specific residues and their spatial arrangement. In particular, it would be interesting to correlate the conservation patterns of key residues with their potential to make specific DNA binding interactions. Unfortunately, the experimental evidence (Taylor et al., 2000) is not sufficient to unambiguously define the details of how a DNA double helix might be bound. Moreover, several distinct modes of DNA binding are known for proteins of the winged-helix superfamily, of which the APSES domain is a member.

In this assignment you will construct a molecular model of the APSES domain from the Mbp1 RBM orthologue in MYSPE.

For the following, please remember the following terminology:

Target
The protein that you are planning to model.
Template
The protein whose structure you are using as a guide to build the model.
Model
The structure that results from the modelling process. It has the Target sequence and is similar to the Template structure.


 

The basic idea - a Point Mutation

To illustrate how force fields modify protein structure in principle, let's consider changing the sequence of a single amino acid, based on a structural template and minimize the structure's energy.

Such minimal changes to structure models can be done directly in Chimera. Let us consider the residue A 42 of the 1BM8 structure. It is oriented towards the core of the protein, but as the MSA shows, most other Mbp1 orthologs have a larger amino acid in this position: V, or even I.

Task:

  1. Open 1BM8 in Chimera, hide the ribbons and show all protein atoms as a stick model.
  2. Color the protein white.
  3. Open the sequence window and select A 42. Color it red. Choose Actions → Set pivot. Then study how nicely the alanine sidechain fits into the cavity formed by its surrounding residues.
  4. To emphasize this better, select the protein atoms and display them as a sphere model to better appreciate the packing, i.e. the Van der Waals contacts. Use the Favorites → Side view panel to move the clipping plane and see a section through the protein. Study the packing, in particular, note that the additional methyl groups of a valine or isoleucine would not have enough space in the structure. Then restore the clipping planes so you can see the whole molecule.
  5. Lets simplify the view: choose Actions → Atoms/Bonds → show and Actions → Atoms/Bonds → backbone only → chain trace. Then select A 42 again in the sequence window and choose Actions → Atoms/Bonds → show.
  6. Add the surrounding residues: choose Select → Zone.... In the window, see that the box is checked that selects all atoms at a distance of less then 5Å to the current selection, and check the lower box to select the whole residue of any atom that matches the distance cutoff criterion. Click OK and choose Actions → Atoms/Bonds → show. You now have a very clear scene of the alanine residue in red, the surrounding side chains, and the rest of the structure as a C-alpha trace. You also see three water molecules. Spend a bit of time again, to get a sense for the spatial context[1].
  7. Select A 42 again: left-click (control click) on any atom of the alanine to select the atom, then up-arrow to select the entire residue. Now let's mutate this residue to isoleucine.
  8. Choose Tools → Structure Editing → Rotamers and select ILE as the rotamer type. Click OK, a window will pop up that shows you the possible rotamers for isoleucine together with their database-derived probabilities; you can select them in the window and cycle through them with your arrow keys. But note that the probabilities are very different - and thus show you high-energy and low-energy rotamers to choose from. Therefore, unless you have compelling reasons to do otherwise, try to find the highest-probability rotamer that may fit. This is where your stereo viewing practice becomes important, if not essential. It is really, really hard to do this reasonably in a 2D image! It becomes quite obvious in 3D.
  9. I find that the first rotamer is actually not such a bad fit. The CD atom comes close to the sidechains of I 25 and L 96. But we can assume that these are somewhat mobile and can accommodate a denser packing, because - as you can easily verify in your MSA - it is NOT the case that sequences that have I 42, have a smaller residue in position 25 and/or 96. So let's accept the most frequent ILE rotamer by selecting it in the rotamer window and clicking OK (while existing side chain(s): replace is selected).
  10. Done.

If you want to go over this in more detail, check the video tutorial on YouTube published by the NIAID bioinformatics group here. I would also encourage you to go over Part 2 of the video tutorial that discusses how to check for and resolve (by energy minimization) steric clashes. But do remember that it is not clear whether energy minimization will make your structure more correct in the sense of a smaller overall RMSD with the real, mutated protein.

Incidentally: What we have done here with one residue is exactly the way homology modeling works with entire sequences. The homology modelling program simply changes all amino acids to the residues of the target sequence, based on the template structure. Let's now build a homology model for MYSPE Mbp1.


 

Preparation

  • We need to define our Target sequence;
  • find a suitable structural Template; and
  • build a Model.


Target sequence

We have encountered the PDB 1BM8 structure before, the APSES domain of saccharomyces cerevisiae Mbp1. This is a useful template to model the DNA binding domain of your RBM match. You have defined the sequence in the BIN-ALI-Optimal_sequence_alignment unit. Let's retrieve it. Open RStudio and load the project.


 
library(msa)

# Recreate the database
source("makeProteinDB.R")

# A: Define your TARGET sequence.
#      You have defined a feature annotation for the MYSPE APSES domain in
#      the BIN-ALI-Optimal_sequence_alignment unit's R code. Retrieve it's
#      sequence from the feature annotation to get the TARGET sequence.
#

(targetName <- sprintf("MBP1_%s", biCode(MYSPE)))

# Get the protein IDs.
(sel <- which(myDB$protein$name == targetName))
(proID <- myDB$protein$ID[sel])

# Find the feature ID in the feature table
(ftrID <- myDB$feature$ID[myDB$feature$name == "APSES fold"])

# Get the annotation ID.
(fanID <- myDB$annotation$ID[myDB$annotation$proteinID == proID &
                             myDB$annotation$featureID == ftrID])

# Get the feature start and end:
(start <- myDB$annotation$start[fanID])
(end   <- myDB$annotation$end[fanID])

# Extract the feature from the sequence
targetSeq <- substring(myDB$protein$sequence[sel], first = start, last = end)

# Name it
names(targetSeq) <- targetName

targetSeq


 

Template choice and template sequence

The SWISS-MODEL server provides several different options for constructing homology models. The easiest option requires only a target sequence as input. In this mode the program will automatically choose suitable templates and create an input alignment. I think that is not the best way to use such a service: template choice and alignment both may be significantly influenced by biochemical reasoning, and an automated algorithm cannot make the necessary decisions. Should you use a structure of reduced resolution that however has a ligand bound? Should you move an indel from an active site to a loop region even though the sequence similarity score might be less? Questions like that may have answers that are different from the best choices an automated algorithm could make. But Swiss Model is flexible and allows us to upload an explicit alignment between target and template. Please note: the model you will produce is "easy" - the sequence similarity is high and there are no significant indels to consider, the automated mode would have done just as well. But the strategy we pursue here is also suitable for much more difficult problems. The automated strategy maybe not. More control over the process is a good thing.

Template choice is the first step. Often more than one related structure can be found in the PDB. The degree of sequence identity is the most important criterion, but there are many other factors to consider. Please refer to the template choice principles page on this Wiki where I discuss more details and alternatives. To find related structures, you can search the PDB itself through its Advanced Search interface; for example one can search for sequence similarity with a BLAST search, or search for structural similarity by accessing structures according to their CATH or SCOP classification. But the BLAST search is probably the method of choice: after all, the most important measure of the probability of success for homology modelling is sequence similarity.


Defining a template means finding a PDB coordinate set that has sufficient sequence similarity to your target that you can build a model based on that template. To find suitable PDB structures, we will perform a BLAST search at the PDB.




Task:

  1. Navigate to the PDB.
  2. Click on Advanced Search to enter the advanced search interface.
  3. Open the menu to Choose a Query Type:
  4. Find the Sequence features section and choose Sequence (BLAST...)
  5. Copy the targetSequence from the R console and paste it into the Sequence field, select BLAST as the search tool, select not to mask low-complexity regions and Submit Query. Since the E-value is set rather high by default, you will get a number of low-confidence hits as well as the actual homologs, these have very low E-values.

All hits that are homologs are potentially suitable templates, but some are more suitable than others. Consider how the coordinate sets differ and which features would make each more or less suitable for creating a homology model: you should consider ...

  • sequence similarity to your target
  • size of expected model (= length of alignment)
  • presence or absence of ligands
  • experimental method and quality of the data set

Sequence similarity is the most important, but we can have the PDB tabulate the other features concisely for this task.

  1. There is a menu to create Reports: - select customizable table.
  2. Select (at least) the following information items:
Structure Summary
  • Experimental Method
Sequence
  • Chain Length
Ligands
  • Ligand Name
Biological details
  • Macromolecule Name
refinement Details
  • Resolution
  • R Work
  • R free
  1. click: Create report.

Unfortunately you don't get the E-values into the report, and those should strongly influence your final decision. As of October 2017, you should find four reasonable candidate structures from 2 species, three of which are from the same species. Some of the yeast sequences have a longer chain-lengths ... but those are only disordered residues (otherwise these would be better suited templates; regrettably, you'd need to check that in the real world, there is no automatic tool to evaluate disorder and its effects on template choice). Depending on MYSPE, your ideal template will be either be 1BM8 or 4UX5. Let's consider both.

Finally
Click on the ID to navigate to the structure page for those templates and save the FASTA sequences to your project directory. Name one 1BM8_A.fa and the other 4UX5_A.fa (save only chain A for 4UX5). These are template sequence.


 


The input alignment

The sequence alignment between target and template is the single most important factor that determines the quality of your model. No comparative modeling process will repair an incorrect alignment; it is useful to consider a homology model rather like a three-dimensional map of a sequence alignment rather than a structure in its own right. In a homology modeling project, typically the largest amount of time should be spent on preparing the best possible alignment. Even though automated servers like the SwissModel server will align sequences and select template structures for you, it would be unwise to use these just because they are convenient. You should take advantage of the much more sophisticated alignment methods available. Analysis of wrong models can't be expected to produce right results.

The best possible alignment is constructed from a multiple sequence alignment that includes at least the target and template sequence and other related sequences as well. The additional sequences are an important aid in identifying the correct placement of insertions and deletions. Your alignment should have been carefully reviewed by you and wherever required, manually adjusted to move insertions or deletions between target and template out of the secondary structure elements of the template structure.

In most of the Mbp1 orthologues, we do not observe indels in the APSES domain regions, but in some we do. Evolutionary pressure on the APSES domains has selected against indels in the more than 600 million years these sequences have evolved independently in their respective species. To obtain an alignment between the template sequence and the target sequence from your species we fetch the Mbp1 sequences from our database, add the template sequences, and convert them to an AAStringSet.


 

Task:
Here's how we do this in R:

# Get all MBP1 Sequences
sel <- grep("^MBP1_", myDB$protein$name)

# Extract the sequences
MBP1Set <- myDB$protein$sequence[sel]

# Name the sequences
names(MBP1Set) <- myDB$protein$name[sel]

# Read the template sequences
seq1BM8 <- dbSanitizeSequence(readLines("1BM8_A.fa"))
names(seq1BM8) <- "1BM8_A"
seq4UX5 <- dbSanitizeSequence(readLines("4UX5_A.fa"))
names(seq4UX5) <- "4UX5_A"

# Add the template sequences to the MBP1set
MBP1Set <- c(MBP1Set, seq1BM8, seq4UX5)

# Turn it into an AAStringSet
(MBP1Set <- AAStringSet(MBP1Set))   # You should have 13 sequences.

# Calculate an msa
(MBP1msa <- msaMuscle(MBP1Set))

# Inspect the msa
writeALN(fetchMSAmotif(MBP1msa, seq1BM8)) # and ...
writeALN(fetchMSAmotif(MBP1msa, seq4UX5))

You need to decide which of the templates you will use. Choose either 1BM8 or 4UX5 - depending on which template has higher sequence similarity to the target. Next, extract aligned target and template sequences, while masking gaps that are not needed for the aligned pair.

# Write the alignments to file, we will need it later. Depending on which
# template you have decided on, execute ...
writeMFA(fetchMSAmotif(MBP1msa, seq1BM8), myCon = "APSES-MBP1.fa") # or ...
writeMFA(fetchMSAmotif(MBP1msa, seq4UX5), myCon = "APSES-MBP1.fa")

# We extract the TARGET and TEMPLATE sequence, and remove any hyphens that
# they both share. Remember: the TARGET is the MYSPE sequence in this alignment,
# the TEMPLATE is either 1BM8_A or 4UX5_A. You need to edit this code so it
# identifies the correct sequences for your situation:

myT <- seq1BM8 # either ...
myT <- seq4UX5 # ... or .

targetSeq   <- as.character(fetchMSAmotif(MBP1msa, myT)[targetName])
templateSeq <- as.character(fetchMSAmotif(MBP1msa, myT)[names(myT)])

# Drop positions in which both sequences have hyphens.
targetSeq   <- unlist(strsplit(targetSeq,   ""))
templateSeq <- unlist(strsplit(templateSeq, ""))
gapMask <- ! ((targetSeq == "-") & (templateSeq == "-"))
targetSeq   <- paste0(targetSeq[gapMask], collapse = "")
templateSeq <- paste0(templateSeq[gapMask], collapse = "")

# Assemble sequences into a set
TTset <- character()
TTset[1] <- targetSeq
TTset[2] <- templateSeq
names(TTset) <- c(targetName, names(myT))

writeMFA(TTset)  # write output to multi FASTA format


The result should be a two sequence alignment in multi-FASTA format, that was constructed from a number of supporting sequences and that contains your aligned target and template sequence. This is your input alignment for the homology modeling server. For MBP1_CRYNE aligned to 4UX5 the result looks like this:

>MBP1_CRYNE
MGKKVIASGGDNGPNTIYKATYSGVPVYEMVCR-DVAVMRRRSDAYLNATQILKVAGFDKPQRTRVLEREVQKGEHEKVQ
GGYGKYQGTWIPIERGLALAKQYGVEDILRPIIDYVPTSVSPPPAPKHSVAPPSKARRDK
>4UX5_A
MVKAAAAAASAPTGPGIYSATYSGIPVYEYQFGLKEHVMRRRVDDWINATHILKAAGFDKPARTRILEREVQKDQHEKVQ
GGYGKYQGTWIPLEAGEALAHRNNIFDRLRPIFEFSPGPDSPPPAPRH----TSKPKQPK


 

Homology model

The alignment defines the residue by residue relationship between target and template sequence. All we need to do now is to change every residue of the template to the target sequence - that's what the homology modelling server will do.


SwissModel

Access the Swissmodel server at https://swissmodel.expasy.org and click on the Start Modelling button. Under the Supported Inputs, choose Target-Template Alignment.

Task:

  • Paste the aligned sequences of the MYSPE target and the template into the form field. SwissModel will analyse the sequences and ask you to identify target and template. The MYSPE sequence is your target. The 1BM8 or 4UX5 sequence is the template.
  • Click Build Model to start the modeling process. This will take about a minute or so.
  • The resulting page returns information about the resulting model and its quality. You can rotate the model in the window on the right with the mouse. Regions that have a reddish hue have lower quality scores, i.e. they were harder to model or could not be modelled well with good geometry. Hovering the mouse over parts of the structure highlights the respective region of the sequence alignment.
  • Mouse over the Model 01 dropdown menu (under the icon of the template structure), and choose the PDB file. Note that the B-factor column of the coordinate section contains the QMEAN scores (between 0 and 1) that the server has calculated. Higher is better. Save the PDB file in your project directory call it MBP1_MYSPE-APSES.pdb.
  • Open the SwissModel documentation in a new tab. Read about the modelling process. there are a number of important technical details that help to understand what the computed coordinates of your model mean, you should pay special attention to the GQME and QMEAN quality scores.
  • Also save:
    • The output page as pdf (for reference)
    • The modeling report (as pdf)


Model interpretation

We have spent a significant amount of time to prepare data for the analysis and in practice it usually seems to turn out that way, that the preparation of data occupies the greatest part of our efforts. The actual computational analysis is generally quite fast. And, unfortunately, the interpretation of results is often somewhat neglected. Don't be that way. Data does not explain itself. The interpretation of your computational results in a bilological context is the most important part.


 


The PDB file

Task:
Open your model coordinates PDB file in RStudio (which is an excellent plain-text editor) and consider the following questions:

  • What is the residue number of the first residue in the model? What should it be, based on the alignment? If you read about a sequence number such as "residue 45" in a manuscript, which residues of your model correspond to that number?

That's not easy to tell. But it should be.


Renumbering the model

As you can see from the coordinate file, SwissModel numbers the first residue "1" in the 1BM8-derived structure, and 14 in the 4UX5 structure: it does not keep the numbering of the template. We should renumber the model so we can compare the model and the template with the same residue numbers and thus interpret our model with reference to sequence numbers we find in the manuscript describing the template structure. (An alternative renumbering would renumber the model correspond to the sequence it came from. Remember that we have only excised a domain from the full-length sequence.) Carefully doing this by hand will take you a bit less than an hour. Fortunately we can do this with bio3d.

Task:

  1. Explore and execute the following R script. It assumes that your model is in your is project directory and the file is called MBP1_MYSPE-APSES.pdb.
if (! require(bio3d, quietly=TRUE)) {
  install.packages("bio3d")
  library(bio3d)
}
# Package information:
#  library(help = bio3d)       # basic information
#  browseVignettes("bio3d")    # available vignettes
#  data(package = "bio3d")     # available datasets

PDB_INFILE      <- "MBP1_MYSPE-APSES.pdb"
PDB_OUTFILE     <- "MBP1_MYSPE-APSESrenum.pdb"


iFirst <-  4  # residue number for the first residue if your template was 1BM8
iFirst <- 14  # residue number for the first residue if your template was 4UX5


# == Read the MYSPE pdb file
MYSPEmodel <- read.pdb(PDB_INFILE) # read the PDB file into a list

MYSPEmodel           # examine the information
MYSPEmodel$atom[1,]  # get information for the first atom

# Explore ?read.pdb and study the examples.

# == Modify residue numbers for each atom
resNum <- as.numeric(MYSPEmodel $atom[,"resno"])
resNum
resNum <- resNum - resNum[1] + iFirst  # add offset
MYSPEmodel $atom[ , "resno"] <- resNum   # replace old numbers with new

# check result
MYSPEmodel $atom[ , "resno"]
MYSPEmodel $atom[1, ]

# == Write output to file
write.pdb(pdb = MYSPEmodel, file=PDBout)

# Done. Open the renumbered PDB file in the RStudio editor
# and confirm that this has worked.


 

First visualization - colouring the model by energy

 


SwissModel calculates energies for each residue of the model with a molecular mechanics forcefield. The SwissModel modeling summary page contains a plot of these energies as a function of sequence number like. The values - between 0.0 and 1.0 - are stored in the PDB file's B-factor field.

Task:

  1. Start Chimera and load the model coordinates that you have just renumbered.
  2. Select all, hide Ribbons and show Atoms, bonds to view the entire model structure.
  3. Choose Tools → Depiction → Render by attribute and select attributes of atoms, Attribute: bfactor, check color atoms and click Apply. Note that you can change the way the spectrum is mapped to the values by moving the blue, white and red bars over the histogram with your mouse.
  4. Study the result: It seems that residues in the core of the protein have better energies (higher values) than residues at the surface, i.e. Swiss-Model was more confident in the predicted conformationstes. Why could that be the case?

Study the options of this window a bit, rendering by attribute is a powerful way to store and depict all manners of information with the molecule. You can simply write a little R script that uses bio3D to replace the B-factor or occupancy values with any value you might be interested in: energies, conservation scores, information ... whatever. Then render this property to map it on the 3D structure of your molecule...


 

Modelling DNA binding

One of the really interesting questions we can discuss with reference to our homology model is how sequence variation might result in changed DNA recognition sites, and then lead to changed cognate DNA binding sequences. In order to address this, we would need to generate a plausible structural model for how DNA is bound to APSES domains.

Since there is currently no software available that would reliably model such a complex from first principles[2], we will base a model of a bound complex on homology modelling as well. This means we need to find a similar structure for which the position of bound DNA is known, then superimpose that structure with our model. This places the DNA molecule into the spatial context of the model we are studying. As a result of the PDB BLAST search we found 4UX5, from the Magnaporthe oryzae Mbp1 orhologue PCG2[3]: this is a protein-DNA complex structure.


 

A homologous protein/DNA complex structure

Task:

The PCG2 / DNA complex
  • Open Chimera.
  • load the 4UX5 structure. Spend some time exploring it. There are two domains of identical sequence bound to one DNA molecule.
    • If your homology model was based on 4UX5, Swiss-Model has already made two copies, and their orientation is the same as the template, so no superposition is required.
    • If your homology model was based on 1BM8: make a second copy of your model. Open the ToolsGeneralModel Panel and use the copy/combine button to create a copy of your model. Then superimpose one copy on chain A of 4UX5, and the other copy on chain B: open a MatchMaker dialogue window with ToolsStructure comparisonMatchMaker. Choose the radio button two match two specific chains and select 4UX5 chain A as the Reference chain, and one of your models as the Chain to match. Click Apply. Similarly superimpose the other copy of the model on chain B.
  • Color the 4UX5 protein chains grey.
  • Color the 4UX5 nucleic acid chains "by element", hide ribbons, show Atoms/Bonds and set nucleotide objects offf.
  • Now color your model by conservation score:
    • In the Multalign Viewer window choose PreferencesHeaders, and in the Headers window choose the Headers tab and select Conservation styleAL2CO[4]. Click OK.
    • In the Multalign Viewer window choose StructureRender by Conservation to open the "Render/Select by Attribute" Window. Select your Model. Select mavConservation as the "Attribute" to render. Note that you can move the blue white and red coloured bars to adjust the way the colour scale is applied to the values. Click on the blue, white and red bar in turn and then on the colour swatch to change the colour. Choose a bright orange red for the low value threshold (high diversity), a dark red for the midpoint, and a dark greay for the high conservation values. Click on Apply. Are all residues that make protein-DNA interactions in the complex conserved between target and template? Are they conserved across the entire family?
  • Do the two molecules bind to the same DNA motif - the CGCG core of the "MCB-box"? Do the chains have protein:DNA interfaces with the cognate sequence, or are one (or both) proteins non-specific complexes? The conditions under which proteins crystallize can be harsh, and physiological function under these conditions is not guaranteed.[5] Indeed, Liu et al. (2015) report that at low concentrations a 1:1 complex is formed and the 2:1 Protein:DNA complex only forms at high concentrations. Figure 3. of their paper shows that the detailed contacts between protein and DNA are in fact not identical.
  • Select one of the residues of that loop in chain A by <control>-clicking on it and use ActionSet pivot to set the centre of rotation to that residue: this makes it easier to visualize the binding situation when you make the molecules larger.
  • Study the situation. Focus on Gly 84.A, especially the interaction of its carbonyl oxygen, which hydrogen bonds to the N2 atom of G8.D chain. Gln 89.A hydrogen bonds to the N2 atom of G8.C chain. Gly 84 and Gln 82 thus recognize a G:C C:G pair. In the B chain, Gly 84.B does not contact the DNA well, since it contacts residues of chain A, especially Gln 82.A. The carbonyl atom of Gly 84.B hydrogen bonds to Gln 89.B. and therefore Gln89.B is not available to contact nucleotide bases. What do you think?


 

In summary: superimposing our homology model with a protein:DNA complex has allowed us to consider how our target sequence might perform its function. This is supported by considering variations in structure between chain A and B of the protein DNA complex that may point to different binding modes, and it is further supported by being able to map structural conservation onto our model, to understand which residues play a structural or functional role that is shared within the entire family.


 


 


Further reading, links and resources

 


Notes

  1. Chimera uses a default distance to screen that is too close and that exaggerates the depth of the scene to a degree that it is difficult to fuse the stereo pairs. Choose Tools → Viewing controls → Camera and set the distance to screen to 50 cm. This will make stereo viewing easier and will also give a better sense of distance estimates in all three dimensions.
  2. Rosetta may get the structure approximately right, Autodock may get the complex approximately right, but the coordinate changes involved in induced fit makes the result unreliable - and we have no good way to validate whether the predicted complex is correct.
  3. Liu et al. (2015) Structural basis of DNA recognition by PCG2 reveals a novel DNA binding mode for winged helix-turn-helix domains. Nucleic Acids Res 43:1231-40. (pmid: 25550425)

    PubMed ] [ DOI ] The MBP1 family proteins are the DNA binding subunits of MBF cell-cycle transcription factor complexes and contain an N terminal winged helix-turn-helix (wHTH) DNA binding domain (DBD). Although the DNA binding mechanism of MBP1 from Saccharomyces cerevisiae has been extensively studied, the structural framework and the DNA binding mode of other MBP1 family proteins remains to be disclosed. Here, we determined the crystal structure of the DBD of PCG2, the Magnaporthe oryzae orthologue of MBP1, bound to MCB-DNA. The structure revealed that the wing, the 20-loop, helix A and helix B in PCG2-DBD are important elements for DNA binding. Unlike previously characterized wHTH proteins, PCG2-DBD utilizes the wing and helix-B to bind the minor groove and the major groove of the MCB-DNA whilst the 20-loop and helix A interact non-specifically with DNA. Notably, two glutamines Q89 and Q82 within the wing were found to recognize the MCB core CGCG sequence through making hydrogen bond interactions. Further in vitro assays confirmed essential roles of Q89 and Q82 in the DNA binding. These data together indicate that the MBP1 homologue PCG2 employs an unusual mode of binding to target DNA and demonstrate the versatility of wHTH domains.

  4. Pei & Grishin (2001) AL2CO: calculation of positional conservation in a protein sequence alignment. Bioinformatics 17:700-12. (pmid: 11524371)

    PubMed ] [ DOI ] MOTIVATION: Amino acid sequence alignments are widely used in the analysis of protein structure, function and evolutionary relationships. Proteins within a superfamily usually share the same fold and possess related functions. These structural and functional constraints are reflected in the alignment conservation patterns. Positions of functional and/or structural importance tend to be more conserved. Conserved positions are usually clustered in distinct motifs surrounded by sequence segments of low conservation. Poorly conserved regions might also arise from the imperfections in multiple alignment algorithms and thus indicate possible alignment errors. Quantification of conservation by attributing a conservation index to each aligned position makes motif detection more convenient. Mapping these conservation indices onto a protein spatial structure helps to visualize spatial conservation features of the molecule and to predict functionally and/or structurally important sites. Analysis of conservation indices could be a useful tool in detection of potentially misaligned regions and will aid in improvement of multiple alignments. RESULTS: We developed a program to calculate a conservation index at each position in a multiple sequence alignment using several methods. Namely, amino acid frequencies at each position are estimated and the conservation index is calculated from these frequencies. We utilize both unweighted frequencies and frequencies weighted using two different strategies. Three conceptually different approaches (entropy-based, variance-based and matrix score-based) are implemented in the algorithm to define the conservation index. Calculating conservation indices for 35522 positions in 284 alignments from SMART database we demonstrate that different methods result in highly correlated (correlation coefficient more than 0.85) conservation indices. Conservation indices show statistically significant correlation between sequentially adjacent positions i and i + j, where j < 13, and averaging of the indices over the window of three positions is optimal for motif detection. Positions with gaps display substantially lower conservation properties. We compare conservation properties of the SMART alignments or FSSP structural alignments to those of the ClustalW alignments. The results suggest that conservation indices should be a valuable tool of alignment quality assessment and might be used as an objective function for refinement of multiple alignments. AVAILABILITY: The C code of the AL2CO program and its pre-compiled versions for several platforms as well as the details of the analysis are freely available at ftp://iole.swmed.edu/pub/al2co/.

  5. This particular crystal structure however was crystallized from a Tris-buffer with 50mM NaCl at pH 8.0 - comparatively gentle conditions actually.


 


Self-evaluation

 



 




 

If in doubt, ask! If anything about this learning unit is not clear to you, do not proceed blindly but ask for clarification. Post your question on the course mailing list: others are likely to have similar problems. Or send an email to your instructor.



 

About ...
 
Author:

Boris Steipe <boris.steipe@utoronto.ca>

Created:

2017-08-05

Modified:

2017-10-30

Version:

1.0

Version history:

  • 1.0 First live version
  • 0.1 First stub

CreativeCommonsBy.png This copyrighted material is licensed under a Creative Commons Attribution 4.0 International License. Follow the link to learn more.