%% BioMed_Central_Tex_Template_v1.06
%%                                      %
%  bmc_article.tex            ver: 1.06 %
%                                       %

%%IMPORTANT: do not delete the first line of this template
%%It must be present to enable the BMC Submission system to
%%recognise this template!!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                     %%
%%  LaTeX template for BioMed Central  %%
%%     journal article submissions     %%
%%                                     %%
%%          <8 June 2012>              %%
%%                                     %%
%%                                     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                 %%
%% For instructions on how to fill out this Tex template           %%
%% document please refer to Readme.html and the instructions for   %%
%% authors page on the biomed central website                      %%
%% http://www.biomedcentral.com/info/authors/                      %%
%%                                                                 %%
%% Please do not use \input{...} to include other tex files.       %%
%% Submit your LaTeX manuscript as one .tex document.              %%
%%                                                                 %%
%% All additional figures and files should be attached             %%
%% separately and not embedded in the \TeX\ document itself.       %%
%%                                                                 %%
%% BioMed Central currently use the MikTex distribution of         %%
%% TeX for Windows) of TeX and LaTeX.  This is available from      %%
%% http://www.miktex.org                                           %%
%%                                                                 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% additional documentclass options:
%  [doublespacing]
%  [linenumbers]   - put the line numbers on margins

%%% loading packages, author definitions

%\documentclass[twocolumn]{bmcart}% uncomment this for twocolumn layout and comment line below
\documentclass{bmcart}

%%% Load packages
%\usepackage{amsthm,amsmath}
%\RequirePackage{natbib}
%\RequirePackage[authoryear]{natbib}% uncomment this for author-year bibliography
%\RequirePackage{hyperref}
\usepackage[utf8]{inputenc} %unicode support
\usepackage{cite}
\usepackage{amsmath,amssymb,amsfonts}
\RequirePackage{natbib}
\usepackage{algorithmic}
\usepackage{graphicx}
\usepackage{textcomp}
\usepackage{xcolor}
\usepackage{epstopdf}
\usepackage{subfigure}
\usepackage{threeparttable}
\usepackage{multirow}
\usepackage{booktabs}
\usepackage{framed}
%\usepackage[applemac]{inputenc} %applemac support if unicode package fails
%\usepackage[latin1]{inputenc} %UNIX support if unicode package fails


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                             %%
%%  If you wish to display your graphics for   %%
%%  your own use using includegraphic or       %%
%%  includegraphics, then comment out the      %%
%%  following two lines of code.               %%
%%  NB: These line *must* be included when     %%
%%  submitting to BMC.                         %%
%%  All figure files must be submitted as      %%
%%  separate graphics through the BMC          %%
%%  submission process, not included in the    %%
%%  submitted article.                         %%
%%                                             %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\def\includegraphic{}
\def\includegraphics{}



%%% Put your definitions there:
\startlocaldefs
\endlocaldefs


%%% Begin ...
\begin{document}

%%% Start of article front matter
\begin{frontmatter}

\begin{fmbox}
\dochead{METHODOLOGY ARTICLE}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the title of your article here     %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{A Sequence Embedding Method For Enzyme Optimal Condition Analysis}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors here                   %%
%%                                          %%
%% Specify information, if available,       %%
%% in the form:                             %%
%%   <key>={<id1>,<id2>}                    %%
%%   <key>=                                 %%
%% Comment or delete the keys which are     %%
%% not used. Repeat \author command as much %%
%% as required.                             %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\author[
   addressref={aff1},
   noteref={n1},
   email={lxj35570747@gmail.com}
]{\inits{XJL}\fnm{Xiangjun} \snm{Li}}
\author[
   addressref={aff2},
   noteref={n1},
   email={bj600800@gmail.com}
]{\inits{ZXD}\fnm{Zhixin} \snm{Dou}}
\author[
   addressref={aff1},
   corref={aff1},
   email={sun\_yuqing@sdu.edu.cn}
]{\inits{YQS}\fnm{Yuqing} \snm{Sun}}
\author[
   addressref={aff2},
   email={lswang@sdu.edu.cn}
]{\inits{LSW}\fnm{Lushan} \snm{Wang}}
\author[
   addressref={aff1},
   email={gb@sdu.edu.cn}
]{\inits{BG}\fnm{Bin} \snm{Gong}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors' addresses here        %%
%%                                          %%
%% Repeat \address commands as much as      %%
%% required.                                %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\address[id=aff1]{%                           % unique id
  \orgname{School of Software, Shandong University}, % university, etc
  \street{Shunhua Road},                     %
  \postcode{250101}                                % post or zip code
  \city{Jinan},                              % city
  \cny{China}                                    % country
}

\address[id=aff2]{
  \orgname{State Key Laboratory of Microbial Technology, Shandong University},
  \street{Binhai Road},
  \postcode{266237}
  \city{Qingdao},
  \cny{China}
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter short notes here                   %%
%%                                          %%
%% Short notes will be after addresses      %%
%% on first page.                           %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{artnotes}
%\note{Sample of title note}     % note to the article
\note[id=n1]{Equal contributor} % note, connected to author
\end{artnotes}

\end{fmbox}% comment this for two column layout

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Abstract begins here                 %%
%%                                          %%
%% Please refer to the Instructions for     %%
%% authors on http://www.biomedcentral.com  %%
%% and include the section headings         %%
%% accordingly for your article type.       %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{abstractbox}

\begin{abstract} % abstract
\parttitle{Background} %if any
An enzyme activity is influenced by the external environment condition. It is important to have an enzyme remain high activity in a specific condition. A usual way is to first determine the optimal condition of an enzyme by either the gradient test or by tertiary structure, and then to use protein engineering to mutate a wild type enzyme for a higher activity in an expected condition. 
\parttitle{Results} %if any
In this paper, we investigate the optimal condition of an enzyme by directly analyzing the sequence. We propose an embedding method to represent the amino acids and the construct information as vectors in the latent space. These vectors contain information about the correlations between amino acids and sites in the aligned amino acid sequences, as well as the correlations with the optimal conditions. We crawled and processed the amino acid sequence in glycoside hydrolase GH11 family, and got 125 amino acid sequences with optimal pH condition. We used probabilistic approximation method to implement the embedding learning method on these samples. % with a machine learning method that relies on a large amount of data.
Based on these embedding vectors, we design a computational score to determine the optimal condition for an enzyme and achieves the accuracy 80\% on the test proteins in the same family. We also  give the mutation suggestion such that it has a higher activity in the expected environment, which is consistent with the professional wet experiments and analysis.  
%We can use the method to ?nd the amino acid sequence fragments that affect the optimal pH condition, and to give the direction of mutation in protein rational design to get a protein that meets expectations.
\parttitle{Conclusion} %if any
A new computational method is proposed for the sequence based enzyme optimal condition analysis. Compared with the traditional process that involves a lot of wet experiments and  requires multiple mutations, this method can get the desired protein for an expected condition in an efficient and effective way. 
\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The keywords begin here                  %%
%%                                          %%
%% Put each keyword in separate \kwd{}.     %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{keyword}
\kwd{Protein sequence analysis}
\kwd{Embedding}
\kwd{Bioinformatics}
\end{keyword}

% MSC classifications codes, if any
%\begin{keyword}[class=AMS]
%\kwd[Primary ]{}
%\kwd{}
%\kwd[; secondary ]{}
%\end{keyword}

\end{abstractbox}
%
%\end{fmbox}% uncomment this for twcolumn layout

\end{frontmatter}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Main Body begins here                %%
%%                                          %%
%% Please refer to the instructions for     %%
%% authors on:                              %%
%% http://www.biomedcentral.com/info/authors%%
%% and include the section headings         %%
%% accordingly for your article type.       %%
%%                                          %%
%% See the Results and Discussion section   %%
%% for details on how to create sub-sections%%
%%                                          %%
%% use \cite{...} to cite references        %%
%%  \cite{koon} and                         %%
%%  \cite{oreg,khar,zvai,xjon,schn,pond}    %%
%%  \nocite{smith,marg,hunn,advi,koha,mouse}%%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%% start of article main body
% <put your article body there>

%%%%%%%%%%%%%%%%
%% Background %%
%%
\section*{Background}
Proteins are made up of hundreds of monomers called amino acids that are attached to one another by peptide bonds, forming a long chain defined as primary structure. It further constructs the second structure and the tertiary structure of the protein, which finally determines the protein function and the optimal condition of its activity \cite{marks2012protein}, such as the acid resistance, or the salt tolerance.

Enzymes are a kind of catalytic proteins.  The enzymes with the same function are referred as a family, which have some identical or conserved amino acid fragments in  sequences that are not easy to be mutated. %As the results of natural selection, these parts are not easy to be mutated.
Although these enzymes have the same biological function, some of them are more active  than others at the same alkaline  or temperature condition, which are caused by the different parts in the sequences, defined as the non-conserved fragments. For example,  the enzymes in GH11 family can degradate the heteroxylans that constitute the lignocellulosic plant cell wall, but they have different optimal temperatures or pH \cite{paes2012gh11}. %There are also some different parts between these sequences that induce the different optimal conditions, defined as the non-conserved fragments.
In many practical applications, we need to find the most active enzyme in a family under a given  condition or to mutate an enzyme for a higher activity under an expected condition than the existing enzymes. For these purposes, biological researchers generally perform a series of gradient tests to measure the optimal condition of each enzyme in a family. Then they adopt the protein engineering to produce an enzyme with an expected  optimal condition \cite{wu2018enhancement}. But the above biological methods analyze only an enzyme at each wet experiment such that they cost much  time,  power and resources to get the expected enzyme.

Recently, many works adopt machine learning methods to predict the optimal temperature or pH for enzymes \cite{lin2011prediction}. For example, Dijk et al. use the ratio between the number of residues inside and outside the folding structure of protein as the measure for its hydrophobicity so as to induce the temperature dependence of the protein \cite{van2015hydrophobic}. Pucci et al. use temperature-dependent statistical potentials to compute the fold free energies at three different temperatures, and use Gibbs-Helmholtz equation to predict the protein thermal stability curve \cite{pucci2014stability}. The SCooP method predicts the full T-dependent stability curve by protein structures and host organisms \cite{pucci2017scoop}. However, these methods require the tertiary structure \cite{kc2016recent}, free energy \cite{liu2012importance}, etc., which are not easy to obtained in many cases.
Another representative method is based on the correlations between the tertiary structure and free energy \footnote{A physics concept to describe the force
that causes chemical reactions.} of an enzyme. It judges the stability of the mutant by analyzing the change of free energy in this process \cite{malinka2015prediction}. For example, Dehouck et al. use a linear combination of statistical potentials whose coefficients depend on the solvent accessibility of the mutated residue to predict the stability changes caused by single site mutations in proteins \cite{dehouck2011popmusic}. Wijma et al. create a library of potentially stable mutations by calculating the change in the free energy of mutations, reduce the size of the library by eliminating false positive predictions and choose the most stable combination of mutations to mutate \cite{wijma2014computationally}. These methods are based on molecular dynamics to simulate the state of proteins such that they perform a single protein in one experiment. But it can not analyze a group of proteins at the same time.
%Considering the tertiary structural of an enzyme is determined by its primary structure, it is feasible to  learn the  optimal condition directly from an amino acid sequence.

In this paper, we propose an embedding method to predicate the optimal condition of an enzyme directly stemming from the amino acid sequence. Amino acids and the construct information in sequence are represented as vectors in the same latent space, which are learned by the compatibility between the sites and expected condition. We propose the compatibility objectives on both a single-site with an amino acid and on multiple sites with different amino acids. Using these vectors, we can predict the optimal condition for a new amino acid sequence. We select the enzymes in the GH11 family as the practical usage, whose optimal pH are already determined. Since from the view of machine learning the quantity of samples is small, we use the statistics to approximate the probability distribution on amino acids on each site of the aligned sequences and generate more samples for the embedding training.

Based on these embeddings, we analyzed the non-conserved segment and the conserved segment of these enzymes. These results are consistent with the professional results on WebLogo\footnote{A sequence logo is created from a collection of aligned sequences and is widely adopted to depict the consensus sequence and diversity of the sequences, such as protein-binding sites in DNA or functional units in proteins. http://weblogo.threeplusone.com} \cite{crooks2004weblogo, schneider1990sequence}. We then design two experiments for biological purposes. One is to compare two enzymes which has a higher optimal pH or which has a higher activity in the same pH environment. Another is to quantify the probability on whether one or more  sites mutation on a given enzyme would result it still active in an expected optimal condition. The experimental results are consistent with the biological results. Comparing with other methods,  the embedding method is more efficient and effective. % and searching the enzyme for an optimal condition.

%The rest of this paper is organized as follows. We present the practical data used in this work and the preprocess on data in the Dataset section. The Results section discusses the application of the model and presents the experimental details. % We deliberate how to mutate amino acid sequence in the Discussion section. 
%After we make a conclusion of the full text,  we present our proposed sequence embedding method in the Method section at the end.

\section*{Notions and Dataset}\label{data_process}

\subsection*{Notions}
There are 20 kinds of amino acids known in nature. Let $aa$ denote  the set of amino acids in bioinformatics, i.e. $aa=\{Ala,Arg,Asp,Cys,$ $Gln,Glu,His, Ile,Gly,$ $Asn,Leu,Lys,Met,Phe,Pro,Ser,Thr,$ $Trp,Tyr,Val\}$.
%
Since the amino acid sequences in a family are often with different lengths, they are aligned before analysis. After alignment, the identical or similar segments are in the successive columns, which are helpful to find the functional or evolutionary structure. We denote the gap between amino acids as  a symbol $'-'$ and the set of elements consisting an aligned amino acid sequence  are replaced by \emph{$A=\{-\}\cup aa$}. Given a family of amino acid sequences, denoted by $\mathbf{F}$, where $\mathbf{F}$ represents the familay of amino acid sequences, the length after alignment is denoted by $l$. For an amino acid sequence, we use \emph{k} to indicate the $k_{th}$ site of the sequence, $c_k$ to represent the amino acid at the $k_{th}$ site, and $c_k^a$ for the amino acid of the specific amino acid sequence \emph{a} at site \emph{k}.

In this paper, We define the event $(k,a)$ as  amino acid \emph{$c_k^a\in A$} being on the $k^{th}, k\in [1..l]$ site of an aligned sequence $\emph{a}\in \mathbf{F}$, and introduce the metric of \emph{capability score} $s(k,a)$ to evaluate how this event influences the optimal condition. Taking alkali resistance as an example, if the \emph{capability score} is higher when amino acid $c_k$ appears at site \emph{k}, it indicates that the occurrence of this event leads to an increase in the optimal pH of the sequence. 


We also consider the correlations between the optimal condition and multiple amino acids appearing on different sites, namely how multiple elements together affect the optimal condition. We define the event $(n_i^a,n_j^a)$ as $c_i^a$ appears on  site  $i$ and $c_j^a$ appears on site $j$ in the sequence $a\in \mathbf{F}$. We also introduce the \emph{suitability score} $s(n_i^a,n_j^a)$ to denote how much this event induces a better optimal condition.
 Different with the single site analysis, the site $k$ and the amino acid $c_k^a$ on it are considered together by concatenating their vectors. 
%\begin{equation}
  %s(n_i^a,n_j^a)=\bf n_i^a\cdot n_j^a=(v_i\oplus c_i^a )\cdot (v_j\oplus c_j^a)
  %\label{score_double}
%\end{equation}
A higher score means that the occurrence of event $(n_i^a,n_j^a)$ induces a better optimal condition and the probability is formalized by  equation \ref{score_probability_double}. 

To combine these two objectives together, a hyper parameter $\alpha\in(0..1)$ is introduced to balance them in the embedding learning process. 
Let $\mathcal{L}_T(\theta)$ and $\cdot\mathcal{L}_C(\theta)$ denote the objective functions for above two scores in the embedding learning process, respectively.  The whole objective function is then defined as 
$\mathcal{L}(\theta)=\alpha\cdot\mathcal{L}_T(\theta)+(1-\alpha)\cdot\mathcal{L}_C(\theta)$.  A larger $\alpha$ induces a higher bias on $s(k,a)$.
Conversely, a smaller $\alpha$ considers more on the correlations between the optimal condition and multiple amino acids appearing on different sites.  
 We would discuss how to learn the embeddings in the Method section. 
 

\subsection*{Dataset}
We crawl the amino acid sequences of the GH11 family on the CAZy website\footnote{The CAZy database describes the families of structurally-related catalytic and carbohydrate-binding modules (or functional domains) of enzymes that degrade, modify, or create glycosidic bonds. http://www.cazy.org/}. In order to extract the required optimal pH of each sequence, we investigate the papers related to the proteins in the GH11 family on the $Web\ of\ Science$ website. There are 272 amino acid sequences in the GH11 family and the length of each sequence is within 128 and 335. After alignment, the length of amino acid sequences is $l$=380.

Considering the conserved and non-conserved segments in the aligned sequences, we quantify the importance of each site against the optimal pH, namely which sites the amino acids being on highly influence the alkali resistance.  We adopt the information gain as a metric. Give a sequence $x$, the information entropy $H(y)$ quantify the uncertainty on whether $x's$ optimal pH $y$ is either higher or lower than the average, i.e. $H(y)=-\sum_{i\in\{l,h\}}{p(x\in S_{i})\cdot \log\! p(x\in S_{i})}$, where  $S_l\subset F$ is the collection of sequences whose optimal pH are lower than the average, $S_h\subset F$ is the collection of sequences whose optimal pH are higher than the average, and $p(x\in S_{i})=\frac{|S_i|}{|\mathbf{F}|}$. Then we quantify how much a concrete site reduces this uncertainty. For site $k$, the conditional entropy $H(y|c_k)$ quantities the uncertainty on the optimal pH after we know the amino acid $c_k$, i.e $H(y|c_k)=-\sum_{i\in\{l,h\}}{p(x\in S_{i}|c_k)\cdot \log\! p(x\in S_{i}|c_k)}$. The conditional entropy on site $k$ is the expectation on different amino acids, i.e. $H(y|k)=\sum_{c_k\in A}\frac{|F_{c_k}|}{|\mathbf{F}|}\cdot H(y|c_k)$, where $F_{c_k}$ denotes the set of sequences with the amino acid $c_k$ on site $k$. The information gain on  site $k$ is $Gain(k)=H(y)-H(y|k)$. The higher this score, the more probability different amino acids on site $k$ influences the optimal condition. 
The results on GH11 family are scatter plotted in Figure \ref{fig:IG}, where the X-axis denotes sites of aligned amino acid sequences and Y-axis denotes the information gain on a site. This point would be considered together with a concrete amino acid for mutation, discussed in next section.




\subsection*{Preprocess}\label{sampling}
We choose the higher alkali-resistance as the expected condition. Then the proteins with the optimal pH higher than 7 are classified into $S_h$ and others into $S_l$. We generate samples in the following three steps: 1) Sample a protein $a\in S_h$ and a protein $b\in S_l$. 2) Select the  site $k$ where the amino acids in $a$ and $b$ are different, i.e. $c_k^a \ne c_k^b$. 3) Put $(k,c_k^a)$ into the positive sample set and $(k,c_k^b)$ into the negative sample set. When learning the embeddings, we  select one or more samples from the two sets for each time.

%\subsection{An approximation sample method}
Since the size of a family dataset  is small, we adopt the approximation method to generate more samples by the statistics on amino acids  in the training dataset. To understand how this method works, we compare the probability distribution on the training samples and all samples by the KL distance, an often adopted metric to evaluate two probability distributions \cite{zhang2015model},
\begin{equation}
  \emph{KL-D}(P_i\parallel Q_i)=\sum_{c\in A}P_i(c)\cdot \log{\frac{P_i(c)}{Q_i(c)}}
  \label{KL}
\end{equation}
where $P_i$ and $Q_i$ denote the distribution on the amino acids at site \emph{i} of an aligned sequence in two datasets, respectively. The smaller, the better.

We set the proportion 50\%, 60\%, 70\%, 80\%, and 90\%, respectively, to randomly select samples as the training part. The comparison results against the proportion are listed in Figure \ref{fig:KL}, where  X-axis represents the site and the Y-axis represents the KL distance. We can see that the difference decreases with an increasing size of training set. We would further analyze how this difference influences the computationally biological results in next section.


\section*{Results}\label{sec:experiments}
In this section, we discuss how to use the learnt embeddings for biological purposes so as to help biologists for efficient analysis without wet experiments.
%
We consider the important task on  an enzyme mutant for an expected optimal condition, which is better than others in a family. A usual way is first to verify the optimal condition of each enzyme in a family by wet experiments and select the one most approaching the expected  condition. Then they perform some mutation on the non-conserved segments and wet experiments to verify the activity of the mutant under the expected  environment. These biological purposes can be transferred as the following two computational problems: the pH predication and mutation suggestion. %, which would be discussed in this section and the next section.

\subsection*{Sequence based optimal pH predication}\label{compare}
 Given two amino acid sequences in a  family, we  want to determine by computation  which one has a better optimal condition, such as a higher temperature resistance. This  result can be used to select an expected enzyme from  a family by comparing  each pair of them.


For two new amino acid sequences $a$ and $b$, we calculate the \emph{quantitative scores} of them, which is defined as the combination of  \emph{capability score} and \emph{suitability score}, i.e. $Score(a)=\alpha\cdot\sum_{k=1}^lS(k,a_k)+(1-\alpha)\cdot \frac {\sum_{i,j=1}^lS(n_i^a,n_j^a)}{l}$, and would be also discussed in the Method Section.  
If $Score(a)>Score(b)$, we judge that $a$ is more alkali-resistance than $b$ or $a$ has a higher activity than $b$ in the same alkaline environment, and vice versa. Since the sites have different effects on the alkali-resistance of proteins,  we select the site with the information gain larger than a threshold $\omega$, which would be discussed in details in the experiments.

In addition to comparing the optimal pH between different amino acid sequences from the same family, we also use the embeddings to predict the optimal pH of amino acid sequences. We try different prediction models, and finally choose the Support Vector Regression(SVR) \cite{chang2011libsvm} method, where the kernel function is \emph{rbf}, punishment term is 10, and the parameter gamma of \emph{rbf} is 0.0001, as the prediction model. The input of the SVR model is the feature vector of the amino acid sequence and the output is the optimal pH.

\subsection*{Comparison methods and evaluation metric}
We use the accuracy as the evaluation criteria. For two amino acid sequences selected from the test set, we compare their optimal pH.
The accuracy is computed by the ratio of correct judgments against the number of independent experiments.

Since the current related works of predicting the optimal condition are based on the tertiary structure and/or free energy, there is not any closely related works on the sequence based predication. Therefore, we select several machine leaning algorithms as the comparison methods.

\begin{itemize}
  \item SVR \cite{chang2011libsvm}. This is the traditional predication method. We try several kernel functions, such as $rbf$, $poly$, $linear$ and $sigmoid$, respectively, and adjust the parameters to verify the effects of the SVR. We select the best results from several variants. The corresponding parameters include: $poly$ as the kernel function, 10 as the punishment term C, 3 as the dimensions of poly functions.
   \item Neural Network(NN) \cite{goodfellow2016deep}. There are many factors affect the results of the Neural Network method, such as the activation function, optimization function, the number of hidden layer and the number of nodes on the neural network model.We explore the effects on different settings and select the best one as: the number of nodes in hidden layer to 1000 in one hidden layer, \emph{tanh} as the activation function and $adam$ as the optimization function.
\end{itemize}

As a comparison to embeddings, we use the one-hot vector on the set $aa$ of amino acids as the sample features to feed these methods.
Namely, each site $k$ of a sequence maps to a vector $x\in \{0,1\}^{|aa|}$ and $x_i=1$ if and only if the amino acid $c_k$ is at rank $i$ in the lexicographical order of the set of amino acids. The gap maps to zero in a vector. For a given sequence $a\in F$,  we predict the value of optimal pH $ \hat{y}_a $ and calculate the residual by $\hat{y}_a-y_a$ with its real optimal pH $y_a$, which are shown in Figure \ref{fig:svrnn}.


From the results, we can see that the residuals in embedding model obeys normal distribution, and the variance is the smallest. Our method is better than the comparison methods in predicting the optimal pH.

Then we compare our method with other methods and show the results in Figure \ref{accuracy}, together with several settings on the ratio between the training set and the test set, i.e. 5:5, 6:4, 7:3, 8:2 and 9:1. The threshold of the information gain $\omega$ is set 0.0478, The dimension of embeddings is set 30 and the preference of the combinative objectives $\alpha$ is set 0.4. We use the 10-fold cross-validation method to verify the accuracy.

Overall, the embedding method gives the best results in most cases comparing with other methods. When the ratio of training set to test set is 8:2, the accuracy is the highest by our method. Comparatively, NN is better than SVR, especially in predicting the optimal pH value for mutants which would be discussed in next subsection.

Considering the ratio of the training set and the test set,  a larger ratio often induces a higher accuracy. This is because when the ratio is small, the training set can not reflect the probability distribution of the whole dataset. But a dominate ratio, i.e. 90\%, may induce the number of test set too small such that the results are with more randomness and can not reflect representative results.

\subsection*{Experiment on the influence of parameters}
In this subsection, we discuss the influence by the super parameter, including the dimensions of embeddings $d$, the preference of the combinative objectives $\alpha$ for embedding learning and the threshold of the information gain $\omega$ for selecting the sites on computation of score.

First, we discuss the influence of different dimensions of embeddings on the accuracy. We try several settings on the dimension of embeddings, i.e. d=5, 10, 30, 50, 90 and 130, respectively and list the result in Table \ref{tab:accuracy}

\begin{table*}[h!]
  \centering
  \fontsize{7}{10}\selectfont
  \begin{threeparttable}
  \caption{Summary on the accuracy against different ratios on training data, the dimensions of embedding, and comparison methods. }
  \label{tab:accuracy}
    \begin{tabular}{c cccccc c c}
    \toprule
    \multirow{2}{*}{train:test}&
    \multicolumn{6}{c}{ Dimension of embedding vector }\cr
    \cmidrule(lr){2-7}

    &5&10&30&50&90&130 \cr
    \midrule
    \hline
    5:5 &0.734&0.753&0.706&0.704&0.694&0.685 \cr
    6:4 &0.755&0.752&0.744&0.751&0.749&0.746\cr
    7:3 &0.743&0.77&0.787&0.792&0.794&0.794 \cr
    8:2 &0.74&0.779&\textbf{0.807}&0.800&0.789&0.787\cr
    9:1 &0.745&0.773&0.787&0.794&0.799&0.799\cr
    \bottomrule
    \end{tabular}
    \end{threeparttable}
\end{table*}

Considering the influence by the dimension of embeddings, we verify the accuracy against different settings. We can see that after the dimension is larger than a threshold 30, there is not obvious difference on accuracy. This illustrates that a small dimension cannot hold  much implicit information on amino acids and the sites of an aligned sequence. However, a too large dimension may induce the sparse problem and is difficult for convergence due to a small quantity of samples.

The embedding dimension in following experiments is set 30, and the ratio between training set to test set is 8:2. We choose the threshold $\omega$ as  0.111, 0.0478, 0.026, 0.0105 and 0.002, and $\alpha$ as 0, 0.1, 0.2, 0.3,..., 1, respectively.

The accuracy results against different parameters, i.e. the balance factor $\alpha$ and the threshold $\omega$ on information gain, are shown in Figure \ref{fig:alpha}, where the X-axis represents the values of $\alpha$, the Y-axis represents the accuracy of our method. It achieves the highest accuracy at $\omega=0.0478$ and $\alpha=0.4$.

Considering the threshold of the information gain $\omega$ for selecting the sites on computing the score, a suitable $\omega$ is required. The best setting is at $\omega=0.0478$, colored red. A large  threshold on the information gain $\omega $ induces fewer sites being selected and some important non-conserved sites being ignored. Conversely a small $\omega$ may introduce more conserved sites in calculating the score as Eq.\ref{score} such that it reduces the accuracy.

Considering the preference on $\alpha$, both the type of amino acids on single site and interactions between multiple sites contribute to accuracy.
When $\alpha$ is at 0.4, i.e. the type of amino acids on single site, our method achieves the best. When $\alpha$ is less than 0.4, our optimization objectives would consider more interaction of different sites. The special case $\alpha=0$ learnt less about the semantics on amino acids, which only considers the correlations between the optimal condition and multiple amino acids appearing on different sites. However a larger $\alpha>0.4$ does not influence much on predication.

\subsection*{Bioinformatics verification}
In this subsection, we verify the predication results with the biological wet experiments.   We take the enzyme $Xylanase\ A$ from the GH11 family  as an example, whose optimal pH is 6.0. Ruller, Alponti and Deliberto et al. generate 5 mutants on $Xylanase\ A$ which can survive in alkaline environment by multiple mutation experiments and determine their activity at 5.5 pH environment by wet experiments \cite{ruller2014concommitant}. We use Eq.\ref{sum_score} to calculate the activity scores of $Xylanase\ A$ and the 5 mutants. The higher the score, the higher the activity of the enzyme at the same pH environment.

\begin{table}[h!]
  \centering
  \fontsize{6.3}{10}\selectfont
  \begin{threeparttable}
  \caption{The activity at 5.5pH and the activity score on $Xylanase\ A$ and Mutant.}
  \label{tab:mutant}
    \begin{tabular}{c|c|c|c}
    \toprule
     Wild and Mutants& Mutation Sites& Activity at 5.5 pH& Score\cr
    \midrule
    \hline
    Wild Xylanase A & N/A &2.27U/ml &4022.35\cr
    Mutant1 &S22P &\textbf{3.58U/ml} $\uparrow$ &\textbf{4024.73} $\uparrow$ \cr
    Mutant2 &G13R &1.89U/ml $\downarrow$ &4018.68 $\downarrow$ \cr
    Mutant3 &Q7H/G13R/S22P &\textbf{2.45U/ml} $\uparrow$ &\textbf{4022.88}$\uparrow$ \cr
    Mutant4 &S22P/H156L/S179C &\textbf{3.01U/ml } $\uparrow$ &\textbf{4024.73}$ \uparrow $\cr
    Mutant5 &S31Y &2.20U/ml $\downarrow$ &4022.35\cr
    \bottomrule
    \end{tabular}
    \end{threeparttable}
\end{table}

Both the activity at 5.5 pH \cite{ruller2014concommitant} and the scores of these enzymes calculated by our method are together listed in Table \ref{tab:mutant}. By comparing the changes in activity and scores between mutants and wild-type $Xylanase\ A$, we can see that the trend of activity are consistent with the scores. These are positive correlations.

A special case is mutant5, which may be caused by  the mutation on site 31. Since
 the information gain on this site is only 0.039 and is lower than the threshold 0.0478, it is not included in calculating the score. But considering the little change in activity after mutation,  it is not important to select this mutant for an expected optimal condition. % also proves that the site 31 is not important for changing optimal pH.


\subsection*{Embeddings based mutation suggestion}\label{sec:mutation}
%\subsection*{Mutation Suggestion} 
In this subsection, we discuss how to guide mutations on enzymes so as to help biologists wet experiments. %For a given wild-type enzyme, we select one or more sites for generating a candidate mutate by embeddings so as to have the optimal pH of mutant higher than wild-type enzyme.
%
 Given a family of enzymes and an expected optimal condition like alkali-resistance, we design a new amino acid sequence such that it has a better optimal condition than the wild-type enzymes. %This can help biological researchers to mutate the enzymes for an expected direction, which reduces too much wet experiments.


For a given wild-type enzyme, we propose two rules for a single-site mutation: 1) Do not mutate the site that is filled by gap $'-'$ so as to remain the original length of the amino acid sequence. 2) Do not mutate the sites with information gain smaller than 0.0478.
Then we calculate the scores of different single-site mutations and select the mutant with the highest score. For multiple-sites mutation, we combine single-site mutations to get multi-sites mutants on a wild-type enzyme, and retain the mutant with the highest score. %For multi-sites mutation, we mutate any single-site of multi-sites for wild-type enzyme will make the optimal pH higher than the wild-type enzyme.

We perform single-site and multi-sites mutations on the known wild-type amino acid sequences in GH11 family. From a biological point of view, the more mutation sites, the greater the uncertainty of the mutation and the higher the cost of mutation. So we only do at most 3-sites mutation for amino acid sequences. We use Neural Network model to predict the optimal pH of mutants, and calculate the change of optimal pH after mutation. The results are shown in Figure \ref{deta}, where the X-axis represents the change of optimal pH and the Y-axis represents the number of amino acid sequence under the change.


For both mutations, the optimal pH of the 88\% amino acid sequence increases. From the prediction results, both mutation suggestions can increase the optimal pH of amino acid sequences and the effect of multi-sites mutation suggestion is better than that of single-site mutation suggestion. Specially, The optimal pH of the 5\% amino acid sequence in single-site mutation and 10\% amino acid sequence in multi-sites mutation increase more than 1.0.


\section*{Discussion}\label{SEC:Semantic analysis on embeddings}
In this proposed method, the learnt embeddings are expected to include the biological semantics for an expected optimal condition. Since the embedding vectors are in the latent space, we adopt the Distributed Stochastic Neighbor Embedding (\emph{t-SNE}) method \cite{maaten2008visualizing} to reduce the high-dimension vectors to 2d coordinates for illustration.

 We first verify whether the capability of amino acids with the sites of a sequence is embedded in the vectors. The learned \emph{capability scores} $s(k,a)$  are illustrated as the thermodynamic chart in Figure \ref{fig:weight_color}, where the X-axis denotes the sites and Y-axis denote the amino acids.  Red color maps to high scores, while green color maps to low scores. The obviously green areas 1-20, 157-166, 215-238, 250-254 and 366-380 indicate amino acids on these sites  influence less the optimal pH, while the red areas indicate that different amino acids on these sites highly influence the optimal pH.
Compared with the results by biologist analysis, these areas are mapped to the conserved segments and  the potential non-conserved segments, which is the basic task for analyzing the enzymes in a family. This convinces the semantics of embeddings.


Then we verify the biological differences between amino acids.  We select Top10 information gain sites and two alkali-resistance amino acids \emph{Lys} and \emph{Arg}, as well as two acid-resistances amino acids \emph{Asp} and \emph{Glu}. We combine the embeddings of each pair of them as Eq.\ref{add} and color them as points in Figure \ref{fig:t-SNE-a} by $t-SNE$. The results show that the points that are combined with the same amino acid are clustered together rather than others. Moreover, the points combined with the  alkali-resistance amino acids, colored green and pink are clustered more closely than with the acid-resistance amino acid, colored blue and red.

At last we compare the sites with different information gain, the top2 sites 278, 80 and the bottom2 sites 1 and 380. The combinative embeddings by sites and amino acids, as Eq.\ref{add}, are colored differently in Figure \ref{fig:t-SNE-b}, where they are obviously clustered into two groups: sites with high information gain and sites with low information gain. We zoom in the details on site 80 and 278, and the results in Figure \ref{fig:t-SNE-c} show that the pink and purple points co-exist in pairs. The points related to alkali-resistance amino acids are clustered together, while the points related acid-resistance amino acids are clustered together. Comparatively, the points on alkali-resistance amino acids are closer than the acid-resistance amino acids. This illustrates that the embeddings have learned much information on distinguishing the amino acids and their influences on alkali-resistance functions.


As for future work, we plan to analyze the optimal conditions of proteins in different families. Since they have different protein structures and aligned sequences, it would be more complex and challenging to analyze their conserved and non-conserved sequences and learning the embeddings of sites. Another interesting direction is to analyze multiple optimal conditions together. For different optimal conditions, the effectiveness of sites and amino acids in a protein may be different. The multi-objective analysis would help us find the relationships between different optimal conditions. This requires an elaborate design on the embeddings method so that they can contain more information to distinguish the features of different optimal conditions.


\section*{Conclusion} \label{sec:conclusion}
In this paper, we proposed an embedding method to represent the amino acids and the construct information as vectors in the latent space. These vectors contain information about the correlations between amino acids and sites in aligned amino acid sequences. Based on these embeddings, we then predict the optimal pH of a new amino acid sequence in the same family and design a method to suggest the site and direction of a mutation for an expected condition by embeddings.
 We adopt the amino acid sequences in glycoside hydrolase GH11 family for the verification of our method. And we design two computational experiments and verify the results by wet experiments, which are to predict the optimal pH of amino acid sequences and to give a mutation suggestion. Compared with the traditional method, this method does not require the tertiary structure of a protein or the situation of free energy of the amino acid sequence, which is more efficient and effective.

These advantages attribute to two aspects. One is that we take into account both the correlations between the amino acid at a single site and the interaction between multiple-sites. Another is that the vectors we have learnt reveal the information about the optimal condition implied in amino acids sequence.


\section*{Method}\label{sec:method}
In this section, we present the details on the embedding method and the learning process.  
\subsection*{The embedding objectives}
The optimal condition of an amino acid sequence is determined by which amino acids consist the enzyme and how they construct together in the sequence. Motivated by this point, we introduce two objective functions to learn the correlations between the optimal condition and the sequence:  1) the influence by one amino acid on each site. 2) the mutual influences by two amino acids on different sites.
We propose an embedding method to combine these two objectives together by modeling the sites and the amino acids of a sequence into the vectors in the same latent space, denoted by ${\bf v},{\bf c}\in R^d$, respectively.


The first  objective learns  the correlations between the optimal condition and the type of amino acid on each site.
We define the event $(k,a)$ as  amino acid \emph{$c_k^a\in A$} being on the $k^{th}, k\in [1..l]$ site of an aligned sequence $\emph{a}\in \mathbf{F}$, and introduce the metric of \emph{capability score} $s(k,a)$ to evaluate how this event influences the optimal condition, formalized by:
\begin{center}
\begin{equation}
    s(k,a)=\bf{v_k \cdot c_k^a}
\label{score}
\end{equation}
\end{center}

where $\textbf{v}_k,\textbf{c}_k^a\in R^d$ are the embeddings of site $k$ and the amino acid $c_k^a$, respectively. A higher score means $c_k^a$ on site $k$  resulting the sequence a better activity under a given condition.



Let \textbf{E} be the set of  all occurrences of different events, namely the combination of amino acids and sites. We adopt the \emph{softmax} function to model the probability of such an event.
\begin{equation}
    p_\theta(k,a)=\frac{\exp\!(s(k,a ))}{\sum_{(k',a')\in \textbf{E}} \exp\!(s(k',a'))}
    \label{score_un}
\end{equation}
where $\theta$ denotes the parameters in this model. To maximize the likelihood of the occurrences, the objective loss function is formalized by
\begin{equation}
  \ \mathcal{L}_T(\theta)=-\sum_{(k,a)\in \textbf{E}}\log\!p_{\theta}(k,a)
\label{target}
\end{equation}

Since the normalization part in the denominator in Eq.\ref{score_un} cost much computation, we use the Noise Contrastive Estimation (NCE) method, proposed by Gutmann et al. \cite{Gutmann2010Noise}, to estimate the optimal parameters $\theta^*$. It treats the normalization part as an additional  parameter, denoted by \emph{C}. The Eq.\ref{score_un} is then re-written as:
\begin{equation}
p_\theta(k,a)=\exp\!(s(k,a)+C)
\label{score_probability}
\end{equation}

According to NCE, we add artificially generated noise data to the training set. The parameters in probability density function and normalization part can be estimated by discriminating the original data and noise data. Let $p(D=1|(k,a),\theta)$ denote the probability that the optimal condition gets higher when  $c_k^a\in A$ appears on site \emph{k} in sequence $a \in \mathbf{F}$. Let  $p(D=0|(k,a),\theta)$ denote the probability that the optimal condition gets lower when $c_k^a$ appears on site \emph{k}.
\begin{equation}
\begin{split}
p(D=1|(k,a),\theta) = \frac{p_\theta(k,a)}{p_\theta(k,a)+p_n(k,a)}\\
=\sigma(\log\!p_\theta(k,a)-\log\!p_n(k,a))
\end{split}
\end{equation}

\begin{equation}
\begin{split}
p(D=0|(k,a),\theta) =1-p(D=1|(k,a),\theta)\\
=1-\sigma(\log\!p_\theta(k,a)-\log\!p_n(k,a))
\end{split}
\end{equation}
where $\sigma(x)=\frac{1}{1+\exp\!(-x)}$ is the \emph{sigmoid} function and $p_n(k,a)$ is the artificial noise distribution.

We fit the model by maximizing the expectation of log-posterior probability over the mixture of observed samples and noise samples. The expectation and loss function are formulated by Eq.\ref{equ:pro} and Eq.\ref{equ:loss}:
\begin{equation}\label{equ:pro}
E_{p_\theta}[\log\!p(D=1|(k,a),\theta)]+E_{p_n}[\log\!p(D=0|(k,a),\theta)]
\end{equation}
\begin{equation}\label{equ:loss}
\begin{split}
  \mathcal{L}_T(\theta)=-[\log\!\sigma(\log\!p_\theta(k,a)-\log\!p_n(k,a))\\
  +\log\!(1-\sigma(\log\!p_\theta(k',a')-\log\!p_n(k',a')))]
\end{split}
\end{equation}

%\subsubsection{Embedding objective the correlations between the optimal condition and multiple amino acids appearing on different sites}
The second embedding objective learns the correlations between the optimal condition and multiple amino acids appearing on different sites, namely how multiple elements together affect the optimal condition. %objective determines the non-conserved segments that highly influence  an enzyme's optimal condition.
%
 Different with the single site analysis, the site $k$ and the amino acid $c_k^a$ on it are considered together by concatenating their vectors,
\begin{equation}
    \bf{n_k^a}=\bf{v_k}\oplus \bf{c_k^a}
    \label{add}
\end{equation}
where $\textbf{v}_k,\textbf{c}_k^a\in R^d$ are the embeddings of site $k$ and the amino acid $c_k^a$, respectively, and $\textbf{n}_k^a \in R^{2d}$ is their joint embedding.
%
We define the event $(n_i^a,n_j^a)$ as $c_i^a$ appears on  site  $i$ and $c_j^a$ appears on site $j$ in the sequence $a\in \mathbf{F}$. We also introduce the \emph{suitability score} $s(n_i^a,n_j^a)$ to denote how much this event induces a better optimal condition.
\begin{equation}
  s(n_i^a,n_j^a)=\bf n_i^a\cdot n_j^a=(v_i\oplus c_i^a )\cdot (v_j\oplus c_j^a)
  \label{score_double}
\end{equation}
A higher score means that the occurrence of event $(n_i^a,n_j^a)$ induces a better optimal condition and the probability is formalized by  equation \ref{score_probability_double}.
\begin{equation}
\begin{aligned}
  p_\theta(n_i^a,n_j^a)&= \frac{\exp\!(s(n_i^a,n_j^a))}{ \sum_{(i',j')\in \textbf{E}} \exp\!(s({n_{i'}^a},{n_{j'}^a}))}\\
\label{score_probability_double}
\end{aligned}
\end{equation}


Similar to the discussion of Eq.\ref{score_probability} to Eq.\ref{equ:loss}, the loss function is formulated as follow:
\begin{equation}
\begin{aligned}
&\mathcal{L}_C(\theta)=-[\log\!\sigma(\log\!p_\theta(n_i^a,n_j^a)-\log\!p_n(n_i^a,n_j^a))\\
&+\log\!(1-\sigma(\log\!p_\theta({n_i^a}',{n_j^a}')-\log\!p_n({n_i^a}',{n_j^a}')))]
\end{aligned}
\label{loss_double}
\end{equation}

We model the loss function as a linear combination of  the above two loss functions, where $\alpha\in[0,1]$ is a preference parameter.
\begin{equation}\label{equ:objectives}
\mathcal{L}(\theta)=\alpha\cdot\mathcal{L}_T(\theta)+(1-\alpha)\cdot\mathcal{L}_C(\theta)
\end{equation}

Through the above introduction, we define two score functions $s(k,a)$ and $s(n_i^a,n_j^a)$. Based on two score functions, we introduce a \emph{quantitative score} as a specific optimal condition for a specific sequence $a\in F$. It considers two parts: the compatibility of the amino acid on each site and the influence induced by multiple amino acids on different sites, which are combined by the preference parameter $\alpha$, defined in Eq.\ref{equ:objectives}.
\begin{equation}\label{sum_score}
  Score(a)=\alpha\cdot\sum_{k=1}^lS(k,a_k)+(1-\alpha)\cdot \frac {\sum_{i,j=1}^lS(n_i^a,n_j^a)}{l}
\end{equation}


\subsection*{The embedding learning process}
We adopt the stochastic gradient descent algorithm  to optimize the parameters, proposed by \emph{Adam} \cite{Kingma2014Adam}. Given the set of amino acid sequences in a family, we partition them into two sets: one for training and another for test. For the expected optimal condition, the training set is further classified into the high optimal condition set $S_h$ and the low optimal condition set $S_l$. For example, if a higher pH is expected as the optimal condition, the amino acid sequences in $S_h$ have higher pH than those in $S_l$. We select a sequence $a$ from $S_h$ and the amino acid $c_k^a$ on site $i$ of $a$  as a positive sample, i.e. the pair $(i,a)$. A negative sample $(j,b)$ is similarly chosen by site $j$ on $b\in S_l$ so as to calculate the gradient of parameters. Taking Eq.\ref{equ:loss} as an example, the gradient of the embedding $\bf{v}_i$ of site \emph{i} in the loss function is calculated by:
\begin{equation}
\begin{split}
&\frac{\partial\mathcal{L}_{T}}{\partial \bf{v}_i}=-[\sigma(\log\!p_\theta(i,a)-\log\!p_n(i,a))-1]\bf{c}_i^a
\end{split}
\label{gradient}
\end{equation}


To calculate Eq.\ref{loss_double}, we need two positive samples $(i,a)$, $(j,a)$ to compose event $(n_i^a,n_j^a)$, and two negative samples $(i',b)$, $(j',b)$ to compose event $(n_{i'}^b,n_{j'}^b)$.
\begin{equation}
    \begin{split}
        \frac{\partial\mathcal{L}_{C}}{\partial \bf{v}_{i}}=-[\sigma(\log\!p_\theta(n_i^a,n_j^a)-\log\!p_n(n_i^a,n_j^a))-1]\bf n_j^a
    \end{split}
\label{gradient}
\end{equation}

The learning process would stop until the objective function is converged. %Due to the limitation of space, we omit the equations on parameter updates one by one.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Backmatter begins here                   %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{backmatter}

\section*{Abbreviations}
GH11: Glycoside Hydrolase Family 11; CAZy: Carbohydrate-Active enzymes; NCE: Noise Contrastive Estimation; KL Distance: Kullback Leibler Distance; t-SNE: Distributed Stochastic Neighbor Embedding; SVR: Support Vector Regression; NN: Neural Network
\section*{Ethics approval and consent to participate}
Not applicable
\section*{Consent for publication}
Not applicable
\section*{Availability of data and material}
All data used in the experiments are from CAZy database.
\section*{Competing interests}
  The authors declare that they have no competing interests.
\section*{Acknowledgements}
The scientific calculations in this paper have been done on the HPC Cloud Platform of Shandong University.
This work was supported by the National Key R\&D Program of China (2018YFB0204000,2018YFC0831401), the Major Project of NSF Shandong Province (ZR2018ZB0420), and National Natural Science Foundation of China(91646119,31770054), and the Key Research and Development Program of Shandong Province (2017GGX10114).

\section*{Author's contributions}
YS is responsible for model building, paper structure editing, and proofreading. LW is responsible for the work guidance and planning of the biological direction, explaining the biological significance of the experiment. BG is responsible for the experimental code structure planning and proofreading work. XL is responsible for experimental work, data analysis, paper writing, and proofreading. ZD prepared raw data, wrote and edited the manuscript, and interpreted biological meaning of experimental results. All authors read and approved the final manuscript.
\section*{Funding}
This work was supported by the National Key R\&D Program of China (2018YFB0204000,2018YFC0831401), the Major Project of NSF Shandong Province (ZR2018ZB0420), and National Natural Science Foundation of China(91646119,31770054), and the Key Research and Development Program of Shandong Province (2017GGX10114).

%This work was supported by the National Key R\&D Program of China (2018YFC0831401), the National Natural Science Foundation of China (91646119), the Major Project of NSF Shandong Province (ZR2018ZB0420), the National Major Science and Technology Projects of China (91646119), and National Natural Science Foundation of China(31770054), and the Key Research and Development Program of Shandong Province (2017GGX10114).


The funding body played no role in the design of the study, the collection, analysis, and interpretation of the data or in writing of the manuscript.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                  The Bibliography                       %%
%%                                                         %%
%%  Bmc_mathpys.bst  will be used to                       %%
%%  create a .BBL file for submission.                     %%
%%  After submission of the .TEX file,                     %%
%%  you will be prompted to submit your .BBL file.         %%
%%                                                         %%
%%                                                         %%
%%  Note that the displayed Bibliography will not          %%
%%  necessarily be rendered by Latex exactly as specified  %%
%%  in the online Instructions for Authors.                %%
%%                                                         %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% if your bibliography is in bibtex format, use those commands:
\bibliographystyle{bmc-mathphys} % Style BST file (bmc-mathphys, vancouver, spbasic).
\bibliography{bmc_article}      % Bibliography file (usually '*.bib' )

\section*{Figure Legends}

\begin{figure}[h!]
  \caption{\csentence{Information gain for the sites in aligned amino acid sequences in GH11 family.}}
  \label{fig:IG}
\end{figure}

\begin{figure}[h!]
\caption{\csentence{Evaluation on the approximate method}}
\label{fig:KL}
\end{figure}

\begin{figure}[h!]
\center
\caption{\csentence{The thermodynamic chart for the capability scores on PH for the GH11 family.} The X-axis represents the aligned 380 sites of sequences and the Y-axis represents amino acids and the gap in set A. }
\label{fig:weight_color}
\end{figure}

\begin{figure}[h!]
\center
\caption{\csentence{The $t-SNE$ coordinates of the embeddings.}Each point denotes a combinative embedding by a site and an amino acid. The selected amino acids  \emph{Ala}, \emph{Glu}, \emph{Arg}, \emph{Asp} and \emph{Lys} are colored by orange, red, green, blue and pink, respectively.}
\label{fig:t-SNE-a}
\end{figure}

\begin{figure}[h!]
\center
\caption{\csentence{The $t-SNE$ coordinates of the embeddings.}The sites 1, 80, 278 and 380 are colored brown, pink, purple and green, respectively.}
\label{fig:t-SNE-b}
\end{figure}

\begin{figure}[h!]
\center
\caption{\csentence{The $t-SNE$ coordinates of the embeddings.}The sites 80 and 278 with different amino acids are labeled pink and purple, respectively.}
\label{fig:t-SNE-c}
\end{figure}

\begin{figure}[h!]
\caption{\csentence{Probability distribution on the predication residuals $\hat{y} - y$.} The $\hat{y}$ is the predicted value of optimal pH for a given sequence, and the $y$ is the real optimal pH.}
\label{fig:svrnn}
\end{figure}

\begin{figure}[h!]
    \caption{\csentence{The accuracy of the embedding method and comparison methods in different ratios on training data.}}
\label{accuracy}
\end{figure}

\begin{figure}[h!]
 \caption{\csentence{The accuracy against different parameter settings.} Each curve represents the accuracy under the given threshold $\omega$.}
 \label{fig:alpha}
\end{figure}

\begin{figure}[h!]
    \centering
    \caption{\csentence{Statistical of the change of optimal pH.} Each cylinder represents that how many mutants change their optimal pH. X-axis represents the change in optimal pH. }
\label{deta}
\end{figure}

\end{backmatter}
\end{document}
