% -*- mode: noweb; noweb-default-code-mode: R-mode; -*-
%\VignetteIndexEntry{hapFabia: Manual for the R package}
%\VignetteDepends{hapFabia}
%\VignettePackage{hapFabia}
%\VignetteKeywords{hapFabia,models,multivariate,cluster,hplot,file,genetics,haplotype,identity by descent,bicluster,next generation sequencing,genotype,single nucleotide polymorphism,single nucleotide variation,rare variants,rare SNPs, rare SNVs,rare haplotypes,short haplotypes}
%setwd("c:/sepp/work/hapfabia/hapFabia_10/hapFabia/vignettes")
%Sweave("hapFabia.Rnw")
%tools::texi2dvi("hapFabia.tex",pdf=TRUE)
%Stangle("hapFabia.Rnw")
%source("hapFabia.R")


\documentclass[article]{bioinf}

\usepackage{amsmath,amssymb}
\usepackage{bm}
\usepackage{natbib}
\usepackage{hyperref}
\usepackage{booktabs}
\usepackage{afterpage}
\usepackage{float}

\hypersetup{colorlinks=false,
   pdfborder=0 0 0,
   pdftitle={Identification of rare and short haplotype clusters in large sequencing data},
   pdfauthor={Sepp Hochreiter}}

\renewcommand{\dblfloatpagefraction}{0.99}
\renewcommand{\topfraction}{0.99}
\renewcommand{\bottomfraction}{0.99}
\renewcommand{\textfraction}{0.01}
\setcounter{dbltopnumber}{2}
\setcounter{bottomnumber}{2}


\def\Up{\bm{\Upsilon}}
\def\pp{\bm{\Psi}}
\def\epp{\bm{\epsilon}}
\def\Ncal{\mathcal{N}}
\def\X{\bm{X}}
\def\x{\bm{x}}
\def\xii{\bm{\xi}}
\def\Xii{\bm{\Xi}}
\def\oo{\mathrm{old}}
\def\nn{\mathrm{new}}
\def\sign{\mathrm{sign}}

\def\La{\bm{\Lambda}}
\def\la{\bm{\lambda}}
\def\z{\bm{z}}
\def\Z{\bm{Z}}
\def\Rb{\mathbb{R}}
\def\E{\mathbf{\mathrm{E}}}

\newcommand{\R}{\textrm{R} }
\newcommand{\Rfunction}[1]{{\texttt{#1}}}
\newcommand{\Robject}[1]{{\texttt{#1}}}
\newcommand{\Rpackage}[1]{{\texttt{#1}}}




\title{hapFabia: Identification of rare and short haplotype clusters in large sequencing data\\ \textit{---
    Manual for the R package ---}}

\author{Sepp Hochreiter}
\affiliation{Institute of Bioinformatics, Johannes Kepler University
Linz\\Altenberger Str. 69, 4040 Linz, Austria\\
\email{hochreit@bioinf.jku.at}}

\usepackage[noae]{Sweave}


\SweaveOpts{eps=FALSE}

\begin{document}
<<echo=FALSE,results=hide>>=
options(width=60)
set.seed(0)
library(hapFabia)
hapFabiaVer<-packageDescription("hapFabia")$Version
@
%$

\newcommand{\hapFabiaVer}{\Sexpr{hapFabiaVer}}

\manualtitlepage[Version \hapFabiaVer, \today]

\newlength{\auxparskip}
\setlength{\auxparskip}{\parskip}
\setlength{\parskip}{0pt}
\tableofcontents
\clearpage
\setlength{\parskip}{\auxparskip}


\section{Introduction}

This package \Rpackage{hapFabia} provides software
for the method HapFABIA which identifies rare and short haplotype
clusters in
large sequencing data with a focus on rare variants.

Individuals that inherited a particular DNA segment
from the same founder constitute a haplotype cluster
by sharing minor alleles of variants that tag this segment.
Knowledge of haplotype clusters are relevant for
genetic and genomic studies and for population genetics where they
shed light on the evolutionary history of humans.

\begin{figure*}[htb!]
\includegraphics[width=\textwidth]{./figures/pedi}
\caption{The DNA segment marked in yellow decended from a founder to
  different individuals which constitute a haplotype
  cluster by possessing the DNA segment.\label{fig:pedi}}
\end{figure*}
%\end{figure}

Rare haplotype clusters are identified by biclustering that is
clustering individuals only on a subset of single
nucleotide variants (SNVs), the tagSNVs (these SNVs tag the haplotype cluster).
Biclustering combines linkage disequilibrium (LD) information across individuals via
correlations of SNVs and identity by descent (IBD) information along the chromosome
via contiguous identical alleles.
The LD information enters HapFABIA through
FABIA implemented in the package \Rpackage{fabia}
which detects correlations of multiple SNVs by a factor
analysis model. The IBD information enters HapFABIA through
identification of haplotype clusters as local agglomerations of correlated SNVs.



%\begin{figure}[htb!]
%\includegraphics[width=\columnwidth]{./figures/plotACv}
\begin{figure*}[htb!]
\includegraphics[width=\textwidth]{./figures/plotAC}
\caption{Different approaches to detect haplotype clusters. The $y$-axis
enumerates individuals and the $x$-axis SNVs. Yellow and violet bars
indicate major and minor alleles, respectively. Haplotype cluster
segements are marked in gold.
Panel {\bf A}: Linkage disequilibrium (LD) is computed for two SNVs
across all samples. Identity by descent (IBD) is detected for two
individuals for a long haplotype region.
Panel {\bf B}: A haplotype block is depicted as a region of
increased LD comprising different haplotypes.
Biclustering clusters individuals together that have the same founder DNA segment
indicated by identical minor alleles of tagSNVs. \label{fig:bicluster}}
\end{figure*}
%\end{figure}

Fig.~\ref{fig:bicluster} shows in panel A
linkage disequilibrium (LD) and identity by descent (IBD)
and in panel B a haplotype block and a bicluster.
Individuals that belong to a haplotype cluster are similar to
each other because they share the same minor alleles in an inherited
founder DNA segment.
This similarity define a bicluster according to
the FABIA model \cite{Hochreiter:10}.
Consequently, a haplotype cluster constitutes a bicluster to which SNVs
belong if they tag the haplotype cluster and to which individuals belong if
they possess the founder DNA segment (see the bicluster in panel B of
Fig.~\ref{fig:bicluster}).

The multiplicative biclusters of FABIA are able to represent
homozygous regions that is two occurrences of a
haplotype in one diploid individual via their factors.
Overlapping haplotype clusters in one diploid individual
are represented through the additivity of biclusters in the FABIA model.
Examples of haplotype clusters found by hapFabia are given in
Fig.~\ref{fig:example1} and Fig.~\ref{fig:example2}.

\begin{figure*}[p]
\includegraphics[width=0.9\textwidth]{./figures/example1}
\caption{Example of a haplotype cluster. The $y$-axis
gives the chromosomes and the $x$-axis consecutive SNVs/Indels/SVs. Yellow
indicates major alleles, violet minor alleles of tagSNVs, and blue
minor alleles of other SNVs.
``model L'' indicates tagSNVs
identified by hapFabia in violet.   \label{fig:example1}}
\end{figure*}

\begin{figure*}[p]
\includegraphics[width=0.9\textwidth]{./figures/example8}
\caption{Another example of a haplotype cluster. See Fig.~\ref{fig:example1}
 for a description.\label{fig:example2}}
\end{figure*}

\clearpage

\section{Getting Started}

\subsection{Typical Analysis Pipeline}

First, we briefly describe a typical analysis pipeline.
Assume we have the genotype data of chromosome 1 in the file \texttt{filename.vcf.gz} in
compressed \texttt{vcf} format.
To prepare the data for \Rpackage{hapFabia} we have to perform
preprocessing steps.
First \texttt{filename.vcf.gz}
must be 1.\ uncompressed, then 2.\ converted to the sparse matrix
format, and then 3.\ split into segments. The following command line
commands perform these steps:
\begin{enumerate}
\item \texttt{gunzip filename.vcf.gz}
\item \texttt{./vcftoFABIA filename ./}
\item \texttt{./split\_sparse\_matrix filename \_mat.txt 10000 5000 1}
\end{enumerate}
The command line tools for step 2.\ and 3.\ are provided by the package
\Rpackage{hapFabia} in \texttt{inst/commandline/arch/}.
However step 2.\ and 3.\ can be performed in \R as well (see below).
The commandline parameters for \texttt{vcftoFABIA} are
1) filename without \texttt{.vcf} and 2) path to the file.
The commandline parameters for \texttt{split\_sparse\_matrix} are
1) filename without \texttt{.vcf} 2) extension (default\texttt{\_mat.txt})
3) segment length 4) shift size 5) indicator whether annotation is
present (is generated by  \texttt{vcftoFABIA} as default).
The data is split into segments of 10,000 SNVs where the distance
between adjacent segments is 5,000 thus they overlap by 5,000 SNVs.

After providing the file \texttt{filename.vcf} the following steps
constitute a typical analysis pipeline in \R:
\begin{Sinput}
R> #####define segments, overlap, filename #######
R> shiftSize <- 5000
R> segmentSize <- 10000
R> fileName="filename" # without type
R>
R> #####load libraries#######
R> library(hapFabia)
R> library(fabia)
R>
R> #####convert from .vcf to _mat.txt: step 2. above#######
R> vcftoFABIA(fileName=fileName)
R>
R> #####split/ generate segments: step 3. above#######
R> split_sparse_matrix(fileName=fileName,segmentSize=segmentSize,
+  shiftSize=shiftSize,annotation=TRUE)
R>
R> #####compute how many segments we have#######
R> ina <- as.numeric(readLines(paste(fileName,"_mat.txt",sep=""),n=2))
R> snps <- ina[2]
R> over <- segmentSize%/%shiftSize
R> N1 <- snps%/%shiftSize
R> endRunA <- (N1-over+2)
R>
R> #####analyze each segment#######
R> #####may be done by parallel runs#######
R> iterateSegments(startRun=1,endRun=endRunA,shift=shiftSize,
+  segmentSize=segmentSize,fileName=fileName,samples=0,
+  upperBP=0.05,p=10,iter=40,alpha=0.03,cyc=50,IBDlength=50,
+  Lt = 0.1,Zt = 0.2,thresCount=1e-5,minSNPsFactor=3/4,
+  pMAF=0.035,haplotypes=TRUE)
R>
R> #####identify duplicates#######
R> identifyDuplicates(fileName=fileName,startRun=1,endRun=endRunA,
+  shift=shiftSize,segmentSize=segmentSize)
R>
R> #####analyze results; parallel#######
R> anaRes <- analyzeIBDs(fileName=fileName,startRun=1,endRun=endRunA,
+  shift=shiftSize,segmentSize=segmentSize)
R> print("Number haplotype clusters:")
R> print(anaRes$noIBD)
R> print("Haplotye cluster length:")
R> print(anaRes$avibdLengthS)
R> print("Number of individuals belonging to a haplotye cluster:")
R> print(anaRes$avnoSampS)
R> print("Number of tagSNVs of a haplotye cluster:")
R> print(anaRes$avnoSnpsS)
R> print("MAF of tagSNVs of a haplotype cluster:")
R> print(anaRes$avnoFreqS)
R> print("MAF within the group of tagSNVs of a haplotye cluster:")
R> print(anaRes$avnoGroupFreqS)
R> print("Number of changes between major and minor allele frequency:")
R> print(anaRes$avnosnpChangeS)
R> print("Number of tagSNVs per individual:")
R> print(anaRes$avnosnpsPerSampleS)
R> print("Number of individuals having a tagSNV:")
R> print(anaRes$avnosamplePerSnpS)
R>
R> #####load result for segment 50#######
R> posAll <- 50 # (50-1)*5000 = 245000: segment 245000 to 255000
R> start <- (posAll-1)*shiftSize
R> end <- start + segmentSize
R> pRange <- paste("_",format(start,scientific=FALSE),"_",
+  format(end,scientific=FALSE),sep="")
R> load(file=paste(fileName,pRange,"_resAnno",".Rda",sep=""))
R> IBD <- resHapFabia$mergedIBD # $
R>
R> #####plot result for segment 50#######
R> plotIBDs(IBD,filename=paste(fileName,pRange,"_mat",sep=""))
R>    ##attention: filename without type ".txt"
R>
R> #####plot specific haplotype clusters (the first) in segment 50#######
R> plotOneIBD(IBD,filename=paste(fileName,pRange,"_mat",sep=""),
+  ibdC=1)
R>    ##attention: filename without type ".txt"
\end{Sinput}
First the packages \Rpackage{hapFabia} and \Rpackage{fabia} are loaded.
Then \texttt{filename.vcf} is converted by \Rfunction{vcftoFABIA} to sparse matrix format
\texttt{filename\_mat.txt} where the annotation file
\texttt{filename\_annot.txt} and the sample labels \texttt{filename\_samples.txt}
are generated.
The function \Rfunction{split\_sparse\_matrix} splits
chromosome 1 into segments of length 10,000 SNVs with a distance of
5,000 SNVs between the segments.
This results in 640 segments because we had more than 3,200,000 and
less than 3,210,000 SNVs in the genotype data.
Using the function \Rfunction{iterateSegments} haplotype clusters are
identified in these segments and the results stored in a EXCEL like
\texttt{.csv} format and as \R data object.
The function \Rfunction{identifyDuplicates} marks and memorizes
duplicates of haplotype custers which occur because the segments overlap.
Next the function \Rfunction{analyzeIBDs} analyzes the results where
duplicates as marked in previous step are not considered.
Results are listed by \texttt{anaRes}.

The next example shows how to view all haplotype clusters of a segment, for
which we chose segment 50 which corresponds to chromosome 1 interval
from 245,000 to 255,000 ($(50-1)*5000 = 245000$).
Then we plot a specific haplotype cluster, in this case the first (\texttt{ibdC=1}),
which can also be used to store a \texttt{.pdf} or a \texttt{.fig} for
editing with Xfig.
Examples of this plot function are given in
Fig.~\ref{fig:example1} and Fig.~\ref{fig:example2}.


An \R source file  \Rfunction{pipeline.R} of above pipeline
can be created and executed as follows:
\begin{Sinput}
R> makePipelineFile("filename",shiftSize=5000,segmentSize=10000)
R>
R> source("pipeline.R")
\end{Sinput}
NOTE: sourcing may take a while for large datasets.


\subsection{Examples}

Next we show an example how to use  \Rpackage{hapFabia}.

First the packages are loaded and then the data \texttt{simu}.
The data \texttt{simu} is written into three files:
\texttt{dataSim1fabia\_samples.txt} (sample names),
\texttt{dataSim1fabia\_annot.txt} (SNV annotation information), and
\texttt{dataSim1fabia\_mat.txt} (the data in sparse matrix format).
These are files which are in the standard pipeline produced by
\texttt{vcftoFABIA} and by \texttt{split\_sparse\_matrix}.



<<createData>>=
library(hapFabia)
library(fabia)
data(simu)
namesL <- simu[[1]]
haploN <- simu[[2]]
snps <- simu[[3]]
annot <- simu[[4]]
alleleIimp <- simu[[5]]
write.table(namesL,file="dataSim1fabia_samples.txt",
    quote = FALSE,row.names = FALSE,col.names = FALSE)
write(as.integer(haploN),file="dataSim1fabia_annot.txt",
    ncolumns=100)
write(as.integer(snps),file="dataSim1fabia_annot.txt",
    append=TRUE,ncolumns=100)
write.table(annot,file="dataSim1fabia_annot.txt", sep = " ",
    quote = FALSE,row.names = FALSE,col.names = FALSE,append=TRUE)
write(as.integer(haploN),file="dataSim1fabia_mat.txt",ncolumns=100)
write(as.integer(snps),file="dataSim1fabia_mat.txt",
    append=TRUE,ncolumns=100)

for (i in 1:haploN) {

  a1 <- which(alleleIimp[i,]>0.01)

  al <- length(a1)
  b1 <- alleleIimp[i,a1]

  a1 <- a1 - 1
  dim(a1) <- c(1,al)
  b1 <- format(as.double(b1),nsmall=1)
  dim(b1) <- c(1,al)

  write.table(al,file="dataSim1fabia_mat.txt", sep = " ",
     quote = FALSE,row.names = FALSE,col.names = FALSE,append=TRUE)
  write.table(a1,file="dataSim1fabia_mat.txt", sep = " ",
     quote = FALSE,row.names = FALSE,col.names = FALSE,append=TRUE)
  write.table(b1,file="dataSim1fabia_mat.txt", sep = " ",
     quote = FALSE,row.names = FALSE,col.names = FALSE,append=TRUE)

}
@

Now the haplotype clusters can be extracted from the data:
<<callHapFabia,results=hide>>=
hapRes <- hapFabia(fileName="dataSim1fabia",prefixPath="",
     upperBP=0.15,p=10,iter=1,IBDlength=10,pMAF=0.1)
@

To view the results the first haplotype cluster is plotted:
<<assignHaplotypeList>>=
mergedIBD <- hapRes$mergedIBD # $
@
<<plotHaplotype,fig=TRUE>>=
plotOneIBD(mergedIBD,filename="dataSim1fabia_mat",ibdC=1)
@

Finally we remove the temporary data files:
<<removeData>>=
unlink("dataSim1fabia_mat.txt")
unlink("dataSim1fabia_annot.txt")
unlink("dataSim1fabia_samples.txt")
unlink("dataSim1fabia_resIBD.Rda")
unlink("dataSim1fabia_mat_res_L.txt")
unlink("dataSim1fabia_mat_res_lapla.txt")
unlink("dataSim1fabia_mat_res_Psi.txt")
unlink("dataSim1fabia_mat_res_Z.txt")
unlink("dataSim1fabia_mat_res_Z_full.txt")
@


Here another example with random data:
<<secondExampleData>>=
simulateForFabia(minruns=2,maxruns=2)
@
<<secondExampleCallHapFabia,results=hide>>=
hapRes <- hapFabia(fileName="dataSim2fabia",prefixPath="",
      upperBP=0.15,p=10,iter=1,IBDlength=10,pMAF=0.1)
@
<<secondExampleAssignHaplotype>>=
mergedIBD <- hapRes$mergedIBD # $
@
<<secondExamplePlot,fig=TRUE>>=
plotOneIBD(mergedIBD,filename="dataSim2fabia_mat",ibdC=1)
@
<<secondExampleRemove>>=
unlink("dataSim2fabia_mat.txt")
unlink("dataSim2fabia_annot.txt")
unlink("dataSim2fabia_samples.txt")
unlink("dataSim2fabia_resIBD.Rda")
unlink("dataSim2fabia_mat_res_L.txt")
unlink("dataSim2fabia_mat_res_lapla.txt")
unlink("dataSim2fabia_mat_res_Psi.txt")
unlink("dataSim2fabia_mat_res_Z.txt")
unlink("dataSim2fabia_mat_res_Z_full.txt")
@



\section{hapFabia Method}
\label{sec:methods}

We introduce our novel method HapFABIA for extracting rare and short
haplotype clusters from large sequencing data.
HapFABIA first applies FABIA biclustering to phased or unphased
genotype data. The first step, biclustering, extracts similarities between
individuals based on a subset of SNVs but does not
consider that haplotype clusters consist of contiguous
nucleotides.
In the second step, HapFABIA therefore extracts haplotype clusters
from FABIA models by considering local tagSNV accumulations and pruning
spurious correlations of SNVs with a haplotype cluster.
These two HapFABIA steps are described in the next two subsections.


\subsection*{FABIA for genotype data}
\label{sec:fabia}

We propose identifying similarities between individuals by biclustering.
A bicluster corresponds to a haplotype cluster where
individuals are similar to each other by virtue of being identical at minor
alleles of tagSNVs.
Panel B in
Fig.~\ref{fig:bicluster} shows such a bicluster.


For biclustering we use the
``Factor Analysis for Bicluster Acquisition'' (FABIA)
biclustering model \cite{Hochreiter:10}.
In contrast to other biclustering methods such as BIMAX
\cite{Prelic:06} and QUBIC \cite{Li:09qubic}, FABIA can represent
homozygous regions and overlapping haplotype clusters because of its
multiplicative bicluster model and a data model that is additive in
its biclusters (see below).
FABIA can be applied to discrete genotype data but also to real values
that give the minor allele
likelihood or the minor allele dosage.
We use FABIA not only because it is well suited for genotyping data, but also
because it outperformed other biclustering methods
in extensive comparisons on different data sets \cite{Hochreiter:10}.

\subsubsection*{FABIA models genotype data by haplotype clusters}


FABIA describes genotype data $\X$
by an outer product $\z \ \la^T$ of two vectors $\la$ and $\z$.
The vector $\la$ corresponds to a haplotype cluster that
contains ones for haplotype cluster tagSNVs and zeros otherwise.
Vector $\z$ indicates the number of segments of an
individual that belong to the haplotype cluster (multiploidy).
FABIA can represent a
homozygous region, that is, two occurrences of a
haplotype in one diploid individual, by $z_i=2$.
Fig.\ \ref{fig:vectors} visualizes this representation of a genotype
matrix by an outer product for one haplotype cluster.


A diploid
individual may participate at two
haplotype clusters at a particular locus.
In this case genotyping sums the occurrences
of minor alleles, which is reflected by the additive FABIA model.
If we assume genotyping errors,
the FABIA model for genotype data $\X$ is
\begin{align}
\label{eq:expgen}
\X \ &= \ \sum_{i=1}^{p} \z_i \ \la_i^T \ + \ \Up \ = \  \Z \ \La \ +
\ \Up \ ,
\end{align}
where $\X \in \Rb^{l \times n}$ is the genotyping data; $\Z \in \Rb^{l
  \times p}$ is the matrix that gives for each individual the number of segments that
belong to a haplotype cluster; $\La \in \Rb^{p
  \times n}$ is the haplotype cluster tagSNV matrix;
$\Up \in \Rb^{l \times n}$ is additive noise; $n$ is the number
of SNVs; $l$ is the number of individuals (or chromosomes for phased
genotypes); $p$ is the number of haplotype clusters;
$\la_i  \in \Rb^n$ is the $i$th haplotype cluster tagSNV vector (the $i$th
row of $\La$);
and $\z_i  \in \Rb^l$  gives for each individual the number of segments that
belong to the $i$-th haplotype cluster (the $i$th column of $\Z$).
The additive noise not only covers genotyping errors but also
genotypes which cannot be explained by haplotype clusters.
Such unexplained genotypes may arise from recently acquired SNVs or
haplotype cluster segments that were fragmented by recombination events.

The FABIA model in Eq.~\eqref{eq:expgen} allows a generative
interpretation by a factor analysis model with $p$ factors:
\begin{align}
\label{eq:factormodel}
\x \ = \ \sum_{i=1}^{p} \la_i \  z_i  \ + \ \epp \ = \  \La^{T} \ \z \
+ \ \epp \ ,
\end{align}
where $\x \in \Rb^n$
is an observed genotype (a row of $\X$ written as
column vector), the vector of factors $\z \in \Rb^p$ gives the number of
segments belonging to each haplotype cluster for genotype $\x$,
$z_i$ is the number of segments that
belong to the $i$-th haplotype cluster (a
component of both $\z_i$ and $\z$),
and $\epp \in \Rb^{n}$ is the additive noise
(a row of $\Up$ written as
column vector).



As illustrated in Fig.\ \ref{fig:vectors},
both the vector $\z_i$
and the vector $\la_i$ of haplotype cluster tagSNVs should be sparse.
Sparse $\z_i$ means that only few individuals belong to
the haplotype cluster, that is,
the haplotype cluster is rare.
Sparse $\la_i$ means that only few SNVs are tagSNVs, which implies
short haplotype clusters.
Sparse $\z_i$ can be achieved if all components $z_i$ are sparse,
that is, if  the vector of factors $\z$ is sparse.
Note, that $\z_i \in \Rb^l$ in Eq.~\eqref{eq:expgen} gives for each
of the $l$ individuals the number of segments that
belong to the $i$-th haplotype cluster, while
$\z \in \Rb^p$ gives the number of
segments belonging to each of the $p$ haplotype clusters for one genotype $\x$.
In contrast to standard factor analysis, FABIA's model selection is
tailored to sparse factors and sparse
loadings, which are essential for haplotype cluster detection.
Sparseness in the FABIA model is obtained by a
component-wise independent Laplace
distribution both for the prior on the
parameters $\la_i$ and the
distribution of the factors $\z$ \cite{Hyvarinen:99a}:
\begin{align} \label{eq:laplaceprior1}
p(\z) \ = \ \left( {\textstyle\frac{1}{\sqrt{2}}} \right)^p
\prod_{i=1}^{p}
    e^{-\ \sqrt{2} \ |z_i|} \\ \label{eq:laplaceprior2}
p(\la_i) \ = \ \left( {\textstyle\frac{1}{\sqrt{2}}} \right)^n
\prod_{k=1}^{n} e^{-\ \sqrt{2} \
|\lambda_{ki}|}
\end{align}


The Laplace distribution of the factors leads to the analytically
intractable likelihood:
\begin{align}  \label{eq:likelihood}
p( \x \mid \La, \pp ) \ = \ \int p( \x \mid \z, \La, \pp) \ p(
\z) \ d\z \ .
\end{align}
Therefore, FABIA model selection is performed by means of
variational expectation maximization, which is
a variational optimization in
the expectation maximization (EM) framework
to maximize the posterior of the
parameters \cite{Girolami:01,Palmer:06,Hochreiter:10,Clevert:11,Klambauer:12}.
The idea of the variational approach is to express the prior $p(\z)$ by
the maximum
\begin{align}
 p(\z) \ = \ \underset{\xii}{\max} \ p(\z \mid \xii)
\end{align}
over a model family $p(\z \mid \xii)$
that is parametrized
by the variational parameter $\xii$
or by scale mixtures
\begin{align}
 p(\z) \ = \ \int p(\z \mid \xii) \ d\mu(\xii) \ .
\end{align}
A Laplace distribution can be expressed exactly by the maximum of a
Gaussian family or by Gaussian scale mixtures
\cite{Girolami:01,Palmer:06}.
Therefore for each $\x$, the maximum $\hat \xii$ of the variational parameter
$\xii$ allows to represent the Laplacian prior by a Gaussian:
\begin{align}
 \hat \xii \ = \ \underset{\xii}{\arg\max} \ p(\xii \mid \x) \ .
\end{align}
The maximum $\hat \xii$ can be computed analytically
(see Eq.~\eqref{eq:updateVar} below)
because for each Gaussian
the likelihood Eq.~\eqref{eq:likelihood} can be computed analytically.



If we denote the $j$th genotype by $\x_j \in \Rb^n$
with corresponding factors $\z_j \in \Rb^p$, then
we obtain the following variational E-step \cite{Hochreiter:10}:
\begin{align} \label{eq:updateE1}
\E \big(\z_j   \mid \x_j \big) \ = \  \big( \La^T \
  \pp^{-1} \ \La \
  + \  \Xii_j^{-1} \big)^{-1} \ \La^T \  \pp^{-1} \ \x_j \ , \\ \label{eq:updateE2}
\E \big(\z_j \ \z_j^T  \mid  \x_j \big) \ = \  \big(
\La^T \
  \pp^{-1} \ \La \
  + \  \Xii_j^{-1} \big)^{-1} \ + \\ \nonumber
 \E \big(\z_j   \mid  \x_j \big)  \ \E (\z_j  \mid
\x_j)^T \ ,
\end{align}
where $\Xii_j$ means $\mathrm{diag} \left(\xii_j \right)$. The
update for the variational parameter $\xii_j$ is
\begin{align}
\label{eq:updateVar}
\xii_{j} \ = \  \mathrm{diag}\left( \sqrt{\E (\z_j \
\z_j^T \mid \x_j)} \right)\ .
\end{align}
The variational M-step is \cite{Hochreiter:10}
\begin{align} \label{eq:updateM1}
\La^{\nn} \ = \  \frac{\frac{1}{l} \ \sum_{j=1}^{l} \x_j \ \E
(\z_j \mid  \x_j)^T  \ - \ \frac{\alpha}{l} \  \pp \ \sign
(\La
  )}{\frac{1}{l} \sum_{j=1}^{l} \E (\z_j \ \z_j^T \mid \x_j
  )} \\ \label{eq:updateM2}
\mathrm{diag} \left(\pp^{\nn}  \right) \ = \ \pp^{\mathrm{EM}} \ +
\ \mathrm{diag} \Big(\frac{\alpha}{l} \  \pp \ \sign (\La
  ) ( \La^{\nn} )^T \Big)\ , \\ \label{eq:updateM3}
\pp^{\mathrm{EM}} \  = \ \mathrm{diag} \bigg( \frac{1}{l}
\sum_{j=1}^l \x_j  \ \x_j^T \ - \ \La^{\nn}  \frac{1}{l} \sum_{j=1}^l
 \E \left(  \z_j   \mid  \x_j \right) \ \x_{j}^{T} \bigg) .
\end{align}
The parameter $\alpha$ controls the degree of sparseness (an
expectation of how rare the haplotype clusters are) and can be
introduced as a parameter of
the Laplacian prior of the factors \cite{Hochreiter:10}.

Note, that the number of bicluster must not be determined a priori if
$p$ is chosen large enough. The sparseness constraint will remove
spurious biclusters by setting $\la$ to a zero vector. In this way
FABIA automatically determines the number of biclusters.


\subsubsection*{Adaptation of FABIA for haplotype cluster detection}

Since an entry in the genotype matrix $\X$
reports how often the minor
allele is present, FABIA must
explain occurrences of minor alleles through haplotype clusters.
Because both the counts of minor alleles and the  the occurrences of
segments in haplotype clusters are non-negative,
we modified FABIA to enforce non-negative
loadings $\la_i$ by projecting negative components of $\la_i$ to zero.
If both $\x_j$ and $\La$ are non-negative, the posterior mean
$\E \big(\z_j   \mid \x_j \big)$ of $\z_j$ is non-negative, too.
Therefore also $\La^{\nn}$ is non-negative, because the prior term in the
update rule is not allowed to change the sign of the loadings.
Therefore it is sufficient to initialize $\La$ by positive values
to enforce non-negative factors and loadings.
$\Z$ is estimated by the prior mean $\E \big(\z_j   \mid \x_j \big)$
and used to identify individuals belonging to a haplotype cluster.

Further we developed a sparse matrix algebra which represents only
non-zero values and indices for
FABIA biclustering in order to efficiently analyze large genotyping data.
In  Eq.~\eqref{eq:updateE1} to Eq.~\eqref{eq:updateM3} the values of
the genotype vectors $\x_j$ and the values of the loading vector
$\la_i$ and the loading matrix $\La$ are sparse. Therefore we
introduced a sparse matrix algebra for multiplying both two sparse
vectors and a sparse vector by a dense vector.

To further speed up the computation, we developed an iterative version
of FABIA. Each FABIA iteration detects $p$ biclusters.
These $p$ biclusters are removed from the genotype matrix $\X$ before
starting the next iteration.


The vectors $\la_i$ and $\z_i$ acquire a new
interpretation when detecting rare haplotype clusters.
Components of $\la_i$ correspond to
tagSNVs and indicate to what degree they belong to the $i$-th
haplotype cluster. These components in particular indicate how many bicluster
individuals possess the minor allele of the corresponding tagSNV.
Components of $\z_i$ correspond to individuals and
indicate both the number of segments that belong to the
haplotype cluster and to what degree the individual's genotype
matches the haplotype cluster.




\subsection*{Haplotype cluster extraction from FABIA models}
\label{sec:extract}

FABIA biclustering extracts correlated SNVs but
does not consider that a haplotype cluster consists of contiguous
nucleotides which lead to local accumulations of a haplotype cluster's tagSNVs.
This information can be used to separate
unwanted spurious correlations between SNVs from desired correlations between tagSNVs.
Further, haplotype clusters may overlap at ancient, preexisting SNVs
and therefore may be joined in a FABIA model.
Note that not all ancient SNVs are common and can be filtered out in a
preprocessing step.
FABIA may also join haplotype clusters if they contain the same
individuals.


FABIA biclustering does not regard the order of SNVs or individuals,
thus random shuffling of SNVs does not change its result.
Therefore randomly correlated SNVs that are found by FABIA would be
uniformly distributed along the chromosome.
However SNVs that are correlated because they are tagSNVs of a
haplotype cluster agglomerate locally as they originate from an
ancient segment.
Deviations from the null hypothesis of uniformly distributed SNVs can
be detected by a binomial test for the number of expected SNVs within an
interval if the SNV frequency is given.
A low $p$-value hints at local agglomerations of bicluster SNVs
stemming from a haplotype cluster.


We propose a four-step procedure to extract haplotype clusters from FABIA
models:
\begin{enumerate}
\item
identifying local agglomerations of correlated SNVs based on a
binomial test,
\item
disentangling haplotype clusters and re-assigning individuals or chromosomes
to haplotype clusters,
\item
pruning haplotype clusters from SNVs with spurious correlations based
on an exponential test,
\item
merging similar haplotype clusters, and joining the parts of large haplotype clusters
that were divided by the bins from the first step.
\end{enumerate}

{\bf Step 1:}
FABIA model selection is independent of the
order of both the individuals/chromosomes and the SNVs.
Therefore, spurious correlated SNVs are unlikely to
agglomerate at a DNA locus, whereas local SNV agglomerations
hint at a haplotype cluster.
To detect agglomerations, we compute histogram counts of the largest
values of the FABIA model parameters $\la_i$.
The threshold ``Lt''  (fixed to 0.1)
gives the percentage of largest $\la_i$ values that are
used for the histogram.
The HapFABIA parameter ``IBDlength''
gives the largest haplotype cluster (in kbp or centiMorgans) to be considered.
The histogram bin size in number of SNVs is computed from ``IBDlength''
using the average genomic distance between adjacent SNVs.
To account for haplotype clusters that exceed bin borders, the histogram is
computed a second time with bins shifted by half the bin size.

The histogram bins with more counts than expected
are assumed to contain haplotype clusters.
We rank bin counts that exceed the expected value by a binomial test in
order to select bins.
We need to compute how many SNVs are expected in a bin, that
is, the probability of observing $k$ or more bin counts.
Without errors in the factor analysis model,
the minor allele of each tagSNV is present in
each of the $t$ individuals belonging to the haplotype cluster.
Let $p$ be the probability of a random
minor allele match between $t$ individuals.
The probability of observing $k$ or more matches for $n$ SNVs in
a bin is given by one minus the binomial distribution $F(k;n,p)$:
\begin{align}
\label{eq:binomial}
1 \ - \ F(k-1;n,p) \ &= \ \Pr(K \geq k) \ = \ \sum_{i=k}^n \ {n\choose i} \
p^i \ (1-p)^{n-i}  \ .
\end{align}
If $q$ is the minor allele frequency (MAF) for one SNV,
the probability $p$ of observing the minor allele of this SNV in
all $t$ individuals is $p=q^t$.
We assume that all SNVs have the same MAF $q$ --- in the
experiments we used the average MAF.
For $b$ bins, the probability of observing $k$ or more counts in at
least one bin is
\begin{align}
\label{eq:counts}
b \ {l\choose t} \ \sum_{i=k}^n \ {n\choose i} \
q^{it} \ \left(1-q^t \right)^{n-i} \ ,
\end{align}
where $l$ is the number of individuals and  ${l\choose t}$ is the
number of possibilities to chose $t$ individuals from the $l$ individuals.
If the probability Eq.~\eqref{eq:counts} is below the threshold ``thresCount'',
the according bin is selected for haplotype cluster extraction.
If $k_{\textrm{min}}$ is the minimum $k$ for which
Eq.~\eqref{eq:counts} is below the threshold ``thresCount'', then all
bins with counts $k \geq k_{\textrm{min}}$ are selected.
In our experiments, we allow haplotype clusters of only two individuals
(IBD runs), and therefore set $t=2$.



If a bin is selected, SNVs and individuals must be assigned to it.
Note that bicluster memberships of FABIA biclusters cannot be used directly,
because they include all bins.
First, those SNVs that contributed to its count are assigned to the selected bin.
Then, individuals or chromosomes are assigned to the selected bin
if they possess a minor allele at one or more SNVs that have been assigned
to the bin.
Individuals are only chosen from the top $z$-values of the FABIA model
to ensure that individuals are indeed similar to each other.
The parameter ``Zt'' (fixed to 0.2)
gives the percentage of top $z$-scores that are considered.





{\bf Step 2:}
In this step, haplotype clusters in a selected bin are disentangled, considering only
SNVs and individuals that have been assigned to the bin.
A haplotype cluster is initialized by an IBD run between
two core individuals that are identical at $m$ or more minor alleles.
The number $m$ is computed as $m = \textrm{minSNPsFactor} \times
k_{\textrm{min}}$.
All individuals that are identical in least $m$ minor alleles
to one of the two core haplotype cluster individuals are
classified as belonging to the haplotype cluster.

The tagSNVs of a haplotype cluster are
SNVs that have their minor allele in at least 2
individuals of the haplotype cluster.
Alternatively, tagSNVs must have their minor
allele in at least a specific percentage of haplotype cluster individuals (but
in at least 2).

Step 2 is repeated after removing the current haplotype cluster by
deleting the cluster's tagSNVs
until no more core individuals are found.


{\bf Step3 :} This step prunes haplotype cluster borders of SNVs that have
spurious correlations to the haplotype cluster.
Spurious correlations may still be present in a bin leading to an
overestimation of the cluster length.
Such SNVs can be identified by deviations of their MAFs
from those of other tagSNVs.
However, this criterion is not reliable for rare SNVs.
Therefore, we identify SNVs with spurious correlations to the haplotype cluster
on the basis of unusually large distances to other tagSNVs.
The deviation from an expected distance
is quantified by means of an exponential distribution with the median distance
between tagSNVs as parameter.
SNVs with distances leading to $p$-values below 1e-3 are removed.
The two furthest upstream and the two furthest downstream tagSNVs
are tested for their distances to other
tagSNVs. If the second-furthest
up- or downstream tagSNV is removed, then
the furthest up- or downstream tagSNV is removed, too.

{\bf Step 4:}
Haplotype clusters which are very similar to each other are merged.
Further, large haplotype cluster that were divided by the bins into smaller
parts are reconstructed in this step.
Thus, Haplotype clusters larger than given by
the bin sizes can be detected.
In order to compute similarities, we
assess how many tagSNVs or individuals of the smaller haplotype cluster are
explained by the larger haplotype cluster. This criterion is expressed by the
``overlap coefficient''
\begin{align}
\label{eq:dmin}
O(A,B) \ &= \ {{|A \cap B|}\over{\min\{|A|,|B|\}}} \ .
\end{align}
Using the overlap coefficient for both tagSNVs and individuals,
we define a distance-like measure between
haplotype clusters $H_1$ and $H_2$ by
\begin{align}
\label{eq:dist}
D(H_1,H_2) \ &= \ 1 \ - \ O(S_{H_1},S_{H_2}) \ O(I_{H_1},I_{H_2}) \ ,
\end{align}
where $S_{H_i}$ and $I_{H_i}$ are the tagSNVs and individuals belonging
to haplotype cluster $H_i$, respectively.
Using measure $D$,
haplotype clusters within overlapping (shifted) or neighboring bins
are clustered by hierarchical clustering using complete linkage.
Haplotype clusters are merged if their segments
are clustered together below a cutting height of 0.8.



\section{Tools to Analyze \Rpackage{fabia} Results}
\label{sec:tools}

To analyze the \Rpackage{fabia} results we provide some
functions. This might be convenient if parameters are optimized
for a specific data set.

Accumulations of \Rpackage{fabia} loadings can be
given as histogram counts to see
locations of accumulations:
<<histograms>>=
data(res)
h1 <- histL(res,n=1,p=0.9,intervv=50,off=0)
print(h1$counts)
h1 <- histL(res,n=1,w=0.5,intervv=50,off=0)
print(h1$counts)
@


\Rpackage{fabia} loadings can be plotted to identify locations of
accumulations:
<<plotLhistP,fig=TRUE>>=
data(res)
plotL(res,1,p=0.95,type="histogram",intervv=50)
@
<<plotLpointsP,fig=TRUE>>=
data(res)
plotL(res,1,p=0.95)
@
<<plotLpointsW,fig=TRUE>>=
data(res)
plotL(res,1,w=0.5)
@
<<plotLsmoothP,fig=TRUE>>=
data(res)
plotL(res,1,p=0.95,type="smooth")
@
<<plotLsmoothW,fig=TRUE>>=
data(res)
plotL(res,1,w=0.5,type="smooth")
@


Finally the largest \Rpackage{fabia} loadings $L$ and factors $Z$
can be listed. The largest values must exceed a threshold either
given by quantile $p$ or a value $w$:
<<topLZ>>=
data(res)

topLZ(res,1,LZ="L",p=0.95)
topLZ(res,1,LZ="L",w=0.95)

topLZ(res,1,LZ="Z",p=0.95)
topLZ(res,1,LZ="Z",w=0.4)

topLZ(res,1,LZ="L",indices=FALSE,p=0.95)
topLZ(res,1,LZ="L",indices=FALSE,w=0.95)

topLZ(res,1,LZ="Z",indices=FALSE,p=0.95)
topLZ(res,1,LZ="Z",indices=FALSE,w=0.4)
@



\bibliographystyle{natbib}

\bibliography{ibd}




\end{document}
