### R code from vignette source 'vignettes/VariantAnnotation/inst/doc/VariantAnnotation.Rnw'

###################################################
### code chunk number 1: options
###################################################
options(width=72)


###################################################
### code chunk number 2: data
###################################################
library(VariantAnnotation)
data(variants)
head(variants)


###################################################
### code chunk number 3: location
###################################################
library(TxDb.Hsapiens.UCSC.hg18.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg18.knownGene 
loc <- locateVariants(variants, txdb)
head(loc)


###################################################
### code chunk number 4: singleSNP
###################################################
singleSNP <- GRanges(seqnames="chr16",
        ranges=IRanges(start=c(23004226, 22316999, 24133766),
               width=c(1,1,1)),
        varAllele=DNAStringSet("T"),
        comments=c("Celera assembly", "HuRef assembly", "reference assembly"))

singleSNP


###################################################
### code chunk number 5: aacoding
###################################################
library(BSgenome.Hsapiens.UCSC.hg18)
cdsByTx <- cdsBy(txdb)
grl <- keepSeqlevels(cdsByTx, c("chr1", "chr2", "chr16"))

aaCoding <- predictCoding(variants, grl, seqSource=Hsapiens,
    varAllele="varAllele")
head(aaCoding)


###################################################
### code chunk number 6: aacodingResultOneVariant
###################################################
aaCoding[width(aaCoding$varAA) == 0,] 


###################################################
### code chunk number 7: createFilters
###################################################
library("SNPlocs.Hsapiens.dbSNP.20090506")
## create filter to identify variants present in dbSNP 
snpFilt <- dbSNPFilter("SNPlocs.Hsapiens.dbSNP.20090506")

## create filter to identify variants present in introns 
regionFilt <- regionFilter(txdb, region="intron")
regionFilt


###################################################
### code chunk number 8: filtering
###################################################
## return the subset of records that passed the filter
snpResult <- snpFilt(variants)
snpResult 
metadata(snpResult)

## return a VAFilterResult object
snpFilt(variants, subset=FALSE)


###################################################
### code chunk number 9: readVCF
###################################################
## compress vcf file , create index
vcfFile <- system.file("extdata", "ex1.vcf", package="VariantAnnotation")
from <- vcfFile
to <- tempfile()
compressVcf <- bgzip(from, to)
idx <- indexTabix(compressVcf, "vcf")
tab <- TabixFile(compressVcf, idx)
head <- headerTabix(tab)

## select ranges with param 
gr <- GRanges(seqnames="1",
              ranges=IRanges(start=0, end=100000))
vcf <- readVcf(tab, param=gr)


###################################################
### code chunk number 10: readVCF_show
###################################################
vcf

exptData(vcf)[["HEADER"]]

exptData(vcf)[["HEADER"]][["META"]]


###################################################
### code chunk number 11: readVCF_rowData
###################################################
exptData(vcf)[["HEADER"]][["FILTER"]]
exptData(vcf)[["HEADER"]][["INFO"]]

rowData(vcf)

## ALT as a DNAStringSetList
head(as.list(values(rowData(vcf))[["ALT"]]), 4)


###################################################
### code chunk number 12: readVCF_assays
###################################################
exptData(vcf)[["HEADER"]][["FORMAT"]]

assays(vcf)

head(assays(vcf)$GT)


###################################################
### code chunk number 13: SIFT_and_PP
###################################################
library(SIFT.Hsapiens.dbSNP132)
library(PolyPhen.Hsapiens.dbSNP131)

## rsids in the package 
head(keys(SIFT.Hsapiens.dbSNP132))

## list available columns
cols(SIFT.Hsapiens.dbSNP132)

rsids <- c("rs2142947", "rs8692231", "rs3026284")
subst <- c("RSID", "PREDICTION", "SCORE", "AACHANGE", "PROTEINID")
select(SIFT.Hsapiens.dbSNP132, keys=rsids, cols=subst)

select(PolyPhen.Hsapiens.dbSNP131, keys="rs3026284", 
    cols=c("TRAININGSET", "PREDICTION", "PPH2PROB", "PPH2FPR", "PPH2FDR"))


###################################################
### code chunk number 14: sessionInfo
###################################################
sessionInfo()