### R code from vignette source 'vignettes/SeqVarTools/inst/doc/SeqVarTools.Rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: SeqVarTools.Rnw:44-50 ################################################### library(SeqVarTools) vcffile <- seqExampleFileName("vcf") gdsfile <- "tmp.gds" seqVCF2GDS(vcffile, gdsfile, verbose=FALSE) gds <- seqOpen(gdsfile) gds ################################################### ### code chunk number 2: SeqVarTools.Rnw:55-59 ################################################### ref(gds) head(refChar(gds)) alt(gds) head(altChar(gds)) ################################################### ### code chunk number 3: SeqVarTools.Rnw:63-64 ################################################### table(nAlleles(gds)) ################################################### ### code chunk number 4: SeqVarTools.Rnw:70-74 ################################################### 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 5: SeqVarTools.Rnw:78-80 ################################################### table(isSNV(gds)) isSNV(gds)[multi.allelic] ################################################### ### code chunk number 6: SeqVarTools.Rnw:85-88 ################################################### head(seqGetData(gds, "chromosome")) head(seqGetData(gds, "position")) granges(gds) ################################################### ### code chunk number 7: SeqVarTools.Rnw:92-94 ################################################### head(seqGetData(gds, "sample.id")) head(seqGetData(gds, "variant.id")) ################################################### ### code chunk number 8: SeqVarTools.Rnw:102-105 ################################################### rsID <- seqGetData(gds, "annotation/id") head(rsID) length(unique(rsID)) == length(rsID) ################################################### ### code chunk number 9: SeqVarTools.Rnw:110-114 ################################################### seqClose(gds) setVariantID(gdsfile, rsID) gds <- seqOpen(gdsfile) head(seqGetData(gds, "variant.id")) ################################################### ### code chunk number 10: SeqVarTools.Rnw:122-125 ################################################### geno <- getGenotype(gds) dim(geno) geno[1:10, 1:5] ################################################### ### code chunk number 11: SeqVarTools.Rnw:129-131 ################################################### geno <- getGenotypeAlleles(gds) geno[1:10, 1:5] ################################################### ### code chunk number 12: SeqVarTools.Rnw:145-148 ################################################### 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 13: SeqVarTools.Rnw:153-157 ################################################### library(GenomicRanges) gr <- GRanges(seqnames="22", IRanges(start=1, end=250000000)) geno <- applyMethod(gds, getGenotype, variant=gr) dim(geno) ################################################### ### code chunk number 14: SeqVarTools.Rnw:167-169 ################################################### titv(gds) head(titv(gds, by.sample=TRUE)) ################################################### ### code chunk number 15: SeqVarTools.Rnw:175-183 ################################################### 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 16: SeqVarTools.Rnw:186-196 ################################################### 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 17: SeqVarTools.Rnw:199-207 ################################################### 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 18: SeqVarTools.Rnw:210-219 ################################################### 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 19: SeqVarTools.Rnw:226-229 ################################################### 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 20: SeqVarTools.Rnw:236-242 ################################################### seqSetFilter(gds, variant.id=filt) het <- heterozygosity(gds, margin="by.sample") homnr <- homozygosity(gds, margin="by.sample", allele="alt") hethom <- het / homnr hist(hethom, main="", xlab="Het/Hom Non-Ref") seqSetFilter(gds) ################################################### ### code chunk number 21: SeqVarTools.Rnw:250-253 ################################################### pc <- pca(gds) names(pc) plot(pc$eigenvect[,1], pc$eigenvect[,2]) ################################################### ### code chunk number 22: SeqVarTools.Rnw:260-268 ################################################### pval <- hwe(gds) pval <- -log10(sort(pval)) n <- length(pval) a <- 1:n b <- a/n x <- -log10(b) plot(x, pval, xlab="-log10(expected P)", ylab="-log10(observed P)") abline(0,1,col="red") ################################################### ### code chunk number 23: SeqVarTools.Rnw:273-275 ################################################### f <- inbreedCoeff(gds, margin="by.variant") hist(f) ################################################### ### code chunk number 24: SeqVarTools.Rnw:279-281 ################################################### ic <- inbreedCoeff(gds, margin="by.sample") range(ic) ################################################### ### code chunk number 25: SeqVarTools.Rnw:289-294 ################################################### data(pedigree) pedigree[pedigree$family == 1463,] err <- mendelErr(gds, pedigree, verbose=FALSE) table(err$by.variant) err$by.trio ################################################### ### code chunk number 26: SeqVarTools.Rnw:299-300 ################################################### seqClose(gds) ################################################### ### code chunk number 27: SeqVarTools.Rnw:303-304 ################################################### unlink(gdsfile)