Diffbind adding my own consensus peaks

Question

[lara.bedeyan]

Hi

I am using Diffbind to determine peaks that are significantly changing between my 2 conditions in my ATAC-Seq experiment. I want to add my own consensus peaklist that I created using bedtools to be able to integrate with my other analyses. Currently my consensus peaklist is in a bed format. looks like this:

chr1    3185152 3185266 +   0

Currently I am importing the list as table

merged_narrowPeak_sorted <- read.table("merged_narrowPeak_sorted.txt",header = FALSE, sep="\t",stringsAsFactors=FALSE, quote="")

and then I am trying to add it as a consensus peaklist using dba.peakset

consensus.peaks <- dba.peakset(NULL,peaks= merged_narrowPeak_sorted, peak.caller = "bed")

I get this error

Error in mergeScores(merged, def, peakset, TRUE) : 
  Not compatible with requested type: [type=character; target=double]

Should I be changing my bed file or should I be using different arguments for the dba.peakset?

I tried playing around with peak caller (tried "narrow") and scoreCol =0. Currently I just want to provide a list of genomic regions I want Diffbind to quantify each of my samples for each of my peaks.

Answer 1

[Rory Stark]

I'm not entirely sure what you are actually trying to do. Are you trying to add this peakset to use as a consensus peakset for counting? If so you can specify it as a parameter to dba.count() by setting peaks=merged_narrowPeak_sorted without having to load it to look like a sample (which it's not). Alternaitvely, if you have a bunch of these files, all with the same merged peakset but different counts pre-computed, you can add them via dba.peakset() using the counts= parameter, with the count values in the fourth columns (e.g. counts=merged_narrowPeak_sorted[,-4]).

If you really do need to create a new DBA object with a consensus peakset as a sample, there are some with dba.peakset(), which I'll sort out soon. In the meantime I can give you some ways to get past this.

First, you should be setting peak.format="bed", since peak.caller only sets a metadata string, nothing more.

However there is a bug right now in doing this. There are a number of alternative workarounds you can try while a bugfix is pending:

• You can pass the bedfile instead of the table: peaks="merged_narrowPeak_sorted.txt"
• You can set scoreCol=5
• You can pass the score in column 4: peaks=merged_narrowPeak_sorted[,-4]

In all of these cases, there is a bug when adding an initial peakset using DBA=NULL. You have to add a line after the first call to dba.peakset(). for example:

consensus.peaks <- dba.peakset(NULL,peaks="merged_narrowPeak_sorted.txt"`, peak.format = "bed")
consensus.peaks <- dba(consensus_peaks)

I'll be addressing the bugs soon in maintenance release.