Clarion is a multi-label problem transformation method for identifying mRNA subcellular localizations

Benchmark dataset

In this study, all subcellular localization annotations and mRNA sequences were collected from the RNALocate database (version 2.0) [32]. The latest version of RNALocate (updated in June 2021) contained more than 210 000 RNA subcellular localization entries, encompassing more than 110 000 RNAs with 171 subcellular localizations across 104 different species. Its version 2.0 provides more accurate localization annotations than the first version, facilitating the construction of a reliable benchmark dataset. More specifically, the benchmark dataset was constructed according to the five following major steps:

• 1) We downloaded all RNA subcellular localization annotation entries from RNALocate (version 2.0) and accordingly collected 84 792 mRNA subcellular localization entries as the initial dataset.

• 2) According to the statistics of the initial dataset, there were 150 different types of annotated subcellular localizations. However, some had minimal and incomplete entries and as such, we removed those subcellular localization types whose corresponding entry numbers were less than 3000. As a result, we obtained nine types of subcellular localizations with 152 887 unique transcripts including exosome, nucleus, cytosol, chromatin, nucleoplasm, ribosome, nucleolus, cytoplasm and membrane.

• 3) Next, we redefined the mapping relationships between mRNAs and subcellular localizations based on multiple localizations of mRNAs in the transcriptome. In particular, an mRNA can be labelled with multiple subcellular localizations instead of being only labelled with one subcellular localization.

• 4) To reduce the effect of sequence redundancy on the performance of the classifier, we applied CD-HIT-EST [33] to remove the redundant sequences with the 80% sequence identity threshold to ensure the similarity between any two nucleotide sequences was less than 80%. Finally, 36 971 mRNAs were obtained and used as the benchmark dataset.

• 5) We analyzed the distribution of sequence length of these 36 971 mRNAs, which varied from 119 nt to 12 000 nt. In view of the computing complexity and limitation of the feature engineering algorithms, the mRNA lengths were adjusted to no longer than 6000 nt. Specifically, for those mRNAs with more than 6000 nt, the first 3000 nt and the last 3000 nt were extracted and merged.

Sequence vectorization

mRNA sequences need to be encoded as numeric vectors prior to the training of machine learning models. The -mer nucleotide composition is one of the widely used sequence encoding methods, which has been successfully applied in a variety of bioinformatics studies [34–38]. Given an mRNA sequence with length nt, , where and represents the nucleotide acid at position , one type of adenine (A), cytosine (C), guanine (G) and uracil/thymine (U/T). Accordingly, when using the -mer () method to encode features, the feature vectors can be calculated as:

where represents the number of -mer type along the sequence. From the equation, we can observe that the -mer vector dimension increases exponentially with the increase of value. These numerous features may contain a great deal of redundancy and noise, which may cause extra training time and even have a negative influence on the model quality. In view of the dimension restriction and training effectiveness, we used 1-mer, 2-mer, 3-mer, 4-mer, 5-mer and 6-mer in this study.

Weighted series

In this study, we proposed a novel problem transformation method named weighted series (WS) to tackle multi-label learning problems, which is specific for the subcellular localization identification of mRNAs. The weighted series method involves two modules of binary classifiers, including a non-label module and a fusion-label module. The non-label module is concerned with training the model only from pre-extracted features. While the fusion-label module incorporates the priori information about the labels in the model training, whose learned priori label distributions could contribute to the model predictions. The final predictions of the two modules are combined by a weight (ranging from 0 to 1) that requires user customization, reflecting the labels' relevance. When there is a strong correlation between the labels, a smaller value of will promote the prediction performance.

Suppose represents the dimensional instance space and represents the dimensional label space. Given a multi-label data set that containing samples, where and . Both non-label module and fusion-label module require training q binary classifiers. When training the -th model () of non-label module, the training label is , and the training input is the feature vector of the training samples, i.e., . With regard to training the -th model () of the fusion-label module, the priori information on other labels is added to the training process by combining features as input. That is to say, the label is still , whereas the training input is the fusion of and . After model training, 2q binary classifiers of weighted series, and , can be used to conduct predictions. The prediction process of weighted series includes three steps: module prediction, module prediction and integration. Given a query mRNA sequence and its feature vector are , firstly, non-label models are used to conduct the prediction and generate the prediction probabilities ; would then be fused with as the input to conduct prediction of fusion-label models and output the probabilities ; finally, and would be integrated with a user-defined weight that ranges from 0 and 1. The detailed training and prediction procedures of the weighted series are illustrated in Figure 1 and outlined in Algorithm 1.

Performance evaluation metrics

The model evaluation for MLL tasks is more complicated than binary classification problems because the predictive performance for all labels should be taken into account. In this study, we employed six widely used MLL evaluation metrics [29, 39, 40] to evaluate the performance of Clarion, including example-based accuracy (), average precision, coverage, one-error, ranking loss and Hamming loss. Let be the learned multi-label classifier and be a multi-label instance, and represent the true and predicted label set for the instance, represents the complementary set of , accordingly the above metrics can be formulated as follows:

where represents the rank of y in Y based on the descending order, represents the cardinality of set while q is the cardinality of , ∆ stands for the symmetric difference of two sets, while counts the times that meet the condition.

Statistical analysis of the dataset

The benchmark dataset used in this study contained a total of 36 971 mRNA sequences, each of which might be localized to a single or multiple subcellular compartments. 12 884 mRNAs had only one localization/compartment, 4060 mRNAs had two localizations, 3442 had three localizations, 3165 mRNAs had four localizations, 3518 mRNAs had five localizations, 4258 mRNAs had six localizations, 4079 mRNAs had seven localizations, 1443 mRNAs had eight localizations and 122 mRNAs had nine localizations, as illustrated in Figure 2A. In addition, we also plotted the distribution of the positive and negative samples in the nine compartments, as shown in Figure 2B. Compared with the nucleus, nucleoplasm, chromatin, nucleolus, cytosol, membrane and ribosome, the data of exosome and cytoplasm were unevenly distributed. For the exosome, the number of the positive samples was much larger than that of the negatives, while for cytoplasm, the number of the negative samples was much larger than that of the positives. In addition, there were 31 448 mRNAs in exosome, 21 439 in nucleus, 14 237 in nucleoplasm, 14 328 in chromatin, 4016 in cytoplasm, 11 124 in nucleolus, 16 312 in cytosol, 6739 in membrane and 8680 in ribosome, respectively. The distribution of the label (i.e., subcellular localization) number in these nine compartments is shown in Figure 2C.

We noticed that those mRNAs in the cytoplasm were mostly single-localized, mRNAs in the nucleoplasm, chromatin, nucleolus, cytosol, membrane and ribosome were mostly multi-localized, while mRNAs in exosome and nucleus appeared to be uniformly localized. In this study, we separated the benchmark dataset into the training_validation dataset and the independent test dataset by random sampling. Accordingly, 33,274 mRNAs were included in the training_validation dataset (90% of the total) and 3697 mRNAs in the independent test set (10% of the total). The former was used to compare the algorithm and optimize the parameters by tenfold cross-validation, while the latter was used to evaluate and validate the model performance. Detailed localization distributions of the training_validation and independent test datasets can be found in the Supplementary Table S1.

Selection of the base classifier

In this section, we performed a preliminary analysis to select a suitable base classifier. To do this, we evaluated seven popular machine learning algorithms, including k-nearest neighbour (KNN), logistic regression (LR), random forest (RF), LightGBM, XGBoost, CatBoost and multilayer perception (MLP) [41, 42] based on the WS framework to determine the optimal base classifier. For the sake of comparison, the weight of WS was assigned to one-quarter, two-quarters and three-quarters of its range, respectively. This can save the running time of the algorithm and accelerate the determination of the base classifier. We employed the default parameters of each algorithm in the Scikit-learn Package and conducted a 10-fold cross-validation test on the training_validation dataset for performance comparison. As can be seen from Figure 3A, in the case of w = 0.25, 0.5 and 0.75, XGBoost secured the best predictive performance among the seven different algorithms in terms of accuracy. In the case of w = 0.25, we found XGBoost achieved a slightly inferior performance to RF in terms of coverage and ranking loss but attained the best performance in terms of the other metrics (refer to Supplementary Table S2 for detailed results). In addition, XGBoost was still the best-performing algorithm for w = 0.5 and 0.75 in terms of all the evaluation metrics. Consequently, the XGBoost algorithm was adopted as the base classifier of the weighted series for developing Clarion. The relationship between the model performance and the weight values will be discussed in more detail in the following section.

Effect of weight w

In this section, we evaluated the effect of weight in weighted series structure on the model performance, which ranged from 0 to 1. In particular, we evaluated the model performance on the training data through a ten-fold cross-validation test with 19 different candidates ranging from 0.05 to 0.95 with a step size of 0.05. Figure 3B illustrates the model performance with different candidate weights in terms of all six performance metrics. There is a clear peak/valley in each metric sub-figure, whereas the corresponding candidate weights are not consistent. As it is much more difficult to directly determine the first-rank value of weighted series for Clarion, we employed a method by assigning the weights score. According to the performance ranking of each evaluation metric, we assigned 5, 4, 3, 2 and 1 points to the top five candidate weights and 0 to those after the sixth, respectively. For instance, the candidate was assigned 2 points on Acc_exam as the accuracy was the fourth highest out of the 19 candidates. Similarly, the candidate was assigned 5 points for the average precision, 0 points for the coverage, 5 points for the One-error, 0 points for the ranking loss and 4 points for the hamming loss, respectively. We then obtained an overall score of 16 points for the candidate by summing up all the above scores. Using this procedure, we calculated the overall scores of the other 18 candidate weights, whose detailed results and score statistics are provided in Supplementary Table S3 and Supplementary Table S4. The candidate reached the best score of 22 points out of the 19 candidates and accordingly, it was adopted as the fixed weight of weighted series for Clarion.

Performance comparison with other problem transformation strategies

To demonstrate the capacity of the weighted series strategy in dealing with the multi-label multi-class mRNA subcellular localization prediction tasks, we used XGBoost with the default parameters as the base classifier to benchmark and compare our proposed method with the other three well-known problem transformation methods, including BR, CC and LP, on the 10-fold cross-validation tests using the training_validation set. Specifically, we performed 10 times of experiments with 10 groups of randomly generated label orders for CC because the input label order could directly affect the quality of the CC model. Among these 10 experiments, the best result was used for the comparison with BR, LP and WS. From the performance comparison results provided in Table 2, we found that the WS strategy displayed the best performance in terms of the evaluation metrics of Acc_exam, average precision, coverage, ranking loss and hamming loss. Although WS achieved a slightly lower performance than LP in terms of one-error, it showed a clear superiority in predicting mRNA subcellular localizations. To further improve the model performance, we optimized the key parameters, including learning_rate, n_estimators and max_depth, for XGBoost based on the 10-fold cross-validation test. The specific hyperparameters are provided in Supplementary Table S5. Subsequently, we retrained the final model named Clarion on the whole training_validation set using the weighted series strategy with the optimized hyperparameters of XGBoost.

Table 2.

Performance comparison of weighted series with binary relevance, classifier chains and label powerset

Strategy	Accexam	Average precision	Coverage	One-error	Ranking loss	Hamming loss
BR	0.600	0.651	6.121	0.706	0.345	0.194
CC	0.456	0.580	7.384	0.838	0.598	0.299
LP	0.558	0.626	6.271	0.601	0.443	0.229
WS ()	0.627	0.670	6.029	0.629	0.344	0.182

Open in a new tab

Performance comparison with existing state-of-the-art tools

In this section, Clarion’s performance was benchmarked and compared with several state-of-the-art approaches for predicting mRNA subcellular localizations. Firstly, we compared Clarion with DM3Loc and Wang’s method, the only two multi-label predictors. Clarion has only five overlapping predictable compartments with these two predictors, including cytosol, exosome, membrane (cytoplasm for Wang’s method), nucleus and ribosome. Therefore, we compared their five-label prediction performance via an independent test. Clarion and DM3Loc achieved the Acc_exam of 0.722/0.441, average precision of 0.769/0.618, coverage of 3.019/3.957, one-error of 0.463/0.869, ranking loss of 0.204/0.533 and hamming loss of 0.150/0.330 on the independent dataset. Similarly, Clarion and Wang’s method achieved Acc_exam of 0.745/0.281, average precision of 0.767/0.679, coverage of 3.127/4.516, one-error of 0.445/0.882, ranking loss of 0.241/0.716 and hamming loss of 0.146/0.375. The above results indicated that Clarion outperformed DM3Loc and Wang’s method in multi-label prediction tasks.

Afterwards, we also compared Clarion’s single-label prediction performance with that of the other state-of-the-art methods. To facilitate the performance comparison, only the methods with accessible webservers were used for prediction and comparison, including iLoc-mRNA [19], mRNALoc [20], mRNALocator [21], DM3Loc [25] and Wang’s method [26]. In particular, the RNA sequences of the independent set were uploaded to their webservers, which then outputted the prediction labels. As shown in Table 3, Clarion clearly outperformed the other methods in predicting cytoplasm, cytosol, exosome, membrane, nucleus and ribosome. We also found that Clarion secured over 80% accuracies in almost all compartment predictions with the only exception of cytosol and nucleus. Notably, the mRNAs of the independent set may appear in the training set of other methods, which may account for Clarion’s slightly lower F1 scores than Wang’s method on the prediction of cytoplasm and ribosome (more details can be found in Supplementary Table S6 of the supplementary file). These comparison results further demonstrated Clarion’s superior prediction power on single-label tasks.

Effect of non-label and fusion-label modules in the weighted series structure

In this section, we further examined the effect of the non-label and fusion-label modules used in the WS framework. A total of 18 binary classifiers were trained in Clarion for predicting nine mRNA subcellular localizations respectively, including nine non-label models and nine fusion-label models. Here, we employed these binary models to predict the RNA sequences in the independent test dataset and accordingly evaluated the performance for each localization. Figure 3C illustrates the predictive performance of non-label models and fusion-label models in terms of accuracy. As a result, we found that the fusion-label models performed better than the non-label counterparts for all nine locations, highlighting the necessity and effectiveness of fusion-label models in the WS framework. However, it is noteworthy that the better performance of the fusion-label models originated from the outputs of non-label models because the fusion-label models used the prediction results of the non-label models as input features for the model training. Therefore, we conclude that the non-label and fusion-label modules are complementary and essential to the WS framework of Clarion.

Model interpretation

Shapley additive explanations (SHAP) is a powerful method based on the cooperative game theory that can interpret machine learning models [31] using the Shapley value, which can be used to rank and evaluate the importance of each feature and explain the predictions. SHAP has been successfully applied in a variety of bioinformatics tasks [43–46]. In this study, we used the Shapley value to assess the -mer fragments of mRNAs that are important for subcellular localization prediction. With the nine binary models of the non-label module, the Shapley value of each -mer feature was calculated and ranked using the SHAP Python package (https://shap.readthedocs.io/en/latest/index.html). Figure 4 and Supplementary Figures S1-S3 show the top 15 important -mer features of Clarion for predicting mRNA subcellular localizations. We found some features are only important for one certain localization, such as ‘TTG’, ‘GGGCGC’ in the nucleus, ‘GACGC’ and ‘GCGGCA’ in chromatin, and ‘GGATCT’ and ‘GGCCG’ in the ribosome. For example, ‘TTG’ was identified as the most important feature for nucleus prediction in terms of the SHAP value but did not appear in the list of the top 15 features for the prediction of the other eight localizations. For ‘TTG’ shown in Figure 4A, each point represents an instance with a value from small to large corresponding to the colour from blue to red. When ‘TTG’ took high Shapley values, it would have an influence on the model to make the positive prediction of nucleus and visa verse. Therefore, the effect of ‘TTG’ can be summarized in a way that its larger value promotes the nucleus localization prediction. In contrast, the larger value of ‘GGGCGC’ promotes the non-nucleus localization prediction. These -mer segments may be part of or related to protein recognition motifs for mRNA specific localization. Interestingly, the repeat motif ‘UUCAC’ has been found to be crucial for localization via binding with Vg1PBP [47–49], which corresponded to ‘TTCACC’ that was ranked thirteenth in exosome prediction.

In addition, it was also found that certain -mer features were important for the prediction of several subcellular localizations. For instance, ‘GCGGC’ was ranked second in nucleus prediction, tenth in chromatin prediction, and tenth in ribosome prediction, respectively, as shown in Figure 4A–C. This could be also observed in Figure 4D, which shows the top 15 features ranked according to the average of the absolute SHAP values and shows the important proportion in different mRNA subcellular localizations. The feature ‘GCGGC’ shown in Figure 4D is more important for the prediction of the nucleus, ribosome, chromatin, as well as nucleoplasm. Moreover, the feature ‘AAAAAA’ was important for almost all compartments of subcellular localizations, whose larger values drove the positive prediction of cytoplasm, but the negative predictions of exosome, chromatin, nucleoplasm, nucleolus, cytosol and ribosome.

Webserver implementation

To facilitate the wider research community to make subcellular localization predictions, we developed a web server for Clarion that is freely available at http://monash.bioweb.cloud.edu.au/Clarion/. The web server is maintained by Nectar Research Cloud and configured on a Linux server equipped with a 4-core CPU, 8-GB memory and 30-GB hard disk. The web page was implemented using PHP and has been tested on several popular web browsers, including Google Chrome, Microsoft Edge, Internet Explorer, Mozilla Firefox and Safari. Users are required to copy and paste query mRNA sequences in the textbox or alternatively upload a sequence file in the FASTA format via the file-selection dialogue box. It should be noted the query sequences will be truncated to ensure that no mRNA is greater than 6000 nt. All the generated prediction results will be saved in a table format containing detailed information regarding the sequences and predicted subcellular localization types. The web server also provides a probability score in the range of 0–1 to indicate the probability of the subcellular localization type. The prediction results can be easily exported to widely used file formats, including CSV, Excel, PDF and plain text. More detailed instructions for using the Clarion webserver can be found on the help page of the webserver.

Fuyi Li received his PhD in Bioinformatics from Monash University, Australia. He is currently a professor in the College of Information Engineering, Northwest A&F University, China. His research interests are bioinformatics, computational biology, machine learning, and data mining.

Xudong Guo received his MEng degree from Ningxia University, China. He is currently a research assistant at the College of Information Engineering, Northwest A&F University, China. His research interests are bioinformatics and data mining.

Zhikang Wang is currently a PhD student in the Department of Biochemistry and Molecular Biology and Biomedicine Discovery Institute, Monash University, Australia. His interests are bioinformatics, computational pathology, pattern recognition and deep learning.

Tong Pan is currently a PhD student in the Department of Biochemistry and Molecular Biology and Biomedicine Discovery Institute, Monash University, Australia. Her interests are bioinformatics, protein function analysis, deep learning and pattern recognition.

Yuming Guo is a professor of Global Environmental Health and Biostatistics & Head of the Monash Climate, Air Quality Research (CARE) Unit. His research focuses on environmental epidemiology, biostatistics, global environmental change, air pollution, climate change, urban design, residential environment, remote sensing modelling and infectious disease modelling.

Geoffrey I. Webb is a professor in the Faculty of Information Technology and a research director of the Monash Data Futures Institute at Monash University. His research interests include machine learning, data mining, computational biology and user modelling.

Jianhua Yao is a group leader in Tencent AI Lab, China. He received his PhD degree from John Hopkins University. His research interests include bioinformatics, computational biology, medical imaging and pattern recognition.

Cangzhi Jia is a professor in the college of science, Dalian Maritime University, China. She obtained her PhD degree in the school of Mathematical Sciences from the Dalian University of Technology in 2007. Her major research interests include mathematical modelling in bioinformatics and machine learning.

Jiangning Song is an associate professor and a group leader in the Monash Biomedicine Discovery Institute, Monash University. He is also affiliated with the Monash Data Futures Institute, Monash University. His research interests include bioinformatics, computational biomedicine, machine learning, data mining and pattern recognition.