Hello all.
I am attempting to use "featureAlignedDistribution" in the "ChIPpeakAnno" package to get the binding curves for some subset of annotated peaks I have stored as GRanges object.
My experimental test function is as follows:
bindingPattern<-function(gr, bwpath, outprefix, ups=2000, downs=2000) {
wid <- width(gr)
gr_recentered <- gr_center <- gr
start(gr_center) <- start(gr) + floor(wid/2)
width(gr_center) <- 1
start(gr_recentered) <- start(gr_center) - ups
end(gr_recentered) <- end(gr_center) + downs
sig <- import(bwpath, format="BigWig", which=gr_recentered, as="RleList")
aligned_signal <- featureAlignedSignal(sig, gr_center, upstream=ups, downstream=downs)
distout <- paste(outprefix, "distribution.png", sep="")
checkMakeFile(distout)
png(filename=distout)
featureAlignedDistribution(aligned_signal, gr_center, upstream=ups, downstream=downs, type="l")
dev.off()
}
(Note that bwpath is the filepath to a bigwig file I have made from a bedgraph file which in turn was generated from MACS2 peakcalling. Both show up fine in IGV)
However, featureAlignedSignal() fails with the error:
Error in deparse(sig) : object 'sig' not found
I thought perhaps sig was null or the wrong type, so after the sig<-import... line I checked it out:
1.class(sig) shows that sig is of type RleSimpleList, which is one of the accepted object types for featureAlignedSignal
2.Printing it out shows: RleList of length 17 $chr1 numeric-Rle of length 258084 with 80751 runs Lengths: 39199 1 1 1 1 ... 1 2 3 1 34199 Values : 0 67.9800796508789 69.5986633300781 67.9800796508789 66.3615036010742 ... 6.474289894104 4.85572004318237 3.23714995384216 1.61856997013092 0 $chr10 numeric-Rle of length 766161 with 310066 runs Lengths: 32592 1 1 2 2 ... 29 4 3 7 2282 Values : 0 51.7943496704102 50.1757698059082 48.557201385498 46.9386291503906 ... 9.71144008636475 8.09286975860596 6.474289894104 4.85572004318237 0
...etc...
I tried looking at the source code for featureAlignedSignal() but I did not see any obvious reason why this error should occur. Intestingly, I found that if I convert sig to a list:
sig <- list(sig)
then featureAlignedSignal() runs and sets aligned_signal as a list object, which prints as:
[1] "Signal_aligned class after featureAlignedSignal: list"
[[1]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,] 90.721035 126.2082518 150.6082455 133.289509 114.9996338 103.952873 90.478250 98.6115808 94.241433 92.01589451 82.9518829 78.5412695 93.391681 118.7628157 124.5087505
[2,] 213.489816 227.0858353 241.6529968 241.491138 249.9886482 254.601581 261.601912 261.1163364 248.774718 227.69279861 217.9408966 208.3103859 225.507725 229.5541584 225.4672607
...etc...
but that leads to a problem in featureAlignedDistribution():
Error in legend("topright", legend = colnames(density), col = col, lty = lty, :
'legend' is of length 0
In addition: Warning message:
In featureAlignedDistribution(aligned_signal, gr_center, upstream = ups, :
cvglists contain NA values. NA value will be omit.
I imagine it's because featureAlignedSignal wasn't meant to take a 'list' object in the first place, as its source code shows that it only wants RleList, SimpleRleList, or CompressedRleList.
I'm stumped. Can anyone help me out? Your time and help is much appreciated!
(Just a little extra info: Running on Ubuntu. R version is 3.3.1-1xenial, chippeakanno version is 3.6.5. Using the latest version of bioconductor)
