Computational analysis of functional SNPs in Alzheimer’s disease-associated endocytosis genes

Dataset retrieval

The information of the genes of interest was retrieved from the Ensembl genome database (release 87) (http://www.ensembl.org/index.html) in March 2017 for computational analysis. Ensembl comprises the variant data from multiple sources, such as dbSNP, COSMIC and DVGa. The positions of the variants were mapped on the human reference genome GRCh38, which is also equivalent to the UCSC (University of California, Santa Cruz) hg38.

Analysis of UTRs variants

Both 3′and 5′ UTR sequences of the genes of interest were collected from the curated database, UTRdb (http://utrdb.ba.itb.cnr.it/). These sequences with different variants were submitted to the UTRScan, which is a pattern matcher to search for UTR motifs collected in UTRSite (Grillo et al., 2010). UTRSite is a collection of functional sequence patterns located in 3′ and 5′ UTR sequences, such as upstream open reading frame (uORF), internal ribosome entry site (IRES), and polyadenylation signal (PAS). If different functional patterns are found between wild type and mutant UTR sequences, the variant is predicted to be functionally significant.

Analysis of splice site variants

The impact of splice site variants was predicted by using MaxEntScan integrated in Ensembl Variant Effect Predictor (VEP) (http://www.ensembl.org/Tools/VEP), which implements the Maximum Entropy Distribution approach. The input for the VEP requires only variant identifier (variant ID). For each variant, a consensus score with value in the range of −20 to +20 is calculated by MaxEntScan. Then, the score difference (in percentage) between wild type and mutant is calculated manually. The consensus score threshold of 3 and score difference threshold of 30% were applied to predict the splice site variants from both positive and negative controls (Desmet et al., 2009). Table S1 shows the possible consequences of splice site variants predicted by MaxEntScan.

Analysis of sSNPs

PredictSNP2 (http://loschmidt.chemi.muni.cz/predictsnp2/) is a web server that predicts the effect of nucleotide substitution in any region of the genome. The server integrates five prediction tools—CADD, DANN, FATHMM, FunSeq2 and GWAVA (Bendl et al., 2016). The consensus score generated by PredictSNP2 ranges from −1 to +1, in which the variant is considered neutral if the consensus score falls between −1 and 0, and is deleterious if the consensus score falls between 0 and +1.

Analysis of nsSNPs

Five prediction software packages that implement different prediction methodologies, namely SIFT (Ng & Henikoff, 2001), PolyPhen-2 (Adzhubei, Jordan & Sunyaev, 2013), Mutation Assessor (Reva, Antipin & Sander, 2011), I-Mutant2.0 (Capriotti, Fariselli & Casadio, 2005) and SNPs&GO (Capriotti et al., 2013), were used to predict the deleteriousness of all the nsSNPs. The nsSNPs that passed all the default cutoff points were considered to be deleterious nsSNPs. Then, the degree of evolutionary conservation of the deleterious nsSNPs was estimated using the ConSurf server (Ashkenazy et al., 2016), and the mutant protein structures were modelled and refined by SPARKS-X (Yang et al., 2011) and ModRefiner (Xu & Zhang, 2011).

SIFT (Sorting Intolerant From Tolerant, version 5.2.2) from http://sift.jcvi.org/ is a sequence conservation-based package that predicts the functional impact of amino acid substitution based on sequence homology and physiochemical diversity (Ng & Henikoff, 2001). The cutoff point was set at 0.05, so variants with a SIFT score lower than 0.05 were considered deleterious.

PolyPhen-2 (Polymorphism phenotype 2, version 2.2.2) from http://genetics.bwh.harvard.edu/pph2/, is a sequence- and structure-based prediction package. It predicts the impact of nsSNPs based on sequential and structural features collected from multiple databases and tools, such as UniProtKB, BLAST, Pfam and DSSP (Adzhubei, Jordan & Sunyaev, 2013). The cutoff point was set at 0.9, so variants with PolyPhen-2 score (known as PSIC score) higher than 0.9 were categorized as probably damaging.

Mutation Assessor (release 3) from http://mutationassessor.org/r3/ performs the sequence conservation-based prediction by assessing the entropy differences of an alignment position. It uses the sequence homology information of both protein families and subfamilies to assess the level of conservation of each residue (Reva, Antipin & Sander, 2011). The functional impact of the variants was categorized into “neutral”, “low”, “medium” or “high” impact categories. In this study, the cutoff level for the deleterious variant was set at “medium” or “high” impact.

I-Mutant2.0 (http://folding.biofold.org/i-mutant/i-mutant2.0.html) is a structure analysis-based prediction software that analyzes the protein stability change for a single point variant (Capriotti, Fariselli & Casadio, 2005). I-Mutant2.0 works by calculating the change in Gibbs free energy between the wild type and mutant protein structures. Variants with energy change value (DDG) less than −1 or greater than 1 indicate that these variants are destabilizing or stabilizing the protein, respectively.

SNPs&GO (Single Nucleotide Polymorphism Database & Gene Ontology), at http://snps.biofold.org/snps-and-go/snps-and-go.html, combines both sequence and structure information to perform prediction. Information from Gene Ontology terms is also taken into consideration in predicting the functional impact of the variants (Capriotti et al., 2013). The cutoff point was set at 0.5, so variants with probability greater than 0.5 were predicted as disease-related.

ConSurf (http://consurf.tau.ac.il/2016/) is a widely used evolutionary conservation analysis tools based on amino acid positions among the homologous protein sequences (Ashkenazy et al., 2016). In ConSurf, the conservation score of each protein residue indicates the level of conservation, with score “1” indicating the residue is highly variable and score “9” indicating the residue is highly conserved. Residues are also predicted to be buried (b), exposed (e), highly conserved and exposed (functional, f), or highly conserved and buried (structural, s).

SPARKS-X fold recognition server (http://sparks-lab.org/yueyang/server/SPARKS-X/index.php) is a template-based protein structure modelling tool (Yang et al., 2011). The modelled protein structures were refined by using ModRefiner, which is a structure refinement server that improves the structure quality, in terms of side-chain positioning and backbone topology (Xu & Zhang, 2011). The structural impacts of the nsSNPs can be predicted by studying the total free energy changes between the wild type and mutant protein structures, as well as the root mean square deviation (RMSD) of the mutant structures from the native protein structure.

Dataset retrieval from Ensembl

The SNP datasets of PICALM, SYNJ1 and SH3KBP1 genes were retrieved from the Ensembl genome database. It contained a total of 1741 SNPs, including 395 UTRs (22%) variants, 188 (11%) splice site variants, 399 (23%) sSNPs and 759 (43%) nsSNPs. Out of a total of 759 nsSNPs in the genes of interest, 106 were classified as “deleterious” in SIFT and “probably damaging” in PolyPhen-2, according to the results in the Ensembl genome database. Figure S1 shows the total number of variants used in this study.

Prediction analysis of UTRs variants

The information of the retrieved UTR sequences and the total number of variants on the corresponding region are shown in Table S2. There were no 5′ UTR variants in SYNJ1 gene in the Ensembl genome database.

In the wild type UTR sequences of PICALM gene, a total of 11 motifs were identified from the UTRScan. Out of a total of 35 UTR variants in PICALM gene, only one 3′  UTR variant, rs796143140, caused an additional motif, namely K-box at position 212–219 of 3′ UTR sequence. K-box is one of the 3′ UTR motifs that mediates the negative post-transcriptional regulation of RNA (Ganguli et al., 2013). K-box motif was not present in the wild type UTR sequence of PICALM gene and the addition of the motif may affect transcript turnover and translational efficiency.

In the wild type 3′ UTR sequence of SYNJ1 gene, a total of 20 motifs were identified from UTRScan. Three 3′ UTR variants, including rs777138831, rs147590143 and rs111516740, caused an addition of upstream open reading frame (uORF) while another four variants, including rs767649429, rs143739621, rs761852713 and rs79652470, resulted in a deletion of uORF in 3′ UTR. Variant rs780084162 caused an additional motif named UNR-binding site (UNR-bs). UNR is one of the RNA-binding proteins that acts as a regulator in the initiation of internal ribosomal entry site (IRES)-mediated translation and cap-dependent translation (Saltel et al., 2017). Therefore, variant rs780084162 may affect the expression of SYNJ1 gene and the actual biological effect of the additional UNR-bs needs to be further experimentally validated.

For SH3KBP1 gene, there were a total of 21 motifs identified in the wild type 3′ and 5′ UTR sequences. Variant rs192424738 caused an addition of uORF; variant rs755895016 caused a deletion of uORF in the 3′ UTR sequence. Variant rs778448080 resulted in a deletion of GY-box at position 2064–2070 of the 3′ UTR sequence. Like K-box and Brd-box, GY-box is one of the 3′ UTR motifs that mediates the negative post-translational regulation of RNA (Lai, Tam & Rubin, 2005). Therefore, deletion of the GY-box caused by variant rs778448080 may affect the translation efficiency of SH3KBP1 gene.

Table 1 lists the variants that change the number of motif matched in UTRSite, as compared with their wild type UTR sequences (see Table S3 for complete prediction analysis results). Both 3′and 5′UTRs are enriched with cis-acting regulatory elements, and both UTRs are important for the transcriptional and translational regulation. In this study, a total of 12 UTR variants were predicted to cause an addition or a deletion of regulatory elements in the UTR sequences. Further analysis on the impacts of these regulatory elements, including the effect of the miRNA binding to the UTR sequence is outside the scope of this study. However, it should be characterized in the future.

Prediction analysis of splice site variants

The prediction analysis of splice site variants was done using MaxEntScan, which is integrated in VEP. Submission to MaxEntScan requires only the variant IDs. The consensus score and score difference between wild type and mutant sequences were obtained after the submission. The consensus score for each variant was calculated based on different protein-coding transcripts and the same variant may have different consensus scores on different transcripts. This allowed users to study the impact of the splicing variant on different transcripts. Table 2 shows the variants with score difference exceeding the defined threshold. Columns “MaxEntScan Ref” and “MaxEntScan Alt” indicate the wild type and mutant consensus scores, respectively. “Score difference” shows the differences between the wild type and mutant consensus scores (as percentages) to predict the variant effects in causing exon extension, intron retention or splice site creation. Table S4 comprises the complete prediction analysis results for all the splice site variants.

A total of 18 variants were predicted to either disrupt or create splice sites, including 3, 13 and 2 variants from PICALM, SYNJ1 and SH3KBP1 genes, respectively. Among these variants, 17 were predicted to disrupt the splice site, in which splicing was omitted or disrupted, resulting in exon extension or intron retention. Variant rs533064963 in SYNJ1 gene was the only variant predicted to create a new splice site in the intron.

Genetic variations are known to affect alternative splicing. For example, some risk alleles in genes including CLU, PICALM and PTK2B exert their biological impacts through alternative splicing in the pathogenesis of AD (Raj et al., 2018). The 18 splice site variants that were predicted to cause alternative splicing in this study should be further investigated for their structural and functional impacts on the proteins.

Prediction analysis of sSNPs

Out of the 399 sSNPs we studied, eight sSNPs were predicted as deleterious sSNPs in PredictSNP2, including five from PICALM gene and three from SYNJ1 gene. Table 3 shows the eight deleterious sSNPs. Tables S5–S7 describe the complete prediction results of sSNPs in PICALM, SYNJ1 and SH3KBP1 genes, respectively.

Among the 99 sSNPs in PICALM gene, five were predicted to have deleterious effects. Four out of those five, including rs779752380, rs777199260, rs775830584 and rs759002840, are located on the PICALM ANTH domain. Among the three deleterious sSNPs in SYNJ1 gene, only rs768565296 is located on the 5-phosphatase domain of the protein. For SH3KBP1 gene, all the 112 sSNPs were predicted to have neutral effects on the gene.

sSNPs affect protein expression and conformation through the regulation of RNA splicing, mRNA stability, translation dynamics, translation rate and co-translational folding of the protein (Sauna & Kimchi-Sarfaty, 2011). The actual biological impacts of the eight deleterious sSNPs predicted in this study should be validated in the future.

Prediction analysis of nsSNPs

All 759 nsSNPs in the three genes were analyzed by five prediction packages. The Ensembl genome database contains prediction results from SIFT and PolyPhen-2, in which a total of 106 nsSNPs were predicted as “damaging” in SIFT and “probably damaging” in PolyPhen-2. Together with the prediction results from Mutation Assessor, I-Mutant2.0 and SNPs&GO, a total of 18 nsSNPs exceeded the cutoff points of all five packages and were identified as deleterious nsSNPs. The prediction results of the deleterious nsSNPs are summarized in Table 4 while the complete prediction results of nsSNPs using five prediction packages are shown in Table S8. All the deleterious nsSNPs were predicted with high SIFT score and most of them have a PSIC score larger than 0.95 in PolyPhen-2. Prediction results from SIFT and PolyPhen-2 showed that all deleterious nsSNPs were highly conserved in the proteins (Table 4).

There were four deleterious nsSNPs in PICALM gene predicted in this study. All of them are located in amino acids 106 to 179 of the PICALM sequence. Positions 106 to 179 are located on the PICALM protein ANTH domain (23 to 283 aa) and ENTH domain (14 to 145 aa), and these two domains were reported to be functionally redundant (Maldonado-Baez et al., 2008; Manna et al., 2015). The ANTH domain of the protein is important in the binding to phosphatidylinositol 4,5-bisphosphate (PIP2). The binding of PICALM to PIP2 is important for the formation of clathrin-coated pit, which is one of the key functions of PICALM in CME (Ishikawa et al., 2015). These deleterious nsSNPs were predicted to cause conformation change and affect PICALM protein function. Fig. 1 shows the sticks representation of the protein structural changes caused by the deleterious nsSNPs in PICALM gene. In Fig. 1, the variant residues are colored yellow while red dashed lines indicate the hydrogen bonds between the residues. Variants rs780443419 (F109S) and rs765338634 (L179P) resulted in an addition or a deletion of hydrogen bond formation between the mutant and neighboring amino acids. Therefore, the substitution of these protein residues could significantly affect the ANTH domain function as well as the overall PICALM protein structure.

For SYNJ1 gene, 13 nsSNPs were predicted as deleterious. Six of them including rs781675993, rs398122403, rs762909719, rs771755243, rs768897710 and rs779479360 are located on the Sac domain while another four including rs756845805, rs751110096, rs147929290 and rs745418083 are located on the 5-phosphatase domain. Both Sac domain and 5-phosphatase domain mediate the crucial SYNJ1 protein function of phosphate hydrolysis, and they were reported to function in a concerted way (see review in Hsu & Mao (2013)). Besides that, a previous study has shown that both Sac and 5-phosphatase domains are important in the synaptic vesicle endocytosis and SYNJ1 D730A mutant completely eliminated the activity of 5-phosphastase domain (Mani et al., 2007). Although the actual impact of the predicted deleterious variants is yet to be determined, these variants may affect the domains function and influence the role of SYNJ1 in endocytosis as well.

Variant rs770229859 is the only predicted deleterious nsSNP in SH3KBP1 gene. This variant changed residue arginine to glutamine at position 648, which is the C-terminal coiled-coil domain of the protein. The coiled-coil domain of SH3KBP1 was previously reported to be important in mediating the interaction between SH3KBP1 and phosphatidic acid (Zhang et al., 2009). Moreover, deletion of the coiled-coil domain directly affected the binding of SH3KBP1 to c-Cbl ubiquitin ligases, in which the interaction is essential for the initiation of CME regulated by SH3KBP1 (Zhang et al., 2009). Interestingly, another study has found that alanine mutation of four basic amino acids (K645, K646, R648 and R650) on the coiled-coil domain can completely disrupt the interaction of SH3KBP1 with phosphatidic acid and c-Cbl (Zheng, Zhang & Liao, 2014). Since the deleterious nsSNP rs770229859 resulted in amino acid substitution at position 648 of SH3KBP1, this variant is predicted to have a high likelihood of influencing the functions of the coiled-coil domain.

The structural and functional importance of the 18 deleterious nsSNPs were further analyzed using ConSurf analysis tools. Evolutionary conservation analysis determines the level of conservation of each protein residue and predicts the potential structural and functional importance of these deleterious variants to the protein. Figure 2 shows that 17 out of 18 (94%) deleterious nsSNPs were analyzed to be “conserved”, with 12 of them (70%) “highly conserved” (score “9”) through homologous sequence alignment. Only one deleterious nsSNP, rs745418083 (L776S) in SYNJ1 gene, was estimated to be “intermediate” in terms of evolutionary conservation. Besides that, 11 of the 18 (61%) deleterious variants were predicted as structural residues and the rest (39%) were functional residues. Figures S2–S4 shows the conservation scores of full length proteins of PICALM, SYNJ1 and SH3KBP1, respectively.

Besides the prediction analysis of the functional and structural importance of the deleterious nsSNPs and their level of conservation on the proteins, the changes of physical and chemical properties between wild type and mutant amino acids were studied. Table S9 shows the hydropathy, polarity and charge differences between the wild type and mutant amino acids of the deleterious nsSNPs.

During protein synthesis, both covalent and non-covalent interactions are driven by the physical and chemical properties of the amino acids. The hydropathy, size and charge of the amino acids all play important roles in the protein synthesis (Biro, 2006). Table S9 shows that the hydropathy in eight deleterious nsSNPs has changed from hydrophobic to hydrophilic, and the polarity of four nsSNPs has changed from non-polar to polar. The substitution of amino acids may affect both covalent and non-covalent interactions among amino acids, subsequently influencing the stability and conformation of protein structure.

To study the role of nsSNPs in affecting the total free energy and the stability of protein structures, full lengths of both wild type and mutant protein structures were modelled with the use of SPARKS-X server (Yang et al., 2011), and the structures were refined by using ModRefiner (Xu & Zhang, 2011). Energy minimization of the protein structures was performed by GROMOS96 force field using Swiss-PdbViewer. Table 5 shows the energy deviation and RMSD between wild type protein structures of PICALM, SYNJ1 and SH3KBP1, and the 18 mutant protein structures of deleterious nsSNPs. RMSD values were calculated to provide an insight into the deviation between the wild type and mutant protein structures. The RMSD values are proportional to the deviation between the wild type and mutant structures. The total free energy of wild type PICALM structure was −24,238.746 kJ/mol, and the total free energy of the four PICALM mutant structures was marginally decreased. The RMSD values also showed a minor deviation ranging from 1.13 to 1.45 Å. For SYNJ1 protein, the total free energy of the wild type structure was -19307.051 kJ/mol. The deviation of total free energy between SYNJ1 wild type and mutant protein structures is significant. Among the 13 SYNJ1 mutant structures, variant rs398122403 (R258Q) has the highest total free energy of 3,057.503 kJ/mol, and variant rs775515863 (G487R) has the lowest total free energy of −28,591.270 kJ/mol. Among all the SYNJ1 variants, rs781675993 (N200K) has the highest RMSD value of 2.34 Å. Figure 3 shows the structural comparison of the wild type SYNJ1 structure (green) and rs781675993 mutant structure (blue). Superposition of the two structures shows the structural deviation caused by the variant. The only SH3KBP1 mutant structure caused by variant rs770229859 (R648Q) showed a decreased total free energy of −1,6167.838 kJ/mol, with a RMSD value of 1.88 Å.

Verification of nsSNP prediction principle

To demonstrate the reliability of nsSNP prediction, we predicted the functional consequences of other nsSNPs that have been previously studied in other benchwork experiments. Ten PSEN1 pathogenic nsSNPs that have been validated experimentally to affect amyloid-beta (Aβ) level were retrieved from ALZFORUM (https://www.alzforum.org/). Besides that, another five PSEN1 nsSNPs, including three non-pathogenic variants and two variants that have never been reported to be deleterious or disease-associated, were selected as negative controls. The prediction results of the total 15 nsSNPs in PSEN1 gene are shown in Table S10. The prediction results show that nine out of these ten pathogenic nsSNPs were predicted deleterious by all five prediction packages used in this paper. The only pathogenic nsSNP that was not predicted as deleterious variant, which is rs63750231 (E280A), has I-Mutant DDG and SNPs&GO score that is lower than the cutoff point. All these five negative controls of PSEN1 nsSNPs were predicted to be non-deleterious to the protein.

Computational prediction of the biological consequences of variants was previously reported to yield a high false positive (FP) and low false negative (FN) rate, suggesting that computational prediction tends to overestimate the deleterious effects of variants (Miosge et al., 2015). To reduce the FP rate, our study implemented meta-prediction method to perform functional prediction of nsSNPs. Despite one FN prediction result (rs63750231) was reported, the overall prediction results of the PSEN1 pathogenic and non-pathogenic nsSNPs indicate that the prediction approaches which were used, reliably predicted the functional consequences of nsSNPs.

NGS data of AD studies

To check if any of the 56 deleterious variants were identified in the next generation sequencing (NGS) studies of AD, genotype calling data (variant calling format (.VCF) files) were retrieved from three studies including Alzheimer’s Disease Sequencing Project (ADSP), Alzheimer’s Disease Neuroimaging Initiative (ADNI) and Mount Sinai Brain Bank (MSBB) studies. ADSP aims to identify genetic variants that increase the risk for or protect against AD. The study consists of two parts, including a whole genome sequencing (WGS) of 584 samples family-based study and a whole exome sequencing (WES) of 10,061 samples case-control study (Beecham et al., 2017). Genotype calling data of both WGS and WES were retrieved from the study. For ADNI, the study is primarily designed to discover novel biomarkers from a total of 800 individuals for the early detection of AD. The 800 individuals consist of 200 elderly controls, 400 mild cognitive impairment (MCI) and 200 AD patients. WGS of the 800 individuals was performed in the following add-on study and the genotype calling data is publicly available (Weiner et al., 2015). For MSBB study, it aims to discover novel molecular mechanism of AD pathogenesis from genomic, transcriptomic and proteomic data consisting of 364 AD case-control human brains (Wang et al., 2018). The raw and processed data of MSBB study are available through the Synapse software platform (https://www.synapse.org/#!Synapse:syn3159438) and the .VCF file of WES of the study was retrieved.

Among the 56 deleterious variants predicted in this study, seven variants were also identified in ADSP, ADNI or MSBB studies. Table S11 shows the deleterious variants that were identified in the NGS studies of AD. Identification of these deleterious variants in these NGS studies further suggests their potential association with AD. In these studies, the genotype data was mapped on the reference genome GRCh37 while the variant data retrieved from Ensembl genome database was mapped on GRCh38. Therefore, the variant positions in both GRCh37 and GRCh38 are provided in Table S11.