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

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


###################################################
### code chunk number 2: roi
###################################################
p53 <- gmapR:::exonsOnTP53Genome("TP53")


###################################################
### code chunk number 3: callVariants
###################################################
library(VariantTools)
bams <- LungCancerLines::LungCancerBamFiles()
bam <- bams$H1993
tally.param <- VariantTallyParam(gmapR::TP53Genome(), 
                                 readlen = 100L,
                                 high_base_quality = 23L,
                                 which = range(p53))
called.variants <- callVariants(bam, tally.param)


###################################################
### code chunk number 4: callVariants-exonic
###################################################
subsetByOverlaps(called.variants, p53, ignore.strand = TRUE)


###################################################
### code chunk number 5: tallyVariants
###################################################
raw.variants <- tallyVariants(bam, tally.param)


###################################################
### code chunk number 6: qaVariants
###################################################
qa.variants <- qaVariants(raw.variants)


###################################################
### code chunk number 7: callVariants
###################################################
called.variants <- callVariants(qa.variants)


###################################################
### code chunk number 8: VariantQAFilters
###################################################
qa.filters <- VariantQAFilters()


###################################################
### code chunk number 9: VariantQAFilters-summary
###################################################
summary(qa.filters, raw.variants)


###################################################
### code chunk number 10: VariantQAFilters-subset
###################################################
qa.variants <- subsetByFilter(raw.variants, qa.filters)


###################################################
### code chunk number 11: VariantQAFilters-fisherStrandP
###################################################
qa.filters.custom <- VariantQAFilters(fisher.strand.p.value = 1e-4)
summary(qa.filters.custom, raw.variants)


###################################################
### code chunk number 12: VariantQAFilters-compare
###################################################
fs.original <- eval(qa.filters["fisherStrand"], raw.variants) 
fs.custom <- eval(qa.filters.custom["fisherStrand"], raw.variants)
raw.variants[fs.original != fs.custom]


###################################################
### code chunk number 13: VariantCallingFilters
###################################################
calling.filters <- VariantCallingFilters()
summary(calling.filters, qa.variants)


###################################################
### code chunk number 14: variantGR2VCF
###################################################
vcf <- variantGR2Vcf(called.variants, sample.id = "H1993",
                     project = "VariantTools_Vignette")


###################################################
### code chunk number 15: writeVcf (eval = FALSE)
###################################################
## writeVcf(vcf, "H1993.vcf", index = TRUE)