We retrieved the relevant publications in the PubMed database up to March 31, 2021, to obtain the complete and detailed information of genes and variants related to HL. The keywords used for searching were “(hearing loss OR hearing impairment) AND (mutation OR variant) (Title/Abstract).” For SHL, common diseases with HL, including Alport Syndrome (Kandai et al., 2019; Nozu et al., 2019), Branchio-Oto-Renal Syndrome (Kochhar et al., 2007), CHARGE Syndrome (Hsu et al., 2014), Jervell & Lange-Nielsen Syndrome (Cusimano et al., 1991), Norrie Disease (Sowden et al., 2020), Pendred Syndrome (Reardon et al., 1996), Perrault Syndrome (Pan et al., 2020), Stickler Syndrome (David 2002), Treacher Collins Syndrome (Dixon 1996), Usher Syndrome (Mathur et al., 2019), Waardenburg Syndrome (Pingault et al., 2010) were searched according to Hereditary Hearing Loss Homepage (Van et al., 2021). Next, all HL-related publications were reviewed manually. The annotation of genes and variants, detailed information of clinical data of patients (ethnicity, hearing, and other symptoms) were extracted through in-depth reading of the full text of each publication (Figure 1).

Gene4HL integrated the information of 62 public databases (Supplementary Table S1) and annotated six aspects from the gene level. The meaningful annotations for each gene were collected from public databases, including basic information, phenotype and disease, protein sequences, protein-protein interaction, gene expression pattern, and gene functions, drug-gene interactions, and the gene-drug abilities for precision medicine. Each extracted variation was annotated in detail using ANNOVAR (Wang et al., 2010), such as the functional effects and functional prediction of variants as well as allele frequencies of variants in different populations based on the definitions of transcripts from RefSeq. The annotations of each variation included the basic information (location, PubMed ID, variant type, genetic pattern, the allele frequencies of different populations, etc.), predictive scores and pathogenicity, variant and related diseases or phenotype information, the total genetic testing samples, the positive samples, etc.

Based on the functional effects and ReVe (Rare Exome and Variant Effect Scoring Tool) score, which were previously proposed by our group (Li et al., 2018), variants were divided into five classes: 1) loss-of-function (LoF) variants, including stop-gain, stop-loss and splicing variants, 2) damaging missense variants (ReVe score >0.7), 3) tolerant missense variants (ReVe score ≤0.7), 4) synonymous variants, 5) non-frameshift indels variants, including non-frameshift deletion, non-frameshift insertion, non-frameshift substitution, and 6) noncoding variants.

Next, we extended our previously developed scoring system in other genetic diseases (Li et al., 2021; Zhao et al., 2021) to prioritize genes associated with HL to quantify the contribution of the different classes of variants and integrate them into the entire prioritization of deafness genes (Supplementary Table S2). Briefly, the LoF variants were assigned an evidence score of five; the damaging missense variants and tolerant missense variants were assigned an evidence score of 3 and 2, respectively; the non-frameshift indels and noncoding variants were assigned an evidence score of 1. An integrated summing up evidence score for each gene was calculated by adding the evidence scores of each type of variant from all integrated studies. All genes integrated in this study were classified into three classes: high confidence (score ≥20), strongly associated (score of 10–20), and suggestively associated (score of 5–10). We performed a permutation test to evaluate the interconnectivity and functional correlation among the HL-related genes (score ≥5). Specifically, we randomly simulated 1,000,000 permutation tests to evaluate the interconnectivity among the HL-related genes.

The spatiotemporal transcriptase of the mouse inner-ear organs (cochlea and utricle) at four developmental stages (E16, P0, P4, and P7) were sourced from a previous study (Shen et al., 2015) (https://shield.hms.harvard.edu/). We applied signed hybrid-weighted gene co-expression network analysis (WGCNA) (Langfelder et al., 2008) to analyze all 16 samples using the standard method with a power of eight to cluster the spatiotemporal-expression patterns and prenatal laminar-expression profiles of a given gene set.

A permutation test was performed to investigate the interconnectivity and functional correlation of HL-related genes (evidence score ≥5) using the STRING v11.0 database (Szklarczyk et al., 2019). In addition, we constructed a protein-protein interaction network using the STRING online analysis platform (https://string-db.org/) with a confidence score >0.4. Moreover, the functional networks were clustered by multiple biological processes of Gene Ontology (GO) (The Gene Ontology 2017) (http://www.geneontology.org/).

Gene4HL was developed, supported by versatile browsing and searching functionalities, similar to the Gene4Denovo database (Zhao et al., 2020). All data were stored in a MySQL database. Users could access the genetic data or extended analysis results freely through this web interface. The web interface of Gene4HL contained Search, Analysis, Browse, Upload, and Download modules.

A total of 1,608 publications met the inclusion criteria and were used for extracting further information. Gene4HL database (http://www.genemed.tech/gene4hl/) integrated 326 HL-related genes, including 170 NSHL-related genes and 156 SHL-related genes, involving 3,872 genetic variations. The corresponding detailed genetic and clinical information was also integrated into the Gene4HL database. The database was comprehensively annotated at the variant-level and gene-level, including 1) functional effects of variants (nonsense, nonsynonymous, and frameshift, etc.), 2) disease and phenotype-related information for variant- and gene-level implications, 3) functional consequences of variants through 24 in silico predictive algorithms, 4) allele frequency in different populations of public databases, 5) meaningful gene-level information, such as protein sequences, protein-protein interactions, the gene expression patterns in human tissues, and gene functions, etc., and 6) drug-gene interactions and the gene-drug abilities for precision medicine (Figure 1).

Notably, we found that 7 of the 10 most frequently occurring variants of Gene4HL were located on GJB2, the most common HL gene associated with DFNB1A and DFNA3A (Supplementary Figure S1A). Next, we classified variants based on the DVD and the American College of Medical Genetics and Genomics (ACMG) guidelines for the interpretation of variants in the context of Gene4HL. The variants were classified as pathogenic (P), likely pathogenic (LP), variant of uncertain significance (VUS), likely benign (LB), or benign (B). Of these variants, 2,115 (54.7%) were classified as P, 748 (19.3%) were LP, 127 (3.3%) were LB, 474 (12.2%) were B, and 406 (10.5%) were VUS (Supplementary Figure S1B). Of the 3,873 variants listed in the Gene4HL, 2,343 were found in the DVD.

The query interface contains panels for quick and advanced searching. The quick searching function is the main tool for quick access to detailed information regarding genes or variations, which can be found on the homepage. The quick search automatically identifies seven key terms, such as gene symbol, genomic region, cytoband, transcript accession, the nucleic acid change in a certain gene or transcript, the genomic coordinate of a variant, and the Gene4HL ID (Figure 2).

Advanced searching supports batch search and allows users to specify annotated datasets (Figure 2). The advanced searching provides options for primary information, prediction algorithms for nonsynonymous variations, allele frequencies in different populations, and disease-related and phenotype-related information. The advanced searching also has the same seven input forms like quick searching. The advanced searching query form and the corresponding result set have been displayed on the same page to improve the users’ experience. Also, users can specify datasets and obtain various information.

Both quick and advanced searching results provide detailed annotation data of genes and variations. The searching results are contained in two tables, a summary of the variation of each gene and detailed annotation at the variant level. The summary table presents the number of loss of function, damaging missense, tolerance missense, synonymous, non-frameshift, and noncoding variants. The detailed table contains the specific information of each variant: the Gene4HL ID, gene symbol, location, functional variant, amino acids change, the pathogenicity, inheritance pattern, and a number of reports. In this searching result interface, users can access annotation on gene-level by clicking the gene symbol in the table. Gene4HL integrates the information of the main public databases and annotates six aspects from the gene level: 1) Basic information, including the official name, location, gene type, gene summary (UniProt 2018), the genic intolerance score (Forge et al., 2003; Itan et al., 2015; Petrovski et al., 2015; Aggarwala et al., 2016; Lek et al., 2016; Fadista et al., 2017; Han et al., 2018); 2) Gene function, including the molecular function, gene ontology terms (The Gene Ontology 2017), domain information (Finn et al., 2017), protein–protein interactions (Li et al., 2017), and biological pathway information (Geer et al., 2010); 3) Phenotype and disease, including phenotype data (Amberger et al., 2015), clinical variation data (Landrum et al., 2016), mammalian phenotype (Eppig et al., 2017), and human phenotype ontology (Kohler et al., 2017); 4) Gene expression, including spatio-temporal expression profiles (Miller et al., 2014), gene expression data in tissues (GTEx Consortium, et al., 2017), and the subcellular location (Uhlen et al., 2015); 5) Variants in different populations; 6) Drug–gene interactions (Wagner et al., 2016).

Gene4HL provides users with the interface to freely analyze their own genetic data. In this section, users can upload data of genetic variants (VCF4 format) to identify the co-segregated variants associated with HL and Gene4HL, prioritizing the HL-risk genes in this panel. This upload process involves four simple steps: 1) enter your e-mail address, 2) choose options Trio or Non-trio, 3) upload your file, and 4) fill up trio information or genotype information. Then, datasets of your interest need to be specified and submitted. The main annotation data sets that can be selected are: primary information (gene system, splicing threshold, whether to send all co-segregated variants, reference and Gene4HL variant-level annotation), identify co-segregated variants, perform comprehensive annotations (basic information annotation, pathogenicity prediction of missense variants, allele frequency in various populations, related clinical databases), rare damaging variants (in silico predictive algorithms for nonsynonymous variants, allele frequency in different populations). Concurrently, users can specify their cut-off values of quality control, the annotation data sources, and the methods to detect rare damaging variants. The results are sent to the users’ e-mail once the annotation completes.

Gene4HL also contains other five useful sections: 1) Browse, to facilitate users to browse the HL-related genes in Gene4HL; 2) Upload, to encourage the users to upload their variants list to Gene4HL database; 3) Download, to allow users to freely access all released datasets in Gene4HL and export and download as Excel files; 4) Data Source included in Gene4HL; 5) Tutorial, which provides a further description of Gene4HL.

We developed a weighted scoring system to prioritize HL-related genes by combining all the above genetic evidence. We prioritized 71 high confidence genes (score ≥20), 42 strong associated genes (20 > score ≥10), 52 suggestive associated genes (10 > score ≥5), total of 164 genes (Table 1).

Furthermore, we performed a permutation test based on random resampling, which was used to assess whether functional interconnections of 164 HL-related genes were more than random expectation (evidence score ≥5, Table 1). The permutation test took gene length into consideration. We observed 135 of 164 HL-related genes (p < 10^–6) that interacted with each other and had 809 interconnections (p < 10^–6), which were significantly higher than the random expectation (Supplementary Figure S2).

Next, we performed WGCNA in the tissue samples from FACS and identified three independent modules (M1-M3) comprising 109 genes to characterize the spatiotemporal-expression patterns of the 164 HL-related genes in convergent networks during inner-ear development (Figure 3 and Supplementary Table S3). Compared with the surrounding cells, we found that majority of genes within M1 (n = 44) were expressed lower in the hair cells of cochlea and utricle during E16, P0, P4, and P7 periods. In contrast, the expression of M2 genes (n = 35) was elevated in the hair cells of the cochlea and utricle. The pattern of gene expression in the surrounding cells of the cochlea and utricle within the M3 module was similar to the M2 module (n = 30). However, the gene expression in the hair cells of the cochlea and utricle within M2 was lower in the hair cells of utricle during the P0, P4, and P7 periods.

Furthermore, we performed PPI base on the STRING v11.0 database and identified three independent networks (N1 red, N2 green, and N3 blue) to characterize the protein-protein interaction (PPI) patterns of the 164 genes in convergent networks. Specifically, the interacted functional network had 163 of the 164 HL-related genes with 930 connections of each other at protein-level (PPI enrichment P-value: < 1.0e-16) (Figure 4A and Supplementary Table S4).

GO enrichment analysis of the 163 genes identified several pathways associated with HL (Figure 4B and Supplementary Table S5). Most of these GO terms, such as sensory perception of sound (GO:0007605, FDR = 1.62E-82), sensory perception of mechanical stimulus (GO:0050954, FDR = 8.43E-80), sensory perception (GO:0007600, FDR = 2.53E-49), nervous system process (GO:0050877, FDR = 1.64E-45), inner ear development (GO:0048839, FDR = 8.70E-43), ear development (GO:0043583, FDR = 2.71E-42), system process (GO:0003008, FDR = 3.83E-41), sensory organ development (GO:0007423, FDR = 4.58E-33), ear morphogenesis (GO:0042471, FDR = 1.61E-27), inner ear morphogenesis (GO:0042472, FDR = 2.06E-26), multicellular organismal process (GO:0032501, FDR = 6.36E-25), inner ear receptor cell differentiation (GO:0060113, FDR = 8.91E-22), were regarded as critical functional signaling pathways associated with HL. The functional network suggested that the prioritized HL-related genes shared a common signaling mechanism and were functionally correlated.

In addition, we inspected the contribution of genes in three independent networks N1, N2, and N3, to the functional signaling pathways analyzed by GO enrichment analysis (Figure 4B). The contribution of the number of genes in N1 accounted for more than 50% for the following pathways: sensory perception of sound (GO:0007605), sensory perception of mechanical stimulus (GO:0050954), sensory perception (GO:0007600), inner ear development (GO:0048839), ear development (GO:0043583), system process (GO:0003008), ear morphogenesis (GO:0042471), inner ear morphogenesis (GO:0042472), inner ear receptor cell differentiation (GO:0060113).