4
$\begingroup$

When using GOseq analysis on RNA-seq data, I often find many 'false positives'. What I mean with that is that some genes in certain categories are not really involved in the process, but are only involved due to 'Electronic Annotation' (IEA). Most of the genes are IEA. A list of all evidence codes is here.

Is it possible to (easily) do GOseq analysis without this evidence group?

I mean, I know how to do it the hard way (make a custom GO-genes annotation in e.g., biomaRt removing all IEA evidence code genes, and perform GOseq with custom genome). But is there also a easy way in GOseq available?

Improve this question
$\endgroup$
10
  • $\begingroup$ For what it is worth, there have been papers that say that the quality of the IEA is not so bad (as bad as other categories), also the GO says that the evidence codes are not designed to be used as quality filters. $\endgroup$
    – llrs
    Mar 28, 2018 at 9:33
  • $\begingroup$ Thanks @Llopis, I mostly agree with their statement about IEA quality. However, I was thinking of ways to prevent getting enriched GO terms such as 'spermatogenesis' in female mice experiments... Another motivation is to get 'real' target genes for further investigation (often you get false positives when further investigating). $\endgroup$
    – benn
    Mar 28, 2018 at 9:50
  • $\begingroup$ Then I recommend to give a look at the topGO package, it should filter most of them out, via taking into account the DAG structure of GO. $\endgroup$
    – llrs
    Mar 28, 2018 at 9:59
  • $\begingroup$ @Llopis, do you mean with elim function? That reduces the number of significant GO terms, but won't eliminate false positive genes (targets) for further research. $\endgroup$
    – benn
    Mar 28, 2018 at 10:10
  • 1
    $\begingroup$ I think the answer to the quesiton is that there is not easy way to do it. But as a comment i'd say that I would never trust a GO analysis to provide "answers", only suggestions for further analysis, and as such any genes in significant categories are always going to have to be examined by a domain expert. $\endgroup$
    – Ian Sudbery
    Mar 29, 2018 at 10:14

1 Answer 1

Reset to default
4
$\begingroup$

A bit weird to answer my own question, but for who's interested I can show the 'hard' way using biomaRt and then use goseq with gene2cat defined.

So first you need biomaRt to get annotation of the genes. I have ensembl genes from mm10.

library(biomaRt)

table <- read.table(file, header = T, stringsAsFactors = F)

ensembl <- useMart("ENSEMBL_MART_ENSEMBL")

ensembl <- useDataset("mmusculus_gene_ensembl", mart = ensembl)

allGenes <- table$Geneid

GOs <- getBM(attributes=c('ensembl_gene_id', 'go_id', 'name_1006', 'go_linkage_type', 'namespace_1003'), 
      filters = 'ensembl_gene_id', 
      values = allGenes, 
      mart = ensembl)

dim(GOs)
[1] 477082      5

table(GOs$go_linkage_type)
          EXP    IBA     IC    IDA    IEA    IEP    IGI    IMP    IPI    ISA    ISM    ISO    ISS    NAS     ND    RCA    TAS 
 36931     36  27918    525  34038 210114    630   6891  30219   6205   4234     20  77365  24673    609  12821     99   3754

# I only want BP and no IEA evidence code
GOs_no_IEA <- GOs[GOs$go_linkage_type != "IEA" & GOs$namespace_1003 == "biological_process",]

dim(GOs_no_IEA)
[1] 111944      5

Then you can continue with goseq.

library(goseq)

All <- read.table(input_All, sep = "\t", header = T)

All_genes <- as.character(All$Geneid)

TopTable <- read.table(input, sep = "\t", header = T)

sig_genes <- as.character(TopTable$Ensembl)

# Make gene vector

gene.vector <- as.integer(All_genes %in% sig_genes)

names(gene.vector) <- All_genes

gene.lengths <- All$Length

pwf <- nullp(gene.vector, "mm10", "ensGene", bias.data = gene.lengths)
head(pwf)

GO.wall <- goseq(pwf, "mm10", "ensGene", gene2cat = GOs_no_IEA[, c(1:2)], use_genes_without_cat = F, test.cats = "GO:BP")
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.