\documentclass[12pt,a4paper]{article}

\title{Manual for ProbABEL v0.4.3}
\author{\emph{Current Programmers:} Lennart Karssen$^{1,2}$, Maarten
  Kooyman$^2$, \\
  Yurii Aulchenko$^{1,3}$ \\
  \emph{Former Programmers:} Maksim Struchalin
  \\
  \\
  $^{1}${\small YuriiA Consulting} \\
  $^{2}${\small Erasmus MC, Rotterdam}\\
  $^{3}${\small Institute of Cytology and Genetics SD RAS, Novosibirsk}
}
\date{April 1, 2014}


\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{textcomp}

\usepackage[svgnames]{xcolor}
\definecolor{webgreen}{rgb}{0,.5,0}

\usepackage{verbatim}

\usepackage{listings}
\lstloadlanguages{Bash}
\definecolor{lstbgcolor}{rgb}{0.9,0.9,0.9}
\lstset{
  tabsize=4,
  rulecolor=,
  basicstyle=\ttfamily,
  upquote=true,
  columns=fixed,
  showstringspaces=false,
  extendedchars=true,
  breaklines=true,
  breakatwhitespace,
  prebreak = \raisebox{0ex}[0ex][0ex]{\ensuremath{\hookleftarrow}},
  frame=single,
  showtabs=false,
  showspaces=false,
  showstringspaces=false,
  keywordstyle=\color[rgb]{0,0,1},
  commentstyle=\color[rgb]{0,0.4,0},
  stringstyle=\color[rgb]{0.5,0,1},
  basicstyle=\footnotesize\ttfamily,
  backgroundcolor=\color{lstbgcolor},
}

\usepackage{titleref}
\usepackage{amsmath}
\usepackage{makeidx}
\usepackage[pdftex,hyperfootnotes=false,pdfpagelabels]{hyperref}
\hypersetup{%
  linktocpage=false, % If true the page numbers in the toc are links
  % instead of the section headings.
  pdfstartview=FitH,% pdfstartpage=3,
  breaklinks=true, pageanchor=true, %
  pdfpagemode=UseOutlines, plainpages=false, bookmarksnumbered, %
  bookmarksopen=true, bookmarksopenlevel=1, hypertexnames=true, %
  pdfhighlight=/O, %hyperfootnotes=true,%nesting=true,%frenchlinks,%
  pdfauthor={\textcopyright\ Y.~Aulchenko, M.~Struchalin, L.C.~Karssen},
  pdfsubject={ProbABEL manual},
  colorlinks=true, urlcolor=blue, linkcolor=blue, citecolor=webgreen %
}
% get the links to the figures and tables right:
\usepackage[all]{hypcap} % to be loaded after hyperref package

% lowercase letters as footnote numerals (to avoid confusion with powers).
\renewcommand{\thefootnote}{\alph{footnote}}

\DeclareMathOperator{\var}{\mathbf{var}}

\makeindex

\newcommand{\PA}{\texttt{ProbABEL}}
\newcommand{\GA}{\texttt{GenABEL}}
\newcommand{\DA}{\texttt{DatABEL}}

\begin{document}
\maketitle
\tableofcontents

%\begin{abstract}
%This document describes ProbABEL package for analysis
%of genome-wide imputed data \ldots
%\end{abstract}

\section{Motivation}

Many statistical and experimental techniques, such as imputations and
high-throughput sequencing, generate data which are informative for
genome-wide association analysis and are probabilistic in the nature.

When we work with directly genotyped markers using such techniques as
SNP or microsatellite typing, we would normally know the genotype of
a particular person at a particular locus with very high degree of
confidence, and, in case of biallelic marker, can state whether
the genotype is $AA$, $AB$ or $BB$.

On the other hand, when dealing with imputed or high-throughput
sequencing data, the genotypic status of the person is known with a
much lower confidence. Instead of dealing with
known genotypes we work with a probability distribution that is based
on observed information, and we have estimates that the true underlying
genotype is either $AA$, $AB$ or $BB$. The degree of confidence about
the real status is measured with the probability distribution
$\{P(AA), P(AB), P(BB)\}$.

Several techniques may be applied to analyse such data. The most
simplistic approach would be to pick up the genotype with highest
probability, i.e.~$\max_g[P(g=AA), P(g=AB), P(g=BB)]$ and then
analyse the data as if directly typed markers were used. The
disadvantage of this approach is that it does not take into
account the probability distribution -- i.e.~the uncertainty
about the true genotypic status. Such
analysis is statistically wrong: the estimates of association
parameters (regression coefficients, odds or hazard ratios, etc.)
are biased, and the bias becomes more pronounced with greater
probability distribution uncertainty (entropy).

One of the solutions that generate unbiased estimates
of association parameters and takes the
probability distribution into account is achieved by
performing association analysis by means of regression of the
outcome of interest onto estimated genotypic probabilities.

The \PA{} package was designed to perform such regression
in a fast, memory-efficient and, consequently, genome-wide feasible manner.
Currently, \PA{} implements linear and logistic regression,
as well as the Cox proportional hazards model. The corresponding analysis
programs are called \texttt{palinear},  \texttt{palogist},
and \texttt{pacoxph}.

For more information, please have a look at the GenABEL project
website at \url{http://www.genabel.org}. The \PA{}-specific home
page is located at
\url{http://www.genabel.org/packages/ProbABEL}. For user support
questions, please turn to our forum at
\url{http://forum.genabel.org}. Bugs in \PA{} can be reported in the
GenABEL project bug tracker at
\url{https://r-forge.r-project.org/tracker/index.php?group_id=505&atid=2058}.

\section{Obtaining and installing \PA}
\label{sec:obtaininstall}
\PA{} is a tool that is mostly used on computers running the Linux
operating system. We try to publish binary packages for Windows as
well, but these aren't tested. We strongly suggest using \PA{} on
Linux.

\subsection{Precompiled packages}
\PA{} can be obtained in several ways:
\begin{itemize}
\item If you are using Ubuntu Linux and have administrative rights on
  the machine you can add the GenABEL PPA to your APT configuration
  and install it from there. The PPA can be found at
  \url{https://launchpad.net/~l.c.karssen/+archive/genabel-ppa}. Instructions
  on how to add the PPA can also be found there.
\item If your computer runs Debian Linux\footnote{At the moment \PA{}
    is only available in Debian testing and unstable.} (and you have
  administrative rights on it), you can install ProbABEL like this:
  \begin{lstlisting}
user@server:~$ apt-get install probabel
  \end{lstlisting}
\item Zip files with pre-compiled binaries (if available) can be found
  on the ProbABEL web page
  (\url{http://www.genabel.org/packages/ProbABEL}).
\item If you don't fall in any of the aforementioned
  categories\footnote{We know that many people have use Red Hat Linux,
    CentOS, Scientific Linux or any other Red Hat
    derivative. Unfortunately we haven't got \texttt{rpm} files
    yet. Any help in creating those will be highly appreciated}, you
  can install \PA{} manually by downloading the source code of the
  latest version from the website and compiling it yourself. This will
  be explained in section~\ref{sec:obtain}.
\end{itemize}


\subsection{Obtaining the source code and compiling it yourself}
\label{sec:obtain}
If you can't use any of the aforementioned pre-compiled packages, you
can download the source code of \PA{} yourself, compile it and run it
from your own home directory. This section details the steps you need
to take. More information can be found in the \texttt{doc/INSTALL}.

On the \href{http://www.genabel.org/packages/probabel}{\PA{}} website
you can find the link to the latest version of the source code of \PA{}
in a \texttt{tar.gz} file\footnote{The \texttt{tar.gz} file archive
  format is the most commonly used format for distributing source code
  on Linux/UNIX systems. These are compressed files, similar to
  \texttt{zip} files.}. A \texttt{.asc} file with the same base name
as the source code archive is also provided. This file contains a
so-called GPG signature of the \texttt{tar.gz} file. Using this file
and the \texttt{gpg} tool you can verify the authenticity of the
source code by typing this command on the command line of a Linux
shell\footnote{The \$ sign indicates the end of the command line
  prompt. You don't need to type it.}:
\begin{lstlisting}[]
user@server:~$ gpg --verify probabel-0.4.3.tar.gz.asc
gpg: Signature made Thu Jan  2 02:38:25 2014 CET using DSA key ID DA9CD509
gpg: Good signature from "L.C. Karssen (GPG key for personal stuff) <lennart@karssen.org>"
gpg:                 aka "L.C. Karssen (My GMail address) <l.c.karssen@gmail.com>"
\end{lstlisting}
Notice the ``Good signature'' message and the fact that the package was
signed by Lennart Karssen, the ProbABEL maintainer. If a malicious
hacker would have replaced the source code file (for example with one
including a virus), he won't be able to sign the package using the
same key (with key ID DA9CD509). If, for some reason, the
\texttt{tar.gz} file has changed (e.g.~by such a hacker or because
the file didn't get downloaded correctly) you will see output like
this (notice the ``BAD signature'' message):
\begin{lstlisting}[]
user@server:~$ gpg --verify probabel-0.4.2.tar.gz.asc
gpg: Signature made Thu Jan  2 02:38:25 2014 CET using DSA key ID DA9CD509
gpg: BAD signature from "L.C. Karssen (GPG key for personal stuff) <lennart@karssen.org>"
user@server:~$
\end{lstlisting}

Before continuing, it is important to mention that \PA{} can make use
of the EIGEN library\footnote{EIGEN is a library for fast matrix
  multiplication.}. We strongly recommend compiling \PA with EIGEN as
it will speed up your analyses considerably. Moreover, we plan to
remove the non-EIGEN part of the code in a future release. So, go to
\url{http://eigen.tuxfamily.org} and download the \texttt{tar.gz} file
of the latest version of EIGEN (3.2.1 at the time of writing). Extract
the files:
\begin{lstlisting}
user@server:~$ tar -xzf 3.2.1.tar.gz
\end{lstlisting}
This will create a directory called \texttt{eigen-eigen} followed by a
series of letters and digits. For simplicity we rename it to EIGEN
\begin{lstlisting}
user@server:~$ mv eigen-eigen-6b38706d90a9 EIGEN
\end{lstlisting}

Now it's time to extract the \PA{} source code and move into the
directory that is created:
\begin{lstlisting}
user@server:~$ tar -xzf probabel-0.4.3.tar.gz
user@server:~$ cd probabel-0.4.3
\end{lstlisting}
With the following command we will indicate where the EIGEN files can
be found and where we want to install \PA{}. Let's install in a
subdirectory of your home directory,
e.g.~\texttt{/home/yourusername/ProbABEL}:
\begin{lstlisting}
user@server:~$ ./configure \
   --prefix=/home/yourusername/ProbABEL/ \
   --with-eigen-include-path=/home/yourusername/EIGEN
\end{lstlisting}
This will be followed by a series of checks to see if all tools
required for compilation and installation are present on your
system. If you don't see any warnings you can continue to
compile\footnote{Compilation is the process of converting the source
  files containing human readable program code to a files with machine
  readable instructions.} the code using the \texttt{make}
command\footnote{If you work on a machine with multiple processors (or
  processor cores), which should be the case on modern servers, but
  also on most PCs, you can speed up the process by adding this number
  to the \texttt{-j} option. For example for four cores run
  \texttt{make -j4}.} The next step will check the compiled code,
after wich you install the program, documentation and examples to the
directory you specified previously with the \texttt{--prefix} argument
to the \texttt{./configure} command.
\begin{lstlisting}
user@server:~$ make
user@server:~$ make check
user@server:~$ make install
\end{lstlisting}
Note that each of these steps will scroll a lot of output on the
screen. Please watch it for any warnings or errors. Please ask any
questions on \href{http://forum.genabel.org/}{our support forum}.

If all went well you will find the executable programs
(\texttt{palinear}, \texttt{palogist}, and \texttt{pacoxph}) in the
directory \texttt{/home/yourusername/ProbABEL/bin/}. You are now ready
to analyse your data!

\section{Input files}
\PA{} takes three files as input: a file containing SNP
information (e.g.~the MLINFO file of MaCH), a file with genome- or
chromosome-wide predictor information (e.g.~the MLDOSE or MLPROB file
of MaCH or \texttt{minimac}),
and a file containing the phenotype of interest and covariates.

Optionally, the map information can be supplied (e.g.~the "legend"
files of HapMap).

The dose/probability file may be supplied in filevector format
in which case \PA{} will operate much faster, and
in low-RAM mode (approx.~128 MB). See the R libraries \GA{} and
\DA{} on how to convert MaCH and IMPUTE2 files to
filevector format (functions: \texttt{mach2databel()} and
\texttt{impute2databel()}, respectively).

\subsection{SNP information file}
\label{ssec:infoin}
In the simplest scenario, the SNP information file is an MLINFO
file generated by MaCH/\texttt{minimac}. This must be a space or tab-delimited file
containing SNP name, coding for allele 1 and 2 (e.g.~A, T, G or C),
frequency of allele 1, minor allele frequency and two quality
metrics (``Quality'', the average maximum posterior probability and
``Rsq'', the proportion of variance decrease after imputations).

Actually, for \PA{}, it (almost) does not matter what is written in
this file -- this information is simply copied to the output. However,
\textbf{it is critical} that the number of columns is
seven\footnote{This means that for \texttt{minimac} output files the number of
  columns needs to be reduced. This can be done using e.g.~GAWK or
  \texttt{cut}.} and the number of lines in the file is equal to the
number of SNPs in the corresponding DOSE file (plus one for the header
line). Also make sure that the ``Rsq'' column contains values $>0$,
otherwise you will end up with $\beta$'s set to \texttt{nan}.

The example of SNP information file content follows here (also
to be found in \texttt{ProbABEL/examples/test.mlinfo})

\verbatiminput{test.mlinfo}
Note that a header line is present in the file. The file describes
five SNPs.

\subsection{Genomic predictor file}
\label{ssec:dosein}

Again, in the simplest scenario this is an MLDOSE or MLPROB file
generated by MaCH and/or \texttt{minimac}.  Such file starts with two
special columns plus, for each of the SNPs under consideration, a
column containing the estimated allele 1 dose (MLDOSE).  In an MLPROB
file, two columns for each SNP correspond to posterior probability
that person has two ($P_{A_1A_1}$) or one ($P_{A_1A_2}$) copies of
allele 1.  The first ``special'' column is made of the sequential id,
followed by an arrow followed by study ID (the one specified in the
MaCH input files). The second column contains the method keyword
(e.g.~``MLDOSE'').

An example of the few first lines of an MLDOSE file for
five SNPs described in SNP information file follows here (also
to be found in the file \texttt{ProbABEL/examples/test.mldose})

\verbatiminput{short_test.mldose}
%\immediate\write18{head -n 10 INSTALL > tmp.txt}


\textbf{The order of SNPs in the SNP information file and DOSE or PROB
  file must be the same}. This should be the case if you just used
MaCH/\texttt{minimac} outputs.
Consequently, the number of columns in the genomic predictor file
must be the same as the number of lines in the SNP information file
plus one in the case of a DOSE file. Similarly, for a PROB file the
number of columns must be equal to two times the number of SNPs plus
1.

The dose/probability file may be supplied in filevector format
(\texttt{.fvi} and \texttt{.fvd} files) in which case
\texttt{ProbABEL} will operate much faster, and in low-RAM mode
(approx.~128 MB). On the command line simply specify the \texttt{.fvi}
file as argument for the \texttt{--dose} option
(cf.~section~\ref{sec:runanalysis} for more information on the options
accepted by \texttt{ProbABEL}). See the R libraries \GA{} and
\DA{} on how to convert MaCH and IMPUTE files to
filevector format (functions: \texttt{mach2databel()} and
\texttt{impute2databel()}, respectively).


\subsection{Phenotypic file}
\label{ssec:phenoin}

The phenotypic data file contains phenotypic data, but also specifies
the analysis model. There is a header line, specifying the variable
names.  The first column should contain personal study IDs. It is
assumed that \textbf{both the total number and the order of these IDs
  are exactly the same as in the genomic predictor (DOSE/PROB) file
  described in previous section}. This is not difficult to arrange
using e.g.~\texttt{R}; an example is given in the \texttt{examples}
directory.

\textbf{Missing data should be coded with 'NA', 'N' or 'NaN' codes.} Any
other coding will be converted to some number which will be used in
analysis! E.g.~coding missing as '-999.9' will result in an analysis which
will consider -999.9 as indeed a true measurements of the trait/covariates.

In the case of linear or logistic regression (programs
\texttt{palinear} and \texttt{palogist}, respectively), the second
column specifies the trait under analysis, while the third, fourth,
etc.~provide information on covariates to be included into analysis.
As an example, a few lines of a phenotypic information file designed
for linear regression analysis follow here (also to be found in
\texttt{examples/height.txt})

\verbatiminput{short_height.txt}
Note again that the order of IDs is the same in the MLDOSE file
and the phenotypic data file. The model specified by this file is
\begin{equation*}
\textrm{height} \sim \mu + \textrm{sex} + \textrm{age},
\end{equation*}
where $\mu$ is the intercept.

Clearly, you can for example include \texttt{sex $\times$ age} interaction terms by
specifying another column having a product of sex and age here.

For logistic regression, it is assumed that in the second column cases are
coded as ``1'' and controls as ``0''. An couple of example lines of a phenotypic
information file designed for logistic regression analysis follow here (also
to be found in \texttt{examples/logist\_data.txt})

\verbatiminput{short_logist_data.txt}

You can see that in the first 10 people, there are three cases, as indicated
by ''chd'' equal to one. The model specified by this file
is
\begin{equation*}
\textrm{chd} \sim \mu + \textrm{sex} + \textrm{age} +
\textrm{other cov}.
\end{equation*}

In case of the Cox proportional hazards model, the composition of the
phenotypic input file is a bit different. In the second column and
third column, you need to specify the outcome in terms of follow-up
time (column two) and event (column three, ``1'' if an event occurred
and zero if censored). Columns starting from four (inclusive) specify
covariates to be included into the analysis. An example few lines of
a phenotypic information file designed for the Cox proportional hazards model
analysis follow here (also to be found in
\texttt{examples/coxph\_data.txt})

\verbatiminput{short_coxph_data.txt}
You can see that for the first ten people, the event occurs for three of
them, while for the other seven there is no event during the follow-up
time, as indicated by the ``chd'' column. Follow-up time is specified in the preceding
column. The covariates included into the model are age (presumably
at baseline), sex and ``othercov''; thus the model, in terms of
\texttt{R/survival} is
\begin{equation*}
\textrm{Surv}(\textrm{fuptime\_chd}, \textrm{chd})
\sim \textrm{sex} + \textrm{age} + \textrm{other cov}.
\end{equation*}

\subsection{Optional map file}
If you would like map information (e.g.~base pair position) to
be included in your outputs, you can supply a map file. These follow
HapMap "legend" file format. For example, for the five SNPs we considered
the map-file may look like (example can be found in
\texttt{examples/test.map})

\verbatiminput{test.map}

The order of the SNPs in the map file should follow that in the SNP information
file. Only information from the second column -- the SNP location -- is
actually used to generate the output.

\section{Running an analysis}
\label{sec:runanalysis}
To run linear regression, you should use the program called
\texttt{palinear}; for logistic analysis use \texttt{palogist}, and
for the Cox proportional hazards model use \texttt{pacoxph} (all are
found in the \texttt{bin/} directory after you have compiled the
program).

There are in total 11 command line options you can specify to the
\PA{} analysis functions \texttt{palinear} or \texttt{palogist}. If
you run either program with the \texttt{--help} option, you will get a
short explanation to the command line options:
\begin{verbatim}
user@server:~$ palogist --help
probabel v. 0.4.3
(C) Yurii Aulchenko, Lennart C. Karssen, Maksim Struchalin, EMCR

Using EIGEN version 3.1.2 for matrix operations

Usage: src/palogist options
Options:
         --pheno   : phenotype file name
         --info    : information (e.g. MLINFO) file name
         --dose    : predictor (e.g. MLDOSE/MLPROB) file name
         --map     : [optional] map file name
         --nids    : [optional] number of people to analyse
         --chrom   : [optional] chromosome (to be passed to output)
         --out     : [optional] output file name (default is regression.out.txt)
         --skipd   : [optional] how many columns to skip in the predictor
                     (dose/prob) file (default 2)
         --ntraits : [optional] how many traits are analysed (default 1)
         --ngpreds : [optional] how many predictor columns per marker
                      (default 1 = MLDOSE; else use 2 for MLPROB)
         --separat : [optional] character to separate fields (default is space)
         --score   : use score test
         --no-head : do not report header line
         --allcov  : report estimates for all covariates (large outputs!)
         --interaction: Which covariate to use for interaction with SNP analysis (default is no interaction, 0)
         --interaction_only: like previous but without covariate acting in interaction with SNP (default is no interaction, 0)
         --mmscore : score test in samples of related individuals. File with inverse of variance-covariance matrix (for palinear) or inverse correlation (for palogist) as input parameter
         --robust  : report robust (aka sandwich, aka Hubert-White) standard errors
         --help    : print help
\end{verbatim}
More information on the options can also be found in the manual page
of each of the programs. For example, to access the manual page for
\texttt{palinear} on a Linux system run\footnote{Use the \texttt{q}
  key to quit the man page viewer.}:
\begin{verbatim}
user@server$ man palinear
\end{verbatim}


\subsection{Basic analysis options}
However, for a simple run only three options are mandatory, which
specify the necessary files needed to run the regression analysis.

These options are \texttt{--dose} (or \texttt{-d}), specifying the
genomic predictor/MLDOSE file described in section \ref{ssec:dosein};
\texttt{--pheno} (or \texttt{-p}), specifying the phenotypic data file
described in section \ref{ssec:phenoin}; and \texttt{--info} (or
\texttt{-i}), specifying the SNP information file described in section
\ref{ssec:infoin}.

If you change to the \texttt{examples} directory you can run
an analysis of height by running
\begin{verbatim}
palinear -p height.txt -d gtdata/test.mldose -i gtdata/test.mlinfo
\end{verbatim}
Output from the analysis will be stored in the
\texttt{regression.out.csv} file.
The analysis of a binary trait (e.g.~chd) can be run with
\begin{verbatim}
palogist -p logist_data.txt -d gtdata/test.mldose \
    -i gtdata/test.mlinfo
\end{verbatim}
To run a Cox proportional hazards model, try
\begin{verbatim}
pacoxph -p coxph_data.txt -d gtdata/test.mldose \
    -i gtdata/test.mlinfo
\end{verbatim}

Please have a look at the shell script files \texttt{example\_qt.sh},
\texttt{example\_bt.sh} and \texttt{example\_all.sh} to have
a better overview of the analysis options.

To run an analysis with MLPROB files, you need specify the MLPROB file
with the \texttt{-d} option and also specify that there are two
genetic predictors per SNP, e.g.~you can run a linear model with
\begin{verbatim}
palinear -p height.txt -d gtdata/test.mlprob -i gtdata/test.mlinfo \
    --ngpreds=2
\end{verbatim}

When using genomic predictor files (dosages or probabilities) stored
in filevector (a.k.a.~DatABEL) format (i.e.~a combination of
\texttt{.fvi} and \texttt{.fvd} files) you can specify these like you
would with ordinary text files. This is how the previous example would
change:
\begin{verbatim}
palinear -p height.txt -d gtdata/test.mlprob.fvi -i gtdata/test.mlinfo \
    --ngpreds=2
\end{verbatim}

\subsection{Advanced analysis options}
The option \texttt{--interaction} allows you to include interaction
between SNPs and any covariate. If for example your model is
\begin{equation*}
  \textrm{trait} \sim \textrm{sex} + \textrm{age} + \textrm{SNP},
\end{equation*}
running the program with the option \texttt{--interaction=2} will model
\begin{equation*}
  \textrm{trait} \sim \textrm{sex} + \textrm{age} + \textrm{SNP} +
  \textrm{age} \times \mathrm{SNP}.
\end{equation*}
The option \texttt{--robust} allows you to compute so-called
``robust'' (a.k.a.~``sandwich'', a.k.a.~Hubert-White) standard errors
(cf.~section \ref{sec:methodology} ``Methodology'' for details).

With the option \texttt{--mmscore} a score test for association
between a trait and genetic polymorphisms in samples of related
individuals is performed. For quantitative traits (\texttt{palinear})
a file with the inverse of the variance-covariance matrix goes as input
parameter with that option, e.g.~\texttt{--mmscore <filename>}. The
file has to contain the first column with id names exactly like in
phenotype file, BUT OMITTING people with no measured phenotype. The
rest is a matrix. The phenotype file in case of using the
\texttt{--mmscore} argument may contain any amount of covariates (this
is different from previous versions). The first column contains id
names, the second the trait. The others are covariates.

For binary traits (\texttt{palogist}) the file should contain the
inverse of the correlation matrix. \textbf{Note: this is an
  experimental feature!} This matrix can be obtained through (in
GenABEL notation):
\begin{verbatim}
  h2.obj$InvSigma * h2.obj$h2an$estimate[length(h2.obj$h2an$estimate)]
\end{verbatim}
where \texttt{h2.obj} is the object returned by GenABEL's
\texttt{polygenic()}. The GenABEL documentation explains this
procedure in more detail.

An example of how a polygenic object estimated by \GA{} can be used
with ProbABEL is provided in \texttt{examples/mmscore.R}

Though technically \texttt{--mmscore} allows for inclusion of multiple
covariates, these should be kept to minimum as this is a score test. We
suggest that any covariates explaining an essential proportion of
variance should be fit as part of \GA{}'s
\texttt{polygenic} procedure.

%Option \texttt{--interaction\_only} is like \texttt{--interaction} but does not
%include in the model the main effect of the covariate, which is acting in
%interaction with SNP. This option is useful when running \texttt{--mmscore},
%in which case the main effect should normally estimated in the polygenic
%model and only the interaction term in the \PA{} analysis.

\subsection{Running multiple analyses at once: \texttt{probabel}}
The \texttt{bin/probabel} script is a handy wraper for the \PA{}
functions. To start using it the configuration file
\texttt{etc/probabel\_config.cfg.example} needs to be edited and
renamed to \texttt{etc/probabel\_config.cfg}. The configuration file
consists of five columns, separated by commas. Each column except the
first is a pattern for files produced by MaCH or \texttt{minimac}
(imputation tools). The column named ``cohort'' is an identifying name
of a population (``STUDY\_1'' in the example), the column
``info\_path'' is the full path to ``info'' files, including a pattern
where the chromosome number has been replaced by
\texttt{\_.\_chr\_.\_}. In case the imputations were run on chunks of
chromosomes, the pattern \texttt{\_.\_chunk\_.\_} will be replaced
with the corresponding chunk number. Chunk numbers should start at 1
for each chromosome. The columns ``dose\_path'', ``prob\_path'' and
``legend\_path'' are paths and patterns for ``dose'', ``prob'' and
``legend'' files, respectively. These also need to include the pattern
for the chromosome as used in the column for the ``info'' files.
Empty lines and lines starting with a \texttt{\#} are ignored.

The \texttt{make install} installation procedure should have set all
paths in the \texttt{probabel} script correctly. If that is not the
case you will have to change the variable \texttt{\$config} in the
script to point to the full path of the configuration file and the
variables \texttt{\$base\_path} and \texttt{@anprog} to point the full
path to the \PA{} scripts.


\section{Output file format}
Let us consider what comes out of the linear regression analysis
described in the previous section. After the analysis has run, in
the output file you will find something like
\begin{small}
\verbatiminput{short_height_base_add.out.txt.save}
\end{small}

Here, only the first three lines of output have been shown. Note that lines
starting with \texttt{+>} are actually the ones continuing the
previous line -- they have just been wrapped so we can see
these long lines.

The header provides a short description of what can be found in a
specific column. The first column provides the SNP name and
next six are descriptions which were taken directly from the
SNP information file. Therefore, these describe allele frequencies and
the quality in your total imputations, not necessarily in the data under
analysis.

In contrast, starting with the next column, named \texttt{n},
the output concerns the data analysed. Column 8 (\texttt{n}) tells the
number of subjects for whom complete phenotypic information was
available. At this point, unless you have complete measurements on all
subjects, you should feel alarmed if the number here is exactly the
number of people in the file -- this may indicate you did not code
missing values according to \PA{} format ('NA', 'NaN', or 'N').

The next column, nine (``Mean\_predictor\_allele''), gives the estimated
frequency of the predictor allele (\texttt{A1}) in subjects with complete
phenotypic data.

If the \texttt{--chrom} option was used, in the next column you will
find the value specified by this option. If \texttt{--map} option was
used, in the subsequent column you will find map location taken from
the map-file. The subsequent columns provide coefficients of
regression of the phenotype onto the genotype ($\beta$), corresponding
standard errors ($\text{SE}_\beta$), and the $\chi^2$ test value based
on the likelihood ratio test. Note that for the additive, recessive,
dominant and overdominant genetic models this is a $\chi^2$ of 1
degree of freedom (df), whereas for the genotypic model this is a
$\chi^2$ of 2df. If the \texttt{--mmscore} option is used, the
$\chi^2$ values are calculated from the Wald statistic
(1df)\footnote{For the genotypic (2df) model the $\chi^2$ values can't
  be simply calculated using the Wald statistic, consequently the
  $\chi^2$ values are set to \texttt{nan}. A fix for this needs to be
  implemented still.}.

\section{Preparing input files}
After installing \PA{} you can find the \texttt{prepare\_data.R} file
in the \texttt{scripts} directory. It is an R script that arranges
phenotypic data in the right format. Please read this script for
details.

\section{Memory use and performance}
Maximum likelihood regression is implemented in
\PA{}. With 6,000 people and 2.5 million SNPs, a
genome-wide scan is completed in less that an hour for a linear model
with 1-2 covariates and overnight for logistic regression or the Cox
proportional hazards model (figures for a PC bought back in 2007).

Memory may be an issue with \PA{} if you use MaCH/\texttt{minimac}
text dose/probability files, e.g. for large chromosomes, such as
chromosome one consumed up to 5 GB of RAM with 6,000 people.

We suggest that the genotype dosage/probability file is to be supplied
in filevector format in which case \PA{} will operate about 2-3 times
faster, and in low-RAM mode (approx.~128 MB). See the R libraries
\GA{} and \DA{} on how to convert MaCH and IMPUTE files to filevector
format (functions: \texttt{mach2databel()} and
\texttt{impute2databel()}, respectively).

When the \texttt{--mmscore} option is used, the analysis takes
more time.

\section{Methodology}
\label{sec:methodology}
\subsection{Analysis of population-based data}
\subsubsection{Linear regression assuming normal distribution}
Standard linear regression theory is used to estimate coefficients of
regression and their standard errors. We assume a linear model with
expectation
\begin{equation}
  E[\mathbf{Y}] = \mathbf{X}\, \boldsymbol{\beta}
\label{eq:expectation}
\end{equation}
and variance-covariance matrix
$$
\mathbf{V} = \sigma^2 \mathbf{I},
$$
where $\mathbf{Y}$ is the vector of phenotypes of interest,
$\mathbf{X}$ is the design matrix, $\boldsymbol{\beta}$ is the vector
of regression parameters, $\sigma^2$ is the variance and $\mathbf{I}$
is the identity matrix.

The maximum likelihood estimates (MLEs) for the regression parameters
are given by
\begin{equation}
  \hat{\beta} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{Y}
\end{equation}
and the MLE of the residual variance is
\begin{equation}
  \hat{\sigma}^2 = \frac{(\mathbf{Y} - \mathbf{X}\hat{\beta})^T
    (\mathbf{Y} - \mathbf{X}\hat{\beta})} {N-r_X},
\end{equation}
where $N$ is the number of observations and $r_X$ is the rank of
$\mathbf{X}$ (i.e.~the number of columns of the design matrix).

The variance-covariance matrix for the parameter estimates under
alternative hypothesis can be computed as
\begin{equation}
\var_{\hat{\beta}} = \hat{\sigma}^2 (\mathbf{X}^T\mathbf{X})^{-1}.
\end{equation}

For the $j$-the element $\hat{\beta}(j)$ of the vector of estimates the standard
error under the alternative hypothesis is given by the square root of the
corresponding diagonal element of the above matrix, $\var_{\hat{\beta}}(jj)$,
and the Wald test can be computed with
$$
T^2(j) = \frac{ \hat{\beta}(j)^2 }{ \var_{\hat{\beta}}(jj) },
$$
which asymptotically follows the $\chi^2$ distribution with one degree of
freedom under the null hypothesis.

When testing significance for more than one parameter simultaneously,
several alternatives are available. Let us first partition the vector of
parameters into two components, $\beta = (\beta_g,\beta_x)$, and our
interest is testing the parameters contained in $\beta_g$ (SNP effects),
while $\beta_x$ (e.g. effects of sex, age, etc.) are considered nuisance
parameters. Let us define the vector of the parameters of interest
which are fixed to certain values under the null hypothesis as
$\beta_{g,0}$.

Firstly, the likelihood ratio test can be obtained with
$$
LRT = 2 \left(\mathrm{logLik}(\hat{\beta}_g,\hat{\beta}_x) -
\mathrm{logLik}(\beta_{g,0},\hat{\beta}_x) \right),
$$
which under the null hypothesis is asymptotically distributed as
$\chi^2$ with the
number of degrees of freedom equal to the number of parameters specified
by $\beta_g$. Assuming the normal distribution, the log-likelihood of a
model specified by the vector of parameters $\beta$ and residual variance
$\sigma^2$ can be computed as
$$
\mathrm{logLik}(\beta,\sigma^2) =
-\frac{1}{2} \left( N \cdot \log_e \sigma^2 +
(\mathbf{Y} - \beta \mathbf{X})^T (\mathbf{I}/\sigma^2) (\mathbf{Y} -
\beta \mathbf{X}) \right).
$$

Secondly, the Wald test can be used; for that the inverse
variance-covariance matrix of $\hat{\beta}_g$ should be computed as
$$
\var_{\hat{\beta}_g}^{-1} = \var_{\hat{\beta}}^{-1}(g,g) -
\var_{\hat{\beta}}^{-1}(g,x) \left(
  \var_{\hat{\beta}}^{-1}(x,x) \right)^{-1}
  \var_{\hat{\beta}}^{-1}(x,g),
$$
where $\var_{\hat{\beta}}^{-1}(a,b)$ correspond to sub-matrices of the
inverse of the variance-covariance matrix of $\hat{\beta}$, involving
either only parameters of interest $(g,g)$, nuisance parameters
$(x,x)$ or combination of these $(x,g)$, $(g,x)$.

The Wald test statistics is then computed as
$$
W^2 = (\hat{\beta}_g - \beta_{g,0})^T \,
\var_{\hat{\beta}_g}^{-1} (\hat{\beta}_g - \beta_{g,0}),
$$
which asymptotically follows the $\chi^2$ distribution with the number
of degrees of freedom equal to the number of parameters specified by
$\beta_g$. The Wald test generally is computationally easier than the
LRT, because it avoids estimation of the model specified by the
parameter's vector $(\beta_{g,0},\hat{\beta}_x)$.

Lastly, similar to the Wald test, the score test can be performed by use
of $\var_{(\beta_{g,0},\hat{\beta}_x)}$ instead of $\var_{\hat{\beta}}$.

%% Comparative advantages of these testing approaches in the context of GWAS will be discussed in
%% ''\titleref{implementation}'' section (\ref{implementation}).

\subsubsection{Logistic regression}
For logistic regression, the procedure to obtain
parameters estimates, their variance-covariance matrix, and tests are
similar to these outlined above with several modifications.

The expectation of the binary trait is defined as the expected
probability of the event as defined by the logistic
function
$$
E[\mathbf{Y}] = \pi = \frac{ 1 }{ 1 + e^{-(\mathbf{X}\beta)} }.
$$
The estimates of the parameters are obtained not in one
step, as is the case of the linear model, but using an iterative
procedure (iteratively re-weighted least squares). This
procedure is not described here for the sake of brevity.

The log-likelihood of the data is computed using the
binomial probability formula:
$$
\mathrm{logLik}(\beta) = \mathbf{Y}^T \log_e \pi + (\mathbf{1} -
\mathbf{Y})^T \log_e (\mathbf{1}-\pi),
$$
where $\log_e \pi$ is a vector obtained by taking the natural
logarithm of every value contained in the vector $\pi$.

\subsubsection{Robust variance-covariance matrix of parameter estimates}
For a linear model, these are computed using the equation
$$
\var_r = (\mathbf{X}^T\mathbf{X})^{-1} (\mathbf{X}^T\mathbf{R}\mathbf{X})
(\mathbf{X}^T\mathbf{X})^{-1},
$$
where $\mathbf{R}$ is a diagonal matrix containing squares of residuals
of $\mathbf{Y}$. The
same equation may be used for ``standard'' analysis, in which case
the elements of the $\mathbf{R}$ matrix are constant, namely the mean
residual sum of squares (the estimate of $\sigma^2$).

Similar to that, the robust matrix is computed for logistic regression with
$$
\var_r = (\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1} (\mathbf{X}^T\mathbf{R}\mathbf{X})
(\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1},
$$
where $\mathbf{1}$ is the vector of ones and $\mathbf{W}$ is the diagonal matrix
of ''weights'' used in logistic regression.


\subsubsection{Cox proportional hazards model}
The implementation of the Cox proportional hazard model used
\index{Cox proportional hazards model}
\index{proportional hazards model}
\index{regression!Cox proportional hazards}
in \PA{} is entirely based on the code of \texttt{R}
library \texttt{survival} developed by Thomas Lumley
(function \texttt{coxfit2}), and is therefore not described here.

Many thanks to Thomas for making his code available under GNU GPL!

\subsection{Analysis of pedigree data}
The framework for analysis of pedigree data follows the two-step logic
developed in the works of Aulchenko \emph{et al}. (2007) and Chen and
Abecasis (2007). The general analysis model is a linear mixed model
where the expectation of the trait is defined as
$$
E[\mathbf{Y}] = \mathbf{X} \mathbf{\beta},
$$
identical to that defined for a linear model
(cf.~Eq.~\ref{eq:expectation}). To account for correlations between
the phenotypes of relatives which may be induced by family relations
the variance-covariance matrix is defined to be proportional to the
linear combination of the identity matrix $\mathbf{I}$ and the
relationship matrix $\mathbf{\Phi}$:
$$
\mathbf{V}_{\sigma^2,h^2} = \sigma^2 \left( 2 h^2 \mathbf{\Phi} + (1-h^2)
\mathbf{I} \right),
$$
where $h^2$ is the heritability of the trait. The relationship matrix
$\mathbf{\Phi}$ is twice the matrix containing the coefficients of
kinship between all pairs of individuals under consideration; its
estimation is discussed separately in section \ref{kinship}.

Estimation of a model defined in such a way is possible by numerical
maximization of the likelihood function, however, the estimation of
this model for large pedigrees is laborious, and is not
computationally feasible for hundreds of thousands to millions of SNPs
to be tested in the context of GWAS, as we have demonstrated
previously (Aulchenko \emph{et al}., 2007).

\subsubsection{Two-step score test for association}
A two-step score test approach is therefore used to decrease the
computational burden. Let us first re-define the expectation of the
trait by splitting the design matrix in two parts, the ``base'' part
$\mathbf{X}_x$, which includes all terms not changing across all SNP
models fit in GWAS (e.g.\ effects of sex, age, etc.), and the part
including SNP information, $\mathbf{X_g}$:
$$
E[\mathbf{Y}] = \mathbf{X}_x \mathbf{\beta}_x +
\mathbf{X}_g \mathbf{\beta}_g.
$$
Note that the latter design matrix may include not only the main SNP
effect, but e.g.\ SNP by environment interaction terms.

In the first step, a linear mixed model not including SNP effects
$$
E[\mathbf{Y}] = \mathbf{X}_x \mathbf{\beta}_x
$$
is fitted. The maximum likelihood estimates (MLEs) of the model
parameters (regression coefficients for the fixed effects
$\hat{\mathbf{\beta}}_x$, the residual variance $\hat{\sigma}^2_x$ and
the heritability $\hat{h}^2_x$) can be obtained by numerical
maximization of the likelihood function
$$
\mathrm{logLik}(\beta_x,h^2,\sigma^2) = -\frac{1}{2} \left(
  \log_e|\mathbf{V}_{\sigma^2,h^2}| + (\mathbf{Y} - \beta_x
  \mathbf{X}_x)^T \, \mathbf{V}_{\sigma^2,h^2}^{-1} \, (\mathbf{Y} -
  \beta_x \mathbf{X}_x) \right ),
$$
where $\mathbf{V}_{\sigma^2,h^2}^{-1}$ is the inverse and
$|\mathbf{V}_{\sigma^2,h^2}|$ is the determinant of the
variance-covariance matrix.

In the second step, the unbiased estimates of the fixed effects of the
terms involving SNP are obtained with
$$
\hat{\beta}_g = (\mathbf{X}^T_g
\mathbf{V}^{-1}_{\hat{\sigma}^2,\hat{h}^2}
\mathbf{X}_g)^{-1}
\mathbf{X}^T_g \mathbf{V}^{-1}_{\hat{\sigma}^2,\hat{h}^2}
\mathbf{R}_{\hat{\beta}_x},
$$
where $\mathbf{V}^{-1}_{\hat{\sigma}^2,\hat{h}^2}$ is the inverse
variance-covariance matrix at the point of the MLE estimates of
$\hat{h}^2_x$ and $\hat{\sigma}^2_x$, and $\mathbf{R}_{\hat{\beta}_x}
= \mathbf{Y} - \hat{\beta}_x \mathbf{X}_x$ is the vector of residuals
obtained from the base regression model. Under the null model, the
inverse variance-covariance matrix of the parameter's estimates is
defined as
$$
\var_{\hat{\beta}_g} = \hat{\sigma}^2_x (\mathbf{X}^T_g \mathbf{V}^{-1}_{\hat{\sigma}^2,\hat{h}^2} \mathbf{X}_g)^{-1}.
$$
Thus the score test for joint significance of the terms involving SNP
can be obtained with
$$
T^2 = (\hat{\beta}_g - \beta_{g,0})^T \,
\var_{\hat{\beta}_g}^{-1} (\hat{\beta}_g - \beta_{g,0}),
$$
where $\beta_{g,0}$ are the values of parameters fixed under the null
model.  Under the null hypothesis this test statistic asymptotically
follows the $\chi^2$ distribution with the number of degrees of
freedom equal to the number of parameters tested.  The significance of
an individual $j$-th element of the vector $\hat{\beta}_g$ can be
tested with
$$
T^2_j = \hat{\beta}_{g}^2(j) \var_{\hat{\beta}_g}^{-1}(jj),
$$
where $\hat{\beta}_{g}^2(j)$ is the square of the $j$-th element of
the vector of estimates $\hat{\beta}_{g}$, and
$\var_{\hat{\beta}_g}^{-1}(jj)$ corresponds to the $j$-th diagonal
element of $\var_{\hat{\beta}_g}^{-1}$. The latter statistic
asymptotically follows $\chi^2_1$.

\subsubsection{Estimation of the kinship matrix}
\label{kinship}

The relationship matrix $\mathbf{\Phi}$ used in estimation of the
linear mixed model for pedigree data is twice the matrix containing
the coefficients of kinship between all pairs of individuals under consideration.
This coefficient is defined as the probability that two gametes randomly sampled
from each member of the pair are identical-by-descent (IBD), that is they are copies
of exactly the same ancestral allele. The expectation of kinship
can be estimated from pedigree data using standard methods, for example the
kinship for two outbred sibs is $1/4$, for grandchild-grandparent is $1/8$, etc.
For an outbred person, the kinship coefficient is $1/2$ -- that is two gametes
sampled from this person at random are IBD only if the same gamete is
sampled. However, if the person is inbred, there is a chance that a maternal
and paternal chromosomes are also IBD. The probability of this is characterized
by kinship between individual's parents, which is defined as the individual's
inbreeding coefficient, $F$. In this case, the kinship coefficient for the
individual is $F + 1/2$. Similar logic applies to computation of the kinship
coefficient for other types of pairs in inbred pedigrees.

The kinship matrix can be computed using the pedigree data using standard methods.
However, in many cases, pedigree information may be absent, incomplete, or not
reliable. Moreover, the estimates obtained using pedigree data reflect the
expectation of the kinship, while the true realization of kinship may vary
around this expectation. In presence of genomic data it may therefore be
desirable to estimate the kinship coefficient from these, and not from pedigree.
It can be demonstrated that unbiased and positive semi-definite estimator
of the kinship matrix can be obtained (Astle and Balding, 2010; Amin \emph{et al.}, 2007)
by computing the kinship coefficients between individuals $i$ and $j$ with
$$
\hat{K}_{ij} = \frac{1}{L} \sum_{l=1}^L \frac{ (g_{l,i} - p_l) (g_{l,j} - p_l) }{ p_l (1-p_l) }
$$
where $L$ is the number of loci, $p_l$ is the allelic frequency at $l$-th locus
and $g_{l,j}$ is the genotype of $j$-th person at the $l$-th locus, coded
as $0$, $1/2$, and $1$, corresponding to the homozygous, heterozygous, and
other type of homozygous genotype. The frequency is computed for the allele
which, when homozygous, corresponds to the genotype coded as ``1''.


\section{How to cite}

If you used \PA{} for your analysis please give a link to the
\texttt{GenABEL project} home page
\begin{quote}
\url{http://www.genabel.org/}
\end{quote}
and cite the \PA{} paper to give us some credit:
\begin{quote}
Aulchenko YS, Struchalin MV, van Duijn CM.\\
\emph{ProbABEL package for genome-wide association analysis of imputed data.}\\
BMC Bioinformatics. 2010, 11:134.
\end{quote}
A proper reference may look like
\begin{quote}
For the analysis of imputed data, we used \PA{} v.0.4.3
from the \texttt{GenABEL} suite of programs (Aulchenko \emph{et al.}, 2010).
\end{quote}

If you have used the Cox proportional hazard model, please mention the
\texttt{R} package \texttt{survival} by Thomas Lumley. Additionally
to the above citation, please include
\begin{quote}
The Cox proportional hazards model implemented in \PA{}
makes use of the source code of the \texttt{R} package ''\texttt{survival}''
as implemented by T. Lumley.
\end{quote}


%\subsection{Generalized linear models}

%Let us define linear predictor $\eta$ as
%$$
%\eta = \mathbf{X} \beta
%$$
%For regular linear regression model described in previous section
%this linear predictor models the expectation of the trait of interest,
%while

%\begin{table}
%\begin{center}
%\caption{Link .. functions}
%\begin{tabular}{lccc}
%\hline
%Distribution & $\mu$ & $\tau^2$ & $v(\mu)$ \\
%\hline
%Normal  &            $\eta$ & $\sigma^2$ & \\
%Binomial & $1/(1+e^{-\eta})$ & & \\
%Poisson  &        $e^{\eta}$ & & \\
%\hline
%\end{tabular}
%\end{center}
%\end{table}

%Iteratively re-weighted least squares are used to obtain parameters'
%estimates.

\printindex

\end{document}
