inst/doc/scoreInvHap.R

In scoreInvHap: Get inversion status in predefined regions

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