### R code from vignette source 'SeqVarTools.Rnw'

###################################################
### code chunk number 1: style-Sweave
###################################################
BiocStyle::latex()


###################################################
### code chunk number 2: SeqVarTools.Rnw:37-43
###################################################
library(SeqVarTools)
vcffile <- seqExampleFileName("vcf")
gdsfile <- "tmp.gds"
seqVCF2GDS(vcffile, gdsfile, verbose=FALSE)
gds <- seqOpen(gdsfile)
gds


###################################################
### code chunk number 3: SeqVarTools.Rnw:48-52
###################################################
ref(gds)
head(refChar(gds))
alt(gds)
head(altChar(gds))


###################################################
### code chunk number 4: SeqVarTools.Rnw:56-57
###################################################
table(nAlleles(gds))


###################################################
### code chunk number 5: SeqVarTools.Rnw:63-67
###################################################
multi.allelic <- which(nAlleles(gds) > 2)
altChar(gds)[multi.allelic]
altChar(gds, n=1)[multi.allelic]
altChar(gds, n=2)[multi.allelic]


###################################################
### code chunk number 6: SeqVarTools.Rnw:71-73
###################################################
table(isSNV(gds))
isSNV(gds)[multi.allelic]


###################################################
### code chunk number 7: SeqVarTools.Rnw:78-81
###################################################
head(seqGetData(gds, "chromosome"))
head(seqGetData(gds, "position"))
granges(gds)


###################################################
### code chunk number 8: SeqVarTools.Rnw:85-87
###################################################
head(seqGetData(gds, "sample.id"))
head(seqGetData(gds, "variant.id"))


###################################################
### code chunk number 9: SeqVarTools.Rnw:95-98
###################################################
rsID <- seqGetData(gds, "annotation/id")
head(rsID)
length(unique(rsID)) == length(rsID)


###################################################
### code chunk number 10: SeqVarTools.Rnw:103-107
###################################################
seqClose(gds)
setVariantID(gdsfile, rsID)
gds <- seqOpen(gdsfile)
head(seqGetData(gds, "variant.id"))


###################################################
### code chunk number 11: SeqVarTools.Rnw:115-118
###################################################
geno <- getGenotype(gds)
dim(geno)
geno[1:10, 1:5]


###################################################
### code chunk number 12: SeqVarTools.Rnw:122-124
###################################################
geno <- getGenotypeAlleles(gds)
geno[1:10, 1:5]


###################################################
### code chunk number 13: SeqVarTools.Rnw:138-141
###################################################
samp.id <- seqGetData(gds, "sample.id")[1:10]
var.id <- seqGetData(gds, "variant.id")[1:5]
applyMethod(gds, getGenotype, variant=var.id, sample=samp.id)


###################################################
### code chunk number 14: SeqVarTools.Rnw:146-150
###################################################
library(GenomicRanges)
gr <- GRanges(seqnames="22", IRanges(start=1, end=250000000))
geno <- applyMethod(gds, getGenotype, variant=gr)
dim(geno)


###################################################
### code chunk number 15: SeqVarTools.Rnw:160-162
###################################################
titv(gds)
head(titv(gds, by.sample=TRUE))


###################################################
### code chunk number 16: SeqVarTools.Rnw:168-176
###################################################
binVar <- function(var, names, breaks) {
  names(var) <- names
  var <- sort(var)
  mids <- breaks[1:length(breaks)-1] + 
    (breaks[2:length(breaks)] - breaks[1:length(breaks)-1])/2
  bins <- cut(var, breaks, labels=mids, right=FALSE)
  split(names(var), bins)
}


###################################################
### code chunk number 17: SeqVarTools.Rnw:179-189
###################################################
variant.id <- seqGetData(gds, "variant.id")
afreq <- alleleFrequency(gds)
maf <- pmin(afreq, 1-afreq)
maf.bins <- binVar(maf, variant.id, seq(0,0.5,0.02))
nbins <- length(maf.bins)
titv.maf <- rep(NA, nbins)
for (i in 1:nbins) {
  capture.output(titv.maf[i] <- applyMethod(gds, titv, variant=maf.bins[[i]]))
}
plot(as.numeric(names(maf.bins)), titv.maf, xlab="MAF", ylab="TiTv")


###################################################
### code chunk number 18: SeqVarTools.Rnw:192-200
###################################################
miss <- missingGenotypeRate(gds)
miss.bins <- binVar(miss, variant.id, c(seq(0,0.5,0.05),1))
nbins <- length(miss.bins)
titv.miss <- rep(NA, nbins)
for (i in 1:nbins) {
  capture.output(titv.miss[i] <- applyMethod(gds, titv, variant=miss.bins[[i]]))
}
plot(as.numeric(names(miss.bins)), titv.miss, xlab="missing rate", ylab="TiTv")


###################################################
### code chunk number 19: SeqVarTools.Rnw:203-212
###################################################
depth <- seqApply(gds, "annotation/format/DP", mean, margin="by.variant", 
                  as.is="double", na.rm=TRUE)
depth.bins <- binVar(depth, variant.id, seq(0,200,20))
nbins <- length(depth.bins)
titv.depth <- rep(NA, nbins)
for (i in 1:nbins) {
  capture.output(titv.depth[i] <- applyMethod(gds, titv, variant=depth.bins[[i]]))
}
plot(as.numeric(names(depth.bins)), titv.depth, xlab="mean depth", ylab="TiTv")


###################################################
### code chunk number 20: SeqVarTools.Rnw:219-222
###################################################
miss.var <- missingGenotypeRate(gds, margin="by.variant")
het.var <- heterozygosity(gds, margin="by.variant")
filt <- seqGetData(gds, "variant.id")[miss.var <= 0.1 & het.var <= 0.6]


###################################################
### code chunk number 21: SeqVarTools.Rnw:227-231
###################################################
seqSetFilter(gds, variant.id=filt)
hethom <- hethom(gds)
hist(hethom, main="", xlab="Het/Hom Non-Ref")
seqSetFilter(gds)


###################################################
### code chunk number 22: SeqVarTools.Rnw:239-242
###################################################
pc <- pca(gds)
names(pc)
plot(pc$eigenvect[,1], pc$eigenvect[,2])


###################################################
### code chunk number 23: SeqVarTools.Rnw:253-259
###################################################
hw <- hwe(gds)
pval <- -log10(sort(hw$p))
hw.perm <- hwe(gds, permute=TRUE)
x <- -log10(sort(hw.perm$p))
plot(x, pval, xlab="-log10(expected P)", ylab="-log10(observed P)")	
abline(0,1,col="red")


###################################################
### code chunk number 24: SeqVarTools.Rnw:264-265
###################################################
hist(hw$f)


###################################################
### code chunk number 25: SeqVarTools.Rnw:269-271
###################################################
ic <- inbreedCoeff(gds, margin="by.sample")
range(ic)


###################################################
### code chunk number 26: SeqVarTools.Rnw:279-284
###################################################
data(pedigree)
pedigree[pedigree$family == 1463,]
err <- mendelErr(gds, pedigree, verbose=FALSE)
table(err$by.variant)
err$by.trio


###################################################
### code chunk number 27: SeqVarTools.Rnw:296-315
###################################################
library(Biobase)
sample.id <- seqGetData(gds, "sample.id")
pedigree <- pedigree[match(sample.id, pedigree$sample.id),]
n <- length(sample.id)
pedigree$phenotype <- rnorm(n, mean=10)
pedigree$case.status <- rbinom(n, 1, 0.3)
sample.data <- AnnotatedDataFrame(pedigree)

seqData <- SeqVarData(gds, sample.data)

## continuous phenotype
assoc <- regression(seqData, outcome="phenotype", covar="sex",
                    model.type="linear")
head(assoc)
pval <- -log10(sort(assoc$Wald.Pval))
n <- length(pval)
x <- -log10((1:n)/n)
plot(x, pval, xlab="-log10(expected P)", ylab="-log10(observed P)")
abline(0, 1, col = "red")


###################################################
### code chunk number 28: SeqVarTools.Rnw:325-333
###################################################
assoc <- regression(seqData, outcome="case.status", covar="sex",
                    model.type="firth")
head(assoc)
pval <- -log10(sort(assoc$PPL.Pval))
n <- length(pval)
x <- -log10((1:n)/n)
plot(x, pval, xlab="-log10(expected P)", ylab="-log10(observed P)")
abline(0, 1, col = "red")


###################################################
### code chunk number 29: SeqVarTools.Rnw:338-339
###################################################
seqClose(gds)


###################################################
### code chunk number 30: SeqVarTools.Rnw:342-343
###################################################
unlink(gdsfile)