Hi

I am using DiffBind for analyses of my Chip-seq results. I was wondering if there is any way I can get the original peak fold enrichment scores (from MACS2 output) for the consensus peak set that DiffBind generates.

Thanks

Hi, two quick comments:

The problem isn't uniquely defined, because more than one peak in the original peak file may correspond with a consensus peak.  (For example, in sample A, peaks 1 and 2 are near each other; in sample B, peak 3 overlaps both 1 and 2: the merged peak will include all 3 peaks.)

Having said that, DiffBind doesn't actually store the fold change, only the score, so notwithstanding the above, there won't be a direct way.  The following code snippet gives you the peaks in the original peak file corresponding with the consensus peaks, with the fold change, which is not what you want, but might be of some small help:

library(IRanges)
library(GenomicRanges)
library(GenomicAlignments)
library(DiffBind)

macs2granges <- function(filename) {
  df <- read.table(filename,blank.lines.skip=TRUE,header=TRUE,comment.char="#",
                   sep="\t",stringsAsFactors=FALSE)
  gr <- GRanges(seqnames=Rle(df$chr),ranges=IRanges(df$start,df$end),
                strand=Rle("*",nrow(df)),score=df$fold_enrichment)
  return(gr)
}

dba2granges <- function(d) {
  gr <- GRanges(seqnames=Rle(d$chrmap[d$vectors$CHR]),
                ranges=IRanges(d$vectors$START,d$vectors$END),
                strand=Rle("*",nrow(d$vectors)))
  return(gr)
}

orig <- macs2granges('your_peaks.xls')
dba_obj <- dba(sampleSheet='your_sampleSheet.csv')

overlap <- subsetByOverlaps(orig,dba2granges(dba_obj))

You can get the full matrix of scores, as a data frame,  as follows:

> scoreMatrix <- dba.peakset(DBA, bRetrieve=TRUE, DataType=DBA_DATA_FRAME)

Before counting, the scores are derived from the peak scores in the peak files (in your case MACS2), normalized to a 0..1 scale. After counting via dba.count(), the scores are computed  The default scoring method, DBA_SCORE_TMM_MINUS_FULL, uses the TMM normalization method from the edgeR package, on read counts minus the control read counts (set to a minimum value of 1), and based on the full library size (number of reads in the bam file). Other scores include raw reads, RPKM, and summit heights, as described in the man page for dba.count for the score parameter.

You can re-set the scores without re-counting (and re-export the matrix) by calling dba.count with peaks=NULL, e.g.

> DBA <- dba.count(DBA, peaks=NULL, score=DBA_SCORE_RPKM)
> readCounts <- dba.peakset(DBA, bRetrieve=TRUE, DataType=DBA_DATA_FRAME)

Hope this helps-

Rory