\documentclass[shortnames,nojss,notitle,article]{jss}
\usepackage{booktabs}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{tablefootnote}
%\VignetteIndexEntry{protr: R package for generating various numerical representation schemes of protein sequence}
%\VignetteKeywords{R, protr, protein sequence, amino acid, descriptor calculation, feature extraction, parallel computation, structural similarity, functional similarity, sequence alignment, Gene Ontology, GO}
%\VignettePackage{protr}

% header and footer
\usepackage{fancyhdr}
\pagestyle{fancy}
\fancyhead[LE]{The \pkg{protr} Manual}
\fancyhead[RO]{The \pkg{protr} Manual}
\fancyhead[RE]{}
\fancyhead[LO]{}
\cfoot{\thepage}

\author{Nan Xiao\\Central South University \And
        Qing-Song Xu\\Central South University \And
        Dong-Sheng Cao\\Central South University}
\Plainauthor{Nan Xiao, Qing-Song Xu, Dong-Sheng Cao}

\title{\pkg{protr}: \proglang{R} package for generating various numerical representation schemes of protein sequence}
\Plaintitle{protr: R package for generating various numerical representation schemes of protein sequence}

\Volume{xx}
\Issue{x}
\Month{Xxxxx}
\Year{2014}
\Submitdate{2014-xx-xx}
\Acceptdate{2014-xx-xx}

\Address{
  Nan Xiao \\
  School of Mathematics and Statistics\\
  Central South University\\
  Changsha, Hunan, China\\
  E-mail: \email{me@nanx.me}\\
  URL: \url{http://nanx.me}\\

  Qing-Song Xu\\
  School of Mathematics and Statistics\\
  Central South University\\
  Changsha, Hunan, China

  Dong-Sheng Cao\\
  School of Pharmaceutical Sciences\\
  Central South University\\
  Changsha, Hunan, China

}

%
<<prelim,echo=FALSE,print=FALSE>>=
protr.version = '1.2-0'
now.date = strftime(Sys.Date(), "%Y-%m-%d")
@
%

\begin{document}
\SweaveOpts{concordance=TRUE}

    \begin{center}
    \vspace*{\baselineskip}
    \rule{\textwidth}{1.6pt}\vspace*{-\baselineskip}\vspace*{2pt}
    \rule{\textwidth}{0.4pt}\\[2\baselineskip]
    {\LARGE \pkg{protr}\textbf{: \proglang{R} package for generating various numerical representation schemes of protein sequence}}\\[1.2\baselineskip]
    \rule{\textwidth}{0.4pt}\vspace*{-\baselineskip}\vspace{3.2pt}
    \rule{\textwidth}{1.6pt}\\[2\baselineskip]
    {\Large Nan Xiao, Qing-Song Xu, Dong-Sheng Cao}\\[\baselineskip]
    {\large Package Version: \Sexpr{protr.version}}\par
    {\large \Sexpr{now.date}}\par
    \vfill
    \begin{figure}[h]
    \centering
    \includegraphics[width=0.6\textwidth]{fig/logo-panel-text.pdf}
    \end{figure}

%    {\scshape Computational Biology and Drug Design Group} \\
%    {\scshape Central South University, China}\par
    \end{center}

\thispagestyle{empty}

\clearpage

\tableofcontents

\thispagestyle{empty}

\clearpage

\setcounter{page}{1}

\begin{center}
\textbf{Abstract}
\end{center}

The \pkg{protr} package offers a unique and comprehensive
toolkit for generating various numerical representation schemes of protein
sequence. The descriptors included are extensively
utilized in Bioinformatics and Chemogenomics research. The commonly used
descriptors listed in \pkg{protr} include amino acid composition,
autocorrelation, CTD, conjoint traid, quasi-sequence order, pseudo amino
acid composition, and profile-based descriptors derived by
Position-Specific Scoring Matrix (PSSM). The descriptors for
proteochemometric (PCM) modeling, includes the scales-based
descriptors derived by principal components analysis, factor analysis,
multidimensional scaling, amino acid properties (AAindex), 20+ classes of
2D and 3D molecular descriptors (Topological, WHIM, VHSE, etc.), and
BLOSUM/PAM matrix-derived descriptors. The \pkg{protr} package also
integrates the function of parallelized similarity computation derived by
pairwise protein sequence alignment and Gene Ontology (GO) semantic
similarity measures. \pkg{ProtrWeb}, the web server built on \pkg{protr}, is located
at: \url{http://protr.org}.

\textbf{Keywords:} protein sequence, amino acid, descriptor, structural similarity, functional similarity, sequence alignment, Gene Ontology

\section{The Full Workflow Using protr}\addcontentsline{toc}{section}{1. The Full Workflow Using protr}

Here we use the subcellular localization dataset of human proteins presented in the study of \cite{chou2008cell} to demonstrate the full workflow when using protr.

The complete dataset includes 3134 protein sequences (2750 different proteins), classified into 14 human subcellular locations. We select two classes of proteins as our benchmark dataset. Class 1 contains 325 \emph{extracell} proteins, and class 2 includes 307 \emph{mitochondrion} proteins.

First, we load the \pkg{protr} package, then read the protein sequences stored in two separated FASTA files with \code{readFASTA()}:

\begin{CodeInput}
require(protr)

# load FASTA files
extracell = readFASTA(system.file('protseq/extracell.fasta',
                                   package = 'protr'))
mitonchon = readFASTA(system.file('protseq/mitochondrion.fasta',
                                   package = 'protr'))
\end{CodeInput}

To read protein sequences stored in PDB format files, use \code{readPDB()} instead. The loaded sequences will be stored as two lists in \proglang{R}, and each component in the list is a character string representing one protein sequence. In this case, there are 325 \emph{extracell} protein sequences and 306 \emph{mitonchon} protein sequences:

\begin{CodeInput}
length(extracell)
\end{CodeInput}

\begin{CodeOutput}
## [1] 325
\end{CodeOutput}

\begin{CodeInput}
length(mitonchon)
\end{CodeInput}

\begin{CodeOutput}
## [1] 306
\end{CodeOutput}

To assure that the protein sequences only have the twenty standard amino acid types which is required for the descriptor computation, we use the \code{protcheck()} function in \pkg{protr} to do the amino acid type sanity checking and remove the \emph{non-standard} sequences:

\begin{CodeInput}
extracell = extracell[(sapply(extracell, protcheck))]
mitonchon = mitonchon[(sapply(mitonchon, protcheck))]
\end{CodeInput}

\begin{CodeInput}
length(extracell)
\end{CodeInput}

\begin{CodeOutput}
## [1] 323
\end{CodeOutput}

\begin{CodeInput}
length(mitonchon)
\end{CodeInput}

\begin{CodeOutput}
## [1] 304
\end{CodeOutput}

Two protein sequences were removed from each class. For the remaining sequences, we calculate the Type II PseAAC descriptor, i.e., the amphiphilic pseudo amino acid composition (APAAC) descriptor \citep{chouapaac} and make class labels for classification modeling.

\begin{CodeInput}
# calculate APAAC descriptors
x1 = t(sapply(extracell, extractAPAAC))
x2 = t(sapply(mitonchon, extractAPAAC))
x  = rbind(x1, x2)

# make class labels
labels = as.factor(c(rep(0, length(extracell)), rep(1, length(mitonchon))))
\end{CodeInput}

In \pkg{protr}, the functions of commonly used descriptors for protein sequences and proteochemometric (PCM) modeling descriptors are named after \code{extract...()}.

Next, we will split the data into a $75\%$ training set and a $25\%$ test set.

\begin{CodeInput}
# split training and test set
set.seed(1001)
tr.idx = c(sample(1:nrow(x1), round(nrow(x1) * 0.75)),
           sample(nrow(x1) + 1:nrow(x2), round(nrow(x2) * 0.75)))
te.idx = setdiff(1:nrow(x), tr.idx)
x.tr   = x[tr.idx, ]
x.te   = x[te.idx, ]
y.tr   = labels[tr.idx]
y.te   = labels[te.idx]
\end{CodeInput}

We will train a random forest classification model on the training set with 5-fold cross-validation, using the \pkg{randomForest} package.

\begin{CodeInput}
require(randomForest)
rf.fit = randomForest(x.tr, y.tr, cv.fold = 5)
print(rf.fit)
\end{CodeInput}

The training result is:

\begin{CodeOutput}
## Call:
##  randomForest(x = x.tr, y = y.tr, cv.fold = 5)
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 8
##
##         OOB estimate of  error rate: 25.11%
## Confusion matrix:
##     0   1 class.error
## 0 196  46   0.1900826
## 1  72 156   0.3157895
\end{CodeOutput}

With the model trained on the training set, we predict on the test set and plot the ROC curve with the \pkg{pROC} package, as is shown in figure \ref{fig:roc}.

\begin{CodeInput}
# predict on test set
rf.pred = predict(rf.fit, newdata = x.te, type = 'prob')[, 1]

# plot ROC curve
require(pROC)
plot.roc(y.te, rf.pred, col = '#0080ff', grid = TRUE, print.auc = TRUE)
\end{CodeInput}

The area under the ROC curve (AUC) is:

\begin{CodeOutput}
## Call:
## plot.roc.default(x = y.te, predictor = rf.pred, col = "#0080ff",
##                  grid = TRUE, print.auc = TRUE)
##
## Data: rf.pred in 81 controls (y.te 0) > 76 cases (y.te 1).
## Area under the curve: 0.8697
\end{CodeOutput}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.6\textwidth]{fig/roc.pdf}
\caption{ROC curve for the test set of protein subcellular localization data} \label{fig:roc}
\end{figure}

\section{Package Overview}\label{sec:intro}\addcontentsline{toc}{section}{2. Package Overview}

The \pkg{protr} package \citep{CRAN:protr} implemented most of the state-of-the-art protein sequence descriptors with \proglang{R}. The \pkg{protr} package is freely available from the Comprehensive \proglang{R} Archive Network (\url{http://CRAN.R-project.org/package=protr}). This vignette corresponds to \pkg{protr} version \Sexpr{protr.version} and was typeset on \Sexpr{now.date}.

Generally, each type of the descriptors (features) could be calculated with a function named \code{extractX()} in the \pkg{protr} package, where \code{X} stands for the abbrevation of the descriptor name. The descriptors and the function names implemented are listed below:

\begin{itemize}
  \item Amino acid composition
  \begin{itemize}
  \item \code{extractAAC()} - Amino acid composition
  \item \code{extractDC()} - Dipeptide composition
  \item \code{extractTC()} - Tripeptide composition
  \end{itemize}
  \item Autocorrelation
  \begin{itemize}
  \item \code{extractMoreauBroto()} - Normalized Moreau-Broto autocorrelation
  \item \code{extractMoran()} - Moran autocorrelation
  \item \code{extractGeary()} - Geary autocorrelation
  \end{itemize}
  \item CTD
  \begin{itemize}
  \item \code{extractCTDC()} - Composition
  \item \code{extractCTDT()} - Transition
  \item \code{extractCTDD()} - Distribution
  \end{itemize}
  \item Conjoint triad descriptors
  \begin{itemize}
  \item \code{extractCTriad()} - Conjoint triad descriptors
  \end{itemize}
  \item Quasi-sequence-order descriptors
  \begin{itemize}
  \item \code{extractSOCN()} - Sequence-order-coupling number
  \item \code{extractQSO()} - Quasi-sequence-order descriptors
  \end{itemize}
  \item Pseudo-amino acid composition
  \begin{itemize}
  \item \code{extractPAAC()} - Pseudo-amino acid composition
  \item \code{extractAPAAC()} - Amphiphilic pseudo-amino acid composition
  \end{itemize}
  \item Profile-based descriptors
  \begin{itemize}
  \item \code{extractPSSM()}
  \item \code{extractPSSMAcc()}
  \item \code{extractPSSMFeature()}
  \end{itemize}
\end{itemize}

The descriptors commonly used in Proteochemometric Modeling (PCM) implemented in \pkg{protr} include:

\begin{itemize}
  \item \code{extractScales()} and \code{extractScalesGap()} - Scales-based descriptors derived by Principal Components Analysis
  \begin{itemize}
  \item \code{extractProtFP()} and \code{extractProtFPGap()} - Scales-based descriptors derived by amino acid properties from AAindex (a.k.a. Protein Fingerprint)
  \item \code{extractDescScales()} - Scales-based descriptors derived by 20+ classes of 2D and 3D molecular descriptors (Topological, WHIM, VHSE, etc.)
  \end{itemize}
  \item \code{extractFAScales()} - Scales-based descriptors derived by Factor Analysis
  \item \code{extractMDSScales()} - Scales-based descriptors derived by Multidimensional Scaling
  \item \code{extractBLOSUM()} - BLOSUM and PAM matrix-derived descriptors
\end{itemize}

The \pkg{protr} package integrates the function of parallelized similarity score computation derived by local or global protein sequence alignment between a list of protein sequences, the sequence alignment computation is provided by \pkg{Biostrings}, the corresponding functions listed in the \pkg{protr} package include:

\begin{itemize}
  \item \code{twoSeqSim()} - Similarity calculation derived by sequence alignment between two protein sequences
  \item \code{parSeqSim()} - Parallelized pairwise similarity calculation with a list of protein sequences
\end{itemize}

The \pkg{protr} package also integrates the function of parallelized similarity score computation derived by Gene Ontology (GO) semantic similarity measures between a list of GO terms / Entrez Gene IDs, the GO similarity computation is provided by \pkg{GOSemSim}, the corresponding functions listed in the \pkg{protr} package include:

\begin{itemize}
  \item \code{twoGOSim()} - Similarity calculation derived by GO-terms semantic similarity measures between two GO terms / Entrez Gene IDs
  \item \code{parGOSim()} - Pairwise similarity calculation with a list of GO terms / Entrez Gene IDs
\end{itemize}

To use the \code{parSeqSim()} function, we suggest the users to install the packages \pkg{foreach} and \pkg{doParallel} first, in order to make the parallelized pairwise similarity computation available.

In the next sections, we'll introduce the descriptors and function usage in this order.

\section{Commonly Used Descriptors}\label{sec:descriptors}\addcontentsline{toc}{section}{3. Commonly Used Descriptors}

\textbf{Disclaimer.} Users of the \pkg{protr} package need to intelligently evaluate the underlying details of the descriptors provided, instead of using protr with their data blindly, especially for the descriptor types with more flexibility. It would be wise for the users to use some negative and positive control comparisons where relevant to help guide interpretation of the results.

A protein or peptide sequence with $N$ amino acid residues could be generally represented as $\{\,R_1, R_2, \ldots, R_n\,\}$, where $R_i$ represents the residue at the $i$-th position in the sequence. The labels $i$ and $j$ are used to index amino acid position in a sequence, and $r$, $s$, $t$ are used to represent the amino acid type. The computed descriptors are roughly divided into 4 groups according to their known applications described in the literature.

A protein sequence could be divided equally into segments and the methods, described as follows for the global sequence, could be applied to each segment.

\subsection{Amino Acid Composition (AAC)}\label{subsec:aac}

The Amino Acid Composition (AAC) is the fraction of each amino acid type within a protein. The fractions of all 20 natural amino acids are calculated as:

$$
f(r) = \frac{N_r}{N} \quad r = 1, 2, \ldots, 20.
$$

where $N_r$ is the number of the amino acid type $r$ and $N$ is the length of the sequence.

As was described above, we could use the function \code{extractAAC()} to extract the descriptors (features) from protein sequences:

%
<<extractAAC>>=
require(protr)
x = readFASTA(system.file('protseq/P00750.fasta', package = 'protr'))[[1]]
extractAAC(x)
@
%

Here with the function \code{readFASTA()} we loaded a single protein sequence (P00750, Tissue-type plasminogen activator) from a FASTA format file. Then extracted the AAC descriptors with \code{extractAAC()}. The result returned is a named vector, whose elements are tagged with the name of each amino acid.

\subsection{Dipeptide Composition (DC)}\label{subsec:dc}

The Dipeptide Composition (DC) gives 400 descriptors, defined as:

$$
f(r, s) = \frac{N_{rs}}{N - 1} \quad r, s = 1, 2, \ldots, 20.
$$

where $N_{rs}$ is the number of dipeptide represented by amino acid type $r$ and type $s$. Similar to \code{extractAAC()}, here we use \code{extractDC()} to compute the descriptors:

%
<<extractDC>>=
dc = extractDC(x)
head(dc, n = 30L)
@
%

Here we only showed the first 30 elements of the result vector and omitted the rest of the result. The element names of the returned vector are self-explanatory as before.

\subsection{Tripeptide Composition (TC)}\label{subsec:tc}

The Tripeptide Composition (TC) gives 8000 descriptors, defined as:

$$
f(r, s, t) = \frac{N_{rst}}{N - 2} \quad r, s, t = 1, 2, \ldots, 20
$$

where $N_{rst}$ is the number of tripeptides represented by amino acid type $r$, $s$ and $t$. With function \code{extractTC()}, we could easily obtain the length-8000 descriptor, to save some space, here we also omitted the tedious outputs:

%
<<extractTC>>=
tc = extractTC(x)
head(tc, n = 36L)
@
%

\subsection{Autocorrelation Descriptors}\label{subsec:autocorr}

Autocorrelation descriptors are defined based on the distribution of amino acid properties along the sequence. The amino acid properties used here are various types of amino acids index (Retrieved from AAindex Database: \url{http://www.genome.jp/dbget/aaindex.html}, see \citet{aaindex1999}, \citet{aaindex2000}, and \citet{aaindex2008}, see Figure \ref{fig:AAindex} for an illustrated example). Three types of autocorrelation descriptors are defined here and described below.

All the amino acid indices are centralized and standardized before the calculation, i.e.

$$
P_r = \frac{P_r - \bar{P}}{\sigma}
$$

where $\bar{P}$ is the average of the property of the 20 amino acids:

$$
\bar{P} = \frac{\sum_{r=1}^{20} P_r}{20} \quad \textrm{and} \quad \sigma = \sqrt{\frac{1}{2} \sum_{r=1}^{20} (P_r - \bar{P})^2}
$$

\begin{figure}[htbp]
\centering
\includegraphics[width=0.8\textwidth]{fig/AAindex.pdf}
\caption{An illustrated example in the AAIndex database} \label{fig:AAindex}
\end{figure}

\textbf{\subsubsection{Normalized Moreau-Broto Autocorrelation Descriptors}}\label{sssec:moreaubroto}

Moreau-Broto autocorrelation descriptors application to protein sequences could be defined as:

$$
AC(d) = \sum_{i=1}^{N-d} P_i P_{i + d} \quad d = 1, 2, \ldots, \textrm{nlag}
$$

where $d$ is called the lag of the autocorrelation and $P_i$ and $P_{i+d}$ are the properties of the amino acids at position $i$ and $i+d$, respectively. $\textrm{nlag}$ is the maximum value of the lag.

The normalized Moreau-Broto autocorrelation descriptors are defined as:

$$
ATS(d) = \frac{AC(d)}{N-d} \quad d = 1, 2, \ldots, \textrm{nlag}
$$

The corresponding function for this descriptor is \code{extractMoreauBroto()}. A typical call could be:

%
<<extractMoreau1>>=
moreau = extractMoreauBroto(x)
head(moreau, n = 36L)
@
%

The 8 default properties used here are:

\begin{itemize}
\item \textbf{AccNo. CIDH920105} --- {Normalized Average Hydrophobicity Scales}
\item \textbf{AccNo. BHAR880101} --- {Average Flexibility Indices}
\item \textbf{AccNo. CHAM820101} --- {Polarizability Parameter}
\item \textbf{AccNo. CHAM820102} --- {Free Energy of Solution in Water, kcal/mole}
\item \textbf{AccNo. CHOC760101} --- {Residue Accessible Surface Area in Tripeptide}
\item \textbf{AccNo. BIGC670101} --- {Residue Volume}
\item \textbf{AccNo. CHAM810101} --- {Steric Parameter}
\item \textbf{AccNo. DAYM780201} --- {Relative Mutability}
\end{itemize}

Users could change the property names of AAindex database with the argument \code{props}.
 The AAindex data shipped with \pkg{protr} could be loaded by \code{data(AAindex)}, which has the detailed information of each property. With the argument \code{customprops} and \code{nlag}, users could specify their own properties and lag value to calculate with. For illustration, we could use:

%
<<extractMoreau2>>=
# Define 3 custom properties
myprops = data.frame(AccNo = c("MyProp1", "MyProp2", "MyProp3"),
                     A = c(0.62,  -0.5, 15),  R = c(-2.53,   3, 101),
                     N = c(-0.78,  0.2, 58),  D = c(-0.9,    3, 59),
                     C = c(0.29,    -1, 47),  E = c(-0.74,   3, 73),
                     Q = c(-0.85,  0.2, 72),  G = c(0.48,    0, 1),
                     H = c(-0.4,  -0.5, 82),  I = c(1.38, -1.8, 57),
                     L = c(1.06,  -1.8, 57),  K = c(-1.5,    3, 73),
                     M = c(0.64,  -1.3, 75),  F = c(1.19, -2.5, 91),
                     P = c(0.12,     0, 42),  S = c(-0.18, 0.3, 31),
                     T = c(-0.05, -0.4, 45),  W = c(0.81, -3.4, 130),
                     Y = c(0.26,  -2.3, 107), V = c(1.08, -1.5, 43))

# Use 4 properties in the AAindex database, and 3 cutomized properties
moreau2 = extractMoreauBroto(x, customprops = myprops,
                             props = c('CIDH920105', 'BHAR880101',
                                       'CHAM820101', 'CHAM820102',
                                       'MyProp1', 'MyProp2', 'MyProp3'))
head(moreau2, n = 36L)
@
%

About the standard input format of \code{props} and \code{customprops}, see \code{?extractMoreauBroto} for details.

\textbf{\subsubsection{Moran Autocorrelation Descriptors}}\label{sssec:moran}

Moran autocorrelation descriptors application to protein sequence may be defined as:

$$
I(d) = \frac{\frac{1}{N-d} \sum_{i=1}^{N-d} (P_i - \bar{P}') (P_{i+d} - \bar{P}')}{\frac{1}{N} \sum_{i=1}^{N} (P_i - \bar{P}')^2} \quad d = 1, 2, \ldots, 30
$$

where $d$ and $P_i$ and $P_{i+d}$ are defined in the same way as in the first place, and $\bar{P}'$ is the considered property $P$ along the sequence, i.e.,

$$
\bar{P}' = \frac{\sum_{i=1}^N P_i}{N}
$$

$d$, $P$, $P_i$ and $P_{i+d}$, $\textrm{nlag}$ have the same meaning as above.

With \code{extractMoran()}, which has exactly the same arguments with \code{extractMoreauBroto()}, we could compute the Moran autocorrelation descriptors (only output the first 36 elements of the result):

%
<<extractMoran>>=
# Use the 3 custom properties defined before
# and 4 properties in the AAindex database
moran = extractMoran(x, customprops = myprops,
                     props = c('CIDH920105', 'BHAR880101',
                               'CHAM820101', 'CHAM820102',
                               'MyProp1', 'MyProp2', 'MyProp3'))
head(moran, n = 36L)
@
%

\textbf{\subsubsection{Geary Autocorrelation Descriptors}}\label{sssec:geary}

Geary autocorrelation descriptors for protein sequence could be defined as:

$$
C(d) = \frac{\frac{1}{2(N-d)} \sum_{i=1}^{N-d} (P_i - P_{i+d})^2}{\frac{1}{N-1} \sum_{i=1}^{N} (P_i - \bar{P}')^2} \quad d = 1, 2, \ldots, 30
$$

where $d$, $P$, $P_i$ and $P_{i+d}$, $\textrm{nlag}$ have the same meaning as above.

For each amino acid index, there will be $3 \times \textrm{nlag}$ autocorrelation descriptors. The usage of \code{extractGeary()} is exactly the same with \code{extractMoreauBroto()} and \code{extractMoran()}:

%
<<extractGeary>>=
# Use the 3 custom properties defined before
# and 4 properties in the AAindex database
geary = extractGeary(x, customprops = myprops,
                     props = c('CIDH920105', 'BHAR880101',
                               'CHAM820101', 'CHAM820102',
                               'MyProp1', 'MyProp2', 'MyProp3'))
head(geary, n = 36L)
@
%

\subsection{Composition / Transition / Distribution}\label{subsec:ctd}

These descriptors are developed and described by \citet{dubchak1} and \citet{dubchak2}.

\begin{figure}[htbp]
\centering
\includegraphics[width=\textwidth]{fig/CTD.pdf}
\caption{The sequence of a hypothetic protein indicating the construction of composition, transition and distribution descriptors of a protein. Sequence index indicates the position of an amino acid in the sequence. The index for each type of amino acids in the sequence (`1', `2' or `3') indicates the position of the first, second, third, ... of that type of amino acid. 1/2 transition indicates the position of `12' or `21' pairs in the sequence (1/3 and 2/3 are defined in the same way.).} \label{fig:ctd}
\end{figure}

\textbf{Step 1: Sequence Encoding}

The amino acids are divided in three classes according to its attribute and each amino acid is encoded by one of the indices 1, 2, 3 according to which class it belonged. The attributes used here include hydrophobicity, normalized van der Waals volume polarity, and polarizability, as in the references. The corresponding division is in the table 1.

\begin{table}[htbp]
  \centering
  \caption{Amino acid attributes and the division of the amino acids into three groups for each attribute}
  \scriptsize{
    \begin{tabular}{cccc}
    \addlinespace
    \toprule
          & \textbf{Group 1} & \textbf{Group 2} & \textbf{Group 3} \\
    \midrule
    Hydrophobicity & Polar & Neutral & Hydrophobicity \\
          & R, K, E, D, Q, N & G, A, S, T, P, H, Y & C, L, V, I, M, F, W \\
          \midrule
    Normalized van der  & 0-2.78  & 2.95-4.0  & 4.03-8.08 \\
    Waals Volume    & G, A, S, T, P, D, C & N, V, E, Q, I, L  & M, H, K, F, R, Y, W \\
    \midrule
    Polarity & 4.9-6.2  & 8.0-9.2  & 10.4-13.0 \\
          & L, I, F, W, C, M, V, Y  & P, A, T, G, S  & H, Q, R, K, N, E, D \\
          \midrule
    Polarizability & 0-1.08  & 0.128-0.186  & 0.219-0.409 \\
          & G, A, S, D, T & C, P, N, V, E, Q, I, L  & K, M, H, F, R, Y, W \\
          \midrule
    Charge & Positive  & Neutral  & Negative \\
          & K, R  & \tiny{A, N, C, Q, G, H, I, L, M, F, P, S, T, W, Y, V}  & D, E \\
          \midrule
    Secondary & Helix  & Strand  & Coil \\
    Structure & E, A, L, M, Q, K, R, H  & V, I, Y, C, W, F, T  & G, N, P, S, D \\
    \midrule
    Solvent  & Buried  & Exposed  & Intermediate \\
    Accessibility & A, L, F, C, G, I, V, W  & R, K, Q, E, N, D  & M, S, P, T, H, Y \\
    \bottomrule
    \end{tabular}
  }
  \label{tab:ctd}
\end{table}

For example, for a given sequence ``MTEITAAMVKELRESTGAGA'', it will be encoded as ``32132223311311222222'' according to its hydrophobicity division.

\textbf{Step 2: Compute Composition, Transition and Distribution Descriptors}

Three descriptors, \emph{Composition} (\emph{C}), \emph{Transition} (\emph{T}), and \emph{Distribution} (\emph{D}) were calculated for a given attribute as follows.

\textbf{\subsubsection{Composition}}

It is the global percent for each encoded class in the sequence. In the above example using hydrophobicity division, the numbers for encoded classes ``1'', ``2'', ``3'' are 5, 10, 5 respectively, so the compositions for them are $5/20=25\%$, $10/20=10\%$, and $5/20=25\%$ respectively, where 20 is the length of the protein sequence. Composition can be defined as

$$
C_r = \frac{n_r}{n} \quad r = 1, 2, 3
$$

where $n_r$ is the number of amino acid type $r$ in the encoded sequence and $N$ is the length of the sequence. An example for \code{extractCTDC()} could be:

%
<<extractCTDC>>=
extractCTDC(x)
@
%

The result shows the elements whose names are \code{PropertyNumber.GroupNumber} in the returned vector.

\textbf{\subsubsection{Transition}}

A transition from class 1 to 2 is the percent frequency with which 1 is followed by 2 or 2 is followed by 1 in the encoded sequence. Transition descriptor can be calculated as

$$
T_{rs} = \frac{n_{rs} + n_{sr}}{N - 1} \quad rs = \textrm{`12'}, \textrm{`13'}, \textrm{`23'}
$$

where $n_{rs}$, $n_{sr}$ is the numbers of dipeptide encoded as ``rs'' and ``sr'' respectively in the sequence and $N$ is the length of the sequence. An example for \code{extractCTDT()} could be:

%
<<extractCTDT>>=
extractCTDT(x)
@
%

\textbf{\subsubsection{Distribution}}

The ``distribution'' descriptor describes the distribution of each attribute in the sequence.

There are five ``distribution'' descriptors for each attribute and they are the position percents in the whole sequence for the first residue, 25\% residues, 50\% residues, 75\% residues and 100\% residues, respectively, for a specified encoded class. For example, there are 10 residues encoded as ``2'' in the above example, the positions for the first residue ``2'', the 2th residue ``2'' (25\%*10=2), the 5th ``2'' residue (50\%*10=5), the 7th ``2'' (75\%*10=7) and the 10th residue ``2'' (100\%*10) in the encoded sequence are 2, 5, 15, 17, 20 respectively, so the distribution descriptors for ``2'' are: 10.0 (2/20*100), 25.0 (5/20*100), 75.0 (15/20*100), 85.0 (17/20*100) , 100.0 (20/20*100), respectively.

Finally, an example for \code{extractCTDD()} could be:

%
<<extractCTDD>>=
extractCTDD(x)
@
%

\subsection{Conjoint Triad Descriptors}\label{ssec:ctraid}

Conjoint triad descriptors are proposed by \citet{shenjw}. These conjoint triad descriptors abstracts the features of protein pairs based on the classification of amino acids. In this approach, each protein sequence is represented by a vector space consisting of descriptors of amino acids. To reduce the dimensions of vector space, the 20 amino acids were clustered into several classes according to their dipoles and volumes of the side chains. The conjoint triad descriptors are calculated as follows:

\textbf{Step 1: Classification of Amino Acids}

Electrostatic and hydrophobic interactions dominate protein-protein interactions. These two kinds of interactions may be reflected by the dipoles and volumes of the side chains of amino acids, respectively. Accordingly, these two parameters were calculated, respectively, by using the density-functional theory method B3LYP/6-31G and molecular modeling approach. Based on the dipoles and volumes of the side chains, the 20 amino acids could be clustered into seven classes (See Table 2). Amino acids within the same class likely involve synonymous mutations because of their similar characteristics.

\begin{table}[htbp]
  \centering
  \caption{Classification of amino acids based on dipoles and volumes of the side chains}
    \begin{tabular}{cccc}
    \addlinespace
    \toprule
    No.   & Dipole Scale\tablefootnote{Dipole Scale (Debye): $-$, Dipole < 1.0; $+$, 1.0 < Dipole < 2.0; $++$, 2.0 < Dipole < 3.0; $+++$, Dipole > 3.0; $+'+'+'$, Dipole > 3.0 with opposite orientation.} & Volume Scale\tablefootnote{Volume Scale (\AA$^3$): $-$, Volume < 50; $+$, Volume > 50.} & Class \\
    \midrule
    1     & $-$      & $-$     & Ala, Gly, Val \\
    2     & $-$      & $+$     & Ile, Leu, Phe, Pro \\
    3     & $+$      & $+$     & Tyr, Met, Thr, Ser \\
    4     & $++$     & $+$     & His, Asn, Gln, Tpr \\
    5     & $+++$    & $+$     & Arg, Lys \\
    6     & $+'+'+'$ & $+$     & Asp, Glu \\
    7     & $+$\tablefootnote{Cys is separated from class 3 because of its ability to form disulfide bonds.} & $+$     & Cys \\
    \bottomrule
    \end{tabular}
  \label{tab:ctriad}
\end{table}

\textbf{Step 2: Conjoint Triad Calculation}

The conjoint triad descriptors considered the properties of one amino acid and its vicinal amino acids and regarded any three continuous amino acids as a unit. Thus, the triads can be differentiated according to the classes of amino acids, i.e., triads composed by three amino acids belonging to the same classes, such as ART and VKS, could be treated identically. To conveniently represent a protein, we first use a binary space  $(\mathbf{V}, \mathbf{F})$ to represent a protein sequence. Here, $\mathbf{V}$ is the vector space of the sequence features, and each feature $v_i$ represents a sort of triad type; $\mathbf{F}$ is the frequency vector corresponding to $\mathbf{V}$, and the value of the $i$-th dimension of $\mathbf{F} (f_i)$ is the frequency of type $v_i$ appearing in the protein sequence. For the amino acids that have been catogorized into seven classes, the size of $\mathbf{V}$ should be $7 \times 7 \times 7$; thus $i = 1, 2, \ldots, 343$. The detailed description for ($\mathbf{V}$, $\mathbf{F}$) is illustrated in Figure \ref{fig:ctriad}.

\begin{figure}[htbp]
\centering
\includegraphics[width=\textwidth]{fig/ctriad.pdf}
\caption{Schematic diagram for constructing the vector space ($\mathbf{V}$, $\mathbf{F}$) of protein sequence. $\mathbf{V}$ is the vector space of the sequence features; each feature ($v_i$) represents a triad composed of three consecutive amino acids; $\mathbf{F}$ is the frequency vector corresponding to $\mathbf{V}$, and the value of the $i$-th dimension of $\mathbf{F} (f_i)$ is the frequency that $v_i$ triad appeared in the protein sequence.} \label{fig:ctriad}
\end{figure}

Clearly, each protein correlates to the length (number of amino acids) of protein. In general, a long protein would have a large value of $f_i$, which complicates the comparison between two heterogeneous proteins. Thus, we defined a new parameter, $d_i$, by normalizing $f_i$ with the following equation:

$$
d_i = \frac{f_i - \min\{\,f_1, f_2 , \ldots, f_{343}\,\}}{\max\{\,f_1, f_2, \ldots, f_{343}\,\}}
$$

The numerical value of $d_i$ of each protein ranges from 0 to 1, which thereby enables the comparison between proteins. Accordingly, we obtain another vector space (designated $\mathbf{D}$) consisting of $d_i$ to represent protein.

To compute conjoint triads of protein sequences, we could simply use:

%
<<extractCTriad>>=
ctriad = extractCTriad(x)
head(ctriad, n = 65L)
@
%

by which we only outputted the first 65 of total 343 dimension to save space.

\subsection{Quasi-sequence-order Descriptors}\label{subsec:qsod}

The quasi-sequence-order descriptors are proposed by \citet{chouqsoe}. They are derived from the distance matrix between the 20 amino acids.

\textbf{\subsubsection{Sequence-order-coupling Number}}\label{sssec:socn}

The $d$-th rank sequence-order-coupling number is defined as:

$$
\tau_d = \sum_{i=1}^{N-d} (d_{i, i+d})^2 \quad d = 1, 2, \ldots, \textrm{maxlag}
$$

where $d_{i, i+d}$ is the distance between the two amino acids at position $i$ and $i+d$.

\textbf{Note}: maxlag is the maximum lag and the length of the protein must be not less than $\textrm{maxlag}$.


The function \code{extractSOCN(x)} is used for computing the sequence-order-coupling numbers:

%
<<extractSOCN>>=
extractSOCN(x)
@
%

Users could also specify the maximum lag value with the \code{nlag} argument.

\textbf{Note}: In addition to Schneider-Wrede physicochemical distance matrix \citep{wrede} used by Kuo-Chen Chou, another chemical distance matrix by \citet{grantham} is also used here. So the descriptors dimension will be \code{nlag * 2}. The quasi-sequence-order descriptors described next also utilized the two matrices.

\textbf{\subsubsection{Quasi-sequence-order Descriptors}}\label{sssec:qso}

For each amino acid type, a quasi-sequence-order descriptor can be defined as:

$$
X_r = \frac{f_r}{\sum_{r=1}^{20} f_r + w \sum_{d=1}^{\textrm{maxlag}} \tau_d} \quad r = 1, 2, \ldots, 20
$$

where $f_r$ is the normalized occurrence for amino acid type $i$ and $w$ is a weighting factor ($w=0.1$). These are the first 20 quasi-sequence-order descriptors. The other 30 quasi-sequence-order are defined as:

$$
X_d = \frac{w \tau_{d-20}}{\sum_{r=1}^{20} f_r + w \sum_{d=1}^{\textrm{maxlag}} \tau_d} \quad d = 21, 22, \ldots, 20 + \textrm{maxlag}
$$

\begin{figure}[htbp]
\centering
\includegraphics[width=0.75\textwidth]{fig/QSO.png}
\caption{A schematic drawing to show (a) the 1st-rank, (b) the 2nd-rank, and (3) the 3rd-rank sequence-order-coupling mode along a protein sequence. (a) Reflects the coupling mode between all the most contiguous residues, (b) that between all the 2nd most contiguous residues, and (c) that between all the 3rd most contiguous residues. This figure is from \citet{chouqsoe}.} \label{fig:qso}
\end{figure}

An minimal example for \code{extractQSO()} could be:

%
<<extractQSO>>=
extractQSO(x)
@
%

where users could also specify the maximum lag with argument \code{nlag} and the weighting factor with argument \code{w}.

\subsection{Pseudo-Amino Acid Composition (PAAC)}\label{subsec:paac}

This groups of descriptors are proposed in \citet{choupaac}. PAAC descriptors are also called the \emph{type 1 pseudo-amino acid composition}. Let $H_1^o (i)$ , $H_2^o (i)$, $M^o (i)$ ($i=1, 2, 3, \ldots, 20$) be the original hydrophobicity values, the original hydrophilicity values and the original side chain masses of the 20 natural amino acids, respectively. They are converted to following qualities by a standard conversion:

$$
H_1 (i) = \frac{H_1^o (i) - \frac{1}{20} \sum_{i=1}^{20} H_1^o (i)}{\sqrt{\frac{\sum_{i=1}^{20} [H_1^o (i) - \frac{1}{20} \sum_{i=1}^{20} H_1^o (i) ]^2}{20}}}
$$

$H_2^o (i)$ and $M^o (i)$ are normalized as $H_2 (i)$ and $M (i)$ in the same way.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.75\textwidth]{fig/PAAC.png}
\caption{A schematic drawing to show (a) the first-tier, (b) the second-tier, and (3) the third-tiersequence order correlation mode along a protein sequence. Panel (a) reflects the correlation mode between all the most contiguous residues, panel (b) that between all the second-most contiguous residues, and panel (c) that between all the third-most contiguous residues. This figure is from \citet{choupaac}.} \label{fig:paac}
\end{figure}

Then, a correlation function could be defines as

$$
\Theta (R_i, R_j) = \frac{1}{3} \bigg\{ [ H_1 (R_i) - H_1 (R_j) ]^2 + [ H_2 (R_i) - H_2 (R_j) ]^2 + [ M (R_i) - M (R_j) ]^2 \bigg\}
$$

This correlation function is actually an averaged value for the three amino acid properties: hydrophobicity value, hydrophilicity value and side chain mass. Therefore we can extend this definition of correlation function for one amino acid property or for a set of n amino acid properties.

For one amino acid property, the correlation can be defined as:

$$
\Theta (R_i, R_j) = [H_1 (R_i) - H_1(R_j)]^2
$$

where $H (R_i)$ is the amino acid property of amino acid $R_i$ after standardization.

For a set of n amino acid properties, it can be defined as: where $H_k (R_i)$ is the $k$-th property in the amino acid property set for amino acid $R_i$.

$$
\Theta (R_i, R_j) = \frac{1}{n} \sum_{k=1}^{n} [H_k (R_i) - H_k (R_j)]^2
$$

where $H_k(R_i)$ is the $k$-th property in the amino acid property set for amino acid $R_i$.

A set of descriptors called sequence order-correlated factors are defined as:

\begin{align*}
\theta_1  & = \frac{1}{N-1} \sum_{i=1}^{N-1} \Theta (R_i, R_{i+1})\\
\theta_2  & = \frac{1}{N-2} \sum_{i=1}^{N-2} \Theta (R_i, R_{i+2})\\
\theta_3  & = \frac{1}{N-3} \sum_{i=1}^{N-3} \Theta (R_i, R_{i+3})\\
          & \ldots \\
\theta_\lambda  & = \frac{1}{N-\lambda} \sum_{i=1}^{N-\lambda} \Theta (R_i, R_{i+\lambda})
\end{align*}

$\lambda$ ($\lambda < L$) is a parameter to be chosen. Let $f_i$ be the normalized occurrence frequency of the 20 amino acids in the protein sequence, a set of $20 + \lambda$ descriptors called the pseudo-amino acid composition for a protein sequence can be defines as:

$$
X_c = \frac{f_c}{\sum_{r=1}^{20} f_r + w \sum_{j=1}^{\lambda} \theta_j} \quad (1 < c < 20)
$$

$$
X_c = \frac{w \theta_{c-20}}{\sum_{r=1}^{20} f_r + w \sum_{j=1}^{\lambda} \theta_j} \quad (21 < c < 20 + \lambda)
$$

where $w$ is the weighting factor for the sequence-order effect and is set as $w = 0.05$ in \pkg{protr} as suggested by Kuo-Chen Chou.

With \code{extractPAAC()}, we could compute the PAAC descriptors:

%
<<extractPAAC>>=
extractPAAC(x)
@
%

The \code{extractPAAC()} fucntion also provides the \code{props} and \code{customprops} arguments, which is similar to the functions for Moreau-Broto/Moran/Geary autocorrelation descriptors. For minor differences, see \code{?extracPAAC}. Users could specify the lambda parameter and the weighting factor with arguments \code{lambda} and \code{w}.

\textbf{Note}: In the work of Kuo-Chen Chou, the definition for ``normalized occurrence frequency'' was not given. In this work, we define it as the occurrence frequency of amino acid in the sequence normalized to $100\%$ and hence our calculated values are not the same as values by them.

\subsection{Amphiphilic Pseudo-Amino Acid Composition (APAAC)}\label{subsec:apaac}

Amphiphilic Pseudo-Amino Acid Composition (APAAC) was proposed in \citet{choupaac}. APAAC is also recognized as the \emph{type 2 pseudo-amino acid composition}. The definitions of these qualities are similar to the PAAC descriptors. From $H_1 (i)$ and $H_2 (j)$ defined before, the hydrophobicity and hydrophilicity correlation functions are defined respectively as:

\begin{align*}
H_{i, j}^1  & = H_1 (i) H_1 (j)\\
H_{i, j}^2  & = H_2 (i) H_2 (j)
\end{align*}

From these qualities, sequence order factors can be defines as:

\begin{align*}
\tau_1  & = \frac{1}{N-1} \sum_{i=1}^{N-1} H_{i, i+1}^1\\
\tau_2  & = \frac{1}{N-1} \sum_{i=1}^{N-1} H_{i, i+1}^2\\
\tau_3  & = \frac{1}{N-2} \sum_{i=1}^{N-2} H_{i, i+2}^1\\
\tau_4  & = \frac{1}{N-2} \sum_{i=1}^{N-2} H_{i, i+2}^2\\
        & \ldots \\
\tau_{2 \lambda - 1} & = \frac{1}{N-\lambda} \sum_{i=1}^{N-\lambda} H_{i, i+\lambda}^1\\
\tau_{2 \lambda}     & = \frac{1}{N-\lambda} \sum_{i=1}^{N-\lambda} H_{i, i+\lambda}^2
\end{align*}

\begin{figure}[htbp]
\centering
\includegraphics[width=\textwidth]{fig/APAAC.pdf}
\caption{A schematic diagram to show (\textbf{a1}/\textbf{a2}) the first-rank, (\textbf{b1}/\textbf{b2}) the second-rank and (\textbf{c1}/\textbf{c2}) the third-rank sequence-order-coupling mode along a protein sequence through a hydrophobicity/hydrophilicity correlation function, where $H_{i, j}^1$ and $H_{i, j}^2$ are given by Equation (3). Panel (a1/a2) reflects the coupling mode between all the most contiguous residues, panel (b1/b2) that between all the second-most contiguous residues and panel (c1/c2) that between all the third-most contiguous residues.  This figure is from \citet{chouapaac}.} \label{fig:apaac}
\end{figure}

Then a set of descriptors called \emph{Amphiphilic Pseudo-Amino Acid Composition} (\emph{APAAC}) are defined as:

$$
P_c = \frac{f_c}{\sum_{r=1}^{20} f_r + w \sum_{j=1}^{2 \lambda} \tau_j} \quad (1 < c < 20)
$$

$$
P_c = \frac{w \tau_u}{\sum_{r=1}^{20} f_r + w \sum_{j=1}^{2 \lambda} \tau_j} \quad (21 < u < 20 + 2 \lambda)
$$

where $w$ is the weighting factor and is taken as $w = 0.5$ in \pkg{protr} as in the work of Chou KC.

A minimal example for \code{extracAPAAC()} is:

%
<<extractAPAAC>>=
extractAPAAC(x)
@
%

This function has the same arguments as \code{extractPAAC()}.

\subsection{Profile-based Descriptors}\label{subsec:profile}

The profile-based descriptors for protein sequences are available in the \pkg{protr} package. The feature vectors of profile-based methods were based on the PSSM by running PSI-BLAST, and often show good performance. See \cite{ye2011assessment} and \cite{rangwala2005profile} for details. The functions \code{extractPSSM()}, \code{extractPSSMAcc()} and \code{extractPSSMFeature()} are used to generate these descriptors. Users need to install the NCBI-BLAST+ software package first to make the functions fully functional.

\section{Descriptors for Proteochemometric Modeling}\label{sec:pcm}\addcontentsline{toc}{section}{4. Descriptors for Proteochemometric Modeling}

Proteochemometric (PCM) modeling utilizes statistical modeling techniques to model ligand-target interaction space. The below descriptors implemented in \pkg{protr} are extensively used in Proteochemometric modeling.

\begin{itemize}
\item Scales-based descriptors derived by Principal Components Analysis
\begin{itemize}
\item Scales-based descriptors derived by Amino Acid Properties from AAindex (Protein Fingerprint)
\item Scales-based descriptors derived by 20+ classes of 2D and 3D molecular descriptors (Topological, WHIM, VHSE, etc.)
\end{itemize}
\item Scales-based descriptors derived by Factor Analysis
\item Scales-based descriptors derived by Multidimensional Scaling
\item BLOSUM and PAM matrix-derived descriptors
\end{itemize}

Note that each of the scales-based descriptor functions are freely to combine with the more than 20 classes of 2D and 3D molecular descriptors to construct highly customized scales-based descriptors. Of course, these functions are designed to be flexible enough that users could provide totally self-defined property matrices to construct scales-based descriptors.

For example, to compute the ``topological scales'' derived by PCA (using the first 5 principal components), one could use \code{extractDescScales()}:

%
<<extractDescScales>>=
x = readFASTA(system.file('protseq/P00750.fasta', package = 'protr'))[[1]]
descscales = extractDescScales(x, propmat = 'AATopo',
                               index = c(37:41, 43:47),
                               pc = 5, lag = 7, silent = FALSE)
@
%

the argument \code{propmat} involkes the \code{AATopo} dataset shipped with \pkg{protr} package, and the argument \code{index} selects the 37 to 41 and the 43 to 47 columns (molecular descriptors) in the \code{AATopo} dataset to use, the parameter \code{lag} was set for the Auto Cross Covariance (ACC) for generating scales-based descriptors of the same length. At last, we printed the summary of the first 5 principal components (standard deviation, proportion of variance, cumulative proportion of variance).

The result is a length 175 named vector, which is consistent with the descriptors before:

%
<<extractDescScales2>>=
length(descscales)
head(descscales, 15)
@
%

For another example, to compute the descriptors derived by BLOSUM62 matrix and use the first 5 scales, one could use:

%
<<extractBLOSUM>>=
x = readFASTA(system.file('protseq/P00750.fasta', package = 'protr'))[[1]]
blosum = extractBLOSUM(x, submat = 'AABLOSUM62',
                       k = 5, lag = 7, scale = TRUE, silent = FALSE)
@
%

The result is a length 175 named vector:

%
<<extractBLOSUM2>>=
length(blosum)
head(blosum, 15)
@
%

\textbf{Dealing with gaps.} In proteochemometrics, (sequence alignment) gaps can be very useful, since a gap in a certain position contains information. The \pkg{protr} package has built-in support for such gaps. We deal with the gaps by using a dummy descriptor to code for the 21\textsuperscript{st} type of amino acid. The function \code{extractScalesGap()} and \code{extractProtFPGap()} could be used to deal with such gaps. See \code{?extractScalesGap} and \code{?extractProtFPGap} for details.

\section{Similarity Calculation by Sequence Alignment}\addcontentsline{toc}{section}{5. Similarity Calculation by Sequence Alignment}

Similarity computation derived by local or global protein sequence alignment between a list of protein sequences is great need in the protein related research and applications. However, this sort of pairwise similarity computation often computationally intensive, especially when there exists many protein sequences. Luckily, this process is also highly parallelizable, the \pkg{protr} package integrates the function of parallelized similarity computation derived by local or global protein sequence alignment between a list of protein sequences.

The function \code{twoSeqSim()} calculates the alignment result between two protein sequences, and the function \code{parSeqSim()} calculates the pairwise similarity calculation with a list of protein sequences in parallel:

\begin{CodeInput}
> s1 = readFASTA(system.file('protseq/P00750.fasta', package = 'protr'))[[1]]
> s2 = readFASTA(system.file('protseq/P08218.fasta', package = 'protr'))[[1]]
> s3 = readFASTA(system.file('protseq/P10323.fasta', package = 'protr'))[[1]]
> s4 = readFASTA(system.file('protseq/P20160.fasta', package = 'protr'))[[1]]
> s5 = readFASTA(system.file('protseq/Q9NZP8.fasta', package = 'protr'))[[1]]
> plist   = list(s1, s2, s3, s4, s5)
> psimmat = parSeqSim(plist, cores = 4, type = 'local', submat = 'BLOSUM62')
> print(psimmat)
\end{CodeInput}

\begin{CodeOutput}
           [,1]       [,2]       [,3]       [,4]       [,5]
[1,] 1.00000000 0.11825938 0.10236985 0.04921696 0.03943488
[2,] 0.11825938 1.00000000 0.18858241 0.12124217 0.06391103
[3,] 0.10236985 0.18858241 1.00000000 0.05819984 0.06175942
[4,] 0.04921696 0.12124217 0.05819984 1.00000000 0.05714638
[5,] 0.03943488 0.06391103 0.06175942 0.05714638 1.00000000
\end{CodeOutput}

It should be noted that for a small number of proteins, calculating their pairwise similarity scores derived by sequence alignment in parallel may not significantly reduce the overall computation time, since each of the task only requires a relatively small time to finish, thus, computational overheads may exist and affect the performance. In testing, we used about 1,000 protein sequences on 64 CPU cores, and observed significant performance improvement comparing to the sequential computation.

Users should install the packages \pkg{foreach} and \pkg{doParallel} before using \code{parSeqSim()}, according to their operation system. The \pkg{protr} package will automatically decide which backend to use.

\section{Similarity Calculation by GO Semantic Similarity Measures}\addcontentsline{toc}{section}{6. Similarity Calculation by GO Semantic Similarity Measures}

The \pkg{protr} package also integrates the function of similarity score computation derived by Gene Ontology (GO) semantic similarity measures between a list of GO terms / Entrez Gene IDs.

The function \code{twoGOSim()} calculates the similarity derived by GO-terms semantic similarity measures between two GO terms / Entrez Gene IDs, and the function \code{parGOSim()} calculates the pairwise similarity with a list of GO terms / Entrez Gene IDs:

\begin{CodeInput}
# by GO Terms
> go1 = c('GO:0005215', 'GO:0005488', 'GO:0005515',
+         'GO:0005625', 'GO:0005802', 'GO:0005905')  # AP4B1
> go2 = c('GO:0005515', 'GO:0005634', 'GO:0005681',
+         'GO:0008380', 'GO:0031202')                # BCAS2
> go3 = c('GO:0003735', 'GO:0005622', 'GO:0005840',
+         'GO:0006412')                              # PDE4DIP
> glist = list(go1, go2, go3)
> gsimmat1 = parGOSim(glist, type = 'go', ont = 'CC')
> print(gsimmat1)
\end{CodeInput}

\begin{CodeOutput}
      [,1]  [,2]  [,3]
[1,] 1.000 0.077 0.055
[2,] 0.077 1.000 0.220
[3,] 0.055 0.220 1.000
\end{CodeOutput}

\begin{CodeInput}
# by Entrez gene id
> genelist = list(c('150', '151', '152', '1814', '1815', '1816'))
> gsimmat2 = parGOSim(genelist, type = 'gene')
> print(gsimmat2)
\end{CodeInput}

\begin{CodeOutput}
       150   151   152  1814  1815  1816
150  0.689 0.335 0.487 0.133 0.169 0.160
151  0.335 0.605 0.441 0.171 0.198 0.274
152  0.487 0.441 0.591 0.151 0.178 0.198
1814 0.133 0.171 0.151 0.512 0.401 0.411
1815 0.169 0.198 0.178 0.401 0.619 0.481
1816 0.160 0.274 0.198 0.411 0.481 0.819
\end{CodeOutput}

\section{ProtrWeb}\addcontentsline{toc}{section}{7. ProtrWeb}

The web service built on \pkg{protr}, namely \pkg{ProtrWeb}, is located at:

\begin{center}
\url{http://protr.org}
\end{center}

\pkg{ProtrWeb} (Figure \ref{fig:protrweb}) does not require any knowledge of programming for the users, it is a user-friendly and one-click-to-go online platform for computing the protein descriptors presented in the \pkg{protr} package.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.9\textwidth]{fig/protrweb.png}
\caption{A screenshot of the web server \pkg{ProtrWeb}} \label{fig:protrweb}
\end{figure}

A step-by-step instruction on how to use \pkg{ProtrWeb} could be found at:

\begin{center}
\url{http://protr.org/intro.html}
\end{center}

\section{Miscellaneous Tools}\label{sec:misc}\addcontentsline{toc}{section}{8. Miscellaneous Tools}

In this section, we will briefly introduce some useful tools provided by the \pkg{protr} package.

\subsection{Retrieve Protein Sequences from UniProt}

This function \code{getUniProt()} gets protein sequences from uniprot.org by protein ID(s). The input \code{ID} is a character vector specifying the protein ID(s). The returned sequences are stored in a list:

%
\begin{CodeChunk}

\begin{CodeInput}
> ids = c('P00750', 'P00751', 'P00752')
> prots = getUniProt(ids)
> print(prots)
\end{CodeInput}

\begin{CodeOutput}
[[1]]
[1] "MDAMKRGLCCVLLLCGAVFVSPSQEIHARFRRGARSYQVICRDEKTQMIYQQHQSWLRPVLRSNRVEYCWCN
SGRAQCHSVPVKSCSEPRCFNGGTCQQALYFSDFVCQCPEGFAGKCCEIDTRATCYEDQGISYRGTWSTAESGAECT
NWNSSALAQKPYSGRRPDAIRLGLGNHNYCRNPDRDSKPWCYVFKAGKYSSEFCSTPACSEGNSDCYFGNGSAYRGT
HSLTESGASCLPWNSMILIGKVYTAQNPSAQALGLGKHNYCRNPDGDAKPWCHVLKNRRLTWEYCDVPSCSTCGLRQ
YSQPQFRIKGGLFADIASHPWQAAIFAKHRRSPGERFLCGGILISSCWILSAAHCFQERFPPHHLTVILGRTYRVVP
GEEEQKFEVEKYIVHKEFDDDTYDNDIALLQLKSDSSRCAQESSVVRTVCLPPADLQLPDWTECELSGYGKHEALSP
FYSERLKEAHVRLYPSSRCTSQHLLNRTVTDNMLCAGDTRSGGPQANLHDACQGDSGGPLVCLNDGRMTLVGIISWG
LGCGQKDVPGVYTKVTNYLDWIRDNMRP"

[[2]]
[1] "MGSNLSPQLCLMPFILGLLSGGVTTTPWSLARPQGSCSLEGVEIKGGSFRLLQEGQALEYVCPSGFYPYPVQ
TRTCRSTGSWSTLKTQDQKTVRKAECRAIHCPRPHDFENGEYWPRSPYYNVSDEISFHCYDGYTLRGSANRTCQVNG
RWSGQTAICDNGAGYCSNPGIPIGTRKVGSQYRLEDSVTYHCSRGLTLRGSQRRTCQEGGSWSGTEPSCQDSFMYDT
PQEVAEAFLSSLTETIEGVDAEDGHGPGEQQKRKIVLDPSGSMNIYLVLDGSDSIGASNFTGAKKCLVNLIEKVASY
GVKPRYGLVTYATYPKIWVKVSEADSSNADWVTKQLNEINYEDHKLKSGTNTKKALQAVYSMMSWPDDVPPEGWNRT
RHVIILMTDGLHNMGGDPITVIDEIRDLLYIGKDRKNPREDYLDVYVFGVGPLVNQVNINALASKKDNEQHVFKVKD
MENLEDVFYQMIDESQSLSLCGMVWEHRKGTDYHKQPWQAKISVIRPSKGHESCMGAVVSEYFVLTAAHCFTVDDKE
HSIKVSVGGEKRDLEIEVVLFHPNYNINGKKEAGIPEFYDYDVALIKLKNKLKYGQTIRPICLPCTEGTTRALRLPP
TTTCQQQKEELLPAQDIKALFVSEEEKKLTRKEVYIKNGDKKGSCERDAQYAPGYDKVKDISEVVTPRFLCTGGVSP
YADPNTCRGDSGGPLIVHKRSRFIQVGVISWGVVDVCKNQKRQKQVPAHARDFHINLFQVLPWLKEKLQDEDLGFL"

[[3]]
[1] "APPIQSRIIGGRECEKNSHPWQVAIYHYSSFQCGGVLVNPKWVLTAAHCKNDNYEVWLGRHNLFENENTAQF
FGVTADFPHPGFNLSLLKXHTKADGKDYSHDLMLLRLQSPAKITDAVKVLELPTQEPELGSTCEASGWGSIEPGPDB
FEFPDEIQCVQLTLLQNTFCABAHPBKVTESMLCAGYLPGGKDTCMGDSGGPLICNGMWQGITSWGHTPCGSANKPS
IYTKLIFYLDWINDTITENP"
\end{CodeOutput}
\end{CodeChunk}
%

\subsection{Read FASTA Format files}

The \code{readFASTA()} function provides a convenient way to read protein sequences stored in FASTA format files. See \code{?readFASTA} for details. The returned sequences are stored in a named list, whose components are named with the protein sequences' names.

\subsection{Read PDB Format files}

The Protein Data Bank (pdb) file format is a textual file format describing the three dimensional structures of protein. The \code{readPDB()} function provides the function to read protein sequences stored in PDB format files. See \code{?readPDB} for details.

\subsection{Sanity Check of the Amino Acid Types}

The \code{protcheck()} function checks if the protein sequence's amino acid types are in the 20 default types, which returns a \code{TRUE} if all the amino acids in the sequence belongs to the 20 default types:

%
<<protcheck>>=
x = readFASTA(system.file('protseq/P00750.fasta', package = 'protr'))[[1]]
# A real sequence
protcheck(x)
# An artificial sequence
protcheck(paste(x, 'Z', sep = ''))
@
%

\subsection{Protein Sequence Partition}

The \code{protseg()} function partitions the protein sequences to create sliding windows. This is usually required when creating feature vectors for machine learning tasks. Users could specify a sequence \code{x}, and a character \code{aa}, one of the 20 amino acid types, and a positive integer \code{k}, which controls the window size (half of the window).

This function returns a named list, each component contains one of the segmentations (a character string), names of the list components are the positions of the specified amino acid in the sequence. See the example below:

%
<<protseg>>=
protseg(x, aa = 'M', k = 5)
@
%

\subsection{Auto Cross Covariance (ACC) Computation}

Auto Cross Covariance (ACC) is extensively used in the scales-based descriptors computation, this approach calculates the auto covariance and auto cross covariance for generating scale-based descriptors of the same length. Users could write their own scales-based descriptor functions with the help of \code{acc()} function in the \pkg{protr} package.

\subsection{Pre-computed 2D and 3D Descriptor Sets for the 20 Amino Acids}

The \pkg{protr} package ships with more than 20 pre-computed 2D and 3D descriptor sets for the 20 amino acids to use with the scales-based descriptors, see \code{data(package = 'protr')} for all the datasets included in the \pkg{protr} package.

\subsection{BLOSUM and PAM Matrices for the 20 Amino Acids}

The BLOSUM and PAM matrices for the 20 amino acids could be used to calculate BLOSUM and PAM matrix-derived descriptors with function \code{extractBLOSUM()}, the datasets are named in \code{AABLOSUM45}, \code{AABLOSUM50}, \code{AABLOSUM62}, \code{AABLOSUM80}, \code{AABLOSUM100}, \code{AAPAM30}, \code{AAPAM40}, \code{AAPAM70}, \code{AAPAM120}, and \code{AAPAM250}.

\subsection{Meta Information of the 20 Amino Acids}

As the reference, the \code{AAMetaInfo} dataset includes the meta information of the 20 amino acids used for the 2D and 3D descriptor calculation in the \pkg{protr} package. This dataset include each amino acid's name, one-letter representation, three-letter representation, SMILE representation, PubChem CID and PubChem link. See \code{data(AAMetaInfo)} for details.

\section{Summary}\addcontentsline{toc}{section}{9. Summary}

The summary of the descriptors in the \pkg{protr} package is listed in table \ref{tab:all}.

\begin{table}[htbp]
  \centering
  \caption{List of commonly used descriptors in \pkg{protr}}
  \scriptsize{
    \begin{tabular}{p{130pt}p{130pt}p{80pt}p{80pt}}
    \addlinespace
    \toprule
    Descriptor Group & Descriptor Name & Descriptor Dimension & Function Name\\
    \midrule
    Amino Acid Composition & Amino Acid Composition & 20 & \code{extractAAC()} \\
          & Dipeptide Composition & 400 & \code{extractDC()} \\
          & Tripeptide Composition & 8000 & \code{extractTC()} \\
    Autocorrelation & Normalized Moreau-Broto Autocorrelation & 240\footnotemark[1]  & \code{extractMoreauBroto()} \\
          & Moran Autocorrelation & 240\footnotemark[1] & \code{extractMoran()} \\
          & Geary Autocorrelation & 240\footnotemark[1] & \code{extractGeary()} \\
    CTD   & Composition & 21  & \code{extractCTDC()}, \code{extractCTDCClass()} \\
          & Transition & 21 & \code{extractCTDT()}, \code{extractCTDTClass()} \\
          & Distribution & 105 & \code{extractCTDD()}, \code{extractCTDDClass()} \\
    Conjoint Triad & Conjoint Triad & 343 & \code{extractCTriad()}, \code{extractCTriadClass()} \\
    Quasi-Sequence-Order & Sequence-Order-Coupling Number & 60\footnotemark[2] & \code{extractSOCN()} \\
          & Quasi-Sequence-Order Descriptors & 100\footnotemark[2] & \code{extractQSO()} \\
    Pseudo-Amino Acid Composition & Pseudo-Amino Acid Composition & 50\footnotemark[3] & \code{extractPAAC()} \\
          & Amphiphilic Pseudo-Amino Acid Composition & 80\footnotemark[4] & \code{extractAPAAC()} \\
    \bottomrule
    \end{tabular}
  }
  \label{tab:all}
\end{table}

\footnotetext[1]{The number depends on the choice of the number of properties of amino acids and the choice of the maximum values of the \code{lag}. The default is use 8 types of properties and \code{lag} = 30.}
\footnotetext[2]{The number depends on the maximum value of \code{lag}. By default \code{lag} = 30. And two distance matrices were used, so the descriptor dimension is $30 \times 2 = 60$ and $(20 + 30) \times 2 = 100$.}
\footnotetext[3]{The number depends on the choice of the number of the set of amino acid properties and the choice of the $\lambda$ value. The default is use 3 types of properties proposed by Kuo-Chen Chou and $\lambda = 30$.}
\footnotetext[4]{The number depends on the choice of the $\lambda$ vlaue. The default is that $\lambda = 30$.}

The summary of the scales-based PCM descriptors in the \pkg{protr} package is listed in table \ref{tab:all2}.

\begin{table}[htbp]
  \centering
  \caption{List of PCM descriptors in \pkg{protr}}
  \footnotesize{
    \begin{tabular}{p{130pt}p{200pt}p{90pt}}
    \addlinespace
    \toprule
    Derived by & Descriptor Class & Function Name \\
    \midrule
    Principal Components Analysis & Scales-based descriptors derived by Principal Components Analysis & \code{extractScales()}, \code{extractScalesGap()} \\
          & Scales-based descriptors derived by amino acid properties from AAindex  (a.k.a. Protein Fingerprint) & \code{extractProtFP()}, \code{extractProtFPGap()}\\
          & Scales-based descriptors derived by 2D and 3D molecular descriptors (Topological, WHIM, VHSE, etc.) & \code{extractDescScales()} \\
    Factor Analysis & Scales-based descriptors derived by Factor Analysis & \code{extractFAScales()} \\
    Multidimensional Scaling & Scales-based descriptors derived by Multidimensional Scaling & \code{extractMDSScales()} \\
    Substitution Matrix & BLOSUM and PAM matrix-derived descriptors & \code{extractBLOSUM()} \\
    \bottomrule
    \end{tabular}
  }
  \label{tab:all2}
\end{table}

The summary of the amino acid descriptor sets used by scales-based descriptors provided in the \pkg{protr} package is listed in table \ref{tab:all3}. Note that the non-informative descriptors (like the descriptors have only one value across all the 20 amino acids) in these datasets have already been filtered out.

\begin{table}[htbp]
  \centering
  \caption{List of the pre-calculated descriptor sets of the 20 amino acids in \pkg{protr}}
  \footnotesize{
    \begin{tabular}{lccc}
    \addlinespace
    \toprule
    Dataset Name & Descriptor Set Name & Dimensionality & Calculated by \\
    \midrule
    \code{AA2DACOR} & 2D Autocorrelations Descriptors & 92 & Dragon \\
    \code{AA3DMoRSE} & 3D-MoRSE Descriptors & 160 & Dragon \\
    \code{AAACF} & Atom-Centred Fragments Descriptors & 6 & Dragon \\
    \code{AABurden} & Burden Eigenvalues Descriptors & 62 & Dragon \\
    \code{AAConn} & Connectivity Indices Descriptors & 33 & Dragon \\
    \code{AAConst} & Constitutional Descriptors & 23 & Dragon \\
    \code{AAEdgeAdj} & Edge Adjacency Indices Descriptors & 97 & Dragon \\
    \code{AAEigIdx} & Eigenvalue-Based Indices Descriptors & 44 & Dragon \\
    \code{AAFGC} & Functional Group Counts Descriptors & 5 & Dragon \\
    \code{AAGeom} & Geometrical Descriptors & 41 & Dragon \\
    \code{AAGETAWAY} & GETAWAY Descriptors & 194 & Dragon \\
    \code{AAInfo} & Information Indices Descriptors & 47 & Dragon \\
    \code{AAMolProp} & Molecular Properties Descriptors & 12 & Dragon \\
    \code{AARandic} & Randic Molecular Profiles Descriptors & 41 & Dragon \\
    \code{AARDF} & RDF Descriptors & 82 & Dragon \\
    \code{AATopo} & Topological Descriptors & 78 & Dragon \\
    \code{AATopoChg} & Topological Charge Indices Descriptors & 15 & Dragon \\
    \code{AAWalk} & Walk and Path Counts Descriptors & 40 & Dragon \\
    \code{AAWHIM} & WHIM Descriptors & 99 & Dragon \\
    \code{AACPSA} & CPSA Descriptors & 41 & Accelrys Discovery Studio \\
    \code{AADescAll} & All the 2D Descriptors Calculated by Dragon & 1171 & Dragon \\
    \code{AAMOE2D} & All the 2D Descriptors Calculated by MOE & 148 & MOE \\
    \code{AAMOE3D} & All the 3D Descriptors Calculated by MOE & 143 & MOE \\
    \bottomrule
    \end{tabular}
  }
  \label{tab:all3}
\end{table}

In this manual, we discussed the functions of the \pkg{protr} package, which is trying to offer a comprehensive and unique toolkit for protein sequence descriptor calculation and similarity computation.

\section*{Acknowledgments}

This work is financially supported by the National Natural Science Foundation of China (Grant No. 11271374) and the Postdoctoral Science Foundation of Central South University. The studies meet with the approval of the university's review board.

\clearpage

\nocite{*} % list uncited bibs
\bibliography{protr}

\vspace*{-0.35cm}

\end{document}
