A comprehensive in silico analysis of the functional and structural consequences of the deleterious missense nonsynonymous SNPs in human GABRA6 gene

doi:10.21203/rs.3.rs-3776390/v1

Download PDF

Research Article

A comprehensive in silico analysis of the functional and structural consequences of the deleterious missense nonsynonymous SNPs in human GABRA6 gene

https://doi.org/10.21203/rs.3.rs-3776390/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Epilepsy, a prevalent neurological disorder, affects more than 50 million individuals worldwide and is characterized by recurring seizures. Nonsynonymous single nucleotide polymorphisms (nsSNPs) found within coding regions of epilepsy-related genes are believed to have significant impacts on protein function. This is due to their tendency to cause mutations in the encoded amino acids, which can subsequently lead to pathogenic alterations in protein structure and function. Consequently, nsSNPs have the potential to serve as diagnostic markers for epilepsy and other neuropsychiatric conditions. The primary objective of this study is to evaluate the deleterious effects of missense nsSNP mutations on the GABRA6 gene. The GABRA6 gene encodes the alpha-6 subunit of the \({\text{G}\text{A}\text{B}\text{A}}_{\text{A}}\) receptor and is associated with genes implicated in epilepsy. To achieve this, we employed various computational tools, including SIFT, PolyPhen-2, PROVEAN, Condel, SNPs&GO, PMut, SNAP2, MutPred2, and SNPeffect4.0, for predicting missense nsSNPs. Additionally, we used I-Mutant3.0 and MUpro to analyze protein stability, ConSurf to assess evolutionary conservation, FTSite and COACH to predict ligand binding sites, SOPMA and PSIPRED to analyze protein secondary structures, project HOPE to predict structural changes, and I-TASSER to model the 3D structure. Furthermore, structural validation was conducted using the PROCHECK and ERRAT servers, while molecular dynamics simulations were performed using GROMACS to gain a better understanding of the effects of mutations on protein structure. Among the 451 missense nsSNPs identified within the GABRA6 gene, three were found to have pathogenic effects on the structure and function of the protein, potentially contributing to the development of epilepsy or other neuropsychiatric disorders. Notably, two of these nsSNPs (W87S and W112R) were located within the ligand-binding domain, while the third (C310R) was situated in the transmembrane domain. These identified missense nsSNPs hold promise as potential diagnostic markers and drug targets.

Epilepsy

GABRA6 gene

Single nucleotide polymorphism

In silico analysis

Genetic marker

Epilepsy is a neurological disorder characterized by recurring seizures (Fisher et al. 2014). While approximately 20–30% of epilepsy cases are associated with tumors, strokes, or head injuries, the remaining 70–80% are believed to have genetic factors (Myers and Mefford 2015). Genetic epilepsy refers to cases where seizures are a primary symptom of the disorder and are directly caused by presumed or known mutations. It is important to note that the underlying mutations responsible for genetic epilepsies are still largely unknown (Scheffer et al. 2017).

In a study conducted by Wang et al. in 2017, 977 genes associated with epilepsy were identified through a comprehensive search of databases such as OMIM, HGMD, EpilepsyGene, and PubMed publications. These genes were classified into four groups based on the phenotypes of epilepsy manifestation. Among the 977 genes, 84 were designated as epilepsy genes, 73 as neurodevelopment-related genes, 536 as epilepsy-related genes, and 284 as putatively epilepsy-associated genes, which require further investigation (Wang et al. 2017). The focus of this study is to evaluate the GABRA6 gene, which falls under the putatively epilepsy-associated genes category.

\({\text{G}\text{A}\text{B}\text{A}}_{\text{A}}\) receptors, which belong to the Cys-loop ion channel superfamily, play a crucial role as the major inhibitory mediators in the mammalian central nervous system (CNS) (Brejc et al. 2001). These receptors are composed of pentameric subunits, with 19 different subtypes identified (ρ1–ρ3, γ1–γ3, β1–β3, α1–α6, θ, π, ε, and δ) (Macdonald and Olsen 1994). The diverse molecular composition and expression of \({\text{G}\text{A}\text{B}\text{A}}_{\text{A}}\) receptors' subunits contribute to the variations in their properties, including agonist binding affinity, conductance, kinetics, and distribution within the brain (Brooks-Kayal and Russek 2012). These differences may potentially contribute to the pathogenesis of epilepsy and chronic seizure recurrence (Fritschy 2008).

The GABRA6 gene encodes the alpha-6 subunit of the \({\text{G}\text{A}\text{B}\text{A}}_{\text{A}}\)receptor. Previous research has demonstrated that variations in the GABRA6 gene are associated with genetic epilepsy disorders such as childhood absence epilepsy (CAE), autosomal dominant juvenile myoclonic epilepsy (JME), and infantile spasms (ISs) caused by substitution mutations (Hirose 2014). For instance, Dibbens et al. reported a mutation (Arg46Trp) in the GABRA6 gene in a patient with CAE. This mutation is located in a homologous region to the γ2 subunit mutation (Arg82Gln) in the gamma-2 subunit of the \({\text{G}\text{A}\text{B}\text{A}}_{\text{A}}\)receptor, which is known to be associated with febrile seizures and CAE in humans. The γ2 subunit mutation results in reduced surface receptor levels and receptor currents, indicating impaired assembly and functionality of \({\text{G}\text{A}\text{B}\text{A}}_{\text{A}}\) receptors (Dibbens et al. 2009). Furthermore, Hernandez et al. found that the CAE-associated Arg46Trp mutation in the GABRA6 gene can cause neuronal disinhibition in the \({\text{G}\text{A}\text{B}\text{A}}_{\text{A}}\) receptor, leading to an increased predisposition to generalized seizures due to reduced function and expression of αβγ and αβδ receptors. This mechanism shares similarities with the γ2 subunit mutation (Arg82Gln) linked to CAE in humans (Hernandez et al. 2011).

When comparing genomic DNA sequences from different individuals, it is observed that certain positions can contain two or more bases. These variations, known as single nucleotide polymorphisms (SNPs), are highly prevalent, occurring at a frequency of approximately one out of every 0.3–1 kilobase in the human genome. Additionally, SNPs can have distinct effects on phenotypic levels depending on their location (Syvänen 2001). SNPs with a minimum frequency of 1% in the general population are believed to account for approximately 90% of inter-individual variability and nearly 100,000 amino acid differences. Due to their low mutation rate, SNPs serve as valuable markers for investigating complex genetic characteristics. Consequently, SNP markers are increasingly utilized in genetic studies due to their numerous advantages. In general, SNPs can be categorized as either intronic or exonic. Intronic SNPs are located in non-coding regions of the genome and do not impact protein products when translated. On the other hand, nsSNPs are exonic SNPs found in coding regions that result in changes to protein length or amino acid variants (Mah and Chia 2007). Among the various types of SNPs, nsSNPs in coding regions are considered to have the most significant effect on protein functions, often leading to mutations in encoded amino acids that can detrimentally affect protein structure and function (George Priya Doss et al. 2008).

The emergence of large, complex biological data sets has necessitated rapid advancements in the field of bioinformatics or computational biology. While not a new field, bioinformatics has become increasingly indispensable as biology transitions towards high-throughput methods for analyzing entire genomes. Bioinformatics encompasses the scientific discipline involved in the acquisition, storage, distribution, processing, interpretation, and analysis of biological information. Consequently, numerous databases and software applications have expanded to facilitate the analysis of genetic and physical data in order to infer the functions of model organisms and individuals. Computational analysis provides effective and alternative methods for generating new hypotheses, designing appropriate investigations, and interpreting vast amounts of data derived from genome-scale research. Computational or in silico analyses supplement traditional experimental biology techniques by significantly enhancing predictive capabilities (Fielden et al. 2002). The term "in silico" refers to experiments conducted using computers and is analogous to the more familiar biological terms "in vitro" and "in vivo" (Ekins, Mestres, and Testa 2007). Given the exponential number of SNPs, it is impractical to determine the biological significance of each SNP through wet laboratory techniques alone. However, computational tools can be employed to screen potentially deleterious SNPs that may alter protein function and structure before proceeding to wet laboratory analysis (Mah, Low, and Lee 2011).

The objective of this study was to identify and validate the most deleterious missense nsSNPs associated with the GABRA6 gene and assess their impact on protein function and structure. This approach involved investigating the underlying molecular mechanisms using various in silico methods. Predicted missense nsSNPs, along with the native form, were further analyzed through molecular dynamics simulations to gain a better understanding of their effects on protein structure. This investigation aimed to shed light on how these missense nsSNPs may contribute to epilepsy and other neuropsychiatric disorders by influencing the expression of the protein encoded by the GABRA6 gene. The findings of this study could aid in the discovery of biomarkers for SNPs in epilepsy diagnosis resulting from these variations, as well as contribute to the development of targeted therapeutic interventions in the future.

The deleterious nature of missense nsSNPs in the function and structure of the human GABRA6 gene was predicted by several computational algorithms. To ensure that the results have high accuracy, various in silico tools were utilized for prediction. The approach followed to analyze missense nsSNPs in the present study is shown in Fig. 1.

Figure 1. Schematic representation of in silico workflow to analyze deleterious missense nsSNPs in human GABRA6 gene.

2.1. Data mining

The missense nsSNPs for the human GABRA6 gene were retrieved from the National Center for Biotechnology Information (NCBI) dbSNP database (Sherry et al. 2001) (https://www.ncbi.nlm.nih.gov/snp/). Moreover, the protein sequence encoded by the human GABRA6 gene was obtained from the UniProt database (Bateman et al. 2021) (https://www.uniprot.org/).

2.2. Predicting deleterious missense nsSNPs

To predict deleterious missense nsSNPs, in silico tools including SIFT, PolyPhen-2, PROVEAN, and Condel were utilized.

SIFT (Sorting Intolerant from Tolerant): An algorithm that uses sequence homology and the physical properties of amino acids to predict an amino acid substitution effect on protein function. The algorithm considered evolutionarily conserved regions to be relatively intolerant to mutations, and therefore, amino acid substitution in these regions is more likely to be expected to affect protein function (Sim et al. 2012). In this study, a faster version of SIFT called SIFT 4G (SIFT for genomes) (Vaser et al. 2016) is used (The package is available at https://github.com/rvaser/sift4g). SIFT takes a protein sequence and a list of substitutions which is searched against a protein database for homology searching. UniRef90 (16 February 2021 release) (Suzek et al. 2007) is used as the protein database. SIFT score ranges from 0 to 1. The amino acid substitution is considered deleterious if the score < 0.05, while the score ≥ 0.05 is considered tolerated.

PolyPhen-2 (Polymorphism Phenotyping v2): A tool that predicts the structural and functional impacts of an amino acid substitution in human proteins. PolyPhen-2 (The package is available at http://genetics.bwh.harvard.edu/pph2/dokuwiki/downloads) uses a Naı̈ve Bayes classifier which is trained by supervised learning algorithms. HumDiv and HumVar datasets from the UniProt database were used to train PolyPhen-2 prediction models. The classifier which used HumDiv should be used for rare alleles assessment, natural selection analysis, and dense mapping of regions determined through genome-wide association studies. The classifier which used HumVar should be used for Mendelian disease recognition which requires recognizing mutations with severe effects from all the remaining human variation, including many mildly deleterious alleles (Adzhubei, Jordan, and Sunyaev 2013). PolyPhen-2 takes protein UniProt ID and substitutions as an input file and identifies homologs of the protein sequences by Basic Local Alignment Search Tool (BLAST) (Altschul et al. 1997) in the UniRef100 (27 February 2021 release) database. The PolyPhen-2 score ranges from 0 (neutral) to 1 (deleterious), and functional prediction is categorized into probably damaging (damaging with high confidence), possibly damaging (damaging with low confidence), and benign (benign with high confidence).

PROVEAN (Protein Variation Effect Analyzer): A tool based on an alignment scoring method which predicts functional effects of single amino acid substitutions, in-frame insertions, deletions, and multiple substitutions (Choi et al. 2012). PROVEAN (http://provean.jcvi.org/index.php) uses a delta alignment score (or delta score) to measure the effect of a variation. The protein sequence and a list of substitutions were submitted as inputs to the algorithm. PROVEAN’s default threshold score is -2.5. A score < -2.5 is considered deleterious, while a score > -2.5 is considered neutral (Choi and Chan 2015).

Condel (CONsensus DELeteriousness): A method to predict the outcome of non-synonymous SNPs, which is based on a weighted average of the normalized scores (WAS) of various tools. Condel (https://bbglab.irbbarcelona.org/fannsdb/) uses a consensus deleteriousness score that integrates the output of five tools (Logre (Clifford et al. 2004), MAPP (Stone and Sidow 2005), Mutation Assessor (Reva, Antipin, and Sander 2007), Polyphen2, and SIFT) into a unified classification (González-Pérez and López-Bigas 2011). Protein UniProt ID and substitutions were submitted as an input file to the algorithm. Condel score ranges from 0 to 1, which a higher score predicting SNPs as deleterious.

2.3. Predicting disease-related variants

To disease-related prediction of selected missense nsSNPs, SNPs&GO and PMut tools were used in this step.

SNPs&GO: A support vector machine (SVM) based tool for predicting disease-related SNPs using Gene Ontology (GO) terms (http://snps.biofold.org/snps-and-go) (Emidio Capriotti et al. 2013). Protein SwissProt code and substitutions were submitted as inputs to the algorithm. When the input sequence is a SwissProt code, GO terms used for the prediction will be retrieved automatically. The result consists of three different algorithms; SNPs&GO, PhD-SNP (E. Capriotti, Calabrese, and Casadio 2006), and PANTHER (Thomas and Kejariwal 2004). SNPs&GO score ranges from 0 to 1, if the score > 0.5 is considered a disease-related variation while a score ≤ 0.5 is considered benign.

PMut: A tool that predicts pathological properties of single amino acid substitutions. PMut (http://mmb.irbbarcelona.org/PMut) uses a neural network-based method that is trained using a manually curated database SwissProt (October 2016 release), includes 27,203 deleterious and 38,078 benign mutations for 12,141 proteins (López-Ferrando et al. 2017). Protein UniProt ID and substitutions were submitted as inputs to the algorithm. PMut score ranges from 0 to 1. A score > 0.5 is classified as pathological variation and a score < 0.5 is classified as neutral.

2.4. Predicting the molecular and phenotypic consequences of variants

Further investigation for the prediction of molecular and phenotypic consequences of selected missense nsSNPs was done through SNAP2, MutPred2, and SNPeffect4.0 tools.

SNAP2: A neural network-based classifier tool that predicts the functional effects of non-synonymous SNPs (https://rostlab.org/services/snap2web/). SNAP2 uses evolutionary information and structural features like solvent accessibility and predicted secondary structure to predict variants’ effects (Hecht, Bromberg, and Rost 2015). The protein sequence was submitted as input to the algorithm. SNAP2 score ranges from − 100 (strongly neutral) to + 100 (strongly effective).

MutPred2: A machine learning-based method which predicts the molecular and phenotypic consequences caused by amino acid substitutions as pathogenic or benign (The package is available at http://mutpred.mutdb.org/). It was trained using 53,180 deleterious and 206,946 putatively benign mutations from the Human Gene Mutation Database (HGMD) (Stenson et al. 2020), SwissVar database (Mottaz et al. 2010), dbSNP database, and inter-species pairwise alignments. MutPred2 takes protein UniProt ID and substitutions as an input file. MutPred2 gives two output scores, general score (g) and property score (pr). The g score demonstrates the pathogenicity of the substitution and ranges between 0 and 1, a higher score indicates a greater probability to be pathogenic. The pr score is the posterior probability of losing or gaining a certain property because of the substitution and it also ranges between 0 and 1, a higher score indicates more alteration of the property in the molecular mechanism of the disease (Pejaver et al. 2020).

SNPeffect 4.0: A database that uses 4 tools to predict molecular and structural phenotypic consequences of human protein-coding SNVs (https://snpeffect.switchlab.org/). SNPeffect 4.0 uses TANGO (Fernandez-Escamilla et al. 2004) for aggregation-prone regions prediction, WALTZ (Maurer-Stroh et al. 2010) for amyloidogenic regions prediction, LIMBO (Van Durme et al. 2009) for hsp70 chaperone-binding sites prediction, and FoldX (Schymkowitz et al. 2005) for analyzing the possible impact of protein stability (De Baets et al. 2012). Protein UniProt ID and substitutions were submitted as inputs to the algorithm. Mutations can increase (dTANGO, dWALTZ, and dWALTZ > 50), decrease (dTANGO, dWALTZ, and dWALTZ < -50), or have no effect (dTANGO, dWALTZ, and dWALTZ between − 50 and 50), aggregation propensity, amyloid propensity, and chaperone binding.

2.5. Predicting the effect of variants on protein stability

Analyzing missense nsSNPs effects on protein stability was accomplished by I-Mutant3.0 and MUpro tools, through Gibbs free energy change calculation of mutant protein and wild-type protein.

I-Mutant: A support vector machine (SVM) based web server tool, which predicts protein stability changes upon single amino acid substitutions (https://folding.biofold.org/i-mutant/). I-Mutant estimates protein stability by calculating Gibbs free energy change (ΔΔG, kcal/mol), of the mutant protein and wild-type protein (ΔΔG = ΔG (mutant type) - ΔG (wild-type)) (Emidio Capriotti et al. 2005). This server can evaluate the protein stability change by the protein sequence or the protein structure. In this study, protein sequence was used. The protein sequence and substitutions were submitted as inputs to the algorithm. If 0 < ΔΔG considered to increase protein stability, while ΔΔG < 0 is considered decreases protein stability. I-Mutant 3.0 predicts which a single substitution has a low effect (− 0.5 ≤ ΔΔG ≤ 0.5 kcal/mol) or extensively increases protein stability (ΔΔG > 0.5 Kcal/mol) or extensively decreases protein stability (ΔΔG < -0.5 Kcal/mol).

MUpro: A tool based on support vector machine (SVM) and neural network methods, which is used for protein stability prediction resulting from single amino acid substitutions (http://mupro.proteomics.ics.uci.edu/) (Cheng, Randall, and Baldi 2006). MUpro also likes I-Mutant can predict the stability changes using protein sequence or protein structure. The protein sequence and substitutions were submitted as inputs to the algorithm. When the energy change (ΔΔG) is positive the mutation increases protein stability. Instead, if ΔΔG is negative the mutation decreases protein stability.

2.6. Evolutionary conservation analysis

For the identification of selected missense nsSNPs position in evolutionarily conserved regions, the ConSurf tool was utilized.

ConSurf: A tool used for evolutionarily conservation analysis of amino acids in a protein, DNA, or RNA to reveal important functional and/or structural areas (https://consurf.tau.ac.il/). ConSurf uses an empirical Bayesian method that estimates the conservation score from protein sequence (Ashkenazy et al. 2016). The protein sequence was submitted as input to the algorithm. The conservation score is grouped into 9 grades, where 1 indicates most rapidly evolving areas, 5 indicates mildly evolving areas, and 9 indicates most evolutionary conserved areas.

2.7. Predicting protein secondary structure

Secondary structures of selected missense nsSNPs were predicted by SOPMA and PSIPRED tools.

SOPMA: A web version of the self-optimized prediction method (SOPM) (Geourjon and Deléage 1995), which use for secondary structure prediction of amino acids (https://npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl?page=/NPSA/npsa_sopma.html). The joint prediction using SOPMA, and a neural network method (PHD) correctly predicts 82.2% of residues for 74% of co-predicted amino acids. The protein sequence was submitted as input to the algorithm.

PSIPRED: A highly accurate method that predicts the secondary structure of proteins (http://bioinf.cs.ucl.ac.uk/psipred/). PSIPRED (McGuffin, Bryson, and Jones 2000) includes two feed-forward neural networks that perform an analysis of the output obtained from PSI-BLAST (Altschul et al. 1997). The protein sequence was submitted as input to the algorithm.

2.8. Predicting structural effects of variants

Prediction of structural changes such as size, charge, and hydrophobicity-value between mutant and wild residues was done by the HOPE project.

HOPE (Have (y)Our Protein Explained): A web tool to predict the structural effects of single amino acids substitutions on proteins (https://www3.cmbi.umcn.nl/hope/). HOPE uses 3D structures of the proteins that are available in the UniProt database and Distributed Annotation (DAS) servers to simulate the structural features of mutations on the wild-type protein (Venselaar et al. 2010). However, if necessary, the HOPE server can build homology models independently. HOPE sever also predicts structural changes between mutant and wild residues. The protein sequence was submitted as input to the algorithm.

2.9. 3D structure modeling and visualizing of natives and variants structures

The 3D structure of the protein is crucial to annotate the protein function of native and mutant structures. Identification of domain regions and locating missense nsSNPs position in domains were done through the Pfam database (Mistry et al. 2021) search (http://pfam.xfam.org/). Then, generating 3D structures of each native domain using I-TASSER (Iterative Threading ASSEmbly Refinement) was accomplished separately (https://zhanglab.dcmb.med.umich.edu/I-TASSER/). I-TASSER is a method for predicting protein structure that consists of 3 steps. In the first step, after the domain sequence was submitted to the I-TASSER server as input, a meta-threading program called LOMETS (Wu and Zhang 2007) retrieve structural templates from the PDB library (Berman et al. 2000). In the second step, full-length models are constructed by assembling well-aligned continuous fragments excised from the PDB templates, with unaligned regions structures built by ab-initio modeling based on replica-exchange Monte Carlo simulations (Zhang, Kolinski, and Skolnick 2003). I-TASSER simulations generate thousands of conformations named decoys. I-TASSER uses the SPICKER program (a clustering algorithm) (Zhang and Skolnick 2004) to cluster all decoys, through clustering, SPICKER identifies low free-energy states. Eventually, the functional annotations are obtained by matching the predicted structure models with known proteins in the BioLiP function library (Yang, Roy, and Zhang 2013a). Up to the top five models will predict by I-TASSER. The quality of predicted models is estimated by a confidence score or C-score. C-score is typically in the range of [-5,2], wherein a higher value signifies a better-quality model; when C score > − 1.5, indicates the model has a correct fold. I-TASSER takes protein sequence as input, and up to five models are predicted as output. In cases in which the output final models are less than five indicating the similarity of the top templates identified via LOMETS, and cause converging of I-TASSER simulations. In such cases, the C-score is usually high and predicts a high-quality model (Roy, Kucukural, and Zhang 2010; Yang et al. 2014; Yang and Zhang 2015).

Structural validation for the predicted models was carried out by PROCHECK and ERRAT. PROCHECK by Psi/Phi Ramachandran plot analysis assess the stereochemical quality of a protein structure (Laskowski et al. 1993). ERRAT is another protein structural verification tool, which demonstrates the reliability of the model. For the best quality of the predicted model, the ERRAT value must be over 80% (Colovos and Yeates 1993). PROCHECK and ERRAT are available at Structural Analysis and Verification Server (SAVES) (https://saves.mbi.ucla.edu/). PDB file was submitted as input to the server for both tools. Modeling of mutated protein structures was performed using the mutagenesis feature in the PyMol (DeLano 2002) visualization tool (2.3.0 version) using the predicted native type by I-TASSER as a reference.

2.10. Predicting ligand binding sites

The FTSite and COACH tools were used to predict the presence of selected missense nsSNPs in ligand-binding sites.

FTSite: An accurate computational method for ligand binding site prediction, which is based on experimental evidence with experimental accuracy of 94% (https://ftsite.bu.edu/) (Ngan et al. 2012). PDB file was submitted as input to the algorithm. FTSite identifies three potential binding sites in proteins.

COACH (COnsensus approACH): A tool for ligand binding site prediction, which is based on two comparative methods, TM-SITE (binding-specific substructure comparison) and S-SITE (sequence profile alignment) which identify ligand binding templates from the BioLiP library. These two methods in combination with COFACTOR (Roy, Yang, and Zhang 2012), ConCavity (Capra et al. 2009), and FINDSITE (Brylinski and Skolnick 2008) methods generate final predictions (Yang, Roy, and Zhang 2013b). After the protein sequence was submitted to the I-TASSER for generating 3D structures, then it will be fed into the COACH pipeline for predicting ligand binding sites.

2.11. Molecular dynamics simulation in water

Molecular dynamics simulations were performed using the GROMACS (2021.2 version) to reveal changes at the atomic level in different time scales for wild and mutant types (Abraham et al. 2015). The Amber ff99SB protein force field was used for simulations. A cubic box was formed by extending 1 nm on each side of the protein. The system was implemented by adding TIP3P water model molecules and neutralizing them with Cl − ions. The energy minimization step was performed using the steepest descent algorithm before the simulation. Long-range electrostatic interactions were modeled by the particle mesh Ewald (PME) method. And short-range electrostatic and van der Waals interactions cut-off radius was set to 1 nm. Periodic boundary conditions were maintained to eliminate surface effects. Leap-frog algorithm was used for simulations to integrate Newton’s equations of motion. A time step interval of 2 fs was used for all simulations. To constrain bonds involving hydrogens used the Lincs algorithm. Minimized systems should be equilibrated to reach the desired temperature and pressure before the real dynamics start. Equilibration is conducted in two ensembles, NVT (constant Number of particles, Volume, and Temperature) and NPT (constant Number of particles, Pressure, and Temperature). The system temperature and pressure were coupled to a v-rescale thermostat with a time constant of 0.1 ps at 300 K, and a Parrinello-Rahman barostat with a time constant of 2 ps at 1 bar. We first conducted a 200-ps NVT equilibration and then a 1000-ps NPT equilibration. The molecular dynamics simulation of each equilibrated system of the wild and mutant types was run for 100 ns. At last, the MDS trajectory files analyses were calculated by GROMACS build-in programs to get RMSD (root mean square deviations), RMSF (root mean square fluctuations), Rg (radius of gyration), SASA (solvent accessible surface area), and DSSP (definition secondary structure of the protein) analyses. RMSD, RMSF, Rg, and SASA of Cα atoms in the trajectories were evaluated. The changes in the secondary structure of the protein during molecular dynamics simulations were also calculated.

3.1. Missense nsSNPs datasets

For the biological consequences study of missense nsSNPs, we retrieved 451 missense nsSNPs with 386 rsIDs in the human GABRA6 gene mapped to NM_000811.3 RefSeq (S1 Table) from the dbSNP database. And the protein sequence of the GABRA6 gene (with Q16445 ID) was retrieved from the UniProt database.

3.2. Prediction of most deleterious missense nsSNPs

Four in silico tools including SIFT, PolyPhen-2, PROVEAN, and Condel were used to predict deleterious missense nsSNPs in this study. Out of 451 missense nsSNPs, 303 were predicted as deleterious and 148 as tolerated by SIFT, 306 were predicted as deleterious and 145 as neutral by HumDiv PolyPhen-2 (143 benign, 49 possibly damaging, and 259 probably damaging), 290 were predicted as deleterious and 161 as neutral by HumVar PolyPhen-2 (157 benign, 55 possibly damaging, and 239 probably damaging), 225 were predicted as deleterious and 226 as neutral by PROVEAN, and 298 were predicted as deleterious and 153 as neutral by Condel (S2 Table and Fig. 2).

Figure 2. Distribution of predicted deleterious and neutral missense nsSNPs by 4 tools for human GABRA6 gene. Deleterious nsSNPs shown in dark blue and neutral nsSNPs in light blue.

Deleterious missense nsSNPs were predicted with default scores for each tool. To obtain the high confidence deleterious missense nsSNPs, high scores across the tools are considered. SIFT score = 0, PolyPhen-2 score > 0.99, PROVEAN score < -9, and Condel score > 0.7 are used. After integration of the scores, 3 missense nsSNPs in the GABRA6 gene (Table 1) were obtained which are used for further analysis.

Table 1

The most deleterious missense nsSNPs by SIFT, PolyPhen-2, PROVEAN, and CONDEL tools in the *GABRA6* gene.
rs ID	Transcript change	Substitution	SIFT4G Score	PolyPhen-2 Score		PROVEAN Score	CONDEL Score
rs ID	Transcript change	Substitution	SIFT4G Score	HumDiv probability	HumVar probability	PROVEAN Score	CONDEL Score
rs1428649051	c.260G > C	W87S	0	1	1	-12.906	0.709363974
rs1199782347	c.334T > C	W112R	0	1	1	-12.837	0.735214831
rs1317373536	c.928T > C	C310R	0	1	1	-10.964	0.764688222

3.3. Prediction of disease-related variants

To get more accurate results selected missense nsSNPs were analyzed by SNPs&GO and PMut to predict disease-related nsSNPs. SNPs&GO's result consists of three different algorithms; SNPs&GO, PhD-SNP, and PANTHER. Three selected missense nsSNPs in the human GABRA6 gene were submitted to SNPs&GO and PMut to analyze disease-related nsSNPs. The output of SNPs&GO, PhD-SNP, PANTHER, and PMut predicted that all three selected missense nsSNPs as disease-related (Table 2).

Table 2

List of disease-related missense nsSNPs of *GABRA6* by SNPs&GO and PMut tools.
Substitution	SNPs&GO									PMut
	SNPs&GO			PhD-SNP			PANTHER			PMut
	Prediction	RI	Probability	Prediction	RI	Probability	Prediction	RI	Probability	Prediction	Score
W87S	Disease	8	0.919	Disease	8	0.922	Disease	10	0.999	Disease	0.8214
W112R	Disease	9	0.93	Disease	9	0.937	Disease	10	0.979	Disease	0.8214
C310R	Disease	9	0.949	Disease	9	0.97	Disease	9	0.974	Disease	0.8627
RI: Reliability Index

3.4. Prediction of the molecular and phenotypic consequences of variants

To predict the molecular and phenotypic consequences of selected missense nsSNPs, we further investigated nsSNPs through SNAP2, MutPred2, and SNPeffect4.0 tools. The SNPeffect tool predicts SNPs’ consequence on aggregation-prone regions by TANGO, amyloid-forming regions by WALTZ, and hsp70 chaperone binding sites by LIMBO. The output of SNAP2 of three selected missense nsSNPs in the GABRA6 gene predicted that all three selected missense nsSNPs have damaging effects on protein structure. Since the score > 0.5 suggests pathogenicity in MutPred2, the output of MutPred2 showed that all three selected missense nsSNPs in the GABRA6 gene have damaging functional and structural effects (Detailed information in S3 Table). According to the SNPeffect results, TANGO and WALTZ analysis revealed that the C310R variant (dTANGO equals 302.70 and dWALTZ equals − 553.68) increases the aggregation tendency of the protein and decreases the amyloid propensity of the protein. However, none of the variations does affect the chaperone binding tendency of the protein (Table 3).

Table 3

List of analyzed missense nsSNPs of the *GABRA6* gene by SNAP2, MutPred2, and SNPeffect4.0 tools.
Substitution	SNAP2		MutPred2 score	SNPeffect4.0
Substitution	Prediction	Score	MutPred2 score	dTANGO	dWALTZ	dLIMBO
W87S	effect	87	0.953	0.00	-0.07	0.00
W112R	effect	87	0.948	0.00	-0.12	0.00
C310R	effect	87	0.964	302.70	-553.68	0.00

(d in dTANGO, dWALTZ, and dLIMBO is the score change between the mutant protein and wild-type protein)

3.5. Prediction of the effect of variants on protein stability

The effects of selected missense nsSNPs on protein stability were analyzed by I-Mutant3.0 and MUpro tools by calculating Gibbs free energy change of mutant protein and wild-type protein. I-Mutant3.0 analyses showed W87S and W112R variants decreased protein stability, and the C310R variant increased protein stability. While all three variants by Mupro prediction showed decreasing protein stability. W87S and W112R variants have ΔΔG values < -1 kcal/mol in both tools which are expected to alter the function and structure of the protein (Table 4).

Table 4

List of missense nsSNPs of the *GABRA6* gene which was analyzed for protein stability by I-Mutant3.0 and MUpro tools.
Substitution	I-Mutant3.0			Mupro
Substitution	Prediction	ΔΔG	RI	Prediction	ΔΔG
W87S	Decrease	-1.69	9	Decrease	-1.0174971
W112R	Decrease	-1.44	9	Decrease	-1.2580621
C310R	Increase	-0.32	2	Decrease	-0.85748451
RI: Reliability Index

3.6. Evolutionary conservation analysis

To analyze the evolutionary conservation of selected missense nsSNPs, the ConSurf tool was used, which grouped amino acids based on conservation scores in 9 grades. As shown in Fig. 3, ConSurf results indicated that W87 and W112 residues located in highly conserved regions with conservation scores of 9 and predicted as buried residues and have a structural impact on the protein, and C310 predicted as buried residue with conservation scores of 8.

Figure 3. ConSurf analysis of GABRA6 gene residues. The black boxes indicate the most deleterious missense nsSNPs.

3.7. Prediction of protein secondary structure

Amino acids’ secondary structure corresponding to the proteins was predicted by SOPMA and PSIPRED tools. SOPMA prediction showed the distributions of alpha helix, extended strand, beta-turn, and random coil in proteins. PSIPRED predicted the distributions of strands, helices, and coils in proteins and validated the secondary structure of the proteins. SOPMA results showed 135 alpha helices (29.80%), 104 extended strands (22.96%), 11 beta turns (2.43%), and 203 random coils (44.81%) in the predicted secondary structure. In the 3 amino acid residues that correspond to selected missense nsSNPs, W87 and W112 are located in random coils, and C310 in alpha helices (Fig. 4). PSIPRED results indicated that in the 3 amino acid residues, W87 is located in strands, Y186 in coils, and C310 in helices (S1 Fig).

Figure 4. SOPMA analysis of GABRA6 gene residues. The black boxes indicate the most deleterious missense nsSNPs.

3.8. Prediction of structural effects of variants

HOPE was used to predict structural changes including size, charge, and hydrophobicity value between mutant and wild residues. Out of 3 selected missense nsSNPs in the GABRA6 gene, Ser and Arg residues, respectively in W87S and W112R mutants, are smaller than Trp residue in wild-types forms, while Arg residue in C310R mutant is bigger than the Cys residue in wild-type form. The charge of Trp and Cys residues respectively in W112R and C310R mutants were neutral, then turn to Arg residues with a positive charge. There is not any significant charge change in the W87S mutant. The wild-type residues in all three nsSNPs are more hydrophobic than the mutant residues. (S4 Table). According to the project HOPE results, when the mutant residue is smaller than the wild-type residue, it might lead to loss of interactions in protein structure, and if the mutant residue is bigger than the wild-type residue, it might lead to bumps in protein structure. When the charge of the wild-type residue is lost, it might cause a loss of interactions with other molecules or residues. If the mutation introduces a charge, it might cause the repulsion of ligands or other residues with the same charge. When the hydrophobicity of the wild-type residue is lost or decreased, the hydrophobic interactions will be lost either in the core of the protein or on the surface. If the mutation introduces a more hydrophobic residue, it might cause a loss of hydrogen bonds and/or disturb correct folding.

3.9. 3D structure modeling and visualizing of wild and mutant structures

The Pfam server was used to identify domain regions in the GABRA6 gene and locate selected nsSNPs positions in different domains. In this study, we just need domains that involved deleterious missense nsSNPs (Table 5). Generating 3D structure models was performed by the I-TASSER server. Predicted models’ validation was done by PROCHECK and ERRAT servers. PROCHECK predicted the quality of wild models by Ramachandran plot analysis. The ERRAT indicated the overall quality factor of the predicted models. Visualization of wild and mutant protein structures was performed by PyMol software.

Table 5

List of domains that involved selected deleterious missense nsSNPs in the *GABRA6* gene.
Protein name	Domain name	Start residue	Stop residue	Domain length	Involved nsSNPs
GBRA6_HUMAN	Neurotransmitter-gated ion-channel ligand-binding domain	32	240	209 aa	W87S and W112R
GBRA6_HUMAN	Neurotransmitter-gated ion-channel transmembrane domain	247	399	153 aa	C310R

Pfam reported two domains in the GABRA6 gene, including the neurotransmitter-gated ion-channel ligand-binding domain (32–240), and a neurotransmitter-gated ion-channel transmembrane domain (247–399). Out of three selected missense nsSNPs for further analysis, two (W87S and W112R) are located in the ligand-binding domain, and one (C310R) is located in the transmembrane domain.

3D structure prediction of wild types for the ligand-binding domain, and the transmembrane domain in the GABRA6 gene was modeled by I-TASSER. For the ligand-binding domain, one model predicted a C-score of 1.16, which was a high-quality model. And in the transmembrane domain, out of 5 predicted models, we selected the first model with the C-score of -1.62 (S2-a Fig and S3-a Fig). Validation of predicted models’ quality was done by PROCHECK and ERRAT servers. The result of PROCHECK for the ligand-binding domain model showed, 79.7% of residues in most favored regions, 18.7% in additional allowed regions, 1.6% in generously allowed regions, and 0% in disallowed regions (S2-b Fig). The ERRAT result showed that the overall quality factor for the predicted model was 87.940 (S2-c Fig). PROCHECK’s result for the transmembrane domain model showed, 77.3% of residues in most favored regions, 12.8% in additional allowed regions, 7.1% in generously allowed regions, and 2.8% in disallowed regions (S3-b Fig). The result of ERRAT showed that the overall quality factor for the predicted model was 83.916 (S3-c Fig).

Finally, the predicted models were used to model mutant types by using the mutagenesis feature in PyMol. Structural models for wild-types and deleterious missense nsSNPs in the ligand-binding domain (W87S and W112R) and the transmembrane domain (C310R), are shown in Fig. 5 to Fig. 7.

Figure 5. (W87S)

The amino acid Tryptophan (green) changed to Serine (red) at position 87 in the ligand-binding domain. Visualization was done by PyMol software and HOPE result.

Figure 6. (W112R)

The amino acid Tryptophan (green) changed to Arginine (red) at position 112 in the ligand-binding domain. Visualization was done by PyMol software and HOPE result.

Figure 7. (C310R)

The amino acid Cysteine (green) changed to Arginine (red) at position 310 in the transmembrane domain. Visualization was done by PyMol software and HOPE result.

3.10. Prediction of ligand binding sites

To determine the presence of the selected deleterious missense nsSNPs in protein-binding regions of the GABRA6 gene, we employed the FTSite and COACH algorithms, which are protein-ligand docking tools. FTSite identified three potential binding sites in the ligand-binding domain of the protein.

FTSite predicted three potential binding sites in the ligand-binding domain. Site 1 contains 15 binding residues, site 2 contains 4 binding residues, and site 3 contains 12 binding residues (Fig. 8-a and S5 Table). Additionally, COACH predicted 20 amino acid residues as potential binding sites (S6 Table).

Similarly, in the transmembrane domain, FTSite predicted three potential binding sites. Site 1 contains 9 binding residues, site 2 contains 16 binding residues, and site 3 contains 6 binding residues (Fig. 8-b and S5 Table). COACH predicted 44 amino acid residues as potential binding sites (S6 Table).

Figure 8. Illustration of GABRA6, ligand-binding domain (A), and transmembrane domain (B), with ligand-binding site predictions. FTSite were annotated the ligand-binding pockets. Site 1 illustrated in red, while site 2 and 3 illustrated in green and blue, respectively.

3.11. Molecular dynamics simulation of WT and mutant types

To comparatively study the conformational changes of the WT and mutant types in physiological environments, we performed 100 ns molecular dynamic simulations for each domain. Various parameters, such as root mean square deviations (RMSD), root mean square fluctuations (RMSF), the radius of gyration (Rg), and solvent accessible surface area (SASA), were analyzed using the time-dependent function of molecular dynamic simulation. These parameters were calculated for Cα atoms during the molecular dynamics simulations, with reference to their WT structures.

The RMSD analysis of the WT and mutant types revealed significant deviations in their structural stability. The W87S and W112R mutants in the ligand-binding domain exhibited similar RMSD values. Compared to the WT structure, the mutant types in the ligand-binding domain showed higher fluctuation, as depicted in Fig. 9-a. The C310R mutant in the transmembrane domain also deviated from the WT structure, with the mutant type displaying lower fluctuation than the WT, as shown in Fig. 9-b. The average RMSD values for the WT, W87S, and W112R are 1.88 Å, 2.13 Å, and 2.38 Å, respectively. Furthermore, the WT and C310R mutant have average RMSD values of 6.18 Å and 4.89 Å, respectively.

Figure 9. Analysis of RMSD values of the domain backbone of WT and mutant structures, ligand-binding domain (A), and transmembrane domain (B) at 100 ns simulation. The ordinate is RMSD (Å), and the abscissa is time (ns).

The RMSF analysis of each residue illustrates the effect of mutations on their dynamics. The analysis of RMSF values revealed significant differences in fluctuation between the WT and mutant structures in the N-terminal region of the W87S mutant and the N-terminal region as well as positions 90–99 in the W112R mutant in the ligand-binding domain after 100 ns of molecular dynamic simulation (Fig. 10-a). Similarly, in the transmembrane domain, the analysis of RMSF values indicated significant differences in fluctuation between the WT and mutant structures in positions 38–57, 77–92, and the C-terminal region of the C310R mutant (Fig. 10-b). The RMSF plots show that residues in positions 1–22 and 65–100 of the W112R mutant in the ligand-binding domain, as well as residues in positions 39–53 and 77–99 of the C310R mutant in the transmembrane domain, exhibit a relatively flexible region compared to other residues. Additionally, the highest residual fluctuation is observed at positions 1 (8.01 Å) and 2 (5.96 Å) in the W87S mutant, positions 1 (6.90 Å) and 2 (4.84 Å) in the W112R mutant, and positions 85 (7.95 Å) and 86 (7.71 Å) in the C310R mutant, when compared to their respective WT structures.

Figure 10. Analysis of RMSF values of the domain backbone of WT and mutant structures, ligand-binding domain (A), and transmembrane domain (B) over the entire simulation. The ordinate is RMSF (Å), and the abscissa is residue.

Based on the Rg analysis of the WT and mutant structures, it is evident that the W112R mutant exhibits a higher average Rg value (19.83 Å) in the ligand-binding domain compared to the WT (19.48 Å) and W87S mutant (19.41 Å), as depicted in Fig. 11-a. The WT and W87S mutant structures in the ligand-binding domain display similar average Rg values. Conversely, in the transmembrane domain, the C310R mutant demonstrates a significantly lower average Rg value (16.31 Å) compared to its WT structure (17.03 Å), as illustrated in Fig. 11-b. This suggests a potential decrease in the flexibility of the C310R mutant, and interestingly, the C310R mutant appears to deviate from its Rg value after 30 ns.

Figure 11. Analysis of Rg of the domain backbone of WT and mutant structures, ligand-binding domain (A), and transmembrane domain (B) over the entire simulation. The ordinate is Rg (Å), and the abscissa is time (ns).

Furthermore, the SASA analysis reveals that in the ligand-binding domain, the W112R mutant exhibits a higher average SASA value (11663.50 Å²) than the WT (11468.27 Å), while the W87S mutant displays a lower average SASA value (11370.25 Å²) than the WT, as shown in Fig. 12-a. In the transmembrane domain, the C310R mutant exhibits a lower average SASA value (8796.62 Å²) compared to its WT (8936.70 Å²), as depicted in Fig. 12-b. Since a higher SASA value indicates protein expansion, it can be inferred that in the ligand-binding domain, the WT and W87S mutant are more stable than the W112R mutant, and the W87S mutant is more stable than the WT. In the transmembrane domain, the C310R mutant is more stable than its WT structure.

Figure 12. Analysis of SASA of the domain backbone of WT and mutant structures, ligand-binding domain (A), and transmembrane domain (B) over the entire simulation. The ordinate is SASA (Å²), and the abscissa is time (ns).

To further investigate the structural changes resulting from substitutions in the GABRA6 gene in each mutant type, the number of different secondary structures present in the mutated types was calculated and compared with the wild types, as presented in Supplementary Figs. 4 to 10. Additionally, the contribution of different secondary structures in the protein structure during the simulation is summarized in Table 6. Moreover, to understand the secondary structural changes of the mutant types, the DSSP parameter was calculated during the simulation, as shown in Figs. 13 and 14.

Figure 13. Variation in secondary structure elements for WT (A), W87S (B), and W112R (C) structures with respect to time, in the ligand-binding domain.

Figure 14. Variation in secondary structure elements for WT (A), and C310R (B) structures with respect to time, in the transmembrane domain.

Table 6

Contribution of different secondary structures in protein structure during the simulation.
Protein	Structure	Coil	β-sheet	β-bridge	Bend	Turn	α-helix	5-helix	3₁₀-helix
WT ligand-binding domain	0.63	0.22	0.48	0.00	0.13	0.12	0.02	-	0.02
W87S mutant	0.62	0.23	0.46	0.01	0.12	0.14	0.01	-	0.03
W112R mutant	0.60	0.26	0.45	0.00	0.13	0.13	0.02	-	0.02
WT transmembrane domain	0.63	0.15	0.00	0.00	0.13	0.21	0.42	0.00	0.08
C310R mutant	0.58	0.17	0.00	0.00	0.15	0.17	0.41	0.02	0.08
(Structure = α-helix + β-sheet + β-bridge + turn)

The average coil percentage in the ligand-binding domain was determined to be 22.09%, 22.86%, and 25.85% for the WT, W87S mutant, and W112R mutant structures, respectively (S4-a Fig). In the transmembrane domain, the average coil percentage for the WT and C310R mutant structures was found to be 15.49% and 17.28%, respectively (S4-b Fig). Secondary structure analysis in the ligand-binding domain reveals a significant difference in fluctuation between the WT and W112R mutant structures in the coil region residues. Similarly, in the transmembrane domain, a notable difference in fluctuation in the coil region can be observed between the WT and C310R mutant structures after approximately 50 ns of molecular dynamics simulation.

In the ligand-binding domain, the average percentage of β-sheet was found to be 48.32%, 46%, and 44.54% for the WT, W87S mutant, and W112R mutant structures, respectively (S5-a Fig). However, in the transmembrane domain, the average β-sheet percentage was practically 0 (S5-b Fig). The secondary structure analysis in the ligand-binding domain demonstrates a considerable difference in fluctuations between the WT and mutant structures, particularly in the residues of the β-sheet region.

The average β-bridge percentages in the ligand-binding domain were determined to be 0.35% for the WT structure, 0.83% for the W87S mutant structure, and 0.28% for the W112R mutant structure (S6-a Fig). In the transmembrane domain, the β-bridge percentages were found to be 0.03% for the WT structure and 0.04% for the C310R mutant structure (S6-b Fig).

The average percentages of bend were found to be similar among the WT, W87S mutant, and W112R mutant structures in the ligand-binding domain, with values of 12.69%, 12.11%, and 12.65% respectively (S7-a Fig). In the transmembrane domain, the average bend percentages were 12.85% for the WT structure and 14.75% for the C310R mutant structure (S7-b Fig). Notably, the C310R mutant structure exhibited a significant increase in fluctuation in the bend region residues during the first 30 ns compared to the WT structure in the ligand-binding domain.

Regarding the average turn percentages, the ligand-binding domain showed values of 12.10% for the WT structure, 14.06% for the W87S mutant structure, and 13.43% for the W112R mutant structure (S8-a Fig). In the transmembrane domain, the turn percentages were determined to be 20.67% for the WT structure and 16.60% for the C310R mutant structure (S8-b Fig). Notably, there were significant differences in fluctuation between the WT and three mutant structures in both domains within the turn region residues.

In terms of average α-helix percentages, the ligand-binding domain exhibited values of 2% for the WT structure, 1.18% for the W87S mutant structure, and 1.62% for the W112R mutant structure (S9-a Fig). In the transmembrane domain, the average α-helix percentages were 42.47% for the WT structure and 41.37% for the C310R mutant structure (S9-b Fig). Notably, there was a considerable difference in fluctuation in the α-helix region residues between the WT and C310R mutant structures throughout the molecular dynamics simulation times in the transmembrane domain.

No 5-helix secondary structure was observed in the ligand-binding domain. However, in the transmembrane domain, the average percentages of 5-helix were 0.21% for the WT structure and 1.95% for the C310R mutant structure.

Lastly, the average 3₁₀-helix percentages in the ligand-binding domain were found to be 2.41% for the WT structure, 2.93% for the W87S mutant structure, and 1.60% for the W112R mutant structure (S10-a Fig). In the transmembrane domain, the average 3₁₀-helix percentages were 8.24% for the WT structure and 7.98% for the C310R mutant structure (S10-b Fig). Notably, a decrease in fluctuation was observed in the W112R mutant structure compared to the WT structure in the ligand-binding domain. Conversely, an increase in fluctuation was visible in the transmembrane domain of the C310R mutant structure compared to the WT structure.

The objective of this study was to identify the most detrimental missense nsSNPs associated with the human GABRA6 gene, which encodes the alpha-6 subunit of the \({\text{G}\text{A}\text{B}\text{A}}_{\text{A}}\) receptor, and to investigate their molecular, functional, and structural effects on the protein. It is widely observed that regions that are evolutionarily conserved tend to contain disease-related SNPs. According to the neutral theory of evolution proposed by Kimura in 1983, amino acids that differ between species or SNPs that do not occur in coding regions of a protein are usually not influenced by natural selection or are under low selection pressure (Kimura 1983). This suggests that such amino acid changes are generally tolerable and have minimal impact on protein function. Consequently, these amino acids are considered less mutable and more stable. However, nsSNPs occurring in coding regions of human genes can have phenotypic effects and may undergo natural selection. Amino acid substitution patterns can provide insights into selection patterns in protein-coding genes. An interesting example illustrating this is the distribution pattern of disease-related and non-pathogenic mutations in human genes. For instance, studies utilizing disease-related mutation data and phylogenetic analysis of multiple species have demonstrated that disease-associated substitutions (DAS) occur more frequently in evolutionarily conserved regions compared to regions subject to variation. Conversely, silent and polymorphic mutations exhibit the opposite trends, being randomly distributed. These patterns support the notion that protein conserved regions are subjected to evolutionary pressure as these amino acids are crucial for the normal functioning of a particular gene. Conversely, due to their random distribution, silent mutations have minimal impact on organisms and may not be influenced by natural selection (Shastry 2007). Therefore, non-synonymous SNPs in coding regions have been selected for analysis in this study.

Experimental investigation through wet laboratory techniques can be time-consuming, expensive, and resource-intensive. However, prior utilization of bioinformatics tools can determine the molecular, functional, and structural effects of SNPs on proteins, thereby reducing the time frame and economic costs of the preclinical stage, before further in vivo and in vitro analysis. To predict missense nsSNPs, a combination of SIFT, PolyPhen-2, PROVEAN, Condel, SNPs&GO, PMut, SNAP2, MutPred2, and SNPeffect4.0 tools were employed. Additionally, I-Mutant3.0 and MUpro tools were used for protein stability analysis, ConSurf for evolutionary conservation analysis, and SOPMA and PSIPRED for protein secondary structure analysis. Project HOPE was utilized to predict structural changes, such as size, charge, and hydrophobicity value, between mutant and wild residues. For 3D structure modeling, the Pfam server was employed for domain regions identification, the I-TASSER server for 3D structure model generation, and the PROCHECK and ERRAT servers for structural validation. Ligand binding sites prediction was conducted using the FTSite and COACH servers. Finally, to investigate conformational changes in physiological environments, molecular dynamic simulations were performed for each WT and mutant structure using GROMACS, with RMSD, RMSF, Rg, SASA, and DSSP values calculated.

4.1. Missense nsSNPs prediction

In this study, 3 most deleterious missense nsSNPs in the human GABRA6 gene coding region (Table 1), after analyzing by SIFT, PolyPhen-2, PROVEAN, and Condel tools were identified with high confidence scores, which are high-risk missense nsSNPs.

These selected missense nsSNPs were further analyzed using five tools to study the effects of mutations on the protein. The reasons why each computational method has different results are due to predict based on different databases and the use of different algorithms. SNPs&GO and PMut were used to predict disease-related nsSNPs; SNPs&GO's result consists of three different algorithms; SNPs&GO, PhD-SNP, and PANTHER. All 3 missense nsSNPs in the GABRA6 gene (W87S, W112R, and C310R, Table 2), were predicted as deleterious in all 4 algorithms. SNAP2, MutPred2, and SNPeffect4.0 tools were used to predict the molecular and phenotypic consequences of selected missense nsSNPs; SNPeffect tool predicts SNPs consequence by TANGO (predict aggregation-prone regions), WALTZ (predict amyloid forming regions), and LIMBO (predict hsp70 chaperone binding sites). Results of SNAP2 and MutPred2 of 3 selected missense nsSNPs in the GABRA6 gene showed that all 3 nsSNPs have damaging effects on protein structure and function. TANGO analysis showed that the C310R variant increases the aggregation tendency of the protein and WALTZ analysis showed that this variant decreases the amyloid propensity of the protein (Table 3).

4.2. Protein stability analysis

Protein stability is critical for the function, activity, and regulation of biomolecules. Missense mutations often impair the affected polypeptide propensity for the functional conformation folding and decrease the functional conformation stability; that leads to an increase in the proportion of mutant polypeptides that exist in a non-functional conformation, and the non-functional conformation is more prone to aggregation or degradation than the functional conformation. Any decreased stability and incorrect folding are the major consequences of deleterious missense variants (Bross et al. 1999). The effects of selected missense nsSNPs on protein stability were estimated by I-Mutant3.0 and MUpro tools by calculating Gibbs free energy change (ΔΔG) of the mutant protein and wild-type protein. The free energy of folding (ΔG) is used for measuring the thermodynamic stability of a protein and corresponds to the difference in free energy between the folded and unfolded protein state; both mutant protein and wild-type protein have their own ΔG values (Ferrer-Costa, Orozco, and De La Cruz 2002). When 0 < ΔΔG indicates increased stability in the mutant type, while ΔΔG < 0 shows decreased stability. Selected missense nsSNPs in the GABRA6 gene showed decreasing protein stability in the mutant types by both tools; except the C310R variant, which showed increasing protein stability in the mutant type by I-Mutant3.0 analysis. W87S and W112R variants in the GABRA6 gene have ΔΔG values < -1 kcal/mol in both tools which are expected can alter the function and structure of the protein.

4.3. Evolutionary conservation analysis and ligand binding sites prediction

Proteins conformational changes are essential for the function of proteins. In recent years, some important investigations have proposed the structure of complexes that exhibit conformational changes upon binding, and further understanding of the functions of these proteins. As more experimental work is done for characterizing the dynamics of protein interactions, it has become increasingly clear that proteins can exist in different conformational states. The flexibility within the protein regions allows it to adopt a new conformation and then bind ligands with different structures. The proteins' ability to adopt several structures allows functional diversity without relying on the sequence diversity evolution, which can greatly facilitate the potential to rapidly evolve new functions and new structures (Goh, Milburn, and Gerstein 2004). Many properties of proteins derive from two opposing factors, flexibility, and rigidity of their conformation. Protein molecules need flexibility and rigidity in order to function. Flexibility is required for the regulation and function, and rigidity is essential for the maintaining globular structure. In some proteins, there are flexible local sites that are highly conserved, this indicates the relevance of flexibility to the overall stability of proteins. So, flexibility and rigidity in protein conformation must be exquisitely balanced (Vihinen 1987).

ConSurf was used to identify selected missense nsSNPs position in evolutionarily conserved regions; and, FTSite and COACH algorithms were used to predict whether these selected deleterious missense nsSNPs are present in protein binding regions or not. ConSurf results for the human GABRA6 gene showed that two residues (W87 and W112) located in highly conserved regions (with scores 9) and predicted to have structural impacts on the protein, and another residue (C310) predicted as buried residue (with score 8). Further analysis using FTSite and COACH algorithms did not reveal any significant deleteriousness in selected missense nsSNPs located at binding sites in its two protein domains. All the missense nsSNPs which were predicted as highly conserved by the ConSurf server could be used in pathogenic SNPs screening, because neutral nsSNPs are more frequent in variable regions, whilst the deleterious nsSNPs are more common in conserved regions (Bond et al. 2010).

4.4. Protein secondary structure analysis

According to secondary structure predicting tools of proteins used in this study, 74.61% of amino acids sites of GABRA6 protein, are located in alpha helices and random coils (Fig. 4); which corresponds to the evidence from the literature that sites of pathogenic and polymorphic mutations are commonly found to be in alpha helices and random coils regions, and not frequently in beta-strand regions (Kucukkal et al. 2015). Out of 3 amino acid residues corresponding to GABRA6 selected missense nsSNPs, SOPMA results indicated one residue (C310) located in alpha helices, and other two residues (W87 and W112) in random coils; and PSIPRED results indicated W87 residue located in strands, W112 residue in coils, and C310 residue in helices (S1 Fig).

4.5. Project HOPE analysis

Mutations in or near specific amino acids that contribute to the conformation of functional spatial are very likely to cause pathology. Missense mutations can cause amino acids to be substituted, resulting in changes in amino acid size, charge, and hydrophobicity, which can affect protein folding and interactions. According to the project HOPE analysis, mutation-related changes either lead to loss of interaction or structural disorder, especially in the transmembrane domain. Furthermore, the introducing or loss of hydrophobicity or charge could cause misfolding, loss of interactions, or repulsion. As Ser and Arg residues in W87S and W112R mutants are smaller than Trp residue in wild-type forms, it may lead to loss of interactions in protein structure. Also, Arg residue in the C310R mutant is bigger than the Cys residue in the wild-type form which may lead to bumps in protein structure. The charge of Trp and Cys residues in W112R and C310R mutants were neutral, then turned to Arg residues with a positive charge, whenever the mutations introduce any charges, it might cause repulsion of ligands or other residues with the same charge. In the W87S mutant there is not any charge change. The hydrophobicity in these three nsSNPs decreased than the wild-type forms, which might cause hydrophobic interactions lost in the core of the protein or on the surface. (S4 Table).

4.6. 3D structure modeling

Pfam server was utilized for domain region identification in GABRA6 protein; and, then selected missense nsSNPs positions were located in different domains (Table 5). After generating 3D structure models by I-TASSER server and quality validation of predicted models by PROCHECK and ERRAT servers, PyMol software was used to visualize the 3D structural changes in the wild and mutant protein structures; the differences in the structures between the wild type and mutant type which visualized by PyMol software and 2D schematic structures from project HOPE proved the pathogenic impact of these deleterious missense nsSNPs on the protein (Figs. 5_7).

4.7. Substitution effects on structure via molecular dynamics simulations

We conducted molecular dynamics simulations to investigate the potential deleterious effects of selected mutations on the protein. Specifically, we performed 100 ns simulations to examine the impact of these mutations on the structural dynamics of the protein.

In the case of the W87S mutant, the RMSD value was consistently higher during most of the 100 ns simulation, particularly in the first 50 ns. This suggests that this mutant is likely to destabilize the protein structure compared to the WT. Although the RMSF value remained relatively constant in most regions, it was higher in certain areas, indicating a potential loss of stability in the protein structure. Additionally, the Rg value was slightly lower throughout the simulation, which may contribute to an increase in the stability of the protein structure. Furthermore, the SASA value was lower than both the WT and W112R mutant structures, indicating that the W87S mutant exhibits the highest level of stability among these variants.

For the W112R mutant, the RMSD value consistently remained higher throughout the entire simulation, especially in the first 20 ns and the last 30 ns. This suggests that this mutant is likely to destabilize the protein structure compared to the WT. Similarly, the RMSF value was higher in specific regions, indicating a potential loss of stability in the protein structure. Moreover, the Rg value was higher throughout the 100 ns simulation, suggesting a potential decrease in the stability of the protein structure. The SASA value of the W112R mutant was the highest among the WT and mutant structures, indicating that this mutant is likely to be the least stable among them.

In the case of the C310R mutant, the RMSD value consistently remained much lower throughout the simulation, suggesting that this mutant is more likely to stabilize the protein structure compared to the WT. The RMSF value showed variations, both higher and lower than the WT, indicating the possibility of both stability loss and gain in the protein structure. The Rg value was lower throughout the simulation, particularly after 30 ns, suggesting a potential gain in the stability of the protein structure. Additionally, the SASA value was lower than the WT structure, suggesting that this mutant may cause a loss of stability in the protein structure.

To understand the disruption in secondary structures during the simulation, we conducted secondary structure analyses. We observed significant differences in coil, β-sheet, bend, turn, α-helix, and 3₁₀-helix regions for the mutant types compared to the WT structures, as shown in Figs. 13 and 14. The results revealed that the total number of secondary structures involved in the formation of secondary structures is higher in the WT compared to all mutant types.

Specifically, compared to the WT, the W87S and W112R mutants exhibited a decrease of 2.32% and 3.78% in β-sheets and turns, respectively, and an increase of 1.96% and 1.33% in β-sheets and turns, respectively. The W112R mutant also showed a reduction of 3.76% in coils. On the other hand, the C310R mutant displayed an increase of 1.79% in coils, 1.9% in bends, and 1.74% in 5-helices, along with a decrease of 4.07% in turns, compared to the WT. These alterations in secondary structures in the mutant types can influence the tertiary structure of the protein and potentially lead to changes in its three-dimensional structure or instability.

Furthermore, the analysis revealed that in the W87S and W112R mutants, the amino acids removed from the secondary structure (Try) resulted in a reduction of β-sheets and conversion to turn and coil structures. Conversely, in the C310R mutant, there was no significant preference for conversion to any specific structure, as turns were reduced and converted to 5-helices, coils, and bends. Overall, among the secondary structures, the turn structure exhibited the most significant changes in the mutated types.

In conclusion, our molecular dynamics simulations suggest that the selected mutations have the potential to affect the stability and secondary structure of the protein. The W87S and W112R mutants are likely to destabilize the protein structure, while the C310R mutant may contribute to its stabilization. These findings provide valuable insights into the impact of these mutations on the protein's structural dynamics and can aid in further investigations and understanding of their functional implications.

This study aimed to investigate the functional and structural effects of missense nsSNPs in the human GABRA6 gene using advanced in-silico prediction tools. While experimental methods are more reliable for discerning deleterious nsSNPs, conducting repeated experiments on all pathogenic mutations is time-consuming. The approaches employed in this study provide evidence of the diverse effects of mutations, which can aid in characterizing their pathogenicity. Notably, three high-risk missense mutations in the GABRA6 gene were identified, exhibiting the most detrimental impacts on protein structure and function. These nsSNPs hold potential as diagnostic markers or drug targets for epilepsy and other neuropsychiatric disorders. Consequently, future experimental designs can incorporate this in-silico data to analyze the biological context. To explore avenues for mitigating the effects of these three high-risk mutations in the mutated protein, it is suggested to conduct research, such as targeted mutagenesis in vitro or animal models, to modify the functional or structural consequences of the mutations.

Author Contribution

TM collected, analyzed, and interpreted the data and also drafted the manuscript. RM designed and revised the work. Both authors read and approved the final manuscript.

Abraham M, James, et al. Gromacs: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX. 2015;1–2:19–25.
Adzhubei I, Jordan DM, Sunyaev SR. 2013. Current Protocols in Human Genetics Predicting Functional Effect of Human Missense Mutations Using PolyPhen-2.
Altschul SF, et al. Gapped BLAST and PSI-BLAST: A New Generation of Protein Database Search Programs. Nucleic Acids Res. 1997;25(17):3389–402. https://doi.org/10.1093/nar/25.17.3389.
Ashkenazy H, et al. ConSurf 2016: An Improved Methodology to Estimate and Visualize Evolutionary Conservation in Macromolecules. Nucleic Acids Res. 2016;44(W1):W344–50.
De Baets, Greet, et al. SNPeffect 4.0: On-Line Prediction of Molecular and Structural Effects of Protein-Coding Variants. Nucleic Acids Res. 2012;40(D1):935–39.
Bateman A, et al. UniProt: The Universal Protein Knowledgebase in 2021. Nucleic Acids Res. 2021;49(D1):D480–89.
Berman HM, et al. The Protein Data Bank. Nucleic Acids Res. 2000;28(1):235–42. https://doi.org/10.1093/nar/28.1.235.
Bond LM, et al. Gene Repression by Minimal Lac Loops in Vivo. Nucleic Acids Res. 2010;38(22):8072–82.
Brejc Katjuša, et al. Crystal Structure of an ACh-Binding Protein Reveals the Ligand-Binding Domain of Nicotinic Receptors. Nature. 2001;411(6835):269–76.
Brooks-Kayal AR, Shelley JR. 2012. Regulation of GABAA Receptor Gene Expression and Epilepsy. Jasper’s Basic Mechanisms of the Epilepsies [Internet]. 4th edition.
Bross P, et al. Protein Misfolding and Degradation in Genetic Diseases. Hum Mutat. 1999;14(3):186–98.
Brylinski M. and Jeffrey Skolnick. 2008. A Threading-Based Method (FINDSITE) for Ligand-Binding Site Prediction and Functional Annotation. Proceedings of the National Academy of Sciences 105(1): 129–34. https://www.pnas.org/content/105/1/129.
Capra JA et al. 2009. Predicting Protein Ligand Binding Sites by Combining Evolutionary Sequence Conservation and 3D Structure. PLoS Computational Biology 5(12).
Capriotti E, Calabrese R, Casadio R. Predicting the Insurgence of Human Genetic Diseases Associated to Single Point Protein Mutations with Support Vector Machines and Evolutionary Information. Bioinformatics. 2006;22(22):2729–34.
Capriotti E et al. 2013. WS-SNPs&GO: A Web Server for Predicting the Deleterious Effect of Human Protein Variants Using Functional Annotation. BMC Genomics 14 Suppl 3(Suppl 3).
Capriotti E, Fariselli P, Calabrese R. and Rita Casadio. 2005. Predicting Protein Stability Changes from Sequences Using Support Vector Machines. Bioinformatics 21(SUPPL. 2): 54–58.
Cheng J, Randall A, and Pierre Baldi. Prediction of Protein Stability Changes for Single-Site Mutations Using Support Vector Machines. Proteins: Struct Function Genet. 2006;62(4):1125–32.
Choi Y et al. 2012. Predicting the Functional Effect of Amino Acid Substitutions and Indels. PLoS ONE 7(10).
Choi Y, Chan AP. PROVEAN Web Server: A Tool to Predict the Functional Effect of Amino Acid Substitutions and Indels. Bioinformatics. 2015;31(16):2745–47.
Clifford RJ, Edmonson MN, Nguyen C, Buetow KH. Large-Scale Analysis of Non-Synonymous Coding Region Single Nucleotide Polymorphisms. Bioinformatics. 2004;20(7):1006–14.
Colovos C, and Todd O. Yeates. Verification of Protein Structures: Patterns of Nonbonded Atomic Interactions. Protein Sci. 1993;2(9):1511–19.
DeLano WL. Pymol: An Open-Source Molecular Graphics Tool. CCP4 Newsl protein Crystallogr. 2002;40:82–92. http://www.ccp4.ac.uk/newsletters/newsletter40/11_pymol.pdf.
Dibbens LM, et al. The Role of Neuronal GABAA Receptor Subunit Mutations in Idiopathic Generalized Epilepsies. Neurosci Lett. 2009;453(3):162–65.
Van Durme J et al. 2009. Accurate Prediction of DnaK-Peptide Binding via Homology Modelling and Experimental Data. PLoS Comput Biol 5(8).
Ekins S, Mestres J, Testa B. In Silico Pharmacology for Drug Discovery: Methods for Virtual Ligand Screening and Profiling. Br J Pharmacol. 2007;152(1):9–20.
Fernandez-Escamilla A, Maria F, Rousseau J, Schymkowitz, and Luis Serrano. Prediction of Sequence-Dependent and Mutational Effects on the Aggregation of Peptides and Proteins. Nat Biotechnol. 2004;22(10):1302–6.
Ferrer-Costa C, Orozco M, De La X Cruz. Characterization of Disease-Associated Single Amino Acid Polymorphisms in Terms of Sequence and Structure Properties. J Mol Biol. 2002;315(4):771–86.
Fielden MR et al. 2002. 32 Critical Reviews in Toxicology In Silico Approaches to Mechanistic and Predictive Toxicology: An Introduction to Bioinformatics for Toxicologists.
Fisher RS, et al. ILAE Official Report: A Practical Clinical Definition of Epilepsy. Epilepsia. 2014;55(4):475–82.
Fritschy JM. 2008. Epilepsy, E/I Balance and GABAA Receptor Plasticity. Frontiers in Molecular Neuroscience 1(MAR).
Doss GP, C., et al. A Novel Computational and Structural Analysis of NsSNPs in CFTR Gene. Genom Med. 2008;2(1–2):23–32.
Geourjon C, Deléage G. Sopma: Significant Improvements in Protein Secondary Structure Prediction by Consensus Prediction from Multiple Alignments. Bioinformatics. 1995;11(6):681–84.
Goh C, Sing D, Milburn, and Mark Gerstein. Conformational Changes Associated with Protein-Protein Interactions. Curr Opin Struct Biol. 2004;14(1):104–9.
González-Pérez A, Nuria, López-Bigas. Improving the Assessment of the Outcome of Nonsynonymous SNVs with a Consensus Deleteriousness Score, Condel. Am J Hum Genet. 2011;88(4):440–49.
Hecht M, Bromberg Y, and Burkhard Rost. Better Prediction of Functional Effects for Sequence Variants. BMC Genomics. 2015;16(8):1–12.
Hernandez CC, Katharine N, Gurba N, Hu, Macdonald RL. The GABRA6 Mutation, R46W, Associated with Childhood Absence Epilepsy, Alters Α6β2γ2 and Α6β2δ GABA A Receptor Channel Gating and Expression. J Physiol. 2011;589(23):5857–78.
Hirose S. 2014. 213 Progress in Brain Research Mutant GABAA Receptor Subunits in Genetic (Idiopathic) Epilepsy. 1st ed. Elsevier B.V. http://dx.doi.org/10.1016/B978-0-444-63326-2.00003-X.
Kimura M. The Neutral Theory of Molecular Evolution. Cambridge University Press; 1983.
Kucukkal TG, Petukh M, Li L, and Emil Alexov. Structural and Physico-Chemical Effects of Disease and Non-Disease NsSNPs on Proteins. Curr Opin Struct Biol. 2015;32:18–24. http://dx.doi.org/10.1016/j.sbi.2015.01.003.
Laskowski RA, MacArthur MW, Moss DS, Thornton JM. PROCHECK: A Program to Check the Stereochemical Quality of Protein Structures. J Appl Crystallogr. 1993;26(2):283–91.
López-Ferrando Víctor, et al. PMut: A Web-Based Tool for the Annotation of Pathological Variants on Proteins, 2017 Update. Nucleic Acids Res. 2017;45(W1):W222–28.
Macdonald RL, Richard W, Olsen. GABAA Receptor Channels. Annu Rev Neurosci. 1994;17(1):569–602.
Mah JTL, Chia KS. A Gentle Introduction to SNP Analysis: Resources and Tools. J Bioinform Comput Biol. 2007;5(5):1123–38.
Mah JTL, Esther SH, Low, and Edmund Lee. In Silico SNP Analysis and Bioinformatics Tools: A Review of the State of the Art to Aid Drug Discovery. Drug Discovery Today. 2011;16(17–18):800–9. http://dx.doi.org/10.1016/j.drudis.2011.07.005.
Maurer-Stroh S, et al. Exploring the Sequence Determinants of Amyloid Structure Using Position-Specific Scoring Matrices. Nat Methods. 2010;7(3):237–42.
McGuffin LJ, Bryson K, and David T. Jones. The PSIPRED Protein Structure Prediction Server. Bioinformatics. 2000;16(4):404–5.
Mistry J, et al. Pfam: The Protein Families Database in 2021. Nucleic Acids Res. 2021;49(D1):D412–19.
Mottaz Anaïs, Fabrice PA, David AL, Veuthey, Yip YL. Easy Retrieval of Single Amino-Acid Polymorphisms and Phenotype Information Using SwissVar. Bioinformatics. 2010;26(6):851–52.
Myers CT, Mefford HC. Advancing Epilepsy Genetics in the Genomic Era. Genome Med. 2015;7(1):1–11. http://dx.doi.org/10.1186/s13073-015-0214-7.
Ngan CH, et al. FtSite: High Accuracy Detection of Ligand Binding Sites on Unbound Protein Structures. Bioinformatics. 2012;28(2):286–87.
Pejaver V et al. 2020. Inferring the Molecular and Phenotypic Impact of Amino Acid Variants with MutPred2. Nature Communications 11(1). http://dx.doi.org/10.1038/s41467-020-19669-x.
Reva B, Antipin Y. and Chris Sander. 2007. Determinants of Protein Function Revealed by Combinatorial Entropy Optimization. Genome Biol 8(11).
Roy A, Alper Kucukural, and, Zhang Y. I-TASSER: A Unified Platform for Automated Protein Structure and Function Prediction. Nat Protoc. 2010;5(4):725–38.
Roy A, Yang J, Zhang Y. COFACTOR: An Accurate Comparative Algorithm for Structure-Based Protein Function Annotation. Nucleic Acids Res. 2012;40(W1):471–77.
Scheffer IE, et al. ILAE Classification of the Epilepsies: Position Paper of the ILAE Commission for Classification and Terminology. Epilepsia. 2017;58(4):512–21.
Schymkowitz J et al. 2005. The FoldX Web Server: An Online Force Field. Nucleic Acids Research 33(SUPPL. 2): 382–88.
Shastry BS. SNPs in Disease Gene Mapping, Medicinal Drug Development and Evolution. J Hum Genet. 2007;52(11):871–80.
Sherry ST, et al. DbSNP: The NCBI Database of Genetic Variation. Nucleic Acids Res. 2001;29(1):308–11.
Sim N, Leng, et al. SIFT Web Server: Predicting Effects of Amino Acid Substitutions on Proteins. Nucleic Acids Res. 2012;40(W1):452–57.
Stenson PD, et al. The Human Gene Mutation Database (HGMD®): Optimizing Its Use in a Clinical Diagnostic or Research Setting. Hum Genet. 2020;139(10):1197–207. https://doi.org/10.1007/s00439-020-02199-3.
Stone EA, and Arend Sidow. Physicochemical Constraint Violation by Missense Substitutions Mediates Impairment of Protein Function and Disease Severity. Genome Res. 2005;15(7):978–86.
Suzek BE, et al. UniRef: Comprehensive and Non-Redundant UniProt Reference Clusters. Bioinformatics. 2007;23(10):1282–88.
Syvänen A-C. Accessing Genetic Variation: Genotyping Single Nucleotide Polymorphisms. Nat Rev Genet. 2001;2(12):930–42.
Thomas PD, and Anish Kejariwal. Coding Single-Nucleotide Polymorphisms Associated with Complex vs. Mendelian Disease: Evolutionary Evidence for Differences in Molecular Effects. Proc Natl Acad Sci USA. 2004;101(43):15398–403.
Vaser R, et al. SIFT Missense Predictions for Genomes. Nat Protoc. 2016;11(1):1–9. http://dx.doi.org/10.1038/nprot.2015-123.
Venselaar H et al. 2010. Protein Structure Analysis of Mutations Causing Inheritable Diseases. An e-Science Approach with Life Scientist Friendly Interfaces. BMC Bioinformatics 11.
Vihinen M. Relationship of Protein Flexibility to Thermostability. Protein Eng. 1987;1(6):477–80.
Wang J, et al. Epilepsy-Associated Genes. 2017;44:11–20.
Wu S, and Yang Zhang. LOMETS: A Local Meta-Threading-Server for Protein Structure Prediction. Nucleic Acids Res. 2007;35(10):3375–82.
Yang J, et al. The I-TASSER Suite: Protein Structure and Function Prediction. Nat Methods. 2014;12(1):7–8. http://dx.doi.org/10.1038/nmeth.3213.
Yang J, Roy A, and Yang Zhang. BioLiP: A Semi-Manually Curated Database for Biologically Relevant Ligand-Protein Interactions. Nucleic Acids Res. 2013a;41(D1):1096–103.
Yang J, Roy A, and Yang Zhang. Protein-Ligand Binding Site Recognition Using Complementary Binding-Specific Substructure Comparison and Sequence Profile Alignment. Bioinformatics. 2013b;29(20):2588–95.
Yang J, and Yang Zhang. I-TASSER Server: New Development for Protein Structure and Function Predictions. Nucleic Acids Res. 2015;43(W1):W174–81.
Zhang Y, Kolinski A, and Jeffrey Skolnick. TOUCHSTONE II: A New Approach to Ab Initio Protein Structure Prediction. Biophys J. 2003;85(2):1145–64. http://dx.doi.org/10.1016/S0006-3495(03)74551-2.
Zhang Y, and Jeffrey Skolnick. SPICKER: A Clustering Approach to Identify near-Native Protein Folds. J Comput Chem. 2004;25(6):865–71.

No competing interests reported.

SupplementaryFigures.docx
S1 Fig. PSIPRED analysis of GABRA6 gene residues. The black boxes indicate the most deleterious missense nsSNPs. S2 Fig. The first predicted 3D structure model of wild-type (A), the result of PROCHECK server (B), and the result of ERRAT server (C) for the ligand-binding domain. S3 Fig. The first predicted 3D structure model of wild-type (A), the result of PROCHECK server (B), and the result of ERRAT server (C) for the transmembrane domain. S4 Fig. coil secondary structure comparison between WT and mutant structures. S5 Fig. β-sheet secondary structure comparison between WT and mutant structures. S6 Fig. β-bridge secondary structure comparison between WT and mutant structures. S7 Fig. bend secondary structure comparison between WT and mutant structures. S8 Fig. turn secondary structure comparison between WT and mutant structures. S9 Fig. α-helix secondary structure comparison between WT and mutant structures. S10 Fig. 3₁₀-helix secondary structure comparison between WT and mutant structures.
SupplementaryTables.zip
S1 Table. All the missense nsSNPs in the GABRA6 gene. S2 Table. Prediction of deleterious missense nsSNPs by four tools (SIFT, PolyPhen-2, PROVEAN, and CONDEL) in the GABRA6 gene. S3 Table. Detailed output of MutPred2 tool for the three most deleterious nsSNPs in the GABRA6 gene. S4 Table. Prediction of structural effects of the three most deleterious nsSNPs in the GABRA6 gene by project HOPE. S5 Table. Predicted ligand-binding sites by FTSite. S6 Table. Predicted ligand-binding sites by COACH.

Download PDF

Version 1

posted

You are reading this latest preprint version

A comprehensive in silico analysis of the functional and structural consequences of the deleterious missense nonsynonymous SNPs in human GABRA6 gene

Status:

Version 1

Abstract

Figures

1. Introduction

2. Materials and methods

2.1. Data mining

2.2. Predicting deleterious missense nsSNPs

2.3. Predicting disease-related variants

2.4. Predicting the molecular and phenotypic consequences of variants

2.5. Predicting the effect of variants on protein stability

2.6. Evolutionary conservation analysis

2.7. Predicting protein secondary structure

2.8. Predicting structural effects of variants

2.9. 3D structure modeling and visualizing of natives and variants structures

2.10. Predicting ligand binding sites

2.11. Molecular dynamics simulation in water

3. Results

3.1. Missense nsSNPs datasets

3.2. Prediction of most deleterious missense nsSNPs

3.3. Prediction of disease-related variants

3.4. Prediction of the molecular and phenotypic consequences of variants

3.5. Prediction of the effect of variants on protein stability

3.6. Evolutionary conservation analysis

3.7. Prediction of protein secondary structure

3.8. Prediction of structural effects of variants

3.9. 3D structure modeling and visualizing of wild and mutant structures

3.10. Prediction of ligand binding sites

3.11. Molecular dynamics simulation of WT and mutant types

4. Discussion

4.1. Missense nsSNPs prediction

4.2. Protein stability analysis

4.3. Evolutionary conservation analysis and ligand binding sites prediction

4.4. Protein secondary structure analysis

4.5. Project HOPE analysis

4.6. 3D structure modeling

4.7. Substitution effects on structure via molecular dynamics simulations

5. Conclusion

Declarations

Author Contribution

References

Additional Declarations

Supplementary Files

Status:

Version 1