Consensus Peakset DIffBind

Question

[g.persic1991]

Hi everyone, I'm analyzing H3K4me3 ChIP-seq and i want to find different bound sites using DiffBind. I have Control samples (7 replicates) and Treated (5 replicates)

After peak calling and load a samplesheet, i have tried to make consensus peak for each condition.

consensus <- dba.peakset(dbObj, consensus = DBA_CONDITION, minOverlap=0.3)

dba.show(consensus, consensus$masks$Consensus)

ID Factor Condition Replicate Caller Intervals

13 Control_k4 K4me3 Control_k4 1-2-3-4-5-6-7 raw 46408

14 Severe_k4 K4me3 Severe_k4 1-2-3-4-5 raw 35976

This allow to discard region enriched just in few samples (1 samples in this case) and take in account only "real enriched regions"....right??

Now the question is...how to count reads in this regions using dba.count???

it's conceptually correct or it's wrong?

i'm expected that the regions investigated are 82384 (or less considering that some overlap)...right?

Thanks in advance

Giuseppe

Answer 1

[Rory Stark]

With minOverlap=.3, you will keep peaks called in at least 3 Control samples (peaks only called in 1 or 2 Control samples will not be included), and at least 2 of the Treated samples (peaks only called in 1 Treated sample will not be included). If you only ewant to exclude peaks called in a single replicate for both Control and Treated, use minOverlap=2.

Forming a binding matrix with read counts using this consensus peakset takes three steps:

> consensus <- dba(consensus, consensus$masks$Consensus, minOverlap=1)
> consensus.peaks <- dba.peakset(consensus, bRetrieve=TRUE)
> counts <- dba.count(dbObj, peaks=consensus.peaks)

The first step forms a DBA object containing all of the consensus peaks (both Control and Treated). The second step retrieves the consensus peakset. The third step counts reads in the consensus peakset, using the sample information in the original DBA object.