
==== Front
Plant Commun
Plant Commun
Plant Communications
2590-3462
Elsevier

S2590-3462(24)00319-5
10.1016/j.xplc.2024.101002
101002
Resource Article
Dual-extraction modeling: A multi-modal deep-learning architecture for phenotypic prediction and functional gene mining of complex traits
Ren Yanlin 14
Wu Chenhua 14
Zhou He 14
Hu Xiaona huxiaona109@nwafu.edu.cn
2∗
Miao Zhenyan miaozhenyan@nwafu.edu.cn
13∗∗
1 State Key Laboratory for Crop Stress Resistance and High-Efficiency Production, Center of Bioinformatics, College of Life Sciences, Northwest A&F University, Yangling, Shaanxi 712100, China
2 College of Chemistry & Pharmacy, Northwest A&F University, Yangling, Shaanxi 712100, China
3 Key Laboratory of Biology and Genetics Improvement of Maize in Arid Area of Northwest Region, Ministry of Agriculture, Northwest A&F University, Yangling, Shaanxi 712100, China
∗ Corresponding author huxiaona109@nwafu.edu.cn
∗∗ Corresponding author miaozhenyan@nwafu.edu.cn
4 These authors have contributed equally to this article.

13 6 2024
09 9 2024
13 6 2024
5 9 10100225 1 2024
27 5 2024
11 6 2024
© 2024 The Author(s)
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Despite considerable advances in extracting crucial insights from bio-omics data to unravel the intricate mechanisms underlying complex traits, the absence of a universal multi-modal computational tool with robust interpretability for accurate phenotype prediction and identification of trait-associated genes remains a challenge. This study introduces the dual-extraction modeling (DEM) approach, a multi-modal deep-learning architecture designed to extract representative features from heterogeneous omics datasets, enabling the prediction of complex trait phenotypes. Through comprehensive benchmarking experiments, we demonstrate the efficacy of DEM in classification and regression prediction of complex traits. DEM consistently exhibits superior accuracy, robustness, generalizability, and flexibility. Notably, we establish its effectiveness in predicting pleiotropic genes that influence both flowering time and rosette leaf number, underscoring its commendable interpretability. In addition, we have developed user-friendly software to facilitate seamless utilization of DEM’s functions. In summary, this study presents a state-of-the-art approach with the ability to effectively predict qualitative and quantitative traits and identify functional genes, confirming its potential as a valuable tool for exploring the genetic basis of complex traits.

This study introduces the dual-extraction modeling (DEM) architecture for prediction of phenotypes and identification of functional genes through analysis of multi-modal omics data. DEM combines new methods, high performance, interpretability, and easy-to-use software, making it a valuable tool for understanding complex traits in the postgenomic era.

Key words

deep learning
complex traits
multi-omics
phenotypic prediction
gene mining
Published: June 13, 2024
==== Body
pmcIntroduction

Phenotypic variations in complex traits are intricately linked to genetic and epigenetic variations in multiple molecular processes. These processes involve single-nucleotide polymorphisms (SNPs), DNA methylation, RNA transcription, protein expression, and more (Gallusci et al., 2017; Springer et al., 2019). Leveraging the complementary information embedded in diverse molecular layers can significantly enhance the accuracy of phenotypic prediction and expedite the identification of pertinent biomarkers. In the postgenomic era, high-throughput sequencing technologies have generated copious multi-omics data, encompassing genomic, epigenomic, transcriptomic, and proteomic datasets. Despite offering comprehensive insights into the intricate relationships among molecular layers, the integration and translation of multi-modal omics data into meaningful biological insights remain significant challenges in complex trait research.

Machine learning (ML), a subset of artificial intelligence technology, has demonstrated substantial capacity, efficiency, and flexibility in encoding and learning hidden biological knowledge from intricate multi-omics data, addressing phenotypic prediction challenges (Xu and Jackson, 2019; Yan and Wang, 2023). Some traditional ML approaches such as random forest (RF), support vector machine, and LightGBM have enabled improvements in prediction accuracy across various crop traits (Long et al., 2011; Holliday et al., 2012; Yan et al., 2021). However, these algorithms require complex data preprocessing before modeling, and an abundance of noisy features can significantly diminish their predictive performance.

Deep learning (DL), as an ML paradigm, has demonstrated enhanced predictive performance in supervised learning tasks, providing an effective means for exploring the relationship between omics data and phenotypic traits (Ma et al., 2018; Liu et al., 2020; Kang et al., 2022; Wang et al., 2023). For instance, the DeepGS approach, based on a convolutional neural network (CNN) model, accurately predicted eight wheat phenotypes using 33 709 wheat genomic markers (Ma et al., 2018). Li and collaborators developed a CNN multi-omics prediction approach, DNNGP, which uses multi-omics data to achieve accurate phenotypic predictions across different species (Wang et al., 2023). However, owing to the limited receptive fields of CNN models, these approaches have struggled to build relationships between omics features and phenotypes from a global perspective. In some human complex-trait classification tasks, phenotypic prediction models using multi-omics data, such as CustOmics (Benkirane et al., 2023) and MOGAT (Tanvir et al., 2024), outperformed modeling based on single-omics data alone. Nevertheless, these approaches relied on the data distribution of input data and suffered from over-smoothing issues, thereby limiting the depth and generalizability of the network.

The self-attention mechanism, stacked with multiple encoder modules, updates sample feature weights globally, facilitating the learning of dependencies between distant nodes and yielding greater prediction accuracy. Although the self-attention mechanism holds promise for translating multi-omics data features into meaningful biological insights for human disease classification (Jiang et al., 2022; Kaczmarek et al., 2022), challenges persist in terms of interpretability and generalizability when implementing it for the analysis of multi-modal omics data across diverse plant species. In addition, no computational tools have demonstrated superior performance in both classification and regression tasks.

In this study, we present a dual-extraction modeling (DEM) architecture that incorporates a multi-head self-attention mechanism and a fully connected feedforward neural network (FFN) for phenotypic prediction and functional gene mining of complex traits, achieving superior performance and interpretability. Systematic benchmarks using multi-modal omics datasets of qualitative and quantitative traits demonstrate the effectiveness of our method in various challenging tasks across different species.

Results

Overview of DEM architecture

DEM functions as a multi-modal supervised DL method designed for predicting phenotypes and pinpointing functional genes associated with complex traits from heterogeneous multi-omics datasets. The architecture comprises four pivotal functional components (Figure 1).(1) Data preprocessing. Initially, a neural network is constructed to transform sparse and discrete SNP features into dense and continuous values (Supplemental Figure 1). Subsequently, raw omics features with a high proportion of nulls and low variances are removed. The filtered features are then imputed and scaled. For classification tasks, omics features that demonstrate distinguishable differences among phenotypic classes are selected for subsequent DEM. For regression tasks, phenotypic values are scaled, and the RF algorithm is used to select informative feature candidates for subsequent DEM.

(2) Dual-extraction modeling. The modeling process unfolds in two stages. Initially, individual single-omics feature matrices, along with their corresponding phenotypic labels, serve as independent inputs to the multi-head self-attention network, which is followed by an FFN. This supervised learning process is executed separately for each omics in the high-dimensional feature space. Concurrently, all single-omics features are concatenated into a multi-omics joint matrix, serving as an independent input for training the multi-omics joint model. In the second stage, latent spatial information from each independent modeling process in the high-dimensional feature space is combined as input to the multi-head self-attention network, aiming to filter out noise and transfer omics features to a low-dimensional space, improving separability and representation. Upon minimizing the total losses of phenotypic prediction in the system, the optimal DEM model is established for subsequent phenotypic prediction and functional gene mining.

(3) Phenotypic prediction. The phenotypic prediction results are derived from the dot product operation between the weights of the DEM model and preprocessed omics features.

(4) Functional gene mining. A post hoc interpretation method is devised for estimating feature importance scores (refer to methods). It compares the phenotypic prediction performance of the DEM model using actual feature values against that based on randomly shuffled feature values. High-ranking features are subsequently screened out and accumulated with their corresponding genes to represent trait-related genes.

Figure 1 The architecture of DEM.

DEM is a multi-modal deep-learning method that integrates the multi-head self-attention mechanism and fully connected feedforward neural networks. It consists of four essential functional components: data preprocessing, feature extraction modeling, phenotype prediction, and functional gene mining. The primary objective of DEM is to extract representative features from multi-omics data, which are then used to predict phenotypes associated with complex traits and to mine functional genes. Versatile in nature, DEM is capable of accommodating both classification and regression tasks.

The end-to-end DEM architecture excels in harnessing latent space information derived from individual omics models, thereby complementing omics-specific details absent from the hidden layer of the multi-omics joint model. This capability facilitates the extraction of multi-omics features, implying both cross-omics correlations and omics-specific knowledge. Such capabilities prove advantageous for interpreting phenotypic predictions and mining functional genes, as demonstrated in subsequent benchmarking experiments.

Evaluating the effectiveness of DEM and its variations through ablation experiments

To assess the impact of the DL algorithms within the DEM architecture on both regression and classification tasks, we conducted extensive ablation studies using two datasets: the flowering time (FT) dataset of Arabidopsis from the 1001 Genomes Project for regression tasks and the breast invasive carcinoma (BRCA) dataset obtained from The Cancer Genome Atlas (TCGA) Program for classification tasks.

Initially, we used five DL algorithms (CNN, Graph Convolutional Network [GCN], autoencoder [AE], long short-term memory [LSTM], and multi-layer perceptron [MLP]) to replace the multi-head self-attention mechanism in both the high-dimensional and low-dimensional feature spaces of the DEM architecture. Hyperparameters for each DL algorithm underwent fine-tuning, and the two best-performing sets were selected for each DL algorithm to ensure a fair comparison (Supplemental Figures 2 and 3). We then compared the performance of our proposed DEM model with that of the other variations. As depicted in Figure 2A–2C, DEM exhibited a higher Pearson correlation coefficient (PCC) and the lowest mean absolute error (MAE) and mean squared error (MSE) compared with the other variant models in terms of the phenotypic regression task. In addition, we evaluated the performance of the classification task on the BRCA dataset using the mean accuracy (ACC), weighted F1 score (F1-weighted), and macro F1 score (F1-macro). The results indicated that the DEM still achieved the highest mean ACC, F1-weighted, and F1-macro (Supplemental Figure 4).Figure 2 Evaluating the performance of DEM through ablation experiments.

(A–C) Performance comparisons between the proposed DEM model and other variant models, evaluating PCC (A), MSE (B), and MAE (C) on the Arabidopsis FT datasets. Five deep-learning algorithms (CNN, GCN, AE, LSTM, and MLP) were used to replace the multi-head self-attention mechanisms in both the high-dimensional and low-dimensional feature spaces of the DEM architecture to obtain the DEM variations. The values in the grids represent the performance metrics (PCC, MSE, and MAE) of DEM and the other variant models.

(D and E) Performance comparisons of DEM with two SEM models on the FT and BRCA datasets. The performance was quantified by mean PCC, MSE, and MAE (D) and ACC, F1-weighted, and F1-macro (E) values (n = 10 repeats with the outer loop of nested cross-validation).

To scrutinize the necessity of the dual-extraction process in the DEM architecture, we compared the performance of DEM with two single-extraction modeling (SEM) architectures (Supplemental Figure 5). As illustrated in Figure 2D and 2E, DEM outperformed the two SEM architectures in both regression and classification tasks, highlighting the enhanced capability of DEM in extracting representative features for phenotypic prediction.

DEM exhibits superior performance in phenotypic prediction

We began by evaluating the performance of the DEM architecture in quantitative plant trait regression tasks. Considering the continuous nature of phenotypic variations in quantitative plant traits influenced by genetic variations and gene expression regulation across multiple molecular layers, we used multi-omics data, including population genome, transcriptome, and DNA methylome data, for two Arabidopsis phenotypic datasets (FT and rosette leaf number [RLN]) and multi-omics data with four phenotypic datasets (days to tasseling [DTT], plant height [PH], kernel number per ear [KNPE], and kernel weight per ear [KWPE]) for maize.

Initially, we transformed raw sparse and discrete SNP features into dense and continuous features representing genomic variation features (Supplemental Figure 1) and preprocessed each omics dataset to filter out noisy features. After assessing the impact of feature number on predictive performance, we selected the most effective feature set with an optimal number of features from each omics dataset for phenotypic prediction modeling (Supplemental Figures 6 and 7). Three hyperparameters (number of encoders, dropout probability, and learning rate) were fine-tuned using nested cross-validation to achieve optimal quantitative trait regression prediction performance (Supplemental Figure 8). We then compared the performance of DEM with that of two published phenotypic prediction methods (DNNGP and Bayesian ridge regression [BRR] [Perez et al., 2010]) using PCC, MAE, and MSE metrics. DEM consistently yielded the highest mean PCC and the lowest mean MAE and MSE for both the Arabidopsis and maize datasets (Figure 3A–3F), demonstrating its superiority in regression prediction tasks.Figure 3 Superior performance of DEM models in phenotypic prediction.

(A–F) Comparative analysis of regression performance using three models (BRR, DNNGP, and DEM) on the Arabidopsis FT (A) and RLN (B) datasets, as well as four maize datasets (KWPE [C], KNPE [D], DTT [E], and PH [F]). Performance metrics (PCC, MAE, and MSE) were calculated as the average of 10 repetitions with the outer loop of nested cross-validation.

(G and H) Comparative analysis of classification performance using four models (CustOmics, AttentionMOI, MOGAT, and DEM) on the two rice datasets BR (G) and SS (H). Performance metrics (ACC, F1-weighted, and F1-macro) were calculated as the average of 10 repetitions with the outer loop of nested cross-validation.

For qualitative plant trait classification tasks, we used two multi-classification datasets (straight head susceptibility [SS] and blast resistance [BR]) from rice (Zhao et al., 2011). After determining optimal model configurations for each dataset through hyperparameter fine-tuning (Supplemental Figure 9), we compared the performance of DEM with that of three published classification prediction methods: AttentionMOI (Pang et al., 2023), CustOmics (Benkirane et al., 2023), and MOGAT (Tanvir et al., 2024). DEM exhibited the best performance in multi-classification tasks, as evaluated by ACC, F1-weighted, and F1-macro metrics (Figure 3G and 3H). Considering that the other three methods have been reported to perform well on tasks related to human disease classification, we further evaluated the performance of DEM using human disease datasets (BRCA and brain lower grade glioma [LGG]), on which it demonstrated superior performance in human disease classification prediction (Supplemental Figure 10A and 10B).

We also performed independent external validation experiments to assess the generalizability of DEM. For regression tasks, FT and RLN datasets from Arabidopsis grown in Sweden and DTT, PH, KNPE, and KWPE datasets from maize grown in Beijing, China, were used as external validation datasets. For classification tasks, SS and BR datasets from rice grown in the United States, and BRCA and LGG datasets generated from formalin-fixed paraffin-embedded tissues were used as external validation datasets. Compared with five other methods (DNNGP, BRR, AttentionMOI, CustOmics, and MOGAT), DEM consistently achieved the best phenotypic prediction performance on these independent validation datasets (Supplemental Figures 10C, 10D, and 11). These results underscore the superiority of DEM over other methods for both qualitative and quantitative trait prediction tasks.

DEM can identify trait-related functional genes

Following the post hoc interpretation approach introduced in the methods section, we quantified and ranked the feature importance levels for phenotypic prediction in Arabidopsis, maize, rice, and human disease (Supplemental Table 1). On the basis of the rate of change in feature importance scores, 137, 126, 101, 173, 151, 166, 181, 118, 133, and 169 gene-level features were identified as trait-related candidates from the FT, RLN, DTT, PH, KNPE, KWPE, SS, BR, BRCA, and LGG datasets, respectively (Figure 4A–4H; Supplemental Figure 12A and 12B). Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed to interpret the biological functions of these functional genes, and the enriched GO terms and KEGG pathways demonstrated strong associations with the corresponding traits (Figure 4I and Supplemental Figure 12C).Figure 4 Identification of trait-related features through feature importance analysis.

(A–H) Omics features identified as trait-related candidates from datasets of eight plant traits. Features were ranked on the basis of their feature importance scores, and those with mean rates of change >0.05 were selected as trait-related candidates.

(I) Significantly enriched GO terms and KEGG pathways associated with trait-related candidates in datasets of eight traits. Significance was set to a Q value of <0.05 for enriched terms.

Among the top 30 important candidates from each dataset, some genes have been reported to play roles in the corresponding traits (Supplemental Tables 2–11). For example, FLOWERING LOCUS C (FLC) encodes a key regulatory protein that acts as a dose-dependent repressor of Arabidopsis flowering (Michaels and Amasino, 1999). SUPPRESSOR OF OVEREXPRESSION OF CONSTANS 1 (SOC1) encodes a multifunctional protein that regulates not only FT but also floral patterning and floral meristem determinacy (Melzer et al., 2008). Suppression of NADPH dehydrogenase B1 (NDB1) in Arabidopsis affects central metabolism and rosette growth (Wallstrom et al., 2014). Natriuretic peptide receptor 1 (NPR1) is a key positive factor in disease resistance and plays an essential role in resistance to blast fungus disease (Sugano et al., 2010). SQUAMOSA promoter-binding protein-like 4 (SPL4) positively regulates resistance to blast disease, enhancing immune responses in rice and delaying infection by Magnaporthe oryzae (Zhang et al., 2022).

Together, these findings highlight the strong generalizability and interpretability of DEM for phenotypic predictions of qualitative and quantitative traits involving different species.

Flexibility of DEM in the context of omics data variations

Capitalizing on the self-attention mechanism, which adeptly assigns weights to the most relevant features in a flexible manner, the DEM architecture exhibits the potential to handle noise information and capture essential features in diverse types of omics data. To assess the robustness of DEM in the presence of low-quality omics data, Gaussian noise conforming to the standard normal distribution was introduced into each training set. The results from the phenotypic regression tasks revealed that the performance of DEM on noisy data remained comparable to that on raw data (Figure 5A–5F and Supplemental Figure 13A–13F). The introduction of noise into the Arabidopsis datasets resulted in a reduction of about 4% in DEM performance. Even in the maize datasets, the impact of noise injection did not exceed 5%. Similar robustness was observed for phenotypic classification tasks, as the performance of DEM was not significantly affected by noise injection into the rice and human disease datasets (Figure 5G and 5H; Supplemental Figures 13G, 13H, 14A, and 14B).Figure 5 Stability and flexibility of DEM in the context of omics data variations.

(A–H) Evaluation of the stability of DEM performance on the two Arabidopsis datasets FT (A) and RLN (B), the four maize datasets DTT (C), PH (D), KNPE (E), and KWPE (F), and the two rice datasets SS (G) and BR (H). Performance was measured in terms of PCC (for regression tasks) and ACC (for classification tasks). The reduction ratio (RR) of the performance was calculated as the mean RR of the metric values resulting from random addition of Gaussian noise 10 times during modeling.

(I–P) Comparative analysis of predictive performance among six methods (DEM, DNNGP, BRR, AttentionMOI, CustOmics, and MOGAT) across varying noise ratios. Performance was measured in terms of PCC and ACC. The bands represent mean metric values ± SD (n = 10 repeats with the outer loop of nested cross-validation).

(Q and R) Comparative analysis of DEM performance across diverse combinations of omics modalities in regression tasks, evaluated using the PCC. Boxplots depict the distribution of PCC values, with centerlines representing the median values and triangles indicating the mean values. The experiments were performed with 10 repeats using the outer loop of nested cross-validation.

(S and T) Comparative analysis of regression performance between DEM and two alternative methods (DNNGP and BRR) across variations in omics modalities. The bands illustrate the mean PCC values along with the standard deviation for 10 repeats using the outer loop of nested cross-validation.

We next compared the performance of DEM, trained on raw data and tested on noisy data, with that of five published phenotypic prediction methods (DNNGP and BRR for regression tasks; AttentionMOI, CustOmics, and MOGAT for classification tasks). The results indicated that although the performance of DEM was influenced by noise injection into the testing data, it remained more robust to noise, consistently outperforming the other five methods under different noise ratios (Figure 5I–5P; Supplemental Figures 14C, 14D, and 15).

In evaluating the ability of DEM to handle different omics data types, we compared its performance when trained using combinations of various omics modalities (single-omics modeling, double-omics modeling, and triple-omics modeling). As illustrated in Figure 5Q and 5R and Supplemental Figure 16, in regression tasks, with an increase in the number of omics types, the triple-omics model achieved the best performance, and the performance of double-omics models was notably better than that of single-omics models. In addition, in classification tasks, DEM trained with multi-omics data types outperformed DEM trained with single-omics data types (Supplemental Figure 17).

Next, we compared the performance of DEM with that of five published phenotypic prediction methods (AttentionMOI, CustOmics, MOGAT, DNNGP, and BRR) concerning different omics modalities. The results consistently demonstrated the superiority of DEM relative to the other five methods (Figure 5S and 5T; Supplemental Figure 18) across all combinations of omics data types, underscoring its advanced capability for omics-specific and integrative multi-omics learning.

DEM enhances the prediction of pleiotropic genes

Pleiotropy, often referred to as a genetic tradeoff, is a prevalent phenomenon whereby one gene influences multiple traits, contributing to phenotypic correlations (Solovieff et al., 2013). The developmental processes of flowers and leaves are intricately intertwined, exemplifying their close relationship (Mendez-Vigo et al., 2010). Some genes involved in flowering have pleiotropic effects on multiple traits, especially leaf development (Auge et al., 2019). In a case study, we used the DEM architecture, updated with a modified loss function for simultaneous double-trait prediction, to integrate Arabidopsis FT and RLN datasets for identification of pleiotropic genes that influence both flower and leaf development (see Figure 6A and methods). Overlapping samples among the FT and RLN data were first selected for double-trait modeling. The loss of each single-trait prediction was defined as the difference between the predicted and true phenotypic values. The total loss of the double-trait model, defined as the sum of the loss values derived from two single-trait predictions, was iteratively updated until the best-performing model with a minimum loss value was obtained for feature importance analysis and pleiotropic gene prediction.Figure 6 Effective prediction of pleiotropic genes.

(A) Updated DEM workflow featuring a modified loss function for simultaneous double-trait prediction.

(B and C) Performance comparisons between double-trait and single-trait models in terms of FT (B) and RLN (C) prediction. Modeling performance was assessed using mean values of PCC, MAE, and MSE (n = 10 repeats with the outer loop of nested cross-validation).

(D) Comparisons of importance scores predicted by double-trait and single-trait models for pleiotropic candidates. Violin plots depict median values (centerlines) and mean values (triangles). Statistical analysis was performed using Student’s t-test. ∗p < 0.05, ∗∗p < 0.01.

(E) Phenotypic analysis of wild-type (WT) Arabidopsis Columbia-0 and mutants of newly identified pleiotropic genes. Early (upper) and late (bottom) flowering phenotype images were captured for 23-day-old and 27-day-old plants grown in soil in a glasshouse under long-day (16-h light/8-h dark) conditions. Graphs represent the time (in days) to the opening of the first flower and the number of rosette leaves at the start of the flowering period. Data are presented as mean ± SE (three independent biological replicates, each replicate comprising 15 plants per genotype). ∗p < 0.05, ∗∗p < 0.01, and ∗∗∗p < 0.001 (Student’s t-test).

Because a single feature may have various levels of relevance to different traits, the proportion of irrelevant features for double-trait prediction is theoretically higher than that for each single-trait prediction case, potentially leading to degradation in the performance of the double-trait prediction model. To retain as many representative features as possible for double-trait modeling, we estimated the relevance degree of each feature for double-trait prediction and selected the most effective feature set with an optimal number of features from each type of omics data for training the double-trait model (Supplemental Figure 19). To optimize the performance of the double-trait model, three hyperparameters were fine-tuned (Supplemental Figure 20). As a result, the double-trait model yielded a 0.697 mean PCC, a 0.482 mean MAE, and a 0.363 mean MSE for FT prediction and yielded a 0.608 mean PCC, a 0.603 mean MAE, and a 0.492 mean MSE for RLN prediction, which were comparable to those of the single-trait models (Figure 6B and 6C).

By analyzing the feature importance levels of the double-trait model, 24 protein-coding genes were predicted to be putative pleiotropic genes associated with both FT and RLN (Supplemental Table 12). Notably, the importance levels of these pleiotropic candidates in double-trait modeling were significantly higher than those in single-trait modeling (Figure 6D), indicating that the influences of these pleiotropic candidates on individual traits are likely not dominant. Nevertheless, some of these genes have previously been reported to participate in regulating flowering and leaf development. For example, ULTRAPETALA1 (ULT1) acts as a putative transcriptional regulator and is involved in floral meristem differentiation and determinacy (Engelhorn et al., 2014). The ult1 mutant flowers late and produces more rosette leaves than wild-type Arabidopsis (Pu et al., 2013). SENSITIVITY TO RED LIGHT REDUCED 1 (SRR1) acts as an integrator between photoperiodic regulation and other pathways to maintain the repression of flowering in Arabidopsis (Johansson and Staiger, 2014), and loss-of-function srr1 mutation caused a delay in leaf initiation rate (Johansson et al., 2023).

In addition, several pleiotropic genes that affected both FT and RLN were newly identified. As shown in Figure 6E, single T-DNA insertion mutants of at3g12020 (Kinesin-like), at4g11040 (Reduced dormancy 5), and at5g45890 (Senescence-associated gene 12) exhibited phenotypes of early flowering and decreased RLN. Single mutants of at2g34320, at5g35390 (Pollen receptor like kinase 1), and at5g27830 (Folate receptor gene) displayed late flowering and increased RLN phenotypes. Interestingly, at5g27830 encodes a folate binding protein that helps to increase the stabilization of folates in plant foliage (Puthusseri et al., 2018). Loss-of-function at5g27830 mutation caused late flowering and an increase in RLN, implying a potential role of folates in the transition from vegetative to reproductive development.

Software for DEM applications

To facilitate the utilization of DEM architecture within the research community, we have developed a comprehensive software solution, aiming to streamline the prediction of phenotypes and the mining of functional genes and providing researchers with a user-friendly and versatile tool (Figure 7). The DEM software consists of four functional modules that can be launched separately, enabling users to make selective use of specific functionalities according to their analytical needs. The software supports both classification and regression tasks, accepting multi-modal omics datasets and phenotypic data as inputs. This flexibility enables users to customize dataset splits and adjust parameters according to their specific experimental designs. Whether initiating the modeling process from scratch or utilizing a pretrained model, the software offers the option to obtain predictions, enhancing its adaptability to diverse research scenarios. The DEM software and its user manual are freely accessible to users at https://github.com/cma2015/DEM/. We welcome user feedback to further refine and improve the DEM software for the benefit of the research community.Figure 7 Functional modules of DEM software.

DEM software comprises four functional modules. (1) Data preprocessing: this module focuses on preparing and refining input data for subsequent analysis, ensuring data compatibility with the DEM architecture. (2) Dual-extraction modeling: the core module of the software, which uses the dual-extraction architecture to extract representative features from multi-omics data for enhanced modeling. (3) Phenotypic prediction: a module dedicated to predicting phenotypes using the optimal DEM model and the processed multi-omics data, facilitating accurate and efficient predictions. (4) Functional gene mining: this module is designed for the identification of functional genes associated with complex traits, enhancing the understanding of molecular mechanisms. These modules can be launched independently, providing users with flexibility in the use of specific functionalities.

Discussion

In the postgenomic era, leveraging bio-omics knowledge is crucial for phenotypic prediction and gaining insights into the molecular mechanisms of complex traits. DL has emerged as a powerful tool for integrating multi-modal omics data, offering advanced techniques for modeling complex traits in plant breeding. However, the lack of universally applicable analysis methods has hindered progress, necessitating the development of an effective approach with superior performance for phenotypic prediction and functional gene mining. The innovations and limitations of this study are discussed below.

Appropriate architecture design for representative feature extraction

Existing multi-omics integration methods for supervised learning can be categorized into concatenation-based, model-based, and transformation-based methods (Reel et al., 2021). Concatenation-based methods construct a joint matrix from combined multi-omics data, but they often lack the ability to learn omics-specific knowledge. Model- and transformation-based methods, exemplified by MOGONET (Wang et al., 2021), DeepTTA (Jiang et al., 2022), and MOLI (Sharifi-Noghabi et al., 2019), create intermediate extraction models for each omics dataset and then build a final integrated model. However, they may overlook hidden cross-omics correlations in higher-dimensional feature spaces.

This study introduces the DEM architecture, designed to extract hidden omics-specific and cross-omics correlation information from high-dimensional feature space and enable representative feature extraction in a low-dimensional latent space with improved separability. Comparative assessments of various DL algorithms in high- and low-dimensional feature spaces demonstrate the superior performance of the DEM architecture.

Superior performance of DEM in complex trait studies

Accurate extraction of representative features from multi-omics data is pivotal for phenotypic interpretation in complex trait studies. The DEM architecture addresses this challenge by incorporating a dual-extraction modeling approach and the multi-head self-attention mechanism. Benchmarked against multiple published phenotypic prediction methods, DEM exhibited superior accuracy and effectively identified potential functional genes strongly associated with target traits, demonstrating its efficacy in discerning representative features in multi-modal omics data.

Generalizability is another critical concern in multi-omics integration analysis. Although many existing self-attention-related methods are tailored to specific species, prediction tasks, or omics modalities, DEM showcases superior performance across various tasks, including quantitative trait prediction in Arabidopsis and maize and phenotype classification in rice and human disease. Its flexibility in handling diverse omics data types, including the ability to handle noise information effectively, sets DEM apart from existing methods.

The interpretability of DEM is highlighted by its ability to identify pleiotropic genes that influence multiple traits, such as FT and RLN in plants. This aspect enhances its potential for unraveling complex biological relationships. The subsequent validation of predicted genes through experimental studies adds credibility to the predictions made by DEM, emphasizing its practical impact in real-world biological investigations.

Limitations of the study

First, multi-omics data collection focuses primarily on the upper part of the central dogma (DNA to mRNA) rather than the lower part (mRNA to protein), owing to constraints on data availability. Incorporating genomic, transcriptomic, and proteomic data into our multi-omics analysis could provide a more comprehensive understanding of biological interactions and potentially improve DEM predictions. Second, the number of features in the dataset used for model training exceeds the number of samples. The “curse of dimensionality” may impact the optimal predictive performance of our DEM model. Collecting more unpaired multi-omics datasets and utilizing reliable multi-omics imputation algorithms may help address this limitation. Third, while aiming for excellent performance, the complexity of the DEM model may compromise computational efficiency. Although DEM achieved the best predictive performance across all datasets with various batch sizes, a comparison of computational efficiency among the methods revealed that DEM required longer training times and more GPU resources than the other methods (Supplemental Figures 21 and 22).

In conclusion, this study not only introduces the innovative DEM architecture but also makes a significant contribution to the field by addressing key challenges in multi-modal omics data analysis for phenotypic prediction and functional gene mining. The combination of methodological advances, superior performance, interpretability, and user-friendly software positions DEM as a valuable tool for advancing our understanding of complex traits in the postgenomic era.

Methods

Data description

In this study, we used datasets derived from four species: Arabidopsis thaliana, Oryza sativa, Zea mays, and Homo sapiens. These datasets exhibited diversity in terms of species, traits, omics types, and population sizes, allowing for a comprehensive analysis. The Arabidopsis omics datasets, covering mRNA expression and DNA methylation, were sourced from the Gene Expression Omnibus database with accession numbers GEO: GSE43857 and GSE80744. The Arabidopsis SNP and phenotype datasets were obtained from the 1001 Genomes Project Data Center (https://1001genomes.org/data/GMI-MPI/releases/v3.1/). The maize RNA-sequencing data were downloaded from the CNCB–NGDC Genome Sequence Archive (https://ngdc.cncb.ac.cn/gsa/) with accession number CRA001241. The maize genomic variation and phenotype data were acquired from CropGS-Hub (Chen et al., 2024) with the ID GSTP004. The rice genomic variations and phenotypes were obtained from CropGS-Hub with the ID GSTP009. The human disease datasets (BRCA and LGG), including mRNA expression, miRNA expression, and DNA methylation, were obtained from TCGA (http://gdac.broadinstitute.org/). Details of these datasets are provided in Supplemental Table 13.

The reference genomes and annotations used were TAIR10 (downloaded from https://www.arabidopsis.org), MSU Rice Genome Annotation Project Release 7 (downloaded from http://rice.uga.edu/pub/data/Eukaryotic_Projects/o_sativa/), Zm-B73-REFERENCE-NAM-5.0 (downloaded from https://download.maizegdb.org/Zm-B73-REFERENCE-NAM-5.0/), and B73 RefGen_v4 (downloaded from https://ngdc.cncb.ac.cn/gwh/ncbi_assembly/55201/). The maize assembly chain files were downloaded from Ensembl Genomes (https://ftp.ebi.ac.uk/ensemblgenomes/pub/release-55/plants/assembly_chain/zea_mays/).

Genotype file processing

Genotype files in BED (Browser Extensible Data) format were converted to VCF (Variant Call Format) files using PLINK (version 1.9) (Chang et al., 2015). The genotypes of the required samples were retained, and INDELs were removed using VCFtools (version 0.1.16) (Danecek et al., 2011), along with quality control for missing data and minor allele frequency (vcftools –maf 0.05 –remove-indels –max-missing 0.9 –recode –recode-INFO-all). The raw VCF file of maize genotypes was converted to the assembly B73 RefGen_v4 from AGPv3 by CrossMap (version 0.7.0) (Zhao et al., 2014) and then converted to the genome coordinates of Zm-B73-REFERENCE-NAM-5.0.

RNA-sequencing data analysis

For maize RNA data processing, clean reads were aligned to the Zm-B73-REFERENCE-NAM-5.0 reference genome using HISAT2 (version 2.2.1) (hisat2 -k 5 –max-intronlen) (Kim et al., 2019). SAMtools (version 1.7) was used to extract high-quality reads and organize the output BAM files (Li et al., 2009). Quantitative analysis of transcript expression was performed using StringTie (version 2.1.7) (Kovaka et al., 2019) and featureCounts (version 2.0.1) (featureCounts -p -t exon -g gene_id) (Liao et al., 2014).

Omics feature preprocessing

To enhance the quality of the multi-omics datasets and mitigate the influence of redundant features and sequencing errors, we executed a feature preprocessing step in high-dimensional feature spaces. Initially, omics features with more than 25% null values were excluded (Zhang et al., 2018). Subsequently, we used 10 outer and 5 inner nested cross-validations (Cawley and Talbot, 2010) to partition the dataset. The following procedures were applied to different partitions: (1) null values in the training set, inner validation set, outer test set, and independent external validation set were imputed with the mean of the training set features; (2) the training set, inner validation set, outer test set, and independent external validation set were linearly transformed to the [0, 1] range using the maximum and minimum feature values.

Two distinct preprocessing pipelines were then implemented, tailored to the nature of the task, distinguishing between regression and classification tasks.

For regression tasks

The phenotypic values of the training set, inner validation set, outer test set, and independent external validation dataset were scaled using Z scores based on the mean and standard deviation of the phenotypic values. By comparing the modeling performance of DEM on regression tasks under different variance thresholds, the omics features with variances below 0.01 were filtered out (Supplemental Figure 23). To address challenges arising from the sparsity of the raw SNP dataset, a neural network comprising a locally connected layer and a fully connected layer was used. This network aimed to transform the sparse and discrete SNP features into dense and continuous features. To preserve the actual genotype information for model training, we applied 10-bit one-hot encoding to the 10 possible dinucleotide combinations of allelic variants (A/A, C/C, G/G, T/T, A/C, A/G, A/T, C/G, C/T, G/T) to represent the variant status of each SNP on the gene. Missing genotypes were encoded as a 10-bit zero, making them free of imputation, as the gradient descent algorithm was used in SNP-to-gene modeling. In the locally connected layer, the gene variation degree (Vgene) was calculated using Equation 1. The fully connected layer defined the predicted phenotypic value (Yˆ) and genomic variation features (L) using Equation 2. The loss function (LossSNP) of SNP transform modeling was defined as in Equation 3, considering the true and predicted phenotypic values.

The RF algorithm was then used to construct a phenotypic regression model to derive feature importance scores and select informative feature candidates from the training set. On the basis of the results of feature selection from the training set (Supplemental Figures 6 and 7), the consistent features in the validation set were used to validate the performance of the model. The model underwent nested cross-validation (10 outer, 5 inner). The phenotypic value (Yˆk) predicted by the RF model and its corresponding loss function (Lossk) were defined using Equations 4 and 5, respectively. The relevant score of the qth feature of the kth omics was then calculated as in Equation 6 on the basis of out-of-bag errors in the RF model. By assessing the impact of the number of top features on model performance, the most effective feature set with an optimal number of features was selected from each omics dataset for subsequent modeling.(Equation 1) Vgene=wSNP+b

Here, SNP denotes the SNP feature matrix consisting of 10-bit one-hot encoding, w denotes the weight of the SNP learned by the neural network, and b denotes the random bias of the locally connected layer.(Equation 2) Y,ˆL=wVgene+b

Here, w represents the gene variation weight and b represents the random bias of the fully connected layer.(Equation 3) LossSNP=∑m=1n|Yˆm−Ym|

Here, Ym denotes the true phenotypic value of the mth sample, Yˆm denotes the predicted phenotypic value of the mth sample, and n is the total number of samples.(Equation 4) Yˆk=RFR(Xk)

Here, RFR(·) represents the objective function of the RF model, Yˆk∈Rn×c denotes the predicted phenotypic values of n samples belonging to the kth omics type, and Xk∈Rn×dk denotes the dataset of the kth omics type.(Equation 5) Lossk=|Yˆk−Yk|

Here, Yk∈Rn×c denotes the true phenotypic labels of the n samples for the kth omics type.(Equation 6) xkq=1n∑i=1n(errOOBiq−errOOB′iq)

Here, n is the number of decision trees in RFR(Xk), set to 5000 in this study; errOOBiq denotes the qth out-of-bag error in the ith decision tree; and errOOB′iq denotes the out-of-bag error of the noisy qth feature in the ith decision tree.

For classification tasks

To filter out the omics features with no signals or low variances, we evaluated the modeling performance of DEM on classification tasks under different variance thresholds and selected a variance threshold of 0.05 for mRNA expression, DNA methylation, and miRNA expression in the SS, BR, BRCA, and LGG datasets (Supplemental Figure 23). To transform the sparse and discrete SNP features into dense and continuous features in classification datasets, we use a method akin to that used in regression tasks, with modifications made solely to the loss function (LossSNPC); it was defined as(Equation 7) LossSNPC=∑c=1n|Yˆc−Yc|

Here, Yˆc∈Rc denotes the predicted category probability of the mth sample, and Yc∈Rc denotes the true category probability of the mth sample.

The ANOVA F value of each feature was then calculated to assess whether a feature significantly differed across different phenotype classes. Features with the lowest F values were progressively removed until the first principal component of the remaining features explained less than 50% of the variance in the training set (Wang et al., 2021). The features consistent with the processed training set were then retained for the inner validation set, outer test set, and independent external validation dataset.

Feature extraction modeling

The modeling processes for both classification and regression tasks follow a two-step pipeline, involving the first extraction modeling in the high-dimensional feature space and the second extraction modeling in the low-dimensional feature space.

For classification tasks

In the initial step, individual omics feature matrices (Ok) and a multi-omics joint matrix served as separate inputs for the multi-head self-attention network in the high-dimensional feature space. The joint matrix was constructed by concatenating individual omics matrices, as expressed in Equation 8. Attention vectors for each omics feature were calculated using Equation 9, incorporating query (Q), key (K), and value (V) vectors. The multi-head attention weight (Wa) for each omics feature was computed according to Equation 10 using the softmax activation function (g(·)) and scaling coefficient (d). The fully connected FFN iteratively updated model weights, producing a final output (F) defined in Equation 11. Predicted phenotypes (Yˆ) and features (h) extracted from the last hidden layer were then defined in Equation 12, where a nonlinear activation function (σ(·)) was applied. The loss function (LossHC) for classification modeling in the high-dimensional feature space was defined in Equation 13, considering predicted and true phenotypic labels.

Features extracted from each high-dimensional modeling process were combined for input into the multi-head self-attention network in the low-dimensional feature space. The combined feature matrix (Xint) was formulated as in Equation 14. In the low-dimensional feature space, the multi-head self-attention modeling process was similar to that in the high-dimensional feature space, with a different loss function (LossLC) defined in Equation 15. The total loss function (Losstotal) for DEM was a combination of high-dimensional and low-dimensional components, as specified in Equation 16.(Equation 8) jointmatrix=Concat(O1,O2,…,Ok)

Here, Ok∈Rn×dk denotes the input of the kth omics, n is the number of omics samples, and dk is the number of input features of the kth omics.(Equation 9) Q,K,V=X·(WQ,WK,WV)

Here, Q denotes the query vector, K denotes the key vector, V denotes the value vector, X denotes omics features, and WQ, WK, and WV denote the weights of the query, key, and value vectors, respectively.(Equation 10) Wa=Concat(g(QKTd)VW1h,g(QKTd)VW2h,…,g(QKTd)VWih)

Here, g(·) denotes the softmax nonlinear activation function, Wih denotes the weight of the ith head, and d denotes the scaling coefficient of the attention weight.(Equation 11) F=Circlee(ReLU(Wa)P+b)

Here, P denotes the parameter weight of the FFN, b denotes random bias, Circlee(·) denotes looping the feedforward process and updating the model weights e times, and e denotes the number of encoders.(Equation 12) Yˆ,h=σ(F)P′+b

Here, σ(·) denotes a nonlinear activation function, P′ is the last hidden layer parameter weight of the FFN, and b denotes random bias.(Equation 13) LossHC=−∑m=1nc(∑k=1zlog(eyˆmk·ymc∑deyˆm,dk)+log(eyˆmjoint·ymc∑deyˆm,djoint))

Here, nc is the total number of samples, c is the number of phenotypic classes, z is the number of omics types, yˆmk∈Rc denotes the predicted category probability of the mth sample from the kth omics, yˆmjoint∈Rc denotes the predicted category probability of the mth sample from jointmatrix, ymc∈Rc denotes the true phenotypic label of the mth sample, and yˆm,dk and yˆm,djoint denote the dth element in the vectors yˆmk and yˆmjoint, respectively.(Equation 14) Xint=Concat(h1,h2,…,hk,hjoint)

Here, hk and hjoint denote the feature values calculated with Equations 9, 10, 11, and 12 using Ok and jointmatrix, respectively.(Equation 15) LossLC=∑m=1nc−log(eyˆm·ymc∑deyˆm,d)

Here, nc is the total number of samples, c is the number of phenotype classes, yˆm∈Rc denotes the predicted category probability of the mth sample from the feature integration matrix, ymc∈Rc denotes the true phenotypic label of the mth sample, and yˆm,d denotes the dth element in the vector yˆm.(Equation 16) Losstotal=LossHC+αLossLC

Here, α denotes the tradeoff parameter between the high-dimensional extraction component and the low-dimensional extraction component.

For regression tasks

The DEM process for regression tasks mirrored that for classification tasks, with distinctions in the loss functions. The loss function (LossHR) for regression modeling in the high-dimensional feature space was defined in Equation 17, considering predicted and true phenotypic values for individual omics and the joint matrix. In the low-dimensional feature space, the loss function (LossLR) for regression modeling was defined by Equation 18, focusing on predicted and true phenotypic values.(Equation 17) LossHR=∑m=1n(∑k=1z|yˆmk−ym|+|yˆmjoint−ym|)

Here, ym denotes the true phenotypic value of the mth sample, n is the total number of samples, z is the number of omics types, yˆmk denotes the predicted phenotypic value of the mth sample from the kth omics, and yˆmjoint denotes the predicted phenotypic value of the mth training sample from jointmatrix.(Equation 18) LossLR=∑m=1n|yˆm−ym|

Here, ym denotes the true phenotypic value of the mth sample, yˆm denotes the predicted phenotypic value of the mth sample, and n is the total number of samples.

In both classification and regression cases, the optimal DEM model was established by minimizing the total losses, consisting of a tradeoff between high-dimensional and low-dimensional extraction components (parameterized by α in Equation 16). For this study, α was set to 1.

Phenotypic prediction

The collaborative execution of phenotypic prediction relies on the integration of optimal DEM model weights and omics features.

The predicted phenotype, denoted as yDEM, is formulated as(Equation 19) yDEM=WDEM⨀(O1,O2,…,Ok)

Here, ⨀ represents the dot product operation, WDEM signifies the optimal DEM model weights, Ok∈Rn×dk denotes the input of the kth omics, n is the number of omics samples, and dk is the number of input features of the kth omics.

Post hoc interpretation method for functional gene mining

To evaluate the significance of omics features in the context of phenotypic prediction, we devised a post hoc interpretation method aiming to elucidate the functional gene landscape. This was achieved by analyzing the performance of DEM with raw feature values versus randomly shuffled ones.

The importance score (Scorej) for the jth feature was determined by evaluating the model’s predictive performance under shuffled conditions across multiple iterations. The calculation is expressed as(Equation 20) Scorej=E−1N∑n=1NEn,j

Here, N represents the number of times the jth feature is shuffled, E signifies the predictive performance of the DEM model with raw feature values, and En,j denotes the model performance achieved by shuffling the jth feature at the nth time. Model performance was assessed using cross-entropy for the classification task and MAE for the regression task.

Subsequently, an adjustment was introduced to contextualize the feature importance scores within each individual dataset. The adjusted importance score (Ascorej,k) for the jth feature in the kth dataset is defined as(Equation 21) Ascorej,k=Scorej,kEkNkQk

Here, Nk represents the number of samples in the kth dataset, Qk is the number of features in the kth dataset, and Ek denotes the cross-entropy (for classification) or MAE (for regression) of the model using the kth dataset.

Inspired by the mathematical concept of the rate of change, our objective was to identify trait-related features by observing shifts in feature ranking relative to their importance scores. When a change in ranking corresponds to a significant alteration in the importance score, this suggests that the feature is relevant for phenotypic prediction. Therefore, after sorting all features on the basis of their importance scores, we calculate the average rate of change between the feature importance and its ranking within each sliding window range. When the average rate of change for features within the window is less than 0.05, this suggests that, starting from this window, the rate of change for subsequent features is also less than 0.05. In other words, from the first feature within this window onward (inclusive of the first feature), the subsequent ranked features have minimal impact on the overall feature importance. Thus, when the first occurrence of an average rate of change below 0.05 is observed for a sliding window, the features before this window are considered candidate genes related to the phenotype.

The rate of change in the importance score (RCj,k) for the jth feature in the kth dataset is calculated as follows:(Equation 22) RCj,k=α·∑p=1m∑j=1n−1mn(Ascore(j+1),kp−Ascorej,kpRank(j+1),kp−Rankj,kp)

In this equation, n is the size of the sliding window containing n features, m is the number of sliding steps used by the feature window, Ascorej,kp denotes the jth adjusted feature importance score in the pth sliding window of the kth dataset, Rankj,kp denotes the rank of the jth adjusted feature importance score in the pth sliding window of the kth dataset, and α is an adjustment factor for the rate-of-change value. In this study, m and n were set to 10, and α was set to 100. Features with rates of change exceeding 0.05 in their importance scores were then identified and accumulated along with their corresponding genes.

Construction of the multi-trait prediction model

First, omics features from the FT and RLN datasets with variances below 0.01 were excluded. RF regression was then applied separately to select features highly relevant to FT and RLN phenotypes for subsequent DEM modeling.

In the high-dimensional feature space, the attention vector for each feature was obtained using Equation 9. The multi-head attention weight was computed as per Equation 10, and the FFN output was generated following Equation 11.

The predicted phenotypes of FT (YˆFT) and RLN (YˆRLN), along with the relevant features extracted from the hidden layer in the high-dimensional space (h), are defined as(Equation 23) YˆFT,YˆRLN,h=σ(F)P+b

Here, σ(·) denotes a nonlinear activation function, F denotes the output of the FFN, P denotes the parameter weight of the FFN, and b denotes the random bias of feature extraction modeling.

The loss function for multi-trait feature extraction modeling in the high-dimensional feature space (LossHigh) is defined as(Equation 24) LossHigh=∑m=1n(∑k=1z(|yˆFTk,m−yFTm|+|yˆRLNk,m−yRLNm|)+|yˆFTJm−yFTm|+|yˆRLNJm−yRLNm|)

Here, n is the total number of samples, z is the number of omics types, and the variables are defined as before.

In the low-dimensional feature space, the feature integration matrix was concatenated following Equation 14. The loss function of the modeling process in the low-dimensional feature space (LossLow) is defined as(Equation 25) LossLow=∑m=1n(|yˆFTm−yFTm|+|yˆRLNm−yRLNm|)

Here, n is the total number of samples, and the variables are defined as before.

Ultimately, the total loss function of multi-trait modeling (Losstotal) is defined as(Equation 26) Losstotal=LossHigh+γLossLow

Here, LossHigh represents the multi-trait prediction loss in the high-dimensional feature space, LossLow represents the multi-trait prediction loss in the low-dimensional feature space, and γ denotes the tradeoff parameter between the feature extraction and feature integration components in the multi-trait prediction model. In this study, γ was set to 1.

To identify features related to both FT and RLN, the feature importance scores of the FT and RLN models were calculated using Equation 20. For each trait model, the features were ranked by their adjusted importance scores obtained via Equation 21. The overlapping features among the top 100 ranked features of the two traits were identified as pleiotropic candidates.

Systematic benchmarks

The performance of DEM was benchmarked against five published methods: DNNGP and BRR for regression tasks, and AttentionMOI, CustOmics, and MOGAT for classification tasks. The DNNGP method was implemented using the Python software package available at https://github.com/AIBreeding/DNNGP/. The BRR methods were executed via BayesianRidge of the sklearn v.1.4.2 Python packages (https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.BayesianRidge.html). The AttentionMOI script was executed using the Python software package available at https://github.com/BioAI-kits/AttentionMOI/. The CustOmics method, accessed at https://github.com/HakimBenkirane/CustOmics, was implemented using its Python script. The MOGAT method was executed using the Python script provided at https://github.com/MezbahJUCSE39/MOGAT. To prevent the model training from cheating, each dataset was divided into a training set, an inner validation set, an outer test set, and an independent external validation set. Hyperparameters of each model were fine-tuned using nested cross-validation to achieve optimal prediction performance (Supplemental Figures 24 and 25).

Prediction of miRNA target genes

Three software programs, TargetScan (McGeary et al., 2019), PicTar (Krek et al., 2005), and miRTarBase (Huang et al., 2022), were used for miRNA target gene prediction with default parameters. The final miRNA targets were determined by identifying overlapping results across the three programs.

GO and KEGG enrichment analysis

The DAVID online tool (Sherman et al., 2022) was used to perform GO and KEGG pathway enrichment analyses. Significance was set at a Q value of <0.05 for enriched terms.

Plant materials, growth conditions, and phenotype measurements

A. thaliana genotypes, including the wild type (Columbia-0) and six T-DNA insertion mutants (Supplemental Figure 26), were sourced from AraShare (https://www.arashare.cn/index). Primers used for mutant identification are detailed in Supplemental Table 14. Arabidopsis seeds were germinated on half-strength Murashige and Skoog (1/2 MS) nutrient agar plates (PhytoTechnology Laboratories, Lenexa, KS, USA) and underwent a 2-day incubation at 4°C. Seedlings were grown in a glasshouse at 22°C under long-day conditions (16-h light/8-h dark), then transferred to soil for blooming under the same conditions. FT was recorded as the days from germination to the first flower opening, and the number of rosette leaves was measured at the onset of flowering.

Data and code availability

The DEM script codes and preprocessed datasets used for modeling are provided on GitHub (https://github.com/cma2015/DEM) and PyPI (https://pypi.org/project/biodem).

Funding

This work was supported by the 10.13039/501100001809 National Natural Science Foundation of China (32370723 , 32000410 ).

Author contributions

Z.M. and X.H. conceived the study and supervised the research. Y.R. and C.W. developed the computational framework and carried out benchmarking. H.Z. cultivated plant materials and executed phenotypic experiments. Z.M. drafted the manuscript with contributions from Y.R., C.W., and H.Z. All authors contributed to the interpretation of the manuscript.

Supplemental information

Document S1. Supplemental Figures 1‒26

Supplemental Tables 1‒14

Document S2. Article plus supplemental information

Acknowledgments

We thank the High-Performance Computing of Northwest A&F University for providing computing resources. No conflict of interest is declared.

Published by the Plant Communications Shanghai Editorial Office in association with Cell Press, an imprint of Elsevier Inc., on behalf of CSPB and CEMPS, CAS.

Supplemental information is available at Plant Communications Online.
==== Refs
References

Auge G.A. Penfield S. Donohue K. Pleiotropy in developmental regulation by flowering-pathway genes: is it an evolutionary constraint? New Phytol. 224 2019 55 70 10.1111/nph.15901 31074008
Benkirane H. Pradat Y. Michiels S. Cournède P.H. CustOmics: A versatile deep-learning based strategy for multi-omics integration PLoS Comput. Biol. 19 2023 e1010921 10.1371/journal.pcbi.1010921 36877736
Cawley G.C. Talbot N.L. On over-fitting in model selection and subsequent selection bias in performance evaluation J. Mach. Learn. Res. 11 2010 2079 2107
Chang C.C. Chow C.C. Tellier L.C. Vattikuti S. Purcell S.M. Lee J.J. Second-generation PLINK: rising to the challenge of larger and richer datasets GigaScience 4 2015 7 10.1186/s13742-015-0047-8 25722852
Chen J. Tan C. Zhu M. Zhang C. Wang Z. Ni X. Liu Y. Wei T. Wei X. Fang X. CropGS-Hub: a comprehensive database of genotype and phenotype resources for genomic prediction in major crops Nucleic Acids Res. 52 2024 D1519 D1529 10.1093/nar/gkad1062 38000385
Danecek P. Auton A. Abecasis G. Albers C.A. Banks E. DePristo M.A. Handsaker R.E. Lunter G. Marth G.T. Sherry S.T. The variant call format and VCFtools Bioinformatics 27 2011 2156 2158 10.1093/bioinformatics/btr330 21653522
Engelhorn J. Moreau F. Fletcher J.C. Carles C.C. ULTRAPETALA1 and LEAFY pathways function independently in specifying identity and determinacy at the Arabidopsis floral meristem Ann. Bot. 114 2014 1497 1505 10.1093/aob/mcu185 25288633
Gallusci P. Dai Z. Génard M. Gauffretau A. Leblanc-Fournier N. Richard-Molard C. Vile D. Brunel-Muguet S. Epigenetics for plant improvement: current knowledge and modeling avenues Trends Plant Sci. 22 2017 610 623 10.1016/j.tplants.2017.04.009 28587758
Holliday J.A. Wang T. Aitken S. Predicting adaptive phenotypes from multilocus genotypes in Sitka spruce (Picea sitchensis) using random forest G3 (Bethesda) 2 2012 1085 1093 10.1534/g3.112.002733 22973546
Huang H.Y. Lin Y.C.D. Cui S. Huang Y. Tang Y. Xu J. Bao J. Li Y. Wen J. Zuo H. miRTarBase update 2022: an informative resource for experimentally validated miRNA-target interactions Nucleic Acids Res. 50 2022 D222 D230 10.1093/nar/gkab1079 34850920
Jiang L. Jiang C. Yu X. Fu R. Jin S. Liu X. DeepTTA: a transformer-based model for predicting cancer drug response Brief. Bioinform. 23 2022 bbac100 10.1093/bib/bbac100 35348595
Johansson M. Staiger D. SRR1 is essential to repress flowering in non-inductive conditions in Arabidopsis thaliana J. Exp. Bot. 65 2014 5811 5822 10.1093/jxb/eru317 25129129
Johansson M. Steffen A. Lewinski M. Kobi N. Staiger D. HDF1, a novel flowering time regulator identified in a mutant suppressing sensitivity to red light reduced 1 early flowering Sci. Rep. 13 2023 1404 10.1038/s41598-023-28049-6 36697433
Kaczmarek E. Jamzad A. Imtiaz T. Nanayakkara J. Renwick N. Mousavi P. Multi-omic graph transformers for cancer classification and interpretation Pac. Symp. Biocomput. 27 2022 373 384 34890164
Kang M. Ko E. Mersha T.B. A roadmap for multi-omics data integration using deep learning Brief. Bioinform. 23 2022 bbab454 10.1093/bib/bbab454 34791014
Kim D. Paggi J.M. Park C. Bennett C. Salzberg S.L. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype Nat. Biotechnol. 37 2019 907 915 10.1038/s41587-019-0201-4 31375807
Kovaka S. Zimin A.V. Pertea G.M. Razaghi R. Salzberg S.L. Pertea M. Transcriptome assembly from long-read RNA-seq alignments with StringTie2 Genome Biol. 20 2019 278 10.1186/s13059-019-1910-1 31842956
Krek A. Grün D. Poy M.N. Wolf R. Rosenberg L. Epstein E.J. MacMenamin P. da Piedade I. Gunsalus K.C. Stoffel M. Rajewsky N. Combinatorial microRNA target predictions Nat. Genet. 37 2005 495 500 10.1038/ng1536 15806104
Li H. Handsaker B. Wysoker A. Fennell T. Ruan J. Homer N. Marth G. Abecasis G. Durbin R. 1000 Genome Project Data Processing Subgroup The Sequence Alignment/Map format and SAMtools Bioinformatics 25 2009 2078 2079 10.1093/bioinformatics/btp352 19505943
Liao Y. Smyth G.K. Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features Bioinformatics 30 2014 923 930 10.1093/bioinformatics/btt656 24227677
Liu J. Li J. Wang H. Yan J. Application of deep learning in genomics Life Sci. 63 2020 1860 1878 10.1007/s11427-020-1804-5
Long N. Gianola D. Rosa G.J.M. Weigel K.A. Application of support vector regression to genome-assisted prediction of quantitative traits Theor. Appl. Genet. 123 2011 1065 1074 10.1007/s00122-011-1648-y 21739137
Ma W. Qiu Z. Song J. Li J. Cheng Q. Zhai J. Ma C. A deep convolutional neural network approach for predicting phenotypes from genotypes Planta 248 2018 1307 1318 10.1007/s00425-018-2976-9 30101399
McGeary S.E. Lin K.S. Shi C.Y. Pham T.M. Bisaria N. Kelley G.M. Bartel D.P. The biochemical basis of microRNA targeting efficacy Science 36610.1126/science.aav1741 366 2019 eaav1741
Melzer S. Lens F. Gennen J. Vanneste S. Rohde A. Beeckman T. Flowering-time genes modulate meristem determinacy and growth form in Arabidopsis thaliana Nat. Genet. 40 2008 1489 1492 10.1038/ng.253 18997783
Mendez-Vigo B. de Andres M.T. Ramiro M. Martinez-Zapater J.M. Alonso-Blanco C. Temporal analysis of natural variation for the rate of leaf production and its relationship with flowering initiation in Arabidopsis thaliana J. Exp. Bot. 61 2010 1611 1623 10.1093/jxb/erq032 20190039
Michaels S.D. Amasino R.M. FLOWERING LOCUS C encodes a novel MADS domain protein that acts as a repressor of flowering Plant Cell 11 1999 949 956 10.1105/tpc.11.5.949 10330478
Pang J. Liang B. Ding R. Yan Q. Chen R. Xu J. A denoised multi-omics integration framework for cancer subtype classification and survival prediction Briefings Bioinf. 24 2023 bbad304 10.1093/bib/bbad304
Perez P. de Los Campos G. Crossa J. Gianola D. Genomic-Enabled Prediction Based on Molecular Markers and Pedigree Using the Bayesian Linear Regression Package in R Plant Genome-Us 3 2010 106 116 10.3835/plantgenome2010.04.0005
Pu L. Liu M.S. Kim S.Y. Chen L.F.O. Fletcher J.C. Sung Z.R. EMBRYONIC FLOWER1 and ULTRAPETALA1 act antagonistically on Arabidopsis development and stress response Plant Physiol. 162 2013 812 830 10.1104/pp.112.213223 23632855
Puthusseri B. Divya P. Veeresh L. Kumar G. Neelwarne B. Evaluation of folate-binding proteins and stability of folates in plant foliages Food Chem. 242 2018 555 559 10.1016/j.foodchem.2017.09.049 29037729
Reel P.S. Reel S. Pearson E. Trucco E. Jefferson E. Using machine learning approaches for multi-omics data analysis: A review Biotechnol. Adv. 49 2021 107739 10.1016/j.biotechadv.2021.107739 33794304
Sharifi-Noghabi H. Zolotareva O. Collins C.C. Ester M. MOLI: multi-omics late integration with deep neural networks for drug response prediction Bioinformatics 35 2019 i501 i509 10.1093/bioinformatics/btz318 31510700
Sherman B.T. Hao M. Qiu J. Jiao X. Baseler M.W. Lane H.C. Imamichi T. Chang W. DAVID: a web server for functional enrichment analysis and functional annotation of gene lists (2021 update) Nucleic Acids Res. 50 2022 W216 W221 10.1093/nar/gkac194 35325185
Solovieff N. Cotsapas C. Lee P.H. Purcell S.M. Smoller J.W. Pleiotropy in complex traits: challenges and strategies Nat. Rev. Genet. 14 2013 483 495 10.1038/nrg3461 23752797
Springer N. de León N. Grotewold E. Challenges of translating gene regulatory information into agronomic improvements Trends Plant Sci. 24 2019 1075 1082 10.1016/j.tplants.2019.07.004 31377174
Sugano S. Jiang C.J. Miyazawa S.I. Masumoto C. Yazawa K. Hayashi N. Shimono M. Nakayama A. Miyao M. Takatsuji H. Role of OsNPR1 in rice defense program as revealed by genome-wide expression analysis Plant Mol. Biol. 74 2010 549 562 10.1007/s11103-010-9695-3 20924648
Tanvir R.B. Islam M.M. Sobhan M. Luo D. Mondal A.M. MOGAT: A Multi-Omics Integration Framework Using Graph Attention Networks for Cancer Subtype Prediction Int. J. Mol. Sci. 25 2024 2788 10.3390/ijms25052788 38474033
Wallstrom S.V. Florez-Sarasa I. Araujo W.L. Aidemark M. Fernandez-Fernandez M. Fernie A.R. Ribas-Carbo M. Rasmusson A.G. Suppression of the external mitochondrial NADPH dehydrogenase, NDB1, in Arabidopsis thaliana affects central metabolism and vegetative growth Mol. Plant 7 2014 356 368 10.1093/mp/sst115 23939432
Wang K. Abid M.A. Rasheed A. Crossa J. Hearne S. Li H. DNNGP, a deep neural network-based method for genomic prediction using multi-omics data in plants Mol. Plant 16 2023 279 293 10.1016/j.molp.2022.11.004 36366781
Wang T. Shao W. Huang Z. Tang H. Zhang J. Ding Z. Huang K. MOGONET integrates multi-omics data using graph convolutional networks allowing patient classification and biomarker identification Nat. Commun. 12 2021 3445 10.1038/s41467-021-23774-w 34103512
Xu C. Jackson S.A. Machine learning and complex biological data Genome Biol. 20 2019 76 10.1186/s13059-019-1689-0 30992073
Yan J. Wang X. Machine learning bridges omics sciences and plant breeding Trends Plant Sci. 28 2023 199 210 10.1016/j.tplants.2022.08.018 36153276
Yan J. Xu Y. Cheng Q. Jiang S. Wang Q. Xiao Y. Ma C. Yan J. Wang X. LightGBM: accelerated genomically designed crop breeding through ensemble learning Genome Biol. 22 2021 271 10.1186/s13059-021-02492-y 34544450
Zhang H. Xie P. Xing E. Missing value imputation based on deep generative models Preprint at arXiv preprint 2018 10.48550/arXiv.1808.01684
Zhang L.L. Huang Y.Y. Zheng Y.P. Liu X.X. Zhou S.X. Yang X.M. Liu S.L. Li Y. Li J.L. Zhao S.L. Osa-miR535 targets SQUAMOSA promoter binding protein-like 4 to regulate blast disease resistance in rice Plant J. 110 2022 166 178 10.1111/tpj.15663 34997660
Zhao H. Sun Z. Wang J. Huang H. Kocher J.P. Wang L. CrossMap: a versatile tool for coordinate conversion between genome assemblies Bioinformatics 30 2014 1006 1007 10.1093/bioinformatics/btt730 24351709
Zhao K. Tung C.W. Eizenga G.C. Wright M.H. Ali M.L. Price A.H. Norton G.J. Islam M.R. Reynolds A. Mezey J. Genome-wide association mapping reveals a rich genetic architecture of complex traits in Oryza sativa Nat. Commun. 2 2011 467 10.1038/ncomms1467 21915109
