%\VignetteIndexEntry{GWAS Data Cleaning}
%\VignetteDepends{GWASTools, GWASdata, gdsfmt, SNPRelate}

\documentclass[11pt]{article}
\usepackage{fullpage}
\usepackage{Sweave}
\usepackage{amsmath}  
\usepackage{graphicx}
\usepackage[utf8]{inputenc}
\usepackage[pdftex,plainpages=false, letterpaper, bookmarks, bookmarksnumbered, colorlinks, linkcolor=blue, citecolor=blue, filecolor=blue, urlcolor=blue]{hyperref}

% Different font in captions
\newcommand{\captionfonts}{\footnotesize}

\begin{document}

\title{GWAS Data Cleaning}
\author{GENEVA Coordinating Center\\
  Department of Biostatistics\\
  University of Washington}

\maketitle
\tableofcontents

\SweaveOpts{keep.source=TRUE, eps=FALSE} 

\section{Overview}
This vignette takes a user through the data cleaning steps developed and used for genome wide association data as part of the Gene Environment Association
 studies (GENEVA) project. This project (\href{http://www.genevastudy.org}{http://www.genevastudy.org}) is a collection of whole-genome studies supported by the NIH-wide 
 Gene-Environment Initiative.  
The methods used in these vignettes have been published in 
Laurie et al.\ (2010).\footnote{Laurie, Cathy C., et
  al. Quality Control and Quality Assurance in Genotypic Data for
  Genome-Wide Association Studies. \textit{Genetic Epidemiology} \textbf{34}, 591-602 (August 2010).}

 For replication purposes the data used here are taken from the HapMap project.
 These data were kindly provided by 
 the Center for Inherited Disease Research (CIDR) at Johns Hopkins University and the Broad Institute of MIT and Harvard University (Broad). The data 
 are in the same format as these centers use in providing data to investigators: the content and format of these data are a little different from those
  for processed data available at the HapMap project site. The data supplied here should not be used for any purpose other than this tutorial.

\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Preparing Data}

\subsection{Data formats used in GWASTools}

The GWASTools package provides containers for storing annotation data
called {\tt SNPAnnotationDataFrame} and {\tt ScanAnnotationDataFrame} (derived
from the {\tt AnnotatedDataFrame} class in the Biobase package).  The name
``scan'' refers to a single genotyping instance.  Some subjects in a
study are usually genotyped multiple times for quality control
purposes, so these subjects will have duplicate scans.  Throughout
this tutorial, ``scan'' or ``sample'' refers to a unique genotyping instance.

The AnnotationDataFrame classes provide a way to store metadata about
an annotation variable in the same R object as the variable itself.
When a new column is added to an AnnotationDataFrame, we also add a
column to the metadata describing what that data means.
The SNP and scan AnnotationDataFrame objects are stored in R data objects (.RData files) which can be directly loaded into R.

The raw and called genotype data are stored in the Network Common Data Format (NetCDF). 
NetCDF is a set of software libraries and machine-independent data 
formats, designed specifically for large sets of array-oriented scientific data, 
and its use is widespread in those sciences that require very large data 
sets. It is maintained by the \href{http://www.unidata.ucar.edu/software/netcdf/}{Unidata program 
at the University Corporation for Atmospheric Research} (UCAR).
This portable binary file format is convenient since it allows for efficient multi-dimensional arrayed data
(although we only use up to two dimensions). 

In the GWASTools package, access to the NetCDF files is provided by
the {\tt NcdfGenotypeReader} and {\tt NcdfIntensityReader} classes.
These classes are built on top of the ncdf package and provide access
to a standard set of variables defined for GWAS data.  Additionally,
the NetCDF files and SNP and scan annotation can be linked through
the {\tt GenotypeData} and {\tt IntensityData} classes, which have
slots for a {\tt NcdfGenotypeReader} (or {\tt NcdfIntensityReader})
object, a {\tt SnpAnnotationDataFrame} object, and a {\tt
  ScanAnnotationDataFrame} object.  When an object of one of these classes is
created, it performs checks to ensure that the annotation matches the
data stored in the NetCDF file and all required information is
present.  The majority of the functions in the GWASTools package take
{\tt GenotypeData} or {\tt IntensityData} objects as arguments.

  

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

\subsection{Creating the SNP Annotation Data Object}
\label{create_snpann}

All of the functions in GWASTools require a minimum set of variables in the SNP annotation data object.  The minimum required variables are

\begin{itemize}
\item{{\tt snpID}, a unique integer identifier for each SNP}
\item{{\tt chromosome}, an integer mapping for each chromosome, with values 1-27, mapped in order from 1-22, 23=X, 24=XY (the pseudoautosomal region), 25=Y, 26=M (the mitochondrial probes), and 27=U (probes with unknown positions)}
\item{{\tt position}, the base position of each SNP on the chromosome.}
\end{itemize}

We create the integer chromosome mapping for a few reasons.  The
chromosome is stored as an integer in the NetCDF files, so in order to
link the SNP annotation with the NetCDF file, we use the integer
values in the annotation as well.
For convenience when using GWASTools functions,
the chromosome variable is most times assumed to be an integer value.  Thus, for the sex chromosomes, we can simply
use the {\tt chromosome} values.  For presentation of results, it is important to have the mapping of the
integer values back to the standard designations for the chromosome
names, thus the {\tt getChromosome()} functions in the GWASTools
objects have a {\tt char=TRUE} option to return the characters 1-22,
X, XY, Y, M, U.
The position variable should hold all numeric values of the physical position of a probe.  {\em The SNP annotation file is assumed to list the probes in order of chromosome and 
position within chromosome.}

<<>>=
library(GWASTools)
library(GWASdata)
# Load the SNP annotation (simple data frame)
data(affy_snp_annot)
# Create a SnpAnnotationDataFrame
snpAnnot <- SnpAnnotationDataFrame(affy_snp_annot)
# names of columns
varLabels(snpAnnot)
# data
head(pData(snpAnnot))
# Add metadata to describe the columns
meta <- varMetadata(snpAnnot)
meta[c("snpID", "chromosome", "position", "rsID", "probeID"), 
  "labelDescription"] <- c("unique integer ID for SNPs",                      
  paste("integer code for chromosome: 1:22=autosomes,",
   "23=X, 24=pseudoautosomal, 25=Y, 26=Mitochondrial, 27=Unknown"), 
  "base pair position on chromosome (build 36)",
  "RS identifier",
  "unique ID from Affymetrix")
varMetadata(snpAnnot) <- meta
@ 



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

\subsection{Creating the Scan Annotation Data Object}
\label{create_sampann}

The scan annotation file holds attributes for each genotyping scan that are relevant to genotypic data
cleaning.  These data include processing variables such as tissue type, DNA extraction method, and genotype processing
batch.  They also include individual characteristics such as gender and race.  The initial sample
annotation file is created from the raw data supplied by the genotyping center and/or study investigator, providing a 
mapping from the raw data file(s) for each sample scan to other sample information such as sex,
coded as {\tt M} and {\tt F}, ethnicity,  
unique scan identifier, called {\tt scanID}, and unique subject identifier.  
Since a single subject may have been genotyped multiple times as a
quality control measure, it is important to distinguish between the
scanID (unique genotyping instance) and subjectID (person
providing a DNA sample).

<<>>=
# Load the scan annotation (simple data frame)
data(affy_scan_annot)
# Create a ScanAnnotationDataFrame
scanAnnot <- ScanAnnotationDataFrame(affy_scan_annot)
# names of columns
varLabels(scanAnnot)
# data
head(pData(scanAnnot))
# Add metadata to describe the columns
meta <- varMetadata(scanAnnot)
meta[c("scanID", "subjectID", "family", "father", "mother",
  "CoriellID", "race", "sex", "status", "genoRunID", "plate",
  "alleleFile", "chpFile"), "labelDescription"] <- 
   c("unique integer ID for scans", 
  "subject identifier (may have multiple scans)",
  "family identifier", 
  "father identifier as subjectID",
  "mother identifier as subjectID",
  "Coriell subject identifier",
  "HapMap population group",
  "sex coded as M=male and F=female",
  "simulated case/control status" ,
  "genotyping instance identifier",
  "plate containing samples processed together for genotyping chemistry",
  "data file with intensities",
  "data file with genotypes and quality scores")
varMetadata(scanAnnot) <- meta
@ 


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


\subsection{Creating the NetCDF Files}


The data for genotype calls, allelic intensities and other variables such as BAlleleFrequency are stored as NetCDF files. 
More information on the NetCDF file format in general can be found at \\ \href{http://www.unidata.ucar.edu/software/netcdf/}{http://www.unidata.ucar.edu/software/netcdf/}.
The GWASTools package depends on the ncdf library.
Documentation for the ncdf R library can be found at \\
\href{http://cran.r-project.org/web/packages/ncdf/ncdf.pdf}{http://cran.r-project.org/web/packages/ncdf/ncdf.pdf}.
The ncdf library provides a convenient R interface to create, populate and extract data from NetCDF files. 

For each study, three different NetCDF files are  
created to be used in subsequent cleaning and analysis steps.   
This format is used for the ease with which multi-dimensional arrays of data can be stored and accessed.

All NetCDF files created have two dimensions, one called {\tt snp} and one titled {\tt sample}.  The {\tt snp}
dimension is of the same length as the number of probes that were released from the genotyping center and listed in the
SNP annotation file.  The 
{\tt sample} dimension is of length equal to the number of genotyping scans released as listed in the sample annotation file.
Further, all NetCDF files have three variables in common: {\tt sampleID, chromosome} and {\tt position}.
The {\tt sampleID} is used for indexing the columns of the two dimensional values stored in the NetCDF files
(genotype calls, for example).  The {\tt sampleID} ordering must match the {\tt scanID} values as listed
in the sample annotation file, see Section~\ref{create_sampann}.
The index to the SNP probes in the NetCDF file is the {\tt snpID}, which is stored as values of the SNP dimension.  Since {\tt snpID} is in chromosome and position order, these variables also provide a check on ordering and are often used to select subsets of SNPs for analysis.  Analogous to the sample ordering, these values must match the {\tt snpID} values listed
in the SNP annotation file, see \ref{create_snpann}.
To prevent errors in ordering samples or SNPs, the functions in the GWASTools package take as arguments R objects
which will return an error on creation if the sample and SNP
annotation does not match the NetCDF file.
We recommend always checking the ordering of these variables before writing new versions of the SNP or sample annotation
data files.

\subsubsection*{Genotype NetCDF Files}

The genotype NetCDF files store genotypic data in 0, 1, 2 format
indicating the number of ``A'' alleles in the genotype 
(i.e. AA=2, AB=1, BB=0 and missing=-1).  The conversion from AB format and forward strand (or other)
allele formats can be stored in the SNP annotation file.

The genotypic data are stored as a two-dimensional array, where rows are SNPs and columns are samples.
To store the genotype data, the raw data files are opened and checked to ensure the sample identifier from
the sample annotation file and the genotype data file match.  If no discrepencies exist, the probes
listed in the file are checked against the expected list of probes, then ordered and written to the NetCDF file.
This process iterates over each file (sample).  Diagnostics are stored as the process continues so that after the
data are written one can ensure the function performed as expected.   

\subsubsection*{Creating the Genotype NetCDF file for Broad}
                                     
\begin{itemize}
  
\item{The first step in creating a NetCDF file is to create a `shell' NetCDF file to which the 
data will be written.  The function {\tt ncdfCreate} creates the two dimensions and 
three common variables as described above. It also assigns values of the SNP dimension as {\tt snpID} values and populates the {\tt chromosome} and {\tt position} variables.

In this case, we also want the {\tt genotype} variable to be created, so the {\tt vars} argument must be set to 
{\tt "genotype"}. In the interest of computational feasibility, only 3 samples will be written in this step.}

<<>>=
geno.nc.file <- "tmp.geno.nc"
# get snp annotation data frame for function
snp <- affy_snp_annot[,c("snpID", "chromosome", "position")]
ncdfCreate(ncdf.filename = geno.nc.file,
       snp.annotation = snp,
       variables = c("genotype"),
       array.name = "AffyGenomeWideSNP_6",
       genome.build = "36",
       n.samples = 3,
       precision = "single")
@ 

\item{ Now the shell has been created so we can call the {\tt ncdfAddData} function to 
populate the NetCDF with the genotype data stored in the raw data files.  The data are written to the NetCDF file one sample at a time and, simultaneously, the corresponding sample identifier {\tt scanID} is written to the {\tt sampleID} variable. The {\tt chpFile}
variable from the sample annotation file
holds the name of the CHP file for each sample scan; 
these are the files we must read in to get genotype data for each sample.}

<<>>=
# first 3 samples only
scan.annotation <- affy_scan_annot[1:3, c("scanID", "genoRunID", "chpFile")]
names(scan.annotation) <- c("scanID", "scanName", "file")
# indicate which column of SNP annotation is referenced in data files
snp.annotation <- affy_snp_annot[,c("snpID", "probeID")]
names(snp.annotation) <- c("snpID", "snpName")
@ 

\item{The arguments to {\tt ncdfAddData} must be created. We need {\tt
      col.nums}, which is an integer vector indicating which columns
    of the raw text file contain variables for input.}

<<>>=
col.nums <- as.integer(c(2,3,5,6))
names(col.nums) <- c("snp","geno","a1","a2")
# Define a path to the raw data CHP text files which are read by 
#   ncdfAddData to access the raw genotypic data
path <- system.file("extdata", "affy_raw_data",
  package="GWASdata")
@ 

\item{All required arguments have been defined so we are ready to call {\tt ncdfAddData}.
This function will take the previously created genotype NetCDF file and populate it with the genotype
data from the CHP text files and with the sample identifier {\tt scanID} corresponding to that CHP file in the sample annotation data.frame.
A set of diagnostic values are written and stored in {\tt diag.geno}, so we must look at those to
ensure no errors occurred.
For more information on the function {\tt ncdfAddData}, such as argument specification, diagnostics recorded and
examples, please see the function documentation.}

<<>>=
diag.geno.file <- "diag.geno.RData"
diag.geno <- ncdfAddData(path = path,
  ncdf.filename = geno.nc.file, 
  snp.annotation = snp.annotation,
  scan.annotation = scan.annotation, 
  sep.type = "\t",
  skip.num = 1,
  col.total = 6,
  col.nums = col.nums,
  scan.name.in.file = -1,
  scan.start.index = 1,
  diagnostics.filename = diag.geno.file,
  verbose = FALSE)
# Look at the values included in the "diag.geno" object which holds 
#   all output from the function call
names(diag.geno)
# `read.file' is a vector indicating whether (1) or not (0) each file 
#   specified in the `files' argument was read successfully
table(diag.geno$read.file)
# `row.num' is a vector of the number of rows read from each file
table(diag.geno$row.num) 
# `sample.match' is a vector indicating whether (1) or not (0) 
#   the sample name inside the raw text file matches that in the 
#   sample annotation data.frame
table(diag.geno$sample.match)
# `snp.chk' is a vector indicating whether (1) or not (0) 
#   the raw text file has the expected set of SNP names
table(diag.geno$snp.chk)
# `chk' is a vector indicating whether (1) or not (0) all previous
#   checks were successful and the data were written to the NetCDF file
table(diag.geno$chk)
@ 

\item{Although the diagnostic values indicated no issues, we will open
    the NetCDF file and extract a small sampling of data.  This
    illustrates the use of the {\tt NcdfGenotypeReader} class for
    retrieving data from a NetCDF file.
  }
  
<<>>=
(genofile <- NcdfGenotypeReader(geno.nc.file))
# Take out genotype data for the first 3 samples and 
#   the first 5 SNPs
(genos <- getGenotype(genofile, snp=c(1,5), scan=c(1,3)))
# Close the NetCDF file
close(genofile)
@ 

\item{Run the function {\tt ncdfCheckGenotype} to check that the
    NetCDF file contains the same data as the raw data files.}

<<>>=
check.geno.file <- "check.geno.RData"
check.geno <- ncdfCheckGenotype(path = path,
  ncdf.filename = geno.nc.file, 
  snp.annotation = snp.annotation,
  scan.annotation = scan.annotation, 
  sep.type = "\t",
  skip.num = 1,
  col.total = 6,
  col.nums = col.nums,
  scan.name.in.file = -1,
  check.scan.index = 1:3,
  n.scans.loaded = 3,
  diagnostics.filename = check.geno.file,
  verbose = FALSE)
# Look at the values included in the "check.geno" object which holds 
#   all output from the function call
names(check.geno)
# 'snp.order' is a vector indicating whether (1) or not (0) the snp ids 
#   are in the same order in each file.
table(check.geno$snp.order)
# 'geno.chk' is a vector indicating whether (1) or not (0) the genotypes 
#   in the netCDF match the text file
table(check.geno$geno.chk)
@ 

\end{itemize}

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

\subsubsection*{Intensity NetCDF Files}


The intensity NetCDF files store quality scores and allelic intensity data for each SNP. 
The normalized X and Y intensities as well as the confidence scores are written to the NetCDF for all
samples, for all SNPs. (A separate NetCDF file will store the BAlleleFreq and LogRRatio data.)

Aside from the three variables held in common with all NetCDF files ({\tt sampleID, chromosome} and {\tt position}),
the intensity and quality data are written to the intensity NetCDF in a two dimensional format, with SNPs corresponding
to rows and samples corresponding to columns.
To write the intensity data, the raw data files are opened and the intensities and quality score
are read.  Like with the genotype data, if all sample and probe identifiers match between the data files
and the annotation files, the data are populated in the NetCDF and diagnostics are written.

Affymetrix data are provided in two files per genotyping scan.  The CHP file holds the genotype calls, used to create the 
genotype NetCDF file, as well as the confidence score, which is written to the {\tt quality} variable in
the intensity NetCDF file.
The normalized X and Y intensity data are stored in the {\tt allele\_summary} files in the format of two rows
per SNP, one for each allelic probe. 
A separate function {\tt ncdfAddIntensity} reads these data and writes them to the NetCDF file. 
Thus, when writing the intensity NetCDF file using Affymetrix data, there are two function calls
needed, each of which opens and reads from the two sets of files.

Illumina data are provided in one file per genotyping scan.  The confidence score, genotype call, normalized intensities
as well as the BAlleleFrequency and LogRRatio values are all stored in one file.
Because of this, we do not need to invoke the {\tt ncdfAddIntensity} function when reading in Illumina
data to the intensity NetCDF.

\subsubsection*{Creating the Intensity NetCDF file for Broad}
                       
\begin{itemize}
\item{The first step in creating the intensity NetCDF is to create a `shell' NetCDF
file to which the data will be written.  We can use the same function used for creating the shell
genotype NetCDF, {\tt ncdfCreate}.  The two dimensions and three common variables will be created
along with the X, Y and quality variables; the {\tt vars} argument is set to write these.
For a reasonable computation time, we will create the file to hold data for 3 samples only.}

<<>>=
qxy.nc.file <- "tmp.qxy.nc"
# get snp annotation data frame for function
snp <- affy_snp_annot[,c("snpID", "chromosome", "position")]
ncdfCreate(ncdf.filename = qxy.nc.file,
       snp.annotation = snp,
       variables = c("quality","X","Y"),
       array.name = "AffyGenomeWideSNP_6",
       genome.build = "36",
       n.samples = 3,
       precision = "single")
@ 

		
\item{We next call the {\tt ncdfAddData} function to populate the intensity NetCDF with the 
X and Y intensities and the quality score stored in the raw CHP text files. The sample identifier {\tt scanID} is also populated at the same time.  The {\tt chpFile}
variable from the sample annotation file holds the name of the CHP file for each sample scan.  
From this file we can get all needed data for the intensity NetCDF file.}

<<>>=
# first 3 samples only
scan.annotation <- affy_scan_annot[1:3, c("scanID", "genoRunID", "chpFile")]
names(scan.annotation) <- c("scanID", "scanName", "file")
# indicate which column of SNP annotation is referenced in data files
snp.annotation <- affy_snp_annot[,c("snpID", "probeID")]
names(snp.annotation) <- c("snpID", "snpName")
@ 

\item{The arguments to {\tt ncdfAddData} must be created. We need {\tt
      col.nums}, which is an integer vector indicating which columns
    of the raw text file contain variables for input.}

<<>>=
col.nums <- as.integer(c(2,4))
names(col.nums) <- c("snp","qs")
# Define a path to the raw data CHP text files which are read by 
#   ncdfAddData to access the raw genotypic data
path <- system.file("extdata", "affy_raw_data",
  package="GWASdata")
@ 

\item{All required arguments have been defined so we are ready to call {\tt ncdfAddData}.
This function will take the previously created intensity NetCDF file and populate it with 
the quality score from the CHP text files.  Recall the intensity values are stored in a different set of files 
so those will be populated in the next step.
A set of diagnostic values are written and stored in {\tt diag.qual}.}

<<>>=
diag.qual.file <- "diag.qual.RData"
diag.qual <- ncdfAddData(path = path,
  ncdf.filename = qxy.nc.file, 
  snp.annotation = snp.annotation,
  scan.annotation = scan.annotation, 
  sep.type = "\t",
  skip.num = 1,
  col.total = 6,
  col.nums = col.nums,
  scan.name.in.file = -1,
  scan.start.index = 1,
  diagnostics.filename = diag.qual.file,
  verbose = FALSE)
@ 

\item{As alluded to above, the normalized X and Y intensity values are stored in the {\tt allele files}.  
To write the intensities, we will call {\tt ncdfAddIntensity} with similar arguments to the {\tt ncdfAddData}
function.  For further explanation of this function, please refer to the function help documentation.}

<<>>=
scan.annotation <- affy_scan_annot[1:3, c("scanID", "genoRunID", "alleleFile")]
names(scan.annotation) <- c("scanID", "scanName", "file")
diag.xy.file <- "diag.xy.RData"
diag.xy <- ncdfAddIntensity(path = path,
  ncdf.filename = qxy.nc.file, 
  snp.annotation = snp.annotation, 
  scan.annotation = scan.annotation, 
  scan.start.index = 1, 
  n.consecutive.scans = 3,
  diagnostics.filename = diag.xy.file,
  verbose = FALSE)
@ 

\item{We will open
    the NetCDF file and extract a small sampling of data.  In this
    case we use the {\tt NcdfIntensityReader} class.}
 
<<>>=
# Open the NetCDF file we just created
(intenfile <- NcdfIntensityReader(qxy.nc.file))
# Take out the normalized X intensity values for the first 
#     5 SNPs for the first 3 samples
(xinten <- getX(intenfile, snp=c(1,5), scan=c(1,3)))
# Close the NetCDF file
close(intenfile)
@

\item{Run the function {\tt ncdfCheckIntensity} to check that the
    NetCDF file contains the same data as the raw data files.}

<<>>=
scan.annotation <- affy_scan_annot[1:5,
  c("scanID", "genoRunID", "chpFile", "alleleFile")]
names(scan.annotation) <- c("scanID", "scanName", "file", "inten.file")
check.qxy.file <- "check.qxy.RData"
check.qxy <- ncdfCheckIntensity(path = path,
  intenpath = path,
  ncdf.filename = qxy.nc.file, 
  snp.annotation = snp.annotation,
  scan.annotation = scan.annotation, 
  sep.type = "\t",
  skip.num = 1,
  col.total = 6,
  col.nums = col.nums,
  scan.name.in.file = -1,
  check.scan.index = 1:3,
  n.scans.loaded = 3,
  affy.inten = TRUE,
  diagnostics.filename = check.qxy.file,
  verbose = FALSE)
@ 

\end{itemize}


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


\subsubsection*{BAlleleFrequency and LogRRatio NetCDF Files}


The BAlleleFrequency and LogRRatio NetCDF file stores these values for every sample by SNP.
For Affymetrix data, these values must be calculated, but for Illumina data these values are 
calculated by the BeadStudio software and may be provided by the genotyping center.  
For the purposes of this tutorial, we will be calculating
the BAlleleFrequency and LogRRatio values from the Broad data as explained in the following manner.
For a thorough explanation and presentation of an application of these values, please refer to 
Peiffer, Daniel A., et al. (2006).\footnote{Peiffer, Daniel A., et al. High-resolution genomic profiling of chromosomal aberrations using
Infinium whole-genome genotyping. \textit{Genome Research} \textbf{16}, 1136-1148 (September 2006).}
          
                                                                                      
For a given sample and SNP, R and $\theta$ are calculated using the X and Y intensities, where
\begin{equation}
\label{rthe}
R=X+Y 
\end{equation}
\begin{equation*}
\theta = \frac{2\arctan(Y/X)}{\pi}
\end{equation*}

Illumina data may come with $R$ and $\theta$ values for each sample, where 
$\theta$ corresponds to the polar coordinate angle and R is the sum of the normalized X and Y intensities (not, as one might assume, the magnitude of the polar coordinate vector).
AffyMetrix data uses X and Y values from a Cartesian system.  
To transform to the R and $\theta$ polar coordinates, we use equation \ref{rthe}.

Regardless of which platform our data come from, we are able to find the R and $\theta$ values for 
every sample and every SNP.  It is from these values that we calculate the LogRRatio and BAlleleFrequency.
The LogRRatio is given below.  The expected value of R is derived from a plot of $\theta$ versus R for a given SNP.  
It is the predicted value of R derived from a line connecting the centers of the two nearest genotype clusters.  
\begin{equation}
 {\rm LogRRatio} =\log \left( \frac{R_{\text{observed values}}}{R_{\text{expected values}}}\right)
\end{equation} 
Variation in the LogRRatio across a single chromosome indicates possible duplication or deletion, and is an indication
of overall sample quality. 
 
The BAlleleFrequency is the frequency of the B allele in the population of cells from which the DNA is extracted.  
Each sample and SNP combination has a BAlleleFrequency value.  Note the BAlleleFrequency values vary for a subject
with each DNA extraction and tissue used. 
After all SNPs have been read and all samples have been clustered for a probe, 
the mean $\theta$ ``cluster'' value is calculated for each probe, for each of the three genotype clusters, 
resulting in $\theta_{AA}, \theta_{AB} \text{ and }\theta_{BB}$ for every probe.  
Then the $\theta$ value for each sample, call it $\theta_{n}$, 
is compared to $\theta_{AA}, \theta_{AB}\text{ and }\theta_{BB}$.
The BAlleleFrequency is calculated
\begin{equation*}
 \text{BAlleleFrequency}= 
\begin{cases}
0 & \text{if } \theta_{n} < \theta_{AA}\\
 \mbox{}\\
\frac{{\textstyle (1/2)(\theta_{n}-\theta_{AA})}}{{\textstyle \theta_{AB}-\theta_{AA}}} & \text{if } \theta_{AA} \leq \theta_{n} < \theta_{AB}\\
\mbox{}\\
\frac{{\textstyle 1}}{{\textstyle 2}}+\frac{{\textstyle (1/2)(\theta_{n}-\theta_{AB})}}
  {{\textstyle \theta_{BB}-\theta_{AB}}} & \text{if } \theta_{AB} \leq \theta_{n} < \theta_{BB}\\
 \mbox{}\\
1 & \text{if } \theta_{n} \geq \theta_{BB}
\end{cases}
\end{equation*}
A $\theta_{n}$ value of 0 or 1 corresponds to a homozygote genotype for sample $n$ at that
particular probe, and a $\theta_{n}$ value of 1/2 indicates a heterozygote genotype.
Thus, $BAlleleFrequency \in[0,1]$ for each probe.  Across a chromosome, three bands are expected, 
one hovering around 0, one around 1 and one around 0.5, and any deviation from this is considered aberrant.

We use the BAlleleFrequency and LogRRatio values to detect mixed samples or samples of low quality, 
as well as chromosomal duplications and deletions.  Samples that have a significantly large 
(partial or full chromosome) aberration for a particular chromosome as detected from the BAlleleFrequency 
values are recommended to be filtered out, for the genotype data are not reliable in these situations.  
Because of these applications, the BAlleleFrequency and LogRRatio values are a salient part of the data cleaning steps.

In addition to the three variables held in common with all NetCDF files ({\tt sampleID, chromosome} and {\tt position}),
the BAlleleFrequency and LogRRatio values are calculated and written to this NetCDF in a two dimensional format, 
with SNPs corresponding to rows and samples corresponding to columns.
Because we have already completed the creation of both the genotype and intensity NetCDF files,
we simply use those files to access the data.
The BAlleleFrequency and LogRRatio values are calculated in subsets for efficiency and written to the 
corresponding subset indices in the NetCDF file.

\subsubsection*{Creating the BAlleleFrequency and LogRRatio NetCDF file for Broad}

\begin{itemize}

\item{The first step in creating the BAlleleFrequency NetCDF is to create a `shell' NetCDF
file to which the data will be written.  We can use the same function used for creating the other shell
NetCDF files, {\tt ncdfCreate}.  The two dimensions and three common variables will be created
along with the BAlleleFrequency and LogRRatio variables.
In order to allow for a reasonable computation time, we will create the file to hold data for 3 samples only.}

<<>>=
bl.nc.file <- "tmp.bl.nc"
# get snp annotation data frame for function
snp <- affy_snp_annot[,c("snpID", "chromosome", "position")]
ncdfCreate(ncdf.filename = bl.nc.file,
       snp.annotation = snp,
       variables = c("BAlleleFreq","LogRRatio"),
       array.name = "AffyGenomeWideSNP_6",
       genome.build = "36",
       n.samples = 3,
       precision = "single")
@ 

\item{We now will calculate the BAlleleFrequency and LogRRatio values for each sample by SNP and
write these values to the NetCDF by calling the function {\tt BAFfromGenotypes}.
We will also select ``by.study'' as the call method, so all 3 samples
have their genotype clusters called together.  In normal usage, we
recommend calling Affymetrix genotypes ``by.plate'' (in which case the {\tt
  plate.name} argument is passed to the function).
For more detail regarding the {\tt BAFfromGenotypes} function, please see the function documentation.
After the function is complete, we will look at a few values to ensure the file was created successfully.}

<<>>=
xyNC <- NcdfIntensityReader(qxy.nc.file)
genoNC <- NcdfGenotypeReader(geno.nc.file)
BAFfromGenotypes(xyNC, genoNC,
           bl.ncdf.filename = bl.nc.file, 
           min.n.genotypes = 0,
           call.method ="by.study")
close(xyNC)
close(genoNC)
# Open the NetCDF file we just created
(blfile <- NcdfIntensityReader(bl.nc.file))
# Look at the BAlleleFrequency values for the first 5 SNPs
(baf <- getBAlleleFreq(blfile, snp=c(1,5), scan=c(1,3)))
# Close the NetCDF file
close(blfile)
@ 

\end{itemize}

<<echo=FALSE, results=hide>>=
file.remove(geno.nc.file, qxy.nc.file, bl.nc.file)
file.remove(diag.geno.file, diag.qual.file, diag.xy.file)
file.remove(check.geno.file, check.qxy.file)
@ 

\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Batch Quality Checks}

The overall goal of this step is to check the quality of the sample batches. Substantial quality control is done by
the genotyping centers prior to releasing the genotype data; however it is our experience that despite the stringent 
quality controls it is still possible for batches with lower than desired quality to pass the pre-release data quality 
checks. If a lower quality batch is detected then it may be necessary to re-run the genotyping for that batch. We check
the batch quality by comparing the missing call rates between batches and looking for significant allele frequency 
differences between batches.

\subsection{Calculate Missing Call Rate for Samples and SNPs}

The first step is to calculate the missing call rates for each SNP and for each sample. A high missing call rate for 
a sample is often indicative of a poorly performing sample. It has been seen that samples from DNA that has undergone 
whole-genome amplification (WGA) have a relatively higher missing call rate. Similarly a high missing call rate for a 
SNP is indicative of a problem SNP. Experience from the GENEVA studies has shown that there seem to be a subset of 
SNPs from which genotype calls are more difficult to make than others. 
We calculate the missing call rates in a two step process: first the missing call rates over all samples and SNPs 
are calculated, then the missing call rates are calculated again, filtering out SNPs and samples that have an 
initial missing call rate greater than 0.05. The initial SNP missing call rate over all samples is saved in the 
SNP annotation data file as {\tt missing.n1}.  The analogous idea is applied to the samples: {\tt missing.e1} 
is saved in the sample annotation file and corresponds to the missing call rate per sample over all SNPs, 
excluding those SNPs with all calls missing. The {\tt missing.n2} is calculated as the call rate per SNP over 
all samples whose {\tt missing.e1} is less than 0.05. Again, similarly for the samples, {\tt missing.e2} is calculated 
for each sample over all SNPs with {\tt missing.n2} values less than 0.05. It is important to remember that the Y 
chromosome values should be calculated for males only, since we expect females to have no genotype values for the 
Y chromosome.

\subsubsection*{Calculate {\tt missing.n1} }
\label{sec:missing.n1}
                                         
\begin{itemize}

\item{This step calculates and examines {\tt missing.n1}, the missing call rate per SNP over all samples by calling the 
function {\tt missingGenotypeBySnpSex}.  
This function takes a {\tt GenotypeData} object as an argument, and
requires that the scan annotation of this object contains a ``sex''
column.
There is also an option to send a vector of SNPs to exclude from the
calculation, which is what we will use later to find {\tt missing.n2}.  For now, we will use all SNPs for each sample, 
being sure to calculate by gender.  The function returns a list, with one element that holds the missing counts for each
SNP and one element that holds the sex counts; we will save these results to be used later.}

<<>>=
library(GWASTools)
library(GWASdata)
# sex is required for this function, so we load the scan annotation
data(affy_scan_annot)
scanAnnot <- ScanAnnotationDataFrame(affy_scan_annot)
# open the netCDF file and create a GenotypeData object
ncfile <- system.file("extdata", "affy_geno.nc",
  package="GWASdata")
nc <- NcdfGenotypeReader(ncfile)
genoData <-  GenotypeData(nc, scanAnnot=scanAnnot)
# Calculate the number of missing calls for each snp over all samples 
#     for each gender separately
miss <- missingGenotypeBySnpSex(genoData)
# Examine the results
length(miss)
names(miss)
head(miss$missing.counts)
miss$scans.per.sex
head(miss$missing.fraction)
@ 

\item{The Y chromosome should be missing for all females, but an occasional probe on the Y chromosome is called in 
a female. \texttt{missingGenotypeBySnpSex} excludes females when
calculating the missing rate for Y chromosome SNPs. Note this may need to be changed 
later if there are some gender mis-annotations because the Y chromosome SNP missing call rates may change.
We add the missing call rates to the SNP annotation table and save a new version. }

<<>>=
# get snp annotation
data(affy_snp_annot)
snpAnnot <- SnpAnnotationDataFrame(affy_snp_annot)
snpAnnot$missing.n1 <- miss$missing.fraction
varMetadata(snpAnnot)["missing.n1", "labelDescription"] <- paste(
  "fraction of genotype calls missing over all samples",
  "except that females are excluded for Y chr SNPs")
summary(snpAnnot$missing.n1)
@ 

\item{We will create some plots of the missing call rates so we can easily identify any outliers.
Plots will be shown of Y chromosome SNP missing call rates for males, and overall SNP missing call rates.
We also find the number of SNPs with 100\% missing, and the fraction of SNPs with missing call rate less than
0.05 for each chromosome type. }

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
hist(snpAnnot$missing.n1[snpAnnot$chromosome == 25],
     xlab="Fraction of Y chr SNP missing calls", 
     main="Y Chromosome SNP Missing Call Rate for Males")
@ 
\caption{Missing call rate over all probes on the Y chromosome for the male samples only.  The missing call 
rate for all females over the Y chromosome is identically 1. Although it may seem as if three samples have a higher missing
call rate, one must note the small scale on the x-axis.  There are no samples here that have an unacceptable missing call rate.}
\label{fig:Miss_call_rate_Y}
\end{figure}


\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
hist(snpAnnot$missing.n1, ylim=c(0,100),
     xlab="SNP missing call rate", 
     main="Missing Call Rate for All Probes")
@ 
\caption{Histogram of overall missing call rate per SNP.  Note the y-axis has been truncated in order to show the higher 
values of the missing call rate on the x-axis.}
\label{fig:Miss_call_rate_snps}
\end{figure}

<<>>=
# Find the number of SNPs with every call missing
length(snpAnnot$missing.n1[snpAnnot$missing.n1 == 1])
# Fraction of autosomal SNPs with missing call rate < 0.05
x <- snpAnnot$missing.n1[snpAnnot$chromosome < 23]
length(x[x < 0.05]) / length(x)
# Fraction of X chromosome SNPs with missing call rate < 0.05
x <- snpAnnot$missing.n1[snpAnnot$chromosome == 23]
length(x[x < 0.05]) / length(x)
# Fraction of Y chromosome SNPs with missing call rate < 0.05
x <- snpAnnot$missing.n1[snpAnnot$chromosome == 25]
length(x[x < 0.05]) / length(x)
@ 

\end{itemize}


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

\subsubsection*{Calculate {\tt missing.e1} }
                         
\begin{itemize}
  
\item{This step calculates {\tt missing.e1}, which is the missing call rate per sample over all SNPs, by chromosome.
We read in the new SNP annotation file which holds the {\tt missing.n1} variable.  
Since we have both sample and SNP annotation loaded, we create a {\tt
  GenotypeData} object including both these annotations to ensure consistency.
For those SNPs with a {\tt missing.n1}
value less than one, we call the {\tt missingGenotypeByScanChrom} function that returns a list with one element holding
the missing counts per sample by chromosome and one element holding the number of SNPs per chromosome.  We check these values
to ensure they are as expected.} 

<<>>=
# Want to exclude all SNP probes with 100% missing call rate
# Check on how many SNPs to exclude
sum(abs(snpAnnot$missing.n1 - 1) < 1e-6)
sum(snpAnnot$missing.n1 == 1)
# Create a variable that contains the IDs of these SNPs to exclude
snpexcl <- snpAnnot$snpID[abs(snpAnnot$missing.n1 - 1) < 1e-6 | 
                          snpAnnot$missing.n1 == 1]
length(snpexcl)
@ 
<<>>=
# Use the missingGenotypeByScanChrom function
miss <- missingGenotypeByScanChrom(genoData, snp.exclude=snpexcl)
length(miss)
names(miss)
head(miss$missing.counts)
head(miss$snps.per.chr)
# Check to make sure that the correct number of SNPs were excluded 
sum(miss$snps.per.chr)
nrow(snpAnnot) - sum(miss$snps.per.chr)
@ 

\item{{\tt missingGenotypeByScanChrom} calculates the missing call
    rate for each sample over all SNPs.  For females, the missing call rate does not include the probes on the Y chromosome.
The values for {\tt missing.e1} are added to the sample annotation file.}

<<>>=
head(miss$missing.fraction)
# Check the ordering matches the sample annotation file
allequal(names(miss$missing.fraction), scanAnnot$scanID)
# Add the missing call rates vector to the sample annotation file
scanAnnot$missing.e1 <- miss$missing.fraction
varMetadata(scanAnnot)["missing.e1", "labelDescription"] <- paste(
  "fraction of genotype calls missing over all snps with missing.n1<1",
  "except that Y chr SNPs are excluded for females")
summary(scanAnnot$missing.e1)
@ 

\item{We will create a histogram of the overall missing call rate per sample in order to identify any samples with
a relatively larger missing call rate.  It is known that genotype data taken from DNA that has been through whole-genome
amplification (WGA) has an overall higher missing call rate; this is something that we would see at this step if any 
samples are of WGA origin. We also look at the summary of the missing call rate for females and males separately to 
ensure there are no large gender differences. Finally, we calculate the number of samples with a missing call rate 
greater than 0.05.  In this case, there are no samples but in other data this may not be the case.
If any samples have a high missing rate, we recommend further investigation of what may be causing the missing calls;
the samples with a missing call rate greater than 0.05 should be filtered out due to low sample quality. }


\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
hist(scanAnnot$missing.e1, 
     xlab="Fraction of missing calls over all probes", 
     main="Histogram of Sample Missing Call Rate for all Samples")
@ 
\caption{Histogram of missing call rate for all samples over all probes.}
\label{fig:Miss_calls_over_probes}
\end{figure}

<<>>=
# Look at missing.e1 for males
summary(scanAnnot$missing.e1[scanAnnot$sex == "M"])
# Look at missing.e1 for females
summary(scanAnnot$missing.e1[scanAnnot$sex == "F"])
# Number of samples with missing call rate > 5%
sum(scanAnnot$missing.e1 > 0.05)
@ 


\item{For some analyses we require the missing call rate for autosomes and the X chromosome to be separated.
We calculate these values here and add them to the sample annotation file.  
Also, we will create a logical {\tt duplicated} variable.  We can identify the duplicated scans in the sample
annotation file by identifying the subject ids that occur more than once.  Among samples with the same subject id, the 
one with the lowest {\tt missing.e1} value will have the variable {\tt duplicated} set to {\tt FALSE}.
We then write a new version of the sample annotation with these added variables. }

<<>>=
auto <- colnames(miss$missing.counts) %in% 1:22
missa <- rowSums(miss$missing.counts[,auto]) / sum(miss$snps.per.chr[auto])
summary(missa)
missx <- miss$missing.counts[,"X"] / miss$snps.per.chr["X"]
summary(missx)
# check they match sample annotation file
allequal(names(missa), scanAnnot$scanID)
allequal(names(missx), scanAnnot$scanID)
# Add these separate sample missing call rates to the sample 
# annotation
scanAnnot$miss.e1.auto <- missa
scanAnnot$miss.e1.xchr <- missx
# Order scanAnnot by missing.e1 so duplicate subjectIDs 
# with a higher missing rate are marked as duplicates
scanAnnot <- scanAnnot[order(scanAnnot$subjectID, scanAnnot$missing.e1),]
scanAnnot$duplicated <- duplicated(scanAnnot$subjectID)
table(scanAnnot$duplicated, exclude=NULL)
# Put scanAnnot back in scanID order; this is very important!!
# Save the updated sample annotation
scanAnnot <- scanAnnot[order(scanAnnot$scanID),]
allequal(scanAnnot$scanID, sort(scanAnnot$scanID))
@ 

\end{itemize}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection*{Calculate {\tt missing.n2} }

\begin{itemize}
  
\item{This step calculates {\tt missing.n2}, which is the missing call rate per SNPs with {\tt missing.e1} less 
than 0.05 over all samples.  This calculation is done separately for each sex as in {\tt missing.n1} in 
Section~\ref{sec:missing.n1}. We read in the sample annotation file which holds the {\tt missing.e1} variable.  
In most cases, there will be samples with missing call rate greater than 0.05.  However, because of the high 
quality of the HapMap data, there are no samples in this case.  We will continue with the steps as if there are 
samples we must exclude from the {\tt missing.n2} calculation. We call the {\tt missingGenotypeBySnpSex} 
function just as we did to calculate for {\tt missing.n1}, but this time we include the list of sample numbers to 
exclude from the calculation (although here that list is empty).  We review the values returned from the function 
and save the results. }

<<>>=
sum(scanAnnot$missing.e1 > 0.05)
# Even though there are no samples with missing.e1 > 0.05, 
# we will go through the full process for the tutorial
# Find the samples with missing.e1 > .05 and make a vector of 
# scanID to exclude from the calculation
scan.exclude <- scanAnnot$scanID[scanAnnot$missing.e1 > 0.05]
# Call missingGenotypeBySnpSex and save the output
miss <- missingGenotypeBySnpSex(genoData, scan.exclude=scan.exclude)
length(miss)
names(miss) # Check that all SNPs are listed
dim(miss$missing.counts)
miss$scans.per.sex
dim(miss$missing.fraction)
# Make sure ordering matches snp annotation
allequal(snpAnnot$snpID, as.numeric(names(miss$missing.fraction)))
snpAnnot$missing.n2 <- miss$missing.fraction
varMetadata(snpAnnot)["missing.n2", "labelDescription"] <- paste(
  "fraction of genotype calls missing over all samples with missing.e1<0.05",
  "except that females are excluded for Y chr SNPs")
summary(snpAnnot$missing.n2)
@ 

\item{As we did in the step for calculating {\tt missing.n1}, we will create some plots of the missing call rates so 
we can easily identify any outliers. Plots will be shown of Y chromosome SNP missing call rates for males, and overall 
SNP missing call rates. We also find the number of SNPs with 100\% missing, and the fraction of SNPs with missing call 
rate less than 0.05 for each chromosome type. }

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
hist(snpAnnot$missing.n2[snpAnnot$chromosome == 25], 
     xlab="Fraction of Y chr SNP missing calls", 
     main="Y Chromosome SNP Missing Call Rate for Males")
@ 
\caption{Missing call rate over all probes on the Y chromosome for the male samples excluding all samples with a
missing call rate greater than 0.05.  The missing call rate for all females over the Y chromosome is identically 1.
Although it may seem as if three samples have a higher missing call rate, one must note the small scale on the x-axis. 
There are no samples here that have an unacceptable missing call rate.}
\label{fig:miss_call_snp_Y}
\end{figure}

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
hist(snpAnnot$missing.n2, ylim=c(0,10000), 
     xlab="SNP missing call rate", main="Histogram of SNP Missing Call Rate")
@ 
\caption{Histogram of overall missing call rate per SNP excluding samples with missing call rate greater than 0.05.
Note the y-axis has been truncated in order to show the higher values of the missing call rate on the x-axis.}
\label{fig:SNP_miss_call_rate}
\end{figure}

<<>>=
# Find the number of SNPs with every call missing
length(snpAnnot$missing.n2[snpAnnot$missing.n2 == 1])
# Fraction of autosomal SNPs with missing call rate < .05
x <- snpAnnot$missing.n2[snpAnnot$chromosome < 23]
length(x[x < 0.05]) / length(x)
# Fraction of X chromosome SNPs with missing call rate <.05
x <- snpAnnot$missing.n2[snpAnnot$chromosome == 23]
length(x[x < 0.05]) / length(x)
# Fraction of Y chromosome SNPs with missing call rate <.05
x <- snpAnnot$missing.n2[snpAnnot$chromosome == 25]
length(x[x < 0.05]) / length(x)
@ 

\end{itemize}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection*{Calculate {\tt missing.e2} }

\begin{itemize}

\item{This step calculates {\tt missing.e2}, which is the missing call rate per sample over all SNPs 
with {\tt missing.n2} less than 0.05. We read in the new SNP annotation file which holds the {\tt missing.n2} variable.  
Because we must exclude those SNPs with a sufficiently high missing call rate, we will create a vector of SNPs that 
have {\tt missing.n2} greater than 0.05.  Then, for the non-excluded SNPs we call the {\tt missingGenotypeByScanChrom}
function, which returns a list with one element holding the missing counts per sample by chromosome and one element 
holding the number of SNPs per chromosome.  We check these values to ensure they are as expected.} 

<<>>=
# Create a vector of the SNPs to exclude.
snpexcl <- snpAnnot$snpID[snpAnnot$missing.n2 >= 0.05]
length(snpexcl)
@ 
<<>>=
# Use the missingGenotypeByScanChrom function
miss <- missingGenotypeByScanChrom(genoData, snp.exclude=snpexcl)
length(miss)
names(miss)
head(miss$missing.counts)
head(miss$snps.per.chr)
# Check to make sure that the correct number of SNPs were excluded 
sum(miss$snps.per.chr)
nrow(snpAnnot) - sum(miss$snps.per.chr)
head(miss$missing.fraction)
# Check the ordering matches the sample annotation file
allequal(names(miss$missing.fraction), scanAnnot$scanID)
# Add the missing call rates vector to the sample annotation file
scanAnnot$missing.e2 <- miss$missing.fraction
varMetadata(scanAnnot)["missing.e2", "labelDescription"] <- paste(
  "fraction of genotype calls missing over all snps with missing.n2<0.05",
  "except that Y chr SNPs are excluded for females")
summary(scanAnnot$missing.e2)
# Make sure the samples are still in order
allequal(scanAnnot$scanID, sort(scanAnnot$scanID))
@ 

\item{We will create a histogram of the overall missing call rate per sample in order to identify any samples with
a relatively larger missing call rate.  It is known that genotype data taken from DNA that has been through whole-genome
amplification (WGA) has an overall higher missing call rate; we could identify that here.
We also look at the summary of the missing call rate for females and males separately to ensure there are no large gender
differences. Finally, we plot the samples sorted by missing call rate.  }

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=> 
hist(scanAnnot$missing.e2, xlab="Fraction of missing calls over all probes 
     with missing call rate < 0.05", 
     main="Histogram of Sample Missing Call Rate for all Samples")
@ 
\caption{Histogram of missing call rate for all samples over all probes with missing call rate less than 0.05.  
There are 896,250 SNPs below this threshold.  Note the highest sample has an overall missing call rate of 2 percent.}
\label{fig:missing_e2}
\end{figure}

<<>>=
# Look at missing.e2 for males
summary(scanAnnot$missing.e2[scanAnnot$sex == "M"])
# Look at missing.e2 for females
summary(scanAnnot$missing.e2[scanAnnot$sex == "F"])
@ 

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=> 
plot(sort(scanAnnot$missing.e2), rank(sort(scanAnnot$missing.e2)), 
     xlab="Fraction of missing calls over probes with MCR < 0.05",
     ylab="rank", cex.lab=0.75,
     main="Sorted Fraction of Missing Call Rate \nOver All Probes with MCR < 0.05")
@ 
\caption{In increasing order, a plot of the fraction of missing calls for each sample over all probes with missing call
rate less than 0.05.}
\label{fig:missing_e2_sort}\end{figure}

\end{itemize}


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

\subsection{Calculate Missing Call Rates by Batch} 

\begin{itemize} 

\item{Next, the missing call rate by batch is calculated to check that there are no batches with comparatively lower
call rates.  We find the distribution of number of samples per batch from the sample annotation file variable {\tt geno.batch}.
We calculate the mean missing call rate for all samples in each of the batches by simply taking the mean 
of all samples in a given batch of the 
missing call rate per sample value, stored in the sample annotation file as {\tt missing.all}. }

<<>>=
varLabels(scanAnnot)
# Check how many batches exist and how many samples are in each batch
sum(is.na(scanAnnot$plate))
length(unique(scanAnnot$plate))
table(scanAnnot$plate, exclude=NULL)
@ 

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
# Plot the distribution of the number of samples per batch.
barplot(table(scanAnnot$plate), 
        ylab="Number of Samples", xlab="Batch", 
        main="Distribution of Samples per Batch", 
        cex.axis=0.85, cex.names=0.75, cex.lab=0.9)
@ 
\caption{The distribution of the samples among the batches.}
\label{fig:samples_per_batch}
\end{figure}

<<>>=
# Examine the mean missing call rate per batch for all SNPs
batches <- unique(scanAnnot$plate)
bmiss <- rep(NA,length(batches)); names(bmiss) <- batches
bn <- rep(NA,length(batches)); names(bn) <- batches
for(i in 1:length(batches)) {
  x <- scanAnnot$missing.e1[is.element(scanAnnot$plate, batches[i])]
  bmiss[i] <- mean(x)
  bn[i] <- length(x)
}
@ 

\item{To find the slope of the regression line from the mean missing call rate per batch regressed on the number
of samples per batch, we will take the results from ANOVA.  Then we can plot the mean missing call rate against the 
number of samples in the batch with the regression line.  For studies with more batches, this test can identify any
batch outliers with regard to missing call rate for samples in a given batch.
We can do the same analysis using the mean missing call rate for autosomal SNPs, or SNPs on the X chromosome in the exact
same way, substituting {\tt missing.e1} with either {\tt miss.e1.auto} or {\tt miss.e1.xchr}.  Because the results are nearly
identical, we will not show them here.}

<<>>=
y <- lm(bmiss ~ bn)
anova(y)
@ 

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
plot(bn, bmiss, 
 xlab="Number of samples per batch", ylab="Mean missing call rate", 
 main="Mean Missing Call Rate vs\nSamples per Batch", cex.main=0.8)
abline(y$coefficients)
@ 
\caption{The number of samples per batch plotted against the mean missing call rate.  The 
regression line has a nominal slope taken from the ANOVA results for the mean missing call 
rate per batch regressed on the number of samples per batch.}
\label{fig:samples_by_miss_rate}
\end{figure}

\end{itemize}

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

\subsection{Chi-Square Test of Allelic Frequency Differences in Batches}

\begin{itemize}

\item{In this step, the chi-square test for differences in allelic frequency is performed between each batch 
individually and a pool of all the other batches in the study.  We then look at the mean $\chi^2$ statistic over all SNPs for each batch as a function of the ethnic composition of samples in a batch. We use the {\tt plate} variable 
in the scan annotation  to identify the samples in each batch,
so we must include the scan annotation in the {\tt GenotypeData} object.
Then we call the function {\tt batchChisqTest} 
which calculates the $\chi^2$ values from 2x2 tables for each SNP,
comparing each plate with the other plates.  This function returns the
genomic inflation factors for each batch, as well as matrix of $\chi^2$ values for
each SNP.}

<<>>=
# Perform the Chi-square test using the batchChisqTest function
res <- batchChisqTest(genoData, batchVar="plate", return.by.snp=TRUE)
close(genoData)
# chi-square values for each SNP
dim(res$chisq)
# genomic inflation factor
res$lambda
# Examine the average chi-square test statistics for each of the batches
res$mean.chisq
@ 

\item{Next we test for association between batches and population groups, using a $\chi^2$ contingency test.  Then we look at the relationship between the ethnic composition of each batch and the previously calculated $\chi^2$ test of allelic frequency between each batch and a pool of the other batches.  The point is to look for batches that differ from others of similar ethnic composition, which might indicate a batch effect due to genotyping artifact.  In this experiment, there are only a few batches and wide variations in ethnicity among batches, so it is difficult to interpret the results.  In larger GWAS experiments, we generally observe a U-shaped curve of allelic frequency test statistic as a function of ethnic composition.}

<<>>=
x <- table(scanAnnot$race, exclude=NULL)
x
x[1] / sum(x)
x[2] / sum(x)
x <- table(scanAnnot$race, scanAnnot$plate)
x
# Run an approximate chi-square test to see if there are ethnic effects
chisq <- chisq.test(x)
chisq$p.value
# Check if ethnicity CEU has a batch effect
# Calculate the fraction of samples in each batch that are CEU ethnicity
batches <- unique(scanAnnot$plate)
length(batches)
eth <- rep(NA,length(batches)); names(eth) <- sort(batches)
for(i in 1:length(batches)){
  x <- scanAnnot$race[is.element(scanAnnot$plate, batches[i])]
  xl <- length(x[x == "CEU"])
  eth[i] <- xl / length(x)
}
allequal(names(eth), names(res$mean.chisq))
@ 

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
# Plot the average Chi-Square test statistic against the 
#     fraction of samples that are CEU
plot(eth, res$mean.chisq, xlab="Fraction of CEU Samples per Batch", 
  ylab="Average Chi-square Test Statistic",
  main="Fraction of CEU Samples per Batch 
  vs Average Chi-square Test Statistic")
abline(v=mean(eth), lty=2, col="red")
@ 
\caption{A plot of the fraction of samples in each batch from the CEU population group versus the average chi-square statistic
for the particular batch.  The red line is the average fraction of samples that are of European ancestry in a batch.}
\label{fig:eth_vs_chi_sq}
\end{figure}

\end{itemize}


\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Sample Quality Checks}

In this step we examine sample quality using three methods.  We check for outliers in genotype quality score;
we check for anomalous sample-chromosome pairs using BAlleleFreq variance analysis; lastly, we check 
sample missingness and heterozygosities.  We use Illumina data, since
the B Allele Frequency and Log R Ratio are included in the dataset.

\subsection{Sample genotype quality scores}


Genotype calling algorithms report quality scores and classify genotypes with insufficient confidence as missing.  
This code calculates the mean and median genotype quality score for each sample.

\begin{itemize}

\item{Calculate quality scores by sample.  The {\tt
      qualityScoreByScan} function requires both an {\tt
      IntensityData} object, to read the quality scores, and a {\tt
      GenotypeData} object, to determine which scans have missing
    genotypes and should be omitted from the calculation.}
  
<<>>=
library(GWASTools)
library(GWASdata)
data(illumina_scan_annot)
scanAnnot <- ScanAnnotationDataFrame(illumina_scan_annot)
qxyfile <- system.file("extdata", "illumina_qxy.nc", 
  package="GWASdata")
qualNC <- NcdfIntensityReader(qxyfile)
qualData <- IntensityData(qualNC, scanAnnot=scanAnnot)
genofile <- system.file("extdata", "illumina_geno.nc", 
  package="GWASdata")
genoNC <- NcdfGenotypeReader(genofile)
genoData <- GenotypeData(genoNC, scanAnnot=scanAnnot)
qual.results <- qualityScoreByScan(qualData, genoData)
close(qualData)
@ 

\begin{figure}[ht!] 
\centering
<<fig=TRUE>>=
hist(qual.results[,"median.quality"], main="Median Genotype Quality Scores 
  of Samples", xlab="Median Quality")
@ 
\caption{Median quality scores for each sample.}
\label{fig:medQualScore}
\end{figure}
\end{itemize}

Figure~\ref{fig:medQualScore} shows the distribution of median quality scores; it is unsurprising that these are all 
good, given that some quality checking happens at the genotyping centers. Clear outliers in this plot would be
cause for concern that the sample(s) in question were of significantly lower quality than the other samples.


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

\subsection{BAlleleFreq variance analysis}

The "BAlleleFreq" was previously calculated for each SNP-sample combination. ``BAlleleFreq'' is a standardized version of the polar coordinate angle. It calculates
the frequency of the B allele within a single sample. Under normal circumstances, the true frequency
is 0, $\frac{1}{2}$, or 1. In cases of allelic imbalance the true frequencies may vary. For example, in a 
population of trisomic cells, the true frequencies would be 0, $\frac{1}{3}$, $\frac{2}{3}$, or 1. Here we calculate the variance of BAF (for SNPs called as heterozygotes) within a sliding window along each chromosome for each sample.  Each chromosome is divided into 12 sections with equal numbers of SNPs and the variance is calculated in a window of two adjacent sections (one-sixth of the chromosome), which slides along the chromosome in increments of one section. Regions (windows)
with very high BAF variance can indicate chromosomal anomalies. 


\subsection*{Calculate the sliding window BAF standard deviation}
This process identifies chromosome-sample pairs that have windows with very high BAlleleFreq
standard deviation, with ``very high'' defined as more than 4 standard deviations from the window's
mean BAlleleFreq standard deviation over all samples. The output is a matrix listing all sample-chromosome 
pairs with high BAlleleFreq standard deviations, the number of windows with high SDs in 
each pair, and the sample's sex. We examine plots of BAlleleFreq by
position for each identified chromosome-sample pair (though only a
subset of plots are shown here).
      
\begin{itemize}

\item{First, run the {\tt meanBAFbyScanChromWindow} function. This
    requires both an {\tt IntensityData} object with BAlleleFreq
    and a {\tt GenotypeData} object. Its output is a list of matrices, with
one matrix for each chromosome containing the standard
deviation of BAlleleFreq at each window in each scan.}
  
<<>>=
blfile <- system.file("extdata", "illumina_bl.nc", 
  package="GWASdata")
blNC <- NcdfIntensityReader(blfile)
blData <- IntensityData(blNC)
nbins <- rep(12, 3)
slidingBAF12 <- sdByScanChromWindow(blData, genoData, nbins=nbins)
length(slidingBAF12)
length(slidingBAF12)
dim(slidingBAF12[[1]])
dim(slidingBAF12[[2]])
@ 

\item{The function {\tt meanBAFSDbyChromWindow} calculates the mean and standard 
deviation of the BAlleleFreq standard deviations in each window in each chromosome over
all samples. For the X chromosome, males and females are calculated separately, and we
save the results split by sex.}

<<>>=
sds.chr <- meanSdByChromWindow(slidingBAF12, 
                               scanAnnot$sex)
sds.chr[["21"]]
sds.chr[["X"]]     
@ 

\item{Next, identify windows within sample-chromosome pairs that have very high BAlleleFreq
standard deviations compared to the same window in other samples. There are 80 sample-chromosome
pairs that are flagged in this analysis, and seem to be evenly distributed across chromosomes.}  

<<>>=
res12bin4sd <- findBAFvariance(sds.chr, slidingBAF12, 
                               scanAnnot$sex, sd.threshold=4)
nrow(res12bin4sd)
head(res12bin4sd)
sum(is.na(res12bin4sd))
table(res12bin4sd[, 2])
all(res12bin4sd[res12bin4sd[, "sex"] == "F", "scanID"] %in% 
    scanAnnot$scanID[scanAnnot$sex == "F"])
@ 

\item{Call {\tt chromIntensityPlot} to plot the BAlleleFreq of all SNPs on the indicated chromosome-sample
pair against position.  This yields many plots that must be individually examined
to distinguish noisy data from chromosomal abnormalities.}  

\begin{figure}[ht!]
\centering
<<fig=TRUE>>= 
scanID <- as.integer(res12bin4sd[, "scanID"])
chrom <- as.integer(res12bin4sd[, "chromosome"])
chrom[res12bin4sd[, "chromosome"] == "X"] <- 23
bincode <- paste("Bin", res12bin4sd[, "bin"], sep = " ")
sexcode <- paste("Sex", res12bin4sd[, "sex"], sep = " ")
code <- paste(bincode, sexcode, sep = ", ")
plotind <- 1
chromIntensityPlot(blData, scanID[plotind], chrom[plotind],
                   code=code[plotind], type="BAF")
@
\label{fig:bafcheck}
\end{figure}

<<>>=
close(blData)
@ 

\end{itemize}

At this stage, we have generated plots of those chromosomes (over all chromosomes and samples) that have unusually 
high BAF standard deviation. The next step in the process is to examine each of these plots to look for anomalies. 
Note that this step will be labor intensive for any reasonably large study. Better methods of detecting BAF anomalies
are an active area of research.

\clearpage

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

\subsection{Missingness and heterozygosity within samples}

This step calculates the percent of missing and heterozygous genotypes in each chromosome of each sample. 
We create boxplots of missingness by individual chromosome, as well as autosomal and X chromosome heterozygosity in each population. 
This allows for identification of samples that may have relatively high heterozygosity for all chromosomes, indicating a possible mixed sample.
Further, we are able to identify any outliers with regard to missingness.
Plotting by chromosome enables visualization of chromosomal artifacts on a particular subset of 
SNPs that lie on a chromosome.

\begin{itemize}

\item{We will call the function {\tt missingGenotypeByScanChrom} to
    calculate the missing call rate.  Since the function returns
    missing counts per chromosome as well as snps per chromosome, we
    divide to find the missing call rate per chromosome.}
  
<<>>=
miss <- missingGenotypeByScanChrom(genoData)
miss.rate <- t(apply(miss$missing.counts, 1, function(x) {
  x / miss$snps.per.chr}))
miss.rate <- as.data.frame(miss.rate)
@ 

\item{First, Figure~\ref{fig:missByChrom} shows a boxplot of missingness in the autosomes, the X chromosome, 
and the pseudoautosomal region.
Second, Figure~\ref{fig:XChromMissSex} shows a boxplot of X chromosome missingness for each sex.}

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
cols <- names(miss.rate) %in% c(1:22, "X", "XY")
boxplot(miss.rate[,cols], main="Missingness by Chromosome", 
  ylab="Proportion Missing", xlab="Chromosome")
@
\caption{Boxplot of missingness by chromosome.}
\label{fig:missByChrom}
\end{figure}

\begin{figure}[ht!]
<<fig=TRUE>>=
boxplot(miss.rate$X ~ scanAnnot$sex, 
  main="X Chromosome Missingness by Sex", 
  ylab="Proportion Missing")
@ 
\centering
\caption{Boxplot of X chromosome missingness by sex.}
\label{fig:XChromMissSex}
\end{figure}
 
\item{We will call the function {\tt hetByScanChrom} to calculate the heterozygosity.
  We want to store the heterozygosity 
calculations in the sample annotation, so after ensuring the ordering is correct we will write the heterozygosity
values and store them in the sample annotation.}
 
<<>>=
# Calculate heterozygosity by scan by chromosome
het.results <- hetByScanChrom(genoData)
close(genoData)
# Ensure heterozygosity results are ordered correctly
allequal(scanAnnot$scanID, rownames(het.results))
# Write autosomal and X chr heterozygosity to sample annot
scanAnnot$het.A <- het.results[,"A"]
scanAnnot$het.X <- het.results[,"X"]
varMetadata(scanAnnot)["het.A", "labelDescription"] <- 
  "fraction of heterozygotes for autosomal SNPs"
varMetadata(scanAnnot)["het.X", "labelDescription"] <- 
  "fraction of heterozygotes for X chromosome SNPs"
@ 

\item{There are two plots for heterozygosity. Figure~\ref{fig:autoMiss} shows a boxplot of heterozygosity over the 
autosomes, subsetted by population. We recommend examining BAF plots for high heterozygosity outliers, to look for evidence of sample contamination (more than 3 bands on all chromosomes).  Examination of low heterozygosity samples may also identify chromosomal anomalies with wide splits in the intermediate BAF band.
Figure~\ref{fig:XChromMissPop} shows a boxplot of female heterozygosity on the X chromosome, subsetted by 
population.}

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
boxplot(scanAnnot$het.A ~ scanAnnot$race, 
  main="Autosomal Heterozygosity")
@ 
\caption{Boxplot of autosomal heterozygosity by continental ancestry.}
\label{fig:autoMiss}
\end{figure}

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
male <- scanAnnot$sex == "M"
boxplot(scanAnnot$het.A[male] ~ scanAnnot$race[male], 
  main="X Chromosome Heterozygosity in Females")
@
\caption{Boxplot of X chromosome heterozygosity by population.}
\label{fig:XChromMissPop}
\end{figure}

\end{itemize}
 

\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Sample Identity Checks}

This step performs a series of identity checks on the samples. First, samples are analyzed to determine if 
there exist any gender discrepancies between the annotated and genetic genders in the sample. 
Next, the relatedness among samples is investigated through IBD estimation. Finally, the samples are checked 
for potential population substructure, which if unidentified can threaten the validity of subsequent analyses.  

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



\subsection{Mis-annotated Gender Check}

This section looks for discrepancies between the annotated and genetic genders. Gender identity is usually inferred 
from X chromosome heterozygosity, but our experience is that this variable can give ambiguous results when used 
alone (for example, in XXY males or due to genotyping artifacts). Plots of the mean allelic intensities of SNPs on the X and Y chromosomes can identify mis-annotated gender as well as sex chromosome aneuploidies. It is important to have accurate 
gender annotation not only for completeness but also for analyses which treat male and female samples separately. 
Any found gender mis-annotations are presented to the investigators in order to resolve discrepancies. If a genetic and recorded gender do not match, 
a collective decision must be made regarding the inclusion of those genetic data. In some cases a recording error explains the discrepancy, but more often the discrepancy is unexplained. These cases are assumed to be a sample mis-identification and these samples are excluded from subsequent analyses.


\begin{itemize}
  
\item{In order to compare the mean X and Y chromosome intensities for all samples, we must calculate the mean intensity
for each sample by chromosome. The function {\tt meanIntensityByScanChrom} calculates for each sample the mean
and standard deviation of the sum of the two allelic intensities for each probe on a given chromosome.  
A matrix with one row per sample and one column per chromosome with entries $[i,j]$ corresponding to either the mean or 
standard deviation of all probe intensities for the $i^{\text th}$ sample and the $j^{\text th}$ chromosome is returned 
from the function.}

<<>>=
library(GWASTools)
library(GWASdata)
qxyfile <- system.file("extdata", "affy_qxy.nc",
  package="GWASdata")
intenNC <- NcdfIntensityReader(qxyfile)
inten.by.chrom <- meanIntensityByScanChrom(intenNC)
close(intenNC)
length(inten.by.chrom)
names(inten.by.chrom)
@ 

\item{Now we will use the calculated mean intensities by sample to identify any gender mis-annotation or sex 
chromosome aneuploidies. For the plots, we will create a color coding corresponding to the annotated gender, with blue 
for males and red for females. We will also load in the SNP annotation file to find the probe counts for the X and Y 
chromosomes; we use these in the plot axis labels.  }

<<>>=
mninten <- inten.by.chrom[[1]]  # mean intensities
dim(mninten)
data(affy_scan_annot)
scanAnnot <- ScanAnnotationDataFrame(affy_scan_annot)
# Check to be sure sample ordering is consistent
allequal(scanAnnot$scanID, rownames(mninten))
# Assign each gender a color
xcol <- rep(NA, nrow(scanAnnot))
xcol[scanAnnot$sex == "M"] <- "blue"
xcol[scanAnnot$sex == "F"] <- "red"
data(affy_snp_annot)
snpAnnot <- SnpAnnotationDataFrame(affy_snp_annot)
nx <- sum(snpAnnot$chromosome == 23)
ny <- sum(snpAnnot$chromosome == 25)
@ 

\item{For two of the plots we will create next, we use the autosome and X chromosome heterozygosity values calculated
in an earlier step and stored in the sample annotation file. It is assumed the sample 
annotation file holds these values. Four plots will now be created: mean X chromosome intensity versus mean Y chromosome 
intensity, mean X chromosome intensity versus X chromosome heterozygosity, mean X chromosome heterozygosity versus 
mean Y chromosome intensity and mean autosomal heterozygosity versus mean X chromosome heterozygosity. The fourth plot 
applies to annotated females only, since males are expected to have zero heterozygosity on the X chromosome.}

<<results=hide>>=
#All intensities
x1 <-mninten[,"X"]; y1 <- mninten[,"Y"]
main1 <- "Mean X vs \nMean Y Chromosome Intensity"
#Het on X vs X intensity
x2 <- mninten[,"X"]; y2 <- scanAnnot$het.X
main2 <- "Mean X Chromosome vs
  Mean X Chromosome Heterozygosity"
# Het on X vs Y intensity
y3 <- mninten[,"Y"]; x3 <- scanAnnot$het.X
main3 <- "Mean X Chromosome Heterozygosity vs
  Mean Y Chromosome Intensity"
# X vs A het
x4 <- scanAnnot$het.A[scanAnnot$sex == "F"]
y4 <- scanAnnot$het.X[scanAnnot$sex == "F"]
main4 <- "Mean Autosomal Heterozygosity vs
  Mean X Chromosome Heterozygosity"
cols <- c("blue","red")
mf <- c("male", "female")
xintenlab <- paste("X intensity (n=", nx, ")", sep="")
yintenlab <- paste("Y intensity (n=", ny, ")", sep="")
pdf("DataCleaning-gender.pdf")
par(mfrow=c(2,2))
plot(x1, y1, xlab=xintenlab, ylab=yintenlab, 
  main=main1, col=xcol, cex.main=0.8)
legend("topright",mf,col=cols,pch=c(1,1),cex=0.8)
plot(x2, y2, col=xcol, xlab=xintenlab, 
  ylab="X heterozygosity", main=main2, cex.main=0.8)
legend("topleft", mf, col=cols, pch=c(1,1),cex=0.8)
plot(x3, y3, col=xcol, ylab=yintenlab, 
  xlab="X heterozygosity", main=main3, cex.main=0.8)
legend("topright", mf, col=cols, pch=c(1,1),cex=0.8)
plot(x4,y4, col="red", xlab="Autosomal heterozygosity", 
  ylab="X heterozygosity", main=main4, cex.main=0.8)
legend("topleft", mf, col=cols, pch=c(1,1),cex=0.8)
dev.off()
@ 

\begin{figure}[ht!]
\centering
\includegraphics[width=5.67in,height=5.67in,keepaspectratio]{DataCleaning-gender}
\caption{Intensity plots for all samples, with blue indicating an annotated male and red
corresponding to annotated females}
\label{fig:inten_plots}
\end{figure}

\end{itemize}



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

\subsection{Relatedness and IBD Estimation}

In most studies, there are discrepancies between pedigrees provided and relatedness inferred from the genotype data.  To infer genetic relatedness, we estimate coefficients of identity by descent (IBD). It is important to identify and record unannotated 
relationships so that analyses assuming all subjects are unrelated can use a filtered subset of samples.
From our experience, it is difficult to accurately estimate low levels of relatedness, but higher levels can be more
reliably determined. Users are encouraged to employ analyses which take into accounts the IBD estimates themselves
rather than discrete relationship coefficients for any relationships. This step calculates identity by descent (IBD) 
estimates for each population group using an expectation-maximization (EM) approach for maximum likelihood 
estimation (MLE) of the IBD probabilities {\tt k1} and {\tt k0} 
\footnote{Weir, B. S., Anderson, A. D. \& Hepler, A. B. Genetic relatedness analysis: modern data and new challenges.
\textit{Nature Reviews Genetics} \textbf{7}, 771-780 (October 2006).}, 
and then graphs {\tt k1} versus {\tt k0}. The calculations are done
using the {\tt gdsfmt} and {\tt SNPRelate} packages developed as part 
of the GENEVA project.   {\tt gdsfmt} is a library that uses a file
format similar to the NetCDF files, and {\tt SNPRelate} performs the
IBD calculations using this data format.
                             

\subsubsection*{IBD Using Maximum Likelihood Estimation (MLE)}                             
      
\begin{itemize}
\item{We assume the sample annotation file has the {\tt duplicated} variable, while the SNP file needs to have 
{\tt chromosome} and {\tt missing.n1} variables. We will convert the genotype NetCDF file to the {\tt gds} format with the function {\tt convertNcdfGds}. For more function information, please refer to the function 
documentation.}

<<>>=
library(gdsfmt)
library(SNPRelate)
ncfile <- system.file("extdata", "affy_geno.nc",
  package="GWASdata")
gdsfile <- "broad.hapmap.geno.gds"
convertNcdfGds(ncfile, gdsfile)
@ 

\item{SNP selection for the IBD estimation analysis has proved to be a salient component of the analysis.  
One may consider choosing SNPs based upon patterns of linkage, minor allele frequency or position along the genome.  
Here we simply use autosomal SNPs with missing call rate less than 0.05 to include in the
MLE IBD estimation.}

<<>>=
snp.auto <- pData(snpAnnot)[snpAnnot$chromosome < 23, ]
snp.auto.nonmiss <- snp.auto[snp.auto$missing.n1 < 0.05, ]
snp.select <- snp.auto.nonmiss$snpID
length(snp.select)
@ 

\item{Next, we perform the IBD estimation for the CEU population group.
The {\tt snpgdsIBDMLE} function gives square matrices with rows and columns equal to the number of
samples where entries are the {\tt k0} and {\tt k1} coefficient estimates, respectively, for each pairwise combination 
of samples. Note for larger studies this is a significant analysis since every possible pair of samples must be 
considered.}

<<>>=
sample.sel <- scanAnnot$scanID[scanAnnot$race == "CEU"]
length(sample.sel)
gdsobj <- openfn.gds(gdsfile)
ibd <- snpgdsIBDMLE(gdsobj, sample.id=sample.sel,
                    snp.id=snp.select, method="EM")
closefn.gds(gdsobj)
dim(ibd$k0)
ibd$k0[1:5,1:5]
@ 

\item{We can now find the expected relationships between samples in the CEU group based on the pedigree data
that is stored in the sample annotation file. We will create a subset of the sample annotation that has one line per 
sample and columns that hold family, father and mother ids, where an entry of {\tt 0} indicates no familial data.  
Then the function {\tt pedigreeClean} is called, which checks for basic pedigree errors such as not all mothers annotated as 
female or not all fathers as male, or that there exist no individual ids that appear in both the mother and father columns.
Please see the function documentation for more information on this function. }

<<>>=
samp <- pData(scanAnnot)[scanAnnot$race == "CEU",
  c("family", "subjectID", "father", "mother", "sex")]
dim(samp)
names(samp) <- c("family", "individ", "father", "mother", "sex")
samp[1:5, ]
pedigreeClean(samp)
@ 

\item{The functions that determine expected relationships require no duplicates in the pedigree. This will be verified by a
function call to {\tt pedigreeFindDuplicates}. We will also call {\tt pedigreeCheck} which determines if there are any 
singleton families, mothers/fathers whose sex does not match,
impossible relationships, or subfamilies.}

<<>>=
dups <- pedigreeFindDuplicates(samp)
dups
uni.samp <- pedigreeDeleteDuplicates(samp, dups$dups.match)
fc <- pedigreeCheck(uni.samp)
@ 

\item{The output is a list consisting of  vectors and a matrix: 
{\tt one.person} is vector of family ids for one-person families. 
{\tt mismatch.sex} consists of family ids where sex of mother and/or father is incorrect.
{\tt impossible.related} consists of family ids where either child is mother of self or 
an individual is both child and mother of same person subfamilies. {\tt ident} is a matrix
with family id (families with subfamilies), subfamily identifier, individual ids of
persons in the subfamily. (In {\tt subfamilies.ident} the individual id's include 
individuals identified as mother or father who may not be in {\tt individ}.) }
 
<<>>=
length(fc$one.person)
length(fc$mismatch.sex)
length(fc$impossible.related)
length(fc$subfamilies.ident)
@ 

\item{There are multiple subfamilies identified, so we will need to
    assign new family IDs to the subfamililes.}
  
<<>>=
subf <- fc$subfamilies.ident
head(subf)
table(subf$family)
subf.ids <- subf$individ[subf$subfamily == 2]
newfam <- uni.samp$individ %in% subf.ids
uni.samp$family[newfam] <- uni.samp$family[newfam] + 10000
table(uni.samp$family)
pedigreeCheck(uni.samp)
@ 

\item{ Now from the verified sample list excluding duplicate samples, we can calculate the expected relationships
among the samples by calling the function {\tt pedigreePairwiseRelatedness}. The relationships looked for as annotated are:
UN = unrelated, PO = parent/offspring, FS = full siblings, HS = half siblings, and FC = first cousins.
Families where mothers and fathers are related are also looked for among the family annotations.
It is this function output with annotated expected relationships that we will compare with the relationship estimates
calculated earlier and saved in the R object {\tt ibd}. }

<<>>=
rels <- pedigreePairwiseRelatedness(uni.samp)
length(rels$inbred.fam)
relprs <- rels$relativeprs
relprs[1:5,]
table(relprs$relation)
@ 

\item{In order to plot the IBD coefficient estimates color coded by expected relationships, some data manipulation must
occur. The samples must be coded in terms of subject id and each pair of samples must be annotated with the expected 
relationship.}

<<>>=
relprs$s12 <- paste(relprs$Individ1, relprs$Individ2)
ibd$subjectID <- scanAnnot$subjectID[is.element(
  scanAnnot$scanID, ibd$sample.id)]
sample1.temp <- rep(ibd$subjectID, length(ibd$subjectID))
sample1.temp <- matrix(sample1.temp, nrow = dim(ibd$k0)[1])
sample2.temp <- t(sample1.temp)
ibd$sample1.subjectID <- as.array(sample1.temp)
ibd$sample2.subjectID <- as.array(sample2.temp)
ibd$s12 <- paste(ibd$sample1.subjectID, ibd$sample2.subjectID)
ibd$s21 <- paste(ibd$sample2.subjectID, ibd$sample1.subjectID)
ibd$rel <- array("U", length(ibd$s12))
ibd$rel[is.element(ibd$s12, relprs$s12[is.element(relprs$relation, 
  "PO")]) | is.element(ibd$s21, relprs$s12[is.element(
  relprs$relation, "PO")])] <- "PO"
ibd$rel[ibd$sample1.subjectID == ibd$sample2.subjectID] <- "Dup"
table(ibd$rel, exclude=NULL)
@ 

\item{ Now the pedigree information is in the proper format for the 
IBD estimates to be plotted for each pair of samples, color coded by expected relationship, see Figure~\ref{fig:ibd_ceu}. 
The orange line segments span the values for the expected {\tt k0}, {\tt k1} values, $\pm$2 standard deviations. }

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
ibdPlot(ibd$k0, ibd$k1, relation=ibd$rel,
  main="HapMap CEU IBD Estimates \nEM Algorithm")
@ 
\caption{IBD estimates for the CEU samples.}
\label{fig:ibd_ceu}
\end{figure}


\end{itemize}


\subsubsection*{IBD Using Method of Moments (MoM)}                             

IBD analysis can also be performed using the PLINK Method of Moments (MoM) approach. For the MoM approach, we use
a selection of autosomal SNPs with missing call rate less than 5
percent and minor allele frequency $>0$ that are spaced at least 15 kb apart.

\begin{itemize}
\item{We calculate the allele frequency in the set of samples that will be
analyzed (in this case, YRI samples).
The function {\tt apartSnpSelection} randomly selects SNPs for which no pair of SNPs is closer than a user specified
threshold; we choose 15 kb here.}

<<>>=
# Allele frequency calculation in YRI samples
# open the NetCDF file and create a GenotypeData object
nc <- NcdfGenotypeReader(ncfile)
genoData <-  GenotypeData(nc, scanAnnot=scanAnnot, snpAnnot=snpAnnot)
# sample selection - unduplicated YRI samples
sample.excl <- scanAnnot$scanID[scanAnnot$race != "YRI" |
                                scanAnnot$duplicated]
length(sample.excl)
afreq <- alleleFrequency(genoData, scan.exclude=sample.excl)
close(genoData)
# Add allele frequency to SNP annotation
afreq <- afreq[,"all"] # males and females combined
snpAnnot$A.freq.yri <- afreq
varMetadata(snpAnnot)["A.freq.yri", "labelDescription"] <- 
  "frequency of allele A in YRI subjects"
# Initial selection: autosomal SNPs with low missing call rate and MAF>0
init.sel <- snpAnnot$missing.n1 < 0.05 & snpAnnot$chromosome < 23 &
  afreq > 0 & afreq < 1
sum(init.sel)
# Call the function to select SNPs no closer than 15000 bases
set.seed(1001)
n <- 15000
rsnp2 <- apartSnpSelection(snpAnnot$chromosome, snpAnnot$position, 
                           min.dist=n, init.sel=init.sel)
sum(rsnp2)
# Add the SNP selection to the SNP annotation file
snpAnnot$r15kb <- rsnp2
varMetadata(snpAnnot)["r15kb", "labelDescription"] <- paste(
  "Logical indicator for autosomal SNPs with missing.n1<0.05",
  "and YRI MAF>0, at least 15 kb apart")
@ 

\item{We now use the function {\tt snpgdsIBDMoM} to perform IBD
    for the YRI samples using the PLINK MoM approach.}
  
<<>>=
sample.sel <- scanAnnot$scanID[scanAnnot$race == "YRI"]
snp.select <- snpAnnot$snpID[snpAnnot$r15kb]
gdsobj <- openfn.gds(gdsfile)
ibd <- snpgdsIBDMoM(gdsobj, sample.id=sample.sel,
                      snp.id=snp.select)
closefn.gds(gdsobj)
dim(ibd$k0)
ibd$k0[1:5,1:5]
@ 

\item{We can now find the expected relationships between samples in the YRI group based on the pedigree data
that is stored in the sample annotation file. 
We will create a subset of the sample annotation that has one line per sample and columns that hold family, father
and mother ids, where an entry of {\tt 0} indicates no familial data.
As before, we use the function {\tt pedigreeClean} to check for
pedigree errors.}

<<>>=
samp <- pData(scanAnnot)[scanAnnot$race == "YRI",
  c("family", "subjectID", "father", "mother", "sex")]
dim(samp)
names(samp) <- c("family", "individ", "father", "mother", "sex")
samp[1:5, ]
pedigreeClean(samp)
@ 

\item{As with the Europeans, we use {\tt pedigreeFindDuplicates} to
    identify duplicates and and {\tt pedigreeCheck} to look for
singleton families, mothers/fathers whose sex does not match,
impossible relationships, or subfamilies.
In studies with family data, sorting pedigree inconsistencies
is a convoluted task and it is important to utilize automated methods to help detect any inconsistencies. }

<<>>=
dups <- pedigreeFindDuplicates(samp)
dups
uni.samp <- pedigreeDeleteDuplicates(samp, dups$dups.match)
fc <- pedigreeCheck(uni.samp)
(subf <- fc$subfamilies.ident)
@ 

\item{There is a subfamily with two unrelated people (likely
    founders), so we remove this family from the pedigree.}

<<>>=
uni.samp <- uni.samp[!(uni.samp$family %in% subf$family),]
table(uni.samp$family)
pedigreeCheck(uni.samp)
@ 

\item{Now from the verified sample list excluding duplicate samples, we can calculate the expected relationships
among the samples by calling the function {\tt pedigreePairwiseRelatedness}.}

<<>>=
rels <- pedigreePairwiseRelatedness(uni.samp)
length(rels$inbred.fam)
relprs <- rels$relativeprs
relprs[1:5,]
table(relprs$relation)
relprs$s12 <- paste(relprs$Individ1, relprs$Individ2)
ibd$subjectID <- scanAnnot$subjectID[is.element(
  scanAnnot$scanID, ibd$sample.id)]
sample1.temp <- rep(ibd$subjectID, length(ibd$subjectID))
sample1.temp <- matrix(sample1.temp, nrow = dim(ibd$k0)[1])
sample2.temp <- t(sample1.temp)
ibd$sample1.subjectID <- as.array(sample1.temp)
ibd$sample2.subjectID <- as.array(sample2.temp)
ibd$s12 <- paste(ibd$sample1.subjectID, ibd$sample2.subjectID)
ibd$s21 <- paste(ibd$sample2.subjectID, ibd$sample1.subjectID)
ibd$rel <- array("U", length(ibd$s12))
ibd$rel[is.element(ibd$s12, relprs$s12[is.element(relprs$relation, 
  "PO")]) | is.element(ibd$s21, relprs$s12[is.element(
  relprs$relation, "PO")])] <- "PO"
ibd$rel[ibd$sample1.subjectID == ibd$sample2.subjectID] <- "Dup"
table(ibd$rel, exclude=NULL)
@ 

\item{ Now the pedigree information is in the proper format for the 
IBD estimates to be plotted for each pair of samples, color coded by expected relationship, see Figure~\ref{fig:ibd_yri}. 
Several samples in the YRI population show higher than expected relatedness.\footnote{A 
second generation human haplotype map of over 3.1 million SNPs, supplementary information. The International HapMap
Consortium, \textit{Nature} \textbf{449}, 851-861 (18 October 2007).} }

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
ibdPlot(ibd$k0, ibd$k1, relation=ibd$rel,
  main="HapMap YRI IBD Estimates \nPLINK MOM Approach")
@ 
\caption{Method of Moment (MoM) estimates for the YRI samples. The yellow bars are the expected positions (+/- 2 SD) for full siblings, half siblings and first cousins.}
\label{fig:ibd_yri}
\end{figure}

\end{itemize}


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

\subsection{Population Structure}

\subsubsection*{Principal Component Analysis on all ethnic groups}
In this section, we perform principal component analysis (PCA) in order to detect any population substructure that 
may exist among samples in a study. After calculating the eigenvectors for the samples, we plot the values for each of the first 4 
eigenvectors in a pairwise fashion for each individual. By color coding the plots by annotated ethnicity, we can 
identify any individuals whose recorded self-identified ethnicity differs from their inferred genetic ancestry. Further, we can use the PCA-identified continental ancestry when stratifying 
samples by population group. It may also be useful to include the values of some eigenvectors as covariates in 
association tests.
     
\begin{itemize}

\item{We must load the SNP and sample annotation files to select autosomal SNPs with values in {\tt missing.n1} less 
than 0.05. Here we use all SNPs that meet these criteria, but it is generally advisable to select a set of SNPs within which each pair has a low level of linkage disequilibrium (e.g.\ $r^2<0.2$). We must also ensure no duplicate samples are used for the principal component calculations. A color code 
based upon sample continental ancestry is assigned for plotting; CEU is coded as black, YRI as red. The {\tt snpgdsPCA} function is called with the SNP and sample subsets to calculate the first 32 
eigenvectors. To read more about this function, please see the function documentation. }
 
<<>>=
snp.sel <- snpAnnot$snpID[snpAnnot$missing.n1 < 0.05 & 
                          snpAnnot$chromosome < 23]
sample.sel <- scanAnnot$scanID[scanAnnot$duplicated == FALSE]
gdsobj <- openfn.gds(gdsfile)
pca <- snpgdsPCA(gdsobj, sample.id=sample.sel, 
                 snp.id=snp.sel)
closefn.gds(gdsobj)
names(pca)
length(pca$eigenval)
dim(pca$eigenvect)
@ 

\item{We will make a pairs plot showing the first four eigenvectors. A simple calculation is made to find the fraction
of variance among the samples as explained by each eigenvector. }

<<>>=
# Calculate the percentage of variance explained 
# by each principal component.  
pc.frac <- pca$eigenval/sum(pca$eigenval)
lbls <- paste("EV", 1:4, "\n", format(pc.frac[1:4], digits=2), sep="")
samp <- pData(scanAnnot)[match(pca$sample.id, scanAnnot$scanID),]
cols <- rep(NA, nrow(samp))
cols[samp$race == "CEU"] <- "black"
cols[samp$race == "YRI"] <- "red"
@ 

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
pairs(pca$eigenvect[,1:4], col=cols, labels=lbls,
  main = "CEU: black, YRI: red")
@ 
\caption{First four eigenvectors calculated from the principal components analysis.}
\label{fig:pca_plot}
\end{figure}


\subsubsection*{Parallel Coordinates Plot}

\item{A handy method of visualizing the effects of eigenvectors on clusters 
for a principal components analysis is the parallel coordinates plot. 
Figure~\ref{fig:par_coord_plot} shows this plot for all samples, colored by ethnicity. 
The genetic diversity in the YRI group is apparent in the later eigenvectors, while
the remaining groups remain in clusters throughout.}

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
par.coord <- pca$eigenvect
rangel <- apply(par.coord, 2, function(x) range(x)[1])
rangeh <- apply(par.coord, 2, function(x) range(x)[2])
std.coord <- par.coord
for (i in 1:14)
  std.coord[,i] <- (par.coord[,i] - rangel[i])/(rangeh[i]-rangel[i]) 
plot(c(0,15), c(0,1), type = 'n', axes = FALSE, ylab = "", xlab = "", 
 main = "Parallel Coordinates Plot
 CEU: black, YRI: red")
for (j in 1:13)
  for (i in sample(1:nrow(std.coord)) )
    lines(c(j,j+1), std.coord[i,c(j,j+1)], col=cols[i], lwd=0.25)
axis(1, at = 1:14, labels = paste("PC",1:14, sep = "."))
@ 
\caption{Parallel coordinates plot of first 14 eigenvectors, colored by ethnicity.}
\label{fig:par_coord_plot}
\end{figure}

\end{itemize}

<<echo=FALSE, results=hide>>=
file.remove(gdsfile)
@ 


\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Case-Control Confounding}

We recommend checking for case-control confounding as part of the data
cleaning process for GWAS.  This involves checking both the principal
components and the missing call rate for a relationship with case status.

\subsection{Principal Components Differences}
\label{sec:pca_case_cntrl}

This step examines differences in principal components according to case-control status.

\begin{itemize}
\item{Collate PCA information with sample number, case-control status, and population group.}

<<>>=
princomp <- as.data.frame(pca$eigenvect)
samples.nodup <- pData(scanAnnot)[!scanAnnot$duplicated,]
princomp$scanID <- as.factor(samples.nodup$scanID)
princomp$case.ctrl.status <- as.factor(samples.nodup$status)
princomp$race <- as.factor(samples.nodup$race)
@ 

\item{The code below gives what percent of variation is accounted for by the principal component for the first
32 PCs.  It is clear to see the first two eigenvectors hold the largest percentage of variance among the population,
although the total variance accounted for is still less the one-quarter of the total.}

<<>>=
pc.percent <- 100 * pca$eigenval[1:32]/sum(pca$eigenval)
pc.percent
lbls <- paste("EV", 1:3, "\n", format(pc.percent[1:3], digits=2), "%", sep="")
table(samples.nodup$status)
cols <- rep(NA, nrow(samples.nodup))
cols[samples.nodup$status == 1] <- "green"
cols[samples.nodup$status == 0] <- "red"
@ 

\item{Plot the principal component pairs for the first three PCs, by
case-control status. }

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
pairs(pca$eigenvect[,1:3], col=cols, labels=lbls,
  main = "First Three EVs by Case-Control Status")
@ 
\caption{The values from the first three eigenvectors plotted for each sample, color coded by case-control status.  
Green denotes case status and red indicates a control.}
\end{figure}

\item{Do boxplots for the first few PCs to show differences between cases and controls, along with
a two-factor ANOVA accounting for case-control status and population group. Since we are using
randomized case-control status, we do not expect to see a significant difference in principal 
components between cases and controls, when considering population group. }

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
boxplot(princomp[, 1] ~ princomp$case.ctrl.status, 
  ylab = "PC1", main = "PC1 vs. Case-control Status")
@ 
\caption{A boxplot of the values of the first eigenvector stratified by case-control status, where 0 indicates a 
control and 1 indicates a case.}
\label{fig:pca1_casectrl}
\end{figure}

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
boxplot(princomp[, 2] ~ princomp$case.ctrl.status, 
  ylab = "PC2", main = "PC2 vs. Case-control Status")
@ 
\caption{A boxplot of the values of the second eigenvector stratified by case-control status, where 0 indicates a 
control and 1 indicates a case.}
\label{fig:pca2_casectrl}
\end{figure}

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
boxplot(princomp[, 3] ~ princomp$case.ctrl.status, 
  ylab = "PC3", main = "PC3 vs. Case-control Status")
@ 
\caption{A boxplot of the values of the third eigenvector stratified by case-control status, where 0 indicates a 
control and 1 indicates a case.}
\label{fig:pca3_casectrl}
\end{figure}

<<>>=
aov.p1 <- aov(princomp[,1] ~ princomp$race * 
  princomp$case.ctrl.status, princomp)
summary(aov.p1)
aov.p2 <- aov(princomp[,2] ~ princomp$race * 
  princomp$case.ctrl.status, princomp)
summary(aov.p2)
aov.p3 <- aov(princomp[,3] ~ princomp$race * 
  princomp$case.ctrl.status, princomp)
summary(aov.p3)
@ 

\end{itemize}


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

\subsection{Missing Call Rate Differences}

This step determines whether there are differences in missing call rates between
cases and controls. As in section~\ref{sec:pca_case_cntrl}, we use simulated case-control status to demonstrate this step,
since the HapMap II data does not contain information
on cases and controls. 

\begin{itemize}
\item{Investigate the difference in mean missing call rate by case-control status, using the sample annotation variable
{\tt missing.e1}. Here, since the case-control status was randomly assigned, we do not expect to see a difference
in any of the missing call rates with respect to case-control status. }

<<>>=
lm.all <- lm(scanAnnot$missing.e1 ~ scanAnnot$status)
summary(aov(lm.all))
@ 

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
boxplot(scanAnnot$missing.e1 ~ scanAnnot$status, ylab = 
  "Mean missing call rate", main="Mean missing call rate by case status")
@ 
\caption{A boxplot showing the mean missing call rate stratified by case-control status, where 0 indicates a 
control and 1 indicates a case.}
\label{fig:miss_call_rate_vs_case_status}
\end{figure}

\end{itemize}



\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Chromosome Anomaly Detection}

This step looks for large chromosomal anomalies that may be filtered
out during the final analysis.

\subsection{B Allele Frequency filtering}

\begin{itemize}

\item{Create an \texttt{IntensityData} object and a \texttt{GenotypeData} object.}

<<>>=
library(GWASTools)
library(GWASdata)
data(illumina_scan_annot)
scanAnnot <- ScanAnnotationDataFrame(illumina_scan_annot)
data(illumina_snp_annot)
snpAnnot <- SnpAnnotationDataFrame(illumina_snp_annot)

blfile <- system.file("extdata", "illumina_bl.nc", package="GWASdata")
blnc <- NcdfIntensityReader(blfile)
blData <-  IntensityData(blnc, scanAnnot=scanAnnot, snpAnnot=snpAnnot)

genofile <- system.file("extdata", "illumina_geno.nc", package="GWASdata")
genonc <- NcdfGenotypeReader(genofile)
genoData <-  GenotypeData(genonc, scanAnnot=scanAnnot, snpAnnot=snpAnnot)
@ 

\item{Identify some low quality samples by looking at the standard deviation of BAF.}

<<>>=
baf.sd <- sdByScanChromWindow(blData, genoData, var="BAlleleFreq")
med.baf.sd <- medianSdOverAutosomes(baf.sd)
low.qual.ids <- med.baf.sd$scanID[med.baf.sd$med.sd > 0.05]
@ 

\item{Decide which SNPs to exclude based on genome build.}

<<>>=
chrom <- snpAnnot$chromosome
pos <- snpAnnot$position
build <- getAttribute(blnc, attname="genome_build")
if (build == 36) {
  data(HLA.hg18)
  hla <- chrom == 6 & pos >= HLA.hg18$start.base & pos <= HLA.hg18$end.base
  data(pseudoautosomal.hg18)
  xtr <- chrom == 23 & pos >= pseudoautosomal.hg18["X.XTR", "start.base"] & 
    pos <= pseudoautosomal.hg18["X.XTR", "end.base"]
  data(centromeres.hg18)
  centromeres <- centromeres.hg18
} else if (build  ==  37) {
  data(HLA.hg19)
  hla <- chrom == 6 & pos >= HLA.hg19$start.base & pos <= HLA.hg19$end.base
  data(pseudoautosomal.hg19)
  xtr <- chrom == 23 & pos >= pseudoautosomal.hg19["X.XTR", "start.base"] & 
    pos <= pseudoautosomal.hg19["X.XTR", "end.base"]
  data(centromeres.hg19)
  centromeres <- centromeres.hg19
}
ignore <- snpAnnot$missing.n1 == 1 #ignore includes intensity-only and failed snps
snp.exclude <- ignore | hla | xtr
snp.ok <- snpAnnot$snpID[!snp.exclude]
@ 

\item{We use circular binary segmentation to find change points in B Allele
Frequency (BAF).}

<<>>=
scan.ids <- scanAnnot$scanID[1:10]
chrom.ids <- 21:23
baf.seg <- anomSegmentBAF(blData, genoData, scan.ids=scan.ids,
  chrom.ids=chrom.ids, snp.ids=snp.ok, verbose=FALSE)
head(baf.seg)
@ 

\item{Filter segments to detect anomalies, treating the low quality samples differently.}

<<>>=
baf.anom <- anomFilterBAF(blData, genoData, segments=baf.seg, 
  snp.ids=snp.ok, centromere=centromeres, low.qual.ids=low.qual.ids,
  verbose=FALSE)
names(baf.anom)
baf.filt <- baf.anom$filtered
head(baf.filt)
@ 

\end{itemize}

\subsection{Loss of Heterozygosity}

\begin{itemize}

\item{We look for Loss of Heterozygosity (LOH) anomalies by identifying 
homozygous runs with change in LogRRatio.
Change points in LogRRatio are found by circular binary segmentation.
Known anomalies from the BAF detection are excluded. }


<<>>=
loh.anom <- anomDetectLOH(blData, genoData, scan.ids=scan.ids,
  chrom.ids=chrom.ids, snp.ids=snp.ok, known.anoms=baf.filt,
  verbose=FALSE)
names(loh.anom)
loh.filt <- loh.anom$filtered
head(loh.filt)
@ 
\end{itemize}

\subsection{Statistics}

\begin{itemize}

\item{Calculate statistics for the anomalous segments found with the BAF
and LOH methods.}

<<>>=
# create required data frame
baf.filt$method <- "BAF"
if (!is.null(loh.filt)) {
  loh.filt$method <- "LOH"
  cols <- intersect(names(baf.filt), names(loh.filt))
  anoms <- rbind(baf.filt[,cols], loh.filt[,cols])
} else {
  anoms <- baf.filt
}
anoms$anom.id <- 1:nrow(anoms)

stats <- anomSegStats(blData, genoData, snp.ids=snp.ok, anom=anoms,
                      centromere=centromeres)
names(stats)
@ 

\item{Plot the anomalies with relevant statistics, one anomaly per plot. Each plot has
two parts: upper part is a graph of LRR and lower part is a graph of BAF.  }

<<fig=TRUE>>=
snp.not.ok <- snpAnnot$snpID[snp.exclude]
anomStatsPlot(blData, genoData, anom.stats=stats[1,],
  snp.ineligible=snp.not.ok, centromere=centromeres)
@ 

\end{itemize}

\subsection{Identify low quality samples}

\begin{itemize}
  
\item{To identify low quality samples, one measure we use is the standard deviation
of BAF and LRR.  BAF results were found previously, now we find
results for LRR.  Unlike for BAF, all genotypes are included.}

<<>>=
lrr.sd <- sdByScanChromWindow(blData, var="LogRRatio", incl.hom=TRUE)
med.lrr.sd <- medianSdOverAutosomes(lrr.sd)
@ 

\item{We also need the number of segments found using 
circular binary segmentation in anomaly detection.}

<<>>=
baf.seg.info <- baf.anom$seg.info
loh.seg.info <- loh.anom$base.info[,c("scanID", "chromosome", "num.segs")]
@

\item{We identify low quality samples separately for BAF and LOH, using
different threshold parameters.  A SnpAnnotationDataFrame with an
``eligible'' column is required.
BAF detected anomalies for low quality BAF samples tend 
to have higher false positive rate. LOH detected anomalies
for low quality LOH samples tend to have higher false positive rate.}

<<>>=
snpAnnot$eligible <- !snp.exclude
baf.low.qual <- anomIdentifyLowQuality(snpAnnot, med.baf.sd, baf.seg.info,
  sd.thresh=0.1, sng.seg.thresh=0.0008, auto.seg.thresh=0.0001)
loh.low.qual <- anomIdentifyLowQuality(snpAnnot, med.lrr.sd, loh.seg.info,
  sd.thresh=0.25, sng.seg.thresh=0.0048, auto.seg.thresh=0.0006)
@ 

\item{Close the \texttt{IntensityData} and \texttt{GenotypeData} objects.}

<<>>=
close(blData)
close(genoData)
@ 

\end{itemize}

\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{SNP Quality Checks}

This step finds SNPs that may not be suitable for use in GWAS studies due to genotyping artifacts. 
Three methods are used to look at the genotyping error rates for each SNP: duplicate sample discordance, Mendelian 
error rates and deviation from Hardy-Weinberg equilibrium.

\subsection{Duplicate Sample Discordance}

This step calculates the discordance of genotype calls between samples that are duplicates.
Genotype discordance is evaluated by comparing the genotypes of samples that were genotyped more than once. 
We can examine the discordance rate with respect to samples or with regard to SNPs.
The discordance rate for a pair of samples is the fraction of genotype calls that differ over all
SNPs for which both calls are non-missing. The discordance rate for a SNP is the number of calls that differ divided by the number of duplicate pairs in which both calls are non-missing.

\begin{itemize}

\item{Keep the samples with a low enough value for the missing call rate, {\tt missing.e1}.  The threshold chosen here
is 0.05.}

<<>>=
library(GWASTools)
library(GWASdata)
data(affy_scan_annot)
scanAnnot <- ScanAnnotationDataFrame(affy_scan_annot)
scan.excl <- scanAnnot$scanID[scanAnnot$missing.e1 >= 0.05]
length(scan.excl)
@ 

\item{We make a vector of SNP {\tt snpIDs} with {\tt missing.n1} = 1 to exclude from the comparison.
We then call the {\tt duplicateDiscordance} function and save the output file.
This function finds subjectIDs for which there is more than one scanID.
To look at the discordance results, we will calculate the percentage value and look at the summary of the values for 
each of the duplicate pairs.
We will plot the rates color coded by continental ancestry, since experience has shown the values often differ based upon the population 
group.}

<<>>=
data(affy_snp_annot)
snpAnnot <- SnpAnnotationDataFrame(affy_snp_annot)
snp.excl <- snpAnnot$snpID[snpAnnot$missing.n1 == 1]
length(snp.excl)
genofile <- system.file("extdata", "affy_geno.nc",
  package="GWASdata")
genoNC <- NcdfGenotypeReader(genofile)
genoData <- GenotypeData(genoNC, snpAnnot=snpAnnot, scanAnnot=scanAnnot)
dupdisc <- duplicateDiscordance(genoData, subjName.col="subjectID",
  scan.exclude=scan.excl, snp.exclude=snp.excl)
names(dupdisc)
head(dupdisc$discordance.by.snp)
length(dupdisc$discordance.by.subject)
dupdisc$discordance.by.subject[[1]]
# each entry is a 2x2 matrix, but only one value of each
# is important since these are all pairs
npair <- length(dupdisc$discordance.by.subject)
disc.subj <- rep(NA, npair)
subjID <- rep(NA, npair)
race <- rep(NA, npair)
for (i in 1:npair) {
  disc.subj[i] <- dupdisc$discordance.by.subject[[i]][1,2]
  subjID[i] <- names(dupdisc$discordance.by.subject)[i]
  race[i] <- scanAnnot$race[scanAnnot$subjectID == subjID[i]][1]
}
dat <- data.frame(subjID=subjID, disc=disc.subj, pop=race,
                  stringsAsFactors=FALSE)
summary(dat$disc)
# Assign colors for the duplicate samples based on population group.
dat$col <- NA
dat$col[dat$pop == "CEU"] <- "black"
dat$col[dat$pop == "YRI"] <- "red"
table(dat$col, exclude=NULL)
dat <- dat[order(dat$disc),]
dat$rank <- 1:npair
@ 

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
# Plot the sample discordance rates color coded by ethnicity.
plot(dat$disc, dat$rank, 
  xlab="Discordance rate between duplicate samples", 
  ylab="rank", col=dat$col, 
  main="Duplicate Sample Discordance by Continental Ancestry")
legend("topleft", unique(dat$pop), 
  pch=rep(1,2), col=unique(dat$col))
@ 
\caption{The discordance for duplicate samples plotted in increasing values, color coded by population group.}
\label{fig:disc_ethn}
\end{figure}

\end{itemize}



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

\subsection{Mendelian Error Checking}

This step calculates and examines the Mendelian error rates. Mendelian errors are
detected in parent-offspring trios or pairs as offspring genotypes that are inconsistent with Mendelian inheritance.
We use the {\tt
  mendelErr} function to calculate a Mendelian error rate per SNP.
Lastly some checks are done on Mendelian error rates per family. 

\begin{itemize}

\item{To call the Mendelian error checking function, we first must
    create a {\tt mendelList} object.  We will call {\tt mendelList} that creates a list of trios, checking for any gender 
inconsistencies among annotated father and mother samples. Then, {\tt
  mendelListAsDataFrame} puts this list into a data frame for easier checking. Finally, we can call the {\tt mendelErr} function to find the 
Mendelian errors for SNPs with {\tt missing.n1} less than 0.05.}  

<<>>=
men.list <- mendelList(scanAnnot$family, scanAnnot$subjectID,
  scanAnnot$father, scanAnnot$mother, scanAnnot$sex,
  scanAnnot$scanID)
res <- mendelListAsDataFrame(men.list)
head(res)
dim(res)
# Only want to use SNPs with missing.n1 < 0.05
snp.excl <- snpAnnot$snpID[snpAnnot$missing.n1 >= 0.05]
length(snp.excl)
mend <- mendelErr(genoData, men.list, snp.exclude=snp.excl)
names(mend)
head(mend$trios)
names(mend$snp)
@ 

\end{itemize}

\subsection*{Mendelian Errors per SNP}

\begin{itemize}
  
\item{The Mendelian error rate is calculated for each SNP by dividing the number of errors per SNP for 
all trios by the
number of trios used in the error checking. We expect to get some {\tt
  NA} values in this calculation for those probes for which
there were no valid trios to check for errors, causing the rate to have a denominator of zero.}
  
<<>>=
# Calculate the error rate
err <- mend$snp$error.cnt / mend$snp$check.cnt
table(err == 0, exclude=NULL)
@ 

\item{We will plot the error rate per SNP as calculated above, sorting the values by increasing rate. }

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
plot(err, rank(err), xlab="Error Rate (fraction)",
  ylab="rank", main="Mendelian Error Rate per SNP, ranked")
@ 
\caption{A plot of the Mendelian error rate per SNP, ranked in increasing order.}
\end{figure}

\item{Next we will look at the Mendelian error rates among the trios we have in the HapMap data. Looking at the summary of the number of families with at least one 
error over all SNPs, we can see that the maximum number of errors per SNP. Next, we can look at subsets of SNPs with 
greater than 0, 1 and 2 errors per SNP. Finally, for those SNPs that have valid trios to detect errors, we get the 
fraction of SNPs with no errors.}

<<>>=
fam <- mend$snp$error.cnt
n <- mend$snp$check.cnt
fam[is.na(err)] <- NA
table(is.na(fam), exclude=NULL)
summary(fam)
# SNPs with errors
length(fam[n > 0 & fam > 0])
# SNPs for which more than one family has an error
length(fam[n > 0 & fam > 1])
length(fam[n > 0 & fam > 2])
# Get the SNPs with valid trios for error detection
val <- length(fam[n > 0])
noerr <- length(fam[n > 0 & fam == 0])
# Divide to get fraction with no errors
noerr / val
@ 

\item{A new version of the SNP annotation table will be saved with the Mendelian error values included. 
The number of families with at least one error per SNP, {\tt mend\$snp\$error.cnt}, gets saved as {\tt mendel.err.count}. 
The number of valid families for checking, {\tt mend\$snp\$check.cnt}, gets saved as {\tt mendel.err.sampsize}.}

<<>>=
# Write the error rates to the SNP table
snp.sel <- match(names(mend$snp$error.cnt), snpAnnot$snpID)
snpAnnot$mendel.err.count[snp.sel] <- mend$snp$error.cnt
snpAnnot$mendel.err.sampsize[snp.sel] <- mend$snp$check.cnt
allequal(snpAnnot$snpID, sort(snpAnnot$snpID))
# The high number of NA values is due to the filtering out of SNPs 
#    before the Mendelian error rate calculation
sum(is.na(snpAnnot$mendel.err.count))
sum(is.na(snpAnnot$mendel.err.sampsize))
# Save the updated SNP annotation table with 
#   the new Mendelian error data
varMetadata(snpAnnot)["mendel.err.count", "labelDescription"] <- 
  paste("number of Mendelian errors detected in trios averaged over",
        "multiple combinations of replicate genotyping instances")
varMetadata(snpAnnot)["mendel.err.sampsize", "labelDescription"] <- 
  "number of opportunities to detect Mendelian error in trios"
@ 

\item{To further investigate SNPs with a high Mendelian error rate, we will make a cluster plot for 3 SNPs
with the highest Mendelian error rate. We expect the 
plots to lack defined genotype clusters, leading to a poor call rate. For this example, we will plot the top 3 SNPs  
using the function {\tt genoClusterPlotByBatch}. }

<<results=hide>>=
# Get a vector of SNPs to check
snp <- pData(snpAnnot)
snp$err.rate <- snp$mendel.err.count / 
  snp$mendel.err.sampsize
snp <- snp[order(snp$err.rate, decreasing=TRUE),]
snp <- snp[1:3,]
xyfile <- system.file("extdata", "affy_qxy.nc",
  package="GWASdata")
xyNC <- NcdfIntensityReader(xyfile)
xyData <- IntensityData(xyNC, snpAnnot=snpAnnot, scanAnnot=scanAnnot)
pdf(file="DataCleaning-mendel.pdf")
par(mfrow = c(3,3))
mtxt <- paste("SNP", snp$rsID, "\nMendelian Error Rate",
   format(snp$err.rate, digits=5))
genoClusterPlotByBatch(xyData, genoData, snpID=snp$snpID, main.txt=mtxt,
                       plot.type="XY", batchVar="plate", cex.main=0.9) 
dev.off()
close(xyData)
@ 

\begin{figure}[ht!]
\centering
\includegraphics[width=5.67in,height=5.67in,keepaspectratio]{DataCleaning-mendel}
\caption{Cluster plots of the 3 SNPs with the highest Mendelian error rate.  Blue indicates a sample
with a ``BB'' genotype, green is an ``AB'' genotype and red is an ``AA'' genotype. The black X marks indicate
a sample with a missing genotype for that SNP.
All of these cluster plots are clearly problematic with regard to the genotype clusters and the called genotypes.  }
\end{figure}

\end{itemize}


\subsection*{Mendelian Errors per Family}

\begin{itemize}

\item{This section does some analyses on the Mendelian Errors for each family (trio). The variable {\tt all.trios} contains 
results of all combinations of duplicate samples. The variable {\tt trios} contains the averages of unique trios 
(averages of duplicates from {\tt all.trios}).}

<<>>=
# Calculate the fraction of SNPs with an error for each trio
trios <- mend$trios
trios$Mend.err <- trios$Men.err.cnt/trios$Men.cnt
summary(trios$Mend.err)
# Start by pulling out the vectors needed from `trios'
tmp <- trios[, c("fam.id", "Mend.err")]; dim(tmp)
# Change fam.id to match the sample annotation column name
names(tmp) <- c("family", "Mend.err.rate.fam")
# Merge the variables into the sample annotation file
scanAnnot$mend.err.rate.fam <- NA
for(i in 1:nrow(tmp))
{
  ind <- which(is.element(scanAnnot$family,tmp$family[i]))
  scanAnnot$mend.err.rate.fam[ind] <- tmp$Mend.err.rate.fam[i]
}
allequal(scanAnnot$scanID, sort(scanAnnot$scanID))
head(scanAnnot$mend.err.rate.fam)
varMetadata(scanAnnot)["mend.err.rate.fam", "labelDescription"] <-
  "Mendelian error rate per family"
@ 

\item{The Mendelian error rate per family, broken up by continental ancestry, could illuminate issues with SNPs that 
may not be accurately called across all ethnicities for the minor allele. We will plot the Mendelian error rate per family, color 
coded by population group. The error rates are higher for the YRI families as a whole, which is expected due to the 
higher level of genetic diversity.}

<<>>=
# Get the families that have non-NA values for the family 
#     Mendelian error rate
fams <- pData(scanAnnot)[!is.na(scanAnnot$mend.err.rate.fam) & 
  !duplicated(scanAnnot$family), c("family", 
  "mend.err.rate.fam", "race")]
dim(fams)
table(fams$race,exclude=NULL)
# Assign colors for the different ethnicities in these families
pcol <- rep(NA, nrow(fams))
pcol[fams$race == "CEU"] <- "black"
pcol[fams$race == "YRI"] <- "red"
@ 

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
plot(fams$mend.err.rate.fam*100, rank(fams$mend.err.rate.fam), 
  main="Mendelian Error rate per Family, ranked",
  xlab="Mendelian error rate per family (percent)", 
  ylab="rank", col=pcol)
legend("bottomright", c("CEU", "YRI"), 
  pch=c(1,1), col=c("black", "red"))
@ 
\caption{A plot of the Mendelian error rate for families, ranked in order of increasing rate.
The black indicates a family with CEU continental ancestry and red indicates a family with YRI ancestry.}
\end{figure}

\end{itemize}

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

\subsection{Hardy-Weinberg Equilibrium Testing}

This section uses Fisher's exact test to examine each SNP for departure from Hardy-Weinberg Equilibrium.  
For each SNP, p-values are obtained; those SNPs with extremely low values will be considered for filtering. 
QQ-plots of the p-values are made for both the autosomes and X chromosome. 

\begin{itemize}

\item{To run the Hardy-Weinberg test, we will filter out duplicates and non-founders. We will run {\tt gwasExactHW} for the samples with European continental ancestry only,
although the process is just the same for all population groups.  This
function filters out the males for the X chromosome and all samples for the
Y and mitochondrial probes.  Note that an optional argument to {\tt
  gwasExactHW} is a a scan by chromosome by filter matrix, to exclude
scan-chromosome combinations with aneuploidy.}

<<>>=
head(pData(scanAnnot)[,c("father", "mother")])
nonfounders <- scanAnnot$father != 0 & 
               scanAnnot$mother != 0
table(nonfounders)
scan.excl <- scanAnnot$scanID[scanAnnot$race != "CEU" |
   nonfounders | scanAnnot$duplicated]
length(scan.excl)
hwe <- gwasExactHW(genoData, scan.exclude=scan.excl)
close(genoData)
@ 

\item{ We will examine the CEU population data for the
purposes of the tutorial.  The remaining population groups and their results are analyzed in exactly the same manner.

We will look at the values calculated from the function call to {\tt gwasExactHW}, which include p-values, minor allele frequency, and genotype counts for each SNP on each of the chromosome 
types.  All the p-values are missing for the Y and M chromosomes. There are minor allele frequencies on the same chromosome types; we look at the first 50 
values. Finally, there are genotype counts for all chromosome types, which are stored as {\tt nAA, nAB} 
and {\tt nBB}.} 

<<>>=
names(hwe)
dim(hwe)
# Check on sample sizes for autosomes and X chromosome
hwe$N <- hwe$nAA + hwe$nAB + hwe$nBB
summary(hwe$N[is.element(hwe$chromosome,1:22)])
summary(hwe$N[is.element(hwe$chromosome,23)])
hwe$p.value[1:10]
sum(is.na(hwe$p.value[hwe$chromosome == 24])) # XY
sum(is.na(hwe$p.value[hwe$chromosome == 23])) # X
hwe$MAF[1:10]
hwe[1:3, c("nAA", "nAB", "nBB")]
@ 

\item{Next we want to estimate the inbreeding coefficient per SNP calculated using the minor allele frequencies 
and the total sample number count.  A histogram shows the distribution is centered around 0, which indicates 
there is most likely no significant population substructure.}

<<>>=
summary(hwe$f)
@ 

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
hist(hwe$f, main="Histogram of the Inbreeding Coefficient 
  For CEU Samples", xlab="Inbreeding Coefficient")
@ 
\caption{A histogram of the inbreeding coefficient values
from running the Hardy-Weinberg analysis for the unrelated CEU samples.}
\end{figure}

<<>>=
# Check the MAF of those SNPs with f=1
chkf <- hwe[!is.na(hwe$f) & hwe$f==1,]; dim(chkf)
summary(chkf$MAF)
@ 

\item{To see at what value the SNPs begin to deviate from the Hardy-Weinberg expected values, we will make QQ-plots 
that exclude SNPs where MAF = 0. Make a plot of the observed p-values vs. the expected p-values for the autosomes and 
X chromosome separately by calling the function {\tt qqPlot}. }

<<>>=
idx <- which(hwe$MAF == 0); length(idx)
hwe.0 <- hwe[-idx,]; dim(hwe.0)
# Check to makes sure that the correct number of SNPs were removed
length(hwe$p.value) - length(hwe.0$p.value)
# Only keep the autosomal SNPs for first plot
pval <- hwe.0$p.value[is.element(hwe.0$chromosome, 1:22)]
length(pval)
pval <- pval[!is.na(pval)]
length(pval)
# X chromosome SNPs for plot 2
pval.x <- hwe.0$p.value[is.element(hwe.0$chromosome, 23)]
length(pval.x)
pval.x <- pval.x[!is.na(pval.x)]
length(pval.x)
@ 
<<results=hide>>=
pdf(file = "DataCleaning-hwe.pdf")
par(mfrow=c(2,2))
qqPlot(pval=pval, truncate = FALSE, 
  main="Autosomes, all")
qqPlot(pval=pval, truncate = TRUE, 
  main="Autosomes, truncated")
qqPlot(pval=pval.x, truncate = FALSE, 
  main="X chromosome, all")
qqPlot(pval=pval.x, truncate = TRUE, 
  main="X chromosome, truncated")
dev.off()
@ 

\begin{figure}[ht!]
\centering
\includegraphics[width=5.67in,height=5.67in,keepaspectratio]{DataCleaning-hwe}
\caption{Four QQ-plots showing the deviation of Hardy-Weinberg p-values for each SNP for the samples from the 
CEU population group.} 
\end{figure}

\item{We plot the p-values against MAF for all SNPs with MAF greater than zero.}

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
plot(hwe.0$MAF, -log10(hwe.0$p.value), 
  xlab="Minor Allele Frequency", ylab="-log(p-value)", 
  main="Minor Allele Frequency vs\nP-value")
@ 
\caption{A plot of the p-value as a function of the minor allele frequency for SNPs for the CEU samples,
excluding the SNPs with a minor allele frequency equal to 0.}
\end{figure}

\end{itemize}


\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Preliminary Association Tests}

The final step in the data cleaning process is to perform preliminary association tests.
This step creates and examines QQ, `Manhattan' signal, regional association and genotype cluster plots.
If significant SNPs appear as a result of the association test, SNP cluster plots must be examined
to determine if the association is driven from a poorly clustering SNP.
Note that HapMap data do not come with phenotypic outcomes, thus, for purposes of the tutorial
we use simulated binary outcomes instead. The tests conducted are logistic regression based tests;
the samples are filtered by quality criteria and only unrelated subjects are included.
In the code below we do not include any covariates in the logistic regression as these data are not part of a case
control study. For data in the GENEVA project and other GWA studies we discuss which variables should be considered for inclusion as covariates in the
preliminary association tests.
The determination is made by analyzing a model including these covariates but without genotypes;
covariates with significant effects are then included in the final model.

\subsection{Association Test}

\begin{itemize}

\item{To run the association test, we call the function {\tt
      assocTestRegression}.  For the tutorial, all samples will be
    kept, but we would use the argument {\tt scan.exclude} for those samples we wish to filter out for the association test. Typically, we do not filter out SNPs for the association test -- we run all SNPs and then filter the results.  The omission of filters may 
cause some SNPs to return significant p-values for association.}

<<>>=
library(GWASTools)
library(GWASdata)
data(affy_scan_annot)
scanAnnot <- ScanAnnotationDataFrame(affy_scan_annot)
genofile <- system.file("extdata", "affy_geno.nc",
  package="GWASdata")
genoNC <- NcdfGenotypeReader(genofile)
genoData <- GenotypeData(genoNC, scanAnnot=scanAnnot)
assoc.file <- "assoc"
assocTestRegression(genoData, outcome="status", covar.list=c(""),
  ivar.list=c(""), model.type="logistic", robust=TRUE,
  chromosome.set=c(24:26), outfile=assoc.file,)
@ 

\end{itemize}

After running the association test on the selected subset of SNPs and samples we must analyze
the results to determine if any probes with significant p-values are spurious or truly associated with the phenotype
of interest.  Quantile-quantile, `Manhattan' and SNP cluster plots will all be used to further understand
those probes with significant p-values.

\subsection{QQ Plots}

\begin{itemize}

\item{To create QQ plots of the ordered p-values from the association tests versus the ordered expected p-values, we first load the saved p-values data.
Then, we will call {\tt qqPlot} to plot the Wald test p-values.
Given that the samples were split randomly between cases and controls, not surprisingly there are no outliers
visible in the QQ plot in Figure~\ref{fig:qq_wald}.}

<<>>=
assoc <- getobj(paste(assoc.file, 
  ".model.1.additive.chr.24_26.RData", sep=""))
names(assoc)
@ 

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
qqPlot(pval=assoc$model.1.additive.pvalue.G,
  truncate=TRUE, main="QQ Plot of Wald Test p-values")
@ 
\caption{QQ plot of Wald statistic p-values for association test.}
\label{fig:qq_wald}
\end{figure}

\end{itemize}
  

\subsection{``Manhattan'' Plots of the P-Values}

\begin{itemize}
  
\item{To create the `Manhattan' plots, we will call the function {\tt manhattanPlot}.  We take the negative log
transformation of the p-values and plot them for each probe.  See the resulting plot in Figure~\ref{fig:manh_wald}. }

<<>>=
data(affy_snp_annot)
snpAnnot <- SnpAnnotationDataFrame(affy_snp_annot)
chrom <- getChromosome(snpAnnot)
chrom.sel <- chrom %in% 24:26
@ 

\begin{figure}[ht!]
\centering
<<fig=TRUE>>=
manhattanPlot(assoc$model.1.additive.pvalue.G,
  chromosome=chrom[chrom.sel], 
  chrom.labels=c("XY","Y","M"))
@ 
\caption{``Manhattan'' plot of Wald statistic p-values for association test.}
\label{fig:manh_wald}
\end{figure}

\end{itemize}
  

\subsection{SNP Cluster Plots}

\begin{itemize}
  
\item{Next, we will create SNP cluster plots for the probes with significant p-values.
It is important to examine cluster plots of all top hits, as poor clusters not picked up by other quality
checking steps may still show up as having low p-values.
We plot SNPs with the 9 most significant p-values from
the Likelihood Ratio test. }

<<results=hide>>=
# Identify SNPs with lowest p-values
snp <- pData(snpAnnot)[chrom.sel, c("snpID", "rsID")]
allequal(snp$snpID, assoc$snpID)
snp$pval <- assoc$model.1.additive.pvalue.G
snp <- snp[order(snp$pval),]
snp <- snp[1:9,]
xyfile <- system.file("extdata", "affy_qxy.nc",
  package="GWASdata")
xyNC <- NcdfIntensityReader(xyfile)
xyData <- IntensityData(xyNC, snpAnnot=snpAnnot, scanAnnot=scanAnnot)
pdf(file="DataCleaning-cluster.pdf")
par(mfrow = c(3,3))
mtxt <- paste("SNP", snp$rsID, "\np =", format(snp$pval, digits=4))
genoClusterPlot(xyData, genoData, snpID=snp$snpID, main.txt=mtxt,
                plot.type="XY") 
dev.off()
close(xyData)
close(genoData)
@ 

\begin{figure}[ht!]
\centering
\includegraphics[width=5.67in,height=5.67in,keepaspectratio]{DataCleaning-cluster}
\caption{Cluster plots of the 9 SNPs with smallest p-values.}
\label{fig:snp_cluster}
\end{figure}

\end{itemize}

The cluster plots in Figure~\ref{fig:snp_cluster} show that some of the SNPs have poor cluster plots, indcating genotyping artifact.   In a real association study, the evidence from such SNPs would be discounted.

<<echo=FALSE, results=hide>>=
unlink(paste(assoc.file, "*.RData", sep=""))
@ 


\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Acknowledgements}
This manual reflects the work of many people. In the first place the
methods described were developed and implemented by a team headed by
Cathy Laurie. The team included David Crosslin, Stephanie Gogarten,
David Levine, Caitlin McHugh, Sarah Nelson, Jess Shen, Bruce Weir, Qi
Zhang and Xiuwen Zheng.  Before any the work started, valuable advice
was provided by Thomas Lumley and Ken Rice. Preparation of the manual
began with a team headed by Ian Painter and Stephanie Gogarten. The team included Marshall Brown, Matthew Conomos, Patrick Danaher, Kevin Rubenstein, Emily Weed and Leila Zelnick. 

The data cleaning activities of the GENEVA Coordinating Center have
been greatly helped by the experience and advice from other
participants in the GENEVA program: the genotyping centers at CIDR and the Broad; the dbGaP group at the National Center for Biotechnology Information (NCBI); and the many study investigators. Particular thanks to Kim Doheny and Elizabeth Pugh at CIDR and Stacey Gabriel and Daniel Mirel at the Broad and Justin Paschall at NCBI.

Funding for the GENEVA project includes HG 004446 (PI: Bruce Weir) for the Coordinating Center, U01 HG 004438 (PI: David Vallee) for CIDR, HG 004424 (PI: Stacey Gabriel) for the Broad.

 
The continuing guidance of Dr. Teri Manolio of NHGRI is deeply appreciated.


\end{document}
