Hello all,

I've attempted to use GOseq, and while I know that, due to the gene length correction, it is more appropriate for RNAseq data than topGO, I like that topGO takes the topology of the GO graph in to account. I'm trying to find a middle ground.

In his great great tutorial (http://www-huber.embl.de/users/klaus/Teaching/DESeq2Predoc2014.html#gene-ontology-enrichment-analysis), Bernd Klaus, uses the genefinder function to match 10 genes of similar expression levels (mean normalized count) to his DE genes of interest. He then uses this set as a background set to proceed with enrichment using topGO. I liked this approach, though I have an RNAseq data set with 4 different treatments with varied gene expression among the four and so I thought that matching based on average gene count might yield some inappropriate matches. I then thought couldn't I just use this matching method to pick genes with similar length I by that control for that bias similarly to GOseq. 

Does this make sense? 

Another question is, in the tutorial we search for 10 gene matches. Is this enough? I have a specific set of 70+ genes I'm interested in GO terms for and figured the background matching would only be 700 genes and that wouldn't be so representative of the whole set. Thoughts?

Any comments or suggestions are greatly appreciated.

-Ben

My code based on the tutorial:

geneLengthMatrix<-as.matrix(geneLengthData)

backG <- genefinder(geneLengthMatrix, rownames(DEgeneSet),numResults = 100, method = "manhattan")

backG <- rownames(geneLengthMatrix)[as.vector(sapply(backG, function(x)x$indices))]

backG <- setdiff(backG, rownames(DEgeneSet))

library(geneplotter)

multidensity( list( 

      all= log2(geneLengthMatrix[,1]),

      foreground =log2(geneLengthMatrix[rownames(DEgeneSet),]), 

      background = log2(geneLengthMatrix[backG,])), 

      xlab="log2 gene length", main = "Matching for enrichment analysis")

geneIDs = rownames(dds)

inUniverse = geneIDs %in% c(rownames(DEgeneSet),  backG) 

inSelection =  geneIDs %in% rownames(DEgeneSet) 

alg <- factor( as.integer( inSelection[inUniverse] ) )

names(alg) <- geneIDs[inUniverse]



tgd<-new("topGOdata", ontology="BP", allGenes=alg, nodeSize=ns, annot= annFUN.gene2GO, gene2GO= CukeGO)

Hi Ben,

thank you for your positive feedback, great that you like my tutorial!

1.) Matching by actual expresssion:

I simply used the basemeans for matching and 10 matches to speed up the whole process and since it gave reasonable results. This workflow was intended for teaching and thus needed to run relatively quickly. 

For the Fisher-test bases GO analysis that topGO performs, it is best to have as many genes in the background as possible in order to have maximal power, having just 700 in the background is probably a bit too few genes.

Also, you can provide an actual expression matrix of your genes to genefinder function, in order to match the genes by their expression pattern across samples, rather than just by base mean.

So you could try to run something like this:

# do rlog to have homoscedastic data

rld_mat <- assay(rlog(dds, blind = FALSE))

de_idx <- which(rownames(rld_mat) %in% de_genes)

genefinder(rld_mat, de_idx, 100, method = "manhattan")

2.) Matching by gene length

if you rather would want to match by gene length, you can adapt the following code to find the background genes

using the MatchIt package, just replace the gene mean count by the gene length.

library(MatchIt)
background <- c()

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

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

## repeat matching multiple times since
## each differentially expressed gene is
## matched by exactly one non-expressed gene
## you can run this as often as you want ...

for( i in 1:3 ){
  match_res <- matchit( sign ~ baseMean, df,
                 method="nearest", distance="mahalanobis")
  background <- c(background, match_res$match.matrix[,1])
  df <- df[which(!rownames(df) %in% background),]
}

background <- unique( na.exclude(background) )
## total number of matched de_genes
length(background )

## no DE genes in Background:
intersect(background, de_genes)

Hope that helps, feel free to ask for clarifications if needed,

Bernd