Gene Onotology: Can't get the goana function to work
1
0
Entering edit mode
Adelyn • 0
@47ceacfa
Last seen 3 months 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 • 172 views
ADD COMMENT
0
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
0
Entering edit mode
@gordon-smyth
Last seen 1 hour 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

Login before adding your answer.

Traffic: 258 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6