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