Limma: Top 10% gene list based on p-values
2
0
Entering edit mode
Voke AO ▴ 760
@voke-ao-4830
Last seen 9.6 years ago
Hi all, Thanks to Dan, I now know how to get genes highlighted in a Volcano plot based on p-values from the results of my limma analysis. I would like to get a list of the genes highlighted. For example, the code below gives me a volcano plot with the top 10% genes based on p-values highlighted. windows() pdf("VolcanoPlot_GSExxxxxx_10perc.pdf"); results$threshold = as.factor(results$P.Value<=quantile(results$P.Value, 0.1)) g = ggplot(data=results, aes(x=logFC, y=-log10(P.Value), colour=threshold)) + geom_point(alpha=0.4, size=1.75) + opts(legend.position = "none") + xlim(c(-5, 5)) + ylim(c(0, 8)) + xlab("log2 fold change") + ylab("-log10 p-value") g dev.off() This code below gives me the exact 10% threshold value: results = topTable(fit2, adjust ="BH", number = nrow(gse25724dat[[1]]), sort.by="P") results2 = quantile(results$P.Value, 0.1) What I would like however is a list/table of the genes that meet this criteria, i.e, the genes highlighted from the first code. Any help will be greatly appreciated. Thanks. -Avoks
limma limma • 1.7k views
ADD COMMENT
0
Entering edit mode
@sean-davis-490
Last seen 3 months ago
United States
On Wed, May 9, 2012 at 8:21 AM, Ovokeraye Achinike-Oduaran < ovokeraye@gmail.com> wrote: > Hi all, > > Thanks to Dan, I now know how to get genes highlighted in a Volcano > plot based on p-values from the results of my limma analysis. I would > like to get a list of the genes highlighted. For example, the code > below gives me a volcano plot with the top 10% genes based on p-values > highlighted. > > windows() > pdf("VolcanoPlot_GSExxxxxx_10perc.pdf"); > results$threshold = as.factor(results$P.Value<=quantile(results$P.Value, > 0.1)) > g = ggplot(data=results, aes(x=logFC, y=-log10(P.Value), > colour=threshold)) + > geom_point(alpha=0.4, size=1.75) + > opts(legend.position = "none") + > xlim(c(-5, 5)) + ylim(c(0, 8)) + > xlab("log2 fold change") + ylab("-log10 p-value") > g > dev.off() > > This code below gives me the exact 10% threshold value: > > results = topTable(fit2, adjust ="BH", number = > nrow(gse25724dat[[1]]), sort.by="P") > results2 = quantile(results$P.Value, 0.1) > > What I would like however is a list/table of the genes that meet this > criteria, i.e, the genes highlighted from the first code. Any help > will be greatly appreciated. > > top10pct = results[results$P.Value<=results2,] As an aside, I noticed that you are working with unadjusted p-values. I would not recommend doing so and would suggest using an adjusted p-value instead in most cases. Sean [[alternative HTML version deleted]]
ADD COMMENT
0
Entering edit mode
Voke AO ▴ 760
@voke-ao-4830
Last seen 9.6 years ago
Thanks Sean. Duly noted. I understand that the unadjusted p-values are not corrected for multiple testing. But thanks a bunch. -Avoks On Wed, May 9, 2012 at 2:21 PM, Ovokeraye Achinike-Oduaran <ovokeraye at="" gmail.com=""> wrote: > Hi all, > > Thanks to Dan, I now know how to get genes highlighted in a Volcano > plot based on p-values from the results of my limma analysis. I would > like to get a list of the genes highlighted. For example, the code > below gives me a volcano plot with the top 10% genes based on p-values > highlighted. > > windows() > pdf("VolcanoPlot_GSExxxxxx_10perc.pdf"); > results$threshold = as.factor(results$P.Value<=quantile(results$P.Value, 0.1)) > g = ggplot(data=results, aes(x=logFC, y=-log10(P.Value), colour=threshold)) + > ?geom_point(alpha=0.4, size=1.75) + > ?opts(legend.position = "none") + > ?xlim(c(-5, 5)) + ylim(c(0, 8)) + > ?xlab("log2 fold change") + ylab("-log10 p-value") > g > dev.off() > > This code below gives me the exact 10% threshold value: > > results = topTable(fit2, adjust ="BH", number = > nrow(gse25724dat[[1]]), sort.by="P") > results2 = quantile(results$P.Value, 0.1) > > What I would like however is a list/table of the genes that meet this > criteria, i.e, the genes highlighted from the first code. Any help > will be greatly appreciated. > > Thanks. > > -Avoks
ADD COMMENT

Login before adding your answer.

Traffic: 451 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