Hi,

I have used GOexpress to analyze my DESeq2 results following the instructions from this GOexpress: use DESeq2 data object as input to get my DESeq data into an ExpressionSet. The custom annotations that I've used are shown below. The GOanalyze step seems to work fine, but when I use the pValue_GO function it gives me an error message. Any ideas why this is occurring?

Thanks,

Brian

> head(GOgenes.Ensembl)
           gene_id      go_id
1 CADAFLAG00000001 GO:0004601
2 CADAFLAG00000001 GO:0020037
3 CADAFLAG00000001 GO:0055114
4 CADAFLAG00000002 GO:0007155
5 CADAFLAG00000002 GO:0005515
6 CADAFLAG00000003 GO:0008483

> head(allgenes.Ensembl)
           gene_id external_gene_id                                                                                description
1 CADAFLAG00000001      AFLA_071780                   Dyp-type peroxidase family protein  [Source:UniProtKB/TrEMBL;Acc:B8MVM1]
2 CADAFLAG00000002      AFLA_071790                          Galactose oxidase, putative  [Source:UniProtKB/TrEMBL;Acc:B8MVM2]
3 CADAFLAG00000003      AFLA_071800 Branched-chain amino acid aminotransferase, putative  [Source:UniProtKB/TrEMBL;Acc:B8MVM3]
4 CADAFLAG00000004      AFLA_071910                             MFS tranporter, putative  [Source:UniProtKB/TrEMBL;Acc:B8MVM4]
5 CADAFLAG00000005      AFLA_071920                 Class III aminotransferase, putative  [Source:UniProtKB/TrEMBL;Acc:B8MVM5]
6 CADAFLAG00000006      AFLA_072030                     Putative uncharacterized protein  [Source:UniProtKB/TrEMBL;Acc:B8MVM6]

> head(goresults_aw$genes)
                  Score Rank external_gene_name
CADAFLAG00000082 0.0040    1        AFLA_072990
CADAFLAG00001182 0.0040    1        AFLA_084020
CADAFLAG00002622 0.0040    1        AFLA_032470
CADAFLAG00002994 0.0040    1        AFLA_036190
CADAFLAG00000583 0.0035    5        AFLA_078010
CADAFLAG00001308 0.0035    5        AFLA_085290
                                                                                            description
CADAFLAG00000082                 Putative uncharacterized protein  [Source:UniProtKB/TrEMBL;Acc:B8MVV2]
CADAFLAG00001182 Homeobox and C2H2 transcription factor, putative  [Source:UniProtKB/TrEMBL;Acc:B8MZB3]
CADAFLAG00002622           Phosphoglycolate phosphatase, putative  [Source:UniProtKB/TrEMBL;Acc:B8N1U3]
CADAFLAG00002994              Transcription factor Rap1, putative  [Source:UniProtKB/TrEMBL;Acc:B8N2D4]
CADAFLAG00000583                      Sugar transporter, putative  [Source:UniProtKB/TrEMBL;Acc:B8MXL4]
CADAFLAG00001308                 Putative uncharacterized protein  [Source:UniProtKB/TrEMBL;Acc:B8MZN9]

> head(minimalSet)
ExpressionSet (storageMode: lockedEnvironment)
assayData: 1 features, 4 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: X30_91 X30_99 X37_91 X37_99
  varLabels: aw temp
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation:  

> goresults_aw=GO_analyse(eSet = minimalSet, f = 'aw', GO_genes = GOgenes.Ensembl, all_GO=allGO.Ensembl, all_genes = allgenes.Ensembl)
Using custom GO_genes mapping ...
8489 features from ExpressionSet found in the mapping table.
Using custom GO terms description ...
Analysis using method randomForest on factor aw for 13485 genes. This may take a few minutes ...
ntree      OOB      1      2
  100:  50.00% 50.00% 50.00%
  200:  50.00% 50.00% 50.00%
  300:  50.00% 50.00% 50.00%
  400:  50.00% 50.00% 50.00%
  500:  50.00% 50.00% 50.00%
  600:  50.00% 50.00% 50.00%
  700:  50.00% 50.00% 50.00%
  800:  50.00% 50.00% 50.00%
  900:  50.00% 50.00% 50.00%
1000:  50.00% 50.00% 50.00%
Using custom gene descriptions ...
Merging score into result table ...
> goresults_aw.pval = pValue_GO(result=goresults_aw)
[==================================================] 1000 of 1000
Error in `[.data.frame`(GO.new, , c("go_id", "ave_rank", "ave_score",  :
  undefined columns selected

Actually, i'm pretty sure I found the typo.

Line 350 of https://github.com/kevinrue/GOexpress/blob/master/R/analysis.R

should be ] <- "namespace_1003" instead of ] <- "namespace_1006"

I suspect your "allGO" annotations have a column "namespace". A quick fix while I correct the code, would be for you to rename the column (in your custom GO annotations) to "namespace_1003".

This is simply because I allowed "namespace" as a shortcut, but to remain consistent with the Ensembl Biomart, I rename it namespace_1003 behind the scenes.

I will introduce this bug fix in both release and devel version in the next hour. A new release 1.4.1 in release and 1.5.2 in devel)  should be available from the Bioconductor within the next 24 hours (need to pass the nightly built to be published).

Apologies for the delay in your analyses, and I hope I guessed correctly the source of your error. Otherwise, please do post here the head of your GO annotations, and the head of goresults_aw$GO.

Sincerely,

Kevin

Hi Brian,

This is strange, if Go_analyse ran correctly.

May I have a peek at your head(goresults_aw$GO) ?

Thank you!

Kevin