Fabry:Homology based structure predictions/Journal

From Bioinformatikpedia

Fabry Disease » Homology based structure predictions » Journal

Dataset preparation

The homology search was performed online and resulted in the two output files hhpred.out and coma.out for the structure search with HHpred and COMA, respectively. In both cases, we used the default values, thresholds and databases. From both the resulting files we tried in each case to create three distinct datasets with the demanded sequence identity to the target protein with the following calls and scripts.

$ perl make_dataset_hhpred.pl hhpred.out 0.0000000000000001

This resulted in the HHpred datasets mentioned in Dataset preparation Table 1 and the corresponding pdb structure files.

$ perl make_dataset_coma.pl coma.out 0.002

This resulted in the Coma datasets mentioned in Dataset preparation Table 2 and the corresponding pdb structure files.

Calculation of models


The following steps resulted in 10 models, which can all be downloaded here.

Default settings

For the standard homology modeling with Modeller two basic scripts were used, which build up on the ones described in the recommended tutorial: 1_align.py and 2_Single_template_modeling.py (in this case example file for 1ktb) which need the appropriate template pdb structure files, as well as the target (AGAL) sequence in pir format as input:

sequence:1R46:::::::0.00: 0.00

The first script creates an alignment only based on the structure of template and target, as well as a 2d-alignment also based on the given structural data. The second script creates the actual model and assesses it based on the DOPE score and the GA341 method (see Melo et al., 2002 and John & Šali, 2003). The following runs were performed on our home computers:

#3HG3 100%
$ mod9.10 1_align_3hg3.py
$ mod9.10 2_Single_template_modeling_3hg3.py

#1KTB 53%
$ mod9.10 1_align.py
$ mod9.10 2_Single_template_modeling.py

#3CC1 25%
$ mod9.10 1_align_3cc1.py
$ mod9.10 2_Single_template_modeling_3cc1.py

The scores were read from the log files with the perl script read_logfiles_modeller.pl.

$ perl read_logfiles_modeller.pl 1_align_3hg3.log 2_Single_template_modeling_3hg3.log
$ perl read_logfiles_modeller.pl 1_align.log 2_Single_template_modeling.log 
$ perl read_logfiles_modeller.pl 1_align_3cc1.log 2_Single_template_modeling_3cc1.log

The vizualisation was done with the help of Pymol and the commands listed in pymol_single.txt

Multiple templates

For multiple templates three scripts are needed, which again base on the tutorial and are here again only exemplified for one case: 3_multiAli_1.py, 4_add1R46_1.py and 5_multipleMod_1.py. In the last script we added the methods for the assessment scores, since otherwise there would have been no quality scores for comparison. The following combinations were performed:

#Multiple 3HG3,1KTB
$ mod9.10 3_multiAli_1.py
$ mod9.10 4_add1R46_1.py
$ mod9.10 5_multipleMod_1.py

#Multiple 3HG3,1KTB,3CC1
$ mod9.10 3_multiAli_2.py
$ mod9.10 4_add1R46_2.py
$ mod9.10 5_multipleMod_2.py

#Multiple 3CC1, 3ZSS, 3A24
$ mod9.10 3_multiAli_3.py
$ mod9.10 4_add1R46_3.py
$ mod9.10 5_multipleMod_3.py

#Multiple 3CC1, 3HG3
$ mod9.10 3_multiAli_4.py
$ mod9.10 4_add1R46_4.py
$ mod9.10 5_multipleMod_4.py

The scores were read from the log file with the perl script read_logfiles_modeller_CHAS.pl.

$ perl read_logfiles_modeller_CHAS.pl 5_multipleMod_1.log MULTI1
$ perl read_logfiles_modeller_CHAS.pl 5_multipleMod_2.log MULTI2
$ perl read_logfiles_modeller_CHAS.pl 5_multipleMod_3.log MULTI3
$ perl read_logfiles_modeller_CHAS.pl 5_multipleMod_4.log MULTI4

The vizualisation was done with the help of Pymol and the commands listed in pymol_multi.txt

Edited Alignment input

For the edited alignment models, only the second Modeller script 2_Single_template_modeling.py was needed. Edited alignments were provided as input. We rearranged them with the program SeaView.

#Active Site shifted right to next D (7 and 1) in -2d.ali
$ mod9.10 2_Single_template_modeling_3hg3_changed_actSite.py

#Active Site shifted right to next D (7 and 1) in both ali files
$ mod9.10 2_Single_template_modeling_3hg3_changed_actSite_2.py

#Active Site shifted right to next D (7 and 1) in both ali files + Substrate binding region (203-207) forced to be consecutive
$ mod9.10 2_Single_template_modeling_3hg3_changed_actSite_3.py

The scores were read from the log file with the perl script read_logfiles_modeller_CHAS.pl.

perl read_logfiles_modeller_CHAS.pl 2_Single_template_modeling_3hg3_changed_actSite.log CHAS
perl read_logfiles_modeller_CHAS.pl 2_Single_template_modeling_3hg3_changed_actSite_2.log CHAS2
perl read_logfiles_modeller_CHAS.pl 2_Single_template_modeling_3hg3_changed_actSite_3.log CHAS3

The vizualisation was done with the help of Pymol and the commands listed in pymol_changed.txt


The SwissModel calculations were performed online.
Pictures were created either with Pymol (see pymol_swiss.txt) or with R:

$ R CMD BATCH localError.R


As input sequence, the already known Alpha-galactosidase sequence was used for I-Tasser.

The models were calculated using Option II: Exclude some templates from I-TASSER template library. to exclude homologous templates. The field Exclude homologous templates was

  • left blank
  • set to 80%
  • set to 30%
  • set to 20%

in the different runs. This results in up to five models for each run, which were downloaded and stored in folders named AGAL_HUMAN*. Eventually, the script calc_itasser_scores.sh was applied on the models to calculate the different measures using the TMscore and sap binaries. It also invokes the scripts tm2wikirow.pl and sap2wikirow.pl to parse the output of these tools.

$ bash calc_itasser_scores.sh


As input sequence, the already known Alpha-galactosidase sequence was used for 3D-Jigsaw and the upload mode was selected. The input model sets are listed in detail at the 3D-Jigsaw section. These sets were generated by the following scripts which make use of repairPDB to fix the input files:

$ bash gen_top5_30.sh
$ bash gen_top5_40_80.sh
$ bash gen_top5_80.sh

When the jobs were done, they were downloaded with

$ bash fetch_jigsaw_results.sh

and the RMSD, TM-Score, etc. were calculated:

$ bash calc_3djigsaw_scores.sh



The scores were computed on a home computer with the TM-score version of 2012/05/07. We used a bash script, which calls the TM-score program for all models, as well as the perl script read_TMoutput.pl that extracts the calculated scores and automatically outputs them in Media-wiki table format into a file. The bash script requires the template pdb file and the directory in which all models are located as input.

$ ./calculate_TMScore.sh ../Modeller/pdb/1R46.pdb ../Modeller/final_models/
$ ./calculate_TMScore.sh ../Modeller/pdb/1R47.pdb ../Modeller/final_models/

$ ./calculate_TMScore.sh ../Modeller/pdb/1R46.pdb ../Swissmodel/final_models/
$ ./calculate_TMScore.sh ../Modeller/pdb/1R47.pdb ../Swissmodel/final_models/


Due to some server errors, we calculated most of the RMSD values online with the SAP Web Tool. From this we obtained the same output as from the command line tool. The output was read by another pair of bash and perl script, to extract the calculated scores and create a table.

$ ./readRMSD.sh

The pngs used for the animated gifs were created with Pymol and the commands listed in pdb_superimpose.txt. Afterwards the pngs were animated with the program GIF Movie Gear.

RMSD around catalytic site

The all atom RMSD of the atoms within a radius of 6 Angstrom around the catalytic centre was calculated with the help of pymol and the commands, that can be found in [ RMSD_6A.txt].


The DOPE profile was read with another script, based on a tutorial on the MODELLER website. It reads the DOPE score for each residue and writes them into a file. These scores can be plotted with the R script Profile.R.

$ mod9.10 6_DOPE_score.py
$ R CMD BATCH Profile.R 


Visual comparisons were performed with Pymol and the commands listed in comparison.txt