\documentclass{article}

\usepackage[utf8]{inputenc}
\usepackage{grffile}
\usepackage{graphicx}
\usepackage{xcolor}
%\usepackage{hyperref}
\usepackage{subfigure}
\usepackage{caption}
%\usepackage[export]{adjustbox}% http://ctan.org/pkg/adjustbox

%\VignetteIndexEntry{compcodeR}

%\VignetteEngine{knitr::knitr}
<<style, eval=TRUE, echo=FALSE, results="asis">>=
BiocStyle::latex()
@
%\renewcommand{\thetitle}{"compcodeR"}

\title{Comparing methods for differential expression analysis of RNAseq data with the \Rpackage{compcodeR} package}
\author{Charlotte Soneson}
\date{\today}

\begin{document}
\maketitle 

\tableofcontents

\section{Introduction}
\Rpackage{compcodeR} is an \R{} package that provides an interface to several popular methods for differential expression analysis of RNAseq data and contains functionality for comparing the outcomes and performances of several differential expression methods applied to the same data set. The package also contains a function for generating synthetic RNAseq counts, using the simulation framework described in more detail in 
\begin{quote}
C Soneson and M Delorenzi: A comparison of methods for differential expression analysis of RNA-seq data. {\em BMC Bioinformatics} 2013, {\bf 14}:91.
\end{quote}
This vignette provides a tutorial on how to use the different functionalities of the \Rpackage{compcodeR} package. 

Currently, the differential expression interfaces provided in the package are restricted to comparisons between two conditions. However, many of the comparison functions are more general and can also be applied to test results from other contrast types, as well as to test results from other data types than RNAseq.

\vskip 0.5cm
\begin{center}
\fbox{
\begin{minipage}{.8\linewidth}
{\color{red}{Important!}} Since compcodeR provides interfaces to differential expression analysis methods implemented in other \R{} packages, take care to cite the appropriate references if you use any of the interface functions (see the reference manual for more information). Also be sure to check the code that was used to run the differential expression analysis, using e.g. the \Rfunction{generateCodeHTMLs} function (see Section~\ref{sec:comparison} for more information) to make sure that parameters etc. agree with your intentions and that there were no errors or serious warnings. 
\end{minipage}
}
\end{center}
\vskip 0.5cm



<<label = settings, echo = FALSE, eval = TRUE>>=
options(width = 55)
@

\section{The \Robject{compData} class}
Within the \Rpackage{compcodeR} package (version $\geq$ 0.2.0), data sets and results are represented as objects of the \Robject{compData} class. The functions in the package are still compatible with the list-based representation used in version 0.1.0, but we strongly encourage users to use the \Robject{compData} class, and all results generated by the package will be given in this format. If you have a data or result object generated with \Rpackage{compcodeR} version 0.1.0, you can convert it to a \Robject{compData} object using the \Rfunction{convertListTocompData} function. 

A \Robject{compData} object has at least three slots, containing the count matrix, sample annotations and a list containing at least an identifying name and a unique ID for the data set. It can also contain variable annotations, such as information regarding genes that are known to be differentially expressed. After performing a differential expression analysis, the \Robject{compData} object contains additional information, such as which method was used to perform the analysis, which settings were used and the gene-wise results from the analysis. More detailed information about the \Robject{compData} class are available in Sections~\ref{sec:dataformat} and \ref{sec:resultformat}.

\section{A sample workflow}
This section contains a sample workflow showing the main functionalities of the \Rpackage{compcodeR} package. We start by generating a synthetic count data set (Section~\ref{sec:simulation}), to which we then apply three different differential expression methods (Section~\ref{sec:diffexp}). Finally, we compare the outcome of the three methods and generate a report summarizing the results (Section~\ref{sec:comparison}).

\subsection{Simulating data}\label{sec:simulation}
The simulations are performed following the description by (Soneson and Delorenzi, 2013). As an example, we use the \Rfunction{generateSyntheticData} function to generate a synthetic count data set containing 12,500 genes and two groups of 5 samples each, where 10\% of the genes are simulated to be differentially expressed between the two groups (equally distributed between up- and downregulated in group 2 compared to group 1). Furthermore, the counts for all genes are simulated from a Negative Binomial distribution with the same dispersion in the two sample groups, and no outlier counts are introduced. We filter the data set by excluding only the genes with zero counts in all samples (i.e., those for which the total count is 0). This simulation setting corresponds to the one denoted $B_{625}^{625}$ in (Soneson and Delorenzi, 2013). The following code creates a \Robject{compData} object containing the simulated data set and saves it to a file named \texttt{"B\_625\_625\_5spc\_repl1.rds"}.

<<label = loadpackage, echo = TRUE, eval = TRUE, results = 'hide', message = FALSE>>=
library(compcodeR)
@
                                 
<<label = simulatedata, echo = TRUE, eval = FALSE, tidy = FALSE>>=
B_625_625 <- generateSyntheticData(dataset = "B_625_625", n.vars = 12500, 
                                   samples.per.cond = 5, n.diffexp = 1250, 
                                   repl.id = 1, seqdepth = 1e7, 
                                   fraction.upregulated = 0.5, 
                                   between.group.diffdisp = FALSE, 
                                   filter.threshold.total = 1, 
                                   filter.threshold.mediancpm = 0, 
                                   fraction.non.overdispersed = 0, 
                                   output.file = "B_625_625_5spc_repl1.rds")
@

The \Rfunction{summarizeSyntheticDataSet} function provides functionality to check some aspects of the simulated data by generating a report summarizing the parameters that were used for the simulation, as well as including some diagnostic plots. The report contains two MA-plots, showing the estimated average expression level and the log-fold change for all genes, indicating either the truly differentially expressed genes or the total number of outliers introduced for each gene. It also shows the log-fold changes estimated from the simulated data versus those underlying the simulation. The input to the \Rfunction{summarizeSyntheticDataSet} function can be either a \Robject{compData} object or the path to a file containing such an object. 
<<label = reportsimulated, echo = TRUE, eval = FALSE, tidy = FALSE>>=
summarizeSyntheticDataSet(data.set = "B_625_625_5spc_repl1.rds", 
                          output.filename = "B_625_625_5spc_repl1_datacheck.html")
@

Figure~\ref{fig:synthetic-data-check} shows two of the figures generated by this function. The left panel shows an MA plot with the genes colored by the true differential expression status. The right panel shows the relationship between the true log-fold changes between the two sample groups underlying the simulation, and the estimated log-fold changes based on the simulated counts. 

\begin{figure}
\subfigure{\resizebox{85mm}{!}{\includegraphics{compcodeR_check_figure/maplot-trueDEstatus.png}}}
\subfigure{\resizebox{85mm}{!}{\includegraphics{compcodeR_check_figure/logfoldchanges.png}}}
\caption{Example figures from the summarization report generated for a simulated data set. The left panel shows an MA plot, with the genes colored by the true differential expression status. The right panel shows the relationship between the true log-fold changes between the two sample groups underlying the simulation, and the estimated log-fold changes based on the simulated counts. Also here, the genes are colored by the true differential expression status.}
\label{fig:synthetic-data-check}
\end{figure}

\subsection{Performing differential expression analysis}\label{sec:diffexp}
We will now apply some of the interfaced differential expression methods to find genes that are differentially expressed between the two conditions in the simulated data set. This is done through the \Rfunction{runDiffExp} function, which is the main interface for performing differential expression analyses in \Rpackage{compcodeR}. The code below applies three differential expression methods to the data set generated above: the \Rfunction{voom} transformation from the \Biocpkg{limma} package (combined with \Rfunction{limma} for differential expression), the exact test from the \Biocpkg{edgeR} package, and a regular t-test applied directly on the count level data. 

<<label = rundiffexp1, echo = TRUE, eval = FALSE, tidy = FALSE>>=
runDiffExp(data.file = "B_625_625_5spc_repl1.rds", 
           result.extent = "voom.limma", Rmdfunction = "voom.limma.createRmd", 
           output.directory = ".", norm.method = "TMM")
runDiffExp(data.file = "B_625_625_5spc_repl1.rds", 
           result.extent = "edgeR.exact", Rmdfunction = "edgeR.exact.createRmd", 
           output.directory = ".", norm.method = "TMM", 
           trend.method = "movingave", disp.type = "tagwise")
runDiffExp(data.file = "B_625_625_5spc_repl1.rds", result.extent = "ttest", 
           Rmdfunction = "ttest.createRmd", 
           output.directory = ".", norm.method = "TMM")
@

The code needed to perform each of the analyses is provided in the \Rfunction{*.createRmd} functions. To obtain a list of all available \Rfunction{*.createRmd} functions (and hence of the available differential expression methods), we can use the \Rfunction{listcreateRmd()} function. Example calls are also provided in the reference manual (see the help pages for the \Rfunction{runDiffExp} function).
<<label = listcreatermd, echo = TRUE, eval = TRUE>>=
listcreateRmd()
@

You can also apply your own differential expression method to the simulated data (see Section~\ref{sec:ownmethods}).

\subsection{Comparing results from several differential expression methods}\label{sec:comparison}
Once we have obtained the results of the differential expression analyses (either by the methods interfaced by \Rpackage{compcodeR} or in other ways, see Section~\ref{sec:ownmethods}), we can compare the results and generate an HTML report summarizing the results from the different methods from many different aspects (see Section~\ref{sec:compmethods} for an overview of the comparison metrics). In \Rpackage{compcodeR}, there are two ways of invoking the comparison functionality; either directly from the command line or via a graphical user interface (GUI). The GUI is mainly included to avoid long function calls and provide a clear overview of the available methods and parameter choices. To use the GUI the R package \Rpackage{rpanel} must be installed (which assumes that BWidget is available). Moreover, the GUI may have rendering problems on certain platforms, particularly on small screens and if many methods are to be compared. Below, we will show how to perform the comparison using both approaches.

\subsubsection{The graphical user interface}
First, we consider the \Rfunction{runComparisonGUI} function, to which we provide a list of directories containing our result files, and the directory where the final report will be generated. Since the three result files above were saved in the current working directory, we can run the following code to perform the comparison:
<<label = runcomparison, echo = TRUE, eval = FALSE, tidy = FALSE>>=
runComparisonGUI(input.directories = ".", 
                 output.directory = ".", recursive = FALSE)
@
This opens a graphical user interface (Figure~\ref{fig:screenshot-gui}) where we can select which of the data files available in the input directories that should be included as a basis for the comparison, and which comparisons to perform (see Section~\ref{sec:compmethods} for an overview). Through this interface, we can also set p-value cutoffs for significance and other parameters that will govern the behaviour of the comparison. 

\vskip 0.5cm
\begin{center}
\fbox{
\begin{minipage}{.8\linewidth}
{\color{red}{Important!}} When you have modified a value in one of the textboxes (p-value cutoffs etc.), press Return on your keyboard to confirm the assignment of the new value to the parameter. Always check in the resulting comparison report that the correct values were recognized and used for the comparisons.
\end{minipage}
}
\end{center}
\vskip 0.5cm

\begin{figure}
\begin{center}
\subfigure{\resizebox{30mm}{!}{\includegraphics{"screenshot-gui-1}}}
\subfigure{\resizebox{120mm}{!}{\includegraphics{"screenshot-gui-2}}}
%\resizebox{120mm}{!}{\includegraphics{"screenshot-gui"}}
\caption{Screenshot of the graphical user interface used to select data set (left) and set parameters (right) for the comparison of differential expression methods. The available choices for the ``Data set'', ``DE methods'', ``Number of samples'' and ``Replicates'' are automatically generated from the \Robject{compData} objects available in the designated input directories. Only one data set can be used for the comparison. In the lower part of the window we can set (adjusted) p-value thresholds for each comparison method separately. For example, we can evaluate the true FDR at one adjusted p-value threshold, and estimate the TPR for another adjusted p-value threshold. We can also set the maximal number of top-ranked variables that will be considered for the false discovery curves. }
\label{fig:screenshot-gui}
\end{center}
\end{figure}

After the selections have been made, the function will perform the comparisons and generate an HTML report, which is saved in the designated output directory and automatically named \texttt{compcodeR\_report\_$<$timestamp$>$.html}. Please note that depending on the number of compared methods and the number of included data sets, this step may take several minutes. \Rpackage{compcodeR} will notify you when the report is ready (with a message \textbf{Done!} in the console). 

\vskip 0.5cm
\begin{center}
\fbox{
\begin{minipage}{.8\linewidth}
{\color{red}{Note!}} Depending on the platform you use to run \R{}, you may see a prompt ($">"$) in the console before the analysis is done. However, \Rpackage{compcodeR} will always notify you when it is finished, by typing \textbf{Done!} in the console.
\end{minipage}
}
\end{center}
\vskip 0.5cm

The comparison function will also generate a subdirectory called \texttt{compcodeR\_code}, where the R code used to perform each of the comparisons is detailed in HTML reports, and another subdirectory called \texttt{compcodeR\_figure}, where the plots generated in the comparison are saved. The code HTML reports can also be generated manually for a given result file (given that the \Robject{code} component is present in the \Robject{compData} object (see Section~\ref{sec:resultformat}), using the function \Rfunction{generateCodeHTMLs}. For example, to generate a report containing the code that was run to perform the t-test above, as well as the output from the \R{} console, we can write:
<<label = generatecode, echo = TRUE, eval = FALSE>>=
generateCodeHTMLs("B_625_625_5spc_repl1_ttest.rds", ".")
@
These reports are useful to check that there were no errors or warnings when running the differential expression analyses. 

\subsubsection{Direct command-line call}
Next, we show how to run the comparison by calling the function \Rfunction{runComparison} directly. In this case, we need to supply the function with a list of result files to use as a basis for the comparison. We can also provide a list of parameters (p-value thresholds, the differential expression methods to include in the comparison, etc.). The default values of these parameters are outlined in the reference manual. The following code provides an example, given the data generated above.
<<label = nogui-comparison, echo = TRUE, eval = FALSE, tidy = FALSE>>=
file.table <- data.frame(input.files = c("B_625_625_5spc_repl1_voom.limma.rds", 
                                         "B_625_625_5spc_repl1_ttest.rds", 
                                         "B_625_625_5spc_repl1_edgeR.exact.rds"), 
                         stringsAsFactors = FALSE)
parameters <- list(incl.nbr.samples = NULL, incl.replicates = NULL, 
                   incl.dataset = "B_625_625", incl.de.methods = NULL, 
                   fdr.threshold = 0.05, tpr.threshold = 0.05, 
                   typeI.threshold = 0.05, ma.threshold = 0.05, 
                   fdc.maxvar = 1500, overlap.threshold = 0.05, 
                   fracsign.threshold = 0.05, 
                   comparisons = c("auc", "fdr", "tpr", "ma", "correlation"))
runComparison(file.table = file.table, parameters = parameters, output.directory = ".")
@
By setting \Robject{incl.nbr.samples}, \Robject{incl.replicates} and \Robject{incl.de.methods} to \Robject{NULL}, we ask \Rpackage{compcodeR} to include all results provided in the \Robject{file.table}. By providing a vector of values for each of these variables, it is possible to limit the selection to a subset of the provided files. Please note that the values given to \Robject{incl.replicates} and \Robject{incl.nbr.samples} are matched with values of the \Rcode{info.parameters\$repl.id} and \Rcode{info.parameters\$samples.per.cond} slots in the data/result objects, and that the values given to \Robject{incl.de.methods} are matched with values stored in the \Rcode{method.names\$full.name} slot of the result objects. If the values do not match, the corresponding result object will not be considered in the comparison. Consult the package manual for the full list of comparison methods available for use with the \Rfunction{runComparison} function. Setting \Rfunction{parameters = NULL} implies that all results provided in the \Robject{file.table} are used, and that all parameter values are set to their defaults (see the reference manual). Note that only one dataset identifier can be provided to the comparison (that is, \Robject{parameters\$incl.dataset} must be a single string).


\section{Using your own data}
The \Rpackage{compcodeR} package provides a straightforward function for simulating count data (\Rfunction{generateSyntheticData}). However, it is easy to apply the interfaced differential expression methods to your own data, given that it is provided in a \Robject{compData} object (see Section~\ref{sec:dataformat} below for a description of the data format). You can use the function \Rfunction{check\_compData} to check that your object satisfies the necessary criteria to be fed into the differential expression methods. 

To create a \Robject{compData} object from a count matrix and a data frame with sample annotations, you can use the function \Rfunction{compData}. The following code provides a minimal example. Note that you need to provide a \Robject{dataset} name (a 'description' of the simulation settings) as well as a unique data set identifier \Robject{uID}, which has to be unique for each \Robject{compData} object (e.g., for each simulation instance, even if the same simulation parameters are used). 
<<create-compData, eval = TRUE, tidy = FALSE>>=
count.matrix <- matrix(round(1000*runif(4000)), 1000, 4)
sample.annot <- data.frame(condition = c(1, 1, 2, 2))
info.parameters <- list(dataset = "mytestdata", uID = "123456")
cpd <- compData(count.matrix = count.matrix, 
                sample.annotations = sample.annot, 
                info.parameters = info.parameters)
check_compData(cpd)
@

\section{Providing your own differential expression code}\label{sec:ownmethods}
The \Rpackage{compcodeR} package provides an interface for calling some of the most commonly used differential expression methods developed for RNAseq data. However, it is easy to incorporate your own favorite method. In principle, this can be done in one of two ways:
\begin{enumerate}
\item Write a \Rfunction{XXX.createRmd} function (where \Robject{XXX} corresponds to your method), similar to the ones provided in the package, which creates a \Robject{.Rmd} file containing the code that is run to perform the differential expression analysis. Then call this function through the \Rfunction{runDiffExp} function. When implementing your function, make sure that the output is a \Robject{compData} object, structured as described in Section~\ref{sec:resultformat}. The \Rfunction{*.createRmd} functions provided in the package take the following input arguments:
\begin{itemize}
\item \Robject{data.path} -- the path to the \Robject{.rds} file containing the \Robject{compData} object to which the differential expression will be applied.
\item \Robject{result.path} -- the path to the \Robject{.rds} file where the resulting \Robject{compData} object will be stored.
\item \Robject{codefile} -- the name of the code (with extension .Rmd) where the code will be stored.
\item any arguments for setting parameters of the differential expression analysis.
\end{itemize}

\item Run the differential expression analysis completely outside the \Rpackage{compcodeR} package, and save the result in a \Robject{compData} object with the slots described in Section~\ref{sec:resultformat}. 
\end{enumerate}

You can use the \Rfunction{check\_compData\_results} function to check if your object satisfies the necessary conditions for being used as the output of a differential expression analysis and compared to results obtained by other methods with \Rpackage{compcodeR}.

\section{The format of the data and result objects}
This section details the format of the data and result objects generated and used by the \Rpackage{compcodeR} package. Both objects are of the \Robject{compData} class. The format guidelines below must be be followed if you apply the functions in the package to a data set of your own, or to differential expression results generated outside the package. Note that for most of the functionality of the package, the objects should be saved separately to files with a \Robject{.rds} extension, and the path of this object is provided to the functions. 

\subsection{The data object}\label{sec:dataformat}
The \Robject{compData} data object used by \Rpackage{compcodeR} is an S4 object with the following slots:
\begin{itemize}
\item \Robject{count.matrix} [\Robject{class matrix}] ({\color{blue}{mandatory}}) -- the count matrix, with rows representing genes and columns representing samples.
\item \Robject{sample.annotations} [\Robject{class data.frame}] ({\color{blue}{mandatory}}) -- sample annotations. Each row corresponds to one sample, and each column to one annotation. The data objects generated by the \Rfunction{generateSyntheticData} function have two annotations:
\begin{itemize}
\item \Robject{condition} [\Robject{class character} or \Robject{numeric}] ({\color{blue}{mandatory}}) -- the class for each sample. Currently the differential expression implementations in the package supports only two-group comparisons, hence the \Robject{condition} should have exactly two unique values.
\item \Robject{depth.factor} [\Robject{class numeric}] -- the depth factor for each sample. This factor, multiplied by the \Robject{sequencing depth}, corresponds to the target library size of the sample when simulating the counts. 
\end{itemize}
\item \Robject{info.parameters} [\Robject{class list}] -- a list detailing the parameter values that have been used for the simulation. 
\begin{itemize}
\item \Robject{dataset} [\Robject{class character}] ({\color{blue}{mandatory}}) -- the name of the data set.
\item \Robject{samples.per.cond} [\Robject{class numeric}] -- the number of samples in each of the two conditions.
\item \Robject{n.diffexp} [\Robject{class numeric}] -- the number of genes that are simulated to be differentially expressed between the two conditions.
\item \Robject{repl.id} [\Robject{class numeric}] -- a replicate number, which can be set to differentiate between different instances generated with the exact same simulation settings.
\item \Robject{seqdepth} [\Robject{class numeric}] -- the "base" sequencing depth that was used for the simulations. For each sample, it is modified by multiplication with a value sampled uniformly between \Robject{minfact} and \Robject{maxfact} to generate the actual sequencing depth for the sample. 
\item \Robject{minfact} [\Robject{class numeric}] -- the lower bound on the values used to multiply the \Robject{seqdepth} to generate the actual sequencing depth for the individual samples.  
\item \Robject{maxfact} [\Robject{class numeric}] -- the upper bound on the values used to multiply the \Robject{seqdepth} to generate the actual sequencing depth for the individual samples.  
\item \Robject{fraction.upregulated} [\Robject{class numeric}] -- the fraction of the differentially expressed genes that were simulated to be upregulated in condition 2 compared to condition 1. Must be in the interval [0, 1]. The remaining differentially expressed genes are simulated to be downregulated in condition 2.
\item \Robject{between.group.diffdisp} [\Robject{class logical}] -- whether the counts from the two conditions were simulated with different dispersion parameters.
\item \Robject{filter.threshold.total} [\Robject{class numeric}] -- the filter threshold that is applied to the total count across all samples. All genes for which the total count is below this number have been excluded from the data set in the simulation process.
\item \Robject{filter.threshold.mediancpm} [\Robject{class numeric}] -- the filter threshold that is applied to the median count per million (cpm) across all samples. All genes for which the median cpm is below this number have been excluded from the data set in the simulation process.
\item \Robject{fraction.non.overdispersed} [\Robject{class numeric}] -- the fraction of the genes in the data set that are simulated according to a model without overdispersion (i.e., a Poisson model).
\item \Robject{random.outlier.high.prob} [\Robject{class numeric}] -- the fraction of extremely high random outlier counts introduced in the simulated data sets. Please consult (Soneson \& Delorenzi 2013) for a detailed description of the different types of outliers. 
\item \Robject{random.outlier.low.prob} [\Robject{class numeric}] -- the fraction of extremely low random outlier counts introduced in the simulated data sets. Please consult (Soneson \& Delorenzi 2013) for a detailed description of the different types of outliers. 
\item \Robject{single.outlier.high.prob} [\Robject{class numeric}] -- the fraction of extremely high single outlier counts introduced in the simulated data sets. Please consult (Soneson \& Delorenzi 2013) for a detailed description of the different types of outliers. 
\item \Robject{single.outlier.low.prob} [\Robject{class numeric}] -- the fraction of extremely low single outlier counts introduced in the simulated data sets. Please consult (Soneson \& Delorenzi 2013) for a detailed description of the different types of outliers. 
\item \Robject{effect.size} [\Robject{class numeric}] -- the degree of differential expression, i.e., a measure of the minimal effect size. Alternatively, a vector of provided effect sizes for each of the variables. 
\item \Robject{uID} [\Robject{class character}] ({\color{blue}{mandatory}}) -- a unique identification number given to each data set. For the data sets simulated with \Rpackage{compcodeR}, it consists of a sequence of 10 randomly generated alpha-numeric characters.
\end{itemize}
\item \Robject{filtering} [\Robject{class character}] -- a summary of the filtering that has been applied to the data.
\item \Robject{variable.annotations} [\Robject{class data.frame}] -- annotations for each of the variables in the data set. Each row corresponds to one variable, and each column to one annotation. No annotation is mandatory, however some of them are necessary if the differential expression results for the data are going to be used for certain comparison tasks. The data sets simulated within \Rpackage{compcodeR} have the following named variable annotations:
\begin{itemize}
\item \Robject{truedispersions.S1} [\Robject{class numeric}] -- the true dispersions used in the simulations of the counts for the samples in condition 1.
\item \Robject{truedispersions.S2} [\Robject{class numeric}] -- the true dispersions used in the simulations of the counts for the samples in condition 2.
\item \Robject{truemeans.S1} [\Robject{class numeric}] -- the true mean values used in the simulations of the counts for the samples in condition 1.
\item \Robject{truemeans.S2} [\Robject{class numeric}] -- the true mean values used in the simulations of the counts for the samples in condition 2.
\item \Robject{n.random.outliers.up.S1} [\Robject{class numeric}] -- the number of extremely high random outliers introduced for each gene in condition 1.
\item \Robject{n.random.outliers.up.S2} [\Robject{class numeric}] -- the number of extremely high random outliers introduced for each gene in condition 2.
\item \Robject{n.random.outliers.down.S1} [\Robject{class numeric}] -- the number of extremely low random outliers introduced for each gene in condition 1.
\item \Robject{n.random.outliers.down.S2} [\Robject{class numeric}] -- the number of extremely low random outliers introduced for each gene in condition 2.
\item \Robject{n.single.outliers.up.S1} [\Robject{class numeric}] -- the number of extremely high single outliers introduced for each gene in condition 1.
\item \Robject{n.single.outliers.up.S2} [\Robject{class numeric}] -- the number of extremely high single outliers introduced for each gene in condition 2.
\item \Robject{n.single.outliers.down.S1} [\Robject{class numeric}] -- the number of extremely low single outliers introduced for each gene in condition 1.
\item \Robject{n.single.outliers.down.S2} [\Robject{class numeric}] -- the number of extremely low single outliers introduced for each gene in condition 2.
\item \Robject{M.value} [\Robject{class numeric}] -- the estimated log2-fold change between conditions 1 and 2 for each gene. These values were estimated using the \Rpackage{edgeR} package.
\item \Robject{A.value} [\Robject{class numeric}] -- the estimated average expression in conditions 1 and 2 for each gene. These values were estimated using the \Rpackage{edgeR} package.
\item \Robject{truelog2foldchanges} [\Robject{class numeric}] -- the ``true'' log2-fold changes between conditions 1 and 2 for each gene, based on the parameters used for simulation.
\item \Robject{upregulation} [\Robject{class numeric}] -- a binary annotation (0/1) indicating which genes are simulated to be upregulated in condition 2 compared to condition 1. (Upregulated genes are indicated with a 1.)
\item \Robject{downregulation} [\Robject{class numeric}] -- a binary annotation (0/1) indicating which genes are simulated to be downregulated in condition 2 compared to condition 1. (Downregulated genes are indicated with a 1.)
\item \Robject{differential.expression} [\Robject{class numeric}] ({\color{blue}{mandatory}} for many comparisons, such as the computation of false discovery rates, true positive rates, ROC curves, false discovery curves etc.)-- a binary annotation (0/1) indicating which genes are simulated to be differentially expressed in condition 2 compared to condition 1. In other words, the sum of \Robject{upregulation} and \Robject{downregulation}.
\end{itemize}
\end{itemize}
To apply the functions of the package to a \Robject{compData} object of the type detailed above, it needs to be saved to a file with extension \Robject{.rds}. To save the object \Robject{cpd} to the file \Robject{saveddata.rds}, simply type 
<<label = savedata, echo = TRUE, eval = FALSE>>=
saveRDS(cpd, "saveddata.rds")
@

\subsection{The result object}\label{sec:resultformat}
When applying one of the differential expression methods interfaced through \Rpackage{compcodeR}, the \Robject{compData} object is extended with some additional slots. These are described below. 
\begin{itemize}
\item \Robject{analysis.date} [\Robject{class character}] -- the date and time when the differential expression analysis was performed.
\item \Robject{package.version} [\Robject{class character}] -- the version of the packages used in the differential expression analysis.
\item \Robject{method.names} [\Robject{class list}] ({\color{blue}{mandatory}}) -- a list containing the names of the differential expression method, which should be used to identify the results in the comparison. Contains two components:
\begin{itemize}
\item \Robject{short.name} [\Robject{class character}] -- a short name, used for convenience.
\item \Robject{full.name} [\Robject{class character}] -- a fully identifying name of the differential expression method, including e.g. version numbers and parameter values. In the method comparisons, the results will be grouped based on the \Robject{full.name}.
\end{itemize}
\item \Robject{code} [\Robject{class character}] -- the code that was used to run the differential expression analysis, in R markdown (.Rmd) format.
\item \Robject{result.table} [\Robject{class data.frame}] ({\color{blue}{mandatory}}) -- a data frame containing the results of the differential expression analysis. Each row corresponds to one gene, and each column to a value generated by the analysis. The precise columns will depend on the method applied. The following columns are used by at least one of the methods interfaced by \Rpackage{compcodeR}:
\begin{itemize}
\item \Robject{pvalue} [\Robject{class numeric}] -- the nominal p-values.
\item \Robject{adjpvalue} [\Robject{class numeric}] -- p-values adjusted for multiple comparisons.
\item \Robject{logFC} [\Robject{class numeric}] -- estimated log-fold changes between the two conditions.
\item \Robject{score} [\Robject{class numeric}] ({\color{blue}{mandatory}}) -- the score that will be used to rank the genes in order of significance. Note that {\color{red}{high scores always signify differential expression, that is, a strong association with the predictor}}. For example, for methods returning a nominal p-value the score is generally obtained as \Robject{1 - pvalue}.
\item \Robject{FDR} [\Robject{class numeric}] -- the false discovery rate estimate.
\item \Robject{posterior.DE} [\Robject{class numeric}] -- the posterior probability of differential expression.
\item \Robject{prob.DE} [\Robject{class numeric}] -- The conditional probability of differential expression.
\item \Robject{lfdr} [\Robject{class numeric}] -- the local false discovery rate.
\item \Robject{statistic} [\Robject{class numeric}] -- a test statistic from the differential expression analysis.
\item \Robject{dispersion.S1} [\Robject{class numeric}] -- dispersion estimate in condition 1.
\item \Robject{dispersion.S2} [\Robject{class numeric}] -- dispersion estimate in condition 2.
\end{itemize}
For many of the comparison methods, the naming of the result columns is important. For example, the p-value column must be named \Robject{pvalue} in order to be recognized by the comparison method computing type I error. Similarly, either an \Robject{adjpvalue} or an \Robject{FDR} column must be present in order to apply the comparison methods requiring adjusted p-value/FDR cutoffs. If both are present, the \Robject{adjpvalue} column takes precedence over the \Robject{FDR} column. 
\end{itemize}

To be used in the comparison function, the result \Robject{compData} object must be saved to a \Robject{.rds} file. 

\section{The evaluation/comparison methods}\label{sec:compmethods}
This section provides an overview of the methods that are implemented in \Rpackage{compcodeR} for comparing differential expression results obtained by different methods. The selection of which methods to apply is made through a graphical user interface that is opened when the \Rfunction{runComparisonGUI} function is called (see Figure~\ref{fig:screenshot-gui}). Alternatively, the selection of methods can be supplied to the \Rfunction{runComparison} function directly, to circumvent the GUI. 

\subsection{ROC (one replicate/all replicates)}\label{sec:roc}
This method computes a ROC curve for either a single representative among the replicates of a given data set, or for all replicates (in separate plots). The ROC curves are generated by plotting the true positive rate (TPR, on the y-axis) versus the false positive rate (FPR, on the x-axis) when varying the cutoff on the \Robject{score} (see the description of the result object above). A good ranking method gives a ROC curve which passes close to the upper left corner of the plot, while a bad method gives a ROC curve closer to the diagonal. Calculation of the ROC curves requires that the differential expression status of each gene is provided (see the description of the data object in Section~\ref{sec:dataformat}). Figure~\ref{fig:roc} shows an example of ROC curves for a single replicate.

\begin{center}
\resizebox{140mm}{!}{\includegraphics[trim=0 0 170mm 0, clip = true]{compcodeR_figure/rocone-1.png}}%
\captionof{figure}{}\label{fig:roc}%
\end{center}

\subsection{AUC}\label{sec:auc}
This comparison method computes the area under the ROC curve (see section~\ref{sec:roc}) and represents the result in boxplots, where each box summarizes the results for one method across all replicates of a data set. A good ranking method gives to a large value of the AUC. This requires that the differential expression status of each gene is provided (see the description of the data object in Section~\ref{sec:dataformat}). Figure~\ref{fig:auc} shows an example.

%\noindent%
%\begin{minipage}{\linewidth}
%\makebox[\linewidth]{%
%  \includegraphics[keepaspectratio=true,scale=0.6]{compcodeR_figure/auc.png}}
%\captionof{figure}{}% only if needed
%\label{fig:auc}
%\end{minipage}

\begin{center}
\resizebox{150mm}{!}{\includegraphics{compcodeR_figure/auc.png}}%
\captionof{figure}{}\label{fig:auc}%
\end{center}

\subsection{Type I error}
This approach computes the type I error (the fraction of the genes that are truly non-differentially expressed that are called significant at a given nominal p-value threshold) for each method and sample size separately, and represent it in boxplots, where each box summarizes the results for one method across all replicates of a data set. This requires that the differential expression status of each gene is provided (see the description of the data object above) and that the \Robject{pvalue} column is present in the \Robject{result.table} of the included result objects (see the description of the result object above). An example is given in Figure~\ref{fig:typeIerror}. The dashed vertical line represents the imposed nominal p-value threshold, and we wish that the observed type I error is lower than this threshold.

\begin{center}
\resizebox{150mm}{!}{\includegraphics{compcodeR_figure/typeIerror.png}}%
\captionof{figure}{}\label{fig:typeIerror}%
\end{center}

\subsection{FDR}
Here, we compute the observed false discovery rate (the fraction of the genes called significant that are truly non-differentially expressed at a given adjusted p-value/FDR threshold) for each method and sample size separately, and represent the estimates in boxplots, where each box summarizes the results for one method across all replicates of a data set. This requires that the differential expression status of each gene is provided (see the description of the data object above) and that the \Robject{adjpvalue} or \Robject{FDR} column is present in the \Robject{result.table} of the included result objects (see the description of the result object above). An example is given in Figure~\ref{fig:fdr}. The dashed vertical line represents the imposed adjusted p-value threshold (that is, the level at which we wish to control the false discovery rate).

\begin{center}
\resizebox{150mm}{!}{\includegraphics{compcodeR_figure/fdr.png}}%
\captionof{figure}{}\label{fig:fdr}%
\end{center}

\subsection{FDR as a function of expression level}
Instead of looking at the overall false discovery rate, this method allows us to study the FDR as a function of the average expression level of the genes. For each data set, the average expression levels are binned into 10 bins of equal size (i.e., each containing 10\% of the genes), and the FDR is computed for each of them. The results are shown by means of boxplots, summarizing the results across all replicates of a data set. Figure~\ref{fig:fdrvsexpr} shows an example.

\begin{center}
\resizebox{150mm}{!}{\includegraphics{compcodeR_figure/fdrvsexpr-1.png}}%
\captionof{figure}{}\label{fig:fdrvsexpr}%
\end{center}

\subsection{TPR}
With this approach, we compute the observed true positive rate (the fraction of the truly differentially expressed genes that are called significant at a given adjusted p-value/FDR threshold) for each method and sample size separately, and represent the estimates in boxplots, where each box summarizes the results for one method across all replicates of a data set. This requires that the differential expression status of each gene is provided (see the description of the data object above) and that the \Robject{adjpvalue} or \Robject{FDR} column is present in the \Robject{result.table} of the included result objects (see the description of the result object above). Figure~\ref{fig:tpr} shows an example. There is often a trade-off between achieving a high TPR (which is desirable) and controlling the number of false positives, and hence the TPR plots should typically be interpreted together with the type I error and/or FDR plots.

\begin{center}
\resizebox{150mm}{!}{\includegraphics{compcodeR_figure/tpr.png}}%
\captionof{figure}{}\label{fig:tpr}%
\end{center}

\subsection{False discovery curves (one replicate/all replicates)}
This choice plots the false discovery curves, depicting the number of false discoveries (i.e., truly non-differentially expressed genes) that are encountered when traversing the ranked list of genes ordered in decreasing order by the \Robject{score} column of the \Robject{result.table} (see the description of the result object above). As for the ROC curves, the plots can be made for a single replicate or for all replicates (each in a separate figure). Figure~\ref{fig:fdc} shows an example for a single replicate. A well-performing method is represented by a slowly rising false discovery curve (in other words, few false positives among the top-ranked genes).

\begin{center}
\resizebox{140mm}{!}{\includegraphics[trim=0 0 170mm 0, clip = true]{compcodeR_figure/fdcone-1.png}}%
\captionof{figure}{}\label{fig:fdc}%
\end{center}

\subsection{Fraction/number of significant genes}
We can also compute the fraction and/or number of the genes that are called significant at a given adjusted p-value/FDR threshold for each method and sample size separately, and represent the estimates in boxplots, where each box summarizes the results for one method across all replicates of a data set. This requires that the \Robject{adjpvalue} or \Robject{FDR} column is present in the \Robject{result.table} of the included result objects (see the description of the result object above). An example is given in Figure~\ref{fig:fracsign}. Since this plot does not incorporate information about the number of the identity of the truly differentially expressed genes, we can not conclude whether a large or small fraction of significant genes is preferable. Thus, the plot is merely an indication of which methods are more liberal (giving more significantly differentially expressed genes) and which are more conservative.

\begin{center}
\resizebox{150mm}{!}{\includegraphics{compcodeR_figure/fracsign.png}}%
\captionof{figure}{}\label{fig:fracsign}%
\end{center}

\subsection{Overlap, one replicate}
For each pair of methods, this approach computes the overlap between the sets of genes called differentially expressed by each of them at a given adjusted p-value/FDR threshold. Only one representative of the replicates of a data set is used. This requires that the \Robject{adjpvalue} or \Robject{FDR} column is present in the \Robject{result.table} of the included result objects (see the description of the result object above). Since the overlap depends on the number of genes called differentially expressed by the different methods, a normalized overlap measure (the Sorensen index) is also provided (see Section~\ref{sec:sorensen} below). The results are represented as a table of the sizes of the overlaps between each pair of methods. 

\subsection{Sorensen index, one replicate}\label{sec:sorensen}
For each pair of methods, this approach computes a normalized pairwise overlap value between the sets of genes called differentially expressed by each of them at a given adjusted p-value/FDR threshold. Only one representative of the replicates of a data set is used. This requires that the \Robject{adjpvalue} or \Robject{FDR} column is present in the \Robject{result.table} of the included result objects (see the description of the result object above). The Sorensen index (also called the Dice coefficient) for two sets $A$ and $B$ is given by $$S(A,B)=\frac{2|A\cap B|}{|A|+|B|},$$where $|\cdot|$ denotes the cardinality of a set (the number of elements in the set). The results are represented in the form of a table, and also as a colored heatmap, where the color represents the degree of overlap between the sets of differentially expressed genes found by different methods. An example is shown in Figure~\ref{fig:sorensen}.

\begin{center}
\resizebox{100mm}{!}{\includegraphics{compcodeR_figure/sorensen-1.png}}%
\captionof{figure}{}\label{fig:sorensen}%
\end{center}

\subsection{MA plot}
With this approach we can construct MA-plots (one for each differential expression method) for one replicate of a data set, depicting the estimated log-fold change (on the y-axis) versus the estimated average expression (on the x-axis). The genes that are called differentially expressed are marked with color. Figure~\ref{fig:maplot} shows an example.

\begin{center}
\resizebox{150mm}{!}{\includegraphics{compcodeR_figure/maplot-1.png}}%
\captionof{figure}{}\label{fig:maplot}%
\end{center}

\subsection{Spearman correlation}
This approach computes the pairwise Spearman correlation between the gene scores obtained by different differential expression methods (the \Robject{score} component of the \Robject{result.table}, see the description of the result object above). The correlations are visualized by means of a color-coded table and used to construct a dissimilarity measure for hierarchical clustering of the differential expression methods. Examples are given in Figures~\ref{fig:correlation-1} and \ref{fig:correlation-2}.

\begin{center}
\resizebox{100mm}{!}{\includegraphics{compcodeR_figure/correlation-11.png}}%
\captionof{figure}{}\label{fig:correlation-1}%
\end{center}

\begin{center}
\resizebox{80mm}{!}{\includegraphics{compcodeR_figure/correlation-12.png}}%
\captionof{figure}{}\label{fig:correlation-2}%
\end{center}

\subsection{Score distribution vs number of outliers}
This method plots the distribution of the gene scores obtained by different differential expression methods (the \Robject{score} component of the \Robject{result.table}, see the description of the result object above) as a function of the number of outliers imposed for the genes. The distributions are visualized by means of violin plots. Figure~\ref{fig:scorevsoutlier} shows an example, where there are no outliers introduced. Recall that a high score corresponds to more significant genes. In the presence of outliers, the score distribution plots can be used to examine whether the presence of an outlier count for a gene shifts the score distribution towards higher or lower values.

\begin{center}
\resizebox{150mm}{!}{\includegraphics{compcodeR_figure/scorevsoutlier-1.png}}%
\captionof{figure}{}\label{fig:scorevsoutlier}%
\end{center}

\subsection{Score distribution vs average expression level}
This method plots the gene scores (the \Robject{score} component of the \Robject{result.table}, see the description of the result object above) as a function of the average expression level of the gene. An example is shown in Figure~\ref{fig:scorevsexpr}. Recall that high scores correspond to more significant genes. The colored line shows the trend in the relationship between the two variables by means of a loess fit. These plots can be used to examine whether highly expressed genes tend to have (e.g.) higher scores than lowly expressed genes, and more generally to provide a useful characteristic of the methods. However, it is not clear what would be the "optimal" behavior. 

\begin{center}
\resizebox{150mm}{!}{\includegraphics{compcodeR_figure/scorevsexpr-1.png}}%
\captionof{figure}{}\label{fig:scorevsexpr}%
\end{center}

\subsection{Score vs 'signal' for genes expressed in only one condition}
This method plots the gene scores (the \Robject{score} component of the \Robject{result.table}, see the description of the result object above) as a function of the 'signal strength' for genes that are detected in only one of the two conditions. The signal strength can be defined either as the average (log-transformed) normalized pseudo-count in the condition where the gene is expressed, or as the signal-to-noise ratio in this condition, that is, as the average log-transformed normalized pseudo-count divided by the standard deviation of the log-transformed normalized pseudo-count. Typically, we would expect that the score increases with the signal strength. Figure~\ref{fig:scorevssignal} shows an example.

\begin{center}
\resizebox{150mm}{!}{\includegraphics{compcodeR_figure/scorevssignal-1.png}}%
\captionof{figure}{}\label{fig:scorevssignal}%
\end{center}

\subsection{Matthew's correlation coefficient}
Matthew's correlation coefficient summarizes the result of a classification task in a single number, incorporating the number of false positives, true positives, false negatives and true negatives. A correlation coefficient of +1 indicates perfect classification, while a correlation coefficient of -1 indicates perfect "anti-classification", i.e., that all objects are misclassified.

\section{Session info}
<<session-info, eval = TRUE, echo = FALSE, results="asis">>=
toLatex(sessionInfo())
@

\end{document}