Hi there,
I have a bed file from MACS, I did the motif-enrichment by MEME-ChIP and found the primary motif. Now, I want to find the secondary motif, but I want to resize the bed file to 500bp centered on the primary motif. Is there an easy way to do this? I know one can use resize function for GenomicRanges object, but it can only center on the middle of the ranges.
Thank you,
Ming
2 answers
If you know the location of the primary motif for the features in your bed file, simply use the starts and ends of the motifs to create a new GRanges object, and resize that with a fixed center.
This function returns an aligned GRanges object based on matchPattern or will also work with matchPWM. It's a little slow because it has to retrieve the DNA sequence for each range but I wanted to make sure there was only one site under each peak.
returnAlignedGR <- function(gr,pattern,width) {
# pattern<-"WWWGWWGAANTTCWWCWWW"
c <-as.matrix("chr0",nrow=length(gr))
s <- as.matrix(rep(0,width(gr[1])))
for (i in 1:length(gr)) {
seq<-as.character(seqnames(gr[i])@values)
c[i]<-seq
start<-ranges(gr[i])@start
width2<-ranges(gr[i])@width
peak <- DNAString(Mmusculus[[seq]],start=start,nchar=width2)
site <- matchPattern(pattern,peak,fixed=F)
if(length(site) == 1) {
s[i]<-start+start(site)+8 #assuming motif is 15mer, 8 is center
} else {
s[i]= 0
}
}
gr_outRD<-RangedData(space=c, IRanges(start=s-width, end=s+width ), strand="*")
return(as(gr_outRD,"GRanges"))
}
Log in to answer this question.