Hello,

I have some ATAC-seq data and I did peak calling with MACS2. Now I want to extract a matrix with raw counts for the called peaks.

My sample sheet looks like this:

SampleID Condition bamReads Peaks PeakCaller mut1 mut ATAC_peaks/mapped/mut1.bam ATAC_peaks/peaks/mut1_peaks.xls macs mut2 mut ATAC_peaks/mapped/mut2.bam ATAC_peaks/peaks/mut2_peaks.xls macs wt1 wt ATAC_peaks/mapped/wt1.bam ATAC_peaks/peaks/wt1_peaks.xls macs wt2 wt ATAC_peaks/mapped/wt2.bam ATAC_peaks/peaks/wt2_peaks.xls macs

The data comes from different samples with different conditions, but sadly it is not possible to have replicates. I know it is not ideal.

Now I try to create a DBA object and then extract the raw counts as follows:

samples <- read.xlsx("ATAC_peaks/metadata_atac.xlsx")
peaks  <- dba(sampleSheet=samples)
black <- dba.blacklist(peaks, blacklist=TRUE, greylist=FALSE)

rawcounts <- dba.count(black, peaks=NULL, score=DBA_SCORE_READS)

counts <- dba.peakset(rawcounts , bRetrieve=TRUE, writeFile = "ATAC_counts.csv")

However when I get to dba.count(black, peaks=NULL, score=DBA_SCORE_READS), I get the following error:

Error in dba.count(black, peaks = NULL, score = DBA_SCORE_READS) : 
  DBA object must contain only counts

However, I have no problem when running the same code without peaks=NULL, score=DBA_SCORE_READS

I am using R version 4.3.1 and DiffBind 3.10.1.

Help is very much appreciated.