inst/doc/waveTiling-vignette.R

In waveTiling: Wavelet-Based Models for Tiling Array Transcriptome Analysis

### R code from vignette source 'waveTiling-vignette.Rnw'

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


###################################################
### code chunk number 2: waveTiling-vignette.Rnw:51-55
###################################################
library(waveTiling)
library(waveTilingData)
library(pd.atdschip.tiling)
data(leafdev)


###################################################
### code chunk number 3: waveTiling-vignette.Rnw:60-65
###################################################
leafdev <- as(leafdev,"WaveTilingFeatureSet")
leafdev <- addPheno(leafdev,noGroups=6,
	groupNames=c("day8","day9","day10","day11","day12","day13"),
	replics=rep(3,6))
leafdev


###################################################
### code chunk number 4: waveTiling-vignette.Rnw:68-70
###################################################
rm(leafdev)
gc()


###################################################
### code chunk number 5: waveTiling-vignette.Rnw:75-83
###################################################
library(BSgenome.Athaliana.TAIR.TAIR9)

# leafdevMapAndFilterTAIR9 <- filterOverlap(leafdev,remap=TRUE,
#	BSgenomeObject=Athaliana,chrId=1:7,
#	strand="both",MM=FALSE)
data(leafdevMapAndFilterTAIR9)

# leafdevBQ <- bgCorrQn(leafdev,useMapFilter=leafdevMapAndFilterTAIR9)


###################################################
### code chunk number 6: waveTiling-vignette.Rnw:91-100
###################################################
data(leafdevBQ)
chromosome <- 1
strand <- "forward"
leafdevFit <- wfm.fit(leafdevBQ,filter.overlap=leafdevMapAndFilterTAIR9,
	design="time",n.levels=10,
	chromosome=chromosome,strand=strand,minPos=22000000,
	maxPos=24000000,var.eps="marg",prior="improper",
	skiplevels=1,save.obs="plot",trace=TRUE)
leafdevFit


###################################################
### code chunk number 7: waveTiling-vignette.Rnw:107-111
###################################################
delta <- log(1.2,2)
leafdevInfCompare <- wfm.inference(leafdevFit,
	contrasts="compare",delta=c("median",delta))
leafdevInfCompare


###################################################
### code chunk number 8: waveTiling-vignette.Rnw:118-121
###################################################
sigGenomeRegionsCompare <- getGenomicRegions(leafdevInfCompare)
sigGenomeRegionsCompare[[2]]
length(sigGenomeRegionsCompare)


###################################################
### code chunk number 9: waveTiling-vignette.Rnw:128-135
###################################################
library(TxDb.Athaliana.BioMart.plantsmart22)
sigGenesCompare <- getSigGenes(fit=leafdevFit,inf=leafdevInfCompare,
	biomartObj=TxDb.Athaliana.BioMart.plantsmart22)
head(sigGenesCompare[[2]])
nonAnnoCompare <- getNonAnnotatedRegions(fit=leafdevFit,inf=leafdevInfCompare,
	biomartObj=TxDb.Athaliana.BioMart.plantsmart22)
head(nonAnnoCompare[[2]])


###################################################
### code chunk number 10: waveTiling-vignette.Rnw:140-145
###################################################
leafdevInfTimeEffect <- wfm.inference(leafdevFit,contrasts="effects",
	delta=c("median",2,0.2,0.2,0.2,0.2))

leafdevInfMeans <- wfm.inference(leafdevFit,contrasts="means",
	delta=4,minRunPos=30,minRunProbe=-1)


###################################################
### code chunk number 11: waveTiling-vignette.Rnw:150-176
###################################################
custDes <- matrix(0,nrow=18,ncol=6)
orderedFactor <- factor(1:6,ordered=TRUE)
desPoly <- lm(rnorm(6)~orderedFactor,x=TRUE)$x
custDes[,1] <- 1
custDes[,2:6] <- apply(desPoly[,2:6],2,rep,getReplics(leafdevBQ))
custDes

leafdevFitCustom <- wfm.fit(leafdevBQ,filter.overlap=leafdevMapAndFilterTAIR9,
	design="custom",design.matrix=custDes,n.levels=10,
	chromosome=chromosome,strand=strand,minPos=22000000,
	maxPos=24000000,var.eps="marg",prior="improper",
	skiplevels=1,save.obs="plot",trace=TRUE)

noGroups <- getNoGroups(leafdevBQ)
myContrastMat <- matrix(0,nrow=noGroups*(noGroups-1)/2,ncol=noGroups)
hlp1 <- rep(2:noGroups,1:(noGroups-1))
hlp2 <- unlist(sapply(1:(noGroups-1),function(x) seq(1:x)))
for (i in 1:nrow(myContrastMat))
{
	myContrastMat[i,hlp1[i]] <- 1
	myContrastMat[i,hlp2[i]] <- -1
}
myContrastMat

leafdevInfCustom <- wfm.inference(leafdevFitCustom,contrast.matrix=myContrastMat,
	delta=c("median",log(1.2,2)))


###################################################
### code chunk number 12: waveTiling-vignette.Rnw:184-192
###################################################
trs <- transcripts(TxDb.Athaliana.BioMart.plantsmart22)
sel <- trs[elementMetadata(trs)$tx_name %in% "AT1G62500.1",]
start <- start(ranges(sel))-2000
end <- end(ranges(sel))+2000
plotWfm(fit=leafdevFit,inf=leafdevInfCompare,
	biomartObj=TxDb.Athaliana.BioMart.plantsmart22,
	minPos=start,maxPos=end,two.strand=TRUE,
	plotData=TRUE,plotMean=FALSE,tracks=c(1,2,6,10,11))


###################################################
### code chunk number 13: waveTiling-vignette.Rnw:195-201
###################################################
start <- start(ranges(sel))-4000
end <- end(ranges(sel))+4000
plotWfm(fit=leafdevFit,inf=leafdevInfTimeEffect,
	biomartObj=TxDb.Athaliana.BioMart.plantsmart22,
	minPos=start,maxPos=end,two.strand=TRUE,
	plotData=TRUE,plotMean=FALSE,tracks=1)


###################################################
### code chunk number 14: waveTiling-vignette.Rnw:204-208
###################################################
plotWfm(fit=leafdevFit,inf=leafdevInfMeans,
	biomartObj=TxDb.Athaliana.BioMart.plantsmart22,
	minPos=start,maxPos=end,two.strand=TRUE,
	plotData=TRUE,plotMean=FALSE,tracks=1:6)


###################################################
### code chunk number 15: waveTiling-vignette.Rnw:216-234
###################################################
getGenomeInfo(leafdevFit)
dataOrigSpace <- getDataOrigSpace(leafdevFit)
dim(dataOrigSpace)
dataOrigSpace[1:8,1:8]
dataWaveletSpace <- getDataWaveletSpace(leafdevFit)
dim(dataWaveletSpace)
dataWaveletSpace[1:8,1:8]
getDesignMatrix(leafdevFit)
probepos <- getProbePosition(leafdevFit)
length(probepos)
head(probepos)

effects <- getEff(leafdevInfCompare)
dim(effects)
effects[1:8,1:8]
fdrs <- getFDR(leafdevInfCompare)
dim(fdrs)
fdrs[1:8,1:8]