An Efficient and Easy-to-Use Network-Based Integrative Method of Multi-Omics Data for Cancer Genes Discovery

Identifying personalized driver genes is essential for discovering critical biomarkers and developing effective personalized therapies of cancers. However, few methods consider weights for different types of mutations and efficiently distinguish driver genes over a larger number of passenger genes. We propose MinNetRank (Minimum used for Network-based Ranking), a new method for prioritizing cancer genes that sets weights for different types of mutations, considers the incoming and outgoing degree of interaction network simultaneously, and uses minimum strategy to integrate multi-omics data. MinNetRank prioritizes cancer genes among multi-omics data for each sample. The sample-specific rankings of genes are then integrated into a population-level ranking. When evaluating the accuracy and robustness of prioritizing driver genes, our method almost always significantly outperforms other methods in terms of precision, F1 score, and partial area under the curve (AUC) on six cancer datasets. Importantly, MinNetRank is efficient in discovering novel driver genes. SP1 is selected as a candidate driver gene only by our method (ranked top three), and SP1 RNA and protein differential expression between tumor and normal samples are statistically significant in liver hepatocellular carcinoma. The top seven genes stratify patients into two subtypes exhibiting statistically significant survival differences in five cancer types. These top seven genes are associated with overall survival, as illustrated by previous researchers. MinNetRank can be very useful for identifying cancer driver genes, and these biologically relevant marker genes are associated with clinical outcome. The R package of MinNetRank is available at https://github.com/weitinging/MinNetRank.


INTRODUCTION
Rapid technological advances in high-throughput sequencing have driven the development of omics field. Omics data types include genomics, transcriptomics, proteomics, epigenomics, and metabolomics (Hasin et al., 2017). However, a single type of "omics" only provides limited insights into the biological mechanisms of diseases. Additionally, the different omics data events are somewhat interdependent. An integrative study of multiomics data contributes to a holistic understanding of the molecular function (Sun and Hu, 2016). An essential question in cancer genomics is distinguishing driver genes, which are causally implicated in oncogenesis, from biologically neutral passenger genes that are immaterial to neoplasia (Greenman et al., 2007). Passenger mutations can become driver mutations (and vice versa) under changing environmental conditions and selection pressures, increasing the complexity of intratumor heterogeneity (Yap et al., 2012). Accumulating evidence suggests that identifying personalized driver genes is essential for the development of effective personalized therapies and realizing the goals of precision medicine (Dagogo-Jack and Shaw, 2018). A critical but challenging step is to incorporate different omics data in a meaningful and efficient way to discover cancer driver genes and elucidate potential causative changes of cancer (Huang et al., 2017). The main approaches for distinguishing driver genes from passenger genes can be divided into frequency-based methods and network-based approaches.
Frequency-based methods estimate the background mutation rate (BMR) representing the rate of random passenger mutations and identify driver genes that harbor significantly more somatic mutations than BMR (Ding et al., 2008;Pon and Marra, 2015). However, accurately estimating BMR is difficult because of the variability among cancer types, among samples of the same cancer type, and between genomes (Pon and Marra, 2015). Subsequent frequency-based methods, such as MuSiC and MutSigCV, have been developed to correct for one or more of these factors (Dees et al., 2012;Lawrence et al., 2013). Somatic mutations are characterized by a small number of frequently mutated genes and many infrequently mutated genes. Moreover, more than 99.9% of the somatic mutations in tumors are passengers (Vogelstein et al., 2013). It is challenging to identify infrequent or rare driver genes by methods based only on mutation frequency.
Network-based approaches have emerged as promising and powerful methods to detect low-frequency and high-frequency mutated driver genes due to their ability to model gene interactions. For network-based approaches, nodes representing genes and edges are links between two genes if there is an interaction between them (Huang et al., 2017). Network-based methods have been successfully applied to many biomedical fields, such as the discovery of mutation subnetwork (Vandin et al., 2011), prediction of drug-target interaction, and cancer gene prioritization (Bashashati et al., 2012;Chen et al., 2012;Yu et al., 2013). HotNet2 uses a network diffusion model to simultaneously assess the frequency of somatic mutation and the local topology of the interaction network and detects significantly mutated subnetworks (Leiserson et al., 2015). Mutations for Functional Impact on Network Neighbors (MUFFIN) is a method for prioritizing cancer genes accounting for mutation frequency of genes and their direct neighbors in functional network (Cho et al., 2016). Both HotNet2 and MUFFIN use mutation data only without integrating other omics data. DawnRank is a single patient approach to rank potential driver genes based on their impact on downstream differential expression genes in the interaction network (Hou and Ma, 2014). NetICS predicts mediator genes affected by proximal upstream-located aberrant genes and proximal downstreamlocated differentially expressed genes (Dimitrakopoulos et al., 2018). Both DawnRank and NetICS consider only incoming degree or outgoing degree of interaction network for single omics. For example, DawnRank only considers incoming degree for expression data. It is desirable to use incoming and outgoing degree simultaneously. Driver_IRW (Driver genes discovery with Improved Random Walk method) assigns different transition probabilities for different genes of the interaction network (Wei et al., 2020). DeepDriver predicts cancer driver genes based on mutation-based features and gene similarity networks using deep convolutional neural networks . None of these methods consider the different weights for the different types of mutations; however, the weighting method is essential for sample-specific study. Furthermore, none of these methods investigate the relationship between the top rankings of genes and overall survival. Therefore, we develop a more meaningful and efficient method that considers different weight coefficients for the various types of mutations, simultaneously considers the incoming and outgoing degree of interaction network for single omics, and uses minimum strategy to integrate multiomics data.
We present a new method called MinNetRank that uses minimum strategy among multi-omics data to prioritize cancer genes (Figure 1). The main steps of MinNetRank include (1) single-omics data analysis: calculating mutation relevance scores and expression relevance scores of genes for each sample using network diffusion based on incoming and outgoing degree. We further consider different weight coefficients for the different types of mutations and propose Weighted_MinNetRank. (2) The integration of multi-omics data: calculating the minimum value of mutation relevance score and expression relevance score as an integrated score for each gene in each sample. A higher minimum value reflects a higher mutation relevance score and expression relevance score simultaneously; (3) prioritizing driver genes: aggregating the sample-specific and integratedscore-based rankings of genes into a robust populationlevel gene ranking.
We apply Weighted_MinNetRank and MinNetRank to analyze five The Cancer Genome Atlas (TCGA) datasets (hepatocellular carcinoma, stomach adenocarcinoma, bladder urothelial carcinoma, lung adenocarcinoma, and skin cutaneous melanoma) and one International Cancer Genome Consortium (ICGC) dataset (hepatocellular carcinoma). We select the top 50 genes of population-level ranking as candidate driver genes. We systematically examine the performance of Weighted_MinNetRank and MinNetRank from three aspects. Firstly, Weighted_MinNetRank and MinNetRank outperform other methods [Mean, Maximum, DawnRank, NetICS, and a commonly used frequency-based method (Freq)] in terms of precision, F1 score, and partial area under the curve (AUC) value of selecting cancer driver genes. Secondly, Weighted_MinNetRank and MinNetRank detect rare and novel candidate driver genes (e.g., SP1 in hepatocellular carcinoma). Finally, the top seven genes can be used as prognostic biomarkers for risk stratification. The survival difference between two FIGURE 1 | Overview of MinNetRank workflow. MinNetRank integrates mutation data and expression data into the interaction network; MinNetRank utilizes minimum strategy to select the candidate driver genes with both high mutation relevance score and high expression relevance score; the sample-specific and integrated-score-based rankings of genes are integrated into the overall rankings. We access the performance in predicting known cancer genes, discovering personalized driver genes, and survival risk stratification of tumor samples.
subtypes (low-risk and high-risk groups) is statistically significant in all six datasets.

RESULTS
We propose a new method (MinNetRank) that uses minimum strategy among multi-omics data to prioritize cancer genes. For comparison, we also add the performance of mean (Mean) and maximum (Maximum) to integrate the mutation data and expression data. All mutations have the same weight for MinNetRank. We further consider different weight coefficients for the different types of mutations (Weighted_MinNetRank). In this study, Weighted_MinNetRank and MinNetRank are compared with other five methods [Mean, Maximum, DawnRank (Hou and Ma, 2014), NetICS (Dimitrakopoulos et al., 2018), and Freq] on five types of cancer (liver hepatocellular carcinoma, stomach adenocarcinoma, lung adenocarcinoma, bladder urothelial carcinoma, and skin cutaneous melanoma). Freq is a simple and common method based only on mutation frequency, which compares the mutation frequency of genes in tumor patient (Dimitrakopoulos et al., 2018;Guo et al., 2018). Weighted_MinNetRank and MinNetRank are an efficient and easy-to-use network-based method for cancer genes discovery by integrating multi-omics data, as shown in the subsequent results.

Overview of MinNetRank
The schematic in Figure 1 illustrates the three-step procedure of our new method MinNetRank. MinNetRank requires three input files: gene mutations, gene expression for tumor and normal samples, and the interaction network.
Step 1: calculating mutation relevance score and expression relevance score using RWR (Random Walker with Restart) algorithm. The n×m matrix S M is the gene mutation status for each sample, where n is the number of genes, and m is the number of samples. S M ik = 1 if gene i is mutated in sample k and S M ik = 0 otherwise. We further consider different weight coefficients for the different types of mutations and supplement a new method (Weighted_MinNetRank). We normalize each column of S M by S M /colSum( S M ). We define the n×m mutation relevance score matrix W M as multiplication between diffused matrix D and S M : The D ij reflects the connectivity between gene i and gene j, and S M ik reflects the mutation status of gene i in sample k. The product W M ik is gene i's mutation relevance score in sample k, defined as the proximity of gene i to mutation genes.
Similarly, the n×m matrix S E is RNA differential expression score (Absolute value of Log2 Fold-Change, ALFC) for each sample. We define the expression relevance score matrix W E as, Step 2: minimum value of mutation relevance score and expression relevance score. To integrate multi-omics data (gene mutation and expression data), the mutation relevance score and expression relevance score are combined to produce a gene minscore for each sample. The min-score is the minimum value of W M ik and W E ik : pmin is R function and returns the minimum of the corresponding elements of the two input vectors. W ik is the minimum value of W M ik and W E ik (i ∈ 1 · · · n, k ∈ 1 · · · m), where n is the number of genes, and m is the number of samples. The high score of W ik means that gene i is proximal to many mutation genes and differentially expressed genes for each k. The minimum value is a meaningful and efficient way to integrate multi-omics data for the following two reasons: Firstly, the minimum strategy reduces extreme values that may be potential outliers in highly skewed distributions. The probability distribution of W M * k (the mutation relevance scores for genes in sample k) and W E * k (the expression relevance scores for genes in sample k) is a positively skewed distribution. This means that some genes have extremely high scores. These high scores may be due to the technical noise of high-throughput sequencing and the incomplete interaction network. For example, as shown in Figure 2, sample TCGA-BC-A10X has three mutated genes in TCGA-LIHC, and only one gene (OR2C3) of these is in the interaction network. The OR2C3 mutation relevance score in TCGA-BC-A10X is evidently high (W M ik = 0.48, i =OR2C3 and k =TCGA-BC-A10X) and is ranked 1st. Meanwhile, the OR2C3 expression relevance score in TCGA-BC-A10X is 3.24-06 and is ranked 8, 221st. Henceforth, the high mutation relevance score needs to be cautiously processed. Lastly, the min-score of OR2C3 mutation relevance score and expression relevance score is ranked 1, 943rd. OR2C3 is an olfactory receptor protein and probably is not a potential driver gene (Malnic et al., 2004;Riessland et al., 2017).
Secondly, the minimum ("double high") strategy is necessary to prioritize cancer genes having a higher biological relevance. If one gene has a relatively high mutation relevance score but low expression relevance score (such as OR2C3 in TCGA-BC-A10X), this gene may not be a potential driver gene since differential gene expression is the downstream events of DNA mutation (Sager, 1997). In the other case, the SI expression relevance score in TCGA-DD-AAE2 is ranked 8th (W E ik = 0.0012, i =SI, and k =TCGA-DD-AAE2), and the mutation relevance score is ranked last. Only MGAM interacts with SI in the interaction network, and TCGA-DD-AAE2 has no SI or MGAM mutation. We hope the candidate driver genes have a high mutation relevance score and high expression relevance score.
MinNetRank used a minimum strategy to integrate multiomics data (mutation data and expression data). We further investigated which data have the greatest effect on the minimum score. We calculated the proportion of mutation relevance score and expression relevance score in minimum scores for the top 50 candidate cancer genes. The proportion of mutation relevance score was 0.657 in all six datasets, and expression relevance score was 0.347. Mutation relevance score affected the minimum score more.
Step 3: integrating sample-specific rankings of genes into a population-level ranking. We transform the min-scores into rankings, since min-scores indicate the relative importance of each sample's genes. To integrate the sample-specific rankings of genes into a robust population-level ranking, we calculate the sum of per-sample ranking. Each step of MinNetRank is based on single sample analysis, such as using the per-sample network diffusion, calculating the minimum value of mutation relevance score and expression relevance score for each gene in each sample, and transforming min-scores into rankings for each sample. We calculate the sum of per-sample ranking as the population-level ranking.
To perform a systematic comparison of seven methods (Weighted_MinNetRank, MinNetRank, Mean, Maximum, DawnRank, NetICS, and Freq), the 576 genes annotated in cancer gene census (CGC) are used as the gold standard cancer driver gene set, and the genes not in CGC are the negative set. The evaluation metrics (precision, F1 score, and partial AUC value) are based on the top 50 genes of six different datasets (five TCGA datasets and one ICGC dataset). The five TCGA datasets are regarding hepatocellular carcinoma (TCGA-LIHC), stomach adenocarcinoma (TCGA-STAD), bladder urothelial carcinoma (TCGA-BLCA), lung adenocarcinoma (TCGA-LUAD), and skin cutaneous melanoma (TCGA-SKCM), respectively. The one ICGC dataset includes hepatocellular carcinoma data from LIRI-JP (Liver Cancer-RIKEN, JP) project (LIRI-LIHC) (Fa et al., 2019). Skin cutaneous melanoma, lung adenocarcinoma, bladder urothelial carcinoma, and stomach adenocarcinoma have a high mutation burden (Martincorena and Campbell, 2015), and LIHC has two different datasets. Both are common cancer types and pose increasing public concerns. The detailed descriptions of six datasets are provided in Table 1. The somatic mutations include non-synonymous simple nucleotide variation (SNV) and insertions and deletions (InDels) in coding regions.

MinNetRank Accurately Predicted Cancer Gene
In general, considering the weights for the different types of mutations (Weighted_MinNetRank) had a better performance than other six methods (MinNetRank, Mean, Maximum, NetICS, DawnRank, and Freq) in all six cancer datasets (TCGA-LIHC, TCGA-STAD, TCGA-BLCA, TCGA-LUAD, TCGA-SKCM, and LIRI-LIHC). Weighting for the different types of mutations was essential for a personalized analysis. As shown in , Weighted_MinNetRank and MinNetRank achieved a higher precision, F1 score, and AUC in all six datasets, namely, Weighted_MinNetRank and MinNetRank could rank the known gold standard cancer driver genes higher. The AUC of Freq was not calculated as the mutation frequency for some genes were the same.

MinNetRank Robustly Predicted Cancer Gene
The Weighted_MinNetRank and MinNetRank also had the advantage of obtaining robust and stable results using the subset Frontiers in Genetics | www.frontiersin.org FIGURE 2 | The heatmap of rankings of mutation relevance scores, expression relevance scores, and min-scores. Each row represents a gene. Ranking differences are shown in different colors. Red means high ranking (high score), and blue means low ranking. The rankings of genes are ordered by the rankings of min-scores. The left is the enlarged drawing of the top 50 genes with both high mutation relevance score and high expression relevance score. OR2C3 has high mutation relevance score (red) and low expression relevance score (blue). The ranking of OR2C3 for mutation relevance score, expression relevance score, and min-score in TCGA-BC-A10X is 1, 8,821, and 1,943, respectively.
of samples with different sample sizes. We calculated the mean and standard deviation (SD) of the precision values P (mean precision of the top 50 genes), F1 scores, and partial AUC values after 10 runs. The precision value was proportional to the area under the precision curve ( Figure 3A). All six methods used the same subset of samples, and the subset of samples was randomly selected from all samples by R. Using the same subset of samples, we compared the results of six methods. The mean of the precision, F1 score, and partial AUC for Weighted_MinNetRank and MinNetRank was higher than other methods, and the SD was smaller [ Furthermore, Weighted_MinNetRank and MinNetRank still performed well, even with a smaller number of samples.
In order to evaluate the contribution of each part of Weighted_MinNetRank and MinNetRank (calculating the relevance score using both incoming and outgoing degree of the interaction network for single omics, using minimum strategy to integrate multi-omics data, and the different weighted methods), we calculated the precision, F1 score, and partial AUC value of the top 50 candidate cancer genes. We also added network metrics (degree centrality, betweenness centrality, and the mean of degree and betweenness centrality). We needed to calculate the baselines of the network only once, and the results were the same for all datasets. As shown in Table 2, Weighted_MinNetRank had a better performance than all other methods in terms of precision, F1 score, FIGURE 3 | Comparison of precision, F1 score, and AUC for different methods in TCGA-LIHC and LIRI-LIHC datasets. (A) The X-axis is the top 50 candidate cancer genes, and the Y -axis is the precision according to known cancer genes (in CGC). (B) The X-axis is the top 50 candidate cancer genes, and the Y -axis is the F1 score according to known cancer genes. (C) The ROC curve of the top 50 candidate cancer genes. The X-axis is the subset of samples, and the Y -axis is the mean and SD of the precision values after 10 runs using the subset of samples. (B) The X-axis is the subset of samples, and the Y -axis is the mean and SD of the F1 score after 10 runs using the subset of samples. (C) The X-axis is the subset of samples, and the Y -axis is the mean and SD of the partial AUC after 10 runs using the subset of samples. and partial AUC in all six datasets. For weighted methods, Weighted_MinNetRank_PrCID had better performance than PrDSM weighted methods (Weighted_MinNetRank_PrDSM and Weighted_MinNetRank_Filter_PrDSM) in all datasets. There was no significant difference between Weighted_MinNetRank_PrCID and Weighted_MinNetRank. There were some possible reasons for this phenomenon. Firstly, there were many synonymous mutations in all datasets (32,381 synonymous mutations on average); however, the percentage of deleterious synonymous mutations was relatively small (9.76% in the study of PrDSM) (Cheng et al., 2019). Many benign synonymous mutations increased noise. We may need to pre-process the scores of synonymous mutations (Weighted_MinNetRank_Filter_PrDSM performed better than Weighted_MinNetRank_PrDSM). Secondly, the number of missense mutations was the largest, and the number of frameshift mutations was small, so Weighted_MinNetRank weighting for missense mutations had almost the same performance as Weighted_MinNetRank_PrCID weighting for missense mutations and frameshift mutations. LIRI-LIHC dataset did not provide the position information of frameshift mutations in cDNA, so Weighted_MinNetRank_PrCID was not available for LIRI-LIHC dataset.

MinNetRank Discovered Rare and Novel Driver Genes
In addition to obtaining the accurate and robust results, one of the main advantages of MinNetRank was to discover rare and personalized cancer genes. Personalized driver genes could contribute to the development of personalized medicine.
As shown in   We further validated SP1 expression on the protein level, and the differential protein expression between tumor and normal samples was significant in LIHC (Wilcoxon Signed Rank test, P = 4.14e-13). Only LIHC had protein expression data from CPTAC (The National Cancer Institute's Clinical Proteomic Tumor Analysis FIGURE 5 | The SP1 differential expression between tumor and normal samples. From left to right: SP1 RNA differential expression in TCGA-LIHC dataset, SP1 RNA differential expression in LIRI-LIHC dataset, and protein differential expression for LIHC from CPTAC dataset. Consortium) dataset. These results suggested that SP1 can be the biomarker of hepatocellular carcinoma.

Top Genes of MinNetRank Were Associated With Clinical Outcome
For each dataset, we selected seven genes with top ranking and high SD as biomarkers for tumor stratification (mentioned in the section "Materials and Methods"). We performed unsupervised K-means clustering using obtained biomarkers to assign each patient into either high-risk or low-risk groups. The Kaplan-Meier survival curves of the two groups are well separated, and the log-rank P-values of the survival difference between two groups are 9.21e-04, 1.23e-05, 2.42e-03, 3.75e-03, 9.21e-04, and In the two liver cancer datasets (TCGA-LIHC and LIRI-LIHC), there were six shared genes (CTNNB1, JUN, PIK3R1, RAC1, SRC, and TP53). All these genes used for tumor stratification are biologically relevant. CTNNB1 regulated cell growth and adhesion and was predictive for recurrence in aggressive fibromatosis (van Broekhoven et al., 2015). JUN (AP-1 Transcription Factor Subunit) participated in regulating a diverse array of cellular processes, including proliferation, apoptosis, differentiation, and survival (Trop-Steinberg and Azar, 2017). PIK3R1 was a prognostic biomarker for breast cancer (Cizkova et al., 2013). RAC1 regulated a wide range of cellular events, including the control of cell growth and the activation of protein kinases (Lou et al., 2018). SRC was prognostic relevant to colon cancer and rectal cancer (Martínez-Pérez et al., 2017). TP53 was one of the most frequent alterations and potential prognostic markers in human cancers (Olivier et al., 2010). GRB2 was the special biomarker for TCGA-LIHC, and MAPK14 was for LIRI-LIHC. GRB2 was evaluated as a prognostic marker for lung adenocarcinoma (Toki et al., 2016). MAPK14 was a member of the MAP kinase family. MAPK pathway regulated cell proliferation, differentiation, and development (Fang and Richardson, 2005). The seven biomarkers are the same in TCGA-STAD and TCGA-BLCA (CTNNB1, GRB2, JUN, RAC1, SP1, SRC, and TP53). These seven genes were reported to be related to prognosis (Hang et al., 2016). For TCGA-LUAD and TCGA-SKCM, there were six shared genes (CTNNB1, JUN, RAC1, SRC, TP53, and GRB2). GNB1 was the special biomarker for TCGA-LUAD, and FYN was for TCGA-SKCM. FYN was tyrosine kinases and was an essential molecule in cancer pathogenesis and drug resistance (Elias and Ditzel, 2015). In summary, the top seven genes were associated with clinical outcome and were biologically relevant in all six datasets. These results suggested that MinNetRank could also be a promising method for tumor stratification.
NetICS and DawnRank did not investigate the prognostic value of top genes in cancer. To evaluate the performance of predicting the clinical outcome for different methods, we used the same criterion to choose the top seven genes for each method in six datasets. Compared with NetICS and DawnRank, only Weighted_MinNetRank and MinNetRank obtained a statistically significant survival risk difference between the high-risk and low-risk groups in all six datasets (Supplementary Table 2).

DISCUSSION
Extensive genetic heterogeneity exists between tumors of different tissues and between individuals with the same tumor type (Burrell et al., 2013). The personalized mutation profile is the key to advance personalized disease diagnosis and therapy in the clinic (Sheng et al., 2015;Olivier et al., 2019). However, few methods could efficiently prioritize driver genes over many passenger genes in a single patient. The critical challenge facing today is to predict rare and even personalized driver genes with higher accuracy. We develop MinNetRank, an efficient and easyto-use method that integrates the mutation data, expression data and interaction network to prioritize each sample's driver genes. Weighted_MinNetRank further considers the different weights for the different types of mutations.
Weighted_MinNetRank and MinNetRank achieve a higher precision, F1 score, and partial AUC value of prioritizing cancer genes in five TCGA datasets (TCGA-LIHC, TCGA-STAD, TCGA-BLCA, TCGA-LUAD, and TCGA-SKCM). We also utilize an additional liver cancer cohort (LIRI-LIHC) to validate the result of TCGA-LIHC. Better performance in all datasets demonstrates the proposed approach's robustness (Figure 3 and Table 2). We use top candidate driver genes for pathway enrichment analysis and find some signaling pathways previously studied in cancer, such as the Ras signaling pathway and ErbB signaling pathway. Furthermore, we first investigate the relationship between the top seven genes and clinical outcome and find the statistically significant survival difference between the low-risk and high-risk groups in all six datasets only for Weighted_MinNetRank and MinNetRank. The top seven genes are biologically relevant and could be used as biomarkers for survival risk stratification. Accurate outcome prediction is important for personalized cancer therapies in clinical practice, for instance, a low-risk patient can be advised to select a less radical therapy.
We demonstrate that MinNetRank can discover rare and novel cancer genes. Personalized driver genes could contribute to developing personalized diagnosis and therapy. SP1 is considered a candidate driver gene only by MinNetRank and is ranked top three in all six datasets. The RNA expression of SP1 is significantly higher in LIHC tumor samples (TCGA-LIHC and LIRI-LIHC datasets) and STAD tumor samples (TCGA-STAD dataset). The differential expression is further validated on the protein level in LIHC. SP1 is the biomarker for tumor stratification in TCGA-STAD and TCGA-BLCA, and SP1 RNA expression is associated with survival outcome in TCGA-STAD dataset (Cox proportional hazards model, P = 0.02). These results are in accordance with the reports in literatures (Shi and Zhang, 2019). Targeting SP1 is highly promising strategy in cancer chemotherapy (Vizcaíno et al., 2015).
Using both the incoming and outgoing degree of interaction network, the minimum strategy and weighting for the different types of mutations all contribute to the accuracy and robustness of prioritizing driver genes. Known cancer genes have a higher incoming and outgoing degree, and simultaneously considering incoming and outgoing degree is rational. MinNetRank adopts a minimum strategy to prioritize cancer genes with a high mutation relevance score and high expression relevance score. These enable our method to select more relevant genes and avoid the potential outliers, which are common in highthroughput sequencing technologies due to the positively skewed distributions of mutation and expression relevance scores. Weighting for different types of mutations is essential for samplespecific study and finding personalized driver genes.
There are some limitations to MinNetRank and similar methods. Firstly, MinNetRank largely depends on the interaction network. Although many interaction sources exist, such as experiment, co-expression, and text mining, the interaction network is still incomplete. If the mutation gene or differentially expressed gene is not in the interaction network, this gene would not be used for network diffusion and not be as a candidate cancer gene. Secondly, MinNetRank uses paired tumor and normal samples to calculate ALFC; however, TCGA datasets have a limited number of normal samples with expression data. Thirdly, MinNetRank only integrates mutation data and expression data into the interaction network. Besides mutation data, other events, such as miRNA differential expression, epigenetic changes, copy number variation, and structure variation, could also contribute to cancer progression. Differential expression data, including RNA expression data and protein expression data, could be combined. We may need to improve MinNetRank from two aspects in the future. On one hand, we could integrate the gene co-expression network with the interaction networks (Hou et al., 2019;Wei et al., 2020). We also need to incorporate additional types of omics data (genomics, transcriptomics, proteomics, epigenomics, and images). On the other hand, Weighted_MinNetRank only considers mutations in coding region. We may need to incorporate non-coding mutations. We also need to give weight coefficients for all mutations through multiple techniques.
Integrating different types of omics data is often used to better elucidate the molecular function. However, sound study designs and solid analytical strategies are needed to advance human disease research further. For example, the mean precision of the top 50 cancer genes is 0.61 (MinNetRank) and 0.56 (NetICS) in TCGA-LIHC and 0.61 (MinNetRank) and 0.54 (NetICS) in TCGA-BLCA. The top 50 candidate cancer genes of NetICS used here are from the published paper (Dimitrakopoulos et al., 2018). In this article, NetICS integrates different types of data that include somatic mutation, copy number variation, methylation, miRNA expression, gene expression, and protein expression. Although MinNetRank only focuses on integrating the mutation data and expression data, the mean precision of MinNetRank is still higher than that of NetICS.

CONCLUSION
This article developed a new method (denoted as MinNetRank) by setting weights for different types of mutations and using the minimum strategy to integrate multi-omics for cancer genes discovery. Minimum strategy reduced the influence of extreme scores in highly skewed distributions and was the "double high" strategy to prioritize cancer genes, having a relatively high mutation score and expression score. Different weight coefficients for the different types of mutations contributed to the better performance. We demonstrated our method's accuracy and robustness in prioritizing driver genes on five TCGA datasets and one ICGC dataset. Besides, MinNetRank has the advantage of discovering rare and personalized cancer genes. The top seven candidate driver genes stratified patients into two subtypes (high-risk and low-risk groups) exhibiting significant survival differences and could be used as prognostic biomarkers for survival. Of course, our method has room for improvement. Gene co-expression network and more types of omics data should be incorporated, and different weight coefficients should be considered.

Dataset
The genes annotated in the CGC can be used to benchmark known cancer genes (Tate et al., 2019). This gold standard known cancer gene set includes 576 genes (July 2019) 1 . Many cancer

Interaction Network
We used the interaction network that has been widely used in the related paper (Hou and Ma, 2014;Guo et al., 2018). The interaction network integrated a variety of resources, including the network used in MEMo as well as the up-to-date information from Reactome (Croft et al., 2011;Ciriello et al., 2012), the NCI-Nature Pathway Interaction Database (Schaefer et al., 2009), and KEGG (Kanehisa et al., 2016). The resulting interaction network consisted of 11,648 genes and 211,794 edges. The average degree centrality of interaction network was 34.20, and the average betweenness centrality was 1.58E-04.

MinNetRank
MinNetRank uses an interaction network that could discover cancer driver genes more efficiently (Leiserson et al., 2015). One of the main reasons for this is the high connectivity (high incoming degree and outgoing degree) of known cancer genes in the interaction network. For example, the mean and median of incoming degree for known cancer genes (in CGC) are 36.06 and 17, which are much higher than those of the genes that are not classified as known cancer genes (17.41 and 3, respectively). Also, the mean and median outgoing degree of known cancer genes are 30.37 and 12, which are much higher than those of the genes that are not in CGC (17.66 and 4,respectively). To a certain extent, this is expected since genes with high connectivity could exert a more significant influence on the biological system (Winter et al., 2012). RWR algorithm models how closely related the two genes are and measures both the direct and indirect neighbors of each gene in the interaction network, making it more sensitive for prioritizing cancer driver genes (Dimitrakopoulos et al., 2018). Unlike NetICS and DawnRank, we consider both incoming and outgoing degree of interaction network for single omics.

Diffused Matrix
Let A be the n×n adjacency matrix of an interaction network where n represents the number of nodes (the number of genes in the interaction network). A is a 0-1 matrix and a ij = 1 if there is a directed edge from node j to node i . A is the transpose of matrix A and a ji = 1 if there is a directed edge from node i to node j . We denote deg out j = N i=1 a ij as the outgoing degree of node j or the number of outgoing edges. While deg in j = N i=1 a ji is the incoming degree of node j . MinNetRank considers both the incoming degree and outcoming degree, so we define the normalized adjacency matrix A norm as, We define the diffused matrix D as, The value of D ij lies between 0 and 1 and reflects the connectivity between nodes j and i. Higher score means that two genes are more closely related. The restart probability of β(0 ≤ β ≤ 1) determines the degree of diffusion, namely, how far the random walker can move in the network. When β=1, there is no diffusion, namely, we do not use the information of the interaction network. When β=0, gene mutation score or differential expression score (see below) diffuses to the whole network. β depends on the interaction network and is independent of any mutation data or expression data. We chose β to balance diffusion and retainment (Leiserson et al., 2015), and β is 0.48 in this study. The diffused matrix D needs to be computed only once for a given interaction network.

ALFC
For each patient k , we calculate the Absolute value of Log2 Fold-Change ( ALFC ) of gene i for the paired tumor and normal samples as a differential expression score. The fold change, or relative difference, is widely used to measure differential gene expression (Love et al., 2014). The absolute value of fold change is taken in order to capture both upregulation and downregulation.

Weighted_MinNetRank
Weighted_MinNetRank uses SIFT scores (between 0 and 1) as the weight coefficients for missense mutations and gives the same weight with 1 to other mutations (stop-gain, stop-loss, frameshift, and non-frameshift) (Ng and Henikoff, 2001). Although synonymous mutations do not alter amino acids, some deleterious synonymous mutations play important roles in cancer (Wen et al., 2016). We further incorporate synonymous mutations and use PrDSM scores as the weights for synonymous mutations (Weighted_MinNetRank_PrDSM). We also use PrDSM scores greater than 0.38 as the weights (Weighted_MinNetRank_Filter_PrDSM). If a PrDSM score is greater than 0.308, the corresponding synonymous mutation is considered as deleterious (Cheng et al., 2019). Besides, we use PredCID scores as the weights for frameshift mutations (Weighted_MinNetRank_PrCID) (Yue et al., 2020).

Assessing the Performance in Predicting Known Cancer Genes
In order to assess the performance in predicting known cancer genes, our method (Weighted_MinNetRank and MinNetRank) was compared with NetICS (Dimitrakopoulos et al., 2018), DawnRank (Hou and Ma, 2014), and Freq. The top 50 genes of the population-level ranking were identified as candidate driver genes and compared with the positive genes in CGC. We used the precision, F1 score, and partial AUC value to evaluate the performance. The precision was defined as expression (7) and can be viewed as the measure of exactness. The recall was the percentage of total known cancer genes correctly predicted by MinNetRank. F1 score combined recall and precision using the harmonic mean. There were many more negative genes than positive genes (positives/negatives = 0.052) and even fewer positive genes when we considered cancer type-specific known cancer genes (positives/negatives ≈ 0.0029). It was more informative to use partial AUC, which considered the number of true positives scored higher than the nth highest scoring negatives, measured for all values from 1 to n (Dimitrakopoulos et al., 2018). Precision, F1 score, and partial AUC were based on the top 50 genes. precison = (CGC genes) ∩ (Top N predicted driver genes) Top N predicted driver genes .
recall = (CGC genes) ∩ (Top N predicted driver genes) CGC genes . (8) F1 Score = 2 × precision × recall precision + recall . (9) where T was the total number of known cancer genes in CGC, and T i was the number of positives scored higher than the ithhighest scoring negatives.

Assessing the Robustness Using the Subset of Samples
In order to further compare these methods, we calculated the precision, F1 score, and partial AUC using the subset of samples with different sample sizes. We experimented with sample sizes of n = 10, 25, 50 * 1, 50 * 2, ... , 50 * N/50 , and N was the total sample size of multi-omics data. For each sample size, we performed 10 random samples. We defined the precision value P= mean(p i ), where p i was the precision of top i candidate cancer gene, i = 1, 2, ... , 50. The mean and SD of precision value, F1 score, and partial AUC value for 10 runs were used to measure the robustness.

Tumor Stratification
Some papers used gene mutation data and expression data to identify genes that were indicators for survival. Using these biomarkers, patients can be stratified into subtypes (Haider et al., 2014). We further investigated the relationship between the top genes of population-level ranking and patients' survival time.
Genes whose expression with a low variation between tumors provided very limited information for tumor stratification (Winter et al., 2012). According to the genes' rankings, we selected the top seven genes with a greater SD of expression than five as biomarkers for each dataset (Winter et al., 2012). Using these seven biomarkers, K-means clustering (unsupervised learning algorithm) assigned each patient to one of the two clusters (high-risk and low-risk groups). The log-rank test was then used to compare the survival differences of the two groups (R survival package).

DATA AVAILABILITY STATEMENT
The mutation data, expression data, and clinical data of the TCGA dataset are available in the TCGA Data Portal (https://portal.gdc.cancer.gov/projects/). Those from LIRI-JP are available in ICGC Data Portal (https://dcc.icgc.org/projects/ LIRI-JP). The LIHC protein expression data are from CPTAC Data Portal (https://proteomics.cancer.gov/data-portal). The detail descriptions of these data are provided in