How to correct this error More fragment sizes than libraries
in DiffBind Package ? What does that means ?
How to correct this error More fragment sizes than libraries
in DiffBind Package ? What does that means ?
This error occurs in a call to dba.count()
if the fragmentSize
parameter is a vector whose length is greater than the number of libraries (unique BAM files). Either you are specifying a fragmentSize
value directly, or the default value has changed. Check:
> myDBA$config$fragmentSize
By default, this is set to 125.
-Rory
You have identified bugs in both dba.mask()
(not building the correct object) and dba()
(not masking the fragmentSize
vector). I have checked in fixes for these to DiffBind_2.4.5
.
When the fixes are available in the next day or two, you don't need to use dba.mask()
at all here. You can do this in one step directly using dba()
:
> dba_chip <- dba(dba_chip_init, mask=dba_chip_init$masks$K27AC)
If dba_chip_init
is a ChIPQexperiment
object, you can get the masked DBA
object in a couple of ways. Either:
> dba_chip <- dba(dba_chip_init, mask=dba_chip_init@DBA$masks$K27AC)
or
> dba_chip_dba <- dba_chip_init@DBA > dba_chip <- dba(dba_chip_dba, mask=dba_chip_dba@DBA$masks$K27AC)
If you need a workaround immediately, before the fixes are available, you can use:
> dba_chip <- dba(dba_chip_init, mask=dba_chip_init@DBA$masks$K27AC)
> dba_chip$config$fragmentSize <- dba_chip$config$fragmentSize[dba_chip_init@DBA$masks$K27AC]
Cheers-
Rory
For the workaround :
dba_chip <- dba(dba_chip_init, mask=dba_chip_init@DBA$masks$K27AC)In dba(dba_chip_init, mask = dba_chip_init@DBA$masks$K27AC) :
Returning new DBA object (not ChIPQCexperiment object)
dba.show(dba_chip, bContrast=T)
[1] "contrast"
Warning message:
In checkQCobj(DBA$ChIPQCobj, res) :
ChIPQexperiment out of sync -- returning new DBA object
There are warnings but it seems to work.
Update : _ctr removed. dba_chip <- dba(dba_chip_init, mask=dba_chip_init@DBA$masks$K27AC) dba_chip$config$fragmentSize <- dba_chip$config$fragmentSize[dba_chip_init@DBA$masks$K27AC] print("contrast") dba_chip <- dba.contrast(dba_chip,categories=DBA_CONDITION,minMembers=2) dba.show(dba_chip, bContrast=T) dba.show(dba_chip)
Warning message:
In dba(dba_chip_init, mask = dba_chip_init@DBA$masks$K27AC) :
Returning new DBA object (not ChIPQCexperiment object)
[1] "contrast"
Warning message:
In checkQCobj(DBA$ChIPQCobj, res) :
ChIPQexperiment out of sync -- returning new DBA object
Group1 Members1 Group2 Members2
1 T1 2 T7 2
2 T1 2 unT7 2
3 T7 2 unT7 2
....
As I am into it I continue to post my trouble in the same thread.
[1] "analyse for 3 contrasts..."
dba_chip <- dba.analyze(dba_chip,bParallel=T) UPDATE : bParallel=F correct the bug but still have warning at the end :
In checkQCobj(DBA$ChIPQCobj, res) :
ChIPQexperiment out of sync -- returning new DBA object
converting counts to integer mode
converting counts to integer mode
converting counts to integer mode
gene-wise dispersion estimates
gene-wise dispersion estimates
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
mean-dispersion relationship
mean-dispersion relationship
final dispersion estimates
final dispersion estimates
Error in sendMaster(try(lapply(X = S, FUN = FUN, ...), silent = TRUE)) :
long vectors not supported yet: fork.c:376
Calls: dba.analyze ... <Anonymous> -> mclapply -> lapply -> FUN -> sendMaster
Error in sendMaster(try(lapply(X = S, FUN = FUN, ...), silent = TRUE)) :
long vectors not supported yet: fork.c:376
Calls: dba.analyze ... <Anonymous> -> mclapply -> lapply -> FUN -> sendMaster
Error in sendMaster(try(lapply(X = S, FUN = FUN, ...), silent = TRUE)) :
long vectors not supported yet: fork.c:376
Calls: dba.analyze ... <Anonymous> -> mclapply -> lapply -> FUN -> sendMaster
Warning messages:
1: In mclapply(arglist, fn, ..., mc.preschedule = TRUE, mc.allow.recursive = TRUE) :
all scheduled cores encountered errors in user code
2: In drec$counts <- NULL : Coercing LHS to a list
3: In drec$counts <- NULL : Coercing LHS to a list
4: In drec$counts <- NULL : Coercing LHS to a list
5: In checkQCobj(DBA$ChIPQCobj, res) :
ChIPQexperiment out of sync -- returning new DBA object
You will need to update to Bioconductor 3.5
to get the fixes I am making, but I may have a workaround for you.
Basically what is happening relates to issues (mostly warnings) around working directly on the DBA
object created within a ChIPQCexperiment
object. The workaround is to extract a "clean" DBA
object and just work on that in parallel. Here's workaround code that should eliminate the warnings:
> dba_chip_dba <- dba_chip_init@DBA > dba_chip_dba$ChIPQCobj <- NULL > dba_chip <- dba(dba_chip_dba, mask=dba_chip_dba$masks$K27AC)
At this point, the dba_chip
object will be a pure DiffBind
object without any link back to the ChIPQCexperiment
object and it should work fine in the version you are using.
-Rory
Thanks but still :
Error: trying to get slot "DBA" from an object (class "DBA") that is not an S4 object Execution halted
I'm not root on this machine and Bioconductor is 3.4 and can't upgrade more only if I upgrade R itself...That's a bit boring...
But anyway I will go with warnings ...
I sent you a private email with some questions not to overload this thread...Hope you ll take 10 minutes to answer . Thanks !
Also notice something weird thing with :
dba.report(dba_chip,contrast=1,bCalled=TRUE,bCounts=TRUE,bCalledDetail=TRUE,ext='tsv',file=filename,initString="")
If you want to write file not in the current directory you can't pass absolute path via file because even if you set initString empty, it will add "_" to your absolute path...and there is no dir option , so you are force to write in directory where Rscript is launched.
The second line should return the first number of the list , no ?
dba.overlap(dba_chip,mode=DBA_OLAP_RATE) [1] 186450 102988 78376 65514 53238 43346 dba.overlap(dba_chip,mode=DBA_OLAP_RATE,1) [1] 94888
UPDATE : It returns the number of intervals for the sample in line 1 from your dba object. In doc it is said that it should return intervals present at least the minimum value you give. no ?
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
From doc :
MODE: Compute overlap rates at different stringency thresholds: dba.overlap(DBA, mask, mode=DBA_OLAP_RATE, minVal)
Finaly minVal is not valid option because I'm using an old version...and in fact it was passing 1 to mask...
Not sure what you are trying to do with the second line? the value of 1
will be assigned to the second argument, mask
, so this will return the number of peaks in the first sample. It is the same as saying:
dba_first <- dba(dba_chip, mask=1) dba.overlap(dba_first,mode=DBA_OLAP_RATE)
If you just want the first value, try:
dba.overlap(dba_chip,mode=DBA_OLAP_RATE)[1]
-R
Why total number of peaks per condition plotted in VennDiagram is not the same as showed in the dba object ?
dba_chip_consensus <- dba.peakset(dba_chip, consensus=c(DBA_CONDITION),minOverlap=0.99) dba.show(dba_chip_consensus) 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 4 L182S18 MCF10 K27AC T7 Tamoxifen 3 macs 87323 5 L15S15 MCF10 K27AC unT7 unTreated 1 macs 106930 6 L162S16 MCF10 K27AC unT7 unTreated 3 macs 109774 7 T1 MCF10 K27AC T1 Tamoxifen 1-3 macs 71453 8 T7 MCF10 K27AC T7 Tamoxifen 1-3 macs 69588 9 unT7 MCF10 K27AC unT7 unTreated 1-3 macs 65616
In the tab the values for each consensus peaks per condition are 71453,69588,65616.
In the plot, when I make the sum, I get 65824 for T1...What is the reason of this difference ?
Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
# Here I want to keep a subset to make the count
#' ChIPQCexperiment Object is returned so you need to coerce it to a new DBA object
print("Use a subset")
# It seems fragmentSize keep all the samples. not only my subset , it's why i'm getting an error.
So I will try to do the count before creating the subset with dba.mask. But I did that to leverage the count process,to go faster , to count only the library I will process, that was the idea behind.