Question

Gene Onotology: Can't get the goana function to work

0

Entering edit mode

Adelyn ▴ 10

@47ceacfa

Last seen 2.1 years ago

United States

Hello,

Using the edgeR manual, I'm trying to follow the examples to run the goana function. I run the glmQLFTest (without any issues) and then use its output "qlf" as the input for goana or the topTags(qlf). I get the following errors.


> qlf <- glmQLFTest(fit,coef=2)
> go <- goana(qlf, species="Mm")
Error in goana.default(de = DEGenes, universe = universe, ...) : 
  No annotated genes found in universe
> 
> toptags <- topTags(qlf, n=50)
> go <- goana(toptags, species="Mm")
Error in goana.default(toptags, species = "Mm") : 
  components of de should be vectors
> TT <- as.vector(toptags)
> go <- goana(TT, species="Mm")
Error in goana.default(TT, species = "Mm") : 
  components of de should be vectors

> toptags
Coefficient:  groups2 
                        logFC    logCPM          F       PValue          FDR
ENSMUSG00000028195 -3.3561168  7.661106 1076.90952 6.348826e-22 8.937877e-18
ENSMUSG00000006445 -1.7990272  7.681729  995.09376 1.651467e-21 1.162468e-17
ENSMUSG00000032531 -2.8981796  6.854250  878.57035 7.424919e-21 2.730134e-17
ENSMUSG00000022178 -1.6372362  6.168699  875.38451 7.757165e-21 2.730134e-17
ENSMUSG00000028967 -1.9122235  6.359604  842.42713 1.231471e-20 3.293641e-17
ENSMUSG00000028525 -2.1673682  5.400544  833.31214 1.403740e-20 3.293641e-17
ENSMUSG00000019970 -1.9011901  7.168098  630.19507 3.991683e-19 8.027845e-16

I try to get them into the list of character vectors, but have difficulty converting them and get the same errors. My geneIDs are of the form ENSMUSG00000026596. Could that be causing my error?

I have followed other entries on the help forum to figure it out without any resolution. I'm fairly new to R and definitely new to this type of analysis; so, any help is greatly appreciated.

Thank you for your time, Adelyn

edgeR GO • 1.5k views

ADD COMMENT • link 2.7 years ago Adelyn ▴ 10

1

Entering edit mode

Finally got it to work! (seems that once you ask for help, you always figure it out)

My issue must have been in the ensembl IDs I adapted the code from the manual 4.1.3 and got the feature to run.

Here is that code for ensembl ID to EntrezIDs for a mouse genome:

>x <- read.delim("subread_counts_R_EnsemblID.txt")
> groups <- factor(c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,7,7,7,8,8,8))

> y <- DGEList(counts=x[,2:25], genes=x[1],group=groups)



> library(org.Mm.eg.db)

>idfound <- y$genes$EnsemblID %in% mappedRkeys(org.Mm.egENSEMBL)
> y <- y[idfound,]
> dim(y)
    [1] 32061    24
 > egENSEMBL <- toTable(org.Mm.egENSEMBL)
> head(egENSEMBL)
> m<- match(y$genes$EID, egENSEMBL$ensembl_id)
> y$genes$EntrezGene <- egENSEMBL$gene_id[m]
> keep <- filterByExpr(y)
> y <- y[keep,,keep.lib.sizes=FALSE]
> y <- calcNormFactors(y)
> design <- model.matrix(~groups)
> y <- estimateDisp(y,design)
> fit <- glmQLFit(y,design)
> qlf <- glmQLFTest(fit,coef=2)
> go <- goana(qlf)
> topGO(go, ont="BP", sort="Up", n=30, truncate=30)
                                     Term Ont   N Up Down         P.Up    P.Down
GO:0006335 DNA replication-dependent n...  BP   9  6    2 0.0008335504 0.4971445
GO:0034723 DNA replication-dependent n...  BP   9  6    2 0.0008335504 0.4971445
GO:0045653 negative regulation of mega...  BP   7  5    1 0.0015465112 0.7482637
GO:0045652 regulation of megakaryocyte...  BP  10  6    1 0.0018069733 0.8607809

It just took me a few reads through until I got it! Hope this comment helps if other people are struggling like I did. Thanks Bioconductors for writing such a thorough manual and an awesome program!!!

ADD REPLY • link 2.7 years ago Adelyn ▴ 10

score 0 · Answer 1 · 2022-03-03

0

Entering edit mode

Gordon Smyth 51k

@gordon-smyth

Last seen 5 hours ago

WEHI, Melbourne, Australia

goana requires genes to be identified by NCBI Gene IDs whereas you have Ensembl IDs. Please see the goana documentation by typing help("goana").

ADD COMMENT • link 2.7 years ago Gordon Smyth 51k