ExIR enables prioritizing driver and biomarker genes from omics data in a reference free manner

Article type: Research Article

Keywords: health sciences, medicine, bioinformatics

Authors: Adrian Salavaty, Alon M. Douek, Jan Kaslin, Mirana Ramialison, Peter D. Currie

Affiliations: Australian Regenerative Medicine Institute, Monash University, Clayton, VIC 3800, Australia; Sir Peter MacCallum Department of Oncology, The University of Melbourne, Parkville, VIC 3010, Australia; Division of Cancer Research, Peter MacCallum Cancer Centre, Melbourne, VIC, Australia; Systems Biology Institute Australia, Monash University, Clayton, VIC 3800, Australia; Institute of Biomedicine, Faculty of Medicine, University of Turku, Turku, Finland; Novo Nordisk Foundation Centre for Stem Cell Medicine, Murdoch Children’s Research Institute, Royal Children’s Hospital, Flemington Road, Parkville, VIC 3052, Australia; Department of Paediatrics, The Royal Children’s Hospital, University of Melbourne, Parkville, VIC 3052, Australia; EMBL Australia, Monash University, Clayton, VIC 3800, Australia

License: © 2026 The Authors CC BY 4.0 This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

Article links: DOI: 10.1016/j.isci.2026.116303  |  PMC: PMC13264208

Relevance: Moderate: mentioned 3+ times in text

Full text: PDF (14.4 MB)

Introduction

The emergence of advanced high-throughput technologies, including bulk and single-cell RNA sequencing as well as sophisticated mass spectrometry-based proteomic and metabolomic methodologies, have greatly enhanced the capacity for generating hypotheses relating to complex biological processes.ref. bib1^,ref. bib2^,ref. bib3^,ref. bib4 Nevertheless, the subsequent experimental validation of these hypotheses invariably necessitates lower-throughput focused analytical approaches, including extensive manual curation and expert review of the generated data in the context of existing literature.ref. bib2^,ref. bib5 Such approaches are heavily time-consuming, and the identification and selection for further investigation of the most relevant candidates from high-throughput datasets frequently represent the most significant bottleneck in progressing exploratory omics-based studies to verified biological findings.ref. bib2^,ref. bib5

To address this bottleneck, a variety of computational models have been proposed to prioritize candidate features from lists of genes or proteins, such as those extracted from differential expression or abundance analyses. These prioritization algorithms employ a diverse range of strategies, including network-based algorithms and machine learning (ML) models, to filter and rank candidate features (e.g., genes and proteins). For example, graph-based frameworks have been developed to integrate multimodal data for enhanced prediction accuracy in related biological contexts such as gene regulatory networks,ref. bib6^,ref. bib7 while diffusion models enable robust inference of cell-specific networks from noisy single-cell data.ref. bib8 Similarly, feature selection techniques using recursive methods with random forests have improved classification of structural classes in low-similarity sequences, demonstrating the value of intrinsic data-driven prioritization.ref. bib9 However, most of these models heavily rely on external—often manually curated and inherently incomplete—sources of evidence, such as gene ontology and pathway annotations, and protein-protein interactions (PPIs).ref. bib10^,ref. bib11^,ref. bib12 Knowledge-based prioritization models manifest three significant shortcomings.

First, reliance on prior knowledge as an immutable input source often fails to account for the molecular heterogeneity inherent to a given biological sample.ref. bib13 Second, the dependency on external information sources, particularly those assembled by text-mining of literature (e.g., co-mentioning of gene names) such as gene ontology (GO) and STRING databases,ref. bib14 introduces bias favoring extensively researched genes,ref. bib15 consequently elevating the risk of introducing both type I and II errors (i.e., false positive and false negative) during feature prioritization. Third, these models cannot be robustly applied in contexts where prior knowledge or supporting evidence sources is lacking; for example, in novel biological processes, or in molecular studies using non-model organisms. By contrast, a robust gene prioritization model solely reliant on the experimental data under investigation could circumvent these deficiencies, thereby augmenting accuracy of candidate feature prioritization. Approaches that detect sequence-level signatures from single genomes without annotated priors exemplify this potential, using statistical divergence measures and multiscale testing to identify atypical regions with high fidelity.ref. bib16^,ref. bib17^,ref. bib18 Complementary methods that incorporate additional modalities like epigenetic signals, such as those that perform enhancer prioritization via ML on chromatin interactions, further highlight how data-intrinsic models can link regulatory elements to disease risk without external biases.ref. bib19

In addition to the ranking of features to streamline experimental validation of high-throughput data, the ability to apply meaningful classifications to these features based on their inferred relevance to the underlying biological processes is of significant utility. While numerous models have emerged for the prioritization of genes, they largely overlook feature classification, including the crucial distinctions between the “driver,” “biomarker,” and “mediator” feature classification paradigms we introduce in this study (Box 1). Recent advances in inferring transcription factor activity from multimodal single-cell data offer promising avenues for such classifications, enabling prediction of drivers in cellular transitions or disease states by integrating regulatory networks with expression profiles.ref. bib20^,ref. bib21^,ref. bib22^,ref. bib23 Here, “drivers” denote data-driven, process-level features inferred to influence biological state transitions, rather than genetically defined oncogenic driver mutations.

To address the above shortcomings, we developed ExIR (experimental data-based integrative ranking), a workflow that combines both supervised and unsupervised ML models as well as graph-based influence analyses to perform data-driven classification and prioritization of candidate features (e.g., genes, proteins, etc.). Uniquely, ExIR performs these tasks using solely the high-throughput experimental data under investigation (i.e., without reliance on external sources of knowledge). Benchmarking using several different transcriptomic and proteomic datasets demonstrates that ExIR outperforms other prioritization methods in specifically and sensitively identifying, classifying and ranking features by their functional importance to the biological process represented in a given dataset. Moreover, we demonstrate the utility of ExIR to decipher complex biological processes by applying the workflow to a transcriptomic dataset of a zebrafish model of the pediatric lysosomal storage disease MPS IIIA,ref. bib26 facilitating identification of several mechanistic drivers of MPS IIIA disease progression as well as potentially clinically useful candidate diagnostic biomarkers.

Results

ExIR is a versatile data-driven framework for feature classification and prioritization

We developed ExIR, a data-centric framework that classifies and prioritizes biological features directly from high-throughput experimental data without reliance on external annotations or curated knowledge sources (Figure 1; STAR Methods). Starting from normalized input data, ExIR constructs feature-feature association networks using correlation-based adjacency matrices that capture global co-variation structure across samples. These adjacency matrices are subsequently used for network reconstruction and influence estimation. Supervised random forest modeling identifies features with predictive relevance to experimental conditions, while unsupervised principal component analysis captures dominant variance contributions independent of class labels. Together, these scores report feature influence and prioritization within the reconstructed network (Figure 1A; STAR Methods). ExIR integrates complementary supervised and unsupervised ML approaches with network-based influence analysis to derive multiple, independent feature-level scores capturing differential behavior, predictive importance, variance structure and network context. These data-derived scores are aggregated to classify features as drivers, biomarkers and mediators, followed by ranking within each class (Figure 1B; Table 1; STAR Methods). In this way, ExIR identifies features that function primarily as process-level drivers, diagnostic biomarkers, or network mediators, while also allowing biologically meaningful overlap between driver and biomarker roles.

Unlike other gene/protein prioritization methods such as Endeavourref. bib28 and ToppGene,ref. bib29 which require a seed or training set to derive correlated features, ExIR relies only on interrogation of self-contained experimental data (Figure 1) for the detection of drivers, biomarkers, and mediators. The random forest algorithm prioritizes genes that would significantly affect inter-group variation between samples; ExIR then applies principal component analysis to weight features with a greater potential to differentiate all samples from each other, regardless of which group they belong to. Accordingly, those genes that have more sensitivity to inter-group variation and to phenotypic distinction within each group (i.e., intra-group variation) are assigned a higher rank.

ExIR outperforms existing feature prioritization methods in driver ranking

To examine ExIR’s performance in ranking influential driver genes/proteins against existing feature prioritization methods, we used a diverse set of published high-throughput transcriptomic and proteomic datasets representing several disease states (prostate cancer,ref. bib30 head and neck squamous cell carcinoma [HNSCC],ref. bib31^,ref. bib32 glioblastoma,ref. bib33 schizophrenia,ref. bib34 breast cancer,ref. bib32^,ref. bib35 thyroid cancer,ref. bib32 lung adenocarcinoma,ref. bib32^,ref. bib36 stomach adenocarcinoma,ref. bib32 lung squamous cell carcinoma,ref. bib32 and hepatocellular carcinomaref. bib32). The use of pathology-associated data facilitates easier examination of model performance, as the transcriptional and proteomic changes associated with these pathological states are well-defined and provide a consistent ground truth. We additionally sought to examine ExIR’s performance outside of a disease context, and thus included a dataset representing a physiological process (oogenesisref. bib37). Collectively, the experimental modalities tested included bulk- and single-cell RNA sequencing, RNA microarray, as well as mass spectrometry-based proteomics datasets. In all cases, and to facilitate comparisons with methods that require training feature sets (i.e., Endeavor and ToppGene), benchmarking was performed by first assigning ground-truth (i.e., true positive and true negative) drivers on the basis of a consistent set of criteria from both gene expression/protein abundance data and ontological annotations (Figures 2A and 2B and STAR Methods).

Benchmarking using receiver operating characteristic (ROC) analyses demonstrated that ExIR consistently achieved the highest driver feature prioritization performance across datasets (Figure 2C). Across all prioritization benchmarks, ExIR attained a mean area under the ROC curve (AUC) of 0.96 (Figure 2D), exceeding the performance of ToppGeneref. bib29 (mean AUC = 0.84), log₂ fold change (0.77), Endeavourref. bib28 (0.73) and GeneMANIAref. bib38 (0.27). ExIR was the top-performing method in all datasets except TCGA THCA, where ToppGene marginally outperformed ExIR. These results indicate that ExIR provides robust and accurate prioritization across diverse experimental contexts while remaining independent of external knowledge sources. The low ROC performance observed for GeneMANIA reflects its reliance on external interaction networks and functional annotations,ref. bib38 which are not optimized for prioritizing experimentally derived genes and proteins in this benchmarking context.

ExIR improves biomarker gene prioritization relative to existing methods

Likewise, for biomarker gene prioritization, true positive and true negative biomarker sets were defined using a uniform set of criteria across both cancer and non-cancer disease datasets (Figures 3A and 3B and STAR Methods). ROC-based benchmarking demonstrated that ExIR consistently achieved superior performance in biomarker feature prioritization across datasets (Figure 3C). In particular, ExIR achieved a mean AUC of 0.87 (Figure 3D), outperforming Student’s t test (0.73), mutual informationref. bib39 (MI) (0.72), and correlation-based approaches (0.53).

In addition to generally outperforming all other feature classification and ranking methods, our application of ExIR to the Cancer Genome Atlas lung adenocarcinoma RNA-sequencing dataset (TCGA LUAD) identified several biomarkers (Table 2), many of which fail to be detected as top-ranked biomarkers using other methods. Additionally, since all of the top five LUAD biomarkers prioritized by ExIR were down-regulated genes, we rationalized that detection of an up-regulated biomarker might represent a more traceable feature in a diagnostic context. We proceeded to assess the top five identified up-regulated biomarkers prioritized by ExIR; except for SFTPC, PYCR1, TEDC2, and TOP2A (each inconsistently prioritized by various methods), none of the other known or ExIR-predicted down-/up-regulated LUAD biomarkers were consistently identified as top-ranked biomarkers by methods benchmarked against ExIR (Table 2).

As the above biomarkers were retrieved from interrogation of an RNA sequencing dataset, and appreciating that quantified gene expression may not always accurately reflect protein abundance,ref. bib40 we next sought to determine whether the ExIR-derived biomarker genes indeed evoke changes at the protein level in LUAD samples. Examination of public immunohistochemical data for top-ranked ExIR biomarkers in LUAD from the Human Protein Atlasref. bib41 revealed differential immunoreactivity of all previously established LUAD biomarkers between normal and cancer samples (Figures 4A–4E). Furthermore, the top five ExIR-predicted biomarkers (all down-regulated at the transcript level) exhibited pronounced immunoreactivity within healthy pneumocytes, yet were largely absent from tumor cells (Figures 4A–4F(i)). Subsequent examination of the top five ExIR-prioritized upregulated LUAD biomarkers likewise revealed differential immunoreactivity of these biomarkers between normal and LUAD samples (Figures 4J–4N). Collectively, this demonstrates ExIR’s capacity for sensitive identification of disease-associated biomarkers.

ExIR identifies disease-associated driver and biomarker genes in a zebrafish model of MPS IIIA

Mucopolysaccharidosis IIIA (MPS IIIA, or Sanfilippo syndrome A; OMIM #252900) is a lysosomal glycosaminoglycan storage disorder causing progressive neurological decline and onset of dementia-like symptoms in children, often leading to death during teenage years. Few effective therapeutic options are presently available for MPS IIIA patients; this is partially due to an incomplete understanding of the complex nature of MPS IIIA pathophysiology despite its monogenic etiology. To decipher the molecular features underlying MPS IIIA, we applied ExIR to a transcriptomic dataset of adult brains from the sgsh^Δex5-6 zebrafish mutantref. bib26 (ZFIN line designation sgsh^mnu301, hereafter referred to as sgsh), a zebrafish MPS IIIA model we previously generated. Interrogating the ExIR rankings for the 50 top-ranked driver and biomarker genes identified 67 distinct features; 33 genes were top-ranked in both classifications, with 17 genes each being top drivers but not biomarkers, and vice versa (Table S1). Examining the relationship between ExIR classification of differentially expressed genes (DEGs) and their differential expression levels demonstrated that identification of top drivers and biomarkers did not necessarily correlate with the DEGs exhibiting the largest fold-changes (Figure 5A). This further demonstrates that the network-based approach to gene classification and ranking is agnostic to basic parameters of differential gene expression analysis. Strikingly, of the 33 genes ranked as both top drivers and biomarkers detected in the transcriptome of sgsh brain, the vast majority of them clustered among the up-regulated DEG cohort (Figure 5A).

To further investigate the nature of ExIR-prioritized drivers and biomarkers, these features were next subdivided on the basis of fold-change directionality. Reflecting the above concordance between top drivers and biomarkers, the top two accelerating drivers—ppargc1b and baiap2b (a transcriptional coactivator reported to positively regulate mitochondrial biogenesisref. bib42 and an adapter protein involved in cytoskeletal reorganization of actin filaments,ref. bib43 respectively)—were also the top two up-regulated biomarkers (Figure 5B), while the subsequent top-ranked accelerating drivers (tasor2, ifrd1, and rnf144ab) were also among the top-ranked up-regulated biomarkers (ranks 9, 6, and 8, respectively). Similarly, rhpn2, mmadhca, and plekhm1 (up-regulated biomarker ranks 3, 4, and 5) were among the top-ranked accelerating drivers (8, 9, and 6, respectively; Table S2). This relationship was also observed, albeit to a lesser extent, among the decelerating drivers and down-regulated biomarkers. btg2 and npas4a (both immediate-early genes [IEGs] associated with acute transcriptional responses to neuronal activityref. bib44) were both the top two decelerator drivers and down-regulated biomarkers (Figure 5C), though subsequent top-ranked decelerator drivers (nr4a1, fosab, and egr4, also well-known IEGs) tended to be lower-ranked down-regulated biomarkers (32, 33, and 34, respectively; Table S2). This is somewhat expected, given the observed bias in this dataset toward representation of up-regulated genes among top drivers and biomarkers (Figure 5A).

To further explore the nature of differential expression of the genes ranked as top drivers and biomarkers by ExIR in the sgsh zebrafish brain, we performed in situ hybridization against transcripts for the accelerating/up-regulated genes ppargc1b and baiap2b, and the decelerating/down-regulated genes btg2 and npas4a in the telencephalon. As anticipated by its ranking as both a top up-regulated biomarker and accelerating driver in the sgsh CNS, ppargc1b was detectable in the sgsh telencephalon, but not that of wild-type siblings (Figure 5D(i and ii)); in the sgsh telencephalon, ppargc1b expression was strongest in the subventricular domain of the ventral subpallium (Figure 5D(ii)), with weaker expression in both the medial subventricular domain of the pallium and the parenchyma of the lateral pallium (Figure 5D(ii)). Similarly, little to no baiap2b expression was detected in the wild-type telencephalon (Figure 5D(iii)), while robust expression was observed in the sgsh telencephalon in the anterior parvocellular preoptic nucleus associated with the anterior-most hypothalamic region (Figure 5D(iv)), as well as the supracommissural nucleus of the ventral telencephalic area (Figure 5D(iv)).

Conversely, expression of down-regulated drivers and biomarkers btg2 and npas4a were robustly detected in wild-type telencephalon, but were either reduced in a spatially distinct manner (in the case of btg2), or absent in expression (as with npas4a), in the sgsh brain. Wild-type expression of btg2 was observed in abundant neurons in both dorsal, medial, and ventral subdivisions of the subventricular domain, with sparser expression in the posterior zone of the dorsal telencephalic area (Figure 5D(v)). Subventricular btg2 expression was largely abolished in the sgsh telencephalon, though expression in the dorsolateral pallium was partially preserved (Figure 5D(vi)). Similarly, while npas4a expression was detected in all major subventricular domains of the wild-type telencephalon (Figure 5D(vii)), expression was effectively completely absent from these regions in homozygous sgsh siblings (Figure 5D(viii)).

Taken together, application of ExIR to a transcriptomic analysis of the brain of a zebrafish MPS IIIA model identified and prioritized features strongly characteristic of the disease state. That the abnormal expression patterns of all tested top-ranked genes were associated with specific neuronal domains suggests that certain neuronal subtypes or neuroanatomical regions may be variably susceptible to functional perturbation in MPS IIIA. By revealing specific neuronal domains susceptible to functional and molecular perturbations in MPS IIIA and prioritizing crucial features of the disease, our application of ExIR to this dataset opens avenues for deeper spatial and phenotypic exploration of the mechanisms underpinning MPS IIIA pathogenesis and progression.

ExIR identifies genes mediating the transduction of driver signals

An additional utility of ExIR is its capacity to assign features to a class we refer to as “mediators.” Mediators are features that are important for the flow of information between prioritized drivers in the ExIR-generated network (Figure 6A). Mediators may either be DE or not, depending on the nature of the experimental design employed (see Box 1). Most basically, mediators are never DE if an experiment contains only two comparator groups; detection of DE mediators requires an experimental design of >2 groups so as to fulfill the requirement for a given feature to be inconsistently DE (i.e., DE in some but not all pairwise comparisons within the dataset) and/or show opposing differential expression between different pairwise comparisons. Altogether, mediators exhibit either non-differential, or inconsistent/opposing differential expression/abundance patterns across multi-time point/multi-condition experiments (Figure 6B).

Mediators are predicted to directly and/or indirectly associate with major drivers of the biological process in question, inferred from their interactions with drivers within the ExIR network. However, they themselves do not necessitate a consistent differential expression/abundance pattern between conditions. Importantly, mediator prioritization reflects driver-contextualized network influence rather than nonspecific high-degree connectivity, distinguishing biologically meaningful signal transducers from dense correlation hubs. To our knowledge, only one other model (machine learning-assisted network inference, MALANI)ref. bib45 has so far been described with the capacity to compute a class of genes akin (albeit not entirely identical) to our definition of mediator features. Additionally, MALANI has so far only been utilized in the context of transcriptomic analyses of cancer.ref. bib45 To benchmark ExIR’s performance in detecting mediator features (which MALANI refers to as a “class II” gene), we applied ExIR to the same TCGA BRCA RNA-sequencing dataset as was previously interrogated with MALANI,ref. bib45 and performed GO- and KEGG-based overrepresentation analyses on identified ExIR mediators/MALANI class II features. ExIR significantly outperformed MALANI in identification of mediators associated with breast cancer-related GO biological processes (GO-BPs) and KEGG pathways. The mediator-biological term association analysis was conducted using data from the Comparative Toxicogenomics Database.ref. bib46 Extending this analysis to mediator identification across additional disease-related datasets revealed that, on average, more than 55% of GO-BP terms and KEGG pathways enriched among ExIR-derived mediators are implicated in the corresponding disease (Figures 6C–6F; Tables S3 and S4). Taken together, identification and classification of mediator nodes in the ExIR network allows for more thorough deconvolution of complex biological and pathological processes from high-throughput datasets.

ExIR identifies potential non-DE transcriptional mediators in neuropathology of the MPS IIIA zebrafish

As all gene classes in ExIR—including mediators—are classified in an unbiased fashion, we sought to explore functional enrichment among the identified mediators in the ExIR network pertaining to the MPS IIIA zebrafish brain transcriptome. The MPS IIIA zebrafish RNA-seq dataset employed in this study comprised only two cohorts, and thus can only detect and rank non-DE mediators (Figure 6B). From the 22,022 genes detected in this RNA-seq dataset, 704 genes were defined as DE (in this case, differential expression is defined solely as Q < 0.05, with no log-fold change cut-off applied) and 21,318 non-DE features. Of the latter, ExIR identified 5,577 features as non-DE mediators based on their association with driver genes in the generated network (Figure 6G).

First, to account for limitations associated with zebrafish Gene Ontology annotationsref. bib47^,ref. bib48^,ref. bib49 (compared to the more thoroughly curated human annotations), all zebrafish mediators were first converted to human orthologues, and the resulting gene set was then interrogated against the Gene Ontology resource.ref. bib50 Strikingly, hierarchical clustering of significantly enriched (Q < 0.05, Benjamini-Hochberg p value adjustment) GO terms across all top-level ontology classes identified several GO term clusters of direct relevance to MPS IIIA pathology; these included terms associated with catabolism of the sugar/glycosaminoglycan substituents uronic acid and hexose. Further, among the significantly enriched terms were several associated with neural development and axonal projection (Figure 6H), likely relevant to the overt neuropathology characteristic of MPS IIIA. Collectively, this analysis demonstrates that ExIR is able to identify mediators which, despite their non-differential expression, are likely important lower-order regulators of identified drivers.

Discussion

ExIR is a versatile model that simultaneously extracts, classifies, and prioritizes candidate features from high-throughput experimental data independent of external, curated knowledge sources. ExIR initially employs ML and graph-based influence evaluation techniques, then performs several filtration steps on the input normalized data to derive several different ranking scores. Combinations of these ranking scores are then used to derive summary scores for feature classification and ranking as drivers, biomarkers, and mediators of biological processes. From a biological perspective, manipulation of driver features prioritized by ExIR (potentially alongside co-manipulation of associated mediators) could have the most prominent impact on the progression of a biological process/disease as well as phenotype manifestation, thereby prioritizing candidates for therapeutic targeting. While highly ranked ExIR biomarkers are anticipated to have the highest sensitivity to a biological condition and the severity of the phenotype, their manipulation may not necessarily affect the progression of the process in the way highly ranked driver features might.

Our comparative analyses show that ExIR robustly outperforms existing tools and algorithms in the unbiased prioritization of ground truth driver, biomarker, and mediator features across diverse datasets representing multiple distinct experimental modalities. The underlying reasons for the superior accuracy and robustness of the ExIR model include (1) co-implementation of both supervised and unsupervised ML techniques; (2) integration of both ML and network-based models; (3) optimization of mathematical operations for score integration; and (4) independence from potentially-confounding external sources of information.

In order to experimentally validate ExIR’s utility in deciphering complex disease states, we employed this workflow to a transcriptomic analysis of the first zebrafish model of the mucopolysaccharidosis Sanfilippo syndrome Type A (MPS IIIA).ref. bib26 The etiology of mucopolysaccharidoses, specifically the deficiency in glycosaminoglycan (GAG) catabolism, particularly heparan sulfate (HS) in the case of MPS IIIA, has long been established. However, the precise mechanism by which this deficiency results in the distinctive progressive functional neurodegeneration observed in individuals with mucopolysaccharidoses remains elusive. Several studies have implicated multiple factors contributing to disease progression, including microgliosis and neuroinflammation,ref. bib26^,ref. bib51^,ref. bib52 CNS atrophy,ref. bib53 and more recently an increasing recognition of the role of functional synaptic impairments.ref. bib54^,ref. bib55^,ref. bib56 However, interventions targeting these factors have generally not been sufficient to ameliorate the severity of MPS IIIA or halt its progression. A systems-level analysis of the MPS IIIA CNS provides an opportunity to perform an unbiased survey of the molecular changes associated with disease onset and progression, which in turn may lead to the rational development of more effective, evidence-guided therapies. However, such analyses inevitably produce large volumes of information that require manual curation and filtration for interpretability; this often represents a significant analytical bottleneck. ExIR, by classifying and ranking features in an unbiased manner, provides an effective means of feature prioritization which may aid in focusing research efforts toward the most important biological changes in a given assay. Indeed, application of ExIR to the MPS IIIA zebrafish transcriptome has provided an unbiased basis for deconvolution of the complex transcriptional signature of the disease state, facilitating refinement of distinct candidate features and pathways for further study of the most critical elements of MPS IIIA neuropathology.

Taken together, ExIR establishes a reference-free framework for the classification and prioritization of genes and proteins directly from experimental data, without reliance on external annotations or curated knowledge bases. By integrating network-based inference with data-driven ranking, ExIR provides a generalizable strategy applicable to bulk and single-cell transcriptomic as well as proteomic datasets. Notably, the introduction of a distinct class of “mediator” genes expands the analytical landscape beyond conventional driver and biomarker identification, enabling more nuanced interrogation of regulatory hierarchies underlying complex biological processes and disease progression.

As high-throughput technologies continue to generate increasingly complex and high-dimensional datasets, approaches that distil interpretable and experimentally actionable signals are essential for making optimal use of these data. By streamlining feature classification and prioritization, ExIR facilitates hypothesis generation and accelerates downstream validation efforts. Its standalone and annotation-independent design positions it for broad applicability across basic, translational and clinical research settings, including precision medicine initiatives and rare disease studies where robust identification of candidate drivers and biomarkers is critical.

Limitations of the study

Despite ExIR’s performance and generality, it has several limitations. As a purely data-driven and unsupervised framework, its prioritization is inherently dependent on the quality, depth, and statistical power of the input data, with smaller sample sizes potentially limiting the robustness of correlation-based network inference. While ExIR incorporates adjustable correlation and mutual rank filtering strategies and integrates multiple independent evidence streams to mitigate spurious associations, unmodeled technical variation or batch effects may still influence results if not adequately controlled upstream. In its current implementation, ExIR operates on a single data modality at a time and therefore does not capture cross-modal regulatory relationships arising from concurrent transcriptomic, proteomic, or chromatin accessibility measurements. Extending ExIR to natively integrate multimodal data, including joint RNA-seq, proteomic, and ATAC-seq profiles, represents an important direction for future development and will enable more comprehensive inference of regulatory drivers and mediators underlying complex biological processes.

Resource availability

Lead contact

Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Peter D. Currie (peter.currie@monash.edu).

Materials availability

This study did not generate new unique reagents.

Data and code availability

• All datasets analyzed in this study are publicly available. The Sanfilippo zebrafish RNA sequencing dataset generated in this study has been deposited in the NCBI Gene Expression Omnibus (GEO) and is publicly available as of the date of publication under accession number GSE304151. Additional datasets used include GBM scRNA-seq (GSE84465), oogenesis scRNA-seq (GSE130212), HNSCC scRNA-seq (GSE103322), schizophrenia microarray data (GSE93577), the Prostate Cell Atlas dataset (EGAS00001005787), TCGA datasets accessed via the Broad GDAC Firehose (https://gdac.broadinstitute.org), as well as proteomics datasets including the LUAD dataset (IPX0001804001, iProX), and the BRCA dataset (PXD002057, PRIDE/ProteomeXchange). Accession numbers are listed in the key resources table.
• All code used in this study is publicly available. The ExIR method is implemented in the influential R and Python packages and is accessible via GitHub (https://github.com/asalavaty/influential; https://github.com/asalavaty/python-influential), CRAN (https://cran.r-project.org/package=influential), and PyPI (https://pypi.org/project/influential/). Documentation and usage examples are available in the package vignettes.
• Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Acknowledgments

The authors would like to thank Lan Nguyen for their constructive feedback on the development of ExIR and proteomic data analysis, and Monash AquaCore, Monash Micro Imaging, and the Monash Genomics and Bioinformatics Platform facilities for their excellent support. The results shown in this study are in part based on data generated by the TCGA Research Network (http://cancergenome.nih.gov/). This work was supported by 10.13039/501100001779Monash University, a 10.13039/501100000925National Health and Medical Research Council (NHMRC) Fellowship (1136567) to P.D.C., an 10.13039/501100000925NHMRC grant (1180905) to M.R., an NHMRC grants (GNT2013305, GNT2037953) to J.K., 10.13039/501100002341Research Council of Finland and 10.13039/501100006306Sigrid Juselius Foundation grants to J.K., an incubator grant from the 10.13039/100013918Sanfilippo Children’s Foundation (Australia) and the Cure Sanfilippo Foundation (US) to J.K., a Translational grant from the 10.13039/100013918Sanfilippo Children’s Foundation (Australia), the Cure Sanfilippo Foundation (US), Fundacja Sanfilippo (Poland), Sanfilippo Initiative (Germany), and H.A.N.D.S (Portugal, Spain and France) to J.K. A.D. and A.S. are supported by 10.13039/100015539Australian Government Research Training Program (RTP) scholarships. The Australian Regenerative Medicine Institute is supported by grants from the 10.13039/100007222State Government of Victoria and the 10.13039/100015539Australian Government.

Author contributions

A.S., A.D., J.K., M.R., and P.D.C. conceptualized the study; A.S. developed ExIR and performed in silico benchmarking and data analysis; A.D. performed all in vivo experiments and analysis; A.S. and A.D. wrote the manuscript with contributions from J.K., M.R., and P.D.C.

Declaration of interests

The authors declare no competing interests.

Declaration of generative AI and AI-assisted technologies in the writing process

During the preparation of this work, the authors used LLMs in order to improve readability and language. After using this tool, the authors reviewed and edited the content as needed and take full responsibility for the content of the published article.

STAR★Methods

Key resources table

Experimental model and study participant details

All protocols and procedures using zebrafish of ages greater than 7 dpf were approved by the Monash University Animal Ethics Committee (ERM14481, ERM22161 and ERM17993). Zebrafish were maintained under standard housing and breeding conditionsref. bib57 in the AquaCore facility, Monash University. Embryos and larvae were maintained in E3 medium, while adult zebrafish were maintained in system water. Adult zebrafish (3-months-old) were randomly selected for inclusion in experiments without consideration of sex, which is not anticipated to be an important covariate in this study. Mutant strains used in this study were sgsh^mnu301 (referred to either as sgsh^Δex5-6 or simply sgsh).ref. bib26

Method details

Overview of the ExIR framework

ExIR is a data-driven computational framework designed to classify and prioritize biological features (for example, genes or proteins) as drivers, biomarkers, or mediators of biological processes directly from high-throughput experimental data. ExIR operates on normalized feature-by-sample matrices and associated differential or regression statistics derived from the same experiment. The framework comprises two conceptual stages: (i) extraction of multiple complementary feature-level scores capturing differential behavior, predictive importance, variance structure, and network influence; and (ii) integration of these scores using biologically-motivated additive and synergistic rules to classify and rank features according to their putative functional roles.

Unless otherwise stated, all computations were performed in R, using user-supplied, pre-normalized data.

Input data and preprocessing

ExIR requires two mandatory inputs. The first is a table of feature-level differential or regression statistics, such as log₂ fold changes, regression coefficients, or associated statistical significance values derived from any standard differential expression or abundance analysis. In ExIR’s current form, these statistics may originate from bulk RNA sequencing, single-cell RNA sequencing, microarray, or proteomic assays. In multi-condition or time-course experiments, this table may contain statistics from multiple pairwise comparisons or regression analyses.

The second required input is a normalized feature-by-sample matrix representing the complete experimental dataset, in which rows correspond to features (genes or proteins) and columns correspond to samples or cells. ExIR does not impose a specific normalization strategy; instead, users are expected to apply an appropriate normalization method upstream, consistent with the experimental modality. An optional log-transformation may be applied within ExIR if the input data are not already log-scaled. ExIR deliberately does not perform data normalization, batch correction or differential analysis internally, as these preprocessing steps are highly modality- and study-specific and typically require careful, context-dependent parameterization and quality control by the user. This design ensures that ExIR remains agnostic to experimental platform while avoiding the introduction of inappropriate preprocessing assumptions.

Optionally, users may provide a third input consisting of a custom list of features of interest. When supplied, this list is used to restrict downstream analyses to a predefined feature subset, provided that all listed features are present in the normalized data matrix.

Additive and synergistic score integration strategy

ExIR integrates heterogeneous sources of evidence using two complementary principles, additive integration and synergistic integration. These reflect distinct biological assumptions:

Additive integration is applied when combining multiple independent observations of the same type of evidence, such as differential statistics across contrasts, regression coefficients across conditions, or neighborhood-level influence scores. In this setting, additive integration captures the cumulative magnitude or consistency of evidence without assuming interaction effects. Examples include the calculation of differential scores across multiple comparisons, regression scores across ordered conditions, and neighborhood driver scores computed as the sum of primitive driver scores of all first-order network neighbors. In experiments with a single comparison (for example, two-condition, single–time point studies), additive integration is omitted where not applicable.

Synergistic integration is applied when combining distinct, biologically orthogonal evidence streams for the same feature, under the assumption that concordant signals jointly strengthen functional relevance. In these cases, scores are combined multiplicatively to prioritize features that are simultaneously significant, influential, and predictive, while penalizing features supported by only a single evidence source. Synergistic integration is used in the calculation of primitive driver scores, final driver scores, biomarker scores, and both differentially expressed and non-differentially expressed mediator scores, as described below.

Calculation of differential, regression, and significance scores

For each feature, ExIR computes summary scores that quantify the magnitude, consistency, and statistical confidence of differential behavior across experimental conditions.

For experiments with two or more conditions, a differential score is calculated for each feature by summing the absolute values of user-supplied differential statistics (for example, log₂ fold changes) across all relevant comparisons:

In multi-condition or time-course experiments, features are further subdivided into consistently altered and inconsistently altered categories based on whether the direction of change is preserved across comparisons. Separate differential scores are computed for each category.

For experiments involving ordered conditions or time points, an optional regression score is calculated for each feature by summing user-supplied regression coefficients or coefficients of determination (R² values).

Statistical confidence is quantified using a significance score, computed as the sum of the negative log-transformed significance values (for example, adjusted P-values or false discovery rates):

As with differential scores, significance scores are computed separately for consistently and inconsistently altered features in multi-condition designs.

Supervised machine learning for feature importance estimation

To estimate predictive importance, ExIR employs random forest classification,ref. bib58 using sample-level labels (for example, disease versus control) as the response variable and feature expression or abundance values as predictors. Each sample constitutes one observation, and each feature constitutes one predictor variable.ref. bib59

Random forest models were implemented using the ranger R package.ref. bib60 The number of trees was set to a large default value (10,000) to ensure stable convergence of impurity-corrected feature importance estimates. Because ExIR relies on relative ranking and permutation-based significance rather than absolute importance magnitudes, prioritization outcomes are robust to reasonable variation in ensemble size once importance estimates have stabilised. The number of variables considered at each split set to the square root of the total number of features. Feature importance was quantified using impurity-corrected importance scores, yielding a single importance value per feature that reflects its contribution to accurate sample classification.

To assess the statistical significance of feature importance scores, ExIR applies a permutation-based approach, in which sample labels are permuted to generate null distributions of importance values.ref. bib61 By default, 100 permutations are used, providing a balance between computational efficiency and robust P-value estimation. P-values are calculated following the method of Altmann et al.,ref. bib62 and transformed into a supervised-learning significance score using negative log₁₀ transformation.

Following supervised learning, features may be further filtered based on statistical significance thresholds, reducing noise in downstream analyses.

Unsupervised machine learning for variance-driven feature prioritization

To capture unsupervised variation in the data, ExIR performs principal component analysis (PCA) on the filtered, normalized feature-by-sample matrix. PCA is conducted using base R functions.

For each feature, ExIR extracts the absolute loading (rotation value) on the first principal component, which captures the largest proportion of variance across samples. Higher-order components were not considered to avoid introducing dataset-specific noise or technical variation. Features with high absolute rotation values are interpreted as strong contributors to sample separation along the dominant variance axis. These rotation values are used to increase the sensitivity of biomarker prioritization by identifying features that explain major sources of biological variation independently of class labels.

Construction of gene–gene association networks

To infer feature–feature associations, ExIR computes pairwise Spearman correlations between all features across samples using a matrix-based linear algebra formulation, resulting in a symmetric feature-by-feature correlation matrix. Correlations are computed across the entire dataset, independent of experimental conditions, to capture global association structure. Both positive and negative correlations were considered by thresholding on the absolute correlation coefficient. By default, associations with |r| ≥ 0.5 were retained, balancing network sparsity and biological interpretability; however, this threshold is user-adjustable and may be relaxed or tightened depending on sample size and expected regulatory density. As an alternative filtering strategy, ExIR supports mutual rank (MR)–based thresholding, in which associations are retained if their mutual rank falls below a user-defined threshold (default MR ≤ 20). MR filtering reduces the influence of indirect or unstable correlations and is particularly advantageous for large or heterogeneous datasets.

From this global matrix, a refined adjacency matrix is constructed by retaining correlations involving features present in the filtered dataset. First-order associated features are identified, and second-order associations are subsequently retrieved by iteratively filtering the global correlation matrix. This stepwise approach yields a final adjacency matrix comprising first- and second-order feature associations, enabling efficient network reconstruction while limiting computational complexity.

Network reconstruction and influence estimation

The final adjacency matrix is used to reconstruct a feature association network using the igraph R package. Within this network, the Integrated Value of Influence (IVI) is computed for each node using the influential R package (https://cran.r-project.org/package=influential).ref. bib27 IVI integrates multiple centrality and connectivity measures into a single score that quantifies the overall influence of each feature within the network.

Driver score calculation

To identify candidate driver features, ExIR computes a primitive driver score for each feature by multiplying its IVI value by the significance score of consistently altered features. This score captures features that are both statistically perturbed and highly influential within the network.

To account for local network context, ExIR further computes a neighborhood driver score by summing the primitive driver scores of all first-order neighboring features. The final driver score is calculated as the product of the primitive and neighborhood driver scores, prioritizing features that are influential, consistently altered, and embedded within clusters of other influential features.ref. bib63 Driver scores quantify the likelihood that a feature functionally contributes to progression between experimental conditions based on integrated statistical and network evidence, and do not imply direct genetic causality. Driver features are further subclassified as accelerator or decelerator drivers based on the direction of differential change.

Mediator score calculation

Mediator features are classified into differentially expressed (DE mediators and non-differentially expressed (non-DE) mediators. Mediator prioritization requires concurrent satisfaction of high IVI and driver-enriched neighborhood scores, ensuring that mediators arise from driver-contextualized network influence rather than from nonspecific high-degree connectivity in dense correlation structures. Also, because ExIR mediators are inferred from network influence and neighborhood context rather than from pairwise causal models, they are ranked by relative mediator scores rather than assigned a single mediation effect size.

DE mediators are features that exhibit inconsistent differential behavior across conditions but maintain high network influence and proximity to driver features. Their scores integrate significance scores of inconsistently altered features with IVI and neighborhood driver scores.

Non-DE mediators are features that do not show statistically significant differential expression but are prioritized based on their IVI values and neighborhood driver scores, reflecting potential regulatory or coordinating roles within the network.

Biomarker score calculation

Biomarker prioritization integrates multiple complementary criteria. For each feature, a biomarker score is computed by multiplying its differential score, regression score (if applicable), significance score, supervised learning importance and significance scores, and unsupervised PCA rotation value. This composite score prioritizes features that are consistently altered, statistically robust, predictive of sample class, and explanatory of major variance in the data. Biomarkers are subclassified as up- or down-regulated based on directionality.

Statistical scaling and multiple-testing correction

All derived scores are standardized using Z-score transformation.ref. bib64 P-values are computed from Z-score distributions and adjusted for multiple testing using the Benjamini–Hochberg procedure to control the false discovery rate.

Data preparation for ExIR evaluation and benchmarking

For evaluation and benchmarking of ExIR, the following criteria were required for selection of datasets:

• The dataset should correspond to a disease/biological process with >40 curated driver genes in the DisGeNET database or Gene Ontology resource;
• The dataset should have >100 samples/cells;
• The dataset should include at least 2 conditions (e.g., diseased vs. unaffected, different time points);
• The dataset should have a comparable number of samples/cells within each condition.

Ground-truth driver and biomarker sets were selected from context-specific, independent annotation resources (GO-BP for normal biological processes; DisGeNET/MutPanning for disease drivers; NCI EDRN and disease-specific biomarker databases for biomarkers), and were used exclusively for downstream evaluation rather than model construction, thereby minimizing annotation leakage and database-specific bias.

Accordingly, the following datasets were used to evaluate ExIR and demonstrate its applicability in the classification and prioritization of features from microarray, bulk and single-cell RNA-seq (scRNA-seq), and proteomics experimental modalities.

Glioblastoma (GBM) dataset

GSE84465 is a scRNA-seq dataset generated using plate-based protocols from four patients with confirmed cases of primary GBM and comprises 3,589 cells.ref. bib33 Tissue samples for this dataset originated from either the tumor core or the peritumoural cortical space. Major classes of cells, including neoplastic and non-neoplastic cells, were identified using immunopanningref. bib65 and further confirmed by comparison with other single-cell and bulk RNA-seq data. The dataset was originally normalized based on the read counts generated by HTSeqref. bib66 and filtration of genes with very low counts. Amongst all cell types previously identified, we selected two major subsets of cells including periphery regular (non-neoplastic) cells (n=1,184 cells) as the ‘normal’ set and tumor neoplastic cells (n=1,029 cells) as the cancer set, and filtered out genes with zero counts across all selected cells for downstream analysis.

Prostate Cell Atlas dataset

EGAS00001005787 is a scRNA-seq datasetref. bib30 generated using the 10x Chromium Single Cell 3’ v2 protocol from 10 patients aged 50-72 undergoing image-guided biopsies for suspicion of prostate cancer and comprises 15,492 cells. Tissue samples for this dataset originated from paired cancer biopsies and adjacent normal prostate tissue. Major classes of cells, including epithelial, stromal, and immune cells, were identified through unbiased clustering and validated by flow cytometry and confocal imaging. The dataset was originally processed using CellRanger for alignment and UMI counting, with QC filtration removing cells with fewer than 200 or more than 2,500 genes, mitochondrial content exceeding 30%, and genes expressed in fewer than 3 cells. Amongst all cell types previously identified, we selected luminal epithelial cells including normal (n=3,078 cells) as the ‘normal’ set and tumor (n=5,088 cells) as the cancer set, and filtered out cells with fewer than 200 genes and genes expressed in fewer than 3 cells for downstream analysis.

HNSCC scRNA-seq dataset

GSE103322 is a scRNA-seq datasetref. bib31 generated using a modified SMART-Seq2 protocol from 18 patients with primary oral cavity HNSCC tumors (including five with matched lymph node metastases) and comprises 5,902 cells. Tissue samples for this dataset originated from primary tumors and matched lymph node metastases in treatment-naïve patients. Major classes of cells, including neoplastic (n=2,539) and non-neoplastic (n=3,363) cells, were identified by inferring large-scale chromosomal copy-number variations (CNVs) from averaged expression across chromosomal intervals, epithelial marker expression scoring, and global expression patterns via clustering.ref. bib31 The dataset was originally processed by aligning reads to the hg19 human genome reference, quantifying expression as log2(TPM/10+1).ref. bib31

Oogenesis dataset

GSE130212 is a scRNA-seq dataset generated using plate-based protocols from fetal mouse ovaries at three developmental stages (E12.5, E14.5 and E16.5), together encompassing 19,144 FACS-sorted high-quality murine female germ cells.ref. bib37 The dataset was originally normalized and DEGs between all time-points were detected in Seurat.ref. bib67

TCGA BRCA

TCGA BRCA is a bulk RNA-seq of breast cancer generated by the TCGA project. Here, only primary tumor (#1,095) and solid tissue normal (#113) samples were retrieved using the TCGAbiolinks R package.ref. bib68 The raw RNA-seq data were pre-processed based on the Array-Array Intensity Correlation (AAIC) method in TCGAbiolinks with default parameters (r > 0.6). Processed data underwent quantile normalization using the default parameters of TCGAanalyze_Filtering function for downstream analyses. All other TCGA datasets used in this work were pre-processed and normalized according to the same methods and parameters applied to this dataset.

TCGA THCA

TCGA THCA is a bulk RNA-seq dataset of thyroid carcinoma generated by the TCGA project. Samples retrieved were primary tumor (#505) and solid tissue normal (#59).

TCGA LUAD

TCGA LUAD is a bulk RNA-seq dataset of lung adenocarcinoma generated by the TCGA project. Samples retrieved were primary tumor (#515) and solid tissue normal (#59).

TCGA LIHC

TCGA LIHC is a bulk RNA-seq dataset of liver hepatocellular carcinoma generated by the TCGA project. Samples retrieved were primary tumor (#371) and solid tissue normal (#50).

TCGA HNSC

TCGA HNSC is a bulk RNA-seq dataset of head & neck squamous cell carcinoma generated by the TCGA project. Samples retrieved were primary tumor (#520) and solid tissue normal (#44).

TCGA STAD

TCGA STAD is a bulk RNA-seq dataset of stomach adenocarcinoma generated by the TCGA project. Samples retrieved were primary tumor (#412) and solid tissue normal (#36).

TCGA LUSC

TCGA LUSC is a bulk RNA-seq dataset of lung squamous cell carcinoma generated by the TCGA project. Samples retrieved were primary tumor (#502) and solid tissue normal (#51).

Schizophrenia dataset

GSE93577 is a total RNA microarray dataset generated from dysfunctional dorsolateral prefrontal cortex layer 3 parvalbumin neurons in 36 matched pairs of schizophrenia and unaffected cases using the Affymetrix Human Genome U219 array.ref. bib34 This dataset contains 141 samples, including 71 healthy and 70 schizophrenia samples. The raw microarray data were retrieved from the GEO database utilizing the GEOquery R packageref. bib69 and log-transformed prior to downstream analyses.

LUAD proteomics dataset

IPX0001804001 is a proteomics dataset generated from primary LUAD samples with paired non-cancerous tumor-adjacent tissues from treatment-naive patients by means of high-performance liquid chromatography-mass spectrometry (HPLC-MS) and label-free quantification.ref. bib36 This dataset contains 206 samples including 103 normal and 103 LUAD samples. The MaxQuant-based pre-processed data was retrieved from the Integrated Proteome Resource.ref. bib70

BRCA proteomics dataset

PXD002057 is a proteomics dataset generated from human breast cancer cell lines SKBR3 and BT474 and their lapatinib-resistant derivative cells by means of nano-scale HPLC-MS and label-free quantification.ref. bib35 This dataset contains 20 samples, including 10 benign and 10 malignant samples. The MaxQuant-based pre-processed data of this dataset was retrieved from the LFQ-Analyst website (https://bioinformatics.erc.monash.edu/apps/LFQ-Analyst/).

Definition of ground truth sets for benchmarking

Ground truth driver genes were obtained from DisGeNET, the Gene Ontology resource, and the MutPanning database, providing experimentally supported annotations. Ground truth non-driver genes were selected from statistically significant features with the least significant adjusted P-values. Ground truth biomarkers were obtained from the NCI Early Detection Research Network and a literature-derived biomarker database, while negative biomarkers were selected from differentially expressed features with minimal fold changes.

Differential expression/abundance analyses in ExIR benchmarking

scRNA-seq datasets

Differential expression analysis (DEA) was performed using the Seurat R package.ref. bib71 scRNA-seq data were prepared and processed within the Seurat framework by creating a Seurat object for each sample via CreateSeuratObject, followed by data normalization using NormalizeData with the “LogNormalize” method and a scale factor of 10,000, where applicable. DEA was conducted via the FindMarkers function to identify differentially expressed genes (DEGs) between cancer and normal cells, employing the default Wilcoxon rank sum test for statistical comparisons. Subsequently, log2-transformed fold change values were used to classify up- and down-regulated genes, with DEGs filtered to retain only those with an adjusted P-value ≤ 0.05.

TCGA bulk RNA-seq datasets

All TCGA datasets used in this work underwent DEA using the TCGAbiolinks R package, which implements functions of edgeR.ref. bib72 Specifically, a common negative binomial dispersion was first estimated across all genes, and a negative binomial log-linear model was fit to the read counts for each gene. Then, pair-wise tests for differential expression between the two groups were performed. All P-values were adjusted, and DEGs with Padj > 0.05 were filtered out.

Schizophrenia dataset

Microarray DEA was performed using the limma R package,ref. bib73 where a linear model was first fit to the expression data of each probe, followed by computing the contrasts of the fitted models with moderated empirical Bayes statistics. DEGs with Padj > 0.1 were filtered out.

Proteomics datasets

Proteomics data were analyzed using the R package DEP.ref. bib74 More precisely, initially each dataset was filtered for proteins that had a maximum of 20 percent missing values in at least one condition within each dataset. Next, the variance of each dataset was normalized followed by a missing value imputation using the “man” algorithm.ref. bib74 Lastly, the differential abundance of proteins was calculated using limma, and P-values were adjusted using the Benjamini-Hochberg procedure. Differentially abundant proteins with Padj > 0.05 were filtered out. In the case of the LUAD proteomics dataset, the differentially abundant proteins with |log₂FC| < 1 were filtered out to maintain the most prominent characteristics of the disease for downstream analyses and benchmarking.

Evaluation of driver prioritization performance in ExIR

To assess the performance of ExIR in driver gene prioritization, sets of curated disease/biological process-associated genes of the above datasets were retrieved from either DisGeNET v7 or the Gene Ontology resourceref. bib50 and considered as the ground truth (Table S5). Additionally, the ground truth drivers of the TCGA datasets were complemented with the driver genes proposed by the MutPanning web server (http://www.cancer-genes.org/).ref. bib75 Next, an intersection analysis was performed to identify common genes between the sets of ground truth driver genes and the sets of significantly up-regulated genes in the selected datasets. As genes with the most statistically significant differential expression are more likely to be driver genes,ref. bib76 sets of significantly differentially expressed genes with the least significant adjusted Padj were selected as true negatives, with each negative set matched in size to its corresponding positive set, thereby ensuring balanced class labels for ROC analysis. Then, these sets of true positive and negative genes were combined as an input into ExIR as the desired lists of genes. Importantly, curated gene sets were used exclusively for benchmarking and did not inform feature prioritization within ExIR, which operates independently of external annotations. The outputs were compared with four commonly used gene prioritization methods; log₂ fold change (log2FC), GeneMANIA,ref. bib38 Endeavour,ref. bib28 and ToppGene.ref. bib29 Owing to the discontinuation and limited availability of the Endeavour platform at the time of manuscript revision, Endeavour could be evaluated for only 8 of the 14 benchmark datasets. The performance of ExIR and comparator methods in prioritizing driver genes was assessed using receiver operating characteristic (ROC) analyses implemented with the plotROC R package.ref. bib77 For each method, all genes in a given dataset were ranked once using the full dataset, without partitioning into training and test sets, as the evaluation focuses on ranking accuracy rather than predictive generalization. Ground-truth driver and non-driver gene sets were used to define true positive and negative labels. The plotROC function systematically varies the ranking threshold to generate ROC curves, enabling quantitative comparison of prioritization performance across methods. To obtain the prioritized driver genes, default parameters were used to run the GeneMANIA and ToppGene models. For Endeavour, the gene ontologies, Reactome pathwaysref. bib78 and STRING PPIsref. bib79 were selected for building the models. Moreover, the training sets required to run the Endeavour and ToppGene models included all genes previously retrieved from the DisGeNET, MutPanning and Gene Ontology databases, except for those genes selected for testing the models (Figures 2A and 2B).

Evaluation of biomarker prioritization performance in ExIR

To evaluate the performance of ExIR in sensitively and specifically identifying biomarkers, lists of biomarkers for both cancer and non-cancer diseases were obtained from the NCI EDRN (https://edrn.nci.nih.gov/biomarkers) and a knowledge-based database of disease-related biomarkersref. bib80 (accessed August 19, 2020), respectively. The EDRN has proposed >20 biomarkers for four out of the nine cancer types corresponding to the above datasets – specifically, lung and breast cancer. The lists of protein/proteomic biomarkers were retrieved from the EDRN (Table S6). Additionally, a list of schizophrenia biomarkers was derived from the database of disease-related biomarkersref. bib80 (Table S6), which employs a knowledge-driven text-mining approach to extract biomarkers for a wide variety of diseases. The common features between these sets and the previously obtained lists of significantly differentially expressed genes/proteins were considered as true positive biomarkers; similarly, the same number of DEGs with the lowest absolute fold changes were selected as true negative biomarkers (Figures 2C and 2D). Importantly, the curated gene sets were used solely for benchmarking purposes and did not influence feature prioritization within ExIR, which functions independently of any external annotations. The combined true positive and negative lists were input into ExIR as the desired lists of genes, and the outputs were compared with four commonly used biomarker prioritization methods: mutual information (MI),ref. bib39 Student’s t-test, the point-biserial correlation coefficient, and the Spearman correlation coefficient. The evaluation and comparison of biomarker prioritization methods were performed based on ROC analyses. MI was calculated between the expression profile of each gene and the binary (0,1) sample labels using the mpmi R package (https://cran.r-project.org/package=mpmi).ref. bib81 Similarly, point-biserial and Spearman correlation coefficients, as well as the Student’s t-test, were computed between gene expression profiles and binary sample labels using the stats R package. Additionally, immunohistochemical data for the specific lung adenocarcinoma (LUAD) subtype analyzed were available from the Human Protein Atlas. We therefore examined protein-level expression of the top five ExIR-prioritized LUAD biomarkers—including both previously established biomarkers from the ground-truth set and ExIR-predicted candidates (considering all genes as well as the subset of upregulated genes)—in unaffected and LUAD tissues using this resourceref. bib41 (http://www.proteinatlas.org). For this analysis, ExIR was applied to the full LUAD dataset without providing predefined true positive or negative gene sets.

Evaluation of mediator identification and prioritization performance in ExIR

As there presently exists no centralized resource containing validated sets of gene mediators of biological processes or diseases, the performance of ExIR in identifying and prioritizing mediators was evaluated based on functional annotation of ExIR outputs. Initially, the performance of ExIR was assessed in comparison to the mediator genes inferred by the MALANI algorithm.ref. bib45 In the context of cancer datasets, MALANI proposes two classes of genes; one being genes frequently differentially expressed or mutated, and another being those genes that are not differentially expressed but may mediate the coordination of oncogenic signals between DE/mutated genes.ref. bib45 A set of breast cancer mediator genes has been proposed based on application of MALANI to the TCGA BRCA dataset; thus, to compare the performance of ExIR against MALANI for mediator detection, the entire TCGA BRCA dataset was input to ExIR without prior provision of any desired gene list. Next, an overrepresentation analysis (ORA) of all ExIR- and MALANI-derived mediators for biological processes and KEGG pathways was performed using the enrichR R package,ref. bib82 and statistically non-significant terms were filtered out (P < 0.05). The association of significant biological processes and KEGG pathways corresponding to ExIR- and MALANI-derived mediators in breast cancer were then separately interrogated using the Comparative Toxicogenomics Database (CTD, accessed October 19, 2021),ref. bib46 a manually-curated repository for literature-based and computationally inferred associations between genes, phenotypes, diseases, etc. Additionally, this benchmarking workflow was then applied beyond breast cancer to all other examined disease datasets.

RNA sequencing

For each experimental replicate, total RNA was extracted from n=3 freshly-dissected 3-month-old zebrafish brains in TRIzol (Invitrogen, 15596). Three independent experimental replicates were used for bulk RNA sequencing, and all samples were assayed for RNA integrity on an Agilent 2100 Bioanalyzer using the Agilent RNA 6000 Nano Kit. 150 bp paired-end sequencing was performed by BGI (Hong Kong) using the DNBseq platform. SOAPnukeref. bib83 was used for adaptor removal and low-quality read filtration, and genome mapping was performed with HISAT2.ref. bib84 Clean reads were mapped to the reference genome using Bowtie2,ref. bib85 and gene-length scaled transcript abundance estimates were calculated with Salmon.ref. bib86 DESeq2ref. bib87 was used for differential gene expression analysis. Downstream visualizations for RNA-seq data were generated in R using EnhancedVolcano,ref. bib88 ExIR (this paper) and clusterProfiler.ref. bib89 DESeq2 differential expression analysis outputs are provided in Table S7, and GO enrichment analyses of RNA datasets are listed in Table S8.

In situ hybridization

Riboprobes for in situ hybridization were generated by cloning transcript-specific PCR products from zebrafish brain cDNA libraries into pGEM-T-Easy (Promega, A1360) using the TA-cloning method. Insert directionality was confirmed by either PCR and/or Sanger sequencing, and plasmids were linearised to facilitate in vitro transcription of antisense riboprobes using either SP6 or T7 RNA polymerase and DIG-RNA labelling mix (Roche, 11277073910). All primers used to clone riboprobe sequences are listed in Table S9.

Adult zebrafish were rapidly euthanized in an ice-water slurry, and exsanguinated on ice by severing major vessels posterior to the anal pore. Whole brains were dissected from the neurocranium in cold 1x phosphate-buffered saline (PBS) pH 7.4, and immediately transferred to 4% paraformaldehyde (PFA, Sigma, 158127) in PBS for overnight fixation at 4°C with gentle rocking. After fixation, brains were cryoprotected in a sucrose-EDTA solution (20% sucrose, 20% 0.5 M EDTA pH 8 in 1x PBS) overnight at 4°C, then cryo-embedded in a mixture of sucrose and fish gelatin as previously described.ref. bib90 Brains were serially cryosectioned at 16 μm thickness using a Leica CS3050S cryostat.

For in situ hybridization, sections were pre-fixed with 4% PFA in PBS pH 7.4 for one hour, then washed twice with 1x PBS with 0.3% Triton X-100 (PBS-Tx 0.3%) for 20 minutes. 10 mg/mL Proteinase K (Roche, 3115879001) was diluted 1:500 in PBS-Tx 0.3%, and sections were digested at RT for five minutes. Sections were then quickly washed with PBS-Tx, and Proteinase K digestion was stopped by incubating sections with 4% PFA at RT for 10 minutes, followed by 2x 10-minute washes with PBS-Tx 0.3%. A hybridization chamber was assembled using a slide box, containing filter paper saturated with a hybridization chamber solution (5 mL 10x Salt solution [1.95 M NaCl, 89 mM Tris-HCl, 11 mM Tris base, 50 mM NaH₂PO₄·2H₂O, 50 mM Na₂HPO₄ and 63.68 mM EDTA], 25 mL formamide and 20 mL ddH₂O, and the hybridization chamber was preheated to 60°C in an incubator. Antisense riboprobes were diluted 1:200 in hybridization buffer (1 mg/mL Torula RNA, 50% formamide, 1x Salt solution, 10% dextran sulfate, 1x Denhardt’s buffer, in ddH₂O), vortexed and denatured at 70°C for 10 minutes prior to addition to sections. Parafilm was used to mitigate probe evaporation during hybridization. Probe hybridization was carried out overnight at 60°C. Unbound/weakly-hybridized probe was then removed by varying stringency through washing with 1x SSC buffer and 50% formamide in ddH₂O 1x 15 minutes, then 2x 30 minutes at 62°C, followed by 2x 30-minute washes at RT in MABT. Sections were blocked in 2% DIG blocking reagent (Roche, 11096176001) in MABT for two hours at RT, then incubated for four hours at RT in AP-conjugated anti-DIG sheep Fab fragments (Roche, 11093274910) diluted 1:2000 in 2% DIG blocking reagent in MABT. Sections were then washed 4x 20 minutes in MABT at RT, and equilibrated in staining buffer (0.1 M NaCl, 0.05 M MgCl₂, 0.1 M Tris pH 9.5, 0.1% Tween-20, all in ddH₂O) for 5 minutes. Chromogenic detection was then performed using NBT-BCIP stock solution (Roche, 11681451001) diluted 1:50 in staining buffer until sufficient signal was observed. Development was terminated by incubating sections in 4% PFA for 30 minutes, then washed 3x 10 minutes in PBS. Sections were mounted with 50% glycerol prior to imaging on a Zeiss Imager.Z2 slide-scanning microscope.

Quantification and statistical analysis

All statistical analyses were performed in R (v4.5.0) as described in the STAR Methods. Differential expression/abundance testing was performed using the Wilcoxon rank-sum test (Seurat), edgeR, limma, or DESeq2, depending on the dataset and experimental modality. Multiple-testing correction was applied where applicable using the Benjamini–Hochberg procedure. Benchmarking of driver, biomarker, and mediator prioritization was evaluated by ROC analyses. All statistical details (test name, test statistic, exact P-values, effect sizes, sample sizes, and number of biological units) are provided in the STAR Methods.

References

1. K. Suhre, G.R. Venkataraman, H. Guturu, A. Halama, N. Stephan, G. Thareja, H. Sarwath, K. Motamedchaboki, M.K.R. Donovan, A. Siddiqui. Nanoparticle enrichment mass-spectrometry proteomics identifies protein-altering variants for precise pQTL mapping. Nat. Commun., 2024. [DOI | PubMed]
2. R. Vitorino. Transforming Clinical Research: The Power of High-Throughput Omics Integration. Proteomes, 2024. [DOI | PubMed]
3. F. Danzi, R. Pacchiana, A. Mafficini, M.T. Scupoli, A. Scarpa, M. Donadelli, A. Fiore. To metabolomics and beyond: a technological portfolio to investigate cancer metabolism. Signal Transduct. Target. Ther., 2023. [DOI | PubMed]
4. A. Manchel, M. Gee, R. Vadigepalli. From sampling to simulating: Single-cell multiomics in systems pathophysiological modeling. iScience, 2024. [DOI | PubMed]
5. A.S. Tang, K.P. Rankin, G. Cerono, S. Miramontes, H. Mills, J. Roger, B. Zeng, C. Nelson, K. Soman, S. Woldemariam. Leveraging electronic health records and knowledge networks for Alzheimer’s disease prediction and sex-specific biological insights. Nat. Aging, 2024. [DOI | PubMed]
6. P. Xiong, H. Chen, J. Zhou, Y. Zeng, Q. Dai. MGRL-DDI: Multiview Graph Representation Learning for Accurate Drug-Drug Interaction Prediction. J. Chem. Inf. Model., 2025. [DOI | PubMed]
7. B. Teji, S. Roy, D. Bhandari, J. Kalita. GT-GRN: a graph transformer framework for enhanced gene regulatory network inference via multimodal embedding of expression data and existing network knowledge. Brief. Bioinform., 2025. [DOI]
8. C. Wang, Z.-P. Liu. Diffusion-based generation of gene regulatory networks from scRNA-seq data with DigNet. Genome Res., 2024. [DOI]
9. Y. Wang, Y. Xu, Z. Yang, X. Liu, Q. Dai. Using Recursive Feature Selection with Random Forest to Improve Protein Structural Class Prediction for Low-Similarity Sequences. Comput. Math. Methods Med., 2021. [DOI]
10. O. Zolotareva, M. Kleine. A Survey of Gene Prioritization Tools for Mendelian and Complex Human Diseases. J. Integr. Bioinform., 2019. [DOI]
11. M.R. Raj, A. Sreeja. Analysis of Computational Gene Prioritization Approaches. Procedia Comput. Sci., 2018. [DOI]
12. E. Seyyedrazzagi, N.J. Navimipour. Disease genes prioritizing mechanisms: a comprehensive and systematic literature review. Network Model. Anal. Health Inf. Bioinf., 2017. [DOI]
13. S.J. Altschuler, L.F. Wu. Cellular heterogeneity: do differences make a difference?. Cell, 2010. [DOI | PubMed]
14. D. Szklarczyk, K. Nastou, M. Koutrouli, R. Kirsch, F. Mehryary, R. Hachilif, D. Hu, M.E. Peluso, Q. Huang, T. Fang. The STRING database in 2025: protein networks with directionality of regulation. Nucleic Acids Res., 2025. [DOI | PubMed]
15. R.M. Piro, F. Di Cunto. Computational approaches to disease-gene prediction: rationale, classification and successes. FEBS J, 2012. [DOI | PubMed]
16. X. Guo, Y. Guo, H. Chen, X. Liu, P. He, W. Li, M.Q. Zhang, Q. Dai. Systematic comparison of genome information processing and boundary recognition tools used for genomic island detection. Comput. Biol. Med., 2023. [DOI | PubMed]
17. R. Kong, X. Xu, X. Liu, P. He, M.Q. Zhang, Q. Dai. 2SigFinder: the combined use of small-scale and large-scale statistical testing for genomic island detection from a single genome. BMC Bioinf., 2020. [DOI]
18. Q. Dai, C. Bao, Y. Hai, S. Ma, T. Zhou, C. Wang, Y. Wang, W. Huo, X. Liu, Y. Yao. MTGIpick allows robust identification of genomic islands from a single genome. Brief. Bioinform., 2016. [DOI]
19. L. Wang, S. Baek, G. Prasad, J. Wildenthal, K. Guo, D. Sturgill, T. Truongvo, E. Char, G. Pegoraro, K. McKinnon. Predictive prioritization of enhancers associated with pancreatic disease risk. Cell Genom., 2026. [DOI | PubMed]
20. T. Włodarczyk, A. Lun, D. Wu, M. Shi, X. Ye, S. Menon, S. Toneyan, K. Seidel, L. Wang, J. Tan. Epiregulon: Single-cell transcription factor activity inference to predict drug response and drivers of cell states. Nat. Commun., 2025. [DOI | PubMed]
21. P.M. Barcelos, I.S. Filgueiras, A.L. Nóbile, J.N. Usuda, A.S. Adri, D.G. de Alburquerque, Y.L.G. Côrrea, F.Y.N. do Vale, I.A.F. Bahia, R.G. Nava. Gene regulatory networks analysis for the discovery of prognostic genes in gliomas. Sci. Rep., 2025. [DOI | PubMed]
22. W.M. Yashar, J. Estabrook, H.D. Holly, J. Somers, O. Nikolova, Ö. Babur, T.P. Braun, E. Demir. Predicting transcription factor activity using prior biological information. iScience, 2024. [DOI | PubMed]
23. T. Wang, Y. Pan, F. Ju, S. Zheng, C. Liu, Y. Min, Q. Jiang, X. Liu, H. Xia, G. Liu. CellNavi predicts genes directing cellular transitions by learning a gene graph-enhanced cell state manifold. Nat. Cell Biol., 2025. [DOI | PubMed]
24. U.D. Akavia, O. Litvin, J. Kim, F. Sanchez-Garcia, D. Kotliar, H.C. Causton, P. Pochanard, E. Mozes, L.A. Garraway, D. Pe’er. An integrated approach to uncover drivers of cancer. Cell, 2010. [DOI | PubMed]
25. D. Rajamani, M.K. Bhasin. Identification of key regulators of pancreatic cancer progression through multidimensional systems-level analysis. Genome Med., 2016. [DOI | PubMed]
26. A.M. Douek, M. Amiri Khabooshan, J. Henry, S.A. Stamatis, F. Kreuder, G. Ramm, M.L. Änkö, D. Wlodkowic, J. Kaslin. An Engineered sgsh Mutant Zebrafish Recapitulates Molecular and Behavioural Pathobiology of Sanfilippo Syndrome A/MPS IIIA. Int. J. Mol. Sci., 2021. [DOI | PubMed]
27. A. Salavaty, M. Ramialison, P.D. Currie. Integrated Value of Influence: An Integrative Method for the Identification of the Most Influential Nodes within Networks. Patterns (N Y), 2020. [DOI]
28. L.C. Tranchevent, A. Ardeshirdavani, S. ElShal, D. Alcaide, J. Aerts, D. Auboeuf, Y. Moreau. Candidate gene prioritization with Endeavour. Nucleic Acids Res., 2016. [DOI | PubMed]
29. J. Chen, E.E. Bardes, B.J. Aronow, A.G. Jegga. ToppGene Suite for gene list enrichment analysis and candidate gene prioritization. Nucleic Acids Res., 2009. [DOI | PubMed]
30. Z.K. Tuong, K.W. Loudon, B. Berry, N. Richoz, J. Jones, X. Tan, Q. Nguyen, A. George, S. Hori, S. Field. Resolving the immune landscape of human prostate at a single-cell level in health and cancer. Cell Rep., 2021. [DOI | PubMed]
31. S.V. Puram, I. Tirosh, A.S. Parikh, A.P. Patel, K. Yizhak, S. Gillespie, C. Rodman, C.L. Luo, E.A. Mroz, K.S. Emerick. Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Head and Neck Cancer. Cell, 2017. [DOI | PubMed]
32. J.N. Weinstein, E.A. Collisson, G.B. Mills, K.R.M. Shaw, B.A. Ozenberger, K. Ellrott, I. Shmulevich, C. Sander, J.M. Stuart, J.M. Stuart. The Cancer Genome Atlas Pan-Cancer analysis project. Nat. Genet., 2013. [DOI | PubMed]
33. S. Darmanis, S.A. Sloan, D. Croote, M. Mignardi, S. Chernikova, P. Samghababi, Y. Zhang, N. Neff, M. Kowarsky, C. Caneda. Single-Cell RNA-Seq Analysis of Infiltrating Neoplastic Cells at the Migrating Front of Human Glioblastoma. Cell Rep., 2017. [DOI | PubMed]
34. J.F. Enwright Iii, Z. Huo, D. Arion, J.P. Corradi, G. Tseng, D.A. Lewis. Transcriptome alterations of prefrontal cortical parvalbumin neurons in schizophrenia. Mol. Psychiatry, 2018. [DOI | PubMed]
35. H. Creedon, L. Gómez-Cuadrado, Ž. Tarnauskaitė, J. Balla, M. Canel, K.G. MacLeod, B. Serrels, C. Fraser, A. Unciti-Broceta, N. Tracey. Identification of novel pathways linking epithelial-to-mesenchymal transition with resistance to HER2-targeted therapy. Oncotarget, 2016. [DOI | PubMed]
36. J.Y. Xu, C. Zhang, X. Wang, L. Zhai, Y. Ma, Y. Mao, K. Qian, C. Sun, Z. Liu, S. Jiang. Integrative Proteomic Characterization of Human Lung Adenocarcinoma. Cell, 2020. [DOI | PubMed]
37. Z.H. Zhao, J.Y. Ma, T.G. Meng, Z.B. Wang, W. Yue, Q. Zhou, S. Li, X. Feng, Y. Hou, H. Schatten. Single-cell RNA sequencing reveals the landscape of early female germ cell development. FASEB J, 2020. [DOI | PubMed]
38. D. Warde-Farley, S.L. Donaldson, O. Comes, K. Zuberi, R. Badrawi, P. Chao, M. Franz, C. Grouios, F. Kazi, C.T. Lopes. The GeneMANIA prediction server: biological network integration for gene prioritization and predicting gene function. Nucleic Acids Res., 2010. [DOI | PubMed]
39. C.E. Shannon. A mathematical theory of communication. Bell Syst. Tech. J., 1948
40. Y. Liu, A. Beyer, R. Aebersold. On the Dependency of Cellular Protein Levels on mRNA Abundance. Cell, 2016. [DOI | PubMed]
41. M. Uhlen, C. Zhang, S. Lee, E. Sjöstedt, L. Fagerberg, G. Bidkhori, R. Benfeitas, M. Arif, Z. Liu, F. Edfors. A pathology atlas of the human cancer transcriptome. Science, 2017. [DOI]
42. D. Shao, Y. Liu, X. Liu, L. Zhu, Y. Cui, A. Cui, A. Qiao, X. Kong, Y. Liu, Q. Chen. PGC-1β-Regulated mitochondrial biogenesis and function in myotubes is mediated by NRF-1 and ERRα. Mitochondrion, 2010. [DOI | PubMed]
43. A. Disanza, S. Mantoani, M. Hertzog, S. Gerboth, E. Frittoli, A. Steffen, K. Berhoerster, H.-J. Kreienkamp, F. Milanesi, P.P.D. Fiore. Regulation of cell shape by Cdc42 is mediated by the synergic actin-bundling activity of the Eps8–IRSp53 complex. Nat. Cell Biol., 2006. [PubMed]
44. K.M. Tyssowski, N.R. DeStefino, J.H. Cho, C.J. Dunn, R.G. Poston, C.E. Carty, R.D. Jones, S.M. Chang, P. Romeo, M.K. Wurzelmann. Different Neuronal Activity Patterns Induce Different Gene Expression Programs. Neuron, 2018. [DOI | PubMed]
45. M. Ghanat Bari, C.Y. Ung, C. Zhang, S. Zhu, H. Li. Machine learning-assisted network inference approach to identify a new class of genes that coordinate the functionality of cancer networks. Sci. Rep., 2017. [PubMed]
46. A.P. Davis, C.J. Grondin, R.J. Johnson, D. Sciaky, J. Wiegers, T.C. Wiegers, C.J. Mattingly. Comparative Toxicogenomics Database (CTD): update 2021. Nucleic Acids Res., 2021. [DOI | PubMed]
47. P. Gaudet, C. Dessimoz. The gene ontology handbook, 2017
48. S. Yon Rhee, V. Wood, K. Dolinski, S. Draghici. Use and misuse of the gene ontology annotations. Nat. Rev. Genet., 2008. [PubMed]
49. A.M. Altenhoff, R.A. Studer, M. Robinson-Rechavi, C. Dessimoz. Resolving the Ortholog Conjecture: Orthologs Tend to Be Weakly, but Significantly, More Similar in Function than Paralogs. PLoS Comput. Biol., 2012
50. M. Ashburner, C.A. Ball, J.A. Blake, D. Botstein, H. Butler, J.M. Cherry, A.P. Davis, K. Dolinski, S.S. Dwight, J.T. Eppig. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat. Genet., 2000. [DOI | PubMed]
51. A. Arfi, M. Richard, C. Gandolphe, D. Bonnefont-Rousselot, P. Thérond, D. Scherman. Neuroinflammatory and oxidative stress phenomena in MPS IIIA mouse model: the positive effect of long-term aspirin treatment. Mol. Genet. Metab., 2011. [DOI | PubMed]
52. L.D. Archer, K.J. Langford-Smith, B.W. Bigger, J.E. Fildes. Mucopolysaccharide diseases: a complex interplay between neuroinflammation, microglial activation and adaptive immunity. J. Inherit. Metab. Dis., 2014. [DOI | PubMed]
53. C. Lee, T.E. Dineen, M. Brack, J.E. Kirsch, V.M. Runge. The mucopolysaccharidoses: characterization by cranial MR imaging. AJNR Am. J. Neuroradiol., 1993. [PubMed]
54. C. Pará, P. Bose, L. Bruno, E. Freemantle, M. Taherzadeh, X. Pan, C. Han, P.S. McPherson, J.C. Lacaille, É. Bonneil. Early defects in mucopolysaccharidosis type IIIC disrupt excitatory synaptic transmission. JCI Insight, 2021. [DOI]
55. I. Sambri, R. D’Alessio, Y. Ezhova, T. Giuliano, N.C. Sorrentino, V. Cacace, M. De Risi, M. Cataldi, L. Annunziato, E. De Leonibus, A. Fraldi. Lysosomal dysfunction disrupts presynaptic maintenance and restoration of presynaptic function prevents neurodegeneration in lysosomal storage diseases. EMBO Mol. Med., 2017. [DOI | PubMed]
56. C.A. Dwyer, S.L. Scudder, Y. Lin, L.E. Dozier, D. Phan, N.J. Allen, G.N. Patrick, J.D. Esko. Neurodevelopmental Changes in Excitatory Synaptic Structure and Function in the Cerebral Cortex of Sanfilippo Syndrome IIIA Mice. Sci. Rep., 2017. [DOI]
57. M. Westerfield. The zebrafish book: a guide for the laboratory use of zebrafish. 2000
58. L. Breiman. Random forests. Mach. Learn., 2001. [DOI]
59. S. Wenric, R. Shemirani. Using Supervised Learning Methods for Gene Selection in RNA-Seq Case-Control Studies. Front. Genet., 2018. [DOI | PubMed]
60. M.N. Wright, A. Ziegler. ranger : A Fast Implementation of Random Forests for High Dimensional Data in C++ and R. J. Stat. Softw., 2017. [DOI]
61. S. Nembrini, I.R. König, M.N. Wright. The revival of the Gini importance?. Bioinformatics, 2018. [DOI | PubMed]
62. A. Altmann, L. Toloşi, O. Sander, T. Lengauer. Permutation importance: a corrected feature importance measure. Bioinformatics, 2010. [DOI | PubMed]
63. D. Nitsch, L.C. Tranchevent, B. Thienpont, L. Thorrez, H. Van Esch, K. Devriendt, Y. Moreau. Network analysis of differential expression for the identification of disease-causing genes. PLoS One, 2009. [DOI]
64. C. Cheadle, M.P. Vawter, W.J. Freed, K.G. Becker. Analysis of microarray data using Z score transformation. J. Mol. Diagn., 2003. [DOI | PubMed]
65. Y. Zhang, S.A. Sloan, L.E. Clarke, C. Caneda, C.A. Plaza, P.D. Blumenthal, H. Vogel, G.K. Steinberg, M.B. Edwards, G. Li. Purification and Characterization of Progenitor and Mature Human Astrocytes Reveals Transcriptional and Functional Differences with Mouse. Neuron, 2016. [DOI | PubMed]
66. S. Anders, P.T. Pyl, W. Huber. HTSeq–a Python framework to work with high-throughput sequencing data. Bioinformatics, 2015. [DOI | PubMed]
67. T. Stuart, A. Butler, P. Hoffman, C. Hafemeister, E. Papalexi, W.M. Mauck, Y. Hao, M. Stoeckius, P. Smibert, R. Satija. Comprehensive Integration of Single-Cell Data. Cell, 2019. [DOI | PubMed]
68. A. Colaprico, T.C. Silva, C. Olsen, L. Garofano, C. Cava, D. Garolini, T.S. Sabedot, T.M. Malta, S.M. Pagnotta, I. Castiglioni. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res., 2016. [DOI]
69. S. Davis, P.S. Meltzer. GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics, 2007. [PubMed]
70. J. Ma, T. Chen, S. Wu, C. Yang, M. Bai, K. Shu, K. Li, G. Zhang, Z. Jin, F. He. iProX: an integrated proteome resource. Nucleic Acids Res., 2019. [DOI | PubMed]
71. Y. Hao, T. Stuart, M.H. Kowalski, S. Choudhary, P. Hoffman, A. Hartman, A. Srivastava, G. Molla, S. Madad, C. Fernandez-Granda, R. Satija. Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nat. Biotechnol., 2024. [DOI | PubMed]
72. M.D. Robinson, D.J. McCarthy, G.K. Smyth. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics, 2010. [DOI | PubMed]
73. M.E. Ritchie, B. Phipson, D. Wu, Y. Hu, C.W. Law, W. Shi, G.K. Smyth. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res., 2015. [DOI | PubMed]
74. X. Zhang, A.H. Smits, G.B. van Tilburg, H. Ovaa, W. Huber, M. Vermeulen. Proteome-wide identification of ubiquitin interactions using UbIA-MS. Nat. Protoc., 2018. [DOI | PubMed]
75. F. Dietlein, D. Weghorn, A. Taylor-Weiner, A. Richters, B. Reardon, D. Liu, E.S. Lander, E.M. Van Allen, S.R. Sunyaev. Identification of cancer driver genes based on nucleotide context. Nat. Genet., 2020. [DOI | PubMed]
76. C.J. Tokheim, N. Papadopoulos, K.W. Kinzler, B. Vogelstein, R. Karchin. Evaluating the evaluation of cancer driver genes. Proc. Natl. Acad. Sci. USA, 2016. [DOI | PubMed]
77. M.C. Sachs. plotROC: A Tool for Plotting ROC Curves. J. Stat. Softw., 2017. [DOI]
78. B. Jassal, L. Matthews, G. Viteri, C. Gong, P. Lorente, A. Fabregat, K. Sidiropoulos, J. Cook, M. Gillespie, R. Haw. The reactome pathway knowledgebase. Nucleic Acids Res., 2019. [DOI]
79. D. Szklarczyk, A.L. Gable, D. Lyon, A. Junge, S. Wyder, J. Huerta-Cepas, M. Simonovic, N.T. Doncheva, J.H. Morris, P. Bork. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res., 2019. [DOI | PubMed]
80. À. Bravo, M. Cases, N. Queralt-Rosinach, F. Sanz, L.I. Furlong. A knowledge-driven approach to extract disease-related biomarkers from the literature. Biomed Res. Int., 2014. [DOI]
81. C. Pardy, S. Galbraith, S.R. Wilson. Integrative Exploration of Large High-Dimensional Datasets. Ann. Appl. Stat., 2018. [DOI]
82. M.V. Kuleshov, M.R. Jones, A.D. Rouillard, N.F. Fernandez, Q. Duan, Z. Wang, S. Koplev, S.L. Jenkins, K.M. Jagodnik, A. Lachmann. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res., 2016. [DOI | PubMed]
83. Y. Chen, Y. Chen, C. Shi, Z. Huang, Y. Zhang, S. Li, Y. Li, J. Ye, C. Yu, Z. Li. SOAPnuke: a MapReduce acceleration-supported software for integrated quality control and preprocessing of high-throughput sequencing data. Gigascience, 2018. [DOI]
84. D. Kim, J.M. Paggi, C. Park, C. Bennett, S.L. Salzberg. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol., 2019. [DOI | PubMed]
85. B. Langmead, S.L. Salzberg. Fast gapped-read alignment with Bowtie 2. Nat. Methods, 2012. [DOI | PubMed]
86. R. Patro, G. Duggal, M.I. Love, R.A. Irizarry, C. Kingsford. Salmon provides fast and bias-aware quantification of transcript expression. Nat. Methods, 2017. [DOI | PubMed]
87. M.I. Love, W. Huber, S. Anders. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol., 2014. [DOI | PubMed]
88. K. Blighe, S. Rana, M. Lewis. EnhancedVolcano: Publication-ready volcano plots with enhanced colouring and labeling. R package version, 2019
89. S. Xu, E. Hu, Y. Cai, Z. Xie, X. Luo, L. Zhan, W. Tang, Q. Wang, B. Liu, R. Wang. Using clusterProfiler to characterize multiomics data. Nat. Protoc., 2024. [DOI | PubMed]
90. B.W. Lindsey, G.E. Aitken, J.K. Tang, M. Khabooshan, A.M. Douek, C. Vandestadt, J. Kaslin. Midbrain tectal stem cells display diverse regenerative capacities in zebrafish. Sci. Rep., 2019. [DOI | PubMed]