## ----message=FALSE------------------------------------------------------------
library(Matrix)
library(SeqArray)
library(SAIGEgds)

# Common markers used by the GRMs: 1,000 samples and 10,000 SNPs
grm.fn <- system.file("extdata", "grm1k_10k_snp.gds", package="SAIGEgds")

# For this compact example, variants to group and test use the same file.
# Production analyses should use a separate exome or whole-genome GDS file.
assoc.fn <- grm.fn

# Binary phenotype and covariates
pheno.fn <- system.file("extdata", "pheno.txt.gz", package="SAIGEgds")
pheno <- read.table(pheno.fn, header=TRUE, as.is=TRUE)
head(pheno)
table(pheno$y)

## ----eval=FALSE---------------------------------------------------------------
# grm.gds <- seqOpen(grm.fn)
# sp.grm <- seqFitSparseGRM(grm.gds, sample.id=pheno$sample.id,
#     nsnp.sub.random=2000, rel.cutoff=0.125, num.thread=2)
# seqClose(grm.gds)
# 
# # Save once and reuse for phenotypes measured in the same cohort
# saveRDS(sp.grm, "sparse_grm.rds")

## -----------------------------------------------------------------------------
sp.grm.fn <- system.file("extdata", "grm1k_10k_sp_grm.rds",
    package="SAIGEgds")
sp.grm <- readRDS(sp.grm.fn)

dim(sp.grm)
class(sp.grm)
nnzero(sp.grm)
nnzero(sp.grm) / prod(dim(sp.grm))
stopifnot(setequal(colnames(sp.grm), pheno$sample.id))

## -----------------------------------------------------------------------------
glmm <- seqFitNullGLMM_SPA(y ~ x1 + x2, pheno,
    gdsfile=grm.fn, grm.mat=sp.grm,
    trait.type="binary", sample.col="sample.id",
    use.cateMAC=TRUE, num.thread=2, verbose=FALSE)

## -----------------------------------------------------------------------------
glmm$converged
glmm$tau
head(glmm$var.ratio)
c(Sigma_inv = !is.null(glmm$Sigma_inv),
    chol_inv_X_Sigma = !is.null(glmm$chol_inv_X_Sigma))

## ----eval=FALSE---------------------------------------------------------------
# glmm <- seqFitNullGLMM_SPA(y ~ x1 + x2, pheno,
#     gdsfile=grm.fn, grm.mat=TRUE,
#     trait.type="binary", sample.col="sample.id",
#     use.cateMAC=TRUE, nsnp.sub.random=2000, rel.cutoff=0.125)

## -----------------------------------------------------------------------------
assoc.gds <- seqOpen(assoc.fn)
units <- seqUnitSlidingWindows(assoc.gds, win.size=5000, win.shift=5000)
units

## -----------------------------------------------------------------------------
burden <- seqAssocGLMM_Burden(assoc.gds, glmm, units,
    maxMAF=c(0.01, 0.005), parallel=1, verbose=FALSE)
burden

## -----------------------------------------------------------------------------
skat <- seqAssocGLMM_SKAT(assoc.gds, glmm, units,
    maxMAF=c(0.01, 0.005), collapse.mac=10, parallel=1, verbose=FALSE)
skat

## -----------------------------------------------------------------------------
acatv <- seqAssocGLMM_ACAT_V(assoc.gds, glmm, units,
    maxMAF=c(0.01, 0.005), collapse.mac=10, parallel=1, verbose=FALSE)
acatv

## -----------------------------------------------------------------------------
acato <- seqAssocGLMM_ACAT_O(assoc.gds, glmm, units,
    maxMAF=c(0.01, 0.005), collapse.mac=10, parallel=1, verbose=FALSE)
acato

## ----echo=FALSE---------------------------------------------------------------
seqClose(assoc.gds)

## -----------------------------------------------------------------------------
sessionInfo()