NCBI was used to retrieve NRAS gene variants using variation viewer and by selecting dbSNP as the source database (https://www.ncbi.nlm.nih.gov/variation/view/). “NRAS” or 4,893 [geneid] was used as an entry. We filtered those retrieved SNPs and selected only the missense variants for further analysis and screening.

We used six different in silico tools for predicting the most deleterious variants on the function of NRAS protein: SIFT (Sorting Intolerant from Tolerant), PolyPhen (Polymorphism Phenotyping), PROVEAN (Protein Variation Effect Analyzer), SNAP, SNP&GO, and PHD-SNP. The SNPs that were predicted to be deleterious by all tools were considered the most deleterious ones. Using six different tools with varying algorithms, an approach that was applied previously (Hossain et al., 2020; Dhar et al., 2022; Rajput and Gahlay, 2021), could increase the confidence and accuracy of our analysis.

The SIFT server relies on a sequence homology approach in addition to the properties of residues to analyze the probability of developing harmful impacts with the missense variants (https://sift.bii.a-star.edu.sg/) (Sim et al., 2012). PolyPhen-2 analyzes the impacts of amino acid substitutions based on physical considerations in addition to comparative methods (http://genetics.bwh.harvard.edu/pph2) (Adzhubei et al., 2010). The PROVEAN tool uses fast methods to compute the scores of pairwise alignment and calculates the required prediction of the impact of the amino acid substitutions (http://provean.jcvi.org/seq_submit.php) (Choi and Chan, 2015). The SNPs&GO tool depends on the functional annotation of the analyzed proteins to identify the deleterious gene variants (https://snps.biofold.org/snpsand-go/snps-and-go.html) (Capriotti et al., 2013). The PHD-SNP server uses support vector machines for analyzing the association between the resulting phenotypes due to amino acid substitution and the development of human diseases (http://snps.biofold.org/phd-snp/phd-snp.html) (Capriotti et al., 2006). SNAP2 has a novel neural network that enables it to identify the effect of SNPs and the neutral SNPs (https://rostlab.org/services/snap/) (Hecht et al., 2015).

SNPs’ impact on NRAS protein stability was analyzed using both the I-Mutant 2.0 server and the Mu-Pro tool. The I-Mutant 2.0 server (https://folding.biofold.org/i-mutant/i-mutant2.0.html) relies on a support vector machine for the prediction of the direction and the value of the free energy change (DDG) (Capriotti et al., 2005). I-Mutant 2.0 was tested on the ProTherm database, which comprised the largest experimental information about the change in free energy related to the stability of proteins with mutations (Bava et al., 2004). The Mu-Pro tool (http://mupro.proteomics.ics.uci.edu/) provides a support vector machines method that was tested and cross-validated with an accuracy of 84% (Cheng et al., 2006).

We used InterPro to determine the locations of the most deleterious variants on the domains of NRAS protein (https://www.ebi.ac.uk/interpro/). InterPro could perform the functional analysis related to the nominated protein and determine its functional domains and important motifs (Blum et al., 2021).

We used the ConSurf bioinformatics tool to analyze the evolutionary conservation of the residues of NRAS protein (https://consurf.tau.ac.il). ConSurf can analyze the phylogenetic relationships that exist between various homologous sequences to perform the needed analysis with the calculation of the conservation score that ranges from one to nine. Furthermore, the structural and functional residues are identified as well (Berezin et al., 2004; Ashkenazy et al., 2016).

We used the SOPMA server (https://npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl?page=/NPSA/npsa_sopma.html) for identifying the predicted secondary structures of NRAS protein and the alignment of the mutant amino acids in these secondary structures. The SOPMA server can predict the secondary structures of a selected protein by analyzing the multiple alignments of this protein’s sequence (Geourjon and Deléage, 1995).

The filtration step of missense SNPs resulted in 17 SNPs. The 3D structures of these filtered SNPs were predicted using SWISS-MODEL, a fully automated protein homology modeling server (Waterhouse et al., 2018). This server performs a fully automated prediction assessment, relying on the continuously deposited protein tertiary structures in the protein data bank. Regarding the NRAS gene, available data on UniProt demonstrated that the site 166–185 represents a hypervariable region, and for this purpose, the sequence 1–165 was uploaded to SWISS-MODEL to select our template for the tertiary structure prediction of 17 mutated proteins with the filtered SNPs. Following structure prediction, the generated models were further validated through Ramachandran plot analysis, which was integrated as an assessment tool in the SWISS-MODEL server, in addition to the ProSA webserver (Wiederstein and Sippl, 2007), which was also utilized for the validation process. Moreover, the TM-Align web server (https://zhanglab.ccmb.med.umich.edu/TM-align/) was employed to estimate the TM-score and root-mean-square deviation (RMSD) values of the parent model selected by SWISS-MODEL and the mutant ones that were predicted. The TM-score gives information about the degree of similarity between submitted structures, and its value ranges between 0 and 1 (a value of 1 represents a complete match between analyzed structures). On the other hand, the RMSD value indicates the average distance between alpha-carbon backbones of analyzed models, and the higher this value is, the more deviation is predicted between analyzed molecules (Zhang and Skolnick, 2005).

To estimate the consequences of filtered SNPs on the cellular behavior of NRAS protein, a molecular docking study between NRAS (the wild-type and the predicted mutated models) and GTP molecule (the ligand of NRAS active state) was proposed. As mentioned in the introduction section, NRAS is a GTP-binding protein that is switched on by incoming signals, leading to a turn-on for genes related to cell growth and differentiation. Therefore, SNPs in NRAS may stimulate a permanently activated form (NRAS-GTP), leading to overactive cellular signaling and ultimately cancer. The NRAS wild tertiary structure (PDB ID: 5UHV) was downloaded in PDB format to act as a control receptor for the current docking study; the deposited molecules in the receptor were first removed and the generated molecule saved to act as a control for its binding score with GTP to validate the docking scores of GTP-binding with other mutated models. The docking was performed through AutoDock Vina (Oleg and Arthur, 2010), the 3D conformation of the docked complexes was visualized by the molecular graphics system PyMOL (Seeliger and De Groot, 2010), and a 2D chart of the interacting residues for the generated complexes was investigated by LIGPLOT (Wallace et al., 1995).

Docked GTPase NRAS models were subjected to molecular dynamics simulations using GROMACS-5.1.4, CHARMM36m, and CHARMM-General forcefields within the TIP3P water solvation model and under periodic boundary conditions (Páll et al., 2015). Protein ionization was set at physiological pH 7.4, and system neutralization was done via chloride and potassium ions. Steepest-descent algorithm-minimization steps were performed at 5 ps, followed by equilibration at the initial NVT ensemble (Berendsen-temp method; 100 ps at 303.15 K) and final NPT ensemble (Parrinello-Rahmann barostat; 100 ps at one atmospheric pressure and 303.15 K) (Elhady et al., 2021). Molecular dynamics were run for 100 ns under NPT ensemble and Particle-Mesh-Ewald algorithms for computing long-range electrostatic interaction. Trajectory analysis was performed using RMSDs (Å), RMS fluctuations (RMSFs; Å), gyration radii (Rg; Å), and solvent-accessible surface area (SASA; nm²). Free binding energies for GTP substrate binding to GTPase NRAS were estimated using Molecular Mechanics/Poisson-Boltzmann (MM-PBSA; kJ/mol) single trajectory calculations (Kumari et al., 2014b). Visual Molecular Dynamics (VMD V.1.9.3) software (Illinois University, Urbana-Champaign, United States) was used for hydrogen bond analysis. Hydrogen bond distance/angle cut-offs were defined at 3.0 Å/20°. Conformational analysis and visualization of the simulated complexes at specified timeframes were performed using the PyMOL software.

The collective flexibilities/motion functions of the constructed GTPase NRAS were assessed using the iMODS online server (http://www.imods.chaconlab.org/) by applying the elastic torsional network model, whose degrees of freedom are the protein backbone torsion angles (López-Blanco et al., 2014). This approach is fast and accurate while being capable of assessing the collective protein’s motion and inherited dynamics based on normal-state analyses of respective internal dihedral angles and torsional coordinates. Furthermore, this server can predict several parameters reflecting structural flexibility/deformation reflecting significant deviation from the normal distribution values obtained from thousands of deposited reference sets.

We used the HOPE server (https://www3.cmbi.umcn.nl/hope/) for analyzing the SNPs’ impact on the NRAS protein 3D structure. The HOPE server depends on several sources for collecting the needed information in addition to building homology models using the YASARA program to perform the required function (Venselaar et al., 2010).

The MusiteDeep server (https://www.musite.net) was used to predict the sites of different types of PTMs. PTMs have a major role in regulating the function of proteins; thus, the identification of PTMs is important in analyzing disease pathogenesis (Wang et al., 2017; Wang et al., 2019; Wang et al., 2020).

We used the GeneMANIA tool (http://www.genemania.org) for generating the gene–gene interaction network of the NRAS gene. The GeneMANIA tool could predict the strongly interacted genes with a selected gene using various sorts of data and resources: co-expression information, physical interaction information, co-localization data, functional relationships, information regarding pathways, genetic interaction data, and data regarding protein domains (Warde-Farley et al., 2010).

The NRAS gene (NCBI Gene ID: 4893) is a protein-coding gene that is located at 1p13.2. It comprises seven exons and has a length of 12,303 nucleotides; it has one transcript (ensemble.org). The membrane protein encoded by this gene shuttles between the plasma membrane and the Golgi apparatus and is characterized by GTPase activity (https://www.ncbi.nlm.nih.gov/gene/4893). The subcellular localization of the NRAS gene is illustrated in Supplementary Figure S2A (Compartments.jensenlab.org/), and its gene ontology is shown in Supplementary Figure S2B (Genecards.org).

In the NRAS gene, 4,542 single nucleotide variations were found (accessed 31 December 2021). Among these SNPs, there were 113 missense SNPs, 62 synonymous SNPs, 73 5′ untranslated region (UTR) variants, 1,405 3′ UTR variants, and 2,435 intron SNPs, in addition to the other downstream and upstream SNPs.

Six different bioinformatics tools (SIFT, PolyPhen-2, PROVEAN, SNAP, SNP&GO, and PHD-SNP) were used to determine the deleterious variants in the NRAS gene with a significant impact on NRAS protein function. Of note, 16 variants were found to be disease-causing and deleterious by all six tools, and the percentage of each SNP type is shown in Supplementary Figure S3. Among these 16 variants, one SNP (rs121434595) was found to have two mutant alleles, and both alleles could result in amino acid substitution (G13C, G13R). Table 1 shows the prediction and scores generated by the six tools for these 16 SNPs.

The impacts of SNPs on NRAS protein stability were analyzed using the I-Mutant 2.0 server and the Mu-Pro tool. The I-Mutant 2.0 server found that 12 SNPs had a decreasing effect on NRAS protein stability, while the Mu-Pro tool showed a decreasing effect with all SNPs except rs1163400692 SNP. The predicted impacts and values are shown in Table 2 for the 16 SNPs.

The InterPro tool analyzed NRAS protein functionally and revealed the presence of an important domain: a small GTP-binding protein domain (InterPro entry: IPR005225). The positions of the 16 deleterious variants were analyzed, and all SNPs were found to be located in this domain, except for rs1658971260, as shown in Table 3.

The evolutionary conservation of NRAS protein amino acids was analyzed as displayed in Supplementary Figure S4. All SNPs were located in highly conserved positions, except rs1658971260 SNP, which was located in an intermediately conserved position, as demonstrated in Table 3. In addition, six SNPs were located on functional and exposed residues (rs1163400692, rs754428086, rs267606920, rs1557982817, rs1465850103, rs1308441238), while one SNP was located on a structural and buried residue (rs757968407).

The predicted secondary structures of NRAS were analyzed via the SOPMA tool as displayed in Supplementary Figure S5. SOPMA revealed that 86 residues were associated with alpha helix (45.5%), 55 with random coil (29.1%), 35 with extended strand (18.52%), and 13 with beta-turn (6.88%). Furthermore, analysis of the SNPs’ alignment in the secondary structures was performed as displayed in Supplementary Figure S6; nine deleterious SNPs existed in the alpha helix secondary structure, and five deleterious SNPs existed in the random coil secondary structure. Meanwhile, no deleterious SNPs were found in the beta-turn or the extended strand. The detailed alignment of the SNPs is shown in Table 3.

To predict the tertiary structure of the models with 17 different mutations, SWISS-MODEL presented the wild-type NRAS 3D model with PDB ID of 5UHV as the best choice for homology modeling with 99.39% identity with all the mutated submitted sequences. Validation of these models was performed by detecting the ProSA Z score and the percentage of residues in the favored region for each molecule (Table 4 and Supplementary Figure S7), where the corresponding scores for all predicted models confirmed the high quality of the predicted models. Regarding the structural similarity between the native and predicted models, the high values of TM (close to one for all predicted models) and the small scores of RMSD (close to zero for all predicted models) demonstrated the high similarity between the native model and the designed ones (Table 4). This was reasonable as SWISS-MODEL in the first place selected a native structure with high similarity with the submitted mutated sequences to build the required models, and this small variation came from the substitution of one amino acid that represents the SNP.

The current docking study was performed to demonstrate the difference in the binding affinity between native and mutated models with the GTP ligand. The binding score between the native NRAS and GTP was −11.3 kcal/mol, and for the mutant predicted structures, the respective binding score is shown in Table 4. Analysis of the scores showed that all models, except models #11 and #15, exhibited a higher affinity to GTP than the native model, where models #1, #2, #3, #6, and #16 were at the top of the list regarding their binding affinity to GTP.

The 3D complexes of the native model, model #1 (the one with the highest binding affinity), and model #11 (the one with the lowest binding affinity) with GTP are visualized by PyMOL, shown in Figure 1, while a presentation of reacting residues of the same models is shown as a 2D chart by LIGPLOT (Supplementary Figure S8). The general topology of the predicted models illustrated the distinctive two lobes (effector lobe: 1–86 residues and allosteric lobe: 87–166 residues) of the GTPase catalytic domain. The active site for GTP-binding is formed via the P-loop, Switch-I, and switch-II with residue ranges of 10–17, 30–40, and 60–76 amino acids, respectively (Johnson et al., 2017). The GTPase-conserved motifs at the effector lobe, P-loop GxxxxGK (residues 10–16), and DxxG/Thr35 (residues 57–60) were shown to support the GTP-α-β-phosphate groups and Mg2+ coordination (Araki et al., 2011). On the other half, the allosteric lobe showed the other GTPase-conserved motifs including NKxD (residues 116–119) responsible for guanine nucleotide specificity, whereas, ExSAK (residues 143–147) together with Phe28 being responsible for stabilized nucleotide binding, respectively (Parker and Mattos, 2015).

Simulated proteins of the GTPase NRAS models showed steady equilibrated RMSD trajectories, respective to initial proteins’ alpha-carbons, across the 100-ns runs (Figure 2A). The steadiest alpha-carbon RMSDs were depicted for simulated model numbers 2, 6, and 16 as compared to other models, with average tone values of 1.98 ± 0.39 Å, 1.95 ± 0.23 Å, and 2.06 ± 0.17 Å, respectively. Notably, all simulated NRAS models were of much-equilibrated RMSD tones as compared to the native NRAS both at its respective holo (GTP bounded: 2.13 ± 0.45 Å) and apo (unliganded: 2.41 ± 0.65 Å) states. The end molecular dynamics timeframe (60–100 ns) showed high RMSD fluctuations for the apo/unliganded NRAS model, reaching up to ∼5.70 Å. On the other hand, the holo native NRAS showed a slight increase in the protein RMSD following the 75-ns timeframe and till the end of the simulation run.

Monitoring the RMSD trajectories of the GTP ligand, in relation to its initial/reference position, showed steady tones for the simulated models (Figure 2B). The RMSDs of the ligand were much lower than those of their respective bounded proteins, the thing that conferred stable GTP-pocket accommodation across the whole 100-ns simulation runs. The simulated models were further investigated by monitoring the protein’s Rg and SASA trajectories (Figures 2C, D). Low average Rg—SASA values, 15.30 ± 0.08 Å//91.64 ± 1.96 nm², 15.24 ± 0.07 Å//89.64 ± 2.25 nm², and 15.19 ± 0.05 Å//86.69 ± 1.74 nm², were assigned for the simulated model numbers 2, 6, and 16, respectively. Again, the native holo NRAS protein showed a slight increase in its Rg and SASA tones, which then proceeded constantly till the end of the simulation runs. Similar to the protein’s RMSDs, the monitored SASA trajectories were of the highest values for the simulated apo native NRAS protein as compared to other simulated congerants.

Dissecting the stability profile of the simulated GTPase NRAS proteins into their respective residue-oriented level was performed by monitoring the alpha-carbon RMSF trajectories. Difference RMSF values for a particular GTP-NRAS model in relation to the native apoprotein (ΔRMSF = RMSF_Apo—RMSF_Holo) was of better representation for the residue-wise protein stability and fluctuation. Residues depicting ΔRMSF values ≥ 0.3 cut-offs were considered significantly rigid, showing negligible mobility (Albuquerque et al., 2020). Notably, the simulated proteins showed typical thermodynamic behavior, where terminal residues were of high negative ΔRMSF values as compared to most of the core amino acids (Figure 3). On the other hand, certain effector lobe residue ranges depicted high immobility profiles being confined with the GTPase-binding site. These top-immobile residue ranges were the P-loop GxxxGK motif (residues 10–16), switch-I, and vicinal amino acids (residues 26–39), as well as the switch-II range (residues 60–70) with the reported Mg⁺²-coordinating DXXG motif (Table 5).

Moving toward the allosteric lobe residues, the two motifs responsible for guanine nucleotide specificity and stabilized nucleotide binding (NKxD and ExSAK GTPase motifs) depicted a significant immobility profile with ΔRMSF above the 0.30 Å cut-off value. Notably, stability residue-wise ranges were more prominent for the effector lobe rather than those at the allosteric half of the GTPase catalytic domain. Concerning the comparative NRAS flexibility profiles, model #3, along with the native holo protein, showed the highest mobility tones with more negative and less positive ΔRMSF values. It is interesting that SNP NRAS models #1, #2, and #3 depicted consistent mobility and flexibility profiles for the 37–50 residue range, showing almost all its respective residues with negative ΔRMSF values. Moreover, selected switch-II residues (Tyr64, Ser65, Ala66) exhibited high ΔRMSF values for the same three Gly13 SNP models. It is worth noting that the native holo protein depicted negative ΔRMSF values for the same previous switch-I and -II residue ranges yet at less prominent negative numbers. On different bases, the three Gly13 SNP residues, G13R, G13C, and G13V, in models #1, #2, and #3, respectively, showed a slightly higher immobility profile (ΔRMSF ∼0.60 Å) than their native Gly13 in holo NRAS protein. The P34R SNP residue in model #6 showed a more preferential stability profile than its native Pro34 congerant. Likewise, the V152F SNP in model #16 was more rigid, reaching a positive ΔRMSF value (0.24 Å), as compared to its negative-valued native residue (Val152).

Conformational analysis for the simulated NRAS was performed to grasp the main secondary structure changes and SNP residue conformational shift throughout the molecular dynamics simulations. Trajectories at the start and end of the simulation runs were extracted, minimized at 0.0001, and finally overlayed. As expected, the docked GTP molecule was attained within the NRAS-binding site at the end of molecular dynamics simulation runs (Figure 4). Most SNP residues showed limited orientation/conformational changes at the end of the simulation runs. The latter was obvious for G13C, G13V, and V152F in models #2, #3, and #16, respectively. Nevertheless, the G13R SNP residue in model #1 exhibited the highest conformational drift (alpha-carbon RMSD = 1.034 Å) as compared to its G13 SNP congerants (alpha-carbon RMSD ∼ 0.800 Å). The G13R residue showed closer proximity to the GTP phosphate groups as compared to its start frame, favoring phosphate group stability and extended hydrogen bonding. Concerning the P34R SNP residue in model #6, a moderate conformational shift was depicted for this residue with a depicted alpha-carbon RMSD of 0.941 Å.

It is worth mentioning that significant secondary structure conformational alterations were also depicted and such observations were differential across each simulated model. Regarding the G13x SNP model, a significant shift of the switch-I loop of the simulated NRAS proteins at the end of the molecular simulation runs. A switch from the closed state-1 to an open state-2 was depicted with models #1–3, being more obvious for model #1. This switch-I drift was accompanied by an outward movement for the α-helix 3 and loop 7 away from the switch-II domain. The latter drift permitted the switch-II to attain the catalytically active R-state conformation rather than the catalytically incompetent T-state at the start of simulation runs. The above-described switch-I and switch-II movements were less prominent with SNP models #6 and #16, where minimal or insignificant conformational changes were depicted. The native holo NRAS complex exhibited moderate switch-I drift toward the open state-1, yet without significant outward movement of its α-helix 3/loop 7 nor its switch-II domain.

Estimating the free binding affinity of GTP toward the constructed NRAS models as well as the native holo state was performed using the trajectory-oriented MM-PBSA calculations. Using the single trajectory approach and SASA-only model (ΔG Total = ΔG Apolar + ΔG Polar + ΔG Molecular Mechanics), the total binding energies for GTP in all SNP NRAS models were superior to that of the native holo protein (Table 6). This SNP-favored binding affinity is consistent with the results of the preliminary molecular docking investigation. The obtained free binding energies were dissected into their contributing energy terms, showing a general trend in several GTPase NRSA complexes. Dominant energy contributions of the van der Waals interactions (ΔG Van der Waals) over that of the Coulomb’s electrostatic potentials (ΔG Electrostatics) were depicted for SNP models #2, #3, and #16 as well as the native holo complex. The latter van der Waal energy preferentiality was most obvious for SNP model #6. On the other hand, both SNP models #1 and #6 depicted almost two-fold higher ΔG Electrostatics contribution over the hydrophobic potentials, the thing that was also remarkably associated with high polar solvation energies (ΔG Solvation; Polar). Moreover, the non-polar solvation energy terms were almost consistent for all simulated NRAS complexes, being around 19.0 kJ/mol.

The obtained total binding-free energies were further decomposed down to the residue-wise binding energy contributions (Figure 5) to identify the key binding residues. The highest negative energy contribution values were assigned to residues of the above-described motifs being confined with GTP-binding and recognition. The P-loop Lys16 of the GxxxxGK motif and Asp119 of the NKxD motif depicted the highest residue-wise energy contributions across all simulated NRAS complexes, with values of −131.2 up to −157.4 kJ/mol and −130.9 up to −187.5 kJ/mol, respectively. Other GTP active site residues showed significant energy term contributions, including Val14, Ala18, Phe28, Asp30, Tyr32, Asp57, Tyr64, Lys117, and/or Ala146.

Most of the latter top-energy-contributing polar residues showed high hydrogen bond % frequencies with GTP throughout the entire simulation runs using the VMD hydrogen–bond tools. SNP model #1: Arg13-side (30.64%), Lys16-side (93.41%), Ala18-main (44.31%), Asp30-main (30.34%), Asp119-side (123.05%); SNP model #2: Cys13-side/main (22.45%), Lys16-side (69.06%), Ala18-main (14.22%), Asp30-side (50.82%), Asp119-side (118.66%); SNP model #3: Lys16-side/main (107.29%), Ala18-main (52.10%), Asp119-side (123.45%); SNP model #6: Gly13-main (30.84%), Lys16-side/main (131.62%), Ala18-main (34.73%), Asp30-main (32.63%), Tyr32-side (53.29%), Arg34-side (38.82%), Asp119-side (123.25%); SNP model #16: Lys16-side/main (132.23%), Tyr32-side (63.97%), Asp119-side (115.87%); and native holo NRAS: Lys16-side (105.39%), Ala18-main (41.72%), Asp30-main (35.53%), Asp119-side (120.86%).

On the contrary, other hydrophilic pocket-lining residues, including Glu37, Arg38, Glu62, Glu63, and/or Asp69, exhibited unfavored positive-valued energies. The latter binding repellent residues were of the highest values for SNP models #1 and #6. Concerning the SNP residues, G13R showed a higher negative energy contribution in model #1 as compared to its SNP or native congerants. Likewise, model #6 showed significant binding energy for its SNP P34R, with a value reaching up to −113.0 kJ/mol. The SNP V152F residue showed relevant energy contribution in model #16, whereas the native congerant illustrated either poor energy contributions (∼−0.45 kJ/mol) or even positive-value energies (∼0.14 kJ/mol).

Collective functional flexibility/mobility of modeled SNP and native holo NRAS proteins was assessed via the iMODS server and represented in Supplementary Figure S9. In particular, the residue range 33–55, corresponding to the switch-I domain and vicinal residues, depicted the highest deformability across all NRAS proteins (Supplementary Figure S9A). This inherited flexibility pattern was also recapitulated through the estimated highest B-factor scores for these residue regions (Supplementary Figure S9B). The B-factor analytical parameter can quantify the atomic displacement amplitudes around a conformational equilibrium (López-Blanco et al., 2014). The estimated eigenvalues (1) of SNP models #1, #2, and #3 complexes were the highest (∼8.8e⁻⁰⁴) as compared to those of other SNP models as well as the native NRAS protein (Supplementary Figure S9C). In addition, the depicted inverse eigenvalue-variance relationship related to each normal mode predicted a significantly higher intrinsic mobility Gly13 SNP models across the collective function motions (Supplementary Figure S9D).

Regarding the obtained covariance matrices (Supplementary Figure S9E), correlated residue-pair motions (reds) were more assigned to the NRAS proteins as compared to either the uncorrelated (white) or anti-correlated (blue) motions. The depicted correlated residue-pair motions were highly correlated to the switch-I residue ranges as compared to other protein domains. Finally, the elastic network model further described the differential flexibility pattern in Supplementary Figure S9F, where this model illustrated atom-wise pairs linked through springs related to their relative stiffness magnitudes (darker gray strings confer stiffness). Scattered non-continuous dark gray bands, around the normal distribution stiffer string, were assigned for NRAS proteins, particularly the switch-I residues.

The five mutations with the highest binding affinity to GTP were further analyzed for their impacts on NRAS protein structure using the HOPE server. Table 7 explains in detail the different structural impacts in NRAS protein associated with these variants. Structure substitution of wild-type amino acids by mutants is illustrated in Figures 6A–D for G13C, G13R, P34R, and V152F, respectively.

The MusiteDeep server was used to identify the predicted sites of PTM; Figure 7 displays the predicted PTM sites in the NRAS protein. In addition, the analysis predicted the mutant residue C164 to be a palmitoylation site. Moreover, the presence of I55R SNP was found to lead to the appearance of a phosphorylation site at the wild residue of Y64. Furthermore, the presence of D119G SNP was found to lead to the loss of a PTM site at the wild residue of S106.

The GeneMANIA tool was used for generating the NRAS gene–gene interaction network and identifying the genes that strongly interact with the NRAS gene. Figure 8 displays the 20 most strongly connected genes to the NRAS gene. Regarding these genes, the SHOC2 gene (leucine-rich repeat scaffold protein gene) showed the strongest relatedness with a rank of 1, followed by the STK19 (serine/threonine kinase 19) gene and the RGS12 (regulator of G protein signaling 12) gene consecutively.