Diffbind: retrieve raw peaks per sample
1
0
Entering edit mode
ZheFrench ▴ 40
@zhefrench-11689
Last seen 15 months 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 • 935 views
ADD COMMENT
0
Entering edit mode
Rory Stark ★ 5.1k
@rory-stark-5741
Last seen 9 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
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
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

Login before adding your answer.

Traffic: 1016 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6