Problem with finding overlapping chips peaks with ChIPpeakAnno
1
0
Entering edit mode
8.6 years ago
Debbie ▴ 10

I have 2 datasets 1. control and other is 2. RA and I am 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 • 2.8k views
ADD COMMENT
0
Entering edit mode
8.6 years ago
seidel 11k

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 COMMENT
0
Entering edit mode

Thanks will try and get back to you.

Thanks,
Debbie

ADD REPLY

Login before adding your answer.

Traffic: 1832 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6