Using eBayes to find P values for individual contrasts
3
0
Entering edit mode
@jason-shoemaker-4357
Last seen 10.2 years ago
Dear all, I have searched the archives but not found any advice on this issue. I am using the LIMMA package to determine differentially expressed genes. I have been using eBayes plus topTable to find the most differentially expressed genes, but now I want to determine the adjusted p values specific for each contrast. Should I simply do as follows (following the example from http://matticklab.com/index.php?title=Single_channel_analysis_of_Agile nt_microarray_data_with_Limma): contrast.matrix <- makeContrasts("Treatment1-Treatment2", "Treatment1-Treatment3", "Treatment2-Treatment1", levels=design); fit2 <- contrasts.fit(fit, contrast.matrix) fit2 <- eBayes(fit2) P.values<-p.adjust(fit2$p.values,methods="fdr"); in doing so- can I make fair comparisons between p values for each contrast? Or more precisely, if a get a p value of 0.01 for "Treatment1-Treatment2" and large value (P>0.1) for the remaining 2 contrasts, is that gene significant only for "Treatment1-Treatment2"? Thank you! Jason
limma limma • 5.0k views
ADD COMMENT
0
Entering edit mode
Mark Cowley ▴ 910
@mark-cowley-2951
Last seen 10.2 years ago
Hi Jason, I think you're in danger of reinventing the wheel. The adj.P.Val column in the topTable is the corrected p value. Don't forget about the coef topTable parameter to control which coefficient to look at. You can control what method to use via the adjust.method parameter. then take a look at the decideTests method to work out which genes are significant for which contrasts. cheers, mark On 16/11/2010, at 7:28 PM, Jason Shoemaker wrote: > Dear all, > > I have searched the archives but not found any advice on this issue. I am using the LIMMA package to determine differentially expressed genes. I have been using eBayes plus topTable to find the most differentially expressed genes, but now I want to determine the adjusted p values specific for each contrast. Should I simply do as follows (following the example from http://matticklab.com/index.php?ti tle=Single_channel_analysis_of_Agilent_microarray_data_with_Limma): > > contrast.matrix <- makeContrasts("Treatment1-Treatment2", "Treatment1-Treatment3", "Treatment2-Treatment1", levels=design); > fit2 <- contrasts.fit(fit, contrast.matrix) > fit2 <- eBayes(fit2) > > P.values<-p.adjust(fit2$p.values,methods="fdr"); > > in doing so- can I make fair comparisons between p values for each contrast? Or more precisely, if a get a p value of 0.01 for "Treatment1-Treatment2" and large value (P>0.1) for the remaining 2 contrasts, is that gene significant only for "Treatment1-Treatment2"? > Thank you! > Jason > > _______________________________________________ > 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 ----------------------------------------------------- Mark Cowley, PhD Peter Wills Bioinformatics Centre Garvan Institute of Medical Research, Sydney, Australia
ADD COMMENT
0
Entering edit mode
Dear Mark, Thank you for the warning! I was worried about asking something silly. So if I may ask, how can I get topTable to display not just a single adjusted p value for one contrast, but the adiusted P values for all contrasts? I don't seem to see this option. Thus I have been applying p.adjust to the raw P values to adjust the values for each contrast of interest. Thank you! Jason On 11/17/2010 10:32 AM, Mark Cowley wrote: > Hi Jason, > I think you're in danger of reinventing the wheel. > > The adj.P.Val column in the topTable is the corrected p value. Don't forget about the coef topTable parameter to control which coefficient to look at. You can control what method to use via the adjust.method parameter. > > then take a look at the decideTests method to work out which genes are significant for which contrasts. > > cheers, > mark > > On 16/11/2010, at 7:28 PM, Jason Shoemaker wrote: > >> Dear all, >> >> I have searched the archives but not found any advice on this issue. I am using the LIMMA package to determine differentially expressed genes. I have been using eBayes plus topTable to find the most differentially expressed genes, but now I want to determine the adjusted p values specific for each contrast. Should I simply do as follows (following the example from http://matticklab.com/index.php?ti tle=Single_channel_analysis_of_Agilent_microarray_data_with_Limma): >> >> contrast.matrix<- makeContrasts("Treatment1-Treatment2", "Treatment1-Treatment3", "Treatment2-Treatment1", levels=design); >> fit2<- contrasts.fit(fit, contrast.matrix) >> fit2<- eBayes(fit2) >> >> P.values<-p.adjust(fit2$p.values,methods="fdr"); >> >> in doing so- can I make fair comparisons between p values for each contrast? Or more precisely, if a get a p value of 0.01 for "Treatment1-Treatment2" and large value (P>0.1) for the remaining 2 contrasts, is that gene significant only for "Treatment1-Treatment2"? >> Thank you! >> Jason >> >> _______________________________________________ >> 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 > > > ----------------------------------------------------- > Mark Cowley, PhD > > Peter Wills Bioinformatics Centre > Garvan Institute of Medical Research, Sydney, Australia > ----------------------------------------------------- > >
ADD REPLY
0
Entering edit mode
Hi Jason, In short, topTable can only give you the adjusted p-values for a single contrast at a time (or a series of contrasts, but it's calculating the overall F-value, not individual t-values). Instead, see write.fit(). This only writes out to a file, but I often just read it back in to R. You could also just do this: indiv.P.values<-apply(fit2$p.values, 2, p.adjust, method="fdr"); Cheers, Jenny At 08:58 PM 11/17/2010, Jason Shoemaker wrote: >Dear Mark, > >Thank you for the warning! I was worried about asking something >silly. So if I may ask, how can I get topTable to display not just a >single adjusted p value for one contrast, but the adiusted P values >for all contrasts? I don't seem to see this option. Thus I have been >applying p.adjust to the raw P values to adjust the values for each >contrast of interest. > >Thank you! >Jason > >On 11/17/2010 10:32 AM, Mark Cowley wrote: >>Hi Jason, >>I think you're in danger of reinventing the wheel. >> >>The adj.P.Val column in the topTable is the corrected p value. >>Don't forget about the coef topTable parameter to control which >>coefficient to look at. You can control what method to use via the >>adjust.method parameter. >> >>then take a look at the decideTests method to work out which genes >>are significant for which contrasts. >> >>cheers, >>mark >> >>On 16/11/2010, at 7:28 PM, Jason Shoemaker wrote: >> >>>Dear all, >>> >>>I have searched the archives but not found any advice on this >>>issue. I am using the LIMMA package to determine differentially >>>expressed genes. I have been using eBayes plus topTable to find >>>the most differentially expressed genes, but now I want to >>>determine the adjusted p values specific for each contrast. Should >>>I simply do as follows (following the example from >>>http://matticklab.com/index.php?title=Single_channel_analysis_of_Ag ilent_microarray_data_with_Limma): >>> >>>contrast.matrix<- makeContrasts("Treatment1-Treatment2", >>>"Treatment1-Treatment3", "Treatment2-Treatment1", levels=design); >>>fit2<- contrasts.fit(fit, contrast.matrix) >>>fit2<- eBayes(fit2) >>> >>>P.values<-p.adjust(fit2$p.values,methods="fdr"); >>> >>>in doing so- can I make fair comparisons between p values for each >>>contrast? Or more precisely, if a get a p value of 0.01 for >>>"Treatment1-Treatment2" and large value (P>0.1) for the remaining >>>2 contrasts, is that gene significant only for "Treatment1-Treatment2"? >>>Thank you! >>>Jason >>> >>>_______________________________________________ >>>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 >> >> >>----------------------------------------------------- >>Mark Cowley, PhD >> >>Peter Wills Bioinformatics Centre >>Garvan Institute of Medical Research, Sydney, Australia >>----------------------------------------------------- >> > >_______________________________________________ >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
ADD REPLY
0
Entering edit mode
Dear all, Great! Thanks for all the advice. I was doing exactly as Jenny recommended, but I've switched the code to simply cycle through the topTable coefficients and stack the results into a large data frame (as Mark recommended). I also played with the decideTests() which proofed useful in constructed a scenario plots (#genes significant for each contrast, recommended by Sean). Thank you all! Jason On 11/19/2010 6:56 AM, Jenny Drnevich wrote: > Hi Jason, > > In short, topTable can only give you the > adjusted p-values for a single contrast at a > time (or a series of contrasts, but it's > calculating the overall F-value, not individual > t-values). Instead, see write.fit(). This only > writes out to a file, but I often just read it > back in to R. You could also just do this: > > indiv.P.values<-apply(fit2$p.values, 2, > p.adjust, method="fdr"); > > Cheers, > Jenny > > At 08:58 PM 11/17/2010, Jason Shoemaker wrote: >> Dear Mark, >> >> Thank you for the warning! I was worried about >> asking something silly. So if I may ask, how >> can I get topTable to display not just a single >> adjusted p value for one contrast, but the >> adiusted P values for all contrasts? I don't >> seem to see this option. Thus I have been >> applying p.adjust to the raw P values to adjust >> the values for each contrast of interest. >> >> Thank you! >> Jason >> On 11/17/2010 10:32 AM, Mark Cowley wrote: >>> Hi Jason, >>> I think you're in danger of reinventing the >>> wheel. >>> >>> The adj.P.Val column in the topTable is the >>> corrected p value. Don't forget about the coef >>> topTable parameter to control which >>> coefficient to look at. You can control what >>> method to use via the adjust.method parameter. >>> >>> then take a look at the decideTests method to >>> work out which genes are significant for >>> which contrasts. >>> >>> cheers, >>> mark >>> >>> On 16/11/2010, at 7:28 PM, Jason Shoemaker wrote: >>> >>>> Dear all, >>>> >>>> I have searched the archives but not found >>>> any advice on this issue. I am using the >>>> LIMMA package to determine differentially >>>> expressed genes. I have been using eBayes >>>> plus topTable to find the most differentially >>>> expressed genes, but now I want to determine >>>> the adjusted p values specific for each >>>> contrast. Should I simply do as follows >>>> (following the example from >>>> http://matticklab.com/index.php?title=Single_channel_analysis_of_ Agilent_microarray_data_with_Limma): >>>> >>>> contrast.matrix<- >>>> makeContrasts("Treatment1-Treatment2", >>>> "Treatment1-Treatment3", >>>> "Treatment2-Treatment1", levels=design); >>>> fit2<- contrasts.fit(fit, contrast.matrix) >>>> fit2<- eBayes(fit2) >>>> >>>> P.values<-p.adjust(fit2$p.values,methods="fdr"); >>>> >>>> in doing so- can I make fair comparisons >>>> between p values for each contrast? Or more >>>> precisely, if a get a p value of 0.01 for >>>> "Treatment1-Treatment2" and large value >>>> (P>0.1) for the remaining 2 contrasts, is >>>> that gene significant only for >>>> "Treatment1-Treatment2"? >>>> Thank you! >>>> Jason >>>> >>>> _______________________________________________ >>>> 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 >>>> >>> >>> >>> ----------------------------------------------------- >>> >>> Mark Cowley, PhD >>> >>> Peter Wills Bioinformatics Centre >>> Garvan Institute of Medical Research, Sydney, >>> Australia >>> ----------------------------------------------------- >>> >>> >> >> _______________________________________________ >> 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 >> > >
ADD REPLY
0
Entering edit mode
Mark Cowley ▴ 910
@mark-cowley-2951
Last seen 10.2 years ago
Hi Jason, you normally generate one topTable per contrast of interest. so if you have 3 coefficients of interest, you would call topTable 3 times, each time choosing coef=1, then coef=2, then coef=3 in 3 I usually just write a for loop for each coef of interest & store the topTables in a list. cheers, Mark On 18/11/2010, at 1:58 PM, Jason Shoemaker wrote: > Dear Mark, > > Thank you for the warning! I was worried about asking something silly. So if I may ask, how can I get topTable to display not just a single adjusted p value for one contrast, but the adiusted P values for all contrasts? I don't seem to see this option. Thus I have been applying p.adjust to the raw P values to adjust the values for each contrast of interest. > > Thank you! > Jason > > On 11/17/2010 10:32 AM, Mark Cowley wrote: >> Hi Jason, >> I think you're in danger of reinventing the wheel. >> >> The adj.P.Val column in the topTable is the corrected p value. Don't forget about the coef topTable parameter to control which coefficient to look at. You can control what method to use via the adjust.method parameter. >> >> then take a look at the decideTests method to work out which genes are significant for which contrasts. >> >> cheers, >> mark >> >> On 16/11/2010, at 7:28 PM, Jason Shoemaker wrote: >> >>> Dear all, >>> >>> I have searched the archives but not found any advice on this issue. I am using the LIMMA package to determine differentially expressed genes. I have been using eBayes plus topTable to find the most differentially expressed genes, but now I want to determine the adjusted p values specific for each contrast. Should I simply do as follows (following the example from http://matticklab.com/index.php?ti tle=Single_channel_analysis_of_Agilent_microarray_data_with_Limma): >>> >>> contrast.matrix<- makeContrasts("Treatment1-Treatment2", "Treatment1-Treatment3", "Treatment2-Treatment1", levels=design); >>> fit2<- contrasts.fit(fit, contrast.matrix) >>> fit2<- eBayes(fit2) >>> >>> P.values<-p.adjust(fit2$p.values,methods="fdr"); >>> >>> in doing so- can I make fair comparisons between p values for each contrast? Or more precisely, if a get a p value of 0.01 for "Treatment1-Treatment2" and large value (P>0.1) for the remaining 2 contrasts, is that gene significant only for "Treatment1-Treatment2"? >>> Thank you! >>> Jason >>> >>> _______________________________________________ >>> 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 >> >> >> ----------------------------------------------------- >> Mark Cowley, PhD >> >> Peter Wills Bioinformatics Centre >> Garvan Institute of Medical Research, Sydney, Australia >> ----------------------------------------------------- >> >> >
ADD COMMENT
0
Entering edit mode
Mark Cowley ▴ 910
@mark-cowley-2951
Last seen 10.2 years ago
Hi Jason, you normally generate one topTable per contrast of interest. so if you have 3 coefficients of interest, you would call topTable 3 times, each time choosing coef=1, then coef=2, then coef=3 in 3 I usually just write a for loop for each coef of interest & store the topTables in a list. cheers, Mark On 18/11/2010, at 1:58 PM, Jason Shoemaker wrote: > Dear Mark, > > Thank you for the warning! I was worried about asking something > silly. So if I may ask, how can I get topTable to display not just a > single adjusted p value for one contrast, but the adiusted P values > for all contrasts? I don't seem to see this option. Thus I have been > applying p.adjust to the raw P values to adjust the values for each > contrast of interest. > > Thank you! > Jason > > On 11/17/2010 10:32 AM, Mark Cowley wrote: >> Hi Jason, >> I think you're in danger of reinventing the wheel. >> >> The adj.P.Val column in the topTable is the corrected p value. >> Don't forget about the coef topTable parameter to control which >> coefficient to look at. You can control what method to use via the >> adjust.method parameter. >> >> then take a look at the decideTests method to work out which genes >> are significant for which contrasts. >> >> cheers, >> mark >> >> On 16/11/2010, at 7:28 PM, Jason Shoemaker wrote: >> >>> Dear all, >>> >>> I have searched the archives but not found any advice on this >>> issue. I am using the LIMMA package to determine differentially >>> expressed genes. I have been using eBayes plus topTable to find >>> the most differentially expressed genes, but now I want to >>> determine the adjusted p values specific for each contrast. Should >>> I simply do as follows (following the example from http://mattickl ab.com/index.php?title=Single_channel_analysis_of_Agilent_microarray_d ata_with_Limma) >>> : >>> >>> contrast.matrix<- makeContrasts("Treatment1-Treatment2", >>> "Treatment1-Treatment3", "Treatment2-Treatment1", levels=design); >>> fit2<- contrasts.fit(fit, contrast.matrix) >>> fit2<- eBayes(fit2) >>> >>> P.values<-p.adjust(fit2$p.values,methods="fdr"); >>> >>> in doing so- can I make fair comparisons between p values for each >>> contrast? Or more precisely, if a get a p value of 0.01 for >>> "Treatment1-Treatment2" and large value (P>0.1) for the remaining >>> 2 contrasts, is that gene significant only for "Treatment1- >>> Treatment2"? >>> Thank you! >>> Jason >>> >>> _______________________________________________ >>> 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 >> >> >> ----------------------------------------------------- >> Mark Cowley, PhD >> >> Peter Wills Bioinformatics Centre >> Garvan Institute of Medical Research, Sydney, Australia >> ----------------------------------------------------- >> >> >
ADD COMMENT
0
Entering edit mode
On Thu, Nov 18, 2010 at 7:07 AM, Mark Cowley <m.cowley@garvan.org.au> wrote: > Hi Jason, > you normally generate one topTable per contrast of interest. so if you have > 3 coefficients of interest, you would call topTable 3 times, each time > choosing coef=1, then coef=2, then coef=3 in 3 > > I usually just write a for loop for each coef of interest & store the > topTables in a list. > > cheers, > Mark > On 18/11/2010, at 1:58 PM, Jason Shoemaker wrote: > > Dear Mark, >> >> Thank you for the warning! I was worried about asking something silly. So >> if I may ask, how can I get topTable to display not just a single adjusted p >> value for one contrast, but the adiusted P values for all contrasts? I don't >> seem to see this option. Thus I have been applying p.adjust to the raw P >> values to adjust the values for each contrast of interest. >> >> Hi, Jason. Be sure to take a look at the decideTests() function that Mark mentioned in one of his previous replies. I can't tell exactly what you are after, but a common use case of needing to decide what contrasts are significant when multiple contrasts are being evaluated simultaneously is handled by decideTests(). Sean > Thank you! >> Jason >> >> On 11/17/2010 10:32 AM, Mark Cowley wrote: >> >>> Hi Jason, >>> I think you're in danger of reinventing the wheel. >>> >>> The adj.P.Val column in the topTable is the corrected p value. Don't >>> forget about the coef topTable parameter to control which coefficient to >>> look at. You can control what method to use via the adjust.method parameter. >>> >>> then take a look at the decideTests method to work out which genes are >>> significant for which contrasts. >>> >>> cheers, >>> mark >>> >>> On 16/11/2010, at 7:28 PM, Jason Shoemaker wrote: >>> >>> Dear all, >>>> >>>> I have searched the archives but not found any advice on this issue. I >>>> am using the LIMMA package to determine differentially expressed genes. I >>>> have been using eBayes plus topTable to find the most differentially >>>> expressed genes, but now I want to determine the adjusted p values specific >>>> for each contrast. Should I simply do as follows (following the example from >>>> >>>> http://matticklab.com/index.php?title=Single_channel_analysis_of_ Agilent_microarray_data_with_Limma >>>> ): >>>> >>>> contrast.matrix<- makeContrasts("Treatment1-Treatment2", >>>> "Treatment1-Treatment3", "Treatment2-Treatment1", levels=design); >>>> fit2<- contrasts.fit(fit, contrast.matrix) >>>> fit2<- eBayes(fit2) >>>> >>>> P.values<-p.adjust(fit2$p.values,methods="fdr"); >>>> >>>> in doing so- can I make fair comparisons between p values for each >>>> contrast? Or more precisely, if a get a p value of 0.01 for >>>> "Treatment1-Treatment2" and large value (P>0.1) for the remaining 2 >>>> contrasts, is that gene significant only for "Treatment1-Treatment2"? >>>> Thank you! >>>> Jason >>>> >>>> _______________________________________________ >>>> 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 >>>> >>> >>> >>> ----------------------------------------------------- >>> Mark Cowley, PhD >>> >>> Peter Wills Bioinformatics Centre >>> Garvan Institute of Medical Research, Sydney, Australia >>> ----------------------------------------------------- >>> >>> >>> >> > _______________________________________________ > 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 REPLY

Login before adding your answer.

Traffic: 1059 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6