text:
\documentclass{article}

\usepackage{arxiv}

\usepackage[utf8]{inputenc} % allow utf-8 input
\usepackage[T1]{fontenc}    % use 8-bit T1 fonts
\usepackage{lmodern}        % https://github.com/rstudio/rticles/issues/343
\usepackage{hyperref}       % hyperlinks
\usepackage{url}            % simple URL typesetting
\usepackage{booktabs}       % professional-quality tables
\usepackage{amsfonts}       % blackboard math symbols
\usepackage{nicefrac}       % compact symbols for 1/2, etc.
\usepackage{microtype}      % microtypography
\usepackage{graphicx}

\title{Replication of SARS-CoV-2 mutation analysis suggests differences
in per-protein mutation characteristics}

\author{
    William Svoboda
    \thanks{website: \url{https://will.cx}, email:
\href{mailto:wcsvoboda@gmail.com}{\nolinkurl{wcsvoboda@gmail.com}}}
   \\
    Computer Science \\
    Princeton University \\
  Princeton, NJ \\
  \texttt{\href{mailto:wsvoboda@princeton.edu}{\nolinkurl{wsvoboda@princeton.edu}}} \\
   \And
    Brendan McManamon
   \\
    Electrical and Computer Engineering \\
    Princeton University \\
  Princeton, NJ \\
  \texttt{\href{mailto:bm18@princeton.edu}{\nolinkurl{bm18@princeton.edu}}} \\
   \And
    Sara Schwartz
   \\
    Computer Science \\
    Princeton University \\
  Princeton, NJ \\
  \texttt{\href{mailto:sarats@Princeton.edu}{\nolinkurl{sarats@Princeton.edu}}} \\
  }

\usepackage{color}
\usepackage{fancyvrb}
\newcommand{\VerbBar}{|}
\newcommand{\VERB}{\Verb[commandchars=\\\{\}]}
\DefineVerbatimEnvironment{Highlighting}{Verbatim}{commandchars=\\\{\}}
% Add ',fontsize=\small' for more characters per line
\usepackage{framed}
\definecolor{shadecolor}{RGB}{248,248,248}
\newenvironment{Shaded}{\begin{snugshade}}{\end{snugshade}}
\newcommand{\AlertTok}[1]{\textcolor[rgb]{0.94,0.16,0.16}{#1}}
\newcommand{\AnnotationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\AttributeTok}[1]{\textcolor[rgb]{0.77,0.63,0.00}{#1}}
\newcommand{\BaseNTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\BuiltInTok}[1]{#1}
\newcommand{\CharTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\CommentTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textit{#1}}}
\newcommand{\CommentVarTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\ConstantTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\ControlFlowTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{\textbf{#1}}}
\newcommand{\DataTypeTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{#1}}
\newcommand{\DecValTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\DocumentationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\ErrorTok}[1]{\textcolor[rgb]{0.64,0.00,0.00}{\textbf{#1}}}
\newcommand{\ExtensionTok}[1]{#1}
\newcommand{\FloatTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\FunctionTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\ImportTok}[1]{#1}
\newcommand{\InformationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\KeywordTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{\textbf{#1}}}
\newcommand{\NormalTok}[1]{#1}
\newcommand{\OperatorTok}[1]{\textcolor[rgb]{0.81,0.36,0.00}{\textbf{#1}}}
\newcommand{\OtherTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{#1}}
\newcommand{\PreprocessorTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textit{#1}}}
\newcommand{\RegionMarkerTok}[1]{#1}
\newcommand{\SpecialCharTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\SpecialStringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\StringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\VariableTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\VerbatimStringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\WarningTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}

% Pandoc citation processing
\newlength{\csllabelwidth}
\setlength{\csllabelwidth}{3em}
\newlength{\cslhangindent}
\setlength{\cslhangindent}{1.5em}
% for Pandoc 2.8 to 2.10.1
\newenvironment{cslreferences}%
  {}%
  {\par}
% For Pandoc 2.11+
\newenvironment{CSLReferences}[2] % #1 hanging-ident, #2 entry spacing
 {% don't indent paragraphs
  \setlength{\parindent}{0pt}
  % turn on hanging indent if param 1 is 1
  \ifodd #1 \everypar{\setlength{\hangindent}{\cslhangindent}}\ignorespaces\fi
  % set entry spacing
  \ifnum #2 > 0
  \setlength{\parskip}{#2\baselineskip}
  \fi
 }%
 {}
\usepackage{calc} % for calculating minipage widths
\newcommand{\CSLBlock}[1]{#1\hfill\break}
\newcommand{\CSLLeftMargin}[1]{\parbox[t]{\csllabelwidth}{#1}}
\newcommand{\CSLRightInline}[1]{\parbox[t]{\linewidth - \csllabelwidth}{#1}\break}
\newcommand{\CSLIndent}[1]{\hspace{\cslhangindent}#1}

\usepackage{fvextra}
\DefineVerbatimEnvironment{Highlighting}{Verbatim}{breaklines,commandchars=\\\{\}}
\hypersetup{colorlinks=true,linkcolor=blue,citecolor=black,urlcolor=black}
\usepackage{booktabs}

\begin{document}
\maketitle

\def\tightlist{}

\begin{abstract}
The increasing spread of COVID-19, caused by the virus SARS-CoV-2,
raises concerns about the extent to which mutations have occurred across
the viral genome. We present a partial replication of an earlier 2021
study by Wang, R. \emph{et al.} that determined the presence of four
substrains and eleven top mutations in the United States. We analyze a
portion of the authors' data set in order to recreate Figure S1 from the
paper, recapitulating the same features observed in the original figure.
We further generate a summary of mutation characteristics for each of
the 26 named proteins and confirm the significance of the spike protein
at roughly 24\% of all recorded mutations. Our analysis suggests that
additional factors may affect per-protein mutation rate besides protein
length.
\end{abstract}

\keywords{
    COVID-19
   \and
    SARS-CoV-2
   \and
    coronavirus
   \and
    mutation
  }

\hypersetup{urlcolor=blue}

\hypertarget{background}{%
\section{Background}\label{background}}

The spread of coronavirus disease 2019 (COVID-19) continues to be of
scientific and public health interest in the context of the ongoing
global pandemic. Of particular concern is whether the underlying virus,
severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has become
more infectious as a result of new mutations.

A 2021 paper by Wang, R. \emph{et
al.}\textsuperscript{\protect\hyperlink{ref-Wang2021}{1}} analyzed
45,494 complete SARS-CoV-2 genome sequences in order to analyze their
mutations. Of the 12,754 sequences from the United States, the authors
determined the presence of four substrains and eleven top mutations.
Furthermore, the study concluded that two of these substrains had become
more infectious given the mutations found on the spike protein.

We present a partial replication of the original study in order to
validate its findings. A subset of the data used in the paper was
processed in order to recreate selected figures (\autoref{fig:fig1}).
Further analysis was also performed in order to summarize mutation
characteristics for each protein across the entire SARS-CoV-2 genome.
Finally, we discuss the implications of these findings on our current
understanding of the SARS-CoV-2 viral makeup.

\hypertarget{results}{%
\section{Results}\label{results}}

\hypertarget{mutation-tracker}{%
\subsection{Mutation Tracker}\label{mutation-tracker}}

\begin{figure}

{\centering \includegraphics[width=\textwidth]{./figures/original_distribution}

}

\caption{Screenshot of the mutation tracker website showing the distribution of mutations across the entire SARS-CoV-2 genome.}\label{fig:fig1}
\end{figure}

\begin{figure}

{\centering \includegraphics[width=\textwidth]{./figures/updated_distribution}

}

\caption{Genome-wide SARS-CoV-2 mutation distribution as of 31 October 2021. As in the original distribution figure, the natural log frequency of mutations on specific nucleotide positions is plotted on the $y$-axis. The positions of certain landmark proteins are also annotated on the $x$-axis. (author W.S.)}\label{fig:fig2}
\end{figure}

As part of the original
study\textsuperscript{\protect\hyperlink{ref-Wang2021}{1}}, a website
called
\href{https://users.math.msu.edu/users/weig/SARS-CoV-2_Mutation_Tracker.html}{Mutation
Tracker} was created to report the distribution of unique single
mutations across the entire SARS-CoV-2 genome. A screen shot of this
tracker is included in the
paper's\textsuperscript{\protect\hyperlink{ref-Wang2021}{1}}
supplementary information (Figure S1, included as \autoref{fig:fig1} in
this report). We downloaded the summary data provided from the website,
itself generated by Wang, R. \emph{et al.} using sequence data from
GISAID (\url{https://www.gisaid.org/}). The summary data is updated
automatically by the mutation website, with the totals from 31 October
2021 being used in this paper. Split among 26 different files for each
named protein, the summary data was reprocessed into a single data frame
representing the whole viral genome.

As \autoref{fig:fig2} shows, we replicated the original distribution
graph using this updated data. It is evident that NSP3 and the spike
protein are both among the largest proteins and have experienced a high
frequency of unique mutations across their domains. We also observe that
ORF7a, ORF7b, and ORF8 all have relatively high mutation frequencies, in
spite of their short length.

\hypertarget{mutation-summary-statistics}{%
\subsection{Mutation Summary
Statistics}\label{mutation-summary-statistics}}

\begin{table}

\caption{\label{tab:tab1}Summary statistics for each SARS-CoV-2 protein (author W.S.)}
\centering
\begin{tabular}[t]{lrrrr}
\toprule
protein & total mutations & \% total mutations & protein length & mutations/nucleotide\\
\midrule
Spike & 17,134,885 & 23.785 & 3,821 & 4,484.398\\
NSP3 & 12,191,479 & 16.923 & 5,834 & 2,089.729\\
Nucleocapsid & 11,760,145 & 16.324 & 1,259 & 9,340.862\\
RdRp & 7,193,224 & 9.985 & 2,794 & 2,574.525\\
NSP4 & 3,612,766 & 5.015 & 1,499 & 2,410.117\\
\addlinespace
NSP6 & 2,735,245 & 3.797 & 869 & 3,147.578\\
ORF3a & 2,490,822 & 3.457 & 827 & 3,011.877\\
NSP2 & 2,392,844 & 3.321 & 1,913 & 1,250.833\\
Helicase & 2,319,519 & 3.220 & 1,802 & 1,287.191\\
ORF8 & 2,072,786 & 2.877 & 365 & 5,678.866\\
\addlinespace
Exonuclease & 1,851,994 & 2.571 & 1,580 & 1,172.148\\
Membrane & 1,606,101 & 2.229 & 668 & 2,404.343\\
ORF7a & 1,462,425 & 2.030 & 365 & 4,006.644\\
2’-O-ribose MTases & 564,099 & 0.783 & 893 & 631.690\\
endoRNAse & 544,308 & 0.756 & 1,037 & 524.887\\
\addlinespace
3CL & 484,768 & 0.673 & 917 & 528.646\\
NSP1 & 433,783 & 0.602 & 539 & 804.792\\
NSP9 & 302,133 & 0.419 & 338 & 893.885\\
NSP8 & 254,185 & 0.353 & 593 & 428.642\\
ORF10 & 163,272 & 0.227 & 116 & 1,407.517\\
\addlinespace
NSP7 & 125,480 & 0.174 & 248 & 505.968\\
NSP10 & 122,023 & 0.169 & 416 & 293.325\\
ORF6 & 109,947 & 0.153 & 185 & 594.308\\
ORF7b & 56,865 & 0.079 & 131 & 434.084\\
Envelope & 53,802 & 0.075 & 227 & 237.013\\
\addlinespace
NSP11 & 2,578 & 0.004 & 38 & 67.842\\
\bottomrule
\end{tabular}
\end{table}

To quantify our interpretations, we computed further statistics to
summarize the mutation distribution. \autoref{tab:tab1} lists this data
for each protein ordered by total mutation frequency. The spike protein,
NSP3, and the nucleocapsid are the top three proteins, with their
respective frequencies (17,134,885, 12,191,479, and 11,760,145
mutations) being a full order of magnitude greater than the rest.
Interestingly, we note that ORF8, ORF7a, and NSP6 have an unusually high
number of mutations with respect to their length. This corroborates our
previous observations of the mutation distribution graph provided by
\autoref{fig:fig2}.

\hypertarget{mutation-percentage-comparison}{%
\subsection{Mutation Percentage
Comparison}\label{mutation-percentage-comparison}}

\begin{figure}

{\centering \includegraphics[width=\textwidth]{./figures/stacked_bar}

}

\caption{Comparison of total mutation percentages for two collected dates (listed in \texttt{MMDDYYYY} format). For each bar, only the 7 proteins with the largest percentages are named. (author W.S.)}\label{fig:fig3}
\end{figure}

From our summary statistics, we generated a graphical breakdown of the
total mutation percentage. \autoref{fig:fig3} displays this data for the
top 7 proteins between 20 October 2021 and 31 October 2021. Following
from our analysis of the summary statistics, the spike protein made up
roughly 24\% of the total recorded mutation frequency. The graph changes
little between the 10 days comprising each time point. Of note is that a
full 13 proteins make up less than 1\% each of the total.

\hypertarget{discussion}{%
\section{Discussion}\label{discussion}}

Our analysis successfully replicated the original plot of the mutation
distribution represented by \autoref{fig:fig1}. Our version,
\autoref{fig:fig2}, makes use of the updated data provided from the
Mutation Tracker website and recapitulates the features observed there.
Along with the summary statistics from \autoref{tab:tab1}, we
corroborated the significance of the spike protein in the current
pattern of SARS-CoV-2 mutation.

Our approach is limited by the scope of the data used in the study. The
Mutation Tracker website only allows downloading summary data from the
most recent time point. While we had preserved an earlier archive from
20 October 2021, a more thorough comparison than that seen in
\autoref{fig:fig3} would use a greater range of dates. Further study
into the change in mutation over time, either by sampling the tracker or
by directly processing the source data from GISAID, could provide more
insight.

The summary statistics from \autoref{tab:tab1} suggest differences in
mutation characteristics between the 26 studied proteins. Given our
observation of relatively more mutations in shorter proteins, a
conformational analysis could further assess the influence of a
protein's size and 3D structure on the number of mutations it
experiences.

There is also interest in determining the extent, if any, of
evolutionary pressure on mutation rate. Some of the proteins comprising
the lowest percentages of the total mutation frequency, such as NSP11
and the envelope protein, are known to have important roles in the viral
life
cycle\textsuperscript{\protect\hyperlink{ref-Zhang2016}{2},\protect\hyperlink{ref-Schoeman2019}{3}}.
Future work could determine if proteins with functions in viral
replication versus internal regulation experience more mutations.

\hypertarget{conclusion}{%
\section{Conclusion}\label{conclusion}}

In this paper, we presented a partial replication of an earlier study by
Wang, R. \emph{et
al.}\textsuperscript{\protect\hyperlink{ref-Wang2021}{1}}. We
recapitulated their visualization of the genome-wide SARS-CoV-2 mutation
distribution. We further analyzed the underlying data in order to
summarize the influence of each protein on the overall mutation
frequencies. The current global pandemic remains one of the most
significant challenges to public health in recent memory. In light of
emerging variants of SARS-CoV-2, it is even more important to understand
how and where the virus mutates.

\hypertarget{methods}{%
\section{Methods}\label{methods}}

\hypertarget{data-collection-and-preprocessing}{%
\subsection{Data Collection and
Preprocessing}\label{data-collection-and-preprocessing}}

The original paper\textsuperscript{\protect\hyperlink{ref-Wang2021}{1}}
used 45,494 complete genome sequences of SARS-CoV-2 strains from
infected invdividuals. These sequences were sourced from the GISAID
(\url{https://www.gisaid.org/}) database up to 11 September 2020. For
this analysis, we used the 31 October 2021 summary data from the
\href{https://users.math.msu.edu/users/weig/SARS-CoV-2_Mutation_Tracker.html}{Mutation
Tracker} website. The tracker is continually updated with the latest
data from GISAID.

Data preprocessing was accomplished using
R\textsuperscript{\protect\hyperlink{ref-R-base}{4}} and the
dplyr\textsuperscript{\protect\hyperlink{ref-R-dplyr}{5}} package. The
26 \texttt{.csv} files for each named protein were parsed and combined
into a single data frame in order of position on the SARS-CoV-2 genome.
Mutation frequencies were mapped to the corresponding protein name using
nucleotide positions sourced from the NIH
(\url{https://www.ncbi.nlm.nih.gov/gene/}).

\hypertarget{plots-and-figures}{%
\subsection{Plots and Figures}\label{plots-and-figures}}

The plot of the genome-wide SARS-CoV-2 mutation distribution
(\autoref{fig:fig2}) was generated using the Plotly package for
R\textsuperscript{\protect\hyperlink{ref-R-plotly}{6}}. The stacked bar
chart comparing per-protein mutation percentages (\autoref{fig:fig3})
was generated using the ggplot2
package\textsuperscript{\protect\hyperlink{ref-R-ggplot2}{7}} using the
same data as the mutation distribution plot (summarized for each
protein).

\hypertarget{declarations}{%
\section*{Declarations}\label{declarations}}
\addcontentsline{toc}{section}{Declarations}

\hypertarget{data-and-code-availability}{%
\subsection*{Data and Code
Availability}\label{data-and-code-availability}}
\addcontentsline{toc}{subsection}{Data and Code Availability}

This paper uses the summary data from the
\href{https://users.math.msu.edu/users/weig/SARS-CoV-2_Mutation_Tracker.html}{Mutation
Tracker} website up to 31 October 2021. This data was originally sourced
from GISAID (\url{https://www.gisaid.org/}).

All code for this paper is attached in the appendix and is also
available on the author's (W.S.) GitHub at
\url{https://github.com/thisstillwill/QCB455-Fall-2021}.

\hypertarget{acknowledgements}{%
\subsection*{Acknowledgements}\label{acknowledgements}}
\addcontentsline{toc}{subsection}{Acknowledgements}

The authors are grateful to Rui Wang and Guo-Wei Wei for their
correspondance providing the color values used in the Mutation Tracker
and describing the approach used to clean the summary data.

\hypertarget{author-contributions}{%
\subsection*{Author Contributions}\label{author-contributions}}
\addcontentsline{toc}{subsection}{Author Contributions}

W.S. downloaded and preprocessed the summary data from the Mutation
Tracker website. W.S. plotted the mutation frequencies from the summary
data to replicate Figure S1 from the original paper. W.S. conducted
further analysis of the dataset to summarize the mutation data for each
protein.

B.M. downloaded and translated reference genomes for related
coronaviruses into individual proteins. B.M. conducted multiple sequence
alignment for each protein and calculated percent similarity in order to
replicate Figures S5, S7, S8, S10, and S13. Additionally, B.M. added a
table summarizing substitution mutation frequencies.

S.S. conducted multiple sequence alignment and single-nucleotide
polymorphism (SNP) calling on the complete genome sequences to create
SNP profiles for each genome, and used a Jaccard distance-based
representation of the SNP variants as input features into a K-means
clustering algorithm. S.S. plotted the distributions of the resulting
clusters on a U.S. map to replicate Figure 1 from the original paper.

All authors have read and agreed to the published version of the
manuscript.

\hypertarget{competing-interests}{%
\subsection*{Competing Interests}\label{competing-interests}}
\addcontentsline{toc}{subsection}{Competing Interests}

The authors declare no competing interests.

\hypertarget{references}{%
\section*{References}\label{references}}
\addcontentsline{toc}{section}{References}

\hypertarget{refs}{}
\begin{CSLReferences}{0}{0}
\leavevmode\vadjust pre{\hypertarget{ref-Wang2021}{}}%
\CSLLeftMargin{1. }
\CSLRightInline{Wang, R. \emph{et al.}
\href{https://doi.org/10.1038/s42003-021-01754-6}{Analysis of SARS-CoV-2
mutations in the united states suggests presence of four substrains and
novel variants}. \emph{Communications Biology} \textbf{4}, 228 (2021).}

\leavevmode\vadjust pre{\hypertarget{ref-Zhang2016}{}}%
\CSLLeftMargin{2. }
\CSLRightInline{Zhang, M. \emph{et al.}
\href{https://doi.org/10.1128/JVI.01309-16}{Structural biology of the
arterivirus nsp11 endoribonucleases}. \emph{Journal of virology}
\textbf{91}, e01309--16 (2016).}

\leavevmode\vadjust pre{\hypertarget{ref-Schoeman2019}{}}%
\CSLLeftMargin{3. }
\CSLRightInline{Schoeman, D. \& Fielding, B. C.
\href{https://doi.org/10.1186/s12985-019-1182-0}{Coronavirus envelope
protein: Current knowledge}. \emph{Virology Journal} \textbf{16}, 69
(2019).}

\leavevmode\vadjust pre{\hypertarget{ref-R-base}{}}%
\CSLLeftMargin{4. }
\CSLRightInline{R Core Team. \emph{\href{https://www.R-project.org/}{R:
A language and environment for statistical computing}}. (R Foundation
for Statistical Computing, 2021).}

\leavevmode\vadjust pre{\hypertarget{ref-R-dplyr}{}}%
\CSLLeftMargin{5. }
\CSLRightInline{Wickham, H., François, R., Henry, L. \& Müller, K.
\emph{\href{https://CRAN.R-project.org/package=dplyr}{Dplyr: A grammar
of data manipulation}}. (2021).}

\leavevmode\vadjust pre{\hypertarget{ref-R-plotly}{}}%
\CSLLeftMargin{6. }
\CSLRightInline{Sievert, C. \emph{et al.}
\emph{\href{https://CRAN.R-project.org/package=plotly}{Plotly: Create
interactive web graphics via plotly.js}}. (2021).}

\leavevmode\vadjust pre{\hypertarget{ref-R-ggplot2}{}}%
\CSLLeftMargin{7. }
\CSLRightInline{Wickham, H. \emph{et al.}
\emph{\href{https://CRAN.R-project.org/package=ggplot2}{ggplot2: Create
elegant data visualisations using the grammar of graphics}}. (2021).}

\leavevmode\vadjust pre{\hypertarget{ref-R-bookdown}{}}%
\CSLLeftMargin{8. }
\CSLRightInline{Xie, Y.
\emph{\href{https://CRAN.R-project.org/package=bookdown}{Bookdown:
Authoring books and technical documents with r markdown}}. (2021).}

\leavevmode\vadjust pre{\hypertarget{ref-R-formatR}{}}%
\CSLLeftMargin{9. }
\CSLRightInline{Xie, Y.
\emph{\href{https://github.com/yihui/formatR}{formatR: Format r code
automatically}}. (2021).}

\leavevmode\vadjust pre{\hypertarget{ref-R-knitr}{}}%
\CSLLeftMargin{10. }
\CSLRightInline{Xie, Y. \emph{\href{https://yihui.org/knitr/}{Knitr: A
general-purpose package for dynamic report generation in r}}. (2021).}

\leavevmode\vadjust pre{\hypertarget{ref-R-reticulate}{}}%
\CSLLeftMargin{11. }
\CSLRightInline{Ushey, K., Allaire, J. \& Tang, Y.
\emph{\href{https://github.com/rstudio/reticulate}{Reticulate: Interface
to python}}. (2021).}

\leavevmode\vadjust pre{\hypertarget{ref-R-rticles}{}}%
\CSLLeftMargin{12. }
\CSLRightInline{Allaire, J. \emph{et al.}
\emph{\href{https://github.com/rstudio/rticles}{Rticles: Article formats
for r markdown}}. (2021).}

\leavevmode\vadjust pre{\hypertarget{ref-bookdown2016}{}}%
\CSLLeftMargin{13. }
\CSLRightInline{Xie, Y.
\emph{\href{https://bookdown.org/yihui/bookdown}{Bookdown: Authoring
books and technical documents with {R} markdown}}. (Chapman; Hall/CRC,
2016).}

\leavevmode\vadjust pre{\hypertarget{ref-ggplot22016}{}}%
\CSLLeftMargin{14. }
\CSLRightInline{Wickham, H.
\emph{\href{https://ggplot2.tidyverse.org}{ggplot2: Elegant graphics for
data analysis}}. (Springer-Verlag New York, 2016).}

\leavevmode\vadjust pre{\hypertarget{ref-knitr2015}{}}%
\CSLLeftMargin{15. }
\CSLRightInline{Xie, Y. \emph{\href{https://yihui.org/knitr/}{Dynamic
documents with {R} and knitr}}. (Chapman; Hall/CRC, 2015).}

\leavevmode\vadjust pre{\hypertarget{ref-knitr2014}{}}%
\CSLLeftMargin{16. }
\CSLRightInline{Xie, Y.
\href{http://www.crcpress.com/product/isbn/9781466561595}{Knitr: A
comprehensive tool for reproducible research in {R}}. in
\emph{Implementing reproducible computational research} (eds. Stodden,
V., Leisch, F. \& Peng, R. D.) (Chapman; Hall/CRC, 2014).}

\leavevmode\vadjust pre{\hypertarget{ref-plotly2020}{}}%
\CSLLeftMargin{17. }
\CSLRightInline{Sievert, C.
\emph{\href{https://plotly-r.com}{Interactive web-based data
visualization with r, plotly, and shiny}}. (Chapman; Hall/CRC, 2020).}

\end{CSLReferences}

\appendix

\hypertarget{all-code}{%
\section{All Code}\label{all-code}}

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{\# Required packages}
\FunctionTok{library}\NormalTok{(knitr, }\AttributeTok{quietly =} \ConstantTok{TRUE}\NormalTok{)}
\FunctionTok{library}\NormalTok{(formatR, }\AttributeTok{quietly =} \ConstantTok{TRUE}\NormalTok{)}
\FunctionTok{library}\NormalTok{(rticles, }\AttributeTok{quietly =} \ConstantTok{TRUE}\NormalTok{)}
\FunctionTok{library}\NormalTok{(dplyr, }\AttributeTok{quietly =} \ConstantTok{TRUE}\NormalTok{)}
\FunctionTok{library}\NormalTok{(reticulate, }\AttributeTok{quietly =} \ConstantTok{TRUE}\NormalTok{)}
\FunctionTok{library}\NormalTok{(plotly, }\AttributeTok{quietly =} \ConstantTok{TRUE}\NormalTok{)}
\FunctionTok{library}\NormalTok{(bookdown, }\AttributeTok{quietly =} \ConstantTok{TRUE}\NormalTok{)}

\CommentTok{\# Wrap R code in final document}
\NormalTok{knitr}\SpecialCharTok{::}\NormalTok{opts\_chunk}\SpecialCharTok{$}\FunctionTok{set}\NormalTok{(}\AttributeTok{cache =} \ConstantTok{TRUE}\NormalTok{, }\AttributeTok{echo =} \ConstantTok{FALSE}\NormalTok{, }\AttributeTok{message =} \ConstantTok{FALSE}\NormalTok{, }\AttributeTok{error =} \ConstantTok{FALSE}\NormalTok{,}
    \AttributeTok{fig.align =} \StringTok{"center"}\NormalTok{, }\AttributeTok{out.width =} \StringTok{"}\SpecialCharTok{\textbackslash{}\textbackslash{}}\StringTok{textwidth"}\NormalTok{, }\AttributeTok{tidy.opts =} \FunctionTok{list}\NormalTok{(}\AttributeTok{width.cutoff =} \DecValTok{80}\NormalTok{),}
    \AttributeTok{tidy =} \ConstantTok{TRUE}\NormalTok{)}
\NormalTok{knitr}\SpecialCharTok{::}\FunctionTok{include\_graphics}\NormalTok{(}\StringTok{"./figures/original\_distribution.png"}\NormalTok{)}
\CommentTok{\# Read latest mutation summary data into a combined data frame Data retrieved}
\CommentTok{\# from https://users.math.msu.edu/users/weig/SARS{-}CoV{-}2\_Mutation\_Tracker.html}
\NormalTok{path }\OtherTok{\textless{}{-}} \StringTok{"./data/MutationSummary\_10312021/"}
\NormalTok{files }\OtherTok{\textless{}{-}} \FunctionTok{dir}\NormalTok{(path)}
\NormalTok{mutation\_summaries }\OtherTok{\textless{}{-}} \FunctionTok{do.call}\NormalTok{(rbind, }\FunctionTok{lapply}\NormalTok{(}\FunctionTok{paste0}\NormalTok{(path, files), read.csv))}

\CommentTok{\# Add column for raw nucleotide position}
\NormalTok{mutation\_summaries }\OtherTok{\textless{}{-}}\NormalTok{ mutation\_summaries }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{mutate}\NormalTok{(}\AttributeTok{position =} \FunctionTok{as.integer}\NormalTok{(}\FunctionTok{gsub}\NormalTok{(}\StringTok{"[\^{}0{-}9]"}\NormalTok{, }\StringTok{""}\NormalTok{, mutation\_site))) }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{arrange}\NormalTok{(position)}

\CommentTok{\# Create new data frame to track the combined frequencies of each unique}
\CommentTok{\# mutation}
\NormalTok{mutation\_summaries\_combined }\OtherTok{\textless{}{-}}\NormalTok{ mutation\_summaries }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{group\_by}\NormalTok{(position) }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{summarise}\NormalTok{(}\AttributeTok{frequency =} \FunctionTok{sum}\NormalTok{(Total\_frequency))}

\CommentTok{\# Use nucleotide locations from https://www.ncbi.nlm.nih.gov/gene/ to clean}
\CommentTok{\# data Protein name mapping adapted from original paper (Python code): def}
\CommentTok{\# mapDict(self): keys = [\textquotesingle{}[266:805]\textquotesingle{}, \textquotesingle{}[806:2719]\textquotesingle{}, \textquotesingle{}[2720:8554]\textquotesingle{},}
\CommentTok{\# \textquotesingle{}[8555:10054]\textquotesingle{}, \textquotesingle{}[10055:10972]\textquotesingle{}, \textquotesingle{}[10973:11842]\textquotesingle{}, \textquotesingle{}[11843:12091]\textquotesingle{},}
\CommentTok{\# \textquotesingle{}[12092:12685]\textquotesingle{}, \textquotesingle{}[12686:13024]\textquotesingle{}, \textquotesingle{}[13025:13441]\textquotesingle{}, \textquotesingle{}[13442:16236]\textquotesingle{},}
\CommentTok{\# \textquotesingle{}[16237:18039]\textquotesingle{}, \textquotesingle{}[18040:19620]\textquotesingle{}, \textquotesingle{}[19621:20658]\textquotesingle{}, \textquotesingle{}[20659:21552]\textquotesingle{},}
\CommentTok{\# \textquotesingle{}[13442:13480]\textquotesingle{}, \textquotesingle{}[21563:25384]\textquotesingle{}, \textquotesingle{}[25393:26220]\textquotesingle{}, \textquotesingle{}[26245:26472]\textquotesingle{},}
\CommentTok{\# \textquotesingle{}[26523:27191]\textquotesingle{}, \textquotesingle{}[27202:27387]\textquotesingle{}, \textquotesingle{}[27394:27759]\textquotesingle{}, \textquotesingle{}[27756:27887]\textquotesingle{},}
\CommentTok{\# \textquotesingle{}[27894:28259]\textquotesingle{}, \textquotesingle{}[28274:29533]\textquotesingle{}, \textquotesingle{}29558:29674\textquotesingle{}] values = [\textquotesingle{}NSP1\textquotesingle{}, \textquotesingle{}NSP2\textquotesingle{},}
\CommentTok{\# \textquotesingle{}NSP3\textquotesingle{}, \textquotesingle{}NSP4\textquotesingle{}, \textquotesingle{}3CL\textquotesingle{}, \textquotesingle{}NSP6\textquotesingle{}, \textquotesingle{}NSP7\textquotesingle{}, \textquotesingle{}NSP8\textquotesingle{}, \textquotesingle{}NSP9\textquotesingle{}, \textquotesingle{}NSP10\textquotesingle{}, \textquotesingle{}RdRp\textquotesingle{},}
\CommentTok{\# \textquotesingle{}Helicase\textquotesingle{}, \textquotesingle{}Exonuclease\textquotesingle{}, \textquotesingle{}endoRNAse\textquotesingle{}, \textquotesingle{}2’{-}O{-}ribose MTases\textquotesingle{},}
\CommentTok{\# \textquotesingle{}NSP11\textquotesingle{},\textquotesingle{}Spike\textquotesingle{}, \textquotesingle{}ORF3a\textquotesingle{}, \textquotesingle{}Envelope\textquotesingle{}, \textquotesingle{}Membrane\textquotesingle{}, \textquotesingle{}ORF6\textquotesingle{}, \textquotesingle{}ORF7a\textquotesingle{}, \textquotesingle{}ORF7b\textquotesingle{},}
\CommentTok{\# \textquotesingle{}ORF8\textquotesingle{}, \textquotesingle{}Nucleocapsid\textquotesingle{}, \textquotesingle{}ORF10\textquotesingle{}] nameMap = dict(zip(keys,values)) return}
\CommentTok{\# nameMap}
\NormalTok{name\_map }\OtherTok{\textless{}{-}} \FunctionTok{rep}\NormalTok{(}\StringTok{"Others"}\NormalTok{, }\DecValTok{30000}\NormalTok{)}
\NormalTok{name\_map[}\DecValTok{266}\SpecialCharTok{:}\DecValTok{805}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"NSP1"}
\NormalTok{name\_map[}\DecValTok{806}\SpecialCharTok{:}\DecValTok{2719}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"NSP2"}
\NormalTok{name\_map[}\DecValTok{2720}\SpecialCharTok{:}\DecValTok{8554}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"NSP3"}
\NormalTok{name\_map[}\DecValTok{8555}\SpecialCharTok{:}\DecValTok{10054}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"NSP4"}
\NormalTok{name\_map[}\DecValTok{10055}\SpecialCharTok{:}\DecValTok{10972}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"3CL"}
\NormalTok{name\_map[}\DecValTok{10973}\SpecialCharTok{:}\DecValTok{11842}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"NSP6"}
\NormalTok{name\_map[}\DecValTok{11843}\SpecialCharTok{:}\DecValTok{12091}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"NSP7"}
\NormalTok{name\_map[}\DecValTok{12092}\SpecialCharTok{:}\DecValTok{12685}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"NSP8"}
\NormalTok{name\_map[}\DecValTok{12686}\SpecialCharTok{:}\DecValTok{13024}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"NSP9"}
\NormalTok{name\_map[}\DecValTok{13025}\SpecialCharTok{:}\DecValTok{13441}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"NSP10"}
\NormalTok{name\_map[}\DecValTok{13442}\SpecialCharTok{:}\DecValTok{16236}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"RdRp"}
\NormalTok{name\_map[}\DecValTok{16237}\SpecialCharTok{:}\DecValTok{18039}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"Helicase"}
\NormalTok{name\_map[}\DecValTok{18040}\SpecialCharTok{:}\DecValTok{19620}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"Exonuclease"}
\NormalTok{name\_map[}\DecValTok{19621}\SpecialCharTok{:}\DecValTok{20658}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"endoRNAse"}
\NormalTok{name\_map[}\DecValTok{20659}\SpecialCharTok{:}\DecValTok{21552}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"2’{-}O{-}ribose MTases"}
\NormalTok{name\_map[}\DecValTok{13442}\SpecialCharTok{:}\DecValTok{13480}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"NSP11"}
\NormalTok{name\_map[}\DecValTok{21563}\SpecialCharTok{:}\DecValTok{25384}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"Spike"}
\NormalTok{name\_map[}\DecValTok{25393}\SpecialCharTok{:}\DecValTok{26220}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"ORF3a"}
\NormalTok{name\_map[}\DecValTok{26245}\SpecialCharTok{:}\DecValTok{26472}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"Envelope"}
\NormalTok{name\_map[}\DecValTok{26523}\SpecialCharTok{:}\DecValTok{27191}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"Membrane"}
\NormalTok{name\_map[}\DecValTok{27202}\SpecialCharTok{:}\DecValTok{27387}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"ORF6"}
\NormalTok{name\_map[}\DecValTok{27394}\SpecialCharTok{:}\DecValTok{27759}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"ORF7a"}
\NormalTok{name\_map[}\DecValTok{27756}\SpecialCharTok{:}\DecValTok{27887}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"ORF7b"}
\NormalTok{name\_map[}\DecValTok{27894}\SpecialCharTok{:}\DecValTok{28259}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"ORF8"}
\NormalTok{name\_map[}\DecValTok{28274}\SpecialCharTok{:}\DecValTok{29533}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"Nucleocapsid"}
\NormalTok{name\_map[}\DecValTok{29558}\SpecialCharTok{:}\DecValTok{29674}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"ORF10"}
\NormalTok{mutation\_summaries\_combined}\SpecialCharTok{$}\NormalTok{protein }\OtherTok{\textless{}{-}}\NormalTok{ name\_map[mutation\_summaries\_combined}\SpecialCharTok{$}\NormalTok{position]}

\CommentTok{\# Remove rows with default protein name}
\NormalTok{mutation\_summaries\_combined }\OtherTok{\textless{}{-}} \FunctionTok{subset}\NormalTok{(mutation\_summaries\_combined, protein }\SpecialCharTok{!=} \StringTok{"Others"}\NormalTok{)}
\CommentTok{\# Color values adapted from settings used in original paper (Python code):}
\CommentTok{\# color\_continuous\_scale=[ \textquotesingle{}rgb(160,200,255)\textquotesingle{},\textquotesingle{}yellow\textquotesingle{},\textquotesingle{}rgba(13,143,129,0.7)\textquotesingle{},}
\CommentTok{\# \textquotesingle{}rgba(119,74,1750,0.8)\textquotesingle{}, \textquotesingle{}rgba(235,74,64,0.8)\textquotesingle{},\textquotesingle{}rgb(62,109,178)\textquotesingle{},}
\CommentTok{\# \textquotesingle{}rgb(160,200,255)\textquotesingle{}, \textquotesingle{}rgb(253,174,107)\textquotesingle{}, \textquotesingle{}rgba(119,74,1750,0.8)\textquotesingle{},}
\CommentTok{\# \textquotesingle{}rgb(160,200,255)\textquotesingle{},\textquotesingle{}black\textquotesingle{},\textquotesingle{}rgba(13,143,129,0.7)\textquotesingle{}, \textquotesingle{}rgba(119,74,1750,0.8)\textquotesingle{},}
\CommentTok{\# \textquotesingle{}rgba(235,74,64,0.8)\textquotesingle{},\textquotesingle{}rgb(62,109,178)\textquotesingle{}, \textquotesingle{}rgb(160,200,255)\textquotesingle{}, \textquotesingle{}black\textquotesingle{},}
\CommentTok{\# \textquotesingle{}rgb(260,200,255)\textquotesingle{}, \textquotesingle{}rgb(160,200,255)\textquotesingle{},\textquotesingle{}yellow\textquotesingle{},\textquotesingle{}rgba(13,143,129,0.7)\textquotesingle{},}
\CommentTok{\# \textquotesingle{}rgba(119,74,1750,0.8)\textquotesingle{}, \textquotesingle{}rgba(235,74,64,0.8)\textquotesingle{},\textquotesingle{}rgb(62,109,178)\textquotesingle{},}
\CommentTok{\# \textquotesingle{}rgb(160,200,255)\textquotesingle{}, \textquotesingle{}rgb(253,174,107)\textquotesingle{}])}
\NormalTok{scatter\_palette }\OtherTok{\textless{}{-}} \FunctionTok{c}\NormalTok{(}\StringTok{"\#a0c8ff"}\NormalTok{, }\StringTok{"yellow"}\NormalTok{, }\StringTok{"\#0d8f81"}\NormalTok{, }\StringTok{"\#774ad6"}\NormalTok{, }\StringTok{"\#eb4a40"}\NormalTok{, }\StringTok{"\#3e6db2"}\NormalTok{,}
    \StringTok{"\#a0c8ff"}\NormalTok{, }\StringTok{"\#fdae6b"}\NormalTok{, }\StringTok{"\#774ad6"}\NormalTok{, }\StringTok{"\#a0c8ff"}\NormalTok{, }\StringTok{"black"}\NormalTok{, }\StringTok{"\#0d8f81"}\NormalTok{, }\StringTok{"\#774ad6"}\NormalTok{, }\StringTok{"\#eb4a40"}\NormalTok{,}
    \StringTok{"\#3e6db2"}\NormalTok{, }\StringTok{"\#a0c8ff"}\NormalTok{, }\StringTok{"black"}\NormalTok{, }\StringTok{"\#04c8ff"}\NormalTok{, }\StringTok{"\#a0c8ff"}\NormalTok{, }\StringTok{"yellow"}\NormalTok{, }\StringTok{"\#0d8f81"}\NormalTok{, }\StringTok{"\#774ad6"}\NormalTok{,}
    \StringTok{"\#eb4a40"}\NormalTok{, }\StringTok{"\#3e6db2"}\NormalTok{, }\StringTok{"\#a0c8ff"}\NormalTok{, }\StringTok{"\#fdae6b"}\NormalTok{)}

\CommentTok{\# Create scatter plot of mutations}
\NormalTok{scatter\_fig }\OtherTok{\textless{}{-}} \FunctionTok{plot\_ly}\NormalTok{(}\AttributeTok{data =}\NormalTok{ mutation\_summaries\_combined, }\AttributeTok{x =} \SpecialCharTok{\textasciitilde{}}\NormalTok{position, }\AttributeTok{y =} \SpecialCharTok{\textasciitilde{}}\FunctionTok{log}\NormalTok{(frequency),}
    \AttributeTok{type =} \StringTok{"scatter"}\NormalTok{, }\AttributeTok{mode =} \StringTok{"markers"}\NormalTok{, }\AttributeTok{color =} \SpecialCharTok{\textasciitilde{}}\NormalTok{protein, }\AttributeTok{colors =}\NormalTok{ scatter\_palette)}
\NormalTok{scatter\_fig }\OtherTok{\textless{}{-}}\NormalTok{ scatter\_fig }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{layout}\NormalTok{(}\AttributeTok{title =} \StringTok{"Mutation distribution across the SARS{-}CoV{-}2 genome"}\NormalTok{, }\AttributeTok{xaxis =} \FunctionTok{list}\NormalTok{(}\AttributeTok{title =} \StringTok{""}\NormalTok{,}
        \AttributeTok{showgrid =} \ConstantTok{FALSE}\NormalTok{, }\AttributeTok{tickangle =} \DecValTok{60}\NormalTok{, }\AttributeTok{tickmode =} \StringTok{"array"}\NormalTok{, }\AttributeTok{nticks =} \DecValTok{11}\NormalTok{, }\AttributeTok{tickvals =} \FunctionTok{c}\NormalTok{(}\DecValTok{536}\NormalTok{,}
            \DecValTok{5638}\NormalTok{, }\DecValTok{9305}\NormalTok{, }\DecValTok{10514}\NormalTok{, }\DecValTok{11408}\NormalTok{, }\DecValTok{14839}\NormalTok{, }\DecValTok{17138}\NormalTok{, }\DecValTok{23474}\NormalTok{, }\DecValTok{25807}\NormalTok{, }\DecValTok{26359}\NormalTok{, }\DecValTok{28904}\NormalTok{),}
        \AttributeTok{ticktext =} \FunctionTok{c}\NormalTok{(}\StringTok{"NSP1"}\NormalTok{, }\StringTok{"NSP3"}\NormalTok{, }\StringTok{"NSP4"}\NormalTok{, }\StringTok{"3CL"}\NormalTok{, }\StringTok{"NSP6"}\NormalTok{, }\StringTok{"RdRp"}\NormalTok{, }\StringTok{"Helicase"}\NormalTok{, }\StringTok{"S"}\NormalTok{,}
            \StringTok{"ORF3a"}\NormalTok{, }\StringTok{"E"}\NormalTok{, }\StringTok{"N"}\NormalTok{)), }\AttributeTok{yaxis =} \FunctionTok{list}\NormalTok{(}\AttributeTok{title =} \StringTok{"ln(Frequency)"}\NormalTok{, }\AttributeTok{showgrid =} \ConstantTok{FALSE}\NormalTok{),}
        \AttributeTok{font =} \FunctionTok{list}\NormalTok{(}\AttributeTok{family =} \StringTok{"Serif"}\NormalTok{, }\AttributeTok{size =} \DecValTok{32}\NormalTok{, }\AttributeTok{color =} \StringTok{"black"}\NormalTok{), }\AttributeTok{showlegend =} \ConstantTok{FALSE}\NormalTok{,}
        \AttributeTok{margin =} \FunctionTok{list}\NormalTok{(}\AttributeTok{t =} \DecValTok{200}\NormalTok{))}
\FunctionTok{save\_image}\NormalTok{(scatter\_fig, }\StringTok{"./figures/updated\_distribution.pdf"}\NormalTok{, }\AttributeTok{width =} \DecValTok{2304}\NormalTok{, }\AttributeTok{height =} \DecValTok{1152}\NormalTok{)}
\NormalTok{knitr}\SpecialCharTok{::}\FunctionTok{include\_graphics}\NormalTok{(}\StringTok{"./figures/updated\_distribution.pdf"}\NormalTok{)}
\CommentTok{\# Generate mutation statistics by protein}
\NormalTok{aggregate\_summaries }\OtherTok{\textless{}{-}}\NormalTok{ mutation\_summaries\_combined }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{group\_by}\NormalTok{(protein) }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{summarise}\NormalTok{(}\AttributeTok{total\_mutations =} \FunctionTok{sum}\NormalTok{(frequency)) }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{mutate}\NormalTok{(}\AttributeTok{raw\_percent =}\NormalTok{ (total\_mutations}\SpecialCharTok{/}\FunctionTok{sum}\NormalTok{(total\_mutations)) }\SpecialCharTok{*} \DecValTok{100}\NormalTok{, }\AttributeTok{nucleotide\_length =} \FunctionTok{sapply}\NormalTok{(protein,}
        \ControlFlowTok{switch}\NormalTok{, }\AttributeTok{NSP1 =}\NormalTok{ \{}
            \DecValTok{805} \SpecialCharTok{{-}} \DecValTok{266}
\NormalTok{        \}, }\AttributeTok{NSP2 =}\NormalTok{ \{}
            \DecValTok{2719} \SpecialCharTok{{-}} \DecValTok{806}
\NormalTok{        \}, }\AttributeTok{NSP3 =}\NormalTok{ \{}
            \DecValTok{8554} \SpecialCharTok{{-}} \DecValTok{2720}
\NormalTok{        \}, }\AttributeTok{NSP4 =}\NormalTok{ \{}
            \DecValTok{10054} \SpecialCharTok{{-}} \DecValTok{8555}
\NormalTok{        \}, }\StringTok{\textasciigrave{}}\AttributeTok{3CL}\StringTok{\textasciigrave{}} \OtherTok{=}\NormalTok{ \{}
            \DecValTok{10972} \SpecialCharTok{{-}} \DecValTok{10055}
\NormalTok{        \}, }\AttributeTok{NSP6 =}\NormalTok{ \{}
            \DecValTok{11842} \SpecialCharTok{{-}} \DecValTok{10973}
\NormalTok{        \}, }\AttributeTok{NSP7 =}\NormalTok{ \{}
            \DecValTok{12091} \SpecialCharTok{{-}} \DecValTok{11843}
\NormalTok{        \}, }\AttributeTok{NSP8 =}\NormalTok{ \{}
            \DecValTok{12685} \SpecialCharTok{{-}} \DecValTok{12092}
\NormalTok{        \}, }\AttributeTok{NSP9 =}\NormalTok{ \{}
            \DecValTok{13024} \SpecialCharTok{{-}} \DecValTok{12686}
\NormalTok{        \}, }\AttributeTok{NSP10 =}\NormalTok{ \{}
            \DecValTok{13441} \SpecialCharTok{{-}} \DecValTok{13025}
\NormalTok{        \}, }\AttributeTok{RdRp =}\NormalTok{ \{}
            \DecValTok{16236} \SpecialCharTok{{-}} \DecValTok{13442}
\NormalTok{        \}, }\AttributeTok{Helicase =}\NormalTok{ \{}
            \DecValTok{18039} \SpecialCharTok{{-}} \DecValTok{16237}
\NormalTok{        \}, }\AttributeTok{Exonuclease =}\NormalTok{ \{}
            \DecValTok{19620} \SpecialCharTok{{-}} \DecValTok{18040}
\NormalTok{        \}, }\AttributeTok{endoRNAse =}\NormalTok{ \{}
            \DecValTok{20658} \SpecialCharTok{{-}} \DecValTok{19621}
\NormalTok{        \}, }\StringTok{\textasciigrave{}}\AttributeTok{2’{-}O{-}ribose MTases}\StringTok{\textasciigrave{}} \OtherTok{=}\NormalTok{ \{}
            \DecValTok{21552} \SpecialCharTok{{-}} \DecValTok{20659}
\NormalTok{        \}, }\AttributeTok{NSP11 =}\NormalTok{ \{}
            \DecValTok{13480} \SpecialCharTok{{-}} \DecValTok{13442}
\NormalTok{        \}, }\AttributeTok{Spike =}\NormalTok{ \{}
            \DecValTok{25384} \SpecialCharTok{{-}} \DecValTok{21563}
\NormalTok{        \}, }\AttributeTok{ORF3a =}\NormalTok{ \{}
            \DecValTok{26220} \SpecialCharTok{{-}} \DecValTok{25393}
\NormalTok{        \}, }\AttributeTok{Envelope =}\NormalTok{ \{}
            \DecValTok{26472} \SpecialCharTok{{-}} \DecValTok{26245}
\NormalTok{        \}, }\AttributeTok{Membrane =}\NormalTok{ \{}
            \DecValTok{27191} \SpecialCharTok{{-}} \DecValTok{26523}
\NormalTok{        \}, }\AttributeTok{ORF6 =}\NormalTok{ \{}
            \DecValTok{27387} \SpecialCharTok{{-}} \DecValTok{27202}
\NormalTok{        \}, }\AttributeTok{ORF7a =}\NormalTok{ \{}
            \DecValTok{27759} \SpecialCharTok{{-}} \DecValTok{27394}
\NormalTok{        \}, }\AttributeTok{ORF7b =}\NormalTok{ \{}
            \DecValTok{27887} \SpecialCharTok{{-}} \DecValTok{27756}
\NormalTok{        \}, }\AttributeTok{ORF8 =}\NormalTok{ \{}
            \DecValTok{28259} \SpecialCharTok{{-}} \DecValTok{27894}
\NormalTok{        \}, }\AttributeTok{Nucleocapsid =}\NormalTok{ \{}
            \DecValTok{29533} \SpecialCharTok{{-}} \DecValTok{28274}
\NormalTok{        \}, }\AttributeTok{ORF10 =}\NormalTok{ \{}
            \DecValTok{29674} \SpecialCharTok{{-}} \DecValTok{29558}
\NormalTok{        \}), }\AttributeTok{mutations\_per\_nucleotide =}\NormalTok{ total\_mutations}\SpecialCharTok{/}\NormalTok{nucleotide\_length) }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{arrange}\NormalTok{(}\FunctionTok{desc}\NormalTok{(total\_mutations))}
\NormalTok{knitr}\SpecialCharTok{::}\FunctionTok{kable}\NormalTok{(aggregate\_summaries, }\AttributeTok{format =} \StringTok{"latex"}\NormalTok{, }\AttributeTok{booktabs =} \ConstantTok{TRUE}\NormalTok{, }\AttributeTok{digits =} \DecValTok{3}\NormalTok{,}
    \AttributeTok{col.names =} \FunctionTok{c}\NormalTok{(}\StringTok{"protein"}\NormalTok{, }\StringTok{"total mutations"}\NormalTok{, }\StringTok{"\% total mutations"}\NormalTok{, }\StringTok{"protein length"}\NormalTok{,}
        \StringTok{"mutations/nucleotide"}\NormalTok{), }\AttributeTok{caption =} \StringTok{"Summary statistics for each SARS{-}CoV{-}2 protein (author W.S.)"}\NormalTok{,}
    \AttributeTok{format.args =} \FunctionTok{list}\NormalTok{(}\AttributeTok{big.mark =} \StringTok{","}\NormalTok{))}
\CommentTok{\# Read mutation summary data for 10202021}
\NormalTok{path }\OtherTok{\textless{}{-}} \StringTok{"./data/MutationSummary\_10202021/"}
\NormalTok{files }\OtherTok{\textless{}{-}} \FunctionTok{dir}\NormalTok{(path)}
\NormalTok{mutation\_summaries\_10202021 }\OtherTok{\textless{}{-}} \FunctionTok{do.call}\NormalTok{(rbind, }\FunctionTok{lapply}\NormalTok{(}\FunctionTok{paste0}\NormalTok{(path, files), read.csv))}

\CommentTok{\# Find percent mutations by protein for 10202021}
\NormalTok{bar\_data\_10202021 }\OtherTok{\textless{}{-}}\NormalTok{ mutation\_summaries\_10202021 }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{mutate}\NormalTok{(}\AttributeTok{position =} \FunctionTok{as.integer}\NormalTok{(}\FunctionTok{gsub}\NormalTok{(}\StringTok{"[\^{}0{-}9]"}\NormalTok{, }\StringTok{""}\NormalTok{, mutation\_site))) }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{group\_by}\NormalTok{(position) }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{summarise}\NormalTok{(}\AttributeTok{frequency =} \FunctionTok{sum}\NormalTok{(Total\_frequency))}
\NormalTok{bar\_data\_10202021}\SpecialCharTok{$}\NormalTok{protein }\OtherTok{\textless{}{-}}\NormalTok{ name\_map[bar\_data\_10202021}\SpecialCharTok{$}\NormalTok{position]}
\NormalTok{bar\_data\_10202021 }\OtherTok{\textless{}{-}} \FunctionTok{subset}\NormalTok{(bar\_data\_10202021, protein }\SpecialCharTok{!=} \StringTok{"Others"}\NormalTok{)}
\NormalTok{bar\_data\_10202021 }\OtherTok{\textless{}{-}}\NormalTok{ bar\_data\_10202021 }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{group\_by}\NormalTok{(protein) }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{summarise}\NormalTok{(}\AttributeTok{total\_mutations =} \FunctionTok{sum}\NormalTok{(frequency)) }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{mutate}\NormalTok{(}\AttributeTok{raw\_percent =}\NormalTok{ (total\_mutations}\SpecialCharTok{/}\FunctionTok{sum}\NormalTok{(total\_mutations)) }\SpecialCharTok{*} \DecValTok{100}\NormalTok{) }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{arrange}\NormalTok{(}\FunctionTok{desc}\NormalTok{(total\_mutations)) }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{select}\NormalTok{(protein, raw\_percent) }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{head}\NormalTok{(}\DecValTok{7}\NormalTok{)}
\NormalTok{others }\OtherTok{\textless{}{-}} \FunctionTok{list}\NormalTok{(}\StringTok{"Others"}\NormalTok{, }\DecValTok{100} \SpecialCharTok{{-}} \FunctionTok{sum}\NormalTok{(bar\_data\_10202021}\SpecialCharTok{$}\NormalTok{raw\_percent))}
\NormalTok{bar\_data\_10202021 }\OtherTok{\textless{}{-}}\NormalTok{ bar\_data\_10202021 }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{rbind}\NormalTok{(others) }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{arrange}\NormalTok{(raw\_percent) }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{mutate}\NormalTok{(}\AttributeTok{date =} \StringTok{"10202021"}\NormalTok{)}

\CommentTok{\# Find percent mutations by protein for latest data (10312021)}
\NormalTok{bar\_data\_10312021 }\OtherTok{\textless{}{-}}\NormalTok{ aggregate\_summaries }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{select}\NormalTok{(protein, raw\_percent) }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{head}\NormalTok{(}\DecValTok{7}\NormalTok{)}
\NormalTok{others }\OtherTok{\textless{}{-}} \FunctionTok{list}\NormalTok{(}\StringTok{"Others"}\NormalTok{, }\DecValTok{100} \SpecialCharTok{{-}} \FunctionTok{sum}\NormalTok{(bar\_data\_10312021}\SpecialCharTok{$}\NormalTok{raw\_percent))}
\NormalTok{bar\_data\_10312021 }\OtherTok{\textless{}{-}}\NormalTok{ bar\_data\_10312021 }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{rbind}\NormalTok{(others) }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{arrange}\NormalTok{(raw\_percent) }\SpecialCharTok{\%\textgreater{}\%}
    \FunctionTok{mutate}\NormalTok{(}\AttributeTok{date =} \StringTok{"10312021"}\NormalTok{)}

\CommentTok{\# Combine data into single data frame for plotting}
\NormalTok{bar\_data }\OtherTok{\textless{}{-}} \FunctionTok{rbind}\NormalTok{(bar\_data\_10202021, bar\_data\_10312021)}

\CommentTok{\# Created stacked bar chart for each date ordered by percent}
\NormalTok{stacked\_bar }\OtherTok{\textless{}{-}} \FunctionTok{ggplot}\NormalTok{(bar\_data, }\FunctionTok{aes}\NormalTok{(}\AttributeTok{x =}\NormalTok{ date, }\AttributeTok{y =}\NormalTok{ raw\_percent, }\AttributeTok{group =}\NormalTok{ raw\_percent)) }\SpecialCharTok{+}
    \FunctionTok{geom\_col}\NormalTok{(}\FunctionTok{aes}\NormalTok{(}\AttributeTok{fill =}\NormalTok{ protein), }\AttributeTok{width =} \FloatTok{0.9}\NormalTok{) }\SpecialCharTok{+} \FunctionTok{scale\_fill\_brewer}\NormalTok{(}\AttributeTok{type =} \StringTok{"qual"}\NormalTok{,}
    \AttributeTok{palette =} \StringTok{"Paired"}\NormalTok{) }\SpecialCharTok{+} \FunctionTok{geom\_text}\NormalTok{(}\FunctionTok{aes}\NormalTok{(}\AttributeTok{label =} \FunctionTok{paste0}\NormalTok{(protein, }\StringTok{" {-} "}\NormalTok{, }\FunctionTok{round}\NormalTok{(raw\_percent,}
    \AttributeTok{digits =} \DecValTok{1}\NormalTok{), }\StringTok{"\%"}\NormalTok{)), }\AttributeTok{position =} \FunctionTok{position\_stack}\NormalTok{(}\AttributeTok{vjust =} \FloatTok{0.5}\NormalTok{)) }\SpecialCharTok{+} \FunctionTok{ggtitle}\NormalTok{(}\StringTok{"Percentage of total mutations by protein"}\NormalTok{) }\SpecialCharTok{+}
    \FunctionTok{xlab}\NormalTok{(}\StringTok{"Date"}\NormalTok{) }\SpecialCharTok{+} \FunctionTok{ylab}\NormalTok{(}\StringTok{"Percent"}\NormalTok{) }\SpecialCharTok{+} \FunctionTok{labs}\NormalTok{(}\AttributeTok{fill =} \StringTok{"Protein"}\NormalTok{) }\SpecialCharTok{+} \FunctionTok{theme\_classic}\NormalTok{(}\AttributeTok{base\_family =} \StringTok{"serif"}\NormalTok{) }\SpecialCharTok{+}
    \FunctionTok{theme}\NormalTok{(}\AttributeTok{plot.title =} \FunctionTok{element\_text}\NormalTok{(}\AttributeTok{hjust =} \FloatTok{0.5}\NormalTok{))}
\FunctionTok{ggsave}\NormalTok{(}\StringTok{"./figures/stacked\_bar.pdf"}\NormalTok{)}
\NormalTok{knitr}\SpecialCharTok{::}\FunctionTok{include\_graphics}\NormalTok{(}\StringTok{"./figures/stacked\_bar.pdf"}\NormalTok{)}
\CommentTok{\# Cite packages}
\NormalTok{knitr}\SpecialCharTok{::}\FunctionTok{write\_bib}\NormalTok{(}\FunctionTok{c}\NormalTok{(}\FunctionTok{.packages}\NormalTok{(), }\StringTok{"bookdown"}\NormalTok{), }\StringTok{"packages.bib"}\NormalTok{)}
\FunctionTok{sessionInfo}\NormalTok{(}\AttributeTok{package =} \ConstantTok{NULL}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\hypertarget{session-info}{%
\section{Session Info}\label{session-info}}

\begin{verbatim}
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base
##
## other attached packages:
## [1] bookdown_0.24   plotly_4.10.0   ggplot2_3.3.5   reticulate_1.22
## [5] dplyr_1.0.7     rticles_0.21    formatR_1.11    knitr_1.36
##
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.1.1   xfun_0.27          purrr_0.3.4        lattice_0.20-45
##  [5] colorspace_2.0-2   vctrs_0.3.8        generics_0.1.0     htmltools_0.5.2
##  [9] viridisLite_0.4.0  yaml_2.2.1         utf8_1.2.2         rlang_0.4.12
## [13] pillar_1.6.4       glue_1.4.2         withr_2.4.2        DBI_1.1.1
## [17] RColorBrewer_1.1-2 lifecycle_1.0.1    stringr_1.4.0      munsell_0.5.0
## [21] gtable_0.3.0       htmlwidgets_1.5.4  codetools_0.2-18   evaluate_0.14
## [25] labeling_0.4.2     fastmap_1.1.0      crosstalk_1.2.0    fansi_0.5.0
## [29] Rcpp_1.0.7         scales_1.1.1       jsonlite_1.7.2     farver_2.1.0
## [33] png_0.1-7          digest_0.6.28      stringi_1.7.5      grid_4.1.1
## [37] rprojroot_2.0.2    here_1.0.1         tools_4.1.1        magrittr_2.0.1
## [41] lazyeval_0.2.2     tibble_3.1.5       crayon_1.4.1       tidyr_1.1.4
## [45] pkgconfig_2.0.3    ellipsis_0.3.2     Matrix_1.3-4       data.table_1.14.2
## [49] assertthat_0.2.1   rmarkdown_2.11     httr_1.4.2         R6_2.5.1
## [53] compiler_4.1.1
\end{verbatim}

\bibliographystyle{unsrt}
\bibliography{packages.bibreferences.bib}

\end{document}


text:
\title{TEMPLATE\\Preliminary draft}
\titlerunning{TEMPLATE}
\toctitle{TEMPLATE}
\author{J. Author\thanks{Supported by the NSF}}
\authorrunning{J. Author}
\institute{Some University, Department of Mathematics}
\maketitle

\begin{abstract}
Here we would write the abstract.
\end{abstract}

\section*{Introduction}
Here we would write the introduction.  See the next chapter for more detailed instructions.

\section{Some interesting math}
Here we could write some interesting math.  We could even index things like
$H^2(X,\mathbf O_X)$\index{$H^i$} with math.

\section{Some interesting computations}
Here we could write some interesting computations.

\subsection{A subsection title}  This is a subsection.

\subsubsection{A subsubsection title}  This is a subsubsection.

\section{Yet more stuff}
For some reason {\tt bibtex} gets upset if we don't have at least one
citation in the paper \cite{MR47:3318}.

% Local Variables:
% mode: latex
% mode: reftex
% tex-main-file: "chapter-wrapper.tex"
% reftex-keep-temporary-buffers: t
% reftex-use-external-file-finders: t
% reftex-external-file-finders: (("tex" . "make FILE=%f find-tex") ("bib" . "make FILE=%f find-bib"))
% End:

ext: tex


text:
\documentclass[10pt,legalpaper]{book}

%Packages Used
\usepackage[utf8]{inputenc}
\usepackage{amsthm}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}

%Meta-Data
\author{Carlos Baez}
\title{Notes on Optimization and Mathematical Programming}

%Formatting
\theoremstyle{definition}
\newtheorem{defn}{Definition}[chapter]
\theoremstyle{plain}
\newtheorem{prop}{Property}
\theoremstyle{remark}
\newtheorem{proprem}{Remark}[prop]

%Document
\makeindex
\begin{document}
\maketitle
\tableofcontents
\part{Mathematical Background}

\chapter{Mathematical Preliminaries}
\section{Set Theory and Logic}
\section{Norms}
\section{Analysis}

\chapter{Linear Algebra}
\section{Fields and Number Systems}
\section{Matrices}
\begin{defn}[Matrix]
A \textbf{matrix} is defined as a rectangular array of numbers written as:
\begin{center}
\begin{equation}
\textbf{A} = ||a_{ij}|| =
	\left[
		\begin{matrix}
			a_{11} & a_{12} & ... & a_{1n}  \\
			a_{21} & a_{22} & ... & a_{2n}  \\
			:      &        &     & :       \\
			a_{m1} & a_{m2} & ... & a_{mn}  \\
		\end{matrix}
	\right]
\end{equation}
\end{center}
where the matrix is of size \textit{m}$\times$\textit{n}.
\end{defn}
\subsection{Properties of Matrices}
Given matrices \textbf{A}, \textbf{B}, \textbf{C}:
\begin{prop}[Equal Matrices]
\textbf{A} = \textbf{B} iff $a_{ij} = b_{ij} \forall$ i,j.
\end{prop}

\begin{prop}[Adding Matrices]
If \textbf{C}=\textbf{A}+\textbf{B}, the elements of
\textbf{C} are given by $c_{ij}=a_{ij}+b_{ij}.$
\begin{proprem}
\textbf{A} and \textbf{B} can only be added if they both have the same number of rows and columns.
\end{proprem}
\end{prop}

\begin{prop}[Commutative Law]
\begin{equation}
\textbf{A}+\textbf{B}=\textbf{B}+\textbf{A}
\end{equation}
\end{prop}

\begin{prop}[Associative Law]
\begin{equation}
\textbf{A}+(\textbf{B}+\textbf{C})=(\textbf{A}+\textbf{B})+\textbf{C}=\textbf{A}+\textbf{B}+\textbf{C}
\end{equation}
\end{prop}

\begin{prop}[Multiplying Matrices]
Given $\lambda \in \mathbb{R}, \lambda$\textbf{A} = $||\lambda a_{ij}||.$ and $\lambda$\textbf{A} =\textbf{A}$\lambda.$
\end{prop}

\begin{prop}[Product of Matrices]
\textbf{A}\textbf{B} is the product of \textbf{A} and \textbf{B}. If \textbf{A} is an m$\times$n matrix and \textbf{B} is an u$\times$v matrix, then \textbf{C}=\textbf{A}\textbf{B} is defined and \textbf{C} is an m$\times$v matrix.
\begin{proprem}
The product of \textbf{A} and \textbf{B} is defined iff the number of columns of \textbf{A} is equal to the number of rows in \textbf{B}.
\end{proprem}
\begin{proprem}
The elements of \textbf{C} can be computed by
\begin{equation}
c_{ij}=\sum\limits_{k=1}^r a_{ik}b_{kj},\;\;\;i=1,...,m; j=1,...,n.
\end{equation}
\end{proprem}
\end{prop}

Hello1
\subsection{Special Types of Matrices}
Why is this
\section{Vector Spaces}
\section{Determinants}
\section{Linear Transformations}
\section{Similarity}
\section{Polynomials and Polynomial Matrices}
\section{Similarity 2}
\section{Matrix Analysis}

%Mathematical Optimization
\part{Mathematical Optimization}
\chapter{Linear Programming}
\chapter{Convex Optimization}
\section{Convex Programming}
\subsection{Semidefinite Programming}
\subsection{Linear-Fractional Programming}

\chapter{Combinatorial Optimization}
\section{Graphs}
\section{Network Flow}
\section{Integer Programming}

\chapter{Nonlinear Optimization}
\section{Nonlinear Programming}
\section{Fractional Programming}

\chapter{Multiple Objective Optimization}
\chapter{Stochastic Programming}
\chapter{Goal Programming}
\chapter{Constraint Programming}

%Algorithms
\part{Algorithms and Heuristics}
\chapter{Computational Complexity}
\chapter{Constrained Nonlinear Algorithms}
\chapter{Dynamic Programming}

%Heuristics
\chapter{Heuristics}
\section{Evolutionary Algorithms}
\section{Simulated Annealing}

\end{document}

max_stars_repo_name: cabaez/MathematicalProgNotes

max_issues_repo_name: cabaez/MathematicalProgNotes

max_forks_repo_name: cabaez/MathematicalProgNotes


text:
%focus on sub-goal
\chapter{Design and Implementation}\label{chapter:prototype_implementation}

Implementing an \gls{ea} requires a number of considerations. A robust fitness assignment method needs to be determined. A mechanism for uniform diversity preservation must be conceptualized, and a well-thought-out solution to address constraints in the problem definition must be designed. There are successful \glspl{ea} that can be used for this; these algorithms have been well researched, and their domain applicability is well defined.

\Gls{nsga}-II \parencite{Jain2013AnOptimization} is an \gls{ea} that uses an elitist non-dominated sorting mechanism to approximate the Pareto front. It is an improved version of an initial \gls{nsga} algorithm developed by Srinivas et al. \parencite{Srinivas1994MuiltiobjectiveAlgorithms} in 1996.
\Gls{nsga}-II uses crowding function for diversity preservation and ranking of the Pareto fronts. However, it is better suited for two-objective problems and does not perform well on multi-objective problems. \Gls{nsga}-III \parencite{Mkaouer2015Many-objectiveNSGA-III} was introduced to address the inability of \gls{nsga}-II to solve multi-objective problems (i.e., three or more objectives).

Other successful \glspl{ea} applicable to our problem can be found in practice. In the following sections, we implement the \gls{nsga} algorithm and its heuristic adaptations to suit our problem definition.

\section{Non-Dominated Sorting Genetic Algorithms}
Designing domination-based \gls{emo} algorithms for multi-objective problems requires additional considerations. For instance, an increase in the number of objectives may cause an increase in the non-dominated population derived from the randomly generated population. Since domination-based \gls{emo} prioritizes non-dominated solutions, the solution space becomes rather large. Therefore, the search process becomes slow, and the algorithm becomes inefficient. With increased solution space also arises an increase in the number of neighbors. Evaluation of diversity measures (i.e., crowding distance) becomes computationally expensive. \Gls{nsga}-III was designed to address these problems.

As defined in \ref{sec:problem_definition}, our problem can be mapped to many objective optimization problems. Each objective function is a function of a \gls{poi} and its characteristics. A chosen algorithm needs to efficiently compare more than three objectives, and this informs the decision to implement and adapt the \gls{nsga}-III to our problem domain. Next, we explore the aspects that comprise the \gls{nsga}-III and constitute our algorithm.

\subsection*{The \Gls{nsga}-III Algorithm}
The \gls{nsga}-III algorithm was developed to address the inefficiency of \gls{nsga}-II in handling multi-objective problems. In its most basic form, \gls{nsga}-III can be seen as the \gls{nsga}-II algorithm with minor changes; that is, most algorithmic steps in \gls{nsga}-III are from the \gls{nsga}-II algorithm. Figure \ref{fig:nsgaII} and \ref{fig:nsgaIII} depict an overview of the major algorithmic steps in \gls{nsga}-II and \gls{nsga}-III, respectively. From these figures can be observed the similarity between the two algorithms. However, a major difference between them is in the last step of next-population generation. Unlike \gls{nsga}-II, which sorts individuals by crowding distance to decide what part of a Pareto front to reject or accept, \gls{nsga}-III uses a sorting technique based on a reference point. Sorting by crowding distance is explained in Section \ref{sec:diversitypreservation}. The reference-point-based sorting technique and other intrinsic details of the \gls{nsga}-III are explained below.

\begin{figure}
    \centering
    \includegraphics[width=9cm]{NSGA-II}
    \caption{\gls{nsga}-II overview}
    \label{fig:nsgaII}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=9cm]{NSGA-III}
    \caption{\gls{nsga}-III overview}
    \label{fig:nsgaIII}
\end{figure}

\subsubsection{Initial Population}
Let us consider $P_t$ as a parent population of the $t^{th}$ generation of the algorithm (i.e., iteration $t$). $P_t$ with $N$ chromosomes is either a population taken from previous generation  $t-1$ or an initial population taken from a defined heuristic.\\
The heuristic used for initial population generation at the start of the algorithm is particularly important for our problem. Simply sampling randomly from the problem space would be logically inefficient. The number of objective problems in practice is typically capped at 15–20. \Gls{nsga}-III was evaluated on 20-objective DTLZ1-DTLZ2. Considering that the objective function in our problem domain comprises a single \gls{poi} and its characteristics, an intelligent pre-selection process must be defined to reduce the problem space.

\subsubsection{Constraint Handling}\label{sec:sub_conshandling}
Figure \ref{fig:unconstrainedsearchspace} illustrates a typical solution search space for an unconstrained problem. For unconstrained problems, all members of the population are regarded as feasible; however, some individuals can become infeasible (i.e., members in the gray shaded area; Figure \ref{fig:constrainedserachspace}), some individuals become infeasible (i.e., members in gray shaded area). Additional techniques are required for handling constraints in \gls{nsga}-III. Two papers address problems with \gls{nsga}-III both with \parencite{Deb2013AnConstraints} and without \parencite{Jain2013AnApproach} constraints.

\begin{figure}[h!]
\centering
\subfloat[Unconstrained solution search space]{
\includegraphics[width=.4\textwidth]{UnconstrainedSolutionSpace}
\label{fig:unconstrainedsearchspace}}\qquad
\subfloat[Constrained solution search space]{\includegraphics[width=.4\textwidth]{ConstrainedSearchSpace}
\label{fig:constrainedserachspace}}
\caption{Solution Search Space}
\label{fig:constrainedvsunconstrained}
\end{figure}

Consider an individual $x$ in the solution space with the following constraints:

\begin{align}
    \tag{constraint 1 }g_j(x) \leq b_j \label{eq:4a}\\
    \tag{constraint 2}h_k(x) = 0 \label{eq:4b}
\end{align}

The authors of \parencite{Jain2013AnApproach} suggest calculating a constraint valuation value $CV(x)$ to measure the degree to which $x$ violates the constraints using the following formula:

\begin{align}
    \tag{1}CV(x) = \sum_{j=1}^J \langle g'_j(x)\rangle + \sum_{k=1}^K \lvert h'_k(x)\rvert \label{eq:4c}
\end{align}

where $g'_j(x)$ is the normalized constraint function for Constraint \ref{eq:4a}, such that $g'_j(x) = g_j(x) / b_j  - 1 \leq 0$ and $h'_k(x)$ is the normalized constraint function for Constraint \ref{eq:4b}. The angle bracket operator $\langle \alpha \rangle$ returns the negative of $\alpha$, if $\alpha < 0$ and returns zero, otherwise. Assuming that Constraint \ref{eq:4a} is a greater or equal constraint (i.e., $g_j(x) \geq b_j$), then $\langle \alpha \rangle$ returns the negative of $\alpha$, if $\alpha > 0$ and returns zero, otherwise.

The notion of feasibility together with the constraint violation measure $CV(x)$, are combined to derive a constraint domination principle. Consider two individuals $x_1$ and $x_2$ in the solution space. An individual $x_1$ is said to constraint dominate $x_2$ if any of the following conditions holds:
\begin{enumerate}
    \item $x_1$ is feasible and $x_2$ is infeasible;
    \item Both $x_1$ and $x_2$ are infeasible and $x_1$ has a lower \gls{cv} value; or
    \item Both $x_1$ and $x_2$ are feasible and $x_1$ dominates $x_2$ ($x_1 \succ x_2$).
\end{enumerate}

\subsubsection{Offspring Production}
Three major aspects in generating the child population $Q_t$ are defined as follows:

\begin{enumerate}
    \item \textbf{Tournament Selection}: To select two individuals that will be eligible to produce offspring for the next generation, a tournament selection procedure is carried out. Two individuals are selected at random from $P_t$, and a binary selection of the better solution is applied, as shown in Algorithm  \ref{alg:tournament_selection}. Individuals $p_1$ and$p_2$ are compared to select the one that does not violate the constraints. If both individuals violate the constraints, the constraint violation values for the two individuals are computed. Then, the algorithm selects the individual with the lesser constraint violation value. Two parents for reproduction (crossover and mutation) are selected using this tournament procedure.
    \item \textbf{Cross over}: A crossover operator is used to generate offspring according to a crossover rate $p_c$. The crossover rate represents the proportion of parents on which a crossover operator will act. A uniform crossover ensures that each offspring will inherit some characteristics of both parents. Either the 1-point crossover or the n-point crossover is commonly used for binary representations. When using the 1-point crossover, a crossover site $k$ is randomly selected. is randomly selected. Then, two offspring are created by exchanging the segments of the parents \parencite{Talbi2009Metaheuristics:Implementation}. Consequently, n-point crossover implies $n$ crossover sites.
    %The most commonly used rates are in the interval [0.45, 0.95]. Adaptive techniques for the crossover rate may also be useful
    \item \textbf{Mutation}: Mutation operators act on a single individual. An element of the representation is mutated (changed) according to some probability $p_m$. \parencite{Talbi2009Metaheuristics:Implementation} recommends small values for $p_m$ ($p_m \in [0.001, 0.01]$). A flip operator is commonly used when using binary representations of individuals. %Some strategies initialize the mutation probability to 1/k where k is the number of decision variables, that is, in average only one variable is mutated.

\end{enumerate}
During reproduction, duplicate offspring are checked and eliminated. A child population $Q_t$ of size $N$ is generated by repeating tournament selection and reproduction. $R_t$ is generated by combining $P_t$ and $Q_t$ (i.e., $R_t = P_t \cup Q_t$). \Gls{nsga}-III and \gls{nsga}-II use an elitist preservation technique (III and IV in figure \ref{fig:nsgaIII}) to preserve elite members of the new population.

\begin{algorithm}[h!]
  \caption{Tournament selection procedure in  constrained \gls{nsga}-III}\label{alg:tournament_selection}
  \SetKwInOut{Input}{Input}
  \SetKwInOut{Output}{Output}
  \Input{$p_1$, $p_2$}
  \Output{$p'$}
    \uIf{feasible($p_1$) and not feasible($p_2$)}
    {
        $p' = p_1$
    }
    \ElseIf{not feasible($p_1$) and feasible($p_2$)}
    {
        $p' = p_2$
    }
    \ElseIf{not feasible($p_1$) and not feasible($p_2$)}
    {
        \If{CV($p_1$)  $>$ CV($p_2$)}
        {
            $p' = p_2$
        }
        \ElseIf{CV($p_1$)  $<$ CV($p_2$)}
        {
            $p' = p_1$
        }
        \Else{
            $p' = random(p_1, p_2)$
        }
    }
    \ElseIf{feasible($p_1$) and feasbile($p_2$)}
    {
        $p' = random(p_1, p_2)$
    }

  \end{algorithm}

\subsubsection{Non-Dominated Sorting} \label{sec:sub:nondominatedsorting}
Non-dominated sorting ranks members of the population $R_t$ into different non-dominated levels ($F_1$, $F_2$, and so on). In optimization problems without constraints, the sorting is performed using individuals fitness values, while in optimization problems, individuals are ranked according to their constraint domination count. We introduce a new symbol $\succ_+$ for constraint domination. For example, let $A$, $B$, and $C$ be selected members of the population  such that $A \succ_+ B \succ_+ C$ (i.e., individual $A$ constraint dominates individual $B$ and $C$, individual $B$ constraint dominates $C$). Figure \ref{fig:paretolevels} depicts $A$, $B$, and $C$ on a graph. $A$ has the highest domination count and hence belongs to domination level $F_1$. $B$ and $C$ belong to level $F_2$.

\begin{figure}
\centering
\includegraphics[width=.5\textwidth]{NonDominatedFronts}
\caption{Ranking of Pareto Frontier according to Constrained Domination}\label{fig:paretolevels}
\end{figure}

Thereafter, the number of feasible solutions $N_f$ in $R_t$ are counted. If $N_f \leq N$, all the feasible solutions are added to $S_t$ to start the next generation $P_{t+1}$. The remaining $K = N - \lvert P_{t+1} \rvert$ population slots are filled with top-level infeasible solutions, that is solutions having smaller \gls{cv} values.\\
If $N_f > N$, the algorithm follows the unconstrained selection procedure described in \parencite{Deb2013AnConstraints}. Each level is selected one at a time until a size $N$ is reached or $N$ is exceeded for the first time. Assuming the $l^{th}$ level is the level that causes the size of $S_t$ to exceed $N$, then solutions in level $l$ are only partially accepted into $S_t$.
%niche-preservation operator that computes the crowding distance for every last level member as the summation of objective-wise normalized distance between two neighboring solutions.

\subsubsection{Reference Point-Based Sorting}
\Gls{nsga}-III uses a pre-defined set of reference points to ensure diversity in solutions. The reference points can either be preferential points supplied by the user or pre-defined in a structured manner \parencite{Deb2013AnConstraints}. \\
In the case of structured reference points, \gls{nsga}-III uses a systematic approach that creates points on a normalized hyperplane. Points are placed on the hyperplane with the same orienteering in all axes. The total number of reference points ($H$) is given by the following:
\begin{align*}
    H &= \begin{pmatrix}
    M + p - 1 \\
    p
    \end{pmatrix}
\end{align*}
where $M$ is the number of objective functions, and $p$ is the number of divisions to consider along each objective. For example, using three objectives ($M = 3$) and four divisions ($p = 4$), the reference points are placed on a triangle with its apex at $(1,0,0), (0,1,0),$ and $(0,0,1)$. $H = \big(\begin{smallmatrix}
  6\\
  4
\end{smallmatrix}\big)$, so 15 reference points are created. The obtained reference points are widely distributed on the entire normalized hyperplane, and the solutions are widely distributed on or near the Pareto-optimal front. Population members associated with the created reference points are emphasized. The population size is recommended to be the smallest multiple of four that is greater than $H$. In the case of user-supplied reference points, the user can mark $H$ points on the normalized hyperplane.

After determining the reference points on a hyperplane, the algorithm adaptively normalizes the population members. First, the ideal point\footnote{The algorithm presented in the paper focused on minimization problems, thus the ideal point was differently defined. Since, we are faced with a maximization problem, our ideal point is defined differently.} of the population $S_t$ is determined by identifying the maximum value ($z_i^{max}$) for each objective function $i = 1, 2,...., M$ in $\cup_{\tau=0}^{t} S_\tau$ and by constructing the ideal point $\Tilde{z} = (z_1^{max}, z_2^{max},..., z_M^{max})$. Using $f'_i(x) = z_i^{max} - f_i(x)$, each objective value of $S_t$ is translated. Thereafter, the extreme point ($z_i^{ext}$) in each objective axis is identified. Given  a weight vector $d_i$, that is close to the $i^{th}$ objective axis, an achievement scalarizing function (ASF) is given by the following:
\begin{align*}
    ASF(x, d_i) = max_{k=1}^M \frac{f'_i(x)}{d_i}, \qquad for\ i = 1, 2,...,M.
\end{align*}
$z_i^{ext}$ is identified by finding a solution ($x \in S_t$) that results in the minimum corresponding ASF. These $M$ extreme vectors are then used to represent an $M$-dimensional hyperplane. The intercept $a_i$ of the $i^{th}$ objective axis and the linear axis can then be calculated by $a_i = z_i^{max} - z_i^{ext}$. These maximum and extreme values are updated at every generation using the current population. Putting this all together, the objective function can be normalized using the following:

\begin{align}
    f_i^n(x) &=  \frac{z_i^{max} - f_i(x)}{z_i^{max} - z_i^{ext}}, \qquad for\ i = 1, 2,...,M. \label{eq:41a}
\end{align}

After normalizing the objective functions, each population member is associated with a reference point. A reference line for each reference point on the hyperplane is defined by joining the reference point with the ideal point. Then, a perpendicular distance between each population member and each reference line is calculated. The reference point that has the closest reference line to a population individual in the normalized objective space is considered to be associated with the population individual.

A niche-preserving procedure to generate the remaining $K = N - \lvert P_{t+1} \rvert $ population slots is performed. population slots is performed. A reference point may have one or more population members associated with it, or there might not be any associated reference point. Initially, the population members from $P_{t+1} = S_t \setminus F_l$ that are associated with the $j^{th}$ reference point are counted as $\rho_j$ (niche count). Then, the reference point set $J_{min} = {j:argmin_j \ \rho_j}$ having minimum $\rho_j$ is identified. In case of multiple such reference points, a random $\Tilde{j} \in J_{min}$ is selected from $J_{min}$.
If there is no associated $P_{t+1 }$ member to the reference point $\Tilde{j}$ ($\rho_{\Tilde{j}} = 0$), then either of the following scenario is possible:
\begin{enumerate}
    \item There exists one or members in front $F_l$ that are associated with the reference point $\Tilde{j}$, or
    \item The front $F_l$ does not have any member associated with the reference point $\Tilde{j}$.
\end{enumerate}
In case 1, the first case, the reference point with the shortest perpendicular distance from the reference line is added to $P_{t+1}$; then, the count $\rho_{\Tilde{j}}$ for the reference point $\Tilde{j}$ is incremented by 1. In case 2, the reference point is excluded from further consideration for the current generation.\\
In the event that $\rho_{\Tilde{j}} \geq 1$, a randomly chosen member (if one exists) from the front  $F_l$ associated with the reference point $\Tilde{j}$ is added to $P_{t+1}$. The count $\rho_{\Tilde{j}}$ for the reference point $\Tilde{j}$ is incremented by 1. The procedure is repeated $K$ times until all vacant slots in $\rho_{\Tilde{j}}$ are filled.

\section{Implementation and Algorithmic Approaches}\label{sec:impl}
In this section, details of the implementation and several approaches to solving the many objective \glspl{op} discussed in Section \ref{sec:problem_definition} are presented. Our approach consists of an optimization process that comprises varying initialization steps and constraint-handling techniques. Algorithm \ref{alg:opt_framework} outlines the general structure of the approach, which addresses both initialization and the constraint-handling technique. It is important to note a fundamental difference between the non-dominated sorting explained in Section \ref{sec:sub:nondominatedsorting} and outlined in Algorithm \ref{alg:opt_framework}. Our approach performs a non-dominated sorting of the fitness values of individuals, as it is performed in the case of \gls{nsga}-III without constraints. Additionally, we have chosen a criterion-based fitness assignment for our approach. Therefore, non-dominated sorting is solely used for approximating the Pareto front, in our implementation.

The implementation part of this thesis is realized using \textit{Python}. \textit{Pandas}  \parencite{McKinney2010DataPython} and \textit{NumPy} \parencite{Harris2020ArrayNumPy} are Python packages that provide functionalities for storing and processing the data. The dataset used for this work is stored in a local MySQL database. NumPy and Pandas are used to transform the dataset after querying stored data from the database. During runtime, the transformed dataset is stored as Pandas tables. An existing framework containing \gls{nsga}-III is selected to reduce implementation time. The evolutionary computation framework  \gls{DEAP} \parencite{Fortin2012DEAP:Easy} is chosen because of the flexibility of use it guarantees via explicit algorithms and transparent data structures. It can also be used with the Python multi-processing module in order to bypass the Python GIL and parallelize the optimization process. Due to the high dimensionality of the considered problems, efficiency through parallelization provides a benefit for the calculations in test cases. Lastly, \gls{DEAP} is an open-source framework with a vast number of examples, which helps better understand the framework.

\begin{algorithm}
  \caption{General outline of optimization process}\label{alg:opt_framework}
 \SetKwInOut{Input}{Input}
  \SetKwInOut{Output}{Output}
  \Input{$In$ Initialization, $CHT$ Constraint Handling Technique, $M$ No of Objectives}
  \Output{Individuals in first non-dominated front}
  Calculate the number of reference point (H) to place on the hyper-plan using $p$ and $M$\\
  Generate the initial population via method defined by $In$\\
  Evaluate fitness of individuals depending on $CHT$\\
  Realize the non-dominated sorting of fitness values\\
    \While {\# generations $\leq$ MAXGEN}{
        Select two parents $P1$ and $P2$ using the tornament method\\
        Apply crossover between P1 and P2 with a probability $P_c$\\
        Evaluate fitness of individuals depending on $CHT$\\
        Realize the non-dominated sorting of fitness values\\
        Normalize the individuals\\
        Associate individuals with the reference points\\
        Apply the niche preservation operation\\
        Keep the niche obtained solutions for the next generation

    }
    Evaluate fitness of individuals depending on $CHT$\\
    Realize the non-dominated sorting of fitness values\\
    return feasible individuals belonging to the first non-dominated front
  \end{algorithm}

The initialization and constraint-handling technique of the optimization process is realized in different ways. The following sections describe several variations that are implemented and evaluated in this thesis.

\subsection{Population Encoding}
Individuals of the population are encoded as binary arrays. An example individual of the population is encoded as $[1, 0, 1, 1]$ where each bit represents a region under consideration. Thus, a population is a nested array of individuals. For example, a population of two individuals is encoded as $[[1, 0, 1, 1], [1,1,1, 1]]$. Each individual is also assigned a so-called \textit{strategy}; an array of index positions in the randomly selected regions that helps in correctly mapping the bits to the appropriate region.

\subsection{Start Population}
The first step in each approach is to initialize the first population from which the \gls{ea} begins. The performance of the algorithm can vary depending on the starting population used. A main criterion is good diversity at the beginning of the optimization process. Four initialization techniques are used in this thesis. The first is a naive method that randomly generates individuals for the start population. We adapt the metric used in \parencite{cbrecsys2014} to create a second method that generates individuals having high scores on one or more of the defined constraints. The third method uses only feasible individuals for the starting population. The last method is a hybrid of similarity-based and feasibility-based initialization.

\subsubsection{Random Initialization}\label{sec_randominit}
The random initialization technique is a naive approach used by the baseline \gls{nsga}-III algorithm. To generate an individual, a binary array of all 1s is generated via NumPy. Then, the flip-bit operator provided by \gls{DEAP} is used to randomly flip the bits of the array. Variant 1 in Figure \ref{fig:feasibleinfeasiblecomp} shows the results of the comparison between the number of feasible and infeasible individuals generated for the initial population using this approach. The preliminary test was designed to generate 92 individuals for the starting population, and it was repeated 100 times. On average, only 0.81\% of individuals were feasible.

\subsubsection{Similarity-Based Initialization}\label{sec_sim}
Similarity-based initialization incorporates the similarity metric used in \parencite{cbrecsys2014} and is performed in two steps:

\begin{enumerate}
    \item Rank the regions according to similarity metric, and
    \item Randomly generate individuals from regions with higher rank
\end{enumerate}
In step 1, a similarity value between the desired score for preference $f_q$ and the actual score of preference $f_c$ is computed as follows:
\begin{align*}
sim_{feature} (f_q, f_c) = 1 - \frac{\lvert f_q - f_c \rvert} {max(f_q , f_c)}
\end{align*}
where $f_c$ is the normalized score in a range of $[0,1]$. The desired score for each input preference is defined as 1 ($f_c=1$). A ranking of all regions computed using the weighted sum average of $sim_{feature}$ is then computed using the below metric:

\begin{align*}
similarity (q, c) = \frac{\sum_{i=1}^{n} w_i * sim_i(q_i, c_i)} {\sum_{i=1}^{n} w_i}
\end{align*}

The weight $w_i$ prioritizes specific preferences over the other. For example, activity preferences selected by the user have higher priority than other influential parameters such as safety. In our implementation, we prioritize user activity preferences and weather for the chosen months. Safety from crime has second priority, hence it has less weight. We do not consider parameters such as budget or minimal stay duration for each region for the similarity metric. After obtaining a ranking of similarities between the regions and preferences, regions having ranks below 0.7 are not considered for the next steps. For more information on the similarity metric, refer to \parencite{cbrecsys2014}.

In Step 2, the random initialization described in Section \ref{sec_randominit} is used to generate the starting population. Variant 2 in Figure \ref{fig:feasibleinfeasiblecomp} shows the results of the comparison between the number of feasible and infeasible individuals generated for the initial population using this approach. Test parameters similar to those in Variant 1 are used. Results show that this approach does not improve the proportion of feasible to infeasible individuals in an initial population (0.5\% average count of feasible individuals).

\begin{figure}[t]
    \centering
    \includegraphics[width=.7\textwidth]{FeasibleInfeasible12}
    \caption{Resulting initial population distribution of 92 randomly initialized individuals (100 trials)}
    \label{fig:feasibleinfeasiblecomp}
\end{figure}

\subsubsection{Feasibility-Based Initialization}\label{sec_feas}
A third approach to generating the initial population is applying a feasibility-based initialization. In this method, a randomly generated individual is first tested for its adherence to all defined constraints. An individual is deemed infeasible if it does not meet all constraint criteria. Feasible individuals are added to the starting population, while infeasible individuals are discarded. This process is repeated until a population size $N$ is attained. An issue with this approach is the large amount of time needed to generate population size $N$. Depending on the population size $N$ and the number of objectives used, more computational time is often needed to obtain $N$ feasible individuals. Hence, algorithm variants using this approach are assigned a smaller population size.

\subsubsection{Hybrid Initialization}
The hybrid approach to initial population generation implemented in this thesis comprises feasibility-based and similarity-based initialization. In this approach, the similarity metric described in \ref{sec_sim} is first applied to generate a ranking of the regions. Regions having a rank below 0.7 are removed from further consideration. Next, possible region combinations are randomly generated, and only feasible individuals of size $N$ are allowed into the starting population.

\subsection{Constraint Handling Technique}
Handling of constraints and evaluation of the fitness of individuals in \gls{DEAP} are tightly coupled due to the interfaces provided in the framework. In this section, two approaches implemented for handling constraints are described. In both approaches, the constrained objective function is transformed into an unconstrained objective function by applying penalty functions.

\subsubsection{Distance-Based Penalties}
In this approach, the fitness value of infeasible individuals is replaced with a penalized fitness value $f^\mathrm{penalty}_i(x)$ using the penalty metric shown below:
\begin{align*}
    f^\mathrm{penalty}_i(\mathbf{x}) = \Delta_i - \sum_{j=1}^{N_{con}} d_{ij}(\mathbf{x})
\end{align*}
where $\Delta_i$ is the smallest possible score any individual of the population could obtain. In our case, $\Delta_i$ is set as $[0,0,0]$. $N_{con}$ represents the total constraints defined in our problem formulation (i.e., $N_{con} = 5$). $d_{ij}$ sums the distance of the individual to the feasible region for each constraint $j$. Depending on the constraint, each region that comprises an individual is penalized independently. Specifically, for the budget (\ref{eq:3_1b}) and duration (\ref{eq:3_1c}) constraints, $d_i^k$ for a specific region in individual $i$ is computed as follows:

\begin{align*}
    d_i^k =  \frac{c_k}{CV_i(x) * q_k}
\end{align*}
$CV_i(x)$ is the sum over the normalized constraints and represents the degree of constraint violation of individual $i$ as described in Equation \ref{eq:4c}. $c_k$ is the actual budget or duration needed for the region $k$, while $q_k$ is the user-input budget or duration.

Regions that in a combination violate Constraint \ref{eq:3_1e} are assigned the value $d_i^k = \frac{c_k}{q_k}$, where $c_k$ is the number of user preferences the region does not fulfill, and $q_k$ is the total number of user preferences. For all other constraints (\ref{eq:3_1f} and \ref{eq:3_1f2}), the region combinations are penalized together with a constant value.

\subsubsection{Quadratic Penalty}
This approach transforms the objective function $z_i(x)$ into an unconstrained objective fitness function of the following form:
\begin{align*}
    f_i(x) = z_i(x) - \sum_{j=1}^{N_{con}}P_j(x)
\end{align*}
The penalty function $P_j$ is a quadratic function of the penalty coefficient for a given constraint. In the implementation, the penalty coefficient was chosen as the constraint violation degree $CV_i(x)$ over all constraints. Thus, the fitness value of an individual$i$ is computed as follows:
\begin{align*}
    f_{i}(x) := \begin{cases}
                            z_i(x), &\text{if}\quad CV_i(x) \leq 0 ;\\
                            [CV_i(x)]^2, &\text{otherwise};
                             \end{cases}
\end{align*}
Preliminary tests show better results using a quadratic increase of the violation degree as opposed to a linear increase.

\section{Discussion}
In this chapter, we presented and described in great detail \gls{nsga}-III as the baseline algorithm chosen for implementing the solution to our optimization problem. We described various heuristic solutions implemented in this thesis to modify the baseline algorithm. The objective function of the initial problem formulation was transformed in order to better handle the defined constraints. An implication of this transformation is that feasible members of the population will dominate infeasible members during non-dominated sorting. The tournament selection procedure described in Algorithm \ref{alg:tournament_selection} was relaxed such that individuals having higher fitness values are prioritized. Comparison between individual degrees of constraint violation is thus eliminated. However, the niche preservation operation is still used to select between individuals having the same Pareto front.

max_stars_repo_path: Item Combination For DestiRec/chapters/04_protype_implementation.tex

max_stars_repo_name: idaShobs/destirec-composite

max_issues_repo_name: idaShobs/destirec-composite

max_forks_repo_name: idaShobs/destirec-composite


text: \section{Summary \& Outlook} %(Alle)

max_stars_repo_name: rmetzler/NoSQL-for-Dependable-Systems

max_issues_repo_name: rmetzler/NoSQL-for-Dependable-Systems

max_forks_repo_name: rmetzler/NoSQL-for-Dependable-Systems


text:
\section{Conclusion}%
\label{sec:Conclusion}%
Insert your text here.%
\lipsum[1]%
\linebreak

max_stars_repo_name: geraldinepyh/reportWriterDemo

max_issues_repo_name: geraldinepyh/reportWriterDemo

max_forks_repo_name: geraldinepyh/reportWriterDemo


text:
\documentclass{tufte-handout}
\usepackage{fetamont}
\usepackage[T1]{fontenc}
\usepackage{caption}

\title{Research Summary}

\author[Viachaslau Bernat]{Viachaslau Bernat}

%\date{28 March 2010} % without \date command, current date is supplied

\hypersetup{colorlinks}

%\geometry{showframe} % display margins for debugging page layout

\usepackage{graphicx} % allow embedded images
  \setkeys{Gin}{width=\linewidth,totalheight=\textheight,keepaspectratio}
  \graphicspath{{graphics/}} % set of paths to search for images
\usepackage{amsmath}  % extended mathematics
\usepackage{booktabs} % book-quality tables
\usepackage{units}   % non-stacked fractions and better unit spacing
\usepackage{multicol} % multiple column layout facilities
\usepackage{fancyvrb} % extended verbatim environments
  \fvset{fontsize=\normalsize}% default font size for fancy-verbatim environments

% Standardize command font styles and environments
\newcommand{\doccmd}[1]{\texttt{\textbackslash#1}}% command name -- adds backslash automatically
\newcommand{\docopt}[1]{\ensuremath{\langle}\textrm{\textit{#1}}\ensuremath{\rangle}}% optional command argument
\newcommand{\docarg}[1]{\textrm{\textit{#1}}}% (required) command argument
\newcommand{\docenv}[1]{\textsf{#1}}% environment name
\newcommand{\docpkg}[1]{\texttt{#1}}% package name
\newcommand{\doccls}[1]{\texttt{#1}}% document class name
\newcommand{\docclsopt}[1]{\texttt{#1}}% document class option name
\newenvironment{docspec}{\begin{quote}\noindent}{\end{quote}}% command specification environment

\begin{document}

\maketitle% this prints the handout title, author, and date

\begin{abstract}
\begin{fullwidth}
\noindent
My research interests are centered around application of organic synthesis to designing molecules with useful properties for biomedical applications. This broad definition led me to pursue research in the range of subjects from medicinal chemistry to molecular pharmacology and chemical biology of membrane proteins and non-coding RNA, and development of diagnostic fluorescent sensors.
\end{fullwidth}
\end{abstract}

\section{Biased modulators of GPCRs}\label{sec:bias-gpcr}
%\subsection{Structure-based Design of $\mu$-Opioid Receptor Agonists}\label{sec:mor}
\vspace{-0.5cm}

\newthought{Biased $\mu$-opioid receptor agonists.} As a part of multi-center
collaboration I synthesized newly designed biased agonists of $\mu$-opioid receptor. My main contribution
was synthesis and preparative separation of individual enantiomers of the optimized hits. The stereochemistry
turned out to be the key parameter in SAR that defined signaling bias and affinity. The most
promising molecule -- PZM21 -- did not recruit $\beta$-2-arrestin, while being single-digit nanomolar
activator of G\textsubscript{i} protein signaling. Animal testing demonstrated strong analgesia, lack of side effects, and improved
therapeutic index\cite[-1.5cm]{Manglik2016}.
\begin{marginfigure}[-7.5cm]
	\includegraphics{PZM21.png}
	\caption{Structure of PZM21}
\end{marginfigure}

\begin{marginfigure}\label{fig:cxcrHM}
	\includegraphics{cxcr3hpp_small.png}
	\caption{Homology model of CXCR3 with docked small molecule allosteric modulator.}
\end{marginfigure}
\newthought{Biased negative allosteric modulators of CXCR3.} Using homology modeling
I constructed the structural model of CXCR3 based on X-ray structure of CXCR4 and
docked known allosteric modulators (ChEMBL) into the model. This helped to identify amino acid
residues that regulate small molecule binding and downstream signaling. Attempts to engage
a nucleophilic subpocket via reversible covalent interaction with strategically placed electrophilic boronic acid led
to discovery of the first biased negative allosteric modulator or CXCR3 (\textbf{14}, Fig. \ref{fig:cxcr3}).\cite{Bernat2014}

Optimization of physico-chemical properties (cLogP and PSA) via introduction of flexible side
chain led to identificaiton of the negative allosteric modulator
with the opposite signaling bias (\textbf{1b}, Fig. \ref{fig:cxcr3}).
These biased modulators exhibited opposite probe dependence -- they preferentially
downregulated activation of CXCR3 by CXCL10 or CXCL11 chemokines.\cite{Bernat2015a}
Affinities of these small-molecule allosteric modulators were similar (double-digit nM) as
measured in radioligand RAMX3 displacement assay.\cite{Bernat2012}

Structural model of CXCR3 that I prepared was subsequently used \textit{as is} in virtual screening for identification
of completely novel chemotypes of specific and dual CXCR3 and CXCR4 allosteric modulators.\cite{Schmidt2015}

\newpage
\begin{fullwidth}
\phantomsection
\centering
\begin{figure*}
	\includegraphics[height = 40mm]{cxcr3-allosteric-ligands.png}
	\captionof{figure}{Small molecule allosteric modulators of CXCR3}
	\label{fig:cxcr3}
\end{figure*}
\end{fullwidth}

\section{Chemical biology of non-coding RNA}
Non-coding RNA folds into 3D-structures that can be targeted by small molecules.\cite{Bernat2015b}
My research in chemical biology of RNA was about identification of new chemotypes targeting
pathogenic RNA structures and development of chemical tools for studying their role in
cellular markers of disease.
\begin{marginfigure}\label{fig:1a-TOq}
	\includegraphics{1a-TOq.png}
	\caption{Fluorescent reporter for r(G4C2)\textsuperscript{exp} binding assay}
\end{marginfigure}
\newthought{Ellipticines as ligands for $r$(G4C2)\textsuperscript{ext}.} Building upon the previous work in the Disney lab, I developed fluorescent reporter 1a-TO\textsubscript{Q} which served as a tool for quick SAR interrogation \textit{in vitro} and led to discovery of better RNA binders from the same chemotype\cite{Wang2019}. Key insight was that permanent positive charge was not required for potent binding to RNA and activity in cellular assays. This dramatically improved blood-brain barrier penetration in mice.

\newthought{Chemical biology of $r$(CUG)\textsuperscript{exp} and $r$(CCUG)\textsuperscript{exp}.} To study the role of pathogenic non-coding RNAs in myotonic dystrophy I developed series of conjugates that disrupted interaction of r(CUG)\textsuperscript{exp} with proteins (GC-2H-CG, Fig. \ref{fig:kanamycin-dimers}\textcolor{DarkRed}{A}) or targeted the pathogenic r(CCUG)\textsuperscript{exp} for degradation by latent cellular nuclease RNAse L (Fig. \ref{fig:kanamycin-dimers}\textcolor{DarkRed}{B}). Further research in Disney lab built upon this work to develop microRNA degrading tool compounds.

\section{NIR fluorescent glucose sensors}

\begin{marginfigure}\label{fig:profusa}
	\includegraphics{v1dye.png}
	\caption{Example of glucose-sensing polymerizable NIR dye.}
\end{marginfigure}

At Profusa, I developed near-infrared (NIR) fluorescent sensors for continuous sensing of glucose \textit{in vivo}. The design of glucose sensors was based on known bis-boronic anthracene module that was incorporated into biocompatible hydrogel matrix. From a library of 120 dyes that I prepared\cite{Gamsey2018}, I could explore structure-activity relationships. Based on that I was able to rationally improve the performance of the company's prototype dye up to 5-fold. Finally, I scaled the synthesis of the lead dye to produce 1.8 g of material for clinical trials in humans.

\begin{figure*}
	\phantomsection
	A)
	\includegraphics[width=80mm]{GC-2H-CG.png}
	B)
	\includegraphics[width=80mm]{2Kan-p2-5A4.png}
	\caption{A) Structure of peptoid-PNA-Hoechst conjugate for disruption of r(CUG)\textsuperscript{exp} - MBNL1 interaction; B) Structure of peptoid-[2-5]A\textsubscript{4}-Kanamycin conjugate for targeting r(CCUG)\textsuperscript{exp} for degradation by RNAse L}
	\label{fig:kanamycin-dimers}
\end{figure*}

\section{DNA bioconjugation for next generation sequencing}

At Quantapore, I established workflow for immobilization of DNA on flat surfaces and nanoparticles. Via rational optimization of conjugation conditions and surface chemistry I have enabled proof of principle for Quantapore's proprietary sequencing platform.

\bibliography{mypubs}
\bibliographystyle{plainnat}

\end{document}

max_stars_repo_name: slavabernat/researchSummary

max_issues_repo_name: slavabernat/researchSummary

max_forks_repo_name: slavabernat/researchSummary