Question: Adapting old Workshop to recent changes of the Limma package
gravatar for giroudpaul
2.7 years ago by
giroudpaul30 wrote:


I am a Student currently trying to learn and to understand how to analyse microarray data with R and Bioconductor.

I am following this course :

I am currently in the second part, but I already noticed some problems when doing the following :

tab = topTable(,coef=2,number=2000,adjust.method="BH")
topups = tab[tab[, "logFC"] > 1, ]
IDsup = topups$ID

The last line does not work since ID column has been remplaced by the ID in the rownames, so I worked this out alone.

But now I have trouble with this part : =,contrast.matrix) = eBayes(
results = decideTests(,method='global',adjust.method="BH",p.value=0.05,lfc=1)
res.IDsExvsCup =$genes[results[,1]==1,1]
res.IDsExvsCdown =$genes[results[,1]==-1,1]

Because I get this :

> res.IDsExvsCup
> res.IDsExvsCdown

mostly because "genes" doesn't exist in the MArrayLM :

> names(
 [1] "coefficients"     "rank"             "assign"           "qr"              
 [5] "df.residual"      "sigma"            "cov.coefficients" "stdev.unscaled"  
 [9] "Amean"            "method"           "design"           "contrasts"       
[13] "df.prior"         "s2.prior"         "var.prior"        "proportion"      
[17] ""          "t"                ""         "p.value"         
[21] "lods"             "F"                "F.p.value"       

I don't see what they want to refer to here ? Is it again to be replaced by rownames() ? but how ?

Please help :)


ADD COMMENTlink modified 2.7 years ago by Aaron Lun21k • written 2.7 years ago by giroudpaul30
gravatar for Aaron Lun
2.7 years ago by
Aaron Lun21k
Cambridge, United Kingdom
Aaron Lun21k wrote:

It's probably best if you contact the person who manages the course, as it's their code and they'll be better placed to explain what they were trying to do. Nonetheless, here's some general comments:

  • We don't recommend filtering on the log-fold change directly. If you want to find genes that are above a particular log-fold change, we would suggest using treat.
  • The genes element is designed to hold a data frame containing information for each gene (e.g., gene symbols, descriptions) in each row. If you want downstream functions to make use of this, you should assign an appropriate data frame to$genes.
  • That said, if all you want is the identity of the DE genes, you can just do:
res.IDsExvsCup <- rownames([results[,1]==1]
ADD COMMENTlink modified 2.7 years ago • written 2.7 years ago by Aaron Lun21k

The treat method was presented in the "Compare two groups of sample" tutorial (as alternative to eBayes if I understood well). However, in the "Comparing Multiple groups" part, he just use eBayes. Here, my first example was just to show that the topups$ID did'nt worked. Afterward, he filter first on the adj.P.val and then on the LogFC.

Do you have any link or explanation for the gene data frame and how to assign it ?

All I want is indeed identify the DE genes, and I tried that solution, however, I get a dim error :

> rownames([results[,1]==1,1]
Error in rownames([results[, 1] == 1, 1] : 
  incorrect number of dimensions
> dim(
[1] 45101     3
> dim(results)
[1] 45101     3

Any idea ? Thank for your help

ADD REPLYlink written 2.7 years ago by giroudpaul30

Read my code carefully. There is no second index in the subsetting, as the rownames produces a vector, not a matrix. Also, if you want to assign something to genes, just do:$genes <- data.frame(ID=my.ids, Symbol=my.symbols)

... where you supply a vector of gene IDs and (optionally) symbols in my.ids and my.symbols, respectively. 

ADD REPLYlink modified 2.7 years ago • written 2.7 years ago by Aaron Lun21k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 326 users visited in the last hour