Search
Question: Diffbind: retrieve raw peaks per sample
0
gravatar for ZheFrench
5 months ago by
ZheFrench10
ZheFrench10 wrote:

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

 

ADD COMMENTlink modified 5 months ago by Rory Stark2.2k • written 5 months ago by ZheFrench10
0
gravatar for Rory Stark
5 months ago by
Rory Stark2.2k
CRUK, Cambridge, UK
Rory Stark2.2k wrote:

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 COMMENTlink written 5 months ago by Rory Stark2.2k

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 REPLYlink modified 5 months ago • written 5 months ago by ZheFrench10

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

-R

ADD REPLYlink written 5 months ago by Rory Stark2.2k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 328 users visited in the last hour