Answers.Body: <p>A bit weird to answer my own question, but for who's interested I can show the 'hard' way using <code>biomaRt</code> and then use <code>goseq</code> with <code>gene2cat</code> defined.</p>
  
  <p>So first you need biomaRt to get annotation of the genes. I have ensembl genes from mm10.</p>
  
  <pre><code>library(biomaRt)
  
  table &lt;- read.table(file, header = T, stringsAsFactors = F)
  
  ensembl &lt;- useMart("ENSEMBL_MART_ENSEMBL")
  
  ensembl &lt;- useDataset("mmusculus_gene_ensembl", mart = ensembl)
  
  allGenes &lt;- table$Geneid
  
  GOs &lt;- 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 &lt;- GOs[GOs$go_linkage_type != "IEA" &amp; GOs$namespace_1003 == "biological_process",]
  
  dim(GOs_no_IEA)
  [1] 111944      5
  </code></pre>
  
  <p>Then you can continue with goseq.</p>
  
  <pre><code>library(goseq)
  
  All &lt;- read.table(input_All, sep = "\t", header = T)
  
  All_genes &lt;- as.character(All$Geneid)
  
  TopTable &lt;- read.table(input, sep = "\t", header = T)
  
  sig_genes &lt;- as.character(TopTable$Ensembl)
  
  # Make gene vector
  
  gene.vector &lt;- as.integer(All_genes %in% sig_genes)
  
  names(gene.vector) &lt;- All_genes
  
  gene.lengths &lt;- All$Length
  
  pwf &lt;- nullp(gene.vector, "mm10", "ensGene", bias.data = gene.lengths)
  head(pwf)
  
  GO.wall &lt;- goseq(pwf, "mm10", "ensGene", gene2cat = GOs_no_IEA[, c(1:2)], use_genes_without_cat = F, test.cats = "GO:BP")
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-03-29T09:22:05.343Z
Answers.LastActivityDate: 2018-03-29T09:22:05.343Z
Answers.UserDisplayName: benn
Body: <p>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 <a href="http://www.geneontology.org/page/guide-go-evidence-codes" rel="nofollow noreferrer">here</a>.</p>
  
  <p>Is it possible to (easily) do GOseq analysis without this evidence group?</p>
  
  <p>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? </p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-03-28T09:33:32.303Z
Comments.Text: 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.
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-03-28T09:50:23.117Z
Comments.Text: 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).
Comments.UserDisplayName: benn
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-03-28T09:59:51.833Z
Comments.Text: 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.
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-03-28T10:10:53.613Z
Comments.Text: @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.
Comments.UserDisplayName: benn
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-03-28T11:06:34.823Z
Comments.Text: Yes, I mean the `elim` and the `weight01` methods. I agree that they won't eliminate false positive genes but it will let you focus on genes  on categories that are really relevant. Anyway, let's see if someone came up with a solution, even if it is the long one.
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-03-29T10:14:32.250Z
Comments.Text: 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.
Comments.UserDisplayName: Ian Sudbery
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-03-29T14:30:12.520Z
Comments.Text: @Llopis, can you direct me to the papers about the quality of IEA and about not using evidence codes for filtering. I am interested in their motivation, thanks.
Comments.UserDisplayName: benn
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-03-30T10:01:54.913Z
Comments.Text: @b.nota sorry I though I had the references but I don't. I am looking for those articles, if I find them I'll post in a new question and let you know
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-05-22T13:56:11.337Z
Comments.Text: @b.nota I found the reference! I found that [Gene Ontology semantic similarity tools: survey on features and challenges for biological knowledge discovery](https://academic.oup.com/bib/article-abstract/18/5/886/2562801?redirectedFrom=fulltext) recommends (under the section "Influence of GO evidence codes") using all the data and evidence codes, because the tools (that make these associations IEA) are getting better to assign a term (and it provides references to semantic similarities with IEA and without IEA). Edit: Now perhaps I should write a question on itself...
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-05-22T14:01:50.693Z
Comments.Text: @Llopis, thanks for getting back with the reference.
Comments.UserDisplayName: benn
Tags: <rna-seq><go-enrichment><go>
Title: GOseq analysis with evidence code filter
ThreadText: # GOseq analysis with evidence code filter
  
  (Asked by: username_0 on 2018-03-28)
  
  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](http://www.geneontology.org/page/guide-go-evidence-codes).
  
  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?
  
  **Tags:** rna-seq, go-enrichment, go
  
  > **Comments:**
  > **username_1**: 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.  
  > **username_0**: 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).  
  > **username_1**: 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.  
  > **username_0**: @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.  
  > **username_1**: Yes, I mean the `elim` and the `weight01` methods. I agree that they won't eliminate false positive genes but it will let you focus on genes  on categories that are really relevant. Anyway, let's see if someone came up with a solution, even if it is the long one.  
  > **username_2**: 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.  
  > **username_0**: @Llopis, can you direct me to the papers about the quality of IEA and about not using evidence codes for filtering. I am interested in their motivation, thanks.  
  > **username_1**: @b.nota sorry I though I had the references but I don't. I am looking for those articles, if I find them I'll post in a new question and let you know  
  > **username_1**: @b.nota I found the reference! I found that [Gene Ontology semantic similarity tools: survey on features and challenges for biological knowledge discovery](https://academic.oup.com/bib/article-abstract/18/5/886/2562801?redirectedFrom=fulltext) recommends (under the section "Influence of GO evidence codes") using all the data and evidence codes, because the tools (that make these associations IEA) are getting better to assign a term (and it provides references to semantic similarities with IEA and without IEA). Edit: Now perhaps I should write a question on itself...  
  > **username_0**: @Llopis, thanks for getting back with the reference.  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_0 on 2018-03-29. Score: 4
  
  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")
  
  ```
  
  
  ---
  
  