Issue with BiSeq for DMR analysis

[shocker8786]

I'm hoping someone has experience using the BiSeq package in R. I'm going through the "An Introduction to BiSeq" pdf available at the following website:

http://www.bioconductor.org/packages/release/bioc/html/BiSeq.html

So far I have successfully imported my data into R using the readBismark command suggested in the tutorial. I have been able to work through sections 2 and 3 without error. However, when I get to section 4.1, definition of CpG clusters, I run into an error. I am using the same code they do, but changing $group to $treatment, as that is what my column in colData is named. Below is the code, and resulting error.

rrbs.clust<-clusterSites(object=rrbs,groups=colData(rrbs)$treatment,perc.samples=4/5,min.sites=20,mx.dist=100)
Error in .filterBySharedRegions_BSraw(object, groups, perc.samples, minCov = 1) :
Length of no.samples should be 1 or the same as the number of group levels.

There seems to be an issue with the no.samples, but I'm not familiar enough with R to figure out how to fix the problem. Any help would be appreciated. Thank you very much!

[dpryan]

Are you literally following the vignette or is this your own data? Also, what's the output of:

sessionInfo()

It would also be informative if you posted the output of the following:

colData(rrbs)$treatment
levels(colData(rrbs)$treatment)

This sort of error might happen if you originally had more groups in your design and then removed a few (so they wouldn't then appear in colData(rrbs)$treatment but there would, nevertheless, be a factor level for them).

[shocker8786]

Yeah sorry, I guess that wasn't clear. I've gone through the vignette with the test dataset, and am now using my own dataset. Here are the results from the requested commands.

sessionInfo()
R version 3.0.3 (2014-03-06)
Platform: x86_64-unknown-linux-gnu (64-bit)

locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] parallel stats graphics grDevices utils datasets methods
[8] base

other attached packages:
[1] BiSeq_1.2.5 Formula_1.1-1 GenomicRanges_1.14.4
[4] XVector_0.2.0 IRanges_1.20.7 BiocGenerics_0.8.0

loaded via a namespace (and not attached):
[1] annotate_1.40.1 AnnotationDbi_1.24.0 betareg_3.0-4
[4] Biobase_2.22.0 Biostrings_2.30.1 bitops_1.0-6
[7] BSgenome_1.30.0 DBI_0.2-7 flexmix_2.3-11
[10] globaltest_5.16.0 grid_3.0.3 lattice_0.20-29
[13] lmtest_0.9-33 lokern_1.1-5 modeltools_0.2-21
[16] nnet_7.3-8 RCurl_1.95-4.1 Rsamtools_1.14.3
[19] RSQLite_0.11.4 rtracklayer_1.22.7 sandwich_2.3-0
[22] sfsmisc_1.0-25 splines_3.0.3 stats4_3.0.3
[25] survival_2.37-7 tools_3.0.3 XML_3.98-1.1
[28] xtable_1.7-3 zlibbioc_1.8.0 zoo_1.7-11

colData(rrbs)$treatment
[1] "control" "treeatment" "control" "control" "treatment"
[6] "control" "treatment"

> levels(colData(rrbs)$treatment)
NULL

This is the full dataset, I haven't removed any groups or anything. I have a feeling I'm not supposed to receive NULL for the last command.

[dpryan]

Your gut isn't leading you astray on this one. I'm guessing that colData(rrbs)$treatment is a character vector rather than a factor, so try just:

colData(rrbs)$treatment <- factor(colData(rrbs)$treatment)

and see if that fixes the problem (it likely will).

[shocker8786]

Thanks, I think we're getting closer. I received a different error message this time.

> rrbs.clust<-clusterSites(object=rrbs,groups=colData(rrbs)$treatment,perc.samples=4/5,min.sites=20,max.dist=100)
Error in rowSums(part >= minCov, na.rm = TRUE) :
'x' must be an array of at least two dimensions

[dpryan]

I would guess that this is due to you misspelling "treatment" as "treeatment" for one of your samples above (the second sample according to your post above).

[shocker8786]

Yes, that seems to have solved the problem. Thank you very much!

[jeni]

bismark output comparision of treatment and control sample

Hello,

I have imported the bismark file but I have 2 files one is treatment and the other is for control in the bismark2bedGraph output format. But I did not found any command to create the BSraw object which contain all the information from treatment and control file. So, How to create the comparable BSraw object from bismark output file?

Thanks in advance.

[mgavery]

Hi,

I know this is an old post, but I thought I'd try and ask a question here.. I am trying to run the clusterSites command with my data but it is crashing because it's using excessive memory. Did you run into that problem?

[dpryan]

How much are you considering excessive? I recall that the method needed for this is rather RAM demanding.

[mgavery]

The cluster I was running it on was hitting 80 GB of memory then crashing. This was with a subset of 4 samples. I have 20 total samples in this dataset. How much memory did you have to work with when you did this analysis?

[dpryan]

I would think that 80 gigs would suffice, though the ones I was using had 250. My only suggestion would be to do one chromosome at a time.

[mgavery]

Thanks for the suggestion!

[jianl]

It seems to be a relatively old thread, but I'm having the exactly same problem as the original post, just my problem cannot be resolved as it.

I had the same error
"Length of no.samples should be 1 or the same as the number of group levels"
as OP, and it can be fixed by transforming the "colData(rrbs)$treatment" into a factor. However, the second error
"Error in rowSums(part >= minCov, na.rm = TRUE) :
'x' must be an array of at least two dimensions" cannot be fixed by the suggested solution since I have only two samples assigned with "test" and "control", so there's no mis-spelling. I'm wondering what can I possibly do to correct that?

Here is my script:

library(BiSeq)
library("GenomicRanges")
rrbs=readBismark(c(“controlMerged.bismark.cov.gz”,”testMerged.bismark.cov.gz"),colData= DataFrame(group=factor(c(“control","test”)),row.names=c(“C1”,”T1")))

rrbs.clust.unlim <- clusterSites(object = rrbs, groups = colData(rrbs)$group, perc.samples = 4/5 ,min.sites = 2,max.dist = 100)

>sessionInfo()
R version 3.2.5 (2016-04-14)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 7 (wheezy)

locale:
[1] C

attached base packages:
[1] parallel stats4 stats graphics grDevices utils datasets
[8] methods base

other attached packages:
[1] BiSeq_1.10.0 Formula_1.2-1
[3] SummarizedExperiment_1.0.2 Biobase_2.30.0
[5] GenomicRanges_1.22.4 GenomeInfoDb_1.6.3
[7] IRanges_2.4.8 S4Vectors_0.8.11
[9] BiocGenerics_0.16.1

loaded via a namespace (and not attached):
[1] futile.logger_1.4.3 XVector_0.10.0 bitops_1.0-6
[4] futile.options_1.0.0 tools_3.2.5 zlibbioc_1.16.0
[7] annotate_1.48.0 RSQLite_1.0.0 lattice_0.20-34
[10] Matrix_1.2-7.1 DBI_0.5-1 rtracklayer_1.30.4
[13] Biostrings_2.38.4 lmtest_0.9-34 grid_3.2.5
[16] nnet_7.3-12 globaltest_5.24.0 flexmix_2.3-13
[19] AnnotationDbi_1.32.3 survival_2.39-5 XML_3.98-1.4
[22] BiocParallel_1.4.3 lokern_1.1-6 lambda.r_1.1.9
[25] splines_3.2.5 Rsamtools_1.22.0 modeltools_0.2-21
[28] sfsmisc_1.1-0 GenomicAlignments_1.6.3 xtable_1.8-2
[31] betareg_3.1-0 sandwich_2.3-4 RCurl_1.95-4.8
[34] zoo_1.7-13

> colData(rrbs)$group
[1] control test
Levels: control test

> levels(colData(rrbs)$group)
[1] "control" "test"


I also tried to increase the sample numbers, i.e. two tests and two controls, and to test things quickly I only restricted to 10000 CpG sites. However, the command "clusterSites" then failed with error:
"Error in mcfork(detached) :
unable to fork, possible reason: Cannot allocate memory" (I guess that's why it now doesn't complain "'x' must be an array of at least two dimensions")

I don't really think 4 samples with 10000CpGs can require to much memory ... But any advice will be highly appreciated. Thanks!