[BioC] Detection of differential expression using limma
christianeisen at alice-dsl.de
Mon Oct 13 19:18:20 CEST 2008
Hi Sean and everyone,
thank you for your reply. I was aware that for log intensities the fold
change is calculated via substraction.
I have posted on example of intensity values for a ample gene in the two
groups I am interessted in.
I have listed all the respective values for that particular gene after
transformation with the VSN or quantile method
as well as simple log2 transformation of the raw data.
10.4093909361377 12.4491486453754 (log2 transformed data)
10.0910975543547 10.5456558587892 (VSN transformed data)
10.1711963483519 10.4027789253209 (quantile transformed data)
1360 5592 (raw data)
Like mentioned previously I am using the limma package.
I read the tab-deliminted files containing the raw data as advised by
Gordon Smyth and Peter White in a
previous post handling the same topic.
mrawObj <- read.maimages(targets$FileName,
columns = list(G = "gMedianSignal", Gb = "gBGUsed", R =
"gProcessedSignal",Rb = "gBGMedianSignal"),
annotation= c("ProbeName", "GeneName", "SystematicName","ControlType"))
mnormObj <- mrawObj
Then I perform normalization only on the green channel
mnormObj$G <- log2(mnormObj$G)
After that I extract the transformed intensity data for the green
channel into a new
object and name the rows after the genes. As I have data for 16 arrays,
created object contains 16 columns, one for each array, and around
one for each individual gene.
For identification of differentially expressed genes I proceed as described
in the limma userguide by defining a model.matrix and a cont.matrix to
comparisons that I want the program to compute.
After this I fit a linear model to my data.
fit <- lmFit(preproc_data, design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)
and display the results for the different comparisons listed in the
Below a sample output from the VSN transformed data
> topTable(fit2, coef=1, adjust="BH")
ID logFC t
P.Value adj.P.Val B
ebv-miR-BART18-5p -2.788856 -3.925838 0.001719277
hsa-miR-671-5p 2.378257 3.875517 0.001891590
hsa-miR-663 2.308951 3.872746 0.001901578
hsa-miR-583 -1.509068 -3.758990
0.002361635 0.9999662 -4.506259
hsa-miR-671-5p 2.097918 3.661009 0.002848245
hsa-miR-671-5p 2.055704 3.576170 0.003351448
hsa-miR-663 2.354609 3.513237
0.003782239 0.9999662 -4.513012
hsa-miR-671-5p 2.136401 3.480265 0.004029907
hsa-miR-145* -3.167342 -3.452580
0.004250482 0.9999662 -4.514748
hsa-miR-373 1.602810 3.381800
0.004871484 0.9999662 -4.516809
so the strange thing is, when looking at the entry for the top hit
indetified by limma anylsis you get something like this
5.20945336562895 5.4093909361377 (log2 transformation)
-0.971205336133305 1.69436870845910 (VSN transformed data)
5.19967234483636 5.31174831500496 (quantile transformation)
37 42.5 (raw data)
Now you may say, sure this is clearly a relict of VSN normalization. No
doubt about that
but for the other methods the problem is still similar just the names of
the genes are different
and the fold change and p-Values.
So as you see I don't understand why, for example the example above
let-7a doesn't show up
in the top differentially expressed gene list, while others, which a
clearly not differentially expressed
Thanks for any kind of advice helping me out on this!
More information about the Bioconductor