%\VignetteIndexEntry{Summary Tables}
%\VignetteDepends{GeneticsBase}
%\VignetteKeywords{Genetics}
%\VignettePackage{GeneticsBase}

\documentclass[10pt]{article}
\usepackage{url}
\usepackage{amsmath}
\usepackage{natbib}
\usepackage{isorot}
\usepackage{fullpage} % standard 1 inch margins 

\newcommand{\Robject}[1]{{\texttt{#1}}}
\newcommand{\Rfunction}[1]{{\texttt{#1}}}
\newcommand{\Rpackage}[1]{{\textit{#1}}}
\newcommand{\Rclass}[1]{{\textit{#1}}}
\newcommand{\Rmethod}[1]{{\textit{#1}}}
\newcommand{\code}[1]{\texttt{#1}}
\newcommand{\email}[1]{\texttt{#1}}

%%% Hyperlinks for ``PDF Latex'' :
\ifx\pdfoutput\undefined%%--- usual ``latex'' :
  %% Stuff w/out hyperref
\else%%---------------------- `` pdflatex '' : -- still gives funny errors
  \RequirePackage{hyperref}
  %% The following is R's share/texmf/hyperref.cfg :
  %% Stuff __with__ hyperref :
  \hypersetup{%
    %default: hyperindex,%
    colorlinks,%
    %default: pagebackref,%
    linktocpage,%
    %%plainpages=false,%
    linkcolor=Green,%
    citecolor=Blue,%
    urlcolor=Red,%
    pdfstartview=Fit,%
    pdfview={XYZ null null null}%
    }
  \RequirePackage{color}
  \definecolor{Blue}{rgb}{0,0,0.8}
  \definecolor{Green}{rgb}{0.1,0.75,0.1}
  \definecolor{Red}{rgb}{0.7,0,0}
  %% ESS JCGS v2 :
  %%\hypersetup{backref,colorlinks=true,pagebackref=true,hyperindex=true}
  %%\hypersetup{backref,colorlinks=false,pagebackref=true,hyperindex=true}
\fi


\usepackage{Sweave}
\begin{document}

\title{Tutorial: Using \Rpackage{GeneticsBase}}
\author{
        Gregory Warnes \\
        \email{gregory\_warnes\@urmc.rochester.edu} \\
        University of Rochester \\
        \\
        Ross Lazarus \\
        \email{ross.lazarus@channing.harvard.edu}\\
        Channing Laboratory
        }       

\maketitle

%%
%% Set up some defaults
%%
<<show=F>>=
options(width=90)
@ 


\section{Introduction}

This vignette was created as a tutorial for the 2007 BioConductor
User's Conference held in Seattle, WA, USA during August 2007, and was
presented by Dr. Warnes and Dr. Lazarus.  The material is structured
as a tutorial with a small example data set (8184 Markers x 180
Subjects belonging to 50 Families) .


\section{Outline}
\begin{enumerate}
\item Preliminaries
  \begin{enumerate}
  \item  Install the necessary packages
  \item  Load the libraries
  \end{enumerate}
  
\item Loading Data
  \begin{enumerate}  
  \item Read data from files
  \item  Error check the loaded data
  \end{enumerate}
  
\item Descriptive Statistics
  \begin{enumerate}  
  \item  Summary of allele/genotyope frequency
  \item  HWE test
  \item  Visualize Disequilbrium
  \end{enumerate}
  
\item Hypothesis Testing 
  \begin{enumerate}  
  \item  Armitage test
  \item  Logistic with synthetic covariates   
  \item  GLM with synthetic covariates \& outcome
  \end{enumerate}
  
\item Subsetting Results and Formatting Output
  \begin{enumerate}  
  \item  Construct some output tables for top 100 significant markers
  \end{enumerate}      

\item Study planning tools (power, sample size)

\end{enumerate}

\section{Preliminaries}

\subsection{Install RGenetics packages and depenedencies}

For MS-Windows it is necessary to manually install dependencies from CRAN
<<eval=FALSE>>=
install.packages( c("xtable","combinat","gdata","gplots","mvtnorm"), dep=TRUE)
@ 

Now to install the necessary packages:
<<eval=FALSE>>=
repos <- c("http://www.warnes.net/bioc2007/", "http://cran.fhcrc.org")
install.packages(c("GeneticsBase", 
                   #"GeneticsQC", 
                   "GeneticsDesign", "fbat"),
                 repos=repos, type="source", dep=TRUE)
@ 

\subsection{Load the libraries}

<<>>=
library(GeneticsBase)
library(GeneticsDesign)
library(fbat)
@ 

\section{Loading Data}

\subsection{Read data from files}

Load the full data:
<<>>=
hm.a<-readGenes(gfile="hmCEU_YRI_chr22_ALLfbat.ped", gformat="fbat")
print(hm.a)
@ 

We'll also need just the founders later:
<<>>=
hm.f<-readGenes(gfile="hmCEU_YRI_chr22_Foundersfbat.ped", gformat="fbat")
@ 

For the purpose of speeding execution of examples, we'll also create a
smaller subset of 100 markers from the original file.   
<<>>=
hm.a2<-hm.a[1:100,]
@ 

\subsection{Error check the loaded data}

Count frequencies of missing genotypes  (requires 26 seconds on my MacBook Pro)

Number of missing genotypes per subject:
<<>>=
mG<-missGFreq(hm.a2, founderOnly=TRUE, quiet=FALSE)
head(mG$nMissSubjects)
@ 
Column headers:
\begin{description}
  \item{00}  missing both alleles 
  \item{0*} 1st allele missing while the 2nd allele is not missing
  \item{*0} 1st allele is not missing while the 2nd allele is missing
\end{description}

Number of missing genotypes per marker:
<<>>=
head(mG$nMissMarkers)
@ 

Plot counts of missing genotypes:
<<fig=TRUE>>=
par(mfrow=c(2,1))
hist(mG$nMissSubjects[,1], main="",
     xlab="number of missing genotypes")
plot(1:nrow(mG$nMissSubjects),mG$nMissSubjects[,1],
     xlab="subject ID",
     ylab="number of missing genotypes",
     type="h")
title("Counts of missing genotypes")
par(mfrow=c(1,1))
@ 

\section{Descriptive Statistics}
  \begin{enumerate}
  \item  Summary of allele/genotyope frequency
  \item  HWE test
  \item  Visualize Disequilbrium
  \end{enumerate}

Basic data quality checks for markers.  Column headings are:
\begin{description}
\item{ObsHET} observed proportion of heterozygous genotypes per marker
\item{PredHET} predicted proportion of heterozygous genotypes per marker
\item{HWpval} pvalues of Hardy-Weinberg test per marker
\item{pGeno} percentage of non-missing genotypes for markes
\item{MAF} minor allele frequencies. missing allele are excluded from calculation
\item{Rating} 1 if passes HW test; 0 if failed HW test.
\end{description}
<<>>=
################
cM<-checkMarkers(hm.a2)
head(cM)
@ 

Check Mendelian errors:
<<>>=
cMend<-checkMendelian(hm.a2, quiet = FALSE)
@ 

Number of Mendelian errors per marker:
<<>>=
head(cMend$nMerrMarker)
@ 

Number of Mendelian errors per family:
<<>>=
head(cMend$nMerrFamily)
@ 

Plot counts of Mendelian errors
<<fig=TRUE>>=
par(mfrow=c(2,1))
hist(cMend$nMerrMarker, main="",
     xlab="number of Mendelian errors")
plot(1:length(cMend$nMerrFamily),cMend$nMerrFamily, 
     xlab="family ID", ylab="number of Mendian Errors")
par(mfrow=c(1,1))
@ 

\subsection{Summary of allele/genotype frequency based on GeneticsBase functions}

Allele summary, including allele counts, allele frequencies, 95\% CI of allele frequencies:
<<>>=
t1<-alleleSummary(hm.a2[1:10])
t1
@ 

Same for genotypes:
<<>>=
t2<-genotypeSummary(hm.a2[1:10], founderOnly=TRUE)
t2
@ 

HWE  test:
<<>>=
hwe<-HWE(hm.a2[1:10])
hwe
@ 

Defatult graphical display of LD:
<<fig=TRUE>>=
ld.small<-LD(hm.a2[1:20])
plot(ld.small)
@ 

Alternative graphical display of LD:
<<fig=TRUE>>=
ld <- LD(hm.a2)
LDView(ld@"X^2")
@ 

\section{Hypothesis Testing}

\subsection{Armitage test}

For the following examples, suppose 'A' is the minor allele, and 'a'
is the major allele.

Armitage test using additive model to code genotype:
\begin{tabular}{cc}
  genotype &  coding \\
  AA       &  2      \\
  Aa       &  1      \\
  aa       &  0      \\
\end{tabular}
<<>>=
res.A<-Armitage(hm.a2, method="A")
head(res.A)
@ 

Armitage test using recessive model to code genotype:
\begin{tabular}{cc}
  genotype &  coding \\
  AA       &  1 \\
  Aa       &  0 \\
  aa       &  0 \\
\end{tabular}
<<>>=  
res.R<-Armitage(hm.a2, method="R")
head(res.R)
@ 

Armitage test using dominant model to code genotype:
\begin{tabular}{cc}
  genotype &  coding \\
  AA       &  1 \\
  Aa       &  1 \\
  aa       &  0 \\
\end{tabular}
<<>>=
res.D<-Armitage(hm.a2, method="D")
head(res.D)
@

\subsection{Logistic regression}

First, we need to construct some synthetic covariates on the founders

<<>>=
sampleInfo(hm.f)$race <- sampleInfo(hm.f)$affected

raceval <- sampleInfo(hm.f)$race-1

sampleInfo(hm.f)$Norm0.0 <- rnorm(nobs(hm.f), mean=0.0*raceval)
sampleInfo(hm.f)$Norm0.5 <- rnorm(nobs(hm.f), mean=0.5*raceval)
sampleInfo(hm.f)$Norm1.0 <- rnorm(nobs(hm.f), mean=1.0*raceval)
sampleInfo(hm.f)$Norm1.5 <- rnorm(nobs(hm.f), mean=1.5*sampleInfo(hm.f)$race)

doSample <- function(raceval, mult)
  {
    prob <-  c( 0.33-raceval*mult, 0.33+(raceval*mult)/2, 0.33+(raceval*mult)/2)
    factor(sample( x=c("Red","Green","Blue"), size=1, p=prob, rep=T))
  }
    
sampleInfo(hm.f)$Cat0.0  <- sapply( raceval, doSample, mult=0.0 )
sampleInfo(hm.f)$Cat0.1  <- sapply( raceval, doSample, mult=0.1 )
sampleInfo(hm.f)$Cat0.2  <- sapply( raceval, doSample, mult=0.2 )
sampleInfo(hm.f)$Cat0.3  <- sapply( raceval, doSample, mult=0.3 )
@ 

Now, construct a function to fit the regression model and return the
parameters and statistics that are of interest.

%%
%% Put model function here both in executable form and as pure 
%% verbatim so comments are preserved. 
%% 

\begin{verbatim}
model <- function( markerName )
  {
    # extract requested genetic marker
    genotype <- genotypes(hm.f,marker=markerName)

    # get data frame to use for fitting the model
    mframe <- model.frame(race ~ sex + Norm0.0 + Norm0.5 + Norm1.0 + Norm1.5 + genotype,
                          data=sampleInfo(hm.f) )

    # To test significance of a term, best method is to do anova of
    # the full model against a submodel omitting the particular term.
    # This avoids issues with changes in names of factor levels,
    # presence or absence of covaraiates, etc.
    result <- try(
                  {
  fit.with    <- glm( race==1 ~ sex + Norm0.0 + Norm0.5 + Norm1.0 + Norm1.5 + as.factor(genotype),
                     data=mframe, family="binomial")
  fit.without <- glm( race==1 ~ sex + Norm0.0 + Norm0.5 + Norm1.0 + Norm1.5,
                     data=mframe, family="binomial")
                    anova(fit.with, fit.without, test="Chisq")$"P(>|Chi|)"[2]
                  }
                  )

    if(class(result)=="try-error")
      return(NA)  # or return(result) to see the error messages
    else

      result # full result.  Usually we want to specify exactly which
             # parameters and stats get returned so the format is consistent
             # across all markers.
  }
\end{verbatim}

<<show=FALSE>>=
model <- function( markerName )
  {
    # extract requested genetic marker
    genotype <- genotypes(hm.f,marker=markerName)

    # get data frame to use for fitting the model
    mframe <- model.frame(race ~ sex + Norm0.0 + Norm0.5 + Norm1.0 + Norm1.5 + genotype,
                          data=sampleInfo(hm.f) )

    # To test significance of a term, best method is to do anova of
    # the full model against a submodel omitting the particular term.
    # This avoids issues with changes in names of factor levels,
    # presence or absence of covaraiates, etc.
    result <- try(
                  {
  fit.with    <- glm( race==1 ~ sex + Norm0.0 + Norm0.5 + Norm1.0 + Norm1.5 + as.factor(genotype),
                     data=mframe, family="binomial")
  fit.without <- glm( race==1 ~ sex + Norm0.0 + Norm0.5 + Norm1.0 + Norm1.5,
                     data=mframe, family="binomial")
                    anova(fit.with, fit.without, test="Chisq")$"P(>|Chi|)"[2]
                  }
                  )

    if(class(result)=="try-error")
      return(NA)  # or return(result) to see the error messages
    else

      result # full result.  Usually we want to specify exactly which
             # parameters and stats get returned so the format is consistent
             # across all markers.
  }
@ 

Fit to a subset of 50, then 100 and use this information to compute the expected run time for all markers:
<<>>=
t1 <- unix.time(
fits <- sapply( markerNames(hm.f)[1:50], model )
)[3]
t1
@ 
(This takes 5.365 seconds on my MacBook Pro, R 2.4.1)

<<>>=
t2 <- unix.time(
fits <- sapply( markerNames(hm.f)[1:100], model )
)[3]
t2
@ 
(This takes 11.878 seconds on my MacBook Pro, R 2.4.1)



Estimate total time to complete, in minutes
<<>>=
t1 + (t2 - t1)/50 * (nmarker(hm.f)-50) / 60
@ 
(This yields 23.02 minutes on my MacBook Pro, R 2.4.1)


Plot the p-values for the first 100 markers
<<fig=TRUE>>=
fits.sorted <- sort(fits)
plot(-log(fits),  type="h", col="blue")
labels <- c(0.05, 0.01, 0.001, 0.0001, 0.00001, 0.00001)
abline(h=-log(labels), lty=2, col="red")
mtext(text=as.character(labels), side=4, at=-log(labels), col="red", las=1)
title("Per-marker statistical significance")
@ 

\subsection{Family-Based Associaton Test ('fbat')}

Do the fbat calculations:
<<>>=
f<-fbat(hm.a2)
@ 

Show the p-values:
<<>>=
summaryPvalue(f)
@ 

Look at the fit details for a specific marker:
<<>>=
viewstat(f, "22_14884399_rs4911642")
@ 

\section{Study planning tools (GeneticsDesign package)}

\subsection{Power to detect a low-frequencty allele}

Compute the probability of missing an allele with
frequency 0.15 when 20 genotypes are sampled:
<<>>=
gregorius(freq=0.15, N=20)
@ 

Determine what sample size is required to observe all alleles with true
frequency 0.15 with probability 0.95
<<>>=
gregorius(freq=0.15, missprob=1-0.95)
@ 

\subsection{Power for a genetics study using a quantitative outcome}
Calculate power for genetics study using a quantitative outcome:
<<>>=
GeneticPower.Quantitative.Numeric(
                                  N=50,
                                  freq=0.1,
                                  minh="recessive",
                                  alpha=0.05
                                  )

GeneticPower.Quantitative.Factor(
                                  N=50,
                                  freq=0.1,
                                  minh="recessive",
                                  alpha=0.05
                                 )
@ 

For a range of sample sizes:
<<>>=
power.range <- function( N, ... )
  {
    sapply(
           N,
           function(n)
           GeneticPower.Quantitative.Numeric(
                                             N=n,
                                             ...
                                             )
           )
  }

power.range(
            N=c(25,50,100,200,500), 
            freq=0.1, 
            minh="recessive",
            alpha=0.05 
            )

@ 

Create a power table:
<<>>=
fun <- function(p)
  power.range(freq=p,
              N=seq(100,1000,by=100), 
              alpha=5e-2, 
              minh='recessive')

m <- sapply( X=seq(0.1,0.9, by=0.1), fun )
colnames(m) <- seq(0.1,0.9, by=0.1)
rownames(m) <- seq(100,1000,by=100)

print(round(m,2))
@ 

\subsection{Power calculator for genetic linear trend association studies.}

The power is for the test that disease is associated with a
marker, given high risk allele frequency ('A'), disease
prevalence, genotype relative risk ('Aa'), genotype relative risk
('AA'), LD measure ($D'$ or $R^2$), marker allele frequency ('B'),
number of cases, control:case ratio, and probability of the Type I
error. The linear trend test (Cochran 1954; Armitage 1955) is
used.

Using $R^2$ as the measuer of LD:
<<>>=     
res1<-GPC(pA=0.05, pD=0.1, RRAa=1.414, RRAA=2, r2=0.9, pB=0.06, 
                        nCase=500, ratio=1, alpha=0.05, quiet=FALSE)
@ 

Using $D'$ as the measure of LD:
<<>>=
res2<-GPC.default(pA=0.05, pD=0.1, RRAa=1.414, RRAA=2, Dprime=0.9, pB=0.06, 
                  nCase=500, ratio=1, alpha=0.05, quiet=FALSE)

@ 

\begin{center}
  \textbf{
    The End.
    }
\end{center}
  

%\clearpage
%
%\bibliographystyle{biometrics}
%\begin{thebibliography}{}
%
%\item [~\hspace{0.125in}~Warnes~GR.]{ ``The Genetics
%    Package: Utilities for handling genetic data'' \code{R News},
%    Volume 3, Issue 1, June 2003.}
%
%\item [~\hspace{0.125in}~Warnes~GR.] ``genetics'', a package for
%  handling marker-based genetic data within the open-source
%  statistical package ``R''.  The package includes function to compute
%  allele frequencies, use genetic markers in statistical models,
%  estimate disequilibrium, and test for departure from Hardy-Weinberg
%  equilibrium. \\
%  \verb+http://cran.us.r-project.org/src/contrib/PACKAGES.html#genetics+, 
%  2002-
%
%\end{thebibliography}


\end{document}

