[BioC] How to create a design matrix for a three-way design?

James W. MacDonald jmacdon at med.umich.edu
Fri Mar 26 14:42:26 CET 2010


Hi January,

January Weiner wrote:
> Hello,
> 
> I would be grateful for any kind of hint -- redirection to other
> posts, documents, examples etc. I have trouble getting a clear picture
> how I should create the design matrix and extract the contrasts.
> 
> I have the following situation: I am comparing two strains of cells,
> wild type and mutant. The cells are subjected to treatment (T) or kept
> as control (C), and RNA is collected at two different times after
> experiment setup (t1, t2). What I would like to know is whether the
> response to treatment is different in the wild type and the mutant at
> any point in time. Two-color arrays are used with dye-swaps, and in
> each the treatment is compared to the appropriate control.

Since you only have two time points, what you are looking for is called 
the 'interaction term', which measures the difference in response 
between the two sample types at the two different times.

> 
> My targets look like this:
> 
> Cy3       Cy5
> wt_t1_T  wt_t1_C
> wt_t1_C  wt_t1_T
> wt_t2_T  wt_t2_C
> wt_t2_C  wt_t2_T
> mu_t1_T mu_t1_C
> mu_t1_C mu_t1_T
> mu_t2_T mu_t2_C
> mu_t2_C mu_t2_T
> 
> What I did (and what was surely not correct) was to call any "C" as
> "control", and use this as a reference:

Actually, it was correct.

> 
> Cy3   Cy5
> wt_t1  control
> control wt_t1
> wt_t2  control
> control  wt_t2
> mu_t1  control
> control mu_t1
> mu_t2  control
> control mu_t2
> 
> design <- modelMatrix( target, ref="control" )
> fit <- lmFit( ma, design )
> contrast <- makeContrast( foo=(wt_t1-mu_t1)+(wt_t2-mu_t2), levels=design)
> 
> contrast:
>> contrast
>           Contrasts
> Levels      foo
>   mu_t1    -1
>   mu_t2    -1
>   wt_t1      1
>   wt_t2      1
> 
> fit2 <- contrasts.fit( fit, contrast )
> fit2 <- eBayes( fit )
> 
> etc. This seems to work... somewhat. If I run separate contrasts to
> compare wt and mu at t1, and wt and mu at t2, I get different results,
> and also way more significant. I am confused.

Yes you are. ;-D

But you have actually done OK. The code you show here will test for any 
genes that react differently in the wt vs mu samples at the different 
time points. Some examples:

A gene is unchanged in wt at the two time points but goes up (or down) 
in mu.

The converse of the above (i.e., unchanged in mu, but changes in wt).

A gene goes up in wt from time 1 to time 2, but goes down in mu from 
time 1 to time 2.

The converse of the above (i.e., down in wt and up in mu).

Note here that all of the situations above are dependent on what happens 
with both sample types, and what happens has to be different. This is a 
special situation, so you should expect fewer genes to fulfill these 
criteria.

A simpler situation is when we only look at the differences between two 
samples at one time (e.g., wt vs mu at time 1, or time1 vs time2 for mu 
samples). As an example, consider the following scenario:

	T1	T2
MU	-3.4	-3.6
WT	4.9	4.8

You can see here that the gene isn't affected differently between mu and 
wt at the different time points, so won't result in a significant 
interaction term. However, mu and wt are very different from each other 
at each time point, so a simple comparison of wt and mu at e.g., t1 will 
be highly significant.

It's way more likely that you will find genes that are differentially 
expressed between mu and wt at a single time point as compared to the 
more complicated requirements of the interaction, so it is not 
surprising that you find more genes when you compare wt vs mu at single 
time points.

Does that help?

Best,

Jim


> 
> best regards,
> j.
> 

-- 
James W. MacDonald, M.S.
Biostatistician
Douglas Lab
University of Michigan
Department of Human Genetics
5912 Buhl
1241 E. Catherine St.
Ann Arbor MI 48109-5618
734-615-7826
**********************************************************
Electronic Mail is not secure, may not be read every day, and should not be used for urgent or sensitive issues 



More information about the Bioconductor mailing list