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
                          Exploring the Dynamics of the Tumor Microenvironment
                          by seqadmin




                          The complexity of cancer is clearly demonstrated in the diverse ecosystem of the tumor microenvironment (TME). The TME is made up of numerous cell types and its development begins with the changes that happen during oncogenesis. “Genomic mutations, copy number changes, epigenetic alterations, and alternative gene expression occur to varying degrees within the affected tumor cells,” explained Andrea O’Hara, Ph.D., Strategic Technical Specialist at Azenta. “As...
                          07-08-2024, 03:19 PM
                        • seqadmin
                          Exploring Human Diversity Through Large-Scale Omics
                          by seqadmin


                          In 2003, researchers from the Human Genome Project (HGP) announced the most comprehensive genome to date1. Although the genome wasn’t fully completed until nearly 20 years later2, numerous large-scale projects, such as the International HapMap Project and 1000 Genomes Project, continued the HGP's work, capturing extensive variation and genomic diversity within humans. Recently, newer initiatives have significantly increased in scale and expanded beyond genomics, offering a more detailed...
                          06-25-2024, 06:43 AM

                        ad_right_rmr

                        Collapse

                        News

                        Collapse

                        Topics Statistics Last Post
                        Started by seqadmin, 07-16-2024, 05:49 AM
                        0 responses
                        28 views
                        0 likes
                        Last Post seqadmin  
                        Started by seqadmin, 07-15-2024, 06:53 AM
                        0 responses
                        34 views
                        0 likes
                        Last Post seqadmin  
                        Started by seqadmin, 07-10-2024, 07:30 AM
                        0 responses
                        40 views
                        0 likes
                        Last Post seqadmin  
                        Started by seqadmin, 07-03-2024, 09:45 AM
                        0 responses
                        205 views
                        0 likes
                        Last Post seqadmin  
                        Working...
                        X