Question

goana limma- extract list of DE genes and genes in the enriched GO terms?

2

Entering edit mode

carl.mann ▴ 20

@carlmann-8378

Last seen 8.8 years ago

France

From the S3 object returned by the goana() function in limma, is it possible to extract the list of DE genes that are enriched for each GO term (gene lists corresponding to N and DE in the returned data.frame)?

limma goana • 6.2k views

ADD COMMENT • link updated 8.8 years ago by Gordon Smyth 50k • written 8.8 years ago by carl.mann ▴ 20

score 6 · Accepted Answer · 2015-07-11

6

Entering edit mode

Aaron Lun ★ 28k

@alun

Last seen 7 hours ago

The city by the bay

No, I'm pretty sure that information isn't stored in the output. To get it, you'll have to do something like this:

# Basic set-up:
DB <- paste("org", species, "eg", "db", sep = ".")
require(DB, character.only = TRUE)
GO2ALLEGS <- paste("org", species, "egGO2ALLEGS", sep = ".")
EG.GO <- AnnotationDbi::toTable(get(GO2ALLEGS))
d <- duplicated(EG.GO[, c("gene_id", "go_id", "Ontology")])
EG.GO <- EG.GO[!d, ]

# Now choosing DE genes for each GO ID/Ontology combination:
de.by.go <- split(EG.GO$gene_id, paste(EG.GO$go_id, EG.GO$Ontology, sep="."))
de.by.go <- lapply(de.by.go, FUN=function(x) { x[x %in% de] })

I'm assuming you have some species identifier in species, and a vector of Entrez Gene IDs in de. I'm also assuming that you haven't set universe in your original call to goana.

ADD COMMENT • link 8.8 years ago Aaron Lun ★ 28k

0

Entering edit mode

Hello Aaron- yes, species = H.s., vector of Entrez Gene IDs for de, and universe not set in call to goana. I'll give your code a try. Thanks for your prompt help! Carl.

ADD REPLY • link 8.8 years ago carl.mann ▴ 20

0

Entering edit mode

Hi Aaron,

I am also trying to extract a list of DE genes from my GO term analysis. However, I am quite new to using R and I don't completely follow your code/ cannot make it work. Would it be possible to break it down a bit more?

Really hope you can help me out!

Thanks, Lotte

ADD REPLY • link 5.8 years ago Lotte • 0

0

Entering edit mode

This post is 3 years old. Start a new post if you have a new question.
Try reading the documentation. ?duplicated, ?toTable, ?split, etc.
I'm not psychic, I can't know what you don't understand if you don't mention it.

ADD REPLY • link 5.8 years ago Aaron Lun ★ 28k

0

Entering edit mode

Hi Aaron,

Thank you for your very kind answer.

ADD REPLY • link 5.8 years ago Lotte • 0

score 5 · Accepted Answer · 2015-07-14

Here's another way if you're only interested in a few GO Terms. Suppose you're interested in GO:0001516. Then

> library(org.Hs.eg.db)
> x <- org.Hs.egGO2ALLEGS
> Rkeys(x) <- "GO:0001516"
> EG <- mappedLkeys(x)
> EG
 [1] "301"    "551"    "552"    "873"    "972"    "1906"   "1907"   "3248"  
 [9] "4282"   "5601"   "5730"   "5740"   "5742"   "5743"   "6916"   "8644"  
[17] "8877"   "9536"   "10728"  "22949"  "23411"  "27306"  "50640"  "80142" 
[25] "127281" "255189"

gives you the Entrez Gene IDs corresponding to N and

> intersect(EG, your.de.list)

gives you Entrez Gene IDs corresponding to DE.