[Wavetiling-commits] r42 - in pkg/waveTiling: . R inst inst/doc inst/scripts man src vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon May 21 14:54:45 CEST 2012
Author: kdbeuf
Date: 2012-05-21 14:54:45 +0200 (Mon, 21 May 2012)
New Revision: 42
Added:
pkg/waveTiling/CHANGELOG
pkg/waveTiling/DESCRIPTION
pkg/waveTiling/NAMESPACE
pkg/waveTiling/R/
pkg/waveTiling/R/allClasses.R
pkg/waveTiling/R/allGenerics.R
pkg/waveTiling/R/helperFunctions.R
pkg/waveTiling/R/initialize-methods.R
pkg/waveTiling/R/methods-MapFilterProbe.R
pkg/waveTiling/R/methods-WaveTilingFeatureSet.R
pkg/waveTiling/R/methods-WfmFit.R
pkg/waveTiling/R/methods-WfmInf.R
pkg/waveTiling/R/show-methods.R
pkg/waveTiling/inst/
pkg/waveTiling/inst/NEWS.Rd
pkg/waveTiling/inst/doc/
pkg/waveTiling/inst/doc/Makefile
pkg/waveTiling/inst/doc/waveTiling-vignette.Rnw
pkg/waveTiling/inst/doc/waveTiling-vignette.pdf
pkg/waveTiling/inst/scripts/
pkg/waveTiling/inst/scripts/Sweave.sty
pkg/waveTiling/inst/scripts/waveTiling-vignette.Rnw
pkg/waveTiling/inst/scripts/waveTiling-vignette.pdf
pkg/waveTiling/man/
pkg/waveTiling/man/GenomeInfo-class.Rd
pkg/waveTiling/man/MapFilterProbe-class.Rd
pkg/waveTiling/man/WaveTilingFeatureSet-class.Rd
pkg/waveTiling/man/WfmFit-class.Rd
pkg/waveTiling/man/WfmFitCircadian-class.Rd
pkg/waveTiling/man/WfmFitCustom-class.Rd
pkg/waveTiling/man/WfmFitFactor-class.Rd
pkg/waveTiling/man/WfmFitTime-class.Rd
pkg/waveTiling/man/WfmInf-class.Rd
pkg/waveTiling/man/WfmInfCompare-class.Rd
pkg/waveTiling/man/WfmInfCustom-class.Rd
pkg/waveTiling/man/WfmInfEffects-class.Rd
pkg/waveTiling/man/WfmInfMeans-class.Rd
pkg/waveTiling/man/WfmInfOverallMean-class.Rd
pkg/waveTiling/man/addPheno.Rd
pkg/waveTiling/man/bgCorrQn.Rd
pkg/waveTiling/man/cel2TilingFeatureSet.Rd
pkg/waveTiling/man/filterOverlap.Rd
pkg/waveTiling/man/getNonAnnotatedRegions.Rd
pkg/waveTiling/man/getSigGenes.Rd
pkg/waveTiling/man/makeContrasts.Rd
pkg/waveTiling/man/makeDesign.Rd
pkg/waveTiling/man/plotWfm.Rd
pkg/waveTiling/man/selectProbesFromFilterOverlap.Rd
pkg/waveTiling/man/selectProbesFromTilingFeatureSet.Rd
pkg/waveTiling/man/wfm.fit.Rd
pkg/waveTiling/man/wfm.inference.Rd
pkg/waveTiling/src/
pkg/waveTiling/src/waveTiling.c
pkg/waveTiling/vignettes/
pkg/waveTiling/vignettes/waveTiling-vignette.pdf
Log:
added data and annotation packages
Added: pkg/waveTiling/CHANGELOG
===================================================================
--- pkg/waveTiling/CHANGELOG (rev 0)
+++ pkg/waveTiling/CHANGELOG 2012-05-21 12:54:45 UTC (rev 42)
@@ -0,0 +1,3 @@
+2011-12-04 Kristof De Beuf <kristof.debeuf at ugent.be> - committed version 0.1.1
+
+* First version
Added: pkg/waveTiling/DESCRIPTION
===================================================================
--- pkg/waveTiling/DESCRIPTION (rev 0)
+++ pkg/waveTiling/DESCRIPTION 2012-05-21 12:54:45 UTC (rev 42)
@@ -0,0 +1,19 @@
+Package: waveTiling
+Version: 0.99.0
+Date: 2012-05-14
+License: GPL (>=2)
+Title: Wavelet-Based Models for Tiling Array Transcriptome Analysis
+Author: Kristof De Beuf <kristof.debeuf at UGent.be>, Peter Pipelers <peter.pipelers at ugent.be> and Lieven Clement <lieven.clement at gmail.com>
+Maintainer: Kristof De Beuf <kristof.debeuf at UGent.be>
+Depends: oligo, oligoClasses, Biobase, Biostrings, GenomeGraphs
+Imports: methods, affy, preprocessCore, GenomicRanges, waveslim, IRanges
+Suggests: BSgenome, BSgenome.Athaliana.TAIR.TAIR9, waveTilingData, pd.atdschip.tiling
+Description: This package is designed to conduct transcriptome analysis for tiling arrays based on fast wavelet-based functional models.
+Collate: allClasses.R allGenerics.R helperFunctions.R
+ initialize-methods.R methods-MapFilterProbe.R
+ methods-WaveTilingFeatureSet.R methods-WfmFit.R methods-WfmInf.R show-methods.R
+URL: https://r-forge.r-project.org/projects/wavetiling/
+LazyLoad: yes
+LazyData:
+biocViews: MicroArray, DifferentialExpression, TimeCourse, GeneExpression
+Packaged: 2011-12-05 20:18:51 UTC; kdebeuf
Added: pkg/waveTiling/NAMESPACE
===================================================================
--- pkg/waveTiling/NAMESPACE (rev 0)
+++ pkg/waveTiling/NAMESPACE 2012-05-21 12:54:45 UTC (rev 42)
@@ -0,0 +1,74 @@
+useDynLib("waveTiling")
+
+## Import Classes
+importClassesFrom(Biobase, AnnotatedDataFrame, AssayData, MIAxE,
+ Versions)
+
+importClassesFrom(methods, array, character, data.frame, list, matrix,
+ numeric, vector)
+
+importClassesFrom(oligoClasses, TilingFeatureSet)
+
+## Import Methods
+importMethodsFrom(affy, pm, pmindex)
+
+importMethodsFrom(Biobase, exprs, "exprs<-", pData, "phenoData<-")
+
+importMethodsFrom(Biostrings, complement, matchPDict, PDict,
+ reverseComplement, startIndex)
+
+importMethodsFrom(GenomicRanges, seqnames)
+
+importMethodsFrom(IRanges, as.matrix, cbind, findOverlaps, ifelse, intersect,
+ lapply, mad, Map, mean, median, ncol, nrow, order, paste,
+ pintersect, pmin,queryHits, rbind, reduce, Reduce, Rle,
+ rownames,"rownames<-", sapply, stack, start, subjectHits,
+ subseq, t, table, tapply, unique, unlist, values,
+ "values<-", which, width)
+
+importMethodsFrom(methods, initialize, show)
+
+importMethodsFrom(oligo, pmChr, pmPosition, pmSequence, pmStrand)
+
+## Import
+importFrom(affy, bg.adjust, list.celfiles)
+
+importFrom(preprocessCore, normalize.quantiles)
+
+#importFrom(GenomeGraphs, DisplayPars, gdPlot, makeAnnotationTrack,
+# makeGenericArray, makeGenomeAxis, makeRectangleOverlay,
+# makeTextOverlay)
+
+importFrom(GenomicRanges, GRanges, GRangesList)
+
+importFrom(IRanges, IRanges)
+
+importFrom(methods, "@<-", callNextMethod, new)
+
+importFrom(oligo, read.celfiles)
+
+importFrom(stats, contr.helmert, contr.treatment, lm, pnorm, qnorm,
+ rnorm)
+
+importFrom(waveslim, wave.filter)
+
+## Export Classes
+#exportClasses(WaveTilingFeatureSet, MapFilterProbe, GenomeInfo, WfmFit, WfmInf, WfmFitTime, WfmFitFactor, WfmFitCircadian, WfmFitCustom, WfmInfCompare, WfmInfEffects )
+exportClasses(WaveTilingFeatureSet, MapFilterProbe, GenomeInfo, WfmFit, WfmInf)
+
+## Export Methods
+exportMethods(show)
+
+
+## WaveTilingFeatureSet
+exportMethods(addPheno, getNoGroups, getGroupNames, getReplics, filterOverlap, selectProbesFromTilingFeatureSet,bgCorrQn,wfm.fit)
+
+## MapFilterProbe
+exportMethods(getFilteredIndices, getPosition, selectProbesFromFilterOverlap)
+
+## WfmFit and WfmInf
+exportMethods(getProbePosition, getNoProbes, getBetaWav, getVarBetaWav, getSmoothPar, getVarEps, getGenomeInfo, getChromosome, getStrand, getMinPos, getMaxPos, getNoLevels, getDesignMatrix, getPhenoInfo, getDataOrigSpace, getDataWaveletSpace, getWaveletFilter, getKj,getPrior, getAlpha, getDelta, getTwoSided, getSigProbes,getRegions, getGenomicRegions, getFDR, getF, getVarF, getEff,getVarEff, wfm.inference, getSigGenes, getNonAnnotatedRegions, plotWfm)
+
+## Other
+export(cel2TilingFeatureSet, makeContrasts, makeDesign)
+
Added: pkg/waveTiling/R/allClasses.R
===================================================================
--- pkg/waveTiling/R/allClasses.R (rev 0)
+++ pkg/waveTiling/R/allClasses.R 2012-05-21 12:54:45 UTC (rev 42)
@@ -0,0 +1,45 @@
+#MapFilterProbe
+setClass("MapFilterProbe",representation(filteredIndices="vector",chromosome="vector",position="vector",strand="vector"))
+
+#GenomeInfo
+setClass("GenomeInfo",representation(chromosome="vector",strand="character",minPos="numeric",maxPos="numeric"))
+
+#WaveTilingFeatureSet
+setClass("WaveTilingFeatureSet",contains="TilingFeatureSet")
+
+#WfmFit
+setClass(Class="WfmFit",
+ representation = representation (betaWav="matrix",varbetaWav="matrix",smoothPar="matrix",varEps="numeric",dataOrigSpace="matrix",dataWaveletSpace="matrix",design.matrix="matrix",phenoData="data.frame",genome.info="GenomeInfo",n.levels="numeric",probePosition="vector",wave.filt="character",Kj="numeric",prior="character",F="matrix",varF="matrix",P="numeric",Z="matrix",noGroups="numeric",replics="numeric")
+)
+
+### Factor
+setClass(Class="WfmFitFactor", contains="WfmFit")
+
+### Time
+setClass(Class="WfmFitTime", contains="WfmFit")
+
+## Circadian
+setClass(Class="WfmFitCircadian", contains="WfmFit")
+
+## Custom
+setClass(Class="WfmFitCustom", contains="WfmFit")
+
+#WfmInf
+setClass(Class="WfmInf",
+ representation = representation (alpha="numeric",delta="numeric",two.sided="numeric",sigProbes="list",regions="list",GlocRegions="list",FDR="matrix",CI="array",eff="matrix",varEff="matrix")
+)
+
+## Compare
+setClass(Class="WfmInfCompare", contains="WfmInf")
+
+## Effects
+setClass(Class="WfmInfEffects", contains="WfmInf")
+
+## Means
+setClass(Class="WfmInfMeans", contains="WfmInf")
+
+## Overall Mean
+setClass(Class="WfmInfOverallMean", contains="WfmInf")
+
+## Custom
+setClass(Class="WfmInfCustom", contains="WfmInf")
Added: pkg/waveTiling/R/allGenerics.R
===================================================================
--- pkg/waveTiling/R/allGenerics.R (rev 0)
+++ pkg/waveTiling/R/allGenerics.R 2012-05-21 12:54:45 UTC (rev 42)
@@ -0,0 +1,318 @@
+# temporary(?) pmStrand
+#
+# setGeneric("pmStrand",function(object)
+# {
+# standardGeneric("pmStrand")
+# }
+# )
+
+# method data extraction
+
+setGeneric("addPheno",function(object, noGroups, groupNames, replics, ...)
+{
+ standardGeneric("addPheno")
+}
+)
+
+
+setGeneric("filterOverlap",function(object, remap=TRUE, BSgenomeObject, chrId, strand=c("forward","reverse","both"), MM=FALSE, ...)
+{
+ standardGeneric("filterOverlap")
+}
+)
+
+setGeneric("bgCorrQn",function(object, useMapFilter=NULL)
+{
+ standardGeneric("bgCorrQn")
+}
+)
+
+setGeneric("selectProbesFromTilingFeatureSet",function(object, chromosome, strand=c("forward","reverse"), minPos, maxPos)
+{
+ standardGeneric("selectProbesFromTilingFeatureSet")
+}
+)
+
+setGeneric("selectProbesFromFilterOverlap",function(object, chromosome, strand=c("forward","reverse"), minPos=min(getPosition(object)), maxPos=max(getPosition(object)))
+{
+ standardGeneric("selectProbesFromFilterOverlap")
+}
+)
+
+#fit method
+# setGeneric("makeDesign",function(object, method=c("twoGroup","compareGroupsTime","compareGroupsFactor","circadian","twoFactors","meansByGroupTime","meansByGroupFactor","effectsTime"), factor.levels=NULL)
+# {
+# standardGeneric("makeDesign")
+# }
+# )
+
+setGeneric("wfm.fit",function(object, filter.overlap=NULL, design=c("time","circadian","group","factorial","custom"), n.levels, factor.levels=NULL, chromosome, strand, minPos, maxPos, design.matrix=NULL, var.eps=c("margLik","mad"), prior=c("normal","improper"), eqsmooth=FALSE, max.it=20, wave.filt="haar", skiplevels=NULL, trace=FALSE, save.obs=c("plot","regions","all"))
+{
+ standardGeneric("wfm.fit")
+}
+)
+
+setGeneric("wfm.inference",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) {
+ standardGeneric("wfm.inference")
+}
+)
+
+#accessors mapFilterProbe
+setGeneric("getFilteredIndices",function(object)
+{
+ standardGeneric("getFilteredIndices")
+}
+)
+
+setGeneric("getChromosome",function(object)
+{
+ standardGeneric("getChromosome")
+}
+)
+
+setGeneric("getPosition",function(object)
+{
+ standardGeneric("getPosition")
+}
+)
+
+setGeneric("getStrand",function(object)
+{
+ standardGeneric("getStrand")
+}
+)
+
+#accessors TilingFeatureSet
+setGeneric("getNoGroups",function(object)
+{
+ standardGeneric("getNoGroups")
+}
+)
+
+setGeneric("getGroupNames",function(object)
+{
+ standardGeneric("getGroupNames")
+}
+)
+
+setGeneric("getReplics",function(object)
+{
+ standardGeneric("getReplics")
+}
+)
+
+######
+#accessors WfmFit
+setGeneric("getProbePosition",function(object)
+{
+ standardGeneric("getProbePosition")
+}
+)
+
+setGeneric("getNoProbes",function(object)
+{
+ standardGeneric("getNoProbes")
+}
+)
+
+setGeneric("getBetaWav",function(object)
+{
+ standardGeneric("getBetaWav")
+}
+)
+
+setGeneric("getVarBetaWav",function(object)
+{
+ standardGeneric("getVarBetaWav")
+}
+)
+
+setGeneric("getSmoothPar",function(object)
+{
+ standardGeneric("getSmoothPar")
+}
+)
+
+setGeneric("getVarEps",function(object)
+{
+ standardGeneric("getVarEps")
+}
+)
+
+setGeneric("getGenomeInfo",function(object)
+{
+ standardGeneric("getGenomeInfo")
+}
+)
+
+setGeneric("getMinPos",function(object)
+{
+ standardGeneric("getMinPos")
+}
+)
+
+setGeneric("getMaxPos",function(object)
+{
+ standardGeneric("getMaxPos")
+}
+)
+
+setGeneric("getNoLevels",function(object)
+{
+ standardGeneric("getNoLevels")
+}
+)
+
+setGeneric("getWfmMethod",function(object)
+{
+ standardGeneric("getWfmMethod")
+}
+)
+setGeneric("getDesignMatrix",function(object)
+{
+ standardGeneric("getDesignMatrix")
+}
+)
+
+setGeneric("getWfmDesign",function(object)
+{
+ standardGeneric("getWfmDesign")
+}
+)
+setGeneric("getPhenoInfo",function(object)
+{
+ standardGeneric("getPhenoInfo")
+}
+)
+
+setGeneric("getDataOrigSpace",function(object)
+{
+ standardGeneric("getDataOrigSpace")
+}
+)
+setGeneric("getDataWaveletSpace",function(object)
+{
+ standardGeneric("getDataWaveletSpace")
+}
+)
+
+setGeneric("getWaveletFilter",function(object)
+{
+ standardGeneric("getWaveletFilter")
+}
+)
+setGeneric("getKj",function(object)
+{
+ standardGeneric("getKj")
+}
+)
+
+setGeneric("getPrior",function(object)
+{
+ standardGeneric("getPrior")
+}
+)
+
+setGeneric("getAlpha",function(object)
+{
+ standardGeneric("getAlpha")
+}
+)
+
+setGeneric("getDelta",function(object)
+{
+ standardGeneric("getDelta")
+}
+)
+
+setGeneric("getTwoSided",function(object)
+{
+ standardGeneric("getTwoSided")
+}
+)
+
+setGeneric("getSigProbes",function(object)
+{
+ standardGeneric("getSigProbes")
+}
+)
+
+setGeneric("getRegions",function(object)
+{
+ standardGeneric("getRegions")
+}
+)
+
+setGeneric("getGenomicRegions",function(object)
+{
+ standardGeneric("getGenomicRegions")
+}
+)
+
+setGeneric("getFDR",function(object)
+{
+ standardGeneric("getFDR")
+}
+)
+
+setGeneric("getCI",function(object)
+{
+ standardGeneric("getCI")
+}
+)
+
+setGeneric("getF",function(object)
+{
+ standardGeneric("getF")
+}
+)
+
+setGeneric("getVarF",function(object)
+{
+ standardGeneric("getVarF")
+}
+)
+
+setGeneric("getEff",function(object)
+{
+ standardGeneric("getEff")
+}
+)
+
+setGeneric("getVarEff",function(object)
+{
+ standardGeneric("getVarEff")
+}
+)
+
+setGeneric("getDesignMatrix",function(object)
+{
+ standardGeneric("getDesignMatrix")
+}
+)
+
+setGeneric("plotWfm",function(fit, inf, annoFile, minPos, maxPos, trackFeature="exon", overlayFeature=c("gene","transposable_element_gene"), two.strand=TRUE, plotData=TRUE, plotMean=TRUE, tracks=0)
+{
+ standardGeneric("plotWfm")
+}
+)
+
+setGeneric("getSigGenes",function(fit, inf, annoFile)
+{
+ standardGeneric("getSigGenes")
+}
+)
+
+setGeneric("getNonAnnotatedRegions",function(fit, inf, annoFile)
+{
+ standardGeneric("getNonAnnotatedRegions")
+}
+)
+
+#if (!isGeneric("plot")) {
+# setGeneric("plot",function(fit,inf,...){
+# standardGeneric("plot")
+# }
+#)
+#}
+
Added: pkg/waveTiling/R/helperFunctions.R
===================================================================
--- pkg/waveTiling/R/helperFunctions.R (rev 0)
+++ pkg/waveTiling/R/helperFunctions.R 2012-05-21 12:54:45 UTC (rev 42)
@@ -0,0 +1,376 @@
+
+## functions for data input
+
+cel2TilingFeatureSet <- function(dataPath,annotationPackage)
+{
+ expFiles <- paste(dataPath,list.celfiles(dataPath),sep="/")
+ tfs <- read.celfiles(expFiles,pkgname=annotationPackage)
+ return(tfs)
+}
+
+
+extractChrFromFasta <- function(char,fasPath)
+{
+ splt <- unlist(strsplit(char,fasPath))
+ hlp <- splt[seq(2,length(splt),2)]
+ splt2 <- unlist(strsplit(hlp,".fas"))
+ return(splt2)
+}
+
+pm2mm <- function(z)
+{
+ z <- complement(subseq(z,start=13,end=13))
+ return(z)
+}
+
+## functions for model fitting
+
+normVec <- function(vec)
+{
+ return(vec/sqrt(sum(vec^2)))
+}
+
+wave.transform <- function(x,n.levels=log(length(x),2),filt="haar")
+{
+ N <- length(x)
+ J <- n.levels
+ if (N/2^J != trunc(N/2^J)) stop("Sample size is not divisible by 2^J")
+ dict <- wave.filter(filt)
+ L <- dict$length
+ storage.mode(L) <- "integer"
+ h <- dict$hpf
+ storage.mode(h) <- "double"
+ g <- dict$lpf
+ storage.mode(g) <- "double"
+ start <- 1
+ stop <- N/2
+ Y <- rep(0,N)
+ for (j in 1:J)
+ {
+ W <- V <- numeric(N/2^j)
+ out <- .C("dwt", as.double(x), as.integer(N/2^(j - 1)),L, h, g, W = as.double(W), V = as.double(V), PACKAGE = "waveslim")[6:7]
+ Y[start:stop] <- out$W
+ x <- out$V
+ start <- start+N/(2^j)
+ stop <- stop+N/2^(j+1)
+ }
+ stop <- N
+ Y[start:stop] <- x
+ return(Y)
+}
+
+wave.backtransform <- function(D,n.levels,filt="haar")
+{
+ p <- length(D)
+ if ((p%%2^n.levels)!=0) stop("Sample size is not divisible by 2^n.levels")
+ dict <- wave.filter(filt)
+ L <- dict$length
+ storage.mode(L) <- "integer"
+ h <- dict$hpf
+ storage.mode(h) <- "double"
+ g <- dict$lpf
+ storage.mode(g) <- "double"
+ stop <- p
+ start <- p-p/2^n.levels+1
+ X <- D[start:stop]
+ for (j in n.levels:1)
+ {
+ stop <- start-1
+ start <- stop-p/2^j+1
+ N <- length(X)
+ XX <- numeric(2*p/2^j)
+ X <- .C("idwt", as.double(D[start:stop]), as.double(X), as.integer(N),L, h, g, out = as.double(XX), PACKAGE = "waveslim")$out
+ }
+ return(X)
+}
+
+wave.backtransformK <- function(Cum,n.levels,order=1,filt="haar")
+{
+ p <- length(Cum)
+ if ((p%%2^n.levels)!=0) stop("Sample size is not divisible by 2^n.levels")
+ dict <- wave.filter(filt)
+ L <- dict$length
+ storage.mode(L) <- "integer"
+ h <- dict$hpf^order
+ storage.mode(h) <- "double"
+ g <- dict$lpf^order
+ storage.mode(g) <- "double"
+ stop <- p
+ start <- p-p/2^n.levels+1
+ X <- Cum[start:stop]
+ for (j in n.levels:1)
+ {
+ stop <- start-1
+ start <- stop-p/2^j+1
+ N <- length(X)
+ XX <- numeric(2*p/2^j)
+ X <- .C("idwt", as.double(Cum[start:stop]), as.double(X), as.integer(N),L, h, g, out = as.double(XX), PACKAGE = "waveslim")$out
+ }
+ return(X)
+}
+
+MapMar <- function(D,X,vareps,J=0,eqsmooth=TRUE)
+{
+ N <- nrow(X)
+ q <- ncol(X)
+ K <- ncol(D)
+ if (J==0)
+ {
+ ends <- K
+ } else
+ {
+ ends <- cumsum(K/2^c(1:J,J))
+ }
+ B <- matrix(0,nrow=q,ncol=K)
+ varB <- B
+ phi <- B
+ if (eqsmooth)
+ {
+ out <- .C("MAPMARGEQSMOOTH",as.double(D),as.integer(K),as.double(vareps),as.double(B),as.double(varB),as.double(phi),as.double(X),as.double(diag(t(X)%*%X)),as.integer(q),as.integer(N),as.integer(ends), PACKAGE = "waveTiling")[4:6]
+ } else
+ {
+ out <- .C("MAPMARG",as.double(D),as.integer(K),as.double(vareps),as.double(B),as.double(varB),as.double(phi),as.double(X),as.double(diag(t(X)%*%X)),as.integer(q),as.integer(N),as.integer(ends), PACKAGE = "waveTiling")[4:6]
+ }
+ names(out) <- c("beta_MAP","varbeta_MAP","phi")
+ dim(out[[1]]) <- c(K,q)
+ out[[1]] <- t(out[[1]])
+ dim(out[[2]]) <- c(K,q)
+ out[[2]] <- t(out[[2]])
+ dim(out[[3]]) <- c(K,q)
+ out[[3]] <- t(out[[3]])
+ return(out)
+}
+
+MapMarImp <- function(D,X,vareps,J=0,eqsmooth=TRUE)
+{
+ N <- nrow(X)
+ q <- ncol(X)
+ K <- ncol(D)
+ if (J==0)
+ {
+ ends <- K
+ } else
+ {
+ ends <- cumsum(K/2^c(1:J,J))
+ }
+ B <- matrix(0,nrow=q,ncol=K)
+ varB <- B
+ phi <- B
+ if (eqsmooth)
+ {
+ out <- .C("MAPMARGIMPEQSMOOTH",as.double(D),as.integer(K),as.double(vareps),as.double(B),as.double(varB),as.double(phi),as.double(X),as.double(diag(t(X)%*%X)),as.integer(q),as.integer(N),as.integer(ends), PACKAGE = "waveTiling")[4:6]
+ } else
+ {
+ out <- .C("MAPMARGIMP",as.double(D),as.integer(K),as.double(vareps),as.double(B),as.double(varB),as.double(phi),as.double(X),as.double(diag(t(X)%*%X)),as.integer(q),as.integer(N),as.integer(ends), PACKAGE = "waveTiling")[4:6]
+ }
+ names(out) <- c("beta_MAP","varbeta_MAP","phi")
+ dim(out[[1]]) <- c(K,q)
+ out[[1]] <- t(out[[1]])
+ dim(out[[2]]) <- c(K,q)
+ out[[2]] <- t(out[[2]])
+ dim(out[[3]]) <- c(K,q)
+ out[[3]] <- t(out[[3]])
+ return(out)
+}
+
+
+WaveMarEstVarJ <- function(Y,X,n.levels,wave.filt,prior=c("normal","improper"),eqsmooth=TRUE,saveall=FALSE,D,var.eps,max.it,tol=1e-6,trace=FALSE)
+{
+ N <- nrow(D)
+ K <- ncol(D)
+ Kj <- K/2^c(1:n.levels,n.levels)
+ ends <- cumsum(Kj)
+ starts <- c(1,ends[-(n.levels+1)]+1)
+ varEps <- rep(mad(D[,starts[1]:ends[1]])^2,n.levels+1)
+ if (var.eps=="mad")
+ {
+ max.it <- -1
+ saveall <- TRUE
+ }
+ if (max.it<0)
+ {
+ for (i in 2:(n.levels+1))
+ {
+ varEps[i] <- ifelse(i<(n.levels+1),mad(D[,starts[i]:ends[i]])^2,varEps[i-1])
+ }
+ }
+ crit1 <- sum(D^2)
+ if (max.it<1)
+ {
+ if (prior=="normal")
+ {
+ WaveFit <- MapMar(D,X,varEps,n.levels,eqsmooth)
+ }
+ if (prior=="improper")
+ {
+ WaveFit <- MapMarImp(D,X,varEps,n.levels,eqsmooth)
+ }
+ } else
+ {
+ if (trace==TRUE) message("Gauss Seidel algorithm...")
+ for (j in 1:max.it)
+ {
+ crit0 <- crit1
+ if (prior=="normal")
+ {
+ WaveFit <- MapMar(D,X,varEps,n.levels,eqsmooth)
+ }
+ if (prior=="improper")
+ {
+ WaveFit <- MapMarImp(D,X,varEps,n.levels,eqsmooth)
+ }
+ for (i in 1:(n.levels+1))
+ {
+ VinvD <- D[,starts[i]:ends[i]]-X%*%(t(X)%*%D[,starts[i]:ends[i]]/(diag(t(X)%*%X)+1/WaveFit$phi[,starts[i]:ends[i]]))
+ varEps[i] <- sum(D[,starts[i]:ends[i]]*VinvD)/(ends[i]-starts[i]+1)/N
+ }
+ crit1 <- -1/2*sum(log(t(X)%*%X%*%WaveFit$phi+rep(1,ncol(X))))-1/2*(log(varEps)%*%Kj)-N*K*log(2*pi)/2-N*K/2
+ if (trace==TRUE) message("iteration ", j," of ",max.it)
+ if ((abs((crit0-crit1)/crit0))<tol)
+ {
+ if (trace==TRUE) message("Gauss Seidel algorithm converged")
+ break
+ }
+ }
+ if ((abs(crit0-crit1)/crit0)>tol) message("Warning: The maximum number of iterations reached without convergence")
+ }
+ if (!saveall)
+ {
+ WaveFit$beta_MAP <- matrix()
+ WaveFit$varbeta_MAP <- matrix()
+ WaveFit$phi <- matrix()
+ }
+ WaveFit$varEps <- varEps
+ WaveFit$n.levels <- n.levels
+ WaveFit$K <- K
+ WaveFit$N <- N
+ WaveFit$Kj <- Kj
+ WaveFit$X <- X
+ WaveFit$wave.filt <- wave.filt
+ return(WaveFit)
+}
+
+## plot
+# annoFIle needs the following columns: "chromosome", "strand", "feature", "ID", "start", "end"
+makeNewAnnotationTrack <- function(annoFile,chromosome,strand,minBase,maxBase,feature="exon",dp=NULL)
+{
+ if (is.null(dp))
+ {
+ dp <- DisplayPars(exon="lightblue")
+ }
+ return(makeAnnotationTrack(regions=annoFile[(annoFile$chromosome==chromosome)&(annoFile$strand==strand)&(annoFile$start<maxBase)&(annoFile$end>minBase)&(annoFile$feature==feature),c("start","end","feature","group","ID")],dp=dp))
+}
+
+
+makeNewAnnotationTextOverlay <- function(annoFile,chromosome,strand,minBase,maxBase,region,feature=c("gene","transposable_element_gene"),y=.5,dp=NULL)
+{
+ if (is.null(dp))
+ {
+ dp=DisplayPars(cex=1)
+ }
+ annohlp <- annoFile[(annoFile$chromosome==chromosome)&(annoFile$strand==strand)&(annoFile$start<maxBase)&(annoFile$end>minBase)&(is.element(annoFile$feature,feature)),c("ID","start","end")]
+ return(makeTextOverlay(annohlp$ID,xpos=(annohlp$start+annohlp$end)/2,ypos=y,region=region,dp=dp))
+}
+
+makeNewTranscriptRectangleOverlay <- function(sigRegions,locations,start,end,region=NULL,dp=NULL)
+{
+ if (is.null(dp))
+ {
+ dp=DisplayPars(color="black")
+ }
+ detectedRegionsSelect <- matrix(sigRegions[sigRegions[,1]<end&sigRegions[,2]>start,],ncol=2)
+ detectedRegionsSelect[detectedRegionsSelect[,1]<start,1] <- start
+ detectedRegionsSelect[detectedRegionsSelect[,2]>end,2] <- end
+ return(makeRectangleOverlay(start=locations[detectedRegionsSelect[,1]],end =locations[detectedRegionsSelect[,2]],region=region,dp=dp))
+}
+
+#### makeDesign
+makeDesign<-function(design=c("time","circadian","group","factorial"),replics, noGroups, factor.levels=NULL) {
+# if (method=="twoGroup")
+# {
+# Xorig <- matrix(0,nrow=dim(pData(object))[1],ncol=2)
+# Xorig[,1] <- 1
+# Xorig[,2] <- rep(c(1,-1),replics) ## contr.helmert
+# }
+ if (design=="time") #time
+ {
+# Xorig <- matrix(0,nrow=dim(pData(object))[1],ncol=noGroups)
+ Xorig <- matrix(0,nrow=noGroups*replics,ncol=noGroups)
+ orderedFactor <- factor(1:noGroups,ordered=TRUE)
+ desPoly <- lm(rnorm(noGroups)~orderedFactor,x=TRUE)$x
+ Xorig[,1] <- 1
+ Xorig[,2:noGroups] <- apply(desPoly[,2:noGroups],2,rep,replics)
+
+ }
+ else if (design=="circadian") # circadian
+ {
+# Xorig <- matrix(0,nrow=dim(pData(object))[1],ncol=3)
+ Xorig <- matrix(0,nrow=noGroups*replics,ncol=3)
+ Xorig[,1] <- 1
+ Xorig[,2] <- rep(sin(seq(0,2*pi-(pi/noGroups),2*pi/noGroups)),replics)
+ Xorig[,3] <- rep(cos(seq(0,2*pi-(pi/noGroups),2*pi/noGroups)),replics)
+ }
+ else if (design=="group") # group
+ {
+# Xorig <- matrix(0,nrow=dim(pData(object))[1],ncol=noGroups)
+ Xorig <- matrix(0,nrow=noGroups*replics,ncol=noGroups)
+ desHelmert <- contr.helmert(noGroups)
+ Xorig[,1] <- 1
+ Xorig[,2:noGroups] <- apply(desHelmert[,1:(noGroups-1)],2,rep,replics)
+ }
+ else if (design=="factorial")
+ {
+ if (is.null(factor.levels))
+ {
+ stop("No proper factor levels given.")
+ }
+ Xorig <- matrix(0,nrow=sum(replics),ncol=prod(factor.levels))
+ Xorig[,1] <- 1
+ desTrt1 <- contr.treatment(factor.levels[1])
+ desTrt2 <- contr.treatment(factor.levels[2])
+ #desHelmert1 <- contr.helmert(factor.levels[1])
+ #desHelmert2 <- contr.helmert(factor.levels[2])
+ replicsMat <- matrix(replics,nrow=factor.levels[1],ncol=factor.levels[2],byrow=TRUE)
+ desFac1 <- apply(desTrt1,2,function(x) unlist(sapply(1:nrow(replicsMat),function(y) rep(x[y],sum(replicsMat[y,])))))
+ Xorig[,2:(factor.levels[1])] <- desFac1
+ desFac2 <- apply(desTrt2,2,function(x) unlist(sapply(1:nrow(replicsMat),function(y) rep(x,replicsMat[y,]))))
+ Xorig[,(factor.levels[1]+1):(factor.levels[1]+factor.levels[2]-1)] <- desFac2
+ # include interactions
+ colno <- sum(factor.levels)
+ for (i in 1:ncol(desFac1))
+ {
+ for (j in 1:ncol(desFac2))
+ {
+ Xorig[,colno] <- desFac1[,i]*desFac2[,j]
+ colno <- colno + 1
+ }
+ }
+ }
+ else {
+ stop("No proper design specified!")
+ }
+ out <- list()
+ X.qr <- qr(Xorig)
+ out$Xorthnorm <- qr.Q(X.qr)
+ out$Xorig <- Xorig
+ return (out);
+}
+
+#### makeContrasts
+makeContrasts <- function(contrasts, nlevels) {
+ if (contrasts=="compare") {
+ q <- nlevels*(nlevels-1)/2
+ contr <- matrix(0,nrow=q,ncol=nlevels)
+ hlp1 <- rep(2:nlevels,1:(nlevels-1))
+ hlp2 <- unlist(sapply(1:(nlevels-1),function(x) seq(1:x)))
+ for (i in 1:nrow(contr))
+ {
+ contr[i,hlp1[i]] <- 1
+ contr[i,hlp2[i]] <- -1
+ }
+ }
+
+ return (contr);
+}
+
+
Added: pkg/waveTiling/R/initialize-methods.R
===================================================================
--- pkg/waveTiling/R/initialize-methods.R (rev 0)
+++ pkg/waveTiling/R/initialize-methods.R 2012-05-21 12:54:45 UTC (rev 42)
@@ -0,0 +1,128 @@
+#MapFilterProbe
+setMethod("initialize","MapFilterProbe",function(.Object,filteredIndices,chromosome,position,strand)
+{
+ .Object at filteredIndices <- filteredIndices
+ .Object at chromosome <- chromosome
+ .Object at position <- position
+ .Object at strand <- strand
+ return(.Object)
+})
+
+#WaveTilingFeatureSet
+setMethod("initialize","WaveTilingFeatureSet",function (.Object)
+{
+ callNextMethod(.Object);
+}
+)
+
+
+#GenomeInfo
+setMethod("initialize","GenomeInfo",function(.Object,chromosome,strand,minPos,maxPos)
+{
+ .Object at chromosome <- chromosome
+ .Object at strand <- strand
+ .Object at minPos <- minPos
+ .Object at maxPos <- maxPos
+ return(.Object)
+})
+
+#WfmFit
+setMethod("initialize","WfmFit",function(.Object,betaWav,varbetaWav,smoothPar,varEps,dataOrigSpace,dataWaveletSpace,design.matrix,phenoData,genome.info,n.levels,probePosition,wave.filt,Kj,prior,F,varF,P,Z,noGroups,replics)
+{
+ .Object at betaWav <- betaWav
+ .Object at varbetaWav <- varbetaWav
+ .Object at smoothPar <- smoothPar
+ .Object at varEps <- varEps
+ .Object at dataOrigSpace <- dataOrigSpace
+ .Object at dataWaveletSpace <- dataWaveletSpace
+ .Object at design.matrix <- design.matrix
+ .Object at phenoData <- phenoData
+# .Object at design <- design
+ .Object at genome.info <- genome.info
+ .Object at n.levels <- n.levels
+ .Object at probePosition <- probePosition
+ .Object at wave.filt <- wave.filt
+ .Object at Kj <- Kj
+ .Object at prior <- prior
+# .Object at rescale <- rescale
+ .Object at F <- F
+ .Object at varF <- varF
+ .Object at P <- P
+ .Object at Z <- Z
+ .Object at noGroups <- noGroups
+ .Object at replics <- replics
+ return(.Object)
+})
+
+
+# WfmInf
+setMethod("initialize","WfmInf",function(.Object, alpha, delta, two.sided, sigProbes, regions, GlocRegions, FDR, CI, eff, varEff)
+{
+ .Object at alpha <- alpha
+ .Object at delta <- delta
+ .Object at two.sided <- two.sided
+ .Object at sigProbes <- sigProbes
+ .Object at regions <- regions
+ .Object at GlocRegions <- GlocRegions
+ .Object at FDR <- FDR
+ .Object at CI <- CI
+ .Object at eff <- eff
+ .Object at varEff <- varEff
+ return(.Object)
+})
+
+
+## Subclasses
+setMethod("initialize","WfmFitTime",function (.Object, ...)
+{
+ callNextMethod(.Object, ...);
+}
+)
+
+setMethod("initialize","WfmFitCircadian",function (.Object, ...)
+{
+ callNextMethod(.Object, ...);
+}
+)
+
+setMethod("initialize","WfmFitFactor",function (.Object, ...)
+{
+ callNextMethod(.Object, ...);
+}
+)
+
+setMethod("initialize","WfmFitCustom",function (.Object, ...)
+{
+ callNextMethod(.Object, ...);
+}
+)
+
+setMethod("initialize","WfmInfCustom",function (.Object, ...)
+{
+ callNextMethod(.Object, ...);
+}
+)
+
+setMethod("initialize","WfmInfCompare",function (.Object, ...)
+{
+ callNextMethod(.Object, ...);
+}
+)
+
+setMethod("initialize","WfmInfEffects",function (.Object, ...)
+{
+ callNextMethod(.Object, ...);
+}
+)
+
+setMethod("initialize","WfmInfMeans",function (.Object, ...)
+{
+ callNextMethod(.Object, ...);
+}
+)
+
+setMethod("initialize","WfmInfOverallMean",function (.Object, ...)
+{
+ callNextMethod(.Object, ...);
+}
+)
Added: pkg/waveTiling/R/methods-MapFilterProbe.R
===================================================================
--- pkg/waveTiling/R/methods-MapFilterProbe.R (rev 0)
+++ pkg/waveTiling/R/methods-MapFilterProbe.R 2012-05-21 12:54:45 UTC (rev 42)
@@ -0,0 +1,48 @@
+setMethod("getFilteredIndices",signature("MapFilterProbe"),function(object)
+{
+ return(object at filteredIndices)
+}
+)
+
+setMethod("getChromosome",signature("MapFilterProbe"),function(object)
+{
+ return(object at chromosome)
+}
+)
+
+setMethod("getPosition",signature("MapFilterProbe"),function(object)
+{
+ return(object at position)
+}
+)
+
+setMethod("getStrand",signature("MapFilterProbe"),function(object)
+{
+ return(object at strand)
+}
+)
+
+setMethod("selectProbesFromFilterOverlap",signature("MapFilterProbe"),function(object,chromosome,strand=c("forward","reverse"),minPos=min(getPosition(object)),maxPos=max(getPosition(object)))
+{
+ if (class(object)!="MapFilterProbe")
+ {
+ stop("class of object is not MapFilterProbe. Use 'filterOverlap()' to create such an object.")
+ }
+ if ((length(grep("chr",chromosome))>0) | (length(grep("Chr",chromosome))>0))
+ {
+ stop("give only the number (or letter) in the chromosome argument.")
+ }
+ if (minPos > maxPos)
+ {
+ stop("minPos is greater than maxPos")
+ }
+ selChrom <- (1:length(getFilteredIndices(object)))[getChromosome(object)==paste("chr",as.character(chromosome),sep="") | getChromosome(object)==paste("Chr",as.character(chromosome),sep="")]
+ selStrand <- (1:length(getFilteredIndices(object)))[getStrand(object)==strand]
+ selHlp <- intersect(selChrom,selStrand)
+ selPos <- (1:length(getFilteredIndices(object)))[(getPosition(object)>=minPos)&(getPosition(object)<=maxPos)]
+ selectionInit <- intersect(selHlp,selPos)
+ selection <- getFilteredIndices(object)[selectionInit]
+ return(list(selection=selection,selectionInit=selectionInit))
+}
+)
+
Added: pkg/waveTiling/R/methods-WaveTilingFeatureSet.R
===================================================================
--- pkg/waveTiling/R/methods-WaveTilingFeatureSet.R (rev 0)
+++ pkg/waveTiling/R/methods-WaveTilingFeatureSet.R 2012-05-21 12:54:45 UTC (rev 42)
@@ -0,0 +1,445 @@
+# temporary (?)
+# I can only use this after inserting strand information to the database when parsing the bpmap and Cel file in the (updated) package PDInfoBuilder
+# setMethod("pmStrand","WaveTilingFeatureSet",function(object)
+# {
+# conn <- db(object)
+# sql <- paste("SELECT fid, strand", "FROM pmfeature", "INNER JOIN chrom_dict", "USING(chrom)")
+# tmp <- dbGetQuery(conn, sql)
+# tmp <- tmp[order(tmp[["fid"]]),]
+# return(tmp[["strand"]])
+# }
+# )
+
+
+setMethod("addPheno",signature("WaveTilingFeatureSet"),function(object,noGroups,groupNames,replics,...)
+{
+ if (dim(object)[2] != sum(replics))
+ {
+ stop("Number of samples in object differs from total number of replicates.")
+ }
+ if (length(groupNames) != noGroups)
+ {
+ stop("Number of group names differs from number of groups")
+ }
+ if (length(replics)==1)
+ {
+ warning("Only one number of replicates given. This number is used for all groups")
+ replics <- rep(replics,noGroups)
+ }
+ if ((length(replics)!=1) & (length(replics)!=noGroups))
+ {
+ stop("Length of replicates vector differs from number of groups")
+ }
+ pData <- data.frame(matrix(NA,sum(replics),2,dimnames=list(NULL,c("group","replicate"))))
+ pData[1:sum(replics),1] <- rep(groupNames,replics)
+ pData[1:sum(replics),2] <- sequence(replics)
+ rownames(pData) <- paste(pData[,1],pData[,2],sep=".")
+ metadata <- data.frame(labelDescription=c("sample group","replicate number within sample group"))
+ phenoData <- new("AnnotatedDataFrame",data=pData,varMetadata=metadata)
+ phenoData(object) <- phenoData
+ return(object)
+}
+)
+
+setMethod("getNoGroups",signature("WaveTilingFeatureSet"),function(object)
+{
+ if ((names(pData(object))[1]!="group")|((names(pData(object))[2]!="replicate")))
+ {
+ stop("phenoData of WaveTilingFeatureSet object is not correct. Use 'addPheno()' first.")
+ }
+ noGroups <- length(table(pData(object)$group))
+ return(noGroups)
+}
+)
+
+setMethod("getGroupNames",signature("WaveTilingFeatureSet"),function(object)
+{
+ if ((names(pData(object))[1]!="group")|((names(pData(object))[2]!="replicate")))
+ {
+ stop("phenoData of TilingFeatureSet object is not correct. Use 'addPheno()' first.")
+ }
+ groupNames <- unique(pData(object)$group)
+ return(groupNames)
+}
+)
+
+setMethod("getReplics",signature("WaveTilingFeatureSet"),function(object)
+{
+ if ((names(pData(object))[1]!="group")|((names(pData(object))[2]!="replicate")))
+ {
+ stop("phenoData of WaveTilingFeatureSet object is not correct. Use 'addPheno()' first.")
+ }
+ # keep the order safe
+ orderHlp <- order(unique(pData(object)$group))
+ replcsTemp <- as.numeric(table(pData(object)$group))
+ replcs <- replcsTemp[orderHlp]
+ return(replcs)
+}
+)
+
+setMethod("filterOverlap",signature("WaveTilingFeatureSet"),function(object,remap=TRUE,BSgenomeObject,chrId,strand=c("forward","reverse","both"),MM=FALSE)
+{
+ if (!inherits(object,"WaveTilingFeatureSet")) #class(object)!="WaveTilingFeatureSet")
+ {
+ stop("class of object is not WaveTilingFeatureSet.")
+ }
+ pmIndex <- oligo::pmindex(object)
+ dataPMSeq <- pmSequence(object)
+ if (remap)
+ {
+ setTrustBand <- min(unique(width(dataPMSeq)))
+ if (setTrustBand < 15)
+ {
+ warning("There are probes smaller than 15 bp")
+ }
+ minTrustBand <- 15
+ trBand <- max(setTrustBand,minTrustBand)
+ if ((strand=="forward") | (strand=="both"))
+ {
+ message("begin remapping forward strand...")
+ message("extract probe sequences")
+ pmSeqDict <- PDict(dataPMSeq,tb.start=1,tb.end=trBand)
+ message("match probe sequences to DNA sequence")
+ chrSeqList <- list()
+ for (i in chrId)
+ {
+ chrSeqList[[i]] <- BSgenomeObject[[i]]
+ }
+ names(chrSeqList) <- seqnames(BSgenomeObject)[chrId]
+ startPM <- lapply(chrSeqList,function(x) startIndex(matchPDict(pmSeqDict,x)))
+ nposChrPM <- lapply(startPM,function(x) sapply(x,length))
+ pmMatch <- Reduce("+",nposChrPM)==1
+ pmMatchIndex <- (1:length(pmMatch))[pmMatch]
+ index <- lapply(nposChrPM,function(x)
+ {
+ indexhlp <- x==1
+ index <- indexhlp[pmMatchIndex]
+ })
+ mp <- stack(Map(which,index))
+ chrInit <- c()
+ chrInit[mp[,1]] <- as.character(mp[,2])
+ posInit <- c()
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/wavetiling -r 42
More information about the Wavetiling-commits
mailing list