Dear Carlos,

The code looks fine except the following line. You would want to use MyPeakList instead of the first 6 records as MyPeakList[1:6,].

AnnotatedPeak = annotatePeakInBatch(MyPeakList[1:6,], AnnotationData=annotation.cds)

Please let me know if this resolves the issue. Thanks!

Best regards,

Julie


On 3/1/12 12:13 PM, "Carlos Araya" <claraya@stanford.edu> wrote:

Dear Julie,

I'm ready to perform GO enrichment analysis on ~292 ChIP-seq datasets and would like to use your ChIPpeakAnno package. However, I've encountered problems that I believe are due to differences in the gene/feature names used. Here is an excerpt of my code:

library(ChIPpeakAnno)
library(org.Ce.eg.db)

# load peak calls and genomic annotations, filter for gene (CDS) annotations:
peaks.bed = read.table("YL474_HPL-2_L1_yale_stn_1x2_peaks.bed", header=FALSE)
annotation = read.table("in2shape_wbComp.bed", header=TRUE)
annotation.cds = annotation[annotation$feature.type == "cds",]
summary(annotation.pct)
chrm          start               end                     feature      score     strand         name                   gene
 chrI  :4593   Min.   :      58   Min.   :     779   WBGene00000001:    0   .:30163   -:14976   2L52.1 :    0   WBGene00000001:    0
 chrII :5195   1st Qu.: 4687094   1st Qu.: 4688288   WBGene00000002:    0             .:    0   2RSSE.1:    0   WBGene00000002:    0
 chrIII:4345   Median : 8229719   Median : 8235458   WBGene00000003:    0             +:15187   2RSSE.2:    0   WBGene00000003:    0
 chrIV :5126   Mean   : 8415951   Mean   : 8419438   WBGene00000004:    0                       4R79.2a:    0   WBGene00000004:    0
 chrM  :  12   3rd Qu.:11924503   3rd Qu.:11929206   WBGene00000005:    0                       4R79.2b:    0   WBGene00000005:    0
 chrV  :6763   Max.   :20915455   Max.   :20922707   (Other)       :    0                       (Other):    0   (Other)       :    0
 chrX  :4129                                         NA's          :30163                       NA's   :30163   NA's          :30163

# convert bed tables into ranged data:
MyPeakList = BED2RangedData(peaks.bed)
annotation.cds = BED2RangedData(annotation.cds)

# annotate peaks:
AnnotatedPeak = annotatePeakInBatch(MyPeakList[1:6,], AnnotationData=annotation.cds)
as.data.frame(AnnotatedPeak)
 space  start    end width             names        peak strand feature start_position end_position insideFeature distancetoFeature shortestDistance
1     I  16921  17165   245 MACS_peak_2 00006 MACS_peak_2      +   00006          11618        16804    downstream              5303              117
2     I 110632 110850   219 MACS_peak_4 00018 MACS_peak_4      +   00018         111036       112316      upstream              -404              186
3     I 536767 536967   201 MACS_peak_6 00070 MACS_peak_6      +   00070         537125       541634      upstream              -358              158
4     I   3680   4064   385 MACS_peak_1 02345 MACS_peak_1      -   02345           4221        10230    downstream              6550              157
5     I 108389 108741   353 MACS_peak_3 02354 MACS_peak_3      -   02354         101213       108161      upstream              -228              228
6     I 315156 315412   257 MACS_peak_5 02372 MACS_peak_5      -   02372         310981       314992      upstream              -164              164
  fromOverlappingOrNearest
1             NearestStart
2             NearestStart
3             NearestStart
4             NearestStart
5             NearestStart
6             NearestStart

# execute enrichment analysis:
enrichedGO = getEnrichedGO(AnnotatedPeak, orgAnn = "org.Ce.eg.db", maxP = 0.01, multiAdj = TRUE, minGOterm = 10, multiAdjMethod = "BH" )
Error in if (class(go.ids) != "matrix" | dim(go.ids)[2] < 4) { :
  argument is of length zero

I've highlighted in blue the outputs of distinct commands and in red the error that I obtain upon attempting GO analysis. Any ideas as to how I should annotate my features (which codes I should use) or how to prepare the objects in R for the enrichment analysis? Please note that in my annotation BED file the WormBase identifiers for each gene are also available (under the "gene" column). I noticed there is a "org.Ce.egWORMBASE" object in the "org.Ce.eg.db" library which may be useful, but I do not understand how to use it convert between gene IDs.

Kind regards,

Carlos L. Araya, PhD
Stanford University | Department of Genetics
300 Pasteur Dr. #M348 | Stanford, CA 94305
m: 206.661.3648 | e: claraya@stanford.edu



	[[alternative HTML version deleted]]

