diff options
Diffstat (limited to 'deliverable/main.tex')
| -rw-r--r-- | deliverable/main.tex | 1291 |
1 files changed, 1291 insertions, 0 deletions
diff --git a/deliverable/main.tex b/deliverable/main.tex new file mode 100644 index 0000000..521be50 --- /dev/null +++ b/deliverable/main.tex @@ -0,0 +1,1291 @@ +% hello +\input{../preamble.tex} + +\makeglossaries +\input{../bussiness_glossary.tex} +\input{../data_mining_glossary.tex} +\input{../acronyms.tex} + +\begin{document} +\MyTitle{Change in immune cell signaling upon repeat vaccination: a data exploration using the FluPrint database} +\tableofcontents +\printglossary[type=bus] +\printglossary[type=dm] +\printglossary[type=\acronymtype] + + +\section{background} + +Influenza viruses are enveloped \gls{bu:rnaVirus} (\acrshort{rna} virus(es)) and +are divided into three types on the basis of \gls{bu:antigen}ic differences of internal +structural proteins \citep{fdaGuidanceIndustryClinical2007}. + +Two influenza virus types, Type A and B, cause yearly epidemic outbreaks in humans +and are further classified based on the structure of two major external +\gls{bu:glycoprotein}s, hemagglutinin (\acrshort{ha}) and neuraminidase (\acrshort{na}) +\citep{fdaGuidanceIndustryClinical2007}. + +Type B viruses, which are largely restricted to the human host, have a single +\acrshort{ha} and \acrshort{na} subtype. In contrast, numerous \acrshort{ha} +and \acrshort{na} Type A influenza subtypes have been identified to date. Type +A and B influenza variant strains emerge as a result of frequent +\gls{bu:antigen}ic change, principally from \gls{bu:mutation}s in the \acrshort{ha} +and \acrshort{na} \gls{bu:glycoprotein}s \citep{fdaGuidanceIndustryClinical2007}. + +Since 1977, influenza A virus subtypes H1N1 and H3N2, and influenza B viruses +have been in global circulation in humans. The current U.S. licensed +\gls{bu:tiv} are formulated to prevent influenza illness +caused by these influenza viruses. Because of the frequent emergence of new +influenza variant strains, the \gls{bu:antigen}ic composition of influenza vaccines +needs to be evaluated yearly, and the \gls{bu:tiv} are reformulated almost every +year. + +Currently, even with full production, manufacturing capacity would not produce +enough seasonal influenza vaccine to vaccinate all those for whom the vaccine +is now recommended \citep{fdaGuidanceIndustryClinical2007}. + +\subsection{Influenza mortality estimation models} + +Numerous works apply regression models to describe seasonal population +influenza mortality \citep{zhouHospitalizationsAssociatedInfluenza2012, +greenMortalityAttributableInfluenza2013, iulianoEstimatesGlobalSeasonal2018}. +Reported are varying age-specific influenza burdens during different seasonal +epidemics for different regions, but in general young children an elderly are +found to be more susceptible to influenza and are adviced to vaccinated +annually \citep{zhouHospitalizationsAssociatedInfluenza2012}. + +Specifically, within the US based work of +\cite{zhouHospitalizationsAssociatedInfluenza2012}, the highest hospitalization +rates for influenza were among persons aged $>=$65 years and those aged $<$1 +year. And, age-standardized annual rates per 100000 person-years varied +substantially for influenza. A similar pattern is in +\cite{greenMortalityAttributableInfluenza2013}, where an age shift in Wales and +England seasonal influenza burden was observed following the 2009 swine flue +pandemic. It is also estimated that globally 291.243–645.832 influenza associated +seasonal deaths occur annually \citep{iulianoEstimatesGlobalSeasonal2018} These +varying demographic statistics and the volume of influenza patients can confound +decision making on national and international public health policies. +Knowledge on vaccine efficacy and implementation can be a valuable asset for +fighting future seasonal influenza outbreaks. + +\subsection{Vaccine success criteria} + +Due to the volume and vulnerability of population groups most at risk for +influenze, the young and the elderly, a placebo controlled vaccine efficacy +study is extremely costly \citep{zhouHospitalizationsAssociatedInfluenza2012}. +Instead the haemagglutination-inhibiting (HAI) antibody test for influenza +virus antibody is used to assess vaccine protection +\citep{dejongHaemagglutinationinhibitingAntibodyInfluenza2003}. The policy for +a succesful vaccine is an 4-fold increase in HAI antibody titre after +vaccination and a geometric mean HAI titer of $\geq$ 40. The last is predicted +to reduce influenza risk by 50\% +\cite{dejongHaemagglutinationinhibitingAntibodyInfluenza2003}. + +\subsection{Finding immunological factors predicting high vaccine response using machine learning} + +It is known that pre-existing T cell populations are correlated with a HAI +antibody response after vaccination. But, the role of T cells in mediating that +response is uncertain. In one work it was found that under certain +circumstances CD8+ T cells specific to conserved viral epitopes correlated with +protection against symptomatic influenza +\citep{sridharCellularImmuneCorrelates2013}.In other work, populations of CD4+ +T cells that associated with protective antibody responses after seasonal +influenza vaccinations were found \citep{bentebibelInductionICOSCXCR3}. +\cite{trieuLongtermMaintenanceInfluenzaSpecific2017} reports a stable CD8+ T +cell populations and an increased CD4+ T cell populatin after vaccination. It was +also reported that repeat vaccinations are an important factor in maintaining +CD4+ T cell population \citep{trieuLongtermMaintenanceInfluenzaSpecific2017}. +How exactly these T cell populations factor into protective influenza immunity +and vaccination reponse is not well understood. + +Machine learning has been applied to clinical datasets to find influenza +protection markers, such as the described T cell populations and titers of +related molecules \citep{furmanApoptosisOtherImmune2013, +sobolevAdjuvantedInfluenzaH1N1Vaccination2016, tsangGlobalAnalysesHuman2014}. +These type of studies suffer from data quality issues, such as: inconsistencies +between findings depending on the epidemic season, only focussing on one type +of biological assay to get data, and a low amount of patients/samples. A +succesful vaccination is also often not well defined. + +\subsection{Bussiness objectives} + +Due to the high volume population that needs vaccines, it is important to study +immune correlates to vaccine response. For example, repeat vaccination might +not be necessary if the response is low, or a different vaccine is desired on a +person to person basis depending on immune correlates. Moreover, identifying +patterns between vaccine response and immune correlates furthers the +understanding of the underlying immunological mechanism of influenza +protection. + +This work uses the \flup database, which aims to solve data quality issues +and low dimensionality of prior studies using clinical datasets comprised of +viurs, cell and serum sample assays. It does so by incorporating eigth clinical +studies conducted between 2007 to 2015 using in total 740 patients, including +different types of assays and normalizing their values, and by providing a +binary classification of high- and low-responder to a vaccine. + +The objectives of this work are to answer: +\begin{itemize} + \item What kind of studies can be done using the \flup database? + \item What immunological factors correlate to a vaccine responses? + \item What is the effect of repeat vaccination? +\end{itemize} + +Since this work is an independent study performed for an assignment, the +success criteria for these objective will be loosely defined as providing a +statistical description or to provide insigth in the questions posed in the +objectives. + +The rationale for these questions and succes criteria are based on the scope +of the 3EC project as part of the Applied data science profile and the data +available. The paper of \cite{tomicFluPRINTDatasetMultidimensional2019} on +which this work is mostly based on provides these questions as interesting +directions for further analysis, but does not directly provide the data +necessary to answer them, only the MySQL database containing a great volume of +data. + +\section{Assess situation} + +\subsection{data and knowledge sources} + +The sole source of data used in the project is provided by +\cite{tomicFluPRINTDatasetMultidimensional2019} (this work is reffered to +using: "the \flup paper" from now on). \Dpaper describes the MySQL database for +which the installation is described in the +\href{https://github.com/LogIN-/fluprint}{FluPrint Github Repository}. A +template query is provided on the +\href{https://github.com/LogIN-/simon-manuscript}{github page} belonging to an +unpublished follow-up study by the same authors of \dpaper +\autoref{lst:QueryTemplate}. According to the authors, this data is the most +interesting for the bussiness objective of finding repeat vaccination effects +and will be used in this work too (this +unpublished follow-up study is referred to using: "\spaper"). The authors give this brief +description of the data: + +\begin{displayquote} + \textit{"The influenza datasets were obtained from the Stanford Data Miner maintained by + the Human Immune Monitoring Center at Stanford University. This included + total of 177 csv files, which were automatically imported to the MySQL + database to facilitate further analysis. The database, named \flup and + its source code, including the installation tutorial are freely available + here and on project's website. Following database installation, you can + obtain data used in the SIMON publication by following MySQL database + query \autoref{lst:QueryTemplate}"}. +\end{displayquote} + +\subsection{Tools and techniques} + +Installation of the \flup database will require an installation on a +unix operating system of \href{https://www.mysql.com/}{MySQL}, +\href{https://www.php.net/manual/en/install.php}{PHP}. More details are at the +\href{https://github.com/LogIN-/fluprint}{FluPrint Github Repository}. + +Database querying was done using a \href{https://neovim.io/}{neovim} based toolset, +personal configuration can be found +\href{https://github.com/Vinkage/mike_neovim/tree/feature}{here}. + +Since in \dpaper R is used, it is also used here. Especially crucial is the +\href{https://cran.r-project.org/web/packages/mulset/index.html}{R package +mulset}, which was made by the authors of \spaper. This package is used to deal +with missing data between different clinical studies and years, and thus will +be used to generate complete data tables in this paper too. All scripts in this +work were written using \href{https://www.tidyverse.org/}{tidyverse} packages +and make heavy use of the \href{https://dplyr.tidyverse.org/}{dplyr} package +for data wrangling. Additionally the following packages were used: +\href{https://cran.r-project.org/web/packages/ggpubr/index.html}{ggpubr} for +making publication quality figures, the kable function from +\href{knitr}{https://www.r-project.org/nosvn/pandoc/knitr.html} to generate +latex tables, \href{https://topepo.github.io/caret/}{caret} and +\href{https://cran.r-project.org/web/packages/MLeval/index.html}{MLeval} to +streamline model training and evaluation, +\href{https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html}{corrplot} +to visualise correlation between features, and other packages that were used +only once. + +\subsection{Requirements of the project} + +Requirements of this work are to show ability in using data science methods. +As such, most of the insights will inevitably be a replication of the work done +by the authors of the FluPrint database \cite{tomicSIMONAutomatedMachine2019}, +but all the scripts and analysis done are original work and are supplied +together with the final deliverable. + +Since the data type used here is a database this makes it more complicated for +an examinator to reproduce all code, especially since installing the database +requires a unix operating system. This is not considered problematic +since the queried tables from the database will be included in the final +deliverable. + +Reporting of the project follows the CRISP-DM methodology, where at each +stage of the project a separate report is written during the analysis work. In +the end the most important information is kept and incorporated in a final +report that is assumed to be graded in conjunction with the code. + +\subsection{Assumptions of the project} + +This work assumes that the focus point of the evaluation lies on the +methodology used, and the ability to apply the basic data science methods +learned in the Applied Data Science profile. The answer to business objectives +is assumed to be subjective, and it is assumed that the methods used and +clarity of insights into the data gained are more important. + +It is also assumed that the FluPrint database and other methods used by the +authors \cite{tomicFluPRINTDatasetMultidimensional2019, +tomicSIMONAutomatedMachine2019} are of high quality, and that this is +appropriate for this work. Out of the scope of this work is investigating +whether the preprocessing done for the data in the database is valid, since we +are not domain experts. A method for querying, cleaning, and generating +complete data tables has been provided by the authors and will also be used in +this work. It is assumed that the SQL and R methods (in particular the mulset R +package) in question are allowed to be used as a starting point in this +assignment. + +\subsection{Constraints of the project} + +This work is an unsupervised assignment, and only personal hardware were +available. This put constraints on dataset size and computational requirements +of analyses. The work was done on a Macbook air (2017) with the OSX big-sur +operating system. This means that unix tools were available and there were no +technical constraints. The filetypes are only csv files generated by the SQL +server. + +\section{Data mining goals} + +\subsection{Translating the problem in data mining terms} +All bussiness objectives described involve querying data from the FluPrint +database. The goal of the authors of the FluPrint database was to provide a +unqiue opportunity to study immune correlates of high vaccine responders across +different years and clinical studies. The authors also provide a binary +classification for donors. In this work we first and foremost explore the +database, and lastly we apply feature selection methods and classification +models on the most interesting dataset. + +The bussiness objectives can be translated in data mining terminology like so: +\begin{itemize} + \item Explore and describe the database and corresponding tables. + \item Apply wrapper feature selection to the most interesting datasets. + \item Explore features identified by the models trained in the wrapper feature selection. +\end{itemize} + +In data mining terms, the problem type is a combination of exploratory data +analysis and classification. Since this work is for a 2-weeks/3EC assignment +for the Applied Data Science profile, success criteria for all goals are +subjective. For the classification type goals we follow +the model evaluation procedure used by the authors +\cite{tomicSIMONAutomatedMachine2019}, models were evaluated by the AUROC +metric, and accuracy, specificity and sensitivity were also reported. Insights +produced by this work were benchmarked against the work of the original +authors. + +\subsection{Project plan} + +\f{v2_desc_exploration} +{Project plan for the SQL related data mining goal.} +{plan:sql} + +The first part of the project involved querying the database, and collecting +and describing the available data \autoref{plan:sql}. The first goal is to +understand the tables in the SQL database, their key relations, and to describe +the attributes within the tables. Valuable info on this part is already +provided in the original publication of the database +\cite{tomicFluPRINTDatasetMultidimensional2019}, but it was also investigated +in this work. The tools that will be used are SQL for querying and R for +statistical descriptions. + +% The second phase of this plan was an iterative process of finding suitable data +% to answer the modeling and visualisation data mining goals. This is a more +% involved process since it requires exploration of the database to answer the +% questions, and therefore was estimated to take time. + +% \f{model_and_vis_plan} +% {Project plan for the modeling and visualisation data mining goals.} +% {plan:vis} +% +% Relations between attributes in the generated datasets are visualised and +% modelled to see if there exist a pattern in the data that is relevant for the +% business objectives \autoref{plan:vis}. A critical point in this plan is +% deciding whether an objective cannot be answered with the available data. In +% that case the goal was revised and the second phase of the SQL query plan was +% reiterated. When deciding if the exploratory analysis was of sufficient +% quality, the work by the authors of the database used in this work was used as +% a subjective benchmark \cite{tomicSIMONAutomatedMachine2019, +% tomicFluPRINTDatasetMultidimensional2019}. + +\f{feature_selection_classification} +{Project plan for the classification and feature selection data mining goal.} +{plan:cls} + +For the modeling data mining goals the plan was to find the immune correlates +of high immune responders using a wrapper based feature selection strategy +\autoref{plan:cls} + + +\section{Data description} + +\subsection{Volumetric analysis} + +In the work of \cite{tomicFluPRINTDatasetMultidimensional2019} data on +indiviuals enrolled in influenza vaccine studies at the Stanford-LPCH Vaccine +Program was collected, the data was archived at the Stanford Data Miner. This +archive was filtered by assays used in influenza studies, resulting in data +from 740 healthy donors, enrolled in influenza vaccine studies conducted by the +Stanford-LPCH Vaccine Program from 2007 to 2015. These studies are described in +the table accompanying the online publication of the fluprint dataset +\autoref{tbl:studiesDesc}. From those 740 donors a vaccine response +classification was only given for 372 donors \autoref{fig:demoGraph}, by a +method that will be described in the section describing the data table +containing this attribute. Overall there was no major difference in demographic +statistics when stratisfying the data in high or low responder classification +\autoref{fig:demoGraph}. + +Importantly, it is reported that in all studies the donors are only vaccinated +once, except in the study SLVP015 \autoref{tbl:studiesDesc} +\citep{tomicFluPRINTDatasetMultidimensional2019}. However, in later work of the +same authors it is claimed that vaccines are administered as specified by the study +\citep{tomicSIMONAutomatedMachine2019}. + +The donors for which a vaccine respone classification was available from all +clinical studies together span a wide age range \autoref{fig:demoGraph}A from 1 +- 50 \autoref{tbl:demoStats}, in the original work the demographic statistics +include the donors for which no vaccine response classification is given, +therefore they report a greater range of 1-90. Stratisfying the donors on +vaccine response does not affect the demographic attribute distribution, but +the maximum age is lowered in the high responders group +\autoref{fig:demoGraph}B. + +\begin{figure}[htpb] + \includegraphics[width=\textwidth]{demographic} + \caption{\textbf{A.} percentage of donors with factor property within high + and low responder groups. Included are sex, race, and CMV status + information. \textbf{B.} Age distribution of donors with a known response + classification.}\label{fig:demoGraph} +\end{figure} + + + +\begin{table}[htpb] +\centering +\begin{tabular}{ll} +\toprule +\textbf{Age (y)} & \\ +\midrule +Mean $\pm$ SD & 21.02 $\pm$ 8.66\\ +Median (min. to max. range) & 22.5 ( 1 - 50 )\\ +\addlinespace + \textbf{Gender} & \\ +\midrule +Male (\%) & 154 ( 41.4 )\\ +Female & 218 ( 58.6 )\\ +\addlinespace + \textbf{Ethnicity} & \\ +\midrule +Caucasian (\%) & 187 ( 50.3 )\\ +African American (Black) (\%) & 4 ( 1.1 )\\ +Asian (\%) & 53 ( 14.2 )\\ +Hispanic/Latino (\%) & 5 ( 1.3 )\\ +Other (\%) & 121 ( 32.5 )\\ +Unknown (\%) & 2 ( 0.5 )\\ +\bottomrule{} +\end{tabular} +\caption{\textbf{Demographic statistics of donors with known vaccine response classification.}}\label{tbl:demoStats} +\end{table} + +The data from the clinical studies consisted of 121 CSV files that were +imported into the FluPrint database. The data was used to build four tables +which will be described in the next sections, but we will not discuss technical +validation of the database construction, refer to the original work for that +\citep{tomicFluPRINTDatasetMultidimensional2019}. The relation between the +tables is best visualised in the original work of +\citep{tomicFluPRINTDatasetMultidimensional2019}, it describes the MySql +attribute types and columns in the tables \autoref{fig:tablesFluprint} +(copied). The volume of the data is also given in the original work, per table +the number of rows and columns is reported \autoref{tbl:volumeTables}. + +\begin{figure}[htpb] + \includegraphics[width=\textwidth]{tablesFluprint} + \caption{ + \textbf{(taken from original paper)} The FluPRINT database model. The diagram shows a schema of the FluPRINT + database. Core tables, donors (red), donor\_visits (yellow), + experimental\_data (blue) and medical\_history (green) are interconnected. + Tables experimental\_data and medical\_history are connected to the core + table donor\_visits. The data fields for each table are listed, including + the name and the type of the data. CHAR and VARCHAR, string data as + characters; INT, numeric data as integers; FLOAT, approximate numeric data + values; DECIMAL, exact numeric data values; DATETIME, temporal data values; + TINYINT, numeric data as integers (range 0–255); BOOLEAN, numeric data with + Boolean values (zero/one). Maximal number of characters allowed in the data + fields is denoted as number in parenthesis. + }\label{fig:tablesFluprint} +\end{figure} + +\begin{table}[htpb] + \centering + \begin{tabular}{lll} + \toprule{} + \textbf{Table name} & \textbf{Rows} & \textbf{Columns} \\ + \midrule{} + \textit{donors} & 740 & 6 \\ + \textit{donor\_visits} & 2,937 & 18 \\ + \textit{experimental\_data} & 371,260 & 9 \\ + \textit{Medical history} & 740 & 18 \\ + \bottomrule{} + \end{tabular} + \caption{Volume of tables in the Fluprint database.}\label{tbl:volumeTables} +\end{table} + +\subsection{Attribute types and values} + +Because of the great number of attributes in the database, we discuss them by +table starting with the donors \autoref{fig:tablesFluprint}. + +\subsubsection{donors table} + +The \textit{donors.id} attribute is simply an enumeration of unique donors, +importantly, it is used as a key to get attributes from other tables. The +column \textit{study\_donor\_id} is an encrypted identification number. Each +donor belongs to the study identified by the \textit{study\_id}, these are the +last two digist of the name code (those starting with SLVP0 \(\cdot\cdot\)) in +the reference table \autoref{tbl:studiesDesc}, the \textit{study\_internal\_id} +is either the digit or a string containing the digit in \textit{study\_id}. The +\textit{gender} and \textit{race} attribute contain the values used in +\autoref{fig:demoGraph}, a minor note is that in the original paper "American +Indian or Alaska Native" is listed as one of the \textit{race} values but is +not used in the database. There are 5 donors whose race is "NULL", which are +mapped to unkown \autoref{fig:demoGraph}. + +\begin{table}[htpb] + \begin{tabular}{rlrlll} +\toprule{} +id & study\_donor\_id & study\_id & study\_internal\_id & gender & race\\ +\midrule{} +1 & e27ad74ff9a5f2f32d8e852533f054c0 & 30 & 30 & Female & Asian\\ +2 & 4a89ac4d3f4dc869e5c8e8cf862cffda & 30 & 30 & Male & Other\\ +3 & a2cde6e54dec92422b0427dd49244350 & 30 & 30 & Female & Caucasian\\ +4 & 0f7d8d1c13e876017ea465f99d25581f & 30 & 30 & Male & Other\\ +5 & 1ed2f6409584b7b4e9720b28d794fe91 & 30 & 30 & Female & Caucasian\\ +\addlinespace +6 & a575678405e9615bfb87eccfa031f7fc & 30 & 30 & Male & Other\\ +\bottomrule{} +\end{tabular} + \caption{Head of the donors table.}\label{tbl:donorsHead} +\end{table} + +\subsubsection{donor\_visits table} + +The donor visits table is the core table of the database, it contains donor +attributes at visit times during enrolment in clinical studies in rows that are +uniquely identified by an \textit{id} integer. Each +row also includes the \textit{donor\_id} identify the donor that visitted. + +The database combines different clinical studies accross years and the data +from these studies is incomplete leading to an incomplete and hetergenous +database \autoref{tbl:visitsDesc}. For example some donors might miss their +second visit to determine their antibody levels, or the number of parameters +measured by an assay changed in the timespan of a clinical study. Unifying +these clinical studies in one database resulted in normalised but incomplete +data and heterogenous data. More specifically, every attribute in the core +table has missing value, which complicates dataset selection. One examples of +visit data of a donor is discussed to highlight important attributes and +problems in the data: that the number of visits is variable, that all columns +are incomplete, and that classification is sometimes based on single visits or +inconsistent \autoref{tbl:visit166} \autoref{tbl:visitsDesc}. + +\begin{table}[htpb] +\addtolength{\leftskip} {-2cm} % increase (absolute) value if needed +\addtolength{\rightskip} {-2cm} % increase (absolute) value if needed +\begin{tabular}{lrrrrrrrrr} +\toprule{} +stat & age & cmv\_status & ebv\_status & bmi & vaccine & geo\_mean & d\_geo\_mean & vaccine\_resp & total\_data\\ +\midrule{} +n & 2937.0 & 1081.0 & 548.0 & 516.0 & 2794.0 & 984.0 & 1260.0 & 1206.0 & 2937.0\\ +na & 0.0 & 1856.0 & 2389.0 & 2421.0 & 143.0 & 1953.0 & 1677.0 & 1731.0 & 0.0\\ +mean & 47.3 & 0.4 & 0.8 & 24.8 & 3.7 & 87.6 & 8.9 & 0.3 & 126.4\\ +sd & 27.0 & 0.5 & 0.4 & 5.6 & 1.0 & 101.7 & 30.9 & 0.4 & 368.4\\ +se\_mean & 0.5 & 0.0 & 0.0 & 0.2 & 0.0 & 3.2 & 0.9 & 0.0 & 6.8\\ +\addlinespace +IQR & 50.2 & 1.0 & 0.0 & 6.7 & 0.0 & 105.4 & 4.0 & 1.0 & 19.0\\ +skewness & 0.2 & 0.3 & -1.4 & 1.0 & -1.7 & 3.6 & 9.9 & 1.1 & 7.1\\ +kurtosis & -1.5 & -1.9 & -0.1 & 2.1 & 3.0 & 26.6 & 114.9 & -0.9 & 49.7\\ +\bottomrule{} +\end{tabular} +\caption{Descriptive stats of relevant numeric or binary factor columns in the + donor visits table. For geo\_mean 0 is considered as missing data.}\label{tbl:visitsDesc} +\end{table} + +Per donor all visits are enumerated in chronological order by +\textit{visit\_id} \autoref{tbl:visit166}. Further visit info includes: +\textit{visit\_internal\_id} which is a number that indicates the visit order +within an influenza season but this differs per clinical study (e.g. some use +1-2-3, orther use 0-7-28), the \textit{vist\_year} is the influenza season of +the visit, the \textit{visit\_day} is the number of days relative to the date +of vaccination, \textit{age} and \textit{age\_round} indicate the donor's age +at time of the visit, and \textit{bmi} gives the donor bmi at visit time, and +lastly \textit{visit\_type\_hai} is the intent of the visit which is either +"pre", "post", or "other", + +During the "pre" visit a virological assay is performed to determine the CMV +and Epstein-Barr virus (EBV) status of the donor, which are indicated by the +binary variables \textit{cmv\_status} and \textit{ebv\_status}. + +To measure vaccine response to a vaccine which is indicated by an id +\autoref{tbl:remapVaccine} in \textit{vaccine}, the hemagglutination inhibition +assay (HAI assay) is used. The procedure measures the influenza antibody titers +before vaccination during the \textit{visit\_type\_hai} "pre" visit of a +participant, and 28 days after vaccination during a "post" visit. The geometric +mean titer (GMT) at each visit is calculated, and a fold change in GMT is +calculated as the ratio of the GMT at day 28 (post) and during the first visit +(pre). These values are \textit{geo\_mean} and \textit{d\_geo\_mean}, +\textit{d\_single} is the antibody titer fold-change per strain of virus used +in the vaccine, it is unclear how this value is aggregated over different +strains and is left out of further analysis. This data was used to classify +donors in high or low responders according to FDA guidelines \cite{}, +individuals are high-responders if they seroconverted (4-fold or greater rise +in HAI titer) and were seroprotected (GMT HAI \(\ge\) 40) after vaccination. +The seasonal vaccine response classifications are given by the binary variable +\textit{vaccine\_resp}. + +The assays performed to get a serological/immunlogical profile of the donor +before vaccination are described later in the section of the experimental data +table, all assays are listed in the original work +\cite{tomicFluPRINTDatasetMultidimensional2019} and are summarised here +\autoref{tbl:assays}, the total rows of assay data is given by +\textit{total\_data}. + +\begin{table}[htpb] +\addtolength{\leftskip} {-2cm} % increase (absolute) value if needed +\addtolength{\rightskip} {-2cm} % increase (absolute) value if needed +\begin{tabular}{rrrlrrrlrrrrr} +\toprule{} +visit\_id & year & day & type & age & cmv & ebv & bmi & vaccine & geo\_mean & d\_geo\_mean & response & assay\_data\_rows\\ +\midrule{} +1 & 2011 & 0 & pre & 20 & 1 & 1 & 30.31 & 4 & 25.20 & 6 & 0 & 343\\ +2 & 2011 & 7 & other & 20 & 1 & 1 & NULL & 4 & 0.00 & 6 & 0 & 51\\ +3 & 2011 & 28 & post & 20 & 1 & 1 & NULL & 4 & 160.00 & 6 & 0 & 51\\ +4 & 2012 & 0 & pre & 21 & 1 & 1 & 30.31 & 4 & 9.28 & 4 & 0 & 292\\ +6 & 2013 & 0 & pre & 22 & 1 & 1 & 30.31 & 4 & 15.91 & 2 & 0 & 2877\\ +\addlinespace +7 & 2013 & 7 & other & 22 & 1 & 1 & NULL & 4 & 0.00 & 2 & 0 & 63\\ +8 & 2013 & 28 & post & 22 & 1 & 1 & NULL & 4 & 26.75 & 2 & 0 & 82\\ +\bottomrule{} +\end{tabular} +\caption{Visit data of donor 166 from study SLVP021 \autoref{tbl:studiesDesc}, +where participants are only vaccinated once. +Number of visits and data collected at visit varies, classification is +inconsistent with \( \geq 40\) and 4-fold increase +rule in 2011.}\label{tbl:visit166} +\end{table} + +The most important data related to the visits of donor 166 is shown in Table +\ref{tbl:visit166}. The vaccine response classification is calculated based on +the GMT in the "pre" and "post" visits. This classification is done per +influenza season, but the HAI assay requires a "pre" visit and a "post" visit +28 days later to measure the difference in GMT. However, sometimes a +classification is given when there is only one visit record in a season, like +in 2012 for donor 166 \autoref{tbl:visit166}. + +\begin{figure}[htpb] + \includegraphics[width=\textwidth]{season_classification} + \caption{}\label{fig:classInconsistent} +\end{figure} + +The example of donor 166 contains an inconsistency in the classification, in +2011 the GMT \textit{geo\_mean} increases from 25.20 to 160.00, and the +\textit{d\_geo\_mean} is 6, but in this season the donor is wrongly classified +as a low responder \autoref{tbl:visit166}. Because of this the seasonal +classification of donors was investigated using the seroprotection and +seroconversion criteria \ref{fig:seasonalClasses}, records of incorrectly +labelled donors are also saved as a spreadsheet. This data is inconsistent in +the database, but the most likely explanation is that antibody titer for one +strain of virus did not meet the high response classification criteria. In this +work it is considered as inconsistent because individual strain titer data is +not in the database, but classification is therefore not necessarily incorrect. +Hence the classification will be used in this work without further selection. + +\subsubsection{Experimental data table} + +\begin{table}[htpb] + \begin{tabularx}{\textwidth}{Xp{0.5\textwidth}X} +\toprule{} + \textbf{Name} & \textbf{Description} & \textbf{id} (\textit{experimental\_data.assay})\\ +\midrule{} + (Multiplex) cytokine assays & Multiplex ELISA using Luminex polysterene + bead or magnetic bead kits. Measures serum cytokine/hormone level in + z.log2 units using fluorescent antibodies. & 3, 6, 15, 16\\ + \addlinespace + Flow and mass cytometry assays & uses labeled antibodies to detect antigens on + a cell surface to identify a subset of a cell population, units are in + percentage of parent population. & 4, 9, 13, 17 \\ + \addlinespace + Phosphorylation cytometry assays & Uses antibodies to measure + phosphorylation of specific proteins stimulated by an immune system event + belonging to cell population subsets. Units are a fold change between + stimulated and un-stimulated cells, for mass cytometry arcsin readout difference, + fold-change of 90th percentile readout values otherwise. & 7, 10 (mass cytometry) (flow cytometry)\\ + \addlinespace + complete blood count (CBCD) & Different cells are counted using flow + cytometry Units are usually in Count/$\mu$L & 11 \\ + \addlinespace + meso scale discovery assays (MSD) & A setup where serum cytokines or hormones + are captured with antibodies, and then detected by using a detection + antibody. Units are arbitrary intensity & 2, 12, 14 \\ +\bottomrule{} +\end{tabularx} + \caption{assays table}\label{tbl:assays} +\end{table} + + +Assays performed in visits are remapped, but the values in the +database do not correspond to the reported table \autoref{tbl:remapVaccine}. +Actual assay type, data units, and id in the database are reported here +\autoref{tbl:assays}. + +\fpfig{exp_data_numbers}{.7} +{Feature count per individual assay id, assay type, stratisfied in either response status or study} +{caption} +{fig:featureNumbers} + +In total there are data from 14 different assays, not counting the virological +and HAI antibody assays \autoref{tbl:assays}. The virological assays include +the cmv virus status and ebv status, and is not used in this work because it is +done in a smaller subset of studies. Those 14 assays have been aggregated in +this work to 5 different types of experiments: the multiplex assays measure +serum molecules such as cytokines and other signaling molecules, flow and mass +cell cytometry measure the phenotype of specific immune related cells, +phosphorylation flow and mass cytometry measures the phosphorylation signaling +pathway activation after an immune stimulation, the blood count measures the +count of cells in the blood, and meso scale discovery (MSD) measures hormones +or cytokines from the blood. + +\begin{figure}[htpb] + \includegraphics[width=\textwidth]{assay_value_distributions} + \caption{noise in 90th \%tile}\label{fig:assayDistr} +\end{figure} + +The experimental data table contains all features recorded for a donor visit. +The number of features collected for each visit is large and varies greatly +(mean at 126 , \(\pm \)368 SD) \autoref{tbl:visitsDesc}, and in total there are +3285 different features measured across all clinical studies. However, not +every assay is done in every clinical study \autoref{fig:featureNumbers} and +over the years the data generated by assays has changed, so a table with all +features as columns and all donors as rows would be extremely sparse (and +crashes R due to RAM limitations). Describing the 3285 different features in +this sparse table would be impossible, but assay value distributions across +studies are shown to follow normal or power distributions +\autoref{fig:assayDistr}. The features included 102 blood-derived immune cell +subsets analyzed by mass cytometry. It also included the signaling capacity of +over 30 immune cells subsets stimulated with seven conditions, which were +evaluated by measuring the phosphorylation of nine proteins. Additionally, up +to 50 serum analytes were evaluated using Luminex bead arrays +\citep{tomicSIMONAutomatedMachine2019}. + +No correlation analysis was done, since this is complicated by the great number +of features and sparseness in the data. + +\begin{figure}[htpb] + \includegraphics[width=\textwidth]{repeat_visits_per_study} + \caption{the number of donors that visited per number of influenza seasons + they visited (years), per study. The color indicates the number of visits for which a + classification was available, counted within the groups of donors that + visited the same amount of times.}\label{fig:repeatVisits} +\end{figure} + +What further complicates selecting data is repeat visits of donors, and missing +visits. The problem of repeat visits over a span of multiple influenza seasons +is that not the same assay types are done, and that repeat visits are only a +small portion of the database. The data is also not suitable right away for +studying the effect of repeat vaccination on high versus low vaccine reponse, +since the classification in the longitudal study (SLVP015) is mostly not +available \autoref{fig:repeatVisits}. + +For example exploring the effect repeat vaccination has an response rate would +first require manual labelling of high and low responses, at least for the +cases where it is possible based on the GMT data. Those cases are when +classification is set to a null value even though GMT data is available. The +reason for this null value assignment is reported, but the pattern seems to set +the vaccine response to null if there is not enough assay data measured. + +\subsection{Data quality} + +The database has issues that are inherent to combining multiple studies and the +classification is inconsistent in some cases \autoref{fig:classInconsistent}, +or often missing completely because no HAI antibody assay data was available or +the classification was set to a null value by the database authors because +possibly the antibody titer for a single strain of virus in the vaccine was too +low (this data is not in the database) \autoref{fig:repeatVisits}. The main value of +the database is the assay data that is fully represented in all studies and +across all years, but this information is hard to access since all studies do +not use overlapping assays \autoref{fig:featureNumbers}, resulting in high +sparsity data. Further, the sample size that can be used for further studies is +limitted, since the high versus low vaccine response is only available for a +small subset of the data. + +Specific attributes that have great amounts of missing values are the +virological and HAI assay data, the last is used for the vaccine response +classifcation. Potential for Studying the correlation of these values with +vaccine response is thus limitted. Nevertheless assay data is often available +and could be used to identify immunological factors that correlate with other +data, such as repeat vaccination, the exploration of this effect is outside the +scope of this work due to the data sparsity issues. + +\section{Data preparation} + +\subsection{Data selection} + +The data selection used in this work is based the query used in \spaper +\autoref{lst:QueryTemplate}. Using this query generates a subset of \flup +comprised data from 5 clinical studies, most importantly the longitudal study +SLVP015 \autoref{tbl:studiesDesc}. Presumably, the authors of \spaper included +only the first visit of donors because the classification is the most complete +in this dataset \autoref{fig:dataRepeatvisits}. In this work we use this query +to generate initial first visit datasets. However, to explore repeat vaccination, +we select a subset of this data that includes donors with a repeat visit in a +second influenza season. + +\begin{figure}[htpb] + \includegraphics[width=\textwidth]{data_selection} + \caption{caption}\label{fig:dataRepeatvisits} +\end{figure} + +The initial query used in \spaper generated a long table for in total 3285 +different features recorded at the first visit of 195 donors in different +studies and years (referred to as \firstvis). The observable pattern in this +data is that low responding donors are overrepresented and that the +classification is consistent with the log2 \gmt change +\autorefsub{fig:repeatVisits}{A}. Unfortunately, the number of donors in the +\firstvis that returned in other influenza seasons decreases quickly, limiting +possibilities of comparing models built on \firstvis and subsequent visit data +\autorefsub{fig:repeatVisits}{B}. Nevertheless, we selected the \secondvis to +explore repeat vaccinations. The \secondvis has an exacerbated class imbalance +that precludes training any models, therefore we use the \secondvis to explore the +knowledge gained by models trained on the \firstvis +\autorefsub{fig:repeatVisits}{C}. + +\begin{minipage}{\linewidth} + +\begin{lstlisting}[caption=Applying the mulset algorithm and preparing the data, label={lst:mulsetStep}] +generate intersection datasets suitable for analysis +for {each donor in data} do: + Calculate intersection between donor and all other donors using mulset algorithm + Skip sets that have less than 5 features and less than 15 donors in common +end for; + +for {each set in generated datasets} do: + Partition data in training (75\%) and test (25\%) split + Skip sets that have less than 10 donors in the test set +end for; + +for {each set in prepared datasets} do: + calculate number of donors that visited a second influenza season + skip dataset if it is not in the top3 of datasets containing the highest number of second visitors +end for; +\end{lstlisting} +\end{minipage} + +The \firstvis had a total of 640575 cells of which 596736 values were +missing (sparsity of 93\%) because of the heterogeneity in clinical studies and +years where data was collected. In \spaper missing data is not imputed because +there is not enough prior knowledge. And, since every donor had a missing +feature, dropping all rows/donors was not an option either. A solution used in +\spaper was generating complete tables comprising subsets of donors that had +all features in common using the mulset algorithm \autoref{lst:mulsetStep} +\autoref{fig:mulsetAlg}. + +In this work the procedure in \spaper was replicated and extended to generate +usable datasets \autoref{lst:mulsetStep}. Firstly, there were duplicate +measurements of features in the \firstvis, these were aggregated to unique +feature records using the mean. Second, the mulset R package was used to +generate 47 complete datasets. These datasets were then reduced to 36 by +selecting those that had at least 5 features and 15 donors. Finally, the +datasets were split into train (75\%) and test (25\%) sets, and datasets with +less than 10 donors in the test set were discarded reducing the number of +datasets to 20 \autoref{tbl:mulsetDatasets}. + +\begin{table}[htpb] + \begin{tabularx}{\textwidth}{XXXXX} +\toprule{} + dataset & Rows x Cols (Donors x Features) & total (low / high (low \%)) & train (low / high) & test (low / high)\\ +\midrule{} +1 & 61 x 78 & 43 / 18 ( 0.7 ) & 33 / 14 & 10 / 4\\ +2 & 105 x 101 & 62 / 43 ( 0.59 ) & 47 / 33 & 15 / 10\\ +3 & 140 x 50 & 94 / 46 ( 0.67 ) & 71 / 35 & 23 / 11\\ +4 & 63 x 269 & 38 / 25 ( 0.6 ) & 29 / 19 & 9 / 6\\ +5 & 62 x 293 & 38 / 24 ( 0.61 ) & 29 / 18 & 9 / 6\\ +\addlinespace +6 & 68 x 237 & 42 / 26 ( 0.62 ) & 32 / 20 & 10 / 6\\ +7 & 67 x 44 & 47 / 20 ( 0.7 ) & 36 / 15 & 11 / 5\\ +8 & 111 x 93 & 66 / 45 ( 0.59 ) & 50 / 34 & 16 / 11\\ +9 & 73 x 54 & 58 / 15 ( 0.79 ) & 44 / 12 & 14 / 3\\ +10 & 40 x 105 & 28 / 12 ( 0.7 ) & 21 / 9 & 7 / 3\\ +\addlinespace +11 & 46 x 97 & 32 / 14 ( 0.7 ) & 24 / 11 & 8 / 3\\ +12 & 137 x 53 & 78 / 59 ( 0.57 ) & 59 / 45 & 19 / 14\\ +13 & 48 x 42 & 35 / 13 ( 0.73 ) & 27 / 10 & 8 / 3\\ +\textbf{14} & \textbf{91 x 38} & \textbf{62 / 29 ( 0.68 )} & \textbf{47 / 22} & \textbf{15 / 7}\\ +15 & 42 x 37 & 36 / 6 ( 0.86 ) & 27 / 5 & 9 / 1\\ +\addlinespace +\textbf{16} & \textbf{92 x 26} & \textbf{62 / 30 ( 0.67 )} & \textbf{47 / 23} & \textbf{15 / 7}\\ +17 & 88 x 6 & 68 / 20 ( 0.77 ) & 51 / 15 & 17 / 5\\ +18 & 82 x 87 & 56 / 26 ( 0.68 ) & 42 / 20 & 14 / 6\\ +\textbf{19} & \textbf{151 x 51} & \textbf{92 / 59 ( 0.61 )} & \textbf{69 / 45} & \textbf{23 / 14}\\ +20 & 83 x 75 & 56 / 27 ( 0.67 ) & 42 / 21 & 14 / 6\\ +\bottomrule{} +\end{tabularx} + \caption{Datasets generated by applying the mulset algorithm on the \simon + \firstvis, and the balanced train test split that was performed.}\label{tbl:mulsetDatasets} +\end{table} + +A significant number of datasets contained more predictors than samples +\autoref{tbl:mulsetDatasets}. However, we consider this as inevitible and not an +absolute obstacle since the purpose of the models is not to discriminate +vaccine responders with the highest accuracy, but to identify features that +correlate with a vaccine response from the great number of features. + +In this work we select the datasets that best fit to the business objectives of +exploring repeat vaccination effects, as well as finding features that +correlate with vaccine responses. Accordingly, we calculated the number of +donors that visited a second influenza season per dataset and chose the top3 +datasets. The remaining datasets were 14, 16, and 19 +\autorefsub{tbl:mulsetDatasets}{\textbf{bold rows}}. These datasets had, +respectively, 27 out of 91, 27 out of 92, and 21 out of 151 donors that +returned for a second vaccination. Alarmingly, this was less than half of the +dataset in all cases, and for other datasets this number was even smaller. + +Within all three datasets 82 of the donors are shared indicating that using +both dataset 14 and 16 might add little additional information. Furthermore, 26 +of the measured features are shared between dataset 14 and 16, meaning all +features of dataset 16 are in dataset 14. Further, all features are also +phospho flow assay data \autoref{tbl:assays}. Nevertheless, in this work we +include both datasets for modeling. + +The \secondvis corresonding to the selected \firstvis in the chosen datasets +were retrieved during exploration of repeat vaccinations using the modeling +results. + +\subsection{Data cleaning} + +In this work features and rows were not changed for the chosen datasets, since +this would result in a lower number of rows when already limitted data is +suitable for modeling. Furthermore, the objective of this work is not +obtaining optimal models but exploring repeat vaccination and vaccine responses. + +\subsection{Data formatting} + +The final format of the datasets were complete tibbles containing the outcome +and features as columns, and donors as rows. + +\section{Modelling} + +\subsection{Choice of modeling technique} + +In this work a form of wrapper feature selection is used, since we are training +models on different subsets of features and chose those that discriminate the +best between low and high vaccine responders +\citep{hiraReviewFeatureSelection2015}. Altough, technically the aim is to +train an at least fair discriminator on any suitable dataset to then use that +model to identify new knowledge about vaccine response and repeat vaccination. + +Four models were chosen for this task: the naive bayes classifier (nb), the +random forest model (rf), the regularised logistic regression model (reglog), +and regularised linear discriminant analysis (rrlda). In \spaper an automatic +machine learning pipeline is used where 2400 models are trained on all 20 +datasets, and the best models are then used to explore important features that +correlate with a high vaccine response. This approach is out of scope for this +work, and instead we change the objective to specifically identifying repeat +vaccination effects. Additionally, the datasets chosen in this work are not +discussed in \spaper. + +\subsection{Test design} + +The three selected datasets were already split in test and training sets. The +training set was used for training models using 2 times repeated 10 fold +cross-validation where the accuracy was computed on every fold. The models that +had the best cross-validated accuracy were compared using the training and test +area under the curve measure, since we are interested in general discriminative +ability. Using these measures the best discriminator is chosen for further +exploration of repeat vaccination and vaccine response features. + +\subsection{Model parameters and assessment} + +\begin{table}[htpb] +\addtolength{\leftskip} {-2cm} % increase (absolute) value if needed +\addtolength{\rightskip} {-2cm} % increase (absolute) value if needed +\begin{tabular}{llrrrrrrrrrrrrr} +\toprule{} +dataset &model & SENS & SPEC & MCC & PREC & NPV & FPR & F1 & TP & FP & TN & FN & train AUC & test AUC\\ +\midrule{} +14 & rrlda & 0.091 & 0.915 & 0.010 & 0.333 & 0.683 & 0.085 & 0.143 & 2 & 4 & 43 & 20 & 0.50 & 0.62\\ + & nb & 0.636 & 0.702 & 0.321 & 0.500 & 0.805 & 0.298 & 0.560 & 14 & 14 & 33 & 8 & 0.67 & 0.59\\ + & rf & 0.364 & 0.851 & 0.243 & 0.533 & 0.741 & 0.149 & 0.432 & 8 & 7 & 40 & 14 & 0.65 & 0.61\\ + & reglog & 0.227 & 0.766 & -0.007 & 0.312 & 0.679 & 0.234 & 0.263 & 5 & 11 & 36 & 17 & 0.49 & 0.48\\ +\addlinespace +16 & rrlda & 0.000 & 1.000 & NaN & NaN & 0.671 & 0.000 & 0.000 & 0 & 0 & 47 & 23 & 0.48 & 0.61\\ + & nb & 0.652 & 0.617 & 0.253 & 0.455 & 0.784 & 0.383 & 0.536 & 15 & 18 & 29 & 8 & 0.68 & 0.55\\ + & rf & 0.261 & 0.851 & 0.135 & 0.462 & 0.702 & 0.149 & 0.333 & 6 & 7 & 40 & 17 & 0.65 & 0.69\\ + & reglog & 0.391 & 0.723 & 0.116 & 0.409 & 0.708 & 0.277 & 0.400 & 9 & 13 & 34 & 14 & 0.64 & 0.47\\ +\addlinespace +19 & rrlda & 0.533 & 0.391 & -0.075 & 0.364 & 0.562 & 0.609 & 0.432 & 24 & 42 & 27 & 21 & 0.47 & 0.41\\ + & nb & 0.489 & 0.565 & 0.053 & 0.423 & 0.629 & 0.435 & 0.454 & 22 & 30 & 39 & 23 & 0.54 & 0.48\\ + & rf & 0.244 & 0.739 & -0.018 & 0.379 & 0.600 & 0.261 & 0.297 & 11 & 18 & 51 & 34 & 0.54 & 0.52\\ + & reglog & 0.267 & 0.754 & 0.023 & 0.414 & 0.612 & 0.246 & 0.324 & 12 & 17 & 52 & 33 & 0.51 & 0.32\\ +\bottomrule{} +\end{tabular} +\caption{Model evaluation measures on the three chosen datasets. + SENS=sensitivity: proportion of true positives, SPEC=specificity: + proportion of true negatives, MCC=mathews correlation coefficient: + correlation prediction with true labels, PREC=precision: true positive over + predicted positive ratio, NPV=negative predictive value: true negative over + predicted negative ratio, F1=f1-score: harmonic mean precision and + accuracy, TP: true positives, FP: false positives, TN: true negatives, FN: + false negatives, AUC: area under the receiver operator curve.}\label{tbl:modelEval} +\end{table} + +On all three datasets the model with the highest train and test AUC metric was +the naive bayes classifier \autoref{tbl:modelEval}. On dataset 14 and 16 the +naive bayes model reached a training AUC of 0.67-0.68, which could reflect the fact +that these dataset share a large part of donors and features. An AUC value in +this range is considered to be a (somewhat) fair discriminator. Although, +ideally discriminators would have training and test AUC values in the range 0.7 +and up, anything below is considered a weak discriminator +\citep{L_demann_2006}. On dataset 19 all models failed to produce good +discriminators, hence we discard this dataset from further analysis. + +On dataset 14 and 16 the random forest model had similar performance compared +with the naive bayes model, the training AUC score was only slightly lower and +the model performed better on unseen test data. This could indicate that the +random forest model is overfitting the training data less than the naive bayes +model, and would therefore be the preferred choice when choosing a +discriminator to be used for new data. Despite this, in this work we consider +the naive bayes model the best on dataset 14 and 16. Further, we continue the +exploration of vaccine responses and repeat vaccination only using the naive +bayes models on dataset 14 and 16. This is motivated by the fact that we are +not interested in the best model and the random forest model tends to predict +false negatives (sensitivity of 0.364) \autoref{tbl:modelEval}, this last fact +is the most problematic since the negative class is overrepresented in our data. + +The parameters for both naive bayes models were laplace = 0 and usekernel = +TRUE and adjust = 1. + +\section{Exploration of modeling results} + +Using the models built on dataset 14 and 16 our goal was to identify the +features relevant to the generation of antibodies in response to vaccination. +The procedure is the same as in \spaper, we calculate the feature importance +for the classifier model and rank them based on their contribution to the model +from 0 to 100. The top three features with the highest score are explored in +more detail. Furthermore, for these features we look at the measurements of +these features in the \secondvis to explore the effect of repeat vaccination. +Lastly, we also calculated the correlation of all features in dataset 14 and 16 +to identify feature groups related to the top three most important features. + +\subsection{Identifying phospho flow cytometry cell signaling features correlated with vaccine response} + +\begin{figure}[htpb] + \centering + \includegraphics[width=\textwidth]{dataset1_nb_feature_exploration} + \caption{dataset1-nb-feature-exploration} + \label{fig:dataset1-nb-feature-exploration} +\end{figure} + +Firstly, the top ranked feature in dataset 14 was the phosphorylated STAT5 transcription factor in unstimulated B cells \autorefsub{fig:dataset1-nb-feature-exploration}{A}. However, the difference in the value of this feature between the high and low vaccine responders was not found to be significant at FDR $<$ 0.01 \autorefsub{fig:dataset2-nb-feature-exploration}{B}. In contrast, the other two features, IFNg stimulated B-cell phosphorylated STAT5 and CD4 T cell phosphorylated STAT5, were found to be significantly greater in the high responder group (FDR $<$ 0.01). A correlation analysis of all features showed that different STAT protein formed positively correlated clusters as expected \autoref{fig:cor-dataset1} (p < 0.0001). Further, the most important feature had slight negative correlations (pearson's r from -0.2 to -0.5) to a set of stimulated STAT1 cell responses (p < 0.0001 after BH adjustment). The second most important feature has similar correlations as the first, likely since they are both B-cell STAT5 features. Lastly, the unstimulated CD4 positive phenotype T-cells STAT5 phosphorylation also belonged in the same cluster as the previous B-cell features. These correlations might indicate an interaction between the STAT5 and STAT1 phosphorylation in response to a vaccine. + +\begin{figure}[htpb] + \centering + \includegraphics[width=\textwidth]{dataset2_nb_feature_exploration} + \caption{dataset2-nb-feature-exploration} + \label{fig:dataset2-nb-feature-exploration} +\end{figure} + +Secondly, in dataset 16 there were only four features that had a variable importance score greater than 50 \autorefsub{fig:dataset2-nb-feature-exploration}{A}. The top two features were phospohorylated STAT1 in unstimulated B-cells and phosphorylated STAT1 in unstimulated CD8 T-cells. However, only the B-cell feature was found to be significantly greater in the positive class (FDR \(< 0.01\)) \autorefsub{fig:dataset2-nb-feature-exploration}{B}. The B-cell STAT1 feature correlated positively with both unstimulated CD8 and CD4 STAT1 phosphorylation (pearson's r= 0.7 and 0.4, p \(< 0.001\)), and there were mild negative correlations with interferon gamma stimulated monocyte STAT3 and STAT5 phosphorylation (pearson's r= 0.3 and 0.2, p \(< 0.001\)) \autoref{fig:cor-dataset2}. + +\subsection{Repeat vaccination effect on identified features} + +Firstly, in the \secondvis of both datasets there were outliers (donors had a value greater than 1000) and negative values. In this work these were left out, since outliers made the pattern unclear and the negative values were considered as nonsensical values. + +\begin{figure}[htpb] + \centering + \includegraphics[width=\textwidth]{second_visit_change1} + \caption{second-visit-change1} + \label{fig:second-visit-change1} +\end{figure} + +To see how a repeat vaccination affects immune cell signaling, the distribution of the top three features of dataset 14 were compared to their distribution when measured in a subsequent influenza season \autorefsub{fig:dataset1-nb-feature-exploration}{C}. +In the 21 donors that had a second measurement of the features in another influenza season that were not left out (outliers and nonsensical values) there was the consistent pattern that the high responders were classified as low responders in their second visit \autorefsub{fig:dataset1-nb-feature-exploration}{C}. +Although, overall the feature values were consistently greater in the \secondvis \autorefsub{fig:dataset1-nb-feature-exploration}{C, enlarged diamonds}. +Thus, vaccination might increase activity in general signaling pathways of PBMC in subsequent influenza seasons, but the classification does not reflect this as increasing influenza antibody response. +One possibility is that the donor was classified as low responder due to a lack of response to one strain of virus in the vaccine administered in the repeat visit, not necessarily to all strains \autoref{fig:classInconsistent}. + +To explore the overall change in the features of dataset 14 between the first and subsequent in influenza seasons the distribution of changes for donors were visualised and ordered by mean of log2 change (negative values were removed) \autoref{fig:second-visit-change1}. +The overall trend that appeared was that the unstimulated PBMCs had higher values upon a repeated visit. +And, in general STAT5 features increased in value. The values that contributed the most to the model discriminating between high and low responders in the \firstvis also increased the most in a repeat visit. +Although, there are outliers that increased a lot in the subsequent influenza season \autoref{fig:second-visit-change1}. + +On dataset 16 two of the top three features had similar distributions to the \firstvis \autorefsub{fig:dataset2-nb-feature-exploration}{C}. +In contrast, unstimulated monocyte cells had higher STAT5 phoshporylation in the subsequent influenza season \autorefsub{fig:dataset2-nb-feature-exploration}{C}. +Further, the same three donors that were classified as high responders in the \firstvis and as low responders in the \secondvis as in dataset 14 \autorefsub{fig:dataset1-nb-feature-exploration}{C} had increased monocyte cell STAT5 phosphorylation \autorefsub{fig:dataset2-nb-feature-exploration}{C, enlarged diamonds}. +Lastly, the top three features of the model trained on dataset 14 also belonged to those that increased the most between the \firstvis and \secondvis \autoref{fig:second-visit-change1}. + +\section{Discussion and conclusion} + +In this work we gave a brief introduction into influenza vaccination and how vaccine responses are measured, described the \flup database, and applied a similar data mining method as in \spaper and additionally explored the available repeat vaccination data. +The \flup database made it possible to study vaccine responses by providing a classification of donors into high or low responders based on measured antibody level before and after vaccination. +Further, it combined and preprocessed data from multiple clinical studies in an accessible database format. +This resulted in a wide variety of data on immune cell populations, serum signaling molecules, and cell signaling activity that is suitable for studying immune correlates to vaccine responses using data mining method. +We applied a procedure as described by the authors of \flup in \spaper, wrapper feature selection using multiple models trained on interesting data subsets of \flup. Using this procedure we then explored selected features and how they changed in subsequent influenza seasons. +It was found that STAT5 related signaling features correlated with a vaccine response and increased the greatest amount in subsequent influenza seasons. + +Initially, the idea was to focus on building accurate predictors of vaccine response by training models including constructed features based on repeat vaccination. However, during the data understanding phase of this project it became clear that \flup contains only complete classifications in the \firstvis. +Instead, the objective was revised to explore the available data on repeat vaccination using models trained on \firstvis data from a selection of clinical studies that received the same vaccine, as done in \spaper. +Overall, during the data understanding phase it became clear that \flup is not suitable for predicting vaccine response with high accuracy, since data is combined from multiple studies and years. +This means using donors/rows from different years and studies creates highly sparse predictors. +Consequently, using \flup data requires selecting small datasets without missing values, this only slightly increases the available example measurements of features by combining data from different studies. +Further, the available data on repeat vaccinations is limited to mostly one clinical study, and in repeat visits there is often no classification making it impossible to train models using repeat vaccination data. + +During the data understanding phase we also found that classification is missing in a lot of cases. +Further, we identified an inconsistency in the classification data presented in \flup. +However, this is likely due to the fact that the before and after antibody titer against individual influenza strains in the vaccine is not completely available in the database and not because the classification is incorrect. +Thus to check the classification quality it is necessary to study the raw data and scripts used to generate the database, which is considered out of the scope of this work. + +The data preparation and modeling phases included selecting the data that was most suitable for training models and studying repeat vaccinations. +We started with the initial data used in \spaper and also collected repeat vaccination data for the donors in this dataset. +To deal with the sparse data the mulset algorithm was applied to generate twenty small but complete datasets, the three datasets that had the highest amount of donors that received a repeat vaccination were then chosen for modeling and further analysis. +Four models were built all three datasets, but models with fair discriminative ability were built only on dataset 14 and 16. + +The features in dataset 14 and 16 were all from the phospho-flow cytometry phosphorylation assay, from them we used the models to identify features correlated with a high vaccine response. +We found that STAT5 phosphorylation in immune cells from different lineages was associated with a high vaccine response and was increased in subsequent influenza seasons. +However, further study of this result is considered out of the scope of this work where the focus lies on the application of data science tools. +Instead, we show here that data mining methods described in \spaper can be replicated to answer research questions using complex clinical datasets. + +The objectives defined before selecting the data and starting the data preparation and modeling phase of the project were: +\begin{itemize} + \item What kind of studies can be done using the \flup database? + \item What immunological factors correlate to a vaccine responses? + \item What is the effect of repeat vaccination? +\end{itemize} + +In summary, we provided insight into which studies can be done using the \flup database by describing the experimental data tables of \flup. +It became clear that \flup is suitable for correlating immunological features with a vaccine response by selecting small complete datasets, but that the possibility of combining large data across years and different studies is limited in \flup. +Additionally, we found that classifications are not available in a great amount of data points limiting the sample size for classification studies. +Further, we identified a group of immune cells from different lineages that had increased phosphorylation activity correlated to vaccine response and found that this increase was present in subsequent influenza seasons. + +\section{Materials and methods} + +\subsection{Data collection} + +By following the guide on the \href{https://github.com/LogIN-/fluprint}{FluPrint Github Repository} the MySQL +server was set up. +All file paths mentioned refer to the github repository of this project which can be found below. + +In this work the FluPrint github was first added as a submodule. +This module provides the php scripts to import raw data csv's into the MySQL database. +The operating system and versions of php and MySQL used in this work were OSX "Big Sur" (on Mac Book air 2017), php 7.3.24 (built-in mac version), and MySQL 8.0.23 (homebrew). + +In the \href{https://github.com/LogIN-/fluprint}{guide} the dependencies to run +the php import script were installed first. This was also done in this work, +except that the hash-file verification step was skipped. + +After the php dependencies were installed the MySQL server was started. By +default homebrew recommends to use the \lstinline{homebrew services [option] [SERVICE]} command to start the MySQL server. However, in this work the server +is started using \lstinline{mysql.server start} which provides a socket that +was symlinked using \lstinline{sudo ln -s /tmp/mysql.sock /var/mysql/mysql.sock}. This was done to prevent an error +(\href{https://stackoverflow.com/questions/15016376/cant-connect-to-local-mysql-server-through-socket-homebrew/18090173}{StackOverflow: cant connect to local mysql server through socket homebrew}) thrown +by the php import scripts. Before the import scripts were run a user was added to the +MySQL server and a database was created \ref{lst:addUser}, the password type had to be \lstinline{mysql_native_password} +(\href{https://stackoverflow.com/questions/62873680/how-to-resolve-sqlstatehy000-2054-the-server-requested-authentication-metho}{how to resolve [SQLSTATEHY000] 2054 the server requested authentication method.}). + +\begin{lstlisting}[language=sql, caption=Adding user and database to sql server, label={lst:addUser}] +mysql> CREATE USER 'mike'@'localhost' IDENTIFIED BY ';lkj'; +mysql> GRANT ALL PRIVILEGES ON * . * TO 'mike'@'localhost'; +mysql> ALTER USER 'mike'@'localhost' IDENTIFIED WITH mysql_native_password BY 'mike'; +mysql> CREATE DATABASE fluprint; +\end{lstlisting} + +The databasename, the username, and password were added to the +\lstinline{config/configuration.json} of the FlruPrint github module. At this +point the configuration for the php import scripts was finished, and the raw +data downloaded in \lstinline{data/upload} were imported in the MySQL server +using \lstinline{php bin/import.php}. + +\subsection{Statistical methods} + +\subsubsection{Data selection} + +In this work immunological features correlating to a vaccine response were identified using wrapper based feature selection on data from the \flup SQL database. +Suitable datasets without missing values were generated using the \href{https://cran.r-project.org/web/packages/mulset/index.html}{R package mulset}, as described in the data preparation section. +These datasets were split into training and test splits using the createDataPartition function from the R package \href{https://topepo.github.io/caret/}{caret}. +As described in the data preparation and selection sections, datasets were not considered if the test set had less than 10 donors. +Lastly, from the generated datasets the number of donors in the \secondvis was used to choose datasets for further analysis. The \secondvis data was obtained from the database by a query that is avalaible in the github repository of this project. + +\subsubsection{Model training, evaluation, exploration} + +Standard procedure were used for model training, models were trained only on the training datasets using 10-fold cross-validation that was repeated two times. +The test data was used only as an independent dataset to estimate how much the model overfits on the training data. +Model training itself was done using the \href{https://topepo.github.io/caret/}{caret} R package function train. +Additionally, parameters were chosen based on the highest cross-validated accuracy automatically train function. + +Variable importance of the models generated by the \href{https://topepo.github.io/caret/}{caret} package was generated by the function varImp from the same package. +This uses model specific feature contribution statistics and ranks them from most important to not important on a scale from 0 to 100, for example for the naive bayes model it uses the class conditional probabilities of features. + +Confusion matrix metrics were generated using the \href{https://cran.r-project.org/web/packages/MLeval/index.html}{MLeval} R package which accepts caret objects and computes metrics in a table format as shown in the model evaluation section. +Additionally, the test AUC was calculated with another R packages called \href{https://cran.r-project.org/web/packages/pROC/pROC.pdf}{pROC}. + +Correlation plots of the features from the selected datasets 14 and 16 were made using the R package \href{https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html}{corrplot}. + +\subsubsection{Significance tests} + +To see if features identified by the best classifier trained on datasets 14 and 16 had different distribution in between the two classes the significance analysis of micro arrays (SAM) at a FDR \(<\) 0.01 was used in R. +P values for all correlations shown in the correlation plots below were calculated using an R package, and only correlations with a p-value less than 0.001 were shown. + +\subsection{Code and data availability} + +The code and data belonging to this project can be found in the \href{}{github repository}. + +\printbibliography + +\begin{appendices} + + \section{Correlation plots} + +\begin{figure}[htpb] + \centering + \includegraphics[width=\textwidth]{cor_dataset1} + \caption{cor-dataset1} + \label{fig:cor-dataset1} +\end{figure} + +\begin{figure}[htpb] + \centering + \includegraphics[width=\textwidth]{cor_dataset2} + \caption{cor-dataset2} + \label{fig:cor-dataset2} +\end{figure} + + \section{mulset algorithm} + +\begin{figure}[htpb] + \includegraphics[width=\textwidth]{F2.large} + \caption{\textbf{taken from original work}}\label{fig:mulsetAlg} +\end{figure} + +\begin{table}[htpb] +\addtolength{\leftskip} {-2cm} % increase (absolute) value if needed +\addtolength{\rightskip} {-2cm} % increase (absolute) value if needed +\begin{tabular}{rrrrrlrllrrl} +\toprule{} +donor\_id & study & age & outcome & year & type & hai\_response & name & data\_name & assay & data & dup\\ +\midrule{} +285 & 18 & 9.47 & 0 & 2009 & pre & 1 & CD4+ T cells & CD4\_pos\_T\_cells & 13 & 33.8 & TRUE\\ +285 & 18 & 9.47 & 0 & 2009 & pre & 1 & CD4+ T cells & CD4\_pos\_T\_cells & 13 & 34.1 & TRUE\\ +285 & 18 & 9.47 & 0 & 2009 & pre & 1 & CD4+ T cells & CD4\_pos\_T\_cells & 13 & 34.3 & TRUE\\ +285 & 18 & 9.47 & 0 & 2009 & pre & 1 & CD4+ T cells & CD4\_pos\_T\_cells & 13 & 33.0 & TRUE\\ +\bottomrule{} +\end{tabular} + \caption{}\label{tbl:exampleDuplicate} +\end{table} + + + + \section{Query that generates initial \simon data} +\begin{lstlisting}[language=sql, caption=Query of initial SIMON data, label={lst:QueryTemplate}] +SELECT donors.id AS donor_id, + donor_visits.age AS age, + donor_visits.vaccine_resp AS outcome, + experimental_data.name_formatted AS data_name, + experimental_data.data AS data +FROM donors + LEFT JOIN donor_visits + ON donors.id = donor_visits.donor_id + AND donor_visits.visit_id = 1 + INNER JOIN experimental_data + ON donor_visits.id = experimental_data.donor_visits_id + AND experimental_data.donor_id = donor_visits.donor_id +WHERE donors.gender IS NOT NULL + AND donor_visits.vaccine_resp IS NOT NULL + AND donor_visits.vaccine = 4 +ORDER BY donors.study_donor_id DESC +\end{lstlisting} + + \section{Full description of FluPrint clinical studies} +\fptable{studies_table}{.7} +{Reference table of clinical studies} +{Clinical study ID used (but remapped) in the database, age information, +vaccine type information, and assay data types of clinical studies are in the +rest of the columns.} +{tbl:studiesDesc} + + + \section{Remaps used in the \flup} + + \begin{table}[htpb] + \begin{tabular}{lll} + \toprule{} + Vaccine received & Vaccine type ID & Vaccine type name \\ + \midrule{} + FluMist IIV4 0.2 mL intranasal spray & 1 & Flumist \\ + FluMist Intranasal spray & 1 & Flumist \\ + FluMist Intranasal Spray 2009–2010 & 1 & Flumist \\ + FluMist Intranasal Spray & 1 & Flumist \\ + Flumist & 1 & Flumist \\ + Fluzone Intradermal-IIV3 & 2 & Fluzone Intradermal \\ + Fluzone Intradermal & 2 & Fluzone Intradermal \\ + GSK Fluarix IIV3 single-dose syringe & 3 & Fluarix \\ + Fluzone 0.5 mL IIV4 SD syringe & 4 & Fluzone \\ + Fluzone 0.25 mL IIV4 SD syringe & 5 & Paediatric Fluzone \\ + Fluzone IIV3 multi-dose vial & 4 & Fluzone \\ + Fluzone single-dose syringe & 4 & Fluzone \\ + Fluzone multi-dose vial & 4 & Fluzone \\ + Fluzone single-dose syringe 2009–2010 & 4 & Fluzone \\ + Fluzone high-dose syringe & 6 & High Dose Fluzone \\ + Fluzone 0.5 mL single-dose syringe & 4 & Fluzone \\ + Fluzone 0.25 mL single-dose syringe & 5 & Paediatric Fluzone \\ + Fluzone IIV3 High-Dose SDS & 6 & High Dose Fluzone \\ + Fluzone IIV4 single-dose syringe & 4 & Fluzone \\ + Fluzone High-Dose syringe & 6 & High Dose Fluzone \\ + \bottomrule{} + \end{tabular} + \caption{Remaps of vaccine type relevant to to the clinical studies + reference table \autoref{tbl:studiesDesc}, and the section on the donor + visits table.}\label{tbl:remapVaccine} + \end{table} + + \begin{table}[htpb] + \begin{tabular}{ll} + \toprule{} + Original & Remapped \\ + \midrule{} + No& 0 \\ + Yes& 1 \\ + IIV injection/im& 2 \\ + Doesn’t know/doesn’t remember/na/does not remember& 3 \\ + LAIV4 intranasal/laiv\_std\_intranasal/laiv\_std\_ intranasal/nasal/intranasal& 4 \\ + \bottomrule{} + \end{tabular} + \caption{caption}\label{tbl:remapHistory} + \end{table} + + \begin{table}[htpb] + \begin{tabular}{ll} + \toprule{} + Original & Remapped \\ + \midrule{} + CMV EBV & 1 \\ + Other immunoassay & 2 \\ + Human Luminex 62–63 plex & 3 \\ + CyTOF phenotyping & 4 \\ + HAI & 5 \\ + Human Luminex 51 plex & 6 \\ + Phospho-flow cytokine stim (PBMC) & 7 \\ + pCyTOF (whole blood) pheno & 9 \\ + pCyTOF (whole blood) phospho & 10 \\ + CBCD & 11 \\ + Human MSD 4 plex & 12 \\ + Lyoplate 1 & 13 \\ + Human MSD 9 plex & 14 \\ + Human Luminex 50 plex & 15 \\ + Other Luminex & 16 \\ + \bottomrule{} + \end{tabular} + \caption{caption}\label{tbl:remapAssays} + \end{table} +\end{appendices} + +\end{document} |