Hi Ryan and all the list,

here I am again. I have another question regarding my factorial design. As
Ryan suggested me I estimated the Model 2, i.e. I defined each treatment
combination as a group. This allow me to make a big number of comparisons
by defining a contrast for every comparison I can be interested in. For
example, I can find differentially expressed genes between GENOTYPE F with
LOAD High at TIME H and GENOTYPE F with LOAD High at TIME PH; or I can
compare GENOTYPE F with LOAD Medium at TIME H and GENOTYPE F with LOAD
Medium at TIME PH, and so on.
I found a total of 45 contrasts, so I built the Contrasts Matrix:

contrasts<-makeContrasts(
F.High.HvsPH = F.High.H - F.High.PH,
F.Medium.HvsPH = F.Medium.H - F.Medium.PH,
F.Low.HvsPH = F.Low.H - F.Low.PH,

F.HighvsMedium.H = F.High.H - F.Medium.H,
F.HighvsMedium.PH = F.High.PH - F.Medium.PH,
  F.HighvsLow.H = F.High.H - F.Low.H,
F.HighvsLow.PH = F.High.PH - F.Low.PH,
F.MediumvsLow.H = F.Medium.H - F.Low.H,
F.MediumvsLow.PH = F.Medium.PH - F.Low.PH,

  G.High.HvsPH = G.High.H - G.High.PH,
G.Medium.HvsPH = G.Medium.H - G.Medium.PH,
G.Low.HvsPH = G.Low.H - G.Low.PH,

G.HighvsMedium.H = G.High.H - G.Medium.H,
G.HighvsMedium.PH = G.High.PH - G.Medium.PH,
G.HighvsLow.H = G.High.H - G.Low.H,
G.HighvsLow.PH = G.High.PH - G.Low.PH,
G.MediumvsLow.H = G.Medium.H - G.Low.H,
G.MediumvsLow.PH = G.Medium.PH - G.Low.PH,

P.High.HvsPH = P.High.H - P.High.PH,
P.Medium.HvsPH = P.Medium.H - P.Medium.PH,
P.Low.HvsPH = P.Low.H - P.Low.PH,

P.HighvsMedium.H = P.High.H - P.Medium.H,
P.HighvsMedium.PH = P.High.PH - P.Medium.PH,
P.HighvsLow.H = P.High.H - P.Low.H,
P.HighvsLow.PH = P.High.PH - P.Low.PH,
P.MediumvsLow.H = P.Medium.H - P.Low.H,
P.MediumvsLow.PH = P.Medium.PH - P.Low.PH,


FvsG.High.H = F.High.H - G.High.H,
FvsG.High.PH = F.High.PH - G.High.PH,
FvsG.Medium.H = F.Medium.H - G.Medium.H,
FvsG.Medium.PH = F.Medium.PH - G.Medium.PH,
FvsG.Low.H = F.Low.H - G.Low.H,
FvsG.Low.PH = F.Low.PH - G.Low.PH,

FvsP.High.H = F.High.H - P.High.H,
FvsP.High.PH = F.High.PH - P.High.PH,
FvsP.Medium.H = F.Medium.H - P.Medium.H,
FvsP.Medium.PH = F.Medium.PH - P.Medium.PH,
FvsP.Low.H = F.Low.H - P.Low.H,
FvsP.Low.PH = F.Low.PH - P.Low.PH,

GvsP.High.H = G.High.H - P.High.H,
GvsP.High.PH = G.High.PH - P.High.PH,
GvsP.Medium.H = G.Medium.H - P.Medium.H,
GvsP.Medium.PH = G.Medium.PH - P.Medium.PH,
GvsP.Low.H = G.Low.H - P.Low.H,
GvsP.Low.PH = G.Low.PH - P.Low.PH,
levels=mydesign)


My problem is that the contrasts I made allow me to compare only one factor
considering the others fixed, i.e. I can compare two levels of one factor
fixing the other two factors at the same level. So I thik in this way I
can't test for interactions. I would like to be able to find also for which
genes the interaction between two factors is significative. For example, I
would like to find which genes change significantly their expression, for
GENOTYPE G, in a different way from High and Medium CHARGE between TIMEs H
and PH, i.e for which genes of GENOTYPE F the interaction between CHARGE
and TIME is significative. Can I do so by defining a contrast like this?

  (F.High.H - F.High.PH) - (F.Medium.H - F.Medium.PH)


And what if I am interested in obtaining a list of DEG for a comparison of
two levels of a factor indipendently from the levels of the other two
factors? Can I define a contrast like this?

  FvsG = F.High.H + F.High.PH + F.Medium.H + F.Medium.PH + F.Low.H +
F.Low.PH
       - G.High.H - G.High.PH - G.Medium.H - G.Medium.PH - G.Low.H -
G.Low.PH

and then do a LRT to find the differentially expressed genes?
Or does it make more sense to obtain a list of DEG for every contrast I
made before and then obtaining the FvsG list of DEG doing the union of all
the lists of differentially expressed genes like this:

  # Let DE.FvsG.High.H, DE.FvsG.High.PH, DE.FvsG.Medium.H, DE.FvsG.Medium.PH
,
  # DE.FvsG.Low.H, DE.FvsG.Low.PH be the lists of differentially expressed
genes
  # for the related contrasts of the Contrast Matrix.

  DE.FvsG<-unique(c(DE.FvsG.High.H, DE.FvsG.High.PH, DE.FvsG.Medium.H,
      DE.FvsG.Medium.PH, DE.FvsG.Low.H, DE.FvsG.Low.PH))

I'm so sorry for my gaps in models with interactions and more than two
factors of interst but it's my first time with a so complex design so I'm
trying to understand.
Thanks in advance for any replies!

Marco



> sessionInfo()
R version 3.0.2 (2013-09-25)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=Italian_Italy.1252  LC_CTYPE=Italian_Italy.1252
[3] LC_MONETARY=Italian_Italy.1252 LC_NUMERIC=C
[5] LC_TIME=Italian_Italy.1252

attached base packages:
[1] splines   parallel  stats     graphics  grDevices utils     datasets
[8] methods   base

other attached packages:
 [1] EDASeq_1.8.0         ShortRead_1.20.0     Rsamtools_1.14.2
 [4] lattice_0.20-24      Biostrings_2.30.1    GenomicRanges_1.14.4
 [7] XVector_0.2.0        IRanges_1.20.6       Biobase_2.22.0
[10] BiocGenerics_0.8.0   edgeR_3.4.2          limma_3.18.10
[13] aroma.light_1.32.0   matrixStats_0.8.14   yeastRNASeq_0.0.10

loaded via a namespace (and not attached):
 [1] annotate_1.40.0      AnnotationDbi_1.24.0 bitops_1.0-6
 [4] DBI_0.2-7            DESeq_1.14.0         genefilter_1.44.0
 [7] geneplotter_1.40.0   grid_3.0.2           hwriter_1.3
[10] latticeExtra_0.6-26  R.methodsS3_1.6.1    R.oo_1.17.0
[13] RColorBrewer_1.0-5   RSQLite_0.11.4       stats4_3.0.2
[16] survival_2.37-7      XML_3.98-1.1         xtable_1.7-1
[19] zlibbioc_1.8.0





2014-01-30 Marco Pegoraro <marco.pegoraro88@gmail.com>:

> Dear Ryan,
>
> thank you very much for your reply. I agree with you about the easier
> interpretation of the parametrization of Model 2.
> About the differences between Models 1 and 2, you are right, the
> dispersion parameters estimates are nearly identical but the p-values are
> very different between the tests. Probably, as you told me, the reason is I
> make some mistakes in the definition of the contrast.
> Thank you again. Have a good day,
> Marco
>
>
> 2014-01-29 Ryan <rct@thompsonclan.org>
>
> Hi Marco,
>>
>> Note that your models 2 & 3 are exactly the same model. The expression in
>> your formula for model 3:
>>
>> Genotype:Load:Time
>>
>> is more or less equivalent to:
>>
>> factor(paste(design$Genotype, design$Load, design$Time, sep=":"))
>>
>> except that the factor levels might possibly be in a different order. As
>> for your question about doing the same test with Models 1 & 2, I believe
>> that both models are equivalent parametrizations of the same design, so you
>> should get (nearly) identical dispersion estimates, and performing
>> equivalent tests with the two models should give you the same p-values. If
>> you are not getting the same p-values, then you have not properly set up
>> the contrast for Model 1. This isn't surprising, because a 3-way
>> interaction model is really confusing to work with, and I can't even tell
>> you what the correct contrast would be. I would recommend that you stick
>> with Model 2 (or 3), where each coefficient has a direct biological
>> interpretation as the estimated expression level of a group. Designing
>> contrasts for this parametrization is much easier, in my opinion. I think
>> you would agree, since you got this contrast correct in your example code.
>>
>> -Ryan
>>
>>
>> On 1/29/14, 7:11 AM, Marco [guest] wrote:
>>
>>> Hi all!
>>>
>>> I have a question about the way to define the model matrix for a
>>> multi-factorial design in which I want to model all the interactions. In
>>> particular, I have a 3x3x2 factorial design (3 factors) and three
>>> biological replicates, so I have 3x3x2x3 = 54 sample.
>>> The three factors are:
>>> - GENOTYPE with the three levels F, G and P
>>> - LOAD with levels HIGH, MEDIUM and LOW
>>> - TIME with levels H and PH
>>>
>>> So this is the design matrix:
>>>
>>>  design
>>>>
>>>         Time   Load Genotype
>>> FHH_1     H   High        F
>>> FHPH_1   PH   High        F
>>> FHH_2     H   High        F
>>> FHPH_2   PH   High        F
>>> FHH_3     H   High        F
>>> FHPH_3   PH   High        F
>>> FMH_1     H Medium        F
>>> FMPH_1   PH Medium        F
>>> FMH_2     H Medium        F
>>> FMPH_2   PH Medium        F
>>> FMH_3     H Medium        F
>>> FMPH_3   PH Medium        F
>>> FLH_1     H    Low        F
>>> FLPH_1   PH    Low        F
>>> FLH_2     H    Low        F
>>> FLPH_2   PH    Low        F
>>> FLH_3     H    Low        F
>>> FLPH_3   PH    Low        F
>>> GHH_1     H   High        G
>>> GHPH_1   PH   High        G
>>> GHH_2     H   High        G
>>> GHPH_2   PH   High        G
>>> GHH_3     H   High        G
>>> GHPH_3   PH   High        G
>>> GMH_1     H Medium        G
>>> GMPH_1   PH Medium        G
>>> GMH_2     H Medium        G
>>> GMPH_2   PH Medium        G
>>> GMH_3     H Medium        G
>>> GMPH_3   PH Medium        G
>>> GLH_1     H    Low        G
>>> GLPH_1   PH    Low        G
>>> GLH_2     H    Low        G
>>> GLPH_2   PH    Low        G
>>> GLH_3     H    Low        G
>>> GLPH_3   PH    Low        G
>>> PHH_1     H   High        P
>>> PHPH_1   PH   High        P
>>> PHH_2     H   High        P
>>> PHPH_2   PH   High        P
>>> PHH_3     H   High        P
>>> PHPH_3   PH   High        P
>>> PMH_1     H Medium        P
>>> PMPH_1   PH Medium        P
>>> PMH_2     H Medium        P
>>> PMPH_2   PH Medium        P
>>> PMH_3     H Medium        P
>>> PMPH_3   PH Medium        P
>>> PLH_1     H    Low        P
>>> PLPH_1   PH    Low        P
>>> PLH_2     H    Low        P
>>> PLPH_2   PH    Low        P
>>> PLH_3     H    Low        P
>>> PLPH_3   PH    Low        P
>>>
>>>
>>> I am interested in finding all possible interactions between the three
>>> factors, so I want to model also the three way interactions
>>> GENOTYPE:LOAD:TIME. I have already read the edgeR vignette so I was
>>> wondering if the best way to define the model matrix is by defining it with
>>> a formula or by defining each treatment combination as a group.
>>> I'm trying to understand which are the differences between the two
>>> alternatives, so I tried to estimate both the models:
>>>
>>> Model 1)
>>>    mydesign1<-model.matrix(~Genotype * Load * Time, data=design)
>>>    D1 <- estimateGLMCommonDisp(D, mydesign1) # D is the DGEList of the
>>> normalized counts
>>>         D1 <- estimateGLMTrendedDisp(D1, mydesign1)
>>>         fit1 <- glmFit(D1, mydesign1)
>>>
>>> The estimated coefficients are:
>>>         > colnames(fit1)
>>>          [1] "(Intercept)"                 "GenotypeG"
>>>          [3] "GenotypeP"                   "LoadLow"
>>>          [5] "LoadMedium"                  "TimePH"
>>>          [7] "GenotypeG:LoadLow"           "GenotypeP:LoadLow"
>>>          [9] "GenotypeG:LoadMedium"        "GenotypeP:LoadMedium"
>>>         [11] "GenotypeG:TimePH"            "GenotypeP:TimePH"
>>>         [13] "LoadLow:TimePH"              "LoadMedium:TimePH"
>>>         [15] "GenotypeG:LoadLow:TimePH"    "GenotypeP:LoadLow:TimePH"
>>>         [17] "GenotypeG:LoadMedium:TimePH" "GenotypeP:LoadMedium:TimePH"
>>>       Model 2)
>>>    Group <- factor(paste(design$Genotype, design$Load, design$Time,
>>> sep="."))
>>>    # I have 18 unique levels in Group (54 samples / 3 biological
>>> replicates)
>>>
>>>    mydesign2<-model.matrix(~0+Group)
>>>    colnames(mydesign2) <- levels(Group)
>>>
>>>         D2 <- estimateGLMCommonDisp(D, mydesign2)
>>>         D2 <- estimateGLMTrendedDisp(D2, mydesign2)
>>>
>>>         fit2<-glmFit(D2, mydesign2)
>>>         > colnames(fit2)
>>>          [1] "H.High.F"    "H.High.G"    "H.High.P"    "H.Low.F"
>>> "H.Low.G"
>>>          [6] "H.Low.P"     "H.Medium.F"  "H.Medium.G"  "H.Medium.P"
>>>  "PH.High.F"
>>>         [11] "PH.High.G"   "PH.High.P"   "PH.Low.F"    "PH.Low.G"
>>>  "PH.Low.P"
>>>         [16] "PH.Medium.F" "PH.Medium.G" "PH.Medium.P"
>>>    Now, to understand the differences between the two models, my
>>> question was this:
>>> if I define the contrast:
>>>
>>>    contrast1<-c(rep(0,14), 1, 0, -1, 0) # "GenotypeG:LoadLow:TimePH" vs
>>> "GenotypeG:LoadMedium:TimePH"
>>>
>>> on the first model and I make a LR test, will I obtain the same result
>>> as making a LR test on the contrast:
>>>
>>>    contrast2<-makeContrasts(G.LowvsMedium.PH = PH.Low.G - PH.Medium.G,
>>> levels=mydesign2)
>>>
>>> on the second model?
>>> This is the code I used to make the tests:
>>>
>>>    lrt1.G.LowvsMedium.PH<-glmLRT(fit1, contrast=contrast1)
>>>    summary(decideTestsDGE(lrt1.G.LowvsMedium.PH))
>>>    and
>>>
>>>    lrt2.G.LowvsMedium.PH<-glmLRT(fit2, contrast=contrast2[,"
>>> G.LowvsMedium.PH"])
>>>    summary(decideTestsDGE(lrt2.G.LowvsMedium.PH))
>>>
>>> The answer to the question is no, because I've tried to do it on my data
>>> and the results are different between the two tests. But why? Am I doing
>>> two completely different tests? The two parameters of the two models aren't
>>> saying me the same thing?
>>>
>>> In a moment of madness :) I tried to estimate another model, with only
>>> the three-level interactions between the factors.
>>>
>>> Model 3)
>>>    mydesign3<-model.matrix(~0+Genotype:Load:Time, data=design)
>>>
>>>    D3 <- estimateGLMCommonDisp(D, mydesign3)
>>>    D3 <- estimateGLMTrendedDisp(D3, mydesign3)
>>>
>>>    fit3<-glmFit(D3, mydesign3)
>>>    > colnames(fit3)
>>>    [1] "GenotypeF:LoadHigh:TimeH"    "GenotypeG:LoadHigh:TimeH"
>>>    [3] "GenotypeP:LoadHigh:TimeH"    "GenotypeF:LoadLow:TimeH"
>>>    [5] "GenotypeG:LoadLow:TimeH"     "GenotypeP:LoadLow:TimeH"
>>>    [7] "GenotypeF:LoadMedium:TimeH"  "GenotypeG:LoadMedium:TimeH"
>>>    [9] "GenotypeP:LoadMedium:TimeH"  "GenotypeF:LoadHigh:TimePH"
>>>    [11] "GenotypeG:LoadHigh:TimePH"   "GenotypeP:LoadHigh:TimePH"
>>>    [13] "GenotypeF:LoadLow:TimePH"    "GenotypeG:LoadLow:TimePH"
>>>    [15] "GenotypeP:LoadLow:TimePH"    "GenotypeF:LoadMedium:TimePH"
>>>    [17] "GenotypeG:LoadMedium:TimePH" "GenotypeP:LoadMedium:TimePH"
>>>    Well, this model gave me the same parameter estimates as the model 2
>>> (each treatment combination as a group) and, obviously, the same results in
>>> terms of DEG. But is it statistically correct to estimate a model with
>>> three-level interactions without any principal effect?
>>>
>>> Thanks in advance for any replies!
>>>
>>> Marco
>>>
>>>   -- output of sessionInfo():
>>>
>>> R version 3.0.2 (2013-09-25)
>>> Platform: x86_64-w64-mingw32/x64 (64-bit)
>>>
>>> locale:
>>> [1] LC_COLLATE=Italian_Italy.1252  LC_CTYPE=Italian_Italy.1252
>>> [3] LC_MONETARY=Italian_Italy.1252 LC_NUMERIC=C
>>> [5] LC_TIME=Italian_Italy.1252
>>>
>>> attached base packages:
>>> [1] splines   parallel  stats     graphics  grDevices utils     datasets
>>> [8] methods   base
>>>
>>> other attached packages:
>>>   [1] EDASeq_1.8.0         ShortRead_1.20.0     Rsamtools_1.14.2
>>>   [4] lattice_0.20-24      Biostrings_2.30.1    GenomicRanges_1.14.4
>>>   [7] XVector_0.2.0        IRanges_1.20.6       Biobase_2.22.0
>>> [10] BiocGenerics_0.8.0   edgeR_3.4.2          limma_3.18.9
>>> [13] aroma.light_1.32.0   matrixStats_0.8.14   yeastRNASeq_0.0.10
>>>
>>> loaded via a namespace (and not attached):
>>>   [1] annotate_1.40.0      AnnotationDbi_1.24.0 bitops_1.0-6
>>>   [4] DBI_0.2-7            DESeq_1.14.0         genefilter_1.44.0
>>>   [7] geneplotter_1.40.0   grid_3.0.2           hwriter_1.3
>>> [10] latticeExtra_0.6-26  R.methodsS3_1.6.1    R.oo_1.17.0
>>> [13] RColorBrewer_1.0-5   RSQLite_0.11.4       stats4_3.0.2
>>> [16] survival_2.37-6      tools_3.0.2          XML_3.98-1.1
>>> [19] xtable_1.7-1         zlibbioc_1.8.0
>>>
>>> --
>>> Sent via the guest posting facility at bioconductor.org.
>>>
>>> _______________________________________________
>>> Bioconductor mailing list
>>> Bioconductor@r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>>> Search the archives: http://news.gmane.org/gmane.
>>> science.biology.informatics.conductor
>>>
>>
>>
>

	[[alternative HTML version deleted]]

