[BioC] LIMMA lmscFit: difference between linear models

[Koen Bossers]

Francois (and other bioconductor users),

Thanks for your explanation. My thinking was in that direction also, and
you made it very clear to me.

I do have some questions though. The second model is more powerful,
because we can now model the technical variation first, and then the
patient variation. This then gives me far more significant genes.

But the thing is, I think that the technical variation is much smaller
than the biological variation to begin with (high quality Agilent arrays
vs human samples)! For example, the distance between all pat1 channels in
a hierarchical clustering is much smaller than between individuals.
Doesn't that imply that "removing" the technical variance shouldn't have
such a big impact on my analysis results? So is the linear model doing
something else that I'm missing here?

Second, and this is something that is more related to some limitations in
the LIMMA package, the technical replicates are now treated as biological
replicates. So probably the number of significant genes is overestimated.
But to what extend? Is there a way the estimate this overestimation? Or
shoul I let the qPCR validation do the talking?

Does anybody have any comments/suggestions on the best method to detect
differential gene expression in an experiment with low technical variance
and high biological variance?

Thank you,

Koen

-- 
Koen Bossers, PhD student
Netherlands Institute for Brain Research
Meibergdreef 33
1105 AZ Amsterdam, The Netherlands
Phone: +31-20-5665512
Email: k.bossers at nih.knaw.nl

> Hi Koen,
>
> Both models are "valid", but that doesn't mean that they're equally good.
>
> Here is the main difference between the two models (please correct me if
I'm wrong):
>
> In the first one, you are trying to fit a single coefficients for each
gene in patient and control. This interprets all variations between
patients and controls as error and gives you a rather high variance. So,
for each gene you have a 2 value (for the patients and for the controls)
which you won't trust much because of the variance.
>
> In the second case, you're fitting each patient and control separately
first. Assuming you have decent arrays, then the variance between the
replicates is pretty low. So now you have 8 values (4 for the patients,
4 for the controls) which you'll trust more, because they have lower
variances. Then those are then used for the contrast.
>
> So the second method is more powerful, because it can model the
> technical variation separately from the patient variation. This is the
intuition why the second method can call more genes as being
> differentially expressed.
>
> Francois
>
> On Mon, 2005-26-09 at 10:46 +0200, Koen Bossers wrote:
>> Dear all,
>>
>> I'm currently studying gene expression in human brain samples using
Agilent arrays (4 patients, 4 controls). I am analyzing the data using
a single channel approach (lmscFit), which I think valid for my
dataset, as all channels from each individual nicely cluster together,
apart from all other channels.
>>
>> The hybridization setup is as follows:
>>
>> --------------
>>                                    Cy3       Cy5
>> US12302316_251182110152_S01_A01   pat2      ctrl3
>> US12302316_251182110153_S01_A01   pat3      ctrl2
>> US12302316_251182110154_S01_A01   pat1      ctrl4
>> US12302316_251182110155_S01_A01   ctrl4     pat4
>> US12302316_251182110156_S01_A01   ctrl1     pat1
>> US12302316_251182110157_S01_A01   pat4      ctrl1
>> US12302316_251182110158_S01_A01   ctrl2     pat4
>> US12302316_251182110159_S01_A01   pat2      ctrl1
>> US12302316_251182110160_S01_A01   ctrl3     pat1
>> US12302316_251182110176_S01_A01   pat3      ctrl4
>> US12302316_251182110177_S01_A01   ctrl2     pat2
>> US12302316_251182110178_S01_A01   ctrl3     pat3
>>
>> --------------
>>
>> The first analysis I tried was the following: I replaced all individual
labels with a generic one (thus: pat1 becomes pat, ctrl2 becomes ctrl),
and calculated a contrast between pat&ctrl:
>>
>> --------------
>>
>> targets2 <- targetsA2C(targets)
>> u <- unique(targets2$Target)
>> f <- factor(targets2$Target, levels=u)
>> design <- model.matrix(~0+f)
>> colnames(design) <- u
>>
>> corfit <- intraspotCorrelation(MA, design)
>> fit <- lmscFit(MA, design, correlation=corfit$consensus)
>>
>> cont.matrix <- makeContrasts(pat-ctrl,levels=design)
>> fit2 <- contrasts.fit(fit, cont.matrix)
>> fit2 <- eBayes(fit2)
>>
>> --------------
>>
>> This analysis does not yield many significantly regulated genes (5 or
so), which is likely due to the small number of biological replicates
and the large diversity in the human population.
>>
>> I also tried another approach, leaving the individual labels intact,
and fitting a linear model in the following manner:
>>
>> --------------
>>
>> cont.matrix <-
>> makeContrasts((pat1+pat2+pat3+pat4-ctrl1-ctrl2-ctrl3-ctrl4)/4
>>                 ,levels=design)
>>
>>
>> --------------
>>
>> Is this linear model valid?
>> This analysis yields loads of significantly regulated genes (hundreds)!
Neither the MA plot or the M values in fit2 look suspicious, so I do
not have a reason to distrust this data.
>>
>> I do not really understand why there is such a large discrepancy
between the two analysis methods. Is this due to the way replication is
handled?
>>
>> Could anybody comment on the validity of these two analyses, taking
into account individual variation in the human population, and the way
replication is handled in LIMMA?
>>
>> Thank you very much,
>>
>> Koen Bossers
>>
>>
>