Revealing metabolite biomarkers for acupuncture treatment by linear programming based feature ...

16 Pages·2012·2.52 MB·English

Wangetal.BMCSystemsBiology2012,6(Suppl1):S15
http://www.biomedcentral.com/1752-0509/6/S1/S15
RESEARCH Open Access
Revealing metabolite biomarkers for acupuncture
treatment by linear programming based feature
selection
Yong Wang1,4*†, Qiao-Feng Wu2†, Chen Chen1, Ling-Yun Wu1,4, Xian-Zhong Yan3, Shu-Guang Yu2,
Xiang-Sun Zhang1,4*, Fan-Rong Liang2*
From The 5th IEEE International Conference on Computational Systems Biology (ISB 2011)
Zhuhai, China. 02-04 September 2011
Abstract
Background: Acupuncture has been practiced in China for thousands of years as part of the Traditional Chinese
Medicine (TCM) and has gradually accepted in western countries as an alternative or complementary treatment.
However, the underlying mechanism of acupuncture, especially whether there exists any difference between varies
acupoints, remains largely unknown, which hinders its widespread use.
Results: In this study, we develop a novel Linear Programming based Feature Selection method (LPFS) to
understand the mechanism of acupuncture effect, at molecular level, by revealing the metabolite biomarkers for
acupuncture treatment. Specifically, we generate and investigate the high-throughput metabolic profiles of
acupuncture treatment at several acupoints in human. To select the subsets of metabolites that best characterize
the acupuncture effect for each meridian point, an optimization model is proposed to identify biomarkers from
high-dimensional metabolic data from case and control samples. Importantly, we use nearest centroid as the
prototype to simultaneously minimize the number of selected features and the leave-one-out cross validation error
of classifier. We compared the performance of LPFS to several state-of-the-art methods, such as SVM recursive
feature elimination (SVM-RFE) and sparse multinomial logistic regression approach (SMLR). We find that our LPFS
method tends to reveal a small set of metabolites with small standard deviation and large shifts, which exactly
serves our requirement for good biomarker. Biologically, several metabolite biomarkers for acupuncture treatment
are revealed and serve as the candidates for further mechanism investigation. Also biomakers derived from five
meridian points, Zusanli (ST36), Liangmen (ST21), Juliao (ST3), Yanglingquan (GB34), and Weizhong (BL40), are
compared for their similarity and difference, which provide evidence for the specificity of acupoints.
Conclusions: Our result demonstrates that metabolic profiling might be a promising method to investigate the
molecular mechanism of acupuncture. Comparing with other existing methods, LPFS shows better performance to
select a small set of key molecules. In addition, LPFS is a general methodology and can be applied to other high-
dimensional data analysis, for example cancer genomics.
*Correspondence:[email protected];[email protected];[email protected]
†Contributedequally
1AcademyofMathematicsandSystemsScience,ChineseAcademyof
Sciences,Beijing100190,China
2AcupunctureandMoxibustionCollege,ChengduUniversityofTraditional
ChineseMedicine,Chengdu610075,China
Fulllistofauthorinformationisavailableattheendofthearticle
©2012Wangetal.;licenseeBioMedCentralLtd.ThisisanOpenAccessarticledistributedunderthetermsoftheCreativeCommons
AttributionLicense(http://creativecommons.org/licenses/by/2.0),whichpermitsunrestricteduse,distribution,andreproductionin
anymedium,providedtheoriginalworkisproperlycited.
Wangetal.BMCSystemsBiology2012,6(Suppl1):S15 Page2of16
http://www.biomedcentral.com/1752-0509/6/S1/S15
Background molecules largely contribute towards the classification.
Acupuncture, an important therapeutic method in Tra- Ideally we can improve the generalization performance
ditional Chinese Medicine (TCM), has been used to of classifier by identifying only the molecules that are
treat various diseases for thousand years in China. significantly contribute to the classifier. This effect is
Recently it has been gradually accepted in western coun- attributable to the overcoming of the curse of dimen-
tries as an alternative or complementary treatment. sionality. For example, if it is possible to identify a small
However, how the acupuncture works remains an open set of metabolites that is indeed capable of providing
question though acupuncture exists as one of the oldest complete discriminatory information, inexpensive diag-
continuous systems of medicine dating back 4,000 years. nostic assays for only a few metabolites might be devel-
Extensive studies have been conducted on the mechan- oped and be widely deployed in clinical settings.
ism of acupuncture to explain the effects of acupuncture Knowledge of a small set of diagnostically relevant
on various systems and symptoms [1-3]. Compared to metabolites may provide important insights into the
the relatively widespread use of acupuncture, systems mechanisms responsible for acupuncture treatment
biology is a new term to describe the recent trends in itself. Those molecules are usually termed as biomar-
biology research. It emphasizes the high-throughput kers. The procedure to reveal them is referred as feature
measurement of biological systems and focuses on the selection, biomarker identification, or feature ranking.
complex interactions in biological systems [4,5]. We Feature selection is known to be NP-hard [7] and the
highly expect that systems biology, a biology-based search becomes quickly computationally intractable.
inter-disciplinary study field, will provide tremendous Suppose we treat the feature selection task in a brute
opportunities for revealing acupuncture mechanism at force way. Given n features, we need to select m fea-
the molecular level. tures which can get the best classification accuracy (m
In this paper, we use systems biology method to study <<n) regarding to a predefined cost function. Usually in
the acupuncture treatment effect by identifying a subset classification or prediction problem, the cost function is
of important molecules from high-throughput metabolic selected as the accuracy of the prediction. The exhaus-
data. Specifically, we separate the acupuncture from tive search method goes through all the possible combi-
moxibustion and only study the effect of acupuncture nations, with the computation complexity O(nm). Thus,
on normal people by investigating the difference this method is not practical for realistic applications.
between acupuncture at particular acupoint and without Existing feature selection strategies can be roughly
acupuncture. Towards this aim, we utilize 1H nuclear categorized into three types [8]. Considering the partial
magnetic resonance (1H NMR) to investigate the effects ordering properties of the subset space, we can either
of acupuncture at several meridian points on plasma start with an empty set and successively add features, or
metabolites. Then metabolite profiles (vectors) are gen- start with the set of all features and successively filter
erated from a collection of case samples(with acupunc- them. The former type is referred to as forward selec-
ture at meridian point) and control samples (without tion while the latter is referred to as backward elimina-
acupuncture). These high-dimensional profile data is tion. The third type is the combination of the two
very similar to SNP (sequence data), gene expression approaches. However, all the above methods relies on
(transcriptome), mass spectrum (proteome), and small the greedy strategy. As an example of forward feature
molecules (metabolome) data in different levels. Then selection, we might first look for the single most discri-
the straightforward task is to identify differentially minative feature using any classifier. Then we could
expressed molecules and further classify and predict the search the single additional feature that gives the best
diagonostic category of a sample, based on its metabolite class discrimination when considered along with the
profile [6]. first feature. Keeping augmenting the feature set itera-
Generally speaking, there are two difficulties in analyz- tively in this greedy fashion, we stop until cross-valida-
ing these high-dimensional profile data. First, a large tion error estimates are minimized. As a result, the
number of features (metabolites in our case) are avail- global optimal solution usually cannot be guaranteed.
able to predict classes for a relatively small number of In this paper, we proposed a novel linear program-
samples. The presence of a significant number of irrele- ming (LP) model to address this important problem.
vant features that are unrelated to the case status makes Feature selection problem is cast into an optimization
such analysis prone to the curse of dimensionality. Sec- problem with two objectives, one is to minimize the
ond, predictive accuracy is not the only goal and further number of chosen features and the other is to maximize
biological validation and mechanism understanding call the predictive accuracy based on the centroid classifica-
for explanatory power other than black box predictive tion framework. In other words, our feature selection
results. Thus it is especially important to know which method simultaneously improves classification accuracy
Wangetal.BMCSystemsBiology2012,6(Suppl1):S15 Page3of16
http://www.biomedcentral.com/1752-0509/6/S1/S15
and selects features. Comparing with several state-of- Overview of the linear programming based feature
the-art feature selection methods, our Linear Program- selection
ming based Feature Selection (LPFS) method can select To investigate the high-dimensional data for acupunc-
a small set of features by applying strong regularization ture treatment effect, we develop a novel method, LPFS,
while keeping high accuracy. We then apply our method to select a small set of metabolites to characterize acu-
to analyze the metabolite profile data generated for acu- puncture treatment effect. The schematic illustration of
puncture treatment. We identify important molecules LPFS is shown in Figure 2. LPFS performs feature selec-
(biomarkers) related to the acupuncture treatment for tion based on the nearest centroid classifier. On one
several meridian points. Further characterization of the hand, we want to attain the best classification accuracy
biomarkers and the common and difference among sev- by minimizing the loss function. On the other hand, fea-
eral meridian points provide biological insights for acu- ture selection algorithms should be robust to noise and
puncture mechanisms at molecular level. Preliminary outliers in the data by applying strong regularization.
results in this paper were presented in our conference Here we use the parsimony principle (also known as
paper [9]. In this extended paper, we further provide the Occam’s razor) by minimizing the number of selected
detailed derivation of the method and comparison with features. Then the feature selection problem is formu-
existing works. In addition, we performed more analysis lated as a multi-objective programming. The next step is
and descriptions for the biological insights of the acu- to convert this multi-objective programming into a sin-
puncture biomarkers. gle-objective linear programming by applying the ε
method. After solving the linear programming model in
Method an efficient way, the optimal features can be selected.
Analytic workflow Centroid classification prototype
In this paper, the acupuncture treatment effect is inves- A fast and simple algorithm for classification is the cen-
tigated in the framework of systems biology. The basic troid method [6,10]. This algorithm assumes that the
analytic workflow is shown in Figure 1. As the first step, target classes correspond to individual (single) clusters
metabolite profiles are originally generated by 1H NMR and uses the cluster means (or centroids) to determine
from control samples and the samples with acupuncture the class of a new sample point (see Figure 2). A proto-
treatment in meridian points. Then we develop a linear type pattern for class Cj is defined as the arithmetic
programming based feature selection method to com- mean:
pare the two groups of metabolite profiles. Finally, a (cid:3)
1
samcuaplulnsecttuoref mtreeatatmboelnitteesffaercet.selected as biomarkers for μCj = (cid:2)(cid:2)Cj(cid:2)(cid:2)si∈Cjxi (1)
Figure1Ouranalyticworkflowtoidentifybiomarkersforacupuncturetreatmenteffect.Metaboliteprofilesareoriginallygeneratedby1H
NMRfromcontrolsamplesandthesampleswithacupuncturetreatment.Thenwedevelopalinearprogrammingbasedfeatureselection
methodtocomparethetwogroupsofmetaboliteprofiles.Finallyasmallsetofmetabolitesareselectedasbiomarkersforacupunctureeffect.
Wangetal.BMCSystemsBiology2012,6(Suppl1):S15 Page4of16
http://www.biomedcentral.com/1752-0509/6/S1/S15
Figure2TheschematicillustrationofLPFS.Underthenearestcentroidframework,LPFSrequirestominimizetheclassificationerrorandthe
numberofselectedfeatures,whichleadstoamulti-objectiveprogramming.Toensurethecomputationalefficiency,alinearprogrammingis
solvedtoidentifybiomarkers.
where si is the i-th training sample labeled as class L1(s,μ)=(cid:3)s−μ(cid:3)1 (4)
Cj. Recall that the training sample is a metabolite spec-
tra represented as a multi-dimensional vector (denoted with ║y║1 = ∑i |y(i)|, and y(i) being the i-th element of
in bold). In a similar fashion, we can obtain a prototy- vector y. The value L1(s, μ) has a linear cost in the num-
pical vector for all the other classes. During classifica- ber of features. In this study, data sets contain two
tion, the class label of an unknown sample s is classes and hence the number of calls to the distance
determined as: metric is also two. Therefore, the centroid classifier, at
run-time, is linear in the number of features. During
C(s)=argmindis(μC,s) (2) training, two prototypes are computed and the cost of
Cj j computing each prototype is O(mN), where N is the
where dis(x, y) is a distance function or: number of features and m is the number of training
samples which belong to a given class. Note that m only
C(s)=argmaxsim(μ ,s)
Cj Cj (3) vtuarreiessebleecttwioenenpdroactaessseetss. Tanhdusn,owtedcuarningvietrwainminagsoarcfoena--
where sim(x, y) is a similarity metric. This simple clas- stant and the centroid classifier has O(N) cost in the
sifier will form the basis of our LPFS method. The training phase.
advantages is that it works with any number of features. Multi-objective optimization model
And its run-time complexity is proportional to the num- Suppose we have two groups in the training dataset, the
ber of features and the complexity of the distance or case group and the control group as the gold-standard
similarity metric used. According to the experiments in data to classify new samples. We denote them set T and
[11], we select L distance metric, which is robust to F respectively. Supposing |T| = m , |F| = m , and the
1 1 2
outliers and most appropriate for the centroid classifica- computed centroids are μT and μF respectively. A simple
tion algorithm. It is defined by: classification scheme is as follows. Given a new sample
Wangetal.BMCSystemsBiology2012,6(Suppl1):S15 Page5of16
http://www.biomedcentral.com/1752-0509/6/S1/S15
s, we want to decide which group it belongs to. The L wrongly classified. Otherwise this sample should be cor-
1
discrepancy between the sample s and the groups T and rectly classified.
F can be calculated as ║s - μT║1 and ║s - μT║1. Thus a Given a case sample sl = (sl1, sl2, ...,s ln), l Î {1, 2 ...,
simple rule is l }, we should have the following constraint considering
1
the tolerable error
s∈T if (cid:3)s−μT(cid:3)1 <(cid:3)s−μF(cid:3)1 (5) (cid:2) (cid:2) (cid:2) (cid:2)
(cid:2) (cid:2) (cid:2) (cid:2)
s∈F if (cid:3)s−μT(cid:3)1 >(cid:3)s−μF(cid:3)1 (6) (cid:3)n (cid:2)(cid:2)(cid:2)(cid:2)sli−(cid:3)m1 tji/m1(cid:2)(cid:2)(cid:2)(cid:2)xi−yl <(cid:3)n (cid:2)(cid:2)(cid:2)(cid:2)sli−(cid:3)m2 fji/m2(cid:2)(cid:2)(cid:2)(cid:2)xi(11)
i=1 j=1 i=1 j=1
Let the feature number be n. We introduce the vari-
ables for feature selection as x = (x1, x2, ..., xn), where xi Similarly given a control sample sl = (sl1, sl2, ..., sln), l
= 0,1. When xi = 1, it means feature i is selected in the Î {l1 + 1, l1 + 2 ..., l1 + l2}, we should have the following
biomarker set. Otherwise it is not selected. constraint considering the tolerable error
Suppose the test dataset is U. And it is composed by (cid:2)(cid:2) (cid:2)(cid:2) (cid:2)(cid:2) (cid:2)(cid:2)
athnedc|aUseT|gr=olu1p, |UUTF|an=dl 2c.oGntirvoelngraocuapseUsFa.mUp=leUslT=∪(UslF1,, (cid:3)n (cid:2)(cid:2)(cid:2)(cid:2)sli−(cid:3)m1 tji/m1(cid:2)(cid:2)(cid:2)(cid:2)xi+yl >(cid:3)n (cid:2)(cid:2)(cid:2)(cid:2)sli−(cid:3)m2 fji/m2(cid:2)(cid:2)(cid:2)(cid:2)xi(12)
sl2, ..., sln), l Î {1, 2 ..., l1}, if it is classified correctly, we i=1 j=1 i=1 j=1
should have
Thus the objective function composes two parts, i. e.,
(cid:2)(cid:2) (cid:2)(cid:2) (cid:2)(cid:2) (cid:2)(cid:2) (cid:3)n
(cid:3)n (cid:2)(cid:2)(cid:2)(cid:2)sli−(cid:3)m1 tji/m1(cid:2)(cid:2)(cid:2)(cid:2)xi <(cid:3)n (cid:2)(cid:2)(cid:2)(cid:2)sli−(cid:3)m2 fji/m2(cid:2)(cid:2)(cid:2)(cid:2)xi (7) we want to choose as few as features minx1,x2,···,xn i=1 xi
i=1 j=1 i=1 j=1 and at the same time we want to reduce the classifica-
(cid:3)
=W(fkh1,efrke2,t.k..,=fk(nt)k1Î, tFk2,,k...=, t1kn,)2Î, ...T,,mk2.= 1, 2, ...,m 1 and fk ttihoenreerirsoar (tlroasdsef-uonffctrieolna)tioy1n,ys2m,h···iip,nyl1b+l2etwelie1=+n1l2tyhie. Icnlagsesnifeicraal-,
Similarly given a control sample sl = (sl1, sl2, ..., sln), l tion error and the number of features. Hence, the fea-
Î {l + 1, l + 2 ..., l + l }, if it is classified correctly, we
1 1 1 2 ture selection problem can be formulated as a multi-
should have
objective optimization problem with discrete variables x
(cid:2) (cid:2) (cid:2) (cid:2)
(cid:3)n (cid:2)(cid:2)(cid:2)(cid:2)sli−(cid:3)m1 tji/m1(cid:2)(cid:2)(cid:2)(cid:2)xi >(cid:3)n (cid:2)(cid:2)(cid:2)(cid:2)sli−(cid:3)m2 fji/m2(cid:2)(cid:2)(cid:2)(cid:2)xi (8) =y=((yx11,,y2,x·2·,· ,y.l.1.+,l2)xans)shoawndn incoEnqutiantiuoonu(s13)v.ariables
(cid:2) (cid:2) (cid:2) (cid:2)
(cid:4) (cid:5)
i=1 j=1 i=1 j=1 (cid:3) (cid:3)
vector-minimize(x,y) ni=1xi, li1=+1l2yi ,
With the above constraints for variable x = (x1, x2, ..., subjectto (11)(12) withxi∈{0,1},i∈{1,2,···,n}, (13)
xn), the objective function is to choose as few as fea- yi≥0,i∈{1,2,···,l1+l2}
tures, i. e.,
The first term of objective function in Equation (13) is
(cid:3)n to minimize the number of chosen features, and the sec-
x1,mx2,i··n·,xn i=1 xi (9) oMnixdedonineteisgetrolmineinairmpirzoegrtahmemtoitnagl classification error.
The optimal solutions of the two-objective optimization
Thus the feature selection problem is formulated as an
problem consist of a Pareto set, which can be solved by
integer linear programming problem in Equation (10).
transforming the two objectives of (13) into a single
(cid:3)
x1,mx2,i··n·,xn(cid:3) (cid:2)(cid:2) (cid:3) (cid:2)(cid:2) ni=1(cid:3)xi (cid:2)(cid:2) (cid:3) (cid:2)(cid:2) oalbtjeercntiavtee.sOanpeotsyiptiivcealstceaclhanriqpuareaimsethteer(cid:1)l-mteothoobdt,aiwnhtihche
s.t. ni=1(cid:2)(cid:2)pli− mj=11tji/m1(cid:2)(cid:2)xi< ni=1(cid:2)(cid:2)pli− mj=21fji/m2(cid:2)(cid:2)xi, Pareto set, with the formulation in Equation (14).
pl(cid:2)=(pl1,sl2,···,pln(cid:2))∈UT,l∈{(cid:2)1,2,···,l1}, (cid:2) (10) (cid:3) (cid:3)
(cid:3)ni=1(cid:2)(cid:2)(cid:2)pli−(cid:3)mj=11tji/m1(cid:2)(cid:2)(cid:2)xi>(cid:3)ni=1(cid:2)(cid:2)(cid:2)pli−(cid:3)mj=21fji/m2(cid:2)(cid:2)(cid:2)xi mx,iyn (cid:3) (cid:2)(cid:2) (cid:3) ni=1(cid:2)(cid:2)xi+λ li(cid:3)1=+1l2yi(cid:2)(cid:2) (cid:3) (cid:2)(cid:2)
pl=(pl1,sl2,···,pln)∈UF,l∈{1,2,···,l2}, s.t. ni=1(cid:2)(cid:2)pli− mj=11tji/m1(cid:2)(cid:2)xi−yl< ni=1(cid:2)(cid:2)pli− mj=21fji/m2(cid:2)(cid:2)xi
xi=0,1, i∈{1,2,···,n} p(cid:2)l=(pl1,sl2,···,p(cid:2)ln)∈UT,l∈{1,(cid:2)2,···,l1}, (cid:2) (14)
(cid:3) (cid:2) (cid:3) (cid:2) (cid:3) (cid:2) (cid:3) (cid:2)
When we consider the noise in the measured data, not n (cid:2)(cid:2)pli− m1tji/m1(cid:2)(cid:2)xi+yl> n (cid:2)(cid:2)pli− m2fji/m2(cid:2)(cid:2)xi
all the test samples can be classified exactly. We intro- i=1 j=1 i=1 j=1
dsaumceplteheintoUleTra∪bleUeFr.rAorndy =yi{≥y1,0,y2i,ε...,{1y,l21+,l2..}.,fol1r+evle2r}y. xi=0,p1l,=(ipl∈1,{s1l2,,2·,····,·p,lnn)},∈yjU≥F,0l,∈{1j,∈2{,1··,2·,,·l2·}·,,l1+l2}
When yi is not equal to zero, it means sample i is
Wangetal.BMCSystemsBiology2012,6(Suppl1):S15 Page6of16
http://www.biomedcentral.com/1752-0509/6/S1/S15
tiv(e14f)uniscatiomnixiend(1in4t)egiser(cid:3)prongrxaim+mλin(cid:3)g (l1I+Ll2Py).i.TThheeoobrjeetci-- ms.xt,i.yn (cid:3)ni=1(cid:2)(cid:2)(cid:2)(cid:2)pli−(cid:3)jm=11−1t(cid:3)ji/(mni=11x−i+1)λ(cid:2)(cid:2)(cid:2)(cid:2)(cid:3)xi−li1=+1yll2<yi(cid:3)ni=1(cid:2)(cid:2)(cid:2)(cid:2)pli−(cid:3)mj=21fji/m2(cid:2)(cid:2)(cid:2)(cid:2)xi
cthaellyP,awreetocasnetobbytacihnaanlglinogptiti=hm1ealpasroalmuteioti=en1rslbefolorntghiengsinto- ppll==((ppll11,,ssll22,,····(cid:3)··,,ppni=ll1nn))(cid:2)(cid:2)(cid:2)(cid:2)p∈∈liTF−,,l(cid:3)l∈∈{mj{=1111,,t22ji,,/·m···1··(cid:2)(cid:2)(cid:2)(cid:2),,xlli21+}},,yftkkl>==(((cid:3)ftkk11,,ni=ft1kk22(cid:2)(cid:2)(cid:2)(cid:2),,p··li··−··,,(cid:3)ftkknn))jm=∈∈21−FT1\\f{j{ip/p(ll}}m,,k2k∈−∈{{111),,(cid:2)(cid:2)(cid:2)(cid:2)22x,,i······,,mm21}}\\{{ll}} (16)
gle-objective optimization problem (14). Clearly, l trans- xi=0,1, i∈{1,2,···,n},yj≥0, j∈{1,2,···,l1+l2}
forms the number of chosen features into equivalent We adopt leave-one-out experiment since this particu-
classification error in (14), and controls the balance lar form of cross validation is an unbiased estimator of
between them. the generalization performance of classifier. It makes the
By solving the proposed mixed integer linear program- best use of the available data and involves no random
ming model (14), we can get solutions for the feature subsampling. Every time we pick out one sample (l = 1
1
selection variables xi, i Î {1, 2, · · ·, n}, and classification or l2 = 1) from the training data and try to classify it
error variables yj, j Î {1, 2, ...,l 1 + l 2}. Checking if xi is correctly. And by doing m1 + m2 times test we add m1
equal to 1, we can know if the corresponding feature + m constraints.
2
should be selected in the classifier. Meanwhile checking Linear programming approximation
the values of all the yj, we can estimate the classification In general, mixed integer linear programming is difficult
accuracy. For example, suppose the number of all j such to solve. To ensure the computational efficiency, mixed
that yj = 0 is N1 and the number of all j such that yj > 0 ILP in Equation (16) can be relaxed to the correspond-
is N . We can simply calculate the classification accu- ing linear programming (LP). Linear programming is the
2
racy by N /l +l and N /l +l . simplest type of mathematical programming and has
1 1 2 2 1 2
Leave-one-out cross validation been widely used in systems biology study [5,13,14].
The above model (14) is based on the general idea of Therefore, sophisticated algorithm can be adopted to
cross validation, thus it depends on the choice of T and efficiently solve this LP. In terms of computational com-
F. We noticed that there are different ways to do cross plexity, the proposed approach makes the computation
validation in feature selection [12]. One way is that the of biomarker tractable. Finally we construct the LPFS
feature selection is done with all the samples and the model in Equation (17).
cross-validation is only done for the classification proce-
devuarleu.aTtihoins ikninfdavoofrcorfostshevasltiuddatieiodnmmeathyosdevdeureelytobi“ainsftohre- ms.xt,i.yn (cid:3)ni=1(cid:2)(cid:2)(cid:2)(cid:2)pli−(cid:3)jm=11−1t(cid:3)ji/(mni=11x−i+1)λ(cid:2)(cid:2)(cid:2)(cid:2)(cid:3)xi−li1=+1yll2<yi(cid:3)ni=1(cid:2)(cid:2)(cid:2)(cid:2)pli−(cid:3)mj=21fji/m2(cid:2)(cid:2)(cid:2)(cid:2)xi
imsattoioinnclleuadke”tihnetfheeatfueraetusreelecsetiloenctipornocsetdepur.eAinnotthheercrwoasys pl=(pl1,sl2,··(cid:3)·,pni=l1n)(cid:2)(cid:2)(cid:2)(cid:2)p∈liT−,(cid:3)l∈mj{=111,t2ji,/m··1·(cid:2)(cid:2)(cid:2)(cid:2),xli1+},ytkl>=((cid:3)tk1,ni=t1k2(cid:2)(cid:2)(cid:2)(cid:2),p·li·−·,(cid:3)tkn)jm=∈21−T1\fj{i/p(l}m,2k−∈{11),(cid:2)(cid:2)(cid:2)(cid:2)2x,i···,m1}\{l} (17)
pl=(pl1,sl2,···,pln)∈F,l∈{1,2,···,l2},fk=(fk1,fk2,···,fkn)∈F\{pl},k∈{1,2,···,m2}\{l}
validation, i.e., to leave the test sample(s) out from the xi≥0, i∈{1,2,···,n},yj≥0, j∈{1,2,···,l1+l2}
training set before undergoing any feature selection.
After relaxing to continuous value, the value of the
Our feature selection model allows the freedom to optimal solution of xi (LPFS score) indicates the impor-
choose the suitable cross validation procedure according tance of feature i in the nearest centroid classifier. It
to the practical need. If more information is preferred to
should be noted that we can use other distance defini-
select biomarkers due to the scarcity of samples, we can tions instead of L in our model to achieve the non-
1
use all the samples to estimate the cross validation error linear classification effect. The parameter l can be
in the following resubstitution validation in Equation
determined by checking the output leave-one-out pre-
(15).
dictive accuracy. We also notice that LPFS model can
(cid:3) (cid:3) be extended to multi-classification task and n-fold cross
mx,iyn (cid:2) ni=1(cid:2)xi+λ li1=+1l2yi(cid:2) (cid:2) validation.
(cid:3) (cid:2) (cid:3) (cid:2) (cid:3) (cid:2) (cid:3) (cid:2)
s.t. n (cid:2)(cid:2)pli− m1tji/m1(cid:2)(cid:2)xi−yl< n (cid:2)(cid:2)pli− m2fji/m2(cid:2)(cid:2)xi
i=1 j=1 i=1 j=1 Metabonomics profiling by 1H NMR spectra
(cid:3) p(cid:2)(cid:2)l=(p(cid:3)l1,sl2,···,p(cid:2)(cid:2)ln)∈T,l∈(cid:3){1,2(cid:2)(cid:2),···,m(cid:3)1}, (cid:2)(cid:2) (15) Venous blood (3ml) was collected into a heparin sodium
n (cid:2)(cid:2)pli− m1tji/m1(cid:2)(cid:2)xi+yl> n (cid:2)(cid:2)pli− m2fji/m2(cid:2)(cid:2)xi tube and the plasma was collected by centrifugation at
i=1 j=1 i=1 j=1
pl=(pl1,sl2,···,pln)∈F,l∈{1,2,···,l2}, 1000× g at 4°C for 10 minutes. An aliquot of 300 μl
xi=0,1, i∈{1,2,···,n},yj≥0, j∈{1,2,···,l1+l2} plasma was mixed with 250 μl D2O and 50 μl TSP (3-
trimethylsilyl-2H -propionic acid) in D O (1 mg/ml) in 5
Resubstitution error rate indicates only how good are 4 2
mm NMR tube. The D O provided a field-frequency
ourbiomarkersonthetrainingdata.However,thismodel 2
lock solvent for the NMR spectrometer and the TSP
has“informationleak”andwillunderestimatetheclassifi-
served as an internal reference of chemical shift. 1H
cation error.In ourimplement,we choose themodel for
NMR spectra of the plasma samples were acquired on a
leave-one-outcrossvalidationinEquation(16).
Wangetal.BMCSystemsBiology2012,6(Suppl1):S15 Page7of16
http://www.biomedcentral.com/1752-0509/6/S1/S15
Varian INOVA 600 MHz NMR spectrometer at 27°C by measured for their concentration in plasma, and mathe-
using Carr-Purcell-Meiboom-Gill (CPMG) spin-echo matically every sample is represented by a vector in 400
pulse sequence. with a total spin-spin relaxation delay dimensional space.
(2nτ) of 320 ms. The free induction decays (FIDs) were Fifty healthy young males were randomly allocated to
collected into 32K data points with a spectral width of Zusanli(ST36),Liangmen(ST21),Juliao(ST3),Yangling-
8000 Hz and 64 scans. The FIDs were zero-filled to quan (GB34), and Weizhong (BL40) groups (The loca-
double size and multiplied by an exponential line-broad- tions of the meridian points are shown in Figure 3a.
ening factor of 0.5 Hz prior to Fourier transformation Amongthefivepoints,Zusanli,Liangmen,andJuliaoare
(FT). In addition, diffusion-edited experiments were also onthesamemeridian.).Eachgroupcontains10persons.
carried out with BPP-LED (bipolar pulse pair longitudi- Inside each group the corresponding meridian points
nal eddy current delay) pulse sequence [15,16]. The gra- were separately acupunctured for 5consecutive days. In
dient amplitude was set at 35.0 G cm-1, with a diffusion addition,twentyhealthyyoungmalesarerecruitedasthe
delay of 100 ms. A total of 128 transients and 16k data blankcontrolgroups.Allthetwentypeoplearemeasured
points were collected with a spectral width of 8000 Hz. beforethestartof5consecutivedaysandadditionallyten
A line-broadening factor of 1 Hz was applied to FIDs of them are measured after 5 consecutive days. Fasting
before Fourier transformation. venousbloodwastakeninallthesubjects.Plasmameta-
All plasma 1H NMR spectra were manually phased bolites were measured by 1H NMR to derive metabolic
and baseline corrected using VNMR 6.1C software (Var- profiles (see details in Method Section). Furthermore to
ian, Inc.). For CPMG spectra, each spectrum over the exclude possible noises, all the seventy malesare strictly
range of 0.4 to 4.4 was data-reduced into integrated trained to make sure their metabolic profiles are mea-
regions of equal width (0.01 ppm). For BPP-LED data, suredinverysimilarconditions.Thedetailedexperimen-
each spectrum over the range of 0.1 to 6.0 was segmen- talmethodcanbefoundin[20].Insummary,wehave80
ted into regions of equal width (0.01 ppm). The regions samples grouped into Zusanli (10 samples, acupuncture
containing the resonance from residual water (4.6-5.1) point ST36), Liangmen (10 samples, acupuncture point
were excluded. The integral values of each spectrum ST21),Juliao (10samples,acupuncturepoint ST3),Yan-
was normalized to constant sum of all integrals in a glinquan(10 samples, acupuncture point GB34), Weiz-
spectrum in order to reduce any significant concentra- hong (10 samples, acupuncture point BL40), Control I
tion differences between samples [17,18]. Identification (10samples,normalpeoplemeasured after theconsecu-
of metabolites in spectra was accomplished based on lit- tive 5 days), and Control II (20 samples, normal people
eratures and the Chenomx NMR Suite 5.0 (Chenomx, measuredbeforetheconsecutive5days).
Calgary, Canada).
Classification experiments design
Results With the data, we design experiments to identify bio-
Metabonomics data generation markers for the acupuncture treatment of each meridian
To investigate the acupuncture treatment effects, we ori- point. The overall design of biomarker identification
ginally obtained metabonomics data of plasma metabo- experiments is shown in Figure 3b. We categorize eighty
lites in healthy males at five meridian points using samples into 8 groups shown as the circles in Figure 3b.
Proton NMR. Proton NMR (also named as Hydrogen-1 ST36, ST21, ST3, GB34, and BL40 each has 10 samples.
NMR, or 1H NMR) applies nuclear magnetic resonance The 20 samples in Control II are naturally decomposed
in NMR spectroscopy with respect to hydrogen-1 nuclei into two groups with equal size, Pre1 (10 samples with
within the molecules of a substance, in order to deter- follow-up measurement after 5 days) and Pre2 (10 sam-
mine the structure of the molecules [19]. ples without follow-up measurement). Treating the Pre1
As a result, most organic compounds are character- as the common control set, we have seven classification
ized by chemical shift values, which are usually tasks (Exp1 to Exp7) shown as the lines in Figure 3b.
expressed in parts per million (ppm) by frequency and For example, task Exp1 tries to identify a subset of
are in the range +14 to -4 ppm. Chemical shift values metabolites to classify Pre1 as the control and ST36 as
are not precise, but typically they are regarded mainly as the case. In this way, Exp1 to Exp5 aim to identify the
orientational. The exact value of chemical shift depends biomarkers for acupuncture treatment on ST36, ST21,
on molecular structure and the solvent in which the ST3, GB34, and BL40 respectively. While Exp6 tries to
spectrum is being recorded. These chemical shift values capture the metabolite change by 5 consecutive days
can be mapped to eight metabolic subsets (amino acids, without acupuncture. And Exp7 tries to test if there are
carbohydrates, energy, glycans, lipids, nucleotides, sec- significant metabolite change for the people under simi-
ondary metabolites/xenobiotics, vitamins, and cofactors). lar condition. Exp6 and Exp7 serve as the control stu-
In our experiment, 400 chemical shift values are dies to guarantee the significance of our result.
Wangetal.BMCSystemsBiology2012,6(Suppl1):S15 Page8of16
http://www.biomedcentral.com/1752-0509/6/S1/S15
Figure3Overalldesignofthebiomarkeridentificationexperiments.a)Metaboliteprofilesareoriginallygeneratedby1HNMRfromfive
acupuncturepoints(Zusanli(ST36),Liangmen(ST21),Juliao(ST3),Yanglingquan(GB34),andWeizhong(BL40)).b)Themetaboliteprofilesare
groupedinto7setsandthebiomarkeridentificationproblemisdesignedas7binaryclassificationexperiments.
Global characterization of the data 80 samples should be clustered by with or without acu-
We firstperformhierarchical clusteringonthe 80 meta- puncture treatment and then by their meridian points.
bolic profiles. The results are shown in Figure 4a. If the However,allthesamplelabelsaremixedintheclustering
samplescanbeclearlydiscriminatedbyglobalpattern,the result(Figure4a)andwecannotseeclearlyboundaries.
Wangetal.BMCSystemsBiology2012,6(Suppl1):S15 Page9of16
http://www.biomedcentral.com/1752-0509/6/S1/S15
Figure 4 Global characterization of the metabolite profiles. a) Hierarchical clustering of the metabolic profiles of the 80 samples. b)
Centroidsforthesevendatasets.Thehorizontalunitsareexpressionvaluesforthemetabolites.Themetabolitesaresortedbytheirchemical
shiftvalues.
Furthermore, we calculate the centroids for the seven respectively. We define the ordinary two-sample t-statis-
groups of samples in Figure 4b by averaging the 10 sam- tic [21] as
ples for their metabolite expression values. These cen-
x¯ −y¯
troids are plotted side by side in Figure 4b, which shows T i i (18)
i
that these centroids are very similar and it’s very diffi- si
cult to detect the global difference. Where x¯i,y¯i, and si are the mean of case, mean of
The above results together demonstrate that global
control, and the standard deviation of the samples for
pattern in metabolic profiles cannot discriminate the metabolite i. From t-statistic Ti we can easily calculate a
Zusanli, Yanglingquan, Liangmen, Juliao, Weizhong,
p-value. Usually a feature is selected if its corresponding
Pre1, Pre2, and Control I groups. Thus it is necessary to
p-value is smaller than a predefined threshold 0.05.
find the local pattern in the profile data. Our strategy is
The standard definition of the fold-change [21] for
to find a subset of metabolites as biomarkers to achieve metabolite i is
clear discrimination.
(cid:5)x
Comparison with other approaches FCi = (cid:5)yi (19)
Before we conduct the acupuncture biomarker identifi- i
cations, we benchmark our LPFS method by comparing Where x´ij and ýij are the raw expression values of
with several existing state-of-the-art methods. There are metabolite i in sample j in the case and control, respec-
many existing feature selection methods and they can
tively. In our implementation, we computes the differ-
be roughly categorized into three types, filter, wrapper,
ence of means (i.e. the fold-change for log-transformed
and embedded methods. To make the comparison sim-
data) and then rank the metabolites by their absolute
ple and comprehensive, we pick out some representative
values. We choose a cutoff 2 to select the significant
methods in each type to compare in the same dataset.
ones.
Filter methods select features as a preprocessing step
On the other hand, wrapper method ranks features
and feature selection part is independent of a machine
based on their effects on classification accuracy. It takes
learning algorithm (classifier). This is computationally
dependencies of the feature subset on the learning algo-
efficient. Fold change and t-test are the simplest and
rithm into account and is computationally more
popular methods to identify biomarker. They are usually
demanding. Support Vector Machine-Recursive Feature
the representative methods for filter methods.
Elimination (SVM-RFE) is one of the most successful
Let xij and yij denote the log expression values of wrapper method based algorithm in the feature selection
metabolite i in sample j in the case and control,
[22]. SVM-RFE conducts feature selection in a
Wangetal.BMCSystemsBiology2012,6(Suppl1):S15 Page10of16
http://www.biomedcentral.com/1752-0509/6/S1/S15
sequential backward elimination manner, which starts derivation versus the difference of mean expression
with all the features and discards one feature at a time value. This is a scatter-plot used in [27], which plots sig-
by checking the SVM accuracy. It has been widely used nificance versus fold-change on the y- and x-axes,
and extended in high-dimensional data analysis [23]. We respectively. Conceptually, it is very similar to the vol-
compare our LPFS method with SVM-RFE in metabolite cano plot in statistics [28]. Plotting points in this way
data. results in two regions of interest in the plot: those
Since our LPFS method is an embedded method and points that are found towards the bottom of the plot
simultaneously optimize classification accuracy and the and far to either the left- or the right-hand side. These
number of selected features, we specifically choose to represent values that display large magnitude fold
compare with an existing method with similar strategy, changes as well as high statistical significance with small
called sparse multinomial logistic regression approach standard derivation.
(SMLR). It was developed to jointly and simultaneously The t-test based method identifies 172 metabolites if
identify the optimal nonlinear classifier, and select the we choose a cutoff 1.73 (corresponding p-value 0.05). A
optimal set of features via the optimization of a single strict threshold will still select 84 metabolites with cutoff
posterior objective function (see [24] and [25]). SMLR 2.84 (corresponding p-value 0.005). We list the top 10 in
has been extensively applied in problems in systems Table 1 and plot them in Figure 5a. We find that the
biology [26]. SMLR is freely available at http://www.cs. ordinary t-statistic selects metabolites with low standard
duke.edu/~amink/software/smlr/ and we take the default deviations.
values for the parameters in our calculation. The fold change based method identifies 159 metabo-
Withoutlossofgenerality,wetakeExp1inFigure3bas lites if we choose a commonly used cutoff 2 [28]. And
an example to identify a subset of metabolites to discri- there are still 97 metabolites by choosing a cutoff 4.
minate ST36 and Pre1. We select biomarkers from the Again the top 10 are listed in Table 2 and plotted in
metabolicprofilesandcomparetheresultsintwoways. Figure 5b. It’s clear that the fold-changes select metabo-
Firstly, we compare different methods in Figure 5, by lites with large shifts between control and treatment.
assessing the quality of the selected biomarkers. We While SMLR selects 37 features to achieve the 100%
simply plot all the metabolites by their standard leave-one-out predictive accuracy. These 37 metabolites
Figure5Methodscomparisonviavolcanolikeplots.ComparisonofourLPFSmethodwithexistingmethodsregardingtotheidentified
biomarkers.Allthe400metabolitesareplottedintoatwodimensionalplane.Theselectedbiomarkersarehighlightedinred.Thex-axisdenotes
thedifferenceofmeansandthey-axisdenotesthestandardderivation.Goodbiomakersshouldlocateeitherintheleftbottomcornerorinthe
rightbottomcorner.a)volcanolikeplotoft-testmethod.Thetop10featuresareinred.b)volcanolikeplotoffoldchangemethod.Thetop10
featuresareinred.c)volcanolikeplotofSMLRmethod.d)volcanolikeplotofSVM-RFEmethod.e)volcanolikeplotofourLPFSmethod.