%% Eduardo Martinez Pedroza.
%%
%% Note from the author:  My dissertation document was created with this file. It passed the general formatting guidelines for PhD dissertations
%% of the Graduate College at the University of Oklahoma. I accept no responsibility in connection with the use of this file. Use it under your own risk.
%%

\documentclass[12pt]{report}

\usepackage[graphicx]{realboxes}
\usepackage{pdflscape}
\usepackage{graphicx}
\usepackage{rotating}
\usepackage{amsmath,amssymb}
\usepackage{color}
% \newtheorem{thm}{Theorem}%[section]
%\newtheorem{cor}[thm]{Corollary}
% \newtheorem{cor}{Corollary}%[section]%[chapter]
% \newtheorem{prop}{Proposition}%[section]%[chapter]
% \newtheorem{lem}{Lemma}%[section]%[chapter]
% \newtheorem{fact}{Fact}%[section]%[chapter]
% \theoremstyle{definition}
% \newtheorem{rem}{Remark}%[section]
% \newtheorem{property}{Property}
% \newtheorem{definition}{Definition}%[section]%[chapter]
% \newtheorem{exam}{Example}%[section]%[chapter]
% \numberwithin{equation}{section}
% \numberwithin{definition}{section}
% \numberwithin{lem}{section}
% \numberwithin{thm}{section}
% \numberwithin{rem}{section}
% \numberwithin{figure}{section}
\usepackage{algorithm}
\usepackage{algorithmic}
\usepackage{enumerate}
% \usepackage{slashbox}
\usepackage{dblfloatfix}
\usepackage{cite}
\usepackage{subfig}
\usepackage{multirow}

\renewcommand*\contentsname{Table of contents}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Standard AMS packages
\usepackage{amssymb, amsmath, amsthm}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  For mathematical diagrams
\usepackage{xypic}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% For inserting graphic files
\usepackage{graphicx} %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Hyper-references
%\usepackage[pdftex, plainpages=false]{hyperref}   %To be use with PDF-TEX

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Graduate college required margins
\usepackage[left=1.6in, top=1in, right=1in, bottom=1in]{geometry}
%\usepackage[left=1.7in, top=1.1in, right=1.1in, bottom=1.1in]{geometry}
%\usepackage[left=1.8in, top=1.2in, right=1.2in, bottom=1.2in]{geometry}
%\usepackage[left=1.9in, top=1.3in, right=1.3in, bottom=1.3in]{geometry}
%\usepackage[left=2in, top=1.5in, right=1.5in, bottom=1.5in]{geometry}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% line spacing
\usepackage{setspace}
    %\singlespacing
    %\onehalfspacing
    \doublespacing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% To set the chapter sizes
\usepackage{titlesec}
 
\titleformat
{\chapter} % command
[display] % shape
{\bfseries\large} % format
{\centering CHAPTER \ \thechapter} % label
{0.5ex} % sep
{
%     \rule{\textwidth}{1pt}
%     \vspace{1ex}
    \centering
} % before-code
[
%  \vspace{-0.5ex}%
% \rule{\textwidth}{0.3pt}
] % after-code
 
 
 \titleformat{\section}
 {\normalfont\large\bfseries}
 { \thesection}{.5cm}{}
%  
% \titlespacing{\section}{12pc}{1.5ex plus .1ex minus .2ex}{1pc}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Personal Environments and Commands
\newtheorem{theorem}{Theorem}[chapter]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}
\theoremstyle{definition}
% \newtheorem{definition}[theorem]{Definition}
\theoremstyle{remark}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{example}{Example}
\numberwithin{equation}{chapter}
\numberwithin{figure}{chapter}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% My commands
\newcommand{\reals}{\mathbb{R}}   %Real numbers' symbol
\newcommand{\hyperbolic}{\mathbb{H}}  %Hyperbolic symbol
\newcommand{\complex}{\mathbb{C}} %Complex numbers' symbol
\newcommand{\integers}{\mathbb{Z}}   %Integer numbers' symbol
\newcommand{\presentation}[2]{\langle \, #1  \, | \, #2 \,  \rangle}
\newcommand{\GammaH}{Cayley (G, X\cup \widetilde{\mathcal{H}}) }
\newcommand{\Hs}{\{ H_i \}_{i=1}^m}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% THESIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

\title{MEDICAL SIGNALS ALIGNMENT AND PRIVACY PROTECTION USING BELIEF PROPAGATION AND COMPRESSED SENSING}
\author{Aminmohammad Roozgard}
\date{November 2014}

\pagenumbering{roman}
{%% PAGE NUMBERING ROMAN

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% TITLE PAGE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\singlespacing

\newpage
\thispagestyle{empty}
\begin{center}
{ %\large
UNIVERSITY OF OKLAHOMA
\par
\vspace{0.16in}
GRADUATE COLLEGE
\par
\vspace{1.2in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
MEDICAL SIGNALS ALIGNMENT AND PRIVACY PROTECTION USING BELIEF PROPAGATION AND COMPRESSED SENSING
\par
% \vspace{0.17in}
% RELATIVELY HYPERBOLIC GROUPS
\par
\vspace{1.2in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A DISSERTATION
\par
\vspace{0.17in}
SUBMITTED TO THE GRADUATE FACULTY
\par
\vspace{0.17in}
in partial fulfillment of the requirements for the
\par
\vspace{0.17in}
Degree of
\par
\vspace{0.17in}
DOCTOR OF PHILOSOPHY
\par
\vfill
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
By
\par
\vspace{0.17in}
AMINMOHAMMAD ROOZGARD
\par
%\vspace{0.05in}
Norman, Oklahoma
\par
%\vspace{0.05in}
2014
}
\end{center}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SIGNATURE PAGE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\thispagestyle{empty}
\begin{center}
{%\large
MEDICAL SIGNALS ALIGNMENT AND PRIVACY PROTECTION USING BELIEF PROPAGATION AND COMPRESSED SENSING
\par
%\vspace{0.17in}
% RELATIVELY HYPERBOLIC GROUPS
% \par
\vspace{0.9in}
A DISSERTATION APPROVED FOR THE
\par
%\vspace{0.17in}
SCHOOL OF ELECTRICAL AND COMPUTER ENGINEERING
\par
\vspace{0.9in}
BY
\par
\vfill
\begin{flushright}
\begin{tabular}{cr}
\hline
\  \  \  \  \  \  \  & Dr. Samuel Cheng, Chair\\
 \\ &  \\
  \\ &  \\
\hline
  & Dr. Pramode Verma, Co-chair \\
\\ & \\
 \\ &  \\
\hline
 & Dr. James Sluss \\
\\   & \\
 \\ &  \\
\hline
 & Dr.  Hong Liu\\
 \\ &  \\
  \\ &  \\
\hline
   & Dr. William Ray \\
\end{tabular}
\end{flushright}
}
\end{center}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% COPYRIGHT PAGE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\thispagestyle{empty}
\   \
\par
\vfill
\begin{center}
\copyright \    Copyright by AMINMOHAMMAD ROOZGARD \   2014

All rights reserved.
\end{center}
}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ACKNOWLEDGEMENTS PAGE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter*{Acknowledgements}
I would never have been able to finish my dissertation without the guidance of my committee members and my wife. 

I would like to express the deepest appreciation to my advisor Dr. Samuel Cheng for his guidance, caring, patience, and providing me with an excellent atmosphere for doing research. I would like to thank Dr. Pramode Verma for his advice and financially supporting my research. His guidance has made this a thoughtful and rewarding journey. I would like to thank my committee members, Dr. Hong Liu, Dr. William Ray and Dr. James Sluss for their support to finish this dissertation in the last five years.

I would like to thank all the students and staff of the school of electrical and computer engineering, Tulsa campus, in particular, Mrs. Renee Wagenblatt for her support and help and patiently reading and correcting my writings.  

Words can not express my thanks to my parents and my brother. They always encouraged me to follow my dreams in my life, advised me to make the right decisions, supported me in my journey of learning and prayed for my success. I also would like to thank my family in law for their support and belief in me.

Finally, I would like to thank my lovely wife, Nafise. She worked with me, even harder than I did everyday and was always close to me to cheer me up and stand by me through every event of my life, both good ones and bad ones. 

\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% TABLE OF CONTENTS PAGE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

{%\singlespacing
  \thispagestyle{empty}
\tableofcontents
\listoftables
\listoffigures
}


\newpage

} %% END PAGE NUMBERING ROMAN

\chapter*{Abstract}
The advance in human genome sequencing technology has significantly reduced the cost of data generation and overwhelms the computing capability of sequence analysis. Efficiency, efficacy and scalability remain challenging in sequence alignment, which is an important and foundational operation for genome data analysis. 
In this dissemination, I propose a two stage approach to tackle this problem. In the preprocessing step, I match blocks of reference and target genome sequences based on the similarities between their empirical transition probability distributions using belief propagation. I then conduct a refined match using our recently published SCoBeP technique. I extract features from neighbors of an input nucleotide (a genome sequence of neighboring nucleotides that the input nucleotide is its middle nucleotide) and leverage sparse coding to find a set of candidate nucleotides, followed by using Belief Propagation (BP) to rank these candidates. 
Our experimental results demonstrated robustness in nucleotide sequence alignment and our results are competitive to those of the SOAP aligner %\cite{li2009soap2} 
and the BWA algorithm %\cite{li2009fast}
.


In addition, Most genomic datasets are not publicly accessible, due to privacy concerns. Patients’ genomic data contains identifiable markers and can be used to determine the presence of an individual in a dataset. Prior research shows that the re-identification can be possible when a very small set of genomic data is released. To protect patients, the data owners impose an application and evaluation procedure which often takes months to complete and limits the researchers. One solution to the problem is to let each data owner publish a set of pilot data to help data users choose the right datasets based on their needs. The data owners release these pilot data with the noise parameters and the mechanism that they used. A data user can run any kind of association tests and compare the outcomes with the other datasets outputs to get an idea
which datasets can be useful. I present a privacy preserving genomic data dissemination algorithm based on the compressed sensing. In my proposed method, I am adding the noise into the sparse representation of the input vector to make it differentially private. It means I find the sparse representation using using the SubSpace Pursuit and then disturb it with sufficient Laplasian noise. I compare my method with state-of-the-art compressed sensing privacy protection method.%  \cite{li2011compressive}.
%{\color{red} Moreover, compared to SCoBeP alignment, the proposed technique can handle sequences of much longer lengths.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 %\pagenumbering{roman}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INTRODUCTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{INTRODUCTION}
\pagenumbering{arabic}

Advances in DNA information extraction techniques have led to huge sequenced genomes from organisms spanning the tree of life. This increasing volume of genomic information requires Algorithms that can accurately compare multiple genome sequences to aid in the study of populations, pan-genomes, and genome evolution \cite{li2009soap2,li2009fast}. For a particular research, many individual genomes may
be sequenced to investigate genetic diversity. For example, the Cancer Genome Atlas \cite{hampton2006cancer} and 1000 Genomes Project \cite{siva20081000} will generate genome sequences from several thousand people. The complete bacterial genomes in public databases are already over one thousand. 

In this dissertation, I propose a novel nucleotide sequence Indexing and alignment method based on empirical transitional probability, sparse coding and belief propagation to compare the similarity of the nucleotide sequences. The alignment method of this dissemination is inspired by my recent works,
3D-SCoBeP described in chapter \ref{3Dmed}. Thanks to the sparse representation, my mechanism can handle long sequences with reduced memory footprint. I also leverage belief propagation to combine local and neighboring information of candidate nucleotides into consideration and generate matching scores to determine the
best match. First, I index the reference and the read genome sequence using empirical transitional probability and pick the top score indexes from the reference genome sequence to build an over--complete dictionary. I then find a set of candidate nucleotide for each nucleotide of the test sequence using sparse coding from the constructed dictionary. 

In addition, Most genomic datasets are not publicly accessible, due to privacy concerns. Prior research shows that the re-identification can be possible when a very small set of genomic data is released. To protect patients, the data owners release these pilot data with the noise parameters and the mechanism that they used. I present a privacy preserving genomic data dissemination algorithm based on the compressed sensing. In my proposed method, I am adding the noise into the sparse representation of the input vector to make it differentially private. In this chapter, I will briefly review four core techniques need in this dissemination: Sparse Coding(SC), Compressed Sensing(CS), Belief Propagation(BP) and Differential Privacy(DP).

\section{Belief Propagation (BP)}\label{BeliefPropagation}
Belief Propagation (BP) is an efficient inference method used on graphical models such as factor graphs \cite{kschischang2001factor}, Bayesian networks \cite{jensen1996introduction} and Markov random fields \cite{rozanov1982markov}. It was performed by passing messages through the factor graph of my problem. A factor graph is a graph that represent a function of multiple variables factors into a product of multiple functions with few variables. For example, a function $h(x_1,x_2,x_3,x_4)$ can be written as 
\begin{equation}
h(x_1,x_2,x_3,x_4)=f_a(x_2)f_b(x_1,x_2,x_3)f_c(x_2,x_4)f_d(x_4)
\end{equation}
where $f_a$, $f_b$, $f_c$ ,and $f_d$ are factor functions and the corresponding graph is shown in Fig. \ref{FactorGraphSample}.

\begin{figure}[t!]
\centering
\subfloat{\includegraphics[width=10cm]{Pics/SampleFG}}
\caption{{Corresponding factor graph of the equation (1). $x_1$, $x_2$, $x_3$, and $x_4$ are the variable nodes and $f_a$, $f_b$, $f_c$, and $f_d$ are the factor nodes.}}
\label{FactorGraphSample}
\end{figure}


Define $N(i)$ and $N(a)$ as two sets of neighbors of a variable node $i$ and a factor node $a$, respectively, and denote $m_{i\rightarrow a}$ and $m_{a\rightarrow i}$ as the forward and backward messages from node $i$ to node $a$, respectively. 
A message itself is a vector containing current beliefs of a node mapping to all candidate pixels in the reference image. For example, $m_{a \rightarrow i}(g_i)$ can be interpreted as the belief of node $a$ of how probable that the pixel of node $i$ in the test image should map to location $x_i$ in the reference image. 
Message updates for $m_{i\rightarrow a}$ and $m_{a\rightarrow i}$ are based on the messages received by the incoming messages towards nodes $i$ and $a$, respectively. More precisely, in my factor graph, the message update rules are given by  \cite{kschischang2001factor}
\begin{equation}
m_{i \rightarrow a }(x_i) = \prod_{b \in N(i) \backslash {a}}m_{b \rightarrow i}(x_i),
\label{eq:FG6}
\end{equation}
\begin{equation}
m_{a \rightarrow i }(x_i) =\sum_{x_a\backslash x_i } f(x_a) m_{j \rightarrow a }(x_j),
\label{eq:FG7}
\end{equation}
where $N(a)\backslash i$ means all neighbors of node $a$ excluding node $i$; the factor node $x_a$ is located between variable nodes $x_i$ and $x_j$. Also, I model $f(x_a)$ as follows:
\begin{equation}
f(x_a)=\tilde f(x_i,x_j)=e^{-\frac{||\mathcal{L}_i-\mathcal{L}_j||_2}{\sigma^2}}
\label{eq:FG10}
\end{equation}	
where $\sigma^2$ is a parameter to control the relative strength of the geometric constraint imposed by a neighboring node. If I increase the value of $\sigma^2$, the belief of each variable node will have less effect on its neighbors.

\section{Sparse Coding}\label{SparseCoding}

Consider a signal $y \in \mathbb{R}^M$ and a fat matrix $D \in \mathbb{R}^{M \times N}$, where I say the matrix is ``\textbf{fat}" since $M \ll N$. I am interested in representing $y$ with the column space of $D \in \mathbb{R}^{M \times N}$, i.e., finding $\alpha \in \mathbb{R}^N$  such that $y = D \alpha $. Since $D$ is fat, $\alpha$ is not unique. However, if I also restrict $\alpha$ to be the sparsest vector to satisfy $y = D \alpha$ (i.e., $\alpha$ that has fewest number of non-zero elements), then in theory there is a unique solution. Sparse coding precisely considers the aforementioned problem of finding a sparse $\alpha$ such that $y = D \alpha$ is satisfied.


Mathematically, I can write the problem as 

\begin{equation}
\hat{\alpha} = arg~min ||\alpha||_0  \; \; subject~~to \; \; y = D\alpha.
\label{eq:L0}
\end{equation}	

However, this $l^0$ optimization problem is NP-complete \cite{baraniuk2007compressive} and thus several alternative methods have been proposed to solve it \cite{pang2010robust}. For example, when a sufficiently sparse solution actually exists, substituting the $l^1$ norm for the $l^0$ pseudo-norm in (\ref{eq:L0})  as below

\begin{equation}
\hat{\alpha} = arg~min ||\alpha||_1  \; \; subject~~to \; \; y = D\alpha
\label{eq:L1}
\end{equation}	
will still result in the same solution \cite{baraniuk2007compressive} . Moreover, solving this modified problem is much easier since it can be readily transformed into a linear programming problem.
Besides linear programming, many other suboptimal techniques have been proposed to solve (\ref{eq:L1}), including orthogonal matching pursuit \cite{yang2010fast}, gradient projection \cite{maleh2007sparse} and subspace pursuit \cite{dai2009subspace}.



\section{Compressed Sensing}\label{CompressedSensing}

In this section, I briefly review the theory of the compressed sensing and its major processes and elements\footnote{for more information, please read \cite{candes2008restricted,barzigar2013video}}.  Consider an input vector $d \in \mathbb{R}^n$ that I want to represent it by a vector $x \in \mathbb{R}^n$ using an orthonormal basis (a transformation matrix) $\Psi \in \mathbb{R}^{n \times n}$  where $d=\Psi x \in \mathbb{R}^n$ and $x$ is a $s$-sparse vector ($s <n$) which means $x$ has at most  $s$ nonzero entries (see Figure \ref{FIG1}). 

\begin{figure}[t!]
\centering
\includegraphics[width=10cm]{Pics/Psi2x.png}
\caption{The sparse representation of the natural signals; $d$ is a one dimensional signal that can be represented in a transformation domain by a sparse vector $x$ where the transformation basis are the columns of the matrix $\Psi \in \mathbb{R}^{n \times n}$. Note that a vector is sparse if most of its elements are equal to zero. Here, the white squares are representative of the zero elements. }
\label{FIG1}
\end{figure}

Note that if $x$ is not a sparse vector, by zeroing the very small coefficients of $x$, I can make it sparse and  this new vector still keeps the most amount of information of original vector \cite{candes2008restricted,roozgard2013nucleotide,barzigar2013scobep}. Furthermore, the orthonormal basis $\Psi$ can be a standard transformation basis like wavelet basis or discrete cosine transform basis. The vector $d$ is the input of the compressed sensing method. 



The compressed sensing is divide into two processes: a ``\textbf{sampling process}'' and a ``\textbf{reconstruction process}''. The sampling process of the compressed sensing is a probabilistic compassion process using a random matrix $\Phi \in \mathbb{R}^{k \times n}$ that reduces the size of input from $n$ to $k=O\big(s \log (n/s)\big)$. This step is modeled as a linear mapping of the input vector $d$ into its random projection $y$ ($y=\Phi d \in \mathbb{R}^k$). The random matrix $\Phi$ is formed by independent and identically distributed (i.i.d.) entries from a symmetric Bernoulli distribution (see Figure \ref{FIG2}).

\begin{figure}[t!] 
\centering
\includegraphics[width=10cm]{Pics/Phi2d_2.png}
\caption{The sampling process of the compressed sensing method; $d$ is an input signal and a random matrix $\Phi \in \mathbb{R}^{k \times n}$ maps $d$ to a measurement vector $y \in \mathbb{R}^k$ which reduces the size of the input signal from $n$ to $k=O\big(s \log (n/s)\big)$.}
\label{FIG2}
\end{figure}

The reconstruction process of the compressed sensing exactly or approximately reconstructs the original data from the compressed samples. In this step, the vector $d$ is recovered from its random projection $y$ using the sparse representation $\hat{x}$. I consider a matrix $A=\Phi\Psi$ which both $\Phi$ (random matrix) and $\Psi$ (Orthogonal transformation matrix) are known from the sampling process. Using a $l_1$-Norm minimization method like Orthogonal Matching Pursuit (OMP)\cite{pati1993orthogonal}  or SubSpace Pursuit (SSP)\cite{dai2009subspace}, the recovered answer $\hat{x}$ is close enough to the original $x$ even in the presence of noise (see Figure \ref{FIG3}). 

Mathematically, the $l_1$-Norm minimizer  selects the smallest set of columns from $A$ such that 
\begin{equation}
\hat{x}=\arg \ \min\ \| {x}' \|_0 ,\mbox{ subject to } \|y-Ax'\|_2 < \epsilon  
\end{equation}
where $\|x\|_0:=|\{i:x_i \neq 0\}|$\footnote{ $\|x\|_1:=\sum_{i=1}^{n}{|x_i|}$ and $\|x\|_p:=\Big( {\sum_{i=1}^{n}x_i^p }\Big)^{1/p}$ where $p>1$  } and $\epsilon$ is a error threshold that determines how close is the $\hat{x}$ to the original vector $x$; The smaller $\epsilon$ forces the $\hat{x}$ to be closer to the $x$.  If I want to have a $s$-sparse $\hat{x}$, then the $l_1$-Norm minimization method solves the following problem:
\begin{equation}
\hat{x}=\arg \ \min_{ \| {x}' \|_0 \leq s } \|y-Ax'\|_2.
\end{equation}

\begin{figure}[t!]
\centering
\includegraphics[width=10cm]{Pics/A2x.png}
\caption{The reconstruction process of the compressed sensing method; $y$ is a measurement vector and $A=\Phi \Psi \in \mathbb{R}^{k \times n}$. The $l_1$-minimizer select the smallest set of $A$'s columns as the solution which is sparse and its error is less than the threshold $\epsilon$. Note that the columns marked by black are the selected ones.}
\label{FIG3}
\end{figure}



\section{SubSpace Pursuit}\label{SubSpacePursuit}

Subspace pursuit (SSP) is a $l_1$-Norm minimization method which has a reconstruction capability compared to the Linear Programming (LP) methods, and has very low reconstruction complexity of matching pursuit techniques for very sparse signals. For any sampling matrix $A$ satisfying the restricted isometry property (RIP) \cite{baraniuk2008simple} with a constant parameter independent of $K$, the Subspace pursuit algorithm can recover arbitrary K-sparse signals exactly from its noiseless measurements. 

When the measurements are inaccurate and/or the signal is not exactly sparse, the reconstruction distortion is of order a constant multiple of the measurement and/or signal perturbation energy. More precisely, for very sparse signals with $K = O( \sqrt{N})$ where $N$ the number of columns of $A$, which, for instance, the computational complexity of the Subspace pursuit algorithm is upper bounded by $O (mN K)$, but can be further reduced to $O (mN log K)$ when the nonzero entries of the sparse signal decay slowly \footnote{For more information about the details of the Subspace pursuit, please read \cite{dai2009subspace} }.

\section{Differential Privacy}

The current privacy protection techniques attempt to add noise to the allele frequencies of the case group in a way that the absent or present of any individual in the output result data be impossible.

A randomized algorithm $f$, called $\epsilon$--Differential Private, if for all adjacent datasets $\mathcal{D}$ and $\mathcal{D}^\prime$, and any possible output $\mathcal{D}$ in the output space of $f$: 
\begin{equation}
 \frac{Pr[f(\mathcal{D})=\hat{\mathcal{D}}]}{Pr[f(\mathcal{D}^\prime)=\hat{\mathcal{D}}]} \leq e ^ \epsilon.
\end{equation}

Note that $\hat{\mathcal{D}}$ is any dataset or a numerical value depends. The Laplacian mechanism \cite{dwork2006calibrating} is commonly used in data disturbing methods to achieve differential privacy, which adds noises generated from a Laplacian distribution, Laplace(0, $\Delta f/ \epsilon$), to the output of a computation on the dataset. The amount of noises will be calculated based on the sensitivity of the computed data. The sensitivity represents the maximum change of the output when a single modification happens to a dataset. Note that for any $f : \mathcal{D} \to  \mathbb{R}^d$, and all adjacent datasets $\mathcal{D}$ and $\mathcal{D}^\prime$, the sensitivity of $f$ can be calculated as follows: 
\begin{equation}
 \Delta f = \max_{\mathcal{D},\mathcal{D}^\prime} \| f(\mathcal{D})-f(\mathcal{D}) \|_1
\end{equation}

Another popular differential privacy mechanism is called Exponential mechanism \cite{mcsherry2007mechanism}. This mechanism will one output $t \in T$ that has optimum utility function and preserving the differential privacy. The inputs of the exponential mechanism are a data set $\mathcal{D}$, a range $T$, a privacy parameter $\epsilon$ , and a utility function $u(\cdot)$ where  $ u: (\mathcal{D} \times T) \to \mathbb{R}$. Then it that assigns a real value number to the output $t \in T$, where the higher value number shows the better utility. The mechanism induces a probability distribution over the range of $T$ and 
\begin{equation}
t \propto \exp\bigg(\frac{\epsilon u(\mathcal{D},t)}{2\Delta(u)}\bigg),
\end{equation}
where $\Delta u = \max_{\mathcal{D},\mathcal{D}^\prime} \| u(\mathcal{D})-u(\mathcal{D}) \|_1$ is the sensitivity of the utility function $u(\cdot)$.  


\chapter{3D MEDICAL IMAGE REGISTRATION USING SPARSE CODING AND BELIEF PROPAGATION}\label{3Dmed}

There are various medical imaging methods which have been used broadly in clinical and medical research. Consequently, the interests in registering and finding similarities of different images for diagnosis, treatment, and the sake of basic science are increasing. As images are typically captured at different times, angles, and often by different modalities, registering (or aligning) one image with another is challenging. In general, The accuracy of registration techniques will affect the performance and robustness of all subsequent analysis. I propose an efficient 3D medical image registration method based on sparse coding and belief propagation for Computed Tomography (CT) and Magnetic Resonance (MR) imaging. I used 3D image blocks as the input features and then I employed sparse coding with a dictionary of the features to find a set of the candidate voxels. To select optimum matches, belief propagation was subsequently applied on a factor graph of voxels generated by these candidate voxels. The outcome of belief propagation was interpreted as a probabilistic map of aligning the candidate voxels to the source voxels. I compared my proposed method (3D-SCoBeP) with the state-of-the-art medical image registration, MIRT \cite{myronenko2010intensity} and GP-Registration algorithm \cite{periaswamy2003elastic}. My objective results based on Root Mean Square Error (RMSE) are smaller than those from MIRT and GP-Registration. My results prove the effectiveness of my algorithm in registering the reference image to the source image.

\section{Introduction}\label{3D:introduction}


Image registration refers to the process of aligning two or more images obtained from different capturing modules and/or angles, and/or at different times into the same coordination system \cite{zitova2003image}. Registration is essential in many clinical applications including diagnosis \cite{roozgard2012Malignant}, simulating and surgical planning \cite{bro1996medical}. For example, registration techniques have been used to align a Magnetic Resonance (MR) image to a Computer Tomography (CT) image \cite{schnabel2010generic,hill1991registration}. In surgery, radiotherapy, or radiological intervention, preoperative medical data are used to diagnose, plan, simulate, guide, or otherwise assist a surgeon, or possibly a robot  \cite{lanfranco2004robotic}. While the surgical procedure is performed in the coordinate system relative to the patient, the surgical plan is constructed in the coordinate system relative to the preoperative data. The spatial transformation between the plan and the preoperative data  is 
formed by registration.
 Registration as a central step of processing images in the treatments, allows any voxel defined in the preoperative image to be precisely located in the patient coordinate system. This can aid the surgeon by delineating the position of the surgical instruments relative to the ultimate target.


% % % % % % 3D/2D Registration of CT and MR to X-Ray Images************************

Images of similar or different modalities need to be aligned for navigation, detection, data-fusion and visualization in medical applications \cite{maintz1998survey}. Medical image registration still presents many challenges. For example, finding a one-to-one correspondence between several scans of the patient is difficult, because the body of the patient can be subject to sudden changes or the modality of the scans can be different. The first one makes the transformations between scans highly non-rigid and the last one creates significantly different images in overall appearance and resolution\cite{periaswamy2006medical}. 

%  Several notable difficulties are  the nonlinear transformation between images and different appearance or resolution of images.

Many medical image registration methods have been developed in the last two decades \cite{myronenko2010intensity,periaswamy2003elastic,shum2000systems,glocker2008dense,zitova2003image,schnabel2010generic,maintz1998survey,roozgard2012Malignant,periaswamy2006medical,bro1996medical,hill1991registration,brown2007automatic,roche1999unifying,musse2001topology,cain2001projection,lanfranco2004robotic }. These can be divided into two major categories, namely, direct and feature-based matching. Direct methods use all available image data, and they result in very accurate registration if initialization points are close to target points at the beginning of the registration procedure \cite{shum2000systems}. For instance, in \cite{periaswamy2003elastic}, a general-purpose registration algorithm for the medical images has been developed which incorporates both geometric and intensity transformation. The authors modeled the transformation with a local affine model and a global smoothness constraint. Intensity variations 
are also modeled with local changes in brightness, contrast and a global smoothness constraint. Moreover, Myronenko and Song used the definition of the similarity measure to propose a registration method \cite{myronenko2010intensity}. They derived the similarity measure by analytically solving for the intensity correction field and its adaptive regularization. The final measure was interpreted as one that favors a registration with minimum compression complexity of the residual image between the two registered images. 

% feature based

Feature-based registration methods, utilize invariant features (specially around Harris corners) to ensure reliable matching. As a result, feature-based methods are independent from an initialization point \cite{brown2007automatic}. Also, feature-based registration methods obtain the transformation parameters from the set of extracted features. For example, Glocker \textit{et al.} \cite{glocker2008dense} used different levels of smoothness in modeling medical images and then used Markov Random Fields (MRFs) to formulate image deformations. Liu \textit{et al.} \cite{liu2002local,collignon19953d} and Elbakary \textit{et al.} \cite{elbakary2005accurate} registered multi-modal medical images using banks of local Gabor and Gaussian filters to evaluate the frequencies. The number and characteristics of filters in those works were selected empirically. Staring \textit{et al.} \cite{staring2009registration} incorporated multiple image features including the intensity gradients and Hessians. They combined parametric 
cubic B-splines, and an iterative stochastic gradient ascent optimization \cite{klein2007evaluation, klein2009adaptive} to solve the registration problem.



%% rigid & nonrigid

Image-registration techniques based on type of deformation have been divided into two categories: ``rigid'' and ``non-rigid''. In the rigid techniques, like \cite{roche1999unifying}, images are assumed to have rotation and translation only but in the non-rigid techniques, like \cite{likar2001hierarchical}, images can have restricted localized stretching. For example, in brain image registration with different modalities, a rigid body approximation is sufficient due to relatively little changes in brain shape over a short period between scans. In \cite{roche1999unifying}, authors formulated the rigid registration problem based on general image acquisition model and cast the problem of finding a similarity measure into their maximum likelihood problem. Then, they derived similarity measures for different modeling assumptions. Their experimental results concentrated on the multi-modal images of the brain. Sabuncu \textit{et al.}, in \cite{sabuncu2008using} introduced an entropy-based algorithm for registering 
rigid multi-modal images that incorporates spatial information. Spatial feature vectors obtained from the images and a minimum spanning-tree approach were used to estimate the conditional higher-dimensional entropy. They minimized the Jensen-Renyi divergence between the learned and new joint intensity distributions with a gradient descent method.
 
As an example for non-rigid registration techniques, Likar \textit{et al.} \cite{likar2001hierarchical} proposed a hierarchical image subdivision strategy to perform a non-rigid registration method based on mutual information. The non-rigid matching problem was decomposed into a Thin-Plate-Spline-based (TPS) elastic interpolation of multiple local rigid registrations of sub-images. One of the sub-categories of the non-rigid image registration is topology-preserving registration. In these kind of methods, the existing structures are kept, no new structures are allowed to be added, and neighborhood relationships between the structures are preserved. For example, Musse \textit{et al.} \cite{musse2001topology} proposed a parametric topology-preserving deformable image registration using the Gauss--Seidel optimization method. The Jacobian of the mapping was controlled over the domain of the transformation to ensure topology preservation. The authors derived the necessary and sufficient conditions for the 
determinant of the Jacobian of such transformations to be continuously positive everywhere and applied their method to the 2D images.

Based on the dimensions of input data, the registration method can be categorized as 3D or 2D registration techniques. The 3D registration method normally applies to the registration of two tomographic datasets but the 2D registration may apply to the separate slices from tomographic data. Also, a 3D to 2D registration may help to transfer the acquired 3D data to the 2D data, to facilitate treatment planning. In \cite{penney2001validation,hipwell2003intensity,byrne2004assessment}, the authors developed automated intensity-based algorithms for updating a 3D position of an interventional instrument using a single-plane angiogram registered to a 3D volume. In \cite{penney2001validation}, Penney \textit{et al.} aligned preoperative CT and intraoperative fluoroscopy images where the surface-target registration errors were of the order of $1–2$ $mm$. In \cite{hipwell2003intensity}, Hipwell \textit{et al.} expanded the former method to registering 3D cerebral Magnetic Resonance Angiography (MRA) with 2D X-Ray 
angiograms \cite{hipwell2003intensity} where their RMSE were $1.5 \pm 0.9$ $mm$ for 85\% of the clinical images. Byrne \textit{et al.}, in \cite{byrne2004assessment}, extended Penney \textit{et al.} work and registered 3D X-ray Digital Subtraction Angiography (3D-DSA) images. Their registration method accuracy was $1.3 \pm 0.6$ $mm$ in the clinical study of the two images with the same modality. 

Many researchers incorporate smoothness (or spatial coherence) conditions by reformulating matching into an optimization problem \cite{liu2011sift,tang2007non}. For example, Tang and Chung \cite{tang2007non} assigned a vector displacement label indicating the position in the test image to each pixel in the reference image. They used a smoothness constraint based on the first derivative to penalize sharp changes in displacement labels across pixels. Then they employed a graph-cuts method to solve that labeling problem. Moreover, Liu \textit{et al.} \cite{liu2011sift} used belief propagation to optimize cost function incorporated with smoothness constraints which encourage similar displacements of near-by pixels.



In this chapter, I propose a dense registration technique by aligning two CT or MR images using sparse coding and belief propagation. First, I build an overcomplete dictionary out of all 3D features of a reference image \cite{mairal2010online}. Note that since the dictionary is constructed by padding the features directly, I only need to normalize each columns. I then find a set of the candidate voxels for each voxel of the source image using sparse coding out of the constructed dictionary. The match score of each candidate voxel will be evaluated taking both local and neighboring information into account using belief propagation \cite{kschischang2001factor}. The best match will be selected as the candidate with the highest score. For those voxels with belief less than the threshold $\theta$, I use graph-cuts algorithm \cite{boykov2001fast} to find the proper matches. In comparing to the state-of-the-art belief propagation based registration methods, the key innovation of the proposed approach (3D-SCoBeP)  is the inclusion of a preprocessing step to preselect good candidate registration points for each voxel. Belief propagation is very powerful optimization technique, but if the size of the problem increases, it is more difficult to obtain a good local optimum. This restricts the size of the search range for each voxel. In prior approaches such as SIFT-flow \cite{liu2011sift}, the search range is simply chosen as a patch containing neighboring voxels around each target voxel. In contrast, a preprocessing step is used to carefully preselect candidate registration points for each voxel in the 3D-SCoBeP. Since these candidate points are selected from any voxel in the image, the search range of the 3D-SCoBeP is much larger than prior approaches and essentially covers the entire image. This is a main reason for the improvement of the 3D-SCoBeP over the prior works.

 {\color{black}{%Furthermore, the 3D-SCoBeP with this  preprocessing step could answer the question ``How  can I select match points?''.
A na\"{\i}ve approach computes the Mean Square Error (MSE) of the input patch with each possible patch of the reference image and selects patches that have the smallest MSEs. These kind of approaches have poor diversity which means the candidate patches are concentrated in a small region.
In this case, a small shift from the most similar patch generally does not decrease similarities sharply except for very high frequency patches. Consequently, this approach results in patches with very low diversity. 
Furthermore, the na\"{\i}ve approach may fail to find the true corresponding match points. It is possible that the na\"{\i}ve approach returns a set of candidates where all of them concentrate around a wrong point. Instead, I propose to find candidate match points using sparse coding. The intuition is that if these candidate patches are similar enough to the source patch, I should be able to construct a source patch out of good candidate patches (so they correspond to a sparse coding solution). With my technique, the sparse coding outputs the patches that can reconstruct the original patch through a linear combination. The resulting patches of sparse coding are likely to be complementary to each other and so provide a better diversity than the na\"{\i}ve solution.

}}

{\color{black}{The proposed method described here is inspired by my recent works, SCoBeP \cite{barzigar2012scobep} to answer the the current challenges in the medical data alignment which are:

1) the patient body movements while capturing the data, 

2) the patient body tissue changes due to progress of disease or treatment and 

3) different sampling rates of the data because of different sampling rates along the directions. 

A preliminary version of this work has been reported in \cite{roozgard2012medical}. Since then, much research has been done and the novel components  in addition to the aforementioned works are summarized as follows: 
\begin{itemize}
 \item Using 3D feature of input data which makes the proposed method more accurate than the other state-of-the-art methods;
 \item Using 3D feature of input data which makes the proposed method more accurate than the other state-of-the-art methods;
 \item Using a template  and interpolating the voxels of input data to map onto the template where the captured voxel coordination systems are different which increases the robustness of registration method;
 \item Employing a graph-cuts method as a post-processing method which further refines the matches obtained from the belief propagation step.
 
\end{itemize}
% 
% (i) 
% (ii)   
% %(iii)  
% %{\color{red}{because of different sampling rate and the patient movements and body changes}}
% (iii) 
}}


The rest of the chapter is organized as follows. In the next section, I will introduce the concept of my 3D-SCoBeP and the inference algorithm and then, in Section \ref{3D:EXPERIMENTAL}, I will show my simulation results, followed by a brief conclusion in Section \ref{nafiss:conclusion}.

% {\color{red}{The authors should briefly summarize belief propagation since this is in the title of the paper, but the authors only give a reference to what this is.}}

% \section{Background} \label{BackGround}


% \subsection{Belief Propagation and message passing}\label{backBP}
% 
% Belief Propagation (BP) \cite{kschischang2001factor} is an approximate inference method used on graphical models such as factor graphs. It is performed by passing messages through the factor graph of a given problem. 
% Define $N(i)$ and $N(a)$ as two sets of neighbors of a variable node ``$i$'' and a factor node ``$a$'', respectively, and denote $m_{i\rightarrow a}$ and $m_{a\rightarrow i}$ as the forward and backward messages from node ``$i$'' to node ``$a$''. The message itself is a vector containing current beliefs of a node mapping to all possible values. For example, $m_{a \rightarrow i}(x_i)$ can be interpreted as the belief of node ``$a$'' that how probable is the pixel of node $i$ to map to location $x_i$ in the factor graph. Message updates for $m_{i\rightarrow a}$ and $m_{a\rightarrow i}$ will be based on the messages received by the incoming messages towards nodes $i$ and $a$, respectively. More precisely, they are given by  \cite{kschischang2001factor}
% \begin{equation}
% m_{i \rightarrow a }(x_i) = \prod_{b \in N(i) \backslash {a}}m_{b \rightarrow i}(x_i),
% \label{eq:FG6}
% \end{equation}
% \begin{equation}
% m_{a \rightarrow i }(x_i) =\sum_{x_a\backslash x_i } f(x_a) \prod_{j \in N(a)\backslash {i}} m_{j \rightarrow a }(x_j),
% \label{eq:FG7}
% \end{equation}
% where I use $N(a)\backslash i$ to denote the neighbor of node $a$ excluding node $i$.

\section{3D--SCoBeP}\label{3D:Proposed_method}

As mentioned in Section \ref{3D:introduction}, in the medical image applications I need dense registration so that for each point of the \textit{source data} a corresponding match point will be found in the \textbf{reference data}. This section describes the implementation details of my proposed registration method for the 3D  medical data which is based on sparse coding and belief propagation. 


The proposed method described here is inspired by my recent works, SCoBeP \cite{barzigar2012scobep}.
First, I extract the features from the 3D reference data $\left \{ \mathcal{X}_s \right \}_{s=1}^{k_x} \in \mathbb{R}^{M \times N \times K}$ and the 3D source data $\left \{ \mathcal{Y}_s \right \}_{s=1}^{k_y} \in \mathbb{R}^{M \times N \times K}$ where $k_x$ and $k_y$ are the numbers of the reference and the source image slices, respectively. I focus on only using 3D block features even though the proposed approach can generally be applied to other features (such as SIFT-features \cite{lowe1999object} or Gabor-features \cite{feichtinger1998gabor}). The 3D feature block is a rectangular cube neighbor around each voxel of the 3D data which I reorder as a 1D vector. Thus, each feature considered here is essentially a vectorized 3D block centered around a voxel in a 3D data. 

Second, I create a dictionary $\mathcal{D}$ which contains all extracted feature vectors of the  reference data to match 
% {\color{red} {the extracted features of the reference CT image}} 
to the corresponding extracted features of the source data. The dictionary includes all vectorized 1D features as its columns where all of them have been normalized. I then  apply sparse coding to each extracted feature of the source data. Sparse coding will reconstruct a 3D source patch at voxel $[i,j,k]$ as a linear combination of the reference 3D patches. Denote $\alpha_{qlv}$ as the weight vector where each element corresponds to a coefficient in this combination. 
Note that the representation coefficients $\alpha_{ijk}$ should be sparse, i.e., it should be 0 for most coefficients.  To select the $n$ candidate voxels, I simply pick those corresponding to $n$ largest coefficients in the sparse coefficient vector. I denote a set as an $n \times 3$ matrix storing the locations of these candidate voxels and a probability vector $\rho$ as a length-$n$ vector storing the corresponding probabilities of the sparse coefficient vector. Each coefficient in the probability vector $\rho_{ijk}$ serves as a prior probability of matching the 3D source patch centered around voxel  $[i,j,k]$ to a 3D  patch of the reference data. This probability vector is taking only local characteristics into account but ignoring geometric characteristics of the matches. 

Finally, to incorporate these geometric characteristics, I model the problem by a factor graph and apply BP to identify the best matches similar to \cite{cheng2009improved}. I consider a 3D lattice factor graph as follows: For each voxel in the source 3D data, one variable node was assigned and then I connect each variable node to its six neighbors by a factor node (see Fig. \ref{3DFG}). Also, I consider one extra factor node for each variable node to take care of prior probabilities of the candidate points. In my model, the factor function $f(x_{i},x_{j})$ which can be interpreted as the local belief mapping nodes $i$ and $j$ to $x_i$ and $x_j$ can be used to impose the geometric constraint described earlier. Intuitively, since $x_i$ and $x_j$ are the corresponding mapped match points in the reference image of two neighboring voxels in the source image, I expect the probability of getting $x_i$ and $x_j$ to decrease as their distance apart increases. Therefore, I model the function of the factor node between two particular variable nodes $x_i$ and $x_j$ with a Gaussian kernal as \cite{kschischang2001factor}
\begin{equation}
f(x_i,x_j)=e^{-\frac{||x_i-x_j||_2}{\sigma^2}}
\label{eq:FG11}
\end{equation}	
where $\sigma^2$ is a parameter to control the relative strength of the geometric constraint imposed by a neighboring node. If I increase the value of $\sigma^2$, the belief of each variable node will have less effect on its neighbors. 


\begin{figure}
 \centering
 \subfloat{\includegraphics[width=14cm]{Pics/3DFG}}
 \caption{{Three dimensional factor graph of medical data used in Belief Propagation: for each voxel in the source 3D data, one variable node was assigned to incorporate these geometric characteristics. I connect each variable node to its six neighbors by a factor node and incorporate one extra factor node to store initial probabilities. A part of two slices of medical data corresponding factor graph is shown. This factor graph can be extended on X-axis, Y-axis and Z-axis.}}
 \label{3DFG} 
\end{figure}





To synthesize the source image, I replaced each voxel of the source image with the selected candidate voxel from the reference image where its probability is more than the threshold $\theta$. If the final maximum belief of the selected point was less than the threshold $\theta$, I employ graph-cuts algorithm  \cite{boykov2001fast,fulkerson2009class} to find the correspond voxel. {\color{black}{First, for voxels with maximum belief more than $\theta$, I calculate the movement of each voxel of reference data in comparison to the source data and create the displacement matrix $\beta$. Then, I feed the displacement matrix $\beta$ to graph-cuts algorithm to estimate the disparity of voxels with unsatisfactory beliefs. In the graph-cuts algorithm, I initialize the label of each voxel by its displacement value if its input belief is small then threshold $\theta$ and by average displacement value of all voxels otherwise. The data term  is defined by a quadratic function of the distance between the current 
label and the desired label, and the smoothness term is defined by a linear function of the distance between the current label and its neighborhood label. The neighborhoods are the same as which I used in the BP step and the swap algorithm of graph-cuts is applied to label voxel with beliefs less than the threshold $\theta$.}} 


\subsection{Implementation}
 This section describes the implementation's details of my proposed registration method. The main procedure for my proposed registration method is summarized in Algorithm \ref{nafiss:Denoisingalgorithm}.


\begin{algorithm} [t]
\caption{3D-SCoBeP for the medical image registration- estimate version of the registered image $\mathcal{Z}$ }
\begin{algorithmic} 
\renewcommand{\algorithmicrequire}{\textbf{Inputs:}}
% \renewcommand{\algorithmicensure}{\textbf{Output:}}
\REQUIRE {a reference data $\left \{ \mathcal{X}_s \right \}_{s=1}^{k_x} \in \mathbb{R}^{M \times N \times K}$, a source data $\left \{ \mathcal{Y}_s \right \}_{s=1}^{k_y} \in \mathbb{R}^{M \times N \times K}$, a threshold $\theta$
%, reference CT image slices number $k_x$, source CT image slices number $k_y$
, the number of the candidate points $n$} 


\renewcommand{\algorithmicrequire}{\textbf{Extract 3D dense feature and construct dictionary:}}
\REQUIRE {\ }

\begin{itemize}
\item $\textbf{Y} = ExtractDenseFeature \big(\left \{ \mathcal{Y}_s \right \}_{s=1}^{k_y}\big)$
\item $\textbf{X} = ExtractDenseFeature \big(\left \{ \mathcal{X}_s \right \}_{s=1}^{k_x}\big)$
\item $\mathcal{D} = MakeDic \big(\textbf{X}\big)$
\end{itemize}

\renewcommand{\algorithmicrequire}{\textbf{Find the initial estimate of the candidate voxels:}}
\REQUIRE {For each vector $y_{i,j,k}\in \textbf{Y}$ perform:}

\begin{itemize}%[(a)]
\item $\hat{\alpha}_{i,j,k} = FindSCV\big(\mathcal{D},y_{i,j,k}\big)$ 
\item $[\mathcal{L}_{i,j,k}, \rho_{i,j,k}]= FindTopSCV\big(n,\hat{\alpha}_{i,j,k}\big)$
\end{itemize}

\renewcommand{\algorithmicrequire}{\textbf{Refine the candidate voxels:}}
\REQUIRE {\ }
\begin{itemize}
\item $\hat{\rho}= BP\big(\mathcal{L},\rho\big)$
\end{itemize}


\renewcommand{\algorithmicrequire}{\textbf{Find  the correspond voxels:}}
\REQUIRE {\ }
\begin{itemize}
\item $[\mathcal{Z},\beta]= Warp\big(\textbf{X},\hat{\rho},\mathcal{L},\theta \big)$
\item if there is a voxel with probability less than $\theta$ then $\mathcal{Z}= Graphcuts \big( \textbf{X} ,\beta,\theta \big)$
\end{itemize}

\renewcommand{\algorithmicrequire}{\textbf{Output:}}
\REQUIRE {the estimated version of the registered image $\mathcal{Z}$}
\end{algorithmic}
\label{nafiss:Denoisingalgorithm}
\end{algorithm}
 
\clearpage
 

\textbf{Implementation Details:} \label{nafiss:Notation1}
\begin{itemize}

\item $\textbf{Y} = ExtractDenseFeature \big(\left \{ \mathcal{Y}_s \right \}_{s=1}^{k_y}\big)$ presents  a 3D block extractor algorithm using $\left \{ \mathcal{Y}_s \right \}_{s=1}^{k_y}$ as a source data. More precisely, I consider a 3D block of size $S = (2a + 1) \times (2b + 1) \times (2c + 1)$ containing neighboring voxels around each voxel on a 3D data, where $a$, $b$ and $c$ are positive integers. For each voxel $p_{i,j,k}$ in the source data $\left \{ \mathcal{Y}_s \right \}_{s=1}^{k_y}$, I vectorized the 3D block centered around the voxel $p_{i,j,k}$ to a feature vector $y_{i,j,k} \in \mathbb{R}^{S \times 1}$. A source feature $\textbf{Y} \in \mathbb{R}^{M \times N \times K \times S}$ is then constructed from $y_{i,j,k,t}$ as follows:

\begin{eqnarray}\nonumber
\textbf{Y} = \{ y_{i,j,k,t} \ |  &\ 1\leq i\leq M , &1\leq j \leq N ,\\ 
&1\leq k\leq K , &t \in S  \}. 
\label{eq11}
\end{eqnarray}
Note that $\textbf{X}$ is created in the same manner as $\textbf{Y}$ but instead from the reference data $\left \{ \mathcal{X}_s \right \}_{s=1}^{k_y}$. 

% % % % % % % % % % % % % % % % % % % % % 
\item  $\mathcal{D} = MakeDic \big(\textbf{X}\big)$ creates a dictionary $\mathcal{D}$ using the vectors of $\textbf{X}$. Later, the dictionary $\mathcal{D}$ is used to match the extracted features of the source data to the corresponding  extracted features of the reference data.
% , a dictionary which contains feature vectors of $\textbf{X}$ is constructed. Thus, 
I can write $\mathcal{D}$ as 
\begin{equation}
\mathcal{D} = \left [ x_{1,1,1}\ ... \ x_{1,1,K}\ x_{1,2,K}...\ x_{1,N,K} ... \ x_{M,N,K}\right ],
\label{eq2}
\end{equation}
where $x_{i,j,k}$ is a feature vector of  $\textbf{X}$. Note that I normalize dictionary $\mathcal{D}$ to guarantee the norm of each feature vector to be $1$.

% % % % % % % % % % % % % % % % % % % % % 


\item $\hat{\alpha}_{i,j,k} = FindSCV\big(\mathcal{D},y_{i,j,k}\big)$ finds the candidate match voxels using the sparse coding algorithm, where $\hat{\alpha}_{i,j,k}$ is a sparse vector.  Mathematically, I try to solve the following sparse coding problem to find the most sparse coefficient vector $\hat{\alpha}_{i,j,k}$ (see Fig. \ref{3D:SCpicture}) such that
\begin {equation}
 y_{i,j,k} = \mathcal{D} \hat{\alpha}_{i,j,k}.
\label{eq333}
\end {equation}



Although there are several methods to solve (\ref{eq333}) \cite{yang2010fast,dai2009subspace,maleh2007sparse},  in my work, I employ Subspace Pursuit (SP) \cite{dai2009subspace} because of its computational efficiency.

\begin{figure}[t]
\centering
\subfloat{\includegraphics[width=10cm]{Pics/yDalpha}}
\caption{Sparse representation of a feature vector $y_{i,j,k}$ with a dictionary $\mathcal{D}$ : $\hat{\alpha}_{i,j,k}$ as a sparse vector constructs the feature vector $y_{i,j,k}$ using a few columns (highlighted in gray) of dictionary $\mathcal{D}$.}
\label{3D:SCpicture}
\end{figure}

% % % % % % % % % % % % % % % % % % % % % % %  
\item $[\mathcal{L}_{i,j,k},\rho_{i,j,k}] = FindTopSCV\big(n,\hat{\alpha}_{i,j,k}\big)$ picks up the $n$ largest coefficients of $\hat{\alpha}_{i,j,k}$ as $n$ candidates. $\mathcal{L}_{i,j,k}$ as an $n \times 3$ matrix stores the locations of these candidate voxels and $\rho_{i,j,k}$ as a length-$n$ vector stores the corresponding values of $\mathcal{L}_{i,j,k}$. Each coefficient in $\rho_{i,j,k}$  serves as a prior probability of matching the source patch at $[i,j,k]$ to a patch centered around the voxel $x_{i,j,k}$. After finding the candidate locations $\mathcal{L}_{i,j,k}$ and their initial probabilities $\rho_{i,j,k}$ for each voxel, I concatenate the results and construct following matrices:
\begin{equation}
 \mathcal{L}=[\mathcal{L}_{1,1,1}\ ... \ \mathcal{L}_{1,1,K}\ \mathcal{L}_{1,2,K}...\ \mathcal{L}_{1,N,K} ... \ \mathcal{L}_{M,N,K}], \end{equation}
 \begin{equation}
 \rho=[\rho_{1,1,1}\ ... \ \rho_{1,1,K}\ \rho_{1,2,K}...\ \rho_{1,N,K} ... \ \rho_{M,N,K}]
\end{equation}
which I will use to apply belief propagation at the next step.


% % % % % % % % % % % % % % % % % % % % 

\item $\hat{\rho}= BP\big(\mathcal{L},\rho\big)$ models the problem by a factor graph and applies belief propagation \cite{kschischang2001factor} to update probability ${\rho}$. The updated probability $\hat{\rho}$  can be used to register the reference data onto the source data. In my case, I assign a variable node for each voxel in the source data and connect each pair of neighboring voxels with a factor node. Also, I introduce one extra factor node to take care of the prior knowledge obtained in the sparse coding step for each voxel of the source data (for more details, see \cite{barzigar2012scobep}).

\item $[\mathcal{Z},\beta]= Warp\big(\textbf{X},\hat{\rho},\mathcal{L},\theta \big)$ returns the registered image $\mathcal{Z}$ and a displacement matrix $\beta$ which contains the movement of each voxel of reference data. This matrix can be used to refine the result of the voxels with a probability less than $\theta$.

\item $\mathcal{Z}= Graphcuts \big( \textbf{X} ,\beta,\theta \big)$ {\color{black}{applies the graph-cuts method to find the displacement of the voxels with a probability less than $\theta$ and returns the registered image $\mathcal{Z}$. The displacement matrix $\beta$ keeps the movement of the voxels and marks the area with a probability less than $\theta$. I feed the matrix $\beta$ to graph-cuts algorithm and initialize the label of each voxel by the average displacement value of all voxels if it is marked in the matrix $\beta$ and by its displacement value otherwise. I use a quadratic function of the distance between the current label and the desired label as the data term and a linear function of the distance between the current label and the neighboring label as the smoothness term. The graph-cuts algorithm updates the matrix $\beta$ and I use the new displacement matrix to find the corresponding voxels.}}





% % % % % % % % % % % % % % % 

\item $\hat{\mathcal{Z}} = Warp\big(\left \{ \mathcal{X}_s \right \}_{s=1}^{k_y} \ \ \textbf{,} \ \ \hat{\rho} \ \ \textbf{,} \ \ \mathcal{L}     \big)$ displays the registered CT image $\hat{\mathcal{Z}}$ using the updated probabilities $\hat{\rho}$, the candidate voxels location  $\mathcal{L}$, and the reference CT image $\left \{ \mathcal{X}_s \right \}_{s=1}^{k_x}$. In my work, I select the most probable point after the BP step as the best match point. I assume that my registration method successfully finds a match for an input point if the most probable candidate has belief larger than a threshold $\theta$. Otherwise, I assume no best match is found. 





% % % % % % % % % % % % % % % % % % % % % 


\end{itemize}









\section{Experimental Results}\label{3D:EXPERIMENTAL}

 





The utility and novelty of my medical image registration algorithm lies in the fact that it can handle images captured not just from a single plane but also from different planes. Hence, in this section, I will study the performance of my method for both cases and compare it with the different registration methods.

In a brief statement, I will present two experiments in this section: \textit{the 3D CT image registration taken along a same direction} in Section \ref{ex:samedirection} where I consider the problem of registering two lung CT images of one person from two different times, and \textit{3D MR image registration taken along different directions} in Section \ref{3DMRSCOBEP} where two brain MR images were captured along the X--Z (sagittal) and the X--Y (transverse) planes. I implemented the 3D-SCoBeP algorithm in Matlab and tested it on a Pentium 3 GHz (11-GB RAM) machine. 

 
  

\subsection{3D CT Image Registration Taken In Same Directions}\label{ex:samedirection}

\begin{landscape}

\begin{figure}[t]
 \centering
% \subfloat[]{\includegraphics[width=10cm]{Pics/3D/3d1}}
% \subfloat[]{\includegraphics[width=10cm]{Pics/3D/3d2}}
\subfloat[]{\includegraphics[width=5cm]{Pics/3D/S_Lung_40angle}}
\subfloat[]{\includegraphics[width=5cm]{Pics/3D/R_Lung_40angle}}
\subfloat[]{\includegraphics[width=5cm]{Pics/3D/ScoBep_Lung_40angle}}

\subfloat[]{\includegraphics[width=5cm]{Pics/3D/3d_R_vs_S}}
\subfloat[]{\includegraphics[width=5cm]{Pics/3D/3d_R_vs_MIRT2}}
\subfloat[]{\includegraphics[width=5cm]{Pics/3D/3d_R_vs_SCo}}

% \subfloat[]{\includegraphics[width=5cm]{Pics/3D/3d_S_vs_MIRT}}
% \subfloat[]{\includegraphics[width=5cm]{Pics/3D/3d_R_vs_MIRT}}
% \subfloat[]{\includegraphics[width=5cm]{Pics/3D/3d_R_vs_MIRT1}}

\caption{Result of 3D-SCoBeP on Lung CT images. (a) The reference CT image; (b) The source CT image; (c) The 3D-SCoBeP result; (d)  The comparison between corespondent voxel between the source and the reference; (e) The comparison between corespondent voxel between the source and the MIRT result; (f) The comparison between corespondent voxel between the source and the 3D-SCoBeP result; In (d)-(e) I used a RGB image where the first channel of the image was assigned to the source image intensity and the second channel to the reference, MIRT, 3D-SCoBep results, respectively.} %Note that the third channel in (d), (e) and (f) are equal to 255. 
\label{lung_3d}
% \end{center}
\end{figure}

\end{landscape}

To evaluate the performance of my approach, I conducted tests on the data sets LIDC-IDRI \cite{website:LIDCIDRI} where the size of each slice of the CT images is $512 \times 512$ voxels. Throughout the experiments, the following parameters were used: the number of the candidate voxels $n$ is set to be $4$, $a = b = 3$ and $c=2$. To synthesize the source image, I replaced each voxel of the source CT image with the selected candidate voxel from the reference CT image. In other words, I map the reference CT image onto the source image using the updated probabilities and  the candidate voxels location. In my work, I select the most probable voxel after the BP step as the best match voxel. I assume that my registration method successfully finds a match for an input voxel if the most probable candidate has belief larger than a threshold $\theta = 0.3$. Otherwise, I assume no ``best match'' is found.
The threshold $\theta$ can be chosen empirically which was the way that I chose in Fig. \ref{nafiss:fig42-65} and \ref{nafiss:fig63-86} to express the results. 

\begin{figure}[htb]
 \centering
\subfloat[]{\includegraphics[width=3.5cm]{Pics/SourceImage.pdf}}\hspace{0pt}
\subfloat[]{\includegraphics[width=3.5cm]{Pics/ReferenceImage.pdf}}\hspace{0pt}

\subfloat[]{\includegraphics[width=3.5cm]{Pics/RegisterdMIRT.pdf}}\hspace{0pt}
\subfloat[]{\includegraphics[width=3.5cm]{Pics/RegisterdSource_qr.pdf}}\hspace{0pt}
\subfloat[]{\includegraphics[width=3.5cm]{Pics/SCoBeP-3D-42-65.pdf}}\hspace{0pt}

\subfloat[]{\includegraphics[width=3.5cm]{Pics/RefImage-65.pdf}}\hspace{0pt}
\subfloat[]{\includegraphics[width=3.5cm]{Pics/Source-42.pdf}}\hspace{0pt}

\subfloat[]{\includegraphics[width=3.5cm]{Pics/RegisterdSource_qr_ZoomIn.pdf}}\hspace{0pt}
\subfloat[]{\includegraphics[width=3.5cm]{Pics/RegisterdMIRT_ZoomIn.pdf}}\hspace{0pt}
\subfloat[]{\includegraphics[width=3.5cm]{Pics/SCoBeP-3D.pdf}}\hspace{0pt}
\caption{{Registration result of the lung CT images that were captured with six months gap. (a) Source image; (b) Reference image; (c) MIRT \cite{myronenko2010intensity} [RMSE: 24.31]; (d) GP-Registration \cite{periaswamy2003elastic} [RMSE: 28.78]; (e) 3D-SCoBeP [RMSE: 21.73]; (f) Source image (zoom in); (g) Reference image (zoom in); (h) MIRT \cite{myronenko2010intensity} (zoom in);  (i) GP-Registration \cite{periaswamy2003elastic} (zoom in); (j) 3D-SCoBeP (zoom in).}}
\label{nafiss:fig42-65}
% \end{center}
\end{figure}









\begin{figure}[htb]
 \centering
\subfloat[]{\includegraphics[width=3.5cm]{Pics/ReferenceImage-86.pdf}}\hspace{0pt}
\subfloat[]{\includegraphics[width=3.5cm]{Pics/SourceImage-63.pdf}}\hspace{0pt}

\subfloat[]{\includegraphics[width=3.5cm]{Pics/RegisterdMIRT-63-88.pdf}}\hspace{0pt}
\subfloat[]{\includegraphics[width=3.5cm]{Pics/RegisterdSource_qr-63-88.pdf}}\hspace{0pt}
\subfloat[]{\includegraphics[width=3.5cm]{Pics/SCoBeP-3D_63-88.pdf}}\hspace{0pt}

\subfloat[]{\includegraphics[width=3.5cm]{Pics/RefImage-86-Zoom.pdf}}\hspace{0pt}
\subfloat[]{\includegraphics[width=3.5cm]{Pics/SourceImage-63_Zoom.pdf}}\hspace{0pt}

\subfloat[]{\includegraphics[width=3.5cm]{Pics/RegisterdMIRT-63-88_Zoom.pdf}}\hspace{0pt}
\subfloat[]{\includegraphics[width=3.5cm]{Pics/RegisterdSource_qr-63-88_Zoomin.pdf}}\hspace{0pt}
\subfloat[]{\includegraphics[width=3.5cm]{Pics/SCoBeP-3D_63-88-Zoom.pdf}}\hspace{0pt}
\caption{{Registration result of the lung CT images that were captured with three months gap. (a) Source image; (b) Reference image; (c) MIRT \cite{myronenko2010intensity} [RMSE: 7.71]; (d) GP-Registration \cite{periaswamy2003elastic} [RMSE: 7.38]; (e) 3D-SCoBeP [RMSE: 4.18]; (f) Source image (zoom in); (g) Reference image (zoom in); (h) MIRT \cite{myronenko2010intensity} (zoom in);  (i) GP-Registration \cite{periaswamy2003elastic} (zoom in); (j) 3D-SCoBeP (zoom in).}}
\label{nafiss:fig63-86}
% \end{center}
\end{figure}


Fig. \ref{lung_3d} shows the result of the 3D-SCoBeP and MIRT\cite{myronenko2010intensity} with a 3D perspective. In this figure, I decided to show only a part of the CT images because the inside details of the lung are  more important than the tissue around it. Figs. \ref{lung_3d}(a) and \ref{lung_3d}(b) are the reference and the source CT image, respectively. Fig. \ref{lung_3d}(c) shows the result of 3D-SCoBeP where I used the voxel of the reference data to synthesize the source data. I created one RGB image where its first channel was assigned to the intensity of source CT image and its third channel was equal to 255. I assigned the reference CT image intensity, the MIRT result and the 3D-SCoBeP result to the second channel, respectively. Therefore, Fig. \ref{lung_3d}(d) corresponds to the initial state and Fig. \ref{lung_3d}(e) and \ref{lung_3d}(f) are final state of the MIRT and the 3D-SCoBeP. Note that in Fig. \ref{lung_3d}(c), I display a pure result of the proposed method which 
only the voxels with a probability more than the threshold $\theta$ was shown, therefore there are some dark voxels in this figure. For those voxels with probability less than $\theta$, the analysis of motion fields of neighborhood voxels could be used to estimate their motions. It means, one can extract the motions of voxels which have probability more than $\theta$ in each direction and apply the graph-cuts method \cite{boykov2001fast}, median filter, or moving average to estimate the motion of voxels with probability less than threshold $\theta$. 


I now proceed to compare the 3D-SCoBeP with other approaches; Figs. \ref{nafiss:fig42-65} and \ref{nafiss:fig63-86} show the output of my proposed method compared to two of the state-of-the-art methods: the MIRT \cite{myronenko2010intensity} and GP-Registration \cite{periaswamy2003elastic}. In these figures, I select only one slice of CT image to show the weaknesses and strengths of each technique. Figs. \ref{nafiss:fig42-65}(a) and \ref{nafiss:fig63-86}(a) correspond to the reference CT image and Figs. \ref{nafiss:fig42-65}(b) and \ref{nafiss:fig63-86}(b) correspond to the source CT image. Figs. \ref{nafiss:fig42-65}(c)-(e) and \ref{nafiss:fig63-86}(c)-(e) show results using MIRT \cite{myronenko2010intensity}, GP-Registration \cite{periaswamy2003elastic} and 3D-SCoBeP.
The warped images using MIRT with highlighted artifacts are shown in Figs. \ref{nafiss:fig42-65}(h) and \ref{nafiss:fig63-86}(h) and the warped images using GP-Registration with highlighted artifacts are shown in Figs. \ref{nafiss:fig42-65}(i) and \ref{nafiss:fig63-86}(i). The estimated images generated from the 3D-SCoBeP with highlighted areas are shown in \ref{nafiss:fig42-65}(j) and \ref{nafiss:fig63-86}(j). 








To quantify my  registration performance, I used the root mean square error (RMSE) measure between the true and estimated transformations:

\begin{equation}
 \varepsilon_{RMSE} = \sqrt{(1/N)\sum \left \| \tau - \hat\tau \right \|^2},
\end{equation}
where $N$ is the number of voxels in the reference and $\tau$ and $\hat\tau$ are the source image and the estimated transformation respectively.However, the RMSE is insufficient to qualify the accuracy of the registration methods. It can only give a rough estimation of similarity between the estimated image and the source image \cite{li1995contour,haque2012low}. 
In the term of the RMSE, I compare the source image with the output of the 3D-SCoBeP, MIRT and GP-Registration and the results are shown under (c), (d) and (e) of Figs. \ref{nafiss:fig42-65} and \ref{nafiss:fig63-86}. However, in Fig. \ref{nafiss:fig42-65}, the RMSE values for the MIRT and GP-Registration are 24.31 and 28.78, respectively. Although, The MIRT generates the result with the lower RMSE value, the GP-Registration preserves the structures better. My proposed technique shows an improved qualitative and quantitative result. The 3D-SCoBeP preserves the structure as precise as GP-Registration and also has the less RMSE value, 21.73, in comparison to two other methods. In Fig. \ref{nafiss:fig63-86}, all three methods show the same structure preservation but the RMSE value of the 3D-SCoBeP is much less than the RMSE value of MIRT and GP-Registration.





\subsection{3D MR Image Registration Taken In Different Directions}\label{3DMRSCOBEP}

{\color{black}{In this section, I am trying to align the medical data shown in Figs. \ref{brain_3d_Vertical} and \ref{brain_3d_Horizontal} where the brain MR images captured in parallel to the X--Z (sagittal) and the X--Y (transverse) planes, respectively. The X-axis is from left to right along the column direction in Fig. \ref{brain_3d_Vertical} and is from anterior to posterior along each slice in Fig. \ref{brain_3d_Horizontal}. The Z-axis and Y-axis are from the first slice to the last slice along the plane direction in  Figs. \ref{brain_3d_Vertical} and \ref{brain_3d_Horizontal}, respectively.}}


\begin{figure}[h]
 \centering
\includegraphics[width=10cm]{Pics/vertical}
% \vspace{-25pt}
\caption{ Brain MR images captured in parallel to X--Y (transverse) plane. }
\label{brain_3d_Vertical}
% \vspace{-15pt}
\end{figure}



\begin{landscape}
\begin{figure}[htb]
 \centering

\includegraphics[height=9cm]{Pics/Horizontal}

\caption{ Brain MR images captured in parallel to X--Z (sagittal) plane. }
\label{brain_3d_Horizontal}
% \vspace{-15pt}
\end{figure}
\end{landscape}


Here, the captured data are slices of brain image, where each one has $320 \times 320$ voxels. In this case, I use 200 slices which were captured in parallel to the X--Y (sagittal) plane as the reference data and 180 slices in parallel to the X--Z (transverse) plane as the source data. For registering this kind of 3D data, I first map each reference and source to a template data in size of $320 \times 320 \times 320$ using B-Spline interpolation \cite{lehmann2001addendum}.  {\color{black}{Since resampling could lead to an image to a new set of coordinates often provides a loss in image quality, the interpolating should be implemented with great care. Although there will be some expenses in computing time, the image quality can be improved by resampling using the B-Spline interpolation function.}}



Note that I set $a = b = c = 7$ and keep all the parameters as the same as previous section. Then, I apply the 3D-SCoBeP on the interpolated data. Fig. \ref{brain_3d}(a) and \ref{brain_3d}(b) show the reference and the source data, respectively. Fig. \ref{brain_3d}(c) shows the output of 3D-SCoBeP. In Fig. \ref{brain_3d}(d), I present the motion field in one of the slices in parallel to the X--Y (sagittal) plane which the darker points have less movement (minimum three pixels movement) and the lighter points are the area with more displacement (maximum sixteen pixels movement). In the other words, the selected area has three pixels translation in compared to the reference data. Also, the distance between the slices in Y direction is not uniform. These slices are closer to each other in the area that is darker in Fig. \ref{brain_3d}(d). Note that all selected voxels' probabilities are bigger than the threshold $\theta$.      


% {\color{red}{B-Spline comparison fig}}

% {\color{red}{3D MR image registration taken in different directions comparison with the other methods fig}}
%(sagittal)(transverse)

\begin{figure*}[t]
 \centering
%  \subfloat[]{\includegraphics[width=10cm]{Pics/3D/3d1}}
%  \subfloat[]{\includegraphics[width=10cm]{Pics/3D/3d2}}
\subfloat[]{\includegraphics[width=7cm]{Pics/3D/s40angle}}
\subfloat[]{\includegraphics[width=7cm]{Pics/3D/40angle}}

\subfloat[]{\includegraphics[width=7cm]{Pics/3D/output_40angle}}
\subfloat[]{\includegraphics[width=5.5cm]{Pics/3D/ydisparity}}
\caption{Result of 3D-SCoBeP on Brain CT images. (a) The reference CT image where the CT slices were taken in parallel to X--Z (sagittal) plane; (b) The source CT image where the CT slices were taken in parallel to X--Y (transverse) plane; (c) The 3D-SCoBeP result; (d) motion field in the X--Y (transverse) plane for one selected slice. }
\label{brain_3d}
% \end{center}
\end{figure*}


{\color{black}{Fig. \ref{SplineNot} shows the output of 3D-SCoBeP on brain data. While in Fig. \ref{SplineNot}(a), I applied 3D-SCoBeP directly on the original reference and source data with no interpolation approach, Fig. \ref{SplineNot}(b) shows the result of 3D-SCoBeP when I map each reference and source data to a template data using B-Spline interpolation. As seen in this figure, the registration with no interpolation performs poorly with significant misalignment where the interpolation approach brings the accurate performance in registration images.}}

\begin{figure}[t]
 \centering
\subfloat[]{\includegraphics[width=7cm]{Pics/3D/output_40angle}}
\subfloat[]{\includegraphics[width=6cm]{Pics/3D/Non_interpolated2}}
\caption{Result of 3D-SCoBeP on Brain CT images. (a) The 3D-SCoBeP result with B-Spline interpolation; (b) The 3D-SCoBeP result without B-Spline interpolation; }
\label{SplineNot}
% \end{center}
\end{figure}



\section{Discussion and Conclusion}\label{nafiss:conclusion}
In conclusion, I have proposed an effective registration method based on a sparse coding and belief propagation. The proposed method can be used for both rigid and non-rigid registration. My technique executes registration by first running sparse coding over an overcomplete dictionary out of 3D features of the reference image to gather possible match candidates.  Belief propagation is then applied to eliminate bad candidates and to select optimum matches.  The experimental result illustrates that my proposed algorithm compares with the high accuracy MIRT method by Myronenko and Song \cite{myronenko2010intensity} and the state-of-the-art GP-Registration by Periaswamy and Farid \cite{periaswamy2003elastic} over the CT and MR images.     
{\color{black}{The key advantage of my proposed method lies in the fact that it is applicable to both global and local registration \cite{shum2000systems}. Since the voxels act independently in my proposed method, by changing the size of the dictionary and using a part of data to make it, one can employ 3D-SCoBeP for medical image registration with higher speed locally.
In addition, I used a 3D factor graph for the entire 3D data instead of using a 2D factor graph per slice. This new strategy along with using graph-cut method that refines the disparity of each voxel enable us to register 3D data captured parallel to the sagittal plane into 3D data captured parallel to the transverse plane.}}







\chapter{GENOME SEQUENCE ALIGNMENT USING EMPIRICAL TRANSITION PROBABILITY, SPARSE CODING AND BELIEF PROPAGATION}

The latest sequencing technologies have generated numerous sequenced genomes for various species. This increasing volume of data requires tools that can accurately compare multiple genome sequences to aid in the study of populations, pan-genomes, and genome evolution \cite{li2009soap2,li2009fast}. For a particular study, many individual genomes may be sequenced to investigate genetic diversity. For example, the Cancer Genome Atlas \cite{hampton2006cancer} and 1000 Genomes Project \cite{siva20081000} will generate genome sequences from several thousand people. The complete bacterial genomes in public databases are already over one thousand. To better utilize this huge amount of sequenced genome information, many tools have been developed that are capable of efficiently finding similar sequences from whole genomes. %Whole genome sequence alignment are used for studying genome evolution and genetic diversity \cite{batzoglou2005many, dewey2006evolution}. For example, Blanchette \textit{et al.}, defined a Threaded Blockset Aligner (TBA) and built a threaded blockset under the assumption that all matching segments occur in the same order and orientation in a given sequence \cite{blanchette2004aligning}. TBA was designed for aligning megabase-sized regions of multiple mammalian genomes. Darling \textit{et al.} \cite{darling2004mauve} implemented a method for identification and alignment of conserved genomic DNA in the presence of rearrangements and horizontal transfer called Mauve. Mauve has been applied to align nine enterobacterial genomes and to determine global rearrangement structure in three mammalian genomes. There are other whole-genome alignment tools that can align multiple whole genomes such as \cite{dubchak2009multiple,hohl2002efficient,paten2008enredo}.

\section{Introduction}\label{introduction}

In bioinformatics, sequence alignment is an important way to identify similar regions that might be associated with similar functional and structural relationship between sequences. With the quick growth of genomic data, it is important to develop effective sequence alignment techniques that are scalable. The past decade has witnessed the development of many sequence alignment technologies. Cancers are caused by the collection of genomic sequence changes \cite{meyerson2010advances}. Therefore, alignment and analyses of cancer genome sequences provide basics to understand cancer biology, diagnosis and therapy.

In general, pairwise sequence alignment methods can be classified into local and global approaches. The global alignment attempts to find the best match between two strings with similar lengths through global optimization. In contrast, the local alignment is usually used to identify regions of similarity between a short query and a longer sequence. Global alignments \cite{needleman1970general,morgenstern1999dialign,bray2003avid,brudno2003lagan} are less prone
to demonstrating false homology as each letter of one sequence is constrained to being aligned to only one letter of the other. Local alignments \cite{smith1981comparison,altschul1990basic,brudno2002fast,schwartz2003human}, on the other hand, can cope with rearrangements between non-syntenic, orthologous sequences by identifying similar regions in sequences; this, however, comes at the expense of a higher false positive rate due to the inability of local aligners to take into account overall conservation maps \cite{brudno2003glocal}. 

A lot of efforts have been made to improve the efficiency and efficacy of sequence alignments. The ClustalW program proposed by Thompson and Larkin \cite{thompson1994clustal,larkin2007clustal} uses a multi-stage mechanism to weight and to align sub-sequences based on sequence divergences. In addition, sequence annealing technique incrementally builds sequence alignment one at a time by checking whether a single match is consistent with a partial multiple alignments \cite{schwartz2007multiple}. 
Darling {\em et al.} proposed a hidden Markov model that uses a sum-of-pairs breakpoint score to facilitate the detection of rearrangement breakpoints, when genomes have unequal gene content \cite{darling2010progressivemauve}. 
Mummer is a highly efficient suffix tree based matching tool for whole genome alignment, as well as incomplete genomes \cite{kurtz2004versatile}. 

Researchers also proposed heuristics to accelerate sequence alignment. For example, the bounded sparse dynamic programming (BSDP) is used to support rapid approximation of exhaustive alignment in \cite{slater2005automated}. Another heuristic-driven approach, namely FastTree, is a tree-based method that stores profiles of internal nodes in a tree, such that candidate joins can be quickly identified. FastTree is also scalable for handling alignments over 10,000 sequences \cite{price2009fasttree, price2010fasttree}. 

Maximum-likelihood based approaches like PhyML and RAxML-VI-HPC have been developed as well. PhyML \cite{guindon2003simple} used a hill-climbing algorithm that adjusts tree topology and branch length at each tree modification iteration. RAxML-VI-HPC \cite{stamatakis2006raxml}, which stands for randomized accelerated maximum likelihood for high performance computing, takes advantages of a parallel program to support large-scale genome alignment. 




% The latest sequencing technologies have generated numerous sequenced genomes for various species. This increasing volume of data requires tools that can accurately compare multiple genome sequences to aid in the study of populations, pan-genomes, and genome evolution \cite{li2009soap2,li2009fast}. For a particular study, many individual genomes may be sequenced to investigate genetic diversity. For example, the Cancer Genome Atlas \cite{hampton2006cancer} and 1000 Genomes Project \cite{siva20081000} will generate genome sequences from several thousand people. The complete bacterial genomes in public databases are already over one thousand. To better utilize this huge amount of sequenced genome information, many tools have been developed that are capable of efficiently finding similar sequences from whole genomes.

Whole genome sequence alignment are used for studying genome evolution and genetic diversity \cite{batzoglou2005many, dewey2006evolution}. For example, Blanchette \textit{et al.}, defined a Threaded Blockset Aligner (TBA) and built a threaded blockset under the assumption that all matching segments occur in the same order and orientation in a given sequence \cite{blanchette2004aligning}. TBA was designed for aligning megabase-sized regions of multiple mammalian genomes. Darling \textit{et al.} \cite{darling2004mauve} implemented a method for identification and alignment of conserved genomic DNA in the presence of rearrangements and horizontal transfer called Mauve. Mauve has been applied to align nine enterobacterial genomes and to determine global rearrangement structure in three mammalian genomes. There are other whole-genome alignment tools that can align multiple whole genomes such as \cite{dubchak2009multiple,hohl2002efficient,paten2008enredo}. 


Whole-genome alignment tools are classified from collinear multiple sequence alignment tools, such as tools in \cite{bradley2009fast,edgar2004muscle,thompson1994clustal} where they can align very long sequences and detect the presence of rearrangements, duplications, and large-scale sequence gains and losses. For example, Bradley \textit{et al.}, in \cite{bradley2009fast} proposed a program for the alignment of multiple biological sequences that is statistically motivated and fast for practical size problems. It was based on pair hidden Markov models which approximate an insertion/deletion process on a tree and used a sequence annealing algorithm to combine the posterior probabilities estimated from these models into a multiple alignment. Edgar \textit{et al.} proposed another alignment tool named MUSCLE which is a program for creating multiple alignments of protein sequences \cite{edgar2004muscle}. Elements of the algorithm include fast distance estimation using \textit{k-mer} counting, progressive alignment using a log-expectation score, and refinement using tree-dependent restricted partitioning. In spite of collinear alignment technologies, non-collinear alignment such as \cite{bartovs2012intraspecific,song2003gene} contains the elements that are arranged in some non-linear order (see Fig. \ref{Collinear}).  

\begin{figure}[tbh]
\centering
\subfloat[Collinear]{\includegraphics[height=2cm]{Pics/Colinear2}}

\subfloat[Non-collinear]{\includegraphics[height=2cm]{Pics/Colinear1}}
\caption{{Collinear vs. Non--collinear nucleotide sequence alignment.}}
\label{Collinear}
 
\end{figure} 

In this chapter, I propose a novel alignment method that uses sparse coding \cite{olshausen1997sparse} and empirical transition probability to tackle the scalability challenge. Thanks to the sparse representation, my mechanism can handle long sequences with reduced memory footprint. I also leverage belief propagation to combine local and neighboring information of candidate nucleotides into consideration and generate matching scores to determine the best match. First, I index the reference and the read genome sequence using empirical transional probability and pick the top score indexes from the reference genome sequence to build an over-complete dictionary. I then find a set of candidate nucleotide for each nucleotide of the test sequence using sparse coding from the constructed dictionary. The match score of each candidate nucleotide will be evaluated taking both local and neighboring information into account using belief propagation. 
The rest of this chapter is structured as follows. Section \ref{SR_PropsedMethod} introduces my proposed method. Section \ref{EXPERIMENTAL} presents my results, including the comparison against SOAP aligner \cite{li2009soap2} and BWA \cite{li2009fast}. Finally, I draw my conclusions in Section \ref{conclusion}.

% In this paper, I propose a novel alignment technique using sparse coding and belief propagation. First, I build an overcomplete dictionary out of extracted nucleotide features of a reference sequence. I then find a set of candidate nucleotide for each nucleotide of the test sequence using sparse coding from the constructed dictionary. The match score of each candidate nucleotide will be evaluated taking both local and neighboring information into account using belief propagation. The best match will be selected as the candidate with the highest score. The rest of this paper is structured as follows. Section \ref{SR_PropsedMethod} introduces my proposed method. In Section \ref{EXPERIMENTAL}, I show and discuss my simulation results and compare them with SOAP aligner \cite{li2009soap2} and BWA \cite{li2009fast}, followed by a brief conclusion in Section \ref{conclusion}.

% {\color{red}{The current indexing techniques:
% \begin{enumerate}
%  \item Suffix array based
%  \item Burrows–Wheeler transform based
%  \item Anchor based 
%  \item Statistical based
%
%
% \item distance-based \cite{schwartz2007posterior,schwartz2007multiple}
%
%
% Comparative genomics analyses, such as genome annotation, phylogenetic inference, protein structure prediction and protein function prediction, rely heavily on multiple sequence alignment (MSA), which is essentially a homology prediction problem at the character level over multiple sequences. Alignments are also crucial in other settings, including natural language processing, as seen in word alignment in machine translation, or in the problem of citation entity alignment that is introduced in this dissertation. Fundamental as it is, MSA is a hard computational problem, and producing accurate MSAs remains challenging. In this thesis I develop new posterior decoding based methods for multiple alignments of discrete characters. I begin by discussing biological sequence alignment, and describe a new alignment accuracy measure (AMA), which is based on a metric for the space of alignments. In the case of pairwise sequence alignment, it is possible to ﬁnd an alignment that maximizes the expected AMA value using posterior decoding. I then describe an extension of the pairwise alignment method to MSA using 1a novel approach that is based on a poset representation of multiple alignments, and an efficient online topological ordering algorithm. This leads to a sequence annealing algorithm, which is an incremental method for building MSAs one match at a time. 
% The sequence annealing algorithm outperforms all existing algorithms on benchmark test sets of protein sequence alignments. It simultaneously achieves high recall and high precision, dramatically reducing the number of incorrectly aligned residues in comparison to other programs. The method can adjust the recall/precision trade off, and can reliably identify homologous regions among protein sequences. Finally, I demonstrate that similar posterior decoding methods can be applied to the problem of multiple word alignment of citation sentences (citances) in the bioscience literature. I extend the deﬁnition of AMA to many-to-many word alignments, and develop a posterior decoding algorithm for multiple citance alignment.
% %
% %
% Motivation: I introduce a novel approach to multiple alignment that is based on an algorithm for rapidly checking whether single matches are consistent with a partial multiple alignment. This leads to a sequence annealing algorithm, which is an incremental method for building multiple sequence alignments one match at a time. My approach improves significantly on the standard progressive alignment approach to multiple alignment.
% %
% Results: The sequence annealing algorithm performs well on benchmark test sets of protein sequences. It is not only sensitive, but also specific, drastically reducing the number of incorrectly aligned residues in comparison to other programs. The method allows for adjustment of the sensitivity/specificity tradeoff and can be used to reliably identify homologous regions among protein sequences.
% %
% Availability: An implementation of the sequence annealing algorithm is available at http://bio.math.berkeley.edu/amap/
%
% Contact: sariel@cs.berkeley.edu
% %
% %
% \item sequence annealing technique \cite{schwartz2007multiple}
% 
% Motivation: I introduce a novel approach to multiple alignment that is based on an algorithm for rapidly checking whether single matches are consistent with a partial multiple alignment. This leads to a sequence annealing algorithm, which is an incremental method for building multiple sequence alignments one match at a time. My approach improves significantly on the standard progressive alignment approach to multiple alignment.
% 
% Results: The sequence annealing algorithm performs well on benchmark test sets of protein sequences. It is not only sensitive, but also specific, drastically reducing the number of incorrectly aligned residues in comparison to other programs. The method allows for adjustment of the sensitivity/specificity tradeoff and can be used to reliably identify homologous regions among protein sequences.
% 
% Availability: An implementation of the sequence annealing algorithm is available at http://bio.math.berkeley.edu/amap/
% 
% 
% \item MUMmer suffix trie matching tool \cite{kurtz2004versatile}
% 
% The newest version of MUMmer easily handles comparisons of large eukaryotic genomes at varying evolutionary distances, as demonstrated by applications to multiple genomes. Two new graphical viewing tools provide alternative ways to analyze genome alignments. The new system is the first version of MUMmer to be released as open-source software. This allows other developers to contribute to the code base and freely redistribute the code. The MUMmer sources are available at http://www.tigr.org/software/mummer webcite.
% 
% \item the exonerate homology-search program \cite{slater2005automated}
% 
% Background
% Exhaustive methods of sequence alignment are accurate but slow, whereas heuristic approaches run quickly, but their complexity makes them more difficult to implement. I introduce bounded sparse dynamic programming (BSDP) to allow rapid approximation to exhaustive alignment. This is used within a framework whereby the alignment algorithms are described in terms of their underlying model, to allow automated development of efficient heuristic implementations which may be applied to a general set of sequence comparison problems.
% 
% Results
% The speed and accuracy of this approach compares favourably with existing methods. Examples of its use in the context of genome annotation are given.
% 
% Conclusions
% This system allows rapid implementation of heuristics approximating to many complex alignment models, and has been incorporated into the freely available sequence alignment program, exonerate.
% 


% 
% \item FastTree 1.0 \cite{price2009fasttree}
% 
% Gene families are growing rapidly, but standard methods for inferring phylogenies do not scale to alignments with over 10,000 sequences. I present FastTree, a method for constructing large phylogenies and for estimating their reliability. Instead of storing a distance matrix, FastTree stores sequence profiles of internal nodes in the tree. FastTree uses these profiles to implement Neighbor-Joining and uses heuristics to quickly identify candidate joins. FastTree then uses nearest neighbor interchanges to reduce the length of the tree. For an alignment with N sequences, L sites, and a different characters, a distance matrix requires O($N^2$) space and O($N^2L$) time, but FastTree requires just O($NLa + N\sqrt{N}$) memory and O($N\sqrt{N}\log {N}La$) time. To estimate the tree's reliability, FastTree uses local bootstrapping, which gives another 100-fold speedup over a distance matrix. For example, FastTree computed a tree and support values for 158,022 distinct 16S ribosomal RNAs in 17 h and 2.4 GB of memory. Just computing pairwise Jukes–Cantor distances and storing them, without inferring a tree or bootstrapping, would require 17 h and 50 GB of memory. In simulations, FastTree was slightly more accurate than Neighbor-Joining, BIONJ, or FastME; on genuine alignments, FastTree's topologies had higher likelihoods. FastTree is available at http://microbesonline.org/fasttree.
% 
% 
% \item maximum-likelihood based like PhyML \cite{guindon2003simple} and RAxML-VI-HPC \cite{stamatakis2006raxml}

% The increase in the number of large data sets and the complexity of current probabilistic sequence evolution models necessitates fast and reliable phylogeny reconstruction methods. I describe a new approach, based on the maximumlikelihood principle, which clearly satisﬁes these requirements. The core of this method is a simple hill-climbing algorithm that adjusts tree topology and branch lengths simultaneously. This algorithm starts from an initial tree built by a fast distance-based method and modiﬁes this tree to improve its likelihood at each iteration. Due to this simultaneous adjustment of the topology and branch lengths, only a few iterations are sufﬁcient to reach an optimum. I used extensive and realistic computer simulations to show that the topological accuracy of this new method is at least as high as that of the existing maximum-likelihood programs and much higher than the performance of distance-based and parsimony approaches. The reduction of computing time is dramatic in comparison with other maximum-likelihood packages, while the likelihood maximization ability tends to be higher. For example, only 12 min were required on a standard personal computer to analyze a data set consisting of 500 rbcL sequences with 1,428 base pairs from plant plastids, thus reaching a speed of the same order as some popular distance-based and parsimony algorithms. This new method is implemented in the PHYML program, which is freely available on my web page: http://www.lirmm.fr/w3ifa/MAAS/. [Algorithm; computer simulations; maximum likelihood; phylogeny; rbcL; RDPII project.]
% 
% RAxML-VI-HPC (randomized axelerated maximum likelihood for high performance computing) is a sequential and parallel program for inference of large phylogenies with maximum likelihood (ML). Low-level technical optimizations, a modification of the search algorithm, and the use of the GTR+CAT approximation as replacement for GTR+Γ yield a program that is between 2.7 and 52 times faster than the previous version of RAxML. A large-scale performance comparison with GARLI, PHYML, IQPNNI and MrBayes on real data containing 1000 up to 6722 taxa shows that RAxML requires at least 5.6 times less main memory and yields better trees in similar times than the best competing program (GARLI) on datasets up to 2500 taxa. On datasets ≥4000 taxa it also runs 2–3 times faster than GARLI. RAxML has been parallelized with MPI to conduct parallel multiple bootstraps and inferences on distinct starting trees. The program has been used to compute ML trees on two of the largest alignments to date containing 25 057 (1463 bp) and 2182 (51 089 bp) taxa, respectively.
% 
% 
% 
% \item FastTree 2 \cite{price2010fasttree}
% 
% Background
% 
% I recently described FastTree, a tool for inferring phylogenies for alignments with up to hundreds of thousands of sequences. Here, I describe improvements to FastTree that improve its accuracy without sacrificing scalability.

% Methodology/Principal Findings
% 
% Where FastTree 1 used nearest-neighbor interchanges (NNIs) and the minimum-evolution criterion to improve the tree, FastTree 2 adds minimum-evolution subtree-pruning-regrafting (SPRs) and maximum-likelihood NNIs. FastTree 2 uses heuristics to restrict the search for better trees and estimates a rate of evolution for each site (the “CAT” approximation). Nevertheless, for both simulated and genuine alignments, FastTree 2 is slightly more accurate than a standard implementation of maximum-likelihood NNIs (PhyML 3 with default settings). Although FastTree 2 is not quite as accurate as methods that use maximum-likelihood SPRs, most of the splits that disagree are poorly supported, and for large alignments, FastTree 2 is 100–1,000 times faster. FastTree 2 inferred a topology and likelihood-based local support values for 237,882 distinct 16S ribosomal RNAs on a desktop computer in 22 hours and 5.8 gigabytes of memory.
% 
% 
% \item CAT-Based Branch Lengths \cite{stamatakis2006phylogenetic}
% 
% 
% Inference of phylogenetic trees using the maximum likelihood (ML) method is NP-hard. Furthermore, the computation of the likelihood function for huge trees of more than 1,000 organisms is computationally intensive due to a large amount of ﬂoating point operations and high memory consumption. Within this context, the present paper compares two competing mathematical models that account for evolutionary rate heterogeneity: the Γ and CAT models. The intention of this paper is to show that—from a purely empirical point of view—CAT can be used instead of Γ. The main advantage of CAT over $\Gamma$ consists in signiﬁcantly lower memory consumption and faster inference times. An experimental study using RAxML has been performed on 19 real-world datasets comprising 73 up to 1,663 DNA sequences. Results show that CAT is on average 5.5 times faster than $\Gamma$ and—surprisingly enough—also yields trees with slightly superior $\Gamma$ likelihood values. The usage of the CAT model decreases the amount of average L2 and L3 cache misses by factor 8.55
% 
% \item Global alignment algorithms \cite{needleman1970general,morgenstern1999dialign,bray2003avid,brudno2003lagan}
% 
% Motivation: The performance and time complexity of an improved version of the segment-to-segment approach to multiple sequence alignment is discussed. In this approach, alignments are composed from gap-free segment pairs, and the score of an alignment is defined as the sum of so-called weights of these segment pairs. 
% 
% Results: A modification of the weight function used in the original version of the alignment program DIALIGN has two important advantages: it can be applied to both globally and locally related sequence sets, and the running time of the program is considerably improved. The time complexity of the algorithm is discussed theoretically, and the program running time is reported for various test examples. 
% 
% Availability: The program is available on-line at the Bielefeld University Bioinformatics Server (BiBiServ) http://bibiserv.TechFak.Uni-Bielefeld.DE/dialign/ Contact: burkhard.morgenstern@rp-rorer.co.uk
% 
% 
% In this paper I describe a new global alignment method called AVID. The method is designed to be fast, memory efficient, and practical for sequence alignments of large genomic regions up to megabases long. I present numerous applications of the method, ranging from the comparison of assemblies to alignment of large syntenic genomic regions and whole genome human/mouse alignments. I have also performed a quantitative comparison of AVID with other popular alignment tools. To this end, I have established a format for the representation of alignments and methods for their comparison. These formats and methods should be useful for future studies. The tools I have developed for the alignment comparisons, as well as the AVID program, are publicly available. See Web Site References section for AVID Web address and Web addresses for other programs discussed in this paper.

% To compare entire genomes from different species, biologists increasingly need alignment methods that are efficient enough to handle long sequences, and accurate enough to correctly align the conserved biological features between distant species. I present LAGAN, a system for rapid global alignment of two homologous genomic sequences, and Multi-LAGAN, a system for multiple global alignment of genomic sequences. I tested my systems on a data set consisting of greater than 12 Mb of high-quality sequence from 12 vertebrate species. All the sequence was derived from the genomic region orthologous to an ∼1.5-Mb region on human chromosome 7q31.3. I found that both LAGAN and Multi-LAGAN compare favorably with other leading alignment methods in correctly aligning protein-coding exons, especially between distant homologs such as human and chicken, or human and fugu. Multi-LAGAN produced the most accurate alignments, while requiring just 75 minutes on a personal computer to obtain the multiple alignment of all 12 sequences. Multi-LAGAN is a practical method for generating multiple alignments of long genomic sequences at any evolutionary distance. My systems are publicly available at http://lagan.stanford.edu.
% 



% \item local alignment methods \cite{smith1981comparison,altschul1990basic,brudno2002fast,schwartz2003human}
% 
% A new approach to rapid sequence comparison, basic local alignment search tool (BLAST), directly approximates alignments that optimize ameasure of local similarity, the maximal segment pair (MSP) score. Recent mathmatical results on the stochastic properties of MSP scores allow an analysis  of the  perfornmnce of this method as well as the statistical significance of alignments it generates. The basic algorithm is simple and robust; it can be implemented in a number of ways and applied in a variety of contexts including straightforward DNA and protein sequence database searches, motif searches, gene identification searches, and in the analysis of multiple regions of similarity in long DNA sequences. In addition to its flexibility and tractability to mathematical analysis, BLAST is an order of magnitude faster than existing sequence comparison tools of comparable sensitivity.
% 
% 
% Comparative analysis of syntenic genome sequences can be used to identify functional sites such as exons and regulatory elements. Here, the first step is to align two or several evolutionary related sequences and, in recent years, a number of computer programs have been developed for alignment of large genomic sequences. Some of these programs are extremely fast but often time-efficiency is achieved at the expen,S'e of sensitivity. One way of combining speed and sensitivity is to use an anchored-alignment approach. In a first step, a fast heuristic identifies a chain of strong sequence similarities that serve as anchor points. In a second step, regions between these anchor points are aligned using a slower but more sensitive method.
% I present CHAOS, a novel algorithm for rapid identification of chains of local sequence similarities among large genomic sequences. Similarities identified by CHAOS are used as anchor points to improve the running time of the DIALIGN alignment program. Systematic test runs show that this method can reduce the running time of DIALIGN by more than 93\% while affecting the quality of the resulting alignments by only 1\%. 
% 
% The Mouse Genome Analysis Consortium aligned the human and mouse genome sequences for a variety of purposes, using alignment programs that suited the various needs. For investigating issues regarding genome evolution, a particularly sensitive method was needed to permit alignment of a large proportion of the neutrally evolving regions. I selected a program called BLASTZ, an independent implementation of the Gapped BLAST algorithm specifically designed for aligning two long genomic sequences. BLASTZ was subsequently modified, both to attain efficiency adequate for aligning entire mammalian genomes and to increase its sensitivity. This work describes BLASTZ, its modifications, the hardware environment on which I run it, and several empirical 
% studies to validate its results.
% 
% 
% 
% \item Glocal alignment method \cite{brudno2003glocal}
% 
% Motivation: To compare entire genomes from different species, biologists increasingly need alignment methods that are efﬁcient enough to handle long sequences, and accurate enough to correctly align the conserved biological features between distant species. The two main classes of pairwise alignments are global alignment, where one string is transformed into the other, and local alignment, where all locations of similarity between the two strings are returned. Global alignments are less prone
% 
% to demonstrating false homology as each letter of one sequence is constrained to being aligned to only one letter of the other. Local alignments, on the other hand, can cope with rearrangements between non-syntenic, orthologous sequences by identifying similar regions in sequences; this, however, comes at the expense of a higher false positive rate due to the inability of local aligners to take into account overall conservation maps. 
% 
% Results: In this paper I introduce the notion of glocal alignment, a combination of global and local methods, where one creates a map that transforms one sequence into the other while allowing for rearrangement events. I present Shufﬂe-LAGAN, a glocal alignment algorithm that is based on the CHAOS local alignment algorithm and the LAGAN global aligner, and is able to align long genomic sequences. To test Shufﬂe-LAGAN I split the mouse genome into BAC-sized pieces, and aligned these pieces to the human genome. I demonstrate that ShufﬂeLAGAN compares favorably in terms of sensitivity and speciﬁcity with standard local and global aligners. From the alignments I conclude that about 9\% of human/mouse homology may be attributed to small rearrangements, 63\% of which are duplications.
% 
% 
% \item positional homology genome alignment \cite{darling2004mauve}
% 
% As genomes evolve, they undergo large-scale evolutionary processes that present a challenge to sequence comparison not posed by short sequences. Recombination causes frequent genome rearrangements, horizontal transfer introduces new sequences into bacterial chromosomes, and deletions remove segments of the genome. Consequently, each genome is a mosaic of unique lineage-specific segments, regions shared with a subset of other genomes and segments conserved among all the genomes under consideration. Furthermore, the linear order of these segments may be shuffled among genomes. I present methods for identification and alignment of conserved genomic DNA in the presence of rearrangements and horizontal transfer. My methods have been implemented in a software package called Mauve. Mauve has been applied to align nine enterobacterial genomes and to determine global rearrangement structure in three mammalian genomes. I have evaluated the quality of Mauve alignments and drawn comparison to other methods through extensive simulations of genome evolution.
% 
% \item anchoring heuristics \cite{blanchette2004aligning,bray2003mavid,brudno2003lagan,hohl2002efficient,ogurtsov2002owen,brudno2004chaos}

% I define a “threaded blockset,” which is a novel generalization of the classic notion of a multiple alignment. A new computer program called TBA (for “threaded blockset aligner”) builds a threaded blockset under the assumption that all matching segments occur in the same order and orientation in the given sequences; inversions and duplications are not addressed. TBA is designed to be appropriate for aligning many, but by no means all, megabase-sized regions of multiple mammalian genomes. The output of TBA can be projected onto any genome chosen as a reference, thus guaranteeing that different projections present consistent predictions of which genomic positions are orthologous. This capability is illustrated using a new visualization tool to view TBA-generated alignments of vertebrate Hox clusters from both the mammalian and fish perspectives. Experimental evaluation of alignment quality, using a program that simulates evolutionary change in genomic sequences, indicates that TBA is more accurate than earlier programs. To perform the dynamic-programming alignment step, TBA runs a stand-alone program called MULTIZ, which can be used to align highly rearranged or incompletely sequenced genomes. I describe my use of MULTIZ to produce the whole-genome multiple alignments at the Santa Cruz Genome Browser.
% 
% MAVID is a multiple alignment program suitable for many large genomic regions. The MAVID web server allows biomedical researchers to quickly obtain multiple alignments for genomic sequences and to subsequently analyse the alignments for conserved regions. MAVID has been successfully used for the alignment of closely related species such as primates and also for the alignment of more distant organisms such as human and fugu. The server is fast, capable of aligning hundreds of kilobases in less than a minute. The multiple alignment is used to build a phylogenetic tree for the sequences, which is subsequently used as a basis for identifying conserved regions in the alignment. The server can be accessed at http://baboon.math.berkeley.edu/mavid/.
% 
% 
% Motivation: To allow a direct comparison of the genomic DNA sequences of sufficiently similar organisms, there is an urgent need for software tools that can align more than two genomic sequences. 
% 
% Results: I developed new algorithms and a software tool ‘Multiple Genome Aligner’ (MGA for short) that efficiently computes multiple genome alignments of large, closely related DNA sequences. For example, it can align 85\% percent of the complete genomes of six human adenoviruses (average length 35 305 bp.) in 159 seconds. An alignment of 74\% of the complete genomes
% of three of strains of E. coli (lengths: 5 528 445; 5 498 450;4 639 221 bp.) is produced in 30 minutes.
% 
% Summary: OWEN is an interactive tool for aligning two long DNA sequences that represents similarity between them by a chain of collinear local similarities. OWEN employs several methods for constructing and editing local similarities and for resolving conflicts between them. Alignments of sequences of lengths over 106 can often be produced in minutes. OWEN requires memory below 20 L, where L is the sum of lengths of the compared sequences.
% 
% Cross-species sequence comparison is a powerful approach to analyze functional sites in genomic sequences and many discoveries have been made based on genomic alignments. Herein, I present a WWW-based software system for multiple alignment of large genomic sequences. My server utilizes the previously developed combination of CHAOS and DIALIGN to achieve both speed and alignment accuracy. CHAOS is a fast database search tool that creates a list of local sequence similarities. These are used by DIALIGN as anchor points to speed up the final alignment procedure. The resulting alignment is returned to the user in different formats together with a list of anchor points found by CHAOS. The CHAOS/DIALIGN software is freely available at http://dialign.gobics.de/chaos-dialign-submission.
% 
% \item glocal multiple genome alignments \cite{paten2008enredo,dubchak2009multiple}
% 
% Pairwise whole-genome alignment involves the creation of a homology map, capable of performing a near complete transformation of one genome into another. For multiple genomes this problem is generalized to finding a set of consistent homology maps for converting each genome in the set of aligned genomes into any of the others. The problem can be divided into two principal stages. First, the partitioning of the input genomes into a set of colinear segments, a process which essentially deals with the complex processes of rearrangement. Second, the generation of a base pair level alignment map for each colinear segment. I have developed a new genome-wide segmentation program, Enredo, which produces colinear segments from extant genomes handling rearrangements, including duplications. I have then applied the new alignment program Pecan, which makes the consistency alignment methodology practical at a large scale, to create a new set of genome-wide mammalian alignments. I test both Enredo and Pecan using novel and existing assessment analyses that incorporate both real biological data and simulations, and show that both independently and in combination they outperform existing programs. Alignments from my pipeline are publicly available within the Ensembl genome browser.
% 
% Multiple sequence alignments have become one of the most commonly used resources in genomics research. Most algorithms for multiple alignment of whole genomes rely either on a reference genome, against which all of the other sequences are laid out, or require a one-to-one mapping between the nucleotides of the genomes, preventing the alignment of recently duplicated regions. Both approaches have drawbacks for whole-genome comparisons. In this paper I present a novel symmetric alignment algorithm. The resulting alignments not only represent all of the genomes equally well, but also include all relevant duplications that occurred since the divergence from the last common ancestor. My algorithm, implemented as a part of the VISTA Genome Pipeline (VGP), was used to align seven vertebrate and six Drosophila genomes. The resulting whole-genome alignments demonstrate a higher sensitivity and specificity than the pairwise alignments previously available through the VGP and have higher exon alignment accuracy than comparable public whole-genome alignments. Of the multiple alignment methods tested, mine performed the best at aligning genes from multigene families—perhaps the most challenging test for whole-genome alignments. My whole-genome multiple alignments are available through the VISTA Browser at http://genome.lbl.gov/vista/index.shtml.
% 

% \end{enumerate}


% The current alignment
% \begin{enumerate}
% \item Whole genome sequence alignment
% \item Collinear sequence alignment
% \item Non-Collinear sequence alignment
% \end{enumerate}
% 
% }}


% {\color{red}{review \cite{darling2010progressivemauve,bradley2009fast,price2010fasttree,brudno2003glocal,chiaromonte2002scoring}}}

% Given the difficulties in assessing accuracy, recent development has included methods that are statistically motivated and show improved specificity \cite{bradley2009fast,paten2008enredo}. At shorter evolutionary distances with large fractions of identical sequences, there is less ambiguity in alignment outcomes. Yet, even within a bacterial species, aligning multiple genomes is not a trivial task, especially if the sequences contain rearrangements, duplications and exhibit sequence gain and loss. Also, despite relatively short chromosome lengths for bacteria, typically a few million base pairs, the computational complexity of multiple sequence alignment makes it a formidable challenge. Calculation of multiple alignments with a simple sum of pairs scoring scheme is known to be an NP-hard problem \cite{elias2006settling}, which makes calculation of an exact solution infeasible for large inputs. Multiple genome agenome alignment tools rely on heuristics to achieve reasonable run times. There are numerous methods to compare a single pair of wholegenome sequences \cite{bray2003avid,schwartz2003human}. The Nucmer and MUMmer package is a fast whole-genome alignment method that utilizes a sufﬁx tree to seed an alignment with maximal unique matches (MUMs) \cite{kurtz2004versatile}. The sufﬁx tree implementation of MUMmer is especially efficient and can be both

% \section{Background}


\section{Proposed Method}\label{SR_PropsedMethod}

In this section, I present my genome indexing and alignment framework in detail, where the proposed method includes three steps: indexing, index matching, and sequence matching. In this work, I refer to ``\textit{reference sequence}'' as the base-line sequence and try to align a ``\textit{read sequence}'' against the base-line sequence.

\subsection{Indexing}
The current genome indexing methods generate huge indices before performing the actual alignment to decrease the alignment time \cite{liu2012soap3,vouzis2011gpu}. The indexing process can be very time-consuming. In contrast, my proposed indexing technique provides a faster and light-weight alternative for index generation, which is similar to the big data retrieval systems that were proposed \cite{carson1999blobworld,doermann1998indexing,liu2000robust}.
These indices can reduce the search space and provide an estimation of the read sequence locations in the reference sequence. The proposed genome indexing technique models a nucleotide sequence as a graph by counting the transitions between each pair of nucleotides. To be more specific, as shown in  Fig. \ref{Transition_diag}, I consider a graph with four states according to the different types of nucleotides and sixteen vertices according to all possible transitions between nucleotides. I read the first nucleotide of the sequence and treat it as the initial state. Then, I move from one state to the other state by scanning the next nucleotide repeatedly till the end of the sequence. Afterwards, I calculate the number of nucleotide transitions where I count how many times I pass one vertex in the graph and store them in a $4 \times 4$ matrix. Finally, I normalize the resulting matrix as follows:
\begin{equation}
% \text{State Matrix}
I = \ \ \bordermatrix{ & A & C & G & T\cr
			 A& k_{aa} & k_{ac} & k_{ag} & k_{at}\cr
			 C& k_{ca} & k_{cc} & k_{cg} & k_{ct}\cr
			 G& k_{ga} & k_{gc} & k_{gg} & k_{gt}\cr
			 T& k_{ta} & k_{tc} & k_{tg} & k_{tt}}
			 \times \frac{1}{\sum_{s,w \in \{a,c,g,t\}} k_{sw}} 
\label{eq:ACGT}
\end{equation}
where $k_{sw}$ is the number that has the $S$-type nucleotide immediately before the $W$-type nucleotide.

If the length of a sequence is larger than a given threshold i.e., $h$, I divide it into subsequences with maximum length of $h$, where each subsequence will have $o$ nucleotides overlap with their neighbors. I set $o \geq \frac{h}{2}$ so that each pair of nucleotides can be counted at least twice. For each subsequence $i$, I count the transition of the nucleotides from the start of the subsequence till its end to reveal the number of different nucleotides that reside beside each other. In Fig. \ref{ExampleofIndexing}, an input sequence with $h = 250$ is used to demonstrate the proposed indexing process, where $I_i$ is the calculated index for the input sequence based on the transition graph shown on the left hand side. Finally, I normalize the transition matrix, which will be used to find the approximate location of each subsequence in the next step. 

\begin{figure}[t]

\centering
\subfloat{\includegraphics[width=13cm]{Pics/Diag1}}
\caption{{The transition diagram between nucleotides. $k_{sw}$ is the number of appearance  of the $W$-type nucleotide immediately after the $S$-type nucleotide where $s,w \in \{a,c,g,t\}$.}}
\label{Transition_diag}
\end{figure}

\begin{landscape}
\begin{figure}[t]
\centering
\subfloat{\includegraphics[width=22cm]{Pics/Genom_Index_Example}}
\caption{{An example of the indexing procedure for a small sample subsequence.}}
\label{ExampleofIndexing}
\end{figure}
\end{landscape}







\subsection{Index Matching}
The index matching step is designed to find similar indices based on global information of the sequence. I define a symmetric distance function between two index matrices $I$ and $J$ as follows: $D_{MSE}(I,J)=\left \| I-J  \right \|_f$, where $\left \| \cdot  \right \|_f$ is the Frobenius norm of the matrix. 


After generating the indices of the reference sequence and the read sequence, the $D_{MSE}$ distances to all reference sequence indices are calculated, where the top $t$ most similar indices in terms of $D_{MSE}$ are chosen as candidate indices. To find the best matched index, I resort to belief propagation (BP) on a factor graph.  In this dissemination, I provided a concise review about the BP algorithm on factor graph on Section \ref{BeliefPropagation}. Interested readers can check my earlier publications in \cite{barzigar2013scobep,roozgard2012medical,roozgard2013nucleotide} for more details about the factor graph design and the BP algorithm.

For each index in the read sequence, I consider a variable node and a factor node which connects each two neighbor variable nodes. An extra factor node of each variable node stores the prior probabilities which are calculated based on $D_{MSE}$ distance (see Fig. \ref{FG1}).

\begin{figure}[t]
\centering
\subfloat{\includegraphics[width=14cm]{Pics/FG}}
\caption{{Nucleotides model: One dimensional factor graph used in Belief Propagation.}}
\label{FG1}
\end{figure}

I apply BP to the factor graph of the test sequence with $n$ candidate nucleotides as the prior knowledge. BP updates the probability of candidate nucleotides based on the probabilities of their neighbors.


Then, the candidate index numbers are fed to a factor graph and the corresponding $D_{MSE}$ of each of candidates is employed to calculate the initial probability (prior probability) of each candidate.  Then, message passing (i.e., forward and backward) algorithm is applied to calculate the best match indices. The correspond subsequences of these indices is used in the next step.   

\subsection{Sequence Matching}

The sequence matching step is based on sparse coding and BP algorithm. In this step, I use the subsequences that were selected in the previous step to generate an over-complete dictionary. Then, for each nucleotide in the read sequence, I pick $n$ candidate nucleotides using sparse coding. 
To choose the best candidate, a factor graph is employed where for each nucleotide in the read sequence, a variable node was assigned. The relation between neighborhood variable nodes is taking care of by a factor node. Also one extra factor node for each variable node keeps the prior probability of candidates. 
By applying belief propagation to a factor graph, I can obtain the best match for each nucleotide in the read sequence. A detailed description about the sequence matching can be found in my  recent publications \cite{barzigar2012scobep,roozgard2013nucleotide}. A summary of the main procedure for my proposed indexing method is shown in Algorithm \ref{amin:mainIndex}.

\clearpage
\begin{algorithm} [t!]
\caption{Proposed nucleotide sequence alignment algorithm for estimating the location of the input sequence %$\mathcal{Z}$ 
}
\begin{algorithmic} 
\renewcommand{\algorithmicrequire}{\textbf{Inputs:}}
\REQUIRE {a reference sequence $\mathcal{X} \in \mathbb{R}^{M}$, a test sequence $\mathcal{Y} \in \mathbb{R}^{N}$, number of the candidate state matrix $k$, number of the candidate points $n$} 

\renewcommand{\algorithmicrequire}{\textbf{Initialize:}}
\REQUIRE {a $4 \times 4$  state matrix $I$ storing the numbers of nucleotide states (\ref{eq:ACGT}), nucleotide overlap $v$} 



\renewcommand{\algorithmicrequire}{\textbf{Iterate:}}
\REQUIRE {while the length of the sequence $\mathcal{X}$  is not reached}

\renewcommand{\algorithmicrequire}{\textbf{Fill the reference state matrix $I$:}}
\REQUIRE {For each subsequence $x_{i}\in \mathcal{X}$ with $v$ nucleotide overlap in each direction perform:}

\begin{itemize}
\item $I_i = MakeIndex(x_i)$
\end{itemize}


\renewcommand{\algorithmicrequire}{\textbf{Fill the test state matrix $J$:}}
\REQUIRE {For each subsequence $y_{j}\in \mathcal{Y}$ with $v$ nucleotide overlap in each direction perform:}

\begin{itemize}
\item $J_j = MakeIndex(y_j)$
\item $[c_j, {\rho}_j] = FindCandidates(J_j,I,k)$
\end{itemize}


\renewcommand{\algorithmicrequire}{\textbf{Refine the candidate state matrix:}}
\REQUIRE {\ }
\begin{itemize}
\item $\hat{{\rho}} = BP(c,{\rho})$
\end{itemize}

\renewcommand{\algorithmicrequire}{\textbf{Find the correspond nucleotide in the reference sequence $\mathcal{X}$:}}
\REQUIRE {For each subsequence $y_{j}\in \mathcal{Y}$ with $v$ nucleotide overlap in each direction perform:}

\begin{itemize}
\item $z_j= FindBestSubsequence(\mathcal{X},y_{j},\theta, n)$ (see Algorithm \ref{amin:alignment} %\cite{barzigar2012scobep,roozgard2013nucleotide} for more details
)
\end{itemize}



\renewcommand{\algorithmicrequire}{\textbf{Output:}}
\REQUIRE {the estimated version of aligned sequence $\mathcal{Z}$}
\end{algorithmic}
\label{amin:mainIndex}
\end{algorithm}


\clearpage
\newpage

% \subsection{\color{red}{reviewed till here}}

This step is inspired by my recent work, SCoBeP \cite{barzigar2012scobep}. First, I map the  reference nucleotide sequence $\mathcal{X}$ of size $N$ and the test nucleotide sequence $\mathcal{Y}$ of size $M$ into the two integer sequences, $\textbf{X}$ and $\textbf{Y}$, respectively. I then extract the features from the reference sequence $\textbf{X} \in \mathbb{Z}^{N\times 2}$ and the test sequence $\textbf{Y} \in \mathbb{Z}^{M \times 2}$. 

% {\color{red}{How to choose column size? a statistic analysis}}


Second, I create a dictionary $\mathcal{D}$ containing all the extracted feature vectors of the reference sequence $\textbf{X}$ which match the corresponding extracted features of the test sequence $\textbf{Y}$. The dictionary includes all vectorized one dimensional features as its columns where all of them have been normalized. I then apply sparse coding to each extracted feature of the test sequence. Sparse coding will reconstruct a nucleotide vector at each nucleotide $g_i$ as a linear combination of the reference sequences.
Note that the obtained representation coefficients {$\alpha_{i}$ should be sparse, i.e., it should be 0 for most coefficients. The non-zero coefficients of $\alpha_{i}$ indicate the corresponding nucleotides on the reference sequence (see Fig. \ref{SCpicture}).

\begin{figure}[t]
\centering
\subfloat{\includegraphics[width=10cm]{Pics/yDalpha}}
\caption{{Sparse representation of a feature vector $y_{i}$ with a dictionary $\mathcal{D}$ : $\hat{\alpha}_{i}$ as a sparse vector constructs the feature vector $y_{i}$ using a few columns (highlighted in gray) of dictionary $\mathcal{D}$.}}
\label{SCpicture}
\end{figure}



To select $n$ candidate nucleotides, I simply pick those corresponding to $n$ which have the largest coefficients in the sparse coefficient vector. I store the locations of these candidate nucleotides in a length-$n$ vector $\mathcal{L}_i$ and a probability vector $\rho_i$ as a length-$n$ vector which stores the corresponding probabilities for the sparse coefficient vector. Each coefficient in the probability vector $\rho_{i}$ serves as a prior probability of matching the nucleotide at $i$ to a nucleotide of a reference sequence. This probability vector takes only local characteristics into account but ignores neighborhood characteristics of the matches. 

Finally, I expect that nearby nucleotides in the test sequence should also match the nucleotides that are close to each other in the reference sequence. To incorporate these neighborhood characteristics, I model the problem by a factor  graph and apply the Belief Propagation (BP) algorithm to identify the best matches. I consider a one dimensional factor graph as  follows: for each nucleotide in the test sequence, one variable node was assigned and each variable node was connected  to its two neighbors by a factor node. Also, I consider one extra factor node for each variable node to impose the restriction of prior probabilities for the candidate nucleotides (see Fig. \ref{FG1}).

To align the test nucleotide sequence, I select the nucleotide candidate with highest probability and then calculate the displacement $\beta$ between the current nucleotide and the selected candidate nucleotide. Therefore, for each nucleotide in the test sequence, I have $\mathcal{Z}_i$, $\beta_i$ and $\rho_i$ for each nucleotide $g_i$ which are the most probable match nucleotide, the most probable displacement and the probability of the current match, respectively. A summary of the alignment method is shown in Algorithm \ref{amin:alignment}.



\begin{algorithm} [t!]
\caption{FindBestSubsequence($\mathcal{X}, y, \theta, n$)
% \caption{Proposed nucleotide sequence alignment algorithm - estimate the location of the input sequence %$\mathcal{Z}$ 
}
\begin{algorithmic} 
\renewcommand{\algorithmicrequire}{\textbf{Inputs:}}
\REQUIRE {a reference sequence $\mathcal{X} \in \mathbb{R}^{M}$, a test sequence $\mathcal{Y} \in \mathbb{R}^{N}$, a threshold $\theta$, number of the candidate points $n$} 


\renewcommand{\algorithmicrequire}{\textbf{Convert a string sequence to numeric sequence:}}
\REQUIRE {\ }

\begin{itemize}
\item $\textbf{X} = ConvertData(\mathcal{X})$
\item $\textbf{Y} = ConvertData(\mathcal{Y})$
\end{itemize}

\renewcommand{\algorithmicrequire}{\textbf{Extract feature and construct dictionary:}}
\REQUIRE {\ }


\begin{itemize}
\item $\hat Y = ExtractFeature (\textbf{Y})$
\item $\hat X = ExtractFeature (\textbf{X})$
\item $\mathcal{D} = MakeDic (\hat{X})$
\end{itemize}

\renewcommand{\algorithmicrequire}{\textbf{Find the initial estimate of the match location:}}
\REQUIRE {For each vector $y_{i}\in \hat{Y}$ perform:}

\begin{itemize}%[(a)]
\item $\alpha_{i} = FindSparseVector(\mathcal{D},y_{i})$ 
\item $[\mathcal{L}_{i}, \hat\rho_{i}]= FindTopScoreMatch(n,\alpha_{i})$
\end{itemize}


\renewcommand{\algorithmicrequire}{\textbf{Refine the candidate match location:}}
\REQUIRE {\ %For all vector $y_{i} \in \hat{Y}$ perform:
}
\begin{itemize}
\item ${\rho}= BP(\mathcal{L},\hat\rho)$
\end{itemize}


\renewcommand{\algorithmicrequire}{\textbf{Find the correspond nucleotides:}}
\REQUIRE {\ }
\begin{itemize}
\item $[\mathcal{Z},\beta]= Warp(\hat{X},{\rho},\mathcal{L})$
\end{itemize}


\renewcommand{\algorithmicrequire}{\textbf{Output:}}
\REQUIRE {the estimated version of aligned sequence $\mathcal{Z}$}
\end{algorithmic}
\label{amin:alignment}
\end{algorithm}


\clearpage


% 
% ------------------------------------------ bi havas --------------------------------------
% 






\subsection{Implementation Details} \label{nafiss:Notation2}
\begin{itemize}
\item $I_i = MakeIndex(x_i)$ fills the state matrix $I_i$ using the relationship of nucleotides in the subsequence $x_i$. The subsequence $x_i$ is scanned through all its nucleotides and the corresponding counts will be stored into the state matrix $I_i$. For example, $k_{cg}$ in $I_i$ in (\ref{eq:ACGT1}) shows how many times the nucleotide $C$ will be identified, which is next to the nucleotide $G$ in the subsequence $x_i$. Note that each subsequence $x_i$ has a separate state matrix $I_i$, where $i$ is the subsequence index.

 \begin{equation}
 % \text{State Matrix}
 I_{i} = \ \ \bordermatrix{ & A & C & G & T\cr
 			 A& k_{aa} & k_{ac} & k_{ag} & k_{at}\cr
 			 C& k_{ca} & k_{cc} & k_{cg} & k_{ct}\cr
 			 G& k_{ga} & k_{gc} & k_{gg} & k_{gt}\cr
 			 T& k_{ta} & k_{tc} & k_{tg} & k_{tt}}
 			 \times \frac{1}{\sum_{s,w \in \{a,c,g,t\}} k_{sw}} 
 \label{eq:ACGT1}
 \end{equation}

\item $[c_j, {\rho}_j ]= FindCandidates(J_j,I,k)$ identifies $k$ candidate state matrices that are highly similar to the test state matrix $J_j$ in $I$ and stores their indices in vector $c_j$ and their probabilities in vector ${\rho}_j$. Note that the approach will compute the Mean Square Error (MSE) of the test state matrix $J_j$ with each possible $I_i$ of the reference state matrices and select $I_i$ that has the smallest MSEs. 

\item $\hat{{\rho}}= BP(c,{\rho})$ models the problem by a factor graph and applies belief propagation \cite{kschischang2001factor} to update probability ${\rho}$. The updated probability $\hat{\rho}$  can be used to align the reference state matrix index onto the test state matrix index. In my case, I assign a variable node for each test state matrix index and connect each pair of neighboring state matrix indices with a factor node. Also, I introduce one extra factor node to take care of the prior knowledge obtained in the MSE step for each test state matrix index (for more details, see \cite{barzigar2012scobep}).



\item $\textbf{X} = ConvertData(\mathcal{X})$ maps the input string of the nucleotide sequence to a sequence of integer values for the further processing. The mapping function of four nucleotides $\{A,C,G,T\}$ can be defined in a two-dimension space where $A=(1,1)$, $C=(1,-1)$, $G=(-1,1)$, and $T=(-1,-1)$, respectively. Therefore, $\{A,C,G,T\} \mapsto\{(1,1),(1,-1),(-1,1),(-1,-1)\}$.

\item $\hat Y = ExtractFeature (\textbf{Y})$ presents  a vector extractor algorithm using $\textbf{Y}$ as a source sequence, where the result is a two dimensional matrix containing the vectorized one dimensional sequences. To this end, I consider a vector of size $S = 2\times (2a + 1)$ containing neighboring nucleotides on two sides of a nucleotide, where $a$ is a positive integer and the first ``2'' is corresponding to the dimension of the mapping space). For each nucleotide $g_{i}$ in the test sequence $\textbf{Y}$, I vectorized a sequence centered around the nucleotide $g_{i}$ to a feature vector $y_{i} \in \mathbb{Z}^{S \times 1}$. A two dimensional test feature data $\hat{Y} \in \mathbb{Z}^{N \times S}$ is then constructed from $y_{i}$ as follows:

\begin{equation}
\hat{Y} =  \{ y_{i} \ | \ 1\leq i\leq M \}.
\label{eq1}
\end{equation}
Note that $\hat{X}$ is created in the same manner as $\hat{Y}$ but from the reference sequence $\textbf{X}$ instead. 

% % % % % % % % % % % % % % % % % % % % 
\item  $\mathcal{D} = MakeDic (\hat{X})$ creates a dictionary $\mathcal{D}$ using the vectors of $\hat{X}$. Later, the dictionary $\mathcal{D}$ is used to match the extracted features of the source sequence to corresponding the extracted features of the reference sequence.
, a dictionary which contains feature vectors of $\textbf{X}$ is constructed. Thus, 
I can write $\mathcal{D}= \left [ x_{1}\ x_{2} ... \ x_{N}\right ]$
where $x_{i}$ is a feature vector of  $\hat{X}$. Note that I normalize dictionary $\mathcal{D}$ to guarantee the norm of each feature vector to be $1$.

% % % % % % % % % % % % % % % % % % % % 


\item $\alpha_{i} = FindSparseVector(\mathcal{D},y_{i})$ finds the candidate match nucleotide using the sparse coding algorithm, where $\alpha_{i}$ is a sparse vector.  Mathematically, I try to solve the following sparse coding problem to find the most sparse coefficient vector $\alpha_{i}$ (see Fig. \ref{SCpicture}) such that
\begin {equation}
y_{i} = \mathcal{D} \hat{\alpha}_{i}.
\label{eq3}
\end {equation}



Although there are several methods to solve (\ref{eq3}) \cite{yang2010fast,dai2009subspace,maleh2007sparse},  in my work, I employ Subspace Pursuit (SP) \cite{dai2009subspace} because of its computational efficiency.



% % % % % % % % % % % % % % % % % % % % % %  
\item $[\mathcal{L}_{i},\hat\rho_{i}] = FindTopScoreMatch(n,{\alpha}_{i})$ picks up the $n$ largest coefficients of $\alpha_{i}$ as $n$ candidates. $\mathcal{L}_{i}$ is a $n \times 1$ vector that stores the locations of these candidate nucleotide and $\hat\rho_{i}$ is the length-n vector that stores the corresponding probabilities of $\mathcal{L}_{i}$. Each coefficient in $\hat\rho_{i}$  serves as a prior probability of matching the source sequence at $i$ to a sub-sequence centered around the nucleotide $g_i$.
After finding the candidates and their initial probabilities, I concatenate the result of each nucleotide and construct following matrices:
\begin{equation}
\mathcal{L} = \left [ \mathcal{L}_{1}\ \mathcal{L}_{2} ... \ \mathcal{L}_{N}\right ], \hat{\rho} = \left [ \hat{\rho}_{1}\ \hat{\rho}_{2} ... \ \hat{\rho}_{N}\right ], 
\end{equation}
which I will use to apply belief propagation at the next step.


% % % % % % % % % % % % % % % % % % % 

\item $\rho= BP(\mathcal{L},\hat\rho)$ models the problem by a factor graph (see Fig. \ref{FG}) and applies belief propagation \cite{kschischang2001factor} to update probability ${\rho}$. The updated probability $\rho$  can be used to align the reference sequence onto the test sequence. In my case, I assign a variable node for each nucleotide on the source sequence and connect each pair of neighboring nucleotide with a factor node. Also, I introduce one extra factor node to take care of the prior knowledge obtained in the sparse coding step for each nucleotide of the source sequence (for more details, see \cite{barzigar2012scobep}).


\item $[\mathcal{Z},\beta]= Warp(\textbf{X},\rho,\mathcal{L},\theta )$ returns the aligned sequence $\mathcal{Z}$ and a displacement vector $\beta$ which contains the movement of each nucleotide in the reference sequence. Note that in step $BP(\cdot)$, I calculated the refined probability of each candidate nucleotide match.% This matrix can be used to refine the result of the voxels with a probability less than $\theta$.




% % % % % % % % % % % % % % % 


\item $z_j=FindBestSubsequence(\mathcal{X},y_{j},\theta, n)$ finds the corresponding location for a nucleotide $y_{j} \in \mathcal{Y}$. In this step, the reference nucleotide sequence $\mathcal{X}$ and the test nucleotide sequence $\mathcal{Y}$ are converted into two integer sequences. Then, an over-complete dictionary is built with all subsequences in the $\mathcal{X}$. I then apply sparse coding followed by using Belief Propagation (BP) to identify the best matches. (see \cite{barzigar2012scobep, roozgard2013nucleotide} for more details) Note that I used non-overlapped subsequences to build the dictionary. This change decreases the memory usage and the accuracy of the proposed algorithm in compare to 1D-SCoBeP \cite{roozgard2013nucleotide}, but it increases the speed of my alignment algorithm. 
\end{itemize}




\section{Experimental Results}\label{EXPERIMENTAL}

\begin{landscape}
\begin{figure}[t]
\centering
\subfloat[]{\includegraphics[width=10cm]{Pics/nFig1}}
\subfloat[]{\includegraphics[width=10cm]{Pics/nFig1z}}
\caption{{The results of proposed method for non-collinear nucleotide sequence alignment. a) comparison among alignment results of the ground truth, 1D-SCoBeP \cite{roozgard2013nucleotide} and the proposed method. {\color{black}{b) zoomed of the black square in figure \ref{FG}--a to show the gap between the proposed method, ground truth and 1D-SCoBeP \cite{roozgard2013nucleotide} on the jump point. The x--axis and y--axis are the index numbers of the original genome sequences and the shuffled genome sequences, respectively.}}}}
\label{FG}
\end{figure}

\end{landscape}

\begin{figure}[t]
\centering
\subfloat{\includegraphics[width=15cm]{Pics/indels_error_report}}
\caption{Accuracy of BWA \cite{li2009fast}, SOAP aligner \cite{li2009soap2} and the proposed method in present of different Indel rate, where the testing genome sequences were obtained from \cite{website:NCBI}.}
\label{Indel_error1}
\end{figure}


\begin{figure}[t]
\centering
\subfloat{\includegraphics[width=14cm]{Pics/mixed_report}}
\caption{The percentage of successful alignments in present of 0.5\% to 1.5\% indels. {\color{black}{The \textbf{green line} is the percentage successful alignments where the rate of the indels are changing with step equal 0.05\% between 0.7\% and 1.3\%. The \textbf{blue line} is the percentage successful alignments where the rate of the indels are changing with step equal 0.1\% between 0.5\% and 1.5\% and the \textbf{red line} is the same as the Fig. \ref{Indel_error1}.}} Each point represents $10^5$ random site selection with same indels rate. Note that the genome sequences used in this studies were obtained from \cite{website:CGHub} and  \cite{website:NCBI}.}
\label{Indel_error_Mixed}
\end{figure}


%{\color{blue}{relation between factor function in the BP and number of candidates: if the factor function is sharp, the required number of candidates is less. }}

%{\color{red}{Based on \cite{brudno2003glocal} lease divide the experimental results into: \newline 1) No rearrangements(only having insertion or deletion) \newline 2) Trans-location (These first two parts results are already included in the manuscript. still I can more expand the part one.) \newline 3) Duplication \newline 4) Inversion (if it is working, if not I have to find a solution and change the proposed method to make it possible) }}

%{\color{red}{Insertion and deletion test on long and short sequences (based on work in \cite{darling2010progressivemauve} )}}

I designed my experiments based on work in \cite{darling2010progressivemauve} to evaluate the proposed method for aligning the nucleotide sequences and to compare it with SOAP aligner \cite{li2009soap2}, BWA \cite{li2009fast} and 1D-SCoBeP \cite{roozgard2013nucleotide}. I considered the problem of aligning a sequence of human nucleotides from the National Center for Biotechnology Information (NCBI)\cite{website:NCBI} and Cancer Genomics Hub (CGHub) \cite{website:CGHub}.  

To evaluate the performance of my approach, I conducted two sets of tests on the nucleotide sequences. In the first set, I selected fifty short sub-sequences of human genomes and then used SOAP aligner, BWA, 1D-SCoBeP and the proposed method to find the location of selected sub-sequence nucleotide in the human chromosome. All of four algorithms successfully passed this test. I created twenty shuffled sub-sequences of the reference sequence as follows: for each read sequence $R$, I cut it into five pieces $p_1$, $p_2$, $p_3$, $p_4$ and $p_5$. Then I switched $p_2$ with $p_4$. Therefore, I converted a read sequence $R=[p_1,p_2,p_3,p_4,p_5]$ into a new read sequence $\hat{R}=[p_1,p_4,p_3,p_2,p_5]$.

Fig. \ref{FG} shows the result of the 1D-SCoBeP and the proposed method show a better performance with a gap  of 100 to 120 nucleotides away from the ground truth. Since I was using non-overlapped subsequences for the dictionary generation, the gap between the proposed method and the ground truth was larger than these reported in 1D-SCoBeP \cite{roozgard2013nucleotide}. In my experiments, the following parameters were used: the number of candidate points $n$ is set to be $3$, the sparsity factor $k = 3$ and the dictionary column size $a=200$.

To evaluate the robustness of the proposed method, I generate indices for long human genome sequences (i.e.  $5 \times 10^8$ nucleotides) where $h=10000$ and $o=5000$. Moreover, I synthesized insertion, deletion and mutation (i.e., indel) in these sequences. For indel rate, I picked $10^5$ number of subsequences with size of $10^4$ nucleotides. Then, I randomly modified a certain number of nucleotides (based on the indel rate) and aligned them with the references.  I counted how many times the alignment location and real subsequence location (i.e., ground truth) are matched, where the accuracy is defined as the count of the successfully aligned sequences over total number of the subsequences. Fig. \ref{Indel_error1} shows the accuracy of alignment of the proposed method, BWA and SOAP aligner in the presence of the different indel rates. The proposed method showed similar accuracies even when I increased the indel rate to 3\%. Moreover, the proposed algorithm still showed more than 75\% accuracy even after I modified 5\% of the nucleotides in my selected subsequences. In contrast, the accuracy of the BWA and SOAP aligner decreased sharply as the indel rates increase.   

%In addition, based on the fact that the SNPs exist in only 1\% of the human DNA, 
I investigate the impact of small indel rate in the range from 0.5\% to 1.5\% in Fig. \ref{Indel_error_Mixed}. % over three different datasets. 
In this figure, I showed accuracy of 1\% indels in red for the data set used in \ref{Indel_error1} as reference. {\color{black} To verify my result, I repeat the experiments with different indel steps and different read locations and present} the results in green and blue, respectively. Note that each point in this figure was obtained from the evaluation over $10^5$ read sequences. {\color{black} There are slight variation among the curves due to statistical deviation.} The summary of the indel rate accuracy was shown in Table \ref{Table1}.

The computational complexity of proposed is mainly determined by the following three steps: 1) indexing 2) index matching 3) extracting sub-sequence nucleotides as features and constructing the dictionary, 4) finding candidate nucleotides via sparse coding, and 5) applying BP. Assume the size of the read and reference sequences are $N$ and $M$ nucleotides, respectively. The required time for create indexes is $O\big(M+N\big)$, because I have to scan whole read and reference sequences. The number of reference sequence indexes is $I_M=\frac{M}{h}+\frac{2oM}{h^2}=O\big(\frac{M}{h}\big)$ and similarly, $I_N=O\big(\frac{N}{h}\big)$. Therefore, the time for the index match matching is $O\big(\frac{M}{h}\big)\times O\big(\frac{N}{h}\big)=O\big(\frac{MN}{h^2}\big)$. After index matching step, the size of search space reduces from $M$ to $\bar{M}=I_s \times h$ where $I_s$ is the number selected indexes $I_s$ and $h$ is the size of each index. The required time of feature extraction will be $O\big(a(\bar M+N)\big)$, where $a$ is the size of the vector of extracted features for each nucleotide. The dictionary construction step involves the normalization of each column, which requires  $O(a\bar M)$ amount of time. Thus the total time complexity of the first step is $O\big(a(\bar M+N)\big)$. 
In the next step, the time complexity of Subspace Pursuit (SP) is $O\big(\log(f)a\bar M\big)$ \cite{dai2008subspace}, where $f$ is the number of iterations for searching the sparse vector. Since I have to repeat the process to find candidate points for all $N$ feature vectors, the time complexity of finding candidate points by SP is $O\big(\log(f)a \bar MN\big)$. Then, the time complexity of Belief Propagation in my factor graph is $O(vn^2\bar M)$, where $v$ is the number iterations before converging and $n$ is the number of candidates in each variable node. Finally, the time complexity of proposed method will be $O\Big( MN  + \log(f)a\bar MN + vn^2\bar M  \Big)$. 

\begin{table}[t]
\caption{ {Percentage of successfull alignments}}
\centering

\begin{tabular}{ |p{1.8cm}|p{1.5cm}|p{1.5cm}|p{1.5cm}| }
\hline
\% of the Indels  & Accuracy of \textbf{red} line & Accuracy of \textbf{blue} line & Accuracy of \textbf{green} line \\ \hline
 0.00 & 80.33	& --	& --\\
 0.50 & --	& 81.19	& --\\
 0.60 & --	& 79.38	& --\\
 0.70 & --	& 80.90	& 80.64\\
 0.75 & --	& --	& 80.63\\
 0.80 & --	& 79.82	& 80.29 \\
 0.85 & --	& --	& 78.90\\
 0.90 & --	& 78.09	& 80.63\\
 0.95 & --	& --	& 78.74\\
 1.00 & 79.85	& 79.71	& 78.04\\
 1.05 & --	& --	& 80.43\\
 1.10 & --	& 80.49	& 79.70\\
 1.15 & --	& --	& 78.22\\
 1.20 & --	& 79.81	& 79.78\\
 1.25 & --	& --	& 79.16\\
 1.30 & --	& 80.54	& 78.94\\
 1.40 & --	& 80.70	& --\\
 1.50 & --	& 79.52	& --\\
 2.00 & 79.09	& --	& --\\
 3.00 & 78.90	& --	& --\\
 4.00 & 76.33	& --	& --\\
 5.00 & 75.90	& --	& --\\
 6.00 & 72.09	& --	& --\\
 7.00 & 72.86	& --	& --\\
 8.00 & 69.79	& --	& --\\
 9.00 & 67.87	& --	& --\\
 10.00 & 66.42 	& --	& --\\
\hline
\end{tabular}
\label{Table1}
\end{table}


\section{Conclusion}\label{conclusion}
In this chapter, I proposed a sparse coding and BP based method for indexing and alignment genome sequences. The proposed method builds a transition matrix based on the neighboring nucleotides of an input sequence and then reduces the search space by selecting the top $K$ most similar subsequences based on their distances. The proposed algorithm selects candidate nucleotides by using sparse coding with an over-completed dictionary, which was constructed from the nucleotides of reference sequence in the indexing step. BP algorithm is then applied to select the best matches. Through experimental results, I showed that the proposed algorithm are comparable to SOAP aligner \cite{li2009soap2},  BWA \cite{li2009fast} and 1D-SCoBeP \cite{roozgard2013nucleotide} in terms of the alignment accuracy. In addition, the proposed method is robust to insertions, deletions, and mutations in the genome sequences when comparing with SOAP aligner and BWA. Finally, the proposed method is able to process much longer sequences then our previous 1D-SCoBeP approach.




\chapter{GENOME SEQUENCE PRIVACY PROTECTION USING COMPRESSED SENSING}

In this chapter, I present a privacy preserving genomic data dissemination algorithm based on the compressed sensing. I participated in the challenge at the iDASH on March 24, 2014 in La Jolla, California which the result of the challage are available on $http://www.humangenomeprivacy.org$. In my proposed method, I am adding the noise into the sparse representation of the input vector to make it differentially private. It means I find the sparse representation using using the SubSpace Pursuit and then disturb it with sufficient Laplasian noise. I compare my method with state-of-the-art compressed sensing privacy protection method  \cite{li2011compressive}. 



\section{Introduction}

Most genomic datasets are not publicly accessible, due to privacy concerns. Patients' genomic data contains identifiable markers and can be used to determine the presence of an individual in a dataset. Prior research shows that the re-identification can be possible when: There is a very small set of SNPs, Or The genomic data is aggregated; like releasing the frequencies of different SNPs across a population. To protect patients, the data owners impose an application and evaluation procedure. Then an agreement (like IRB approval) needs to be signed before the use of data is permitted. This process often takes months to complete and limits the researchers.

One solution to the problem can be to let each data owner publish a set of pilot data to help data users choose the right datasets based on their needs. Such pilot data comes from adding noise to the original genomic data to ensure that individuals’ information is protected. The data owners release these pilot data with the noise parameters and the mechanism that they used. A data user can download, run any kind of association tests and compare the outcomes with the other datasets outputs to get an idea which datasets can be useful. With such information, the researchers can approach the owners of the most relevant datasets for further research with proper agreements.

I reviewed the basics of compressed sensing in section \ref{CompressedSensing} that is the inspiration of current research. Now, I explain my proposed genome privacy protection mechanism and conclude by demonstrating the results of proposed method on human genome sequences.

\section{Privacy Protection Based On Compressed Sensing}
The proposed privacy protection method is based on the compressed sensing mechanism. The input sequences to the proposed method are the genomic nucleotide ($A$, $C$, $G$, $T$). First, I process the input genomic sequences to find the location of the Single Nucleotide Polymorphisms (SNPs). I check the same location of the all genome sequences and look for the locations that their nucleotides are different which I mark them as SNPs. Then I keep only the SNPs location as representatives of the original genome sequences and calculate the frequency of changes in the same SNP. The nucleotide which appears less in the SNP location called ``minor'' and the nucleotide with most appearance in the SNP location called ``major''. These majors and minors with their frequencies will be used in the proposed privacy protection method.  

To have a differentially private data, I am using the compressed sensing technique and adding the noise to the genomic data representation in the transformation domain. As it was explained in the section \ref{CompressedSensing}, I need to specify an input vector $y$, a random matrix $\Phi$ and a transformation matrix $\Psi$ to be used in the compressed sensing process. The input vector $y$ is the frequency of the minors SNPs which I calculate it from pre-processing step. The random matrix $\phi$ is a binary independent and identically distributed (i.i.d.) of the Bernoulli distribution matrix which in average, half of the coefficient in each row and each column are ``$1$''s and the rest are equal to ``$0$''s. The transformation that I choose is the Haar wavelet transform (HWT) \cite{heil1989continuous,mulcahy1997image} Because I know from \cite{xiao2011differential,li2011compressive,cormode2011differentially} that if I have an $\epsilon$ budget how to use it on the wavelet domain that keeps the differential privacy.


In my proposed method, I am adding the noise into the sparse representation of the input vector to make it differentially private. It means I calculate the matrix $A=\Phi  \Psi$ and then using the SubSpace Pursuit \cite{dai2009subspace}, I reach to the sparse representation of the $y$ which I will call it $\hat{x}$. To making sure that the amount of the noise that I am adding to sparse representation is enough to protect the differential privacy, I add Laplasian noise $\lambda $ to all element of the sparse representation according to \cite{cormode2011differentially}. Note that the element in the sparse representation are correspond to HWT tree (see Figure \ref{W1}), therefore based on \cite{xiao2011differential}, I double the amount of noise when I move up one step in the wavelet tree from the leaves to the root. The result of this result of the adding Laplasian noise to $\hat{x}$ generates $x^*$ which I use it to generate the differential private data to publish. To get the punishable data $y^*$, I need to apply inverse transformation on $x^*$ as follows: 
\begin{equation}
 y^*=\Psi x^*
\end{equation}
where the $\Psi$ is the Haar wavelet transformation matrix. Algorithm \ref{amin:mainAlgorithm}. shows the summary of the proposed method.


\begin{algorithm} [t]
\caption{: Proposed differentially private Protection using Compressed Sensing algorithm for Genomic Data %$\mathcal{Z}$ 
}
\begin{algorithmic} 
\renewcommand{\algorithmicrequire}{\textbf{Inputs:}}
\REQUIRE {a set genome sequences $\mathcal{X}$ } 

\renewcommand{\algorithmicrequire}{\textbf{Pre-processing:}}
\REQUIRE {Converting the genome sequence to SNP sequences} 



% \renewcommand{\algorithmicrequire}{\textbf{Iterate:}}
% \REQUIRE {while the length of the sequence $\mathcal{X}$  is not reached}

\renewcommand{\algorithmicrequire}{\textbf{Initialization:}}
\REQUIRE {Set the initial parameters:}

\begin{itemize}
\item Find the frequencies of major and minor for each SNP location
\item Generate sampling matrix $\Phi$ from a Bernoulli i.i.d. distribution
\item Generate Haar wavelet transformation matrix $\Psi$ 
\item Calculate matrix $A=\Phi\Psi$ 
\item Consider an array contains the frequencies of minors as input vector $y[c_j, {\rho}_j] = FindCandidates(J_j,I,k)$
\end{itemize}


\renewcommand{\algorithmicrequire}{\textbf{Adding Noise:}}
\REQUIRE { Add Laplasian noise:}

\begin{itemize}
\item Finding $\hat{x}$, a sparse representation of the input vector $y$ using SSP \cite{dai2009subspace} where $y=A\hat{x}$  
\item Add sufficient Laplasian noise according to \cite{xiao2011differential,li2011compressive,cormode2011differentially}
\end{itemize}


\renewcommand{\algorithmicrequire}{\textbf{Post-processing:}}
\REQUIRE {\ }
\begin{itemize}
\item $y^*=\Psi x^*$
\end{itemize}

\renewcommand{\algorithmicrequire}{\textbf{Output:}}
\REQUIRE {the noise-added differentially private version of SNP frequencies $y^*$}
\end{algorithmic}
\label{amin:mainAlgorithm}
\end{algorithm}

\clearpage


\begin{figure}[h!] 
\centering
\includegraphics[width=17cm]{Pics/W1.pdf}
\caption{The Haar wavelet tree structure of  [ $(0,0)$, $(1,0)$, $(1,1)$, $(2,0)$, $(2,1)$, $(2,2)$, $(2,3)$, $(3,0)$, $(3,1)$, $(3,2)$, $(3,3)$, $(3,4)$, $(3,5)$, $(3,6)$, $(3,7)$ ]. Note that each $(x,y)$ shows one entry in the array. The node $(0,0)$ is the root of the wavelet  transformation tree and the nodes $(3,t)$, where $t$ is a number between $0$ and $7$, are the leaves of the wavelet transformation tree.  }
\label{W1}
\end{figure}
 


  

\section{Experimental Results}

To evaluate the proposed method, I participated in the challenge at the iDASH (UCSD--based National Center for Biomedical Computing: NIH U54HL108460) workshop \footnote{http://www.humangenomeprivacy.org} on March 24, 2014 in La Jolla, California and the results are available online at \cite{WIDGET}.
The challenge was about sharing aggregate human genomic data (i.e., allele frequencies) to preserve the privacy of the data and to maximize the utility of the data  for Genome--Wide Association Studies (GWAS). I applied the proposed method on the two datasets from the case and control groups of individuals: the case group includes 411 individuals from the Personal Genome Project \cite{church2005personal} \footnote{http://www.personalgenomes.org/}, and the control data includes 174 participants from the CEU population in HapMap \cite{gibbs2003international} \footnote{http://hapmap.ncbi.nlm.nih.gov/index.html.en}. There are two test sets: the first one consists of 311 SNP sites of the human chromosome 2 and The second one consists of 600 SNP sites of human chromosome 10. Table \ref{table1} shows the result of proposed method on these two test sets in compare to the SNP--Based baseline of the challenge organizers.
\begin{landscape}
\begin{table}[ht]
  
\centering 
\begin{tabular}{|c| c|| c| c|| c| c|| c|} 

\hline

\multicolumn{2}{|c||}{Data Set}&\multicolumn{2}{|c||}{SNP--Based baseline}&\multicolumn{2}{|c|
|}{Proposed Method} & Number of the significant SNPs\\
\hline
\hline
\multirow{5}{*}{Dataset 1}& Power &  \multicolumn{2}{c||}{0.05} & \multicolumn{2}{c||}{0.61} & \\ \cline{2-6}
& Cutoffs & TPR & FPR & TPR & FPR & \\ \cline{2-7}
& $5 \times 10^{-2}$ & 0.864 & 0.844 & 1.0 & 0.941 & 22 \\
& $10^{-3}$ & 0.632 & 0.774 & 1.0 & 0.884 & 19 \\
& $10^{-5}$ & 0.642 & 0.700 & 1.0 & 0.879 & 14\\
\hline
\hline
\multirow{5}{*}{Dataset 2}& Power &  \multicolumn{2}{c||}{0.4} & \multicolumn{2}{c||}{0.005} & \\  \cline{2-6}
& Cutoffs & TPR & FPR & TPR & FPR &  \\ \cline{2-7}
& $5 \times 10^{-2}$ & 0.933 & 0.924& 1.0 & 0.958 & 45 \\
& $10^{-3}$ & 0.800 & 0.862 & 1.0 & 0.909& 15 \\
& $10^{-5}$ & 0.625 & 0.788 & 1.0 & 0.876& 8\\
\hline %inserts single line
\end{tabular}
\caption{Results of the proposed method on the challenge datasets. The Dataset 1 refers to 200 participants with 311 SNPs on chromosome 2 and Dataset 2 refers to 200 participants with 610 SNPs on chromosome 10. The ''Power`` row is the ratio of identifiable individuals using the likelihood ratio test in the case group. The false positive rate (FPR) and true positive rate based on $\chi^2$ test are listed per different cutoff threshold. In addition, the last column corresponds to the number of significant SNPs.} 
\label{table1} % is used to refer this table in the text
\end{table}
\end{landscape}




% A statistic  $\bar{\mathcal{L}}$   is then computed over all SNP sites for a victim $i$, which represents the likelihood of her presence in the case group. The probability of a victim in the case group to be re--identified is high if the LR test statistic from her SNPs is significantly higher than those who are not in the case group (e.g., from a test group that does not overlap with the case or reference groups). 

%I used the LR test to evaluate the privacy risks in the submitted data after perturbation from participating teams: the submitted data were considered anonymized if there are about the same number of individuals in the case group as those in a test group (i.e., individuals not in the case or control groups) with LR statistics above the same threshold. The utility if the perturbed data were measured by using the number of significant SNPs identified by the $\chi^2$ association test \cite{chernoff1954use} over the case and control groups. A data perturbation method was considered to have high utility if a majority of significant (above a p--value threshold, e.g., $10^{-5}$ ) SNPs reported by the $\chi^2$--association test on perturbed data were true, i.e., they are identified as truly significant on the un--perturbed data by the same statistic test, and a majority of truly significant SNPs remained to be identified on the perturbed data.

Figures \ref{R1} and \ref{R2} are the result of the proposed method disturbed data for identification of an individual in the data base. Figures \ref{R3} and \ref{R4} show the utility evaluation of the proposed method results on the online evaluation tool WIDGET \cite{WIDGET}.

\begin{landscape}
\begin{figure}[h!] 
\centering
\includegraphics[width=20cm]{Pics/CH2_01.pdf}
\caption{ The online privacy evaluation \cite{WIDGET} of the proposed method on chromosome 2 with p-value 0.01  }
\label{R1}
\end{figure}
\end{landscape}

\begin{landscape}
\begin{figure}[h!] 
\centering
\includegraphics[width=20cm]{Pics/CH10_01.pdf}
\caption{The online privacy evaluation \cite{WIDGET} of the proposed method on chromosome 10 with p-value 0.01}
\label{R2}
\end{figure}
\end{landscape}


\begin{figure}[h!] 
\centering
\subfloat{\includegraphics[width=13cm]{Pics/dig_ch2_01.pdf}}

\subfloat{\includegraphics[width=13cm]{Pics/data_ch2_01.pdf}}
\caption{The online utility evaluation \cite{WIDGET} of the proposed method on chromosome 2 with p-value 0.01 }
\label{R3}
\end{figure}

\begin{figure}[h!] 
\centering
\subfloat{\includegraphics[width=13cm]{Pics/dig_ch10_01.pdf}}

\subfloat{\includegraphics[width=13cm]{Pics/data_ch10_01.pdf}}
\caption{ The online utility evaluation \cite{WIDGET} of the proposed method on chromosome 10 with p-value 0.01}
\label{R4}
\end{figure}

% \begin{figure*}[h!] 
% \centering
% \includegraphics[width=16cm]{Pics/CH10_01.pdf}
% \caption{ }
% \label{R4}
% \end{figure*}

The power of the likelihood ratio test is a privacy risk measure which the number of case individuals who can be recognized with confidence level more than a threshold. The lower likelihood ratio power level shows that the perturbed data has less risk of reidentification. In addition the utility of the privacy protection method can be measured based on the $\chi^2$ test to detect the significants SNPs with different cutoff p--value. As Table \ref{table1} shows the proposed method has the higher true positive significant detected on the perturbed data.

In addition to the challenge, I select 180 SNPs of the Personal Genome Project \cite{church2005personal} and partition it into two subsets. The first subset used as control group. I modified the frequencies of randomly selected SNPs on second subset and generated 3 new datasets with different levels of utility called A,B,C with 27,9 and 4 significant. After adding noise, the p-value of a SNP reported by a association test, $\chi^2$, of these three datasets can be compared with each other to select the best dataset as ``Best Pick'' which is A. Also, I check the results to find the ``Correct Order'' of datasets A better than B better than C. The results are shown in Table \ref{table2}. The numbers in the table are the percentage of successful detection of best dataset and the correct order of them. As this results shows the proposed method dominates the other methods.



\begin{table}[ht]
  
\centering 
\begin{tabular}{|c|| c| c|} 

\hline
&Correct Order \% &Best Pick \%\\
\hline\hline
SNP--Based baseline&19.34&34.54\\
\hline
Reference \cite{li2011compressive}&24.85&47.52\\
\hline
Proposed Method&25.02&75.12\\



\hline 
\end{tabular}
\caption{Results of the SNP--Based baseline, Reference \cite{li2011compressive} and the proposed method for best pick and correct order.} 
\label{table2} % is used to refer this table in the text
\end{table}



\section{Conclusion}


In this chapter, I present a privacy preserving genomic data dissemination algorithm based on the compressed sensing. The input sequences to the proposed method are the genomic nucleotide and I process the input genomic sequences to find the location of the SNPs. Then I specify an input vector of the SNPs frequencies, a random matrix and the wavelet transformation matrix. To making sure that the amount of the noise that I am adding to sparse representation is enough to protect the differential privacy, I add Laplasian noise $\lambda $ to all element of the sparse representation according to \cite{cormode2011differentially}. My experimental results shows that the proposed method outperform the other methods.  


\chapter{CONCLUSION}

In this dissemination, I propose a novel nucleotide sequence Indexing and alignment method based on empirical transitional probability, sparse coding and belief propagation to compare the similarity of the nucleotide sequences inspired by my recent works, 3D-SCoBeP described in chapter \ref{3Dmed}. The proposed method builds a transition matrix based on the neighboring nucleotides of an input sequence and then reduces the search space by selecting the top $K$ most similar subsequences based on their distances. The proposed algorithm selects candidate nucleotides by using sparse coding with an over-completed dictionary, which was constructed from the nucleotides of reference sequence in the indexing step. BP algorithm is then applied to select the best matches. The proposed method is robust to insertions, deletions, and mutations in the genome sequences when comparing with SOAP aligner and BWA. Finally, the proposed method is able to process much longer sequences then our previous 1D-SCoBeP approach. Through experimental results, I showed that the proposed algorithm are comparable to SOAP aligner \cite{li2009soap2},  BWA \cite{li2009fast} and 1D-SCoBeP \cite{roozgard2013nucleotide} in terms of the alignment accuracy. In addition, I present a privacy preserving genomic data dissemination algorithm based on the compressed sensing. The input sequences to the proposed method are the genomic nucleotide and I process the input genomic sequences to find the location of the SNPs. Then I specify an input vector of the SNPs frequencies, a random matrix and the wavelet transformation matrix. Then I add Laplasian noise $\lambda $ to all element of the sparse representation. My experimental results shows that the proposed method outperform the other methods.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% BIBLIOGRAPHY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

{
\singlespacing
% \bibliographystyle{plain}
\bibliographystyle{IEEEtran}
\bibliography{refs}
}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DEDICATION PAGE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\newpage
\newpage
\pagenumbering{roman}

\thispagestyle{empty}
{\singlespacing
\begin{center}
{
\par
\vspace{1.2in}
%\large
DEDICATION
\par
\vspace{0.17in}
to
\par
\vspace{1.2in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
My father and mother
\par
%\vspace{0.17in}
Mohammad Rouzgard and Sorayya Moradi
\par
\vspace{0.6in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
AND % Esperanza Pedroza Paipa
\par
\vspace{0.6in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
 My wife
\par
%\vspace{0.17in}
Nafise Barzigar
\par
\vspace{0.17in}
% Esperanza Pedroza Paipa
\par
\vspace{1.2in}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%For
\par
\vspace{0.17in}
%ncouraging me to follow my dreams
\par
\vfill
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
}
\end{center}
}

\end{document}
