This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Questions regarding GSEA - clusterProfiler

Good morning everyone,

I recently received some bulk RNA sequencing data : 3 controls and 3 treated (coming from mice). I identified the DEGs using the DESeq package. Out of the >17000 genes detected, approximately 150 genes had a pvalue which was inferior to 0.05. If I look at the number of genes with a padj value inferior to 0.05, it's even less than 5 genes.

However, I am fairly new with these analysis, and even if the DEG analysis seem to indicate that not a lot is happening in the presence of my treatment, I ran clusterProfiler. To perform GSEA, I use the entire gene list (>17000 genes) ranked based on the log2FoldChange. I see that I have several KEGG showing enrichment. For example, tyrosine metabolism is showing a NES of 2 with a pvalue at 3,00542E-05 and a padj at 3,00542E-05. Using gseKEGG gave me 54 results.

I was wondering if there was a way to incorporate the results of the pvalue ? If it was necessary to perform GSEA ? For example, should I only use the 150 genes with a pvalue inferior to 0.05 - but I read several posts saying that it was reducing the power of the analysis...

What am I doing wrong ? What would you advice ?

Thanks !

rna-seq gsea
• 2,953 views • link
updated by Irshad • Irshad • written by CaroH • CaroH

1 answer

One of the best ways to run GSEA is by using the test statistic from DESeq2 as the gene scoring metric. It takes into consideration the statistical strength of differential expression as well as the fold change. It will probably be less noisy as compared to the log2foldchange.

Regarding ORA enrichment with clusterprofiler, there really is no rule that you need to select the ones with FDR<0.05. You could easily set the threshold to FDR<0.1 or p<0.01, or take the top 1000 up and downregulated genes for enrichment. A similar approach has been described previously (here). Using less than 200 genes is not recommended for ORA as the sample is so small you are likely to have a low sensitivity.

• 0 views • link
written by mark.ziemann mark.ziemann

Hello. Thank you for your answer. What do you mean by the test statistic from DESeq2 as the gene scoring metric ?

• 0 views • link
written by CaroH • CaroH

The output DE table from DESeq has the following column headers: "baseMean", "log2FoldChange", "lfcSE", "stat", "pvalue" and "padj". The test statistic is the "stat" column

• 0 views • link
written by mark.ziemann mark.ziemann

Hi Mark, mark.ziemann

I am sorry to jump in here but there is something that is bugging me when I am doing the GSEA in clusterProfiler. I have used five different gene scoring metric for my analysis.

  1. rank_pvalue = log2foldchange * -log10(pvalue)
  2. rank_padj = log2foldchnage * -log10(padj)
  3. rank_sign_pval = sign(log2foldchange) * -log10(pvalue)
  4. rank_log2FoldChange = log2FoldChange column from DESeq
  5. rank_stat = stat column from DESeq as you have described here

The results for the first two doesn't return any GO term. For the third, I only get one GO term. While for the four and fifth I get numerous enriched GO terms using the following parameter.

gsea_res <- GSEA(
  geneList = gene_list,
  exponent = 1,
  minGSSize = 10,
  maxGSSize = 500,
  eps = 1e-10,
  TERM2GENE = term2gene,
  TERM2NAME = term2name,
  pvalueCutoff = 0.05,
  pAdjustMethod = "BH",
  by = "fgsea",
  verbose = TRUE, 
  seed = TRUE,

)

My question: If I use "stat", do you think that will be acceptable to the reviewers because I have seen a lot of folks using the first three methods.

• 0 views • link
written by Irshad • Irshad

Log in to answer this question.

More posts like this