Question

Diffbind: retrieve raw peaks per sample

0

Entering edit mode

ZheFrench ▴ 60

@zhefrench-11689

Last seen 12 weeks ago

France

Trying to retrieve my raw peaks . Here trying to get peaks called by macs2 for sample L13S13

I should have 94888 but I get 269818 . Seems to be the number of consensus when taking the whole object.

show(dba_chip_init)

         ID Tissue  Factor Condition Treatment Replicate Caller Intervals
1      L13S13  MCF10   K27AC        T1 Tamoxifen         1   macs     94888
2     L142S14  MCF10   K27AC        T1 Tamoxifen         3   macs    100680
3      L17S17  MCF10   K27AC        T7 Tamoxifen         1   macs    106629
....
19   T1_INPUT  MCF10 Control        T1 Tamoxifen        c1           269818
20   T7_INPUT  MCF10 Control        T7 Tamoxifen        c2           269818
21 unT7_INPUT  MCF10 Control      unT7 unTreated        c3           269818

raw_peak <- dba.peakset(dba_chip_init,sampID="L13S13",consensus=F,bRetrieve=T,DataType=DBA_DATA_GRANGES,minOverlap=0.99)
export.bed(raw_peak, paste0(c(dir_output,"/peaksPerSampleInBed/",opt$output,".peaks.bed"),collapse=""))

EDIT : Or just how can can get the path of raw files to then retrieve via dba.peasket(dbaObject=NULL,peaks=pathToFile,bRetrieve=T) juste the peaks called with any consensus added.

Thnaks

diffbind • 1.1k views

ADD COMMENT • link updated 7.5 years ago by Rory Stark ★ 5.2k • written 7.5 years ago by ZheFrench ▴ 60

score 0 · Answer 1 · 2017-06-19

0

Entering edit mode

Rory Stark ★ 5.2k

@rory-stark-5741

Last seen 5 weeks ago

Cambridge, UK

For retrieving a specific peakset, you need to specify the peakset number if you want just the one:

raw_peak <- dba.peakset(dba_chip_init,peaks=1,bRetrieve=T,DataType=DBA_DATA_GRANGES)

The sampID parameter is inoperative in this case.

-Rory

ADD COMMENT • link 7.5 years ago Rory Stark ★ 5.2k

0

Entering edit mode

Got the following error :

Error in if (numCols > 0) { : argument is of length zero
Calls: dba.peakset -> pv.writePeakset -> pv.do_peaks2bed

UPDATE : Also tried another strategy :

#' Load dba object (from ChIPQC) : 
dba_chip_init <- dba.load(file=filename,dir= dir_output, pre='', ext='RData')
dba.show(dba_chip_init)
expPeaks = peaks(dba_chip_init)
length(expPeaks)
expPeaks$L13S13

I get :

GRanges object with 269262 ranges and 2 metadata columns:
         seqnames               ranges strand |    Counts bedRangesSummitsTemp
            <Rle>            <IRanges>  <Rle> | <integer>            <numeric>
       1     chr1       [ 9941, 10211]      * |       205                10066
       2     chr1       [13473, 14854]      * |        47                14287
       3     chr1       [21149, 22969]      * |        70                21696
       4     chr1       [23466, 24103]      * |        25                23477
       5     chr1       [25383, 66079]      * |      1052                29093
     ...      ...                  ...    ... .       ...                  ...

whereas I should retrieve 94888 raw peaks for this sample:

           ID Tissue  Factor Condition Treatment Replicate Caller Intervals
1      L13S13  MCF10   K27AC        T1 Tamoxifen         1   macs     94888

ADD REPLY • link 7.5 years ago ZheFrench ▴ 60

0

Entering edit mode

I'm not sure I understand your code -- what is the call to peaks() on the third line?

-R

ADD REPLY • link 7.4 years ago Rory Stark ★ 5.2k