[BioC] Questions about GOseq after DESeq2

[Bernd Klaus]

Dear Amadine,

I guess the most important thing in an enrichment analysis is to choose a
proper "background" set of genes. ("measured genes" in the GOSeq parlance).

In order to choose an appropriate set, you could adapt the code below.
It first extracts the differentially expressed genes and then uses the function
matchit to compile a background set that has a comparable distribution.

This way, you can be sure that your background does not have a dramatically
different distribution from your DE genes.

Best wishes,

Bernd


### extract differentially expressed genes

genes <- rownames(subset(DESeqResults, pval.adj <0.1))

#Now, we generate a background set of genes matched for 
#expression strength, avoiding potential biases.


  library(MatchIt)
  backM <- c()

  ## prepare data frame for matching, sign indicates wheather 
  ## the gene is differentially expressed or not
  df <- data.frame( 
    sign=as.numeric( rownames(DESeqResults)  %in% as.character(genes)), 
    IBSres["baseMean"])


  df$baseMean <- round( df$baseMean, 0)
  

  ## repeat matching multiple times since
  ## each differentially expressed gene is 
  ## matched by exactly one non-expressed gene 

    for( i in 1:3 ){
    mm <- matchit( sign ~ baseMean, df, 
      method="nearest", distance="mahalanobis")
    backM <- c(backM, mm$match.matrix[,1])
    df <- df[which(!rownames(df) %in% backM),]
   }
  
  backM <- unique( na.exclude(backM) )
 ## total number of matched genes
  length(backM )
   
 ## no DE genes in Background:
  intersect(backM, genes)


#----------------------------------------------------------
###### Compare distributions of background and significant genes
#----------------------------------------------------------

#We check whether our background has the same distribution
#of expression strength. 



 pdf("matching.pdf" , height = 5 , width = 5)

  library(geneplotter)
  multidensity( list( 
       all= DESeqResults[,"baseMean"] ,
       fore=DESeqResults[rownames(DESeqResults) %in% genes, "baseMean"], 
       back=DESeqResults[rownames(DESeqResults) %in% backM, "baseMean"]), 
     xlab="mean counts", xlim=c(0, 150))
  dev.off()

################################

On Feb 20, 2014, at 3:03 PM, <amandine.fournier at chu-lyon.fr> wrote:

> 
> Dear Bioconductor users and developers,
> 
> I have 2 questions regarding a gene ontology analysis of RNA-Seq data.
> 
> I first used the DESeq2 package to perform the differential analysis of my data and it resulted in a list of about 500 DEG. Now I would like to perform a gene ontology analysis on this dataset with the GOseq package.
> 
> According to the GOseq user manual, this package requires 2 pieces of information : the differentially expressed genes (= my list of 500 genes) and the measured genes. These are defined as "all genes for which RNA-seq data was gathered for your experiment".
> 
> As you known, DESeq2 filters out low counting genes and outliers. So I wonder if I should consider :
> - all sequenced genes including low expressed genes and outliers, 
> - or the sequenced genes whithout outliers (but including low expressed genes)
> - or only remaining genes after low count filtering (final assayed genes) ?
> 
> I don't know which step filters outliers out. It seems that they are taking into account to estimate library sizes, but I am not sure if they are filtered before or after estimating dispersion values.
> Low expressed genes are filtered out at the end of the analysis, juste before calculating FDR, so they are taking into account for calculating library sizes, estimating dispersion values, fitting the GLM model and testing, but not in the final FDR.
> In brief, I wonder if outliers and low expressed genes can be considered as "assayed" in the DESeq2 analysis.
> 
> Do you have an opinion about this ? Is there a commonly accepted / advised methodology for that ?
> 
> The second question is : should I analyse over-expressed genes and under-expressed genes together or separately ?
> 
> Many thanks in advance.
> Best regards,
> Amandine
> 
> -----
> Amandine Fournier
> Lyon Neuroscience Research Center
> & Lyon Civil Hospital (France)
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor