BACKGROUND
cancer is driven by genetic alterations, including single nucleotide variants , small insertions or deletions, large copy-number variations  and structural aberrations that accumulate during the lifetime.
several international large scale cancer genomics projects, such as the cancer genome atlas , and international cancer genome consortium   <cit> , etc., have produced a large volume of data in recent years  <cit>  and provided us with an unprecedented opportunity to better characterize the molecular signatures of human cancers  <cit> .
however, it is still a challenge to integrate information across the different omics data  <cit>  and distinguish driver mutations which can promote the cancer cell to proliferate infinitely and diffuse from passenger mutations whose changes represent neutral variation that does not influence cancer development .
in response to the large volume of mutations being generated from massively parallel sequencing projects, many growing mathematical and statistical approaches to search for driver genes, driver pathways or core modules based on data integration were proposed.
the most basic approach, eg.
mutsig  <cit>  and music  <cit> , is to identify driver genes based on somatic mutation rates in cancer patient populations, that is, the most commonly occurring mutations are more likely to be drivers  <cit> .
also, computational approaches based on evaluating the functional impact of mutations  <cit>  such as polyphen- <dig>  <cit>  and oncodriverfm  <cit>  were proposed.
however, cancer is more closely related with a group of genes interacting together in a gene-gene interaction network.
with the advent of the whole-genome measurements of somatic mutations and cnvs in the mass of cancer samples, many changes altered at network and pathway levels are found, not simply a point mutation  <cit> .
therefore, network- and pathway-based approaches have become one of the most promising methods to prioritize driver mutations and significantly mutated genes due to their abilities to model gene-gene interactions.
varwalker is a network-assisted method to prioritize potential driver genes  <cit> .
another method, dawnrank prioritizes altered genes on a single patient level using pagerank algorithm  <cit> .
drivernet is an integrated analysis framework to identify likely driver mutations by virtue of their effect on mrna expression networks and reveals the prevalence of rare candidate driver mutations  <cit> .
it has been demonstrated that genes which are relatively long compared to the distribution of all human consensus coding sequences  are more likely to mutate while they may be not driver genes  <cit> .
however, drivernet doesn’t consider the effect of gene length.
also, the scale of the network in drivernet is a little small which may miss some genes and the information between genes.
in this work, we develop a network-based method called length-net-driver  to improve the performance of detecting driver genes based on the rationale of drivernet  <cit> .
our goal is to consider the point mutation genes’ length and construct a new interaction network contained more genes and interactions based on human protein reference database   <cit>  instead of its original gene influence graph in drivernet.
furthermore, we integrate somatic snv data, cnvs data and gene expression data using gene-gene interaction network.
then a greedy algorithm is applied to the integrated data to prioritize candidate genes.
the application on three tcga datasets demonstrated that the performance of our method is good.
methods
the overview of lndriver approach
in lndriver method, the population-based genomic and transcriptomic interrogations of tumor types were integrated to identify driver mutations.
the pipeline is shown in fig.
<dig> fig.
<dig> schematic of the lndriver.
genes in somatic mutations are firstly applied to gam to filter long genes and then they will combine with cnv to construct mutation matrix.
the bipartite graph is constructed based on mutation data, expression data and gene-gene influence network, where the blue nodes on the left bipartite graph represent the mutated gene and the black nodes on the right represent the outlying patient-gene events from the gene expression matrix.
then greedy algorithm is applied to identify candidate driver genes.
finally, enrichment analysis is employed to these candidates to explore their roles in pathways
actually, some studies have indicated that genes with long length have a better chance to harbor mutations   <cit> .
it indicated that gene length-based filtering process is essential to perform.
hence, in this study, the generalized additive model  was used to assign the somatic mutation probabilities of all human genes for each sample.
then a resampling test was performed to filter passenger genes whose occurring frequencies are ≥ 5% at random datasets  <cit> .
after the filtering procedure, cnvs are combined with it to construct a binary mutation matrix.
in addition, in order to enrich the information of the gene-gene interaction network, we constructed a new interaction network using human protein reference database   <cit> .
as for gene expression data, we built a binary outlying matrix by nominating genes whose expression values are outside two standard deviation of the gaussian distribution as outliers  <cit> .
next, we formulated associations between mutation and gene expression data using a bipartite graph where the left partition of nodes represented the mutation status and the right partition of nodes represented the outlying status in each of patients.
after the above process, greedy algorithm was applied on the bipartite graph to select those genes in the left partition which have the highest number of outlying expression events, and then nominated them as putative driver genes.
also, the statistical significance test was assessed using a randomization framework.
finally, pathway enrichment analysis was done using the database for annotation, visualization and integrated discovery  online tools  <cit> .
to demonstrate the advantages of the approach, we analyzed three large-scale publicly available genome-transcriptome datasets in head and neck squamous cell carcinoma , thyroid carcinoma  and kidney renal clear cell carcinoma .
filtering long genes
the length of genes in human are very different and so the mutation probabilities of different genes are in vast difference.
there may be some genes which have mutations only because they are long yet they aren’t driver mutations.
so, for each gene, we adopted the filtering strategies of varwalker and computed a probability weight vector  by fitting a generalized additive model for each sample  <cit> .
denoting the vector x as the gene length of cdna, we can adopt the following model to assess the mutated probability of a gene according to its cdna length,  <dig> logitπ=logπ1−π~fx where π=#mutantgenes#ccdsgenes represents the proportion of mutant genes  in the researched samples, and f represents an unspecified smooth function  <cit> .
after the above fitting process, each gene was assigned a weight value which would be used to select genes in the next resampling procedure.
then a resampling test was applied to random gene sets for each sample.
the number of being selected random gene sets is same with mutant genes in specific sample.
and the probability of each gene to be selected is based on the probability weight calculated in the above fitting procedure.
the test was performed  <dig> times in each sample following pwv.
the mutation frequency was calculated for each gene using formula :  <dig> fre=#selectingthegeneinresamplingtest <dig> where #  indicates the times and fre represents the frequency of the gene being selected across  <dig> times in resampling process.
then we filtered those genes whose frequencies were ≥ 5% that indicates the gene may occur at random unless they are cgc genes.
those genes with fre < 5% which represented the gene was unlikely mutated at random were observed.
greedy algorithm
for detecting the candidate driver genes based on processed mutation data and expression data, they were integrated with the gene-gene interaction network into a bipartite graph .
the elements on the left of bipartite graph represent the mutation status of genes in population level.
and the right partition events indicate outlying expression status of the genes  <cit> .
an edge between g
i and g
j will be drawn if the gene g
i in the left partition is mutated , the right gene g
j is outlying expression gene  and g
i interacts with g
j in the gene-gene interaction network.
given the bipartite framework, the aim is to find the mutation genes on the left partition which cover the most events on the right of bipartite graph.
to this end, the optimization method of a greedy algorithm was used to select the most covered genes: at each step, chose a mutated gene which connected to the most uncovered outlying expression genes on the right of bipartite graph.
when all the connected outlying expression events were covered, the program was terminated.
finally, the mutated genes ranked based on their coverage and the mostly covered mutated genes are considered as the candidate driver genes.
significance test
in order to assess the statistical significance of the candidate driver genes, the random framework was used by permuting n =  <dig> times of the original datasets including mutation matrix, processed outlying expression matrix and the gene-gene interaction network.
then the algorithm was run on the n randomly generated datasets.
finally, the real data results were assessed to see whether they are significantly different from the results on randomized datasets.
the null hypothesis h
<dig> is that the gene mutations have no influence on the occurrence of the cancer, and the alternative hypothesis h
<dig> is that the cancer is related to the mutations of the genes.
the definition of the statistical significance of gene g, whose corresponding node coverage is cov
g, is the fraction times of selecting driver genes that are more than cov
g in n =  <dig> random runs of the method.
the calculation is listed as follows:  <dig> p−valueg=∑i=1n∑j=1siδcovgij>covg∑i=1nsi where s
i is the number of candidate driver genes selected in the ith run of the method  <cit> .
then the benjamini-hochberg method was used to correct the p-values for multiple tests and finally we chose the genes whose p-values were less than  <dig> .
RESULTS
datasets and pre-processing
we applied lndriver to  <dig> thca samples,  <dig> hnsc samples and  <dig> kirc samples .
these three datasets comprise somatic snv data, cnv data and gene expression data collected from the cancer genome atlas  data portal  <cit> .table  <dig> description of datasets
the construction of mutation matrix
firstly, we collected somatic snvs in level  <dig> and cnv data in level  <dig> directly from tcga data portal.
secondly, we removed the genes whose item of “variant_classification” is “silent” or “rna” in somatic snv data and whose length are too long according to generalized additive model and resampling test process.
thirdly, the cnv information was extracted by selecting genes from amplified and deleted segments in cnv data.
finally, we integrated cnv data with filtered somatic snv data by getting intersecting samples and union genes to construct a binary matrix m, whose rows indicate samples and columns indicate genes.
each entry of m
ij refers to the mutation status of gene j in sample i and m
ij =  <dig> represents that there is labeled valid mutation in gene j of sample i. otherwise, m
ij =  <dig> indicates the absence of a mutation in the jth gene of the ith sample.
expression outlier matrix
for gene expression dataset e, the values of it contain not available  values.
these values affect the results of the approach.
we substituted them with the mean of all other genes in the specific samples.
also, we adopted the assumption in drivernet that the expression distribution of every gene across all samples is gaussian distribution  <cit> .
based on the hypothesis, we converted the expression data to a binary patient-outlier matrix e ' where e
' =  <dig> means the expression of gene i is an outlier in patient j. the definition of the outliers is that genes whose expression values are outside the two-standard deviation range of the expression values of gene i across all the patients  <cit> .
gene-gene interaction network and gene annotation data
cancer is a disease related with sets of genes which interact with each other in some molecular networks not only related with single gene.
in order to enrich the information gene-gene interaction network in drivernet, we built an influence graph g using hprd  <cit>   which contains  <dig> proteins to server as our reference network.
the influence graph g in our work is an undirected and unweighted binary network where v represents the nodes of genes and e represents the edges among genes.
when there is a correlation between gene i and gene j, g
ij =  <dig>  otherwise g
ij =  <dig>
we used the consensus coding sequences  genes data which have been allocated complementary dna  length based on their coding sequences from varwalker  <cit>  as a benchmark gene resource to select those genes that have matched ccds symbols.
in order to explore the impact of the gene length, we compared genes with somatic snvs with the distribution of all human ccds gene length to filter long genes.
cancer gene census  genes
the cgc is a database that catalogues genes whose mutations have been causally implicated in cancer, which has been widely served as benchmark in many cancer researches.
in this work, we also utilized it as the standard reference list which was downloaded from cosmic  <cit>  and included total of  <dig> genes .
the analysis of the overall performance
in this study, the performance of lndriver’s ability was evaluated using the number of indentifying known drivers in cgc database compared with other methods.
the benchmarks of the above evaluation were precision, recall and f1score which were based on the top n genes as following:  <dig> precision=#mutatedgenesincgc∩#genesfoundinlndrivers#genesfoundinlndrivers
<dig> recall=#mutatedgenesincgc∩#genesfoundinlndrivers#mutatedgenesincgc
<dig> f1score=2×precision×recallprecision+recall
for the sake of performing the property of our method on identifying cancer related drivers, we compared the result of our method to classical frequency-based method, generank method  <cit> , drivernet method and personal-based method of dawnrank.
the results of the experiment on hnsc, kirc and thca datasets are shown in fig.
<dig> fig.
<dig>
a hnsc precision.
b hnsc recall.
c hnsc f1score.
d kirc precision.
e kirc recall.
f kirc f1score.
g thca precision.
h thca recall.
i thca f1score.
the comparison of precision, recall and f1score for top ranking genes in lndriver and other methods.
the x axis represents the number of top ranking genes and the y axis represents the score of the precision, recall and f1score respectively
hnsc, the sixth most common cancer worldwide  <cit> , was analyzed in our method.
as for the overall performance of its top  <dig> genes, it can be seen in fig.
2a-c that lndriver method remarkably outperforms other four methods.
for the top  <dig> genes, there are  <dig> genes contained in cgc database of our method, while  <dig> of dawnrank and  <dig> of drivernet.
there are  <dig> genes being selected as candidates and  <dig> genes of them with p-values less than  <dig>  in our method .
apart from those common genes like tp <dig>  egfr, cdkn2a and pik3ca, the notch <dig> which functioned as tumor suppressor gene in hnsc was also indentified in our method  <cit> .
in addition, casp <dig>  which is ranked  <dig> in our method while  <dig> in drivernet, has been demonstrated that in human papillomavirus  hnsc, concurrent mutations of casp <dig> with hras can target cell cycle, death, nf-κb and other oncogenic pathways  <cit> .
furthermore, ppfia <dig> gene, which was ranked  <dig> in our method while was not detected in drivernet, acts as an invasion inhibitor in hnsc and is the highest upregulated gene in the 11q <dig> amplicon of hnsc cell lines  <cit> .
for kirc data set, our method always remarkably outperforms generank and frequency-based method .
although the performance of the top several genes in lndriver is slightly worse than drivernet and dawnrank, for latter genes, it has a remarkably better performance than drivernet method.
the curves show that the stability of our method and dawnrank is relatively good since the precision of the two methods are similar.
about top  <dig> genes,  <dig> are found in cgc in our method.
in lndriver,  <dig> genes are indentified as candidates and  <dig> of them with p − value ≤ 5% .
indeed, some well validated genes such as vhl, tp <dig>  egfr, pten and so on are ranked in the top rank in our method.
interestingly, ewsr <dig>  in cgc is not nominated as candidate drivers in drivernet and dawnrank, while it is one of the most commonly involved genes in sarcoma translocations  <cit> .
for thca, although the performances of lndriver on top several genes are same with drivernet, the overall effect is better than drivernet, frequency-based, and generank method .
in middle part of the top  <dig> genes , our method performs poor than dawnrank in this dataset, but the top  <dig> genes are all in cgc.
after the significance test, we chose  <dig> genes whose p-values were less than  <dig>  as the cancer driver genes .
with respect to several top genes, like ptpn <dig>  it encodes the protein-tyrosine phosphatase shp <dig> whose protein expression was significantly increased in human thyroid carcinoma  <cit> .
in addition, there are literatures suggesting that somatic gain-of-function mutations of ptpn <dig> are presented in breast cancer  <cit> , lung adenocarcinomas  <cit>  and etc.
braf is ranked as the second impactful driver gene which is an important event in the development of papillary thyroid cancer  <cit> .
for the ras genes , upon activation they can activate the mapk pathway  <cit>  which plays an essential role in the control of the cell cycle and differentiation  <cit> .
the analysis of identifying rare drivers
lndriver can identify not only frequently mutated driver genes, but also rare significant drivers.
the ‘rare significant drivers’ are defined as genes with p − values <  <dig>  and whose alteration frequencies are less than 2% of the patient cohort in mutation data.
in hnsc, we obtained  <dig> rare genes  in  <dig> candidate drivers with p − values <  <dig> .
four of them  are in cgc.
for example, akt <dig> , identified by lndriver, is a serine/threonine protein kinase and its downstream proteins have been reported to be frequently activated in human cancers  <cit> .
the rb <dig> gene is tumor suppressor gene identified and loss of it is considered an accelerating event in retinoblastoma  <cit> .table  <dig> the rare driver genes in hnsc
akt1
rb1
calm1
mapk1
plcg1
zbtb16
setdb1
ptk2
for kirc,  <dig> rare drivers were identified in our method and  <dig> of which are in cgc .
although some rare genes like egfr, ep <dig> and crebbp are found in drivernet, but the ranked positions are more near to the top in our method.
in addition, the activity of src , although it isn’t contained in cgc, is often associated with disease and might contribute to the development of human malignancy  <cit> .
the src family of protein tyrosine kinases provides us with many important landmarks in understanding oncogenic transformation  <cit> .
furthermore, cdkn2a  and rb <dig>  are hallmarks of lung squamous cell carcinoma  <cit>  and glioblastoma  <cit>  respectively.table  <dig> the rare driver genes in kirc
src
egfr
ep300
chd3
ewsr1
atf7ip
rb1
ncoa3
prkcd
crebbp
ddx20
smad9
kdr
pparg
atxn1
hdac1
plg
cdkn2a
met
eif6
jak2
pcna
arf6
frs2
setdb1
nos1
ppp2r1a
rab5a
svil
for thca, in addition to the frequently mutated genes , the rest of the drivers indentified by our method are rare genes .
for example, ptk2b is a member in pak signaling pathway  <cit> .table  <dig> the rare driver genes in thca
rb1
tp53
prkaca
ptk2b
pik3r1
ep300
ptpn6
casp3
jak2
ywhag
cdkn1a
pten
ctnnb1
actb
pml
atm
hsp90aa1
smad3
flnc
brca1
chd3
chek2
grin2b
nedd4
pias4
rasa1
vav1
acta1
sp1
long genes filtering analysis
in this study, we adopted gam to assign every point mutation gene with a probability weight consequently to filter frequent mutations because of long length.
with respect to ttn gene, the longest gene in human, ranked  <dig> as a driver gene of hnsc by drivernet algorithm.
however, after the step of filtering long genes in our improved method, it just ranked  <dig> and wasn’t nominated as a candidate of driver gene.
and in thca, our method didn’t identify ttn as a candidate while it was detected as the fourth ranked gene in frequency-based method.
enrichment analysis
to test biological functions of these predicted candidate drivers, kegg pathway enrichment and go functional enrichment were performed using david tool .
for hnsc, the important candidates are mainly enriched in pathways in cancer, prostate cancer, glioma, non-small cell lung cancer, melanoma, erbb signaling pathway and so on after kegg pathway enrichment .
with respect to the biological process, regulation of apoptosis, programmed cell death, cell death, nitrogen compound metabolic process, cellular biosynthetic process and etc.
are enriched after the go functional enrichment .
concerning the cellular component, identified candidates are enriched in nuclear lumen, nucleoplasm, intracellular organelle lumen, organelle lumen, membrane-enclosed lumen and cytosol etc.
.
furthermore, with regard to important molecular functions, candidate drivers are enriched in identical protein binding, nitric-oxide synthase regulator activity, structure-specific dna binding, transcription factor binding, enzyme binding and so on .
in kirc, pathways in cancer, cell cycle, melanoma and prostate cancer etc.
are enriched in kegg pathways .
in terms of biological process, positive regulation of nitrogen compound metabolic process, cellular biosynthetic process, biosynthetic process, cell cycle, transcription and gene expression etc.
are significantly enriched in go functional enrichment .
as for cellular component, candidates are enriched in nucleoplasm, nuclear lumen, nucleoplasm part, nuclear periphery, chromosome and so on .
in terms of molecular functions, transcription factor binding, protein tyrosine kinase activity, transcription regulator activity and nucleotide binding etc.
are enriched .
in thca, the pathways after kegg enrichment are prostate cancer, pathways in cancer, chronic myeloid leukemia and glioma etc.
.
in terms of biological process in go functional enrichment, candidate drivers are enriched in response to organic substance, apoptosis, programmed cell death and induction of apoptosis by intracellular signals etc.
.
with respect to cellular component, cytosol, nucleoplasm, nuclear lumen, intracellular organelle lumen and so on are enriched .
as for molecular functions, candidates are enriched in enzyme binding, enzyme binding, protein serine/threonine kinase inhibitor activity and protein kinase binding etc.
.
discussion and
CONCLUSIONS
in this work, we introduced a network-based framework by integrating transcriptome and genomics data into a gene-gene interaction network to identify significant driver gene in cancer.
by virtue of the consideration of gene length, the frequently mutated genes with long length may be filtered.
also, we constructed a network containing more genes and interaction information in order to improve the accuracy of driver genes identifying.
lndriver can identify not only frequently mutations but also rare drivers.
application on hnsc, kirc and thca datasets has demonstrated that the performance of our method is remarkably better than frequency-based, generank and drivernet method.
in addition, our method also outperforms dawnrank method in hnsc dataset.
however, in kirc and thca, dawnrank sometimes have a better performance than our method.
we will explore the causes about this phenomenon in our following work and we hope to find a new method which can have a good performance on kirc and thca.
furthermore, there are also some limitations of our method.
firstly, gene length filtering step was only applied to point mutations not including cnvs because point mutations are more inclined to be affected by gene length.
although this step has ability to filter long genes, it has randomness.
we will seek solutions to improve it and enhance robustness of it.
secondly, the information of gene-gene interaction network are more and more abundant with the development of the field.
so, we will try to integrate more information to a new gene-gene interaction network which may help us to mine more information about cancer driver genes.
moreover, it is now acknowledged that precision medicine and personalized medicine are important for patient diagnosis and treatment, so we will major in proposing new method to identify patient-specific and rare driver genes based on individual mutational and expression profiles in different tumors in the future.
additional files
additional file 1: the candidate drivers with p-values less than  <dig>  in hnsc.
additional file 2: the candidate drivers with p-values less than  <dig>  in kirc.
additional file 3: the candidate drivers with p-values less than  <dig>  in thca.
additional file 4: the results of kegg and go enrichment analysis in hnsc.
additional file 5: the results of kegg and go enrichment analysis in kirc.
additional file 6: the results of kegg and go enrichment analysis in thca.
abbreviations
ccdsconsensus coding sequences
cdnacomplementary dna
cgccancer gene census
cnvscopy-number variations
daviddatabase for annotation, visualization and integrated discovery
gamgeneralized additive model
hnschead and neck squamous cell carcinoma
hprdhuman protein reference database
icgcinternational cancer genome consortium
kirckidney renal clear cell carcinoma
lndriverlength-net-driver
ngsnext-generation sequencing
pwvprobability weight vector
snvssingle nucleotide variants
tcgathe cancer genome atlas
thcathyroid carcinoma

 cancer is driven by genetic alterations, including single nucleotide variants, small insertions or deletions, large copy-number variations and structural aberrations that accumulate during the lifetime. several international large scale cancer genomics projects, such as the cancer genome atlas, and international cancer genome consortium, etc., have produced a large volume of data in recent years and provided us with an unprecedented opportunity to better characterize the molecular signatures of human cancers. however, it is still a challenge to integrate information across the different omics data and distinguish driver mutations which can promote the cancer cell to proliferate infinitely and diffuse from passenger mutations whose changes represent neutral variation that does not influence cancer development. in response to the large volume of mutations being generated from massively parallel sequencing projects, many growing mathematical and statistical approaches to search for driver genes, driver pathways or core modules based on data integration were proposed. the most basic approach, eg. mutsig and music, is to identify driver genes based on somatic mutation rates in cancer patient populations, that is, the most commonly occurring mutations are more likely to be drivers. also, computational approaches based on evaluating the functional impact of mutations such as polyphen- and oncodriverfm were proposed. however, cancer is more closely related with a group of genes interacting together in a gene-gene interaction network. with the advent of the whole-genome measurements of somatic mutations and cnvs in the mass of cancer samples, many changes altered at network and pathway levels are found, not simply a point mutation. therefore, network- and pathway-based approaches have become one of the most promising methods to prioritize driver mutations and significantly mutated genes due to their abilities to model gene-gene interactions. varwalker is a network-assisted method to prioritize potential driver genes. another method, dawnrank prioritizes altered genes on a single patient level using pagerank algorithm. drivernet is an integrated analysis framework to identify likely driver mutations by virtue of their effect on mrna expression networks and reveals the prevalence of rare candidate driver mutations. it has been demonstrated that genes which are relatively long compared to the distribution of all human consensus coding sequences are more likely to mutate while they may be not driver genes. however, drivernet doesn’t consider the effect of gene length. also, the scale of the network in drivernet is a little small which may miss some genes and the information between genes. in this work, we develop a network-based method called length-net-driver to improve the performance of detecting driver genes based on the rationale of drivernet. our goal is to consider the point mutation genes’ length and construct a new interaction network contained more genes and interactions based on human protein reference database instead of its original gene influence graph in drivernet. furthermore, we integrate somatic snv data, cnvs data and gene expression data using gene-gene interaction network. then a greedy algorithm is applied to the integrated data to prioritize candidate genes. the application on three tcga datasets demonstrated that the performance of our method is good. methods the overview of lndriver approach in lndriver method, the population-based genomic and transcriptomic interrogations of tumor types were integrated to identify driver mutations. the pipeline is shown in fig.  fig. schematic of the lndriver. genes in somatic mutations are firstly applied to gam to filter long genes and then they will combine with cnv to construct mutation matrix. the bipartite graph is constructed based on mutation data, expression data and gene-gene influence network, where the blue nodes on the left bipartite graph represent the mutated gene and the black nodes on the right represent the outlying patient-gene events from the gene expression matrix. then greedy algorithm is applied to identify candidate driver genes. finally, enrichment analysis is employed to these candidates to explore their roles in pathways actually, some studies have indicated that genes with long length have a better chance to harbor mutations. it indicated that gene length-based filtering process is essential to perform. hence, in this study, the generalized additive model was used to assign the somatic mutation probabilities of all human genes for each sample. then a resampling test was performed to filter passenger genes whose occurring frequencies are ≥ 5% at random datasets. after the filtering procedure, cnvs are combined with it to construct a binary mutation matrix. in addition, in order to enrich the information of the gene-gene interaction network, we constructed a new interaction network using human protein reference database. as for gene expression data, we built a binary outlying matrix by nominating genes whose expression values are outside two standard deviation of the gaussian distribution as outliers. next, we formulated associations between mutation and gene expression data using a bipartite graph where the left partition of nodes represented the mutation status and the right partition of nodes represented the outlying status in each of patients. after the above process, greedy algorithm was applied on the bipartite graph to select those genes in the left partition which have the highest number of outlying expression events, and then nominated them as putative driver genes. also, the statistical significance test was assessed using a randomization framework. finally, pathway enrichment analysis was done using the database for annotation, visualization and integrated discovery online tools. to demonstrate the advantages of the approach, we analyzed three large-scale publicly available genome-transcriptome datasets in head and neck squamous cell carcinoma, thyroid carcinoma and kidney renal clear cell carcinoma. filtering long genes the length of genes in human are very different and so the mutation probabilities of different genes are in vast difference. there may be some genes which have mutations only because they are long yet they aren’t driver mutations. so, for each gene, we adopted the filtering strategies of varwalker and computed a probability weight vector by fitting a generalized additive model for each sample. denoting the vector x as the gene length of cdna, we can adopt the following model to assess the mutated probability of a gene according to its cdna length, logitπ=logπ1−π~fx where π=#mutantgenes#ccdsgenes represents the proportion of mutant genes in the researched samples, and f represents an unspecified smooth function. after the above fitting process, each gene was assigned a weight value which would be used to select genes in the next resampling procedure. then a resampling test was applied to random gene sets for each sample. the number of being selected random gene sets is same with mutant genes in specific sample. and the probability of each gene to be selected is based on the probability weight calculated in the above fitting procedure. the test was performed times in each sample following pwv. the mutation frequency was calculated for each gene using formula: fre=#selectingthegeneinresamplingtest where #  indicates the times and fre represents the frequency of the gene being selected across times in resampling process. then we filtered those genes whose frequencies were ≥ 5% that indicates the gene may occur at random unless they are cgc genes. those genes with fre < 5% which represented the gene was unlikely mutated at random were observed. greedy algorithm for detecting the candidate driver genes based on processed mutation data and expression data, they were integrated with the gene-gene interaction network into a bipartite graph. the elements on the left of bipartite graph represent the mutation status of genes in population level. and the right partition events indicate outlying expression status of the genes. an edge between g i and g j will be drawn if the gene g i in the left partition is mutated, the right gene g j is outlying expression gene and g i interacts with g j in the gene-gene interaction network. given the bipartite framework, the aim is to find the mutation genes on the left partition which cover the most events on the right of bipartite graph. to this end, the optimization method of a greedy algorithm was used to select the most covered genes: at each step, chose a mutated gene which connected to the most uncovered outlying expression genes on the right of bipartite graph. when all the connected outlying expression events were covered, the program was terminated. finally, the mutated genes ranked based on their coverage and the mostly covered mutated genes are considered as the candidate driver genes. significance test in order to assess the statistical significance of the candidate driver genes, the random framework was used by permuting n =  times of the original datasets including mutation matrix, processed outlying expression matrix and the gene-gene interaction network. then the algorithm was run on the n randomly generated datasets. finally, the real data results were assessed to see whether they are significantly different from the results on randomized datasets. the null hypothesis h is that the gene mutations have no influence on the occurrence of the cancer, and the alternative hypothesis h is that the cancer is related to the mutations of the genes. the definition of the statistical significance of gene g, whose corresponding node coverage is cov g, is the fraction times of selecting driver genes that are more than cov g in n =  random runs of the method. the calculation is listed as follows: p−valueg=∑i=1n∑j=1siδcovgij>covg∑i=1nsi where s i is the number of candidate driver genes selected in the ith run of the method. then the benjamini-hochberg method was used to correct the p-values for multiple tests and finally we chose the genes whose p-values were less than.  datasets and pre-processing we applied lndriver to thca samples, hnsc samples and kirc samples. these three datasets comprise somatic snv data, cnv data and gene expression data collected from the cancer genome atlas data portal.table description of datasets the construction of mutation matrix firstly, we collected somatic snvs in level and cnv data in level directly from tcga data portal. secondly, we removed the genes whose item of “variant_classification” is “silent” or “rna” in somatic snv data and whose length are too long according to generalized additive model and resampling test process. thirdly, the cnv information was extracted by selecting genes from amplified and deleted segments in cnv data. finally, we integrated cnv data with filtered somatic snv data by getting intersecting samples and union genes to construct a binary matrix m, whose rows indicate samples and columns indicate genes. each entry of m ij refers to the mutation status of gene j in sample i and m ij =  represents that there is labeled valid mutation in gene j of sample i. otherwise, m ij =  indicates the absence of a mutation in the jth gene of the ith sample. expression outlier matrix for gene expression dataset e, the values of it contain not available values. these values affect the results of the approach. we substituted them with the mean of all other genes in the specific samples. also, we adopted the assumption in drivernet that the expression distribution of every gene across all samples is gaussian distribution. based on the hypothesis, we converted the expression data to a binary patient-outlier matrix e' where e' =  means the expression of gene i is an outlier in patient j. the definition of the outliers is that genes whose expression values are outside the two-standard deviation range of the expression values of gene i across all the patients. gene-gene interaction network and gene annotation data cancer is a disease related with sets of genes which interact with each other in some molecular networks not only related with single gene. in order to enrich the information gene-gene interaction network in drivernet, we built an influence graph g using hprd which contains proteins to server as our reference network. the influence graph g in our work is an undirected and unweighted binary network where v represents the nodes of genes and e represents the edges among genes. when there is a correlation between gene i and gene j, g ij =  otherwise g ij =  we used the consensus coding sequences genes data which have been allocated complementary dna length based on their coding sequences from varwalker as a benchmark gene resource to select those genes that have matched ccds symbols. in order to explore the impact of the gene length, we compared genes with somatic snvs with the distribution of all human ccds gene length to filter long genes. cancer gene census genes the cgc is a database that catalogues genes whose mutations have been causally implicated in cancer, which has been widely served as benchmark in many cancer researches. in this work, we also utilized it as the standard reference list which was downloaded from cosmic and included total of genes. the analysis of the overall performance in this study, the performance of lndriver’s ability was evaluated using the number of indentifying known drivers in cgc database compared with other methods. the benchmarks of the above evaluation were precision, recall and f1score which were based on the top n genes as following: precision=#mutatedgenesincgc∩#genesfoundinlndrivers#genesfoundinlndrivers recall=#mutatedgenesincgc∩#genesfoundinlndrivers#mutatedgenesincgc f1score=2×precision×recallprecision+recall for the sake of performing the property of our method on identifying cancer related drivers, we compared the result of our method to classical frequency-based method, generank method, drivernet method and personal-based method of dawnrank. the results of the experiment on hnsc, kirc and thca datasets are shown in fig.  fig. a hnsc precision. b hnsc recall. c hnsc f1score. d kirc precision. e kirc recall. f kirc f1score. g thca precision. h thca recall. i thca f1score. the comparison of precision, recall and f1score for top ranking genes in lndriver and other methods. the x axis represents the number of top ranking genes and the y axis represents the score of the precision, recall and f1score respectively hnsc, the sixth most common cancer worldwide, was analyzed in our method. as for the overall performance of its top genes, it can be seen in fig. 2a-c that lndriver method remarkably outperforms other four methods. for the top genes, there are genes contained in cgc database of our method, while of dawnrank and of drivernet. there are genes being selected as candidates and genes of them with p-values less than in our method. apart from those common genes like tp egfr, cdkn2a and pik3ca, the notch which functioned as tumor suppressor gene in hnsc was also indentified in our method. in addition, casp which is ranked in our method while in drivernet, has been demonstrated that in human papillomavirus hnsc, concurrent mutations of casp with hras can target cell cycle, death, nf-κb and other oncogenic pathways. furthermore, ppfia gene, which was ranked in our method while was not detected in drivernet, acts as an invasion inhibitor in hnsc and is the highest upregulated gene in the 11q amplicon of hnsc cell lines. for kirc data set, our method always remarkably outperforms generank and frequency-based method. although the performance of the top several genes in lndriver is slightly worse than drivernet and dawnrank, for latter genes, it has a remarkably better performance than drivernet method. the curves show that the stability of our method and dawnrank is relatively good since the precision of the two methods are similar. about top genes, are found in cgc in our method. in lndriver, genes are indentified as candidates and of them with p − value ≤ 5%. indeed, some well validated genes such as vhl, tp egfr, pten and so on are ranked in the top rank in our method. interestingly, ewsr in cgc is not nominated as candidate drivers in drivernet and dawnrank, while it is one of the most commonly involved genes in sarcoma translocations. for thca, although the performances of lndriver on top several genes are same with drivernet, the overall effect is better than drivernet, frequency-based, and generank method. in middle part of the top genes, our method performs poor than dawnrank in this dataset, but the top genes are all in cgc. after the significance test, we chose genes whose p-values were less than as the cancer driver genes. with respect to several top genes, like ptpn it encodes the protein-tyrosine phosphatase shp whose protein expression was significantly increased in human thyroid carcinoma. in addition, there are literatures suggesting that somatic gain-of-function mutations of ptpn are presented in breast cancer, lung adenocarcinomas and etc. braf is ranked as the second impactful driver gene which is an important event in the development of papillary thyroid cancer. for the ras genes, upon activation they can activate the mapk pathway which plays an essential role in the control of the cell cycle and differentiation. the analysis of identifying rare drivers lndriver can identify not only frequently mutated driver genes, but also rare significant drivers. the ‘rare significant drivers’ are defined as genes with p − values <  and whose alteration frequencies are less than 2% of the patient cohort in mutation data. in hnsc, we obtained rare genes in candidate drivers with p − values < . four of them are in cgc. for example, akt, identified by lndriver, is a serine/threonine protein kinase and its downstream proteins have been reported to be frequently activated in human cancers. the rb gene is tumor suppressor gene identified and loss of it is considered an accelerating event in retinoblastoma.table the rare driver genes in hnsc akt1 rb1 calm1 mapk1 plcg1 zbtb16 setdb1 ptk2 for kirc, rare drivers were identified in our method and of which are in cgc. although some rare genes like egfr, ep and crebbp are found in drivernet, but the ranked positions are more near to the top in our method. in addition, the activity of src, although it isn’t contained in cgc, is often associated with disease and might contribute to the development of human malignancy. the src family of protein tyrosine kinases provides us with many important landmarks in understanding oncogenic transformation. furthermore, cdkn2a and rb are hallmarks of lung squamous cell carcinoma and glioblastoma respectively.table the rare driver genes in kirc src egfr ep300 chd3 ewsr1 atf7ip rb1 ncoa3 prkcd crebbp ddx20 smad9 kdr pparg atxn1 hdac1 plg cdkn2a met eif6 jak2 pcna arf6 frs2 setdb1 nos1 ppp2r1a rab5a svil for thca, in addition to the frequently mutated genes, the rest of the drivers indentified by our method are rare genes. for example, ptk2b is a member in pak signaling pathway.table the rare driver genes in thca rb1 tp53 prkaca ptk2b pik3r1 ep300 ptpn6 casp3 jak2 ywhag cdkn1a pten ctnnb1 actb pml atm hsp90aa1 smad3 flnc brca1 chd3 chek2 grin2b nedd4 pias4 rasa1 vav1 acta1 sp1 long genes filtering analysis in this study, we adopted gam to assign every point mutation gene with a probability weight consequently to filter frequent mutations because of long length. with respect to ttn gene, the longest gene in human, ranked as a driver gene of hnsc by drivernet algorithm. however, after the step of filtering long genes in our improved method, it just ranked and wasn’t nominated as a candidate of driver gene. and in thca, our method didn’t identify ttn as a candidate while it was detected as the fourth ranked gene in frequency-based method. enrichment analysis to test biological functions of these predicted candidate drivers, kegg pathway enrichment and go functional enrichment were performed using david tool. for hnsc, the important candidates are mainly enriched in pathways in cancer, prostate cancer, glioma, non-small cell lung cancer, melanoma, erbb signaling pathway and so on after kegg pathway enrichment. with respect to the biological process, regulation of apoptosis, programmed cell death, cell death, nitrogen compound metabolic process, cellular biosynthetic process and etc. are enriched after the go functional enrichment. concerning the cellular component, identified candidates are enriched in nuclear lumen, nucleoplasm, intracellular organelle lumen, organelle lumen, membrane-enclosed lumen and cytosol etc.. furthermore, with regard to important molecular functions, candidate drivers are enriched in identical protein binding, nitric-oxide synthase regulator activity, structure-specific dna binding, transcription factor binding, enzyme binding and so on. in kirc, pathways in cancer, cell cycle, melanoma and prostate cancer etc. are enriched in kegg pathways. in terms of biological process, positive regulation of nitrogen compound metabolic process, cellular biosynthetic process, biosynthetic process, cell cycle, transcription and gene expression etc. are significantly enriched in go functional enrichment. as for cellular component, candidates are enriched in nucleoplasm, nuclear lumen, nucleoplasm part, nuclear periphery, chromosome and so on. in terms of molecular functions, transcription factor binding, protein tyrosine kinase activity, transcription regulator activity and nucleotide binding etc. are enriched. in thca, the pathways after kegg enrichment are prostate cancer, pathways in cancer, chronic myeloid leukemia and glioma etc.. in terms of biological process in go functional enrichment, candidate drivers are enriched in response to organic substance, apoptosis, programmed cell death and induction of apoptosis by intracellular signals etc.. with respect to cellular component, cytosol, nucleoplasm, nuclear lumen, intracellular organelle lumen and so on are enriched. as for molecular functions, candidates are enriched in enzyme binding, enzyme binding, protein serine/threonine kinase inhibitor activity and protein kinase binding etc.. discussion and  in this work, we introduced a network-based framework by integrating transcriptome and genomics data into a gene-gene interaction network to identify significant driver gene in cancer. by virtue of the consideration of gene length, the frequently mutated genes with long length may be filtered. also, we constructed a network containing more genes and interaction information in order to improve the accuracy of driver genes identifying. lndriver can identify not only frequently mutations but also rare drivers. application on hnsc, kirc and thca datasets has demonstrated that the performance of our method is remarkably better than frequency-based, generank and drivernet method. in addition, our method also outperforms dawnrank method in hnsc dataset. however, in kirc and thca, dawnrank sometimes have a better performance than our method. we will explore the causes about this phenomenon in our following work and we hope to find a new method which can have a good performance on kirc and thca. furthermore, there are also some limitations of our method. firstly, gene length filtering step was only applied to point mutations not including cnvs because point mutations are more inclined to be affected by gene length. although this step has ability to filter long genes, it has randomness. we will seek solutions to improve it and enhance robustness of it. secondly, the information of gene-gene interaction network are more and more abundant with the development of the field. so, we will try to integrate more information to a new gene-gene interaction network which may help us to mine more information about cancer driver genes. moreover, it is now acknowledged that precision medicine and personalized medicine are important for patient diagnosis and treatment, so we will major in proposing new method to identify patient-specific and rare driver genes based on individual mutational and expression profiles in different tumors in the future. additional files additional file 1: the candidate drivers with p-values less than in hnsc. additional file 2: the candidate drivers with p-values less than in kirc. additional file 3: the candidate drivers with p-values less than in thca. additional file 4: the results of kegg and go enrichment analysis in hnsc. additional file 5: the results of kegg and go enrichment analysis in kirc. additional file 6: the results of kegg and go enrichment analysis in thca. abbreviations ccdsconsensus coding sequences cdnacomplementary dna cgccancer gene census cnvscopy-number variations daviddatabase for annotation, visualization and integrated discovery gamgeneralized additive model hnschead and neck squamous cell carcinoma hprdhuman protein reference database icgcinternational cancer genome consortium kirckidney renal clear cell carcinoma lndriverlength-net-driver ngsnext-generation sequencing pwvprobability weight vector snvssingle nucleotide variants tcgathe cancer genome atlas thcathyroid carcinoma 
in this paper, in order to identify driver genes effectively, we applied a generalized additive model to mutation profiles to filter genes with long length and constructed a new gene-gene interaction network. then we integrated the mutation data and expression data into the gene-gene interaction network. experiments on three tcga datasets, i.e., head and neck squamous cell carcinoma, kidney renal clear cell carcinoma and thyroid carcinoma, demonstrated that the proposed method was effective. cancer is a complex disease which is characterized by the accumulation of genetic alterations during the patient’s lifetime. cancerdriver genesmutation dataexpression datainteraction networkthe 27th international conference on genome informatics shanghai, china 3- october http://admis.fudan.edu.cn/giw2016/issue-copyright-statement© the author 2016 lastly, greedy algorithm was used to prioritize candidate driver genes from the integrated data.
 cancer is a complex disease which is characterized by the accumulation of genetic alterations during the patient’s lifetime. with the development of the next-generation sequencing technology, multiple omics data, such as cancer genomic, epigenomic and transcriptomic data etc., can be measured from each individual. correspondingly, one of the key challenges is to pinpoint functional driver mutations or pathways, which contributes to tumorigenesis, from millions of functional neutral passenger mutations.  in this paper, in order to identify driver genes effectively, we applied a generalized additive model to mutation profiles to filter genes with long length and constructed a new gene-gene interaction network. then we integrated the mutation data and expression data into the gene-gene interaction network. lastly, greedy algorithm was used to prioritize candidate driver genes from the integrated data. we named the proposed method length-net-driver.  experiments on three tcga datasets, i.e., head and neck squamous cell carcinoma, kidney renal clear cell carcinoma and thyroid carcinoma, demonstrated that the proposed method was effective. also, it can identify not only frequently mutated drivers, but also rare candidate driver genes. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords cancerdriver genesmutation dataexpression datainteraction networkthe 27th international conference on genome informatics shanghai, china 3- october http://admis.fudan.edu.cn/giw2016/issue-copyright-statement© the author 2016 