## ----global_options, include=FALSE--------------------------------------------
knitr::opts_chunk$set(warning=FALSE, message=FALSE)

## ----load_package-------------------------------------------------------------
library(scoreInvHap)

## ----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

## -----------------------------------------------------------------------------
check <- checkSNPs(vcf)
check
vcf <- check$genos

## ----classify-----------------------------------------------------------------
res <- scoreInvHap(SNPlist = vcf, inv = "inv7_005")
res

## ----classify_par, eval=FALSE-------------------------------------------------
# res <- scoreInvHap(SNPlist = vcf, inv = "inv7_005",
#                    BPPARAM = BiocParallel::MulticoreParam(8))

## ----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))

## -----------------------------------------------------------------------------
## Inversion classification
table(classification(res))
## Haplotype classification
table(classification(res, inversion = FALSE))

## ----classify imputed---------------------------------------------------------
res_imp <- scoreInvHap(SNPlist = vcf, inv = "inv7_005", probs = TRUE)
res_imp

## ----compare classifications--------------------------------------------------
table(PostProbs = classification(res_imp), 
      BestGuess = classification(res))

## -----------------------------------------------------------------------------
data(inversionGR)
inversionGR

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