Transformar anotación GOs de genes a un archivo GAF

Para realizar un análisis de enriquecimiento funcional necesitamos en muchos casos un archivo GAF (por ejemplo en ontologizer). Acá te explico la forma de obtenerlo si cuentas con la anotación GO de tus genes.

1. ¿Qué es un archivo GAF?

El archivo GAF (Gene Association File) es una forma de almacenar la información de anotación funcional de una lista de genes.
Entonces, si tienes una lista de genes (por ejemplo un transcriptoma) y además tienes información de la anotación funcional de cada uno, una forma de almacenar los datos es en el formato GAF el cual consiste en 17 columnas con información definida y una línea por cada anotación funcional para cada gen.
Acá un ejemplo:

DB Object ID

DB Object Symbol

Qualifier

GO ID

DB:Reference

Evidence Code

With

Aspect

DB Object Name

DB Object Synonym

DB Object Type

Taxon

Date

Assigned By

EGGNOG

gen1

GO:0003674

GO_REF:nd

molecular_function

taxon:0

20210712

EGGNOG

gen1

GO:0003824

GO_REF:nd

catalytic activity

taxon:0

20210712

EGGNOG

gen2

GO:0003008

GO_REF:nd

system process

taxon:0

20210712

EGGNOG

gen2

GO:0003674

GO_REF:nd

molecular_function

taxon:0

20210712

EGGNOG

gen2

GO:0005215

GO_REF:nd

transporter activity

taxon:0

20210712

EGGNOG

gen2

GO:0006766

GO_REF:nd

vitamin metabolic process

taxon:0

20210712

EGGNOG

gen3

GO:0003008

GO_REF:nd

system process

taxon:0

20210712

EGGNOG

gen3

GO:0003674

GO_REF:nd

molecular_function

taxon:0

20210712

EGGNOG

gen3

GO:0005215

GO_REF:nd

transporter activity

taxon:0

20210712

El orden de las columnas tiene la siguiente explicación:

Column	Content	Required?	Cardinality	Example
1	DB	required	1	UniProtKB
2	DB Object ID	required	1	P12345
3	DB Object Symbol	required	1	PHO3
4	Qualifier	optional	0 or greater	NOT
5	GO ID	required	1	GO:0003993
6	DB:Reference (\|DB:Reference)	required	1 or greater	PMID:2676709
7	Evidence Code	required	1	IMP
8	With (or) From	optional	0 or greater	GO:0000346
9	Aspect	required	1	F
10	DB Object Name	optional	0 or 1	Toll-like receptor 4
11	DB Object Synonym (\|Synonym)	optional	0 or greater	hToll
12	DB Object Type	required	1	protein
13	Taxon(\|taxon)	required	1 or 2	taxon:9606
14	Date	required	1	20090118
15	Assigned By	required	1	SGD
16	Annotation Extension	optional	0 or greater	part_of(CL:0000576)
17	Gene Product Form ID	optional	0 or 1	UniProtKB:P12345-2

Mas info sobre lo que significa cada columna la encontrarás en este link.

2. El archivo de entrada o input

Para obtener un archivo GAF según las instrucciones que siguen el requisito es contar con un archivo donde los genes de la especie en cuestión han sido anotados (por ejemplo con Eggnog).

Dicho archivo annot.txt lucirá como esto:

query GOs
gen1 GO:0003674,GO:0003824,GO:0004176,GO:0005575,GO:0005622,GO:0005623
gen2 GO:0003008,GO:0003674,GO:0005215,GO:0005575,GO:0005623,GO:0005886,GO:0005887,GO:0006766
gen3 GO:0003008,GO:0003674,GO:0005215,GO:0005575,GO:0005623,GO:0006950,GO:0007275
...

Pero queremos que luzca como lo siguiente:
(ahora, si ya tienes este archivo... puedes saltarte al punto "3".

gen1    GO:0003674
gen1    GO:0003824
gen1    GO:0004176
gen1    GO:0005575
gen1    GO:0005622
gen1    GO:0005623
gen2    GO:0003008
gen2    GO:0003674
gen2    GO:0005215
gen2    GO:0005575
gen2    GO:0005623
gen2    GO:0005886
gen2    GO:0005887
gen2    GO:0006766
gen3    GO:0003008
gen3    GO:0003674
gen3    GO:0005215
gen3    GO:0005575
gen3    GO:0005623
gen3    GO:0006950
gen3    GO:0007275
...

Para ello haremos lo siguiente:

sed -e 's/ /\t/' -e 's/,/\t/g' annot.txt > annot.list.txt
columnas_a_filas.pl annot.list.txt > annot.list2.txt

3. Archivo con datos GO

Necesitamos el nombre de las funciones GO para cada uno de nustros IDs GO.
Para ello descargarémos la base de datos GO y la procesaremos con un script que permite "desglosar" el archivo.

Descargamos el archivo go.obo desde geneontology.org (abajo coloco el comando).

Descargamos el script obo_to_term_functions.py y obo_to_term_tables.py desde mi github (abajo coloco el comando y/o presionar botón derecho del mouse en en "raw" y seleccionar guardar enlace. Estos scripts contienen mis modificaciones a los scripts de sgrote.

wget https://github.com/jomaldon/scripts_bioinfo/raw/master/obo_to_term_tables.py
wget https://github.com/jomaldon/scripts_bioinfo/raw/master/obo_to_term_functions.py
wget http://purl.obolibrary.org/obo/go.obo

Usamos los scripts sobre el archivo go.obo y luego renombramos las categorías principales a código de letras.
Los dos scripts y el archivo go.obo tienen que estar en el mismo directorio.
Los scripts funcionan con Python 2.7 y Python 3.7.

# Obtener "desglose" del archivo go.obo
python obo_to_term_tables.py go.obo .

# Transformamos el nombre de las categorías principales en letras pues es lo que necesita la columna "Aspect" el formato GAF, y guadamos las columnas que nos importan.
awk -F'\t' '{if($3=="biological_process") $3="P"; else if($3=="molecular_function") $3="F"; else if($3=="cellular_component") $3="C"; print $4"\t"$3"\t"$2}' term.txt > term.tab

# Finalmente el archivo term.tab tendrá esta estructura
is_a    relationship    is_a
GO:0000001      P       mitochondrion inheritance
GO:0000002      P       mitochondrial genome maintenance
GO:0000003      P       reproduction
...

4. Limpiamos la lista de GOs

En las anotaciones tipo eggnog se obtienen referencias a GOs que son obsoletos... mi recomendación es eliminarlos antes de crear el archivo GAF.

# Obtengo la lista de GOs de mi archivo
cut -f2  annot.list2.txt | sort | uniq > annot.list2.GO_ids

# Obtengo la lista de GOs del archivo terms.tab
cut -f1 term.tab > term.tab.GO_ids

# Cruzo la lista con la lista de GOs del archivo OBO y me quedo con lo que no hizo match
grep -v -f term.tab.GO_ids annot.list2.GO_ids > annot.list2.txt.borrar

# Lo que no hizo match lo elimino de mi archivo de anotación
grep -v -f annot.list2.txt.borrar annot.list2.txt > annot.list3.txt

5. Obtener el GAF

Cruzamos el archivo annot.list3.txt con el archivo term.tab usando los GO como clave. Además, en el proceso voy generando las columnas que necesita el GAF. Finalmente mprimo todo para obtener el archivo GAF final.

Modificar los siguientes campos según corresponda:
EGGNOG = Base de datos utilziada para la anotación.
20210712 = Fecha de la anotación (es un dato para el registro... no es necesario que sea exacta).
JM = Iniciales o nombre de la persona que realizó la anotación (Assigned By). Es un dato para el registro, puede ser cualquier sigla o nombre).

awk -F'\t' 'BEGIN {OFS = FS} NR==FNR {h[$1] = $2; i[$1] = $3; next} {print "EGGNOG",$1,$1,"",$2,"GO_REF:nd","ND","",h[$2],i[$2],"NA","NA","taxon:0","20210712","JM"}' term.tab annot.list3.txt > annot.list3.GAF

Listo, el archivo annot.list3.GAF tiene el formato adecuado para ontologyzer o el programa que sea que requiera un GAF.

EGGNOG      gen1    gen1            GO:0003674      GO_REF:nd       ND              F       molecular_function      NA      NA      taxon:0 20210712        JM
EGGNOG      gen1    gen1            GO:0003824      GO_REF:nd       ND              F       catalytic activity      NA      NA      taxon:0 20210712        JM
EGGNOG      gen1    gen1            GO:0004176      GO_REF:nd       ND              F       ATP-dependent peptidase activity        NA      NA      taxon:0 20210
712     JM
EGGNON      gen1    gen1            GO:0005575      GO_REF:nd       ND              C       cellular_component      NA      NA      taxon:0 20210712        JM
EGGNOG      gen1    gen1            GO:0005622      GO_REF:nd       ND              C       intracellular anatomical structure      NA      NA      taxon:0 20210
712     JM
...