An in-depth understanding of structure is the basis of studying mutations in a structural context, so we investigated the structure of human RPS19 in detail. Its structures both in free state and in the SSU complex were inspected in order to obtain an accurate description of its secondary structure, hydrophobic core, and interaction with other molecules.

We have curated 51 DBA missense mutations (Supplementary Table 2) and 30 neutral ones (Supplementary Table 3). These data presumably constitute a very complete missense mutation dataset of RPS19 until now. We first checked the conservation features of DBA mutations. According to the calculated Consurf scores (Ashkenazy et al., 2016), the conservation levels of DBA mutation sites are significantly higher than the neutral ones (median: -0.665 vs. 0.0995, p = 8.9E-5, one-tailed Mann-Whitney U test, Figure 3A). Nearly one half of DBA mutations are located at the highly conserved sites (Consurf score < -1.0). These results are consistent with the general knowledge that disease mutations tend to occur at conservative sites (Kumar et al., 2009; Lal et al., 2020), which are indicative of functional or structural importance. Although it shows the potential in discriminating DBA and neutral mutations, the conservation feature cannot provide specific clues concerning the pathogenesis. We then analyzed the residue types and counts of the DBA mutations.

Overall, DBA mutations are more likely to occur at positively charged, hydrophobic or small residues (Figure 3B and Supplementary Table 2). Arginine (Arg) is the most mutated residue, accounting for about 19.6% (10/51) of all DBA mutations. Leucine (Leu) is the second most mutated residue (17.6%, 9/51), and small volume residues such as Glycine (Gly) and Alanine (Ala) are also mutated frequently (Figure 3B). We speculate that these Arg sites may participate in interactions with negatively charged rRNAs, and their mutations may disrupt or weaken these interactions. As for Leu, the speculation is that they often form the hydrophobic core, and mutations to polar residues may disturb the stability of the protein. Mutations from small residues to large ones may also influence the stability due to steric hindrance effect.

One striking feature is that a larger number of residues (29.4%, 15/51) have been mutated into Proline (Pro) (Figure 3B), a well-known helix breaker due to its cyclic side chain (a pyrrolidine ring) (Woolfson and Williams, 1990). Considering that the hydrophobic core of RPS19 is mainly composed of the bundle of the five helices, these mutations will presumably affect protein structure stability by destroying helices. There were also many other residues mutated to Arg (Figure 3B), most of which were mutated from hydrophobic residues, such as Leu and Tryptophan (Trp) (Supplementary Table 2). The hydrophobic effect of these residues would be disrupted by charged Arg, and thus protein folding would be affected.

Based on the analyses concerning the residue types, we can roughly infer that DBA mutations may mainly disrupt RPS19 from two perspectives: protein folding stability and interaction, as other studies have done (Zhang et al., 2010; Peng et al., 2018) or reviews have summarized (Stefl et al., 2013; Yates and Sternberg, 2013). However, more systematic data mining is required to consolidate these inferences, and more detailed mechanisms need to be explored in order to answer how a DBA missense mutation destabilize RPS19 or perturb its interactions. Compared to sequence, structure is more directly related to function. In the following analyses, more in-depth investigations in structural contexts have thus been conducted.

Our efforts of understanding the structural basis of DBA mutations of RPS19 can provide clues to predict the pathogenicity of newly identified RPS19 mutations. General prediction tools, such as PMut (Lopez-Ferrando et al., 2017), MutPred2 (Pejaver et al., 2017), PolyPhen2 (Adzhubei et al., 2010), and SIFT (Sim et al., 2012), do not perform well with high false positive rate (FPR) on RPS19 mutations (Figures 8A,B), possibly because they internally have the tendency to overestimate the pathogenicity of mutations, as discussed previously (Andersen et al., 2017; Holland et al., 2017; Peng et al., 2018). Another reason may come from that these general tools do not incorporate specific features about RPS19. If we adopt the RPS19 mutations as a specific dataset and incorporate as features the clues in the understanding of DBA mutations, it should be promising to build a better pathogenicity predictor by using machine learning methods, as demonstrated in developing the IDR-specific disease mutation predictor previously (Zhou et al., 2020).

We adopted the support vector machine (SVM) to build the RPS19-specific prediction tool. After further curation, 29 DBA mutations (positive samples) and 30 neutral ones (negative samples) were selected from the mutation dataset as training data. We then extracted the features based on our understandings of the structural basis of DBA mutations. In total, 18 candidate features were extracted for each mutation concerning the interaction with rRNA, structural stability, conservation, etc. After feature selection, 8 features were finally selected (BSA, rBSA, HB_Num, Δcharge, ΔΔG, Δhelix, BLOSUM62, and Δdisorder) (Supplementary Table 7).

Based on fivefold cross-validation, we identified the best hyper-parameters (C = 100, γ = 0.01), and 26 of the 29 DBA mutations can be correctly predicted in the cross-validation, similar to other well-known tools (Figure 8A). At the same time, the FPR has been significantly decreased compared with others (Figure 8A). Overall, much better performance has been achieved as measured by ACC, F1 score, and MCC (Figure 8B). We manually inspected the nine false positives of SIFT (Pro2Leu, Val4Phe, Thr5Pro, Lys24Asn, Asp35Gly, Thr60Ala, Tyr79Cys, Gln105Arg, His145Tyr), and found most of them are located at conserved sites according to the Consurf scores (Supplementary Table 3). According to its prediction logic, SIFT will tend to predict these variants as deleterious. As for the cross-validation of SVM training here, several features independent of conservation have been adopted, which may be responsible for its lower FPR. For these variants, most of them have small or even negative ΔΔG values (Supplementary Table 3), indicating only trivial perturbation of structural stability, which may serve as the reason why the SVM classifier predicts them as neutral.

Finally, we re-trained the predictor, namely RPS19-SVM, on all the training data with the best hyper-parameters, and utilized it to predict the pathogenicity for all possible missense mutations of RPS19 (Figure 8C and Supplementary Table 8). The resulted pathogenicity scores can be valuable in inferring the disease-association of newly identified RPS19 mutations, or the deleteriousness of newly designed mutagenesis mutations.

As ExAC has been upgraded to gnomAD later (Karczewski et al., 2020), we retrieved 25 additional neutral missense variants of RPS19 from gnomAD v2. They were not utilized in the training of RPS19-SVM, so we can use them as an independent testing dataset. It turns out that RPS19-SVM can accurately predict 23 of them as “neutral” with a low FPR of 8% (2/25), much better than the other four tools (PMut: 60%, MutPred2: 76%, PolyPhen2: 20%, SIFT: 40%), confirming its superior performance (Supplementary Table 9). When new positive data are available in the future, further validation could be conducted as well.

Moreover, the heatmap based on the pathogenicity scores (Figure 8C) also confirms our understanding of their molecular mechanisms and can provide new insights. First, mutations substituted by positive residues (Arg, Lys, and His) or by negative residues (Asp and Glu) have relatively low or high pathogenicity scores, respectively. Second, mutations substituted by Proline are not tolerated by RPS19, especially in helix regions. These results are consistent with its characteristics of interaction with negatively charged 18S rRNA and hydrophobic core with a bundle of helices. Moreover, we also found that the two β-hairpin regions, overlapping with the second and third putative IDRs, had lower mutation tolerance for substitution of hydrophobic residues, though few disease mutations in these regions were reported previously. They participate in the composition of surface Patch III and IV, and substitutions by hydrophobic residues may thus affect the interactions with 18S rRNA.

The RPS19 structure was extracted from the coordinates of human ribosome SSU complex (PDB ID: 6G5H) (Ameismeier et al., 2018), and was used as initial conformation for performing molecular dynamics (MD) simulation. The simulation was conducted by using the Amber16 package (Case et al., 2016) with the force field of RSFF2C (Kang et al., 2018), and the explicit solvent model of TIP3P (Jorgensen et al., 1983) was adopted with a periodic box whose edges had a minimum distance of 8 Å to any atom originally presented in the solute.

Before production simulation, energy minimization was performed to relax possible atom collisions, and the system was equilibrated for 2 ns with the final temperature reaching 300.0 K. The production simulation was conducted in an NPT ensemble and lasted 1,000 ns. A timestep of 2 fs was used in both equilibration and production simulation.

The conformations were saved in the trajectory file with an interval of 1 ps, and were analyzed using CPPTRAJ (Roe and Cheatham, 2013). Conformation clustering (dbscan) was performed by setting the backbone RMSD cut-off to 2 Å. The representative conformation of the resulted largest cluster was used for the assignment of secondary structure by using DSSP (Kabsch and Sander, 1983). The 8 states of secondary structure definition were then simplified into 3 states: helix (G, H, I), sheet (B, E), and coil (T, S, C). The backbone RMSF (root mean square fluctuation) of RPS19 was calculated by sampling a frame per 1 ps and using the average conformation as reference. The residues in the middle of the protein with RMSF greater than 2 Å were considered as intrinsically disordered.

The interaction interface of RPS19 in the ribosome complex was identified by using PDBePISA (version 1.48¹; Krissinel and Henrick, 2007).

The DBA missense mutations were collected from dbagenes (Boria et al., 2010), ClinVar (Landrum et al., 2016), and manual literature review. The canonical transcript (ENST00000598742 in Ensembl or NM_001022.3 in RefSeq database) (O’Leary et al., 2016; Yates et al., 2020), was adopted as the reference to map all the mutations. During the data review, several mutations, incorrectly recorded in their original sources, were corrected in this work, including c.182C > A (Ala61Glu), c.358G > C (Gly120Arg), and c.380G > A (Gly127Arg).

The missense variants of RPS19 in the ExAC (Karczewski et al., 2017) were also collected. Based on that the prevalence of DBA is 5–7 in one million individuals, and 25% of them are caused by RPS19 mutations, and the inheritance pattern is autosomal dominant, the minor allele frequency (MAF) of pathogenic RPS19 mutations can be estimated as 6.25E-7 to 8.75E-7. The collected ExAC variants here have allele frequency values with at least one order of magnitude higher. Moreover, considering that DBA usually presents within the first year of life (Ulirsch et al., 2019) and that ExAC has attempted to exclude severe pediatric diseases according to its description (Lek et al., 2016), one can assume that most of these ExAC variants would be neutral. After removing those that has been explicitly recorded as DBA mutations (c.68A > G, c.164C > T, c.208G > A, c.301C > T), the remaining ones served as the neutral dataset in this work.

The Consurf score of each site of RPS19 was calculated by using Consurf 2016 server (Ashkenazy et al., 2016). The scores were normalized with mean of 0 and standard deviation of 1. The lower the score, the slower the evolution rate, and the higher the conservation level.

The electrostatic potential on the surface of RPS19 was calculated by using the APBS Electrostatics Plugin integrated in PyMOL (version 2.3.0) with default parameters (Baker et al., 2001; Dolinsky et al., 2004; Schrodinger, 2015).

The Consurf scores and electrostatic potential were rendered at the protein’s solvent-excluded molecular surface with different color schemes by using UCSF Chimera (Pettersen et al., 2004).

The effect of a mutation on protein structural stability was measured by the change of folding free energy (ΔΔG), and was calculated by using the FoldX package (Schymkowitz et al., 2005). In detail, the RPS19 coordinates from Cryo-EM ribosome structure (PDB ID: 6G5H) was first processed iteratively by using RepairPDB command of FoldX (an energy minimization process) until the decrease of calculated folding energy was lower than 1 kcal/mol in at least 5 iterations. Then we used the PositionScan command of FoldX to mutate all sites to all possible residues, and to calculate the ΔΔG values for all the mutations.

The solvent accessible surface area (SASA) of each residue was calculated by DSSP (Kabsch and Sander, 1983). The relative SASA (rSASA) was obtained by dividing SASA with the maximum SASA, which was computed by placing the specific residue between two Gly residues in an extended conformation accordingly (Tien et al., 2013).

The neutral mutations collected from ExAC were treated as negative samples. The positive samples were selected from the DBA mutations if they met any one of these criteria: (1) identified in more than one patient, (2) explicitly annotated as “Pathogenic” or “Likely pathogenic” in ClinVar, and (3) experimentally confirmed to affect the physiological function of RPS19 (expression, nucleolar localization, ribosomal abundance, etc.).

Based on our analyses of DBA mutations, 18 features in four categories were extracted for describing each mutation. The details of these features are listed in Supplementary Table 7.

The feature selection and the SVM hyper-parameter searching were conducted by using the scikit-learn (version: 0.21.3) package in Python. In detail, all the possible feature combinations with at least 5 features were enumerated, and the hyper-parameters C and γ were enumerated from 0.001, 0.01, 0.1, 1, and 100. For each combination of features and hyper-parameters, fivefold cross-validation was run on the training dataset containing the positive and negative samples. When the maximum cross-validation MCC was reached, the optimal feature combination and optimal hyper-parameters were obtained accordingly. The resulted optimal hyper-parameters and features were then utilized to train the RPS19-SVM predictor on all the training data. The predicted results of PMut (Lopez-Ferrando et al., 2017), MutPred2 (Pejaver et al., 2017), PolyPhen2 (Adzhubei et al., 2010), and SIFT (Sim et al., 2012) were obtained by submitting the mutations to their web servers with default settings.