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

# the example GDS file: 1,000 samples and 10,000 SNPs
(fn <- system.file("extdata", "grm1k_10k_snp.gds", package="SAIGEgds"))
gdsfile <- seqOpen(fn)

## -----------------------------------------------------------------------------
phenofn <- system.file("extdata", "pheno.txt.gz", package="SAIGEgds")
pheno <- read.table(phenofn, header=TRUE, as.is=TRUE)

set.seed(1000)
n <- nrow(pheno)
ftime <- rexp(n, rate=exp(0.3*pheno$x1 - 0.2*pheno$x2) * 0.05)  # time to event
ctime <- rexp(n, rate=0.03)                                     # censoring time

pheno$status <- as.integer(ftime <= ctime)  # 0/1 event status, the response
pheno$atime  <- pmin(ftime, ctime)          # observed follow-up time

head(pheno)
table(pheno$status)  # 1: event, 0: censored

## -----------------------------------------------------------------------------
glmm <- seqFitNullGLMM_SPA(status ~ x1 + x2, pheno, gdsfile,
    trait.type="survival", event.time="atime", sample.col="sample.id")

## -----------------------------------------------------------------------------
glmm$coefficients   # log hazard ratios for x1 and x2 (no intercept)
glmm$tau            # (Sigma_E, Sigma_G); Sigma_G is the frailty variance
glmm$converged

## -----------------------------------------------------------------------------
summary(glmm$fitted.values)   # expected numbers of events
sum(glmm$residuals)           # martingale residuals sum to ~0
head(glmm$var.ratio)

## -----------------------------------------------------------------------------
# the same model without a GRM (no random effects)
glmm0 <- seqFitNullGLMM_SPA(status ~ x1 + x2, pheno, trait.type="survival",
    event.time="atime", verbose=FALSE)

cf <- rbind(SAIGEgds.GRM=glmm$coefficients, SAIGEgds.noGRM=glmm0$coefficients)
if (requireNamespace("survival", quietly=TRUE))
{
    cx <- survival::coxph(survival::Surv(atime, status) ~ x1 + x2, pheno,
        ties="breslow")
    cf <- rbind(coxph=coef(cx), cf)
}
cf

## -----------------------------------------------------------------------------
assoc <- seqAssocGLMM_SPA(gdsfile, glmm, mac=10)
head(assoc)

## -----------------------------------------------------------------------------
table(assoc$method)

## ----top-variant-km, fig.width=6, fig.height=4, fig.align='center'------------
top <- assoc[which.min(replace(assoc$pval, !is.finite(assoc$pval), Inf)), ]
top[c("id", "chr", "pos", "ref", "alt", "AF.alt", "pval")]

seqFilterPush(gdsfile)
seqSetFilter(gdsfile, sample.id=glmm$sample.id, variant.id=top$id,
    verbose=FALSE)
sample.id <- seqGetData(gdsfile, "sample.id")
dosage <- drop(seqGetData(gdsfile, "$dosage_alt2"))
seqFilterPop(gdsfile)

ii <- match(sample.id, pheno$sample.id)
plot.data <- data.frame(
    time = pheno$atime[ii],
    status = pheno$status[ii],
    genotype = factor(dosage))
plot.data <- plot.data[is.finite(dosage), ]
levels(plot.data$genotype) <- paste0("ALT dosage = ", levels(plot.data$genotype))

km <- survival::survfit(survival::Surv(time, status) ~ genotype, data=plot.data)
cols <- seq_along(km$strata)
plot(km, col=cols, lwd=2, mark.time=TRUE,
    xlab="Follow-up time", ylab="Survival probability")
legend("bottomleft", sub("^genotype=", "", names(km$strata)),
    col=cols, lwd=2, bty="n")

## ----fig.width=4, fig.height=4, fig.align='center'----------------------------
p <- assoc$pval[is.finite(assoc$pval)]
median(qchisq(p, 1, lower.tail=FALSE)) / qchisq(0.5, 1)   # lambda_GC

# QQ plot
plot(-log10(ppoints(length(p))), -log10(sort(p)), pch=20, cex=0.6,
    xlab=expression(Expected~~-log[10](italic(p))),
    ylab=expression(Observed~~-log[10](italic(p))))
abline(0, 1, col="red")

## -----------------------------------------------------------------------------
units <- seqUnitSlidingWindows(gdsfile, win.size=500, win.shift=250)

burden <- seqAssocGLMM_Burden(gdsfile, glmm, units, verbose=FALSE)
head(burden)

acatv <- seqAssocGLMM_ACAT_V(gdsfile, glmm, units, verbose=FALSE)
head(acatv)

## -----------------------------------------------------------------------------
seqClose(gdsfile)

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