[Bioc-devel] Issue on MEDIPS.CpGenrich function

[Welliton Souza]

I found a problem in MEDIPS.CpGenrich function (development version).
When I run the example of this function I got the following error:

Error in seqlengths(seqinfo(x)) :
  error in evaluating the argument 'x' in selecting a method for
function 'seqlengths': Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘seqinfo’ for
signature ‘"numeric"’

I think the problem is related to the calls of getGRange and
getPairedGRange functions. I tried to resolve this problem by adding
the missing parameter `dataset` to these function calls.

Diff file attached.
-------------- next part --------------
diff --git a/R/MEDIPS.CpGenrich.R b/R/MEDIPS.CpGenrich.R
index 1f52035..25ba867 100755
--- a/R/MEDIPS.CpGenrich.R
+++ b/R/MEDIPS.CpGenrich.R
@@ -14,69 +14,70 @@ MEDIPS.CpGenrich <-function(file=NULL, BSgenome=NULL, extend=0, shift=0, uniq=1e
 
 	## Proof correctness....
 	if(is.null(BSgenome)){stop("Must specify a BSgenome library.")}
-	
-	## Read region file		
+
+	## Read region file
 	fileName=unlist(strsplit(file, "/"))[length(unlist(strsplit(file, "/")))]
-	path=paste(unlist(strsplit(file, "/"))[1:(length(unlist(strsplit(file, "/"))))-1], collapse="/") 
-	if(path==""){path=getwd()}		
-	if(!fileName%in%dir(path)){stop(paste("File", fileName, " not found in", path, sep =" "))}	
-	
-	if(!paired){GRange.Reads = getGRange(fileName, path, extend, shift, chr.select, uniq)}
-	else{GRange.Reads = getPairedGRange(fileName, path, extend, shift, chr.select, uniq)}
-	
+	path=paste(unlist(strsplit(file, "/"))[1:(length(unlist(strsplit(file, "/"))))-1], collapse="/")
+	if(path==""){path=getwd()}
+	if(!fileName%in%dir(path)){stop(paste("File", fileName, " not found in", path, sep =" "))}
+
+	## Get chromosome lengths for all chromosomes within data set.
+	cat(paste("Loading chromosome lengths for ",BSgenome, "...\n", sep=""))
+	dataset=get(ls(paste("package:", BSgenome, sep="")))
+
+	if(!paired){GRange.Reads = getGRange(fileName, path, extend, shift, chr.select, dataset, uniq)}
+	else{GRange.Reads = getPairedGRange(fileName, path, extend, shift, chr.select, dataset, uniq)}
+
 	## Sort chromosomes
 	if(length(unique(seqlevels(GRange.Reads)))>1){chromosomes=mixedsort(unique(seqlevels(GRange.Reads)))}
 	if(length(unique(seqlevels(GRange.Reads)))==1){chromosomes=unique(seqlevels(GRange.Reads))}
-	
-	## Get chromosome lengths for all chromosomes within data set.
-	cat(paste("Loading chromosome lengths for ",BSgenome, "...\n", sep=""))		
-	dataset=get(ls(paste("package:", BSgenome, sep="")))
+
 	#chr_lengths=as.numeric(sapply(chromosomes, function(x){as.numeric(length(dataset[[x]]))}))
 	chr_lengths=as.numeric(seqlengths(dataset)[chromosomes])
 
 	ranges(GRange.Reads) <- restrict(ranges(GRange.Reads),+1)
-	
+
 	##Calculate CpG density for regions
 	total=length(chromosomes)
-	cat("Calculating CpG density for given regions...\n")  
-	seq=matrix(unlist(IRanges::lapply(RangedData(GRange.Reads),function(x){		
+	cat("Calculating CpG density for given regions...\n")
+	seq=matrix(unlist(IRanges::lapply(RangedData(GRange.Reads),function(x){
 		i=which(mixedsort(chromosomes)%in%names(x) )
 		ranges(x)<-restrict(ranges(x),end=chr_lengths[which(chromosomes %in% names(x))])
 		y=DNAStringSet(getSeq(dataset, names=space(x), start=start(x), end=end(x), as.character=TRUE))
 		c(sum(as.numeric(vcountPattern("CG",y))),sum(as.numeric(vcountPattern("C",y))),sum(as.numeric(vcountPattern("G",y))),sum(as.numeric(width(y))),length(y))
 		}
 	),use.names=F),ncol=5,nrow=total,byrow=T)
-  
+
 	Value=colSums(seq)
 	unused=length(GRange.Reads)-Value[5]
 	if ( unused!=0 )cat(unused,"unused sequences, limits out of range\n")
-	
+
 	regions.CG=Value[1]
 	regions.C=Value[2]
 	regions.G=Value[3]
 	all.genomic=Value[4]
-	
+
 	regions.relH=as.numeric(regions.CG)/as.numeric(all.genomic)*100
-	regions.GoGe=(as.numeric(regions.CG)*as.numeric(all.genomic))/(as.numeric(regions.C)*as.numeric(regions.G))  
-	
+	regions.GoGe=(as.numeric(regions.CG)*as.numeric(all.genomic))/(as.numeric(regions.C)*as.numeric(regions.G))
+
 	##Calculate CpG density for reference genome
-	cat(paste("Calculating CpG density for the reference genome", BSgenome, "...\n", sep = " "))	
+	cat(paste("Calculating CpG density for the reference genome", BSgenome, "...\n", sep = " "))
 	CG <- DNAStringSet("CG")
 	pdict0 <- PDict(CG)
 	params <- new("BSParams", X = dataset, FUN = countPDict, simplify = TRUE, exclude = c("rand", "chrUn"))
-	genome.CG=sum(bsapply(params, pdict = pdict0))			
+	genome.CG=sum(bsapply(params, pdict = pdict0))
 	params <- new("BSParams", X = dataset, FUN = alphabetFrequency, exclude = c("rand", "chrUn"), simplify=TRUE)
 	alphabet=bsapply(params)
 	genome.l=sum(as.numeric(alphabet))
- 	
+
 	genome.C=as.numeric(sum(alphabet[2,]))
 	genome.G=as.numeric(sum(alphabet[3,]))
 	genome.relH=genome.CG/genome.l*100
 	genome.GoGe=(genome.CG*genome.l)/(genome.C*genome.G);
-	
-	enrichment.score.relH=regions.relH/genome.relH	
-	enrichment.score.GoGe=regions.GoGe/genome.GoGe	
- 	
+
+	enrichment.score.relH=regions.relH/genome.relH
+	enrichment.score.GoGe=regions.GoGe/genome.GoGe
+
 	gc()
 	return(list(genome=BSgenome, regions.CG=regions.CG, regions.C=regions.C, regions.G=regions.G, regions.relH=regions.relH, regions.GoGe=regions.GoGe, genome.C=genome.C, genome.G=genome.G, genome.CG=genome.CG, genome.relH=genome.relH, genome.GoGe=genome.GoGe, enrichment.score.relH=enrichment.score.relH, enrichment.score.GoGe=enrichment.score.GoGe))
 }