in a recent breast 
 cancer study conducted by the clinical 
 proteomic tumor analysis consortium ( cptac ) , extensive global- and phospho- protein profiling were obtained for 
 a subset of breast cancer samples that have been extensively characterized 
 in the the cancer genome atlas ( tcga ) .
this is so far the first attempt to study protein 
 activities in breast cancer samples using sophisticated protein experiments 
 on a large scale .
then , by leveraging genomic analytical outputs from 
 tcga on these samples , we have the unique 
 opportunity to characterize and compare gene coexpression networks 
 and protein coexpression networks based on the same set of samples . 

this is of great interest because knowledge about the common and different 
 coexpression patterns among rnas and proteins could improve our understanding 
 of complicated gene regulatory mechanisms in tumor samples and could 
 also facilitate the detection of important disease genes and therapeutical 
 targets . in the past decade , numerous statistical and machine 
 learning methods 
 have been proposed to construct gene  gene regulatory networks 
 and protein  protein regulatory networks , including coexpression 
 network methods , bayesian network methods , and gaussian graphical models .
review of these and other methods are available in 
 hecker et al . and lee et al .
; however , all of these proposed models are designed 
 to construct one network at a time .
presumably , we can use rna expression 
 data and protein expression data to construct two networks separately 
 and then compare , but such an approach is certainly less optimal , 
 as gene expressions and protein expressions in one tumor sample are 
 closely related . on one hand ,
protein levels are regulated by their 
 rna levels , so we expect the two coexpression networks share common 
 structures . on the other hand ,
protein activities are subject to a 
 large amount of post - transcriptional modifications , and thus we also 
 expect to observe unique interaction patterns in each network .
this 
 motivates us to perform joint learning of both rna - seq and proteomic 
 networks , so we are able to borrow information across different data 
 sets and better capture the common correlation structure , which shall 
 then improve the overall accuracy of network estimation .
some of them have been specifically designed to estimate time - varying graphical 
 models in the context of time - series data .
proposed 
 likelihood - based methods for the joint estimation of multiple related 
 gaussian graphical models .
their approaches rely on two different 
 regularization schemes penalizing differences across partial correlation 
 structures of different classes . while both papers successfully demonstrate 
 the advantage of jointly modeling multiple networks over estimating 
 individual networks separately , the performance of the methods heavily 
 depends on the multivariate gaussian assumptions of gene expression 
 distributions , which may not hold in real biological systems .
recently , random - forest - based methods have been utilized for the 
 construction of grn .
random forest is a decision - tree - based 
 nonparametric method for building powerful prediction models .
its ensemble structure combined with its nonparametric 
 nature delivers excellent performance even when the sample size is 
 moderate .
the 
 key idea is to model the expression of each target gene as a function 
 of the expression of all other genes via random forest .
the superior 
 performance of genie3 was demonstrated in both dream 4 in - silico multifactorial 
 challenge and dream 5 network inference 
 challenge . in both challenges ,
, genie3 was further modified to estimate grn 
 based on time - series data .
proposed irafnet , another random - forest - based method that 
 combines different types of heterogeneous data , such as protein  protein 
 interactions , transcription factor - dna binding , and gene knock - down 
 to estimate grn .
while these random - forest - based methods greatly advance 
 the field of grn construction , none of them is designed to estimate 
 simultaneously multiple related networks .
we propose jrf ( joint 
 random forest ) , a new algorithm for the simultaneous 
 estimation of protein coexpression and gene coexpression networks . 

the key idea is to borrow information across different data by forcing 
 the class - specific tree ensembles to use the same genes for the splitting 
 rules . in this way
, regulatory relationships playing important roles 
 in multiple classes will be detected with better power .
after testing 
 the performance of our algorithm on several in silico experiments 
 and comparing jrf with different methods including genie3 the joint graphical lasso ( jgl ) , we applied jrf for the simultaneous construction 
 of gene coexpression networks and protein coexpression networks . for 
 our analysis , we considered gene expression from tcga and protein 
 expression from cptac data for 62 breast cancer patients .
we derived 
 the two networks and identified protein - specific hub genes and gene 
 modules , showing the unique contribution of proteomic data to breast 
 cancer research .
random forest 
 is a nonparametric algorithm that models a response variable via a 
 collection of decision trees .
specifically , each decision tree is 
 constructed based on a random subset of training samples , and the 
 splitting variable at each node of a decision tree is chosen from 
 a randomly sampled subset of predictors upon maximizing a certain 
 utility function ( e.g. , decrease in node impurity ) .
recently , 
 huynh - thu et al . introduced genie3 , a 
 random - forest - based model for inferring gene regulatory networks ( grns ) . 
 in genie3 , first
, for each target gene k , its expression 
 is modeled as a function of the expression of all other genes via 
 random forest .
then , regulatory event ( j  k ) is measured by ijk , the importance score of gene j in the random forest model of gene k. 
 specifically , the importance score ijk is defined as the summation 
 of node impurities across all nodes that utilize predictor j for the splitting rule divided by the total number of 
 trees in the random forest model of gene k. we start with some notations .
for each class g  { 1 , ... , g } , denote the gene 
 expression data for p genes and ng individuals as xngpg , and 
 the ith observation of gene j under 
 class g as xijg .
an overview of jrf is shown in figure 1 , where , for simplicity , we consider only 
 two classes .
the key step of jrf is to construct g random forest models simultaneously using data from g classes when predicting the expression of a target gene k based on the expression of all other genes .
we propose 
 to use the same predictor variables for splitting rules in different 
 trees corresponding to different classes .
the goal is to borrow information 
 across different classes , so that regulatory relationships can be 
 better detected if there are coherent signals across different classes . 

specifically , when we grow g decision trees in parallel 
 for g classes , at one node  , we decide the 
 splitting variable based on the following procedure : ( 1 ) randomly 
 sample a set   containing n genes from 
 the entire set of genes except gene k. ( 2 ) 
 for each class , the best splitting rule based on each candidate 
 predictor in set   is derived .
let   be the subset of samples allocated to node 
  under class g. the splitting rule based on 
 the jth predictor is of the formwhere   and   are the subsets of samples in class g allocated to the left and right children of node  
 with a splitting rule based on the jth predictor . 
 for each predictor in set , the best threshold agj*( )
is derived as followswhereand   with x representing the mean of set . in other 
 words , cgj is the decrease in node impurity 
 after splitting node  in 
 the gth tree according to gene j. ( 3 ) finally , gene g * is chosen as the splitting 
 variable 
 of node  for all classes based on the following criteria : because 
 jrf utilizes the same splitting variable 
 for corresponding nodes across all classes , predictors associated 
 with gene k in multiple classes are more likely to 
 be chosen for the splitting rules . for each step in the tree construction , 
 ( 13 )
only apply to classes for which  is not a final 
 node . as in the original random forest model
, each tree grows until either the total number of observations allocated 
 to the final leaves falls below a certain prespecified threshold or 
 the maximum number of possible nodes is reached .
it is worth mentioning 
 that when the number of classes is one , jrf reduces to genie3 , the 
 original random forest model for network inference .
it is worth noting 
 that to implement jrf , data sets need to be standardized to mean zero 
 and unit variance .
importance scores depend on the scale of the data 
 and therefore variables ( genes ) need to be standardized before the 
 random forest models are fitted.figure 1jrf schematic . for simplicity ,
let us assume 
 that there are only 
 two classes and that each data contains the same number of samples . 
 for each target gene , gk
, we run random forest for each class and the two tree ensembles 
 are forced to share some structure as follows .
first , for each node 
  , the same set of genes s is randomly sampled from the entire set of genes ( a ) .
then , the 
 same variable in set s is chosen 
 for the splitting rule .
this is achieved by deriving the splitting 
 rule for any gene contained in set s ( b ) and choosing the best variable to maximize the sum of decrease 
 in node impurity over different classes ( c ) .
let us assume 
 that there are only 
 two classes and that each data contains the same number of samples . 
 for each target gene , gk , we run random forest for each class and the two tree ensembles 
 are forced to share some structure as follows .
first , for each node 
  , the same set of genes s is randomly sampled from the entire set of genes ( a ) .
then , the 
 same variable in set s is chosen 
 for the splitting rule .
this is achieved by deriving the splitting 
 rule for any gene contained in set s ( b ) and choosing the best variable to maximize the sum of decrease 
 in node impurity over different classes ( c ) . on the basis of the previously 
 described procedure , jrf constructs g random forest 
 models simultaneously for each target gene k and 
 returns a ranking of gene  gene interactions based on importance 
 scores . to assess the performance of jrf in predicting the true interactions , 
 receiver operating characteristic curves ( roc ) and precision - recall 
 curves can be computed by setting different thresholds on importance 
 scores . in this paper ,
jrf is compared with jgl , genie3-sep , and genie3-comb on several in silico experiments . 

genie3-sep is genie3 used to estimate networks based on data from g classes separately , while genie3-comb is genie3 used to 
 estimate a unique network based on the union of data from all classes . 

all random - forest - based algorithms ( jrf , genie3-comb , and genie3-sep ) 
 provide a ranking of regulatory relationships based on importance 
 scores , and roc curves were computed by setting different thresholds 
 on these scores . instead , for jgl
, roc curves were constructed by 
 considering different values for the two parameters controlling the 
 level of sparsity ( as shown by figure s1 in the supporting information ) .
another approach to evaluate 
 the performance of jrf is choosing a proper cutoff value for importance 
 scores using permutation techniques .
let ijkg be the importance score 
 associated with the regulatory event ( j  k ) in the gth class - specific tree ensemble . 

specifically , this importance score is defined as   where t is the number 
 of trees and   is the set of nodes which utilize gene j for the splitting rule in the tree ensemble used to predict 
 gene k based on the gth class - specific 
 data . because in this paper we are interested in estimating undirected 
 networks , we derive importance score ij 
kg for every edge ( j  k ) as the average between importance 
 scores ikjg and ijkg , and 
 the final undirected networks can therefore be derived by applying 
 a cutoff on the edge importance scores . to derive a proper cutoff 
 value for importance scores , we utilize the following permutation - based 
 procedure : ( a ) for b  { 1 ,  , b } , with b being the number of permutations : ( a.1 ) for any target gene k , we first permute 
 its sample order within each class and fit g random 
 forest models via jrf to predict the expression of gene k based on the expression of all other genes in g classes , respectively .
( a.2 ) repeat ( a.1 ) for all genes and 
 compute the final importance 
 scores for each edge in each class , which are denoted as { ippg(b)}g=1 g .
( b ) for each threshold  , we 
 computewhere 1( ) is 
 the indicator 
 function , equal to one if event  occurs and zero otherwise .
f( ) can serve as an approximation of the 
 false discovery rate ( fdr ) .
in the following 
 numerical studies , we use 0 = min{ : f( )  0.001 } and declare an edge between j and k in class g if ij 
the computational 
 complexity 
 of jrf is o(ptng=1glog(ng)ng ) , where p is the number of genes , t is the number 
 of trees in each forest , n is the number of variables 
 sampled at each node , g is the number of classes , 
 and ng is the number 
 of samples in each class .
this complexity is in the same order as 
 the complexity of applying genie3 to estimate g networks 
 for g classes separately . in the contrast , the time 
 complexity of ggm - based approaches , such as jgl , often has the order 
 of p. so 
 when p is large , jrf and genie3 enjoy better computational 
 efficiency , which is supported by our numerical investigations ( see supporting information ) . in practice ,
the 
 number of trees ( t ) and the number of potential regulators 
 to be sampled at each node ( n ) are user - specified 
 parameters .
t is usually chosen sufficiently large 
 because the tree ensemble provides more accurate results as t increases . in our numerical studies , we reported results 
 for t = 1000 .
we also evaluated larger t values , such as 2000 and 5000 , but observed similar results ( data 
 not shown ) . for n ,
a common choice is   where ( p  1 ) is 
 the number of predictors for each target gene . in general , large value 
 of n results in predictions with high bias ; while 
 low value results in predictions with high variance .
we implemented 
 an r package jrf using c routines wrapped with an 
 r interface .
similarly 
 to other tree - based algorithms for network construction , jrf can be 
 easily parallelized because the estimation process consists of p - independent subproblems .
a comparison between jrf , jgl , 
 genie3-sep , and genie3-comb in terms of computational time is presented 
 in table s1 . for data generation
, we first considered 
 two networks containing 500 nodes connected by 498 and 249 edges , 
 respectively .
in particular , network 1 is a network with two disjoint 
 components and network 2 consists of only the first component of network 
 1 . in this example
, we focused on two network topologies : power law and star topology ( figure s2 in supporting information ) . for each network topology , we simulated 20 replicates involving n = 50 samples from a gaussian graphical model .
network 
 structures and covariance matrices were simulated in the same way 
 as danaher et al .
( 2014 ) ( further details 
 can be found in the supporting information ) . besides this example
, we also consider a series of other different 
 simulation settings , including ( 1 ) two non - nested networks ; ( 2 ) five 
 non - nested networks ; ( 3 ) two nonrelated networks ; and ( 4 ) data from 
 non - gaussian distributions .
the expressions of 5864 proteins 
 for 62 breast cancer tumors ( proteome-ratio-norm-nona-normal.gct(04-jun-2014 
 version ) ) were downloaded from the cptac data coordinating center 
 ( http://proteomics.cancer.gov/programs/cptacnetwork ) sponsored 
 by the national cancer institute . for three samples having technical 
 duplicates
we further standardized protein expression so that each 
 sample had median 0 and median absolute deviation ( mad ) 1 . finally , 
 we excluded proteins with interquartile range < 70% quantile .
level - three rnaseq 
 data of breast tumor samples were obtained from the tcga web site 
 ( http://tcga-data.nci.nih.gov/tcga/ ) .
first , we log - transformed 
 the data and then replaced all missing values with 0 . then , we standardized 
 each sample to have median 0 and mad 1 .
we focused on the 62 samples 
 contained in the proteomic data and excluded genes with more than 
 10% missing values .
then , we selected the top 10% genes with largest 
 interquartile range across samples .
as final subset , we considered 
 528 genes obtained by overlapping the set of 1464 selected genes and 
 the set of 1759 selected proteins . for each data
in this section , 
 jrf was compared with 
 jgl , genie3-sep , and genie3-comb on several 
 in - silico experiments . in all implementations , we used t = 1000 and .
we evaluated 
 the performances of different 
 algorithms based on different network structures , sample sizes , and 
 number of genes .
figure 2 shows the average roc curves 
 of each method 
 over 20 replicates for network 1 and network 2 separately .
jrf outperforms 
 all other methods under all different settings , while genie3-sep is 
 the second best performer .
figure s3 shows 
 the boxplot of the area under the precision - recall curve ( aupr ) for 
 jrf , genie3-sep , and genie3-comb .
the advantage of jrf over genie3-sep is more 
 evident for network 2 , which contains only shared edges .
this is expected 
 because jrf should be able to detect common edges with better power .
mean of 
 roc curves for network 1 ( first column ) and network 2 ( second 
 column ) over 20 replicates for jrf ( red ) , genie3-comb ( green ) , genie3-sep 
 ( blue ) , and jgl ( black ) .
for each replicate , we sampled 50 samples 
 from a gaussian graphical model based on the power law topology ( first 
 row ) and star topology ( second row ) . to further assess the performance of the proposed model , 
 we included 
 more simulation scenarios in the supporting information .
in particular , in supporting information section 1.3 , we considered the case of non - nested networks sharing 
 some structure , while , in section 1.4 , we considered the case where 
 five different networks were estimated simultaneously .
as shown by figures s4 and s5 , jrf outperformed competitors 
 under all of these simulation scenarios .
in addition , it would be 
 interesting to evaluate the performance of jrf when class - specific 
 networks do not share any structure .
for this purpose , we simulated 
 two different networks with no common edge and applied different methods 
 on the corresponding data ( supporting information section 1.5 ) .
as shown in figure s6 , 
 performance of jrf is very comparable to that of genie3-sep on this 
 simulation example .
this suggests that jrf is sufficiently robust 
 and works in the absence of common structure .
moreover , to assess 
 the performance of jrf in the presence of nonlinearities in the data , 
 we further generated data examples using genenetweaver , an open - source software for the generation of 
 in - silico data from a set of ordinary and stochastic differential 
 equations ( supporting information section 
 1.6 ) .
as shown in the supplemental figure s7 and s8 , jrf again outperformed other algorithms in terms of both 
 area under the roc curve ( auc ) and area under precision - recall curve 
 ( aupr ) .
finally , figure s9 shows results 
 for different algorithms in estimating networks with 1000 nodes ( further 
 details can be found in supporting information section 1.7 ) , in which we can see that jrf outperformed competitors 
 in terms of both auc and aupr .
it is worth mentioning that the increased 
 dimensionality did not negatively affect results , suggesting that 
 jrf can handle problem with high dimensionality . in this section ,
we assess the 
 performance of jrf and other algorithms on data sets involving different 
 number of samples .
figure 3 shows boxplots of auc over 20 replicates resulting from jrf , 
 genie3-sep , and genie3-comb for different sample sizes , that is , n = 50 , 100 , 200 .
the aucs of jrf are significantly larger 
 than that of genie3-sep based on 20 replicates under all sample sizes . 

as shown in the supplemental figure s10 , jrf outperforms competitors also in terms of aupr for different 
 samples sizes .
as shown , the smaller the sample size the bigger the 
 improvement of jrf over genie3-sep is .
this result is expected because 
 as the sample size increases genie3-sep should have sufficient information 
 to estimate the networks separately .
boxplot of auc over 20 replicates for 
 jrf ( red ) , genie3-comb ( green ) , 
 and genie3-sep ( blue ) for different sample sizes , that is , n = 50 , 100 , and 200 .
data : samples simulated from a gaussian 
 graphical model on a power law topology . in this section ,
we compare jrf and 
 genie3-sep regarding their abilities to estimate class - specific edges . 

we considered the 20 replicate data sets ( n = 100 ) , 
 which was simulated based on the power law network topology in the 
 previous section . for both jrf and genie3-sep ,
edges were declared 
 using the permutation procedure illustrated above with b = 500 and an fdr threshold of 0.001 .
table 1a shows the performance of jrf and genie3-sep 
 in estimating network 1 and network 2 . for each network , we reported 
 the minimum and the maximum value of true positive rate ( tpr ) , false - positive rate ( fpr ) , and false discovery 
 rate ( fdr ) across 20 replicates . in particular , these 
 quantities are defined as tpr = tp/(tp + fn ) , fpr = fp/(fp + tn ) , and fdr = fp/(fp + tp ) , where tp is the number 
 of true positives , fp is the number of false positives , tn is the number of true negatives , and fn is the number of false negatives .
for both networks , jrf achieved 
 better power ( tpr ) than genie3-sep , while the resulting fprs were quite comparable .
table 1b shows a comparison of the two algorithms 
 regarding the detection of class - specific ( differential ) edges .
again , 
 we computed the rate of true differential edges and false differential 
 edges and showed the minimum and the maximum value of these quantities 
 over 20 replicates . while the tprs for detecting 
 differential edges were similar for the two methods , the fprs of genie3-sep were significantly larger than that of jrf .
in addition , 
 jrf leads to substantially lower fdrs for detecting 
 differential edges than genie3-sep .
for each quantity 
 ( tpr , fpr , and fdr ) we show the minimum 
 value ( min ) and the maximum value ( max ) across 20 replicates .
note 
 that the true discovery rate ( tdr ) is defined as tdr = ( 1  fdr ) .
we applied jrf to construct coexpression 
 networks based on the cptac 
 global proteomics data set and tcga rnaseq data of 62 breast cancer 
 samples . for simplicity
, we will refer to the 
 coexpression networks constructed by jrf from global proteomics data 
 and rnaseq data as protein - network and rna - network . at fdr cutoff 
 of 0.001
, we detected 687 common edges shared by both networks , 502 
 edges unique to protein - network , and 382 edges unique to rna - network . 

the topologies of the two networks are shown in figure 4 , while the full list of interactions can 
 be found in additional file 1 in the supporting information .
interestingly , a few network modules in protein - network 
 have a high percentage of edges unique to protein - network .
figure 5 highlights genes whose degrees ( i.e. , total 
 number of edges connecting to one gene ) in protein - network are much 
 higher than that in rna - network .
three out of the top five genes with 
 more protein - unique edges , c3 , cfb , and mpo , are related to immune 
 response functions .
it is well known that immune response plays a 
 major role in the patient response to chemotherapy treatments .
some of these genes have been shown to be predictive 
 of survival outcomes of breast cancer patients .
for example , c3 and 
 cp are predictive of chemoresistance in breast cancer patients , and high activity of mpo genotypes can enhance 
 efficacy of chemotherapy for early - stage breast cancer .
. 
 genes belonging to the mmp family have been linked to cancer progression 
 for their ability to degrade the extracellular matrix and may be implicated 
 in the formation of metastasis .
network based on proteomic 
 data ( a ) and rnaseq data ( b ) resulting 
 from jrf .
green edges ( 687 ) are shared between proteomic and rnaseq 
 data , red edges ( 502 ) are unique to proteomic data , while blue edges 
 are unique to rnaseq data ( 382 )
. only names of genes with at least 
 ten connecting edges are shown in the plot .
for both networks , the 
 full list of interactions can be found in additional file 1 .
number of connecting edges for top connected 
 genes in the network 
 based on proteomic data . for each gene , the green bar corresponds 
 to the number of connecting edges shared between proteomic and rnaseq 
 data ; the red bar indicates the number of protein - specific edges ; 
 while the blue bar indicates the number of rna - specific edges . for 
 each gene , we list its biological functions .
we then focused on network modules unique to protein - network . 
 we 
 derived network modules based on edge betweenness using function  edge.betweenness.community  
 available in the r package igraph . in figure 6
, we show the protein - network 
 with gene modules enriched of at least one go category .
the full list 
 of genes contained in each module can be found in additional file 
 1 .
figure 6 shows , 
 for each module , the most enriched go category and the corresponding 
 benjamini adjusted p - value .
two interesting protein - specific modules ( with more protein - specific 
 edges than shared edges ) are c1 and c2 .
fibroblast 
 growth factor receptor activity  , as shown in figure 7 , genes are more tightly correlated 
 in the proteomic data compared with the rnaseq data .
various studies 
 have investigated the role of the fgfr pathway as a predictive marker 
 for breast cancer . as shown in figure 7 , fgfr2 and fgfr3
cerliani et al . reported a strong correlation between fgfr2 and fgfr3 based on protein 
 expression mentioning that there is no previous evidence of correlation 
 between these two proteins for breast cancer patients .
in addition , 
 the interaction between fgfr2 and fgfr3 is contained in the string 
 database ( see table s2 ) .
similarly to 
 module c1 , genes in module c2 are higher correlated in the proteomic 
 data compared with the rnaseq data ( figure 7 ) .
this cluster contains genes such as apob , 
 c3 , orm1 , and cp that have been recently shown to be predictive of 
 chemoresistance in breast cancer patients .
another 
 important gene in this module with a high number of connecting edges 
 is apoa1 ( the fifth highest connected node according to figure 5 ) .
this gene has been recently 
 identified as a potential target for disease progression using whole - genome 
 trancriptomic and whole proteomic .
as 
 shown in table s2 , 11 of the relationships 
 in figure 7 are contained 
 in the string database .
these results suggest that proteomic data 
 can provide complementary information to genomics data and help to 
 reveal important biological mechanisms .
, we show the most 
 enriched go category with corresponding benjamini adjusted p - value ( p ) .
in particular , only modules 
 reporting a significant enriched benjamini adjusted p - value ( p  0.01 ) are highlighted .
heatmap of the absolute value of correlation 
 between genes based 
 on proteomic and rnaseq data for module c1 ( a ) and c2 ( b ) in figure 6 . for each module
, 
 we include a pie plot showing the number of protein - specific interactions 
 ( red ) and interactions shared across proteomic and rnaseq data ( green ) .
because we are interested 
 in assessing the utility of a joint learning , jrf was compared with 
 the standard random - forest algorithm that constructs the two networks 
 separately . for both methodologies , we derived undirected networks 
 using the same fdr cutoff ( 0.001 ) .
the protein - network resulting from 
 genie3-sep can be found in figure s11 in the supporting information .
as shown , 57% of the edges are shared across protein- 
 and rna - networks under jrf , while only 26% are shared under genie3-sep . 

this finding can be explained by the fact that a joint learning can 
 facilitate the detection of common edges . in this section
, 
 we validate our findings using go categories . in particular , we consider 
 596 go terms containing fewer than 200 genes ( the list of go terms 
 can be found in tables s3s6 in the supporting information ) .
figure 8a shows the number of edges contained in at least one go term 
 for different fdr cutoff for both protein - networks .
first , we constructed 
 a network based on go terms where an edge ( a  b ) was drawn between every pair of genes ( a , b ) sharing at least one go term .
then , a contingency 
 table was constructed considering variables  network based 
 on go terms  and  network based on proteomic data  
 with categories
 number of edges contained  and  number 
 of edges non - contained  , and a fisher s exact test was 
 performed to derive enrichment p - values .
figure 8a , b shows that jrf 
 results in a network more overlapping with go terms than the standard 
 random - forest algorithm genie3-sep .
we validate protein - networks in figure 4 using go terms and 
 knockout signatures . 

the first row contains network validation based on go terms . for our 
 analysis , we considered 596 go terms containing fewer than 200 genes . 

( a ) number of edges contained in at least one go term for different 
 fdr cutoff .
the second row contains validation based on esr1 ( c ) and 
 gata3 ( d ) signatures . for a particular neighborhood size of either 
 esr1 or gata3
, we show the total number of knockout signatures contained 
 under jrf (  ) and genie3-sep (  ) . to further validate 
 our findings , we overlapped both networks with 
 esr1 and gata3 knock - down signatures .
we identified sirna knock - down 
 signatures of key transcription factors ( tfs ) of breast carcinoma 
 in mcf7 cells1 by accessing gene expression data available under gse31912 
 in gene expression omnibus ( further details 
 about these signatures can be found in section 2 of supporting information ) . for each network ,
a neighborhood 
 was defined by defining a threshold k for the number 
 of edges connecting a particular gene to the target gene of interest 
 ( either esr1 or gata3 ) . for each network in figure 8
, we considered different values of k and counted the number of knockout signatures contained 
 in the neighborhood .
figure 8c , d shows the total number of signatures contained in k - size 
 neighborhoods for esr1 and gata3 , respectively . as shown , also in 
 this case , jrf results in a more enriched network .
figure s12 contains the same comparison for rna networks . in this paper , we developed jrf , a random - forest - based 
 algorithm for the simultaneous construction of gene coexpression and 
 protein coexpression networks .
jrf is designed to borrow information 
 across different expression data by selecting the same set of predictors 
 ( genes ) as splitting variables in the random forest model corresponding 
 to each data . in this way , jrf is able to detect common relationships 
 ( edges ) with better power and detect differential edges - specific to 
 individual data with fewer false positives .
because in this 
 study we were mainly interested in assessing the unique contribution 
 of proteomic data to breast cancer research , we focused our attention 
 on protein - specific hub genes and modules .
we identified two interesting 
 protein - specific modules containing potential targets for breast cancer 
 treatment .
in addition , we compared our algorithm to the original 
 random - forest algorithm , which constructs the two networks separately , 
 and showed that our algorithm leads to networks more overlapping with 
 available knockout signatures and existing go terms .
the current 
 version of jrf requires the feature ( protein ) space 
 to be the same across different classes .
extension of jrf for handling 
 data involving different sets of variables remains as future work . 

this problem arises in the protein domain because protein abundance 
 can be measured for different phosphorylation - sites that map to a 
 unique protein . as future work
, we will design a model able to jointly 
 estimate networks from phosphorylation - site abundance and rnaseq data , 
 simultaneously . besides estimating gene - regulatory network and 
 protein - regulatory 
 network ,
as shown in the paper , jrf can be used in many other applications . 

for example , we can apply jrf to estimate grns for the three breast 
 cancer subcategories ( luminal , basal , and her2 ) , grn for different 
 tissues , and grn for cancer and normal tissues .
another interesting 
 direction would be assessing the association between mirnas - genes 
 and mirnas - proteins simultaneously .
it is well known that the 
 downregulation of some mirnas may play an important role in the progression 
 of cancer .
therefore , a better understanding 
 of the interactions between mirna , mrna , and protein expression is 
 crucial to cast light on the potential disease mechanisms of cancer . 
 in this context , jrf can be easily utilized to jointly detect which 
 mirnas regulate genes and proteins .
jrf can be easily parallelized because the computation can be divided 
 into p - independent subproblems . on the contrary , 
 for many existing methods such as bayesian networks and ggm - based 
 algorithms ,
the merit of jrf depends on the assumption that networks 
 of different 
 classes share some common structures . to assess the performance of 
 jrf in the absence of common structure , we conducted numerical studies 
 to investigate such cases and concluded that jrf is sufficiently flexible 
 to guarantee good performance even when nonrelated networks are considered . in jrf ,
we propose a permutation - based procedure 
 to derive proper cutoffs on importance scores for selecting confident 
 edges .
specifically , we introduce f( ) as an 
 approximation of the false discovery rate ( fdr ) of the selected edge 
 set based on cutoff  ; however , in the numerical studies , we 
 observe that f( ) tends to underestimate fdr . 

new methods to obtain more accurate estimate of fdr for 
 jrf warrant future research .
 we focus on characterizing 
 common and different coexpression patterns 
 among rnas and proteins in breast cancer tumors . to address this problem , 
 
 we introduce joint random forest ( jrf ) , a novel nonparametric algorithm 
 to simultaneously estimate multiple coexpression networks by effectively 
 borrowing information across protein and gene expression data . 
 the 
 performance of jrf was evaluated through extensive simulation studies 
 using different network topologies and data distribution functions . 
 
 advantages of jrf over other algorithms that estimate class - specific 
 networks separately were observed across all simulation settings . 
 
 jrf also outperformed a competing method based on gaussian graphic 
 models . 
 we then applied jrf to simultaneously construct gene and protein 
 coexpression networks based on protein and rnaseq data from cptac - tcga 
 breast cancer study . 
 we identified interesting common and differential 
 coexpression patterns among genes and proteins . 
 this information can 
 help to cast light on the potential disease mechanisms of breast cancer . 