Bioinformatics, 34, 2018, i323–i332
doi: 10.1093/bioinformatics/bty252
ISMB 2018

Dai Hai Nguyen1,*, Canh Hao Nguyen1 and Hiroshi Mamitsuka1,2
1

Bioinformatics Center, Institute for Chemical Research, Kyoto University, Uji 611-0011, Japan and 2Department ofComputer Science, Alato University, Espoo 02150, Finland*To whom correspondence should be addressed.

Abstract
Motivation: Recent success in metabolite identification from tandem mass spectra has been led bymachine learning, which has two stages: mapping mass spectra to molecular fingerprint vectors
and then retrieving candidate molecules from the database. In the first stage, i.e. fingerprint prediction, spectrum peaks are features and considering their interactions would be reasonable for moreaccurate identification of unknown metabolites. Existing approaches of fingerprint prediction are
based on only individual peaks in the spectra, without explicitly considering the peak interactions.
Also the current cutting-edge method is based on kernels, which are computationally heavy anddifficult to interpret.
Results: We propose two learning models that allow to incorporate peak interactions forfingerprint prediction. First, we extend the state-of-the-art kernel learning method by developing
kernels for peak interactions to combine with kernels for peaks through multiple kernel learning(MKL). Second, we formulate a sparse interaction model for metabolite peaks, which wecall SIMPLE, which is computationally light and interpretable for fingerprint prediction. The formulation of SIMPLE is convex and guarantees global optimization, for which we develop
an alternating direction method of multipliers (ADMM) algorithm. Experiments using the
MassBank dataset show that both models achieved comparative prediction accuracy withthe current top-performance kernel method. Furthermore SIMPLE clearly revealed individual
peaks and peak interactions which contribute to enhancing the performance of fingerprint
prediction.
Availability and implementation: The code will be accessed through http://mamitsukalab.org/tools/
SIMPLE/.
Contact: hai@kuicr.kyoto-u.ac.jp

1 Introduction
Metabolites are small molecules which are used in, or created by,chemical reactions occurring in living organisms (Wishart, 2007).
They play a lot of important functions such as energy transport,
signaling, building block of cells and inhibition/catalysis.
Understanding biochemical characteristics of metabolites is an
essential and important part of metabolomics to enlarge theknowledge of biological systems. It is also key to development ofmany applications and areas such as biotechnology, biomedicine or
pharmaceutical sciences. Mass Spectrometry is a common techniquein analytical chemistry (de Hoffmann and Stroobant, 2007) for metabolite identification. A mass spectrometer analyses a chemical

C The Author(s) 2018. Published by Oxford University Press.
Vsample by fragmenting it and measuring the mass-to-charge ratios(m/z) of its fragments. The resulting mass spectrum (MS) or tandem
mass spectrum (MS/MS) is represented by a graph with m/z on thex-axis and the relative abundance of ions with m/z values on the
y-axis (Fig. 1). Another common representation of mass spectrum is
a list of peaks, each defined by its m/z and intensity value (top-rightcorner of Fig. 1).
Many methods have been proposed for metabolite identification.
A traditional approach is to compare a query MS or MS/MS spectrum of unknown compound against a database of a number of reference MS or MS/MS spectra (Scheubert et al., 2013). The candidate
molecules from the database are ranked based on the similarity, andi323

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/),
which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact
journals.permissions@oup.com

Downloaded from https://academic.oup.com/bioinformatics/article-abstract/34/13/i323/5045791 by Kyoto Daigaku Bungakubu Toshokan user on 17 June 2019SIMPLE: Sparse Interaction Model over Peaks ofmoLEcules for fast, interpretable metabolite
identification from tandem mass spectra

i324

Fig. 1. Example MS spectrum from Human Metabolome Database (Wishart
et al., 2012) for 1-Methylhistidine (HMBD00001), with its corresponding chemical structure (top-left) and peak list (top-right)

Our particular focus is the first step, which predicts the fingerprintthat is a binary (or rarely real-valued) vector, indicating the presence(or absence) of certain substructure. This step has been tackled bymany ML methods, including linear discriminative analysis (LDA)(Imre et al., 2008) and decision tree (Hummel et al., 2010). A notable method is FingerID (Heinonen et al., 2012), which used support
vector machine (SVM) with kernels for pairs of mass spectra, including integral mass kernel and probability product kernel (PPK)(Jebara et al., 2004).
We point out three drawbacks of the existing ML approaches: i)all methods are based on the information from individual peaks inspectra, without explicitly considering peak interactions, whichleads to the limitation of predictive performance. ii) There is aninteresting (not ML-based) software, which outputs, given a spectrum, a so-called fragmentation tree (FT) (Böcker and Rasche, 2008;
Rasche et al., 2011). In a FT, possible fragments corresponding topeaks in spectrum are shown as node labels, where parent–childrelationships in this tree are inclusive relations of fragments in chemical structure. FT is indeed interesting, but this software of converting a given spectrum into a FT is very slow. CSI: FingerID (Dührkop
et al., 2015; Shen et al., 2014) used both mass spectra and FTs as input, thinking that structural information of chemical compounds
can be captured by FTs. Indeed using FTs might be similar to considering peak interactions. However if FTs are used as input features,
spectrum must be converted into a FT not only in training andbut also in prediction, which needs a heavy computational load. iii)
The size/length of a fingerprint is not short, while the number ofpeaks is usually small, meaning that most elements of a fingerprint
vector are zero. So the training data is very sparse, while sparsemodels have not been considered yet. In addition, sparse models areadvantageous in that their results are easily interpretable.
Motivated by these drawbacks, we address the following two
problems: i) incorporation of peak interactions into the learningmodel to improve the predictive performance; ii) introduction ofsparsity into the models for interpretation. For the above i), we

Fig. 2. A general scheme to identify unknown metabolites based on molecular fingerprint vectors. There are two main stages: 1) fingerprint prediction: learning amapping from a molecule to the corresponding binary molecular fingerprint vector by classification methods, given a set of MS/MS spectra and fingerprints; 2)
candidate retrieval: using the predicted fingerprints to retrieve candidate molecules from the databases of known metabolites. Note that the step of constructing
an affinity matrix is optional and is used in L-SIMPLE only

Downloaded from https://academic.oup.com/bioinformatics/article-abstract/34/13/i323/5045791 by Kyoto Daigaku Bungakubu Toshokan user on 17 June 2019the query spectrum and the best matched candidates are returned.
However, the main disadvantage of this approach is that the reference database often contains only a small fraction of molecules inreality, leading to unreliable matching results if the molecule of thequery spectrum is not in the database. Consequently, to mitigate theinsufficiency of the reference databases, alternative approaches formetabolite identification are devised to deal with unavailable measured reference spectra. In general, we can divide computational
methods for metabolite identification into the following three categories: i) searching in spectral libraries; ii) in silico fragmentation;iii) machine learning (ML).
We focus on the ML approach, where the common scheme is,given a set of MS/MS spectra, to learn a map from a MS/MS
spectrum to a molecule (Fig. 2). This has two steps: i) fingerprint
prediction: predict a fingerprint with supervised ML, and ii) candidate retrieval: use the predicted fingerprint to query the database.

D.H.Nguyen et al.

SIMPLE: Sparse Interaction Model for fingerprint prediction

2 Related work
In Figure 2, the standard data preprocessing converts spectra into
high-dimensional feature vectors by dividing m/z range into bins andtaking accumulated intensity within each bin as a feature value.
However the width of bins is hard to determine. In fact, wide bins
can cause noise, and narrow bins can induce alignment errors due tomass error. One idea to overcome this issue is using a kernel, sayprobability product kernel (PPK; Jebara et al., 2004), which can becomputed directly from mass spectra. PPK treats each peak in a
spectrum as a two-dimensional Gaussian distribution and a spectrum as a uniform mixture of these distributions. Then, kernel between two spectra is computed by all-against-all matching between
the component Gaussians. The detail is in Section 3.1.1.
Also we briefly explain the current cutting-edge ML method by
MKL, CSI: FingerID (Dührkop et al., 2015; Shen et al., 2014), where
the input is both MS/MS spectra and FTs. Motivation behind thismethod is to use structural information, which might be captured by
FTs. So they use numerous kernels for FTs to capture any informationof FTs, ranging from simple ones, such as node kernels: node binary(NB), node intensity (NI) and edge kernels: loss binary (LB), loss count(LC), loss intensity (LI) to more complex one like common path counting (CPC) which is one of path based kernels. Subsequently, these kernels are combined by MKL methods, such as centered alignment(Gönen and Alpaydin, 2011), quadratic combination and lp norm
regularized combination (Kloft et al., 2011). The combined kernel isthen used in learning the final model for fingerprint prediction.
We again emphasize that these methods consider mainly only
peaks in MS/MS spectra without explicitly taking peak interactions
into account. Also kernels using FTs have the limitation of processing MS/MS spectra, particularly for prediction, due to the need ofcomputationally heavy conversion of spectra into FTs. More importantly, kernels are difficult to interpret and deal with sparse data,such as MS/MS spectra, where each spectrum (fingerprint vector)
has only a few number of peaks (nonzeros).

3 Materials and methods
3.1 Kernel method
We develop kernel for peak interactions and combine this kernel
with PPK (kernel for peaks) (Heinonen et al., 2012) through theframework of MKL.

3.1.1 Preliminary: kernel for peaks
MS/MS spectra have peak information, i.e. m/z and intensity. The
problem is that peaks are not correctly aligned each other, due tomeasurement errors in mass spectrometry devices (alignment error).
To overcome this problem, PPK (Jebara et al., 2004) was developed
to calculate the similarity between two spectra. The idea is that apeak can be considered as a two-dimensional Gaussian distribution
and the spectrum is a uniform mixture of peaks (or distributions).
The kernel between two spectra is computed by all-against-allmatching between the Gaussians.
More specifically, given a mass spectrum S, being a list of peaks,
i.e. fðm1 ; I1 Þ; ðm2 ; I2 Þ; . . . ; ðmNS ; INS Þg. The kth peak of the spectrum
S is represented by a Gaussian distribution pk centered around thepeak measurement ðmk ; Ik Þ and with covariance shared with all
peaks: R ¼ diag r2m ; r2I , where r2m and r2I are the variances forthe mass and intensity, respectively. Hence, the spectrum S can be
represented as a uniform mixture of the peaks contained in it,
P S
0i.e. pðSÞ ¼ N1S N
k¼1 pk . Likewise, another spectrum S is also represented by a mixture of distributions ql, l ¼ 1; 2; 3; ::; NS0 and
PNS0
pðS0 Þ ¼ N1 0 l¼1
ql .
S
With definitions above, the kernel for peaks, Kpeak , between two
spectra S and S0 is given by:ð
Kpeak ðS; S0 Þ ¼
pS ðxÞqS0 ðxÞdx
(1)
R2

¼

X


1
K pi ; qj
NS NS0 1  i  N(2)S

1  j  NS0

where Kðp; qÞ is the PPK between two component Gaussian distributions p and q, computed by (3).
0
0
2 
2 11
m
m
lIp  lIq CC
1B
B lp  lq
Kðp; qÞ ¼
exp@0:5@
þAA (3)
4prm rI
rm
rI
I
where lmp and lp denote the mass and intensity values of the peaks
corresponding to p.

3.1.2 Kernel for peak interactions
Being rather straight-forward, kernel for peak interactions between
two spectra S and S0 can be defined as follows:X
Kinteraction ðS; S0 Þ ¼
Kðpi ; qk Þ  K pj ; ql ;(4)
1  i  j  NS
1  k  l  NS0
where K is the function defined by (3).
Note that Kinteraction also can overcome the alignment error problem by taking advantage of Kpeak . Intuitively, Kpeak can be considered as a probabilistic version of the number of common peaks.
Similarly Kinteraction can be considered as a probabilistic version of
the number of common edges or interactions between two spectra.
3.1.3 Combining kernels for peaks and peak interactions
We combine kernels for peaks and peak interactions by using a
regular approach: Multiple Kernel Learning (MKL; Gönen and
Alpaydin, 2011), which can combine kernels from different sources.
We use three approaches: the first is the uniform combination of thekernels (UNIMKL: the weights for kernels are equal), which can
produce good results for prediction in many practical applications.
The second and third are two different MKL algorithms: centered

Downloaded from https://academic.oup.com/bioinformatics/article-abstract/34/13/i323/5045791 by Kyoto Daigaku Bungakubu Toshokan user on 17 June 2019

propose a kernel for peak interactions and combine this kernel withother kernels through multiple kernel learning (MKL). For the above
both i) and ii), we propose a sparse, interpretable model, which wecall Sparse Interaction Model over Peaks of moLEcules (SIMPLE).
Additionally, we also propose a FT-induced Laplacian regularization to make SIMPLE more robust. We note that in SIMPLE, FTs
are used for regularization only and not for input, by which we do
not need FTs for prediction, meaning that computationally SIMPLE
is much lighter for prediction than (Shen et al., 2014). We formulate
the model as a convex optimization problem, for which we develop
an alternating direction method of multipliers (ADMM) algorithm.
We evaluated our two proposed models by using real dataobtained from the MassBank dataset (Horai et al., 2010), which
was used in (Shen et al., 2014). We found that incorporation ofpeak interactions can significantly improve the prediction accuracyof not using interactions, resulting in comparable performance withthe current top method. Furthermore, SIMPLE could show the interpretability of results, i.e. peaks and peak interactions which contribute to high predictive performance.

i325

i326

D.H.Nguyen et al.

Ka ¼

m
X

aj Kj

3.2.1 Prediction model
Given a MS/MS spectrum, represented by a feature vector,
x ¼ ½x1 ; x2 ; . . . ; xd T 2 Rd , we formulate the model for individual
peaks and peak interactions as follows:f ðx; w; W Þ ¼ b þ

d
X
i¼1

wi xi þd X
dX

Wij xi xj¼ b þ wT x þ xT Wx

(5)

(11)

i¼1 j¼1(12)

j¼1

A popular approach to MKL, which is called two-stage MKL, separate kernel learning from prediction learning. In the first stage, thecombined kernel is learned through the use of an objective function,such as centered alignment that measures the similarity of two kernels over the training dataset. For the purposes of MKL, alignmentis often measured between the combined kernel in (5) and the target
kernel Ky ¼ yyT . Subsequently, the learned combined kernel is usedin the second stage with kernel based classifiers such as relevance
vector machine (RVM) used in (Burden and Winkler, 2015) or SVM
used in our experiments.
Formally, the centered kernel of a given kernel K 2 Rnn isdefined by:

 

eeT
eeT
Kc ¼ I 
K I
(6)
n
n

where b 2 R; w 2 Rd and W 2 Rdd . The prediction function consists of a bias b and two terms: main effect term parameterized by
the weight vector w and interaction term parameterized by theweight matrix W. Their roles are different. While the former capture
information about the peaks, the latter captures information aboutpeak interactions. Since our task is classification which predicts thepresence or absence of properties in fingerprint vector, the output of
the model can be computed by yðxÞ ¼ signðf ðx; w; W ÞÞ 2 f1; 1g.
In order to predict a binary vector of fingerprint, we predict one
response variable (hereafter a property or a task) at each time with aseparate classifier. Let consider a property and give a set of n training input/output pairs fðxi ; yi Þgni¼1 , where xi 2 Rd represent the ith
spectrum and yi 2 f1; 1g indicates the presence (þ1) or absence(-1) of the property. We can learn the parameters by minimizing thefollowing optimization function:

where I is the identity matrix and e is the vector with all ones.
The centered alignment between two kernels K and K0 is defined by:qðK; K0 Þ ¼< Kc ; K0c >jjKc jjF jjKc jjF(7)

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃwhere < Kc ; K0c > ¼ trace KTc K0c and jjKc jjF ¼ trace KTc Kc (8)Cortes et al. (2012) proposed two MKL algorithms: i) a simple
centered alignment algorithm (ALIGN), which assigns alignmentscores computed in (7) to combination weights in (5), i.e.

aj ¼ q Kj ; KY ; j ¼ 1; 2; 3; . . . ; m. and ii) the alignment maximum
algorithm (ALIGNF) seeks the weights a by maximizing the alignment scores between the combined kernels Ka and target kernel Ky ,resulting in the objective function:< Ka ; Ky >jjKa jjF(9)

subject to jjajj ¼ 1 ; a  0(10)

a ¼ argmax
a

min
b;w;W

3.2 Sparse Interaction Model over Peaks of moLEcules(SIMPLE)
Kernel-based methods are difficult to deal with sparse data and lackof interpretation. Thus we present a more interpretable, fast, sparselearning: Sparse Interaction Model over Peaks of moLEcules(SIMPLE), to incorporate peak interactions explicitly. We first preprocess each MS/MS spectrum to generate a feature vector: we divide m/z range into bins and taking accumulated peak intensities in a
bin as a feature value to obtain high dimensional vector for each
spectrum. These feature vectors are normalized such that all feature
values are in [0, 1].

(13)

i¼1

where lðy; ybÞ is the hinge loss function and computed by
lðy; ybÞ ¼ ð1  y; by Þþ , in which the operator ðzÞþ denotes maxð0; zÞ.
RðwÞ and RðW Þ are the regularization terms for vector w and matrix W, respectively.
Our purpose is to seek a model which is accurate as wellas interpretable. Different from SVM using l2 -norm regularizer,RðwÞ ¼ ajjwjj22 , our model uses sparsity-induced regularizer(Tibshirani, 1994), RðwÞ ¼ ajjwjj1 to yield sparse solution andinterpretation. As for interaction term, it is natural to impose lowrankness on matrix W due to the existence of groups of peaks interacting with each other. Thus, we propose to use trace norm, similarto (Blondel et al., 2015), i.e. RðW Þ ¼ bjjWjj . Furthermore, due tothe symmetry of matrix W, we also impose the positive semidefiniteness on matrix W, i.e. W Ɒ 0. Therefore, putting all above things together, the objective function of SIMPLE becomes:minIn our experiments, we conducted comparative experiments byusing these algorithm for combining kernels for peaks andinteractions.

nXlðyi ; f ðxi ; w; W ÞÞ þ RðwÞ þ RðW Þb;w;Wsubject tonXlðyi ; f ðxi ; w; W ÞÞ þ ajjwjj1 þ bjjWjj
i¼1

(14)

W Ɒ0

where a and b are hyperparameters to control the sparsity of w andlow-rankness of W, respectively. Note that (15) is the convex function which guarantees to find the global optimal solution.
3.2.2 Fragmentation trees (FTs) as prior information
As the number of interactions is large, it would be a good idea toregularize (14) with background knowledge. We propose to regularize the interaction matrix W by FTs of the spectra. However, there
are two questions to deal with: i) How to represent the structural information from FTs; ii) how to incorporate them into the convex objective function (14) while preserving its convexity to guarantee theglobal optimal solutions.

Downloaded from https://academic.oup.com/bioinformatics/article-abstract/34/13/i323/5045791 by Kyoto Daigaku Bungakubu Toshokan user on 17 June 2019alignment (ALIGN) and alignment maximization algorithm(ALIGNF) (Cortes et al., 2012).
These algorithms have the same setting. Given a set of kernels
fKj 2 Rnn ; j ¼ 1; ::; mg, computed from n data points and representing different sources of information. The vector y ¼ f1; 1gn
corresponds to the output or labels of data points. The aim of MKL
is to seek a combination of these kernels, as defined by:

SIMPLE: Sparse Interaction Model for fingerprint prediction

i327

xT Wx ¼d X
vTi vj xi xj ;i;j¼1

where vi and vj are the ith and jth bins of spectrum x and correspondto representation of ith and jth features in the space Rk .
We assume that if ith and jth peaks are strongly interacting witheach other, then their representation in the vector space should beclose, i.e. jjvi  vj jj2 should be small. From this observation, we include the following term, Laplacian smoothness, into the objectivefunction (14):RðV Þ ¼

d
X

Aij jjvi  vj jj22



¼ trace V T LV

(16)



¼ trace VV T L ¼ traceðWLÞ

(17)

b;w;W

i¼1

(18)

subject to W Ɒ 0
where f i ¼ f ðxi ; w; W Þ, c is an additional hyperparameter to control
Laplacian smoothness of the objective function (18). It is noteworthy that, our derived formulation is still convex due to convexityof regularization terms, thus, the solution of (18) is guaranteed to be
the global optimal. Below we will present an optimization method
for efficiently solving (18).
3.2.3 Optimization: ADMM algorithm
Equation (18) is convex and guaranteed that the algorithm converges to the global optimal solution. However, it is challenging tosolve this problem directly, because terms are nondifferentiable.
Our optimization method is based on alternating direction methodof multipliers (ADMM). Due to the nondifferentiability of the hinge
loss, we introduce the auxiliary variable C, where C ¼
ðC1 ; C2 ; . . . ; Cn ÞT and Ci ¼ 1  yi f i , then (18) can be reformulated
into the following equivalent constrained problem:
minb;w;W;Csubject tonXðCi Þþ þ ajjwjj1 þ bjjWjj þ ctraceðWLÞ
i¼1

C ¼ 1  YF and W Ɒ 0

where F ¼ ½f 1 ; f 2 ; . . . ; f n T .

Lðw0 ; w; W; C; uÞ ¼nX
ðCi Þþ þ ajjwjj1 þ bjjWjj þ ctraceðWLÞ

1þuT ð1  YF  CÞ þ jj1  YF  Cjj22
2
(20)
T

where L ¼ D  A, D is the degree matrix of A, i.e.
P
Dii ¼ nj¼1 Aij ; i ¼ 1; ::; nThus, we can formulate the following optimization problem,
which we call L-SIMPLE:n
Xlðyi ; f i Þ þ ajjwjj1 þ bjjWjj þ ctraceðWLÞ;The augmented Lagrange function of (19) is defined by:

i¼1

(15)

i;j¼1

min

Fig. 3. Illustration of constructing affinity matrix A from the set of fragmentation trees. The constructed matrix A is used as prior information for regularizing interaction matrix W(19)

where u ¼ ½u1 ; u2 ; . . . ; un  is a dual variable corresponding to theconstraint C ¼ 1  YF. Note that the constraint WⱰ0 is notincluded in (20) because this property will be imposed automaticallyon W after each iteration of updating W, as explained in the
Appendix section. We solve the problem of finding the saddle pointðb ; w ; W  ; C ; u Þ of the augmented Lagrangian function (20)through an iterative algorithm between the primal and the dual optimization as follows:8btþ1 ; wtþ1 ¼ argminb;w Lðb; w; W t ; Ct ; ut Þ
>
>
>
>


>
>
< W tþ1¼ argminW L btþ1 ; wtþ1 ; W; Ct ; ut
(21)


>
>
¼ argminC L btþ1 ; wtþ1 ; W tþ1 ; C; ut
Ctþ1
>
>
>
>
: tþ1
¼ ut þ 1  YFtþ1  Ctþ1u
where the first three steps update the primal variables based on thecurrent estimate of the dual variable ut and the final step updates
the dual variable based on the current estimate of the primal variables. Note that the efficiency of ADMM for solving (20) depends onwhether the subproblems in (21) can be solved quickly. Algorithm 1summarizes the ADMM steps for solving the optimization problem(20). To avoid confusion, the detailed derivation of the update rules
for subproblems of (21) are in Appendix.

3.3 Model summary
We here summarize the advantageous features of (L-)SIMPLE:
1. Peak interactions: SIMPLE has, in its formulation, the term forpeak interactions explicitly, which has not been considered byexisting methods, particularly kernel-based methods.
2. Sparse interpretability: MS spectra are with many zeros andsparse data. We formulate SIMPLE as a sparse model, by whichpeaks or peak interactions which contribute to improve the predictive performance can be checked and found easily.

Downloaded from https://academic.oup.com/bioinformatics/article-abstract/34/13/i323/5045791 by Kyoto Daigaku Bungakubu Toshokan user on 17 June 2019

To answer the first question, our idea is to construct an affinity (adjacency) matrix A for all peaks in the spectra. Concretely, each spectrum is converted into a FT by algorithms in (Rasche et al., 2011). Our
assumption is that, the frequency of an interaction present in the fragmentation trees reflects how strongly the corresponding features areinteracting with each other. From that, we calculate the co-occurrenceof all peak pairs in the trees to construct the affinity matrix (Fig. 3).
To answer the second question, we impose the constraint ofbeing positive semidefinite on interaction matrix W. By this, W canbe decomposed into the low rank matrix, i.e. W ¼ VV T where
V 2 Rnxk ; k ¼ rankðW Þ. Thus, the interaction term in the prediction
function (12) can be rewritten as following:

i328

D.H.Nguyen et al.

Table 1. Micro-average performance of kernels: PPK (Heinonen
et al., 2012) is used to compute Kpeak
Acc (%)

F1-score (%)

PPKComUNIMKLComALIGN
ComALIGNF

75.74 (68.13)
78.41 (66.82)
78.57 (66.24)
79.03 (67.89)

60.59 (613.75)
65.05 (612.16)
65.34 (611.99)
65.67 (613.02)

Note: ComUNIMKL, ALIGN, ALIGNF are combinations of Kpeak andKinteraction by algorithms UNIMKL, ALIGN, ALIGNF, respectively.

3. Convex formulation: our formulation keeps SIMPLE (L-SIMPLE) a
convex model, which guarantees to find the global optimum. We
have developed an alternating direction method of multipliers(ADMM) algorithm, to realize the detection of the global optimum.
4. No fragmentation trees in prediction: L-SIMPLE uses fragmentation trees for regularization in training, and they are not
inputs, meaning that we do not need fragmentation trees in prediction, which can avoid heavy computational cost of generating
fragmentation trees from spectrum.

4 Experimental evaluation
Our focus is to incorporate peak interactions of mass spectra into
fingerprint prediction, and to build sparse models for model interpretability, and so we conducted experiments to answer the following two questions:••(Q1): Can peak interactions due to sets of correlated peaks inspectra be used to predict fingerprint vectors more accurately?
(Q2): For the purpose of interpretation, how to identify a smaller
subset of predictors (i.e. peaks or peak interactions) that exhibit
the strongest effects on fingerprints?

4.1 Data, preprocessing and evaluation measures
MassBank (Horai et al., 2010) was used. We used the same dataset with
402 compounds as (Shen et al., 2014), and followed the same preprocessing steps as in (Shen et al., 2014): for each compound, peaks recorded
from different energies were merged for each MS/MS spectra. Then, spectra were normalized such that the sum of intensities is up to 100%. In
terms of the output for the learning models, molecular fingerprints, which
are binary vectors of 528 bits in total, were generated using OpenBabel(O’Boyle et al., 2011). Fingerprints have a high class imbalance, i.e.
mostly þ1 or reverse. Thus we only trained models for predicting fingerprints in which the majority class occupies less than 90% of instances.
We used micro-average accuracy, F1 score to evaluate the performance of different methods, computed by taking the average ofaccuracies and F1 scores over all tasks.

4.2 Benefit of incorporating interaction
We show the benefit of incorporating peak interactions using kernel
methods described in Section 3.1. The MKL algorithms, i.e. UNIMKL,ALIGN and ALIGNF, were used to combine two kernels, Kpeak
and Kinteraction (we call their results ComUNIMKL, ComALIGN,ComALIGNF, respectively) and compared with PPK, which has only
kernel for peaks. The combined kernels were coupled with SVM to predict the fingerprint properties. Each property was separately trained bya classifier. Five-fold cross-validation was conducted to seek suitable
margin parameter C where C 2 f23 ; 22 ; . . . ; 26 ; 27 g.
Table 1 shows the results, in which the micro-average accuracy
and F1 score of combined kernels were higher than PPK. This shows

1: Inputs:
T

A set of MS/MS spectra X ¼ xT1 ; xT2 ; . . . ; xTn
andassociated output Y ¼ ½y1 ; y2 ; . . . ; yn T . Laplacian L(only for L-SIMPLE).
2: Outputs:
weight vector w and interaction matrix W.
3: Initialize:
b
0; w
0; W
0; C
0; u
0
4: while not converged do
5:
1. Fix W, C, u and update b, w


b1
6:
Precompute FYð1  C þ uÞ  diag XWXT
7:
while not converged do
8


b1
>
z
w  qXT Xw þ b  F
>
>
>
>
>
<
w
S a ðzÞ
>>1>b1  Xw
>
b
sum F
>
>
:
n
8:
9:
10:
11:

end while
2: Fix b, w, C, u and update Wb2
Yð1  C þ uÞ  b  Xw
precompute F
while not converged do
8


b2
>
R
diag XWXT  F
>
>
>
>
>
T
>
DW L
X RX
>
>
>
>
>
<Z
W  qðDW L þ cLÞ>U; EEVDðZÞ
>
>
>
>
>
>
b
>
E
S b ðEÞ
>
>
>
>
:
b T
W
UEU

12:
end while13:
3: Update C
14:
F ¼ b þ Xw þ diagXWXT ; C
15:
4: Update the dual variable u16:uu þ 1  YF  C
17: end while

T 1 ð1  YF þ uÞthat incorporating peak interactions can improve the performanceon predicting molecular fingerprint properties. Additionally, the performance of ALIGNF algorithm was slightly better than the rest.

4.3 Benefit of (L-)SIMPLE, sparse interaction models
Since kernel methods in Section 3.1 are unable to provide sparse solutions for interpretability, we examined (L-)SIMPLE to gain insight
into the models to predict fingerprints. (L-)SIMPLE needs to construct feature vectors from MS/MS spectra: we divided m/z range
into 500 bins and took accumulated peak intensities in a bin as afeature value to obtain high dimensional vector for each spectrum.
These feature vectors were normalized such that all feature values
are in [0, 1]. (L-)SIMPLE was trained by ADMM (see algorithm 1)
with all variables initialized at zero. The algorithm was iterated until
the relative difference in training errors fell below 0.0001 or thenumber of iterations reaches 100. Five-fold cross-validation wasused to evaluate the generalization of the learning machines.
Specifically, parameters a, b and c, for controlling sparsity,

Downloaded from https://academic.oup.com/bioinformatics/article-abstract/34/13/i323/5045791 by Kyoto Daigaku Bungakubu Toshokan user on 17 June 2019

MethodAlgorithm 1 ADMM algorithm for optimizing the objectivefunction (19).

SIMPLE: Sparse Interaction Model for fingerprint prediction

i329

Table 2. Performance comparison between SIMPLE and L-SIMPLE
Task Id/NameSIMPLE
Acc (%)

F1 (%)

Acc (%)

F1 (%)

71.16
91.24
79.36
80.35
80.35
82.35
83.83
69.38
73.12
71.88
81.35
78.33 6 6.05

69.24
95.19
53.74
54.18
71.39
70.98
81.67
62.00
40.42
67.15
23.49
66.69 6 13.03

73.14
90.29
79.85
81.84
80.61
82.6
83.34
69.65
73.88
72.39
84.08
78.86 6 5.87

70.75
94.82
53.98
56.17
72.03
71.16
82.16
62.15
44.46
67.28
27.75
67.59 6 12.35

Table 3. Micro-average performance and computation time (forprediction) of kernel-based methods in Shen et al. (2014) and proposed methods in this paper
MethodAcc(%)

F1 score(%)

Run. time(ms)

PPK (Peaks)
LB (Loss binary)
LC (Loss count)
LI (Loss intensity)
NB (Node binary)
NI (Node intensity)
CPC (Common pathcount)ComFT (combining
all above)
ComALIGNF (Proposed:
MKL)
SIMPLE (Proposed)
L-SIMPLE (Proposed)

75.74 (66.72)
76.63 (67.03)
75.33 (65.4)
74.54 (68.49)
79.11 (65.02)
78.41 (64.99)
79.02 (67.4)

60.59 (614.54)
61.64 (615.48)
61.25 (613.99)
58.46 (616.01)
67.34 (611.75)
66.87 (612.11)
67.55 (612.93)

52.37
1501.021501.021501.021501.091501.01
1501.1180.98 (66.05)

69.04 (611.98)

1559.2079.03 (67.89)

65.67 (613.02)

471.7178.33 (66.05)
78.86 (65.87)

66.70 (613.03)
67.59 (612.35)

4.57
4.32

low-rankness and Laplacian smoothness were chosen from the lists:{1, 2, 3}, f2; 3; 4; 5g and f0:1; 0:5; 1:0g, respectively.
Also, to evaluate the effects of adding the Laplacian smoothness
term into the objective function, we compared the accuracy and F1score over tasks by performing five-fold cross-validation. Table 2
shows the accuracy and F1 score obtained for the first ten tasks(Note that we used only tasks in which the majority class occupies
less than 90% of instances). The micro-average over all tasks is also
displayed at the bottom. We can see that Laplacian regularization
worked to make SIMPLE more robust, resulting in better predictive
performance of L-SIMPLE.
We further compared (L-)SIMPLE with various kernel, namely
PPK, NB, NI, LB, LC, LI, CPC and their combination, which we call
ComFT, all from (Shen et al., 2014). While these kernels (except for
PPK) are all computed from the fragmentation trees, in which thecost for converting MS/MS spectra to these trees is heavy and timeconsuming, our method uses peaks from spectra only and is efficiently computable in prediction.
Table 3 shows the results of accuracy, F1 scores and computation time for prediction by all compared methods. The prediction
time (in milliseconds) was averaged over all spectrum in the dataset.
The first four methods including PPK achieved around the accuracyof 75%, which is clearly worse than the other methods, which

achieved around 78–80% and are very comparable each other. In
fact ComFT, the current cutting-edge MKL-based method, performed best in both accuracy and F1 score, while the second best
was not clear (NB by accuracy and L-SIMPLE by F1 score). On theother had, about computation time, FT-based methods, i.e. from LB
to ComFT, were clearly slower than the others, because they needconvert spectra into FTs. In fact ComFT needed more than 1500 ms,
which is more than 300 times slower than that of (L-)SIMPLE,
which just spent only less than 5 ms. This is a sizable difference
when we have to process a huge amount of spectra produced by thecurrent high-throughput MS/MS. We stress that the performancedifference between (L-)SIMPLE and ComFT was very slight andstatistically insignificant, while (L-)SIMPLE was exceedingly fasterthan ComFT.

4.4 Model interpretation
One advantage of sparse learning models over kernel based methods
is interpretation. For illustration purposes, Figure 4 shows theweights of main effect terms (w) and the interaction weight matrix(W) obtained by L-SIMPLE for nine fingerprint properties (or tasks):29, 37, 56, 70, 139, 192, 236, 356 and 370 (these are randomly
selected for investigation). As observed, the weight of main effectand interaction terms were different between properties, suggestingthat different properties are strongly affected by different subsets ofa few peaks in spectra.
Table 4 shows case studies to illustrate the effects of peak interactions. We consider four interaction pairs frequently present inthese tasks. w1 and w2 denote the weights for peaks and W12
denotes the weight for their interactions. We can raise the following
three interesting findings:
1. Either w1 or w2 (or both) can be zero but the interaction weights
are often nonzero: For example, W12 of interaction (42, 85) aremostly nonzero, while w1 and w2 are both zero with respect totasks 29, 37, 56, 70. This means that individual peaks are notgood predictive descriptors of properties, while their interactions
are. Thus this result clearly shows the importance of consideringpeak interactions in fingerprint prediction.
2. Despite of negative impacts of individual peaks, their interactions can be positive: For example, interaction (85, 227) with respect to tasks 29, 56 and 366. This means that while individual
peaks indicate the absence of a property, the interactions of two
peaks mean the presence of the property.

Downloaded from https://academic.oup.com/bioinformatics/article-abstract/34/13/i323/5045791 by Kyoto Daigaku Bungakubu Toshokan user on 17 June 2019

3 (Aldehyde)
27 (Hydroxy)
29 (Primary alcohol)
30 (Secondary alcohol)
37 (Ether)
38 (Dialkyl etherEther)
45 (Aryl)
50 (Carboxylic acid)
56 (Primary Carbon)
57 (Secondary Carbon)
60 (Alkene)
Avg 6 Std

L-SIMPLE

i330

D.H.Nguyen et al.

(a)

(b)

Table 4. Case studies of weight vector w and interaction matrix W learned by L-SIMPLE over a set of randomly selected tasks
Tasks/Name(42, 85)
w1(42, 163)w2

29 (Primary Alcohol)
0.0
0.0
37 (Ether)
0.0
0.0
56 (Primary Carbon)
0.0
0.0
70 (Alkylarylether)
0.0
0.0
139 (Thioenol)
0.1575 0.3939
192 (Carbonic acid monoester) 0.1542 0.0
236 (Heteroaromatic)
0.0
0.0
356 (1,5-Tautomerizable)
0.0
0.0
366 (Actinide)0.0
0.0(85, 227)

W12

w1

w2

W12

0.00160.01200.02710.01590.0
0.0
0.01070.01700.0245

0.05450.0
0.04420.02600.0
0.0471
0.5657
0.0
0.00470.0
0.0
0.0551
0.2308 0.4094 0.0
0.0
0.0
0.0
0.0
0.2499 0.0537
0.0
0.0
0.0301
0.1891
0.6065 0.0282(130, 201)

w1

w2

W12

w1

w2

W12

0.40850.0
0.9833
0.1238
0.35890.6945
0.3069
0.5539
0.5373

0.05450.02600.5657
0.0
0.2308
0.0
0.0
0.0
0.1891

0.07650.12640.0271
0.09720.0
0.0
0.10620.06070.0399

0.40850.0
0.02180.0
0.0
0.3389
0.9833
0.0
0.0104
0.12380.0
0.0265
0.3589 0.1126 0.0
0.69450.0
0.0
0.30690.1401 0.0201
0.55390.0
0.0107
0.5373
0.0
0.0153

Note: w1 and w2 denote weights corresponding two mass positions and W denotes the weight of their interactions. Four pairs of mass positions which are frequently present in these tasks, including (42, 85), (42, 163), (85, 227) and (130, 201) are shown.

3. Interaction of peaks can be zero, indicating these peaks are independent of each other even if their corresponding weights are nonzeros: For example, interaction (42, 85) with respect to task 139.
They are just part of numerous findings, but even these examples
show a fact that even individual peaks are not good indicators ofsome fingerprints, their interactions with others may significantly
contribute to the prediction. This again confirms the importance ofpeak interactions for fingerprint prediction.

5 Discussion and conclusion
The goal of this work is to propose machine learning models which
are able to incorporate peak interactions for fingerprint prediction.
Our experiments showed that peak interactions are definitely useful
to improve fingerprint prediction, along with discriminative information about peaks.
Our first model is based on kernel learning, defining two kernels,
one for peaks and the other for peak interactions, which are combined
through MKL. Again we note that Shen et al. (2014) usedfragmentation trees for prior structural information of spectral, andconverting spectra into fragmentation trees is definitely computationally expensive. On the other hand, our model of kernel learning uses
only peaks in the spectrum as input for prediction, indicating that ourmodel is much more efficient. Kernel learning does not have to construct feature vectors for spectra, and instead this is done implicitly bykernels defined, which can avoid any error caused when generating
feature vectors. However, a big issue of kernel learning is interpretability. That is, it is difficult for kernel learning to figure out which
subset of peaks or peak interactions exhibit the strongest effects on
fingerprint prediction, despite that clearly each property depends on avery few number of mass positions in each given spectrum.
Our sparse interaction models, (L-)SIMPLE, have a number ofadvantages, which are summarized in Section 3.3: (L-)SIMPLE isformulated as a sparse, convex optimization model, which cancapture peak interactions and also give interpretable solutions.
We emphasize that next generation fingerprint prediction needs a
ML model, which should learn, from huge but sparse spectra, peaks
as well as peak interactions comprehensibly and predict fingerprints

Downloaded from https://academic.oup.com/bioinformatics/article-abstract/34/13/i323/5045791 by Kyoto Daigaku Bungakubu Toshokan user on 17 June 2019

Fig. 4. (a) Weight vectors (w) of the main effect terms and (b) smooth heat map of weight matrices (W) of the interaction terms learned by L-SIMPLE for properties
or tasks: 29 (Primary alcohol), 37 (Ether), 56 (Primary Carbon), 70 (Alkylarylether), 139 (Thioenol), 192 (Carbonic acid monoester), 236 (Heteroaromatic), 356(1,5-Tautomerizable) and 366 (Actinide)

SIMPLE: Sparse Interaction Model for fingerprint prediction
against again huge spectra highly efficiently. (L-)SIMPLE would bea reasonable solution for this situation.

i331

A1. ADMM for updating b and w
The subproblem for optimizing b, w in (21) is equivalent tobtþ1 ; wtþ1

¼ argmin Lðb; w; W t ; Ct ; ut Þb;w

Funding

Conflict of Interest: none declared.

References
Blondel,M. et al. (2015). Convex factorization machines. In Appice,A. et al.
(eds.) Machine Learning and Knowledge Discovery in Databases. Springer
International Publishing, Cham, pp. 19–35.
Böcker,S. and Rasche,F. (2008) Towards de novo identification of metabolites
by analyzing tandem mass spectra. Bioinformatics, 24, i49–i55.
Burden,F.R. and Winkler,D.A. (2015) Relevance vector machines: sparse classification methods for qsar. J. Chem. Inf. Model., 55, 1529–1534.
Cai,J.-F. et al. (2010) A singular value thresholding algorithm for matrix completion. SIAM J. Optim., 20, 1956–1982.
Cortes,C. et al. (2012) Algorithms for learning kernels based on centeredalignment. J. Mach. Learn. Res., 13, 795–828.
de Hoffmann,E. and Stroobant,V. (2007) Mass Spectrometry, Principles andApplications, 3rd edn. John Wiley & Sons.
Dührkop,K. et al. (2015) Searching molecular structure databases with tandem mass
spectra using csi: fingerid. Proc. Natl. Acad. Sci. USA, 112, 12580–12585.
Gönen,M. and Alpaydin,E. (2011) Multiple kernel learning algorithms.
J. Mach. Learn. Res., 12, 2211–2268.
Heinonen,M. et al. (2012) Metabolite identification and molecular fingerprint
prediction through machine learning. Bioinformatics, 28, 2333–2341.
Horai,H. et al. (2010) Massbank: a public repository for sharing mass spectral
data for life sciences. J. Mass Spectrom., 45, 703–714.
Hummel,J. et al. (2010) Decision tree supported substructure prediction of
metabolites from gc-ms profiles. Metabolomics, 6, 322–333.
Imre,T. et al. (2008) Mass spectrometric and linear discriminant analysis ofn-glycans of human serum alpha-1-acid glycoprotein in cancer patients andhealthy individuals. J. Proteomics, 71, 186–197.
Jebara,T. et al. (2004) Probability product kernels. J. Mach. Learn. Res., 5,
819–844.
Kloft,M. et al. (2011) lp-norm multiple kernel learning. J. Mach. Learn. Res.,12, 953–997.
Ma,S. et al. (2011) Fixed point and bregman iterative methods for matrix rank
minimization. Math. Program., 128, 321–353.
O’Boyle,N.M. et al. (2011) Open babel: an open chemical toolbox.
J. Cheminf., 3, 33.
Rasche,F. et al. (2011) Computing fragmentation trees from tandem mass
spectrometry data. Anal. Chem., 83, 1243–1251. PMID: 21182243.
Scheubert,K. et al. (2013) Computational mass spectrometry for small molecules. J. Cheminf., 5, 12.
Shen,H. et al. (2014) Metabolite identification through multiple kernel learning on fragmentation trees. Bioinformatics, 30, i157–i164.
Tibshirani,R. (1994) Regression shrinkage and selection via the lasso. J. R.
Stat. Soc. Ser. B, 58, 267–288.
Watanabe,T. et al. (2014) Scalable fused lasso svm for connectome-based disease prediction. In: 2014 IEEE International Conference on Acoustics,
Speech and Signal Processing (ICASSP), pp. 5989–5993.
Wishart,D.S. (2007) Current progress in computational metabolomics. Brief.
Bioinf., 8, 279–293.
Wishart,D.S. et al. (2012) Hmdb 3.0 – the human metabolome database in
2013. Nucleic Acids Res., 41, D801–D807.

Appendix
We will derive the updating rules of ADMM for steps in (21).

1¼ argmin ajjwjj1 þ jj1  YFt  C þ ujj22
2
b;w

(22)

1
F 1 jj22¼ argmin ajjwjj1 þ jjXw þ b1  b
2
b;w


b1 ¼ Yð1  C þ uÞ  diag XWXT . We denote the operator
where Fdiagð:Þ to produce a vector composed of diagonal elements of given
square matrix.
In fact, solving (22) can be done easily with proximal gradientdescent through alternately updating estimate of w and b as follows:
8


b1
>
zkþ1 ¼ wk  dw XT Xwk þ bk  F
>
>
>
>
>
>
<


wkþ1 ¼ S a zkþ1
(23)
>
>
>
>


>
>
>
b1  Xwkþ1: bkþ1 ¼ 1 sum Fnwhere dw is the learning rate for the proximal gradient updates and set
to 0.01 in our experiments, k is the index of inner loop for optimizingb, w, S k ðtÞ is the element-wise soft-thresholding operator, defined by:
S k ðtÞ ¼ signðtÞmaxð0; jtj  kÞ

(24)

A2. ADMM for updating W
The subproblem of optimizing W in (22) is equivalent to:W tþ1 ¼ argminL btþ1 ; wtþ1 ; W; Ct ; ut
W

¼ argminbjjWjj þ ctraceðWLÞ
W1
þ jj1  YFt  Ct þ ut jj22
2
¼ argminbjjWjj þ ctraceðWLÞ

(25)

W



1
F 2 jj22
þ jjdiag XWXT  b2
b2 ¼ Yð1  C þ uÞ  b1  Xw. It is obvious that thewhere Fobjective function (25) split in two components: ctraceðWLÞþ

T
1b 2 jj2 , which is convex, differentiable and bjjWjj ,F22 jjdiag XWX
which is also convex with inexpensive proximal operator. It is known
that the solution of (25) is given by the matrix shrinkage operation
which corresponds to a singular value decomposition (SVD) [see, e.g.
Cai et al. (2010) and Ma et al. (2011) for more details]. Hence, proximal gradient descent for updating W again is given as follows:
8


kþ1
b2
>
¼ diag XW k XT  F
>
>R
>
>
>
>
> DW L¼ cL þ XT Rkþ1 X
>
>
>
>
kþ1
>
>
¼ W k  d W DW L
>
<Z


(26)
Ukþ1 ; Ekþ1 ¼ EVD Zkþ1
>
>
>
>
>


>
>
>
b kþ1
>
E
¼ S b Ekþ1
>
>
>
>
>
>
: kþ1
b kþ1 Ukþ1W
¼ Ukþ1 E

Downloaded from https://academic.oup.com/bioinformatics/article-abstract/34/13/i323/5045791 by Kyoto Daigaku Bungakubu Toshokan user on 17 June 2019

This work was partially supported by MEXT KAKENHI Grant Number
16H02868, Grant Number JPMJAC1503 ACCEL JST, FiDiPro Tekes (currently Business Finland) and AIPSE Academy of Finland.

i332

D.H.Nguyen et al.
Proposition 1: the solution T k ðtÞ ¼ argminx2R kðxÞþ þ 12 ðx  tÞ2
has the following form:8t  k if t  k>><
T k ðtÞ ¼ 0
if 0  t  k>>:t
if t < 0(28)

Note that components of C are independent of each other in (27).

A3. ADMM for updating C

By applying Proposition 1 in element-wise, we can derive the update

For the subproblem of optimizing C in (22), it is equivalent to
Ctþ1 ¼ argmin L btþ1 ; wtþ1 ; W tþ1 ; C; ut

for C in the following closed form solution:Ctþ1 ¼ T 1 1  YFtþ1 þ ut

C

n
X
1
¼ argmin ðCi Þþ þ jj1  YFtþ1 þ ut  Cjj2
2
C
i¼1

(27)

In order to solve (27), we use the following proposition in(Watanabe et al., 2014):(29)

Downloaded from https://academic.oup.com/bioinformatics/article-abstract/34/13/i323/5045791 by Kyoto Daigaku Bungakubu Toshokan user on 17 June 2019

where DW L is the derivative of the differentiable component of (26)with respect to W, dW is the learning rate for the proximal gradient
updates and set to 0.05 in our experiments, k is the index of inner
loop for optimizing W. U and E are columns of eigenvectors and diagonal matrix of eigenvalues obtained by applying eigendecomposition (EVD) to Z, by which we can guarantee the semidefiniteness ofweight matrix W after each iteration.