BACKGROUND
micrornas  play important regulatory roles in many fundamental biological processes for disease development and progression.
especially, tremendous researches have demonstrated that mirnas can serve as oncogene or tumor suppressor in various cancer types  <cit> .
during the last decade, benefitting from the development of mirna microarray and small rna-seq techniques, mirna expression data has been widely used in the comparison of diseased samples with control samples, or different subtypes of diseased samples.
the mirnas with most discriminant capacity, regarded as biomarkers, have assisted in diagnosis, prognosis prediction and therapeutic assessment of cancers  <cit> , and sometimes they are even more accurate than coding-gene markers  <cit> .
in order to search biomarkers, the analysis of differential gene expression is performed and genes are ranked according to certain criteria.
the evaluation on the quality of biomarkers is mainly based on statistical or machine learning approaches, whose corresponding measurements are statistical significance and classification accuracy, respectively.
till now, a variety of statistical methods have been applied into the gene expression analysis.
fold change has been used as an initial metric for measuring the significance of change in expression levels between two groups of samples  <cit> , and t-test  <cit>  is a widely-used statistical method to select differentially expressed genes.
besides, researchers have developed many alternatives of t-test, such as anova  <cit> , wilcoxon test  <cit> , sam  <cit> , rvm  <cit> , limma  <cit> , varmixt  <cit>  and smvar  <cit> .
most of the present criteria are for univariate analysis.
as the rapid development of high-throughput experimental technology, gene expression experiments have become more and more specialized and diversified.
especially, tissue-specific and condition-specific researches have largely been emerged.
the single-gene biomarkers are often unreliable or have insufficient ability to distinguish subtypes or different conditions for complex diseases.
in order to increase the sensitivity and specificity of biomarkers, in many studies, the top ranked genes according to some selection metric were put together and used as composite biomarkers.
this method is not guaranteed to find optimal biomarkers, since there may be redundant information among the genes because of correlation.
and, many genes individually do not show good discriminative ability between different groups, while perform well together with other genes.
therefore, multivariate analysis is necessary to examine the performance of multiple genes as a whole.
a common method for multivariate statistical analysis is hotelling’s t-square test  <cit> .
note that in gene expression analysis, the number of samples is often very limited.
as the dimensionality increases, the statistical inference often fails to provide reliable results.
feature selection is a major branch of methods for screening biomarkers  <cit> .
from a machine learning point of view, biomarkers correspond to the features with most discerning power.
a multivariate feature selection method scores feature subsets and rank them, usually by their classification accuracy.
for example, in order to select gene combinations, cui et al.
<cit>  and xu et al.
<cit>  used support vector machines to separate cancer and normal tissues, and assessed classification accuracy for all the k-gene combinations, for k≤ <dig> and k≤ <dig>  respectively.
these multivariate analysis methods can avoid feature redundancy but may run into exponential complexity due to the huge search space.
another issue is about the interpretation of computational results.
too complex classifier  and too many variables/features in the composite biomarkers could be useless, because the results are extremely difficult for biological explanation and validation.
mirna expression analysis usually follow the same procedure and approaches as mrna expression analysis, such as hypothesis test  <cit> , clustering  <cit>  and classification  <cit>  based on machine learning models.
meanwhile, the above mentioned problems also exist in mirna data.
besides, due to the low intensity on expression level and small difference between mirna sequences, the selection of mirna biomarkers becomes more challenging.
in this study, instead of screening single mirnas or large mirna sets, we aim to find the combinatorial biomarkers, i.e., k-mirna combinations, where k is a small number.
to avoid exponential number of combinations, we propose a clustering-based method to reduce the number of candidate combinations and conduct a highly efficient search.
the basic idea is to assess only the combinations consisting of representative members from clusters that are generated based on expression level similarity, rather than all combinations.
in order to further reduce the search space, a proper criterion is needed to rank the mirnas in the clusters, and only the most promising ones can be selected as the representatives of their clusters to form the candidate biomarkers.
clustering approaches have been extensively used to find co-expressed genes.
genes in the same clusters are usually functionally related.
there have been some studies that adopted clustering-based methods for feature selection.
for example, jaeger et al.
proposed to use a fuzzy c-means clustering method to pre-filter genes before ranking genes individually  <cit> .
that is, only one representative gene is selected from each cluster and involved in the ranking procedure.
a similar approach was proposed by hanczar et al.
<cit> , who used k-means clustering to select ‘prototype genes’.
in both of these two methods, the number of clusters needs to be pre-defined.
actually, it is an important issue to find the proper number of clusters.
in order to address this issue, wang et al.
developed a novel hybrid approach  <cit> .
they applied hierarchical clustering on these genes to generate a dendrogram, then the optimal number of clusters was determined by a leave-one-out cross-validation  strategy by trying all of the different clusterings by breaking up the dendrogram.
in all of these methods, there is no defined criterion on how to determine the number of clusters or the proper size of clusters, though wang et al.
conducted an empirical analysis of loocv  <cit> .
moreover, these methods typically used genes which are the closest to centers of their clusters as the representative genes, while whether the center gene is the best choice is questionable.
in another similar research proposed by sahu et al.
<cit> , k-means clustering was adopted, while signal-to-noise ratio  was used to rank genes in every clusters.
our approach has two major differences from the previous approaches: i) there is a new criterion to select the most discriminant member in each cluster, ii) there is a defined criterion to determine whether a cluster should be split.
and the goal of this study is slightly different from the aforementioned literatures in that we aim to develop efficient method for identifying mirna combinatorial biomarkers, instead of large feature subsets which are hard to be interpreted in biology.
we have conducted a series of experiments to compare different criteria for selecting representative genes from clusters and splitting raw clusters.
we also compared the new method with some widely used feature selections methods.
the experimental results demonstrate that our proposed method is very effective in screening genes in the clusters.
the resulting clusters can greatly reduce the number of combinations to be assessed, and obtain high-quality combinations in the mean time.
the selected mirna combinations have not only high discriminative ability, but also enriched pathways closely related with tumorigenesis.
moreover, many frequently present mirnas in these combinations have been validated to be associated with breast cancer development in previous literatures.
methods
the proposed method consists of three major steps.
the first step is a pre-screening to remove uninformative mirnas using welch’s t-test.
the second step is a hierarchical clustering on the remaining genes.
in the last step, representative mirnas are selected from every clusters to form mirna combinations as candidate biomarkers.
both the criteria for assessing the qualities of clustering and selecting representative mirnas within clusters are defined through a linear discriminant method.
the flowchart of the method is shown in fig.
<dig>
fig.
<dig> flowchart of the feature extraction method
fisher linear discriminant analysis
fisher linear discriminant analysis   <cit>  seeks a best linear combinations of features to achieve maximum separation on the projected feature space, by optimizing the object function which is a ratio of inter-class difference to intra-class difference of data.
since fda projects original features onto one-dimensional features, it is used not only for classification but also for dimensionality reduction.
different from principal component analysis , fda works in a supervised manner, thus the projected features are more discriminative with respect to the classification task.
the algorithm of fda is briefly described in the below.
for a binary classification problem, suppose x is the training set which has n samples with p dimensions, i.e., x={x
<dig> x
<dig> x
<dig> ⋯,x
n}, where x
is  are p-dimensional sample vectors belonging to class c
<dig> or c
<dig>  let m
<dig> and m
<dig> be the mean vectors of samples in these two classes, respectively, and w be the optimal projection direction.
according to fda’s object function, w is obtained by eq.
,
<dig> w∝sw−1m0−m <dig>
where s
w is the sum of variance within each class, i.e.,
<dig> sw=∑xi∈c0xi−m0xi−m0t+∑xi∈c1xi−m1xi−m1t
given this optimal direction, all x
is are projected onto w to get the new one-dimensional sample sets y={y
<dig> y
<dig> ⋯,y
n}, where
<dig> yi=wtxifori= <dig> ,⋯,n
as for classification, the definition of threshold  has multiple choices.
normally, the threshold, y
<dig>  can be computed as eq.
,
<dig> y0=n0m0+n1m1n0+n <dig>
where m
<dig> and m
<dig> are means for the two classes in the projected data space, i.e.,
<dig> m0=wtm <dig>
<dig> m1=wtm <dig>
in the test phase, a test sample x is firstly projected onto w, then the resulted value y is compared against y
<dig>  if y is larger than or equal to y
<dig>  it will be assigned to class c
<dig>  otherwise, it is regarded as belonging to class c
<dig>
the criteria for selecting representative members
in order to avoid feature redundancy, a representative member is selected from each cluster of mirnas.
although many methods directly choose the mean or center member, it is necessary to define some criterion to rank the members by their contribution or potential in the separation of groups of samples.
as described in “fisher linear discriminant analysis” section, fda aims to find the projection direction, w, with maximum discriminative capacity.
and, w can be regarded as a vector of weights, indicating the importance of features.
intuitively, those features with the largest weights are the most informative for classification.
in other words, the magnitude of each component of w implies the relevance of the corresponding mirna to classification.
let i be the index set of all mirnas, i.e.
i={ <dig> ,…,p}, and i
c be the index set of the mirnas in the cluster c. the index of the representative member of c, i
c, satisfies eq.
:
<dig> |w|=maxj∈ic|),
where w is a component of w.
the criteria for splitting clusters
besides selection of mirnas, the determination of number or size of clusters also has a great impact on the performance of search algorithms.
too many clusters are more likely to find high quality combinatorial biomarkers but can result into huge computational complexity.
the extreme case is the trivial clustering that each single mirna is a cluster.
on the contrary, too few clusters would miss valuable combinations since only a few representative mirnas are considered.
thus here is a tradeoff between accuracy and efficiency.
instead of explicitly specify the number of clusters, we seek proper criteria for determining whether a given cluster needs to be split into smaller clusters.
here, we define the criterion mainly based on the loss of information caused by projection.
intuitively, if the cluster members are diversified, it would be very hard to find a unified direction for projection, so the data samples would suffer great information loss after projection, which indicates that the cluster needs to be split.
thus, we define a measure called mean squared loss , to estimate average information loss in a cluster.
equation  formulates this measure.
let h be the hyperplane that passes the mean point of the data samples and has normal direction of w , then msl is defined as:
<dig> msl=∑i=1n|ph|2n,
where p
h denotes the projection of a vector onto h, m is the mean vector of samples, x
i is the ith data sample.
since h is perpendicular to w, we regard the projection of the difference between x
i and m on h as an approximative loss caused by fda projection.
furthermore, considering that the samples may differ in data magnitudes, we define another criterion called mean loss rate  as shown in eq.
,
<dig> mlr=∑i=1n|ph||xi|n,
where mlr denotes the averaged loss rate, i.e.
the ratio of the loss  to the norm of sample.
the whole pipeline is described in algorithm ?
?, in which the mlr is used as the selection criterion.
evaluation criteria
the performance of different criteria are evaluated using two measures for the resulted combinations which are ranked top  <dig>   <dig> and  <dig>  respectively.
one is average rank, denoted by avgrank, and the other is the proportion of the true top combinations identified by the method, denoted by hitratio.
these two measures are defined in eqs.
and , respectively.
for top n
k-mirna combinations searched by the method,
<dig> avgrankn=∑1≤i≤nrankin,
where r
a
n
k
i is the true rank of the ith best combination identified among all k-mirna combinations .
all these ranks are determined according to classification accuracy.
<dig> hitration=hitnn,
where h
i
t
n is the number of hits in the n best combinations searched by the method.
a hit means the searched result is truly among the top-n combinations.
apparently, small avgrank and high hitratio of the search results indicate good performance of the algorithm for identifying high-quality biomarker candidates.
in addition, to evaluate the classification performance of the selected mirna combinations, we used three accuracy measures, namely sensitivity, specificity and total accuracy .
RESULTS
data sets
in this study, we used two public mirna data sets from ncbi geo  <cit> , gse <dig>  <cit>  and gse <dig>  <cit> , which were measured by illumina human v <dig> mirna panel and agilent- <dig> human mirna microarray platform, respectively.
both of these two studies aim to explore function of micrornas in breast tumorigenesis, and reveal potential therapeutic targets.
there are a total of  <dig> samples collected from  <dig> breast cancer patients, including  <dig> pairs of matched tumor and adjacent peri-tumoral breast tissues, and  <dig> unmatched tissues in gse  <dig>  and in gse <dig>  there are  <dig> samples from  <dig> breast cancer patients, including  <dig> estrogen receptor -negative tissues, and  <dig> er-positive tissues.
the detailed statistics of patient characteristics are shown in table  <dig>
in order to ensure the data quality, we removed the mirnas whose expression levels were not detected or below the threshold value in more than 30% of the samples.
gse <dig> was classified into two categories according to tumor and peri-tumor status, while gse <dig> was divided into two categories according to er status.
finally, the gse <dig> data set contains  <dig> pairs of tumor and peri-tumor profiles and the gse <dig> data set contains  <dig> samples of er-positive and  <dig> of er-negative.
experimental settings
as a pre-screening step, welch’s t-test was conducted on the two data sets.
mirnas with pvalue greater than  <dig>  were filtered out, and the remaining mirnas were clustered by a hierarchical clustering with average-link method.
next, the hierarchical tree was cut into raw clusters.
in order to find natural cluster divisions in the hierarchical tree, we computed inconsistency coefficient for each link in the tree  <cit> .
this value compares the height of a link in a hierarchical tree with the average height of links below it.
inconsistent links indicate the border of naturally divided clusters.
the inconsistency coefficients range from  <dig> to  <dig>  for both the two data sets.
thus we specified an inconsistency coefficient threshold of  <dig> to partition the two hierarchical trees into raw clusters, resulting in  <dig> and  <dig> clusters, respectively.
further, fda-based criteria were used to determine whether or not those clusters should be split into smaller clusters.
after the final clusters were determined, a representative member was selected from each cluster to form mirna combinations.
the comparison on several criteria for selecting representative members within clusters and for splitting clusters are given in the following two sections.
in the final step, each combination was assessed by classification accuracy.
we evaluated the classification accuracies of all k-combinations  comprised by the selected representative mirnas using libsvm  <cit>  with default parameters via 5-fold cross validation.
the avgrank and hitratio were calculated based on the true ranking lists obtained by exhaustive searches.
comparison of criteria for selecting representatives in clusters
in previous researches, the center gene, i.e.
the gene closest to the cluster center, was usually selected as the representative member of a cluster .
also, some researchers proposed specific ranking criteria, such as the signal-to-noise ratio  proposed by sahu et al.
<cit> .
here, we compared the fda measure with three other methods based on center-gene, snr, and pvalue of t-test, respectively.
all the k-mirna combinations  resulted from these selection criteria were assessed.
in order to evaluate the quality of the search results, we examined top  <dig>   <dig> and  <dig> best combinations identified by these four methods and recorded their avgranks and hitratios obtained on gse <dig> and gse <dig> in tables  <dig> and  <dig>  respectively.
acenter denotes the method using center gene as the representative member
the results show that a proper selection criterion is crucial for searching high-quality mirna combinations.
specifically, fda and t-test based criteria have significant advantage over other two methods, and snr is slightly better than center-gene.
for instance, on gse <dig>  fda and t-test successfully identified the best pair and triple mirnas, and the second-best quadruple, whose accuracy is only  <dig> % lower than the best one.
fda and t-test have much smaller avgranks than snr and center-gene, no matter what the k is and how long the top list considered.
moreover, fda hits 80% of the top  <dig> pairs.
both fda and t-test catch majority of the top-ranked mirna pairs and triples.
as k increases to  <dig>  the hit ratio decreases greatly, which is mainly due to the exponentially expanded search space.
avgrank and hitratio values of the top  <dig> lists obtained by the four methods on gse <dig> data set are depicted in figs.
<dig> and  <dig>  respectively.
fig.
<dig>
avgrank of top  <dig> lists obtained by the four methods for gse22220
fig.
<dig>
hitratio of top  <dig> lists obtained by the four methods for gse22220
generally, these methods have consistent performance on the two data sets.
for gse <dig>  the accuracies of combinatorial mirnas are very high.
even a pair of mirnas can achieve the accuracy as high as  <dig> %, and the highest accuracy of quadruples is  <dig> %, which suggests that the k-mirna combinations  are sufficient for separating the samples from two classes.
the goal of gse <dig> is to discriminate tumor and peri-tumor samples, thus the differential expressed signal may be widespread.
if too many combinations can achieve high accuracies, the real biomarkers may become not that notable.
thus, the results of average rank and hit ratio seem to be worse than those of gse <dig>
comparison of criteria for splitting clusters
in this study, we propose two criteria for determining whether a given cluster should be split, i.e., mean squared loss  and mean loss rate .
considering that different clusters contain different numbers of mirnas, instead of using the original msl, we divide the squared loss by m , and use msl′=mslm in the analysis.
the m
s
l
′s for all raw clusters in gse <dig> sorted in ascending order are shown in fig.
<dig>  it can be observed that a dramatic change occurs a little above  <dig>  on the curve.
thus, we set the threshold as  <dig> , where the steepest ascent locates.
and we found that in gse <dig> the value is very close.
fig.
<dig>
msl curve of the initial clusters of gse40525
obviously, mlr would grow rapidly as the number of mirnas in the clusters increases.
here we set the threshold as 1−1m <dig>  which is a relatively loose criterion.
mlr works as a supplement to msl.
in our experiment, if either of these two criteria is not satisfied , the cluster should be partitioned.
we compared the refined clustering  by using these two criteria and the conventional hierarchical clustering  without further splitting.
the results are shown in table  <dig>  considering that the results obtained after refinement are generally better than those from raw clusters because more clusters make larger search space of mirna combinations, we did not use results of the raw clusters in the rc experiment.
instead, we tried different inconsistency coefficients for hc, which produced close number of clusters as rc did, and selected the best results, while in rc the inconsistency coefficient and thresholds of msl and mlr are fixed as mentioned above.
arc: refined clustering, in which the inconsistency coefficient for raw clusters and thresholds of msl and mlr are fixed
bhc: hierarchical clustering, which performs the best by trying different inconsistency coefficients
generally, rc has a comparable or better performance to the best hc.
for the top  <dig> list, hc shows some advantage, while for top  <dig> and  <dig>  rc performs better.
these results suggest that msl and mlr can help to improve the clustering quality, and save effort on searching good threshold to yield clusters in the hierarchical tree.
basically, both msl and mlr measure the part of information that cannot be expressed by the projected features, i.e.
information loss during the projection.
different from the absolute loss represented by msl, mlr measures the relative loss and plays a part in screening low-quality clusters when the variances of mirnas differ greatly.
comparison with existing feature selection methods
we further compared the new method with some widely used feature selection methods, including the correlation-based feature selection   <cit> , best-first search , consistency-based selection  <cit> , chi-square score  <cit> , information gain   <cit> , random forest  filter  <cit> , t-test  <cit>  and the wilcoxon rank-sum test  <cit> .
among these methods, cfs, bfs and consistency-based methods determine the number of selected features automatically.
for other methods, we chose the subsets consisting of top  <dig>   <dig> and  <dig> features in the assessment.
the r package fselector  <cit>  was adopted to implement these eight methods in the comparison experiments.
most feature selection methods shown in table  <dig> are filtering methods, except that bfs is a wrapper method, and the proposed method can be regarded as a hybrid method, which conducts filtering within clusters and then acts as a wrapper method using svms.
the best mirna combinations identified by the new method achieve the highest accuracies on both data sets, increasing the total accuracies by about  <dig> % on gse <dig> and  <dig> % on gse <dig> compared with the best accuracies obtained by other methods.
this result again demonstrates the validity of clustering-based screening and the fda criteria.
given the representative members selected from clusters, the search space is greatly reduced and the best combinations can be efficiently searched.
hence, the new method achieves a good balance between efficiency and accuracy.
χ
2
clufdab
<dig>
<dig>
athe numbers before and after ‘/’ denotes feature numbers of gse <dig> and gse <dig>  respectively sensitivity = tp/, specificity = tn/ ta: total accuracy
bclufda denotes the clustering-based feature selection using fda method for selecting representative mirnas
functional analysis on the selected mirnas
in order to perform functional enrichment on the mirna combinatorial biomarkers, we analyzed the enriched pathways of their target genes by using mirpath  <cit> .
for gse <dig> data set, triples of mirnas have the best discriminant capacity, and the top  <dig> significant pathways for the best triple are: fatty acid biosynthesis, pi3k-akt signaling pathway, prostate cancer, tgf-beta signaling pathway and p <dig> signaling pathway, all of which have pvalues below 5×10− <dig>  for the gse <dig> data set, the enriched pathways include pi3k-akt signaling pathway, nf-kappa b signaling pathway, focal adhesion, etc.
interestingly, pi3k-akt signaling pathway is significantly enriched in both data sets.
this pathway acts as regulator of cell proliferation, differentiation, apoptosis, and plays important roles in tumorigenesis.
in addition, we found that the top-ranked combinations usually have overlapped members.
for example, all the top  <dig> pairs and triples of gse <dig> contains hsa-mir-139-5p.
and, best quadruples often contain best pairs and triples.
therefore, we recorded the most frequent mirnas in pairs and triples respectively and got their intersection set .
there are  <dig> mirnas and  <dig> mirnas for gse <dig> and gse <dig>  respectively.
furthermore, these mirnas were searched against two mirna-disease relationship databases, namely hmdd v <dig>   <cit>  and mir2disease  <cit> .
among the  <dig> most frequent mirnas,  <dig> mirnas were reported in previous literatures as being involved in the development of breast cancer .
it is worth noting that  <dig> of the mirnas are not covered in the top  <dig> list evaluated by statistical significance of the conventional t-test ranking method, but all of them have supporting literatures.
specifically, hsa-mir- <dig> ranks  <dig>  hsa-mir- <dig> ranks  <dig>  hsa-mir- <dig> ranks  <dig>  and hsa-mir- <dig> ranks  <dig>  both mir- <dig> and mir- <dig> have been demonstrated as inhibitors of tumorigenesis with biological-experimental evidence.
hsa-mir-18a*
hsa-mir-146b-5p
hsa-mir-149
hsa-mir-224
hsa-mir-340
hsa-mir-452*
hsa-mir-100
hsa-mir-141
hsa-mir-365
amirnas that have evidence of association with breast cancer  are in bold
discussions
in this paper, we propose to identify mirna combinatorial biomarkers due to the important role that mirnas play in the development of cancer and also some good properties of combinatorial biomarkers.
the reasons for searching biomarkers of mirna combinations are manifold.
firstly, single-gene biomarkers identified by uni-variate analysis are often unreliable with low specificity for discriminating complex disease properties.
thus, multi-gene biomarkers are in great need.
however, the biomarkers containing too many genes, resulted from feature subset selection, are extremely difficult to be interpreted in biomedicine.
for instance, if we have identified a k-tuple combinatorial biomarker, and we want to validate the overexpress/unexpress rule as well as inter-correlation in this biomarker, the over/under express pattern has a total of  <dig> k possibilities.
moreover, correlation coefficient can only be computed between two genes, and now there have been some studies on the conditionally independent properties in triples .
but there have been no effective means to measure or validate the interconnection among multiple genes yet.
moreover, according to our results, combinations with small k have sufficient capability to separate groups of samples.
we have also examined the accuracy of using all representative members selected from every clusters , which are much lower than the best k-mirna combinations , decreasing by about 4% on gse <dig> and 9% on gse <dig>  this result further demonstrates the usefulness of small combinatorial biomarkers.
alla
<dig>
<dig>
<dig>
<dig>
aall: the full set of representative mirnas selected from clusters sensitivity = tp/, specificity = tn/ ta: total accuracy
CONCLUSIONS
mirna expression files have been widely used in the identification of biomarkers for complex diseases.
due to the low specificity of single-gene biomarker and difficulty in interpretating large gene set, this study aims to develop efficient algorithm for searching mirna combinatorial biomarkers with high discriminability.
we propose a clustering-based method to avoid brute force search, and define two types of criteria for refining clusters and selecting representative members.
the former criterion aims to measure the loss during the feature projection by fisher linear discriminant analysis, and determine whether or not to partition the given clusters.
the latter criterion aims to select the most informative mirnas in the clusters according to their contribution for classification in fda model.
we conducted experiments on two breast cancer mirna expression profiles.
the fda-based selection method performs the best with regard to average rank of the top searched results and hit ratio on the true top list.
the fda-based cluster splitting rule has also been demonstrated to be effective in refining the clustering results.
for the two data sets, k-mirna combinations with k≤ <dig> have sufficient capacity to discriminate the samples .
this method can also be applied to the search of combinations with larger k, and mrna expression data.
the top-ranked mirna combinations are worth further study on their functions as well as interactions of the mirnas.
as an additional computational analysis, the most frequent mirnas occurring in top  <dig> pairs and triples have been searched again mirna-disease database.
among the  <dig> most frequent mirnas,  <dig> mirnas have supporting literatures of their roles in the development of breast cancer.
from the fifteenth asia pacific bioinformatics conference shenzhen, china.
16- <dig> january 2017

 micrornas play important regulatory roles in many fundamental biological processes for disease development and progression. especially, tremendous researches have demonstrated that mirnas can serve as oncogene or tumor suppressor in various cancer types. during the last decade, benefitting from the development of mirna microarray and small rna-seq techniques, mirna expression data has been widely used in the comparison of diseased samples with control samples, or different subtypes of diseased samples. the mirnas with most discriminant capacity, regarded as biomarkers, have assisted in diagnosis, prognosis prediction and therapeutic assessment of cancers, and sometimes they are even more accurate than coding-gene markers. in order to search biomarkers, the analysis of differential gene expression is performed and genes are ranked according to certain criteria. the evaluation on the quality of biomarkers is mainly based on statistical or machine learning approaches, whose corresponding measurements are statistical significance and classification accuracy, respectively. till now, a variety of statistical methods have been applied into the gene expression analysis. fold change has been used as an initial metric for measuring the significance of change in expression levels between two groups of samples, and t-test is a widely-used statistical method to select differentially expressed genes. besides, researchers have developed many alternatives of t-test, such as anova, wilcoxon test, sam, rvm, limma, varmixt and smvar. most of the present criteria are for univariate analysis. as the rapid development of high-throughput experimental technology, gene expression experiments have become more and more specialized and diversified. especially, tissue-specific and condition-specific researches have largely been emerged. the single-gene biomarkers are often unreliable or have insufficient ability to distinguish subtypes or different conditions for complex diseases. in order to increase the sensitivity and specificity of biomarkers, in many studies, the top ranked genes according to some selection metric were put together and used as composite biomarkers. this method is not guaranteed to find optimal biomarkers, since there may be redundant information among the genes because of correlation. and, many genes individually do not show good discriminative ability between different groups, while perform well together with other genes. therefore, multivariate analysis is necessary to examine the performance of multiple genes as a whole. a common method for multivariate statistical analysis is hotelling’s t-square test. note that in gene expression analysis, the number of samples is often very limited. as the dimensionality increases, the statistical inference often fails to provide reliable results. feature selection is a major branch of methods for screening biomarkers. from a machine learning point of view, biomarkers correspond to the features with most discerning power. a multivariate feature selection method scores feature subsets and rank them, usually by their classification accuracy. for example, in order to select gene combinations, cui et al. and xu et al. used support vector machines to separate cancer and normal tissues, and assessed classification accuracy for all the k-gene combinations, for k≤ and k≤ respectively. these multivariate analysis methods can avoid feature redundancy but may run into exponential complexity due to the huge search space. another issue is about the interpretation of computational results. too complex classifier and too many variables/features in the composite biomarkers could be useless, because the results are extremely difficult for biological explanation and validation. mirna expression analysis usually follow the same procedure and approaches as mrna expression analysis, such as hypothesis test, clustering and classification based on machine learning models. meanwhile, the above mentioned problems also exist in mirna data. besides, due to the low intensity on expression level and small difference between mirna sequences, the selection of mirna biomarkers becomes more challenging. in this study, instead of screening single mirnas or large mirna sets, we aim to find the combinatorial biomarkers, i.e., k-mirna combinations, where k is a small number. to avoid exponential number of combinations, we propose a clustering-based method to reduce the number of candidate combinations and conduct a highly efficient search. the basic idea is to assess only the combinations consisting of representative members from clusters that are generated based on expression level similarity, rather than all combinations. in order to further reduce the search space, a proper criterion is needed to rank the mirnas in the clusters, and only the most promising ones can be selected as the representatives of their clusters to form the candidate biomarkers. clustering approaches have been extensively used to find co-expressed genes. genes in the same clusters are usually functionally related. there have been some studies that adopted clustering-based methods for feature selection. for example, jaeger et al. proposed to use a fuzzy c-means clustering method to pre-filter genes before ranking genes individually. that is, only one representative gene is selected from each cluster and involved in the ranking procedure. a similar approach was proposed by hanczar et al., who used k-means clustering to select ‘prototype genes’. in both of these two methods, the number of clusters needs to be pre-defined. actually, it is an important issue to find the proper number of clusters. in order to address this issue, wang et al. developed a novel hybrid approach. they applied hierarchical clustering on these genes to generate a dendrogram, then the optimal number of clusters was determined by a leave-one-out cross-validation strategy by trying all of the different clusterings by breaking up the dendrogram. in all of these methods, there is no defined criterion on how to determine the number of clusters or the proper size of clusters, though wang et al. conducted an empirical analysis of loocv. moreover, these methods typically used genes which are the closest to centers of their clusters as the representative genes, while whether the center gene is the best choice is questionable. in another similar research proposed by sahu et al., k-means clustering was adopted, while signal-to-noise ratio was used to rank genes in every clusters. our approach has two major differences from the previous approaches: i) there is a new criterion to select the most discriminant member in each cluster, ii) there is a defined criterion to determine whether a cluster should be split. and the goal of this study is slightly different from the aforementioned literatures in that we aim to develop efficient method for identifying mirna combinatorial biomarkers, instead of large feature subsets which are hard to be interpreted in biology. we have conducted a series of experiments to compare different criteria for selecting representative genes from clusters and splitting raw clusters. we also compared the new method with some widely used feature selections methods. the experimental results demonstrate that our proposed method is very effective in screening genes in the clusters. the resulting clusters can greatly reduce the number of combinations to be assessed, and obtain high-quality combinations in the mean time. the selected mirna combinations have not only high discriminative ability, but also enriched pathways closely related with tumorigenesis. moreover, many frequently present mirnas in these combinations have been validated to be associated with breast cancer development in previous literatures. methods the proposed method consists of three major steps. the first step is a pre-screening to remove uninformative mirnas using welch’s t-test. the second step is a hierarchical clustering on the remaining genes. in the last step, representative mirnas are selected from every clusters to form mirna combinations as candidate biomarkers. both the criteria for assessing the qualities of clustering and selecting representative mirnas within clusters are defined through a linear discriminant method. the flowchart of the method is shown in fig. fig. flowchart of the feature extraction method fisher linear discriminant analysis fisher linear discriminant analysis seeks a best linear combinations of features to achieve maximum separation on the projected feature space, by optimizing the object function which is a ratio of inter-class difference to intra-class difference of data. since fda projects original features onto one-dimensional features, it is used not only for classification but also for dimensionality reduction. different from principal component analysis, fda works in a supervised manner, thus the projected features are more discriminative with respect to the classification task. the algorithm of fda is briefly described in the below. for a binary classification problem, suppose x is the training set which has n samples with p dimensions, i.e., x={x x x ⋯,x n}, where x is are p-dimensional sample vectors belonging to class c or c let m and m be the mean vectors of samples in these two classes, respectively, and w be the optimal projection direction. according to fda’s object function, w is obtained by eq., w∝sw−1m0−m where s w is the sum of variance within each class, i.e., sw=∑xi∈c0xi−m0xi−m0t+∑xi∈c1xi−m1xi−m1t given this optimal direction, all x is are projected onto w to get the new one-dimensional sample sets y={y y ⋯,y n}, where yi=wtxifori=,⋯,n as for classification, the definition of threshold has multiple choices. normally, the threshold, y can be computed as eq., y0=n0m0+n1m1n0+n where m and m are means for the two classes in the projected data space, i.e., m0=wtm m1=wtm in the test phase, a test sample x is firstly projected onto w, then the resulted value y is compared against y if y is larger than or equal to y it will be assigned to class c otherwise, it is regarded as belonging to class c the criteria for selecting representative members in order to avoid feature redundancy, a representative member is selected from each cluster of mirnas. although many methods directly choose the mean or center member, it is necessary to define some criterion to rank the members by their contribution or potential in the separation of groups of samples. as described in “fisher linear discriminant analysis” section, fda aims to find the projection direction, w, with maximum discriminative capacity. and, w can be regarded as a vector of weights, indicating the importance of features. intuitively, those features with the largest weights are the most informative for classification. in other words, the magnitude of each component of w implies the relevance of the corresponding mirna to classification. let i be the index set of all mirnas, i.e. i={,…,p}, and i c be the index set of the mirnas in the cluster c. the index of the representative member of c, i c, satisfies eq.: |w|=maxj∈ic|), where w is a component of w. the criteria for splitting clusters besides selection of mirnas, the determination of number or size of clusters also has a great impact on the performance of search algorithms. too many clusters are more likely to find high quality combinatorial biomarkers but can result into huge computational complexity. the extreme case is the trivial clustering that each single mirna is a cluster. on the contrary, too few clusters would miss valuable combinations since only a few representative mirnas are considered. thus here is a tradeoff between accuracy and efficiency. instead of explicitly specify the number of clusters, we seek proper criteria for determining whether a given cluster needs to be split into smaller clusters. here, we define the criterion mainly based on the loss of information caused by projection. intuitively, if the cluster members are diversified, it would be very hard to find a unified direction for projection, so the data samples would suffer great information loss after projection, which indicates that the cluster needs to be split. thus, we define a measure called mean squared loss, to estimate average information loss in a cluster. equation formulates this measure. let h be the hyperplane that passes the mean point of the data samples and has normal direction of w, then msl is defined as: msl=∑i=1n|ph|2n, where p h denotes the projection of a vector onto h, m is the mean vector of samples, x i is the ith data sample. since h is perpendicular to w, we regard the projection of the difference between x i and m on h as an approximative loss caused by fda projection. furthermore, considering that the samples may differ in data magnitudes, we define another criterion called mean loss rate as shown in eq., mlr=∑i=1n|ph||xi|n, where mlr denotes the averaged loss rate, i.e. the ratio of the loss to the norm of sample. the whole pipeline is described in algorithm??, in which the mlr is used as the selection criterion. evaluation criteria the performance of different criteria are evaluated using two measures for the resulted combinations which are ranked top and respectively. one is average rank, denoted by avgrank, and the other is the proportion of the true top combinations identified by the method, denoted by hitratio. these two measures are defined in eqs. and, respectively. for top n k-mirna combinations searched by the method, avgrankn=∑1≤i≤nrankin, where r a n k i is the true rank of the ith best combination identified among all k-mirna combinations. all these ranks are determined according to classification accuracy. hitration=hitnn, where h i t n is the number of hits in the n best combinations searched by the method. a hit means the searched result is truly among the top-n combinations. apparently, small avgrank and high hitratio of the search results indicate good performance of the algorithm for identifying high-quality biomarker candidates. in addition, to evaluate the classification performance of the selected mirna combinations, we used three accuracy measures, namely sensitivity, specificity and total accuracy.  data sets in this study, we used two public mirna data sets from ncbi geo, gse and gse, which were measured by illumina human v mirna panel and agilent- human mirna microarray platform, respectively. both of these two studies aim to explore function of micrornas in breast tumorigenesis, and reveal potential therapeutic targets. there are a total of samples collected from breast cancer patients, including pairs of matched tumor and adjacent peri-tumoral breast tissues, and unmatched tissues in gse and in gse there are samples from breast cancer patients, including estrogen receptor -negative tissues, and er-positive tissues. the detailed statistics of patient characteristics are shown in table in order to ensure the data quality, we removed the mirnas whose expression levels were not detected or below the threshold value in more than 30% of the samples. gse was classified into two categories according to tumor and peri-tumor status, while gse was divided into two categories according to er status. finally, the gse data set contains pairs of tumor and peri-tumor profiles and the gse data set contains samples of er-positive and of er-negative. experimental settings as a pre-screening step, welch’s t-test was conducted on the two data sets. mirnas with pvalue greater than were filtered out, and the remaining mirnas were clustered by a hierarchical clustering with average-link method. next, the hierarchical tree was cut into raw clusters. in order to find natural cluster divisions in the hierarchical tree, we computed inconsistency coefficient for each link in the tree. this value compares the height of a link in a hierarchical tree with the average height of links below it. inconsistent links indicate the border of naturally divided clusters. the inconsistency coefficients range from to for both the two data sets. thus we specified an inconsistency coefficient threshold of to partition the two hierarchical trees into raw clusters, resulting in and clusters, respectively. further, fda-based criteria were used to determine whether or not those clusters should be split into smaller clusters. after the final clusters were determined, a representative member was selected from each cluster to form mirna combinations. the comparison on several criteria for selecting representative members within clusters and for splitting clusters are given in the following two sections. in the final step, each combination was assessed by classification accuracy. we evaluated the classification accuracies of all k-combinations comprised by the selected representative mirnas using libsvm with default parameters via 5-fold cross validation. the avgrank and hitratio were calculated based on the true ranking lists obtained by exhaustive searches. comparison of criteria for selecting representatives in clusters in previous researches, the center gene, i.e. the gene closest to the cluster center, was usually selected as the representative member of a cluster. also, some researchers proposed specific ranking criteria, such as the signal-to-noise ratio proposed by sahu et al.. here, we compared the fda measure with three other methods based on center-gene, snr, and pvalue of t-test, respectively. all the k-mirna combinations resulted from these selection criteria were assessed. in order to evaluate the quality of the search results, we examined top and best combinations identified by these four methods and recorded their avgranks and hitratios obtained on gse and gse in tables and respectively. acenter denotes the method using center gene as the representative member the results show that a proper selection criterion is crucial for searching high-quality mirna combinations. specifically, fda and t-test based criteria have significant advantage over other two methods, and snr is slightly better than center-gene. for instance, on gse fda and t-test successfully identified the best pair and triple mirnas, and the second-best quadruple, whose accuracy is only % lower than the best one. fda and t-test have much smaller avgranks than snr and center-gene, no matter what the k is and how long the top list considered. moreover, fda hits 80% of the top pairs. both fda and t-test catch majority of the top-ranked mirna pairs and triples. as k increases to the hit ratio decreases greatly, which is mainly due to the exponentially expanded search space. avgrank and hitratio values of the top lists obtained by the four methods on gse data set are depicted in figs. and respectively. fig. avgrank of top lists obtained by the four methods for gse22220 fig. hitratio of top lists obtained by the four methods for gse22220 generally, these methods have consistent performance on the two data sets. for gse the accuracies of combinatorial mirnas are very high. even a pair of mirnas can achieve the accuracy as high as %, and the highest accuracy of quadruples is %, which suggests that the k-mirna combinations are sufficient for separating the samples from two classes. the goal of gse is to discriminate tumor and peri-tumor samples, thus the differential expressed signal may be widespread. if too many combinations can achieve high accuracies, the real biomarkers may become not that notable. thus, the results of average rank and hit ratio seem to be worse than those of gse comparison of criteria for splitting clusters in this study, we propose two criteria for determining whether a given cluster should be split, i.e., mean squared loss and mean loss rate. considering that different clusters contain different numbers of mirnas, instead of using the original msl, we divide the squared loss by m, and use msl′=mslm in the analysis. the m s l ′s for all raw clusters in gse sorted in ascending order are shown in fig. it can be observed that a dramatic change occurs a little above on the curve. thus, we set the threshold as, where the steepest ascent locates. and we found that in gse the value is very close. fig. msl curve of the initial clusters of gse40525 obviously, mlr would grow rapidly as the number of mirnas in the clusters increases. here we set the threshold as 1−1m which is a relatively loose criterion. mlr works as a supplement to msl. in our experiment, if either of these two criteria is not satisfied, the cluster should be partitioned. we compared the refined clustering by using these two criteria and the conventional hierarchical clustering without further splitting. the results are shown in table considering that the results obtained after refinement are generally better than those from raw clusters because more clusters make larger search space of mirna combinations, we did not use results of the raw clusters in the rc experiment. instead, we tried different inconsistency coefficients for hc, which produced close number of clusters as rc did, and selected the best results, while in rc the inconsistency coefficient and thresholds of msl and mlr are fixed as mentioned above. arc: refined clustering, in which the inconsistency coefficient for raw clusters and thresholds of msl and mlr are fixed bhc: hierarchical clustering, which performs the best by trying different inconsistency coefficients generally, rc has a comparable or better performance to the best hc. for the top list, hc shows some advantage, while for top and rc performs better. these results suggest that msl and mlr can help to improve the clustering quality, and save effort on searching good threshold to yield clusters in the hierarchical tree. basically, both msl and mlr measure the part of information that cannot be expressed by the projected features, i.e. information loss during the projection. different from the absolute loss represented by msl, mlr measures the relative loss and plays a part in screening low-quality clusters when the variances of mirnas differ greatly. comparison with existing feature selection methods we further compared the new method with some widely used feature selection methods, including the correlation-based feature selection, best-first search, consistency-based selection, chi-square score, information gain, random forest filter, t-test and the wilcoxon rank-sum test. among these methods, cfs, bfs and consistency-based methods determine the number of selected features automatically. for other methods, we chose the subsets consisting of top and features in the assessment. the r package fselector was adopted to implement these eight methods in the comparison experiments. most feature selection methods shown in table are filtering methods, except that bfs is a wrapper method, and the proposed method can be regarded as a hybrid method, which conducts filtering within clusters and then acts as a wrapper method using svms. the best mirna combinations identified by the new method achieve the highest accuracies on both data sets, increasing the total accuracies by about % on gse and % on gse compared with the best accuracies obtained by other methods. this result again demonstrates the validity of clustering-based screening and the fda criteria. given the representative members selected from clusters, the search space is greatly reduced and the best combinations can be efficiently searched. hence, the new method achieves a good balance between efficiency and accuracy. χ 2 clufdab athe numbers before and after ‘/’ denotes feature numbers of gse and gse respectively sensitivity = tp/, specificity = tn/ ta: total accuracy bclufda denotes the clustering-based feature selection using fda method for selecting representative mirnas functional analysis on the selected mirnas in order to perform functional enrichment on the mirna combinatorial biomarkers, we analyzed the enriched pathways of their target genes by using mirpath. for gse data set, triples of mirnas have the best discriminant capacity, and the top significant pathways for the best triple are: fatty acid biosynthesis, pi3k-akt signaling pathway, prostate cancer, tgf-beta signaling pathway and p signaling pathway, all of which have pvalues below 5×10− for the gse data set, the enriched pathways include pi3k-akt signaling pathway, nf-kappa b signaling pathway, focal adhesion, etc. interestingly, pi3k-akt signaling pathway is significantly enriched in both data sets. this pathway acts as regulator of cell proliferation, differentiation, apoptosis, and plays important roles in tumorigenesis. in addition, we found that the top-ranked combinations usually have overlapped members. for example, all the top pairs and triples of gse contains hsa-mir-139-5p. and, best quadruples often contain best pairs and triples. therefore, we recorded the most frequent mirnas in pairs and triples respectively and got their intersection set. there are mirnas and mirnas for gse and gse respectively. furthermore, these mirnas were searched against two mirna-disease relationship databases, namely hmdd v and mir2disease. among the most frequent mirnas, mirnas were reported in previous literatures as being involved in the development of breast cancer. it is worth noting that of the mirnas are not covered in the top list evaluated by statistical significance of the conventional t-test ranking method, but all of them have supporting literatures. specifically, hsa-mir- ranks hsa-mir- ranks hsa-mir- ranks and hsa-mir- ranks both mir- and mir- have been demonstrated as inhibitors of tumorigenesis with biological-experimental evidence. hsa-mir-18a* hsa-mir-146b-5p hsa-mir-149 hsa-mir-224 hsa-mir-340 hsa-mir-452* hsa-mir-100 hsa-mir-141 hsa-mir-365 amirnas that have evidence of association with breast cancer are in bold discussions in this paper, we propose to identify mirna combinatorial biomarkers due to the important role that mirnas play in the development of cancer and also some good properties of combinatorial biomarkers. the reasons for searching biomarkers of mirna combinations are manifold. firstly, single-gene biomarkers identified by uni-variate analysis are often unreliable with low specificity for discriminating complex disease properties. thus, multi-gene biomarkers are in great need. however, the biomarkers containing too many genes, resulted from feature subset selection, are extremely difficult to be interpreted in biomedicine. for instance, if we have identified a k-tuple combinatorial biomarker, and we want to validate the overexpress/unexpress rule as well as inter-correlation in this biomarker, the over/under express pattern has a total of k possibilities. moreover, correlation coefficient can only be computed between two genes, and now there have been some studies on the conditionally independent properties in triples. but there have been no effective means to measure or validate the interconnection among multiple genes yet. moreover, according to our results, combinations with small k have sufficient capability to separate groups of samples. we have also examined the accuracy of using all representative members selected from every clusters, which are much lower than the best k-mirna combinations, decreasing by about 4% on gse and 9% on gse this result further demonstrates the usefulness of small combinatorial biomarkers. alla aall: the full set of representative mirnas selected from clusters sensitivity = tp/, specificity = tn/ ta: total accuracy  mirna expression files have been widely used in the identification of biomarkers for complex diseases. due to the low specificity of single-gene biomarker and difficulty in interpretating large gene set, this study aims to develop efficient algorithm for searching mirna combinatorial biomarkers with high discriminability. we propose a clustering-based method to avoid brute force search, and define two types of criteria for refining clusters and selecting representative members. the former criterion aims to measure the loss during the feature projection by fisher linear discriminant analysis, and determine whether or not to partition the given clusters. the latter criterion aims to select the most informative mirnas in the clusters according to their contribution for classification in fda model. we conducted experiments on two breast cancer mirna expression profiles. the fda-based selection method performs the best with regard to average rank of the top searched results and hit ratio on the true top list. the fda-based cluster splitting rule has also been demonstrated to be effective in refining the clustering results. for the two data sets, k-mirna combinations with k≤ have sufficient capacity to discriminate the samples. this method can also be applied to the search of combinations with larger k, and mrna expression data. the top-ranked mirna combinations are worth further study on their functions as well as interactions of the mirnas. as an additional computational analysis, the most frequent mirnas occurring in top pairs and triples have been searched again mirna-disease database. among the most frequent mirnas, mirnas have supporting literatures of their roles in the development of breast cancer. from the fifteenth asia pacific bioinformatics conference shenzhen, china. 16- january 2017 
in order to avoid exhaustive search and redundant information, mirnas are firstly clustered, then the combinations of representative cluster members are assessed as potential biomarkers. both the criteria for the partition of clusters and selection of representative members are based on fisher linear discriminant analysis. this study aims to select combinatorial mirna biomarkers, which have higher sensitivity and specificity than single-gene biomarkers. our experimental results demonstrate that the clustering-based method can identify microrna combinatorial biomarkers with high accuracy and efficiency. as the rapid development of high-throughput experimental technology, gene expression experiments have become more and more specialized and diversified. in the comparison with eight common feature selection methods, this clustering-based method performs the best with regard to the discriminative ability of selected biomarkers. the fda-based criterion has been demonstrated to be superior to three other criteria in selecting representative members, and also good at refining clusters.
 micrornas have great potential serving as tumor biomarkers and therapeutic targets. as the rapid development of high-throughput experimental technology, gene expression experiments have become more and more specialized and diversified. the complex data structure has brought great challenge for the identification of biomarkers. in the meantime, current statistical and machine learning methods for detecting biomarkers have the problem of low reliability and biased criteria.  this study aims to select combinatorial mirna biomarkers, which have higher sensitivity and specificity than single-gene biomarkers. in order to avoid exhaustive search and redundant information, mirnas are firstly clustered, then the combinations of representative cluster members are assessed as potential biomarkers. both the criteria for the partition of clusters and selection of representative members are based on fisher linear discriminant analysis. the fda-based criterion has been demonstrated to be superior to three other criteria in selecting representative members, and also good at refining clusters. in the comparison with eight common feature selection methods, this clustering-based method performs the best with regard to the discriminative ability of selected biomarkers.  our experimental results demonstrate that the clustering-based method can identify microrna combinatorial biomarkers with high accuracy and efficiency. our method and data are available to the public upon request. keywords micrornabiomarkerclusteringthe fifteenth asia pacific bioinformatics conference apbc 2017shenzhen, china 16- january issue-copyright-statement© the author 2017 