%\VignetteIndexEntry{SKAT}
\documentclass{article}

\usepackage{amsmath}
\usepackage{amscd}
\usepackage[tableposition=top]{caption}
\usepackage{ifthen}
\usepackage[utf8]{inputenc}

\begin{document}

\title{SKAT Package}
\author{Seunggeun (Shawn) Lee}
\maketitle

\section{Overview}
SKAT package has functions to 1) test an association between a SNP set
and continuous/binary phenotypes and 2) to compute power/sample size for future sequence association studies. 

\section{Association test}

\subsection{Example}
An example dataset (SKAT.example) has a genotype matrix (Z) of  2000 individuals and 67 SNPs, 
vectors of continuous (y.c) and binary (y.b) phenotypes, and a covariates matrix (X).

<<data>>=
library(SKAT)
data(SKAT.example)
names(SKAT.example)

attach(SKAT.example)
@

To test an association, SKAT\_Null\_Model function should be first used to estimate parameters
and to obtain residuals under the null model of no associations. And then, 
SKAT function can be used to get p-values.  

<<SKAT1>>=
# continuous trait 
obj<-SKAT_Null_Model(y.c ~ X, out_type="C")
SKAT(Z, obj)$p.value

# dichotomous trait 
obj<-SKAT_Null_Model(y.b ~ X, out_type="D")
SKAT(Z, obj)$p.value

@

When the trait is binary and the sample size is small, 
SKAT can produce conservative results. We developed a moment matching adjustment method 
that adjusts the asymptotic null distribution by estimating empirical variance and kurtosis.
By default, SKAT ( >= ver 0.7) will conduct a small sample adjustment 
when the sample size $ < 2000$. In the following code, we use only 200 samples to run SKAT.

<<SKAT11>>=

IDX<-c(1:100,1001:1100)	
# With-adjustment
obj.s<-SKAT_Null_Model(y.b[IDX] ~ X[IDX,],out_type="D")
SKAT(Z[IDX,], obj.s, kernel = "linear.weighted")$p.value

@

If you don't want to use the adjustment, please set Adjustment=FALSE in the SKAT\_Null\_Model function.

<<SKAT12>>=
# Without-adjustment
obj.s<-SKAT_Null_Model(y.b[IDX] ~ X[IDX,],out_type="D", Adjustment=FALSE)
SKAT(Z[IDX,], obj.s, kernel = "linear.weighted")$p.value
@


We recently developed an efficient resampling method to compute p-values for binary traits, 
and the methods are implemented in SKATBinary function. When you use this function, Adjustment=TRUE in SKAT\_Null\_Model is not necessary.
Implemented methods are 1) Efficient resampling (ER); 2) ER with adaptive resampling (ER.A);
3) Quantile adjusted moment matching (QA); 4) Moment matching adjustment (MA);
5) No adjustment (UA);  and 6) Hybrid. 
"Hybrid" (default method) selects a method based on the total minor allele count (MAC), the number of individuals with minor 
alleles (m), and the degree of case-control imbalance. 
  
<<SKAT22>>=

# default hybrid approach 
out<-SKATBinary(Z[IDX,], obj.s, kernel = "linear.weighted")
out$p.value

@


\subsection{Assign weights for each SNP}

It is assumed that rarer variants are more likely to be causal variants with large effect sizes. 
To incorporate this assumption, the linear weighted kernel uses a weighting scheme and is formulated as 
$ Z W W Z'$, where $Z$ is a genotype matrix, and $W = diag \{ w_1, \ldots, w_m \}$ is a weight matrix. 
In the previous examples, we used the default beta(1,25) weight, 
$w_i = dbeta(p_i, 1, 25) $, where $dbeta$ is a beta density function,
and $p_i$ is a minor allele frequency (MAF) of SNP $i$.
You can use different parameters for the beta weight by changing the weights.beta parameter. 
For example, if you want to use Madsen and Browning weight, use weight.beta=c(0.5,0.5). 


<<SKAT3>>=
SKAT(Z, obj, kernel = "linear.weighted", weights.beta=c(0.5,0.5))$p.value
@

You can make your own weight vector and use it for the weighting. 
For the logistic weight, we provide a function to generate the weight. 

<<SKAT4>>=
# Shape of the logistic weight

MAF<-1:1000/1000
W<-Get_Logistic_Weights_MAF(MAF, par1=0.07, par2=150)
par(mfrow=c(1,2))
plot(MAF,W,xlab="MAF",ylab="Weights",type="l")
plot(MAF[1:100],W[1:100],xlab="MAF",ylab="Weights",type="l")
par(mfrow=c(1,2))

# Use logistic weight
weights<-Get_Logistic_Weights(Z, par1=0.07, par2=150)
SKAT(Z, obj, kernel = "linear.weighted", weights=weights)$p.value
@



\subsection{Combined Test of burden test and SKAT}

A test statistic of the combined test is
$$Q_{\rho} = (1-\rho) Q_S + \rho Q_B,$$
where $Q_S$ is a test statistic of SKAT, and $Q_B$ is a score test statistic of the burden test. 
You can specify $\rho$ value using the r.corr parameter (default: , r.corr=0).

<<SKAT41>>=
#rho=0
SKAT(Z, obj, r.corr=0)$p.value

#rho=0.9
SKAT(Z, obj, r.corr=0.9)$p.value

#rho=1, burden test
SKAT(Z, obj, r.corr=1)$p.value
@


If method=``optimal.adj'', SKAT-O method is performed, 
which computes p-values with eight different values of $\rho=(0, 0.1^2, 0.2^2, 0.3^2, 0.4^2, 0.5^2, 0.5, 1)$
and then uses the minimum p-value as a test statistic. 
If you want to use the original implementation of SKAT-O, use method=``optimal''. 
We recommend to use ``optimal.adj'', since it has a better type I error control.

<<SKAT42>>=

#Optimal Test
SKAT(Z, obj, method="optimal.adj")$p.value

@

\subsection{Combined test of rare and common variants}

If you want to test combined effects of common and rare variants,  you can use SKAT\_CommonRare function. 

<<SKAT43>>=
# Combined sum test (SKAT-C and Burden-C)

SKAT_CommonRare(Z, obj)$p.value
SKAT_CommonRare(Z, obj, r.corr.rare=1, r.corr.common=1 )$p.value

# Adaptive test (SKAT-A and Burden-A)

SKAT_CommonRare(Z, obj, method="A")$p.value
SKAT_CommonRare(Z, obj, r.corr.rare=1, r.corr.common=1, method="A" )$p.value

@

The detailed description of each method can be found in the following reference. \\

Ionita-Laza, I., Lee, S., Makarov, V., Buxbaum, J. Lin, X. (2013). 
Sequence kernel association tests for the combined effect of rare and common variants.  
\emph{American Journal of Human Genetics}, 92(6):841-53.  \\

\subsection{Imputing missing genotypes.}

If there are missing genotypes, SKAT automatically imputes them
based on Hardy-Weinberg equilibrium. 
You can choose from ``bestguess'', ``fixed'' or ``random''. 
The ``bestguess'' imputes missing genotypes as most likely values (0,1,2), 
the ``fixed'' imputes missing genotypes by assigning the mean genotype value (2p, p is the MAF)
and the "random" imputes missing genotypes by generating binomial(2,p) random variables. The default imputation 
method for the SKAT function is ``fixed'' and for the SKATBinary function is ``bestguess''.

<<SKAT5>>=
# Assign missing 
Z1<-Z
Z1[1,1:3]<-NA

# bestguess imputation
SKAT(Z1,obj,impute.method = "bestguess")$p.value

# fixed imputation
SKAT(Z1,obj,impute.method = "fixed")$p.value

# random imputation
SKAT(Z1,obj,impute.method = "random")$p.value


@

\subsection{Resampling}

SKAT package provides functions to carry out resampling to
compute empirical p-values and to control family wise error rate. 
Two different resampling methods are implemented.
``bootstrap'' conducts the parametric bootstrap to resample residuals from $H_0$ 
with considering covariates. When there is no covariate, ``bootstrap'' is equivalent to
the permutation. ``perturbation'' perturbs the residuals by multiplying 
standard normal random variables. The default method is ``bootstrap''.
From ver 0.7, we do not provide the ``perturbation'' method.


<<SKAT6>>=
# parametric boostrap.
obj<-SKAT_Null_Model(y.b ~ X, out_type="D", n.Resampling=5000, 
type.Resampling="bootstrap")

# SKAT p-value
re<- SKAT(Z, obj, kernel = "linear.weighted")
re$p.value	# SKAT p-value
Get_Resampling_Pvalue(re)	# get resampling p-value
@

When there are many genes/SNP sets to test, 
resampling methods can be used to control family-wise error rate. 
You can find an example in the next section. 

\subsection{Plink Binary format files}

SKAT package can use plink binary format files for genome-wide data analysis. 
To use plink files, plink bed, bim and fam files, and your own setid file 
that contains information of SNP sets are needed. 
Example files can be found on the SKAT/MetaSKAT google group page.

<<SKAT_B1>>=
# To run this code, first download and unzip example files

##############################################
# 	Generate SSD file

# Create the MW File
File.Bed<-"./Example1.bed"
File.Bim<-"./Example1.bim"
File.Fam<-"./Example1.fam"
File.SetID<-"./Example1.SetID"
File.SSD<-"./Example1.SSD"
File.Info<-"./Example1.SSD.info"

# To use binary ped files, you have to generate SSD file first.
# If you already have a SSD file, you do not need to call this function. 
Generate_SSD_SetID(File.Bed, File.Bim, File.Fam, File.SetID, File.SSD, File.Info)
@

Now you can open SSD and Info file and run SKAT. 

<<SKAT_B2>>=
FAM<-Read_Plink_FAM(File.Fam, Is.binary=FALSE)
y<-FAM$Phenotype

# To use a SSD file, please open it first. After finishing using it, you must close it.
 
SSD.INFO<-Open_SSD(File.SSD, File.Info)

# Number of samples 
SSD.INFO$nSample 

# Number of Sets
SSD.INFO$nSets

obj<-SKAT_Null_Model(y ~ 1, out_type="C")
out<-SKAT.SSD.All(SSD.INFO, obj)
out
@

If you have a plink covariate file, you can use Read\_Plink\_FAM\_Cov file to read both FAM and covariate files. 

<<SKAT_B2Cov>>=
File.Cov<-"./Example1.Cov"
FAM_Cov<-Read_Plink_FAM_Cov(File.Fam, File.Cov, Is.binary=FALSE)

# First 5 rows
FAM_Cov[1:5,]

# Run with covariates
X1 = FAM_Cov$X1
X2 = FAM_Cov$X2
y<-FAM_Cov$Phenotype

obj<-SKAT_Null_Model(y ~ X1 + X2, out_type="C")
out<-SKAT.SSD.All(SSD.INFO, obj)
out
@

To use custom weight, you need to make a weight file and read it using
``Read\_SNP\_WeightFile'' function. The weight file should have two columns, 
SNP ID and weight values. The output object of ``Read\_SNP\_WeightFile''
can be used as a parameter in SKAT.SSD functions

<<SKAT_B2Weight>>=

# Custom weight
# File: Example1_Weight.txt
obj.SNPWeight<-Read_SNP_WeightFile("./Example1_Weight.txt")

out<-SKAT.SSD.All(SSD.INFO, obj, obj.SNPWeight=obj.SNPWeight)
out
@

The output object of SKAT.SSD.All has an output dataframe object ``results''. 
You can save it using write.table function.
<<SKAT_B2Save>>=

output.df = out$results
write.table(output.df, file="./save.txt", col.names=TRUE, row.names=FALSE)

@


If more than one gene/SNP sets are to be tested, 
you should adjust for multiple testing to control for family-wise error rate. 
It can be done bonferroni correction. If gene/SNP sets are correlated, however,
this approach can be conservative. 
Alternatively, you can directly control family wise error rate (FWER) using the resampling method. 

<<SKAT_B3>>==
obj<-SKAT_Null_Model(y ~ 1, out_type="C", n.Resampling=1000, type.Resampling="bootstrap")
out<-SKAT.SSD.All(SSD.INFO, obj)

# No gene is significant with controling FWER = 0.05
Resampling_FWER(out,FWER=0.05)

# 1 gene is significnat with controling FWER = 0.5
Resampling_FWER(out,FWER=0.5)
@

``SKAT.SSD.OneSet'' or ``SKAT.SSD.OneSet\_SetIndex'' functions can be used
to test a single gene/SNP set. Alternatively, 
you can obtain a genotype matrix using ``Get\_Genotypes\_SSD'' function and then run SKAT. 

<<SKAT_B4>>==

obj<-SKAT_Null_Model(y ~ 1, out_type="C")

# test the second gene
id<-2
SetID<-SSD.INFO$SetInfo$SetID[id]
SKAT.SSD.OneSet(SSD.INFO,SetID, obj)$p.value
 
SKAT.SSD.OneSet_SetIndex(SSD.INFO,id, obj)$p.value

# test the second gene with the logistic weight.
Z<-Get_Genotypes_SSD(SSD.INFO, id)
weights = Get_Logistic_Weights(Z, par1=0.07, par2=150)
SKAT(Z, obj, weights=weights)$p.value

@

SKAT\_CommonRare function also can be used with SSD files. 

<<SKAT_B5>>==

# test all genes in SSD file
obj<-SKAT_Null_Model(y ~ X1 + X2, out_type="C")
out<-SKAT_CommonRare.SSD.All(SSD.INFO, obj)
out


@


After finishing, please close the SSD file. 

<<SKAT_B5>>==
Close_SSD()
@


\subsection{Plink Binary format files: SKATBinary}

SKATBinary functions can be used with plink formatted files. This section shows 
example code. Example plink files can be found on the SKAT/MetaSKAT google group page.

<<SKAT_BB1>>=

# File names
File.Bed<-"./SKATBinary.example.bed"
File.Bim<-"./SKATBinary.example.bim"
File.Fam<-"./SKATBinary.example.fam"
File.Cov<-"./SKATBinary.example.cov"
File.SetID<-"./SKATBinary.example.SetID"
File.SSD<-"./SKATBinary.example.SSD"
File.Info<-"./SKATBinary.example.SSD.info"

# Generate SSD file, and read fam and cov files
# If you already have a SSD file, you do not need to call this function. 
Generate_SSD_SetID(File.Bed, File.Bim, File.Fam, File.SetID, File.SSD, File.Info)
FAM<-Read_Plink_FAM_Cov(File.Fam, File.Cov, Is.binary=TRUE, cov_header=FALSE)

# open SSD files
 
SSD.INFO<-Open_SSD(File.SSD, File.Info)

# No adjustment is needed
obj<-SKAT_Null_Model(Phenotype ~ COV1 + COV2, out_type="D", data=FAM, Adjustment=FALSE)


# SKAT
out.skat<-SKATBinary.SSD.All(SSD.INFO, obj, method="SKAT")

# SKAT-O
out.skato<-SKATBinary.SSD.All(SSD.INFO, obj, method="SKATO")

# First 5 variant sets, SKAT
out.skat$results[1:5,]

@

The effective number of tests and QQ plots can be obtained using the minimum achievable p-values (MAP).


<<SKAT_BB2>>=

# Effective number of test is smaller than 30 (number of variant sets)
# Use SKAT results
Get_EffectiveNumberTest(out.skat$results$MAP, alpha=0.05)

# QQ plot
QQPlot_Adj(out.skat$results$P.value, out.skat$results$MAP)

@



\section{Power/Sample Size calculation.}


\subsection{Dataset}
SKAT package provides a haplotype dataset (SKAT.haplotypes) 
which contains a haplotype matrix of 10,000 haplotypes over 200kb region (Haplotype), 
and a dataframe with information on each SNP. 
These haplotypes were simulated using a calibrated coalescent model (cosi) with mimicking 
linkage disequilibrium structure of European ancestry. 
If you don't have any haplotype information, use this dataset to compute power/sample size.

<<data>>=
data(SKAT.haplotypes)
names(SKAT.haplotypes)

attach(SKAT.haplotypes)
@

\subsection{Power/Sample Size calculation}

The following example uses the haplotypes in SKAT.haplotypes with 
the following parameters. 

\begin{enumerate}
\item Subregion length = 3k bp 
\item Causal percent = $20 \%$
\item Negative percent = $20 \%$
\item For continuous traits, $\beta = c |log_{10}(MAF)|$ (BetaType = ``Log'') with $\beta = 2$ at MAF = $10^{-4}$
\item For binary traits, $log(OR) = c |log_{10}(MAF)|$ (OR.Type = ``Log'') with OR $= 2$ at MAF = $10^{-4}$, and $50 \%$ of samples are cases and $50 \% $ of samples are controls
\end{enumerate}

<<SKAT_P1>>==
set.seed(500)
out.c<-Power_Continuous(Haplotype,SNPInfo$CHROM_POS, SubRegion.Length=5000,    
Causal.Percent= 20, N.Sim=10, MaxBeta=2,Negative.Percent=20)
out.b<-Power_Logistic(Haplotype,SNPInfo$CHROM_POS, SubRegion.Length=5000,   
Causal.Percent= 20, N.Sim=10 ,MaxOR=7, Negative.Percent=20)

out.c
out.b

Get_RequiredSampleSize(out.c, Power=0.8)
Get_RequiredSampleSize(out.b, Power=0.8)

@

In this example, N.Sim=10 was used to get results quickly. 
When you do the power calculation, please increase it to more than 100. 
When BetaType = ``Log'' or OR.Type = ``Log'', the effect size of  continuous trait 
and the log odds ratio of binary traits are $c |log_{10}(MAF)|$, 
where $c$ is determined by Max\_Beta or Max\_OR. 
For example, $ c= 2/4 = 0.5$ when the Max\_Beta = 2. 
In this case, a causal variant with MAF=0.01 has $\beta = 1$. 
For binary traits, $c= log(7)/4 = 0.486$ with MAX\_OR=7. 
And thus, a causal variant with MAF=0.01 has log OR = 0.972.

If you consider non-zero  r.corr ($\rho$) values to compute power, 
Power\_Continuous\_R or Power\_Logistic\_R functions can be used instead.
Since they use slightly different method to compute power, 
power estimates from Power\_Continuous\_R and Power\_Logistic\_R can be slightly different from 
estimates from Power\_Continuous and Power\_Logistic although r.corr=0.

If you want to computer the power of SKAT-O by estimating the optimal r.corr, use r.corr=2. 
The estimated optimal r.corr is 
$$
r.corr = p_1^2 ( 2p_2-1)^2,
$$
where $p_1$ is the proportion of nonzero $\beta$s, and $p_2$ is the proportion of negative (or positive) $\beta$s 
among the non-zero $\beta$s. 


<<SKAT_P2>>==
set.seed(500)
out.c<-Power_Continuous_R(Haplotype,SNPInfo$CHROM_POS, SubRegion.Length=5000,    
Causal.Percent= 20, N.Sim=10, MaxBeta=2,Negative.Percent=20, r.corr=2)

out.c

Get_RequiredSampleSize(out.c, Power=0.8)

@


\end{document}

