methods_dev.tex

\chapter{Methods Developed During Thesis}
\label{chap:methods_dev}
%\section{Overview/meta-text}

In this Chapter, I outline the rationale and development of various methods used throughout this thesis to examine \glspl{synthetic lethal} in \gls{gene expression} data, \glslink{graph}{graph} structures, models and simulations. Firstly, the \acrfull{SLIPT}, a \gls{bioinformatics} approach to triage \gls{synthetic lethal} candidate genes, will be described. This is one of the main research outputs of this thesis project and is supported by comparisons to an experimental screen from a related project and evaluation of performance on simulated data. These supporting findings will be covered in further chapters but simulation data is included to support the use and design of \gls{SLIPT}. This includes the construction of a statistical model of \glspl{synthetic lethal} in (continuous multivariate Gaussian) \gls{gene expression} data, which enables testing \gls{SLIPT} upon simulated data with known \gls{synthetic lethal} partners. Another key component of this simulation pipeline is the generation of simulated data from a known \glslink{graph}{graph} structure or simulated biological pathway (as applied in Chapter~\ref{chap:simulation}). The development of this simulation procedure and other statistical treatment of graph and \glslink{graph}{network} structures will also be covered. Various R packages have been developed to support this project, including the \texttt{slipt} package to implement the \gls{SLIPT} methodology. Additional R packages for handling \glslink{graph}{graph} structures, simulations, and custom plotting features will be described as research outputs of this thesis, methods applied throughout, and contributions of open-source software.

\section{A Synthetic Lethal Detection Methodology} \label{methods:SLIPT}
%\subsection{Rationale and Design of Test}
%\subsection{Synthetic Lethal Detection Method}

The \gls{SLIPT} methodology identifies \gls{gene expression} patterns consistent with \gls{synthetic lethal} interactions, between a query gene and a panel of candidate interacting partners. \Gls{gene expression} is scored ``low'', ``medium'', or ``high'', sorting samples by tertiles ($\sfrac{1}{3}$-quantiles) for each gene. Genes with insufficient \glslink{gene expression}{expression} across all samples are excluded by requiring that the first tertile of raw counts is above zero. A $\chi^2$ test is then performed between the query gene and each candidate partner.  The p-values for the $\chi^2$ test are corrected for multiple testing using \acrfull{FDR} error control to reduce false positives \citep{fdr1995}. Significance is called for \gls{FDR} adjusted $p < 0.05$. A \gls{synthetic lethal} interaction is predicted  (as shown in Figure~\ref{fig:SLIPT_Method}) when (i) the $\chi^2$ test is significant; (ii) observed low-query, low-candidate samples are less frequent than expected; and (iii) observed low-query, high-candidate and high-query, low-candidate samples are more frequent than expected.
%The query and candidate genes are swapped to replicate the directional condition. %redundant
%Where \glspl{synthetic lethal} is scored SL-Q if it is predicted in query-low samples and SL-C if it is predicted in candidate-low samples (as shown in Figure~\ref{fig:SLIPT_Method}). \Glspl{synthetic lethal} is only reported in this text if it meets both of these conditions and a significant p-value where it is scored SL-2. %too detailed

\begin{figure*}[!b]
%\begin{mdframed}
\begin{center}
  \resizebox{0.8 \textwidth}{!}{
    \fbox{\input{SL_Method.pdf_tex}}
   }
   \end{center}
   \caption[Framework for \gls{synthetic lethal} prediction]{\small \textbf{Framework for \gls{synthetic lethal} prediction.} \gls{SLIPT} was designed to identify candidate interacting genes from \gls{gene expression} data using the $\chi^2$ test against a query gene. Samples are sorted into low, medium, and high \glslink{gene expression}{expression} quantiles for each gene to test for a directional shift. A sample being low in both genes of a \gls{synthetic lethal} pair is unlikely, since loss of both genes will be deleterious, and is expected to be statistically under-represented in a \gls{gene expression} dataset. \textcolor{black}{It is expected for there to be} a corresponding (symmetric) increase in frequency of sample with low-high gene pairs. \Gls{synthetic lethal} candidate partners of a gene were identified by running this procedure on all possible partner genes, selecting those with an \gls{FDR}-adjusted $\chi^2$-derived $p < 0.05$, and meeting the directional criteria. Since \gls{synthetic lethal} genes are partners of each other, commutatively, the symmetric direction criteria were defined such that detected \gls{synthetic lethal} genes are partners of each other.
}
\label{fig:SLIPT_Method}
%\end{mdframed}
\end{figure*}


The \gls{synthetic lethal} prediction procedure has also been performed with \glslink{somatic}{somatic} \gls{mutation} data for the query gene. This is intended for a \textcolor{black}{known query gene} which is recurrently mutated, with the majority of \glspl{mutation} disrupting gene function (e.g.,  null or frameshift \glspl{mutation}). A \gls{synthetic lethal} interaction is predicted  (as shown in Figure~\ref{fig:SLIPT_Method_mtSL}) when (i) the $\chi^2$ test is significant; (ii) observed \gls{mutant}-query, low-candidate samples are less frequent than expected; and (iii) observed \gls{mutant}-query, high-candidate and \gls{wild-type}-query, low-candidate samples are more frequent than expected. \textcolor{black}{This modified methodology will be referred to as mtSLIPT}. %Unless otherwise specified, computationally predicted \gls{synthetic lethal} gene candidates from \gls{SLIPT} used \glslink{gene expression}{expression} data (exprSL) for both genes (as shown in Figure~\ref{fig:SLIPT_Method}) rather than \gls{mutation} data (mtSL) for the query gene (as shown in Figure~\ref{fig:SLIPT_Method_mtSL}).

\begin{figure*}[!tb]
%\begin{mdframed}
  \begin{center}
  \resizebox{0.8 \textwidth}{!}{
    \fbox{\input{SL_Method_mtSL.pdf_tex}}
   }
   \end{center}
   \caption[Synthetic lethal prediction adapted for \gls{mutation}]{\small \textbf{Synthetic lethal prediction adapted for \gls{mutation}.} \gls{SLIPT} was also adapted to identify candidate interacting genes using (\glslink{somatic}{somatic}) \gls{mutation} data of the query gene in the $\chi^2$ test. Samples are sorted into low, medium, and high \glslink{gene expression}{expression} quantiles for each candidate gene and tested for a directional shift against \gls{mutation} status of the query gene. A sample having low \glslink{gene expression}{expression} or \gls{mutation} for the \gls{synthetic lethal} pair is expected to be unlikely with a corresponding increase in frequency of sample with \gls{mutant}-high or \gls{wild-type}-low gene pairs. \Gls{synthetic lethal} (mtSLIPT) candidate partners of a gene were identified from running this procedure on all possible partner genes, selecting those with an \gls{FDR}-adjusted $\chi^2$-derived $p < 0.05$, and meeting the directional criteria. %Synthetic lethal genes are partners of each other commutatively with \gls{synthetic lethal} genes will predicted to be partners of each other.
}
\label{fig:SLIPT_Method_mtSL}
%\end{mdframed}
\end{figure*}

The \gls{SLIPT} methodology can be performed on \glslink{gene expression}{expression} data, including pathway \glspl{metagene} (as generated in Section~\ref{methods:metagene}). The application of the \gls{SLIPT} methodology on public \gls{gene expression} data will be supported with simulation results (in Section~\ref{chapt2:simulation_2015} and Chapter~\ref{chap:simulation}), including comparison to other statistical methods. \gls{SLIPT} results for \textit{CDH1} were compared \textcolor{black}{to the} experimental screen results in a breast cell line \citep{Telford2015}. Primary screen results are discussed in Section~\ref{chapt3:compare_SL_genes} and secondary screen results are presented in Section~\ref{chapt3:secondary_screen}.

%This methodology was adapted to used pathway \gls{metagene} quantiles rather than \gls{gene expression} as an input for pathway \glspl{synthetic lethal} testing. The p-values for $\chi^2$ tests were also corrected for multiple testing with the false positive rate \citep{fdr1995} across all pathways tested from the same database and with significance defined as a \gls{FDR} adjusted p-values $p < 0.05$ as above.

%mtSLIPT method
%A similar methodology was developed in both cases to test for \glspl{synthetic lethal} where the query gene has an inactivating \gls{mutation} in some patients. Since most \glspl{mutation}, particularly in \gls{tumour suppressor} genes, are deleterious all \glslink{somatic}{somatic} non-synonymous \glspl{mutation} were counted as \gls{mutant} and \glspl{synthetic lethal} was tested with the query gene changed accordingly (as shown in Figure~\ref{fig:mtSLIPT_Method}. To distinguish these methods they are abbreviated to exprSLIPT and \acrshort{mtSLIPT} respectively depending on the molecular property used to define low gene activity of the query gene.

\FloatBarrier

\section{Synthetic Lethal Simulation and Modelling} \label{methods:simulation_SL_expression} 

A statistical model of \glspl{synthetic lethal} was developed to generate simulated data and to evaluate the \gls{SLIPT} procedure. This Section describes the \gls{synthetic lethal} model and the simulation procedure for generating \gls{gene expression} data with known \gls{synthetic lethal} partners. Simulation results, to support usage of the \gls{SLIPT} methodology throughout this thesis, will be presented in Section~\ref{chapt2:simulation_2015}. The simulation procedure will also be applied in Chapter~\ref{chap:simulation}, including in combination with simulations from \glslink{graph}{graph} structures (as described in Section~\ref{methods:graphsim}).

\subsection{A Model of Synthetic Lethality in Expression Data} \label{methods:SL_Model}

A conceptual model of \glspl{synthetic lethal} was devised (as shown in Figure~\ref{fig:SL_Model}), which \textcolor{black}{was} used to build a statistical model of \gls{synthetic lethal} \gls{gene expression} and to simulate \glslink{gene expression}{expression} data for assessing various potential \gls{synthetic lethal} prediction methods, including \gls{SLIPT}. In the model, \glspl{synthetic lethal} occurs between genes with related functions \textcolor{black}{leading to} a cell death phenotype when these functions are \textcolor{black}{inactivated}.

\begin{figure*}[!p]
%\begin{mdframed}
\begin{center}
  \resizebox{0.95 \textwidth}{!}{
    %\input{{{"SL_Model.pdf_tex"}}
    \fbox{
    \includegraphics{{"SL_Model"}}
   }
   }
   \end{center}
   \caption[A model of \gls{synthetic lethal} \gls{gene expression}]{\small \textbf{A model of \gls{synthetic lethal} gene expresion.} A conceptual model of \gls{synthetic lethal} interactions between a Query gene and partner gene ($G_X$). Genes that are \gls{synthetic lethal} may not both be non-functional in the same sample without another gene compensating for the loss of function. This is most likely to be detectable as low \gls{gene expression}, whether they are lost by \gls{mutation}, deletion, \acrshort{DNA} methylation, or suppressing regulatory signals. This could manifest as coexpression, mutual exclusivity, or directional shifts in sample frequency. Thus the alternative hypothesis ($H_{A}$) is that \gls{synthetic lethal} genes will have a reduced frequency of co-loss samples while the null hypothesis ($H_{0}$) is that non-synthetic lethal gene pairs would show no such relationship, even if they may be correlated for other means such as pathway relationships. In this model \gls{synthetic lethal} genes may compensate for the loss of each other but this is not assumed, only that loss of both is unfavourable to cell viability and probability of detecting samples with combined gene loss.
}
\label{fig:SL_Model}
%\end{mdframed}
\end{figure*}


This model suggests that \glspl{synthetic lethal} is detectable in measures of gene inactivation across a sample population, namely \gls{mutation}, \acrshort{DNA} copy number, \acrshort{DNA} \gls{methylation}, and \glslink{gene expression}{expression} levels. While any of these mechanisms of gene inactivation could lead to \glspl{synthetic lethal}, \glslink{gene expression}{expression} data is readily available and changes in other mechanisms are likely to impact on the amount of expressed \acrshort{RNA} that is detectable. Functional relationships between genes could manifest in \glslink{gene expression}{expression} data in several ways, including coexpression, mutual exclusivity and directional shifts. Co-expression is overly simplistic \citep{Lu2015} and has previously performed poorly as a predictor of \glspl{synthetic lethal} \textcolor{black}{with an AUROC of 0.623} \citep{Jerby2014}, although this will still be tested with correlation measures in later simulations. The alternative hypothesis is that \glspl{synthetic lethal} will result in a detectable shift in the number of samples which exhibit low or high \glslink{gene expression}{expression} of either gene. This model does not preclude mutual exclusivity, compensating \glslink{gene expression}{expression}, or co-loss under-representation which may occur between \gls{synthetic lethal} genes \citep{Wappett2016, Lu2015}. 

The first condition of the \gls{synthetic lethal} model is that if there are only two \gls{synthetic lethal} genes (e.g., \textit{CDH1} and one SL partner), then they will not both be non-functional in the same sample (in an ideal model). Gene function is thus determined for each sample in a model of \glspl{synthetic lethal} with the proportion of samples which are functional or non-functional for a gene being arbitrary. Whether a gene is functional can similarly be modelled by an arbitrary threshold of continuous and normally distributed \gls{gene expression} data to define gene function (as shown in Figure~\ref{fig:SL_Model_Expression}). For the purposes of modelling \glspl{synthetic lethal} in cancer \glslink{gene expression}{expression} data, a threshold of the 30\textsuperscript{th} percentile of the \glslink{gene expression}{expression} levels was used because approximately 30\% of samples analysed had \textit{CDH1} inactivation (mutations) in breast cancer \citep{TCGA2012}. This was generalised for a model of the proportion of samples inactivated for each gene. The threshold of the 0.3 quantile was used in simulations \textcolor{black}{derived} from this model throughout this thesis. In this ideal case, no samples lowly expressing both of these genes are expected to be observed. While this was not the case, that is to be expected as it is unlikely that only two genes will have an exclusive \gls{synthetic lethal} partnership.

\begin{figure*}[!tb]
%\begin{mdframed}
  \begin{center}
  \resizebox{0.7125 \textwidth}{!}{
  \fbox{
    \includegraphics{{"SL_Model_Expression"}}
   }
   }
   \end{center}
   \caption[Modelling \gls{synthetic lethal} \gls{gene expression}]{\small \textbf{Modelling \gls{synthetic lethal} \gls{gene expression}.} When modelling \gls{synthetic lethal} interactions between a Query gene and partner genes ($G_X$ and $G_Y$) above,  cellular viability requires that at least of genes is not inactivated. As a model of loss of function, genes are regarded as non-functional with expression below a threshold for the purposes of modelling \glspl{synthetic lethal}. Tumour suppressor genes with loss of function also have cancer specific phenotypes (although these thresholds are not assumed to be the same). Expression is modelled by normally (Gaussian) distributed continuous data, such as (log-scale) data from \acrshort{RNA} (\gls{microarray} or \gls{RNA-Seq}), protein, or pathway \glspl{metagene}. This rationale generalises to several genes on a multivariate normal distribution.
}
\label{fig:SL_Model_Expression} 
%\end{mdframed}
\end{figure*}


\begin{figure*}[!p]
%\begin{mdframed}
  \begin{center}
  \resizebox{0.95 \textwidth}{!}{
    %\input{{{"SL_Model.pdf_tex"}}
    \fbox{
    \includegraphics{{"SL_Model_Higher"}}
   }
   }
   \end{center}
   \caption[Synthetic lethality with multiple genes]{\small \textbf{Synthetic lethality with multiple genes.} Higher order \gls{synthetic lethal} interactions may occur between 3 or more genes, affecting the simulated \glslink{gene expression}{expression} (or \gls{synthetic lethal} predictions) even if undetected when observed pairwise. Consider interactions between a Query gene and two partner genes ($G_X$ and $G_Y$). They may interact with the Query pairwise (inviable when either gene pair is lost) or form a higher-order interaction such as the ``synthetic lethal triplet''  if any of the genes provide an \gls{essential} function (inviable only when all are lost). Either is plausible with the potential \glslink{graph}{pathway} structures. A \gls{synthetic lethal} triple has 8 potential combinations of gene function but one is not expected to be observed (due to inviability), however pairwise inactivation may be observed if additional partner genes are functional. The proportion of these combinations varies depending on the functional threshold.
}
\label{fig:SL_Model_Higher}
%\end{mdframed}
\end{figure*}



A \gls{synthetic lethal} pair of genes is unlikely to act in isolation, therefore higher-order \gls{synthetic lethal} interactions (i.e., 3 or more genes) must be considered in the model as shown in Figure~\ref{fig:SL_Model_Higher}. Even when testing pairwise interactions, it is important to model higher level interactions that may interfere. If there are additional \gls{synthetic lethal} partners, there are two possibilities for adding these: 1) that they are independent partners of the query genes interacting pairwise (and not with each other) or 2) that an additional partner gene interacts with both of the \gls{synthetic lethal} genes already in the system and any of the three (or more) are required to be functional for the cell to survive.

The signal (in terms of \gls{gene expression} data) will be weaker for this latter case and this model has the more stringent assumption that all \gls{synthetic lethal} partner genes interact with each other: that only one of these must be expressed to satisfy the model of \glspl{synthetic lethal}. In this model, any of the \gls{synthetic lethal} genes in a higher-order interaction are able to perform the essential function of the others, allowing for higher-level \gls{synthetic lethal} partners to compensate for loss \textcolor{black}{of} a \gls{synthetic lethal} gene pair. While samples that express low levels of the \gls{synthetic lethal} gene pairs will be under-represented, they may not be completely absent from the dataset, due to these higher-level interactions. In the example of three \gls{synthetic lethal} genes (shown in Figure~\ref{fig:SL_Model_Higher}), only one of the genes involved in the higher-order \gls{synthetic lethal} interaction is required for cell viability. For \gls{synthetic lethal} pairs, only a subset of these samples will be inviable (i.e., removed from simulated data), leading to an under-representation.

Samples were not actually removed from a simulated dataset, rather the \glslink{gene expression}{expression} and function of the query gene is generated across samples separately from the pool of potential partner genes. The query gene data was matched to simulated samples (as shown in Figure~\ref{fig:simulate_add_query}) satisfying the \gls{synthetic lethal} condition with the procedure described in Section~\ref{methods:simulating_SL}. This was performed to maintain \textcolor{black}{a comparable sample size} across simulations and \textcolor{black}{to} preserve the (multivariate) normal distribution of the data. 

\FloatBarrier

\subsection{Simulation Procedure} \label{methods:simulating_SL}

Simulations were developed to generate normal distributions of \glslink{gene expression}{expression} data and define gene function with a threshold cut-off. 
%This is the reverse to the procedure of \gls{SLIPT} to predict \gls{synthetic lethal} partners (although the threshold is assumed to be unknown when testing upon simulated data). 
While gene function was used as an intermediary step in modelling \gls{synthetic lethal} genes in \glslink{gene expression}{expression} data, the normal distribution was sampled for simulated data to represent normalised empirical \gls{gene expression} data for which \gls{SLIPT} (and other methods) will be applicable.

\begin{figure*}[!htbp]
%\begin{mdframed}
%  \resizebox{\textwidth}{!}{
         \begin{center}
%
        \subcaptionbox{Simulated \glslink{gene expression}{expression} matrix}{%
            \label{fig:simulate_function:first}
            %\includegraphics[width=0.5\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_mat_inhibiting".png}}
            \includegraphics[width=0.45\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_mat".png}}
        }%
        \subcaptionbox{Corresponding gene function calls}{%
           \label{fig:simulate_function:second}
           %\includegraphics[width=0.5\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_inhib_disc_mat".png}} %%check if same tree order (sample) as \glslink{gene expression}{expression}
           \includegraphics[width=0.45\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_disc_mat".png}}
        }%
%
    \end{center}
   %\caption[Simulating \gls{gene expression} and function]{\small \textbf{\textbf{Simulating gene function.}} Simulated data with samples (columns) and genes A--I (rows) shows how a simulated dataset is transformed from a continuous dataset (on a blue to red colour scale) to a discrete matrix of gene function (samples with functional gene levels are shaded in black and non-functional in grey).}
   \caption[Simulating gene function]{\small \textbf{\textbf{Simulating gene function.}} A simulated dataset with samples (columns) and genes A--I (rows) was transformed from a continuous (coloured blue--red) scale to a discrete matrix of gene function (black for functional levels and grey for non-functional).}
%}
\label{fig:simulate_function}
%\end{mdframed}

\iffalse
%\begin{mdframed}
%  \resizebox{\textwidth}{!}{
         \begin{center}
%
	\subcaptionbox{Simulated gene function with SL genes}{%
            \label{fig:simulate_add_query:first}
            %\includegraphics[width=0.5\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_inhib_SL_disc_mat".png}}
            %\includegraphics[width=0.5\textwidth,trim=4cm 2cm 0cm 0cm,clip]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_SL_disc_mat".png}}
            \includegraphics[width=0.45\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_SL_disc_mat".png}}
        }%
        \subcaptionbox{Query gene added with SL condition}{%
           \label{fig:simulate_add_query:second}
           %\includegraphics[width=0.5\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_inhib_disc_query_mat_graph".png}} %%check if same tree order (sample) as \glslink{gene expression}{expression}
           %\includegraphics[width=0.5\textwidth,trim=4cm 2cm 0cm 0cm,clip]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_disc_query_mat_graph".png}}
           \includegraphics[width=0.45\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_disc_query_mat_graph".png}}
        }%
%
    \end{center}
   %\caption[Simulating \gls{synthetic lethal} gene function]{\small \textbf{\textbf{Simulating \gls{synthetic lethal} gene function.}} Simulated data with samples (columns) and genes (rows) in a discrete matrix of gene function (shaded in black for sample with functional gene levels). Genes A and I are selected to be \gls{synthetic lethal} partners of a ``Query'' gene, which of these genes will be the true partner in each sample is selected randomly and indicated in green which samples are considered for the purposes of simulating \glspl{synthetic lethal} (shaded in forest green for samples with functional gene levels). Note that samples are ordered such that either the query gene or selected partner are functional in any particular sample.}
   \caption[Simulating \gls{synthetic lethal} gene function]{\small \textbf{\textbf{Simulating \gls{synthetic lethal} gene function.}} In a discrete simulated gene function dataset (black for functional levels and grey for non-functional) with samples (columns) and genes (rows), genes A and I are the SL partners of a ``Query'' gene. A partner  gene is selected randomly (shown in green) in each sample for simulating \glspl{synthetic lethal} (forest green for functional genes). %Note that samples are ordered such that either the query gene or selected partner are functional in any particular sample.
   }
%}
\label{fig:simulate_add_query}
%\end{mdframed}
\fi
%\end{figure*}
%\begin{figure*}[!ht]
%\begin{mdframed}
%  \resizebox{\textwidth}{!}{
         \begin{center}
%
	\subcaptionbox{Simulated gene function with SL genes}{%
            \label{fig:simulate_add_query:first}
            %\includegraphics[width=0.5\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_inhib_SL_disc_mat".png}}
            %\includegraphics[width=0.5\textwidth,trim=4cm 2cm 0cm 0cm,clip]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_SL_disc_mat".png}}
            \includegraphics[width=0.45\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_SL_disc_mat".png}}
        }%
        \subcaptionbox{Query gene added with SL condition}{%
           \label{fig:simulate_add_query:second}
           %\includegraphics[width=0.5\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_inhib_disc_query_mat_graph".png}} %%check if same tree order (sample) as \glslink{gene expression}{expression}
           %\includegraphics[width=0.5\textwidth,trim=4cm 2cm 0cm 0cm,clip]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_disc_query_mat_graph".png}}
           \includegraphics[width=0.45\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_disc_query_mat_graph".png}}
        }%
%
    \end{center}
   %\caption[Simulating \gls{synthetic lethal} gene function]{\small \textbf{\textbf{Simulating \gls{synthetic lethal} gene function.}} Simulated data with samples (columns) and genes (rows) in a discrete matrix of gene function (shaded in black for sample with functional gene levels). Genes A and I are selected to be \gls{synthetic lethal} partners of a ``Query'' gene, which of these genes will be the true partner in each sample is selected randomly and indicated in green which samples are considered for the purposes of simulating \glspl{synthetic lethal} (shaded in forest green for samples with functional gene levels). Note that samples are ordered such that either the query gene or selected partner are functional in any particular sample.}
   \caption[Simulating \gls{synthetic lethal} gene function]{\small \textbf{\textbf{Simulating \gls{synthetic lethal} gene function.}} In a discrete simulated gene function dataset (shaded for functional levels and pale otherwise) with samples (columns) and genes (rows), genes A and I were SL partners of a ``Query'' gene. A partner was selected (highlighted in green) randomly in each sample for simulating \glspl{synthetic lethal}, then ordered such that the query gene or an SL partner were functional in each sample.
   }
%}
\label{fig:simulate_add_query}
%\end{mdframed}
\end{figure*}


Sampling a distribution for \glslink{gene expression}{expression} profiles has the advantage of enabling simulating correlation structures with the multivariate normal distribution, using the \texttt{mvtnorm} R package \citep{Genz2009, mvtnorm}. The parameter $\Sigma$, the covariance matrix \textcolor{black}{($\Sigma_{ij} = \sigma_{X_iX_j}$)}, defines the correlation structure between the simulated genes \textcolor{black}{(rows $i$ and $j$ in data matrix $X$)} being sampled. With a diagonal of one \textcolor{black}{($\Sigma_{ii} = \sigma^{2}_{X_i}= 1$)}, this $\Sigma$ matrix simulates genes with a standard deviation of one and the covariance parameters between them are the correlations between each gene. In Figure~\ref{fig:simulate_function}, an example of such a simulated multivariate normal dataset is shown with the functional threshold applied.

Once a simulated dataset has been generated, the samples were compared by gene function (as derived from a functional threshold). The known underlying \gls{synthetic lethal} partners were selected within the dataset and a query gene was generated by sampling from the normal distribution. These were matched (as shown for two \gls{synthetic lethal} partners in Figure~\ref{fig:simulate_add_query}) such that the \gls{synthetic lethal} condition was met: at least one of the \gls{synthetic lethal} partner genes and the query gene are functional in any particular cell. The samples are ordered by functional data (without assuming correlation of underyling \glslink{gene expression}{expression} values) with the query gene in one direction and the remaining dataset ordered by the selected \gls{synthetic lethal} partner.


\begin{figure*}[!htb]
%\begin{mdframed}
%  \resizebox{\textwidth}{!}{
         \begin{center}
%
        \subcaptionbox{Initial \glslink{gene expression}{expression} matrix}{%
            \label{fig:simulate_SL:first}
            %\includegraphics[width=0.5\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_mat_inhibiting".png}}
            \includegraphics[width=0.45\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_mat".png}}
        }%
        \subcaptionbox{Simulated \gls{synthetic lethal} dataset}{%
           \label{fig:simulate_SL:second}
           %\includegraphics[width=0.5\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_inhib_query_mat_graph".png}} %%check if same tree order (sample) as \glslink{gene expression}{expression}
           \includegraphics[width=0.45\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_query_mat_graph".png}}
        }%
%
    \end{center}
   %\caption[Simulating \gls{synthetic lethal} \gls{gene expression}]{\small \textbf{\textbf{Simulating \gls{synthetic lethal} \gls{gene expression}.}} Simulated data with samples (columns) and genes (rows) showing how a simulated continuous dataset (on a blue to red colour scale) is matched to a query gene such that at least one \gls{synthetic lethal} partner is above a functional threshold when the query gene is below it satisfying the \gls{synthetic lethal} model.}
   \caption[Simulating \gls{synthetic lethal} \gls{gene expression}]{\small \textbf{\textbf{Simulating \gls{synthetic lethal} \gls{gene expression}.}} A simulated continuous \glslink{gene expression}{expression} dataset (blue--red scale) with samples (columns) and genes A--I (rows) was matched to a query gene such that at least one \gls{synthetic lethal} partner was above a functional threshold when the query gene was below it which satisfied the \gls{synthetic lethal} model.}
%}
\label{fig:simulate_SL}
%\end{mdframed}
\end{figure*}

This procedure produces a simulated dataset where samples with a non-functional query gene have at least one functional partner gene. Similarly, the query gene is functional in all samples where all of the \gls{synthetic lethal} partner genes are non-functional. In this procedure, a dataset has been generated with known \gls{synthetic lethal} partners (as shown in Figure~\ref{fig:simulate_SL}) with few assumptions about the relationships between the each \gls{synthetic lethal} pair (allowing compensating functions from higher-order interactions). This procedure has been designed to have the most stringent (least detectable) \gls{synthetic lethal} relationships, where higher-order interactions are possible for the purposes of testing pairwise detection procedures such as \gls{SLIPT}.  


\FloatBarrier

\section{Detecting Simulated Synthetic Lethal Partners} \label{chapt2:simulation_2015}

The \gls{synthetic lethal} detection methodology (\gls{SLIPT}), as described in Section~\ref{methods:SLIPT}, was evaluated with simulated data containing known \gls{synthetic lethal} partners, generated using the procedure described in Section~\ref{methods:simulating_SL}. Simulations were performed to demonstrate the methodology and support its use throughout this thesis. These simulations were performed by sampling from statistical distributions, including the multivariate normal distribution with correlated blocks of genes, generated by $\Sigma$ matrices such as those shown. A more complex multivariate normal sampling procedure based on pathway \glslink{graph}{graph} structures, as described in Section ~\ref{methods:graphsim}, was used for further investigations in Chapter~\ref{chap:simulation}. 

\subsection{Binomial Simulation of Synthetic Lethality} \label{chapt2:simulation_binom}
%[relevant?]

The \gls{synthetic lethal} simulation procedure (described in Section~\ref{methods:simulating_SL}) initially used gene function, sampled directly from a binomial distribution using the binomial probability of observing functional gene levels ($p = 0.3$) in one observation ($n = 1$) for each samples: $$X\sim Binomial(n,p)$$ \textcolor{black}{for a binary matrix $X$ of gene function with $n$ trials and a probability of $p$}. Once a query gene with \gls{synthetic lethal} partners has been added, these functional levels were passed directly into \gls{SLIPT} as ``low'' and ``high'' samples.

The simulation procedure was performed with 20,000 total genes (as occurs in \glslink{gene expression}{expression} datasets) with a variable number of true \gls{synthetic lethal} partners and 500, 1000, 2000, or 5000 samples. Each \gls{ROC} curve was derived from the results of 10,000 replicate simulations. The statistical performance (as shown in Figure~\ref{fig:Binomial_AUC}) of the $\chi^2$-derived p-value declined towards random predictions (an \gls{AUROC} of 0.5) with \textcolor{black}{more} underlying \gls{synthetic lethal} partners to detect. However, increased sample size somewhat mitigated this decline, as expected with a statistical predictor, particularly for moderate numbers of \gls{synthetic lethal} partners. 

\begin{figure*}[!hp]
%\begin{mdframed}
  \begin{center}
  \resizebox{0.66 \textwidth}{!}{
  %\fbox{
    \includegraphics[width=0.66\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/RUN_20150311/SL_Model_Binomial_1K_AUC_samples_prop".png}}
   %}
   }
   \end{center}
   \caption[Performance of binomial simulations]{\small \textbf{Performance of binomial simulations.} Gene function was simulated by binomial sampling and tested for \glspl{synthetic lethal} by \gls{SLIPT}. Statistical performance declined with additional \gls{synthetic lethal} partners \textcolor{black}{and lower} sample sizes.}
\label{fig:Binomial_AUC}
%\end{mdframed}

%\begin{mdframed}
  \begin{center}
  \resizebox{1 \textwidth}{!}{
  %\fbox{
    \includegraphics[width=1\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/RUN_20150410y/SL_Model_Test_Graph_10K_Graph1_ROC_Compare_Binom(Feb)_v_Mvtn(Aprxy)_Full".png}}
   %}
   }
   \end{center}
   \caption[Comparison of statistical performance]{\small \textbf{Comparison of statistical performance.} Binomial simulation of \glspl{synthetic lethal} (in colour) in comparison to multivariate normal simulations (in greyscale), in which \gls{SLIPT} consistently had higher performance across parameters.}
\label{fig:Binomial_Compare}
%\end{mdframed}
\end{figure*}

Simulations using this binomial model of \glspl{synthetic lethal} were simplistic but informed the development of more complex simulations including \glslink{gene expression}{expression} and correlation structures. It did not represent the data that \gls{SLIPT} will be applied to but the binomial simulations demonstrated that \gls{SLIPT} is able to distinguish small numbers of \gls{synthetic lethal} partners in a simulated system with behaviours expected with respect to sample size. This supported further development of the \gls{synthetic lethal} model and simulation pipeline (as described in Section~\ref{methods:simulation_SL_expression}) using the multivariate normal distribution.

The multivariate normal simulation procedure more closely recapitulates the (normalised) \glslink{gene expression}{expression} data that \gls{SLIPT} was intended for and enables the methodology procedure to be tested without requiring modifications (in Section~\ref{chapt2:simulation_mvtnorm}). Sampling continuous \glslink{gene expression}{expression} values from a normal distribution enabled the \glslink{gene expression}{expression} threshold for gene function to differ from the categorical ``low'' and ``high'' \glslink{gene expression}{expression} binning performed by \gls{SLIPT} (as discussed in Section~\ref{methods:SL_Model}). The \gls{SLIPT} procedure does not assume a known threshold for \glslink{gene expression}{expression} and instead uses \glslink{gene expression}{expression} as an estimate of gene function which does not compromise the statistical performance of the \gls{SLIPT} in the multivariate normal simulation. The performance was an improvement over the binomial simulation procedure (shown in Figure~\ref{fig:Binomial_Compare}) across simulation parameters in an equivalent simulation (without correlation structure). This multivariate normal model is also more refined since it defines the \gls{synthetic lethal} condition, to ensure that at least one \gls{synthetic lethal} partner was active in query-deficient samples, without disrupting the proportion of samples with each gene being functional.

\FloatBarrier

\subsection{Multivariate Normal Simulation of Synthetic Lethality} \label{chapt2:simulation_mvtnorm}

The multivariate normal simulation procedure was initially performed using the \texttt{mvtnorm} R package \citep{Genz2009, mvtnorm} (as described in Section~\ref{methods:simulation_SL_expression}) without correlation structure. Expression \textcolor{black}{data ($X$)} was sampled from multivariate normal distribution with a mean ($\mu = 0$), standard deviation ($\sigma = 1$), and no correlation between genes ($r = 0$): $$X\sim \textcolor{black}{Normal}(\bar{\mu},\Sigma).$$ Once a query gene with \gls{synthetic lethal} partners has been added, the simulated \glslink{gene expression}{expression} values were tested by \gls{SLIPT}, as described in Section~\ref{methods:SLIPT}.

\begin{figure*}[!hp]
%\begin{mdframed}
%  \resizebox{\textwidth}{!}{
         \begin{center}
%
        \subcaptionbox{Statistical evaluation \label{fig:simulation_Apr15ROC:Perf}}{%
            \includegraphics[width=0.475\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/RUN_20150410y/SL_Model_Apr15mvnormCor_1K_ROC1_samplesx_prop".png}}
        }%
        \subcaptionbox{\gls{ROC} curves\label{fig:simulation_Apr15ROC:ROC}}{%
            \includegraphics[width=0.475\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/RUN_20150410y/SL_Model_Apr15mvnormCor_1K_ROC2_samplesx_prop".png}}
        }%
        
        \subcaptionbox{Statistical performance (\acrshort{AUROC}) \label{fig:simulation_Apr15ROC:AUC}}{%
           \includegraphics[width=0.65\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/RUN_20150410y/SL_Model_Apr15mvnormCor_1K_AUC_samplesx_prop".png}}
        }%
    \end{center}
   \caption[Performance of multivariate normal simulations]{\small \textbf{Performance of multivariate normal simulations.} Simulation of \glspl{synthetic lethal} was performed by sampling from a multivariate normal distribution (without correlation structure). Performance of \gls{SLIPT} declined with increasing numbers of \gls{synthetic lethal} partners \textcolor{black}{and lower} sample sizes (in darker colours). This occurred as the sensitivity decreased with a greater number of true positives to detect, which lead to a trade-off in accuracy as seen in a trough for false positive rate and the \gls{ROC} curves.}
%}
\label{fig:simulation_Apr15ROC}
%\end{mdframed}
\end{figure*}

The statistical accuracy of \gls{SLIPT} as a binary classifier was high across simulations of a full dataset of 20,000 genes (shown in Figure~\ref{fig:simulation_Apr15ROC:Perf}). Using the $\chi^2$-derived p-value as a threshold for prediction, this was largely due to high specificity: the majority of non-synthetic lethal genes were distinguished from the underlying \gls{synthetic lethal} genes. Thus the \gls{SLIPT} methodology performed better with larger datasets with more expected negatives and the results of simulations of smaller numbers of genes (e.g.,  the \glslink{graph}{graph} structures analysed in Section~\ref{chapt5:graphsim_performance}) can be applied to larger datasets, where they are expected to perform comparably or better with a lower false negative rate (as shown in Sections~\ref{chapt5:graphsim_performance_20K} and~\ref{chapt5:graphsim_performance_20K_pway}). Accordingly, key results will be supported by replication with larger numbers of non-synthetic lethal genes added to the simulations. 

The sensitivity of \gls{SLIPT} as a binary classifier of \glspl{synthetic lethal} (as shown in Figure~\ref{fig:simulation_Apr15ROC:Perf}) declined with higher numbers of \gls{synthetic lethal} genes to detect, although this is somewhat mitigated by higher sample sizes. The minority of true \gls{synthetic lethal} partners are more difficult to distinguish when there are more of them (with a weaker \glslink{gene expression}{expression} signal from each). While a reduction of the false positive rate could be achieved for moderate numbers of underlying \gls{synthetic lethal} partners, the number of partners to be detected in analyses of \glslink{gene expression}{expression} data is unknown. However, this simulation procedure is amenable to assessing the performance of \gls{SLIPT} across simulation parameters, \glslink{graph}{graph} structures, and comparisons to other approaches (as presented in Chapter~\ref{chap:simulation}).

Not all of the genes detected by \gls{SLIPT} were true \gls{synthetic lethal} genes but they were among the strongest candidates and \gls{SLIPT} had higher performance with fewer underlying \gls{synthetic lethal} genes to detect. These results support a focus on pathway analyses, in particular, the selection of pathways for further investigation. Pathway over-representation analysis was performed to detect functional groups recurrently detected by \gls{SLIPT} since individually detected gene candidates were not necessarily \gls{synthetic lethal}. The detection of functionally related genes (in Chapter~\ref{chap:SLIPT}) supports the role of a pathway in \gls{synthetic lethal} relationships. The use of pathway \glspl{metagene} can reduce the number of potential pathways, compared to genes, to help identify \glspl{synthetic lethal}. These approaches were both applied in Chapter~\ref{chap:SLIPT} to identify the \gls{synthetic lethal} pathways of \textit{CDH1}. Pathways are also more likely to replicate across experimental models, as demonstrated by \citet{Dixon2008}.
 
The \gls{ROC} curves showed that \gls{SLIPT} is subject to a near equal trade-off between sensitivity and specificity across threshold values (as shown in Figure~\ref{fig:simulation_Apr15ROC:ROC}). The lower sensitivity and higher specificity with a binary classification (as shown in Figure~\ref{fig:simulation_Apr15ROC:Perf}) results from stringent testing by \gls{SLIPT} with \gls{FDR} adjusted p-values. The area under these curves (\gls{AUROC}) was used to compare statistical performance (as shown in Figure~\ref{fig:simulation_Apr15ROC:AUC}), which had lower performance for more underlying \gls{synthetic lethal} partners, and higher performance for larger sample size in multivariate normal simulations.

\FloatBarrier

\subsubsection{Multivariate Normal Simulation with Correlated Genes} \label{chapt2:simulation_mvtnorm_cor}
%\subsubsection{Simulated Expression Heatmaps}

\begin{figure*}[!hp]
%\begin{mdframed}
%  \resizebox{\textwidth}{!}{
         \begin{center}
%
        \subcaptionbox{Input $\Sigma$ matrix parameter  \label{fig:simulation_May4SL:first}}{%
            \includegraphics[width=0.35\textwidth]{{"SL_Model_May15mvnorm_heatmap_4SL_cor_comp_top(1)".pdf}}
        }%
        \subcaptionbox{Simulated correlation matrix \label{fig:simulation_May4SL:second}}{%
            \includegraphics[width=0.35\textwidth]{{"SL_Model_May15mvnorm_heatmap_4SL_cor_comp_top(2)".pdf}}
        }%
        
        \subcaptionbox{Simulated \gls{gene expression}  \label{fig:simulation_May4SL:third}}{%
           \includegraphics[width=0.35\textwidth]{{"SL_Model_May15mvnorm_heatmap_4SL_cor_comp_top(4)".pdf}}
        }%
	\subcaptionbox{Simulated gene function \label{fig:simulation_May4SL:fifth}}{%
           \includegraphics[width=0.35\textwidth]{{"SL_Model_May15mvnorm_heatmap_4SL_cor_comp_top(3)".pdf}}
        }%
    \end{center}
   \caption[Simulating \glslink{gene expression}{expression} with correlated gene blocks]{\small \textbf{Simulating \glslink{gene expression}{expression} with correlated gene blocks.} A $\Sigma$ matrix (a) \textcolor{black}{was used to generate} 100 genes with a multivariate normal distribution, including correlated blocks of genes ($r = 0.8$) with correlation (b) similar to $\Sigma$, on a red--to--green scale. The annotation for genes gives the $\chi^2$ (in blue for in the direction of \gls{SLIPT} or red otherwise) and the gene category (blue for query, cyan for query-correlated, red for SL, orange for SL-correlated, forest green for non-synthetic lethal-correlated, and green for non-synthetic lethal). The simulated \gls{gene expression} (c) and function (d) generated were ordered by $\chi^2$ showing the functional structure of \gls{synthetic lethal} genes and that they were among the strongest \gls{SLIPT} results.}
%}
\label{fig:simulation_May4SL}
%\end{mdframed}
\end{figure*}

\begin{figure*}[!htp]
%\begin{mdframed}
%  \resizebox{\textwidth}{!}{
         \begin{center}
%
        \subcaptionbox{Input $\Sigma$ matrix parameter}{%
            \label{fig:simulation_May4SL1K:first}
            \includegraphics[width=0.35\textwidth]{{"SL_Model_May15mvnorm_heatmap_4SL_cor_comp_1K_top(1)".pdf}}
        }%
        \subcaptionbox{Simulated correlation matrix}{%
            \label{fig:simulation_May4SL1K:second}
            \includegraphics[width=0.35\textwidth]{{"SL_Model_May15mvnorm_heatmap_4SL_cor_comp_1K_top(2)".pdf}}
        }%
        
        \subcaptionbox{Simulated \gls{gene expression}}{%
           \label{fig:simulation_May4SL1K:third}
           \includegraphics[width=0.35\textwidth]{{"SL_Model_May15mvnorm_heatmap_4SL_cor_comp_1K_top(4)".pdf}}
        }%
	\subcaptionbox{Simulated gene function}{%
           \label{fig:simulation_May4SL1K:fifth}
           \includegraphics[width=0.35\textwidth]{{"SL_Model_May15mvnorm_heatmap_4SL_cor_comp_1K_top(3)".pdf}}
        }%
    \end{center}
   \caption[Simulating \glslink{gene expression}{expression} with correlated gene blocks]{\small \textbf{Simulating \glslink{gene expression}{expression} with correlated gene blocks.} Using the $\Sigma$ matrix~(a), sampling 1000 genes from a multivariate normal distribution produced (b) correlated blocks of genes (correlated by \textcolor{black}{$r=0.8$}) on a red--to--green scale. The simulated \gls{gene expression} (c) and function (d) generated were ordered by $\chi^2$ and \gls{SLIPT} direction show that \gls{synthetic lethal} genes are among the strongest \gls{SLIPT} results with high specificity against many potential false positives. These are annotated for log-$\chi^2$ (on a red--to--green scale) and category (blue for query, cyan for query-correlated, red for SL, orange for SL-correlated, forest green for non-synthetic lethal-correlated, and green for non-synthetic lethal) for each gene.}
%}
\label{fig:simulation_May4SL1K}
%\end{mdframed}
\end{figure*}

Correlation structures were added to the simulation procedure (with the $\Sigma$ matrix, as discussed in Section~\ref{methods:simulation_SL_expression}), using correlated blocks of genes (as shown in Figure~\ref{fig:simulation_May4SL:first}). These correlated blocks represent genes with correlated \glslink{gene expression}{expression}, such as co-regulation or shared membership of biological pathways. The example (as shown in Figure~\ref{fig:simulation_May4SL}) shows four \gls{synthetic lethal} genes (out of 100), each with five correlated genes that are not themselves \gls{synthetic lethal} partners of the query gene. These simulations address whether \gls{synthetic lethal} genes are distinguishable from correlated partners. The $\Sigma$ matrix produced a similar correlation structure (Figure~\ref{fig:simulation_May4SL:second}) and \glslink{gene expression}{expression} profiles (Figure~\ref{fig:simulation_May4SL:third}).  Apart from correlated blocks of genes ($r = 0.8$), the remaining genes had small variations due to random sampling. The structure of the dataset, particularly between \gls{synthetic lethal} genes and the query, was evident in the simulated \gls{gene expression} (Figure~\ref{fig:simulation_May4SL:third}) and function (Figure~\ref{fig:simulation_May4SL:fifth}). When these genes were ordered by the \gls{SLIPT} results, the \gls{synthetic lethal} genes were highly ranked \textcolor{black}{and} the majority of them were distinguishable from highly correlated genes.

The use of correlation structure was applied to larger datasets, such as the 1000 genes shown in Figure~\ref{fig:simulation_May4SL1K}. \Gls{synthetic lethal} genes were highly ranked by \gls{SLIPT} and were often distinguishable from correlated genes. As previously discussed in Section~\ref{chapt2:simulation_mvtnorm}, these \gls{synthetic lethal} genes were still detectable among a larger number of non-synthetic lethal genes, and the \gls{SLIPT} methodology performed better on large datasets.


\begin{figure*}[!tb]
%\begin{mdframed}
%  \resizebox{\textwidth}{!}{
         \begin{raggedleft}
%
        \subcaptionbox{Gene category in simulations}{%
            \label{fig:simulation_May4SLreps:first}
            \includegraphics[width=0.475\textwidth]{{"SL_Model_May15mvnorm_heatmap_10XSL_cor_comp2(1)".pdf}}
        }%
        \subcaptionbox{Corresponding $\chi^2$ values}{%
            \label{fig:simulation_May4SLreps:second}
            \includegraphics[width=0.475\textwidth]{{"SL_Model_May15mvnorm_heatmap_10XSL_cor_comp2(2)".pdf}}
        }%

    \end{raggedleft}
   \caption[Synthetic lethal prediction across simulations]{\small \textbf{Synthetic lethal prediction across simulations.} The gene category (a) ordered by $\chi^2$ and the \gls{SLIPT} directional condition is shown across simulations (blue for query, cyan for query-correlated, red for SL, orange for SL-correlated, forest green for non-synthetic lethal-correlated, and green for non-synthetic lethal). For each number (1--10) of \gls{synthetic lethal} partners, 10 simulations show that the increasing numbers of \gls{synthetic lethal} partners became harder detect (i.e., red cells become interspersed in the columns of (a)). The log-$\chi^2$ values on a blue--red scale (b) showed a threshold for \gls{synthetic lethal} and correlated genes when there are fewer of them, distinguishable from correlated genes in this case.}
%}
\label{fig:simulation_May4SLreps}
%\end{mdframed}
\end{figure*}


%5 correlated genes
 \begin{figure*}[!hp]
%\begin{mdframed}
%  \resizebox{\textwidth}{!}{
         \begin{center}
%
        \subcaptionbox{Statistical evaluation}{%
            \label{fig:simulation_Apr15ROC2:Perf}
            \includegraphics[width=0.475\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/RUN_20150410y/SL_Model_Apr15mvnorm_1K_ROC1_samplesx_prop".png}}
        }%
        \subcaptionbox{\gls{ROC} curves}{%
            \label{fig:simulation_Apr15ROC2:ROC}
            \includegraphics[width=0.475\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/RUN_20150410y/SL_Model_Apr15mvnorm_1K_ROC2_samplesx_prop".png}}
        }%
        
        \subcaptionbox{Statistical performance (\acrshort{AUROC}) }{%
           \label{fig:simulation_Apr15ROC2:AUC}
           \includegraphics[width=0.65\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/RUN_20150410y/SL_Model_Apr15mvnorm_1K_AUC_samplesx_prop".png}}
        }%
    \end{center}
   \caption[Performance with correlations]{\small \textbf{Performance with correlations.} Simulation of \glspl{synthetic lethal} was performed by sampling from a multivariate normal distribution (with correlation structure). Performance of \gls{SLIPT} declines for more \gls{synthetic lethal} partners but this is mitigated by increased sample sizes (darker colours). This generally occurs as the sensitivity decreases for a greater number of true positives to detect, leading to a trade-off in accuracy as seen in a trough for false positive rate and the \gls{ROC} curves.}
%}
\label{fig:simulation_Apr15ROC2}
%\end{mdframed}
\end{figure*}


These plots (Figures~\ref{fig:simulation_May4SL} and~\ref{fig:simulation_May4SL1K}) used similar correlated blocks with a non-synthetic lethal gene (true negative) and the query gene (which is not \gls{synthetic lethal} with itself). Neither of these were \gls{synthetic lethal} but they could potentially affect performance methodology, particularly the specificity, as correlated non-synthetic lethal genes may be distinguishable from \gls{synthetic lethal} genes. The non-synthetic lethal correlated block of genes had no impact on \gls{synthetic lethal} detection but the query correlated genes were important (as shown in Sections~\ref{chapt2:simulation_mvtnorm_query_cor} and and~\ref{chapt5:compare_chisq_query_cor}).

%\subsection{Replication Simulation Heatmap}
The simulations of \gls{gene expression} data (with 100 genes) with correlations structure were used to examine the variation between detection in different samples and varying the number of underlying \gls{synthetic lethal} partners. A small number of simulations (10 for each) are shown to demonstrate the variation between replicate simulations from iterative sampling from the same multivariate normal distribution (as shown in Figure~\ref{fig:simulation_May4SLreps}). These simulations showed that \gls{synthetic lethal} genes were highly ranked by \gls{SLIPT} when there are few of them and these were relatively consistent across replicate simulations. However, they were less consistent for higher numbers of \gls{synthetic lethal} partners to detect and were more difficult to distinguish from other genes, particularly those correlated with them. Similarly, the $\chi^2$ values showed clear thresholds for \gls{synthetic lethal} and correlated genes in simple simulations but these were more gradual for higher numbers of \gls{synthetic lethal} partners.

\begin{figure*}[!htb]
%\begin{mdframed}
  \begin{center}
  \resizebox{1 \textwidth}{!}{
  %\fbox{
    \includegraphics{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/RUN_20150410y/SL_Model_Test_Graph_10K_Graph1_ROC_Compare_Mvtn(Apry)_v_Cor(Aprxy)_Full"}}
   %}
   }
   \end{center}
   \caption[Comparison of statistical performance with correlation structure]{\small \textbf{Comparison of statistical performance with correlation structure.} Multivariate simulation of \glspl{synthetic lethal} with correlation structure (in colour) has comparable performance to simulation without correlations (in greyscale)  with known \gls{synthetic lethal} partners across parameters.}
\label{fig:mvtnorm_cor_compare}
%\end{mdframed}
\end{figure*}

%query correlated
%5 correlated genes
 \begin{figure*}[!htbp]
%\begin{mdframed}
%  \resizebox{\textwidth}{!}{
         \begin{center}
%
        \subcaptionbox{Statistical evaluation \label{fig:simulation_May15ROC:Perf}}{%
            
            \includegraphics[width=0.475\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/RUN_20150507y/SL_Model_May15mvnorm_dir_1K_ROC1_samplesx_prop".png}}
        }%
        \subcaptionbox{\gls{ROC} curves}{%
            \label{fig:simulation_May15ROC:ROC}
            \includegraphics[width=0.475\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/RUN_20150507y/SL_Model_May15mvnorm_dir_1K_ROC2_samplesx_prop".png}}
        }%
        
        \subcaptionbox{Statistical performance (\acrshort{AUROC})}{%
           \label{fig:simulation_May15ROC:AUC}
           \includegraphics[width=0.65\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/RUN_20150507y/SL_Model_May15mvnorm_dir_1K_AUC_samplesx_prop".png}}
        }%
    \end{center}
   \caption[Performance with query correlations]{\small \textbf{Performance with query correlations.} Simulation of \glspl{synthetic lethal} was performed by sampling from a multivariate normal distribution (with correlation structure including correlated genes with non-synthetic lethal and query genes). Performance of \gls{SLIPT} declined for more \gls{synthetic lethal} partners and is mitigated by increased sample sizes (darker colours) but the sensitivity remains higher for a greater number of true positives with corresponding improvements in \gls{ROC} curves.}
%}
\label{fig:simulation_May15ROC2}
%\end{mdframed}
\end{figure*}

While the \gls{synthetic lethal} genes were detected in simple simulations (as shown in Figure~\ref{fig:simulation_May4SLreps}), \gls{ROC} analysis was performed to determine whether they were robustly detectable and to make further comparisons. These results (as shown in Figure~\ref{fig:simulation_Apr15ROC2}) were similar to simulations without correlation structure. As a binary classifier, \gls{SLIPT} had low sensitivity for higher numbers of \gls{synthetic lethal} partners to detect and high specificity with the vast majority of non-synthetic lethal genes (for 20,000 genes). This was reflected in a similar reduction in statistical performance for higher numbers of \gls{synthetic lethal} partners and a higher performance with higher sample size. Overall, the statistical performance was no different to simulations without correlation structure (as shown in Figure~\ref{fig:mvtnorm_cor_compare}).

\gls{SLIPT} was robust across correlation structures and is applicable to gene expression data, with \glslink{graph}{pathway} structures and correlations. These correlation structures were not intended to model specific biological pathways or represent them but showed potential impact of correlation structure on the performance of \gls{SLIPT} using highly correlated ($r = 0.8$) gene blocks. More complex correlation structures, such as genes positively correlated with the query gene and derived from pathway \glslink{graph}{graph} structures (as described in~\ref{methods:graphsim}) were examined further in Sections~\ref{chapt2:simulation_mvtnorm_query_cor} and~\ref{chapt5:graphsim_performance} respectively.

In particular, genes correlated with true \gls{synthetic lethal} genes had little impact on the performance of \gls{SLIPT} detection: \gls{synthetic lethal} genes were as distinguishable from correlated \textcolor{black}{genes as} they are from true negative genes. Genes correlated with \gls{synthetic lethal} partners did not interfere with the detection of true \gls{synthetic lethal} genes, although they were often ranked next below them and may support \gls{synthetic lethal} pathways by having related gene functions.

%\FloatBarrier

\subsubsection{Specificity with Query-Correlated Pathways}  \label{chapt2:simulation_mvtnorm_query_cor}

Correlation structures were also considered for non-synthetic lethal genes that were (positively) correlated genes with the query gene. Specifically, five highly correlated ($r = 0.8$) with the query gene were added (as described in Section~\ref{chapt2:simulation_mvtnorm_cor}). These simulations had similar performance (as shown in Figure~\ref{fig:simulation_May15ROC2}) to those without these correlations with a higher specificity and a lower false positive rate (shown in Figure~\ref{fig:simulation_May15ROC:Perf}).


%both have 5x5 cor genes (dir tested for May)
\iffalse
\begin{figure*}[!htbp]
%\begin{mdframed}
  \begin{center}
  \resizebox{0.5 \textwidth}{!}{
  \fbox{
    \includegraphics{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/RUN_20150507y/SL_Model_Test_Graph_10K_Graph1_ROC_Compare_Mvtn(Apry)_v_Cor(Mayxy)_Full"}}
   }
   }
   \end{center}
   \caption[Comparison of performance for query correlations]{\small \textbf{Comparison of performance for query correlations.} Multivariate simulation of \glspl{synthetic lethal} with correlation structure (in colour) clearly has lower performance than simulation including query correlations (in greyscale) across parameters. The query correlation simulation shows predictive potential for \gls{SLIPT} even with many underlying \gls{synthetic lethal} partners as positively correlated genes are distinguished robustly.}
\label{fig:mvtnorm_query_cor_compare}
%\end{mdframed}
\end{figure*}
\fi

%\FloatBarrier

%\subsubsection{Importance of Directional Testing}
The directional criteria of the \gls{SLIPT} procedure was important in this case, enhancing its performance, particularly in distinguishing positively correlated non-synthetic lethal genes. The multivariate normal simulations were performed, with 20,000 genes, including all of the correlation structures discussed (with synthetic lethal, non-synthetic lethal, and query correlated genes). These simulations were compared for the direction \gls{SLIPT} and  the $\chi^2$ testing. There was a considerably higher statistical performance with \gls{SLIPT}, particularly increased sensitivity and lower false positive rate (as shown in Figure~\ref{fig:mvtnorm_dir_compare}).

\begin{figure*}[!tb]
%\begin{mdframed}
%  \resizebox{\textwidth}{!}{
         \begin{center}
%
        \subcaptionbox{$\chi^2$ testing without direction}{%
            \label{fig:mvtnorm_dir_compare:perf_null}
            %\fbox{
            \includegraphics[width=0.475\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/RUN_20150507y/SL_Model_May15mvnorm_1K_ROC1_samplesx_prop".png}}
            %}
        }%
        \subcaptionbox{\gls{SLIPT} with directional criteria}{%
           \label{fig:mvtnorm_dir_compare:perf_dir}
            %\fbox{
           \includegraphics[width=0.475\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/RUN_20150507y/SL_Model_May15mvnorm_dir_1K_ROC1_samplesx_prop".png}}
           %}
        }%

    \end{center}
   \caption[Statistical evaluation of directional criteria]{\small \textbf{Statistical evaluation of directional criteria.} A simulated multivariate normal dataset of 20,000 genes with correlation structures was tested by \gls{SLIPT} with the directional condition and the $\chi^2$ test. \gls{SLIPT} exhibited a consistently higher sensitivity and lower false positive rate.}
%}
\label{fig:mvtnorm_dir_compare}
%\end{mdframed}
\end{figure*}

\begin{figure*}[!htbp]
%\begin{mdframed}
%  \resizebox{\textwidth}{!}{
         \begin{center}
%
        \subcaptionbox{$\chi^2$ testing without direction}{%
            \label{fig:mvtnorm_dir_compare:ROC_null}
            %\fbox{
            \includegraphics[width=0.475\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/RUN_20150507y/SL_Model_May15mvnorm_1K_ROC2_samplesx_prop".png}}
            %}
        }%
        \subcaptionbox{\gls{SLIPT} with directional criteria}{%
           \label{fig:mvtnorm_dir_compare:ROC_dir}
            %\fbox{
           \includegraphics[width=0.475\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/RUN_20150507y/SL_Model_May15mvnorm_dir_1K_ROC2_samplesx_prop".png}}
           %}
        }%


        \subcaptionbox{Statistical performance (\acrshort{AUROC}) }{%
           \label{fig:mvtnorm_dir_compare:AUC}
           \includegraphics[width=1\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/RUN_20150507y/SL_Model_Test_Graph_10K_Graph1_ROC_Compare_Mvtn(Mayy)_v_Cor(Mayxy)_Full".png}}
        }%
%
    \end{center}
   \caption[Performance of directional criteria]{\small \textbf{Performance with directional criteria.} A simulated multivariate normal dataset of 20,000 genes with correlation structures was tested by \gls{SLIPT} and $\chi^2$ test. \gls{SLIPT} had higher performance across simulation parameters, clearly differing from random (grey diagonal) in \gls{ROC} curves up to 100 SL genes (b). The performance (c) of \gls{SLIPT} (in greyscale) was consistently higher than the $\chi^2$ test (in color).}
%}
\label{fig:mvtnorm_dir_compare2}
%\end{mdframed}
\end{figure*}

These results show that the performance of \gls{SLIPT} is appropriate for the analysis of \glslink{gene expression}{expression} datasets, where positively correlated genes commonly occur, with the directional condition robustly improving the performance of \gls{SLIPT} across simulation parameters (compared to the $\chi^2$ test). Without assuming the underlying number of \gls{synthetic lethal} genes, \gls{SLIPT} will performed than the $\chi^2$ test alone, irrespective of the significance threshold as shown by \gls{ROC} analysis (as shown in Figure~\ref{fig:mvtnorm_dir_compare2}). The directional \gls{SLIPT} methodology outperformed the $\chi^2$ test at detecting \gls{synthetic lethal} partners with even up to 100 \gls{synthetic lethal} genes.

Together these simulation results support the application of the \gls{SLIPT} methodology as it has been performed throughout Chapters~\ref{chap:SLIPT} and~\ref{chap:Pathways}. The methodology and simulation procedure were explored further in Chapter~\ref{chap:simulation}, with comparison to other \gls{synthetic lethal} detection approaches and the inclusion of \glslink{graph}{graph} structures.

\FloatBarrier

\section{Graph Structure Methods}
Graph structures have been used in several ways in this project, including novel approaches to analysis and simulations. Procedures were developed for statistical and network analysis of gene states in \glslink{graph}{pathway} structures. Specifically, the relationships between \gls{siRNA} and \gls{SLIPT} genes were tested within biological pathways in Chapter~\ref{chap:Pathways}. These \glslink{graph}{graph} structures were also used in Chapter~\ref{chap:simulation} to derive correlation structure between simulated \gls{gene expression} profiles to represent biological pathways.


\subsection{Upstream and Downstream Gene Detection} \label{methods:pathway_str} 
Comparison of experimental and computational candidate \gls{synthetic lethal} partner genes within \glslink{graph}{pathway} structures was performed to determine whether these sets of genes were related by \glslink{graph}{pathway} structure. Considering the differences in how these candidates were generated, it was unsurprising that they did not detect some identical genes within the candidate biological pathways. However, they could still be related by being upstream or downstream of each other. 

Using the Reactome version 52 \citep{Reactome}, as described in Section~\ref{methods:graph_data}, genes detected by each \gls{synthetic lethal} discovery approach were mapped to the \glslink{graph}{graph} structure for each candidate pathway identified in Chapter~\ref{chap:SLIPT} (with graphs defined as described in Section~\ref{methods:subgraphs}). To test whether \gls{siRNA} candidate genes were upstream of \gls{SLIPT} candidate genes, \glspl{shortest path} were traced between each pair of these genes in a directed network. The paths where the \gls{siRNA} candidate was upstream (``up'') and downstream (``down'') of a \gls{SLIPT} candidate were scored.  This procedure yielded the total number of \glspl{shortest path} which indicated that \gls{siRNA} genes were upstream or downstream of the \gls{SLIPT} genes and measured the difference between these to determine if there was an imbalance in a particular direction. While this difference was indicative of the number of paths between the gene candidate groups in either direction, it was not sufficient to statistically verify structure or relationships between \gls{siRNA} and \gls{SLIPT} genes. It was \textcolor{black}{then} combined with a permutation resampling procedure (as described in Section~\ref{methods:network_permutation}) to test for directional relationships.

Initially, this procedure excluded genes that were detected by both approaches since they would count in both directions. Upon further consideration, these genes were restored to accounted for \textcolor{black}{that} they may contribute unequally to each gene set if there are unequal numbers of genes above or below them in the \glslink{graph}{pathway} structure.

\subsubsection{Permutation Analysis for Statistical Significance} \label{methods:network_permutation}
A permutation procedure was developed to randomly assign members of the pathway to \gls{siRNA} and/or \gls{SLIPT} groups, with the same number of each candidate partner gene set as observed in the pathway. These permuted genes were measured for \glslink{graph}{pathway} structure between the permuted gene groups as performed for the observed candidates (as performed in Section~\ref{methods:pathway_str}). A distribution of \glslink{graph}{pathway} structure relationships expected by chance was generated by permuting iteratively over these pathways. The resulting null distribution was compared to the observed counts of relationships (in either direction). This procedure yielded a permutation p-value as the proportion of permutations \textcolor{black}{which} had a value greater than the observed value. The null hypothesis was that there was no relationship between these gene groups compared to genes that had been selected at random. Thus both the alternate hypotheses that the \gls{siRNA} genes were either upstream of the \gls{SLIPT} genes or that they are downstream of them were testable.

The permutation procedure does not assume the underlying distribution of the data under the null hypothesis and accounts for the total number of \glslink{vertex}{nodes}, \glspl{edge}, \gls{siRNA}, and \gls{SLIPT} genes in each \gls{graph}{pathway} structure. The number of genes detected by both \gls{siRNA} and \gls{SLIPT} was not accounted for under the initial \gls{shortest path} counts procedure that excluded them. Once they were included, it was ensured that the number of intersecting genes was equal to the number observed to test for \glslink{graph}{pathway} structure without changing the intersection size, the subject of prior analyses.

\iffalse
\subsubsection{Hierarchy Based on Biological Context} \label{methods:pathway_rank}
An alternative approach to \glslink{graph}{pathway} structure was based on the biological context, given that genes at the upstream and downstream ends of a pathway perform different functions, such as a kinase signalling cascade receiving signals from external stimuli and passing these on to ribosomes or the nucleus. Genes were assigned to a hierarchy to determine if genes of either candidate group disproportionately performed upstream or downstream functions.

A network-based approach was used to generate the pathway hierarchy of genes in a computationally rational way when applied to different biological pathways with a directed \glslink{graph}{graph} structure, $G$ (without loops). The diameter of the network (i.e., the length of the longest possible \gls{shortest path} between the most distant genes) was used to identify a gene ($z$) at the downstream end of the pathway (at the end of a diameter spanning \gls{shortest path}), which was assigned a hierarchy of: $$hierarchy(z) = 1 + diameter(G).$$ Having identified the downstream end of the pathway, genes upstream (e.g., gene $i$) of this were assigned a hierarchy by the length of their \gls{shortest path} ($d$) to this gene $z$. $$hierarchy(i) = hierarchy(z) - d_{iz}.$$ The remaining unassigned genes (e.g., gene $j$) gained the hierarchy of the length of the \gls{shortest path} downstream from the nearest assigned gene if possible $$hierarchy(j) = hierarchy(i) + d_{ij}.$$ This process could be performed iteratively to fill in pathway hierarchy but it was not necessary to perform further iterations for the candidate \gls{synthetic lethal} pathways investigated which exhibited strong directional structure and the \gls{small world} property (i.e., had a low diameter). Using this procedure, genes in a pathway \glslink{graph}{graph} structure were assigned to an integer valued hierarchy upstream to downstream by this procedure: $$hierarchy \in \{1, 2, 3, ..., 1 + diameter(G)\}$$
%new paragraph?
This hierarchy of pathway directionality (e.g.,  that shown in Figure~\ref{fig:SL_Pathway_PI3K_Ranking}) was used for comparison with measures of the number of \gls{synthetic lethal} partners detected by either approach.
\fi

\subsection{Simulating Gene Expression from Graph Structures} \label{methods:graphsim}
The simulation procedure was refined to generate \glslink{gene expression}{expression} data with correlation structure from a known \glslink{graph}{graph} structure. %rather than correlated blocks for $\Sigma$ in \texttt{mvtnorm} \citet{mvtnorm}
This enabled modelling of \gls{synthetic lethal} partners within a biological pathway and the investigation of the impact of \glslink{graph}{pathway} structure on \gls{synthetic lethal} prediction. Firstly, a simulated pathway was constructed as a \glslink{graph}{graph} structure, with the \texttt{igraph} R package \citet{igraph}, with the state of the \glspl{edge} (i.e, whether they activate or inhibit downstream pathway members). This simulation procedure was intended for biological pathway members with correlated \gls{gene expression} (higher than the background of genes in other pathways) but it may also be applicable to modelling protein levels (e.g, in a kinase regulation cascade) or substrates and products (e.g., in a metabolic pathway).% as it is to gene regulation as it has been applied here.

The \glslink{graph}{graph} structure was constructed \textcolor{black}{to be used to simulate data by} sampling a multivariate normal distribution using the \texttt{mvtnorm} R package \citep{Genz2009, mvtnorm}. Throughout this Section, the simulation procedure will be demonstrated with the relatively simple constructed \glslink{graph}{graph} structure shown in Figure~\ref{fig:simple_graph}. This \glslink{graph}{graph} structure visualisation was specifically developed for (directed) iGraph objects in R and has been released in the \texttt{plot.igraph} package and \texttt{igraph.extensions} library (in Table~\ref{tab:computers_r_packages_dev} and Section~\ref{methods:igraph_extensions}). The \texttt{plot\_directed} function enabled customisation of plot parameters for each \glslink{vertex}{node} or \glslink{edge}{edge} and mixed (directed) \glslink{edge}{edge} types for indicating activation or inhibition. These inhibition \glslink{edge}{links} (which occur frequently in biological pathways) were demonstrated in Figure~\ref{fig:simple_graph:second}.

\begin{figure*}[!htb]
%\begin{mdframed}
%  \resizebox{\textwidth}{!}{
         \begin{center}
%
        \subcaptionbox{Activating \glslink{graph}{pathway} structure }{%
            \fbox{
            \includegraphics[width=0.45\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/simple_graph".png}}
            }
        }%
        \subcaptionbox{Pathway structure with inhibitions \label{fig:simple_graph:second}}{%
            \fbox{
           \includegraphics[width=0.45\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/simple_graph_inhibiting".png}}
           }
        }%
%
    \end{center}
   \caption[Simulated \glslink{graph}{graph} structures]{\small \textbf{\textbf{Simulated \glslink{graph}{graph} structures.}} A constructed \glslink{graph}{graph} structure used as an example to demonstrate the simulation procedure. Activating \glslink{edge}{links} are denoted by blue arrows and inhibiting \glslink{edge}{links} by red \glspl{edge}. Nodes have been coloured accordingly to the simulation results presented in Section~\ref{chapt5:graphsim_str} and Appendix~\ref{appendix:graphsim_str}.}
%}
\label{fig:simple_graph}
%\end{mdframed}
\end{figure*}

\begin{figure*}[!hp]
%\begin{mdframed}
%  \resizebox{\textwidth}{!}{
         \begin{center}
%
        \subcaptionbox{Activating \glslink{graph}{pathway} structure \label{fig:simulation_activating:first}}{%
            \includegraphics[width=0.35\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/simple_graph".png}}
        }%
        \subcaptionbox{Distance matrix \label{fig:simulation_activating:second}}{%
            \includegraphics[width=0.35\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/dist_mat".png}}
        }%
        
        \subcaptionbox{$\Sigma$ (expected correlation) \label{fig:simulation_activating:third}}{%
           \includegraphics[width=0.35\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/sigma_mat".png}}
        }%
	\subcaptionbox{Simulated correlation\label{fig:simulation_activating:fifth}}{%
           \includegraphics[width=0.35\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_cor_mat".png}}
        }%
        	
	\subcaptionbox{Simulated \glslink{gene expression}{expression} data  \label{fig:simulation_activating:fourth}}{%
            \includegraphics[width=0.35\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_mat".png}}
        }%
        \subcaptionbox{Simulated gene function calls \label{fig:simulation_activating:sixth}}{%
           \includegraphics[width=0.35\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_disc_mat".png}}
        }%
    \end{center}
   \caption[Simulating \glslink{gene expression}{expression} from a \glslink{graph}{graph} structure]{\small \textbf{\textbf{Simulating \glslink{gene expression}{expression} from a \glslink{graph}{graph} structure.}} An example \glslink{graph}{graph} structure that was used to derive a correlation structure from the relative distances between \glslink{vertex}{nodes} and simulate continuous \gls{gene expression} (coloured blue--red) with sampling from the multivariate normal distribution.}
%}
\label{fig:simulation_activating}
%\end{mdframed}
\end{figure*}

\begin{figure*}[!hp]
%\begin{mdframed}
%  \resizebox{\textwidth}{!}{
         \begin{center}
%
        \subcaptionbox{Inhibiting pathway structure\label{fig:simulation_inhibiting:first}}{%
            \includegraphics[width=0.35\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/simple_graph_inhibiting".png}}
        }%
        \subcaptionbox{Distance matrix \label{fig:simulation_inhibiting:second}}{%
            \includegraphics[width=0.35\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/dist_mat".png}}
        }%
        
        \subcaptionbox{$\Sigma$ (expected correlation) \label{fig:simulation_inhibiting:third}}{%
           \includegraphics[width=0.35\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/sigma_mat_inhibiting".png}}
        }%
	\subcaptionbox{Simulated correlation\label{fig:simulation_inhibiting:fifth}}{%
           \includegraphics[width=0.35\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_cor_mat_inhibiting".png}}
        }%
        
        \subcaptionbox{Simulated \glslink{gene expression}{expression} data \label{fig:simulation_inhibiting:fourth}}{%
            
            \includegraphics[width=0.35\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_mat_inhibiting".png}}
        }%
        \subcaptionbox{Simulated gene function calls \label{fig:simulation_inhibiting:sixth}}{%
           \includegraphics[width=0.35\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_inhib_disc_mat".png}}
        }%
    \end{center}
   \caption[Simulating \glslink{gene expression}{expression} from \glslink{graph}{graph} structure with inhibitions]{\small \textbf{\textbf{Simulating \glslink{gene expression}{expression} from \glslink{graph}{graph} structure with inhibitions.}} An example \glslink{graph}{graph} structure that was used to derive a correlation structure from the relative distances between \glslink{vertex}{nodes} and simulate continuous \gls{gene expression} (coloured blue--red) with sampling from the multivariate normal distribution.}
%}
\label{fig:simulation_inhibiting}
%\end{mdframed}
\end{figure*}

The simulation procedure was designed to use such \glslink{graph}{graph} structures to inform development of a ``Sigma'' variance-covariance matrix ($\Sigma$) for sampling from a multivariate normal distribution (using the \texttt{mvtnorm} R package). \textcolor{black}{This procedure uses} a \glslink{graph}{graph} structure (or adjacency matrix \textcolor{black}{$A_{ij}$}), such as Figure~\ref{fig:simulation_activating:first}. \textcolor{black}{
%
\[
A_{ij} = 
\begin{dcases*}
   1                         & if genes $i$ and $j$ are adjacent \\
   0                         & otherwise
\end{dcases*}
\]
%
A} relation matrix \textcolor{black}{$R_{ij}$} was calculated based on distance such that nearer \glslink{vertex}{nodes} are given higher weight than farther \glslink{vertex}{nodes}. Throughout this thesis, a geometrically decreasing (relative) distance weighting was used \textcolor{black}{
%
\[
R_{ij} = 
\begin{dcases*}
   1                            & if genes $i$ and $j$ are adjacent \\
   (\sfrac{1}{2})^{d_{ij}}   & if a path can be found  between genes $i$ and $j$ \\
   0                            & if no path exists between genes $i$ and $j$ 
\end{dcases*}
\]
%
where $d_{ij}$ is the shortest path length between genes $i$ and $j$ in graph $G$. Thus} with each more distant \glslink{vertex}{node} being related by $\sfrac{1}{2}$ compared to the next nearest, as shown in Figure~\ref{fig:simulation_activating:second}.
An arithmetically decreasing (absolute) distance weighting is also supported in the \texttt{graphsim} R package release of this procedure. \textcolor{black}{
%
\[
R_{ij} = 
\begin{dcases*}
   1                            & if genes $i$ and $j$ are adjacent \\
   1-\frac{d_{ij}}{diam(G)}     & if a path can be found  between genes $i$ and $j$ \\
   0                            & if no path exists between genes $i$ and $j$ 
\end{dcases*}
\]
%
}


A $\Sigma$ matrix can be derived from this distance weighting matrix \textcolor{black}{
%
\[
\Sigma_{ij} = 
\begin{dcases*}
   1                            & if $i=j$ so a gene is correlated with itself \\
   correlation$\times R_{ij}$  & otherwise
\end{dcases*}
\]
%
}
%
creating a matrix (with a diagonal of $1$) where each \glslink{vertex}{node} has a variance and standard deviation of 1. Thus covariances between adjacent \glslink{vertex}{nodes} were assigned by a correlation parameter and the remaining matrix based on weighting these correlations by the geometrically weighted distance matrix (or the nearest ``positive definite'' matrix for $\Sigma$ weighted for negatively correlated inhibitions). \textcolor{black}{Computing the nearest positive definite matrix was necessary to ensure that the variance-covariance matrix could be inverted when used as a parameter in multivariate normal simulations. Matrices that could not be inverted occurred rarely with the graph structures used but was \textcolor{black}{included} for robustly applicable functions in the R software release. The nearest positive definite matrix was computed using the \texttt{nearPD} function of the \texttt{Matrix} R package \citep{Matrix} to perform Higham's algorithm \citep{Higham2002}, a procedure designed for variance-covariance matrices.}

Throughout this thesis, the correlation parameter was $0.8$, unless otherwise specified (as used for the example in Figure~\ref{fig:simulation_activating:third}). This $\Sigma$ matrix was then used to sample from a multivariate normal distribution such that each gene had a mean of $0$, standard deviation $1$, and covariance within the range $[0,1]$ such that they are correlations. This procedure generated a simulated (continuous normally distributed) \glslink{gene expression}{expression} profile for each \glslink{vertex}{node} (as shown in Figure~\ref{fig:simulation_activating:fourth}) with corresponding correlation structure (Figure~\ref{fig:simulation_activating:fifth}). The simulated correlation structure closely resembled the expected correlation structure (Sigma in~\ref{fig:simulation_activating:third}) even for the relatively modest sample size ($N=100$) illustrated in~\ref{fig:simulation_activating}. Once a simulated \gls{gene expression} dataset has been generated (as in Figure~\ref{fig:simulation_activating:fourth}), then a discrete matrix of gene function was constructed with a functional threshold quantile to simulate functional relationships of \glspl{synthetic lethal} (as shown in Figure~\ref{fig:SL_Model_Expression}). Throughout this thesis, this threshold is the 0.3 quantile (as discussed in Section~\ref{methods:SL_Model}) which generates functional discrete matrices such as those used for \gls{synthetic lethal} simulation in Section~\ref{methods:simulating_SL} (as shown Figure~\ref{fig:simulation_activating:sixth}).

The simulation procedure (depicted in Figure~\ref{fig:simulation_activating}) can be used for pathways containing inhibition \glslink{edge}{links} (as shown in Figure~\ref{fig:simulation_inhibiting}) with several refinements. With the inhibition \glslink{edge}{links} (as shown in Figure~\ref{fig:simulation_inhibiting:first}), distances were calculated in the same manner as before (Figure~\ref{fig:simulation_inhibiting:second}) with inhibitions accounted for by iteratively multiplying downstream \glslink{vertex}{nodes} by $-1$ to form blocks of negative correlations (as shown in Figures~\ref{fig:simulation_inhibiting:third} and~\ref{fig:simulation_inhibiting:fifth}). A multivariate normal distribution with these negative correlations can be sampled to generate simulated data (as shown in Figures~\ref{fig:simulation_inhibiting:fourth} and~\ref{fig:simulation_inhibiting:sixth}).  

These simulated datasets could then be used for simulating \gls{synthetic lethal} partners of a query gene within a graph network. The query gene was assumed to be separate from the graph network pathway and was added to the dataset using the procedure in Section~\ref{methods:simulating_SL}. Thus I can simulate known \gls{synthetic lethal} partner genes within a \gls{synthetic lethal} partner \glslink{graph}{pathway} structure.

\iffalse

	\includegraphics{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/dist_mat".png}}
	\includegraphics{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/sigma_mat".png}}
		\includegraphics{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_mat".png}}
			\includegraphics{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_cor_mat".png}}
		\includegraphics{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_disc_mat".png}}
		
	\includegraphics{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/state_matrix_inhibiting".png}}
	\includegraphics{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/dist_mat".png}}
		\includegraphics{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/sigma_mat_inhibiting".png}}
		\includegraphics{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_inhib_mat".png}}
			\includegraphics{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_inhib_cor_mat".png}}
		\includegraphics{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/expr_inhib_disc_mat".png}}
	
\fi	

\FloatBarrier

\section{Customised Functions and Packages Developed} \label{methods:r_packages}

Various R packages \citep{R_core} have been developed throughout this thesis using \texttt{devtools} \citep{devtools} and \texttt{roxygen} \citep{roxygen} to enable reproducibility of customised analysis and visualisation. Many of these have been documented, demonstrated in vignettes, and released on GitHub (\url{https://github.com/TomKellyGenetics})
to enable the research community to utilise them in their own analysis.
These are summarised in Table~\ref{tab:computers_r_packages_dev}, with the corresponding urls for their GitHub repository which contains instructions for installation with the \texttt{devtools} R package \citep{devtools} and links the relevant vignette(s).

\subsection{Synthetic Lethal Interaction Prediction Tool}
The statistical methodology for detection of \glspl{synthetic lethal} in \gls{gene expression} data (\gls{SLIPT}) is one of the main novel procedures developed in this thesis, as described in Section~\ref{methods:SLIPT}. The \texttt{slipt} R package has been prepared for release %(\url{https://github.com/TomKellyGenetics/slipt})
to accompany a publication demonstrating the applications of the methodology for identifying candidate interacting genes and pathways with \textit{CDH1} in breast cancer \citep{TCGA2012}.

\gls{SLIPT} can be \textcolor{black}{used for} analysis of any effectively continuous measure of gene activity (e.g., \gls{microarray}, \gls{RNA-Seq}, protein abundance, or pathway \glspl{metagene}). Executing \texttt{slipt} is straightforward: the \texttt{prep\_data\_for\_SL} function scores samples as ``low'', ``medium'', or ``high'' for each gene, then the \texttt{detect\_SL} function tests a given query gene against all potential partners by performing the chi-squared test and directional conditions. This function returns a table summarising the observed and expected sample numbers used for the directional criteria, the $\chi^2$ values, and corresponding p-values including adjusting for multiple comparisons. The \texttt{count\_of\_SL} and \texttt{table\_of\_SL} functions serve to facilitate summary and extraction of the positive \gls{SLIPT} hits, respectively, from the table of predictions of \gls{synthetic lethal} partners.

The \gls{SLIPT} methodology in this package release was used in later analyses rather than the corresponding source R code, including use on remote machines and upon simulated data. In particular, the functions in the package facilitate alterations to parameters, such as the proportion of samples called as exhibiting low or high gene activity (as shown in Section~\ref{chapt5:compare_chisq}). This release supports reproducible research and enables wider use of \gls{SLIPT} in future investigations into other disease genes.

\subsection{Data Visualisation} \label{methods:r_packages_vis}
Customised data visualisations in R \citep{R_core} were developed to present data throughout this thesis. %The \texttt{vioplotx} and \texttt{heatmap.2x} packages are enhancements of the \texttt{vioplot} package \citep{vioplot} and \texttt{heatmap.2} provided by the \texttt{gplots} package \citep{gplots}. 
%
The \texttt{vioplotx} package provides an alternative visualisation (of continuous variables against categories) to the more familiar boxplot, showing variability of the data by the width of the plots. As demonstrated in Figure~\ref{fig:vioplot}, this version enables separate plotting parameters for each violin with vector inputs for colour, shape, and size of various elements of the median point, central boxplot, borders, and fill colour for the violin. Scaling violin width to adjust violin area and splitting data by a second categorical variable is also enabled. This function is intended to be backwards compatible with the \texttt{vioplot} package \citep{vioplot}  (applying scalar inputs across all violins) and \texttt{boxplot} (by enabling formula inputs as an S3 method). Each of these features has been demonstrated with examples in respective vignettes on the package \href{https://github.com/TomKellyGenetics/vioplotx}{GitHub repository} (\url{https://github.com/TomKellyGenetics/vioplotx}).

The \texttt{heatmap.2x} function provides extensions for annotation colour bars for both the rows and columns (as shown in Figure~\ref{fig:heatmap.2x}). Multiple bars are enabled on both axes with matrix inputs (rather than single vector for \texttt{heatmap.2} \citep{gplots}) which facilitates additional plotting of gene and sample characteristics for comparison with correlation matrices, \glslink{gene expression}{expression} profiles, or pathway \glspl{metagene}. The annotation bar inputs correspond to their orientation on the plot, each colour bar is provided as a column for the row annotation on the left of the heatmap and as a row for the column annotation on top of the heatmap. Row and column annotation bars are labelled with the column or row names respectively. Additional parameters enable resizing of these annotation bar labels and control of reordering columns for when samples have been ordered in advance (e.g., ranked by a \gls{metagene} or split into groups clustered separately).  These features were used through this thesis and have been provided in a package \href{https://github.com/TomKellyGenetics/heatmap.2x}{GitHub repository} (\url{https://github.com/TomKellyGenetics/heatmap.2x}).

\begin{figure*}[!thb]
%\begin{mdframed}
%  \resizebox{\textwidth}{!}{
         \begin{center}
%
        \subcaptionbox{Customised violin plot}{%
            \label{fig:vioplot:first}
            \includegraphics[width=0.75\textwidth]{{"vioplot1".png}}
        }
        
        \subcaptionbox{Split violin plot}{%
            \label{fig:vioplot:second}
            \includegraphics[width=0.75\textwidth]{{"vioplot2".png}}
        }%
        \end{center}
   \caption[Demonstration of violin plots with custom features]{\small \textbf{Demonstration of violin plots with custom features.} An example of the \texttt{iris} dataset is plotted to show the custom features of the \texttt{vioplotx} package including (a) individual colour, shape and size parameters of each violin, scaling violin widths by area, and (b) splitting violins to compare subsets of data.}
%}
\label{fig:vioplot}
%\end{mdframed}
\end{figure*}

\begin{figure*}[!thb]
%\begin{mdframed}
%  \resizebox{\textwidth}{!}{
         \begin{center}
            \includegraphics[width=1 \textwidth]{{"heatmap2x".png}}
        \end{center}
   \caption[Demonstration of annotated heatmap]{\small \textbf{Demonstration of annotated heatmap}. The example heatmap depicts the additional row and column annotation bars enabled by \texttt{heatmap.2x}, extending the features of \texttt{gplots} with backwards compatible inputs.}
%}
\label{fig:heatmap.2x} 
%\end{mdframed}
\end{figure*}

%\FloatBarrier

\subsection{Extensions to the iGraph Package} \label{methods:igraph_extensions}
The following features were developed during this thesis using ``iGraph'' data objects, building upon the \texttt{igraph} package \citep{igraph}. These have been released as separate packages for each respective procedure and can be installed together as a collection of extensions to the \texttt{igraph} package (\url{https://github.com/TomKellyGenetics/igraph.extensions}).

\subsubsection{Sampling Simulated Data from Graph Structures}
The \texttt{graphsim} package implements the procedure for simulating \gls{gene expression} from \glslink{graph}{graph} structures (as described in Section~\ref{methods:graphsim}). By default, this derives a matrix with a geometrically decreasing weighting by distance (by \glspl{shortest path}) between each pair of \textcolor{black}{\glslink{vertex}{nodes}}. An absolute decreasing weighting is also available with the option \textcolor{black}{to} derive correlation structures from adjacency matrices or the number of \glslink{edge}{links} common partners (i.e., size of the shared ``neighbourhood'' \citep{Hell1976}) between each pair of \glslink{vertex}{nodes}. Functions to compute these are called directly by passing parameters to them when running the \texttt{generate\_expression} or \texttt{make\_sigma\_mat} commands. This package enables simulating \glslink{gene expression}{expression} data directly from a \glslink{graph}{graph} structure (with the intermediate steps automated) or generating $\Sigma$ parameters for \texttt{mvtnorm} from \glslink{graph}{graph} structures or matrices derived from them. These functions support assignment of activating or inhibiting relationships to each \glslink{edge}{edge} (with a \texttt{state} parameter).

\subsubsection{Plotting Directed Graph Structures}
The \texttt{plot.igraph} package provides the \texttt{plot\_directed} function, specifically developed for directed \glslink{graph}{graph} structures, to plot activating or inhibiting for each \glslink{edge}{edge} (as described in Section~\ref{methods:graphsim}). As shown in Figure~\ref{fig:simple_graph2}, this function supports separate plotting parameters for each \glslink{vertex}{node}, \glslink{vertex}{node} label, and \glslink{edge}{edge}. This includes colours of \glslink{vertex}{node} fill, border, label text, and \glspl{edge} and size of \glslink{vertex}{nodes}, \glslink{edge}{edge} widths, arrowhead lengths, and font size of labels. The  \texttt{state} parameter for assigning activating or inhibiting to each \glslink{edge}{edge} determines whether \glspl{edge} were depicted with 30\textdegree\ or 90\textdegree\ arrowheads. Colours are assigned separately so they may be customised. Vectorised parameters are applied across each \glslink{vertex}{node} or \glslink{edge}{edge}, whereas scalar parameters apply the same plotting parameters across them. The default layout function is \texttt{layout.fruchterman.reingold} but any layout function supported by \texttt{plot} function in \texttt{igraph} \citep{igraph} was compatible, such as \texttt{layout.kamada.kawai} used to implement the Kamada--Kawai algorithm \citep{Kamada1989} for graph plots throughout this thesis.
 

 \begin{figure*}[!htb]
%\begin{mdframed}
%  \resizebox{\textwidth}{!}{
         \begin{center}
         \fbox{
          \includegraphics[width=0.5\textwidth]{{"/home/tomkelly/Documents/PhD Otago Uni/SL_Model/graph_sim_method/simple_graph_inhibiting".png}}
	  }
	  \end{center}
   \caption[Simulating \glslink{graph}{graph} structures]{\small \textbf{\textbf{Simulating \glslink{graph}{graph} structures.}} An example \glslink{graph}{graph} structure which has been used throughout demonstrating the simulation procedure from \glslink{graph}{graph} structures. Activating \glslink{edge}{links} are denoted by blue arrows and inhibiting \glslink{edge}{links} by red \glspl{edge}.}
%}
\label{fig:simple_graph2}
%\end{mdframed}
\end{figure*}

%\FloatBarrier
 
\subsubsection{Computing Information Centrality} 
The \glspl{shortest path} of a network were computed by the \texttt{igraph} package \citep{igraph} which can be extended to calculate the network efficiency but was not provided by the package itself (ss described in Section~\ref{methods:network_metrics}). The ``\gls{information centrality}'' of a \glslink{vertex}{vertex} is computed as the relative change in the network efficiency when the \glslink{vertex}{vertex} is removed. \Gls{information centrality} is calculated iteratively for each \glslink{vertex}{node} and the sum of \gls{information centrality} for each \glslink{vertex}{vertex} is the \gls{information centrality} for the network. These metrics were released in the \texttt{info.centrality} package (\url{https://github.com/TomKellyGenetics/info.centrality}).

\subsubsection{Testing Pathway Structure with Permutation Testing}
A network-based procedure \textcolor{black}{was developed to compare} \gls{siRNA} and \gls{SLIPT} candidate genes in a \glslink{graph}{pathway} structure. Such \glslink{graph}{pathway} structure relationships were tested by computing the number of \glspl{shortest path} between two different groups of \glslink{vertex}{nodes} in either direction within a graph . This pathway relationship metric was implemented in the \texttt{pathway.structure.permutation} package (\url{https://github.com/TomKellyGenetics/pathway.structure.permutation}) with permutation testing (as described in Sections~\ref{methods:pathway_str} and~\ref{methods:network_permutation}). 

\subsubsection{Metapackage to Install iGraph Functions}
These features may be installed together with \texttt{igraph.extensions}, which can be accessed from a \href{https://github.com/TomKellyGenetics/igraph.extensions}{GitHub repository} (\url{https://github.com/TomKellyGenetics/igraph.extensions}). This meta-package installs \texttt{igraph} \citep{igraph} and the packages described in Section~\ref{methods:igraph_extensions} including their dependencies for matrix operations and statistical procedures: \texttt{Matrix}, \texttt{matrixcalc}, and \texttt{mvtnorm} \citep{Matrix, matrixcalc, Genz2009, mvtnorm}.


%simulation 2015 committee meeting
\iffalse
\section{Developing a Synthetic Lethal detection methodology}

\subsection{Testing multivariate normal Simulation of \Glspl{synthetic lethal}}

I have developed a model of \glspl{synthetic lethal} in \gls{gene expression} data based on sampling a multivariate normal distribution.  This enables simulation of statistically testing for \gls{synthetic lethal} genes where the true and false positives are known, discovery of the expected test statistic distributions for different conditions, educated thresholds for public data analysis, and building a complex model with known correlation structure between genes.  Sampling a small number of genes from this model shows, in Figure 4, that \glspl{synthetic lethal} is detectable with in a simple model.

Figure 4.  Chi-Square (upper) and p-values (lower) distributions show that \gls{synthetic lethal} partners (red) are distinguishable from correlated (blue) and other genes (black) in an example simulation of sampling 1000 samples and 100 genes, from a multivariate normal distribution with 1 (left), 2 (centre), and 3 (right) \gls{synthetic lethal} partners respectively, showing that \gls{synthetic lethal} genes become more difficult to detect if there are more true partners.

Figure 5.  Chi-Square (upper), \gls{FDR} adjusted p-values (centre), and Holm adjusted p-values (lower) show that show that \gls{synthetic lethal} partners (red) are distinguishable from correlated (blue) and other genes (black) are distinguishable replicated across 1000 replicate simulated sampling of 1000 samples and 100 genes, from a multivariate normal distribution with 1 (left), 2 (centre), and 3 (right) \gls{synthetic lethal} partners respectively, showing \gls{synthetic lethal} genes become more difficult to detect in with more true partners but adjusting p-values may be too stringent an approach to this.

Having shown that the Chi-Square test is capable of detecting \glspl{synthetic lethal}, Figure 5 shows that detecting \glspl{synthetic lethal} in a simple case is largely robust and reproducible across many replicates with \gls{synthetic lethal} and correlated genes clearly having higher test statistic scores and lower adjusted p-values than the null distribution of non-synthetic lethal genes when there are only 1 or 2 \gls{synthetic lethal} partners.  While it is promising that correlated genes and \gls{synthetic lethal} partners could be distinguished from other genes in a simple case, there is also indication that true \gls{synthetic lethal} partners (candidates as robust drug targets) and their correlated genes (or pathways) could be distinguished by test statistic.

However, such clear evidence of \glspl{synthetic lethal} by co-loss under-representation is rarely detected in public data analyses, indicating cryptic additional \gls{synthetic lethal} genes compensating for the loss of both the query and putative \gls{synthetic lethal} partner.  Therefore higher-order \gls{synthetic lethal} is potentially very common, difficult to detect, and confounding attempts to identify \gls{synthetic lethal} pairs from \gls{gene expression} data.  In Figure 5, more than 3 \gls{synthetic lethal} partners will be difficult to identify directly with a Chi-Square test.  Although deeper understanding of the system could still enable use of the procedure to prioritise small numbers of candidate genes, estimate the number of underlying true \gls{synthetic lethal} partners, and identify the biological pathways interacting with a gene to focus complementary experimental approaches.

With higher number of true \gls{synthetic lethal} genes there is no clear threshold for Chi-Square values (or associated p-values) to detect \glspl{synthetic lethal} and choosing any threshold is a trade-off between sensitivity (ensuring all true positives are detected) and specificity (reducing the number of false positives detected).  Receiver operating characteristic (ROC) curves, as shown in Figures 6 and 7, summarise this trade-off to show the statistical performance of a test where the true \gls{synthetic lethal} genes are known in the simulated data.  Performance of a statistical test is measured as the area under the \gls{ROC} (AUROC) curves, as shown in Figures 8 and 9, to compare performance across simulations for different parameters such as type of model, correlation structure, the total number of genes, sample size and number of true \gls{synthetic lethal} genes.  A random predictor has an \gls{AUROC} of 0.5, whereas an ideal predictor has an \gls{AUROC} of 1.0, so intermediate values are expected.

\subsection{Receiver Operating Characteristic Curves}

Figure 6.   \gls{ROC} curves showing statistical performance (by area under the curve) of a \gls{synthetic lethal} simulation based on sampling a Binomial distribution, with 20,000 genes, averaged over 1000 replicates, sample size (1000, 2000, 5000, or 10,000) and number of \gls{synthetic lethal} genes (up to 100) varies by panel and colour showing better performance with fewer \gls{synthetic lethal} genes or higher sample size.    

Figure 7.   \gls{ROC} curves of a \gls{synthetic lethal} simulation based on sampling a multivariate normal distribution, with 20,000 genes, averaged over 1000 replicates, sample size and number of \gls{synthetic lethal} genes varies by panel and colour showing better performance than a Binomial model and similar performance with correlation structure (upper panes).

Figure 8.  Comparison of Binomial (red) and multivariate normal models with (blue) and without (green) correlation structure by simulation with 1000 samples, 20,000 genes, sample size varied by pane, and number of \gls{synthetic lethal} partners on the x axis where performance on the y axis is measured as the \gls{AUROC} showing better performance in the multivariate normal model than the Binomial model and similar performance in the multivariate normal model with correlation structure added for all simulation parameters.  There was better performance with fewer \gls{synthetic lethal} partners or higher sample size with both multivariate normal models.   

Figure 6 shows performance of an earlier model based on the Binomial distribution for gene function calls, based on similar a Normally distributed model of \gls{gene expression} which called gene function from an arbitrary \glslink{gene expression}{expression} cut-off.  This model is shown for comparison with Multivariate model I have chosen to develop since the Multivariate model, as shown in Figure 7, has better performance, allows the inclusion of correlation structure expected in \gls{gene expression} data for biological pathways, and could have variable gene function cut-offs.  

Figures 7 and 8, show that the Multivariate model which corrects this effect by specifying \gls{synthetic lethal} genes differently performs better in simulations, even with correlation structure expected to disrupt the \gls{synthetic lethal} detection.  There is indication in Figure 8 that correlation structure even improves the performance of simulations.  Although replicated across parameters, the difference in performance of simulations with correlated genes (with each \gls{synthetic lethal} partner) is marginal and the number of correlated genes is still vastly outnumbered by the total number of genes (20,000 modelling a complete mammalian \glspl{genome}).  Simulations with fewer total genes may show the impact of correlated genes more clearly, which is biologically plausible since some co-regulated pathways do involve a substantial proportion of the \glspl{genome}.

As indicated, the models behave as expected when performing better when simulated with higher sample size and fewer true \gls{synthetic lethal} genes.  As summarised in Figure 9, this behaviour occurs in simulation with all of the models discussed above.  The number of \gls{synthetic lethal} partners impacts performance with a sigmoidal decay where­­ higher sample size (albeit approaching the limit of feasible \gls{genomic}-scale projects) markedly delay decay of \gls{AUROC} towards random 0.5.  Therefore a large sample size is crucial for \gls{bioinformatics} \gls{synthetic lethal} discovery.  Only a small number of \gls{synthetic lethal} partners will be detectable with a gene-centric approach motivating pathway-centric approaches and accounting for \glslink{graph}{pathway} structure, which has shown be more reproducible between model organism experiments (Dixon et al. 2009).  However, whether potential false positives are more likely to be correlated genes or occur due to the sheer number of non-synthetic lethal genes (and multiple tests) is unclear.  The impact of correlation structure on the simulated data is explored in detail below in Figures 10-12 and the results of these simulations repeated is shown in Figure 13.    Figure 9.  Summary of effect of sample size and number of \gls{synthetic lethal} partners on performance of simulations for prediction of \glspl{synthetic lethal} by \gls{AUROC} on continuous scale (left) and as a barplot (right) showing that sample size (by colour) and number of \gls{synthetic lethal} partners (x axis) affects performance as expected in which was replicated across all 3 models discussed above.

\subsection{Simulated Expression Heatmaps}

In Figures 10-12 below, simulations are summarised with expected (Sigma) and generated (Correlation) structure of \gls{gene expression} patterns in the top figures.  The following line shows how the \glslink{gene expression}{expression} and gene function calls have been distributed with correlation structure and ordering samples (columns) to ensure a \gls{synthetic lethal} partner or query gene is active in each sample.

Figure 10.  Simulation for 1 SL partner (100 genes, 1000 samples)

Figure 11.  Simulation for 2 SL partners (100 genes, 1000 samples)

Figure 12.  Simulation for 3 SL partners (100 genes, 1000 samples)

As shown in the Figures 10-12, the correlation structure of the simulated \gls{gene expression} data (upper right) largely reflects the expected sigma matrix (upper left) used to specify the variation in the multivariate normal distribution with some variation due to low sampling error.  The Sigma and correlation matrices show blocks of correlated genes with each \gls{synthetic lethal} partner where there are 1, 2, or 3 \gls{synthetic lethal} partners in Figures 10, 11, and 12 respectively.  In the \gls{gene expression} heatmap (lower right) and associated discrete gene function calls based on a threshold of the 30\% quantile (lower left), the sample (column) ordering shows how samples were ordered so at least one \gls{synthetic lethal} gene is active in all query deficient samples.  The row (gene) ordering is based on a Chi-Square test statistic value and odds-ratio sign (with negative genes at the top), apart from Query gene at the top (with positive odds-ratio).  The Chi-Square values are shown on the outer colour bar on a log scale and the inner colour bar annotates the known gene class in the simulation: query (blue), \gls{synthetic lethal} (red), correlated (orange), and other (green).

With 1 \gls{synthetic lethal} partner, in Figure 10, the relationship between \gls{synthetic lethal} (and correlated genes with the Query gene is clear and detectable with Chi-Square test (as shown with the colour bars) as expected.  The relationship is clearer in the true \gls{synthetic lethal} partner showing that it should be distinguishable from confounding correlated genes.  With multiple \gls{synthetic lethal} genes, as shown in Figures 11 and 12, the true \gls{synthetic lethal} partner is less related to the \glslink{gene expression}{expression} profile of the Query gene and the co-loss under-representation is more difficult to detect since the number of co-occuring loss of \gls{synthetic lethal} genes expected (even in Query functional samples is low).  In these examples, the Chi-Square test still correctly identifies \gls{synthetic lethal} genes with the highest test statistic, although with a less well defined cut-off and it may not be reproducible (as discussed above).  This is consistent with more \gls{synthetic lethal} partners being able to recover function and ensure cell survival which is plausible given the evolutionary \glslink{genetic robustness}{robustness} of molecular networks, difficulty detecting individual gene pairs in \gls{gene expression} data, and rates of recurrence or drug resistance in cancer patients.  Therefore I have to consider cryptic \gls{synthetic lethal} genes compensating for Query and candidate \gls{synthetic lethal} partners due to higher-order \glslink{functional redundancy}{genetic redundancy}, cancer \gls{genomic} evolution and cellular heterogeneity.

\subsection{Replication Simulation Heatmap}

The declining performance in \gls{ROC} curves with more \gls{synthetic lethal} genes shows that the ability to robustly distinguish \gls{synthetic lethal} genes from other genes (including their correlated genes) declines as the \gls{synthetic lethal} genes do not consistently have a higher Chi-Square test statistic across replicate sampling simulations.  Although it is noted that increased sample size can compensate for this decline, increasing the number of expected co-loss events and sensitivity of the procedure.  The effect of total gene number, impact of correlation structure, and reproducibility of Chi-Square tests across replicate sampling simulations is explored below.

Figure 13 is composed of columns of side colour bars ordered by Chi-Square and odds-ratio sign (with Query in the corrected position at the bottom) as shown in Figures 10-12 with separate columns for repeated sampling with different parameters.  Figure 13 is an example of this visualisation of simulations for a small number of genes (100) and replicates (10 each for 1 to 10 \gls{synthetic lethal} partners).  Even in this small simulation, I see many of the processes discussed above summarised: the effect of number of \gls{synthetic lethal} genes on Chi-Square tests, power to detect \gls{synthetic lethal} and other correlated genes, decaying reproducibility and variation across replicates, lack of a clear threshold, and importance of directional conditions (e.g., odds-ratio sign) to distinguish \gls{synthetic lethal} and co-expressed genes.  This visualisation is an effective way to capture the simulation process and compare conditions which will be valuable for more complex correlation structure and comparison to public data Chi-Square distributions.
    
Figure 13.  Comparison of simulation across various parameters for sampling a multivariate normal model for 100 genes and 1000 samples with correlation structure with 10 replicates (columns) for each number of \gls{synthetic lethal} partners (1 to 10 in the top colour bar) with genes sorted by chi-squared value (and odds-ratio negative at the top) this shows preferential sorting of \gls{synthetic lethal} partners (red) and correlated (orange) genes near the top (on the left) for lower numbers of \gls{synthetic lethal} partners which becomes less clear or consistent across replicates for higher numbers of \gls{synthetic lethal} partners, reflected in less variation in chi-square values (shown in log-scale on the right) and lack of a clear prediction threshold, however positive odds-ratio genes show no preference except for the query gene associated itself as expected.

This framework may also be useful to compare different analyses of public data and infer the true number of \gls{synthetic lethal} partners from the distribution of test statistic scores.  With an effective visualisation, I can further explore more complex correlation structures (as shown in the supplementary Figures S1 and S2).  This will be important to develop simulated data as similar to empirical data as possible, to test whether \gls{synthetic lethal} and correlated genes are robustly detectable, and discover effective drug targets (which are repeatable across a cohort, tissues or species).  The impact of high-order \glspl{synthetic lethal}, genetic background and variation between replicates indicates that more care has to be taken interpreting experimental model systems and \glspl{genomic} analysis will be valuable to ensure candidate drug targets are suitable for clinical application.  I show below that this visualisation scales up and shows similar effects for number of \gls{synthetic lethal} genes in more replicates (Figure 14), more total genes (Figure 15), and both (Figure 16).
    
Figure 14.  Comparison of simulation across various parameters for sampling a multivariate normal model for 100 genes and 1000 samples with correlation structure with 100 replicates (columns) for each number of \gls{synthetic lethal} partners (1 to 10 in the top colour bar) with genes sorted by chi-squared value (and odds-ratio negative at the top) this shows preferential sorting of \gls{synthetic lethal} partners (red) and correlated (orange) genes near the top (on the left) for lower numbers of \gls{synthetic lethal} partners which becomes less clear or consistent across replicates for higher numbers of \gls{synthetic lethal} partners, reflected in less variation in chi-square values (shown in log-scale on the right) and lack of a clear prediction threshold, however positive odds-ratio genes show no preference except for the query gene associated itself as expected.  
   
Figure 15.  Comparison of simulation across various parameters for sampling a multivariate normal model for 1000 genes and 1000 samples with correlation structure with 10 replicates (columns) for each number of \gls{synthetic lethal} partners (1 to 10 in the top colour bar) with genes sorted by chi-squared value (and odds-ratio negative at the top) this shows preferential sorting of \gls{synthetic lethal} partners (red) and correlated (orange) genes near the top (on the left) for lower numbers of \gls{synthetic lethal} partners which becomes less clear or consistent across replicates for higher numbers of \gls{synthetic lethal} partners, reflected in less variation in chi-square values (shown in log-scale on the right) and lack of a clear prediction threshold, however positive odds-ratio genes show no preference except for the query gene associated itself as expected.
    
Figure 16.  Comparison of simulation across various parameters for sampling a multivariate normal model for 1000 genes and 1000 samples with correlation structure with 100 replicates (columns) for each number of \gls{synthetic lethal} partners (1 to 10 in the top colour bar) with genes sorted by chi-squared value (and odds-ratio negative at the top) this shows preferential sorting of \gls{synthetic lethal} partners (red) and correlated (orange) genes near the top (on the left) for lower numbers of \gls{synthetic lethal} partners which becomes less clear or consistent across replicates for higher numbers of \gls{synthetic lethal} partners, reflected in less variation in chi-square values (shown in log-scale on the right) and lack of a clear prediction threshold, however positive odds-ratio genes show no preference except for the query gene associated itself as expected.
\fi