Thanks for your reply. My understanding for the p.adjust from the R help page is that it is adjusted by number of comparison. I mean, if I filter genes (below certain expression threshold) before I carry out multiple testing, the p.adjust would be decreased. But how do we set the threshold for p.adjust, <=0.05? If so, the MKG.56.34.1_at (probe in our array comparison, de facto decreasing expressed gene) will not be statistically significant, although its fold changes reach about 6. Then how do we interpret this result? Thanks million times. Yisong On Fri, Sep 19, 2008 at 7:23 AM, Sean Davis <sdavis2@mail.nih.gov> wrote: > > > On Fri, Sep 19, 2008 at 10:12 AM, James W. MacDonald < > jmacdon@med.umich.edu> wrote: > >> Hi Yisong, >> >> Yisong Zhen wrote: >> >>> Hello All, >>> I used the following script to find the genes that both have p.adj < 0.05 >>> and logFC < -1 in comparisons (TCdn - TCwt and Tcestw - TCwt). >>> I checked one candidate gene in the normalized expression output >>> (affy_My_data.txt). The first two are TCdn, middle two are TCwt and the >>> last >>> two are TCestw): >>> >>> MKG.56.34.1_at 4.30886639379254 5.33637926862047 11.8385942887500 >>> 10.9736438003031 4.22713527087133 3.28741332008267 >>> >>> It seems that this gene have same logFC in two comparisons (TCdn~ TCwt >>> and >>> TCets~TCwt). But the p.adj is very different: one is 0.053759(TCdn~TCwt) >>> and >>> the other is 0.009286 (TCetsw~TCwt). >>> >>> I am not familiar with LIMMA model and the code is simply from >>> Biocoductor >>> manual, please give some advice and tell me how to interpret this >>> result.(How the p.adj is calculated? Why the same gene in two comparisons >>> have same logFC but the p.adj is so different?) >>> >> >> The reason you get a different p-value with the same log FC is because the >> p-value is based on the t-statistic, not the fold change (although the fold >> change is the numerator of that statistic). >> >> If you look at the t-statistics for these two genes you will see that they >> are different. >> > > And even if the t-statistics were identical, the adjustment method you have > chosen ('fdr') is a function of a LIST of genes, and NOT the gene itself. > For details, you can read the help for p.adjust(). > > Sean > > >> >> >> >> >>> Thanks. >>> >>> Yisong >>> >>> >>> >>> library(limma) >>> library(affy) >>> targets<-readTargets("targets.txt") >>> #targets; >>> data<-ReadAffy(filenames=targets$FileName) >>> eset<-rma(data) >>> write.exprs(eset, file="affy_My_data.txt") >>> design<-model.matrix(~-1+factor(c(1,1,2,2,3,3))) >>> colnames(design)<-c("TCdn","TCwt","TCetsw") >>> fit<-lmFit(eset,design) >>> contrast.matrix<-makeContrasts(TCdn-TCwt,TCetsw-TCwt, levels=design) >>> fit2<-contrasts.fit(fit, contrast.matrix) >>> fit2<-eBayes(fit2) >>> write.table(topTable(fit2, coef=1, adjust="fdr", sort.by="B", >>> number=50000), >>> file="limmadn_wt.xls", row.names=F, sep="\t") >>> write.table(topTable(fit2, coef=2, adjust="fdr", sort.by="B", >>> number=50000), >>> file="limmaestw_wt.xls", row.names=F, sep="\t") >>> >>> [[alternative HTML version deleted]] >>> >>> _______________________________________________ >>> Bioconductor mailing list >>> Bioconductor@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 >> Hildebrandt Lab >> 8220D MSRB III >> 1150 W. Medical Center Drive >> Ann Arbor MI 48109-0646 >> 734-936-8662 >> >> >> _______________________________________________ >> Bioconductor mailing list >> Bioconductor@stat.math.ethz.ch >> https://stat.ethz.ch/mailman/listinfo/bioconductor >> Search the archives: >> http://news.gmane.org/gmane.science.biology.informatics.conductor >> > > -- //-------------------------------------------------------------------- --------------------- // We have a hunger of the mind which asks for knowledge // of all around us, and the more we gain, the more is // our desire; the more we see, the more we are capable // of seeing. //-------------------------------------------------------------------- -------------------- Yisong Zhen, Ph.D (Davidson Lab) Molecular & Cellular Biology Molecular Cardiovascular Research Program University of Arizona 1656 E. Mabel, MRB 317 Tucson, AZ 85724 USA Lab: 520-626-8153 Cell: 520-977-7397 @--------------------------------------------------------------------- ------------------ [[alternative HTML version deleted]]

On Fri, Oct 3, 2008 at 7:58 PM, Yisong Zhen <zhenyisong at="" gmail.com=""> wrote: > Thanks for your reply. > > My understanding for the p.adjust from the R help page is that it is > adjusted by number of comparison. > I mean, if I filter genes (below certain expression threshold) before I > carry out multiple testing, the p.adjust would be decreased. > But how do we set the threshold for p.adjust, <=0.05? p.adjust has an argument for the type of adjustment. You need to specify that. If you specify 'fdr" for false discovery, then it makes sense to use a larger cutoff for significance, depending on your needs. Keep in mind that an fdr is interpreted with respect to a LIST OF GENES. So a false discovery rate of 0.2, for example, means that 20% of the genes that occur in the list up to that point are likely false discoveries. You can see how a cutoff of 0.05 might not make sense in this situation, depending on your needs. > If so, the MKG.56.34.1_at (probe in our array comparison, de facto > decreasing expressed gene) will not be statistically significant, although > its fold changes reach about 6. Then how do we interpret this result? It is possible to have a high fold change without statistical significance if the variability in the measure is also high. Sean > On Fri, Sep 19, 2008 at 7:23 AM, Sean Davis <sdavis2 at="" mail.nih.gov=""> wrote: > >> >> >> On Fri, Sep 19, 2008 at 10:12 AM, James W. MacDonald < >> jmacdon at med.umich.edu> wrote: >> >>> Hi Yisong, >>> >>> Yisong Zhen wrote: >>> >>>> Hello All, >>>> I used the following script to find the genes that both have p.adj < 0.05 >>>> and logFC < -1 in comparisons (TCdn - TCwt and Tcestw - TCwt). >>>> I checked one candidate gene in the normalized expression output >>>> (affy_My_data.txt). The first two are TCdn, middle two are TCwt and the >>>> last >>>> two are TCestw): >>>> >>>> MKG.56.34.1_at 4.30886639379254 5.33637926862047 11.8385942887500 >>>> 10.9736438003031 4.22713527087133 3.28741332008267 >>>> >>>> It seems that this gene have same logFC in two comparisons (TCdn~ TCwt >>>> and >>>> TCets~TCwt). But the p.adj is very different: one is 0.053759(TCdn~TCwt) >>>> and >>>> the other is 0.009286 (TCetsw~TCwt). >>>> >>>> I am not familiar with LIMMA model and the code is simply from >>>> Biocoductor >>>> manual, please give some advice and tell me how to interpret this >>>> result.(How the p.adj is calculated? Why the same gene in two comparisons >>>> have same logFC but the p.adj is so different?) >>>> >>> >>> The reason you get a different p-value with the same log FC is because the >>> p-value is based on the t-statistic, not the fold change (although the fold >>> change is the numerator of that statistic). >>> >>> If you look at the t-statistics for these two genes you will see that they >>> are different. >>> >> >> And even if the t-statistics were identical, the adjustment method you have >> chosen ('fdr') is a function of a LIST of genes, and NOT the gene itself. >> For details, you can read the help for p.adjust(). >> >> Sean >> >> >>> >>> >>> >>> >>>> Thanks. >>>> >>>> Yisong >>>> >>>> >>>> >>>> library(limma) >>>> library(affy) >>>> targets<-readTargets("targets.txt") >>>> #targets; >>>> data<-ReadAffy(filenames=targets$FileName) >>>> eset<-rma(data) >>>> write.exprs(eset, file="affy_My_data.txt") >>>> design<-model.matrix(~-1+factor(c(1,1,2,2,3,3))) >>>> colnames(design)<-c("TCdn","TCwt","TCetsw") >>>> fit<-lmFit(eset,design) >>>> contrast.matrix<-makeContrasts(TCdn-TCwt,TCetsw-TCwt, levels=design) >>>> fit2<-contrasts.fit(fit, contrast.matrix) >>>> fit2<-eBayes(fit2) >>>> write.table(topTable(fit2, coef=1, adjust="fdr", sort.by="B", >>>> number=50000), >>>> file="limmadn_wt.xls", row.names=F, sep="\t") >>>> write.table(topTable(fit2, coef=2, adjust="fdr", sort.by="B", >>>> number=50000), >>>> file="limmaestw_wt.xls", row.names=F, sep="\t") >>>> >>>> [[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 >>> Hildebrandt Lab >>> 8220D MSRB III >>> 1150 W. Medical Center Drive >>> Ann Arbor MI 48109-0646 >>> 734-936-8662 >>> >>> >>> _______________________________________________ >>> 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 >>> >> >> > > > -- > //------------------------------------------------------------------ ----------------------- > // We have a hunger of the mind which asks for knowledge > // of all around us, and the more we gain, the more is > // our desire; the more we see, the more we are capable > // of seeing. > //------------------------------------------------------------------ ---------------------- > > Yisong Zhen, Ph.D (Davidson Lab) > Molecular & Cellular Biology > Molecular Cardiovascular Research Program > University of Arizona > 1656 E. Mabel, MRB 317 > Tucson, AZ 85724 USA > Lab: 520-626-8153 > Cell: 520-977-7397 > > @------------------------------------------------------------------- -------------------- > > [[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 >