[BioC] Unexpected results using limma with numerical factor
Matthew Hannah
Hannah at mpimp-golm.mpg.de
Wed Aug 25 09:40:48 CEST 2004
I've recently asked a similar question but have got no feedback, see
point 2 onwards in the following.
https://www.stat.math.ethz.ch/pipermail/bioconductor/2004-August/005802.
html
As I now understand it (perhaps/probably wrongly) the standard approach
(equilivent to ANOVA?) is to use a non-numeric factor for the fit.
However, lm() in R is also capable of regression fits. Looking (but not
understanding) lm.fit in limma, it appears to be an independent function
(lm bit written in fortran?) and doesn't call lm() from stats. So the
question is really if lm.fit can do numeric regression?
Another consideration is that you wouldn't use a design (factor) but a
numeric vector. My guess is that your design is being taken as a factor,
(and if you look in the user guide) the -ve values may indicate dye
swaps, which could be interesting! I've just tried lmfit with a numeric
vector (but as there are replicates each number appears 3 times) and got
meaningless results - all lods>30 and all tiny p-values 1e-24. So
initially it looks like you can't use limma like this, but I'd like to
hear an expert verdict as if the approach is completely wrong.
Anyway if you're looking for correlations then you could perhaps try
pearson (see my previous post about p-values and R2 from lm().
try-
Test1 <- c(0.58,-2.36,-12.24,-14.84,0.15,-3.23,-11.66,-12.91)
Correl <- esApply(eset, 1, function(x) {
cor(x, Test1)
})
be aware that you might get alot of correlations by chance, particularly
as your scores seem to be in 2 groups, close to zero and < -11. And a
straight line between two groups gives a good pearson as it's sensitive
to outliers.
Perhaps it's best to change your test scores to factors anyway-
Test1.high, Test1.low, Test2.high, Test2.low, and do a conventional
limma analysis. Without a regular distribution of test scores,
correlations are going to be largely meaningless.
HTH, and that someone can answer this more definetely.
Cheers,
Matt
Hi again,
For an experiment I'm analyzing, I do not have a series of factors.
Rather, I
have a pair of test-score (numerical) for each replicate animal. This
is
Affymetrix data, and for my initial pass I tried with only a single
test-score.
The commands I used are below:
> library(gcrma);
Welcome to Bioconductor
Vignettes contain introductory material. To view,
simply type: openVignette()
For details on reading vignettes, see
the openVignette help page.
> library(limma);
> cel.files <- c(
+ 'RAE230_2_060104_LH_IM07T.CEL',
+ 'RAE230_2_060104_LH_IM08T.CEL',
+ 'RAE230_2_060104_LH_IM09T.CEL',
+ 'RAE230_2_060104_LH_IM10T.CEL',
+ 'RAE230_2_060204_LH_IM07T.CEL',
+ 'RAE230_2_060204_LH_IM08T.CEL',
+ 'RAE230_2_060204_LH_IM09T.CEL',
+ 'RAE230_2_060204_LH_IM10T.CEL'
+ );
> eset <- ReadAffy(filenames=cel.files, phenoData="phenodata.txt");
> eset <- ReadAffy(filenames=cel.files, phenoData="phenodata.txt");
> eset;
AffyBatch object
size of arrays=834x834 features (43476 kb)
cdf=Rat230_2 (31099 affyids)
number of samples=8
number of genes=31099
annotation=rat2302
> pData(eset);
TestScore1
RAE230_2_060104_LH_IM07T.CEL 0.58
RAE230_2_060104_LH_IM08T.CEL -2.36
RAE230_2_060104_LH_IM09T.CEL -12.24
RAE230_2_060104_LH_IM10T.CEL -14.84
RAE230_2_060204_LH_IM07T.CEL 0.15
RAE230_2_060204_LH_IM08T.CEL -3.23
RAE230_2_060204_LH_IM09T.CEL -11.66
RAE230_2_060204_LH_IM10T.CEL -12.91
> eset <- rma(eset);
Background correcting
Normalizing
Calculating Expression
> design <- model.matrix(~-1 + TestScore1, pData(eset));
> design;
TestScore1
RAE230_2_060104_LH_IM07T.CEL 0.58
RAE230_2_060104_LH_IM08T.CEL -2.36
RAE230_2_060104_LH_IM09T.CEL -12.24
RAE230_2_060104_LH_IM10T.CEL -14.84
RAE230_2_060204_LH_IM07T.CEL 0.15
RAE230_2_060204_LH_IM08T.CEL -3.23
RAE230_2_060204_LH_IM09T.CEL -11.66
RAE230_2_060204_LH_IM10T.CEL -12.91
attr(,"assign")
[1] 1
> fit1 <- lmFit(eset, design);
> fit3 <- eBayes(fit1);
All proceeds well without any error-messages, so I believed I had
successfully
fit my model. When I extract the data, however, I get some unexpected
results:
> topTable(fit3, coef=1, number=20, adjust="fdr");
ID M A t P.Value
B
104 1367555_at -1.212175 14.77173 -6.808742 3.912627e-07
14.75482
105 1367556_s_at -1.203275 14.67175 -6.762709 3.912627e-07
14.48816
3411 1370862_at -1.185363 14.42457 -6.674599 3.912627e-07
13.98156
2777 1370228_at -1.184768 14.46704 -6.666414 3.912627e-07
13.93475
549 1368000_at -1.175866 14.33987 -6.622929 3.912627e-07
13.68683
2697 1370148_at -1.174858 14.32747 -6.617795 3.912627e-07
13.65764
837 1368288_at -1.170412 14.25957 -6.596453 3.912627e-07
13.53649
710 1368161_a_at -1.169621 14.27702 -6.589664 3.912627e-07
13.49801
420 1367871_at -1.166552 14.09481 -6.588461 3.912627e-07
13.49120
2558 1370009_at -1.162241 14.18252 -6.552347 3.912627e-07
13.28707
147 1367598_at -1.161072 14.19550 -6.543620 3.912627e-07
13.23787
2576 1370027_a_at -1.158955 14.13682 -6.536068 3.912627e-07
13.19533
1136 1368587_at -1.154447 14.04036 -6.517046 3.912627e-07
13.08836
196 1367647_at -1.154999 14.07951 -6.516673 3.912627e-07
13.08627
31081 AFFX-r2-P1-cre-5_at -1.153568 14.05140 -6.510380 3.912627e-07
13.05094
19150 1387082_at -1.153501 14.05431 -6.509674 3.912627e-07
13.04697
2898 1370349_a_at -1.153164 14.07047 -6.505935 3.912627e-07
13.02599
1235 1368686_at -1.152660 14.04966 -6.504796 3.912627e-07
13.01960
8200 1375651_at -1.152557 14.04428 -6.504674 3.912627e-07
13.01892
19359 1387291_at -1.151666 14.03699 -6.499744 3.912627e-07
12.99127
The near-identical M, A, and p-values indicate a problem, and none of
the genes
on this list are very plausible biologically for our test-system. Based
on that
I'm pretty sure I've gone astray somewhere.
Is it possible to use numeric scores in fitting a linear model with
limma? If
so, am I asking the question in the right manner? If not, are there any
BioConductor tools appropriate for this kind of question?
Any help very much appreciated,
Paul
More information about the Bioconductor
mailing list