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

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex()


###################################################
### code chunk number 2: code
###################################################
ps.options(font="serif")
choptxt = function(x, nwpl=12) {
 stopifnot(length(x)==1)
 x = strsplit(x, " ")[[1]]
 nx = length(x)
 chunks = ceiling((1:nx)/nwpl)
 x = split(x, chunks)
 sapply(x, function(y) paste(y, collapse=" ", sep=""))
}

genemodel = function(sym, genome="hg19") {
 stopifnot(genome=="hg19")
 require(TxDb.Hsapiens.UCSC.hg19.knownGene)
 require(org.Hs.eg.db)
 num = get(sym, revmap(org.Hs.egSYMBOL))
 exonsBy(txdb, by="gene")[[num]]
}



###################################################
### code chunk number 3: lka
###################################################
library(harbChIP)
if (!exists("harbChIP")) data(harbChIP)
cat("     ", choptxt(abstract(harbChIP),11),sep="\n   ")


###################################################
### code chunk number 4: lkh (eval = FALSE)
###################################################
## library(harbChIP)
## data(harbChIP)


###################################################
### code chunk number 5: doe
###################################################
harbChIP


###################################################
### code chunk number 6: lkpro
###################################################
sampleNames(harbChIP)[1:10]


###################################################
### code chunk number 7: lkv
###################################################
library(viz14)
ace2 = makebs("ACE2")
ace2
QQnorm(ace2)


###################################################
### code chunk number 8: lk2
###################################################
ace2bnd = boundGenes(ace2)


###################################################
### code chunk number 9: ans
###################################################
facs = c("ACE2", "SWI5", "SWI6", "SWI4", "MBP1", 
   "FKH1", "FKH2", "NDD1", "MCM1")
bg = lapply(facs, function(x) boundGenes(makebs(x)))
names(bg) = facs
sapply(bg,length)


###################################################
### code chunk number 10: lkven (eval = FALSE)
###################################################
## library(VennDiagram)
## v1 = venn.diagram( bg[1:5], filename=NULL )
## grid.draw(v1)


###################################################
### code chunk number 11: lkd
###################################################
data(trigFits)
summary(trigFits[, "dtf"])


###################################################
### code chunk number 12: dolapp
###################################################
facs = c("ACE2", "SWI5", "SWI6", "SWI4", "MBP1", 
   "FKH1", "FKH2", "NDD1", "MCM1")
bg = lapply(facs, function(x) boundGenes(makebs(x)))
bgrps = lapply(bg, function(x) trigFits[
   intersect(x, rownames(trigFits)), "dtf" ] ) 
names(bgrps) = facs
sapply(bgrps, function(x) length(na.omit(x)))


###################################################
### code chunk number 13: dobpbg
###################################################
boxplot(bgrps, las=2)


###################################################
### code chunk number 14: getd
###################################################
library(yeastCC)
data(spYCCES)
spYCCES
experimentData(spYCCES)


###################################################
### code chunk number 15: lkd
###################################################
data(orf800)
orf800[1:4]


###################################################
### code chunk number 16: lkabs
###################################################
cat("     ", choptxt(abstract(spYCCES),11),sep="\n   ")


###################################################
### code chunk number 17: lka
###################################################
alp = spYCCES[ , spYCCES$syncmeth=="alpha"]
alp
table(alp$time)


###################################################
### code chunk number 18: lkg
###################################################
yal040c = exprs(alp)["YAL040C",]
df = data.frame(yal040c, time=alp$time)
plot(yal040c~time, data=df)


###################################################
### code chunk number 19: donl
###################################################
m1 = nls(yal040c~b*sin(d+a*time),data=df,start=list(d=.1,b=1,a=.1))
m1


###################################################
### code chunk number 20: sol (eval = FALSE)
###################################################
## ptime = seq(0,120,.1)
## pex = predict(m1, newdata=list(time=ptime))
## plot(pex~ptime, type="l")
## points(alp$time, yal040c)


###################################################
### code chunk number 21: sol2 (eval = FALSE)
###################################################
## res = resid(m1)
## plot(res~alp$time)


###################################################
### code chunk number 22: sol3 (eval = FALSE)
###################################################
## library(ggplot2)
## prdf = data.frame(pred=pex,time=ptime)
## g1 = ggplot(prdf, aes(x=time,y=pred))
## print(g1 + geom_line() + geom_point(data=df, aes(y=yal040c, x=time)) + 
##    stat_smooth(data=df, aes(y=yal040c, x=time)))


###################################################
### code chunk number 23: newm
###################################################
df$ptime = 2*pi*(df$time %% 64)/64
m2 = lm(yal040c ~ sin(ptime) + cos(ptime) -1, data=df)
summary(m2)
sum(resid(m1)^2)
sum(resid(m2)^2)


###################################################
### code chunk number 24: domod
###################################################
gettrm = function(genename, es, period=64) {
 stopifnot("time" %in% names(pData(es)))
 ex = exprs(es)[genename,]
 et = es$time
 ptime = 2*pi*(df$time %% period)/period
 ndf = data.frame(time=et, ptime=ptime)
 ndf[[tolower(substitute(genename))]] = exprs(es)[genename,]
 fm = as.formula(paste(tolower(substitute(genename)), "~ sin(ptime) + cos(ptime) - 1",
    sep=" "))
 lm(fm, data=ndf)
}


###################################################
### code chunk number 25: lkch
###################################################
gettrm("YAL040C", alp)


###################################################
### code chunk number 26: getmat
###################################################
outs = matrix(NA, nc=3, nr=nrow(alp))
yg = featureNames(alp)
rownames(outs) = yg
colnames(outs) = c("msep", "amp", "phase")
suppressMessages({
for (i in 1:nrow(alp)) {
  curg = force(yg[i])
  m = try( gettrm( curg, alp ) )
  if (inherits(m, "try-error")) next
  cm = coef(m)
  outs[i,"msep"] = mean(resid(m)^2)
  outs[i,"amp"] = sqrt(sum(cm^2))
  outs[i,"phase"] = atan(-cm[1]/cm[2])
  }
})


###################################################
### code chunk number 27: lkalu
###################################################
library(viz14)
data(nestrep)
metadata(nestrep)[[2]]$RDataName


###################################################
### code chunk number 28: dolim
###################################################
soma = paste0("chr", 1:22)
alu = nestrep[ grep("^Alu", nestrep$name) ]
library(GenomicRanges)
alu = alu[ which(as.character(seqnames(alu)) %in% soma) ]
alu14 = alu[ which(seqnames(alu)=="chr14") ]
seqlevels(alu) = soma
alus = alu[ sample(1:length(alu), size=5000) ]
library(ggbio)
autoplot(alus, layout="karyogram")


###################################################
### code chunk number 29: dotx14
###################################################
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
ucg = genes(txdb)
ucg14 = ucg[which(seqnames(ucg)=="chr14")]
egovalu = as.character(ucg14[which(ucg14 %over% alu14)]$gene_id)
library(org.Hs.eg.db)
symovalu = unique(unlist(mget(egovalu, org.Hs.egSYMBOL, ifnotfound=NA)))
allg14 = get("14", revmap(org.Hs.egCHR))
alls14 = unlist(mget(allg14, org.Hs.egSYMBOL))
symnoalu = setdiff(alls14, symovalu)


###################################################
### code chunk number 30: dogwa
###################################################
library(gwascat)
gwrngs  # limit to chr14 :
gwr14 = gwrngs[which(seqnames(gwrngs)=="chr14")]


###################################################
### code chunk number 31: getrate
###################################################
mapped14 = unique(gwr14$Mapped_gene)
table(mapped14)[1:10]


###################################################
### code chunk number 32: esti
###################################################
mean(symovalu %in% mapped14)  # very crude


###################################################
### code chunk number 33: doac
###################################################
library(GenomicFiles)
library(Gviz)

GENENAME = "ADSSL1"    # gene to analyze
gm = genemodel(GENENAME)
rgm = range(gm)
at = AnnotationTrack(gm, genome="hg19", name=GENENAME)

#  now set up the references to RNA-seq data
fn = dir(system.file(  # use Herve's package
  "extdata", package="RNAseqData.HNRNPC.bam.chr14"), full=TRUE,
  patt="bam$")
bfv = BamFileViews(fn)  # lightweight reference

# now ingest the BAM data into AlignmentsTrack instances
STACKTYPE = "hide"      # for Gviz
kd1a = AlignmentsTrack( path(fileList(bfv)[[1]]), isPaired=TRUE, 
    name="KD1a", chromosome="chr14", stacking=STACKTYPE )
wt1a = AlignmentsTrack( path(fileList(bfv)[[5]]), isPaired=TRUE, 
    name="WT1a", chromosome="chr14", stacking=STACKTYPE )
kd2a = AlignmentsTrack( path(fileList(bfv)[[3]]), isPaired=TRUE, 
    name="KD2a", chromosome="chr14", stacking=STACKTYPE )
wt2a = AlignmentsTrack( path(fileList(bfv)[[7]]), isPaired=TRUE, 
    name="WT2a", chromosome="chr14", stacking=STACKTYPE )

gt = GenomeAxisTrack(genome="hg19")
#data(nestrep)
#rep14 = nestrep[which(seqnames(nestrep)=="chr14")] 
#alu14 = rep14[ grep("^Alu", rep14$name) ]
alut = AnnotationTrack(alu14, genome="hg19", name="Alu")
ps.options(font="sans")
plotTracks(list(kd1a, kd2a, wt1a, wt2a, 
    alut, at, gt), from=start(rgm), to=end(rgm))


###################################################
### code chunk number 34: setupGGv (eval = FALSE)
###################################################
## library(ggvis)
## library(GenomicFiles)
## library(viz14)
## data(alu14)
## alu14 = alu14[,1:3]
## fn = dir(system.file(
##   "extdata", package="RNAseqData.HNRNPC.bam.chr14"), full=TRUE,
##   patt="bam$")
## bfv = BamFileViews(fn)
## browseWreads()


###################################################
### code chunk number 35: lkb
###################################################
browseWreads


###################################################
### code chunk number 36: L1 (eval = FALSE)
###################################################
##     base = renderGene(sym, cachedir = cachedir, y = y, y2 = y2, 
##         ylim = ylim, padpct = padpct)


###################################################
### code chunk number 37: lkrn (eval = FALSE)
###################################################
##     rng = rangeOfGene(sym, ggmDefC())
##     wid = end(rng) - start(rng)
##     pad = padpct * wid/100
##     left = start(rng) - pad


###################################################
### code chunk number 38: lkal (eval = FALSE)
###################################################
##         return(base %>% rectsAtY(sym = sym, gr = alugr, y = aluy) %>% 
##             addText(sym = "Alu", x = left, y = aluy) %>% readsAtY(sym = sym, 
##             bamf = bamf, y = ready) %>% addText(sym = "reads", 
##             x = left, y = ready))


###################################################
### code chunk number 39: lkren
###################################################
renderGene


###################################################
### code chunk number 40: lkses
###################################################
sessionInfo()