## ---- loadPackage, message=FALSE, warnings=FALSE------------------------- library(genotypeeval) ## ---- loadVCF, message=FALSE, warnings=FALSE----------------------------- vcffn <- system.file("ext-data", "chr22.GRCh38.vcf.gz", package="genotypeeval") mydir <- paste(dirname(vcffn), "/", sep="") myfile <-basename(vcffn) vcf <- ReadVCFData(mydir, myfile, "GRCh38") ## ---- getCoding, message=FALSE, warnings=FALSE--------------------------- suppressWarnings(library("TxDb.Hsapiens.UCSC.hg38.knownGene")) txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene cds <- cds(txdb) seqlevelsStyle(cds) = "NCBI" reg.chrs <- c(as.character(seq(1:22)), "X", "Y") cds <- keepSeqlevels(cds, reg.chrs) genome(cds) <- "GRCh38" ## ---- loadGold, message=FALSE, warnings=FALSE, message=FALSE------------- suppressWarnings(library("SNPlocs.Hsapiens.dbSNP141.GRCh38")) snps <- SNPlocs.Hsapiens.dbSNP141.GRCh38 dbsnp.ch22 <- snplocs(snps, "ch22", as.GRanges=TRUE) seqlevelsStyle(dbsnp.ch22) = "NCBI" gold.param = GoldDataParam(percent.confirmed=0.8) gold <- GoldDataFromGRanges("GRCh38", dbsnp.ch22, gold.param) ## ---- prepAdmix, message=FALSE, warnings=FALSE--------------------------- admix.var <- getVR(vcf)[getVR(vcf)$GT %in% c("0|1", "1|0", "1|1"),][,1:2] admix.var$EAS_AF <- ifelse(admix.var$GT %in% c("1|1"), 1, .5) admix.var$AFR_AF<- 0 admix.var$EUR_AF<- 0 admix.hom <- getVR(vcf)[getVR(vcf)$GT %in% c("0|0"),][,1:2] admix.hom$EAS_AF<- 0 admix.hom$AFR_AF<- 1 admix.hom$EUR_AF<- 1 admix.ref <- c(admix.var, admix.hom) ## ---- runVCF------------------------------------------------------------- vcfparams <- VCFQAParam(count.limits=c(3e9, Inf), titv.noncoding=c(1.99, 4), titv.coding=c(2.8, 4)) ## ---- evaluatenoref, message=FALSE, warning=FALSE------------------------ ev.noref <- VCFEvaluate(vcf, vcfparams) ## ---- evaluate, message=FALSE, warning=FALSE----------------------------- ev <- VCFEvaluate(vcf, vcfparams, gold.ref=gold, cds.ref=cds, admixture.ref = admix.ref) ## ---- passed------------------------------------------------------------- didSamplePassOverall(ev) ## ---- passedDetails------------------------------------------------------ head(didSamplePass(ev)) ## ---- results------------------------------------------------------------ head(getResults(ev)) ## ---- results2----------------------------------------------------------- getResults(ev)[c("admixture.EAS", "admixture.EUR", "admixture.AFR")] ## ---- plots-------------------------------------------------------------- getPlots(ev)$variant_type ## ---- plots2------------------------------------------------------------- getPlots(ev)$chr