%\VignetteIndexEntry{ASSET Vignette}
%\VignettePackage{ASSET}
%\VigetteDepends{ASSET}

\documentclass[a4paper]{article}


\begin{document}
\SweaveOpts{concordance=TRUE}

\title{TREAT}
\maketitle

\section*{Introduction}
The TREAT package is designed to enhance statistical power for testing non-additive association 
between SNP sets and binary outcome in case-control studies. 
Various environmental and genetic covariates can be properly adjusted, 
as long as they are categorical or can be discretized into factors. 
The statistical test implemented in TREAT is complementary to those commonly used multi-locus tests
 that assume an additive risk model.
TREAT is a tree-based test. The state-of-the-art Boolean operation method is adopted to accelerate the
 tree-building procedure, which makes TREAT applicable to large-scale association studies, 
e.g., genome-wide association studies. 


<<start>>=
library(TREAT)
@

\section*{Testing association between SNP sets and a binary outcome}

For the treat function, the input data must be in a data frame.
The data frame should include a binary outcome (y), covariates to be adjusted (covar) and the SNPs to be tested (SNPs).
 Note that TREAT only accepts categorical covariates. If you have a continuous covariate, e.g., an eigenvector 
computed from the GWAS data,
 we recommend discretizing it into a factor variable X with a few levels. 
To do so, you may fit a simple logistic regression model on factor X and other covariates,
 then collapse those levels with similar effects on the outcome.
 The SNPs data must be coded as 0, 1, and 2.
<<data frame>>=
  set.seed(123)
  y     <- c(rep(0, 1000), rep(1, 1000))
  covar <- sample(0:2, 2000, replace=TRUE, prob=rep(1/3,3))
  SNPs  <- matrix(rbinom(2000*10, 2, c(.36, .48, .16)), nrow=2000)
  SNPs  <- as.data.frame(SNPs)
  colnames(SNPs) <- paste("rs", 1:10, sep="")
  data  <- data.frame(y, covar, SNPs)
@

To test an association, simply call the treat function as below.
You can specify the outcome and covariates in a formula. The data frame is passed into the function. 
The association between the outcome and those SNPs specified in snp.vars is then tested.
The function treat also accepts options (op) to control the tree-building procedure. 
Please refer to manual for more information.
<<TREAT>>=
tr <- treat(y ~ factor(covar), data = data, snp.vars = colnames(SNPs))
@

The function returns a p-value adjusted for multiple-comparisons. 
Some details about how the tree is built based on the observed data are also provided in \$tree.model.forward. 
For example, according to SNP rs2, node 1 (root node) is split as node 2 (if rs2 = 0 or 1) and node 3 (if rs2 = 2).
<<results1>>=
tr$adj.pval
tr$tree.model.forward
@


\section*{PLINK data fromat}

The TREAT package can load genotype data from plink binary files directly. 
This feature is particular useful when you are planning to perform a genome-wide gene-based analysis.  
The binary PED file (*.bed), the pedigree/phenotype information file (*.fam) and the extended MAP file (*.bim) 
are required in function treat.plink.
 Please refer to the website of PLINK for more information about how to create these files from the PED file and MAP file.
The TREAT package provides an example dataset in plink format, which were converted from the example provided 
by PLINK (http://pngu.mgh.harvard.edu/~purcell/plink/dist/example.zip). 

Get the data files from TREAT package
<<plink>>=
  bed <- system.file("sampleData", "data.bed", package="TREAT")
  bim <- system.file("sampleData", "data.bim", package="TREAT")
  fam <- system.file("sampleData", "data.fam", package="TREAT")
@

A data frame containing the binary outcome, covariates and id variables is required by the function treat.plink.
The id variables (e.g., "FAMILY", "SUBJECT") are used to describe pedigree information. Here we load them from the *fam file.
<<data>>=
  x <- read.table(fam, header=0, stringsAsFactors=FALSE)
  colnames(x) <- c("FAMILY", "SUBJECT", "FATHER", "MOTHER", "SEX", "Y")
@

The function treat.plink allows testing all given gene sets sequentially. 
Here we specify two genes by constructing a gene object gene.obj.
 It is a vector of gene names or a matrix with at least 3 columns called "Chr", "Start" and "Stop" 
which give the chromosome, starting 
location (physical position in base units) and ending location for each gene to analyze.
 This matrix can also have a column called "Gene" for the name of genes which is used in creating the names of the output files. 
Note that the format of gene.obj used here is very similar with
 that generate from PLINK (see section "Extract a subset of SNPs: file-list options" in  
http://pngu.mgh.harvard.edu/~purcell/plink/dataman.shtml)
<<gene object>>=
gene.obj <- rbind(c(8, 12799052, 12895289, "Gene1"), 
                    c(8, 12868315, 12989321, "Gene2"))
colnames(gene.obj) <- c("Chr", "Start", "Stop", "Gene")
@

Create the formula object, vector of id variables and call the treat.plink function. Note that the example data shown here has small sample size, we therefore modify the options to grow a larger tree, which will illustrate more details of the output. 
<<treat.plink>>=
formula <- as.formula(Y ~ SEX)
id.vars <- c("FAMILY", "SUBJECT")
op <- list(thr.sample = 10, thr.leaf.size = 5)
tr <- treat.plink(c(bed, bim, fam), formula, x, id.vars, gene.obj, op=op)
@

The function will return a list which contains a character vector of the output files created. 
Each output file is an R object file in directory op\$out.dir with name <gene>.rda containing the returned object from treat. 
The R object can be loaded back into R, from which a list with name obj can be obtained.
<<output>>=
tr$summary
load(tr$saved.files["Gene1"])
obj$adj.pval
obj$tree.model.forward
@

\section*{Session Information}
<<sessionInfo>>=
sessionInfo()
@ 

\end{document}


