Question: Problem with finding overlapping chips peaks with ChIPpeakAnno
0
gravatar for Debbie
3.8 years ago by
Debbie10
United States
Debbie10 wrote:

I have 2 datasets 1. control and other is 2. RA  and I ma trying to find overlapping chip seq peaks in theses 2 datasets with the script below:

#To read peak tables#
macspeaks<- read.table("MYCN_ControlPeaks.txt", header=TRUE, sep="\t", stringsAsFactors=FALSE)
head (macspeaks)
macspeaks.RA<- read.table("MYCN_RAPeaks.txt", header=TRUE, sep="\t", stringsAsFactors=FALSE)
head (macspeaks.RA)
--------------------
# function to convert table into a IRanges#

Control.peaks = RangedData(IRanges(start=macspeaks$Start, end=macspeaks$End, names=paste(macspeaks$Chr,macspeaks$Start,sep=":")),strand="*")
head (Control.peaks)

RA.peaks = RangedData(IRanges(start=macspeaks.RA$Start, end=macspeaks.RA$End, names=paste(macspeaks.RA$Chr,macspeaks.RA$Start,sep=":")),strand="*")
head (RA.peaks)
-----------------------
#to find overlapping peaks in control and RA#

t1 =findOverlapsOfPeaks (Control.peaks, RA.peaks, maxgap=0L, minoverlap=1L, NameOfPeaks="Control","RA", select="all", annotate=1)
r = t1$OverlappingPeaks
pie(table(r$overlapFeature))
as.data.frame(t1$MergedPeaks)

-------------------------

While running it I am getting the following error in the find overlapping peaks

+ t1 =findOverlapsOfPeaks (Control.peaks, RA.peaks, maxgap=0L, minoverlap=1L, NameOfPeaks="Control","RA", select="all", annotate=1)
Error in findOverlapsOfPeaks(Control.peaks, RA.peaks, maxgap = 0L, minoverlap = 1L,  : 
  The length of input peaks list should no more than 5

> r = t1$OverlappingPeaks
Error in t1$OverlappingPeaks : 
  object of type 'closure' is not subsettable
> pie(table(r$overlapFeature))
Error in eval(expr, envir, enclos) : object 'r' not found
> as.data.frame(t1$MergedPeaks)
Error in as.data.frame(t1$MergedPeaks) : 
  error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': Error in t1$MergedPeaks : object of type 'closure' is not subsettable

I am still in the learning stage of R and have very limited knowledge. 

Kindly help to solve the problem.

Thank you

chip-seq R • 1.4k views
ADD COMMENTlink modified 3.8 years ago by seidel6.8k • written 3.8 years ago by Debbie10
0
gravatar for seidel
3.8 years ago by
seidel6.8k
United States
seidel6.8k wrote:

One problem is that you have an error in your function call: NameOfPeaks="Control","RA"

This will cause an error because the comma separates the arguments to the function, so the "RA" string is all alone as an argument. If you want to hand NameOfPeaks more than one thing you have to use c(), as in NameOfPeaks=c("Control","RA"). Either way, I don't know what purpose this argument serves for this function. To debug your problem, two approaches you can take are (1) it can be helpful to create some toy data, to see where things are failing, and (2) simplify your function call to the bare minimum. 

# load the library
library(ChIPpeakAnno)

# create 5 random segments on an interval
# that are at least 50 bases wide
peakStarts <- sample(1:5000,5)
peakEnds <- peakStarts + 50 + sample(1:300,1)

# populate Control with fake data
Control.peaks = RangedData(IRanges(start=peakStarts,
                  end=peakEnds,
                  names=paste(rep("Chr1",ntimes=length(peakStarts)),peakStarts,sep=":")),
                  strand="*")

# get 5 new random locations
peakStarts <- sample(1:5000,5)
peakEnds <- peakStarts + 50 + sample(1:300,1)

# populate the RA with fake data
RA.peaks = RangedData(IRanges(start=peakStarts,
             end=peakEnds,
             names=paste(rep("Chr1",ntimes=length(peakStarts)),peakStarts,sep=":")),
             strand="*")

# check to see that we have a legit data structure
head (RA.peaks)

# RangedData with 5 rows and 1 value column across 1 space
#            space       ranges |      strand
#          <factor>    <IRanges> | <character>
# Chr1:1675        1 [1675, 1802] |           *
# Chr1:4661        1 [4661, 4788] |           *
# Chr1:3516        1 [3516, 3643] |           *
# Chr1:1959        1 [1959, 2086] |           *
# Chr1:3031        1 [3031, 3158] |           *

# find overlaps
t1 <- findOverlapsOfPeaks(Control.peaks, RA.peaks, maxgap=0, minoverlap=1)

You might check the help for ?findOverlapsOfPeaks

ADD COMMENTlink written 3.8 years ago by seidel6.8k

Thanks will try and get back to you.

Thanks Debbie 

ADD REPLYlink written 3.8 years ago by Debbie10
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 906 users visited in the last hour