This is a test version of Biostars. For the public version, visit https://www.biostars.org.
GSEA from broad institute normalization method

Hi freinds

I have HT-Seq raw count data of 19000 coding genes and 300 samples of two groups (treatment and control). I want to do gene set enrichment analysis with GSEA from broad institute. do I need to normalize my data before GSEA? How should I normalize for this purpose?

Thanks

rna-seq
• 6,010 views • link
updated by benformatics benformatics • written by Rob • Rob

1 answer

Hi,

As input to the Broad Institute's GSEA program, you should use any type of expression data that is [properly] normalised such that cross-sample differences can be faithfully gauged. This can mean using any of these:

  • normalised RNA-seq counts via DESeq2's 'geometric' normalisation, EdgeR's TMM method, et cetera.
  • normalised + transformed RNA-seq expression levels, such as variance-stabilised (vst) or regularised log (rlog) expression levels from DESeq2, or log2 CPMs from EdgeR
  • normalised microarray data via RMA, GC-RMA, MAS5, neqc, et cetera

This does not mean raw counts or any of these types of expression levels: FPKM, RPKM, TPM, et cetera

More information here:

Kevin

• 0 views • link
written by Kevin Blighe Kevin Blighe

Hi I am trying to normalize my count data to do GSEA. these error come up: what is the solution?

data <- read.csv("myData.csv")
data <- matrix(data)
vst(data)

Error in vst(data) : less than 'nsub' rows,
  it is recommended to use varianceStabilizingTransformation directly


rlog(data)
Error in DESeqDataSet(se, design = design, ignoreRank) : 
  'list' object cannot be coerced to type 'double'
• 0 views • link
written by Rob • Rob

Assuming that your data that is held in data is raw counts, you should first normalise this via:

dds <- DESeqDataSetFromMatrix(
  countData = data,
  colData = coldata,
  design= ~ Group)
dds <- DESeq(dds)

data does not have to be a data matrix.

coldata should be a data frame that represents the metadata for data, and its rows should be perfectly aligned with the columns of data. You should have at leas one column, in this case 'group', that represents treatment and control

Then we can transform these normalised counts:

varStabilised <- vst(data, blind = FALSE)
regularisedLog <- rlog(data, blind = FALSE)

The normalsied + transformed expression levels will then be accessible via:

assay(varStabilised)
assay(regularisedLog)
• 0 views • link
written by Kevin Blighe Kevin Blighe

Thanks Kevin What I understand from your answer is : first I do DESeq work flow with the raw count data, then the result (dds) as the normalized count will come to vst() or rlog():

varStabilised <- vst(dds, blind = FALSE)

Then this varStabilised will be my normalized data for GSEA. Am I getting this correctly?

• 0 views • link
written by Rob • Rob

Yep, but you can also use the normalised counts, accessible via:

counts(dds, normalized = TRUE)

The distributions of the normalised counts and that of the variance stabilised expression levels differ, but this is not a problem due to the fact that GSEA is based on ranking.

• 0 views • link
written by Kevin Blighe Kevin Blighe

Thanks for explanation When I used this I get error:

counts(dds, normalized = TRUE)
Error in .local(object, ...) : 
  first calculate size factors, add normalizationFactors, or set normalized=FALSE

how can I solve this?

• 0 views • link
written by Rob • Rob

You have evidently not yet run DESeq(dds)

• 0 views • link
written by Kevin Blighe Kevin Blighe

Thanks Kevin. it was helpful.

• 0 views • link
written by Rob • Rob

Log in to answer this question.

More posts like this
  • analyze gene expression of 3 genes in RNA-Seq data
    written by pinheirofabiano • pinheirofabiano
    hello, I want to analyze the gene expression of three genes in a RNA-Seq dataset. Can I perform Mann-Whitney test using the raw counts table …
  • RNA-seq, STAR count, GSEA
    written by Rob • Rob
    Hi friends I have about 300 differentially expressed genes for my data of STAR count. This is for two groups: treatment & control. However, for …
  • GSEA RNASeq
    written by Rob • Rob
    Hi friends For gene set enrichment analysis (GSEA), the software from broad institute does not accept ensemble IDs, I want to do the analysis using …
  • RNA-seq GSEA
    written by Rob • Rob
    Hi friends I have 120 patients and 19000 genes and I have differential expression analysis with DESeq2 in R. I found 213 significant genes. but …
  • Choosing the right pathways among gene set enrichment analysis results
    written by jabbari.parnian • jabbari.parnian
    Hello everyone, I have performed gene set enrichment analysis on my RNA-seq data using the GSEA software by Broad Institute. Depending on the comparisons and …
  • differential expression- coding genes
    written by Rob • Rob
    hello friends, I performed differential expression analysis with DESeq2 package. my data is HT-seq count data and protein coding genes. in the result TUNAR gene …
  • GSEA broad institute
    written by Rob • Rob
    Hello friends, I have been doing GDEA with the broad institute software. At the end, the result does not show any significant hallmark based on …
  • Gene set enrichment analysis
    written by Rob • Rob
    Hi friends I am doing gene set enrichment analysis with GSEA java base software from broad institute. I have 44 patients (2 groups: 22 control …
  • HT-Seq count data coding gene
    written by Rob • Rob
    Hi friends I am using R code to download HT-Seq count data Does anybody know what line should I add to the code to download …
  • Exclude non-coding genes
    written by Rob • Rob
    Hello, I have HT-Seq read count with around 60000 genes (with ENSEMBLE ID). I only want to keep protein coding genes. How can I remove …