Hi Guido, Hooiveld, Guido wrote: > Hi James, > > OK, I got it. :) > > Let me, as biologist, rephrase what you said (correct me if I am wrong): > (email is also for the archive, since I do know some aspects of the > issue I had are also faced by fellow biologists, who, like me, are less > trained in statistics). > > > - key words are 'group-means parameterization'. > - my current design is: design <- model.matrix(~0+TS). The term '~0' > indictates that no intercept is calculated. [Thus I was on the correct > track... ;)] > - I defined these two contrasts: cont.matrix <- > makeContrasts(WTwyvWTc=WT.WY-WT.Con, KOwyvKOc=KO.WY-KO.Con, > levels=design). Doing so, for each of the two strains, I only tested for > the difference between treatment (WY) and control. The coefficients > reflect the difference between the treatment and control, which equals > the log2 of the FC. This is what I had seen before: > >> fit2 > An object of class "MArrayLM" > $coefficients > Contrasts > WTwyvWTc KOwyvKOc > 100009600_at -0.02321959 0.1853882 > 100012_at -0.11743549 -0.2804115 > 100017_at 0.41806549 0.2506005 > 100019_at 0.44521455 0.1060852 > 100034251_at -0.39616238 0.1285417 > 11512 more rows ... Oh yeah, my bad. If you want the group means, you want to look at the object you had before you fit the contrast. > > > - After your reply I redefined my contrast to explicitly state the first > 2 experimental groups: cont.matrix <- makeContrasts(WTWY=WT.WY, > WTCon=WT.Con, WTwyvWTc=WT.WY-WT.Con, KOwyvKOc=KO.WY-KO.Con, > levels=design). Since no intercept has been defined, I now also test > whether the expression of a gene in the WT-WY and WT-control groups is > significanty different compared to....? Compared to what I don't know > (overall mean, or zero??), but in my case this is not relevant. In the > same run I also test for the difference between treatment and control, > like I did before. The nice thing now is that the first two coeffecient > represent the average expression of the WT.WY and WT.Con groups, the > last two coefficients (again) the log2 of the FC between treatment vs > control. > >> fit2 > An object of class "MArrayLM" > $coefficients > Contrasts > WT.WY WT.Con WTwyvWTc KOwyvKOc > 100009600_at 2.418141 2.486590 -0.02321959 0.1853882 > 100012_at 1.992032 2.058794 -0.11743549 -0.2804115 > 100017_at 6.804809 6.888234 0.41806549 0.2506005 > 100019_at 3.028478 3.137888 0.44521455 0.1060852 > 100034251_at 4.182340 6.298939 -0.39616238 0.1285417 > 11512 more rows ... > > - However, regarding the standard deviation I am worried about all > identical values that are observed for 5 different probesets within a > same group. I have the feeling the stdev.unscaled are actually not the > 'real' SDs, but that I have to further manipulate the unscaled stdevs. > Thus, are the identical results below to be expected? If not, how to get > the standard deviation instead? > > $stdev.unscaled > Contrasts > WT.WY WT.Con WTwyvWTc KOwyvKOc > 100009600_at 0.5491619 0.478989 0.8439029 0.7600834 > 100012_at 0.5491619 0.478989 0.8439029 0.7600834 > 100017_at 0.5491619 0.478989 0.8439029 0.7600834 > 100019_at 0.5491619 0.478989 0.8439029 0.7600834 > 100034251_at 0.5491619 0.478989 0.8439029 0.7600834 > 11512 more rows ... Again, my bad. You can get the standard deviation appropriate for a t-statistic by fit2$stdev.unscaled/fit2$sigma, but this won't be the same as what you would get by apply(exprs(x.norm)[,1:n], 1, sd), which is what I think you are looking for. Best, Jim > > > Thanks & regards, > Guido > >> -----Original Message----- >> From: James W. MacDonald [mailto:jmacdon at med.umich.edu] >> Sent: 06 July 2009 21:45 >> To: Hooiveld, Guido >> Cc: bioconductor at stat.math.ethz.ch >> Subject: Re: [BioC] Limma: how to extract mean of groups? >> >> Hi Guido, >> >> Hooiveld, Guido wrote: >>> Hi James, >>> >>> Thanks. I indeed had seen that before, but I had the impression the >>> 'coefficients' actually are the log2 of the fold change...?!? >>> Anywayz, I'll have a closer look at it to make sure I understand it >>> correctly. >> It depends on how you set up the design matrix. If you have >> an intercept, then the coefficients will be the difference >> between a baseline level and whatever level you are looking >> at. However, if you don't have an intercept, then the >> coefficients represent the group means. >> >> Best, >> >> Jim >> >> >>> Best, >>> Guido >>> >>> >>> >>> >>>> -----Original Message----- >>>> From: James W. MacDonald [mailto:jmacdon at med.umich.edu] >>>> Sent: 06 July 2009 18:50 >>>> To: Hooiveld, Guido >>>> Cc: bioconductor at stat.math.ethz.ch >>>> Subject: Re: [BioC] Limma: how to extract mean of groups? >>>> >>>> Hi Guido, >>>> >>>> Hooiveld, Guido wrote: >>>>> Dear list, >>>>> >>>>> Is it possible to extract the group means + SD from >> Limma's output? >>>>> >>>>> I do know that the "Amean" can be extracted (fit2$Amean), >>>> but this is >>>>> the mean across all arrays/groups, and i am also >> interested in the >>>>> means of the experimental groups. >>>> Have you read the relevant help page (hint ?MArrayLM-class)? >>>> >>>> I get >>>> >>>> Slots/Components: >>>> >>>> 'MArrayLM' objects do not contain any slots (apart >> from '.Data') >>>> but they should contain the following list components: >>>> >>>> 'coefficients': 'matrix' containing fitted coefficients or >>>> contrasts >>>> >>>> 'stdev.unscaled': 'matrix' containing unscaled standard >>>> deviations >>>> of the coefficients or contrasts >>>> >>>> Which appears to be what you want. >>>> >>>> Best, >>>> >>>> Jim >>>> >>>> >>>>> Therefore, is it somehow possible to extact the mean (+ >> SD) of the >>>>> groups that are defined by the contrast matrix? Thus in my >>>> particular >>>>> case the mean + SD of the WT.Con, WT.WY. KO.Con and KO.WY groups. >>>>> >>>>> Thanks, >>>>> Guido >>>>> >>>>> >>>>> library(affy) >>>>> library(limma) >>>>> >>>>> targets <- readTargets("targets.txt") data <- >>>>> ReadAffy(filenames=targets$FileName) >>>>> x.norm <- rma(data) >>>>> >>>>> TS <- paste(targets$Strain, targets$Treatment, sep=".") TS <- >>>>> factor(TS, levels=c("WT.Con","WT.WY","KO.Con","KO.WY")) >>>>> design <- model.matrix(~0+TS) >>>>> colnames(design) <- levels(TS) >>>>> fit <- lmFit(eset, design) >>>>> cont.matrix <- makeContrasts(WTwyvWTc=WT.WY-WT.Con, >>>>> KOwyvKOc=KO.WY-KO.Con, levels=design) >>>>> fit2 <- contrasts.fit(fit, cont.matrix) >>>>> fit2 <- eBayes(fit2) >>>>> >>>>> ------------------------------------------------ >>>>> Guido Hooiveld, PhD >>>>> Nutrition, Metabolism & Genomics Group Division of Human >> Nutrition >>>>> Wageningen University Biotechnion, Bomenweg 2 >>>>> NL-6703 HD Wageningen >>>>> the Netherlands >>>>> tel: (+)31 317 485788 >>>>> fax: (+)31 317 483342 >>>>> internet: http://nutrigene.4t.com <http: nutrigene.4t.com=""/> >>>>> email: guido.hooiveld at wur.nl >>>>> >>>>> >>>>> >>>>> [[alternative HTML version deleted]] >>>>> >>>>> _______________________________________________ >>>>> Bioconductor mailing list >>>>> Bioconductor at stat.math.ethz.ch >>>>> https://stat.ethz.ch/mailman/listinfo/bioconductor >>>>> Search the archives: >>>>> http://news.gmane.org/gmane.science.biology.informatics.conductor >>>> -- >>>> 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 >>>> >>>> >> -- >> 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 >> >> > -- 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