Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • Issue with BiSeq for DMR analysis

    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:

    The BiSeq package provides useful classes and functions to handle and analyze targeted bisulfite sequencing (BS) data such as reduced-representation bisulfite sequencing (RRBS) data. In particular, it implements an algorithm to detect differentially methylated regions (DMRs). The package takes already aligned BS data from one or multiple samples.


    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!

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

    Code:
    sessionInfo()
    It would also be informative if you posted the output of the following:

    Code:
    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).

    Comment


    • #3
      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.

      Comment


      • #4
        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:

        Code:
        colData(rrbs)$treatment <- factor(colData(rrbs)$treatment)
        and see if that fixes the problem (it likely will).

        Comment


        • #5
          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

          Comment


          • #6
            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).

            Comment


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

              Comment


              • #8
                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.

                Comment


                • #9
                  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?

                  Comment


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

                    Comment


                    • #11
                      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?

                      Comment


                      • #12
                        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.

                        Comment


                        • #13
                          Thanks for the suggestion!

                          Comment


                          • #14
                            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!

                            Comment

                            Latest Articles

                            Collapse

                            • seqadmin
                              Recent Advances in Sequencing Analysis Tools
                              by seqadmin


                              The sequencing world is rapidly changing due to declining costs, enhanced accuracies, and the advent of newer, cutting-edge instruments. Equally important to these developments are improvements in sequencing analysis, a process that converts vast amounts of raw data into a comprehensible and meaningful form. This complex task requires expertise and the right analysis tools. In this article, we highlight the progress and innovation in sequencing analysis by reviewing several of the...
                              05-06-2024, 07:48 AM
                            • seqadmin
                              Essential Discoveries and Tools in Epitranscriptomics
                              by seqadmin




                              The field of epigenetics has traditionally concentrated more on DNA and how changes like methylation and phosphorylation of histones impact gene expression and regulation. However, our increased understanding of RNA modifications and their importance in cellular processes has led to a rise in epitranscriptomics research. “Epitranscriptomics brings together the concepts of epigenetics and gene expression,” explained Adrien Leger, PhD, Principal Research Scientist...
                              04-22-2024, 07:01 AM

                            ad_right_rmr

                            Collapse

                            News

                            Collapse

                            Topics Statistics Last Post
                            Started by seqadmin, 05-14-2024, 07:03 AM
                            0 responses
                            26 views
                            0 likes
                            Last Post seqadmin  
                            Started by seqadmin, 05-10-2024, 06:35 AM
                            0 responses
                            45 views
                            0 likes
                            Last Post seqadmin  
                            Started by seqadmin, 05-09-2024, 02:46 PM
                            0 responses
                            59 views
                            0 likes
                            Last Post seqadmin  
                            Started by seqadmin, 05-07-2024, 06:57 AM
                            0 responses
                            46 views
                            0 likes
                            Last Post seqadmin  
                            Working...
                            X