Diffbind adding my own consensus peaks
1
1
Entering edit mode
8 months ago
lara.bedeyan ▴ 10

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.

dba.peakset ATAC-Seq Diffbind ChIP-Seq • 584 views
1
Entering edit mode
8 months ago
Rory Stark ★ 1.5k

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.