Entering edit mode
Natasha
▴
440
@natasha-4640
Last seen 10.2 years ago
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.factors.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]]
