<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>

CC BY-SA 3.0
2018-03-29T09:22:05.343Z
2018-03-29T09:22:05.343Z
benn
<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>

CC BY-SA 3.0
2018-03-28T09:33:32.303Z
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.
llrs
CC BY-SA 3.0
2018-03-28T09:50:23.117Z
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).
benn
CC BY-SA 3.0
2018-03-28T09:59:51.833Z
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.
llrs
CC BY-SA 3.0
2018-03-28T10:10:53.613Z
@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.
benn
CC BY-SA 3.0
2018-03-28T11:06:34.823Z
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.
llrs
CC BY-SA 3.0
2018-03-29T10:14:32.250Z
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.
Ian Sudbery
CC BY-SA 3.0
2018-03-29T14:30:12.520Z
@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.
benn
CC BY-SA 3.0
2018-03-30T10:01:54.913Z
@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
llrs
CC BY-SA 4.0
2018-05-22T13:56:11.337Z
@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...
llrs
CC BY-SA 4.0
2018-05-22T14:01:50.693Z
@Llopis, thanks for getting back with the reference.
benn
<rna-seq><go-enrichment><go>
GOseq analysis with evidence code filter
# 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")


---

