## ----global_options, include=FALSE----------------------------------------- knitr::opts_chunk$set(warning=FALSE, message=FALSE) ## ---- load_package--------------------------------------------------------- library(scoreInvHap) ## ---- load data------------------------------------------------------------ data("Refs") names(Refs) ## ---- Refs----------------------------------------------------------------- ref <- Refs$inv8p23.1 class(ref) ref[1:2] ## ---- SNPsR2--------------------------------------------------------------- data("SNPsR2") names(SNPsR2) R2s <- SNPsR2$inv8p23.1 head(R2s) ## ---- hetRefs-------------------------------------------------------------- data("hetRefs") names(hetRefs) hRefs <- hetRefs$inv8p23.1 head(hRefs) ## ---- eval=FALSE----------------------------------------------------------- # library(snpStats) # # ## From a bed # snps <- read.plink("example.bed") # # ## From a pedfile # snps <- read.pedfile("example.ped", snps = "example.map") ## ---- Load SNPs, message=FALSE--------------------------------------------- library(VariantAnnotation) vcf_file <- system.file("extdata", "example.vcf", package = "scoreInvHap") vcf <- readVcf(vcf_file, "hg19") vcf ## ---- classify------------------------------------------------------------- res <- scoreInvHap(SNPlist = vcf, SNPsR2 = SNPsR2$inv7p11.2, hetRefs = hetRefs$inv7p11.2, Refs = Refs$inv7p11.2) res ## ---- scoreInvHap results-------------------------------------------------- # Get classification head(classification(res)) # Get scores head(scores(res)) ## ---- scoreInvHap scores--------------------------------------------------- # Get max score head(maxscores(res)) # Get difference score head(diffscores(res)) ## -------------------------------------------------------------------------- plotScores(res, pch = 16, main = "QC based on scores") ## -------------------------------------------------------------------------- # Get Number of scores used head(numSNPs(res)) # Get call rate head(propSNPs(res)) ## -------------------------------------------------------------------------- plotCallRate(res, main = "Call Rate QC") ## -------------------------------------------------------------------------- ## No filtering length(classification(res)) ## QC filtering length(classification(res, minDiff = 0.1, callRate = 0.9)) ## -------------------------------------------------------------------------- ## No filtering table(classification(res)) ## QC filtering table(classification(res, inversion = TRUE)) ## ---- classify imputed----------------------------------------------------- res_imp <- scoreInvHap(SNPlist = vcf, SNPsR2 = SNPsR2$inv7p11.2, hetRefs = hetRefs$inv7p11.2, Refs = Refs$inv7p11.2, imputed = TRUE) res_imp ## ---- compare classifications---------------------------------------------- table(PostProbs = classification(res_imp), BestGuess = classification(res)) ## -------------------------------------------------------------------------- sessionInfo()