Homology modelling Gaucher Disease

The object of this task was to apply homology modeling for predicting the tertiary structure of glycosylceramidase given its sequence. For this, we first selected different templates which were than used to derive the structure of glycosylceramidase using three different homology modeling tools, namely Modeller<ref name="modeller">A. Sali, T.L. Blundell. (1993) Comparative protein modelling by satisfaction of spatial restraints. J. Mol. Biol.</ref> , SWISS-MODEL <ref name="swissmodel">Arnold K., Bordoli L., Kopp J., and Schwede T. (2006). The SWISS-MODEL Workspace: A web-based environment for protein structure homology modelling. Bioinformatics</ref>, and the I-TASSER<ref name="itasser">Ambrish Roy, Alper Kucukural, Yang Zhang. I-TASSER: a unified platform for automated protein structure and function prediction. Nature Protocols</ref> server. The resulting models were evaluated using both quality assessment scores and the native crystal structure 1ogs for comparison. Technical details are reported in our protocol.

Template selection

We used HHsearch for searching the PDB for homologous templates. <xr id="tab:templates"/> lists some of the top-ranking templates. 2nt0_A is identical to the target 1ogs_A and was therefore excluded. Although all listed hits are homologous to the target (HHsearch probability > 97%), their sequence identity was below 30%. We therefore selected 2wnw_A (blue) as a close homolog, 2y24_A (green) as an intermediate homolog, and 3nco_A (yellow) as a more distant homolog. Note that the latter two templates to not cover the complete target which makes the homology modeling process harder. <figtable id="tab:templates">

Hit Nr	Template	Identity	Query HMM	Prob
> 80% sequence identity
1	2nt0_A	100.0	1-497	100.0
40% - 80% sequence identity
< 30% sequence identity
2	2wnw_A	28.0	36-496	100.0
3	3clw_A	14.0	64-495	100.0
4	2y24_A	18.0	66-495	100.0
5	3kl0_A	19.0	65-495	100.0
6	3zr5_A	17.0	65-494	100.0
7	3ik2_A	14.0	65-495	99.2
22	3nco_A	11.0	113-384	97.7
28	1egz_A	12.0	113-387	97.4

Homologs found by HHsearch. Bold: selected templated used for the following modeling. </figtable>

Modeller

Modeller is a popular tool for building models by the satisfaction of spatial restrains which are derived, for instance, from one or several target-template alignments. The alignments can come from any alignment or homology search tool, or they can be built by modeller itself. We created models for our target protein by (1) using a single template and employing Modeller to compute the alignment, (2) using a single template and the alignment from HHsearch, and (3) using multiple templates. The Model quality was assessed via the DOPE and DOPE z-score reported by Modeller as well as the QMEAN6 score from the SWISS-MODEL workspace. We compared the resulting models to the crystal structure 1ogs_A via the weighted all-atom RMSD score computed by SAP, as well as the TM-score, GDT_TS, and GTD_HA score which we computed by the program TMscore.

Single-template modeling using Modeller alignments

<xr id="tab:modeller-single-models"/> shows the resulting models. 2wnw_A produced the best looking model: all major secondary structure elements coincided with the native structure 1ogs_A (red). Only the target range 1-31 which was not covered by template (cf. <xr id="tab:templates"/>) resulted in some deviating loop regions (at the top right corner). Although 2y24_A is less conserved than 2wnw_A, the model came close to the native structure. 3nco_A shares the same TIM beta/alpha-barrel domain (the tube at the center) than 1ogs_A but is missing the glycosyl hydrolase domain (the sheets at the right side) such that the model was less well structured in this region.

<figtable id="tab:modeller-single-models">

2wnw_A

2y24_A

3nco_A

Dope score per residue.

Models built by Modeller using single templates and alignments computed by Modeller itself. Red: 1ogs_A. </figtable>

<xr id="tab:modeller-single-eval"/> shows the evaluation results. Both the quality assessment scores and the comparison to the native structure 1ogs_A were in line with the observations described above. The model derived from 2wnw_A had a lower energy (was more stable) than the models from 2y24_A and 3nco_A, and fitted better with 1ogs_A. Except for the DOPE z-score which should have been lower for 2y24_A, the assessment scores correlated well with the structure comparison scores. The dope score per residue (cf. <xr id="tab:modeller-single-models"/>) of all three models were correlated and lowest for 2wnw_A.

<figtable id="tab:modeller-single-eval">

Template	DOPE	DOPE z-score	QMEAN6	RMSD	TM-score	GDT_TS	GTD_HA
2wnw_A	-55925	-0.471	0.689	1.006	0.824	0.661	0.479
2y24_A	-47194	0.777	0.376	1.222	0.550	0.294	0.223
3nco_A	-44033	1.224	0.139	2.158	0.252	0.093	0.043

Evaluation of models built by Modeller using single templates and alignments computed by Modeller itself. </figtable>

Single-template modeling using HHsearch alignments

Modeller computes alignments by aligning the target sequence to the known structure of the template. Hence, no predicted features of the target sequence are used. Instead, HHsearch computes an HMM-to-HMM alignment where the target HMM comprises more features than the sequence alone. It also contains the secondary structure predicted by PSIPRED and information about the conservation of all residues derived from a sequence profile. HHsearch alignments are therefore thought to be more accurate than those produced by Modeller which should lead to better models. Thus, we tried to improve the models by using HHsearch alignments.

<figtable id="tab:modeller-hhsearch-alis">

2wnw_A

2y24_A

3nco_A

Comparison of Modeller and HHsearch alignments. </figtable>

<xr id="tab:modeller-hhsearch-alis"/> depicts HHsearch alignments compared to Modeller alignments. The alignments produced by HHsearch were more compact, i.e. they exhibited more gaps at the beginning and at the end of the aligned sequences, whereas gaps were more distributed in case of alignments computed by modeller. In comparing the resulting models of <xr id="tab:modeller-hhsearch-models"/> with those of <xr id="tab:modeller-single-models"/>, we found that the core regions better matches the native structure in case of using HHsearch alignments. However, Modeller could not find the correct topology for the ends of the target sequence which were not covered by the template. These regions just became an unfolded thread pointing away from the core. Surprisingly, their dope score per residue was relatively low.

<figtable id="tab:modeller-hhsearch-models">

2wnw_A

2y24_A

3nco_A

Dope score per residue.

Models built by Modeller using HHsearch alignments covering the complete target. Red: 1ogs_A. </figtable>

The DOPE score of the models shown in <xr id="tab:modeller-hhsearch-models"/> was slightly higher compared to the models shown in <xr id="tab:modeller-single-models"/>, but they were more stable according to the QMEAN6 score. The RMSD score increased due to the unfolded ends of the target sequence. However, the TM-score and the GDT scores improved since the actual core to the target became closer to the native structure.

<figtable id="tab:modeller-hhsearch-eval">

Template	DOPE	DOPE z-score	QMEAN6	RMSD	TM-score	GDT_TS	GTD_HA
2wnw_A	-54695	-0.295	0.726	1.079	0.869	0.732	0.538
2y24_A	-47256	0.765	0.566	1.640	0.722	0.553	0.386
3nco_A	-32577	2.857	0.272	9.757	0.398	0.235	0.131

Evaluation of models built by Modeller using HHsearch alignments covering the complete target. </figtable>

Single-template modeling using local HHsearch alignments

Since the ends of the sequences could not be modelled correctly by Modeller and might impact the model evaluation, we repeated the modeling using the same alignments as in the previous section but with the uncovered ends truncated. This resulted in more dense models without the long unfolded threads (cf. <xr id="tab:modeller-hhsearch-local-models"/>).

<figtable id="tab:modeller-hhsearch-local-models">

2wnw_A

2y24_A

3nco_A

Dope score per residue.

Models built by Modeller using local HHsearch alignments covering the complete target. Red: 1ogs_A. </figtable>

Although the DOPE score became worse since the target sequence was shorter, the QMEAN6 score could be further increased. The RMSD could be reduced, in particular in case of 3nco_A, whereas the TM-score and GDT score remained the same. This demonstrated the sensitivity of the RMSD for deviation loop regions whereas the TM-score and GDT score better evaluate the actual fit of the core region.

<figtable id="tab:modeller-hhsearch-local-eval">

Template	DOPE	DOPE z-score	QMEAN6	RMSD	TM-score	GDT_TS	GTD_HA
2wnw_A	-53593	-0.744	0.749	1.003	0.866	0.731	0.532
2y24_A	-46290	-0.030	0.564	1.162	0.724	0.563	0.381
3nco_A	-25047	0.886	0.451	1.827	0.402	0.245	0.138

Evaluation of models built by Modeller using HHsearch alignments covering the complete target. </figtable>

Multiple-template modeling using Modeller alignments

For testing whether multiple templates could improve the model quality, we prepared three groups of templates: (1) 2wnw_A, 2y24_A, and 3kl0_A as close homologs, (2) 2wnw_A and 3nco_A as a close and a more distant homolog, and (3) 1ogs_A, 3ik2_A, and 3nco_A-1egz as remote homologs. A multiple alignment was created for each group via Modeller to which the target was added afterwards. The resulting alignment was used as input for Modeller.

<figtable id="tab:modeller-mult-models">

2wnw_A-2y24_A-3kl0_A

2wnw_A-3nco_A

3ik2_A-3nco_A-1egz_A

Dope score per residue.

Models built by Modeller using multiple templates. Red: 1ogs_A. </figtable>

<xr id="tab:modeller-mult-models"/> shows the resulting models. None of the template combinations led to better model than using 2wnw_A alone (cf. <xr id="tab:modeller-single-models"/>). We assumed that 2wnw_A was more similar to 1ogs_A than 2y24_A and 3kl0_A were. The latter two templates rather drew the model away from the true structure 1ogs_A and did not contain additional information which might have helped building a better model. Surprisingly, using 2wnw_A in combination with 3nco_A resulted in a better model than using two templates with a higher sequence identity. Taking advantage of more than one template improves the model quality in particular if the templates cover different region, e.g. different domains, of the target sequence. However, this had not been the case for our target since 2wnw_A had already the highest coverage (cf. <xr id="tab:templates"/>).

<figtable id="tab:modeller-mult-eval">

Template	DOPE	DOPE z-score	QMEAN6	RMSD	TM-score	GDT_TS	GTD_HA
2wnw_A-2y24_A-3kl0	-39084	1.930	0.314	1.987	0.465	0.277	0.179
2wnw_A-3nco_A	-46962	0.807	0.556	1.256	0.778	0.593	0.409
3ik2_A-3nco_A-1egz_A	-25391	0.3881	0.101	6.607	0.241	0.062	0.029

Evaluation of models built by Modeller using multiple templates. </figtable>

SWISS-MODEL

The SWISS-MODEL webserver provides two modes for carrying out the modeling: given the target sequence, the model is built completely automatically without user interference in case of the automated mode whereas the user can specify the target-template alignment by its own in case of the alignment mode.

Automated mode

For testing the automated mode, we dictated the SWISS-MODEL server to use the selected templates 2wnw_A, 2y24_A, and 3nco_A, instead of identifying the templates by its own. In this way, we could better compare the methods for aligning the target to the template and building the model with the methods of the modeller package.

<figtable id="tab:swiss-auto-models">

2wnw_A

2y24_A

3nco_A

Dope score per residue.

Models built by SWISS-MODEL using the automated mode. Red: 1ogs_A. </figtable>

<xr id="tab:swiss-auto-models"/> shows the models build by the SWISS-MODEL server. Since the SWISS-MODEL server did not model starting and trailing ends of the target sequence which were not covered by the target, the results have to be compared to #Single-template modeling using local HHsearch alignments. Altogether, quality of the models whose comparable to Modeller (cf. <xr id="tab:swiss-auto-eval"/>, <xr id="tab:modeller-hhsearch-local-eval"/>. Hence, the SWISS-MODEL alignments were in our case approximately as accurate as HHsearch alignments and the fragment assembly procedure worked as reliably as Modeller's satisfaction of spatial restrains.

<figtable id="tab:swiss-auto-eval">

Template	DOPE	DOPE z-score	QMEAN6	RMSD	TM-score	GDT_TS	GTD_HA
2wnw_A	-53863	-0.786	0.530	1.003	0.858	0.723	0.540
2y24_A	-42674	0.565	0.430	1.100	0.712	0.539	0.382
3nco_A	-27101	2.109	0.270	2.158	0.422	0.252	0.138

</figtable>

Alignment mode

For testing the alignment mode, we used the same HHsearch alignments as we did for evaluating Modeller (cf. #Single-template modeling using HHsearch alignments). Modeller failed to built a model in case of 2y24_A an complained about alignment position 84, which was a gap. After removing this alignment column, the SWISS-MODEL server succeeded.

<figtable id="tab:swiss-ali-models">

2wnw_A

2y24_A

3nco_A

Dope score per residue.

Models built by SWISS-MODEL using the alignment mode. Red: 1ogs_A. </figtable>

Except for 2y24_A, the resulting models were only slightly of minor quality compared to the automated mode. This reconfirmed our observation that the SWISS-MODEL alignments were as good as HHsearch alignments. The model for 2y24_A was significantly worse which was probably caused be the deletion of alignment column 84.

<figtable id="tab:swiss-ali-eval">

Template	DOPE	DOPE z-score	QMEAN6	RMSD	TM-score	GDT_TS	GTD_HA
2wnw_A	-53963	-0.788	0.540	1.011	0.859	0.722	0.536
2y24_A	-40131	1.055	0.360	1.222	0.510	0.238	0.067
3nco_A	-30244	0.689	0.240	2.001	0.426	0.247	0.135

</figtable>

I-TASSER

We evaluated I-TASSER in two ways: To make the evaluation consistent with the other methods, we first tried to assign the templates 2wnw_A, 2y24_A, and 3nco_A using Option I: Specify template without alignment. Yet, the I-TASSER server used in each case 2nt0_A (cf. <xr id="tab:templates"/>) for the modeling such that we repeated the modeling using Option II: Exclude homologous templates and assigned my maximum sequence identity of 60%.

Specifying templates without alignment

For being able to directly compare the prediction performance of I-TASSER with the other servers, we specified the selected templates 2wnw_A, 2y24_A, and 3nco_A as templates. The modelling took about 72 hours, compared to less than 5 minutes in case of Modeller and the SWISS-MODEL server.

<figtable id="tab:itasser-models">

2wnw_A C-score: 1.47

2y24_A C-score: 0.12

3nco_A C-score: -0.25

Dope score per residue.

Models built by I-TASSER using different templates. Red: 1ogs_A. </figtable>

<figure id="fig:itasser-2wnw_A-loops">

Loop region 1-35 of 2wnw_A

</figure>

The resulting models turned out to be highly precise (cf. <xr id="tab:itasser-models"/>, <xr id="tab:itasser-eval"/>). Although the DOPE score and RMSD was not better compared to Modeller and the SWISS-MODEL server, the TM-score and GDT score was. In contrast to the other methods, I-TASSER also modelled the loop regions very accurate. <xr id="fig:itasser-2wnw_A-loops"/> demonstrates the the high coincidence of the loop region 1-35 with 1ogs_A. Also the side-chains matched better those of 1ogs_A than the models built by the other methods. This also accounted for the high TM-score and GDT scores. However, we found that I-TASSER used in each case also 2nt0_A (cf. <xr id="tab:templates"/> as template which is identical to the target 1ogs_A. This makes the modeling of course trivial.

<figtable id="tab:itasser-eval">

Template	DOPE	DOPE z-score	QMEAN6	RMSD	TM-score	GDT_TS	GTD_HA
2wnw_A	-54684	-0.293	0.73	1.010	0.943	0.812	0.607
2y24_A	-52919	-0.042	0.603	1.332	0.873	0.673	0.460
3nco_A	-50315	0.329	0.568	2.298	0.785	0.497	0.286

Evaluation of models built by I-TASSER using different templates. </figtable>

Specifying a maximum sequence identity of 60%

In order to prevent I-TASSER from using 2nt0_A as template, we specified a maximum target-template sequence identity of 60%. As a result, 2nt0_A was excluded and I-TASSER correctly identified 2wnw_A with the second highest sequence identity of 24% as template.

<figtable id="tab:itasser-id-models">

2wnw_A C-score: -0.52

Dope score per residue.

Models built by I-TASSER using templates with a target-template sequence identity <60%. Red: 1ogs_A. </figtable>

Despite a remarkable higher computation time of more than 72 hours, I-TASSER could only build a slightly model than Modeller and the SWISS-MODEL server (cf. <xr id="tab:itasser-id-models"/>, <xr id="tab:itasser-id-eval"/>), if the TMscore and GDT score is used for comparison. This showed that the alignment produced by the threading tool MUSTER, which is employed by I-TASSER, was as accurate the the alignment of Modeller and HHsearch. Furthermore, I-TASSER's fragment assembly and loop modeling algorithms did not outperform Modeller and the SWISS-MODEL server.

<figtable id="tab:itasser-id-eval">

Template	DOPE	DOPE z-score	QMEAN6	RMSD	TM-score	GDT_TS	GTD_HA
2wnw_A	-54661	-0.290	0.693	1.062	0.882	0.737	0.538

Evaluation of models built by I-TASSER using templates with a target-template sequence identity <60%. </figtable>

3D-JIGSAW

3D-JIGSAW<ref name="jigsaw">Contreras-Moreira B, Fitzjohn PW, Offman M, Smith GR, Bates PA. (2003) Novel use of a genetic algorithm for protein structure prediction: searching template and sequence alignment space. Proteins.</ref> employs a genetic algorithm for finding a conformation which minimizes the POPLUS energy function: several input models are recombined to obtain new models which a selected depending on their energy. The optimization algorithm terminates after a certain number of iterations. For testing whether 3D-JIGSAW could build a better model from a set of input models, we called 3D-JIGSAW with our best Model obtained by Modeller and the I-TASSER server. Since the input models had to be of the same size as the input sequence, we could not use a model from the SWISS-MODEL server which did not cover the complete target sequence.

<figtable id="tab:jigsaw-models">

2wnw_A

Dope score per residue.

Models built by 3D-JIGSAW using the best model of Modeller and I-TASSER as input. Red: 1ogs_A. </figtable>

The results are shown in <xr id="tab:jigsaw-models"/> and <xr id="tab:itasser-id-eval"/>. Noteworthy was both the low DOPE score and the high QMEAN6 score which suggested that 3D-JIGSAW could in fact find a conformation with a lower energy compared to the two input models. However, the model did not come closer to the native structure 1ogs_A. Due to the high coincide of the input models with 1ogs_A, this might be explained by over-fitting, i.e. 1ogs_A was not the conformation with the lowest energy.

<figtable id="tab:jigsaw-eval">

Template	DOPE	DOPE z-score	QMEAN6	RMSD	TM-score	GDT_TS	GTD_HA
2wnw_A	-58732	-0.870	0.737	1.052	0.871	0.726	0.536

Evaluation of models built by 3D-JIGSAW. </figtable>

Conclusions

• Modeller, SWISS-MODEL, and I-TASSER produced models of similar quality.
• Apart from the DOPE z-score, the model assessment scores were highly correlated with the structure comparison scores.
• Modeller is a very versatile package providing plenty of modelling features for different purposes.
• However, deriving spatial restrains from multiple templates made our model worse. Better results might be achieved by weighting the templates differently, depending on how good the alignment is at a particular region.
• HHsearch alignment improved the model quality in regions which were covered by the template.
• The SWISS-MODEL fragment assembly method worked as well the Modeller's satisfaction of spatial restraints.
• I-TASSER built slightly more accurate models at the expense of a much longer computation time.
• 3D-JIGSAW found a model with a lower energy which, however, deviated more from 1ogs_A than the best I-TASSER model.

Segment matching vs. fragment assembly

• Segment matching is employed by:
  • SMM
• Fragment assembly is employed by:
  • Rosetta
  • SWISS-MODEL
  • I-TASSER

The segment matching approach by Levitt breaks the complete target sequence randomly into short pieces which are matched against a database of segments with non structure. The best matching segments are assembled to get the structure of the protein. This procedure is repeated several times and the best model is selected in the end.

Fragment assembly methods first search the PDB for templates which cover certain fragments of the target. The templates structure is assigned to the matching fragments. Parts of the target which are not covered by any template are solved in a subsequent loop modeling step. The loops are either modeled de novo using some force field, or knowledge based by searching a database of experimentally determined loops which link the fragments in a suitable way. De novo and knowledge-based loop moding approaches are often combined.

References

<references/>