### R code from vignette source 'D6_RegulatoryVariants.Rnw'

###################################################
### code chunk number 1: setup
###################################################
library(Morgan2013)
library(org.Hs.eg.db)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(BSgenome.Hsapiens.UCSC.hg19)
library(seqnames.db)
library(gmapR)
library(VariantTools)
library(lattice)
library(MotifDb)
library(parallel)


###################################################
### code chunk number 2: parallel
###################################################
library(parallel); options(mc.cores=detectCores())


###################################################
### code chunk number 3: annotation-libraries
###################################################
library(org.Hs.eg.db)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(BSgenome.Hsapiens.UCSC.hg19)


###################################################
### code chunk number 4: path-to-bam
###################################################
path <- "~/BigData/TERT/bam"


###################################################
### code chunk number 5: sanity-check
###################################################
stopifnot(length(dir(path, "bam$")) == 11)
stopifnot(length(dir(path, "bai$")) == 11)


###################################################
### code chunk number 6: consistent-names
###################################################
seqnameStyle(BSgenome.Hsapiens.UCSC.hg19) <- 
    seqnameStyle(TxDb.Hsapiens.UCSC.hg19.knownGene) <- "NCBI"


###################################################
### code chunk number 7: anno-symbols
###################################################
bsgenome <- BSgenome.Hsapiens.UCSC.hg19
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
geneid <- select(org.Hs.eg.db, "TERT", "ENTREZID", "SYMBOL")


###################################################
### code chunk number 8: VariantTools-setup
###################################################
VARIANTTOOLS_OK <- local({
    sys <- !Sys.info()["sysname"] %in% c("Darwin", "Windows")
    sys && require(VariantTools) && require(gmapR)
})
stopifnot(VARIANTTOOLS_OK)


###################################################
### code chunk number 9: mini-genome
###################################################
roi <- GRanges("5", IRanges(1, seqlengths(bsgenome)[["5"]]))
chr5seq <- getSeq(bsgenome, roi)
names(chr5seq) <- "5" 


###################################################
### code chunk number 10: make-mini
###################################################
library(gmapR)
genome.5 <- GmapGenome(chr5seq, name="hg19_5", create=TRUE)


###################################################
### code chunk number 11: TERT-promoters
###################################################
rng <- transcriptsBy(txdb, "gene")[[geneid$ENTREZID]]


###################################################
### code chunk number 12: TERT-starts
###################################################
unique(start(rng))


###################################################
### code chunk number 13: TERT-promoter
###################################################
unique(start(rng))
pregion <- promoters(range(rng), upstream=330, downstream=0)
seqlevels(pregion) <- "5"


###################################################
### code chunk number 14: setup-variant-calling
###################################################
library(VariantTools)
vtparam <- VariantTallyParam(genome.5, readlen=101L, 
    which=pregion, indels=TRUE)


###################################################
### code chunk number 15: call-variants
###################################################
fls <- dir(path, "bam$", full=TRUE)
called <- lapply(fls, callVariants, tally.param=vtparam)


###################################################
### code chunk number 16: call-variants-tidy
###################################################
len <- elementLengths(called)
called <- do.call(c, called)
id <- sub(".TERT.bam", "", basename(fls))
called$id <- factor(rep(id, len), levels=unique(id))


###################################################
### code chunk number 17: load-called
###################################################
data(called, package="Morgan2013")


###################################################
### code chunk number 18: explore-variants
###################################################
library(lattice)
altCounts <-
    xtabs(cycleCount.10.91 ~ start(called) + id, mcols(called))
plt <- levelplot(altCounts, xlab="Position", ylab=NULL,
          scales=list(x=list(rot=45)), aspect="fill",
          col.regions=rev(gray.colors(100, 0, 1)))


###################################################
### code chunk number 19: variant-levelplot
###################################################
pdf("variant-levelplot.pdf", width=10, height=3)
print(plt)
xx <- dev.off()


###################################################
### code chunk number 20: snp-def
###################################################
snp <- GRanges("5", IRanges(1295228, width=1))
snp <- flank(snp, 10, both=TRUE)


###################################################
### code chunk number 21: ref-seq
###################################################
refSeq <- getSeq(bsgenome, snp)


###################################################
### code chunk number 22: variantSequences
###################################################
library(Morgan2013)
variantSequences
altSeq <- variantSequences(called, snp, refSeq)


###################################################
### code chunk number 23: alt-seq
###################################################
altConsensus <- DNAStringSet(consensusString(altSeq))


###################################################
### code chunk number 24: complement
###################################################
complement(c(refSeq, altConsensus))


###################################################
### code chunk number 25: pwms
###################################################
library(MotifDb)
idx <- with(mcols(MotifDb), 
            organism=="Hsapiens" & dataSource == "JASPAR_CORE")
jasparHumanPWMs <- MotifDb[idx]


###################################################
### code chunk number 26: pwm-match
###################################################
## matchPWMs: helper in Morgan2013
minScore = "90%"             # high min. matching score: 1 base change
(refHits <- matchPWMs(jasparHumanPWMs, reverseComplement(refSeq)[[1]],
                      minScore))
(altHits <- matchPWMs(jasparHumanPWMs, reverseComplement(altConsensus)[[1]],
                      minScore))