[Wavetiling-commits] r43 - in pkg: waveTiling waveTiling/R waveTiling/inst/scripts waveTiling/man waveTilingData waveTilingData/data waveTilingData/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jun 4 14:28:08 CEST 2012
Author: kdbeuf
Date: 2012-06-04 14:28:07 +0200 (Mon, 04 Jun 2012)
New Revision: 43
Added:
pkg/waveTilingData/data/leafdevBQ.rda
pkg/waveTilingData/man/leafdevBQ.Rd
pkg/waveTilingData/man/waveTilingData.Rd
Removed:
pkg/waveTiling/CHANGELOG
pkg/waveTiling/TODO
pkg/waveTilingData/data/tair9gff.rda
pkg/waveTilingData/man/tair9gff.Rd
Modified:
pkg/waveTiling/DESCRIPTION
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-WfmFit.R
pkg/waveTiling/R/methods-WfmInf.R
pkg/waveTiling/R/show-methods.R
pkg/waveTiling/inst/scripts/waveTiling-vignette.Rnw
pkg/waveTiling/man/MapFilterProbe-class.Rd
pkg/waveTiling/man/WaveTilingFeatureSet-class.Rd
pkg/waveTiling/man/WfmFit-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/getNonAnnotatedRegions.Rd
pkg/waveTiling/man/getSigGenes.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/waveTilingData/DESCRIPTION
pkg/waveTilingData/data/leafdevFit.rda
pkg/waveTilingData/data/leafdevInfCompare.rda
pkg/waveTilingData/man/leafdev.Rd
pkg/waveTilingData/man/leafdevFit.Rd
pkg/waveTilingData/man/leafdevInfCompare.Rd
pkg/waveTilingData/man/leafdevMapAndFilterTAIR9.Rd
Log:
solving several issues
Deleted: pkg/waveTiling/CHANGELOG
===================================================================
--- pkg/waveTiling/CHANGELOG 2012-05-21 12:54:45 UTC (rev 42)
+++ pkg/waveTiling/CHANGELOG 2012-06-04 12:28:07 UTC (rev 43)
@@ -1,3 +0,0 @@
-2011-12-04 Kristof De Beuf <kristof.debeuf at ugent.be> - committed version 0.1.1
-
-* First version
Modified: pkg/waveTiling/DESCRIPTION
===================================================================
--- pkg/waveTiling/DESCRIPTION 2012-05-21 12:54:45 UTC (rev 42)
+++ pkg/waveTiling/DESCRIPTION 2012-06-04 12:28:07 UTC (rev 43)
@@ -1,5 +1,5 @@
Package: waveTiling
-Version: 0.99.0
+Version: 0.99.1
Date: 2012-05-14
License: GPL (>=2)
Title: Wavelet-Based Models for Tiling Array Transcriptome Analysis
@@ -7,7 +7,7 @@
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
+Suggests: BSgenome, BSgenome.Athaliana.TAIR.TAIR9, waveTilingData, pd.atdschip.tiling, TxDb.Athaliana.BioMart.plantsmart12
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
Modified: pkg/waveTiling/R/allClasses.R
===================================================================
--- pkg/waveTiling/R/allClasses.R 2012-05-21 12:54:45 UTC (rev 42)
+++ pkg/waveTiling/R/allClasses.R 2012-06-04 12:28:07 UTC (rev 43)
@@ -26,7 +26,7 @@
#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")
+ representation = representation (alpha="numeric",delta="numeric",two.sided="numeric",sigProbes="list",regions="list",GlocRegions="list",FDR="matrix",CI="array",eff="matrix",varEff="matrix",genome.info="GenomeInfo")
)
## Compare
Modified: pkg/waveTiling/R/allGenerics.R
===================================================================
--- pkg/waveTiling/R/allGenerics.R 2012-05-21 12:54:45 UTC (rev 42)
+++ pkg/waveTiling/R/allGenerics.R 2012-06-04 12:28:07 UTC (rev 43)
@@ -291,19 +291,19 @@
}
)
-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)
+setGeneric("plotWfm",function(fit, inf, biomartObj, minPos, maxPos, trackFeature="exon", two.strand=TRUE, plotData=TRUE, plotMean=TRUE, tracks=0)
{
standardGeneric("plotWfm")
}
)
-setGeneric("getSigGenes",function(fit, inf, annoFile)
+setGeneric("getSigGenes",function(fit, inf, biomartObj)
{
standardGeneric("getSigGenes")
}
)
-setGeneric("getNonAnnotatedRegions",function(fit, inf, annoFile)
+setGeneric("getNonAnnotatedRegions",function(fit, inf, biomartObj)
{
standardGeneric("getNonAnnotatedRegions")
}
Modified: pkg/waveTiling/R/helperFunctions.R
===================================================================
--- pkg/waveTiling/R/helperFunctions.R 2012-05-21 12:54:45 UTC (rev 42)
+++ pkg/waveTiling/R/helperFunctions.R 2012-06-04 12:28:07 UTC (rev 43)
@@ -252,24 +252,50 @@
## plot
# annoFIle needs the following columns: "chromosome", "strand", "feature", "ID", "start", "end"
-makeNewAnnotationTrack <- function(annoFile,chromosome,strand,minBase,maxBase,feature="exon",dp=NULL)
+# 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))
+# }
+
+makeNewAnnotationTrack <- function(biomartObj,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))
+ exn <- exons(biomartObj)
+ exnId <- which((as.character(seqnames(exn))==as.character(chromosome)&(as.character(strand(exn))==strand)&(start(ranges(exn))>minBase)&(end(ranges(exn))<maxBase)))
+ exnSel <- exn[exnId]
+ regionObj <- data.frame(start=start(ranges(exnSel)),end=end(ranges(exnSel)),feature=feature,group=elementMetadata(exnSel)$exon_id,ID=rep(NA,length(exnSel)))
+ return(makeAnnotationTrack(regions=regionObj,dp=dp))
}
-makeNewAnnotationTextOverlay <- function(annoFile,chromosome,strand,minBase,maxBase,region,feature=c("gene","transposable_element_gene"),y=.5,dp=NULL)
+# 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))
+# }
+
+makeNewAnnotationTextOverlay <- function(biomartObj,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))
+ trs <- transcripts(biomartObj)
+ trsId <- which((as.character(seqnames(trs))==as.character(chromosome)&(as.character(strand(trs))==strand)&(start(ranges(trs))>minBase)&(end(ranges(trs))<maxBase)))
+ trsSel <- trs[trsId]
+ annohlp <- data.frame(ID=elementMetadata(trsSel)$tx_name,start=start(ranges(trsSel)),end=end(ranges(trsSel)))
+ return(makeTextOverlay(as.character(annohlp$ID),xpos=(annohlp$start+annohlp$end)/2,ypos=y,region=region,dp=dp))
}
makeNewTranscriptRectangleOverlay <- function(sigRegions,locations,start,end,region=NULL,dp=NULL)
Modified: pkg/waveTiling/R/initialize-methods.R
===================================================================
--- pkg/waveTiling/R/initialize-methods.R 2012-05-21 12:54:45 UTC (rev 42)
+++ pkg/waveTiling/R/initialize-methods.R 2012-06-04 12:28:07 UTC (rev 43)
@@ -56,7 +56,7 @@
# WfmInf
-setMethod("initialize","WfmInf",function(.Object, alpha, delta, two.sided, sigProbes, regions, GlocRegions, FDR, CI, eff, varEff)
+setMethod("initialize","WfmInf",function(.Object, alpha, delta, two.sided, sigProbes, regions, GlocRegions, FDR, CI, eff, varEff, genome.info)
{
.Object at alpha <- alpha
.Object at delta <- delta
@@ -68,6 +68,7 @@
.Object at CI <- CI
.Object at eff <- eff
.Object at varEff <- varEff
+ .Object at genome.info <- genome.info
return(.Object)
})
Modified: pkg/waveTiling/R/methods-WfmFit.R
===================================================================
--- pkg/waveTiling/R/methods-WfmFit.R 2012-05-21 12:54:45 UTC (rev 42)
+++ pkg/waveTiling/R/methods-WfmFit.R 2012-06-04 12:28:07 UTC (rev 43)
@@ -455,19 +455,19 @@
#genomeLoc <- new("genomeInfo",chromosome=chromosome,strand=strand,minPos=min(Gloc),maxPos=max(Gloc))
#infObject <- new("Wfm",method=method,alpha=alpha,delta=delta,two.sided=two.sided,sigProbes=sigProbes,regions=regions,GlocRegions=GlocRegions,FDR=FDR,CI=CI)
if (!is.null(contrast.matrix)) {
- infObject <- new("WfmInfCustom",alpha=alpha,delta=delta,two.sided=two.sided,sigProbes=sigProbes,regions=regions,GlocRegions=GlocRegions,FDR=FDR,CI=CI,eff=eff,varEff=varEff)
+ infObject <- new("WfmInfCustom",alpha=alpha,delta=delta,two.sided=two.sided,sigProbes=sigProbes,regions=regions,GlocRegions=GlocRegions,FDR=FDR,CI=CI,eff=eff,varEff=varEff,genome.info=getGenomeInfo(object))
}
else if (contrasts=="compare" ) {
- infObject <- new("WfmInfCompare",alpha=alpha,delta=delta,two.sided=two.sided,sigProbes=sigProbes,regions=regions,GlocRegions=GlocRegions,FDR=FDR,CI=CI,eff=eff,varEff=varEff)
+ infObject <- new("WfmInfCompare",alpha=alpha,delta=delta,two.sided=two.sided,sigProbes=sigProbes,regions=regions,GlocRegions=GlocRegions,FDR=FDR,CI=CI,eff=eff,varEff=varEff,genome.info=getGenomeInfo(object))
}
else if (contrasts=="effects" ) {
- infObject <- new("WfmInfEffects",alpha=alpha,delta=delta,two.sided=two.sided,sigProbes=sigProbes,regions=regions,GlocRegions=GlocRegions,FDR=FDR,CI=CI,eff=eff,varEff=varEff)
+ infObject <- new("WfmInfEffects",alpha=alpha,delta=delta,two.sided=two.sided,sigProbes=sigProbes,regions=regions,GlocRegions=GlocRegions,FDR=FDR,CI=CI,eff=eff,varEff=varEff,genome.info=getGenomeInfo(object))
}
else if (contrasts=="means" ) {
- infObject <- new("WfmInfMeans",alpha=alpha,delta=delta,two.sided=two.sided,sigProbes=sigProbes,regions=regions,GlocRegions=GlocRegions,FDR=FDR,CI=CI,eff=eff,varEff=varEff)
+ infObject <- new("WfmInfMeans",alpha=alpha,delta=delta,two.sided=two.sided,sigProbes=sigProbes,regions=regions,GlocRegions=GlocRegions,FDR=FDR,CI=CI,eff=eff,varEff=varEff,genome.info=getGenomeInfo(object))
}
else if (contrasts=="overallMean" ) {
- infObject <- new("WfmInfOverallMean",alpha=alpha,delta=delta,two.sided=two.sided,sigProbes=sigProbes,regions=regions,GlocRegions=GlocRegions,FDR=FDR,CI=CI,eff=eff,varEff=varEff)
+ infObject <- new("WfmInfOverallMean",alpha=alpha,delta=delta,two.sided=two.sided,sigProbes=sigProbes,regions=regions,GlocRegions=GlocRegions,FDR=FDR,CI=CI,eff=eff,varEff=varEff,genome.info=getGenomeInfo(object))
}
return (infObject);
})
Modified: pkg/waveTiling/R/methods-WfmInf.R
===================================================================
--- pkg/waveTiling/R/methods-WfmInf.R 2012-05-21 12:54:45 UTC (rev 42)
+++ pkg/waveTiling/R/methods-WfmInf.R 2012-06-04 12:28:07 UTC (rev 43)
@@ -48,19 +48,69 @@
{
return(object at varEff)
})
+
+ setMethod("getGenomeInfo",signature("WfmInf"),function(object)
+{
+ return(object at genome.info)
+})
+# setMethod("getSigGenes",signature(fit="WfmFit",inf="WfmInf"),function(fit,inf,annoFile)
+# {
+# #Gloc <- getProbePosition(object)
+# strand <- getStrand(fit)
+# chromosome <- getChromosome(fit)
+# regions <- getGenomicRegions(inf)
+# 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]
+# message("find overlaps with detected regions...")
+# 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)
+# })
-setMethod("getSigGenes",signature(fit="WfmFit",inf="WfmInf"),function(fit,inf,annoFile)
+
+setMethod("getSigGenes",signature(fit="WfmFit",inf="WfmInf"),function(fit,inf,biomartObj)
{
#Gloc <- getProbePosition(object)
strand <- getStrand(fit)
chromosome <- getChromosome(fit)
regions <- getGenomicRegions(inf)
- 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)
+ trs <- transcripts(biomartObj)
if (strand=="forward")
{
strandAlt <- "+"
@@ -70,11 +120,7 @@
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]
+ trsChr <- trs[seqnames(trs)==chromosome]
message("find overlaps with detected regions...")
nList <- length(regions)
annoOver <- GRangesList()
@@ -82,9 +128,9 @@
{
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)
+ overL <- findOverlaps(regGlocGR,trsChr)
regOver <- regGlocGR[queryHits(overL)]
- annoOverj <- annoChrGeneGR[subjectHits(overL)]
+ annoOverj <- trsChr[subjectHits(overL)]
overInt <- pintersect(regOver,annoOverj)
values(annoOverj)$regNo <- queryHits(overL)
values(annoOverj)$percOverGene <- width(overInt)/width(annoOverj)*100
@@ -100,9 +146,61 @@
})
-setMethod("getNonAnnotatedRegions",signature(fit="WfmFit",inf="WfmInf"),function(fit,inf,annoFile)
+# setMethod("getNonAnnotatedRegions",signature(fit="WfmFit",inf="WfmInf"),function(fit,inf,annoFile)
+# {
+# #Gloc <- getProbePosition(object)
+# strand <- getStrand(fit)
+# chromosome <- getChromosome(fit)
+# regions <- getGenomicRegions(inf)
+# if (strand=="forward")
+# {
+# strandAlt <- "+"
+# strandOpp <- "-"
+# } else
+# {
+# strandAlt <- "-"
+# strandOpp <- "+"
+# }
+# message("get annotated regions...")
+# 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)
+# message("find overlaps with detected regions...")
+# 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% as.matrix(overSense)[,1]))
+# noAnnoOppId <- which(!(1:length(regGlocIR) %in% as.matrix(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("getNonAnnotatedRegions",signature(fit="WfmFit",inf="WfmInf"),function(fit,inf,biomartObj)
{
- #Gloc <- getProbePosition(object)
strand <- getStrand(fit)
chromosome <- getChromosome(fit)
regions <- getGenomicRegions(inf)
@@ -116,8 +214,9 @@
strandOpp <- "+"
}
message("get annotated regions...")
- 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")]
+ exn <- exons(biomartObj)
+ annoExons <- exn[(strand(exn)==strandAlt)&(seqnames(exn)==chromosome),]
+ annoExonsOpp <- exn[(strand(exn)==strandOpp)&(seqnames(exn)==chromosome),]
regGlocNoAnnoSense <- list()
regGlocNoAnnoBoth <- list()
nList <- length(regions)
@@ -125,8 +224,8 @@
for (j in 1:nList)
{
regGlocIR <- regions[[j]]
- annoExonsIR <- IRanges(start=annoExons$start,end=annoExons$end)
- annoExonsOppIR <- IRanges(start=annoExonsOpp$start,end=annoExonsOpp$end)
+ annoExonsIR <- ranges(annoExons)
+ annoExonsOppIR <- ranges(annoExonsOpp)
overSense <- findOverlaps(regGlocIR,annoExonsIR)
overOpp <- findOverlaps(regGlocIR,annoExonsOppIR)
noAnnoSenseId <- which(!(1:length(regGlocIR) %in% as.matrix(overSense)[,1]))
@@ -153,15 +252,23 @@
return(out)
})
+
setMethod("plotWfm",signature=c(fit="WfmFit",inf="WfmInf"),
-function(fit,inf,annoFile,minPos,maxPos,trackFeature="exon",overlayFeature=c("gene","transposable_element_gene"),two.strand=TRUE,plotData=TRUE,plotMean=TRUE,tracks=0)
+function(fit,inf,biomartObj,minPos,maxPos,trackFeature="exon",two.strand=TRUE,plotData=TRUE,plotMean=TRUE,tracks=0)
{
- if (missing(annoFile)) {stop("Annotation File is missing!!")}
+ if (missing(biomartObj)) {stop("Annotation object from Biomart is missing!")}
Gloc <- getProbePosition(fit)
if (missing(minPos)) {minPos<-min(Gloc)}
if (missing(maxPos)) {maxPos<-max(Gloc)}
chromosome <- getChromosome(fit)
strand <- getStrand(fit)
+ if (strand=="forward")
+ {
+ strandAlt <- "+"
+ } else
+ {
+ strandAlt <- "-"
+ }
selID <- (1:length(Gloc))[Gloc>minPos & Gloc<maxPos]
sta <- min(selID)
end <- max(selID)
@@ -180,19 +287,19 @@
}
if (two.strand==TRUE)
{
- trackInfo[[trackCount]] <- makeNewAnnotationTrack(annoFile=annoFile,chromosome=chromosome,minBase=minBase,maxBase=maxBase,strand="forward",feature=trackFeature,dp=NULL)
+ trackInfo[[trackCount]] <- makeNewAnnotationTrack(biomartObj=biomartObj,chromosome=chromosome,minBase=minBase,maxBase=maxBase,strand="+",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)
+ overlayInfo[[trackCount]] <- makeNewAnnotationTextOverlay(biomartObj=biomartObj,chromosome=chromosome,minBase=minBase,maxBase=maxBase,strand="+",region=c(trackCount,trackCount),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)
+ trackInfo[[trackCount]] <- makeNewAnnotationTrack(biomartObj=biomartObj,chromosome=chromosome,minBase=minBase,maxBase=maxBase,strand="-",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)
+ overlayInfo[[trackCount]] <- makeNewAnnotationTextOverlay(biomartObj=biomartObj,chromosome=chromosome,minBase=minBase,maxBase=maxBase,strand="-",region=c(trackCount,trackCount),y=0.5)
trackCount <- trackCount + 1
} else
{
- trackInfo[[trackCount]] <- makeNewAnnotationTrack(annoFile=annoFile,chromosome=chromosome,minBase=minBase,maxBase=maxBase,strand=strand,feature=trackFeature)
+ trackInfo[[trackCount]] <- makeNewAnnotationTrack(biomartObj=biomartObj,chromosome=chromosome,minBase=minBase,maxBase=maxBase,strand=strandAlt,feature=trackFeature)
if (strand=="forward")
{
names(trackInfo)[trackCount] <- "F"
@@ -201,7 +308,7 @@
{
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)
+ overlayInfo[[trackCount]] <- makeNewAnnotationTextOverlay(biomartObj=biomartObj,chromosome=chromosome,minBase=minBase,maxBase=maxBase,strand=strandAlt,region=c(trackCount,trackCount),y=0.5)
trackCount <- trackCount + 1
gAxis <- makeGenomeAxis(add53 = TRUE,add35 = TRUE)
trackCount <- trackCount + 1
Modified: pkg/waveTiling/R/show-methods.R
===================================================================
--- pkg/waveTiling/R/show-methods.R 2012-05-21 12:54:45 UTC (rev 42)
+++ pkg/waveTiling/R/show-methods.R 2012-06-04 12:28:07 UTC (rev 43)
@@ -10,8 +10,8 @@
cat("Genome Info :\n")
cat("\tChromosome:",object at chromosome,"\n")
cat("\tStrand:",object at strand,"\n")
- cat("\tMinimum probe position:",object at minPos,"\n")
- cat("\tMaximum probe position:",object at maxPos,"\n")
+ cat("\tMinimum probe position:",format(object at minPos,scientific=FALSE),"\n")
+ cat("\tMaximum probe position:",format(object at maxPos,scientific=FALSE),"\n")
cat("\n")
}
)
@@ -54,10 +54,40 @@
}
)
-# Needs further implementation :: what to show?
setMethod("show",signature("WfmInf"),function(object)
{
- cat("object of class 'WfmInf'\n")
- cat("\n")
+ show(object at genome.info);
}
-)
\ No newline at end of file
+)
+
+setMethod("show",signature("WfmInfCompare"),function(object)
+{
+ cat("Inference object from wavelet based functional model - Pairwise Comparison\n")
+ callNextMethod();
+})
+
+setMethod("show",signature("WfmInfCustom"),function(object)
+{
+ cat("Inference object from wavelet based functional model - Custom Design\n")
+ callNextMethod();
+})
+
+setMethod("show",signature("WfmInfEffects"),function(object)
+{
+ cat("Inference object from wavelet based functional model - Effects\n")
+ cat("Time Effect if Time Design, Circadian Effect if Circadian Design\n")
+ callNextMethod();
+})
+
+setMethod("show",signature("WfmInfMeans"),function(object)
+{
+ cat("Inference object from wavelet based functional model - Groupwise Means\n")
+ callNextMethod();
+})
+
+setMethod("show",signature("WfmInfOverallMean"),function(object)
+{
+ cat("Inference object from wavelet based functional model - Overall Mean\n")
+ callNextMethod();
+})
+
Deleted: pkg/waveTiling/TODO
===================================================================
--- pkg/waveTiling/TODO 2012-05-21 12:54:45 UTC (rev 42)
+++ pkg/waveTiling/TODO 2012-06-04 12:28:07 UTC (rev 43)
@@ -1,12 +0,0 @@
-
-1) plot function: create argument contrasts to plot specific contrasts
-2) wfm.inference: check calculation of confidence intervals
-3) plot: introduce option to give groupnames to put on Y-axis of plot
-4) getNonAnnotatedRegions: make function more generic wrt how the strands and features in the annotation file are defined
-5) getNonAnnotatedRegions: include option to give max./mean expression / FC per region
-6) getNonAnnotatedRegions: add option to set threshold for minimum density of probes within the regions
-7) getSigGenes: add option to include thresholds (density, outliers,...)
-8) getSigGenes: add option to also detect regions that don't overlap but are in the neighbourhood of annotated regions
-9) getSigGenes: add option to include mean / max expression / FC for each gene
-
-
Modified: pkg/waveTiling/inst/scripts/waveTiling-vignette.Rnw
===================================================================
--- pkg/waveTiling/inst/scripts/waveTiling-vignette.Rnw 2012-05-21 12:54:45 UTC (rev 42)
+++ pkg/waveTiling/inst/scripts/waveTiling-vignette.Rnw 2012-06-04 12:28:07 UTC (rev 43)
@@ -46,14 +46,13 @@
\section{Read in and prepare data for analysis}
-First we have to load the \waveTiling{} package and the \Rpackage{waveTilingData} package. The latter contains an \Rclass{TilingFeatureSet} (\textsf{leafdev}) from the \Rpackage{oligoClasses} package \cite{carvalho} with the expression values for the leaf development experiment, and the TAIR 9 \textit{Arabidopsis thaliana} gene identifier data \textsf{tair9gff}. Make sure to also load the \Rpackage{pd.atdschip.tiling} package which contains the tiling array info to map the probe locations on the array to the exact genomic positions. The \Rpackage{pd.atdschip.tiling} package was created by using the \Rpackage{pdInfoBuilder} package \cite{falcon}, which should also be used to build similar packages for other array designs.
+First we have to load the \waveTiling{} package and the \Rpackage{waveTilingData} package. The latter contains an \Rclass{TilingFeatureSet} (\textsf{leafdev}) from the \Rpackage{oligoClasses} package \cite{carvalho} with the expression values for the leaf development experiment. Make sure to also load the \Rpackage{pd.atdschip.tiling} package which contains the tiling array info to map the probe locations on the array to the exact genomic positions. The \Rpackage{pd.atdschip.tiling} package was created by using the \Rpackage{pdInfoBuilder} package \cite{falcon}, which should also be used to build similar packages for other array designs.
<<>>=
library(waveTiling)
library(waveTilingData)
library(pd.atdschip.tiling)
data(leafdev)
-data(tair9gff)
@
We first change the class to \Rclass{WaveTilingFeatureSet}, which is used as input for the wavelet-based transcriptome analysis, and add the phenotypic data for this experiment.
@@ -66,6 +65,11 @@
leafdev
@
+<<echo=FALSE>>=
+rm(leafdev)
+gc()
+@
+
Before starting the transcriptome analysis, the probes that map to several genomic locations (either PM or MM, or forward and reverse strand) are filtered using \Rfunction{filterOverlap}. This function can also be used if the probes have to be remapped to another version of the genome sequence as the version used for the array design. For instance, the probes on the AGRONOMICS1 array are build based on the TAIR 8 genome, and remapped onto the TAIR 9 sequence. The function needs an argument \Rfunarg{BSgenomeObject} available from loading the appropriate \Rpackage{BSgenome} package \cite{pagesBSgenome}. The output is an object of class \Rclass{mapFilterProbe}. After filtering and/or remapping, the expression data are background-corrected and quantile-normalized (\Rfunction{bgCorrQn}). The \Rclass{mapFilterProbe} \Robject{leafdevMapAndFilterTAIR9} is used to make sure only the filtered probes are used in the background correction and normalization step.
<<>>=
@@ -76,7 +80,7 @@
# strand="both",MM=FALSE)
data(leafdevMapAndFilterTAIR9)
-leafdevBQ <- bgCorrQn(leafdev,useMapFilter=leafdevMapAndFilterTAIR9)
+# leafdevBQ <- bgCorrQn(leafdev,useMapFilter=leafdevMapAndFilterTAIR9)
@
\section{Wavelet-based transcriptome analysis}
@@ -85,13 +89,14 @@
The analysis has to be conducted in a chromosome- and strand-wise manner. First, the wavelet-based model is fitted to the expression data, leading to a \Rclass{WfmFit}-class object \textsf{leafdevFit}.
<<>>=
+data(leafdevBQ)
chromosome <- 1
strand <- "forward"
leafdevFit <- wfm.fit(leafdevBQ,filter.overlap=leafdevMapAndFilterTAIR9,
design="time",n.levels=10,
- chromosome=chromosome,strand=strand,
- var.eps="marg",prior="improper",skiplevels=1,
- save.obs="plot",trace=TRUE)
+ chromosome=chromosome,strand=strand,minPos=22000000,
+ maxPos=24000000,var.eps="marg",prior="improper",
+ skiplevels=1,save.obs="plot",trace=TRUE)
leafdevFit
@
@@ -118,13 +123,15 @@
The \Rfunction{getGenomicRegions} accessor outputs a \Rcode{list} of \Rcode{IRanges} objects \cite{pagesIRanges} denoting the start and end position of each significant region. The first element in the \Rcode{list} always gives the significant regions for the mean expression over all arrays (transcript discovery). Elements 2 to 16 in \Robject{sigGenomeRegions} give the differentially expressed regions between any pair of contrasts between different time points. The order is always 2-1, 3-1, 3-2, 4-1,... Hence, \Rcode{sigGenomeRegions[[2]]} gives the differentially expressed regions between time point 2 and time point 1.
-If an annotation file containing gene identifiers is available, we can extract both significantly affected genes with \Rfunction{getSigGenes}, and the non-annotated regions with \Rfunction{getNonAnnotatedRegions}. Both functions output a \Rclass{list} of \Rclass{GRanges} objects \cite{aboyoun}.
+If information on the annotation of the studied organism is available, we can extract both significantly affected genes with \Rfunction{getSigGenes}, and the non-annotated regions with \Rfunction{getNonAnnotatedRegions}. Both functions output a \Rclass{list} of \Rclass{GRanges} objects \cite{aboyoun}. The annotation info can be obtained from ab appropriate object of class TranscriptDb representing an annotation database generated from BioMart. For the current data we make use of the \Rpackage{TxDb.Athaliana.BioMart.plantsmart12} package \cite{carlson}.
<<>>=
-sigGenesCompare <- getSigGenes(fit=leafdevFit,inf=leafdevInfCompare,annoFile=tair9gff)
+library(TxDb.Athaliana.BioMart.plantsmart12)
+sigGenesCompare <- getSigGenes(fit=leafdevFit,inf=leafdevInfCompare,
+ biomartObj=TxDb.Athaliana.BioMart.plantsmart12)
head(sigGenesCompare[[2]])
nonAnnoCompare <- getNonAnnotatedRegions(fit=leafdevFit,inf=leafdevInfCompare,
- annoFile=tair9gff)
+ biomartObj=TxDb.Athaliana.BioMart.plantsmart12)
head(nonAnnoCompare[[2]])
@
@@ -150,8 +157,9 @@
leafdevFitCustom <- wfm.fit(leafdevBQ,filter.overlap=leafdevMapAndFilterTAIR9,
design="custom",design.matrix=custDes,n.levels=10,
- chromosome=chromosome,strand=strand,var.eps="marg",
- prior="improper",skiplevels=1,save.obs="plot",trace=TRUE)
+ 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)
@@ -174,30 +182,30 @@
Plots can be made very easily using the \Rfunction{plotWfm} function which needs both the \Rclass{WfmFit}- and \Rclass{WfmInf}-class objects as input. It also needs an appropriate annotation file. The plot function makes use of the implementations in the \textit{GenomeGraphs}-package \cite{durinck}.
<<fig=TRUE>>=
-gene1 <- tair9gff[tair9gff$ID %in% "AT1G04350",]
-start <- gene1$start-2500
-end <- gene1$end+2500
+trs <- transcripts(TxDb.Athaliana.BioMart.plantsmart12)
+sel <- trs[elementMetadata(trs)$tx_name %in% "AT1G62500.1",]
+start <- start(ranges(sel))-2000
+end <- end(ranges(sel))+2000
plotWfm(fit=leafdevFit,inf=leafdevInfCompare,
- annoFile=gene1,minPos=start,maxPos=end,
- two.strand=TRUE,plotData=TRUE,
- plotMean=FALSE,tracks=c(1,2,6,10,11))
+ biomartObj=TxDb.Athaliana.BioMart.plantsmart12,
+ minPos=start,maxPos=end,two.strand=TRUE,
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/wavetiling -r 43
More information about the Wavetiling-commits
mailing list