[BioC] Limma; a kind of extended paired analyses with or without treatment

john herbert arraystruggles at gmail.com
Sun Oct 14 12:51:34 CEST 2012


Dear James,
I attempted "design <- model.matrix(~0+treatment*time+sample)"; the
design matrix did not make sense to me, unfortunately. So I stepped
back and did your first suggestion.

patient <- factor(targets$patient)
treat <- factor(targets$treat, levels=c("N","Y"))
time <- factor(targets$time)

treatment <- factor(paste(treat,time,sep="."))
design <- model.matrix(~0+treatment)
colnames(design) <- levels(treatment)

# > head(design)
#   N.0 N.24 N.72 Y.24 Y.72
# 1   0    0    0    1    0
# 2   0    0    0    0    1

corfit <- duplicateCorrelation(y_norm_ave,design,block=targets$patient)

fit <- lmFit(y_norm_ave, design, block=targets$patient,
correlation=corfit$consensus)

cm <- makeContrasts(
Treat24_vs_Normal24 = Y.24-N.24,
Treat72_vs_Normal72 = Y.72-N.72,
Treat24_vs_Normal0 = Y.24-N.0,
Treat72_vs_Normal0 = Y.72-N.0,
Normal24_vs_Normal0 = N.24-N.0,
Normal72_vs_Normal0 = N.72-N.0,
Diff=(Y.72-N.0)-(N.72-N.0),
levels=design)

fit2 <- contrasts.fit(fit, cm)
fit2 <- eBayes(fit2)

topTable(fit2, coef="Treat24_vs_Normal24")

topTable(fit2, coef="Treat72_vs_Normal72").......

A look at some top diff genes, looks reasonably encouraging based on
the biology but take your observation that there is no need for a
random effect model. This style of code and matrix I can make sense
of; is there a way of integrating your second suggestion into this
code, so that the design and contrast matrices make sense to me?

design <- model.matrix(~0+treatment) # (above)

and

design <- model.matrix(~0+treatment*time+sample)

produce very different looking design matrices.

Kind regards,

John.

On Fri, Oct 12, 2012 at 8:23 AM, john herbert <arraystruggles at gmail.com> wrote:
> Thanks James,
> appreciated as you have saved me a lot of time.
>
> John.
>
> On Thu, Oct 11, 2012 at 7:48 PM, James W. MacDonald <jmacdon at uw.edu> wrote:
>> Ugh. Jumped the gun. This does *not* require you to fit a random effects
>> model, as you have done every treatment to cells from each patient. You can
>> just block on Sample and then make your comparisons.
>>
>> In other words, if you add Sample to your design matrix, you will in effect
>> be removing the patient-specific effect. Something like
>>
>> design <- model.matrix(~0+treatment*time+sample)
>>
>> Best,
>>
>> Jim
>>
>>
>>
>>
>> On 10/11/2012 2:27 PM, john herbert wrote:
>>>
>>> Thanks James,
>>> This does not have time course but judging by your answer, I can just
>>> add this in, in place of, say, tissue.
>>>
>>> Kind regards,
>>>
>>> John.
>>>
>>> On Thu, Oct 11, 2012 at 7:23 PM, James W. MacDonald<jmacdon at uw.edu>
>>> wrote:
>>>>
>>>> Hi John,
>>>>
>>>>
>>>> On 10/11/2012 2:15 PM, john herbert wrote:
>>>>>
>>>>> Dear all.
>>>>> I have been pondering about constructing a design matrix based on the
>>>>> Limma user guide, where I combine a time course with a paired
>>>>> analyses. The targets file looks like;
>>>>>
>>>>> Sample  treatment       time
>>>>> 1       control 24
>>>>> 1       control 72
>>>>> 1       control 0
>>>>> 1       treatment       24
>>>>> 1       treatment       72
>>>>> 2       control 24
>>>>> 2       control 72
>>>>> 2       control 0
>>>>> 2       treatment       24
>>>>> 2       treatment       72
>>>>> 3       control 24
>>>>> 3       control 72
>>>>> 3       control 0
>>>>> 3       treatment       24
>>>>> 3       treatment       72
>>>>>
>>>>> Sample number refers to an individuals cancer cells, treatment refers
>>>>> to added drug or not and numbers are in hours (time elapsed). So it is
>>>>> a kind of paired, as patient variability is to be considered. The
>>>>> control sample at 0 is the same as treatment at time 0 as these are
>>>>> the same cells without any time/treatment.
>>>>>
>>>>> Please could someone help me understand how I can construct a design
>>>>> matrix and to understand how I can extract differently expressed genes
>>>>> that come about due to time, due to treatment and interaction of them
>>>>> both.
>>>>>
>>>>> Any pointers appreciated, though I am trying to see if the examples in
>>>>> the manual can be applied to this scenario.
>>>>
>>>>
>>>> See the multi-level experiment example in the user guide, starting on p.
>>>> 47.
>>>>
>>>> Best,
>>>>
>>>> Jim
>>>>
>>>>> Thank you.
>>>>>
>>>>> John.
>>>>>
>>>>> _______________________________________________
>>>>> Bioconductor mailing list
>>>>> Bioconductor at r-project.org
>>>>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>>>>> Search the archives:
>>>>> http://news.gmane.org/gmane.science.biology.informatics.conductor
>>>>
>>>>
>>>> --
>>>> James W. MacDonald, M.S.
>>>> Biostatistician
>>>> University of Washington
>>>> Environmental and Occupational Health Sciences
>>>> 4225 Roosevelt Way NE, # 100
>>>> Seattle WA 98105-6099
>>>>
>>
>> --
>> James W. MacDonald, M.S.
>> Biostatistician
>> University of Washington
>> Environmental and Occupational Health Sciences
>> 4225 Roosevelt Way NE, # 100
>> Seattle WA 98105-6099
>>



More information about the Bioconductor mailing list