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

S2590-3462(24)00293-1
10.1016/j.xplc.2024.100985
100985
Resource Article
DeepCBA: A deep learning framework for gene expression prediction in maize based on DNA sequences and chromatin interactions
Wang Zhenye 1237
Peng Yong 147
Li Jie 1237
Li Jiying 5
Yuan Hao 123
Yang Shangpo 123
Ding Xinru 123
Xie Ao 123
Zhang Jiangling 3
Wang Shouzhe 146
Li Keqin 123
Shi Jiaqi 3
Xing Guangjie 3
Shi Weihan 3
Yan Jianbing 14
Liu Jianxiao liujianxiao@mail.hzau.edu.cn
1234∗
1 National Key Laboratory of Crop Genetic Improvement, Huazhong Agricultural University, Wuhan 430070, China
2 Hubei Key Laboratory of Agricultural Bioinformatics, Huazhong Agricultural University, Wuhan 430070, China
3 College of Informatics, Huazhong Agricultural University, Wuhan 430070, China
4 Hubei Hongshan Laboratory, Wuhan 430070, China
5 Microsoft Corporation, Redmond, WA 98052, USA
6 WIMI Biotechnology Co., Ltd., Changzhou 213000, China
∗ Corresponding author liujianxiao@mail.hzau.edu.cn
7 These authors contributed equally to this article.

10 6 2024
09 9 2024
10 6 2024
5 9 10098528 3 2024
25 5 2024
5 6 2024
© 2024 The Authors
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/).
Chromatin interactions create spatial proximity between distal regulatory elements and target genes in the genome, which has an important impact on gene expression, transcriptional regulation, and phenotypic traits. To date, several methods have been developed for predicting gene expression. However, existing methods do not take into consideration the effect of chromatin interactions on target gene expression, thus potentially reducing the accuracy of gene expression prediction and mining of important regulatory elements. In this study, we developed a highly accurate deep learning-based gene expression prediction model (DeepCBA) based on maize chromatin interaction data. Compared with existing models, DeepCBA exhibits higher accuracy in expression classification and expression value prediction. The average Pearson correlation coefficients (PCCs) for predicting gene expression using gene promoter proximal interactions, proximal-distal interactions, and both proximal and distal interactions were 0.818, 0.625, and 0.929, respectively, representing an increase of 0.357, 0.16, and 0.469 over the PCCs obtained with traditional methods that use only gene proximal sequences. Some important motifs were identified through DeepCBA; they were enriched in open chromatin regions and expression quantitative trait loci and showed clear tissue specificity. Importantly, experimental results for the maize flowering-related gene ZmRap2.7 and the tillering-related gene ZmTb1 demonstrated the feasibility of DeepCBA for exploration of regulatory elements that affect gene expression. Moreover, promoter editing and verification of two reported genes (ZmCLE7 and ZmVTE4) demonstrated the utility of DeepCBA for the precise design of gene expression and even for future intelligent breeding. DeepCBA is available at http://www.deepcba.com/ or http://124.220.197.196/.

This study reports the development of DeepCBA, a highly accurate, deep-learning-based framework for gene expression prediction based on maize chromatin interaction data. DeepCBA exhibits high accuracy in expression classification and expression value prediction, and it can identify important motifs involved in promoter proximal interactions and proximal-distal interactions in maize. Promoter editing and verification of ZmCLE7 and ZmVTE4 demonstrate the potential of DeepCBA for use in precise design of gene expression and future intelligent breeding.

Key words

maize
gene expression prediction
chromatin interactions
deep learning
promoter editing
regulatory elements and motifs
Published: June 10, 2024
==== Body
pmcIntroduction

Gene expression plays an important regulatory role in the development, growth, and reproduction of organisms, and a specific amount of gene product is produced in a particular spatiotemporal manner (Zrimec et al., 2020a, 2020b). Prediction of gene expression enables a better understanding of the mechanism and effect of sequence variation on transcriptional regulation and is complementary to population-based association analysis methods (Avsec et al., 2021). Therefore, developing accurate models for gene expression prediction and mining important variation sites can help to reveal the genetic basis of complex traits.

Gene expression is regulated by key genomic regulatory elements in DNA sequences, including promoters, enhancers, silencers, and insulators. Epigenetic features such as histone modifications, DNA methylation, transcription factors (TFs; Liu et al., 2021a, 2021b), and DNase I hypersensitive sites also play significant roles in the expression of target genes. Early studies mainly used traditional machine learning methods to predict gene expression from various types of data (Beer and Tavazoie, 2004; Cheng et al., 2011, Cheng et al., 2019; Dong et al., 2012; Tasaki et al., 2020), with room for improvement in accuracy. Deep learning methods have good fitting capability for processing of complex nonlinear data and can effectively extract complex features. To date, deep learning methods have been widely used in many fields, such as image processing and natural language processing. In recent years, researchers have also used deep learning for genome sequence analysis tasks, including the prediction of sequence functions (Zhou and Troyanskaya, 2015), TF binding sites (Zhao et al., 2021), chromatin interactions, and methylation status (Karlić et al., 2010; Schmidt et al., 2017). Researchers have also developed deep learning models to predict gene expression through the analysis of genome sequences; these include Basenji (Kelley et al., 2018), ExPecto (Zhou et al., 2018), Enformer (Avsec et al., 2021), and Chromoformer (Lee et al., 2022). ExPecto uses a convolutional neural network (CNN) to integrate 40-kb sequences upstream and downstream of gene promoters and uses spatial feature transformation and a linear regression model to predict human gene expression. Enformer uses the Transformer model to analyze the effect of distal elements as distant as 100 kb on gene expression. Studies have shown that chromatin interactions create spatial proximity between distal regulatory elements and target genes on the genome, which has an important effect on gene expression, transcriptional regulation, and phenotypic traits (Peng et al., 2019; Schoenfelder and Fraser, 2019). However, existing deep learning methods do not consider the effect of chromatin interactions (including promoter proximal regions and distal elements) on target gene expression, resulting in the capture of incomplete sequence information, thus affecting prediction accuracy. Although Enformer can predict the effect of distal elements within 100 kb on target gene expression, this method cannot capture the effect of regulatory elements at the genome-wide scale. In addition, the above methods have been used mainly in humans and mice, and there have been few studies in plants.

Maize has one of the largest crop cultivation areas in the world. It is not only the most important food crop but has also come to be used in industry and agriculture. Using the B73 reference genome sequence, chromatin interaction data, and gene expression data from multiple tissues, our study makes three major contributions.(1) We developed an accurate maize gene expression prediction model called DeepCBA (Deep neural networks of CNN module, BiLSTM and Attention mechanism) based on chromatin interaction data. DeepCBA has higher area under the receiver operating characteristic [ROC] curve (AUC) and Pearson correlation coefficient (PCC) values for gene expression classification and regression prediction tasks, and it increases the PCC for prediction of gene expression values by 46.3% compared with the traditional method. Experimental results showed that PCCs for gene expression prediction using promoter proximal interactions, proximal-distal interactions, and both proximal and distal interactions were 0.801, 0.621, and 0.923, respectively.

(2) The DeepCBA model identified 400–800 motifs that affect gene expression through chromatin interactions in maize shoot and ear tissues. These motifs have clear tissue specificity and are significantly enriched in expression quantitative trait loci (eQTLs) and open chromatin regions. The identified motifs are clustered into 6 groups of core sequences in the ear and shoot, and their distribution patterns can be divided into 5 and 4 categories of promoter proximal region interactions (PPIs) and proximal-distal region interactions (PDIs), respectively.

(3) Experimental results for detection of regulatory elements of ZmRap2.7 and ZmTb1, saturation mutations of the promoter regions of ZmCLE7 and ZmVTE4, and construction of cross-tissue and cross-genotype transfer learning models revealed the utility of DeepCBA for mining of distal regulatory elements, precise design of gene expression, and even future intelligent breeding.

Results

Gene expression prediction with the DeepCBA model

The experimental data used in this study are published maize chromatin interaction and expression data (Li et al., 2019) from shoots and ears. According to the type of element that interacts with the genes (see methods), the chromatin interaction data were divided into two categories, PPIs and PDIs. The average number of PPIs in shoot and ear tissue was 50 198, and the average number of PDIs was 11 198. The number of genes involved in the above interaction datasets was 23 707. To balance the number of genes in each category, we classified genes into unexpressed/expressed/highly expressed according to FPKM∈[0–0.1), FPKM∈[0.1–15), and FPKM∈[15–maximum] (FPKM, fragments per kilobase of transcript per million mapped reads). The gene expression values of different tissues were between 0 and 500, accounting for 99.6% of all genes in the maize genome. We defined the DNA sequence for a specific gene as including the sequences 1 kb upstream and 0.5 kb downstream of the transcription start site (TSS) and 0.5 kb upstream and 1 kb downstream of the transcription termination site (TTS) (Figure 1A) (Washburn et al., 2019).Figure 1 The workflow of DeepCBA.

(A) Two types of chromatin interactions: PPI and PDI. 1.5-kb gene proximal sequence of TSS and TTS.

(B) Five steps of DeepCBA: sequence encoding, feature extraction using CNN, temporal and distal feature extraction using BiLSTM, self-attention mechanism, and gene expression prediction.

(C) PCCs obtained using 3 methods for prediction of gene expression values. CNN_No_PPI: the CNN model using only gene upstream and downstream sequences. DeepCBA_No_PPI: the DeepCBA model using only gene upstream and downstream sequences. DeepCBA_PPI: the DeepCBA model using interaction sequences. Data augmentation denotes considering gene order in PPI mode during model training.

We developed a high-precision maize gene expression prediction model, DeepCBA, to make predictions based on chromatin interactions (Figure 1B). DeepCBA includes three modules, and a CNN is used to extract features of the encoded chromatin sequence and reduce the dimensionality. DNA sequences are usually double stranded, with the two strands connected by hydrogen bonds between bases, known as reverse complements. The bidirectional long short-term memory network (BiLSTM) can capture bidirectional and spatial information, and it has the ability to capture the dependencies between features by accessing long-range context. In this study, we used BiLSTM to capture distal interactions among chromatin sequence features. A self-attention mechanism was used to capture the contribution of key features for the model.

To evaluate the reliability of DeepCBA, we compared the accuracy of gene expression classification prediction for the following models: (1) CNN_No_PPI: a CNN model using the 3-kb sequences upstream and downstream of the gene; (2) CNN_PPI: a CNN model using the 6-kb PPI sequence; and (3) DeepCBA_PPI: the DeepCBA model using the 6-kb PPI sequence. DeepCBA has excellent model generalization ability in predicting gene expression classification. The gene expression classification prediction accuracy of DeepCBA using PPI data, or in PPI mode, was significantly better than that of CNN_No_PPI and CNN_PPI (Supplemental Figure 1). The PCCs of predicted gene expression for DeepCBA_PPI in the ear and shoot datasets were 0.954 and 0.967, and the PCCs of predicted gene expression for CNN_No_PPI in the 2 datasets were 0.309 and 0.52 (Figure 1C). To explore the effect of the order of interacting genes on the results of target gene expression, we considered gene order in PPI mode (data augmentation) and the number of chromatin interactions is twice that of the original PPI sequences (no data augmentation). The experimental results showed that DeepCBA can achieve better prediction results when the regulatory order between genes in PPI mode is taken into consideration (Figure 1C). These results show that DeepCBA has a higher accuracy than traditional methods of predicting gene expression.

DeepCBA identifies a dynamic range of gene expression values in different interaction modes

When only PDI sequences were used to predict gene expression, the PCCs of DeepCBA prediction results for the 2 tissues were 0.6214 and 0.6278, respectively, whereas when only PPI sequences were used to predict gene expression, the PCCs of the prediction results were 0.8060 and 0.8290. When both PPI and PDI sequences were used, the PCCs of DeepCBA for shoot and ear tissue were 0.9314 and 0.9266, respectively. When the order of genes was considered in PPI mode (data augmentation), the PCCs of the 2 tissues were 0.9539 and 0.9672 (Figure 2). These results show that the PCCs of DeepCBA were increased by 35.7%, 20%, and 50.2% when PDI, PPI, and PDI + PPI sequences were used, respectively. We also obtained expression classification and regression predictions for 3 other datasets: shoot (Peng et al., 2019), ear (Peng et al., 2019), and tassel (Sun et al., 2020) (see methods). Compared with other methods, DeepCBA showed an increase in regression prediction accuracy of 23.16%, 20.54%, and 19.73% in the 3 datasets (Supplemental Figure 2). For gene expression classification prediction, DeepCBA had a better average AUC than the other 2 methods (Supplemental Figure 3). These results reveal the significant role of information on chromatin interactive regulatory sequences in predicting gene expression and quantifying the effects of different interactive elements on expression regulation. As the size of the chromatin interaction dataset increased, the predictions of gene expression became more accurate.Figure 2 Performance of DeepCBA for prediction of maize gene expression in different modes.

From left to right are the results of predictions based on PDI, PPI, and PDI + PPI, respectively.

(A–D) The distribution of predicted values and true values of gene expression when the DeepCBA model was used to predict gene expression in shoots on the basis of PDI, PPI, and PDI + PPI.

(E–H) The distribution of predicted values and true values of gene expression when the DeepCBA model was used to predict gene expression in ears on the basis of PDI, PPI, and PDI + PPI.

Interestingly, there were some outliers in the above DeepCBA gene expression prediction results. These outliers all showed a tendency of higher true expression and lower predicted values. Taking the expression prediction with PPI sequences as an example, we compared the predicted gene expression value (Pre_exp) and the real gene expression value (Real_exp). When Pre_exp < 0.5 × Real_exp or Pre_exp > 1.5 × Real_exp, we regarded the gene as a candidate gene with a large prediction deviation. The candidate genes were then sorted according to deviation, and approximately 40% of these genes were found to participate in both PPI and PDI interactions (Supplemental Table 5). Moreover, the genes with biased predictions showed clear tissue specificity (Supplemental Figure 4A). These results suggest that an insufficient amount of chromatin interaction data may be the factor causing the deviation in gene expression predictions, revealing the complex regulatory network in the process of gene expression.

Genome-wide mining of PPI-related motifs that affect gene expression

To identify important motifs that affect gene expression, a saliency map (Chu, 2011) was used to calculate the gradient of chromatin sequence. The results showed that the sequence regions around 750-bp upstream of the gene TSS and 250-bp downstream of the gene TSS make a significant contribution to gene expression prediction (Figure 3A), consistent with previous results (Washburn et al., 2019). TF-MoDIsco (TF motif discovery from importance scores; Shrikumar et al., 2018) was next used to mine motifs (Supplemental Figure 4B and 4C), and 812 and 897 motifs were identified in the ear and shoot, respectively (Figure 3B). To assess the reliability of the predicted motifs, we used PlantTFDB (Jin et al., 2014, 2015, 2016; Tian et al., 2020) as the ground truth and queried these motifs against the database. The results showed that 87.6% and 41.14% of the motifs in the shoot and ear were matched with 653 verified conserved domains of higher plants (E value <0.5) (Supplemental Figure 8). The motifs identified in the ear and shoot were then compared with those bound by 104 maize TFs (Bailey et al., 2015, Grant et al., 2011, Tu et al., 2020), and the results showed that 40.4% and 39.4% of the 104 TFs could be matched in the two tissues (Supplemental Data 1 and 2), respectively. Interestingly, 4 motifs (ATTTAA, CAGGAA, TAATAT, and CACAGA) were confirmed to participate in the regulation of gene expression (Fu et al., 2013; Hufford et al., 2021; Liu et al., 2021a, 2021b) (Supplemental Figure 12A and 12B). These results show that the motif sequences identified with DeepCBA have biological significance. The detected motifs were positionally anchored in the 6-kb (3 kb per gene) sequence of the PPI sequence, and the distribution patterns of the motifs could be divided into 5 categories: (1) highly enriched near 250 bp downstream of the TSS; (2) only significantly enriched at specific sites, such as the TSS and TTS; (3) slightly enriched near 250 bp downstream of the TSS; (4) slightly enriched in the TSS and highly enriched in the TTS; and (5) evenly distributed throughout the entire sequence (Figure 3C). By comparing the number of motifs with different patterns in different tissues, we found that the largest proportion of motifs were found in the first category, consistent with the gradient results in sequences (Washburn et al., 2019). The number of overlapping motifs in the ear and shoot dataset was 95. The distribution of these 95 motifs in the 5 enrichment patterns was basically consistent between the ear and shoot (Supplemental Table 6). To further confirm that the motifs identified by DeepCBA affect gene expression, we combined the motifs in pairs to form a real motif composition (Real motif) and used motif compositions formed by random sequence combinations as controls (Random motif). The two kinds of motif compositions were inserted into 3-kb sequences, which were used to perform gene expression prediction (using N coding mode in the one-hot coding). The results showed that the effect on expression of the composited motifs identified by DeepCBA was significantly higher than that of the control group (Supplemental Figure 7B) (p = 2.06e−13). To further assess the effect of specific motif sequence changes on expression, we randomly selected 2 motifs (CCGCCG and CTCTCTC), mutated these 2 motifs in the test dataset, and predicted the expression of the corresponding genes (Supplemental Figure 7A). The numbers of genes in the ear and shoot were 2608 and 4176, respectively. According to the standard that the expression value changed by more than 50%, the results for the ear dataset showed that the proportions of the 2 motif mutations that affected gene expression were as high as 96.24% and 96.43%, respectively. In the shoot dataset, the proportions of the 2 motif mutations that affected gene expression were as high as 92.12% and 92.21% (Supplemental Figure 7C–7F). These results show that the discovered motifs have important functions in gene expression and regulation.Figure 3 Motifs that influence gene expression can be identified on the basis of PPI sequences.

(A) Effect on expression prediction of 2 interacting sequences input into the DeepCBA model.

(B) Venn diagram showing the motifs identified in shoots and ears by DeepCBA in PPI mode.

(C) Five different distribution patterns of motifs identified in shoots and ears: (1) highly enriched near 250 bp downstream of the TSS, (2) highly enriched at specific positions, (3) poorly enriched near 250 bp downstream of the TSS, (4) poorly enriched near TSSs but highly enriched near TTSs, (5) evenly distributed across the whole sequence.

(D and E) Core motif sequences obtained using MetaLogo on the basis of motifs identified in ears and shoots in PPI mode. Six core sequences were obtained in the 2 tissues.

(F and G) Changes in the expression of expressed and highly expressed genes containing different numbers of motifs in ears and shoots.

To further examine whether the motifs identified by DeepCBA have sequence similarities in different categories, MetaLogo (Chen et al., 2022) was used to cluster the motifs into core sequences. Six groups of motifs with similar core sequences were identified in the ear and shoot (Figure 3D, 3E, and Supplemental Figure 9). On the basis of the clustering results, we compared the expression of genes containing different numbers of motifs in their gene sequences. Gene expression tended to increase as the number of motif compositions increased (Figure 3F and 3G), implying that gene expression is the result of joint regulation by different factors.

The motifs identified by DeepCBA reveal the regulation of gene expression

To analyze the apparent characteristics and biological functions of the motifs identified by DeepCBA in PPI mode, the motifs identified in the ear were positionally anchored onto the original gene sequences. We determined the physical location of each motif on the chromosome (Supplemental Data 5) and matched motifs to published eQTLs (Tian et al., 2023). The following 2 methods were used to obtain control sequences: (1) removing motif sequences from PPI sequences and selecting sequences of equal length from the remaining PPI sequences and (2) removing PPI sequences from the whole DNA genome and selecting sequences from the remaining genome sequences. Compared with the controls, the motifs identified by DeepCBA were significantly enriched at eQTLs (Figure 4A and 4B) (∗∗∗∗p < 0.0001, t test). We performed 100 repeated experiments in the above 2 controls to reduce the error introduced by randomized experiments. In addition, we matched these motifs to the open chromatin regions identified in 26 lines of the Nested Association Mapping (NAM) population (Rodgers-Melnick et al., 2016, Woodhouse et al., 2021). The results showed that the motifs identified by DeepCBA were significantly enriched in open chromatin regions identified in the NAM population (Figure 4C and 4D) (∗∗∗∗p < 0.0001, t test). Similarly, the physical locations of important motifs identified in the shoot were significantly enriched in eQTLs and open chromatin regions (Supplemental Figure 10) (∗∗∗∗p < 0.00001, t test). Through matching with 104 TFs (Tu et al., 2020) and eQTLs, a CATGCA motif was identified in the gene sequence of Zm00001d042609. The motif and the downstream gene Zm00001d042600 can be bound simultaneously by the TF NACTF109 (Supplemental Data 9). Variation in CATGCA in the maize association mapping panel was associated with differences in the expression of Zm00001d042600, which affected drought resistance at the seedling stage (Figure 4E). CATGCA (RY-motif) is a highly conserved motif present in the promoters of many seed-specific genes and plays an important role in seed development (Mönke et al., 2004). The TFs ABI3 and FUS3 play important regulatory roles in the development and maturation of Arabidopsis seeds, and they can bind to the RY-motif to regulate the process of abscisic acid–mediated endosperm maturation (Guerriero et al., 2009). The maize homolog of ABI, Vp1, is also a key factor in the regulation of seed maturation. In addition, Vp1 is expressed in the phloem cells of vegetative tissues under drought stress (Cao et al., 2007). ZmABI19, a TF containing the B3 domain, can also bind to the RY-motif upstream of the grain-filling-specific TF gene Opaque2 (O2) promoter to perform transactivation to regulate endosperm development in maize. In addition, deletion of the RY-motif greatly reduced promoter activity in the regulatory regions of the legumin gene of Vicia faba and napin in Brassica napus (Reidt et al., 2000). The above cases confirm that the motifs identified by DeepCBA have important biological functions in multiple species.Figure 4 Epigenetic features and examples of gene expression regulation of motifs identified by DeepCBA (ears).

(A) The matching number of motifs with different lengths and eQTLs in PPI mode. The motif sequences were removed from the PPI sequences, and sequences with lengths of 6–10 were randomly selected from the remaining PPI sequences as controls (∗∗p < 0.05, ∗∗∗p < 0.01, ∗∗∗∗p < 0.0001; t test).

(B) The matching number of motifs with different lengths and eQTLs in PPI mode. PPI interaction sequences were removed from the whole genome, and sequences with lengths of 6–10 were randomly selected from the remaining sequences as controls (∗∗p < 0.05, ∗∗∗p < 0.01, ∗∗∗∗p < 0.0001; t test).

(C and D) The matching number of motifs identified by DeepCBA in PPI mode and open chromatin regions in the NAM population. Controls were selected as described in (A) and (B).

(E) For the CATGCA motif identified in the Zm00001d042609 sequence in PPI mode, the motif and the downstream gene Zm00001d042600 can be bound simultaneously by the TF NACTF109. Variation in CATGCA in the maize association mapping panel was associated with differences in the expression of Zm00001d042600 and thus affected maize drought resistance at the seedling stage.

Using published articles and database reports, we identified seven important motifs involved in the regulation of gene expression, plant growth and development, and other processes. For instance, the RY-motif (CATGCA) is reported to be involved in the regulation of endosperm development (Yang et al., 2021). The Y-patch (TC-motif) is a core regulatory element that can enhance promoter activity (Jores et al., 2024), and some experiments have shown that EjBZR1 can bind to the BRRE-motif in the EjCYP90A promoter to regulate gene expression and fruit cell enlargement (Su et al., 2021). Descriptions of the functions of some important motifs detected by DeepCBA in PPI mode are provided in Supplemental Table 7.

DeepCBA reveals regulatory motifs for gene expression in PDI mode

The gradient results of motifs obtained in PDI mode show that the region approximately 250 bp downstream of the gene TSS is more likely to affect gene expression (Supplemental Figures 5A and 13A). Supplemental Figure 5B shows the motifs identified by DeepCBA based on PPI and PDI sequences. The motifs identified in the 2 tissue types show clear tissue specificity, and the motifs identified by the 2 modes (PPI and PDI) in the same tissue also differ. These results indicate that there are differences in the factors that regulate gene expression through PPI and PDI (Supplemental Figure 13B). To verify the reliability of the motifs identified in PDI mode, we compared these motifs to the PlantTFDB (Jin et al., 2014, 2015, 2016; Tian et al., 2020) (E <0.5). We found that 51.52% and 48.64% of the motifs identified in shoots and ears in PDI mode matched 653 conserved structural regions of higher plants. In addition, 58% and 56% of 104 maize TFs (Tu et al., 2020) matched with the motifs identified in ears and shoots (Supplemental Data 3 and 4). As observed for motifs identified in PPI mode, the motifs identified by DeepCBA in PDI mode (Supplemental Data 6) matched with eQTLs and open chromatin regions. The following 2 methods were used to generate control motifs: (1) removing motif sequences from PDI sequences and selecting sequences of equal length from the remaining PDI sequences, and (2) removing PDI sequences from the whole DNA genome and selecting sequences from the remaining genome sequences. Compared with the controls, the motifs identified by DeepCBA were significantly enriched at eQTLs (Supplemental Figure 5C and 5D) (∗∗∗∗p < 0.0001, t test). We next matched these motifs with the open chromatin regions identified in 26 NAM population lines. The results showed that the motifs identified by DeepCBA were significantly enriched in the open chromatin regions identified in the NAM population (Supplemental Figure 5E and 5F) (∗∗∗∗p < 0.0001, t test). Similarly, the physical locations of important motifs identified in shoots were significantly enriched in eQTLs and open chromatin regions (Supplemental Figures 10 and 11) (∗∗∗∗p < 0.00001, t test).

The physical locations of the motifs identified in PDI sequences of different tissues were clustered and divided into 4 categories: (1) highly enriched near 250 bp downstream of the TSS, (2) highly enriched at specific sites, (3) poorly enriched near 250 bp downstream of the TSS, and (4) distributed evenly throughout the whole sequence (see Figure 6A). The forms of sequence importance of the first and second categories were complementary, suggesting that there may be differences in the way they work. Similar to the results in PPI mode, the distribution of the 44 overlapping motifs in the 4 enrichment patterns was basically the same in ears and shoots (Supplemental Table 6). MetaLogo (Chen et al., 2022) was used to cluster the core sequences of motifs, and 6 groups of core sequences were identified in ears and shoots (Supplemental Figure 6B). Two motifs, GGCCCA and AAAAAA (Supplemental Figures 6C, 6D, and 13C), were also reported in previous studies (Peng et al., 2019; Woodhouse et al., 2021). Further analysis showed that the conserved region in which the TCP TF binds to DNA is also the GGCCCA motif in maize. Therefore, TFs (e.g., TCP) are likely to play an important role in regulating gene expression through distal elements. GGCCCA, also known as the site II motif, was identified in the promoter regions of various highly expressed genes such as ribosomal and DEAD-box RNA helicase genes. TCP and ASR5 TFs are examples of proteins known to bind to the GGCCCA motif. Moreover, Xu et al. (2011) identified several diurnal-related cis elements in seedlings and leaves, including element II of Arabidopsis proliferating cell nuclear antigen 2 (∗GGCCCA∗ or ∗AGCCCA∗).

DeepCBA also detected 12 important motifs in PDI mode whose functions included enhancement of promoter activity, regulation of gene expression, formation of immune complexes, and resistance to heat or salt stress. In particular, GGCCCA has been found in the promoter regions of various highly expressed genes and is a binding site for TCP and ASR5 TFs (Xu et al., 2011; Oka et al., 2017). In addition, association of the Telo box (AAACCTA) with site II (GGCCCA) or TEF cis-acting elements appears to be involved in ribosome biogenesis (Gaspin et al., 2010). The ERF subfamily has been reported to bind to the GCC box (GCCGCC) in response to biotic stress and can also respond to ethylene by enhancing gene expression (Ishige et al., 1999; Wu et al., 2020). Descriptions of the functions of some important motifs detected by DeepCBA in PDI mode are provided in Supplemental Table 8.

DeepCBA identifies regulatory elements in the maize genes ZmRap2.7 and ZmTb1

On the basis of a series of motifs identified in PDI mode, we performed in-depth analysis of the Vgt1 regulatory site of the maize flowering-related gene ZmRap2.7 (Zhao et al., 2018; Ricci et al., 2019). We extracted the 70-kb sequence upstream of ZmRap2.7, split the sequence into 1.5-kb sub-sequences, and used them as input for the DeepCBA model. A saliency map (Chu, 2011) was used to calculate the sequence gradient, and the results showed that the important motifs identified by DeepCBA were found mainly in open chromatin regions (Figure 5A). We narrowed down the regions to two 500-bp regions (chr8: 135 941 716–135942216 and chr8: 135945716–135946216). After matching the identified motif with 104 TFs, we performed enrichment analysis for the above 2 regions. There were 18 TF-binding sites in the first region and 19 TF-binding sites in the second region (Figure 5B). Importantly, there were 16 common TFs in the 2 regions (Supplemental Table 9 and 10 and Supplemental Data 10). Similarly, we extracted a 3-kb sequence upstream of the TSS of ZmRap2.7, used it as input to the DeepCBA model for training, and calculated the sequence gradient (Figure 5C). The detected motifs were matched with open chromatin regions, and the motifs of 7 identified TFs were found to overlap with the open chromatin regions (Figure 5D). We also performed sequence analysis of the maize tillering-related gene ZmTb1 and its regulatory elements. The 70-kb sequence upstream of ZmTb1 was selected for verification using DeepCBA, and the important motifs identified were also found mainly in open chromatin regions (Supplemental Figure 14A–14C). For the 10-kb (chr1: 270482176–270492176) region with the highest gradient value, we trained the DeepCBA model using the PDI sequence and detected 314 and 514 motifs in ears and shoots, respectively, involving 93 TFs. Similarly, we detected 140 and 262 motifs in the 3-kb (chr1: 270549176–270552176) region, with the second highest gradient value in ears and shoots, involving 93 TFs. There were 59 and 118 overlapping motifs in the above 2 regions, involving 84 TFs (Supplemental Figure 14D and 14E). These results reveal that there are similar TF binding clusters in distal regulatory elements and their regulated genes, thus enabling a deep analysis of gene expression regulation.Figure 5 Identification of regulatory elements in ZmRap2.7.

(A) Distribution of open chromatin regions in 70-kb bins upstream of ZmRap2.7 in different tissues and the sequence gradient values calculated by DeepCBA.

(B) Identified motifs and TFs that can be bound in the 2 regions with the highest gradient values (chr8: 135941716–135942216 and chr8: 135945716–135946216).

(C and D) The motifs and TFs that can be bound in the 3-kb region upstream of the TSS of ZmRap2.7.

①②③④ represent DNase-sequencing, assay for transposase accessible chromatin-sequencing, H3K4me3, and H3K9ac, respectively.

DeepCBA accurately predicts gene expression values through promoter saturation mutations

CRISPR-Cas9 is another method for generating weak alleles by targeting coding regions, and it has been widely used to edit cis-regulatory regions in different plants (Rodríguez-Leal et al., 2017; Liu et al., 2021a, 2021b; O’Connor et al., 2020, Song et al., 2022). However, screening lines with continuous expression-gradient changes in target genes from a large amount of gene-edited material comes with many uncertainties. Therefore, it is crucial to utilize computational tools to predict and screen variations in continuous gradient expression by performing saturation mutations of gene promoters. For validation, we selected the editing results for the promoter region of maize ZmCLE7 (Liu et al., 2021a, 2021b). ZmCLE7 affects yield by changing the ear phenotype, and ear tissue was used in this study to build the deep learning prediction model. The 4-kb upstream region (chr4: 8334400–8338400) of ZmCLE7 was selected as the candidate editing region (Figure 6A). Combined with the published results, the CRISPR-Cas9-edited sequences were used as input for the DeepCBA model to predict gene expression (Figure 6B). The predicted gene expression showed a trend consistent with that measured experimentally, and the correlation between the predicted gene expression after editing the target segment and the qPCR value was 0.51 (Figure 6C). To explore more precisely how the 4-kb target sequence affects the expression of ZmCLE7, we used a sliding-window method (window size = 200 bp, step size = 200 bp) to process the 4-kb sequence and obtained 12 sequences with a length of 3 kb (Figure 6D). We then used the DeepCBA model to predict gene expression for these 12 edited sequences. The results showed that a wider range of expression variation types could be produced than were produced in the biological experiment results (Figure 6E).Figure 6 DeepCBA edits the maize genes ZmCLE7 and ZmVTE4 to achieve accurate expression prediction.

(A) Distribution of 4 histone modifications (H3K27ac, H3K4me3, H3K27me3, H3K9ac) and open chromatin regions within the 4-kb upstream region of ZmCLE7.

(B) Schematic diagram of 6 pieces of editing information for ZmCLE7 in the published literature (Liu et al., 2021a, 2021b).

(C) DeepCBA was used to predict the expression of gene-edited sequences in (B) and compare it with quantitative real-time PCR (qPCR) results from the published literature.

(D) Using sliding windows (window size = 200 bp, step size = 200 bp) to process the 4-kb sequence of ZmCLE7.

(E) Expression of the edited sequences in (D) predicted using DeepCBA.

(F) Distribution of 3 histone modifications (H3K27ac, H3K4me3, and H3K9ac) and open chromatin regions within the 4-kb region upstream of ZmVTE4.

(G) Gene editing events in the 4-kb region upstream of ZmVTE4.

(H) Comparison of ZmVTE4 expression predicted by DeepCBA and ZmVTE4 expression measured by leaf quantitative real-time PCR for different gene editing events in the 4-kb upstream region of ZmVTE4.

To further verify the reliability of the DeepCBA model in gene editing applications, we edited the promoter (chr5: 205820586–205829816) of ZmVTE4, which affects the vitamin E content of maize. Combined with the distribution characteristics of different histone modifications in the target region (Figure 6F), we designed 7 fragment deletion types for CRISPR-Cas9 editing (Figure 6G). These sequences were used as input for the DeepCBA model to predict expression in ear tissue. We also extracted RNA from the individual edited plants and detected the expression of target genes; the correlation between the predicted gene expression after editing of the target segment and the qPCR value was 0.57 (Figure 6H). The predicted gene expression and real gene expression showed similar changes in response to promoter editing. Taken together, these results demonstrate the feasibility of DeepCBA for promoter saturation mutations, confirming that it is a powerful tool for the precise design of desired gene expression.

DeepCBA enables cross-tissue and cross-genotype prediction of maize gene expression

With the development of machine learning technologies, many methods have been used to study different tissues, materials, and species (Kelley et al., 2018). To further explore the feasibility of predicting gene expression between different tissues and materials, we constructed a transfer learning model to enable prediction of gene expression across tissues (ears and shoots) and materials (B73 and SK [Yang et al., 2019]) (Supplemental Figure 16). We also used a new shoot dataset (Peng et al., 2019) to further verify the generalization ability of DeepCBA. The new shoot dataset (Peng et al., 2019) contained 43 865 PPIs involving 20 695 genes. The SK shoot dataset (Yang et al., 2019) contained 7394 PPIs involving 7099 genes. The PCC of DeepCBA increased from 0.7601 to 0.8687 after applying transfer learning to the shoot dataset (Yang et al., 2019) (Supplemental Figure 17A). In addition, we compared the runtime of DeepCBA for gene expression prediction with and without transfer learning using different datasets. The results indicated that transfer learning could improve prediction accuracy while reducing the time required for model training and prediction (Supplemental Figure 17B). For cross-tissue prediction of gene expression, we first predicted gene expression for shoots (Li et al., 2019) and ears (Li et al., 2019) using a prediction model trained with PPI and PDI datasets from shoots (Peng et al., 2019) (Supplemental Figure 18A). The PCCs between predicted expression and true expression were 0.8811 and 0.888. We next predicted gene expression for two shoot datasets (Peng et al., 2019 and Li et al., 2019) (Supplemental Figure 18B) using a prediction model trained with PPI and PDI datasets from ears (Li et al., 2019) and obtained PCCs of 0.8848 and 0.8737. We also predicted gene expression of ears (Li et al., 2019) and shoots (Peng et al., 2019) using a prediction model trained with PPI and PDI datasets from shoots (Li et al., 2019) and obtained PCCs of 0.8931 and 0.8689 (Supplemental Figure 18C). To test cross-genotype and cross-tissue predictions, we predicted the gene expression of SK shoot tissue using a model trained on 3 maize B73 datasets (shoot [Peng et al., 2019], ear [Li et al., 2019], and shoot [Li et al., 2019]); the resulting PCCs were 0.8687, 0.8583, and 0.8623, respectively (Supplemental Figure 18D–18F). In all, the PCCs for gene expression prediction among different tissues and materials through transfer learning typically exceeded 0.85. These findings provide a reference for prediction of gene expression across different tissues and materials and broaden the application scope of the DeepCBA model.

DeepCBA website

To facilitate open access of the DeepCBA model, we developed the DeepCBA website (http://www.deepcba.com/ or http://124.220.197.196/), which enables high-precision gene expression prediction based on chromatin interactions in maize and 3 other crops (rice, cotton, and wheat). Users can select any of the 4 crops and input any interaction sequences that meet the requirements for predicting the expression of related genes and sequences. In addition, the website provides a visualization interface to display the gradient importance of the input sequences (Figure 7).Figure 7 The DeepCBA website.

(A) Functions of the DeepCBA website.

(B) DeepCBA enables high-precision gene expression prediction based on chromatin interactions for 4 crops: maize, rice, cotton, and wheat. Users can freely select relevant models to achieve the prediction tasks.

(C) DeepCBA implements a parallel computing algorithm. The prediction results are sent to users via e-mail, and users can view the results according to the Job_id.

(D) DeepCBA provides a visualization interface to display the gradient importance of the input sequences that affect gene expression.

Discussion

Coding regions account for only a small percentage of the maize genome, and most of the genome consists of noncoding regions. Many functional loci have been identified in the noncoding regions of maize through association analysis, and several kinds of regulatory elements have been identified through epigenetics analyses at the genome-wide level. However, it is still unclear how regulatory elements in noncoding regions accurately regulate gene expression. Using deep learning tools to predict the contributions of different regulatory elements to gene expression has important biological significance. As we know, chromatin interactions have an important effect on target gene expression. To what extent gene expression is determined by the chromatin interactions of DNA sequences is therefore an important question.

In this study, we developed the DeepCBA model for high-precision gene expression prediction based on maize chromatin interactions. A CNN is used to extract local features of the DNA sequence, BiLSTM is innovatively used to capture the relationships between distal features, and the self-attention mechanism is used to capture important features. Compared with existing methods, DeepCBA showed higher accuracy in classification of gene expression and prediction of expression values. DeepCBA predicts gene expression accurately by integrating chromatin interaction (PPI and PDI) data from different maize tissues and quantitatively assesses the effects of different regulatory elements on gene expression. DeepCBA revealed that the average contributions of PPI, PDI, and PPI + PDI to gene expression prediction were 0.817, 0.625, and 0.929, representing increases of 0.357, 0.165, and 0.469 compared with the single-sequence method.

Unraveling the black box of deep learning-based applications remains a challenge in biological research. To interpret the reasoning process of DeepCBA, we used a saliency map to calculate the model gradient through reverse calculation and obtained a significance map of DNA sequences. We identified important motifs in the chromatin interaction sequences, together with other latent features yet unknown for the prediction of gene expression. Verification against known databases and the published literature showed that the detected motifs were highly reliable. In terms of molecular characteristics, the motifs identified in this study were enriched mainly in eQTLs and open chromatin regions. Moreover, gene expression trended upward as the number of motifs in the motif composition increased (Figure 3F and 3G). The detected motifs and gradient results from different maize tissues (ears and shoots) showed clear tissue specificity (Supplemental Figure 5B). The identified motifs were clustered into 6 groups of core sequences in ears and shoots (Supplemental Figure 9). In addition, the distribution patterns of the motifs could be divided into 5 and 4 categories in PPI and PDI mode, respectively (Figure 3C and Supplemental Figure 6A).

Alleles that control important traits (e.g., crop yield, resistance) often alter the expression levels of genes, thereby affecting the phenotype. Editing the promoter region helps to intelligently design the expression of target genes, achieving the goal of improving crop yield and resistance (Rodríguez-Leal et al., 2017; Liu et al., 2021a, 2021b; Song et al., 2022). DeepCBA can detect the gradient effect of regulatory elements in noncoding regions at the single-base level and accurately evaluate the functional loci of regulatory elements. The feasibility of DeepCBA for exploration of regulatory elements that affect gene expression was validated using the maize flowering-related gene ZmRap2.7 and the tillering-related gene ZmTb1 (Figure 5; Supplemental Figure 14). Through saturation mutations in the promoter and regulatory regions of specific genes (ZmCLE7 and ZmVTE4), we achieved de novo gene expression prediction (Figure 6). These results validated the reliability of DeepCBA using real examples in maize, and the model can be used widely for precise design of target gene expression levels in different crops, thereby enabling intelligent design and breeding.

To further explore the feasibility of expression prediction in practical applications, we constructed a DeepCBA transfer learning model to enable expression prediction across tissues (ears and shoots) and materials (Supplemental Figure 18). The PCCs of expression prediction across tissues and genotypes with the transfer learning model exceeded 85%, confirming the wide application scope of DeepCBA. To facilitate use of DeepCBA, we developed a user-friendly website for predicting gene expression and visualizing sequence importance in four crops (maize, rice, cotton, and wheat).

Interestingly, the expression values of some genes predicted by DeepCBA were lower than their actual expression values. This is perhaps related to the small number of chromatin interactions in different tissues and materials. That is to say, there were insufficient data for training DeepCBA and influencing the prediction effect. Prediction accuracy will also be improved by integrating multi-omics data, including data on open chromatin regions, TF binding sites, histone modifications, and DNA methylation.

Artificial intelligence is continuing to penetrate various fields, and machine learning offers advantages for exploring the effects of different regulatory elements on gene expression. With the development of different deep learning algorithms, our understanding of different regulatory elements will gradually become clearer. We will have a deeper understanding of the effect of variation on gene expression when considering tissue and spatiotemporal specificity. This study will also provide a theoretical basis for accurate design of gene expression and optimization of intelligent breeding in the future.

Methods

Data collection and processing

Statistical analysis

Published datasets of maize chromatin interactions and gene expression were used as experimental data in this study (Li et al., 2019; Peng et al., 2019). The data involved 2 tissues, shoots and ears. The chromatin interactions included 2 types, PPI and PDI (Li et al., 2019). The ear dataset (Li et al., 2019) contained 35 332 PPIs involving 20 601 genes, and the shoot (Li et al., 2019) dataset contained 65 064 PPIs involving 23 707 genes. We also used 3 additional datasets from shoots (Peng et al., 2019), ears (Sun et al., 2020), and tassels (Sun et al., 2020) (Supplemental Table 3) to evaluate the performance of DeepCBA in the prediction of gene expression classification and regression (Supplemental Figures 2 and 3).

To balance the number of genes in different categories, we classified genes into unexpressed, expressed, and highly expressed according to the expression ranges of [0–0.1), [0.1–15), and [15–maximum] (Supplemental Table 1). To reduce false positives in the training process, we divided the training and test datasets on the basis of gene family information. The numbers of genes with expression in the ranges of [0–100], [1–100], [0–500], and [0–maximum] can be seen in Supplemental Table 4. More than 99% of genes had expression in the range of [0–500], and we used these genes for the experimental analysis of expression prediction.

PDI data processing

The PDI datasets (Li et al., 2019) for shoots and ears included 11 207 and 11 189 PDIs, respectively. The intergenic distal sequences of different tissues were mainly in the range of 1–2 kb, and we set the length of the distal sequences to 1.5 kb. The promoter proximal sequence length was also set to 1.5 kb, including 1 kb upstream and 0.5 kb downstream of the gene TSS. When the length of the distal sequence was less than 1.5 kb, we added Ns at both ends of the sequence to 1.5 kb. When the length of the distal sequence was greater than 1.5 kb, 750-bp sequences were extracted from the middle of the sequence to both sides to form a 1.5-kb sequence.

Data augmentation

For the chromatin interactions of gene A and gene B in PPI mode, gene A may affect the expression of gene B, or gene B may affect the expression of gene A. We therefore doubled the number of PPI interactions by adjusting the order of gene sequences, a process termed data enrichment.

Quantitative real-time PCR

Quantitative real-time PCR data for ZmCLE7 were obtained from Liu et al. (2021a), (2021b). For quantitative real-time PCR analysis, total RNA was extracted from 0.1 g of plant tissue using the Quick RNA Isolation Kit (Huayueyang Biotechnology, Beijing, China). Analyzed leaf tissue was obtained from ZmVTE4-edited materials. EasyScript One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen Biotech, Beijing, China) was used to remove the genomic DNA from the extracted RNA and synthesize first-strand cDNA. Real-time fluorescence quantitative PCR with SYBR Green Master Mix (Vazyme Biotech, Nanjing, China) was performed on a CFX96 Real-Time System (Bio-Rad, Hercules, CA) to quantify gene expression of edited materials. The primers used for quantitative real-time PCR are listed in Supplemental Data 13 and 14. Concrete information on the target locations involved in ZmVTE4 gene editing is shown in Supplemental Figure 15.

Model elaboration

Model architecture

The DeepCBA model uses one-hot encoding and includes 3 modules: CNN, BiLSTM, and self-attention. The CNN module includes 3 parts, each of which contains 2 convolutional layers. Each convolutional layer connects to a maximum pooling layer to achieve feature dimensionality reduction and feature re-extraction. BiLSTM is used to capture feature information about the proximal and distal chromatin sequences and then mine important motif features that affect gene expression. The self-attention mechanism redistributes the weights of the parameters trained in the model to capture important features. Batch normalization and dropout mechanisms are used to reduce overfitting. In the last layer of DeepCBA, softmax or linear activation functions are used to perform the prediction task (Figure 1B; Supplemental Table 13).

The hyperparameters of DeepCBA are batch size, number of convolutional filters, and size of convolutional kernels. Experimental results showed that batch size has a large effect on the prediction results (Supplemental Figure 19A). The prediction accuracy of the model also improves as the number and size of convolutional kernels increase. By comprehensively considering the running time and memory space, we set the number and size of convolution kernels to 64 and 8, respectively (Supplemental Figure 19B and 19C; Supplemental Table 13).

We assumed that the sequence of the promoter proximal region of target gene G (Pseq_G) consisted of 2 parts: (1) 1 kb upstream and 0.5 kb downstream of the TSS and (2) 0.5 kb upstream and 1 kb downstream of the TTS. The expression level of G is Vag, and Pseq_A denotes the promoter proximal region sequence of gene_A that has a PPI interaction with G. The 4 steps described below—one-hot encoding, CNN layer, BiLSTM layer, and self-attention layer—are used to construct the prediction model for gene G based on PPI chromatin interactions.

One-hot encoding

One-hot encoding is used to process Pseq_G and Pseq_A with a length of 3 kb; that is, A = [1,0,0,0]T, C = [0,1,0,0]T, G = [0,0,1,0]T, T = [0,0,0,1]T, and N = [0,0,0,0]T. Then, Pseq_G is encoded as matrix M1∈R4 × 3000, and Pseq_A is encoded as matrix M2∈R4 × 3000. Then, M1 and M2 are concatenated vertically and input into the model P = Concat (M1, M2), P∈R4 × 6000.

CNN layer

In DeepCBA, the convolution operation in the CNN is used to perform dimensionality reduction and extract important features. For an encoded matrix M∈R(M∗N) and a filter W∈RU∗V, U≪M,V≪N, the convolution operation is shown in Equation 1.(Equation 1) Gij=∑u=1U∑v=1VWuvMi−u+1,j−v+1

For example, the dimension reduction operation on matrices M1 and M2 is shown in Equations 2 and3.(Equation 2) G1[m1,n1]=(M1∗W)[m1,n1]=∑j∑kW[j,k]M1[m1−j,n1−k]

(Equation 3) G2[m2,n2]=(M2∗W)[m2,n2]=∑j∑kW[j,k]M2[m2−j,n2−k]

W represents the filter, and mi and ni represent the number of rows and columns of the matrix after dimensionality reduction. We then obtain the reduced-dimension matrices G1[m1,n1] and G2[m2,n2].

The maximum pooling operation is used for secondary dimensionality reduction to solve the problem of overfitting. On the basis of the feature map Gi[mi,ni]∈RM∗N∗D obtained through the convolution operation, each feature map Gd∈RM∗N (1≤d≤D) can be divided into multiple regions Rm,nd (1≤m≤M′,1≤n≤N′). The maximum value within Rm,nd is then selected to represent the region, as shown in Equation 4.(Equation 4) Ym,nd=maxi∈Rm,ndGi

BiLSTM layer

For the promoter proximal sequence Pseq_G of gene G and Pseq_A, denoting the promoter proximal sequence of gene_A that has a PPI interaction with G, Ym1,m2 and Ym1,m2 are obtained through the CNN model. m1 = m2 = 3, n1 = n2 = 128. According to the self-loop update idea of the BiLSTM input gate, forget gate, and output gate, the update operation at time t (taking Ym1,n1 as an example) is shown in (Equation 5-1), (Equation 5-2), (Equation 5-3), (Equation 5-4), (Equation 5-5).(Equation 5-1) fi[m1,n1](t)=σ(bif+∑jUi,jfYj[m1,n1](t)+∑jWi,jfhj[m1,n1](t−1))

(Equation 5-2) Si[m1,n1](t)=fi[m1,n1](t)Si[m1,n1](t−1)gi[m1,n1](t)σ(bi+∑jUi,jYj[m1,n1](t)+∑jWi,jhj[m1,n1](t−1))

(Equation 5-3) gi[m1,n1](t)=σ(big+∑jUi,jgYj[m1,n1](t)+∑jWi,jghj[m1,n1](t−1))

(Equation 5-4) Oi[m1,n1](t)=σ(bio+∑jUi,joYj[m1,n1](t)+∑jWi,johj[m1,n1](t−1))

(Equation 5-5) hi[m1,n1](t)=tanh(Si[m1,n1](t)Oi[m1,n1](t))

where ht denotes the current hidden layer vector. i,j denote the i-th and j-th neurons, respectively, and ht includes the output of all LSTM cells. bf,Uf,Wf represent the bias value, input weight, and cycle weight of the corresponding threshold units, respectively. Oi[m1,n1](t) represents the output of the i-th neuron at the current time t. BiLSTM fully integrates the temporal and pre-/post-feature information of Ym1,n1 and Ym2,n2, reduces its dimensionality to 3∗64, and obtains Oi[m1,n1] and Oi[m2,n2]. m1 = m2 = 3, n1 = n2 = 64.

Self-attention layer

Om1,n1 and Om2,n2 denote the feature matrices of Pseq_G and Pseq_A, respectively. Self-attention mechanism composites Om1,n1 and Om2,n2 vertically and integrates the attention mechanism to enable redistribution of weights and predict target gene expression. By compositing Om1,n1 and Om2,n2 vertically, Os1,s2 is obtained. s1 = m1 + m2 = 6, s2 = n1 + n2 = 64. To improve the accuracy of feature extraction, an attention mechanism is used after the BiLSTM module to enable the redistribution of weights, as shown in Equation 6 (taking Os1,s2 as an example).(Equation 6) fis1,s2=tanh(WwOs1，s2+bw),ai=exp(biTbw)∑iexp(biTbw),Vs=∑iaifis1,s2

bi represents the implicit representation of feature Os1,s2 in the BiLSTM layer. The importance of feature Os1,s2 is measured by the similarity between bi and the sequence vector bw. Then, tanh is used to normalize the weight of each feature. Each feature is multiplied by its corresponding weight through the attention mechanism, and a summation is performed to obtain the output vector Vs. By setting dropout to prevent overfitting, it uses a linear function to obtain the expression Vag of the target gene G.

Transfer learning

Transfer learning is widely used to reduce model training time and achieve better results with a small amount of data. Here, we use transfer learning to fine-tune and train the DeepCBA model, thus enabling the prediction of gene expression across tissues and genotypes (Supplemental Figures 16 and 17). Specifically, we freeze the first 0–19 layers of the pre-trained model and perform fine-tuning through small-batch training (learning rate = 1e−4, batch size = 64, epochs = 200).

DeepCBA mines important DNA sequence regions and motifs

To mine important motifs that affect gene expression in the PPI and PDI modes, a saliency map (Chu, 2011) is used to calculate the significance map (true positive, TP; true negative, TN), and it generates the model gradient through the reverse calculation. The gradient value is denoted as 4 × 6000 × N, where N represents the sum of the number of TPs and TNs. Finally, TF-MoDIsco (Shrikumar et al., 2018) is used to mine important motifs (Supplemental Figure 4B and 4C). The important identified motifs that affect gene expression are compared with the TOMTOM (Cheng et al., 2019, Gupta et al., 2007) platform and 653 conserved structural regions of higher plants in the PlantTFDB (Supplemental Tables 11 and 12) to evaluate the reliability of the motifs.

Different motif compositions affect gene expression

Background sequences are generated by filling a sequence with length 3 kb based on the N coding mode. The important motifs identified by DeepCBA are combined to form a motif composition (Real motif) and A/C/G/T are randomly combined to form a motif composition of equal length (Random motif). The real motif and random motif are embedded into the background sequence, and gene expression is predicted using the DeepCBA model (Supplemental Figure 7B).

Motif enrichment in open chromatin regions

To mine the characteristics of important motifs identified by DeepCBA in shoots and ears, the true physical location of each motif on the chromosome is identified. The physical locations of the identified motifs are then matched to open chromatin regions in the NAM population. The operation is repeated 100 times to establish the following 2 controls: (1) motif sequences removed from the PPI and PDI sequences and sequences of equal length selected from the remaining PPI and PDI sequences, and (2) PPI and PDI sequences removed from the whole genome and sequences of equal length selected from the remaining genome sequences.

Data and code availability

Methods, including statements of data availability and any associated accession codes and references, are available at https://github.com/Jie-Lii/DeepCBA.

Funding

This work has been supported by the Biological Breeding-Major Projects (2023ZD04076), the 10.13039/501100012166 National Key Research and Development Program of China (2022YFD1201504 ), the 10.13039/501100012226 Fundamental Research Funds for the Central Universities (2662022YLYJ010 , 2021ZKPY018 , 2662021JC008 , and SZYJY2021003 ), the Major Project of Hubei Hongshan Laboratory (2022HSZD031 ), the Major Science and Technology Project of Hubei Province (2021AFB002 ), and the Yingzi Tech & Huazhong Agricultural University Intelligent Research Institute of Food Health (IRIFH202209 ).

Author contributions

Conceptualization, J.Y. and J. Liu. Methodology, Z.W. and J. Li. Writing – original draft, Z.W. Writing – review & editing, Z.W., Y.P., and J. Liu. Funding acquisition, J. Liu.

Supplemental information

Document S1. Supplemental Figures 1‒19, Supplemental Tables 1‒13, and Algorithm 1

Supplemental Data 1. The matching results between the motifs mined by ear-1 and 104TF in PPI mode

Supplemental Data 2. The matching results between the motifs mined by shoot-1 and 104TF in PPI mode

Supplemental Data 3. The matching results between the motifs mined by ear-1 and 104TF in PDI mode

Supplemental Data 4. The matching results between the motifs mined by shoot-1 and 104TF in PDI mode

Supplemental Data 5. Motif information mined by ear-1 in PPI mode

Supplemental Data 6. Motif information mined by shoot-1 in PPI mode

Supplemental Data 7. Motif information mined by ear-1 in PDI mode

Supplemental Data 8. Motif information mined by shoot-1 in PDI mode

Supplemental Data 9. Expanding the upstream 3k sequence of a specific gene to mine regulatory element results

Supplemental Data 10. vgt1 validation

Supplemental Data 11. The predicted results of DeepCBA in different datasets under PPI mode

Supplemental Data 12. The predicted results of DeepCBA in different datasets under PDI mode

Supplemental Data 13. Primer list for real-time quantitative PCR of gene CLE7

Supplemental Data 14. Primer list for real-time quantitative PCR of gene VTE4

Supplemental Data 15. Data used for plotting all the figures in the main text

Supplemental Data 16. Data used for plotting all the figures in the supplemental file

Document S2. Article plus supplemental information

Acknowledgments

We thank the anonymous reviewers for their constructive comments, which helped us to substantially improve the manuscript. We thank the Experimental Teaching Center of the College of Informatics at Huazhong Agricultural University for providing the experimental environment and computational 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

Avsec Ž. Agarwal V. Visentin D. Ledsam J.R. Grabska-Barwinska A. Taylor K.R. Assael Y. Jumper J. Kohli P. Kelley D.R. Effective gene expression prediction from sequence by integrating long-range interactions Nat. Methods 18 2021 1196 1203 34608324
Beer M.A. Tavazoie S. Predicting gene expression from sequence Cell 117 2004 185 198 15084257
Bailey T.L. Johnson J. Grant C.E. Noble W.S. The MEME suite Nucleic Acids Res. 43 2015 W39 W49 25953851
Cheng A. Grant C.E. Noble W.S. Bailey T.L. MoMo: discovery of statistically significant post-translational modification motifs Bioinformatics 35 2019 2774 2782 30596994
Cheng C. Yan K.K. Yip K.Y. Rozowsky J. Alexander R. Shou C. Gerstein M. A statistical framework for modeling gene expression using chromatin features and application to modENCODE datasets Genome Biol. 12 2011 R15 R18 21324173
Chu C. Saliency mapping of figure and ground of motion in Chinese J. Chin. Lang. Teach. Assoc. 46 2011 49 69
Chen Y. He Z. Men Y. Dong G. Hu S. Ying X. MetaLogo: a heterogeneity-aware sequence logo generator and aligner Brief. Bioinform. 23 2022 bbab591 35108357
Cao X. Costa L.M. Biderre-Petit C. Kbhaya B. Dey N. Perez P. McCarty D.R. Gutierrez-Marcos J.F. Becraft P.W. Abscisic acid and stress signals induce Viviparous1 expression in seed and vegetative tissues of maize Plant Physiol. 143 2007 720 731 17208960
Dong X. Greven M.C. Kundaje A. Djebali S. Brown J.B. Cheng C. Gingeras T.R. Gerstein M. Guigó R. Birney E. Modeling gene expression using chromatin features in various cellular contexts Genome Biol. 13 2012 R53 22950368
Fu J. Cheng Y. Linghu J. Yang X. Kang L. Zhang Z. Zhang J. He C. Du X. Peng Z. RNA sequencing reveals the complex regulatory network in the maize kernel Nat. Commun. 4 2013 2832 24343161
Guerriero G. Martin N. Golovko A. Sundström J.F. Rask L. Ezcurra I. The RY/Sph element mediates transcriptional repression of maturation genes from late maturation to early seedling growth New Phytol. 184 2009 552 565 19659659
Grant C.E. Bailey T.L. Noble W.S. FIMO: scanning for occurrences of a given motif Bioinformatics 27 2011 1017 1018 21330290
Gaspin C. Rami J.F. Lescure B. Distribution of short interstitial telomere motifs in two plant genomes: putative origin and function BMC Plant Biol. 10 2010 283 312 21171996
Gupta S. Bailey T.L. Noble W.S. Quantifying similarity between motifs Genome Biol. 8 2007 1 9
Hufford M.B. Seetharam A.S. Woodhouse M.R. Chougule K.M. Ou S. Liu J. Ricci W.A. Guo T. Olson A. Qiu Y. De novo assembly, annotation, and comparative analysis of 26 diverse maize genomes Science 373 2021 655 662 34353948
Ishige F. Takaichi M. Foster R. Chua N. Oeda K. AG-box motif (GCCACGTGCC) tetramer confers high-level constitutive expression in dicot and monocot plants Plant J. 18 1999 443 448
Jin J. Tian F. Yang D.C. PlantTFDB 4.0: toward a central hub for transcription factors and regulatory interactions in plants Nucleic Acids Res. 2016 gkw982
Jin J. He K. Tang X. Li Z. Lv L. Zhao Y. Luo J. Gao G. An Arabidopsis transcriptional regulatory map reveals distinct functional and evolutionary features of novel transcription factors Mol. Biol. Evol. 32 2015 1767 1773 25750178
Jin J. Zhang H. Kong L. Gao G. Luo J. PlantTFDB 3.0: a portal for the functional and evolutionary study of plant transcription factors Nucleic Acids Res. 42 2014 D1182 D1187 24174544
Jores T. Tonnies J. Wrightsman T. Synthetic promoter designs enabled by a comprehensive analysis of plant core promoters Nat. Plants 7 2024 842 855 10.1101/2021.01.07.425784
Karlić R. Chung H.R. Lasserre J. Vlahovicek K. Vingron M. Histone modification levels are predictive for gene expression Proc. Natl. Acad. Sci. USA 107 2010 2926 2931 20133639
Kelley D.R. Reshef Y.A. Bileschi M. Belanger D. McLean C.Y. Snoek J. Sequential regulatory activity prediction across chromosomes with convolutional neural networks Genome Res. 28 2018 739 750 29588361
Lee D. Yang J. Kim S. Learning the histone codes with large genomic windows and three-dimensional chromatin interactions using transformer Nat. Commun. 13 2022 6678 36335101
Li E. Liu H. Huang L. Zhang X. Dong X. Song W. Zhao H. Lai J. Long-range interactions between proximal and distal regulatory regions in maize Nat. Commun. 10 2019 2633 31201330
Liu L. Zhang G. He S. Hu X. TSPTFBS: a docker image for trans-species prediction of transcription factor binding sites in plants Bioinformatics 37 2021 260 262 33416862
Liu L. Gallagher J. Arevalo E.D. Chen R. Skopelitis T. Wu Q. Bartlett M. Jackson D. Enhancing grain-yield-related traits by CRISPR–Cas9 promoter editing of maize CLE genes Nat. Plants 7 2021 287 294 33619356
Mönke G. Altschmied L. Tewes A. Reidt W. Mock H.P. Bäumlein H. Conrad U. Seed-specific transcription factors ABI3 and FUS3: molecular interaction with DNA Planta 219 2004 158 166 14767767
O’Connor T. Grant C.E. Bodén M. Bailey T.L. T-Gene: improved target gene prediction Bioinformatics 36 2020 3902 3904 32246829
Oka R. Zicola J. Weber B. Anderson S.N. Hodgman C. Gent J.I. Wesselink J.J. Springer N.M. Hoefsloot H.C.J. Turck F. Genome-wide mapping of transcriptional enhancer candidates using DNA and chromatin features in maize Genome Biol. 18 2017 137 224 28732548
Peng Y. Xiong D. Zhao L. Ouyang W. Wang S. Sun J. Zhang Q. Guan P. Xie L. Li W. Chromatin interaction maps reveal genetic regulation for quantitative traits in maize Nat. Commun. 10 2019 2632 31201335
Ricci W.A. Lu Z. Ji L. Marand A.P. Ethridge C.L. Murphy N.G. Noshay J.M. Galli M. Mejía-Guerra M.K. Colomé-Tatché M. Widespread long-range cis-regulatory elements in the maize genome Nat. Plants 5 2019 1237 1249 31740773
Rodríguez-Leal D. Lemmon Z.H. Man J. Bartlett M.E. Lippman Z.B. Engineering quantitative trait variation for crop improvement by genome editing Cell 171 2017 470 480.e8 28919077
Rodgers-Melnick E. Vera D.L. Bass H.W. Buckler E.S. Open chromatin reveals the functional maize genome Proc. Natl. Acad. Sci. USA 113 2016 E3177 E3184 27185945
Reidt W. Wohlfarth T. Ellerström M. Czihal A. Tewes A. Ezcurra I. Rask L. Bäumlein H. Gene regulation during late embryogenesis: the RY motif of maturation-specific gene promoters is a direct target of the FUS3 gene product Plant J. 21 2000 401 408 10758492
Schmidt F. Gasparoni N. Gasparoni G. Gianmoena K. Cadenas C. Polansky J.K. Ebert P. Nordström K. Barann M. Sinha A. Combining transcription factor binding affinities with open-chromatin data for accurate gene expression prediction Nucleic Acids Res. 45 2017 54 66 27899623
Schoenfelder S. Fraser P. Long-range enhancer–promoter contacts in gene expression control Nat. Rev. Genet. 20 2019 437 455 31086298
Song X. Meng X. Guo H. Cheng Q. Jing Y. Chen M. Liu G. Wang B. Wang Y. Li J. Targeting a gene regulatory element enhances rice grain yield by decoupling panicle number and size Nat. Biotechnol. 40 2022 1403 1411 35449414
Su W. Shao Z. Wang M. EjBZR1 represses fruit enlargement by binding to the EjCYP90 promoter in loquat Hortic. Res. 8 2021 152 34193858
Sun Y. Dong L. Zhang Y. 3D genome architecture coordinates trans and cis regulation of differentially expressed ear and tassel genes in maize Genome Biol 21 2020 1 25
Shrikumar A. Tian K. Avsec Ž. Technical note on transcription factor motif discovery from importance scores (TF-MoDISco) version 0.5. 6.5 Preprint at arXiv:1811.00416 2018 10.48550/arXiv.1811.00416
Tian F. Yang D.C. Meng Y.Q. Jin J. Gao G. PlantRegMap: charting functional regulatory maps in plants Nucleic Acids Res. 48 2020 D1104 D1113 31701126
Tu X. Mejía-Guerra M.K. Valdes Franco J.A. Tzeng D. Chu P.Y. Shen W. Wei Y. Dai X. Li P. Buckler E.S. Reconstructing the maize leaf regulatory network using ChIP-seq data of 104 transcription factors Nat. Commun. 11 2020 5089 33037196
Tasaki S. Gaiteri C. Mostafavi S. Wang Y. Deep learning decodes the principles of differential gene expression Nat. Mach. Intell. 2 2020 376 386 32671330
Tian T. Wang S. Yang S. Yang Z. Liu S. Wang Y. Gao H. Zhang S. Yang X. Jiang C. Genome assembly and genetic dissection of a prominent drought-resistant maize germplasm Nat. Genet. 55 2023 496 506 36806841
Washburn J.D. Mejia-Guerra M.K. Ramstein G. Kremling K.A. Valluru R. Buckler E.S. Wang H. Evolutionarily informed deep learning methods for predicting relative transcript abundance from DNA sequence Proc. Natl. Acad. Sci. USA 116 2019 5542 5549 30842277
Wu J. Deng Y. Hu J. Jin C. Zhu X. Li D. Genome-wide analyses of direct target genes of an ERF11 transcription factor involved in plant defense against bacterial pathogens Biochem. Biophys. Res. Commun. 532 2020 76 81 32828541
Woodhouse M.R. Cannon E.K. Portwood J.L. Harper L.C. Gardiner J.M. Schaeffer M.L. Andorf C.M. A pan-genomic approach to genome databases using maize as a model system BMC Plant Biol. 21 2021 385 410 34416864
Xu W. Yang R. Li M. Transcriptome phase distribution analysis reveals diurnal regulated biological processes and key pathways in rice flag leaves and seedling leaves Plos One 6 2011 e17613
Yang T. Guo L. Ji C. The B3 domain-containing transcription factor ZmABI19 coordinates expression of key factors required for maize seed development and grain filling The Plant Cell 33 2021 104 128 33751093
Yang N. Liu J. Gao Q. Genome assembly of a tropical maize inbred line provides insights into structural variation and crop improvement Nature Genetics 51 2019 1052 1059 31152161
Zhao H. Tu Z. Liu Y. Zong Z. Li J. Liu H. Xiong F. Zhan J. Hu X. Xie W. PlantDeepSEA, a deep learning-based web service to predict the regulatory effects of genomic variants in plants Nucleic Acids Res. 49 2021 W523 W529 34037796
Zhou J. Theesfeld C.L. Yao K. Chen K.M. Wong A.K. Troyanskaya O.G. Deep learning sequence-based ab initio prediction of variant effects on expression and disease risk Nat. Genet. 50 2018 1171 1179 30013180
Zhao H. Zhang W. Chen L. Wang L. Marand A.P. Wu Y. Jiang J. Proliferation of regulatory DNA elements derived from transposable elements in the maize genome Plant Physiol. 176 2018 2789 2803 29463772
Zrimec J. Börlin C.S. Buric F. Muhammad A.S. Chen R. Siewers V. Verendel V. Nielsen J. Töpel M. Zelezniak A. Deep learning suggests that gene expression is encoded in all parts of a co-evolving interacting gene regulatory structure Nat. Commun. 11 2020 6141 33262328
Zrimec J. Börlin C.S. Buric F. Deep learning suggests that gene expression is encoded in all parts of a co-evolving interacting gene regulatory structure Nat. Commun. 11 2020 6141 33262328
Zhou J. Troyanskaya O.G. Predicting effects of noncoding variants with deep learning–based sequence model Nat. Methods 12 2015 931 934 26301843
