[BioC] Microarray analyses. How to know the direction of the expression?
James W. MacDonald
jmacdon at uw.edu
Thu Sep 12 20:57:14 CEST 2013
Hi Mary,
This is another reason to use a 'cell means' parameterization rather
than a 'factor effects' parameterization. This is just statistician
blahblah for computing a mean for each group rather than a baseline
level and a difference between the other sample and the baseline.
In your situation, since the second coefficient has WT in the name (and
you have an intercept term), the intercept is computing the mean of a
baseline sample (in your case the KO sample), and the GroupWT
coefficient is WT - KO.
You could instead do
design <- model.matrix(~0+Group)
contrast <- matrix(c(1,-1), ncol = 1)
colnames(contrast) <- "KO vs WT"
fit <- lmFit(myEset, design)
fit2 <- contrasts.fit(fit, contrast)
fit2 <- eBayes(fit2)
topTable(fit2, 1)
Or you could use relevel()
Group <- relevel(Group, "WT")
and create your design matrix. But IMO it is always easier to just
compute group means and make comparisons directly, because then you
don't have to figure out how to interpret the coefficients.
Best,
Jim
On Thursday, September 12, 2013 1:13:09 PM, Mary Kindall wrote:
> Given a design matrix, how can I know that the logFC in the output of
> ebayesFit is WT/KO or KO/WT
>
> > Group
> [1] KO WT WT WT KO KO
> Levels: KO WT
> > design<- model.matrix(~Group)
> > design
> (Intercept) GroupWT
> 1 1 0
> 2 1 1
> 3 1 1
> 4 1 1
> 5 1 0
> 6 1 0
> attr(,"assign")
> [1] 0 1
> attr(,"contrasts")
> attr(,"contrasts")$Group
> [1] "contr.treatment"
>
> > fit <-lmFit(myEset,design)
> > fit<-eBayes(fit)
> > fit
> An object of class "MArrayLM"
> $coefficients
> (Intercept) GroupWT
> 1415670_at 10.0 -0.11544
> 1415671_at 10.9 0.41568
> 1415672_at 11.6 0.00933
> 1415673_at 10.4 -0.16966
> 1415674_a_at 10.5 -0.01095
> 18937 more rows ...
>
> $rank
> [1] 2
>
> $assign
> [1] 0 1
>
> $qr
> $qr
> (Intercept) GroupWT
> 1 -2.449 -1.225
> 2 0.408 -1.225
> 3 0.408 0.527
> 4 0.408 0.527
> 5 0.408 -0.290
> 6 0.408 -0.290
> attr(,"assign")
> [1] 0 1
> attr(,"contrasts")
> attr(,"contrasts")$Group
> [1] "contr.treatment"
>
>
> $qraux
> [1] 1.41 1.53
>
> $pivot
> [1] 1 2
>
> $tol
> [1] 1e-07
>
> $rank
> [1] 2
>
>
> $df.residual
> [1] 4 4 4 4 4
> 18937 more elements ...
>
> $sigma
> 1415670_at 1415671_at 1415672_at 1415673_at 1415674_a_at
> 0.0736 0.1407 0.0627 0.3145 0.1100
> 18937 more elements ...
>
> $cov.coefficients
> (Intercept) GroupWT
> (Intercept) 0.333 -0.333
> GroupWT -0.333 0.667
>
> $stdev.unscaled
> (Intercept) GroupWT
> 1415670_at 0.577 0.816
> 1415671_at 0.577 0.816
> 1415672_at 0.577 0.816
> 1415673_at 0.577 0.816
> 1415674_a_at 0.577 0.816
> 18937 more rows ...
>
> $pivot
> [1] 1 2
>
> $Amean
> 1415670_at 1415671_at 1415672_at 1415673_at 1415674_a_at
> 9.95 11.12 11.57 10.34 10.53
> 18937 more elements ...
>
> $method
> [1] "ls"
>
> $design
> (Intercept) GroupWT
> 1 1 0
> 2 1 1
> 3 1 1
> 4 1 1
> 5 1 0
> 6 1 0
> attr(,"assign")
> [1] 0 1
> attr(,"contrasts")
> attr(,"contrasts")$Group
> [1] "contr.treatment"
>
>
> $df.prior
> [1] 2.41
>
> $s2.prior
> [1] 0.0378
>
> $var.prior
> [1] 422.8 49.9
>
> $proportion
> [1] 0.01
>
> $s2.post
> 1415670_at 1415671_at 1415672_at 1415673_at 1415674_a_at
> 0.0176 0.0266 0.0167 0.0759 0.0218
> 18937 more elements ...
>
> $t
> (Intercept) GroupWT
> 1415670_at 130.6 -1.0653
> 1415671_at 115.9 3.1225
> 1415672_at 155.0 0.0885
> 1415673_at 65.5 -0.7540
> 1415674_a_at 123.6 -0.0909
> 18937 more rows ...
>
> $df.total
> [1] 6.41 6.41 6.41 6.41 6.41
> 18937 more elements ...
>
> $p.value
> (Intercept) GroupWT
> 1415670_at 3.16e-12 0.3252
> 1415671_at 6.78e-12 0.0188
> 1415672_at 1.05e-12 0.9322
> 1415673_at 2.62e-10 0.4776
> 1415674_a_at 4.49e-12 0.9304
> 18937 more rows ...
>
> $lods
> (Intercept) GroupWT
> 1415670_at 16.9 -6.16
> 1415671_at 16.6 -3.41
> 1415672_at 17.2 -6.75
> 1415673_at 14.4 -6.45
> 1415674_a_at 16.7 -6.75
> 18937 more rows ...
>
> $F
> [1] 16860 13961 24050 4225 15270
> 18937 more elements ...
>
> $F.p.value
> [1] 1.17e-12 2.15e-12 3.75e-13 9.89e-11 1.61e-12
> 18937 more elements ...
>
> > x <- topTable(fit,n=nrow(exprs(myEset)),adjust="BH",coef=2)
> > head(x)
> logFC AveExpr t P.Value adj.P.Val B
> 1452666_a_at 5.36 4.93 45.0 2.92e-09 5.54e-05 8.50
> 1441054_at 4.82 4.67 33.0 2.12e-08 2.01e-04 7.94
> 1424722_at 4.53 4.75 28.0 5.95e-08 3.62e-04 7.53
> 1435378_at 3.91 4.59 26.5 8.52e-08 3.62e-04 7.37
> 1428391_at 4.78 4.97 26.0 9.54e-08 3.62e-04 7.31
> 1419311_at 4.79 5.94 21.8 2.98e-07 9.42e-04 6.71
>
>
>
> Is all these genes here upregulated as compared to WT?
>
--
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