\documentclass[11pt,a4paper]{article}
\usepackage{Sweave}
\usepackage{url}
\usepackage{a4wide}
\usepackage{amsmath,epsfig}
\usepackage{graphicx}
\usepackage[authoryear,round]{natbib}
\bibliographystyle{plainnat}

%% \VignetteIndexEntry{An overview of gap}
%% \Vignettekeywords{human genetics, linkage analysis, association analysis}
%% \VignettePackage{gap}

% \newcommand{\keywords}[1]{\addvspace{\baselineskip}\noindent{{\bf Keywords.} #1}}

\parindent 0in

\begin{document}
\SweaveOpts{engine=R,echo=FALSE,pdf=TRUE}
\title{A Genetic Analysis Package with R}
\author{Jing Hua Zhao}
\date{}
\maketitle

\begin{center}
MRC Epidemiology Unit, Cambridge, UK \\
\url{http://www.mrc-epid.cam.ac.uk}
\end{center}

\tableofcontents

% library(tools)
% Rnwfile<- file.path("/gap/inst/doc","gap.Rnw")
% Sweave(Rnwfile,pdf=TRUE,eps=TRUE,stylepath=TRUE,driver=RweaveLatex())

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

\section{Introduction}

This package was initiated to integrate some C/Fortran/SAS programs I have written 
or used over the years. As such, it would rather be a long-term project, but an
immediate benefit would be something complementary to other packages currently
available from CRAN, e.g. {\bf genetics}, {\bf hwde}, etc. I hope eventually this will
be part of a bigger effort to fulfill most of the requirements foreseen by many, 
e.g. ~\cite{guo00}, within the portable environment of R for data management, 
analysis, graphics and object-oriented programming. My view has been outlined
more formally in ~\cite{zt06a} and ~\cite{zt06b} in relation to other package
systems. Also reported are ~\cite{zhao05} and ~\cite{zhao06} on package {\bf kinship}.\\

The number of functions are quite limited and experimental, but I already feel
the enormous advantage by shifting to R and would like sooner rather than later
to share my work with others. I will not claim this workas exclusively done by
me, but would like to invite others to join me and enlarge the collections and
improve them.

\section{Implementation}

The following list shows the data and functions currently available.

\begin{tabular}{ll}
\\
BFDP              &      Bayesian false-discovery probability\\
FPRP              &      False-positive report probability\\
SNP               &      Functions for single nucleotide polymorphisms (SNPs)\\
ab                &      Test/Power calculation for mediating effect\\
aldh2             &      ALDH2 markers and alcoholism\\
apoeapoc          &      APOE/APOC1 markers and schizophrenia\\
asplot            &      Regional association plot\\
bt                &      Bradley-Terry model for contingency table\\
b2r               &      Obtain correlation coefficients and their variance-covariances\\
ccsize            &      Power and sample size for case-cohort design\\
chow.test         &      Chow's test for heterogeneity in two regressions\\
cf                &      Cystic Fibrosis data\\
comp.score        &      score statistics for testing genetic linkage of quantitative trait\\
crohn             &      Crohn's disease data\\
ESplot            &      Effect-size plot\\
fa                &      Friedreich ataxia data\\
fbsize            &      Sample size for family-based linkage and association design\\
fsnps             &      A case-control data involving four SNPs with missing genotype\\
gc.em             &      Gene counting for haplotype analysis\\
gcontrol          &      genomic control\\
gcontrol2         &      genomic control based on p values\\
gcp               &      Permutation tests using GENECOUNTING\\
genecounting      &      Gene counting for haplotype analysis\\
gif               &      Kinship coefficient and genetic index of familiality\\
hap               &      Haplotype reconstruction\\
hap.em            &      Gene counting for haplotype analysis\\
hap.score         &      Score statistics for association of traits with haplotypes\\
hla               &      HLA markers and schizophrenia\\
htr               &      Haplotype trend regression\\
h2                &      Heritability estimation according to twin correlations\\
hwe               &      Hardy-Weinberg equilibrium test for a multiallelic marker\\
hwe.cc            &      A likelihood ratio test of population Hardy-Weinberg equilibrium \\
                  &      for case-control studies\\
hwe.hardy         &      Hardy-Weinberg equilibrium test using MCMC\\
kin.morgan        &      kinship matrix for simple pedigree\\
klem              &      Haplotype frequency estimation based on a genotype table\\
                  &      of two multiallelic markers\\
LD22              &      LD statistics for two diallelic markers\\
LDkl              &      LD statistics for two multiallelic markers\\
lukas             &      An example pedigree\\
\end{tabular}

\begin{tabular}{ll}
makeped           &      A function to prepare pedigrees in post-MAKEPED format\\
mao               &      A study of Parkinson's disease and MAO gene\\
masize            &      Sample size calculation for mediation analysis\\
metap             &      Meta-analysis of p values\\
metareg           &      Fixed and random effects model for meta-analysis\\
mhtplot           &      Manhattan plot of p values\\
mia               &      multiple imputation analysis for hap\\
mtdt              &      Transmission/disequilibrium test of a multiallelic marker\\
mtdt2             &      Transmission/disequilibrium test of a multiallelic marker\\
                  &      by Bradley-Terry model\\
muvar             &      Means and variances under 1- and 2- locus (diallelic) QTL model\\
mvmeta            &      Multivariate meta-analysis based on generalized least squares\\
nep499            &      A study of Alzheimer's disease with eight SNPs and APOE\\
pbsize            &      Power for population-based association design\\
pbsize2           &      Power for case-control association design\\
pedtodot          &      Converting pedigree(s) to dot file(s)\\
pfc               &      Probability of familial clustering of disease\\
pfc.sim           &      Probability of familial clustering of disease\\
pgc               &      Preparing weight for GENECOUNTING\\
plot.hap.score    &      Plot haplotype frequencies versus haplotype score statistics\\
print.hap.score   &      Print a hap.score object\\
qqfun             &      Quantile-comparison plots\\
qqunif            &      Q-Q plot for uniformly distributed random variable\\
read.ms.output    &      A utility function to read ms output\\
s2k               &      Statistics for 2 by K table\\
snca              &      A study of Parkinson's disease and SNCA markers\\
tscc              &      Power calculation for two-stage case-control design\\
twinan90          &      Classic twin models\\
whscore           &      Whittemore-Halpern scores for allele-sharing\\
\\
\end{tabular}

Assuming proper installation, you will be able to obtain the list by typing
\texttt{library(help=gap)} or view the list within a web browser via 
\texttt{help.start()}. A PDF version of this file can be viewed with command
\texttt{vignette("gap",package="gap")}.\\

You can cut and paste examples at end of each function's documentation.\\

Both \textit{genecounting} and \textit{hap} are able to handle SNPs and multiallelic
markers, with the former be flexible enough to include features such as X-linked data
and the later being able to handle large number of SNPs. But they are unable to
recode allele labels automatically, so functions \textit{gc.em} and \textit{hap.em}
are in \textit{haplo.em} format and used by a modified function \textit{hap.score} in
association testing.\\

It is notable that multilocus data are handled differently from that in {\bf hwde} and
elegant definitions of basic genetic data can be found in {\bf genetics} package.\\

Incidentally, I found my C mixed-radixed sorting routine as in \cite{zhao03} is much
faster than R's internal function.\\

With exceptions such as function \textit{pfc} which is very computer-intensive, most
functions in the package can easily be adapted for analysis of large datasets involving
either SNPs or multiallelic markers. Some are utility functions, e.g. \textit{muvar}
and \textit{whscore}, which will be part of the other analysis routines in the future.\\

The benefit with R compared to standalone programs is that for users, all functions have
unified format. For developers, it is able to incorporate their C/C++ programs more
easily and avoid repetitive work such as preparing own routines for matrix algebra and
linear models. Further advantage can be taken from packages in {\bf Bioconductor}, which
are designed and written to deal with large number of genes.\\

I have included {ms} code and .xls files to accompany {\em read.ms.output} and {\em FPRP}
and {\em BFDP} functions as with a classic twin example for ACE model in {\bf OpenMx}. The
package can be installed with command,

\begin{Schunk}
\begin{Scode}
source('http://openmx.psyc.virginia.edu/getOpenMx.R')
\end{Scode}
\end{Schunk}

\section{Demos}

You can also try several simple examples via \textit{demo}:

\begin{Schunk}
\begin{Scode}
library(gap)
demo(gap)
\end{Scode}
\end{Schunk}

\section{Examples}

I would like to highlight {\em pbsize}, {\em fbsize} and {\em ccsize} functions used for
power/sample calculations in a genome-wide asssociatoin study as reported in~\cite{zhao07}.\\

\subsection{Study design}
\subsubsection*{Family-based design}
The example involving family-based design is as follows,
<<echo=TRUE>>=
library(gap)
models <- matrix(c(
         4.0, 0.01,
         4.0, 0.10,
         4.0, 0.50, 
         4.0, 0.80,
         2.0, 0.01,
         2.0, 0.10,
         2.0, 0.50,
         2.0, 0.80,
         1.5, 0.01,    
         1.5, 0.10,
         1.5, 0.50,
         1.5, 0.80), ncol=2, byrow=TRUE)
outfile <- "fbsize.txt"
cat("gamma","p","Y","N_asp","P_A","H1","N_tdt","H2","N_asp/tdt","L_o","L_s\n",file=outfile,sep="\t")
for(i in 1:12) {
    g <- models[i,1]
    p <- models[i,2]
    z <- fbsize(g,p)
    cat(z$gamma,z$p,z$y,z$n1,z$pA,z$h1,z$n2,z$h2,z$n3,z$lambdao,z$lambdas,file=outfile,append=TRUE,sep="\t")
    cat("\n",file=outfile,append=TRUE)
}
table1 <- read.table(outfile,header=TRUE,sep="\t")
nc <- c(4,7,9)
table1[,nc] <- ceiling(table1[,nc])
dc <- c(3,5,6,8,10,11)
table1[,dc] <- round(table1[,dc],2)
unlink(outfile)
# APOE-4, Scott WK, Pericak-Vance, MA & Haines JL
# Genetic analysis of complex diseases 1327
g <- 4.5
p <- 0.15
cat("\nAlzheimer's:\n\n")
fbsize(g,p)
table1
@

\subsubsection*{Population-based design}
The example involving population-based design is as follows,
<<echo=TRUE>>=
library(gap)
kp <- c(0.01,0.05,0.10,0.2)
models <- matrix(c(
          4.0, 0.01,
          4.0, 0.10,
          4.0, 0.50, 
          4.0, 0.80,
          2.0, 0.01,
          2.0, 0.10,
          2.0, 0.50,
          2.0, 0.80,
          1.5, 0.01,    
          1.5, 0.10,
          1.5, 0.50,
          1.5, 0.80), ncol=2, byrow=TRUE)
outfile <- "pbsize.txt"
cat("gamma","p","p1","p5","p10","p20\n",sep="\t",file=outfile)
for(i in 1:dim(models)[1])
{
   g <- models[i,1]
   p <- models[i,2]
   n <- vector()
   for(k in kp) n <- c(n,ceiling(pbsize(k,g,p)))
   cat(models[i,1:2],n,sep="\t",file=outfile,append=TRUE)
   cat("\n",file=outfile,append=TRUE)
} 
table5 <- read.table(outfile,header=TRUE,sep="\t")
table5
@

\subsubsection*{Case-cohort design}
For case-cohort design, we obtain results for ARIC and EPIC studies.
<<echo=TRUE>>=
library(gap)
# ARIC study
outfile <- "aric.txt"
n <- 15792
pD <- 0.03
p1 <- 0.25
alpha <- 0.05
theta <- c(1.35,1.40,1.45)
beta1 <- 0.8
s_nb <- c(1463,722,468)
cat("n","pD","p1","hr","q","power","ssize\n",file=outfile,sep="\t")
for(i in 1:3)
{
  q <- s_nb[i]/n
  power <- ccsize(n,q,pD,p1,alpha,log(theta[i]))
  ssize <- ccsize(n,q,pD,p1,alpha,log(theta[i]),beta1)
  cat(n,"\t",pD,"\t",p1,"\t",theta[i],"\t",q,"\t",signif(power,3),"\t",ssize,"\n",
      file=outfile,append=TRUE)
}
read.table(outfile,header=TRUE,sep="\t")
unlink(outfile)
# EPIC study
outfile <- "epic.txt"
n <- 25000
alpha <- 0.00000005
power <- 0.8
s_pD <- c(0.3,0.2,0.1,0.05)
s_p1 <- seq(0.1,0.5,by=0.1)
s_hr <- seq(1.1,1.4,by=0.1)
cat("n","pD","p1","hr","alpha","ssize\n",file=outfile,sep="\t")
# direct calculation
for(pD in s_pD)
{
   for(p1 in s_p1)
   {
      for(hr in s_hr)
      {
         ssize <- ccsize(n,q,pD,p1,alpha,log(hr),power)
         if (ssize>0) cat(n,"\t",pD,"\t",p1,"\t",hr,"\t",alpha,"\t",ssize,"\n",
                          file=outfile,append=TRUE)
      }
   }
}
read.table(outfile,header=TRUE,sep="\t")
unlink(outfile)
@

\subsection{Kinship calculation}
Next, I will provide an example for kinship calculation using {\em kin.morgan}.
It is recommended that individuals in a pedigree are ordered so that parents
always precede their children. In this regard, package {\bf pedigree} can be
used, and package {\bf kinship} can be used to produce pedigree diagram as with
kinship matrix.

\subsubsection*{Pedigree diagram}
<<echo=TRUE>>=
# pedigree diagram
data(lukas,package="gap")
library(kinship)
kmat <- with(lukas,kinship(id,father,mother))
ped <- with(lukas,pedigree(id,father,mother,sex))
png("figures/lukas.png")
par(xpd=TRUE)
with(ped, plot(ped, id=paste("\n",id)))
dev.off()
@

The pedigree diagram is as follows,

\includegraphics[height=14cm,width=15cm]{figures/lukas.png}

\subsubsection*{Kinship calculation}
We then turn to the kinship calculation.

<<echo=TRUE>>=
# unordered individuals
library(gap)
gk1 <- kin.morgan(lukas)
write.table(gk1$kin.matrix,"results/gap_1.txt",quote=FALSE)

library(kinship)
kk1 <- kinship(lukas[,1],lukas[,2],lukas[,3])
write.table(kk1,"results/kinship_1.txt",quote=FALSE)

d <- gk1$kin.matrix-kk1
sum(abs(d))

# order individuals so that parents precede their children
library(pedigree)
op <- orderPed(lukas)
olukas <- lukas[order(op),]
gk2 <- kin.morgan(olukas)

write.table(olukas,"olukas.csv",quote=FALSE)
write.table(gk2$kin.matrix,"results/gap_2.txt",quote=FALSE)

kk2 <- kinship(olukas[,1],olukas[,2],olukas[,3])
write.table(kk2,"results/kinship_2.txt",quote=FALSE)

z <- gk2$kin.matrix-kk2
sum(abs(z))
@

We see that in the second case, the result agrees with {\bf kinship}.


\subsection{Graphics examples}
I now include some figures from the documentation that may be of interest.\\

\subsubsection*{Genome-wide association}
The following code is used to obtain a Q-Q plot via {\em qqunif} function,
<<echo=TRUE>>=
library(gap)
png("figures/qqunif.png")
u_obs <- runif(1000)
r <- qqunif(u_obs,pch=21,bg="blue",bty="n")
u_exp <- r$y
hits <- u_exp >= 2.30103
points(r$x[hits],u_exp[hits],pch=21,bg="green")
dev.off()
@

\includegraphics[height=10cm,width=10cm]{figures/qqunif.png}

The code below obtains a Manhattan plot via the {\em mhtplot} function,
<<echo=TRUE>>=
library(gap)
png("figures/mhtplot.png")
data <- with(mhtdata,cbind(chr,pos,p))
glist <- c("IRS1","SPRY2","FTO","GRIK3","SNED1","HTR1A","MARCH3","WISP3","PPP1R3B",
           "RP1L1","FDFT1","SLC39A14","GFRA1","MC4R")
hdata <- subset(mhtdata,gene%in%glist)[c("chr","pos","p","gene")]
color <- rep(c("lightgray","gray"),11)
glen <- length(glist)
hcolor <- rep("red",glen)  
par(las=2, xpd=TRUE, cex.axis=1.8, cex=0.4)
ops <- mht.control(colors=color,yline=1.5,xline=3)
hops <- hmht.control(data=hdata,colors=hcolor)
mhtplot(data,ops,hops,pch=19)
axis(2,pos=2,at=1:16)
title("Manhattan plot with genes highlighted",cex.main=1.8)
dev.off()
@

\includegraphics[height=14cm,width=15cm]{figures/mhtplot.png}

The code below obtains a regional association plot with the {\em asplot} function,
<<echo=TRUE>>=
library(gap)
png("figures/asplot.png")
asplot("rs10811661", "CDKN2A/CDKN2B region", "9", CDKNlocus, CDKNmap, CDKNgenes, 5.4e-8, c(3,6))
dev.off()
@

\includegraphics[height=14cm,width=14cm]{figures/asplot.png}

\subsubsection*{Effect size plot}
The purpose of this function is simply for illustration.

The code below obtains an effect size plot via the {ESplot} function,
<<echo=TRUE>>=
library(gap)
png("figures/ESplot.png")
options(stringsAsFactors=FALSE)
testdata <- data.frame(models=c("Basic model","Adjusted","Moderately adjusted","Heavily adjusted","Other"),
OR = c(4.5,3.5,2.5,1.5,1),
SElogOR = c(0.2,0.1,0.5,0.5,0.2))
ESplot(testdata,v=1)
title("This is a fictitious plot")
dev.off()
@

\includegraphics[height=10cm,width=10cm]{figures/ESplot.png}

Note that all these can serve as templates to customize features of your own.\\

\section{Known bugs}

Unaware of any bug. However, better memory management is expected.

\section{Bibliographic note}

The main references are ~\cite{chow60}, ~\cite{guo92}, ~\cite{wcn92}, ~\cite{gholamic94},
~\cite{hartung08}, ~\cite{risch96}, ~\cite{spielman96}, ~\cite{risch97}, ~\cite{miller97},
~\cite{sham97}, ~\cite{elston75},
~\cite{sham98}, ~\cite{devlin99}, ~\cite{zhao99}, ~\cite{guo00}, ~\cite{hirotsu01}, 
~\cite{zhao02}, ~\cite{zaykin02}, ~\cite{zhao04}, ~\cite{wacholder04}, ~\cite{wang05},
~\cite{skol06}, ~\cite{wakefield07}.

\bibliography{gap}
\end{document}

% Date last modified: 3 May 2011
