Working with DGEExact objects in a loop

[elizabeth000]

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

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!

[dpryan]

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:

extest = c()

Of course, you'll have to deal with extest after the while loop or it'll get over-written.

[elizabeth000]

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

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?

[elizabeth000]

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????

[dpryan]

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

[elizabeth000]

By declaring extest as a list:

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"

[dpryan]

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.

[elizabeth000]

> 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...

[dpryan]

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.

[elizabeth000]

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...

[elizabeth000]

For the moment this works, but is really not ideal.

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()
      }
   }
}

[elizabeth000]

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

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
      }
   }
}