Difference between revisions of "Werner Syndrome"

From Bioinformatikpedia
Line 1,027: Line 1,027:
 
This was probably the first experimental evidence of a dynamical protein world as predicted by Feynman. Around the decade of 1970, biologists were already familiar with protein sequencing techniques and experimental determination of protein structures. The first protein to be sequenced was the bovine insulin, in 1951 by Frederick Sanger, while the first protein structures were determined by X-ray crystallography less than 10 years after that. These were the first ingredients for the development of an in silico approach for observation of protein dynamics. The first use of MD simulations was described 1957 through the work entitled Phase Transitions for a Hard Sphere System, by B.J. Alder and T. E. Wainwright, who were studying the
 
This was probably the first experimental evidence of a dynamical protein world as predicted by Feynman. Around the decade of 1970, biologists were already familiar with protein sequencing techniques and experimental determination of protein structures. The first protein to be sequenced was the bovine insulin, in 1951 by Frederick Sanger, while the first protein structures were determined by X-ray crystallography less than 10 years after that. These were the first ingredients for the development of an in silico approach for observation of protein dynamics. The first use of MD simulations was described 1957 through the work entitled Phase Transitions for a Hard Sphere System, by B.J. Alder and T. E. Wainwright, who were studying the
 
interactions between solid spheres. Subsequently to the application of the new developed technique in this simple system, there was a gap of almost 20 years before its application to a representation of a real system. In 1974, the MD simulation of a liquid water system was performed, opening new possibilities to the use of MD in biology. Today, MD simulations serve as a powerful tool for elucidating the atomic-level behavior of proteins. Since McCammon published the results of the first MD simulation of a protein, the computational environment available for such calculations has become more and more unlimited. As mentioned, the first simulation of BPTI was only approximately 10 ps long. As addressed by Karplus and McCammon, during the years that followed this first simulation, a wide range of motional phenomena were investigated through MD simulations of proteins. Last year, the same system (the BPTI protein) was the subject of a 100000000 times longer simulation, that could only be accomplished with the use of a specific-purpose supercomputer. Other long simulations were performed in the past years, with special note to a 10 ms long simulation of a fast-folding protein
 
interactions between solid spheres. Subsequently to the application of the new developed technique in this simple system, there was a gap of almost 20 years before its application to a representation of a real system. In 1974, the MD simulation of a liquid water system was performed, opening new possibilities to the use of MD in biology. Today, MD simulations serve as a powerful tool for elucidating the atomic-level behavior of proteins. Since McCammon published the results of the first MD simulation of a protein, the computational environment available for such calculations has become more and more unlimited. As mentioned, the first simulation of BPTI was only approximately 10 ps long. As addressed by Karplus and McCammon, during the years that followed this first simulation, a wide range of motional phenomena were investigated through MD simulations of proteins. Last year, the same system (the BPTI protein) was the subject of a 100000000 times longer simulation, that could only be accomplished with the use of a specific-purpose supercomputer. Other long simulations were performed in the past years, with special note to a 10 ms long simulation of a fast-folding protein
domain, that needed three months to be finished. Although, in general the average length of a simulation ranges around 10 ns nowadayså.
+
domain, that needed three months to be finished. Although, in general the average length of a simulation ranges around 10 ns nowadays.
  +
  +
  +
  +
<font size=4>
  +
=== Molecular Dynamics Pipeline ===
  +
  +
<font size=3>
  +
Molecular dynamics (MD) simulations, energy minimization and trajectory analyses were carried out with GROMACS package [51], using GROMOS96 (G53a6) force field [51]. To fully automatize the protocol, a Perl script was designed based on the available version by Marc Offman. The script has several steps, mainly including: 1) Preliminar solvent treatment, where only the so-called structural water molecules (those with conserved positions within the molecules from the crystal and presumably involved in protein function) are maintained in the system; 2) Coordinates file edition, to avoid problems with residue nomenclature and numbering, chain separation, etc; 3) Side-chain treatment, performed by SCWRL to correct andor refine the position of side-chain atoms; 4) The inclusion of solvent and ion molecules (Cl- and Na+) in the system; 5) The minimization of the solvent, with a fix protein structure; 6) The minimization of the protein side-chain, with a fix backbonr; 7) Equilibration phases consisting of short molecular dynamics simulations under constant volume and constant pressure; 8) A production run (usually varying from 2ns to 20ns in length), which is the molecular dynamics simulation itself with all atoms free to move and no volume or pressure constrains.
  +
Explicit SPC water molecules [54] were used in all simulations, in which a 14 Å layer of water molecules were added around the solute molecules, within a cubic water box, using periodic boundary conditions. Counter ions were inserted for system neutralization. LINCS [55] and SETTLE [56] were applied to constraint solute and solvent bond, respectively. Temperature was kept at 298 K by rescaling velocities with a stochastic term [57] and pressure at 1 atm using the Berendsen approach [58]. Electrostatic interactions were corrected with PME method [59], using non-bonded cutoffs of 1.0 nm for Coulomb and 1.2 nm for van der Waals. MD integration time was 2 fs.
  +
  +
A 3-steps energy minimization protocol was used to avoid artifacts in atomic trajectories due to conversion of potential into kinetic energies: firstly, applying the steepest-descent algorithm: i. 5000 steps with solute heavy atom positions restrained to their initial positions using an harmonic constant of 1 kJ/mol.nm in each Cartesian direction, allowing free water and hydrogen movements; and ii. 5000 steps with all atoms free to move. Subsequently, the conjugated gradient algorithm was applied for further energy minimization until an energy gradient of 42 KJ/mol.nm. A preliminary MD (1 ns), with heavy atom positions restrained, was performed for achieving solvent equilibration and system heating until 298 K. In this step, the initial velocities were generated once for each simulation. Then we performed a 10 ns production MD.
  +
  +
  +
  +
<font size=4>
  +
=== References ===
  +
  +
<font size=2>
  +
51. Oostenbrink C, Villa A, Mark AE, van Gunsteren WF: A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force-field parameter sets 53A5 and 53A6. Journal of computational chemistry 2004, 25:1656-1676

Revision as of 21:06, 26 August 2011

ABOUT THIS WIKI

This Wiki Page was created and it is maintained by Mainá Bitar for the "Protein Structure and Function Analysis Practical" of 2011 (SS).


WERNER SYNDROME

You can Download a PDF version of this information here: Media:WS.pdf.

Introduction

The Werner Syndrome (WS) is an autosomal recessive disorder, also known as Adult Progeria. The syndrome was described for the first time in 1904 by Otto Werner (and therefore, named after him), in his PhD thesis entitled “Über katarakt in Verbindung mit Sklerodermie” (which can be translated to “About cataracts connected to sclerodermia”). In the first 90 years of research concerning WS, over 1000 patients were reported, 75% of which were of Japanese descent (figure 1) [1]. WS is one of the several types of segmental progeroid syndromes, which affect multiple tissues and organs (on the other hand, unimodal syndromes predominantly affect a single organ) [2].

Figure 1: Distribution of Werner Syndrome by nationality as registered from 1904 until 1994 [1]. All countries with at least one patient is shaded.


<center|font size=2>Figure 1: Distribution of Werner Syndrome by nationality as registered from 1904 until 1994 [1]. All countries with at least one patient are shaded.


Symptoms

As one can expect, the most notable symptoms of WS mimic the background of the most general condition called Progeria, with a complex phenotype of accelerated aging. The patients prematurely acquire the appearance of someone several decades older, accompanied by loss or graying of hair, scleroderma-like skin and voice alterations, usually around the second or third decade of life [3]. The phenotype of WS was previously summarized as a “caricature of aging” [1]. In general, the subjects develop normally until adolescence, when there is absence of the common growth spurt. The clinical manifestations usually include atherosclerosis, osteoporosis, diabetes, lenticular cataracts, heart failure, cancer and other age-related conditions that appear during early adulthood, following puberty (figure 2) [1, 2]. The typical cause of death is cancer or cardiovascular disease, often occurring between the fourth and fifth decades of life. While in 1966, the median age of death was 47 years [4], in 1997, Makoto Goto reported a surprising median age of death of 54 years [2]. At the cellular level, a reduction in the replicative rate is often observed (cellular senescence) and genomic instability is present in the form of chromosome breaks and translocations, as well as large deletions at the molecular level. A higher occurrence of somatic mutations is also related to the syndrome [1-4].

Figure 2: Appearance of symptoms in Werner Syndrome patients. The horizontal axis represents the average age at which each clinical manifestation occurs, according to Makoto Goto in his study from 1997 [1].


<center|font size=2>Figure 2: Appearance of symptoms in Werner Syndrome patients. The horizontal axis represents the average age at which each clinical manifestation occurs, according to Makoto Goto in his study from 1997 [1].


Related Gene and Protein

The gene related to WS (WRN) was identified as located in the chromosome 8p12 and cloned for the first time in 1996 [5] and pointed as homologous to members of the DNA helicase family. Later on, the protein was shown to catalyze DNA unwinding, as expected [5]. The protein was predicted to have 1432 amino acids and 35 exons, from which 34 are protein coding. At first, four mutations at this gene were identified in WS patients, all correlated to truncated proteins of no more than 1,250 amino acids. The mutated proteins were shorter, but did not present any effects on the helicase domain itself, which is situated between amino acids 569 and 859. Five additional mutations were identified shortly after, being two nonsense mutations, one mutation at a splice-junction site and a deletion leading to a frameshift [3]. The majority of these mutations directly affect the helicase domain, two of those are located within the domain and other two lead to its loss with truncated proteins.

The WRN gene can be divided into three distinct regions. The N-terminal portion, comprising codons 1-539 is mainly acid and first it was described as presenting no homology to known genes. Currently, this region is known to contain an exonuclease domain, an unusual feature among members of this protein family [6]. Similar high concentration of acidic residues are also observed in other DNA repair-deficiency disorders. The median portion of the gene, from codon 540 to 963 is closely related to other helicases from several organisms, presenting all seven conserved motifs that characterize the protein. The C-terminal end of the WRN gene enclosures a nuclear localization signal (NLS) [7]. Two additional domains are found between the helicase domain and the NLS, namely a RecQ helicase conserved region (RQC) and a helicase RNaseD C-terminal conserved region (HDRC, figure 3) [2]. The WRN protein is likely to act on DNA repair, recombination and replication as well as in the maintenance of telomeres.

Figure 3: Basic structure of the WRN gene, with known domains highlighted according to description made by Friedrich et al. [2].


<center|font size=2>Figure 3: Basic structure of the WRN gene, with known domains highlighted according to description made by Friedrich et al. [2].


The majority of the 60 first described mutations associated to the occurrence of WS are related to the loss of the NLS, causing the protein to accumulate outside the nucleus and therefore be incapable of performing its function [2]. Nevertheless, in the past year, 18 new mutations were reported, yielding a total of 71 WRN disease associated mutations identified in clinically diagnosed patients (table 1) [2].

Table 1: Mutations of the WRN gene associated to the Werner Syndrome, according to the mutation report by Friedrick et al., 2010. The total number of known mutations by class is given [2].


<center|font size=2>Table 1: Mutations of the WRN gene associated to the Werner Syndrome, according to the mutation report by Friedrick et al., 2010. The total number of known mutations by class is given [2].


References

1. Goto M (1997). Hierarchical deterioration of body systems in Werner’s syndrome: Implications for normal ageing. Mechanisms of Ageing and Development, 98:239–254

2. Friedrich K, Lee L, Leistritz DF, Nurnberg G, Saha B, Hisama FM, Eyman DK, Lessel D, Nurnberg P, Li C, Garcia FVMJ, Kets CM, Schmidtke J, Cruz VT, Van den Akker PC, Boak J, Peter D, Compoginis G, Cefle K, Ozturk S, Lopez N, Wessel T, Poot M, Ippel PF, Groff-Kellermann B, Hoehn H, Martin GM, Kubisch C and Oshima J (2010). WRN mutations in Werner syndrome patients: genomic rearrangements, unusual intronic mutations and ethnic-specific alterations. Human Genetics, 128:103-11

3. Yu CE, Oshima J, Wijsman EM, Nakura J, Miki T, Piussan C, Matthews S, Fu YH, Mulligan J, Martin GM, Schellenberg GD and the Werner's Syndrome Collaborative Group (1997). Mutations in the consensus helicase domains of the Werner Syndrome gene. The American Society of Human Genetics, 60:330-341

4. Epstein CJ, Martin GM, Schultz AL, Motulsky AG (1966). Werner's syndrome: a review of its symptomatology, natural history, pathologic features, genetics and relationship to the natural aging process. Medicine, 45:177-221

5. Yu CE, Oshima J, Fu YH, Wijsman EM, Hisama F, Alisch R, Matthews S, Nakura J, Miki T, Ouais S, Martin GM, Mulligan J, Schellenberg GD (1996). Positional cloning of the Werner’s syndrome gene. Science, 272:258–262

6. Huang S, Li B, Gray MD, Oshima J, Mian IS, Campisi J (1998) The premature ageing syndrome protein, WRN, is a 3' → 5' exonuclease. Nature Genetics, 20:114–116

7. Suzuki T, Shiratori M, Furuichi Y, Matsumoto T (2001). Diverged nuclear localization of Werner helicase in human and mouse cells. Oncogene, 20:2551-2558

A. Huang S, Lee L, Hanson NB, Lenaerts C, Hoehn H, Poot M, Rubin CD, Chen DF, Yang CC, Juch H, Dorn T, Spiegel R, Oral EA, Abid M, Battisti C, Lucci-Cordisco E, Neri G, Steed EH, Kidd A, Isley W, Showalter D, Vittone JL, Konstantinow A, Ring J, Meyer P, Wenger SL, von Herbay A, Wollina U, Schuelke M, Huizenga CR, Leistritz DF, Martin GM, Mian IS and Oshima J (2006). The spectrum of WRN mutations in Werner syndrome patients. Human Mutation, 27:558-567


SEQUENCE ALIGNMENTS

17.05


About Sequence Alignments

For the Structural Biologist involved in the study of proteins and interested in structure-function relationships or comparative modeling, an accurate alignment is as important as a significant similarity between the sequences [8]. This justifies that, while studying protein structures, one should be careful and dedicate efforts to the generation and analysis of such alignments. In this context, a sequence alignment can be defined as a method of comparison that takes into account typical conservations, substitutions, insertions and deletions of amino acids and represents those in a graphical format. In general, amino acids substitutions not related to changes in physico-chemical properties are considered as conservatives of protein function and structure. On the other hand, when a substitution leads to physico-chemical changes it is frequently non-conservative (although exceptions can be observed to both rules). For this reason, when evaluating the conservation of the structure and function of a protein from its amino acids substitutions, one should be capable of evaluating the relation between this residues. This can be accomplished, for example, by grouping amino acids together according to their properties in a Venn diagram (figure 5). In this context, an optimal alignment can be defined as the one with highest score, given a certain scoring system that considers these characteristics.


Figure 5: Venn diagram for the 20 most common amino acids. In this diagram it is possible to graphically explicit the structural and physico-chemical relationships among the residues, for example.


<center|font size=2>Figure 5: Venn diagram for the 20 most common amino acids. In this diagram it is possible to graphically explicit the structural and physico-chemical relationships among the residues, for example.


The first computer programs capable of computing biological comparison of sequences in a reasonable amount of time were originated more than 25 years ago. From this point on, the process of inferring homology from sequence similarity became routine and considerably more reliable. However, despite the fact that sequence alignment is a well established method, homology inference from sequence similarity can be rather controversial. For this reason, to infeer homology between two proteins, one should not only consider the similarity between those, but also evaluate how likely is that the proteins are actually biologically related [8]. Usually, a 30% sequence similarity is considered an accurate cut-off above which two proteins can be considered possible homologues. The second step then is to conclude that this similarity is not a result of a convergence from evolutionary unrelated origins. To guide this evaluation one can take advantage of specific features of related sequences that can be recognized by mathematical models and well distinguished from what is observed in randomly generated sequences relationships [9]. In 2005, Pearson and Sterk [8] presented important conclusions that include: i) the scores attributed to alignments of non-related sequences are not distinguishable from scores of randomly generated sequences; ii) if the attributed score for similar sequences cannot be reproduced with a set of randomly generated sequences, then those must be evolutionary related; iii) sequences that share statistically significant similarity are considered homologues. In addition to those general rules, when proteins share an uncommonly high degree of structural similarity, it is suggested that they are non-recognized homologues.

One can therefore conclude that high quality alignments are crucial to the analysis of structure and function of a protein and to the search of homologues and related sequences. From this statement we should not exclude the mentioned cases where the protein structures are more conserved than the related sequences.


References

8. Pearson WR and Sierk ML (2005). The limits of protein sequence comparison? Current opinion in structural biology, 15(3):254–260

9. Pearson WR (1998). Empirical statistical estimates for sequence similarity searches. Journal of molecular biology, 276(1):71–84


Methodology

Reference Sequence: To acquire a reference sequence, a search by keyword was performed at the National Center for Biotechnology Information (NCBI [1]) using the term "WRN Protein" and restricting the organisms group to Homo sapiens. The retrieved sequence (figure 6) is 1432 amino acids long and is registered under the Gi number 3719421. The sequence is as follows:


Figure 6: Amino acid sequence of the human WRN protein (Gi 3719421) as retrieved from NCBI.


<center|font size=2>Figure 6: Amino acid sequence of the human WRN protein (Gi 3719421) as retrieved from NCBI.


For this work the WRN protein was used as Query in a BLAST (Basic Local Alignment Search Tool [2]) search against the Database of non-related sequences (nr). This search retrieved 116 related proteins with e-value score between 0 and 3e-91. The algorithm was able to recognize different conserved domains within the protein sequence, as shown below (figure 7).


Figure 7: Conserved domains identified in the sequence of the WRN protein.


<center|font size=2>Figure 7: Conserved domains identified in the sequence of the WRN protein. <center|font size=2>Obs: Compare this result with the gene structure of the protein, presented in Figure 3.


Among the sequences retrieved, the first 20 will be further used to generate sequence alignments and help describe basic features of WRN proteins and related sequences. These 20 best hits are listed in the figure below (figure 8).


Figure 8: The 20 best results for the BLAST search performed for the human WRN protein.


<center|font size=2>Figure 8: The 20 best results for the BLAST search performed for the human WRN protein.


From these 20 sequences, a curated input file was generated for further multiple sequence alignments. In this file, each representative specie has only one related WRN protein sequence (except for Homo sapiens, since different isoforms are represented). A long insertion was observed in the very beginning of the sequence from Canis familiaris. This region was deleted from the input file, to facilitate the visualization of the alignments (the original sequence was saved for possible further analysis of the insertion). The final version of the input file has 16 sequences, namely 5 from Homo sapiens, one synthetic protein and 10 proteins from different mammalian species. The first feature to be noticed is that the proteins are highly conserved. This is expected from the results of the BLAST search, since the proteins were characterized by a high percentage of identity/similarity and query coverage. Highly conserved proteins are often involved in crucial functions in the cell and therefore can bare only subtle changes in their sequence/structure throughout the course of evolution. This is probably the case of WRN protein, since its function is related to DNA maintainance and replication.

Multiple sequence alignments were performed by 4 different tools: ClustalW, T-Coffee, Promals and Cobalt. Each of those alignments were subsequently analyzed with DNATagger [10], a web tool for Colorizing alignments according to codon properties. In this case, the color scheme takes into account physico-chemical properties of the residues, as listed by Lehninger.

Blue (Positive amino acids) - Histidine, Lysine and Arginine (depending on pH).

Red (Negative amino acids) - Aspartic acid and Glutamic acid (depending on pH).

Green (Polar but non-charged amino acids) - Cysteine, Glycine, Asparagine, Glutamine, Serine, Threonine and Tyrosine.

Yellow (Non-polar but non-charged amino acids) - Alanine, Phenylalanine, Isoleucine, Leucine, Methionine, Proline, Valine and Triptophan.



In the next figure is the alignment performed by Promals (figure 9). And the other results for multiple sequence alignments (colored by DNATagger as explained above) can be downloaded in a PDF format as listed below.

Media:Promals.pdf

Media:Clustal.pdf

Media:TCoffee.pdf



Figure 9: The first part of the multiple sequence alignment generated by Promals for the curated set of WRN protein homologues.


<center|font size=2>Figure 9: The first part of the multiple sequence alignment generated by Promals for the curated set of WRN protein homologues.


As each method has its own algorithm and structure from which it generates the alignments, one can expect that those will run in different time frames. Indeed, while submitting the same input file to the different programs, some were faster than others. The summary of time needed to generate the results is presented below (figure 10).

Figure 10: Summary of time needed to generate multiple sequence alignments with each method, given the same input file.


<center|font size=2>Figure 10: Summary of time needed to generate multiple sequence alignments with each method, given the same input file.


SEQUENCE-BASED ANALYSIS

24.05



Unless otherwise stated, all the predictions completed in this task were performed with only one input file, the protein sequence in FASTA format [3].

Secondary Structure Preditcion

In the protein world, the term secondary structure refers to stable local sub-structures [10]. The basic tridimensional patterns that form the overall protein structure. The most frequent and relevant secondary structure elements are the alpha helix and the beta strand (figure 11). These structural patterns arise from the formation of hydrogen bonds, connecting the peptide groups in the main-chain. Other protein regions are ordered, but do not form any regular secondary structure. Those regions are often called loops and can have important roles in protein function.


Figure 11: Structural representations of the main secondary structure elements. An alpha helix (in the bottom) and a beta strand (in the top) coordinating a zinc ion in a so-called zinc finger structure.


<center|font size=2>Figure 11: Structural representations of the main secondary structure elements. An alpha helix (in the bottom) and a beta strand (in the top) coordinating a zinc ion in a so-called zinc finger structure.


Several Bioinformatics tools are currently available to predict such secondary structure elements. One of the most used and accurate is Psipred [11], which was used to infer the secondary structure of the WRN protein. PsiPred is available within an online version [4] together with other sequence-based analysis methods. The full result from PsiPred can be downloaded as a PDF (below). Another secondary structure prediction tool is Jpred [12], which is also available as an online server [5]. Unfortunately, it is not possible to perform a prediction for proteins longer than 800 amino acids. As the WRN protein has more than 1400 residues, the sequence was broken (nearly on half) and two predictions were made.

Media:PsiPred.pdf Media:JpredFirstHalf.pdf Media:JpredSecondHalf.pdf



Transmembrane Helix Prediction

In the same online server that hosts PsiPred, one can run predictions of Transmembrane Helix regions with MEMSAT, a widely used transmembrane topology prediction method [13]. There are two different versions of this method, namely MEMSAT 3 and MEMSAT SVM (that uses supported vector machines). As discussed before, the WRN protein is nuclear and therefore it should not present a transmembrane region. Indeed, this method points to the presence of one possible transmembrane helix and both versions (3 and SVM) of this algorithm identify this helix in the same region (around the portion between residues 550 and 600). The difference would be the orientation of the protein (figure 12). While MEMSAT 3 classifies the N-terminal region as cytoplasmic, MEMSAT SVM suggest the exact opposite, situating the protein 'upside-down' with the N-terminal in the extracellular compartment. The inconsistency of the results obtained point to a non-transmembranar protein. As should be the case for the WRN protein.

Figure 12: Summary of MEMSAT prediction results. The version (3 or SVM) is indicated and the orientation of the molecule is taken from the indication of the N/C terminal regions. The yellow box represents the predicted transmembrane helix and the numbers inside the box correspond to the first and last residue of the helix.


<center|font size=2>Figure 12: Summary of MEMSAT prediction results. The version (3 or SVM) is indicated and the orientation of the molecule is taken from the indication of the N/C terminal regions. The yellow box represents the predicted transmembrane helix and the numbers inside the box correspond to the first and last residue of the helix.


From the Stockholm Bioinformatics Center there is another method for transmembrane helix prediction. It is available within an online server [6] in two versions, OCTOPUS [14] and SPOCTOPUS [15]. The main difference between those version is that SPOCTOPUS has an integrated signal peptide prediction. One reason for this approach is the weakness of many transmembrane region prediction methods, that confuse these regions with signal peptides, since both are composed mainly by hydrophobic residues. When this two algorithms were used, the result was similar to the one observed from MEMSAT predictions. Both OCTOPUS and SPOCTOPUS predicted a transmembrane helix around the residues 550 and 600. The difference is that, again, the programs predict the protein upside-down when compared to one another.


Figure 13: Summary of OCTOPUS prediction results. The version (OCTOPUS or SPOCTOPUS) is indicated and the orientation of the molecule is taken from the color code.


<center|font size=2>Figure 13: Summary of OCTOPUS prediction results. The version (OCTOPUS or SPOCTOPUS) is indicated and the orientation of the molecule is taken from the color code.


A third approach is to use the Phobius method. There are two versions of this algorithm, both available online [7]. Phobius [16] and PolyPhobius [17] predict transmembrane helix and signal peptides. From the results (figure 14) we can conclude that this approach is as successful as the others in predicting the transmembrane helix. On the other hand, no signal peptide was predicted.

Figure 14: Summary of Phobius prediction results. The version (Phobius or PolyPhobius) is indicated.


<center|font size=2>Figure 14: Summary of Phobius prediction results. The version (Phobius or PolyPhobius) is indicated.


Although all of the former mentioned prediction methods show a transmembrane helix in the region around residues 550 - 600, it is unlikely that this is really a transmembrane region. First, this is a nuclear protein, which means a low probability of a transmembrane portion. Second, the predicted helix is contained within the ATP binding domain, a well-defined protein domain that takes part in protein function. The presence of a transmembrane helix in this protein region would disturb protein function and it does not seem likely.

The last tool is also the most well-established of the available methods for transmembrane helix prediction. According to its authors, in 2001 TMHMM [20] was rated best in an independent comparison of programs for prediction of transmembrane helices [21]. The results from this method point to the absence of transmembrane helices in the WRN protein. This result is expected for a nuclear protein and, together with other evidences, it leads to the conclusion that no transmembrane helices are present in the WRN protein.


As the WRN protein then was shown not to present transmembrane helices, therefore yielding a negative result from TMHMM, known transmembrane proteins were chosen to give examples of a positive result by this method (figure sup. 1). Five different proteins were assigned as examples, as follows:

1. The human Retinol-binding protein 4 (GI:62298174).

2. The human Insulin-like peptide INSL5 (GI:205371762).

3. A Bacteriorhodopsin from Halobacterium (GI:114811).

4. The human Lysosome-associated membrane glycoprotein 1 (GI:206729915).

5. The human Amyloid beta A4 protein (GI:112927).


Figure 15: TMHMM results for the five proteins briefly described above.


Figure 15: TMHMM results for the five proteins briefly described above.


Figure 15: TMHMM results for the five proteins briefly described above.


Figure 15: TMHMM results for the five proteins briefly described above.


Figure 15: TMHMM results for the five proteins briefly described above.


<center|font size=2>Figure Sup. 1: TMHMM results for the five proteins briefly described above (in order of listing). The color scheme can be observed in the first plot, with red standing for the predicted transmembrane regions.


The same analysis was performed with OCTOPUS and SPOCTOPUS and the results are as following (figure sup. 2). The most interesting observation that comes from the analysis of these results is the advantage of SPOCTOPUS prediction method. As it is obvious from the plots below, OCTOPUS erroneously predicts several transmembrane helices that are in fact signal peptides according to SPOCTOPUS. For the first two sequences, the retinol-binding protein and the insulin-like peptide, TMHMM predicts no transmembrane helices, while OCTOPUS predicts both sequences to contain one N-terminal helix each. This hypothesis is then concluded to be an artifact after the analysis of SPOCTOPUS results, that point to the presence of a signal peptide, in both cases. This summary of results can be said to justify the development of methodologies that not only predict transmembrane helices, but also signal peptides, which can be pointed as false positives in this context.


Figure 15: results for the five proteins briefly described above.
Figure 15: results for the five proteins briefly described above.
Figure 15: results for the five proteins briefly described above.
Figure 15: results for the five proteins briefly described above.
Figure 15: results for the five proteins briefly described above.


Figure 15: results for the five proteins briefly described above.
Figure 15: results for the five proteins briefly described above.
Figure 15: results for the five proteins briefly described above.
Figure 15: results for the five proteins briefly described above.
Figure 15: results for the five proteins briefly described above.


<center|font size=2>Figure Sup. 2: OCTOPUS (top 5 plots) and SPOCTOPUS (bottom 5 plots) results for the five proteins briefly described above (in order of listing). The color scheme can be observed in the first plot. It is interesting to compare the OCTOPUS results with the correspondent SPOCTOPUS results and observe the false positive transmembrane helices predicted by OCTOPUS and revealed (as signal peptides) by SPOCTOPUS.



Signal Peptide Prediction

All proteins are synthesized in the cytoplasm. But later these molecules can have different subcellular localizations. To correctly target proteins to a specific organelle inside the cell, or even to the extracellular space, there are signals within the sequence (or structure) of these proteins. To perform the task of searching for such signals in protein sequences, several computational methods were designed. These methods take advantage of specific properties of the localization signals to predict the subcellular localization of a protein, given its sequence. In this context, four different programs were used specifically to predict localization signals within the sequence of the WRN protein. Additionally (as discussed) some transmembrane helix prediction tools also report the presence of such signals as an additional information. As discussed, the WRN protein is located inside the nucleus. Therefore, it has a Nuclear Localization Signal (NLS) within its sequence (figure 3). Although this is known and it was proven by experimental results, none of the programs above (MEMSAT, SPOCTOPUS, Phobius and PolyPhobius) were able to predict the NLS. This is probably due to limitations of the methodologies and to address this, other approaches were tested.

As an alternative, SignalP [18] was used. This is a well established and widely used program. Once again, the NLS was not identified by the method (results not shown).

One final attempt was to use the Predict Protein server [19]. This time, the NLS was found by the server. The short sequence (ERKRRL) starts in position 1401 in the protein sequence (figure 15). One possible explanation for the lack of results with other methods would be that the methods employed before were not specifically able (trained) to identify the NLS, which is a specific type of localization signal. Therefore, although suitable for the identification of several distinct localization signals, these methods can fail to predict the presence of a signal to nuclear transport.


Figure 15: Prediction of a NLS within the WRN protein sequence, made by the Predict Protein server.


<center|font size=2>Figure 15: Prediction of a NLS within the WRN protein sequence, made by the Predict Protein server.


The fact that other predictive tools were not able to address the presence of this signal peptide (which was known to be present) can be associated to limitations of the techniques. Some algorithms are actually not able to predict NLS with the same accuracy they predict other localization signs. This hypothesis can be further confirmed by studying the references of each method.


GO Assignment

Gene Ontology (GO) terms assignment is (as defined in Wikipedia [8]) an initiative to unify the representation of gene and gene product attributes across all species. In other words, it is a classification of genes and proteins according to the biological roles they play. This classification can occur according to three major contexts, namely:

Cellular component, the parts of a cell or its extracellular environment;

Molecular function, the elemental activities of a gene product at the molecular level, such as binding or catalysis;

Biological process, operations or sets of molecular events with a defined beginning and end, pertinent to the functioning of integrated living units.


A set of programs was used to assign GO terms to the WRN protein: Pfam [23], Protfun [24] and GOPET [25]. The predictions are represented in the figures below (figures 16, 17 and 18).


Figure 16: Domain class assignment results as predicted by the Pfam server.


<center|font size=2>Figure 16: Domain class assignment results as predicted by the Pfam server.



Figure 17: Gene Ontology (GO) assignment results as predicted by the ProtFun server.


<center|font size=2>Figure 17: Gene Ontology (GO) assignment results as predicted by the GOPET server.



Figure 18: GO assignment results as predicted by the GOPET server.


<center|font size=2>Figure 18: GO assignment results as predicted by the GOPET server.


Disorder Prediction

One important observation about disorder prediction is that the terminal regions of the protein (N-terminal and C-terminal) are often referred to as disordered due to the intrinsic flexibility of this regions [26]. It is important to address the accuracy of such results and treat the terminal regions with a different perspective when analyzing the results.

The method used by IUPRED is based on the assumption that globular proteins usually present a large number of interactions between residues. This is the counterpart for the loss of entropy associated to the folding process. In contrast, disordered regions have unique sequences, that in general do not form the same level of interresidue interactions. With this distinction and a set of known globular proteins (with no disorder), the algorithm can predict the presence of disordered regions [27].


Figure 19: Disorder prediction of short regions performed by IUPRED.


<center|font size=2>Figure 19: Disorder prediction of short regions performed by IUPRED.




Figure 20: Disorder prediction of long regions performed by IUPRED.


<center|font size=2>Figure 20: Disorder prediction of long regions performed by IUPRED.


The Disopred server was used as another predictive tool for disordered regions (figure 21 and 22). The predictions were approximately consistent with the results produced by IUPRED. This is an indication of the validity of these predictions and points to a disordered region present in the protein.


Figure 21: Disorder prediction performed by Disopred. This partial result shows the predicted disordered region of the WRN protein.


<center|font size=2>Figure 21: Disorder prediction performed by Disopred. This partial result shows the predicted disordered region of the WRN protein.




Figure 22: Disorder prediction performed by Disopred.


<center|font size=2>Figure 22: Disorder prediction performed by Disopred.


As a third approach, another method was used. Poodle is a system that predicts disorder regions based on three strategies: short disorder regions prediction, long disorder regions prediction and unfolded protein prediction [29]. Once again, the same main region is predicted as disordered. Additionally, other shorter regions towards the C-terminal portion are also predicted and are consistent with the results of the other predictors.

Figure 23: Disorder prediction performed by Poodle.


<center|font size=2>Figure 23: Disorder prediction performed by Poodle.



To further confirm this prediction, the protein structure was analyzed and what one can observe is a structural gap precisely in the region predicted as disordered. There is no reliable entry in PDB containing the region between residues 250 and 500 of the WRN protein (figure 24) and this my be due to the presence of a disordered region. The same conclusion can be drawn from the observation of the shorter disordered regions. As the structure related to the best (first) hit in this BLAST search (against the PDB database) was solved by X-ray diffraction, a crystal had to be generated, containing several copies of the WRN protein. The presence of a disordered region can lead to a low resolution region in the diffraction plot. In fact, disordered regions were once defined based on this assumption. According to Keith Dunker, a disordered region can be predicted from the observation of a structural gap in a crystal structure [30].


Figure 24: BLAST result for the WRN protein using the database of sequences with related structure present in PDB.


<center|font size=2>Figure 24: BLAST result for the WRN protein using the database of sequences with related structure present in PDB.



References

10. Wikipedia

11. Jones DT (1999). Protein secondary structure prediction based on position-specific scoring matrices. Journal of Molecular Biology, 292:195-202.

12. Cole C, Barber JD and Barton GJ (2008). The Jpred 3 secondary structure prediction server. Nucleic Acids Res, 35(2):W197-W201

13. Nugent T and Jones DT (2009). Transmembrane protein topology prediction using support vector machines. BMC Bioinformatics. 10:159

14. Viklund H and Elofsson A (2008). Improving topology prediction by two-track ANN-based preference scores and an extended topological grammar. Bioinformatics

15. Viklund H, Bernsel A, Skwark M and Elofsson A (2008). A combined predictor of signal peptides and membrane protein topology.

16. Käll L, Krogh A and Sonnhammer ELL (2004). A Combined Transmembrane Topology and Signal Peptide Prediction Method. Journal of Molecular Biology, 338(5):1027-1036.

17. Käll L, Krogh A and Sonnhammer ELL (2005). An HMM posterior decoder for sequence feature prediction that includes homology information. Bioinformatics, 21(1):i251-i257.

18. Nielsen H, Engelbrecht J, Brunak S and von Heijne G (1997). Identification of prokaryotic and eukaryotic signal peptides and prediction of their cleavage sites. Protein Engineering, 10:1-6.

19. Rost B, Yachdav G and Liu J (2003). The PredictProtein Server. Nucleic Acids Research, 32:W321-W326.

20. Krogh A, Larsson B, von Heijne G, and Sonnhammer ELL (2001). Predicting transmembrane protein topology with a hidden Markov model: Application to complete genomes. Journal of Molecular Biology, 305(3):567-580.

21. Moller S, Croning MDR, Apweiler R (2001). Evaluation of methods for the prediction of membrane spanning regions. Bioinformatics, 17(7):646-653.

22.

23. Sonnhammer ELL, Eddy SR and Durbin R (1997). Pfam: a comprehensive database of protein families based on seed alignments. Proteins, 28:405-420.

24. Jensen LJ, Gupta R, Blom N, Devos D, Tamames J, Kesmir C, Nielsen H, Stærfeldt HH, Rapacki K, Workman C, Andersen CAF, Knudsen S, Krogh A, Valencia A and Brunak S.(2002) Ab initio prediction of human orphan protein function from post-translational modifications and localization features. Journal of Molecular Biology, 319:1257-1265.

25. Vinayagam A, del Val C, Schubert F, Eils R, Glatting KH, Suhai S and König R (2006). GOPET: A tool for automated predictions of Gene Ontology terms. BMC Bioinformatics, 7:161.

26. Christian Schaefer

27. Dosztányi Z, Csizmók V, Tompa P and Simon I (2005). IUPred: web server for the prediction of intrinsically unstructured regions of proteins based on estimated energy content. Bioinformatics, 21:3433-3434.

28. Ward JJ, Sodhi JS, McGuffin LJ, Buxton BF and Jones DT (2004). Prediction and functional analysis of native disorder in proteins from the three kingdoms of life. Journal of Molecular Biology, 337:635-645.

29. Hirose1 S, Shimizu K, Kanai S, Kuroda Y and Noguchi T (2007). POODLE-L: a two-level SVM prediction system for reliably predicting long disordered regions. Bioinformatics, 23(16): 2046-2053.

30. Rani MPR, Obradovic Z, and Dunker AK (1998) Annotation of PDB with respect to "Disordered Regions" in proteins. Genome Informatics 9:240-241.



COMPARATIVE MODELING

07.06


About Comparative Modeling

The functional characterization of a protein from its sequence is one of the most frequent problems in biology. This task usually becomes easier when a tridimensional structure for the protein is available. In the absence of an experimentaly solved structure, comparative modeling techniques can be employed, given that the protein is related to at least one known structure (solved by X-ray or nuclear magnetic resonance). Comparative modeling methods are able to predict, with a certain level of probability, the tridimensional structure of a given protein sequence (target) based mainly on its alignment to one or more proteins of known structure (template). The protocol for comparative modeling algorithms is classically divided in four steps: i) the identification of template sequence(s) with the highest level of similarity to the target protein; ii) the generation of a sequence alignment between those proteins; iii) the generation of a set of candidate structures for the target protein, based on the previously generated alignment; iv) the prediction of errors and the selection of the final tridimensional model from the set of candidate structures.

Several computational approaches are available to deal with the modeling process. These tools can require a higher or lower level of intervention from the user and this interaction is usually inversely proportional to the user experience. Ultimately, manual intervention is still extremely important [31]. In the most simple cases, the requires input data consists of the alignment between the target sequence and its template(s) and the atomic coordinates of the protein structure from the template protein(s). Although some web servers are able to start "from scratch", thus requiring only the protein sequence for the target.


Template Assignment

As mentioned, the first step is the template assignment that could be performed for example by a BLAST search against the PDB database to find related sequences with known structures. Another useful tool for template assignment is HHPRED that can be accessed through its web server version. Both algorithms were used and the results were compared in order to reach an unique set of template structures. The figure below (figure 25) shows the templates assigned by each of the methods. Although the results can present obvious differences, after a deeper analysis, it is possible to establish connections between the assigned templates within the PDB database.



Figure 25: Summary of results for template assignment obtained from the BLAST search against the PDB database and from HHPRED.


<center|font size=2>Figure 25: Summary of results for template assignment obtained from the BLAST search against the PDB database and from HHPRED.


From this schematic view, it is easy to notice that the second result is the most extensive since the template proteins cover more of the target sequence. As stated, the templates assigned by each method (BLAST and HHPRED) are somehow related. For the first template protein, 2FBT according to BLAST or 2E6M according to HHPRED, the informations on PDB state that these are analogous proteins that share 70% of sequence similarity. For the second template, the longest one, the proteins are actually 100% similar. For the last template, the PDB entry 1D8B is also related to the entry 2E1E. Additionally, for the result retrieved from HHPRED, the end portion of the WRN protein could be modeled by the protein with PDB identifier 3IUO, the C-terminal domain of the ATP-dependent DNA helicase from Porphyromonas gingivalis . Although, as stated, the templates suggested by HHPRED cover more of the target sequence, the similarity level between the template sequence and the target sequence is low, in general. For that reason, the templates assigned by the BLAST search against the PDB database were chosen. They are as follows (ans as summarized in figure 25):


2FBT - Covers the region between residues 38 and 236 of the WRN protein with 199 identical residues from the 199 aligned residues (100% identity between the sequences). It corresponds to the exonuclease domain of the human WRN.


1OYY - Covers the region between residues 540 and 1030 of the WRN protein with 37% (183 out of 498) identical residues, 56% positive residues and 7% of open gaps in the alignment with the WRN protein. It is annotated as the catalytic core of the E.coli RecQ helicase.


2E1E - Covers the region between residues 1142 and 1242 of the WRN protein with 100% coverage and 100% identity. It stands for the crystal structure of the HRDC Domain of human WRN protein.


Below there is a summary of the regions predicted from homology searches, disorder predictions and signal peptide searches. A direct consequence of the presence of disordered regions is that these portions of the protein will not be correctly modeled, since there is no suitable template for it.


Figure 26: Regions in the WRN protein.


<center|font size=2>Figure 26: Regions in the WRN protein.



When the templates are assigned, the next step is to generate the sequence alignment between these proteins. Unfortunately, none of the fore mentioned programs was able to correctly align the target sequence to the multiple templates.


Model Generation

The first approach to generate the tridimensional models is to use the iTasser server, which is an automated method for generating models from protein sequences. The server generated five different structures which are shown below (figure 27).


Figure 27: Miniature figure of the five models generated by iTasser.


<center|font size=2>Figure 27: Miniature figure of the five models generated by iTasser.


As expected, the model contains 'unfolded' regions which come as a result of the lack of template structures. The long disordered region cannot be modeled since there is no suitable template structure available for it (since it is not possible to crystallize such flexible region). This issue will set an obstacle for the modeling process. As it is not possible to generate a structure that covers the entire protein, the model will have to be build in parts. The schematic view of the protein sequence and the defined domains it contains (figure 26) shows two portions of the protein, separated by the disordered region. The first portion enclosures an exonuclease domain, which is an unique feature of the WRN protein. Among the five human RecQ DNA helicases, this is the only one presenting a functional exonuclease domain. This domain was characterized in its structure, biochemical activity and participation in DNA end joining [32]. This region is isolated from the rest of the protein by the disordered region. The C-terminal half of the protein contains at least three domains and the signal peptide (NLS). These domains are (in order from the N-terminal) the RecQ domain which is the longest domain whithin the protein sequence, the helicase-and-ribonuclease D/C-terminal (HRDC) domain, of unknown function (recently pointed out to be important for dissolution of double Holliday junctions [33]) and a domain which was only characterized as a C-terminal domain. After these domains (briefly separated from each other by a few residues predicted as secondary disordered regions) is the nuclear localization signal. Therefore, one alternative would be to build two separate models. One covering the exonuclease domain and the other covering the RecQ, HRDC and C-terminal domains of the WRN protein.

To build the candidate structures, several approached were employed, as suggested: iTasser [34], SwissModel [35] and Modeller [36]. These programs can be pointed as the most important methods currently available. Modeller is perhaps the most well-established comparative modeling algorithm, available as a "stand-alone" program. It requires the user to be more familiar to "script" language (Python scripts are used by Modeller) and to the process of comparative modeling itself. SwissModel is the most widely used tool for comparative modeling as a web server. It allows the user to choose between three versions that require different levels of interaction between the user and the program. One can start "from scratch" with only the protein sequence as an input, or in a more advanced version, the alignment between target and template can be generated by the user. The iTasser is also available as a web server and, although it is a relatively recent method, it has performed extremely well within the last CASPs (it was ranked as the best server for protein structure prediction in recent CASP7, CASP8 and CASP9 [37] experiments).


iTasser

As already shown in the last figure above (figure 26), the first results were generated by iTasser. The entire protein sequence was submitted to the server and the generated structures present (as expected) long unfolded regions, due to the lack of template which can be in turn attributed to the presence of disordered regions within the protein. For this reason, the protein sequence was subdivided, as mentioned, in two regions. The first subsequence corresponds to the exonuclease domain, positioned in the N-terminal end of the sequence. The second subsequence comprehends the other three identified domains, namely RecQ, HRDC and C-terminal domains. These subsequences were submitted to the server independently, generating two sets of candidate structures, shown in the figure below (figure 27)


Figure 28: Miniature figure of the five models generated by iTasser for the .


<center|font size=2>Figure 27: Miniature figure of the five models generated by iTasser.


Modeller

As stated, Modeller works as a package, to be run locally. For the most basic variant, the program needs as input only the protein sequence. It automatically builds the sequence alignment between the target and the template(s) and performs all other steps until the generation of the candidate structures. Although this would be the easiest way of using Modeller, in this case, the alignment between WRN and the templates was generated manually, since the programs failed to correctly assign each template to its correspondent region of the target protein. The alignment was therefore built considering four templates, one for each domain, as assigned in the previous step with BLAST and HHPRED. The gaps were inserted according to the prediction tools (BLAST and HHPRED) within the locally aligned sequences. One additional step was required by Modeller and this is a common issue to be assessed before the structures can be generated. As the templates are proteins with known tridimensional structures, they are available in PDB and can be retrieved from the database. In a more detailed analysis of the coordinates file (from PDB), it is easy to notice that frequently some residues are missing within the structure. Obviously, when this is the case, the alignment between the available sequence of the protein and the sequence generated from the PDB structure file will present gaps, since these sequences will not be exactly the same. This is due to technical problems with the method of X-ray diffraction, but Modeller is not ready to automatically deal with this divergence. What happens is an error in the program that forces it to terminate without actually generating the structures. To overcome this problem, the template sequences were retrieved from the PDB file (therefore, only the residues that are certainly present in the structure were in the sequence). The only template protein that presented such divergence between sequence and structure was 1OYY and there were four residues missing in the structure (figure 28). The figure below shows the alignment between the protein sequence (as retrieved from PDB) and the sequence generated from the PDB file containing the coordinates of the molecule. It is obvious that these four residues are not present in the solved structure and they were therefore deleted in the alignment between the templates and the target.


Figure 28: Alignment between the sequence of 1OYY retrieved from PDB and the sequence built with the information from the structural PDB file. The only gap within the alignment represents the four residues that are missing in the crystal structure.


<center|font size=2>Figure 28: Alignment between the sequence of 1OYY retrieved from PDB and the sequence built with the information from the structural PDB file. The only gap within the alignment represents the four residues that are missing in the crystal structure.


After solving the problem and deleting the absent residues from the template sequence, the program was able to generate the set of candidate structures for the target sequence. Modeller requires a Python script to set the parameters for each run. In this case, the script is as ilustrated in the figure below (figure 29). It contains simple instructions that are explained in commented lines (starting with #). The sequence alignment between target and templates was generated manually, as discussed above and 10 structures were generated. Different protocols were used, varying the alignment file, the set of templates and the portion of the WRN protein to be modeled (as discussed, the protein was segmented in two main regions and different structures were generated for each region).


Figure 29: Example of a Python script submitted to Modeller as an input to generate partial structural candidates for the WRN protein.


<center|font size=2>Figure 28: Example of a Python script submitted to Modeller as an input to generate partial structural candidates for the WRN protein.



References

32. Perry JJ, Yannone SM, Holden LG, Hitomi C, Asaithamby A, Han S, Cooper PK, Chen DJ, Tainer, JA (2006) WRN exonuclease structure and molecular mechanism imply an editing role in DNA end processing. Nature Structural Molecular Biology 13:414-422

33. Wu L, Chan KL, Ralf C, Bernstein DA, Garcia PL, Bohr VA, Vindigni A, Janscak P, Keck JL and Hickson ID (2005). The HRDC domain of BLM is required for the dissolution of double Holliday junctions. The EMBO Journal, 24:2679 - 2687

34. Roy A, Kucukural A, Zhang Y (2010). I-TASSER: a unified platform for automated protein structure and function prediction. Nature Protocols, 5:725-738.

35. Arnold K, Bordoli L, Kopp J, and Schwede T (2006). The SWISS-MODEL Workspace: A web-based environment for protein structure homology modelling. Bioinformatics, 22:195-201.

36. Sali A and Blundell TL (1993). Comparative protein modelling by satisfaction of spatial restraints. Journal of Molecular Biology 234:779-815.

37. Complete Issue (2009). Proteins: Structure, Function, and Bioinformatics Volume 77, Issue S9, Pages 1-228




MAPPING SINGLE NUCLEOTIDE POLYMORPHISMS

14.06


OMIM - Online Mendelian Inheritance in Man

OMIM [38 [9]] is a comprehensive, authoritative, and timely compendium of human genes and genetic phenotypes. The full-text, referenced overviews in OMIM contain information on all known mendelian disorders and over 12,000 genes. OMIM focuses on the relationship between phenotype and genotype. It is updated daily, and the entries contain copious links to other genetics resources.

When a search is performed for "Werner Syndrome" within the OMIM database (MIM ID #277700), a rich text is retrieved, with a great amount of data from the literature. Upon this information, the following summary was written, addressing the most important points of the disease:


The features of Werner syndrome are scleroderma-like skin changes, cataract, subcutaneous calcification, premature arteriosclerosis, diabetes mellitus, and a wizened and prematurely aged facies (short stature, slender limbs, stocky trunk and beaked nose). Autosomal recessive inheritance was confirmed by Goto and colleagues, in a 1981 study comprised by 42 Japanese families among which 80 subjects were affected. Malignancy was frequent in the patients and in the families generally. The frequency of Werner syndrome in Japan was estimated to be about 3 per million persons. In 1997, Martin gave a thoughtful review of the question of whether the Werner mutation is a reflection of 'normal aging' mechanisms and 4 years later, Mohaghegh and Hickson reviewed the DNA helicase deficiencies associated with cancer predisposition and premature aging disorders.

As for the cellular phenotype related to the syndrome, 'variegated translocation mosaicism' was the designation proposed by W. W. Nichols (Hoehn et al., 1975) for a phenomenon he and others observed in cells from patients: skin fibroblast cell lines were usually composed of one or several clones, each marked by a distinctive, apparently balanced translocation. Salk (1982) found that somatic cells from Werner syndrome patients reveal a propensity to develop chromosomal aberrations, including translocations, inversions, and deletions. Six years later, Gebhart and colleagues concluded that Werner syndrome cells demonstrate some biochemical differences that distinguishing them from those of other classic chromosome instability syndromes. Subsequent studies suggested an increased number spontaneous chromosome rearrangements and deletions in WRN cells. It was later suggested that the deletion mutator uses deletion mutagenesis pathways that are similar or identical to those used in other human somatic cells. In 2007, Crabbe and colleagues proposed a mechanism in which lack of WRN helicase activity results in dramatic telomere loss, causing a DNA damage and repair response that leads to genomic instability and may culminate in cancer.

Matsumoto et al. (1997) presented evidence that the helicase that is defective in Werner syndrome is missing the nuclear localization signal (NLS) and that this leads to impaired nuclear import as a major contributing factor in the molecular pathology of the disorder. The finding helped to explain the enigma that most Werner syndrome patients have similar clinical phenotypes no matter how different their mutations. New findings from 2000 suggested that one consequence of the Werner syndrome defect is an acceleration of normal telomere-driven replicative senescence. A few years later, though, Baird and colleagues proposed that telomere dynamics at the single-cell level in WRN fibroblasts are not significantly different from those in normal fibroblasts, and suggested that the accelerated replicative decline seen in WRN fibroblasts may not result from accelerated telomere erosion.

In 1996, Yu et al. identified 4 mutations in patients with Werner syndrome. Two of which were splice-junction mutations with the predicted result being the exclusion of exons from the final messenger RNA. One of these mutations resulted in a frameshift and a predicted truncated protein. The other 2 were nonsense mutations. The identification of a mutated putative helicase as the gene product of the WRN gene suggested to the authors that defective DNA metabolism is involved in a complex process of aging in Werner syndrome patients. In the same year, Oshima and colleagues reported 9 new WRN mutations, located at different sites across the coding region (in contrast to previously known mutations, which either create a stop codon or cause frameshifts that lead to premature terminations). They noted that the WRN protein is partially homologous to RecQ helicases which act in a number of molecular processes, including unwinding of DNA during replication, DNA repair, and accurate chromosomal segregation. A year later, Yu et al. identified 5 new WRN mutations, four of which either partially disrupted the helicase domain region or resulted in predicted protein products lacking the entire helicase region. These results confirmed that mutations in the WRN gene are responsible for Werner syndrome. In addition, the location of the mutations indicated that the presence or absence of the helicase domain does not influence the Werner syndrome phenotype, suggesting that this syndrome is the result of complete loss of function of the WRN gene product.

Kyng et al. (2003) found that fibroblasts from 4 patients with Werner syndrome and fibroblasts from 5 older control individuals (average age 90 years) showed transcription alteration of 435 (6.3%) of 6,192 genes examined compared to cells from young adult controls. Of the 435 genes, 91% of the 249 genes with known function had similar transcription changes in both Werner syndrome patients and normal old age controls. The major functional categories of the similarly transcribed genes of known function included DNA/RNA metabolism, cell growth, and stress response. Werner syndrome may therefore be a good model for normal aging, since both processes are linked to altered transcription

Lombard et al. (2000) generated mice bearing a mutation that eliminated expression of the C terminus of the helicase domain of the WRN protein. Mutant mice were born at the expected mendelian frequency and did not show any overt histologic signs of accelerated senescence. The mice were capable of living beyond 2 years of age. Cells from these animals did not show elevated susceptibility to 2 genotoxins. However, mutant fibroblasts aged approximately 1 passage earlier than controls. Importantly, mice that were doubly homozygous for WRN and p53 (191170) deficiencies showed an increased mortality rate relative to animals that were heterozygous for WRN deficiency and homozygous for p53 null. Lombard et al. (2000) considered possible models for the synergy between p53 and WRN mutations for the determination of life span.

The references for the cited articles can be found at the OMIM page related to Werner Syndrome [10].


Werner Syndrome Mutation Database

The Werner Syndrome Mutation Database website [39] is a locus-specific mutational database (LSMD) developed by the Washington University in conjunction with the Human Genome Organization (HUGO) and the Human Genome Variation Society (HGVS). From this website, was retrieved the figure below (30), which summarizes WRN mutations found in Werner syndrome patients, from two major sources (indicated by open and closed symbols), and indicates the position of these mutations and SNP's (asterisks) in the WRN open reading frame and mRNA. Additional information on individual mutations and SNPs can be found in this table Media:Mutations.pdf, in a pdf format.


Figure 30: Summary of known mutations of the WRN protein.


<center|font size=2>Figure 30: Summary of known mutations of the WRN protein.


HGMD - Human Gene Mutation Database

As described at the website, the Human Gene Mutation Database (HGMD) [40] includes the first example of all mutations causing or associated with human inherited disease, plus disease-associated/functional polymorphisms reported in the literature. These data comprise various types of mutation within the coding regions, splicing and regulatory regions of human nuclear genes. HGMD does not usually include mutations lacking obvious phenotypic consequences although a few such variants have been included where they could conceivably have some clinical effect.

For the WRN protein encoding gene, 61 mutations are reported at the HGMD, the majority of those are nonsense or missense (figure 31). As described by the authors, this database does not include neutral mutations (mutations that apparently do not affect the protein function).


Figure 31: Search result for mutations within the WRN protein at the Human Gene Mutation Database.


<center|font size=2>Figure 31: Search result for mutations within the WRN protein at the Human Gene Mutation Database.


dbSNP - Single Nucleotide Polymorphism Database

The Single Nucleotide Polymorphism Database (dbSNP) [41] is a free public archive for genetic variation within and across different species developed and hosted by the National Center for Biotechnology Information (NCBI) in collaboration with the National Human Genome Research Institute (NHGRI). Although the name of the database implies a collection of one class of polymorphisms only (i.e., single nucleotide polymorphisms (SNPs)), it in fact contains a range of molecular variation: (1) SNPs, (2) short deletion and insertion polymorphisms (indels/DIPs), (3) microsatellite markers or short tandem repeats (STRs), (4) multinucleotide polymorphisms (MNPs), (5) heterozygous sequences, and (6) named variants.[2] The dbSNP accepts apparently neutral polymorphisms, polymorphisms corresponding to known phenotypes, and regions of no variation. It was created in September 1998 to supplement GenBank, NCBI’s collection of publicly available nucleic acid and protein sequences [42]

From the dbSNP database, 52 mutations were retrieved for the WRN protein.


WRN Mutations Summary

As described, three major databases were consulted regarding information about mutations at the WRN protein-coding gene. These databases are the Werner Syndrome Mutation Database (WSMD herein) the Human Gene Mutation Database (HGMD) and the Single Nucleotide Polymorphism Database (dbSNP). The total amount of mutations found was 129, divided in several categories: Synonymous, Missense, Nonsense, Frameshift, Short/Long Deletion, Short/Long Insertion and Short Indel and Splicing.

Among the 129 unique mutations, 52 were retrieved from dbSNP and 61 from HGMD. From the WSMD database, only 21 entries were unique and therefore added to the list of unique mutations. The overlap between the data available at dbSNP and HGMD is extremely small. Only 5 mutations were common to both databases (with the exception that one mutation was L --> F at dbSNP and F --> L at HGMD). The figure below (figure 32) summarizes these numbers and the table of all unique mutations can be downloaded as a PDF file Media:WRN_Mutations.pdf.


Figure 32: Summary of origin of unique mutations found for the WRN protein coding gene in three different databases (dbSNP, HGMD and WSMD).


<center|font size=2>Figure 32: Summary of origin of unique mutations found for the WRN protein coding gene in three different databases (dbSNP, HGMD and WSMD).



References

38. Online Mendelian Inheritance in Man, OMIM (TM). McKusick-Nathans Institute of Genetic Medicine, Johns Hopkins University (Baltimore, MD) and National Center for Biotechnology Information, National Library of Medicine (Bethesda, MD)

39. http://www.pathology.washington.edu/research/werner/database/

40. HGMD

41. dbSNP

42. http://en.wikipedia.org/wiki/DbSNP


ABOUT RECQ HELICASES

Helicases

All helicases share the ability to unwind nucleic acid duplexes with a distinct directional polarity; they utilize the free energy from nucleoside triphosphate hydrolysis to fuel their translocation along DNA, unwinding the duplex in the process. The WRN protein has helicase activity in the 3'-5' direction


Exonucleases

Exonucleases are enzymes that work by cleaving nucleotides of a polynucleotide chain through a hydrolyzing reaction that breaks phosphodiester bonds at either the 3' or the 5' end. The 3' to 5' human type endonuclease is known to be essential for the proper processing of histone pre-mRNA, in which U7 snRNP directs the single cleavage process. In both archaebacteria and eukaryotes, one of the main routes of RNA degradation is performed by the multi-protein exosome complex, which consists largely of 3' to 5' exoribonucleases. The RNase D exonuclease domain (Pfam) is responsible for the 3'-5' exonuclease proofreading activity of E. coli DNA polymerase I (figure 33) and other enzymes. It catalyses the hydrolysis of unpaired or mismatched nucleotides. This domain consists of the amino-terminal half of the Klenow fragment in E. coli polI it is also found in the WRN protein and ribonuclease D (RNase D).


Human RecQ Helicases

The RecQ helicase family is widely distributed across the three domains of life. The human genome contains five RecQ genes, RecQ1 (Uniprot Q9H9A7), BLM (P54132), WRN (Q14191), RecQ4L (O94761) and RecQ5 (O94762). Mutations in BLM and RecQ4L cause Bloom syndrome and Rothmund-Thomson syndrome, espectively. Werner, Bloom and Rothmund-Thomson syndromes share a predisposition to cancer, but notable pathological differences suggest that each disease pathway is functionally distinct. Biochemical characterization of the human RecQ helicases has shown ATPase activity and unwinding of partialduplex DNA substrates with 3'-5' polarity. Alternate DNA conformations are preferred over double-stranded DNA (dsDNA), suggesting roles for human RecQ helicases in replication, recombination and repair events. According to the Conserved Domain Search results (CDS) retrieved from NCBI, each of the five human RecQ helicases have unique domains. The majority (4 out of the 5, except the RecQ1 protein) of these proteins present two conserved domains, characteristic of this enzyme superfamily (figure 33). The DEXDc domain defines the DEAD-like helicases superfamily (to which all but RecQ1 human helicases are likely to be pertain). These are proteins involved in ATP-dependent RNA and DNA unwinding and this domain enclosures the ATP binding site. The second conserved region corresponds to the helicase superfamily c-terminal (HELICc, cd00079) domain, which is an integral part of a wide variety of helicases and helicase-related proteins. The WRN protein is a unique DNA helicase possessing exonuclease activity. The DnaQ-like, 3'-5' RNaseD exonuclease domain (cd06129) contains three conserved sequence motifs clustered around the active site containing four conserved acidic residues that serve as ligands for the two metal ions required for catalysis. RNaseD is involved in the 3'-end processing of tRNA precursors.

Figure 33: Domains at the five human helicases, as predicted by the CDS tool from NCBI.


<center|font size=2>Figure 33: Domains at the five human helicases, as predicted by the CDS tool from NCBI.


As this is an unique domain for the WRN protein and as it was structurally characterized (PDB 2FDT), the subsequent steps will consider only this protein region (from residue 38 to 236), unless otherwise stated. This restriction will focus the investigation of mutation effects.


WRN Protein and FEN-1

The cellular defects of WS presumably reflect compromised or aberrant function of a DNA metabolic pathway that under normal circumstances confers stability to the genome. A novel interaction of the WRN gene product with the human 5' flap endonuclease/5'-3' exonuclease (FEN-1) was reported. FEN-1 is a DNA structure-specific nuclease implicated in DNA replication, recombination and repair. WRN protein dramatically stimulates the rate of FEN-1 cleavage of a 5' flap DNA substrate. The WRN-FEN-1 functional interaction is independent of WRN catalytic function and mediated by a 144 amino acid domain of WRN that shares homology with RecQ DNA helicases. A physical interaction between WRN and FEN-1 is demonstrated by their co-immunoprecipitation from HeLa cell lysate and affinity pull-down experiments using a recombinant C-terminal fragment of WRN. The underlying defect of WS is discussed in light of the evidence for the interaction between WRN and FEN-1.


References

Mainly extracted from the following papers: Perry JJP, Yannone SM, Holden LG, Hitomi C, Asaithamby A, Han S, Cooper PK, Chen DJ & Tainer JA (2006). WRN exonuclease structure and molecular mechanism imply an editing role in DNA end processing. Nature Structure and Molecular Biology, 13(5):414-422.

Brosh RM Jr, von Kobbe C, Sommers JA, Karmakar P, Opresko PL, Piotrowski J, Dianova I, Dianov GL, Bohr VA (2001). SourceWerner syndrome protein interacts with human flap endonuclease 1 and stimulates its cleavage activity. EMBO Journal 15;20(20):5791-801.




SEQUENCE-BASED MUTATION ANALYSIS

21.06


Subset of 10 Mutations

For this task, the first step is to create a subset of 10 different mutations to be further analyzed. From the total list of mutations, one should be more interested in the substitutions of the "Missense" type, which lead to an amino acid substitution at the sequence, affecting only one position (point mutations). Among the 129 mutations listed for the WRN protein, 42 are point mutations causing an amino acid alteration at the sequence. Among these 42 mutations, 7 are located within the exonuclease domain. To expand the set, new databases were consulted, SwissProt and PMD. The mutations from these databases were retrieved from the another database, comprising dbSNP, PMD and SwissProt. With the new results, a complete table of known mutations at the exonuclease domain was built and is shown in a table below together with prediction results.


About Some Mutations

(Information below was extracted from Perry JJP, Yannone SM, Holden LG, Hitomi C, Asaithamby A, Han S, Cooper PK, Chen DJ & Tainer JA (2006). WRN exonuclease structure and molecular mechanism imply an editing role in DNA end processing. Nature Structure and Molecular Biology, 13(5):414-422.)


About the Y212 mutation:

WRN Tyr212 is crucial in directing water-mediated attack of the phosphodiester bond. The hydroxyl group of Tyr212 is vital to the WRN-exo active site mechanism, as a Y212F substitution substantially reduced but did not eliminate WRN-exo activity. This suggests that the Tyr212 side chain reorients into the active site upon substrate binding.


About the L88 mutation:

In a closely related structure (KF-exo), a leucine side chain, Leu361, stacks on the opposing side of the terminal base. However, the equivalent residue in WRN, Leu88, has its Ca positioned 8.5 A ° away with its side chain pointing into solution. WRN Leu88 and KF Leu361 reside in the flexible b2-b3 loop region. The WRN b2-b3 loop (residues 85–98) has a 2-residue insertion, and its ‘open’ conformation is probably stabilized by crystal-packing interactions of loop residue Arg91. A structurebased L88A mutation had no effect on WRN-exo catalytic activity, indicating that other residues in this substrate-binding loop may be important. Modeling of bound substrate indicates that other residues within the b2-b3 loop that are strictly conserved among the WRN homologs, including Tyr89 and Lys93, or the well-conserved Asn90 and Arg91 may also have roles in substrate binding.


About the V114I and T172P mutations:

The authors compared the nuclease activities of the WRN-exo domain and its E84A, Y212F and L88A mutants with WRN-exo bearing the two known polymorphisms that alter amino acids, V114I and T172P. They observed that the two SNPs did not directly alter exonuclease activity of the minimal domain, when compared to these WRN-exo constructs. The isoleucine variant is likely to be accommodated into the core without disturbing the overall WRN-exo structure, although the extra methyl group could disturb packing and possibly decrease stability. The T172P polymorphism has been previously reported to cause a small reduction in the enzymatic activity of full-length WRN protein. As the WRN-exo T172P activity is not appreciably perturbed, the change in full-length WRN activity may be due to disruption of WRN exonuclease-helicase coordination, reduced stability or folding in the larger protein, or both. At PMD, these mutations are annotated as normal polymorphisms, possibly not related to disease.


(Information was extracted from the correspondent PMD entry.)


About the D82A and E84A mutations:

The mutant is catalytically defective, but retains ability to stimulate FEN-1 cleavage reaction. Amino-acid substitutions at either position 82 (D82A) or 84 (E84A), are two of the five mutations predicted to be critical for exonuclease activity.



SIFT and SNAP

Two programs were used to classify this mutations according to their functional effects upon the protein. SIFT (Sorting Intolerant From Tolerant) is a program that predicts whether an amino acid substitution affects protein function so that users can prioritize substitutions for further study. We have shown that SIFT can distinguish between functionally neutral and deleterious amino acid changes in mutagenesis studies and on human polymorphisms [43]. SNAP is a neural-network based method that uses in silico derived protein information (e.g. secondary structure, conservation, solvent accessibility, etc.) in order to make predictions regarding functionality of mutated proteins. The network takes protein sequences and lists of mutants as input, returning a score for each substitution [44]. The summary of results from these sources is shown below.


Figure 33: Domains at the five human helicases, as predicted by the CDS tool from NCBI.



BLOSUM, PAM and PSSM

The BLOSUM (BLOcks of Amino Acid SUbstitution Matrix) [46] is a matrix representing a scoring system for point mutations widely used for sequence alignment of proteins (figure 34). It was derived by Henikoff and Henikoff through a subset of the BLOCKS database containing very conserved regions of protein families. The retrieved BLOCKS were in general short, ungapped fragments. The relative frequencies of each amino acid was calculated together with its substitution probability. All BLOSUM matrices are based on observed alignments. Several sets of BLOSUM matrices exist using different alignment databases, named with numbers. BLOSUM matrices with high numbers are designed for comparing closely related sequences, while those with low numbers are designed for comparing distant related sequences. Scores within a BLOSUM are log-odds scores. Every possible identity or substitution is assigned a score based on its observed frequences in the alignment of related proteins. A positive score is given to the more likely substitutions while a negative score is given to the less likely substitutions. Here BLOSUM 62 was used to generate the results in the table below.

Figure 34: BLOSUM62 scoring matrix retrieved from NCBI.


<center|font size=2>Figure 34: BLOSUM62 scoring matrix retrieved from NCBI.


Point (or Percent) Accepted Mutation (PAM) [47] is another set of matrices used to score sequence alignments (figure 35). The PAM matrices were introduced by Margaret Dayhoff in 1978 based on 1572 observed mutations in 71 families of closely related proteins. The value in a given cell represents the probability of a substitution of one amino acid for another (according to the respective experimental sequences in the dataset used to derive the matrix). The calculation uses the ratio of the probability value and the frequency of the original amino acid in known sequences. PAM 250 was hereby employed to generate the results in the table below.


Figure 35: PAM250 scoring matrix retrieved from NCBI.


<center|font size=2>Figure 35: PAM250 scoring matrix retrieved from NCBI.


A Position Specific Scoring Matrices (PSSM) is a matrix of score values that gives a weighted match to any given substring of fixed length. It has one row for each symbol of the alphabet, and one column for each position in the pattern. In other words, a PWM score is the sum of position-specific scores for each symbol in the substring. The PSSM for the first domain of the WRN protein was obtained from the NCBI website, while running a PSI-BLAST. The PSSM information was visualized and saved as a matrix using the visualization tool from NCBI [48]. The PSSM file can be downloaded as a pdf Media:PSSM.pdf.



Physico Chemical Information

In addition to the previous results, for each substitution, the physico-chemical properties of each amino acid are included in the results, to illustrate the molecular differences between the wild-type residue and the correspondent mutation. The properties for each residue were retrieved from the Venn diagram presented in figure 5, from which a copy is shown below. These properties are summarized in the table below.


Figure 5: Venn diagram for the 20 most common amino acids. In this diagram it is possible to graphically explicit the structural and physico-chemical relationships among the residues, for example.


<center|font size=2>Figure 5 (copy): Venn diagram for the 20 most common amino acids. In this diagram it is possible to graphically explicit the structural and physico-chemical relationships among the residues, for example.



References

43. Pauline C and Henikoff S (2003). SIFT: predicting amino acid changes that affect protein function. Nucleic Acids Research, 31(13):3812-3814.

44. Bromberg Y and Rost B (2007). SNAP: predict effect of non-synonymous polymorphisms on function. Nucleic Acids Research, 35(11):3823-3835.

45.

46. Henikoff S Henikoff JG (1992). Amino Acid Substitution Matrices from Protein Blocks. PNAS, 89(22):10915–10919

47. Dayhoff MO, Schwartz R and Orcutt BC (1978). A model of Evolutionary Change in Proteins. Atlas of protein sequence and structure. Nat Biomed Research Foundation 5(3):345–358.

48. http://www.ncbi.nlm.nih.gov/Class/Structure/pssm/pssm_viewer.cgi


STRUCTURE-BASED MUTATION ANALYSIS

28.06

Mutant Protein Structures

The first step for a comparison between wild-type (WT) and mutant proteins is the generation of tridimensional structures for the mutant proteins. For this purpose, a script was used to correctly assign the side-chain of the new residue. This script makes use of SCWRL, which is (according to its authors) currently the most advanced algorithm for predicting sidechain conformation [49]. All previous mentioned mutations were created. In addition, SCWRL was also used to analyze and correctly place the side-chain atoms according to the backbone provided in the PDB file of the WT protein.

The mutations were mapped into the sequence and can be visualized in the figure below (figure 36).




<center|font size=2>Figure 37: Structure of the exonuclease domain from the WRN protein (PDB 2FBT). The mutations listed in the table above were highlighted in orange and labeled according to their position and WT residue.



Energy Comparison

FOLD-X [50] is a computer algorithm developed to provide a fast and quantitative estimation of the importance of the interactions contributing to the stability of proteins and protein complexes. The predictive power of the algorithm has been tested on a very large set (1088) of point mutants spanning most of the structural environments found in proteins. FOLD-X uses a full atomic description of the structure of the proteins. The different energy terms taken into account in FOLD-X have been weighted using empirical data obtained from protein engineering experiments. The present energy function uses a minimum of computational resources and can therefore easily be used in protein design algorithms, and in the field of protein structure and folding pathways prediction where one requires a fast and accurate energy function. FoldX was used to calculate the energy of all mutant structures and the WT protein (PDB 2FBT). The results are as follows:



From the obtained results, it is easy to notice that a subset of mutations (D82A, E84A, L88A and T172P) has a lower energy than the WT protein. This roughly means that the correspondent mutations decrease protein flexibility (leading to a more stable protein). These mutated enzymes can lose the capacity of correctly perform its function, since they are likely more rigid than the WT protein. On the other hand, for the rest of the studied mutations (G92V, V114I, K125N, K135E, A136T, W149H, Y212F and R225Q), the mutated structure has a higher energy than the WT protein. This points to a more flexible structure and perhaps, in some cases, to an unstable protein unable to fold and/or to function as expected.

It is probably worth noticing that the mutations for which structure has lower energy differences (either upper or lower than the WT) were predicted as having no effect upon the protein function according to either SNAP, SIFT or both. This is easily observed in the table below, extracted from the complete table. The structures were ordered according to the difference in energy from the WT protein. Those mutations with no effect on function are highlighted in red, the WT structure is highlighted in green. One asterisk (*) means that SNAP predicted some effect while SIFT predicted no effect on protein function for the correspondent mutation and two asterisks (**) means both SNAP and SIFT predicted no effect on function. Additionally to the predictions, it is known that both mutations L88A and T172P have no significant effect on protein function.




Gromacs is a widely used computer program also capable of analyzing a protein structure according to its energetic features. Here, Gromacs was employed as a tool to minimize protein energy and therefore generate information about the energetic profile of the molecules. There are several steps to be followed in order to successfully carry out an energy minimization protocol. As there were 14 structures to be analyzed (13 mutations plus the WT protein), a perl script was written to allow the protocol to be as automated as possible. This script has all information needed by Gromacs and accessories programs. Self-explanatory comments were added to the script in order to allow any user to fully understand each step of the protocol. Below is the complete script with all comments (figure 37).


<center|font size=2>Figure 37: Perl Script to minimize energy of all mutant structures with Gromacs.


The results are obtained as the energy value for each time step during the minimization. The last value (the last time step) was used to generate a table from which it is easy to compare the different mutated structures in terms of potential energy.



References

49. Wang Q, Canutescu AA & Dunbrack RL (2008). SCWRL and MolIDE: computer programs for side-chain conformation prediction and homology modeling. Nature Protocols, 3:1832-1847

50. Schymkowitz JW, Rousseau F, Martins IC, Ferkinghoff-Borg J, Stricher F, Serrano L (2005). Prediction of water and metal binding sites and their affinities by using the Fold-X force field. PNAS, 102:10147-10152.



MOLECULAR DYNAMICS

For a detailed description about Molecular Dynamics techniques, download the PDF file here: Media:MD.pdf.


About Molecular Dynamics

Molecular Dynamics simulations have been used for more than 3 decades to study the molecular motion in the most diverse systems. In Biology, a major use for this computational technique is the examination of protein dynamics on the atomic level. The first Molecular Dynamics (MD) simulation of a biological macromolecule was carried out in 1977 by McCammon. It was a 9.2 ps long simulation of the Bovine Pancreatic Trypsin Inhibitor (BPTI), a small, highly stable protein structure that has been widely used as a case-study for new developments in MD. This was the first consistent analysis of a protein molecule as a dynamical system, although the concept of proteins as flexible entities was already present. As Richard Feynman addressed in his famous book The Feynman Lectures on Physics more than 10 years before this first MD simulation of a protein, “All things are made of atoms and everything that living things do can be understood in terms of the jigglings and wigglings of atoms”. Following this elucidation, in 1955 Landerstrom-Lang waas able to measure the solvent accessibility of protein residues through his innovative technique called Hydrogen-Deuterium Exchange. This was probably the first experimental evidence of a dynamical protein world as predicted by Feynman. Around the decade of 1970, biologists were already familiar with protein sequencing techniques and experimental determination of protein structures. The first protein to be sequenced was the bovine insulin, in 1951 by Frederick Sanger, while the first protein structures were determined by X-ray crystallography less than 10 years after that. These were the first ingredients for the development of an in silico approach for observation of protein dynamics. The first use of MD simulations was described 1957 through the work entitled Phase Transitions for a Hard Sphere System, by B.J. Alder and T. E. Wainwright, who were studying the interactions between solid spheres. Subsequently to the application of the new developed technique in this simple system, there was a gap of almost 20 years before its application to a representation of a real system. In 1974, the MD simulation of a liquid water system was performed, opening new possibilities to the use of MD in biology. Today, MD simulations serve as a powerful tool for elucidating the atomic-level behavior of proteins. Since McCammon published the results of the first MD simulation of a protein, the computational environment available for such calculations has become more and more unlimited. As mentioned, the first simulation of BPTI was only approximately 10 ps long. As addressed by Karplus and McCammon, during the years that followed this first simulation, a wide range of motional phenomena were investigated through MD simulations of proteins. Last year, the same system (the BPTI protein) was the subject of a 100000000 times longer simulation, that could only be accomplished with the use of a specific-purpose supercomputer. Other long simulations were performed in the past years, with special note to a 10 ms long simulation of a fast-folding protein domain, that needed three months to be finished. Although, in general the average length of a simulation ranges around 10 ns nowadays.


Molecular Dynamics Pipeline

Molecular dynamics (MD) simulations, energy minimization and trajectory analyses were carried out with GROMACS package [51], using GROMOS96 (G53a6) force field [51]. To fully automatize the protocol, a Perl script was designed based on the available version by Marc Offman. The script has several steps, mainly including: 1) Preliminar solvent treatment, where only the so-called structural water molecules (those with conserved positions within the molecules from the crystal and presumably involved in protein function) are maintained in the system; 2) Coordinates file edition, to avoid problems with residue nomenclature and numbering, chain separation, etc; 3) Side-chain treatment, performed by SCWRL to correct andor refine the position of side-chain atoms; 4) The inclusion of solvent and ion molecules (Cl- and Na+) in the system; 5) The minimization of the solvent, with a fix protein structure; 6) The minimization of the protein side-chain, with a fix backbonr; 7) Equilibration phases consisting of short molecular dynamics simulations under constant volume and constant pressure; 8) A production run (usually varying from 2ns to 20ns in length), which is the molecular dynamics simulation itself with all atoms free to move and no volume or pressure constrains. Explicit SPC water molecules [54] were used in all simulations, in which a 14 Å layer of water molecules were added around the solute molecules, within a cubic water box, using periodic boundary conditions. Counter ions were inserted for system neutralization. LINCS [55] and SETTLE [56] were applied to constraint solute and solvent bond, respectively. Temperature was kept at 298 K by rescaling velocities with a stochastic term [57] and pressure at 1 atm using the Berendsen approach [58]. Electrostatic interactions were corrected with PME method [59], using non-bonded cutoffs of 1.0 nm for Coulomb and 1.2 nm for van der Waals. MD integration time was 2 fs.

A 3-steps energy minimization protocol was used to avoid artifacts in atomic trajectories due to conversion of potential into kinetic energies: firstly, applying the steepest-descent algorithm: i. 5000 steps with solute heavy atom positions restrained to their initial positions using an harmonic constant of 1 kJ/mol.nm in each Cartesian direction, allowing free water and hydrogen movements; and ii. 5000 steps with all atoms free to move. Subsequently, the conjugated gradient algorithm was applied for further energy minimization until an energy gradient of 42 KJ/mol.nm. A preliminary MD (1 ns), with heavy atom positions restrained, was performed for achieving solvent equilibration and system heating until 298 K. In this step, the initial velocities were generated once for each simulation. Then we performed a 10 ns production MD.


References

51. Oostenbrink C, Villa A, Mark AE, van Gunsteren WF: A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force-field parameter sets 53A5 and 53A6. Journal of computational chemistry 2004, 25:1656-1676