7
$\begingroup$

This question is based on a question on BioStars posted >2 years ago by user jack.

It describes a very frequent problem of generating GO annotations for non-model organisms. While it is based on some specific format and single application (Ontologizer), it would be useful to have a general description of the pathway to getting to a GAF file.

Note, that the input format is lacking a bit of essential information, like how it was obtained. Therefore, it is har to assign evidence code. Therefore, lets assume that the assignments of GO terms were done automagically.


I want to do the Gene enrichment using Ontologizer without a predefined association file(it's not model organism).

I have parsed a file with two columns for that organism like this :

geneA  GO:0006950,GO:0005737
geneB  GO:0016020,GO:0005524,GO:0006468,GO:0005737,GO:0004674,GO:0006914,GO:0016021,GO:0015031
geneC  GO:0003779,GO:0006941,GO:0005524,GO:0003774,GO:0005516,GO:0005737,GO:0005863
geneD  GO:0005634,GO:0003677,GO:0030154,GO:0006350,GO:0006355,GO:0007275,GO:0030528

I have downloaded the .ob file from Gene ontology file which contain this information (from here) :

!
! GO IDs (primary only) and name text strings
! GO:0000000 [tab] text string [tab] F|P|C
! where F = molecular function, P = biological process, C = cellular component
!
GO:0000001  mitochondrion inheritance   P
GO:0000002  mitochondrial genome maintenance    P
GO:0000003  reproduction    P
GO:0000005  ribosomal chaperone activity    F
GO:0000006  high affinity zinc uptake transmembrane transporter activity    F
GO:0000007  low-affinity zinc ion transmembrane transporter activity    F
GO:0000008  thioredoxin F
GO:0000009  alpha-1,6-mannosyltransferase activity  F
GO:0000010  trans-hexaprenyltranstransferase activity   F
GO:0000011  vacuole inheritance P

What I need as output is .gaf file in the following format (in the format of the files here):

!gaf-version: 2.0

!Project_name: Leishmania major GeneDB

!URL: http://www.genedb.org/leish

!Contact Email: mb4@sanger.ac.uk

 GeneDB_Lmajor    LmjF.36.4770    LmjF.36.4770        GO:0003723    PMID:22396527    ISO    GeneDB:Tb927.10.10130    F    mitochondrial RNA binding complex 1 subunit, putative    LmjF36.4770    gene    taxon:347515    20120910    GeneDB_Lmajor       
 GeneDB_Lmajor    LmjF.36.4770    LmjF.36.4770        GO:0044429    PMID:20660476    ISS        C    mitochondrial RNA binding complex 1 subunit, putative    LmjF36.4770    gene    taxon:347515    20100803 GeneDB_Lmajor             GeneDB_Lmajor    LmjF.36.4770    LmjF.36.4770        GO:0016554    PMID:22396527    ISO    GeneDB:Tb927.10.10130    P    mitochondrial RNA binding complex 1 subunit, putative    LmjF36.4770    gene   taxon:347515    20120910    GeneDB_Lmajor       
 GeneDB_Lmajor    LmjF.36.4770    LmjF.36.4770        GO:0048255    PMID:22396527    ISO    GeneDB:Tb927.10.10130    P    mitochondrial RNA binding complex 1 subunit, putative    LmjF36.4770    gene    taxon:347515    20120910    GeneDB_Lmajor  

How to create your own GO association file (gaf)?

$\endgroup$
  • $\begingroup$ Please add some tags for GO, non-model organism, gaf, overrepresentation analysis etc. I just used a random tag because there are not many sensible tags in the list yet. $\endgroup$ – Michael Jun 8 '17 at 7:59
  • $\begingroup$ I modified a bit the question, added tags and I centered the question in creating the file in that format, if you disagree, just revert my changes $\endgroup$ – llrs Jun 8 '17 at 8:07
  • 2
    $\begingroup$ I added gene-ontology and gaf. I'm not sure whether a over-representation-analysis or a enrichment-analysis would be more appropriate. $\endgroup$ – Devon Ryan Jun 8 '17 at 8:08
  • $\begingroup$ @Llopis: Hopefully I didn't just accidentally overwrite your changes! $\endgroup$ – Devon Ryan Jun 8 '17 at 8:09
  • $\begingroup$ No luck, my edit was rejected (or at least doesn't appear) :\ $\endgroup$ – llrs Jun 8 '17 at 8:11
8
$\begingroup$

Here's a Perl script that can do this:

#!/usr/bin/env perl 
use strict;
use warnings;

## Change this to whatever taxon you are working with
my $taxon = 'taxon:1000';
chomp(my $date = `date +%Y%M%d`);

my (%aspect, %gos);
## Read the GO.terms_and_ids file to get the aspect (sub ontology)
## of each GO term. 
open(my $fh, $ARGV[0]) or die "Need a GO.terms_and_ids file as 1st arg: $!\n";
while (<$fh>) {
    next if /^!/;
    chomp;
    my @fields = split(/\t/);
    ## $aspect{GO:0000001} = 'P'
    $aspect{$fields[0]} = $fields[2];
}
close($fh);

## Read the list of gene annotations
open($fh, $ARGV[1]) or die "Need a list of gene annotattions as 2nd arg: $!\n";
while (<$fh>) {
    chomp;
    my ($gene, @terms) = split(/[\s,]+/);
    ## $gos{geneA} = (go1, go2 ... goN)
    $gos{$gene} = [ @terms ];
}
close($fh);

foreach my $gene (keys(%gos)) {
    foreach my $term (@{$gos{$gene}}) {
        ## Warn and skip if there is no aspect for this term
        if (!$aspect{$term}) {
            print STDERR "Unknown GO term ($term) for gene $gene\n";
            next;
        }
        ## Build a pseudo GAF line 
        my @out = ('DB', $gene, $gene, ' ', $term, 'PMID:foo', 'TAS', ' ', $aspect{$term},
                             $gene, ' ', 'protein', $taxon, $date, 'DB', ' ', ' ');
        print join("\t", @out). "\n";
    }
}

Make it executable and run it with the GO.terms_and_ids file as the 1st argument and the list of gene annotations as the second. Using the current GO.terms_and_ids and the example annotations in the question, I get:

$ foo.pl GO.terms_and_ids file.gos 
DB  geneD   geneD       GO:0005634  PMID:foo    TAS     C   geneD       protein taxon:1000  20170308    DB       
DB  geneD   geneD       GO:0003677  PMID:foo    TAS     F   geneD       protein taxon:1000  20170308    DB       
DB  geneD   geneD       GO:0030154  PMID:foo    TAS     P   geneD       protein taxon:1000  20170308    DB       
Unknown GO term (GO:0006350) for gene geneD
DB  geneD   geneD       GO:0006355  PMID:foo    TAS     P   geneD       protein taxon:1000  20170308    DB       
DB  geneD   geneD       GO:0007275  PMID:foo    TAS     P   geneD       protein taxon:1000  20170308    DB       
DB  geneD   geneD       GO:0030528  PMID:foo    TAS     F   geneD       protein taxon:1000  20170308    DB       
DB  geneB   geneB       GO:0016020  PMID:foo    TAS     C   geneB       protein taxon:1000  20170308    DB       
DB  geneB   geneB       GO:0005524  PMID:foo    TAS     F   geneB       protein taxon:1000  20170308    DB       
DB  geneB   geneB       GO:0006468  PMID:foo    TAS     P   geneB       protein taxon:1000  20170308    DB       
DB  geneB   geneB       GO:0005737  PMID:foo    TAS     C   geneB       protein taxon:1000  20170308    DB       
DB  geneB   geneB       GO:0004674  PMID:foo    TAS     F   geneB       protein taxon:1000  20170308    DB       
DB  geneB   geneB       GO:0006914  PMID:foo    TAS     P   geneB       protein taxon:1000  20170308    DB       
DB  geneB   geneB       GO:0016021  PMID:foo    TAS     C   geneB       protein taxon:1000  20170308    DB       
DB  geneB   geneB       GO:0015031  PMID:foo    TAS     P   geneB       protein taxon:1000  20170308    DB       
DB  geneA   geneA       GO:0006950  PMID:foo    TAS     P   geneA       protein taxon:1000  20170308    DB       
DB  geneA   geneA       GO:0005737  PMID:foo    TAS     C   geneA       protein taxon:1000  20170308    DB       
DB  geneC   geneC       GO:0003779  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB       
DB  geneC   geneC       GO:0006941  PMID:foo    TAS     P   geneC       protein taxon:1000  20170308    DB       
DB  geneC   geneC       GO:0005524  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB       
DB  geneC   geneC       GO:0003774  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB       
DB  geneC   geneC       GO:0005516  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB       
DB  geneC   geneC       GO:0005737  PMID:foo    TAS     C   geneC       protein taxon:1000  20170308    DB       
DB  geneC   geneC       GO:0005863  PMID:foo    TAS     C   geneC       protein taxon:1000  20170308    DB       

Note that this is very much a pseudo-GAF file since most of the fields apart from the gene name, GO term and sub-ontology are fake. It should still work for what you need, however.

$\endgroup$
  • $\begingroup$ I think the evidence code should be changed to IEA, because it is unlikely one would have manually annotated GO's for a non-model organism. $\endgroup$ – Michael Jun 8 '17 at 9:47
  • 1
    $\begingroup$ @Michael, it doesn't really make any difference so I chose TAS since these are annotations provided by the OP. But sure, IEA would also work. The field will be ignored in any case when doing an enrichment analysis so you can put whatever you prefer. $\endgroup$ – terdon Jun 8 '17 at 10:09
  • $\begingroup$ Sure, IEA might be even discarded by some tools. $\endgroup$ – Michael Jun 8 '17 at 10:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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