[Wavetiling-commits] r19 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Feb 16 08:00:31 CET 2012
Author: ppipeler
Date: 2012-02-16 08:00:30 +0100 (Thu, 16 Feb 2012)
New Revision: 19
Added:
pkg/R/methods-WfmFit.R
pkg/R/methods-WfmInf.R
Removed:
pkg/R/methods-Wfm.R
Log:
Deleted: pkg/R/methods-Wfm.R
===================================================================
--- pkg/R/methods-Wfm.R 2012-02-15 16:55:39 UTC (rev 18)
+++ pkg/R/methods-Wfm.R 2012-02-16 07:00:30 UTC (rev 19)
@@ -1,540 +0,0 @@
-setMethod("getProbePosition",signature("Wfm"),function(object)
-{
- return(object at probePosition)
-})
-
-setMethod("getNoProbes",signature("Wfm"),function(object)
-{
- return(length(object at probePosition))
-})
-
-setMethod("getBetaMAP",signature("Wfm"),function(object)
-{
- return(object at betaMAP)
-})
-
-setMethod("getVarBetaMAP",signature("Wfm"),function(object)
-{
- return(object at varbetaMAP)
-})
-
-setMethod("getSmoothPar",signature("Wfm"),function(object)
-{
- return(object at smoothPar)
-})
-
-setMethod("getVarEps",signature("Wfm"),function(object)
-{
- return(object at varEps)
-})
-
-
-setMethod("getGenomeInfo",signature("Wfm"),function(object)
-{
- return(object at genome.info)
-})
-
-setMethod("getChromosome",signature("Wfm"),function(object)
-{
- return(object at genome.info@chromosome)
-})
-
-setMethod("getStrand",signature("Wfm"),function(object)
-{
- return(object at genome.info@strand)
-})
-
-setMethod("getMinPos",signature("Wfm"),function(object)
-{
- return(object at genome.info@minPos)
-})
-
-setMethod("getMaxPos",signature("Wfm"),function(object)
-{
- return(object at genome.info@maxPos)
-})
-
-setMethod("getNoLevels",signature("Wfm"),function(object)
-{
- return(object at n.levels)
-})
-
-setMethod("getWfmMethod",signature("Wfm"),function(object)
-{
- return(object at method)
-})
-
-setMethod("getDesign",signature("Wfm"),function(object)
-{
- return(object at design)
-})
-
-setMethod("getPhenoInfo",signature("Wfm"),function(object)
-{
- return(object at phenoData)
-})
-
-setMethod("getDataOrigSpace",signature("Wfm"),function(object)
-{
- return(object at dataOrigSpace)
-})
-
-setMethod("getDataWaveletSpace",signature("Wfm"),function(object)
-{
- return(object at dataWaveletSpace)
-})
-
-setMethod("getWaveletFilter",signature("Wfm"),function(object)
-{
- return(object at wave.filt)
-})
-
-setMethod("getKj",signature("Wfm"),function(object)
-{
- return(object at Kj)
-})
-
-setMethod("getPrior",signature("Wfm"),function(object)
-{
- return(object at prior)
-})
-
-setMethod("getAlpha",signature("Wfm"),function(object)
-{
- return(object at alpha)
-})
-
-setMethod("getDelta",signature("Wfm"),function(object)
-{
- return(object at delta)
-})
-
-setMethod("getTwoSided",signature("Wfm"),function(object)
-{
- return(object at two.sided)
-})
-
-setMethod("getRescale",signature("Wfm"),function(object)
-{
- return(object at rescale)
-})
-
-setMethod("getSigProbes",signature("Wfm"),function(object)
-{
- return(object at sigProbes)
-})
-
-setMethod("getRegions",signature("Wfm"),function(object)
-{
- return(object at regions)
-})
-
-setMethod("getGenomicRegions",signature("Wfm"),function(object)
-{
- return(object at GlocRegions)
-})
-
-setMethod("getFDR",signature("Wfm"),function(object)
-{
- return(object at FDR)
-})
-
-setMethod("getCI",signature("Wfm"),function(object)
-{
- return(object at CI)
-})
-
-setMethod("getF",signature("Wfm"),function(object)
-{
- return(object at F)
-})
-
-setMethod("getVarF",signature("Wfm"),function(object)
-{
- return(object at varF)
-})
-
-setMethod("getEff",signature("Wfm"),function(object)
-{
- return(object at eff)
-})
-
-setMethod("getVarEff",signature("Wfm"),function(object)
-{
- return(object at vareff)
-})
-
-setMethod("plotWfm",signature("Wfm"),function(object,annoFile,minPos,maxPos,trackFeature="exon",overlayFeature=c("gene","transposable_element_gene"),two.strand=TRUE,plotData=TRUE,plotMean=TRUE,tracks=0)
-# complete tracks option for different methods
-{
- Gloc <- getProbePosition(object)
- chromosome <- getChromosome(object)
- strand <- getStrand(object)
- selID <- (1:length(Gloc))[Gloc>minPos & Gloc<maxPos]
- sta <- min(selID)
- end <- max(selID)
- minBase <- Gloc[sta]
- maxBase <- Gloc[end]
- trackCount <- 1
- trackInfo <- list()
- overlayInfo <- list()
- if (plotData==TRUE)
- {
- Y <- getDataOrigSpace(object)
- replcs <- as.numeric(table(getPhenoInfo(object)$group))
- trackInfo[[trackCount]] <- makeGenericArray(intensity=as.matrix(t(Y[,sta:end])),probeStart=Gloc[sta:end],dp=DisplayPars(color=rep(1:length(replcs),replcs),ylim=range(Y[,sta:end]),pointSize=.3,pch=sequence(replcs)))
- names(trackInfo)[trackCount] <- "Data"
- trackCount <- trackCount + 1
- }
- if (two.strand==TRUE)
- {
- trackInfo[[trackCount]] <- makeNewAnnotationTrack(annoFile=annoFile,chromosome=chromosome,minBase=minBase,maxBase=maxBase,strand="forward",feature=trackFeature,dp=NULL)
- names(trackInfo)[trackCount] <- "F"
- overlayInfo[[trackCount]] <- makeNewAnnotationTextOverlay(annoFile=annoFile,chromosome=chromosome,minBase=minBase,maxBase=maxBase,strand="forward",region=c(trackCount,trackCount),feature=overlayFeature,y=0.5)
- trackCount <- trackCount + 1
- trackInfo[[trackCount]] <- makeGenomeAxis(add53 = TRUE,add35 = TRUE)
- trackCount <- trackCount + 1
- trackInfo[[trackCount]] <- makeNewAnnotationTrack(annoFile=annoFile,chromosome=chromosome,minBase=minBase,maxBase=maxBase,strand="reverse",feature=trackFeature)
- names(trackInfo)[trackCount] <- "R"
- overlayInfo[[trackCount]] <- makeNewAnnotationTextOverlay(annoFile=annoFile,chromosome=chromosome,minBase=minBase,maxBase=maxBase,strand="reverse",region=c(trackCount,trackCount),feature=overlayFeature,y=0.5)
- trackCount <- trackCount + 1
- } else
- {
- trackInfo[[trackCount]] <- makeNewAnnotationTrack(annoFile=annoFile,chromosome=chromosome,minBase=minBase,maxBase=maxBase,strand=strand,feature=trackFeature)
- if (strand=="forward")
- {
- names(trackInfo)[trackCount] <- "F"
- }
- if (strand=="reverse")
- {
- names(trackInfo)[trackCount] <- "R"
- }
- overlayInfo[[trackCount]] <- makeNewAnnotationTextOverlay(annoFile=annoFile,chromosome=chromosome,minBase=minBase,maxBase=maxBase,strand=strand,region=c(trackcount,trackCount),feature=overlayFeature,y=0.5)
- trackCount <- trackCount + 1
- gAxis <- makeGenomeAxis(add53 = TRUE,add35 = TRUE)
- trackCount <- trackCount + 1
- }
- effects <- getEff(object)
- regions <- getRegions(object)
- noGroups <- length(table(getPhenoInfo(object)$group))
- if (getWfmMethod(object)=="meansByGroupTime" | getWfmMethod(object)=="meansByGroupFactor")
- {
- plotMean <- FALSE
- }
- if (plotMean==TRUE)
- {
- trackInfo[[trackCount]] <- makeGenericArray(intensity=as.matrix(effects[1,sta:end]),probeStart=Gloc[sta:end],dp=DisplayPars(color="black",ylim=range(effects[1,sta:end])+c(-0.4,0.4),pointSize=.3,pch=1,lwd=1,type="line"))
- names(trackInfo)[trackCount] <- "Mean"
- overlayInfo[[trackCount]] <- makeNewTranscriptRectangleOverlay(sigRegions=as.matrix(data.frame(start(regions[[1]]),end(regions[[1]]))),location=Gloc,start=sta,end=end,region=c(trackCount,trackCount),dp=DisplayPars(color="darkgrey",alpha=.1))
- trackCount <- trackCount + 1
- }
- if (getWfmMethod(object)=="twoGroup" | getWfmMethod(object)=="circadian")
- {
- if (tracks==0)
- {
- trackInfo[[trackCount]] <- makeGenericArray(intensity=as.matrix(effects[2,sta:end]),probeStart=Gloc[sta:end],dp=DisplayPars(color="black",ylim=c(min(-1.3,range(effects[2,sta:end])[1]),max(1.3,range(effects[2,sta:end])[2]))+c(-0.2,0.2),pointSize=.3,pch=1,lwd=1,type="line"))
- if (getWfmMethod(object)=="twoGroup")
- {
- trackInfo[[trackCount]] <- makeGenericArray(intensity=as.matrix(effects[2,sta:end]),probeStart=Gloc[sta:end],dp=DisplayPars(color="black",ylim=c(min(-1.3,range(effects[2,sta:end])[1]),max(1.3,range(effects[2,sta:end])[2])+c(-0.2,0.2)),pointSize=.3,pch=1,lwd=1,type="line"))
- names(trackInfo)[trackCount] <- "FC"
- } else
- {
- trackInfo[[trackCount]] <- makeGenericArray(intensity=as.matrix(sqrt(effects[2,sta:end]^2+effects[3,sta:end]^2)),probeStart=Gloc[sta:end],dp=DisplayPars(color="black",ylim=c(min(-1.3,range(sqrt(effects[2,sta:end]^2+effects[3,sta:end]^2))[1]),max(1.2,range(sqrt(effects[2,sta:end]^2+effects[3,sta:end]^2))))+c(-0.2,0.2),pointSize=.3,pch=1,lwd=1,type="line"))
- names(trackInfo)[trackCount] <- "Ampl"
- }
- overlayInfo[[trackCount]] <- makeNewTranscriptRectangleOverlay(sigRegions=as.matrix(data.frame(start(regions[[2]]),end(regions[[2]]))),location=Gloc,start=sta,end=end,region=c(trackCount,trackCount),dp=DisplayPars(color="darkgrey",alpha=.1))
- hlp <- which(!unlist(lapply(overlayInfo,is.null)))
- if (two.strand==TRUE)
- {
- hlpAnno <- hlp[1:2]
- } else
- {
- hlpAnno <- hlp[1]
- }
- hlpAnnoIndex <- unlist(lapply(hlpAnno,function(x) length(overlayInfo[[x]]@xpos)>0))
- hlpEff <- hlp[!(hlp%in%hlpAnno)]
- hlpEffIndex <- unlist(lapply(hlpEff,function(x) length(overlayInfo[[x]]@start)>0))
- hlpIndex <- c(hlpAnnoIndex,hlpEffIndex)
- overlayInfoCorr <- lapply(hlp,function(x) overlayInfo[[x]])
- gdPlot(trackInfo,minBase=minBase,maxBase=maxBase,overlay=overlayInfoCorr[hlpIndex])
- }
- }
- if (getWfmMethod(object)=="compareGroupsTime" | getWfmMethod(object)=="compareGroupsFactor")
- {
- if (length(tracks)==1 & tracks[1]==0)
- {
-
- effectsToPlot <- rep(0,noGroups)
- effectId <- c(1,3,6,10,15)
- effectNames <- c("2-1","3-2","4-3","5-4","6-5")
- effectNo <- 1
- while (effectNo<length(effectsToPlot))
- {
- trackInfo[[trackCount]] <- makeGenericArray(intensity=as.matrix(effects[effectId[effectNo]+1,sta:end]),probeStart=Gloc[sta:end],dp=DisplayPars(color="black",ylim=c(min(-1.3,range(effects[effectId[effectNo]+1,sta:end])[1]),max(1.3,range(effects[effectId[effectNo]+1,sta:end])[2]))+c(-0.2,0.2),pointSize=.3,pch=1,lwd=1,type="line"))
- names(trackInfo)[trackCount] <- effectNames[effectNo]
- overlayInfo[[trackCount]] <- makeNewTranscriptRectangleOverlay(sigRegions=as.matrix(data.frame(start(regions[[effectId[effectNo]+1]]),end(regions[[effectId[effectNo]+1]]))),location=Gloc,start=sta,end=end,region=c(trackCount,trackCount),dp=DisplayPars(color="darkgrey",alpha=.1))
- effectNo <- effectNo + 1
- trackCount <- trackCount + 1
- }
- trackInfo[[trackCount]] <- makeGenericArray(intensity=as.matrix(effects[effectId[length(effectsToPlot)-2]+2,sta:end]),probeStart=Gloc[sta:end],dp=DisplayPars(color="black",ylim=c(min(-1.3,range(effects[effectId[length(effectsToPlot)-2]+2,sta:end])[1]),max(1.3,range(effects[effectId[length(effectsToPlot)-2]+2,sta:end])[2]))+c(-0.2,0.2),pointSize=.3,pch=1,lwd=1,type="line"))
- names(trackInfo)[trackCount] <- "Last-First"
- overlayInfo[[trackCount]] <- makeNewTranscriptRectangleOverlay(sigRegions=as.matrix(data.frame(start(regions[[effectId[length(effectsToPlot)-2]+2]]),end(regions[[effectId[length(effectsToPlot)-2]+2]]))),location=Gloc,start=sta,end=end,region=c(trackCount,trackCount),dp=DisplayPars(color="darkgrey",alpha=.1))
- hlp <- which(!unlist(lapply(overlayInfo,is.null)))
- if (two.strand==TRUE)
- {
- hlpAnno <- hlp[1:2]
- } else
- {
- hlpAnno <- hlp[1]
- }
- hlpAnnoIndex <- unlist(lapply(hlpAnno,function(x) length(overlayInfo[[x]]@xpos)>0))
- hlpEff <- hlp[!(hlp%in%hlpAnno)]
- hlpEffIndex <- unlist(lapply(hlpEff,function(x) length(overlayInfo[[x]]@start)>0))
- hlpIndex <- c(hlpAnnoIndex,hlpEffIndex)
- overlayInfoCorr <- lapply(hlp,function(x) overlayInfo[[x]])
- gdPlot(trackInfo,minBase=minBase,maxBase=maxBase,overlay=overlayInfoCorr[hlpIndex])
-
- }
- if (length(tracks)>1 | tracks[1]!=0)
- {
- if (length(tracks)>8)
- {
- stop("Too many tracks specified")
- }
- groupNames <- sapply(1:noGroups,function(x) paste("Gr",x,sep=""))
- #groupNames <- sapply(8:13,function(x) paste("D",x,sep=""))
- firstId <- rep(2:noGroups,1:(noGroups-1))
- lastId <- sequence(1:(noGroups-1))
- effectNames <- sapply(1:(noGroups*(noGroups-1)/2),function(x) paste(groupNames[firstId[x]],"-",groupNames[lastId[x]],sep=""))
- for (i in tracks)
- {
- #trackInfo[[trackCount]] <- makeGenericArray(intensity=as.matrix(effects[i+1,sta:end]),probeStart=Gloc[sta:end],dp=DisplayPars(color="black",ylim=c(min(-1.3,range(effects[i+1,sta:end])[1]),max(1.3,range(effects[i+1,sta:end])[2]))+c(-0.2,0.2),pointSize=.3,pch=1,lwd=1,type="line"))
- trackInfo[[trackCount]] <- makeGenericArray(intensity=as.matrix(effects[i+1,sta:end]),probeStart=Gloc[sta:end],dp=DisplayPars(color="black",ylim=range(effects[,sta:end])+c(-0.2,0.2),pointSize=.3,pch=1,lwd=1,type="line"))
- names(trackInfo)[trackCount] <- effectNames[i]
- overlayInfo[[trackCount]] <- makeNewTranscriptRectangleOverlay(sigRegions=as.matrix(data.frame(start(regions[[i+1]]),end(regions[[i+1]]))),location=Gloc,start=sta,end=end,region=c(trackCount,trackCount),dp=DisplayPars(color="darkgrey",alpha=.1))
- trackCount <- trackCount + 1
- }
- hlp <- which(!unlist(lapply(overlayInfo,is.null)))
- if (two.strand==TRUE)
- {
- hlpAnno <- hlp[1:2]
- } else
- {
- hlpAnno <- hlp[1]
- }
- hlpAnnoIndex <- unlist(lapply(hlpAnno,function(x) length(overlayInfo[[x]]@xpos)>0))
- hlpEff <- hlp[!(hlp%in%hlpAnno)]
- hlpEffIndex <- unlist(lapply(hlpEff,function(x) length(overlayInfo[[x]]@start)>0))
- hlpIndex <- c(hlpAnnoIndex,hlpEffIndex)
- overlayInfoCorr <- lapply(hlp,function(x) overlayInfo[[x]])
- gdPlot(trackInfo,minBase=minBase,maxBase=maxBase,overlay=overlayInfoCorr[hlpIndex])
- }
-
- }
- if (getWfmMethod(object)=="meansByGroupTime" | getWfmMethod(object)=="meansByGroupFactor")
- {
- if (length(tracks)==1 & tracks[1]==0)
- {
- effectNames <- sapply(1:noGroups,function(x) paste("Gr",x,sep=""))
- for (i in 1:noGroups)
- {
- trackInfo[[trackCount]] <- makeGenericArray(intensity=as.matrix(effects[i,sta:end]),probeStart=Gloc[sta:end],dp=DisplayPars(color="black",ylim=c(min(-1.3,range(effects[i,sta:end])[1]),max(1.3,range(effects[i,sta:end])[2]))+c(-0.2,0.2),pointSize=.3,pch=1,lwd=1,type="line"))
- names(trackInfo)[trackCount] <- effectNames[i]
- overlayInfo[[trackCount]] <- makeNewTranscriptRectangleOverlay(sigRegions=as.matrix(data.frame(start(regions[[i]]),end(regions[[i]]))),location=Gloc,start=sta,end=end,region=c(trackCount,trackCount),dp=DisplayPars(color="darkgrey",alpha=.1))
- trackCount <- trackCount + 1
- }
- hlp <- which(!unlist(lapply(overlayInfo,is.null)))
- if (two.strand==TRUE)
- {
- hlpAnno <- hlp[1:2]
- } else
- {
- hlpAnno <- hlp[1]
- }
- hlpAnnoIndex <- unlist(lapply(hlpAnno,function(x) length(overlayInfo[[x]]@xpos)>0))
- hlpEff <- hlp[!(hlp%in%hlpAnno)]
- hlpEffIndex <- unlist(lapply(hlpEff,function(x) length(overlayInfo[[x]]@start)>0))
- hlpIndex <- c(hlpAnnoIndex,hlpEffIndex)
- overlayInfoCorr <- lapply(hlp,function(x) overlayInfo[[x]])
- gdPlot(trackInfo,minBase=minBase,maxBase=maxBase,overlay=overlayInfoCorr[hlpIndex])
- }
- if (length(tracks)>1 | tracks[1]!=0)
- {
- if (length(tracks)>8)
- {
- stop("Too many tracks specified")
- }
- #effectNames <- sapply(1:noGroups,function(x) paste("Gr",x,sep=""))
- effectNames <- paste("day",8:13,sep="")
- for (i in tracks)
- {
- trackInfo[[trackCount]] <- makeGenericArray(intensity=as.matrix(effects[i,sta:end]),probeStart=Gloc[sta:end],dp=DisplayPars(color="black",ylim=c(min(-1.3,range(effects[i,sta:end])[1]),max(1.3,range(effects[i,sta:end])[2]))+c(-0.2,0.2),pointSize=.3,pch=1,lwd=1,type="line"))
- #trackInfo[[trackCount]] <- makeGenericArray(intensity=as.matrix(effects[i,sta:end]),probeStart=Gloc[sta:end],dp=DisplayPars(color="black",ylim=range(effects[,sta:end])+c(-0.2,0.2),pointSize=.3,pch=1,lwd=1,type="line"))
- names(trackInfo)[trackCount] <- effectNames[i]
- overlayInfo[[trackCount]] <- makeNewTranscriptRectangleOverlay(sigRegions=as.matrix(data.frame(start(regions[[i]]),end(regions[[i]]))),location=Gloc,start=sta,end=end,region=c(trackCount,trackCount),dp=DisplayPars(color="darkgrey",alpha=.1))
- trackCount <- trackCount + 1
- }
- hlp <- which(!unlist(lapply(overlayInfo,is.null)))
- if (two.strand==TRUE)
- {
- hlpAnno <- hlp[1:2]
- } else
- {
- hlpAnno <- hlp[1]
- }
- hlpAnnoIndex <- unlist(lapply(hlpAnno,function(x) length(overlayInfo[[x]]@xpos)>0))
- hlpEff <- hlp[!(hlp%in%hlpAnno)]
- hlpEffIndex <- unlist(lapply(hlpEff,function(x) length(overlayInfo[[x]]@start)>0))
- hlpIndex <- c(hlpAnnoIndex,hlpEffIndex)
- overlayInfoCorr <- lapply(hlp,function(x) overlayInfo[[x]])
- gdPlot(trackInfo,minBase=minBase,maxBase=maxBase,overlay=overlayInfoCorr[hlpIndex])
- }
- }
- if (getWfmMethod(object)=="effectsTime")
- {
- fct <- 0.2390457
- if (tracks==0 | (length(tracks)==2 & tracks[1]==1 & tracks[2]==2) | tracks==1 | tracks==2)
- {
- if (tracks==0)
- {
- tracks <- c(1,2)
- }
- effectNames <- c("Linear","Quadratic")
- for (i in tracks)
- {
- #trackInfo[[trackCount]] <- makeGenericArray(intensity=as.matrix(effects[i+1,sta:end]),probeStart=Gloc[sta:end],dp=DisplayPars(color="black",ylim=range(effects[i+1,sta:end])+c(-0.5,0.5),pointSize=.3,pch=1,lwd=1,type="line"))
- trackInfo[[trackCount]] <- makeGenericArray(intensity=as.matrix(effects[i+1,sta:end]*fct),probeStart=Gloc[sta:end],dp=DisplayPars(color="black",ylim=range(effects[i+1,sta:end]*fct)+c(-0.2,0.2),pointSize=.3,pch=1,lwd=1,type="line"))
- names(trackInfo)[trackCount] <- effectNames[i]
- overlayInfo[[trackCount]] <- makeNewTranscriptRectangleOverlay(sigRegions=as.matrix(data.frame(start(regions[[i+1]]),end(regions[[i+1]]))),location=Gloc,start=sta,end=end,region=c(trackCount,trackCount),dp=DisplayPars(color="darkgrey",alpha=.1))
- trackCount <- trackCount + 1
- }
- hlp <- which(!unlist(lapply(overlayInfo,is.null)))
- if (two.strand==TRUE)
- {
- hlpAnno <- hlp[1:2]
- } else
- {
- hlpAnno <- hlp[1]
- }
- hlpAnnoIndex <- unlist(lapply(hlpAnno,function(x) length(overlayInfo[[x]]@xpos)>0))
- hlpEff <- hlp[!(hlp%in%hlpAnno)]
- hlpEffIndex <- unlist(lapply(hlpEff,function(x) length(overlayInfo[[x]]@start)>0))
- hlpIndex <- c(hlpAnnoIndex,hlpEffIndex)
- overlayInfoCorr <- lapply(hlp,function(x) overlayInfo[[x]])
- gdPlot(trackInfo,minBase=minBase,maxBase=maxBase,overlay=overlayInfoCorr[hlpIndex])
- } else
- {
- stop("tracks argument has been misspecified")
- }
- }
-})
-
-
-setMethod("getNonAnnotatedRegions",signature("Wfm"),function(object,annoFile)
-{
- #Gloc <- getProbePosition(object)
- strand <- getStrand(object)
- chromosome <- getChromosome(object)
- regions <- getGenomicRegions(object)
- if (strand=="forward")
- {
- strandAlt <- "+"
- strandOpp <- "-"
- } else
- {
- strandAlt <- "-"
- strandOpp <- "+"
- }
- cat("get annotated regions...\n")
- annoExons <- annoFile[(annoFile$strand==strandAlt)&(annoFile$chromosome==chromosome)&((annoFile$feature=="exon")|(annoFile$feature=="pseudogenic_exon")),c("chromosome","strand","feature","ID","start","end")]
- annoExonsOpp <- annoFile[(annoFile$strand==strandOpp)&(annoFile$chromosome==chromosome)&((annoFile$feature=="exon")|(annoFile$feature=="pseudogenic_exon")),c("chromosome","strand","feature","ID","start","end")]
- regGlocNoAnnoSense <- list()
- regGlocNoAnnoBoth <- list()
- nList <- length(regions)
- cat("find overlaps with detected regions...\n")
- for (j in 1:nList)
- {
- regGlocIR <- regions[[j]]
- annoExonsIR <- IRanges(start=annoExons$start,end=annoExons$end)
- annoExonsOppIR <- IRanges(start=annoExonsOpp$start,end=annoExonsOpp$end)
- overSense <- findOverlaps(regGlocIR,annoExonsIR)
- overOpp <- findOverlaps(regGlocIR,annoExonsOppIR)
- noAnnoSenseId <- which(!(1:length(regGlocIR) %in% matchMatrix(overSense)[,1]))
- noAnnoOppId <- which(!(1:length(regGlocIR) %in% matchMatrix(overOpp)[,1]))
- noAnnoBothId <- which((1:length(regGlocIR) %in% noAnnoSenseId) & (1:length(regGlocIR) %in% noAnnoOppId))
- regGlocNoAnnoSense[[j]] <- regGlocIR[noAnnoSenseId]
- regGlocNoAnnoBoth[[j]] <- regGlocIR[noAnnoBothId]
- }
- noAnnoSenseIR <- regGlocNoAnnoSense[[1]]
- noAnnoBothIR <- regGlocNoAnnoBoth[[1]]
- for (i in 2:nList)
- {
- noAnnoSenseIRi <- regGlocNoAnnoSense[[i]]
- noAnnoSenseIR <- c(noAnnoSenseIR,noAnnoSenseIRi)
- noAnnoBothIRi <- regGlocNoAnnoBoth[[i]]
- noAnnoBothIR <- c(noAnnoBothIR,noAnnoBothIRi)
- }
- noAnnoSenseIRAll <- reduce(noAnnoSenseIR)
- noAnnoBothIRAll <- reduce(noAnnoBothIR)
- ## TO DO: include option to give maximum expression / FC per region
- out <- NULL
- out$noAnnoBoth <- GRanges(seqnames=Rle(rep(chromosome,length(noAnnoBothIRAll))),strand=Rle(rep(strandAlt,length(noAnnoBothIRAll))),ranges=noAnnoBothIRAll)
- out$noAnnoSense <- GRanges(seqnames=Rle(rep(chromosome,length(noAnnoSenseIRAll))),strand=Rle(rep(strandAlt,length(noAnnoSenseIRAll))),ranges=noAnnoSenseIRAll)
- return(out)
-})
-
-
-setMethod("getSigGenes",signature("Wfm"),function(object,annoFile)
-{
- #Gloc <- getProbePosition(object)
- strand <- getStrand(object)
- chromosome <- getChromosome(object)
- regions <- getGenomicRegions(object)
- annoFile$strand[annoFile$strand=="forward"] <- "+"
- annoFile$strand[annoFile$strand=="reverse"] <- "-"
- annoFile$strand[!(annoFile$strand %in% c("+","-"))] <- "*"
- annoFileGR <- GRanges(seqnames=Rle(annoFile$chromosome),ranges=IRanges(start=annoFile$start,end=annoFile$end),strand=Rle(annoFile$strand),feature=annoFile$feature,id=annoFile$ID)
- if (strand=="forward")
- {
- strandAlt <- "+"
- strandOpp <- "-"
- } else
- {
- strandAlt <- "-"
- strandOpp <- "+"
- }
- annoChrGR <- annoFileGR[seqnames(annoFileGR)==chromosome]
- geneId <- which(values(annoChrGR)$feature=="gene" | values(annoChrGR)$feature=="transposable_element_gene")
- annoChrGeneGR <- annoChrGR[geneId]
- #annoChrGeneStrandGR <- annoChrGeneGR[strand(annoChrGeneGR)==strandAlt]
- #annoChrGeneStrandOppGR <- annoChrGeneGR[strand(annoChrGeneGR)==strandOpp]
- cat("find overlaps with detected regions...\n")
- nList <- length(regions)
- annoOver <- GRangesList()
- for (j in 1:nList)
- {
- regGlocIR <- regions[[j]]
- regGlocGR <- GRanges(seqnames=rep(chromosome,length(regGlocIR)),ranges=regGlocIR,strand=rep("*",length(regGlocIR)),effectNo=rep(j,length(regGlocIR)))
- overL <- findOverlaps(regGlocGR,annoChrGeneGR)
- regOver <- regGlocGR[queryHits(overL)]
- annoOverj <- annoChrGeneGR[subjectHits(overL)]
- overInt <- pintersect(regOver,annoOverj)
- values(annoOverj)$regNo <- queryHits(overL)
- values(annoOverj)$percOverGene <- width(overInt)/width(annoOverj)*100
- values(annoOverj)$percOverReg <- width(overInt)/width(regOver)*100
- totPercOverGeneHlp <- rep(0,max(subjectHits(overL)))
- totPercOverGeneHlp2 <- tapply(values(annoOverj)$percOverGene,subjectHits(overL),sum)
- totPercOverGeneHlp[as.numeric(names(totPercOverGeneHlp2))] <- totPercOverGeneHlp2
- values(annoOverj)$totPercOverGene <- totPercOverGeneHlp[subjectHits(overL)]
- annoOverj <- GRangesList(annoOverj)
- annoOver <- c(annoOver,annoOverj)
- }
- return(annoOver)
-})
-
-
-
-
-
-
-
Added: pkg/R/methods-WfmFit.R
===================================================================
--- pkg/R/methods-WfmFit.R (rev 0)
+++ pkg/R/methods-WfmFit.R 2012-02-16 07:00:30 UTC (rev 19)
@@ -0,0 +1,431 @@
+setMethod("getProbePosition",signature("WfmFit"),function(object)
+{
+ return(object at probePosition)
+})
+
+setMethod("getNoProbes",signature("WfmFit"),function(object)
+{
+ return(length(object at probePosition))
+})
+
+setMethod("getBetaMAP",signature("WfmFit"),function(object)
+{
+ return(object at betaMAP)
+})
+
+setMethod("getVarBetaMAP",signature("WfmFit"),function(object)
+{
+ return(object at varbetaMAP)
+})
+
+setMethod("getSmoothPar",signature("WfmFit"),function(object)
+{
+ return(object at smoothPar)
+})
+
+setMethod("getVarEps",signature("WfmFit"),function(object)
+{
+ return(object at varEps)
+})
+
+
+setMethod("getGenomeInfo",signature("WfmFit"),function(object)
+{
+ return(object at genome.info)
+})
+
+setMethod("getChromosome",signature("WfmFit"),function(object)
+{
+ return(object at genome.info@chromosome)
+})
+
+setMethod("getStrand",signature("WfmFit"),function(object)
+{
+ return(object at genome.info@strand)
+})
+
+setMethod("getMinPos",signature("WfmFit"),function(object)
+{
+ return(object at genome.info@minPos)
+})
+
+setMethod("getMaxPos",signature("WfmFit"),function(object)
+{
+ return(object at genome.info@maxPos)
+})
+
+setMethod("getNoLevels",signature("WfmFit"),function(object)
+{
+ return(object at n.levels)
+})
+
+
+setMethod("getDesignMatrix",signature("WfmFit"),function(object)
+{
+ return(object at design.matrix)
+})
+
+setMethod("getPhenoInfo",signature("WfmFit"),function(object)
+{
+ return(object at phenoData)
+})
+
+setMethod("getDataOrigSpace",signature("WfmFit"),function(object)
+{
+ return(object at dataOrigSpace)
+})
+
+setMethod("getDataWaveletSpace",signature("WfmFit"),function(object)
+{
+ return(object at dataWaveletSpace)
+})
+
+setMethod("getWaveletFilter",signature("WfmFit"),function(object)
+{
+ return(object at wave.filt)
+})
+
+setMethod("getKj",signature("WfmFit"),function(object)
+{
+ return(object at Kj)
+})
+
+setMethod("getPrior",signature("WfmFit"),function(object)
+{
+ return(object at prior)
+})
+
+
+setMethod("getF",signature("WfmFit"),function(object)
+{
+ return(object at F)
+})
+
+setMethod("getP",signature("WfmFit"),function(object)
+{
+ return(object at P)
+})
+
+setMethod("getZ",signature("WfmFit"),function(object)
+{
+ return(object at Z)
+})
+
+
+setMethod("wfm.inference",signature("WfmFit"),function(object,contrast.matrix=NULL,contrasts=c("compare","means","effects","overallMean"),delta=NULL,two.sided=NULL,minRunPos=90,minRunProbe=1,alpha=0.05,nsim=1000,rescale=NULL)
+{
+
+ sigProbes <- list()
+ regions <- list()
+ GlocRegions <- list()
+ givenDelta <- delta
+ noGroups <- object at noGroups
+ replics<- object at replics
+ F<-object at F
+ varF<-object at varF
+
+ P<-ncol(F) ## so P does not need to be stored!
+
+ Xsel <- cumsum(replics)-replics+1
+ X <- object at design.matrix
+ Xdes <- X[Xsel,]
+
+ Z<-object at Z
+
+ if (!is.null(contrast.matrix)) {
+ ## Given contrast matrix (CustomFit)
+ # Further implementation needed
+ warning("Custom Inference Procedure Not Implemented yet!")
+ }
+ else if (contrasts=="compare") {
+ if (inherits(object,"WfmFitFactor") | inherits(object,"WfmFitTime") | inherits(object,"WfmFitCircadian" | inherits(object,"WfmFitCustom"))) {
+ #q <- noGroups*(noGroups-1)/2
+ q <- noGroups*(noGroups-1)/2
+ contr <- makeContrasts(contrasts=contrasts,nlevels=noGroups);
+ noBetas <- noGroups
+ if (is.null(rescale))
+ {
+ rescale <- contr%*%Xdes
+ rescale <- rbind(c(mean(Xdes[,1]),rep(0,noGroups-1)),rescale)
+ }
+ ## effects!
+ eff <- rescale%*%solve(t(Z)%*%X)%*%F
+ varEff <- (rescale%*%solve(t(Z)%*%X))^2%*%varF
+
+ if (length(alpha)==1)
+ {
+ alpha <- rep(alpha,q+1)
+ }
+ FDR <- matrix(0,nrow=q+1,ncol=P)
+ CI <- rep(0,P*(q+1)*2)
+ dim(CI) <- c(q+1,2,P)
+ if (is.null(givenDelta))
+ {
+ delta <- c(median(getDataOrigSpace(object)),rep(log(1.1,2),q))
+ } else if (length(givenDelta)==1)
+ {
+ delta <- rep(delta,q+1)
+ } else if ((length(givenDelta)==2) & givenDelta[1]=="median")
+ {
+ delta <- rep(0,q+1)
+ delta[1] <- median(getDataOrigSpace(object))
+ delta[2:(q+1)] <- rep(as.numeric(givenDelta[2]),q)
+ } else if ((length(givenDelta)==q+1) & givenDelta[1]=="median")
+ {
+ delta <- rep(0,q+1)
+ delta[1] <- median(getDataOrigSpace(object))
+ delta[2:(q+1)] <- as.numeric(givenDelta[2:(q+1)])
+ }
+
+ if (is.null(two.sided))
+ {
+ two.sided <- c(0,rep(1,q))
+ }
+ for (i in 1:(q+1))
+ {
+ if (two.sided[i]==1)
+ {
+ #FDR[i,] <- pnorm(delta[i],abs(eff[i,]),sqrt(varEff[i,]))
+ FDRUp <- pnorm(delta[i],eff[i,],sqrt(varEff[i,]))
+ FDRDown <- 1-pnorm(-delta[i],eff[i,],sqrt(varEff[i,]))
+ FDR[i,] <- pmin(FDRUp,FDRDown)
+ }
+ if (two.sided[i]==0)
+ {
+ FDR[i,] <- pnorm(delta[i],eff[i,],sqrt(varEff[i,]))
+ }
+ CI[i,1,] <- qnorm(alpha/2,eff[i,],sqrt(varEff[i,]))
+ CI[i,2,] <- qnorm(1-alpha/2,eff[i,],sqrt(varEff[i,]))
+ }
+ } else {
+ stop ("Analysis 'compare' not specified for this design!")
+ }
+ }
+ else if (contrasts=="effects") {
+ if (inherits(object,"WfmFitTime")) {
+ #q <- noGroups-1
+ q <- noGroups-1
+ noBetas <- noGroups
+ rescale <- diag(noBetas)
+ if (length(alpha)==1)
+ {
+ alpha <- rep(alpha,q+1)
+ }
+
+
+ ## effects!
+ eff <- rescale%*%solve(t(Z)%*%X)%*%F
+ varEff <- (rescale%*%solve(t(Z)%*%X))^2%*%varF
+ FDR <- matrix(0,nrow=q+1,ncol=P)
+ CI <- rep(0,P*(q+1)*2)
+ dim(CI) <- c(q+1,2,P)
+ if (is.null(givenDelta))
+ {
+ delta <- c(median(getDataOrigSpace(object)),rep(log(1.1,2),q))
+ } else if (length(givenDelta)==1)
+ {
+ delta <- rep(delta,q+1)
+ } else if ((length(givenDelta)==2) & givenDelta[1]=="median")
+ {
+ delta <- rep(0,q+1)
+ delta[1] <- median(getDataOrigSpace(object))
+ delta[2:(q+1)] <- rep(as.numeric(givenDelta[2]),q)
+ } else if ((length(givenDelta)==q+1) & givenDelta[1]=="median")
+ {
+ delta <- rep(0,q+1)
+ delta[1] <- median(getDataOrigSpace(object))
+ delta[2:(q+1)] <- as.numeric(givenDelta[2:(q+1)])
+ }
+
+ if (is.null(two.sided))
+ {
+ two.sided <- c(0,rep(1,q))
+ }
+ for (i in 1:(q+1))
+ {
+ if (two.sided[i]==1)
+ {
+ #FDR[i,] <- pnorm(delta[i],abs(eff[i,]),sqrt(varEff[i,]))
+ FDRUp <- pnorm(delta[i],eff[i,],sqrt(varEff[i,]))
+ FDRDown <- 1-pnorm(-delta[i],eff[i,],sqrt(varEff[i,]))
+ FDR[i,] <- pmin(FDRUp,FDRDown)
+ }
+ if (two.sided[i]==0)
+ {
+ FDR[i,] <- pnorm(delta[i],eff[i,],sqrt(varEff[i,]))
+ }
+ CI[i,1,] <- qnorm(alpha/2,eff[i,],sqrt(varEff[i,]))
+ CI[i,2,] <- qnorm(1-alpha/2,eff[i,],sqrt(varEff[i,]))
+ }
+ }
+ else if (inherits(object,"WfmFitCircadian")) {
+ #q <- 1
+ noBetas <- 3
+ q <- 1 # check
+ rescale <- diag(noBetas)
+ if (length(alpha)==1)
+ {
+ alpha <- rep(alpha,q+1)
+ }
+
+
+ ## effects!
+ eff <- rescale%*%solve(t(Z)%*%X)%*%F
+ varEff <- (rescale%*%solve(t(Z)%*%X))^2%*%varF
+ FDR <- matrix(0,nrow=q+1,ncol=P)
+ CI <- rep(0,P*(q+1)*2)
+ dim(CI) <- c(q+1,2,P)
+ if (is.null(givenDelta))
+ {
+ delta <- c(median(getDataOrigSpace(object)),rep(log(1.1,2),q))
+ } else if (length(givenDelta)==1)
+ {
+ delta <- rep(delta,q+1)
+ } else if ((length(givenDelta)==2) & givenDelta[1]=="median")
+ {
+ delta <- rep(0,q+1)
+ delta[1] <- median(getDataOrigSpace(object))
+ delta[2:(q+1)] <- rep(as.numeric(givenDelta[2]),q)
+ } else if ((length(givenDelta)==q+1) & givenDelta[1]=="median")
+ {
+ delta <- rep(0,q+1)
+ delta[1] <- median(getDataOrigSpace(object))
+ delta[2:(q+1)] <- as.numeric(givenDelta[2:(q+1)])
+ }
+
+ ## fix me: use of two.sided
+ if (is.null(two.sided))
+ {
+ two.sided <- c(0,rep(1,q))
+ }
+ FDRUp <- pnorm(delta[1],eff[1,],sqrt(varEff[1,]))
+ FDRDown <- 1-pnorm(-delta[1],eff[1,],sqrt(varEff[1,]))
+ FDR[1,] <- pmin(FDRUp,FDRDown)
+ CI[1,1,] <- qnorm(alpha/2,eff[1,],sqrt(varEff[1,]))
+ CI[1,2,] <- qnorm(1-alpha/2,eff[1,],sqrt(varEff[1,]))
+ FDR[2,] <- rep(0,P)
+
+ for (k in 1:nsim)
+ {
+ effSinSampl <- rnorm(P,eff[2,],sqrt(varEff[2,]))
+ effCosSampl <- rnorm(P,eff[3,],sqrt(varEff[3,]))
+ amplSampl <- sqrt(effSinSampl^2 + effCosSampl^2)
+ FDR[2,] <- FDR[2,] + (amplSampl < delta[2])
+ if (k%%100==0) cat(k," ")
+ # calculate CIs: to do
+ }
+ FDR[2,] <- FDR[2,]/nsim
+ }
+
+ else {
+ stop ("Analysis 'effects' not specified for this design!")
+ }
+ }
+ else if (contrasts=="means") {
+ if (inherits(object,"WfmFitFactor") | inherits(object,"WfmFitTime") | inherits(object,"WfmFitCircadian") | inherits(object,"WfmFitCustom")) {
+ #q <- noGroups-1
+ q <- noGroups-1
+ noBetas <- noGroups
+ rescale <- Xdes
+ if (length(alpha)==1)
+ {
+ alpha <- rep(alpha,q+1)
+ }
+
+
+ ## effects!
+ eff <- rescale%*%solve(t(Z)%*%X)%*%F
+ varEff <- (rescale%*%solve(t(Z)%*%X))^2%*%varF
+ FDR <- matrix(0,nrow=q+1,ncol=P)
+ CI <- rep(0,P*(q+1)*2)
+ dim(CI) <- c(q+1,2,P)
+ if (is.null(givenDelta))
+ {
+ delta <- rep(median(getDataOrigSpace(object)),q+1)
+ } else if (length(givenDelta)==1)
+ {
+ delta <- rep(delta,q+1)
+ } else if ((length(givenDelta)==2) & givenDelta[1]=="median")
+ {
+ delta <- rep(0,q+1)
+ delta[1] <- median(getDataOrigSpace(object))
+ delta[2:(q+1)] <- rep(as.numeric(givenDelta[2]),q)
+ } else if ((length(givenDelta)==q+1) & givenDelta[1]=="median")
+ {
+ delta <- rep(0,q+1)
+ delta[1] <- median(getDataOrigSpace(object))
+ delta[2:(q+1)] <- as.numeric(givenDelta[2:(q+1)])
+ }
+
+ if (is.null(two.sided))
+ {
+ two.sided <- rep(1,q+1)
+ }
+ for (i in 1:(q+1))
+ {
+ FDR[i,] <- pnorm(delta[i],eff[i,],sqrt(varEff[i,]))
+ CI[i,1,] <- qnorm(alpha/2,eff[i,],sqrt(varEff[i,]))
+ CI[i,2,] <- qnorm(1-alpha/2,eff[i,],sqrt(varEff[i,]))
+ }
+ }
+ else {
+ stop ("Analysis 'means' not specified for this design!")
+ }
+ }
+ else if (contrasts == "overallMean") {
+ ## overallMean
+ # Further implementation needed
+ warning("Contrast 'overall mean' not yet implemented!")
+ }
+ else {
+ stop ("No contrast matrix of contrast statement specified!")
+ }
+ Gloc <- object at probePosition
+
+ for (i in 1:(q+1))
+ {
+ sigProbes[[i]] <- (1:P)[FDR[i,]<alpha[i]]
+ regions[[i]] <- cbind(sigProbes[[i]][c(TRUE,(sigProbes[[i]][2:length(sigProbes[[i]])]-sigProbes[[i]][1:(length(sigProbes[[i]])-1)])>1)],sigProbes[[i]][c((sigProbes[[i]][2:length(sigProbes[[i]])]-sigProbes[[i]][1:(length(sigProbes[[i]])-1)])>1,TRUE)])
+ if (length(regions[[i]])==0)
+ {
+ regions[[i]] <- matrix(,0,2)
+ }
+ if (length(regions[[i]])==2)
+ {
+ regions[[i]] <- matrix(regions[[i]],nrow=1)
+ }
+ regions[[i]] <- matrix(regions[[i]][(Gloc[regions[[i]][,2]]-Gloc[regions[[i]][,1]])>minRunPos,],ncol=2)
+ regions[[i]] <- matrix(regions[[i]][regions[[i]][,2]-regions[[i]][,1]>minRunProbe,],ncol=2)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/wavetiling -r 19
More information about the Wavetiling-commits
mailing list