Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • Working with DGEExact objects in a loop

    Objects of class DGEExact are the result of the exactTest() function in edgeR. I'm having trouble making loops to assign DGEExact objects to a variable. I'm probably making some very basic mistake in handling R data types...?

    > nlevels(data$samples$group)
    [1] 8

    Code:
    for(i in 1:nlevels(data$samples$group)) {
       index = 1
       while(index < nlevels(data$samples$group)) {
          extest[i] = exactTest(data, pair=c(i,1+index),)
          index = index + 1
       }
    }
    Error in extest[i] = exactTest(data, pair = c(i, 1 + index), ) :
    object 'extest' not found

    Sorry for the basic question, I hope it is easy to answer!

  • #2
    You need to create "extext" beforehand, that's why it's complaining that it can't be found. Just add this above your while loop:
    Code:
    extest = c()
    Of course, you'll have to deal with extest after the while loop or it'll get over-written.

    Comment


    • #3
      Thank you for your fast reply. Yes, this was completely wrong, and I fixed a few problems. Here is the complete loop:

      Code:
      extest = c()
      k = 1
      for(i in 1:nlevels(data$samples$group)) {
         print(c(i,"i"))
         for(j in (i+1):nlevels(data$samples$group)) {
            if (i < nlevels(data$samples$group)) {
               extest[k] = exactTest(data, pair=c(i,j))
               k = k+1
            }
         }
      }
      This is improved because it tests each pair only once, and I prevented extest from being overwritten. But... extest[k] is not structured as expected... for example, I cannot view head(extest[1]), the whole list prints to screen. And:

      > extest[1]$comparison
      NULL

      How can this strange behavior be explained? Is there any way to store DGEExact objects in a list or other R class?

      Comment


      • #4
        I find that extest[1] is a list:

        > class(extest[1])
        [1] "list"

        whereas it should be:

        > class(extest)
        [1] "DGEExact"
        attr(,"package")
        [1] "edgeR"

        So I need to find out how to store the DGEExact objects somewhere while creating them in the loop. There must be a way to do this????

        Comment


        • #5
          You might have to change extest to a list itself, since I'm not sure how/if coercion works in this scenario.

          Comment


          • #6
            By declaring extest as a list:

            Code:
            extest = list()
            k = 1
            for(i in 1:nlevels(data$samples$group)) {
               #print(c(i,"i"))
               for(j in (i+1):nlevels(data$samples$group)) {
                  if (i < nlevels(data$samples$group)) {
                     extest[k] = exactTest(data, pair=c(i,j))
                     #print(et$table)
                     k = k+1
                  }
               }
            }
            I get this warning for each execution of exactTest:
            28: In extest[k] = exactTest(data, pair = c(i, j)) :
            number of items to replace is not a multiple of replacement length

            > class(extest)
            [1] "list"
            > class(extest[1])
            [1] "list"

            Comment


            • #7
              Not sure then, I'd have to play around with it. You might just play around with making a list of two instance by hand.

              Comment


              • #8
                > first = exactTest(data,pair=c(1,2))
                > second = exactTest(data,pair=c(2,3))
                > results = list(first,second)
                > results
                [[1]]
                An object of class "DGEExact"
                $table
                logFC logCPM PValue
                CGI_10000009 -2.3851234 1.493044 0.1171040419
                CGI_10000013 0.1844446 4.147738 0.7403461115
                CGI_10000014 -4.9376158 2.707417 0.0001137864
                CGI_10000015 0.4967844 3.773692 0.4339343285
                CGI_10000028 -0.2890212 2.476181 0.9883519753
                18055 more rows ...

                $comparison
                [1] "28cell" "blastula"

                $genes
                NULL


                [[2]]
                An object of class "DGEExact"
                $table
                logFC logCPM PValue
                CGI_10000009 -2.57188463 1.493044 0.8000000
                CGI_10000013 0.01444976 4.147738 1.0000000
                CGI_10000014 0.41462168 2.707417 0.9686494
                CGI_10000015 -0.36575812 3.773692 0.6137249
                CGI_10000028 0.64099475 2.476181 0.7954756
                18055 more rows ...

                $comparison
                [1] "blastula" "dlarvae"

                $genes
                NULL

                That seems to work... so why doesn't it work in my loop...

                Comment


                • #9
                  Try making the empty list first and then adding to it, since that will match all the coercion processes going on in in your real code.

                  Comment


                  • #10
                    Yes, that's what I was thinking. Instead of filling indices with list objects, I should append the new result to the list. Hmmmm but how do I do that in R?
                    I'm drawing a blank...

                    Comment


                    • #11
                      For the moment this works, but is really not ideal.
                      Code:
                      for(i in 1:nlevels(data$samples$group)) {
                         for(j in (i+1):nlevels(data$samples$group)) {
                            if (i < nlevels(data$samples$group)) {
                               et = exactTest(data, pair=c(i,j))
                               filen = paste("/media/liz/ACERDATA/Elizabeth_DATA/exacttest_dev_cds", i, j, ".txt", sep="")
                               #Open output file
                               sink(filen)
                               print(et$comparison)
                               print(topTags(et, n=100))
                               #Close output file
                               sink()
                            }
                         }
                      }

                      Comment


                      • #12
                        Actually the solution was super easy! You have to use double [[ to index a list! This works:

                        Code:
                        extest = list()
                        k = 1
                        for(i in 1:nlevels(data$samples$group)) {
                           for(j in (i+1):nlevels(data$samples$group)) {
                              if (i < nlevels(data$samples$group)) {
                                 extest[[k]] = exactTest(data, pair=c(i,j))
                                 k = k+1
                              }
                           }
                        }

                        Comment

                        Latest Articles

                        Collapse

                        • seqadmin
                          Current Approaches to Protein Sequencing
                          by seqadmin


                          Proteins are often described as the workhorses of the cell, and identifying their sequences is key to understanding their role in biological processes and disease. Currently, the most common technique used to determine protein sequences is mass spectrometry. While still a valuable tool, mass spectrometry faces several limitations and requires a highly experienced scientist familiar with the equipment to operate it. Additionally, other proteomic methods, like affinity assays, are constrained...
                          04-04-2024, 04:25 PM
                        • seqadmin
                          Strategies for Sequencing Challenging Samples
                          by seqadmin


                          Despite advancements in sequencing platforms and related sample preparation technologies, certain sample types continue to present significant challenges that can compromise sequencing results. Pedro Echave, Senior Manager of the Global Business Segment at Revvity, explained that the success of a sequencing experiment ultimately depends on the amount and integrity of the nucleic acid template (RNA or DNA) obtained from a sample. “The better the quality of the nucleic acid isolated...
                          03-22-2024, 06:39 AM

                        ad_right_rmr

                        Collapse

                        News

                        Collapse

                        Topics Statistics Last Post
                        Started by seqadmin, 04-11-2024, 12:08 PM
                        0 responses
                        25 views
                        0 likes
                        Last Post seqadmin  
                        Started by seqadmin, 04-10-2024, 10:19 PM
                        0 responses
                        28 views
                        0 likes
                        Last Post seqadmin  
                        Started by seqadmin, 04-10-2024, 09:21 AM
                        0 responses
                        24 views
                        0 likes
                        Last Post seqadmin  
                        Started by seqadmin, 04-04-2024, 09:00 AM
                        0 responses
                        52 views
                        0 likes
                        Last Post seqadmin  
                        Working...
                        X