Question: p.adj and logFC
0
gravatar for Yisong Zhen
10.8 years ago by
Yisong Zhen200
Yisong Zhen200 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?) 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]]
limma • 757 views
ADD COMMENTlink modified 10.8 years ago by James W. MacDonald50k • written 10.8 years ago by Yisong Zhen200
Answer: p.adj and logFC
0
gravatar for James W. MacDonald
10.8 years ago by
United States
James W. MacDonald50k 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. Best, Jim > > 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
ADD COMMENTlink written 10.8 years ago by James W. MacDonald50k
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 > [[alternative HTML version deleted]]
ADD REPLYlink written 10.8 years ago by Sean Davis21k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 247 users visited in the last hour