In silico prediction of deleterious non-synonymous SNPs in STAT3

The nsSNP distribution of the STAT3 gene (Accession: NP_003141.2) was collected from the National Center for Biotechnology Information (NCBI) dbSNP database (https://www.ncbi.nlm.nih.gov/projects/SNP) and the protein sequence was retrieved from the UniProt database (UniProtKB – P40763 (STAT3_HUMAN)).

To estimate the functional repercussions of nsSNPs in the coding area acquired from the dbSNP database, 7 best-performing web tools were employed sequentially.

Sorting Intolerant From Tolerant (SIFT; https://sift.bii.a-star.edu.sg/) [10] predicts the amino acid substitution effects and the damaging effect on protein function based on homology to identify the tolerated and deleterious SNPs. The SIFT probability score ≤0.05 indicates deleterious and those >0.05 indicate tolerated. The rsIDs of nsSNPs acquired from NCBI's dbSNP database were submitted as input queries.

Polymorphism Phenotyping v2 (PolyPhen2; http://genetics.bwh.harvard.edu/pph2/) [11] uses the protein sequence and amino acid substitutions in the sequence to predict the structural and functional effects on the protein. The output of the prediction shows the nsSNPs as “PROBABLY DAMAGING” with a score of 0.7–1, “POSSIBLY DAMAGING” with a score of 0.5–0.8, or “BENIGN”.

Based on the alignment score of a protein, Protein Variation Effect Analyzer (PROVEAN; http://provean.jcvi.org/index.php) [12] analyzes the functional effect of amino acid substitutions in the protein. The SNPs with a PROVEAN score of ≤ −2.5 are considered to have a damaging effect and > −2.5 indicates a neutral effect on the protein.

Based on evolutionary relationships, molecular activities, and interactions with other proteins, Protein Analysis through Evolutionary Relationship (PANTHER; https://www.pantherdb.org/tools) [13] predicts that a specific nsSNP will have a functional impact on the protein using position-specific evolutionary conservation scores. The required input query for this prediction was protein sequence and amino acid variants.

SNPs&GO (http://snps.biofold.org/snps-and-go/snps-and-go.html) [14] used the corresponding protein query sequence to predict the use of a reliability index (RI) that would enable the ascertainment of disease-relevant mutations in a protein sequence.

Predictor of human Deleterious Single Nucleotide Polymorphisms (PhD-SNP; http://snps.biofold.org/phd-snp/phdsnp.html) [15] is an online predictor based on the support vector machine (SVM), and is used to classify the nsSNPs as disease-related or neutral polymorphism.

I-Mutant 3.0 (http://gpcr2.biocomp.unibo.it/cgi/predictors/I-Mutant3.0/I-Mutant3.0.cgi) [16] was used to evaluate changes in protein stability due to nsSNPs based on both protein sequence and protein structure. The stability changes and related ΔΔG value were calculated using I-Mutant 3.0. ΔΔG value is the difference in the Gibbs free energy of the mutated protein and the wild-type [17]. A ΔΔG value <0 indicates that the variation has a negative impact on protein stability. A ΔΔG value >0, on the other hand, suggests that the variation improves protein stability.

The ConSurf web server (http://consurf.tau.ac.il) [18] analyzes the evolutionary rate of conservation of the amino acids in the protein sequence. The input was the FASTA format of protein sequence, which calculates conservation scores in the range of 1–9. A score range of 1–4 indicates variable, a 5–6 score indicates intermediate, and a 7–9 score is a conserved one.

NetsurfP-2.0 [19] predicts the accessibility of solvents, secondary structure, disorders, and phi/psi dihedral angles of amino acids in an amino acid sequence. The FASTA format of the STAT3 protein sequence was submitted as an input query for this server.

HOPE (http://www.cmbi.ru.nl/hope/) [20] was used to investigate the impact of point mutations on the protein structure. A protein sequence or an accession code for the protein of interest was used as the input query. The structural data were gathered from the distributed annotation system (DAS) servers as well as the Uniprot database. HOPE also provides information on structural differences between mutant and normal residues.

SOPMA (https://npsa-prabi.ibcp.fr/cgibin/npsa_automat.pl?page=npsa_sopma.html) [21] was used to predict the secondary structure of STAT3 protein using 5 algorithms. The FASTA sequence of the protein was given as input for SOPMA prediction.

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING; https://string-db.org/) database predicts the functional linkages between proteins as well as their associations by combining data from physical interactions and curated biological pathway databases.

Raptor X gives the prediction of the secondary and tertiary structure of the protein, contact and distance maps, solvent accessibility, disordered regions, functional annotation, and binding sites based on a 3D model. Binding site prediction is calculated by the pocket multiplicity, which determines the quality of the projected pocket. The predicted pocket is more accurate when the score is greater than 40.

GeneMANIA [22] (https://genemania.org/) predicts the gene–gene interaction using a large set of functional association data including protein and genetic relationships, pathways, co-expression, co-localization, and protein domain similarity. STAT3 was the input query entered into the software.

AlphaFold (https://alphafold.ebi.ac.uk/) is the most accurate technique for predicting protein structure by neural network–based model. The AlphaFold model includes the coordinates of every heavy atom in a protein as well as “its confidence in the form of a projected LDDT-C score (pLDDT) per residue” [23]. The LDDT score, which ranges from 0 to 100, describes how closely the protein model resembles the reference structure without superposition [24]. The Uniprot accession code of STAT3 was used to acquire the wild-type protein structures from the AlphaFold Protein Structure Database.

The impact of a single point mutation on the stability, conformation, flexibility of proteins, and the visualization of protein dynamics, were evaluated using the DynaMut server (http://biosig.unimelb.edu.au/dynamut/). DynaMut employs normal mode analysis (NMA) to compare the free energy change (ΔΔG) between the wild-type and mutant structures. DynaMut additionally provides structure-based predictions for mCSM [25], SDM [26], and DUET [27], as well as the ΔΔG prediction of an elastic network contact model (ENCoM) based on NMA. Moreover, DynaMut uses ENCoM-based difference in vibrational entropy (ΔΔSVib) to determine whether the mutation will be more or less flexible. The wild-type structure in PDB format and the variant amino acid were given as input.

Various in silico prediction tools such as SIFT, Polyphen-2, PROVEAN, SNP & GO, PANTHER, and PhD–SNP were used to analyze disease-associated SNPs. Initially, all the 417 nsSNPs were loaded to the SIFT server, which predicted 160 nsSNPs as tolerated or deleterious, whereas the remaining SNPs were not found. Out of this, 67 nsSNPs were classified as deleterious with a SIFT score ≤0.05 and the remaining were tolerated. The nsSNPs were then submitted for PolyPhen2 analysis. To increase the accuracy of the prediction, the combined prediction of both SIFT and PolyPhen, such as the nsSNPs with SIFT score ≤0.05 and PolyPhen score >0.90, was selected and 9 nsSNPs were identified as deleterious. The selected 9 nsSNPs were subjected to other in silico tools, namely PROVEAN, SNP & GO, PANTHER, and PhD–SNP, and the results are given in Table 1.

According to PROVEAN results, 6 nsSNPs were predicted as disease-causing and 3 were neutral. Through the PANTHER tool, 6 nsSNPs were predicted as probably damaging, 1 was possibly damaging, and 1 was found to be benign. Moreover, SNP & GO predicted 6 of these 9 nsSNPs as disease-relevant mutations in STAT3 with a RI. Finally, 6 nsSNPs were determined as deleterious by all the above online software. rs145786768, rs193922716, rs193922717, rs193922719, rs1064116, and rs1803125 with a change in amino acid V507F, R335W, E415K, K591M, F561Y, and Q32K were picked as the most deleterious nsSNPs since they showed a mutation in a majority of the in silico prediction tools. In addition, we used the I-Mutant server to analyze the effects of nsSNPs on protein stability. The result showed that 6 nsSNPs have a ΔΔG value < −0.5, which indicates decreasing stability of the STAT3 protein (Table 1).

Through the ConSurf web server, the evolutionary conservancy of 6 nsSNPs of STAT3 protein was analyzed, together with the identification of putative, structural, and functional residues. The results showed that all the 6 nsSNPS, namely V507F, R335W, E415K, K591M, F561Y, and Q32K, were highly conserved and the variants are indicated in black boxes in Figure 1. The result of ConSurf is shown in Table 2.

Figure 1.

The variants with high conservation scores in ConSurf output were assessed for solvent accessibility, stability, and secondary structure prediction by NetsurfP-2.0. The results are displayed in Table 3.

HOPE showed the difference between the wild-type and mutant amino acids in connection with their physical and chemical properties, hydrophobicity, spatial structure, and function. Project HOPE server revealed that the mutant residues V507F, R335W, E415K, F561Y, and Q32K were bigger than the wild-type residues whereas the mutant residue of K591M was smaller than the wild-type. In addition, the mutant residue of R335W, K591M, and F561Y is more hydrophobic than the wild-type. The change in the size and hydrophobicity of the mutant residue can disrupt the H-bonding interaction with neighboring molecules. Moreover, the mutation at 335, 591 positions the wild-type residue charge was lost, whereas at 415^th position the mutant introduced an opposite charge; further, the mutant also introduced a charge at position 32. The variation in charge can cause a loss of interaction with other molecules. Figure 2 illustrates the structural images of 6 deleterious nsSNPs showing both wild-type and mutant residues at the specific protein site.

Figure 2.

SOPMA predicted the secondary structure of STAT3, which explained the distributions of alpha-helix, beta-sheet, and random coil. The secondary structure prediction of STAT3 by SOPMA is shown in Figure 3. SOPMA secondary structure elements calculations showed that 50.39% of sites were in alpha helixes, 33.38% in random coils, 2.86% in beta twists, and 13.38% in extended strands. Out of the 6 high-risk nsSNPs, 2 were found in random coils, 3 in alpha helixes, and 1 in extended strands. The STRING maps were used to depict the protein–protein interaction of STAT3. The STRING results showed the functional interaction pattern of the STAT3 protein with other proteins in a cell, which is useful in interpreting the genotype–phenotype consequences of mutations. Histone acetyltransferase p300 (EP300), E3 SUMO-protein ligase (PIAS3), interleukin-10 receptor (IL10RA), JAK1, JAK2, epidermal growth factor receptor (EGFR), heat shock protein HSP 90-alpha (HSP90AA1), proto-oncogene tyrosine-protein kinase (SRC), and homeobox protein NANOG were all found to have substantial functional connections with STAT3 in the STRING prediction (Figure 4).

Figure 3.

Figure 4.

STAT3 ligand binding sites were predicted using the RaptorX Binding server. A pocket multiplicity number greater than 40, according to the RaptorX Binding server, suggests an accurate prediction. In the STAT3 protein, there were only 1 processed domain and 2 predicted pockets, one with a multiplicity of 11 that binds to R382 L430 I431 S465 N466 Q469 and the other with a multiplicity of 3 that binds to M331 H332 I467. The server predicted that 128 locations were disordered. RaptorX also predicted the secondary structure showing 45% helix, 13% sheets, and 41% coil. Figure 5 shows the structure predicted by the RaptorX Binding server.

Figure 5.

GeneMANIA server predicted the functional gene–gene interaction network of STAT3, which is shown in Figure 6. According to the predicted results, STAT3 interacts physically and genetically with nuclear factor kappa B Subunit 1 (NFKB1Z) and mitogen-activated protein kinase kinase 5 (MAP2K5), EGFR, and STAT1, and is co-expressed mostly with EGFR, protein disulfide-isomerase A3 (PD1A3), interleukin 6 receptor (IL6R), and STAT5B.

Figure 6.

The overall confidence for the whole protein chain is shown by the average pLDDT scores across all residues The Alpha-Fold algorithm produces a pLDDT score for each individual residue that ranges from 0 to 100. In STAT3 3D structure, a very high degree of confidence (pLDDT >90) was obtained and majority of the 3D structural area belongs to α-helical domains (Figure 7).

Figure 7.

DynaMut calculates the effect of point mutations on protein stability and flexibility based on interatomic interactions, as summarized in Table 4. DynaMut identified the ∆∆G prediction of F561Y as destabilizing; however, NMA-based predictions revealed that all 6 deleterious nsSNPs reduce structural stability when compared to the wild-type protein. Figure 8 depicts the difference in vibrational entropy and interatomic interactions between the wild-type and mutant.

Figure 8.