Computational Modeling of complete HOXB13 protein for predicting the functional effect of SNPs and the associated role in hereditary prostate cancer

The human HOXB13 gene encodes a 284 amino acid transcription factor belonging to the homeobox gene family containing a homeobox and a HoxA13 N-terminal domain. It is highly linked to hereditary prostate cancer, the majority of which is manifested as a result of a Single Nucleotide Polymorphism (SNP). In silico analysis of 95 missense SNP’s corresponding to the non-homeobox region of HOXB13 predicted 21 nsSNP’s to be potentially deleterious. Among 123 UTR SNPs analysed by UTRScan, rs543028086, rs550968159, rs563065128 were found to affect the UNR_BS, GY-BOX and MBE UTR signals, respectively. Subsequent analysis by PolymiRTS revealed 23 UTR SNPs altering the miRNA binding site. The complete HOXB13_M26 protein structure was modelled using MODELLER v9.17. Computational analysis of the 21 nsSNP’s mapped into the HOXB13_M26 protein revealed seven nsSNP’s (rs761914407, rs8556, rs138213197, rs772962401, rs778843798, rs770620686 and rs587780165) seriously resulting in a damaging and deleterious effect on the protein. G84E, G135E, and A128V resulted in increased, while, R215C, C66R, Y80C and S122R resulted in decreased protein stability, ultimately predicted to result in the altered binding patterns of HOXB13. While the genotype-phenotype based effects of nsSNP’s were assessed, the exact biological and biochemical mechanism driven by the above predicted SNPs still needs to be extensively evaluated by in vivo and GWAS studies.

Functionally deleterious nsSNP's predicted by SIFT program (Sequence-based homology). From a total of 95 non-homeobox region nsSNP's, 39 nsSNP's were predicted to be functionally deleterious and were marked as "Affects protein structure" by the SIFT server (Table 1). Among those, 18 non-homeobox region nsSNP's were marked as deleterious with a tolerance index of 0.00. The remaining 56 variants were predicted as "Tolerated" by the SIFT program. The detailed analysis of the effect of nsSNP's on the entire non-homeobox region (1-215 and 276-284 AA) of HOXB13 by the SIFT program can be found as Supplementary Fig. S3. The statistical representation of the results is given in Fig. 4. The complete SIFT prediction results can be found as Supplementary Table S1. The functionally damaging nsSNP's predicted by PolyPhen version 2 server (Structure-based homology). The PolyPhen server predicted 47 non-homeobox region nsSNP's to be functionally deleterious to the protein structure. Out of those 47 nsSNP's, 34 nsSNP's were predicted to be "probably damaging" with the score ranging from 0.845 to 1.00 and the remaining 13 nsSNP's were predicted to be "possibly damaging" with the score ranging from 0.537 to 0.851 (Table 2). Interestingly, 22 nsSNP's predicted as deleterious by SIFT were also predicted to be functionally damaging by the PolyPhen server. Therefore, the 22nsSNP's predicted commonly by SIFT and PolyPhen were of functional importance. The statistical representation of the results are given in Fig. 4. The complete PolyPhen prediction results are available as Supplementary Table S2.
The functional validation of deleterious nsSNP's by the PANTHER server. The PANTHER server predicted 25 non-homeobox region nsSNP's to be damaging and the remaining nsSNP's were predicted to be benign. Interestingly, 19 nsSNP's were predicted as deleterious in common among SIFT, PolyPhen, and PANTHER server (Table 3). Additionally, two nsSNP's (G153S, L152M) predicted by PANTHER and  PolyPhen, and one nsSNP (H30Q) predicted by PANTHER and SIFT, were found to be common. The graphical representation of the results are given in Fig. 4. The complete PANTHER prediction results are available as Supplementary Table S3.
The functional impact of deleterious nsSNP's by the PROVEAN server. The PROVEAN server predicted 20 non-homeobox region nsSNP's to be functionally damaging out of the 95 nsSNP's submitted for analysis (Table 4). Among those, 16 nsSNP's were found to be in common as predicted by SIFT, PolyPhen, and PANTHER servers. One nsSNP (D167N) was found to be in common to both predicted by SIFT and PolyPhen. Two nsSNP's (G117E and G117R) were found to be in common with the nsSNP's predicted by SIFT program. The graphical representation of the results are given Fig. 4. The complete PROVEAN prediction results are available as Supplementary Table S4.
The functional impact of deleterious nsSNP's by the nsSNPAnalyzer server. Out of the 95 nsSNP's submitted for analysis, 51 nsSNP's were predicted to be associated with a diseased phenotype. Among those, 20 nsSNP's were common to those predicted by the above four servers (SIFT, PolyPhen, PANTHER, and PROVEAN) ( Table 5). The graphical representation of the results are given in Fig. 4 The functional impact of deleterious nsSNP's by the PhD-SNP server. We used the SVM based method utilizing sequence and profile information algorithm for the analysis of 95 non-homeobox region nsSNP's of the human HOXB13 gene. The server predicted 13 nsSNP's ( Table 6) (Fig. 4) to be functionally associated with the disease and the remaining were considered benign. Among those, ten nsSNP's were common to those predicted by the above-described servers (SIFT, PolyPhen, PANTHER, PROVEAN and nsSNPAnalyzer). The graphical representation of the results are given Fig. 4. The complete PhD-SNP prediction results are available as Supplementary Table S6. Among the 95 HOXB13 non-homeobox region nsSNP's subjected to analysis by SIFT, PolyPhen, PANTHER, PROVEAN, nsSNPAnalyzer and PhD-SNP servers, 21 nsSNP's were found to be functionally significant and causing damaging effects to the HOXB13 protein structure, stability, and function by the servers mentioned above. The list of those 21 nsSNP's are as follows: rs761914407 (R215C), rs779330626 (Q188R), rs570681642 (Q181R), rs777986934 (G177V), rs587780164 (D167N), rs766929278 (G153V), rs575899185 (Q138H), rs770891609 (Y137S), rs769634543 (G135E), rs775273363 (A128V), rs201428095 (R123H), rs8556 (S122R), rs138213197 (G84E), rs772962401 (Y80C), rs778843798 (C66R), rs199813155 (C63Y), rs758166293 (C63G), rs568967699 (K61M), rs770620686 (P59L), rs561048036 (H30Q), rs587780165 (R25Q). For subsequent analysis, these 21 nsSNP's were taken into consideration. The HOXB13 protein structure was available only for the homeobox binding domain (216-275 AA) of the complete protein and since there was no complete structure of HOXB13 in the PDB, the Homology Modelling approach was adopted for simulating the complete protein structure of HOXB13 in silico, so that we could map the above screened 21 nsSNP's into the protein structure and could predict their effects on the protein function, stability, and bioactivity. Functionally significant HOXB13 UTR SNPs predicted by the UTRscan server. Mutations in the untranslated region of the gene were reported very often to be linked with hereditary diseases such as cancer and various immune deficiency syndromes and also plays a key role in mRNA localization, stability, and translational efficiency 45 . Both the 5′ UTR and the 3′ UTR have important functions concerning the stability and expression of the mature mRNA. Mutations in those regions are linked with severe effects on the expression patterns of the gene at the level of mNA processing and translation 46 . The polymorphisms in the 5′ UTR are increasingly related to the altered patterns of ribosomal binding capacity, stability and translational regulation of mRNA, thereby influencing the RNA half-life. Whereas the polymorphisms corresponding to the 3′ UTR are highly involved in altered patterns of polyadenylation, localization, stability, translational efficiency and microRNA (miRNA) binding specificity, thereby rendering a tremendous effect on the gene expression patterns.
The UTRscan server predicts both the 5′ UTR and 3′ UTR SNPs. Among a total of 101 valid 3′ UTR SNPs taken for evaluation, the UTRscan server predicted three SNPs (rs543028086, rs550968159, rs563065128) to be functionally significant to cause a pattern change (Table 7). However, the UTRscan server did not predict any harmful 5′ UTR SNPs. UNR (Upstream of N-ras) is a transcription factor containing five cold shock domains (CSD) that bind to single-stranded DNA and RNA 47 . It controls and plays a major role in transcriptional and post-transcriptional gene expression. UNR is a cytoplasmic protein known to function as an RNA chaperone and is found to be crucial in the control of cell proliferation and death 48,49 . The protein mainly destabilizes the c-fos mRNA and helps in the initiation and activation of cap-independent translation via the IRES for various transcripts, especially the proto-oncogene c-myc, rhinovirus, poliovirus, the cell cycle PISTLRE kinase, and pro-apoptotic factor (Apaf-1) 50,51 . The SNP rs543028086 was predicted to result in the disruption of the UNR Binding Site (UNR_BS) motif in the 3′ UTR of human HOXB13 gene, which results in the deregulation of pro-apoptotic factor (Apaf-1), which might have a negative effect on the control of cell death.     The GY-Box is a conserved motif present in the Notch pathway target genes in Drosophila 52 . It is highly conserved in 3′ UTR regions that have sequence complementarity to the 5′ regions of the miRNA seed region. The result is the formation of RNA duplexes by the interaction of the 3′ UTR end of mRNA and the 5′ end of the miRNA, leading to translational repression 52 . The SNP rs550968159 is present in this region of the 3′ UTR of human HOXB13 gene and leads to the loss of the specific GY-Box pattern, hence voiding the chance of translational repression of HOXB13 mediated by the GY-Box in the 3′ UT region. Thus, an imbalance in the feedback regulation of HOXB13 expression has been predicted to result in the diseased state.
The Mushashi Binding Element (MBE) is an mRNA binding protein, which plays a very important role in the regulation of stem cell renewal process 53,54 by suppressing the translation of all the mRNA coding for the proteins involved in inhibiting cell cycle progression 55,56 . The 3′ UTR SNP rs563065128 results in the loss of the MBE UTR motif in the human HOXB13 gene and is thereby found to lose its natural role of regulating the stem cell renewal process by suppressing the expression of cell cycle progression inhibitors, thereby predicted to result in the loss of stem cell niche. Loss of stem cell niche, in turn, leads to the unavailability of the local stem cell source to replenish the damaged cells of the tissue, thus, leading to the disease state.
These three 3′ UTR SNPs were predicted to have important deleterious effects and functional significance on the expression of human HOXB13 gene. Functionally significant HOXB13 3′ UTR SNPs predicted by the PolymiRTS Database. Out of the 95 nsSNP's under consideration, only 23 nsSNP's were found to have a crucial role in the 3′ UTR region (see Supplementary Table S7). Among those, five nsSNP's (rs8064432, rs79812861, rs148791210, rs184053751 and rs183620920) were found to disrupt only the conserved (ancestral allele with support > = 2) miRNA sites. Two nsSNP's (rs116931900 and rs1056656) were exclusively found to create a new miRNA site. Whereas, the remaining 16 nsSNP's were predicted to be involved in the disruption and creation of a new miRNA site, out of which rs61123825 (disrupting -2 and creating -7) and rs192244427 (disrupting -4 and creating 5) were found to have a maximum of 9 pattern changes.
Modelling of the complete HOXB13 protein using MODELLER v9.17 (Comparative Homology Modelling). The PDB contains only the 3D structure of the Homeobox binding domain of the human HOXB13 protein and not the complete 3D structure. In order to further analyze the effect of the above shortlisted 21 non-homeobox region missense SNPs on the HOXB13 protein structure and function, the complete protein structure was mandatory, since further analysis of those nsSNP's demands mapping the nsSNP's into the protein structure and thereby validating their subsequent effects on structural and functional aspects of the protein in silico. Thus, the complete HOXB13 protein structured was modelled by a technique called comparative homology modelling using the MODELLER v9.17 tool from the Andrej Sali laboratory. The complete modelling procedure and the steps performed were mentioned in the Supplementary Material.
A suitable template structure for developing the model was obtained using psiBLAST by setting PDB as the source database 57 for finding the 3D structure templates. The resulting sequences of at least > 30% similarity and identity were picked for comparative homology modelling. The Supplementary   of psiBLAST. PDB ID 2CRA was found to have 100% identity with the query sequence and PDB ID 2LD5 and 2L7Z showed 78% identity with the query sequence and were chosen as the template for modelling the HOXB13 protein. The respective ".pdb" files of the above-mentioned proteins were downloaded and kept in the same folder where the python script files were located. These three PDB ID structures -2CRA, 2LD5 and 2L7Z were used as the template for modelling the complete HOXB13 protein using comparative homology modelling by MODELLER v9.17. The distance tree of the query sequence and the protein structures from the PDB computed by psiBLAST are available as Supplementary Fig. S4. From the results, we found that 2LD5 and 2L7Z were both structurally and sequentially identical with the same crystallographic resolution of 1.0 Å (as Supplementary Fig. S5). Conversely, the structure 2CRA was found to be diversified from both 2LD5 and 2L7Z with a distance score of 63.5. Hence, the structure 2L7Z was finally selected for modelling of complete HOXB13 because of its high sequence and structural similarity to the query sequence. The alignment of the query sequence (NP_006352) with the template structure 2L7Z was done ( Supplementary Fig. S6). MODELLER v9.17 was instructed to generate 30 similar models of complete HOXB13 protein based on the 2L7Z template structure and "hox13-2l7z.ali" file ( Supplementary Fig. S7). There are several criteria to select the best model among the various models generated by MODELLER v9.17. The most important and widely practised criteria includes selecting the model with the lowest DOPE score 58 and the highest GA341 score 59 . Accordingly, from the summary of the models generated ( Supplementary Fig. S7.), we found that the 26 th model "hoxb13.B99990026.pdb" had the lowest DOPE score of − 10661.23047 and the highest GA341 score of 1.00000. Thus, the selected model "hoxb13.B99990026.pdb" was subjected to further validation of protein structure and folding properties with the help of the Ramachandran Plot and PDBsum.
pdb" was validated and authenticated as the best-generated model and was subjected to analysis for the backbone dihedral angles (phi and psi) of the amino acid residues in the protein structure 44 . For a good protein structure, it is expected that there should be more than 90% of the residues in the core or favoured region of the protein 36,43 . The generated model "hoxb13.B99990026.pdb" was analyzed by RAMPAGE and was found to have 250 residues (88.7%) in the favored region, 18 residues (6.4%) in the allowed region and 14 residues (5.0%) in the outlier region, respectively (Fig. 5.) The model was found to be good and reliable since approximately 89% of the residues fell in the favored region and also because of the low DOPE and high GA341 score.
Model validation by PDBsum. The simulated HOXB13 protein model was further validated with the help of PDBsum for information regarding the motifs, helices, strands, domains, tunnels, angles, positions, error, etc., present in the 3D structure of the proteins 60 . The "hoxb13.B99990026.pdb" was subjected to analysis by PDBsum and was found to have three alpha helices, three helix-helix interactions, 18 beta, and 40 gamma turns, respectively (Fig. 6). The results were in accordance with the features of the homeobox domain of the HOXB13 protein (2CRA). The complete 3D structure of the protein is given in Fig. 6(a). The complete HOXB13 protein, which was modelled using MODELLER v9.17, contained the same features and folding patterns of the homeobox domain of the HOXB13 protein (2CRA), which was clearly evident from Fig. 6(a,b). The detailed protein 3D structure features can be found as Supplementary Fig. S8.
Thus, the protein model "hoxb13.B99990026.pdb" generated by MODELLER v9.17 was found to be the best model based upon the DOPE and GA341 scores and was further validated to be good with the help of the Ramachandra Plot and PDBsum analysis. Hence, this model was taken as the complete human HOXB13 protein structure for further analysis of the corresponding deleterious nsSNP's. The model "hoxb13.B99990026.pdb" was denoted as "HOXB13_M26" protein structure in the subsequent analysis.

Mapping the missense amino acid variation into the protein. Protein template for performing muta-
tion and subsequent analysis. The HOXB13_M26 protein structure was taken as the complete native protein structure for mapping the previously predicted 21 deleterious nsSNP's and also, for further studying their effect on the protein.

Protein mutation and Energy minimization of the native and mutated protein.
The 21 nsSNP's screened to be potentially deleterious by various servers were mapped into the HOXB13_M26 protein using the "mutation" tool in SwissPDBViewer 40,61 . The resulting 21 mutated proteins were denoted as "HOXB13_M26" Mutant.
In order to mimic the in vivo folding conditions and parameters of the protein, energy minimization of both the native (HOXB13_M26) and all the mutant proteins (HOXB13_M26 Mutant) was done with the help of NomadRef Gromacs Server using conjugant gradient force fields 42 . The resulting energy values of all of the native and the mutant structures are given in Table 8.
The total energy of the native protein structure HOXB13_M26 was determined to be − 4505.484 KJ/mol. Among all the 21 mutants, the mutant C66R and S122R were found to have the highest energy of − 4763.567 KJ/mol and − 4745.173 KJ/mol, respectively even after energy minimization, when compared with the native structure. The mutants D167N (− 477.684 KJ/mol), H30Q (− 4675.440 KJ/mol), C63G (− 4663.925 KJ/mol), Q181R (− 4601.128 KJ/mol), C63Y (− 596.328 KJ/mol) and Q188R (− 4588.476 KJ/mol) were found to have considerably higher energy values, whereas the mutants R215C, R123H and Q138H showed very small energy values of − 4273.583 KJ/mol, − 4274.708 KJ/mol, − 4376.856 KJ/mol, respectively, after energy minimization. These nsSNP's with very high and low energies when compared to the native protein structure implied a possible underlying damaging effect on the protein structure, thereby affecting the protein stability and function. The other nsSNP's were found to have near equal energy values as compared with the native structure. The electron cloud density maps of the variants C66R and S122 that had the highest energy are given in Fig. 7(a).
Scientific RepoRts | 7:43830 | DOI: 10.1038/srep43830 RMSD value calculation of the modelled protein. Among the 21 mutants analyzed, the variant H30Q was found to have the highest RMSD value of 2.02 Å, followed by C63Y and P59L, having 1.80 Å each, respectively ( Fig. 7(b)). The mutants D167N, C66R, G177V, R25Q and Y80C, were found to have RMSD values of 1.69 Å, 1.68 Å, 1.56 Å, 1.52 Å and 1.2 Å, respectively ( Table 8). The remaining mutants were found to have RMSD values of less than 1 Å. Among the mutants with high RMSD values, the mutants D167N, C66R, C63Y, P59L and H30Q were found to have both increased energies after energy minimization and RMSD values, which was of critical importance and was taken into further consideration in the subsequent analysis.
Predicting the change in stability of the mutant proteins by I-Mutant Server. Among the 21 nsSNP's submitted, I-Mutant predicted an increase in the stability of 4 mutants, namely rs769634543 (G135E), rs775273363 (A128V), rs8556 (S122R) and rs138213197 (G84E) ( Table 8). The remaining nsSNP's were predicted to be associated with decreased stability. Also, the four variants that were predicted to have increased stability were also found to have low RMSD values. The RMSD values and the I-Mutant results were found to be in conjunction, but the authenticity was yet to be verified by Ramachandran Plot.
Validation of the native and the mutant model using Ramachandran Plot. The energy minimized native (HOXB13_M26) and mutant (HOXB13_M26 Mutant) protein structures in.pdb format were submitted to RAMPAGE for assessment. The native (HOXB13_M26) contained 217 residues (77.2%) in the favored region, 49 residues (17.4%) in the allowed region and 15 residues (5.3%) in the outlier region, respectively (Fig. 8). Interestingly, the mutants G84E, G135E and A128V showed increased positive pattern when compared with the native protein. G84E had 220 (78.3%), 46 (16.4%), and 15 (5.3%) residues in the favored, allowed and outlier regions, respectively (Fig. 8). Three residues from the allowed regions were shifted to the favored region in the G84E mutant and resulted in a better pattern than the native protein. The mutant G135E also showed an Figure 5. Ramachandran Plot for the generated protein model "hoxb13.B99990026.pdb". Almost 89% of the amino acid residues in the modelled protein "hoxb13.B99990026.pdb" occupied the favored region, 6% of the residues occupied the allowed region, and the remaining 5% of the residues occupied the outlier region, respectively.
Scientific RepoRts | 7:43830 | DOI: 10.1038/srep43830 increased and stable amino acid residue pattern with five residues shifting from the allowed region to the favored region with a total of 222 (79%), 43(15.3%) and 16 (5.7%) residues in the favored, allowed, and outlier regions, respectively ( Fig. 8) (Table 8). The mutant A128V also showed a similar increased stabilizing pattern, where two residues from the allowed region became shifted to the favored region. The mutant A128V contained 219 (77.9%), 47 (16.7%) and 15 (5.3%) residues in the three regions, respectively (Fig. 8). The variant S122R, which was predicted to have increased stability by I-Mutant, showed no trace of increased pattern in the Ramachandran plot. It exactly resembled the native protein structure. The mutants Y80C, C66R and P59L, were found to have a destabilizing pattern of amino acid residues. Y80C had (220, 45,16), C66R had (222, 43,16), and P59L had (217, 48,16) residues in the favored, allowed, and outlier regions, respectively. Interestingly, the mutants C66R and P59L were also found to have higher energy and RMSD values, whereas Y80C was reported to have a higher RMSD value. Among all the mutants, the mutant R215C were predicted to have the most destabilizing and damaging combination of amino acid residues with 216 (76.9%), 51 (18.1%), 14 (5.0%) residues in the favored, allowed, and outlier regions respectively. Four residues were shifted from the favored region to the allowed region. The mutant C63G was also found to have a similar destabilizing and damaging pattern to R215C. The remaining mutant models showed near similar or acceptable dihedral angles, which were predicted to confer less damaging effect to the protein when compared to the above-mentioned mutants.

Discussion and Conclusion
The thorough computational analysis of the 21nsSNP's mapped into the HOXB13_M26 protein model, it was predicted that 7 nsSNP's rs761914407 (R215C), rs8556 (S122R), rs138213197 (G84E), rs772962401 (Y80C), rs778843798 (C66R), rs770620686 (P59L) and rs587780165 (R25Q) were found to have seriously damaging and deleterious effects on the HOXB13 with respect to DNA binding and function. Interestingly, the variants G84E, G135E, A128V were found to result in the increased stability of the protein structure. G84E variation was widely reported to be present in large cases of hereditary prostate cancer as epidemiologically reported elsewhere 8,12,17,25,[62][63][64][65] and was in agreement with the results of this study. G135E was also widely reported to be present highly among the Chinese population 66 as reported elsewhere. Thus, the G84E, G135E and A128V variations were predicted to cause some severe structural changes in the protein, which renders it more stable with an increased    half-life. It has also been reported that the gene RFX6 transcribed and regulated by a HOXB13 transcription factor is activated and expressed over a longer period, which results in an imbalance in the feedback mechanism 5,67 that is under the vigilance of HOXB13. It has been scientifically proven that the overexpression of RFX6 helps in the prostate cancer cell migration and disease progression 5,6,12,67 . The variants that cause an increase in the stability of the HOXB13 were found to constitutively express the downstream genes under the influence of HOXB13 and thus it is predicted to be associated with an increase in the risk of prostate cancer, like in the case of RFX6 5 . The mutants R215C, Y80C, C66R and P59L, were found to have highly damaging and deleterious structural and functional properties. This in turn might disturb the role of HOXB13 as a transcriptional factor in activating the genes responsible for cell cycle control and proliferation, eventually leading to the malignancy of the prostate. While G84E, G135E, and A128V were found to increase the stability of the protein structure, the other four nsSNP's, R215C, C66R, Y80C, and S122R, were found to have decreased protein stability. The exact mechanism and the role of those predicted nsSNP's with increased or decreased energy levels and protein stability should further be validated in vitro, since practically either the increased state or the decreased state might possibly involve in the altered patterns of protein structure, function and disease progression. Figure 9 shows the map of the predicted nsSNP's. In addition, three 3′ UTR variations rs543028086, rs550968159, and rs563065128, were found to affect the UNR_BS, GY-BOX and MBE UTR signal, respectively present in the 3′ UTR of the HOXB13 gene, which was predicted to result in the deregulation of the pro-apoptotic factor (Apaf-1). This altered the pattern and regulation of the translational repression of HOXB13 via a feedback mechanism. Thus, it results in the loss of regulating the process of stem cell renewal by blocking or deregulating the cell cycle promotion inhibitors respectively, thereby causing severe damage to the HOXB13 mediated gene expression and function. Out of the 95 nsSNP's subjected to analysis for the variation in the miRNA patterns by the PolymiRTS database five nsSNP's (rs8064432, rs79812861, rs148791210, rs184053751 and rs183620920) were found to disrupt only the conserved (ancestral allele with support > = 2) miRNA sites. Two nsSNP's (rs116931900 and rs1056656) were exclusively found to create a new miRNA site. Conversely, the remaining 16 nsSNP's were predicted to be involved in the disruption and creation of a new miRNA site, out of which rs61123825 (disrupting -2 and creating -7) and rs192244427 (disrupting -4 and creating 5) were found to have a maximum of 9 pattern changes. Thus, the above-mentioned nsSNP's showed up in the progression of prostate tissue malignancy either due to the increase in the stability and half-life of the HOXB13 encoded transcription factor or due to the damaging effects on the protein structure, which resulted in altered binding patterns of the transcription factor, thereby eventually leading to prostate tissue malignancy. The exact mechanism underlying the onset of hereditary prostate cancer by HOXB13 nsSNP's needs to be evaluated and studied extensively with the help of in vivo models and GWAS studies. This study, thus, paves the gateway for future GWAS and clinical studies related to the role of SNPs in hereditary prostate cancer and also has the potential in developing a mechanism for drug targeting and biomarkers for PCa theranostic applications.

Materials and Methods
Datasets. The  Prediction and Screening of deleterious nsSNP's. The highly deleterious missense SNPs associated with the non-homeobox region of HOXB13 gene were predicted using the following in silico servers: The SIFT (Sorting Intolerant From Tolerant) program (http://sift.bii.a-star.edu.sg/www/SIFT_BLink_submit.html) predicts the deleterious or damaging nature of the missense SNPs based upon sequence homology based prediction, physical properties of amino acids and also by calculating the degree of evolutionary conservation of the sequence among various species 68 . The SIFT results were reliable and the scores generated by SIFT program were classified as affecting protein structure (0.00-0.05) and as tolerated (> 0.05). The PolyPhen (Polymorphism and Phenotyping) server (http://genetics.bwh.harvard.edu/pph2), screens and predicts the deleterious nsSNP's based on the observable structural changes induced by the nsSNP's with the help of various proven algorithms 69,70 (THMM, Colis2 program, SignalP program, etc.). These structural changes are in turn known to affect the protein function and stability deleteriously. The relative solvent accessibility and secondary structure details were predicted using DSSP database. The PANTHER (Protein Analysis through Evolutionary Relationships) server (http://pantherdb.org/tools/cSNPscoreForm.jsp?), calculates the duration of a given amino acid that has been evolutionary preserved among various species and predicts the effect of that specific amino acid change on the structural and functional effect on the protein 71 . The longer the amino acid is conserved during the course of evolution, the greater the likelihood of having functional importance in protein structure and function. The PROVEAN (Protein Variation Effect Analyzer) server (http://pantherdb.org/tools/cSNPscoreForm.jsp?) relies upon the data corresponding to the standard properties of the amino acids and protein structure, thereby predicting the effect of the amino acid variations in the protein structure, stability, and function 72 . The nsSNPAnalyzer (http://snpanalyzer.uthsc.edu/) is a tool to predict the phenotypic effect of the missense SNPs based on the data from MSA and three-dimensional protein structure 73 . The PhD-SNP (Predictor of human Deleterious Single Nucleotide Polymorphisms) server (http://SNPs.biofold.org/phd-snp/phd-snp.html ) functions with the help of support vector machines based (SVM-based) and evolutionary information of the sequences 20 . The nsSNP's, which are commonly predicted by more than five servers, were taken into further consideration and analysis. The UTRscan server (http://itbtools.ba.itb.cnr.it/utrscan) is a pattern match identifier that finds the UTR pattern motif match from the protein or nucleotide sequences from the UTRsite databases using UTRblast function 45 . The PolymiRTS (Polymorphism in microRNA and their Target Sites) database (http://compbio.uthsc.edu/ miRSNP/) contains comprehensive data of all the nucleotide variations occurring in the miRNA seed regions and miRNA target sites 74 . The amino acid sequence of the HOXB13 protein (NP_006352) in FASTA format was fed to the server along with its corresponding amino acid variations (ex. G135E).
Modelling the complete HOXB13 protein using MODELLER v9.17. The complete HOXB13 protein was modelled using MODELLER v9.17, which is a computer program used for comparative homology modelling of protein structures. The MODELLER v9.17 was downloaded from the Andrej Sali laboratory website (https:// salilab.org/modeller/). Since MODELLER v9.17 runs on Python scripts, Python was also installed along with MODELLER v9.17 32,34,75 . The Python scripts in MODELLER v9.17 can be executed by the command "mod9.17 script1.py". The basic steps involved in homology modelling using MODELLER are the initial template selection using BLAST, final template selection and alignment of the query sequence with the template structure, building the model based on the final template selected, followed by model evaluation and validation using Ramachandran Plot and PDBsum. The finally validated model waas saved as the HOXB13_M26.
Mapping the missense amino acid variation on the protein structure. Template. The best model was chosen from the various models generated by MODELLER v9.17 and it was used as the template for incorporating the mutations into the protein structure and for subsequent evaluation. This template was taken as the native protein.
Mapping the nsSNP's and Energy minimization of the modelled protein. Each mutant model (21 models) was generated using the "mutation" tool in SwissPDBViewer. Energy minimization of the native and the mutated protein was carried out using NOMAD-Ref Gromacs server (http://lorentz.immstr.pasteur.fr/gromacs/minimiza-tion_submission.php). The NOMAD-Ref Server utilizes Gromacs using conjugant gradient force fields for energy minimization according to the steepest descent, conjugate gradient, or L-BFGS methods. The conjugate gradient method was utilized in this study.

RMSD value calculation of the modelled protein.
The RMSD of the atoms upon superimposing the native and the mutant protein structure was calculated using SwissPDBViewer by the "Calculate RMS" function. The extent of structural deviation between the native and the mutant protein structures was found to have an associated functional effect on the protein, which was predicted by calculating the RMSD by superimposing the native and protein structures. The higher the RMSD value, the higher the structural deviation and associated function of the protein 76,77 .
Predicting the change in stability of the mutant proteins. The stability study of the native and the mutant protein structure was crucial and was carried out with the help of the I-Mutant Server (http://folding.biofold.org/ cgi-bin/i-mutant2.0.cgi). The stability of the protein and its structural changes were predicted by I-Mutant server based on calculating the relative solvent accessibility area, amino acid properties, evolutionary, and structural information of the protein 22 . The server uses the FOLD-X prediction algorithm. The input to the server was the HOXB13 protein sequence (NP_006352), and the amino acid variations were provided manually for each variation.
Validation of the native and the mutant model using Ramachandran Plot and PDBsum. The Ramachandran Plot was used to calculate the dihedral angles of the amino acid residues and to predict the energetically allowed residues based upon their phi and psi dihedral angles, thereby ascertaining the structural and functional properties of the protein structure 36,43,44 . The energy minimized native and the mutant protein models were validated with the online tool RAMPAGE (http://mordred.bioc.cam.ac.uk/~rapper/rampage.php). A percentage of more than 90% residues in the favored region is required for a good protein structure. PDBsum provides the 3D protein structure information regarding the motifs, domains, helices, beta sheets and strands, angles, etc. PDBsum can be accessed online at (http://www.ebi.ac.uk/thornton-srv/databases/cgi-bin/pdbsum/GetPage.pl?pdbcode=index.html).