Multi-task learning for the simultaneous reconstruction of the human and mouse gene regulatory networks

The reconstruction of Gene Regulatory Networks (GRNs) from gene expression data, supported by machine learning approaches, has received increasing attention in recent years. The task at hand is to identify regulatory links between genes in a network. However, existing methods often suffer when the number of labeled examples is low or when no negative examples are available. In this paper we propose a multi-task method that is able to simultaneously reconstruct the human and the mouse GRNs using the similarities between the two. This is done by exploiting, in a transfer learning approach, possible dependencies that may exist among them. Simultaneously, we solve the issues arising from the limited availability of examples of links by relying on a novel clustering-based approach, able to estimate the degree of certainty of unlabeled examples of links, so that they can be exploited during the training together with the labeled examples. Our experiments show that the proposed method can reconstruct both the human and the mouse GRNs more effectively compared to reconstructing each network separately. Moreover, it significantly outperforms three state-of-the-art transfer learning approaches that, analogously to our method, can exploit the knowledge coming from both organisms. Finally, a specific robustness analysis reveals that, even when the number of labeled examples is very low with respect to the number of unlabeled examples, the proposed method is almost always able to outperform its single-task counterpart.

Existing methods usually do not rely on a single theory, but on multiple classes of statistical/mathematical methods and information/machine learning theory. In this context, the Dialogue for Reverse Engineering Assessments and Methods (DREAM) challenges have also contributed to the development of this task. In particular, in follow-up studies, it has been shown that combining multiple approaches 14,15 or multiple sources 16,17 can be beneficial for GRN reconstruction.
Considering Gene Network Reconstruction as a machine learning task, it can be formulated as a link prediction problem via binary classification, where existing relationships among genes (i.e., gene regulation activities that have already been validated in the laboratory) can be considered as the set of positive examples. On the other hand, pairs of genes, for which there is a confirmation about the non-existence of the regulation, can be considered as negative examples. However, validation efforts and resources are usually spent to prove the existence of gene interactions, rather than their non-existence. This means that all the possible gene pairs for which there is no web-lab validation cannot be considered negative examples, but rather unlabeled examples. This context makes the adoption of classical supervised machine learning methods inappropriate or even inapplicable, and requires the design of semi-supervised learning methods, which are also explicitly able to work in the absence of negative examples, i.e., in the positive-unlabeled setting 17 . This is the most challenging setting, especially considering that the number of available positive examples is usually significantly lower than the number of unlabeled examples.
In order to face the above challenges, in this paper we propose a machine learning method for gene network reconstruction, which works in the positive-unlabeled setting and alleviates the issues arising from the limited availability of labeled data. In particular, the method proposed in this paper relies on a transfer learning approach that is able to exploit the knowledge of a source domain D s to improve the result of a task performed on the target domain D t . In our case, the data in each domain represents the expression levels measured for genes of a different organism.
Methodologically, we propose a specific kind of transfer learning approach, namely, a multi-task method 18 , whose main advantage is the ability to simultaneously solve the task on both domains, and possibly exploit dependencies between them that could lead to improved accuracy of reconstruction. In particular, we aim at simultaneously reconstructing the gene regulatory networks of two related organisms, namely, the human and the mouse regulatory networks, by considering a novel instance mapping which is guided by the notion of genetic orthology 19 .
State-of-the-art supervised machine learning methods employ a training set of examples which represents a sample of the population under analysis, described by a feature vector and associated with a known target value. These methods learn a prediction model which is able to assign a target value to unseen examples. This approach is widely proven to be effective if the set of examples is large enough and if the dataset is completely labeled, i.e., each example has a label (target value).
Therefore, classical supervised machine learning algorithms can be naturally applied to solve the task of network reconstruction, where: (i) each example corresponds to a (possible) relationship between two genes; (ii) features correspond to expression data regarding the two genes; (iii) labels can be {Yes, No}, depending on whether the interaction exists or not, or a numerical value representing the degree of certainty of the interaction. In the literature, we can find different approaches to face this challenge, that usually work in the positiveunlabeled learning setting. They can be classified according to three categories 21 : (a) two-step methods, that identify a set of negative examples from the set of unlabeled examples and then, in the second step, exploit off-the-shelf supervised learning methods to build the final predictive model [22][23][24] ; (b) instance-weighting methods, that estimate the reliability of each unlabeled example and exploit it as a weight or a cost while learning the prediction model 25 ; (c) noisy negative methods, that consider the unlabeled set of examples as highly noisy negative examples 22,26 .
The method proposed in this paper partially falls in category (b), that, according to previous studies 15,27 , allows us to avoid the imposition of strong assumptions on the negative examples, made by the methods in categories (a) and (c). However, as we explain in detail in "Methods" section, the estimated weight is used as a target value, rather than as a weight. Moreover, as introduced in "Introduction" section, the proposed method is based on a multi-task approach which simultaneously solves the network reconstruction task for two organisms, namely, human and mouse, possibly exploiting dependencies and similarities among them.
Since the method proposed in this paper solves the gene network reconstruction task by exploiting a transfer learning approach, specifically based on multi-task learning, in the following subsections, we provide some background notions and briefly review existing methods in these fields.
Transfer learning. One possible solution to overcome the scarcity of labeled examples is the adoption of transfer learning approaches 18,28 , that aim at exploiting the knowledge about another related domain D s , called source domain, to improve the quality of the results on the main domain D t , called target domain.
Formally, in a classical supervised learning setting, given (i) the feature space X of training examples, (ii) the output space Y, and (iii) n training examples {(x 1 , y 1 ), (x 2 , y 2 ), . . . , (x n , y n )} , such that x i ∈ X and y i ∈ Y , the goal is to learn a function f: X → Y , that predicts the label/value of unseen, unlabeled examples.
Transfer learning differs from this formulation since it works on two different domains. Formally, given: the goal is to learn a function f t : X t → Y t on the target domain, also exploiting the knowledge acquired by learning a function f s : X s → Y s on the source domain.
In the literature, we can find several transfer learning approaches, which were designed either as generalpurpose frameworks or as specific methods, tailored for solving specific tasks of an application domain. Such approaches can be classified according to two categories (see Fig. 1 (left)): • homogeneous, when the source and the target domains are described according to the same feature space (i.e., X s = X t ); • heterogeneous, where there are no restrictions on the feature spaces (i.e., X s = X t ).
The heterogeneous setting is clearly more difficult to handle, since it is necessary to design a strategy to transform both the feature spaces into a common feature space, or to make them comparable. For example, some heterogeneous transfer learning approaches [29][30][31] assume that the source and the target domains are described with the same number of features and identify a shared feature subspace, where the difference between data distributions is minimized.
Another categorization of transfer learning methods 18 distinguishes among: • instance-based methods, that usually perform a reweighing of the source domain instances, that are then directly used during the training for the target domain (see 16,32,33 ); • parameter-based methods, that aim to transfer the knowledge through some parameters shared by the models learned for the source and the target domains (see 17,34,35 ); • feature-based methods, that perform knowledge transfer by identifying a shared feature space (see [29][30][31]36,37 ).
Focusing on transfer learning approaches proposed in the field of bioinformatics, in the literature we can find a recent work that aims to classify breast cancer tumors, as Estrogen-Receptor-positive (ER-positive) or Estrogen-Receptor-negative (ER-negative), by exploiting two different data sources 38 . A different approach, based on deep learning, has been used for molecular cancer classification 39 , where the feature representation learned while classifying two tumor types also exploits information conveyed, in an unsupervised manner, by other tumor types. Breckels et al. 40 propose to extend a state-of-the-art transfer learning framework to solve the predictive task of mouse protein sub-cellular localization.
To the best of our knowledge, all the cited methods require a fully labeled training set, or assume the presence of some negative examples, following the strategies (a) or (c) described in "Introduction" section. In our previous work 16,17 , we overcame this limitation by designing methods based on strategy (b), i.e., based on instanceweighting. In particular, these methods exploit the knowledge coming from the reconstruction of the mouse GRN for the reconstruction of the human GRN, in a homogeneous transfer learning setting. However, the main limitations of these methods are: (i) their inability to solve the gene network reconstruction task for both organisms simultaneously, and (ii) the homogeneous setting in which they work, that makes them hardly applicable if the gene relationships of the considered organisms are represented in different feature spaces.
The approach we propose in this paper exhibits the advantages of our previous work 16,17 , without their limitations. In particular, the proposed method works in a multi-task learning setting, which aims at solving both gene network reconstruction tasks simultaneously, and which can analyze the considered organisms in either homogeneous or heterogeneous feature spaces. Moreover, according to the second categorization 18 , our method falls in the category of feature-based transfer learning methods, since, as we will describe in "Methods" section, we identify a common feature space by exploiting the concept of genetic orthology.
Multi-task learning. A specific sub-category of transfer learning methods is represented by multi-task learning methods, which aim at simultaneously solving the task for both the source domain D s and the target www.nature.com/scientificreports/ domain D t . Such an advantage is not commonly present in standard transfer learning methods, which usually aim to facilitate or improve the task for the target domain only. On the contrary, multi-task learning approaches are able to optimize both tasks simultaneously, through multiple objective (or loss) functions, or their combination. The simultaneous consideration of the two tasks allows us to take into account possible bidirectional dependencies, which cannot be considered in single-task scenarios, even if a unidirectional transfer learning approach is applied multiple times. Several complex machine learning applications have taken advantage of multi-task approaches, ranging from natural language processing 41 and speech recognition 42 to computer vision 43 and GRN reconstruction 44 . To the best of our knowledge, there is only one multi-task learning method in the literature that is able to work in a positive-unlabeled setting 45 . However, it requires that some of the solved tasks are classical supervised tasks, where the training set also includes negative examples. This makes its application inappropriate in our case, since, in principle, both the gene network reconstruction tasks are posed in the positive-unlabelled setting. This is an important aspect, as well as a strong contribution provided by our method. Indeed, although we can find several multi-task methods that are able to work in the semi-supervised setting (e.g., 46,47 ), they cannot be easily adapted to work in the positive-unlabeled setting for both the considered gene network reconstruction tasks, due to the inherent additional challenges introduced by the absence of negative examples. Therefore, our method simultaneously exhibits all the following characteristics: • it can work with no negative examples, using positive and unlabelled examples of both domains (positiveunlabelled); • it is able to transfer the knowledge acquired in the reconstruction of a GRN of an organism for the reconstruction of the GRN of another organism (transfer learning); • it can simultaneously reconstruct (see Fig. 1 (right)) the GRN of two organisms, i.e., the knowledge is transferred bidirectionally (multi-task learning).
It is noteworthy that multi-task approaches are closely related to multi-target prediction methods. Indeed, multitarget prediction refers to the (possibly simultaneous) prediction of multiple variables of interest for the same units of observation 48 . In our case, we are interested in predicting the existence of relationships between genes of two different organisms. Therefore, considering an output variable for each organism leads the considered task to be conceptually close to a multi-target prediction task (since it is in fact a multi-target link prediction task). This aspect will be clearer in the next section, where we describe how we employ our multi-target prediction approach to solve the network reconstruction task for two organisms simultaneously.

Methods
In this section, we describe our method for simultaneous reconstruction of two GRNs in a multi-task learning setting. In particular, we will focus on the reconstruction of the human and mouse GRNs. To this end, we exploit the Predictive Clustering framework, that has proved its effectiveness in the presence of different forms of autocorrelation 49 , i.e., when objects that are close to each other appear more related than distant ones. This is useful in the context of network data, like the GRNs under study, where genes close in the network are expected to show a similar behavior or to participate in the same biological processes.
In particular, we exploit the Predictive Clustering Tree (PCT) method implemented in the system CLUS 50 . CLUS is a decision tree and rule induction system that unifies unsupervised clustering and predictive modeling. It has been employed in several recent works to solve multi-target prediction tasks for a single domain. The approach proposed in this paper can be considered the first attempt to employ the PCT multi-target prediction method implemented in CLUS to work in a multi-task learning setting, where the variables to be predicted are associated with two different tasks.
Multi-target prediction methods are generally categorized according to whether they build multiple local models, i.e., one model for each target variable, separately, or one global model, i.e., a single predictive model that is able to predict the whole set of target variables simultaneously. Global models are generally more effective than their local counterparts 51 , due to their ability to capture dependencies in both the input and the output spaces. In this paper we exploit the capability of CLUS to learn a global model, and we consider the degrees of existence of each gene interaction in the two organisms as two target variables of a multi-target regression task. This is achieved by representing the same examples of gene interactions, in the two organisms, in a common feature space. In particular, in order to find a match between the genes in the two organisms, we exploit the concept of orthologous genes, that are genes in different species that originated by vertical descent from a single gene of the last common ancestor.
Before describing in detail the proposed multi-task approach, we introduce some useful notions and formally define the problem that we solve. Let: • G h (resp., G m ) be the set of considered genes for the human (resp., mouse) organism; www.nature.com/scientificreports/ • orth hm : G h → G m (resp., orth mh : G m → G h ) be a function that takes a human gene g h ∈ G h (resp., a mouse gene g m ∈ G m ) and returns the corresponding orthologous mouse gene (resp., human gene); • vec h : G h → R r (resp., vec m : G m → R q ) be a function that returns the r-dimensional (resp. q-dimensional) vector of expression levels of a human (resp., mouse) gene; • ⌢ : R n 1 × R n 2 → R n 1 +n 2 be a function that takes as input two vectors in R n 1 and R n 2 , respectively, and returns their concatenation in R n 1 +n 2 .
, be a function that takes as input two human (resp. mouse) genes and returns the concatenation of their vectors of expression levels, representing the features of their interaction. Formally, The task solved by our multi-task learning approach is to find two regression functions, namely: , returns a score representing the degree of certainty of the existence of the interaction between g ′ h and g ′′ h in the human GRN. • f m : R 2q → [0, 1] , that, given a pair of mouse genes g ′ m ∈ G m and g ′′ m ∈ G m represented through the feature vector of their interaction ṽec m (g ′ m , g ′′ m ) , returns a score representing the degree of certainty of the existence of the interaction between g ′ m and g ′′ m in the mouse GRN.
Our goal is to learn both predictive functions simultaneously, by considering the degree of certainty of a given gene pair for the human and the mouse organisms as two different target variables of the same training example. It is noteworthy that this choice allows our method to capture possible dependencies that may exist in the output space (i.e., between the target variables). Specifically, we learn a single regression function f hm that takes as input a pair of genes represented according to the features related to both organisms, and returns the degree of certainty for both organisms. Formally: Note that the construction of all-in-one training examples that can be used to learn a multi-target prediction model needs an additional step, i.e., the identification of a match between human genes and mouse genes. In the following subsections, we describe (i) the details of such a step, (ii) the strategy we adopt to solve the issues of the positive-unlabeled setting, (iii) the construction of the dataset used for learning the multi-target regression function f hm , and (iv) the proposed predictive approach.
Orthologous matching and construction of positive training examples. The first step of our method consists of the identification of possible matches between human and mouse genes. This step is necessary in order to represent each gene pair as a single training example, according to the features (i.e., expression levels) measured for both organisms.
To this aim, we exploit the concept of gene orthology. Ortholog genes are genes of different species that are the result of the speciation of the same originating gene (see Fig. 2). Methodologically, we iterate over the human genes g m ∈ G m and identify the corresponding orthologous gene in the mouse organism (Algorithm 1, Lines 2-27).
At the end of this step, we obtain two new sets of genes: G ho ⊆ G h consisting of the human genes that have orthologs in the set of mouse genes G m , and G mo ⊆ G m consisting of mouse genes that have orthologs in the set of human genes G h .
The subsequent steps of the method work on the orthologous subsets of genes G ho and G mo . From a machine learning viewpoint, the set of genes corresponds to the set of nodes of the GRNs. However, our unit of analysis is a pair of genes, for which we want to estimate/predict the degree of existence of the interaction. Given that a human gene g h is described as a vector of expression levels vec h (g h ) , we represent a pair of genes g ′ h and g ′′ h as the concatenation of their feature vectors, i.e., In particular, given two feature vectors u h (for the human organism) and u m (for the mouse organism) for the same pair of genes, we compute the value of the target variables t h and t m as follows: where cent(c) is the feature vector of the centroid of the cluster c, and sim: R n × R n → [0, 1] is a vector similarity function. In this paper we use sim(a, b) = 1 − √ n i=1 (ai−bi) 2 n , based on the Euclidean distance, after applying a min-max normalization (in the range [0, 1]) to all the features. The identification of the clusters can actually be performed through any centroid-based clustering approach. In this paper we rely on the well-known k-means algorithm. Moreover, in order to optimally identify the number of clusters k h for the human organism and k m for the mouse organism, we use the silhouette cluster analysis 52 . Formally, we define a function sil: P → [1, 2, . . . |P|] , that, given a set of positive examples P ∈ {P h , P m } and the clustering algorithm (in our case, k-means) returns the optimal number of clusters, according to the silhouette analysis. In Algorithm 1, this step is performed at Lines 13-14, whereas the exploitation of the identified clusters for computing t h and t m is performed at Lines 17-23.
After this step, the main issues of the positive-unlabeled setting are solved, since all the examples are associated to a (known or estimated) value for the target variables t h and t m .
Learning the predictive model. The final stage consists of learning the predictive model, in the form of a multi-target regression function, where the two target variables t h and t m represent the degrees of certainty of the existence of the interaction in the human and in the mouse organisms, respectively. With this aim, we build the final training set by concatenating, for each pair of genes for which we identified an ortholog match, (i) the 2r-dimensional feature vector associated to the human organism, (ii) the 2q-dimensional feature vector associated to the mouse organism, (iii) the target variable t h , and (iv) the target variable t m , leading to training examples represented in R 2r+2q , associated to two target variables (see Algorithm 1, Line 24).  www.nature.com/scientificreports/ We learn the predictive model with CLUS 50 , that is based on Predictive Clustering Trees (PCTs). We induce PCTs through a standard approach for the top-down induction of decision/regression trees, that takes as input a set of training examples and returns the induced tree. The heuristics adopted to select the best tests of the internal nodes of the tree is the reduction of variance achieved by partitioning the examples according to such a test. The maximization of the variance reduction leads to maximizing the cluster homogeneity and, therefore, to improving the predictive performance. Therefore, the considered heuristic is formally defined as where Var E (t h ) (resp., Var E (t m ) ) is the variance observed on the target attribute t h (resp., t m ) on the set of examples E falling in a given node of the tree.
This means that the variance reduction, used to identify the best candidate split in the tree construction, is computed as: www.nature.com/scientificreports/

Experiments
In this section, we present the results of our experimental evaluation. All the experiments were performed on a server equipped with a 6-cores CPU @ 3.50Ghz and 128GB RAM. In the following subsections, we first describe the considered competitor systems, the datasets and the experimental setting. Finally, we present and discuss the obtained results.
Competitor approaches. We compared our method with the following competitor approaches: • TJM 29 , that reduces the difference between the two domains by identifying a match between their features and by reweighting the instances to construct a new reduced/shared feature representation; • BDA 31 , that adaptively leverages the importance of the marginal and conditional distribution discrepancies between the two domains; • JGSA 30 , that learns two coupled projections, that are exploited to project the source and the target domain data into low-dimensional subspaces, where the geometrical and the distribution discrepancies are minimized;  www.nature.com/scientificreports/ • no_transfer, that is the single-domain variant of our approach, which reconstructs each single GRN independently.
TJM, BDA, and JGSA are feature-based transfer learning methods that are able to share the knowledge between (also) heterogeneous source and target domains, as long as they are described with the same number of features.
On the other hand, the no_transfer approach can be considered a baseline, that allows us to evaluate the positive contribution of the multi-target approach proposed in this paper or, conversely, to evaluate the possible presence of negative transfer phenomena 53,54 , where the use of knowledge coming from other domains actually compromises the quality of the reconstruction.
Datasets. We built the dataset by downloading a compendium of microarray data of both human (Platform ID: GPL570) and mouse (Platform ID: GPL1261) organisms from Gene Expression Omnibus-GEO (www. ncbi.nlm.nih.gov/geo/), a publicly available web repository hosted by the National Center for Biotechnology Information (NCBI). In total, 174 and 161 raw CEL files related to 54,675 and 45,101 control probesets of 6 different organs were downloaded for human and mouse organisms, respectively (see Supplementary Table S1 for a complete list of accession numbers). More specifically, the 174 CEL files for the human organism are distributed as follows: 17 for bone marrow, 37 for brain, 4 for heart, 7 for liver, 45 for lung, and 64 for skin. On the other hand, the 161 CEL files of the mouse organism are distributed as follows: 14 for bone marrow, 8 for brain, 8 for heart, 124 for liver, 4 for lung, and 3 for skin. We processed the data following the workflow proposed in the DREAM5 challenge 14 (see Fig. 4 for a graphical overview of the followed pipeline). In particular, we performed the Robust Multi-array Average (RMA) 55 normalization through Affymetrix Expression Console Software as one batch per organ. Data were background adjusted, quantile normalized, median polished and log-transformed. The mapping from Affymetrix probeset IDs to gene IDs led to a total of 10,886 human genes and to 11,655 mouse genes.
After the ortholog matching (see "Orthologous Matching and construction of positive training examples" section), we obtained a reduced set of 3, 196 genes for both organisms. The strategy adopted to perform the ortholog matching was based on the gene symbol, that corresponds between homo sapiens and mus musculus organisms, except for differences in the capitalization 56 . Alternative solutions may have been adopted, mostly based on explicit lists of ortholog genes (see, for example, the OMA orthology database 57 ), but the strategy based on the gene symbol provided us 153 additional matches. Note that the set of 3,043 genes was included in the set of 3,196 genes we considered.
The dataset of possible relationships between genes was built by considering all the possible pairs of genes (more than 10 million, excluding the self-links), each associated with the concatenation of the feature vectors of the involved genes (following the strategy described in "Methods" section). This step led to 348-dimensional vectors for human gene pairs and to 322-dimensional vectors for mouse gene pairs. We exploited the database BioGRID 20 as the source of known validated interactions (i.e., to define the sets B h and B m ), while the remaining possible relationships were considered unlabeled. In Table 1 we report a summary of the quantitative characteristics of the considered dataset.
Since some competitor systems, even if they are able to work in the heterogeneous transfer learning setting, require the number of features of all the domains to correspond, we also built a homogeneous version of the dataset. In particular, we aggregated the features associated to each organ (by averaging their value), leading to a homogeneous dataset consisting of 6 features per gene, for both the human and the mouse organisms. Table 3. Summary of settings for which the multi-task approach provided an improvement over the baseline.

Measure
50% 40% 30% 20% 10% 5%    Table 2). It is worth mentioning that the comparison with competitors has been performed on the smallest version of the dataset (i.e., with the ratio 50%), because they were not able to complete the experiments with larger datasets without incurring in RAM exhausting errors on our servers. Indeed, while our approach is based on a top-down induction of regression trees, that is generally efficient, competitor methods are mainly based on matrix computations and easily exhausted the RAM on our servers even with the considered reduced dataset. Results. Figure 6 depicts the improvement obtained by our approach with respect to the baseline no_transfer, in terms of the AUR@K, AUROC and AUPR measures, for both the homogeneous and heterogeneous datasets, and with respect to both organisms.

Heterogeneous human GNR
The charts show that the proposed approach provides a marked improvement over the single-domain counterpart in the reconstruction of the human GRN. Such an advantage is evident for all the variants of the dataset, i.e., for all the considered labeled/unlabeled ratios. On the other hand, the reconstruction of the mouse GRN appears to exploit the knowledge about the human GRN only with higher labeled/unlabeled ratios. This may suggest that the mouse organism can be considered an appropriate model organism for the study of the human GRN, but the contrary may hold to a lesser degree, i.e., only when there is a sufficient amount of biologically validated information.
In general, for both organisms, the higher the labeled/unlabeled ratio, the better the quality of the reconstruction. This is an expected result, since the unlabeled examples could belong to clusters of positive examples that were not properly represented/observed in the set of positive examples, due to their limited availability. Despite such an expected result, the reconstruction performed with our multi-task approach in most of the cases provides an advantage, even with very low labeled/unlabeled ratios. The few cases in which there is not an improvement occur in the reconstruction of the mouse regulatory network. This may indicate that it is easier and more natural to exploit the knowledge coming from a general/simple organism to describe a complex organism, rather than vice versa. In Table 3 we show a summary of the settings where our multi-task approach provided an improvement over the baseline.
Moreover, as explained in "Competitor approaches" section, we compared our results with those obtained by three state-of-the-art approaches. In Fig. 8 we present the results in terms of recall@k, AUR@K, AUROC and AUPR. Observing the figure, it is clear that our approach outperforms all the competitors by a large margin for both the considered organisms.
We also performed a qualitative evaluation of the networks reconstructed by our system. For this specific analysis, we considered the largest version of the dataset, containing 3, 970 positive examples and 75, 430 unlabelled examples (see the variant 5% in Table 2). Since the experiments were performed using 10-fold cross validation, resulting in 10 different rankings (one for each fold), we averaged the scores and analyzed the resulting ranking. We then selected the top 10, 000 ranked interactions for both organisms, we computed some topological measures (see Supplementary Tables S2 and S3 for a detailed overview) and we identified the hub genes, by selecting the top-10% of genes with the highest numbers of regulated genes 59 . Among them, we selected the 352 genes appearing as hubs for both the considered organisms, and we plotted the subnetworks involving each of them, emphasizing the interactions that were present in BioGRID (in black) and those that were predicted by our system (in green). It is noteworthy that the subnetworks of the first 60 hub genes are identical, in both known and predicted interactions, between the two organisms. This confirms the ability of our system in catching crossorganism similarities and in predicting the existence of interactions that appear coherent among the organisms.
In Fig. 9 we depict the first subnetwork that shows some differences between the organisms (i.e., the 61th in the ranking). This is the case of the human gene PIK3R1 (resp., Pik3r1 for the mouse organism). In this case, we can observe 5 (resp., 6) predicted regulated genes for the human (resp., for the mouse) organism and 5 (resp., 5) predicted genes regulating PIK3R1 (resp., Pik3r1). Specifically, it is noteworthy that the interaction Pik3r1 → Alk has been inferred by our method, but is not covered in BioGRID. On the other hand, the interaction Csf1r → Pik3r1 is present in BioGRID for the mouse organism, but our method did not suggest the corresponding Scientific Reports | (2020) 10:22295 | https://doi.org/10.1038/s41598-020-78033-7 www.nature.com/scientificreports/ interaction for the human organism, that is actually absent in BioGRID (preventing, therefore, a possible false positive). This confirms that, although our method exploits the knowledge coming from the simultaneous reconstruction of the regulatory networks of both the organisms, it does not merely mimic the behavior observed on an organism on the other one. On the contrary, it is able to catch possible differences and asymmetries. Finally, we performed an additional analysis regarding the importance of the considered features. In particular, focusing on the homogeneous dataset, we aimed to identify the most relevant organs (i.e., those most relevant to our method during the learning of the multi-target regression model), following the approach of Petković et al. 60 . In Fig. 7 we show the obtained ranking, where we can observe that the features associated to the human skin and heart, together with those associated to the mouse heart and lung, have been detected as the most relevant ones for the gene network reconstruction task. In contrast, it seems that features related to bone marrow (for both organisms) did not provide any relevant contribution. In the middle, we find the features closely related to the brain (for both organisms), the human liver and the human lungs. This behaviour is probably motivated by the fact that some organs show more similar properties between the two organisms, or are better connected through orthologous genes, than others. It is noteworthy that these findings can be profitably exploited to focus future work, where larger sets of samples, related to the organs that provide a higher contribution, can be adopted.

Conclusion
Several computational approaches, mainly based on machine learning methods, can be employed for the reconstruction of GRNs. However, existing gene network reconstruction methods suffer when the number of labeled examples is low, especially when no negative examples are available. In this paper we have proposed a method that overcomes these limitations. Our approach is able to simultaneously reconstruct the GRN of two organisms, by exploiting a multi-target regression approach that, in conjunction with the concept of gene orthology, is able to natively work in a positive-unlabeled learning setting.
The experiments show that our approach is able to really "transfer" knowledge extracted from an organism and profitably use it in another organism. Moreover, the proposed multi-target positive-unlabeled learning algorithm outperforms both its single-task counterpart and three state-of-the-art transfer learning approaches in the reconstruction of both GRNs.
As future work we plan to define a more general approach to map the examples in the considered domains, so that it may be adopted in multiple, even non-biological, applications. Moreover, while in the present paper we presented our novel approach and evaluated its effectiveness, compared with state-of-the-art methods, we plan to extend the experiments to larger datasets, also considering different pipelines.