\documentclass{article}
\usepackage[affil-it]{authblk}
\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}


\usepackage[english]{babel}

\usepackage{nips2013,times}


\usepackage[utf8]{inputenc} % allow utf-8 input
\usepackage[T1]{fontenc}    % use 8-bit T1 fonts
\usepackage{booktabs}       % professional-quality tables
\usepackage{amsfonts}       % blackboard math symbols
\usepackage{nicefrac}       % compact symbols for 1/2, etc.
\usepackage{microtype}      % microtypography


\usepackage[space]{grffile}
\usepackage{natbib}
\usepackage{hyperref}       % hyperlinks
\usepackage{url}            % simple URL typesetting


\bibliographystyle{ieeetr}
\hypersetup{colorlinks=false,pdfborder={0 0 0}}

\nipsfinalcopy

\author{
Alex Rubinsteyn\thanks{\texttt{alex.rubinsteyn@mssm.edu}} \\
Icahn School of Medicine at Mount Sinai \\
\And
Timothy O'Donnell\thanks{\texttt{tim.odonnell@mssm.edu}} \\
Icahn School of Medicine at Mount Sinai \\
\And
Nandita Damaraju\thanks{\texttt{nandita94@gmail.com}} \\
Georgia Institute of Technology \\
\And
Jeff Hammerbacher\thanks{\texttt{jeff.hammerbacher@mssm.edu}} \\
Icahn School of Medicine at Mount Sinai \\
}

\bibliographystyle{ieeetr} 


\begin{document}

\title{Predicting Peptide-MHC Binding Affinities With Imputed Training Data}


 
\author{
Alex Rubinsteyn \\
Icahn School of Medicine at Mount Sinai \\
\And
Timothy O'Donnell \\
Icahn School of Medicine at Mount Sinai \\
\AND
Nandita Damaraju \\
Georgia Institute of Technology \\
\And
Jeff Hammerbacher \\
Icahn School of Medicine at Mount Sinai \\
}


\bibliographystyle{ieeetr}


\date{\today}

\maketitle

\selectlanguage{english}
\begin{abstract}
Predicting the binding affinity between MHC proteins and their peptide ligands is a key problem in computational immunology. State of the art performance is currently achieved by the allele-specific predictor NetMHC and the pan-allele predictor NetMHCpan, both of which are ensembles of shallow neural networks. We explore an intermediate between allele-specific and pan-allele prediction: training allele-specific predictors with synthetic samples generated by imputation of the peptide-MHC affinity matrix. We find that the imputation strategy is useful on alleles with very little training data. We have implemented our predictor as an open-source software package called MHCflurry and show that MHCflurry achieves competitive performance to NetMHC and NetMHCpan.%
\end{abstract}%



\section{Introduction}
%In organisms with adaptive immunity, cells expose fragments of proteins extracted from the cytosol. These peptide fragments, typically but not always 9 amino acids in length, are monitored by cytotoxic T cells, which recognize and kill infected or cancerous cells based on their viral, bacterial, mutant, or unusual peptides~\cite{Anderson_2004}. Presented peptides must bind to a major histocompatability complex I (MHC I) protein, which forms a platform for interaction with T cells. 

In most vertebrates, cytotoxic T-cells enforce multi-cellular order by killing infected or cancerous cells. Each organism possesses a poly-clonal army of T-cells which collectively are able to distinguish unhealthy cells from healthy ones. This amazing feat is achieved through the winnowing and expansion of clonal T-cell populations possessing highly specific T-cell receptors (TCRs)~\cite{Blackman_1990}. Each distinct TCR recognizes a small number of similar peptides bound to an MHC molecule on the surface of a cell~\cite{Huseby_2005}. Though there are many steps in ``antigen processing''~\cite{Cresswell_2005}, it has become apparent that MHC binding is the most restrictive step. Peptide-MHC affinity prediction is the well-studied problem of predicting the binding strength of a given peptide and MHC pair~\cite{Lundegaard_2007}. Early approaches focused on ``sequence motifs''\cite{Sette_1989}, followed by regularized linear models, linear models with interaction terms such as SMM with pairwise features~\cite{Peters_2003}, and more recently the NetMHC family of predictors, a collection of related models based on ensembles of neural networks. Two of these predictors, NetMHC~\cite{Lundegaard_2008} and NetMHCpan~\cite{Nielsen_2007}, have emerged as the methods of choice across multiple fields of study within immunology, including virology~\cite{Lund_2011}, tumor immunology~\cite{Gubin_2015}, and autoimmunity~\cite{Abreu_2012}. 

% Initial approaches to predicting MHC ligands focused on ``sequence motifs''\cite{Sette_1989}, which were quickly replaced by a variety of regularized linear models, which themselves are consistently outperformed by regularized linear models with interaction terms such as SMM~\cite{Peters_2003}. The march toward black box non-linear models reached its local maximum with the NetMHC family of predictors, which are a collection of related models that utilize ensembles of neural networks. T

% Though there are many steps in ``antigen processing''~\cite{Cresswell_2005}, it has become apparent that MHC binding is the most restrictive step and consequently the most important sub-problem of predicting T-cell epitopes.

% (original): In most vertebrates, cytotoxic T-cells enforce multi-cellular order by killing both infected and cancerous cells. Each individual organism possesses a poly-clonal army of T-cells which collectively are able to distinguish rare unhealthy cells from healthy ones. This amazing feat is achieved through the winnowing and expansion of clonal T-cell populations possessing highly specific T-cell receptors (TCRs)~\cite{Blackman_1990}. Each distinct TCR is able to recognize a small number of similar peptides bound to an MHC molecule on the surface of a cell~\cite{Huseby_2005}. Though there are many steps in ``antigen processing''~\cite{Cresswell_2005} (the process by which protein fragments find themselves loaded onto membrane-bound MHCs), it has become apparent that MHC binding is the most restrictive step and consequently the most important sub-problem of predicting T-cell epitopes. 

NetMHC is an {\it allele-specific} method which trains a separate predictor for each allele's binding dataset, whereas NetMHCpan is a {\it pan-allele} method whose inputs are vector encodings of both a peptide and a subsequence of a particular MHC molecule. The conventional wisdom is that NetMHC performs better on alleles with many assayed ligands whereas NetMHCpan is superior for less well-characterized alleles~\cite{Gfeller_2016}.

In this paper we explore the space between {\it allele-specific} and {\it pan-allele} prediction by imputing the unobserved values of peptide-MHC affinities for which we have no measurements and using these imputed values for pre-training of allele-specific binding predictors.

%Paper summary:
%\begin{enumerate}
%\item Use cross-validation over the training data to choose the best ``large'' network configuration for alleles with many samples, as well as a ``smaller'' network which performs best on alleles with few samples.
%\item Use cross-validation over the training data to choose a most accurate imputation/matrix completion algorithm for 
%synthesizing pMHC binding affinities.
%\item Compare performance of the the best ``small'' and ``large'' networks (with and without pre-training on the synthesized samples) on an unseen test set vs. NetMHC, NetMHCpan, and SMM. 
%\end{enumerate}

\section{Data and evaluation metrics}

Two datasets were used from a recent paper studying the relationship between training data and pMHC predictor accuracy\cite{Kim_2014}. The training dataset (BD2009) contained entries from IEDB~\cite{Salimi_2012} up to 2009 and the test dataset (BLIND) contained IEDB entries from between 2010 and 2013 which did not overlap with BD2009 (Table~\ref{tab:datasets}).\selectlanguage{english}
\begin{table}[h!]
\centering
\begin{tabular}{l||cccc}
\toprule
{} & Alleles &  IC50 Measurements & Expanded 9mers \\
\midrule
BD2009 &     106 &                           137,654 &        470,170 \\
BLIND  &      53 &                           27,680 &         83,752 \\
\bottomrule
\end{tabular}
\caption{{Train (BD2009) and test (BLIND) dataset sizes.}}
\label{tab:datasets}
\end{table}

Throughout this paper we will evaluate a pMHC binding predictor using three different metrics:

\begin{itemize}
\item {\bf AUC}: Area under the ROC curve. Estimates the probability that a ``strong binder'' peptide (affinity $\leq 500$nM) will be given a stronger predicted affinity than one whose ground truth affinity is $>500$nM.
\item {\bf F$_1$ score}: Measures trade-off between sensitivity and specificity for predicting ``strong binders'' with affinities $\leq 500$nM. 
\item {\bf Kendall's $\tau$}: Rank correlation across the full spectrum of binding affinities.
\end{itemize}



\section{Comparison of imputation algorithms as predictors}

A dataset of peptide-MHC affinities for $n$ peptides and $a$ alleles may be thought of as a $n \times a$ matrix where peptide/allele pairs without measurements are missing values. The task of predicting values at these positions is known as matrix completion or imputation (depending on the community and data source).  We investigated the performance of several imputation algorithms as a standalone solution to the peptide-MHC affinity prediction problem. The algorithms considered were:

\begin{itemize}
\item {\bf meanFill}: Replace each missing pMHC binding affinity with the mean affinity for that allele. This is a very simple imputation method which serves as a baseline against which other methods can be compared. 

\item {\bf knnImpute}~\cite{Troyanskaya_2001}: Each missing entry $X_{ij}$ is imputed using the values in the $k$ closest columns with observation in row $i$.  Similarity between alleles is computed as $e^{-d_{st}^2}$, where $d_{st}$ is the mean squared difference between observed entries of alleles $s$ and $t$. 

\item {\bf svdImpute}~\cite{Troyanskaya_2001}: Imputation using iterative fixed rank SVD decomposition. 

\item {\bf softImpute}~\cite{Mazumder2010SpectralMatrices}: A singular value thresholding method which iteratively estimates a low-rank matrix completion without forcing the pre-specification of a particular solution rank. Instead, the {\it softImpute} method is parameterized by a shrinkage value $\lambda$ that is subtracted from each singular value. 

\item {\bf MICE}~\cite{Azur_2011}: Average multiple imputations generated using Gibbs sampling from the joint distribution of columns. 
\end{itemize}

We evaluated the performance of these methods using three-fold cross validation on BD2009, only considering peptides which occurred in at least three alleles and excluding alleles with less than five measurements (Table~\ref{tab:imputation}). All imputation methods were implemented in the \textit{fancyimpute} Python library~\cite{fancyimpute-0-0-16}. Since MICE outperformed the other methods on two of the three predictor metrics, we selected it for the subsequent neural network experiments.\selectlanguage{english}
\begin{table}[htbp]
\centering
\begin{tabular}{cl||ccc}
\toprule
Imputation Method & Parameter & AUC & $F_1$ score & Kendall's $\tau$ \\
\midrule 
meanFill & &   0.67665 &  0.04950 &  0.17675 \\
\midrule
knnImpute & $k = 1$ &  0.80907 &  0.57952 &  0.40201 \\
  & $k = 3$  &  0.83189 &  0.57594 &  0.42086 \\
  & $k = 5$ &  0.83103 &  0.56118 &  0.41703 \\
\midrule
MICE & $n = 25$  &  0.85861 &  0.57597 &  0.44978 \\
     & $n = 50$  &  {\bf 0.86127} &  0.56527 &  {\bf 0.45944} \\
\midrule
softImpute & $\lambda=5$ &  0.78981 &  0.39158 &  0.33408 \\
& $\lambda=10$ &  0.83248 &  0.53575 &  0.39763 \\
& $\lambda=20$ &  0.85608 &  {\bf 0.60599} &  0.43754 \\
             
\midrule
svdImpute & rank = 5  &  0.82305 &  0.57040 &  0.39117 \\
& rank = 10  &  0.83667 &  0.58433 &  0.40048 \\
& rank = 20  &  0.82986 &  0.57038 &  0.38817 \\
\bottomrule[1.25pt]
\end{tabular}
\begin{center}
\caption{{Cross-validation performance of imputation algorithms on BD2009 dataset}} \label{tab:imputation}
\end{center}
\end{table}



\section{Neural network architecture}
Each MHCflurry predictor is a feed-forward neural network containing (1) an embedding layer which transforms amino acids to learned vector representations, (2) a single hidden layer with $tanh$ nonlinearity, (3) a sigmoidal scalar output. This network is implemented using Keras~\cite{chollet2015}.  

%% \begin{figure}[h]
%% \centering
%% \includegraphics[scale=0.5]{figures/mhcflurry-gliffy-network.png}
%% \caption{MHCflurry architecture}
%% \label{fig:architecture}
%% \end{figure}

Three-fold cross validation on the training set was used to select the hyper-parameters. The best model had 32 output dimensions for the amino acid vector embedding, a hidden layer size of 64, a dropout rate of 50\%, and 250 training epochs. These hyper-parameters achieved reasonable performance across alleles, but it's likely that performance could be further improved by setting the hyper-parameters separately for each allele.

\section{Data encoding}

Like the NetMHC family of predictors~\cite{lundegaard2008accurate}, MHCflurry uses fixed length 9mer inputs which requires peptides of other lengths to be transformed into multiple 9mers. Shorter peptides are mapped onto 9mer query peptides by introducing a sentinel ``X'' at every possible position, whereas longer peptides are cut into multiple 9mers by removing consecutive stretches of residues. The predicted affinity for a non-9mer peptide is the geometric mean of the predictions for the generated 9mers. When $n$ training samples derive from a single non-9mer sequence then their weights are adjusted to $1/n$. 

We map IC50 concentrations onto a regression targets between 0.0 and 1.0 using the same scheme as NetMHC, $y = 1.0 - max(1.0, log_{50000}(IC50))$.

\section{Training}

For each allele, we train a MHCflurry model using the measured peptide affinities for the allele and the values imputed by MICE based on other alleles in the training set. As training progresses, we place quadratically decreasing weight on the imputed values.

A randomly generated peptide is unlikely to bind a given MHC strongly, but a data acquisition bias toward strong binders in the training set can lead models to assign a high affinity to most peptides. As a form of regularization, we augment the training set at each epoch to include random peptides with affinity set to be maximally weak. The number of random negative peptides is 20\% of the training size (without imputation). At each training epoch, a fresh set of random peptides is generated.


\section{Results}
We evaluated the effect of imputation by drawing subsets of the BD2009 training set for the well-characterized allele HLA-A*02:01. Predictors were trained on a range of simulated training set sizes and tested on the remaining data (Figure~\ref{fig:imputecomparison}). We find that imputation gives a modest improvement up to approximately 100 training samples. With more training data there is no benefit to imputation.

% The figure doesn't show up in the preview but does show up if you export to pdf.\selectlanguage{english}
\begin{figure}[h]
\centering
\includegraphics[scale=0.75]{figures/HLA-A0201-vs-nsamples-hidden-64-activation-tanh-impute-mice-epochs-250-embedding-32-pretrain-quadratic.pdf}
\caption{{MHCflurry performance on down-sampled training data for HLA-A*02:01 with and without imputation}}
\label{fig:imputecomparison} 
\end{figure}

We then compared the performance of MHCflurry against NetMHC, NetMHCpan, and SMM on the blind test data. The MHCflurry ensemble model contains 32 predictors initialized with different random weights. The MHCflurry ensemble is competitive with NetMHC and NetMHCpan.\selectlanguage{english}
\begin{table}[h]
\centering
\begin{tabular}{llll}
\toprule
{} &               AUC &       $F_1$ score &  Kendall's $\tau$ \\
\midrule
MHCflurry (ensemble)        &  0.93260 &           0.78459 &   \textbf{0.58686} \\
MHCflurry (single predictor)    &           0.93225 &           0.78106 &           0.58572 \\
NetMHC                          &           0.93234 &  \textbf{0.80722} &   0.58633 \\
NetMHCpan                       &  \textbf{0.93264} &           0.79957 &           0.58138 \\
SMM-PMBEC                       &           0.92134 &           0.79026 &           0.56488 \\
\bottomrule
\end{tabular}

\caption{{Performance on BLIND dataset}}
\label{tab:measurementweighted}
\end{table}

\section{Discussion}
Imputing training data shows promise in cross-validation as a way to improve performance on alleles with few observations, but only seems to help for very small training sizes ($\leq 100$). Unfortunately, none of the alleles included in the BLIND dataset had fewer than 100 samples in BD2009, and only one had fewer than 200. Thus, additional work is required to assess the accuracy of MHCflurry and other predictors on alleles with scarce training data. Additionally, we need to further investigate the interaction between imputation parameters, the  decay schedule for the weights of imputed samples, and stopping criteria for training individual allele-specific predictors.

% These are generated in the 'paper plots' notebook; do not edit by hand.


\section{Code}
MHCflurry is available at \href{https://github.com/hammerlab/mhcflurry}{https://github.com/hammerlab/mhcflurry}. The data, scripts, and notebooks used to generate the plots and tables in this paper are available at \href{https://github.com/hammerlab/mhcflurry-icml-compbio-2016}{https://github.com/hammerlab/mhcflurry-icml-compbio-2016/}.

\selectlanguage{english}
\FloatBarrier
\bibliographystyle{plain}
\bibliography{bibliography/converted_to_latex.bib%
}

\end{CJK}\end{document}

