Dear List,

I have an experiment with 3 WT (GBT) and 3 KO (GAT) mice, each with 4 time points: 0h,2h,8h,24h.
The PI is interested in the overall effect of KO vs WT (at each time point and overall effect too). From the PCA plot I could see that time has a large effect compared to the groups.

I thus generated differential expression analysis in the following manner:
------
>s.info
   Sample_Group Group Time     GT MousePair Chip Rep
1         GAT0h    KO    0  KO_0h         1    1   1
8         GAT0h    KO    0  KO_0h         4    2   2
13        GAT0h    KO    0  KO_0h         5    3   3
3         GAT2h    KO    2  KO_2h         1    1   1
10        GAT2h    KO    2  KO_2h         4    2   2
19        GAT2h    KO    2  KO_2h         5    4   3
5         GAT8h    KO    8  KO_8h         1    1   1
16        GAT8h    KO    8  KO_8h         4    3   2
21        GAT8h    KO    8  KO_8h         5    4   3
11       GAT24h    KO   24 KO_24h         1    2   1
18       GAT24h    KO   24 KO_24h         4    3   2
23       GAT24h    KO   24 KO_24h         5    4   3
2         GBT0h    WT    0  WT_0h         2    1   1
7         GBT0h    WT    0  WT_0h         3    2   2
14        GBT0h    WT    0  WT_0h         6    3   3
4         GBT2h    WT    2  WT_2h         2    1   1
9         GBT2h    WT    2  WT_2h         3    2   2
20        GBT2h    WT    2  WT_2h         6    4   3
6         GBT8h    WT    8  WT_8h         2    1   1
15        GBT8h    WT    8  WT_8h         3    3   2
22        GBT8h    WT    8  WT_8h         6    4   3
12       GBT24h    WT   24 WT_24h         2    2   1
17       GBT24h    WT   24 WT_24h         3    3   2
24       GBT24h    WT   24 WT_24h         6    4   3

>group = s.info$GT
>design = model.matrix(~0+group)
>colnames(design) =  gsub("group", "", colnames(design))
>fit <- lmFit(sig.norm, design)
>contrasts.KW <- makeContrasts(KW.T0 = KO_0h-WT_0h,
                              KW.T2 = KO_2h-WT_2h,
                              KW.T8 = KO_8h-WT_8h,
                              KW.T24 = KO_24h-WT_24h,
                              IntT2vsT0=(KO_2h-KO_0h)-(WT_2h-WT_0h),
                              IntT8vsT2=(KO_8h-KO_2h)-(WT_8h-WT_2h),
                              IntT24vsT8=(KO_24h-KO_8h)-(WT_24h-WT_8h),
                              levels=design)
>fit2 <- contrasts.fit(fit, contrasts.KW)
>fit2 <- eBayes(fit2)
>colnames(fit2$coefficients)#[1] "KW.T0"      "KW.T2"      "KW.T8"      "KW.T24"     "IntT2vsT0"#[6] "IntT8vsT2"  "IntT24vsT8"

>KW_0h <- topTable(fit2, coef="KW.T0", sort.by="P", adjust='fdr', number=nrow(sig.norm))
>KW_2h <- topTable(fit2, coef="KW.T2", sort.by="P", adjust='fdr', number=nrow(sig.norm))
>KW_8h <- topTable(fit2, coef="KW.T8", sort.by="P", adjust='fdr', number=nrow(sig.norm))
>KW_24h <- topTable(fit2, coef="KW.T24", sort.by="P", adjust='fdr', number=nrow(sig.norm))
>KW.f <- topTableF(fit2, adjust="fdr", number=nrow(sig.norm))

>results <- decideTests(fit2, method="nestedF") # adjust F-test p-values
>summary(results)
#   KW.T0 KW.T2 KW.T8 KW.T24 IntT2vsT0 IntT8vsT2 IntT24vsT8
#-1     4     4     4      5         0         0          0
#0  22564 22564 22564  22563     22568     22568      22568
#1      0     0     0      0         0         0          0
------

However, I think that my F test above is possibly incorrect, because the design for that should have pair as an added covariate (The reason being from each mouse there are 4 time points, thus this would become a paired analysis). Therefore, with this is mind I then do the following:

> pair = as.factor(s.info$MousePair)
>design2 = model.matrix(~0+group+pair)
>colnames(design2) =  gsub("group", "", colnames(design2))

>fit.a <- lmFit(sig.norm, design2)
#Coefficients not estimable: pair6
#Warning message:
#Partial NA coefficients for 22568 probe(s)
>contrasts.KW.b <- makeContrasts(IntT2vsT0=(KO_2h-KO_0h)-(WT_2h-WT_0h),
                              IntT8vsT2=(KO_8h-KO_2h)-(WT_8h-WT_2h),
                              IntT24vsT8=(KO_24h-KO_8h)-(WT_24h-WT_8h),
                              levels=design2)

>fit.b <- contrasts.fit(fit.a, contrasts.KW.b)
>fit.b <- eBayes(fit.b)
>colnames(fit.b$coefficients) #[1] "IntT2vsT0"  "IntT8vsT2"  "IntT24vsT8"

>KW.f2 <- topTableF(fit.b, adjust="fdr", number=nrow(sig.norm))
#Error in topTableF(fit.b, adjust = "fdr", number = nrow(sig.norm)) :
#  F-statistics not available. Try topTable for individual coef instead.

>results2 <- decideTests(fit.b, method="nestedF")
#Warning message:
#In is.na(object$F.p.value) :
#  is.na() applied to non-(list or vector) of type 'NULL'
>summary(results2)
#Warning message:
#In is.na(object$F.p.value) :
#  is.na() applied to non-(list or vector) of type 'NULL'


Now, I am confused as to which of the above is the right way to go about the analysis. Perhaps I have overlooked something basic and misunderstood.
Any advice and suggestion would be much appreciated.

Many Thanks,
Natasha



--


	[[alternative HTML version deleted]]

