\documentclass{article}
\title{Using \textbf{HDTD} to Analyze High-Dimensional Transposable Data: An Application in Genetics}
\author{Anestis Touloumis\footnote{Anestis.Touloumis@cruk.cam.ac.uk}, John C. Marioni and Simon Tavar\'e }
\date{}
\usepackage{Sweave}
\usepackage{amsmath}
\usepackage{graphicx, verbatim}
\setlength{\textwidth}{6.5in} 
\setlength{\textheight}{9in}
\setlength{\oddsidemargin}{0in} 
\setlength{\evensidemargin}{0in}
\setlength{\topmargin}{-1.5cm}
\begin{document}
\SweaveOpts{concordance=TRUE}
%\VignetteIndexEntry{HDTD to Analyze High-Dimensional Transposable Data}
\maketitle
\section{Introduction}
The R/Bioconductor package \textbf{HDTD} is designed to analyze high-dimensional transposable data. The term transposable data implies the following structural information in the dataset: 
\begin{itemize}
  \item the data for each sampling unit (e.g., subject/patient) can be written in a matrix,
  \item the rows and the columns in each matrix correspond to two distinct sets of variables,
  \item dependencies might occur among and/or between the row and column variables.
 \end{itemize} 
The term high-dimensional implies that the sample size (e.g., the number of subjects/patients) is a lot smaller than the total number of row and column variables.\\
Since the statistical methods implemented in \textbf{HDTD} were primarly motivated by studies in genetics, a microarray dataset is utilized to illustrate the functionality of the package. However, the use of \textbf{HDTD} is not limited to gene-expression experiments and we emphasize that it is suitable for analyzing datasets that satisfy the high-dimensional transposable data definition. 
\section{Mouse Dataset}
\textbf{HDTD} includes a subset of the tissue study described in Zahn et al. (2007). This dataset contains expression levels for $40$ mice. For each mouse, the expression levels of $46$ genes that belong to the vascular endothelial growth factor signalling pathway were measured across $9$ tissues (adrenal gland, cerebrum, hippocampus, kidney, liver, muscle, spinal cord, spleen and thymus). The experimental design satisfies the definition of transposable data because: i) the data for each mouse can be written in a matrix form, ii) the rows correspond to genes and the columns to multiple tissues, and iii) we do not expect the gene expression levels across the multiple tissues to vary independently.\\
The dataset is formatted as a single matrix with rows the $46$ genes and columns the $9 \times 40=360$ tissues.
<<>>=
library("HDTD")
data(VEGFmouse)
dim(VEGFmouse)
rownames(VEGFmouse)
@
Further, every 9 consecutive columns belong to the same mouse and the tissues are ordered in the same way for each mouse. For example, we can check the column variables for the first two subjects:
<<>>=
colnames(VEGFmouse)[1:18]
@
It is extremely important to provide datasets in this particular format when using \textbf{HDTD}. To accomplish this, write the data for each subject (mouse) in a matrix form while preserving the order of the row (genes) and column (tissues) variables. The final step is to create a single matrix by stacking column-wise the subject-specific matrices the one after the other.
\section{Mean Relationship of the Genes Across the Tissues}
The user can determine the mean relationship of the genes across tissues by testing and estimating the mean matrix. One intersting hypothesis to be tested is the conservation of the gene expression levels across the tissues, i.e., if the mean gene expression levels vector in the VEGF signaling pathway changes across the $9$ tissues:
<<>>=
meanmat.ts(VEGFmouse,40,group.sizes=9,voi="columns")
@
Since $p$-value $<0.001$, we have strong evidence against the null hypothesis that there is no tissue effect in the gene expression levels. To explore the mean gene expression level pattern across tissues in detail, additional tests can be carried out. We illustrate a more complicated hypothesis that requires data manipulation. Consider testing the hypothesis that the mean gene expression levels vector is constant only across the adrenal glands, the spleen, the kidney and the liver. To do this, we first need to place these $4$ tissues in a successive order in the dataset,   
<<>>=
colnames(VEGFmouse)[1:9]
VEGForder <- orderdata(VEGFmouse,40,order.cols=c(1,4,5,8,2,3,6,7,9)) 
colnames(VEGForder)[1:9]
@
and then to perform the test using the ordered dataset
<<>>=
meanmat.ts(VEGForder,40,group.sizes=c(4,1,1,1,1,1),voi="columns")
@
Note that we included $5$ additional column groups of size one in the \texttt{group.sizes} argument to reflect the fact that the mean gene expression levels vector in each of the remaining $5$ tissues remained unspecified. The null hypothesis is rejected, and hence we may conclude that the mean gene expression levels vector is not constant in the adrenal glands, the spleen, the kidney and the liver.\\ 
Apart from hypothesis testing, it is also important to estimate the mean relationship between the genes and the tissues. In this example the mean matrix seems to be unstructured and thus the mean gene expression levels in the $9$ tissues can be estimated via the sample mean matrix 
<<>>=
sample.mean <- meanmat.hat(VEGFmouse,40)
sample.mean
@
Note that the output preserves the order of the genes and the columns. For example, $0.8399$ is the average $\log_2$ intensity for gene "Akt1" in the adrenal gland based on $40$ mice. The mean matrix for the first 10 genes across the $9$ tissues is
<<>>=
head(round(sample.mean$estmeanmat,4),n=10)
@
\section{Dependence Structure of the Genes and of the Tissues}
The matrix-variate normal distribution is utilized to estimate two covariance matrices, one for the genes (rows) and the other for the multiple tissues (columns). We developed shrinkage estimators for both covariance matrices but we let the user decide if shrinkage is required to both, one or neither of these matrices. In principle, we recommend shrinking both covariance matrices in order to obtain well-defined and invertible covariance matrix estimators. The \texttt{covmat.hat} function provides the corresponding covariance estimators: 
<<>>=
estcovmat <- covmat.hat(datamat = VEGFmouse, N = 40, shrink = "both", centered = FALSE)
estcovmat
@
The output summarizes the results but the user can recover the full covariance matrix estimators. For example, the covariance matrix of the tissues is
<<>>=
round(estcovmat$cols.covmat,4)
@
Moreover, the user can study the gene-wise or tissue-wise correlation by using the \texttt{covmat.ts} function. For example, the identity and sphericity test 
<<>>=
covmat.ts(datamat = VEGFmouse, N = 40, voi = "columns", centered = FALSE)
@
suggest that the tissues might not be uncorrelated since the $p$-values of the sphericity and the identity test are both $<0.001$. 
\section{Citation}
<<>>=
citation("HDTD")
@
\section{Session Info}
<<>>=
sessionInfo()
@
\end{document}