# Anotação funcional com eggNOG-mapper e goatools
## Obtenção das proteínas
O exemplo de análise será realizado utilizando o proteoma da bactéria [*Escherichia coli* strain K12](https://www.uniprot.org/taxonomy/83333) a partir do portal do [UniProt](https://www.uniprot.org).
Link direto para a página do UniProt: [Proteoma de *E. coli* no UniProt](https://www.uniprot.org/proteomes/UP000000625)
Link FTP para download direto das sequências: ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/reference_proteomes/Bacteria/UP000000625_83333.fasta.gz
Referência com mais detalhes sobre o treinamento para obtenção do proteoma a partir do [EBI](https://www.ebi.ac.uk): https://www.ebi.ac.uk/training/online/course/uniprot-exploring-protein-sequence-and-functional/when-use-uniprot-guided-examples/download-p
```bash=
wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/reference_proteomes/Bacteria/UP000000625_83333.fasta.gz
```
### Predição de um proteoma
Se estiver em mãos um transcriptoma recém montado, ainda sem a predição de proteínas e deseja realizar a predição de proteínas, o progrma indicado é o [TransDecoder](https://transdecoder.github.io/).
Para efeito de treinamento, vamos obter o transcriptoma desta mesma *E. coli* a partir do genoma, utilizando a ferramenta [gffread](https://github.com/gpertea/gffread):
Obtendo o genoma a partir do NCBI:
```bash=
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz
gunzip GCF_000005845.2_ASM584v2_genomic.fna.gz
```
Obtendo a anotação do genoma a partir do NCBI:
```bash=
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.gff.gz
gunzip GCF_000005845.2_ASM584v2_genomic.gff.gz
```
Obtendo as sequências do transcriptoma a partir do genoma:
```bash=
gffread GCF_000005845.2_ASM584v2_genomic.gff -g GCF_000005845.2_ASM584v2_genomic.fna -w transcriptome.fa
```
#### Predição
- Extrair as sequências de *Open Reading Frames* (ORFs) mais longas
```bash=
TransDecoder.LongOrfs -t transcriptome.fa
```
::: info
Observe bem os parâmetros que em determinadas situações devem ser alterados!
Em especial:
--genetic_code string (genetic code)
--gene_trans_map string (gene-to-transcript identifier mapping file)
-m int (minimum protein length)
-S strand-specific (only analyzes top strand)
:::
- Obter informações sobre similaridades, a fim de inferir homologia
```bash=
blastp -query transcriptome.fa.transdecoder_dir/longest_orfs.pep \
-db ./uniprot_sprot -max_target_seqs 1 \
-outfmt 6 -evalue 1e-5 -num_threads 10 > blastp.outfmt6
```
::: info
Atenção:
1) O banco de dados para blastp antes deve ser obtido e indexado:
```bash=
wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz
gunzip uniprot_sprot.fasta.gz
makeblastdb -in ./uniprot_sprot.fasta -dbtype prot -out ./uniprot_sprot
```
2) Mude o número de *threads* de acordo com o servidor disponível;
:::
```bash=
hmmscan --cpu 10 --domtblout ./pfam.domtblout ./Pfam-A.hmm ./transcriptome.fa.transdecoder_dir/longest_orfs.pep
```
::: info
Obtenha e prepare o banco de dados do Pfam para execução de buscas HMM:
```bash=
wget ftp://ftp.ebi.ac.uk/pub/databases/Pfam/releases/Pfam31.0/Pfam-A.hmm.gz
gunzip Pfam-A.hmm.gz
hmmpress Pfam-A.hmm
```
:::
- Predição de proteínas
Obtenção de uma única predição para cada transcrito (--single_best_only):
```bash=
TransDecoder.Predict -t transcriptome.fa \
--retain_pfam_hits ./pfam.domtblout \
--retain_blastp_hits ./blastp.outfmt6
```
::: info
Observe bem os parâmetros que em determinadas situações devem ser alterados!
Em especial:
--genetic_code string (genetic code)
--single_best_only (Retain only the single best orf per transcript, prioritized by homology then orf length)
:::
## Predição de proteínas ortólogas
### eggNOG-mapper
O programa [eggNOG-mapper](https://github.com/eggnogdb/eggnog-mapper) realiza uma anotação rápida utilizando como referência o banco de dados eggNOG.
1. Buscas por Homologia
- Dividir as sequências em 4 arquivos:
```bash=
faSplit sequence ./UP000000625_83333.fasta 4 input_file_
```
:::info
Para utilizar o programa faSplit, deve compilar as ferramentas do [UCSC-GenomeBrowser](https://github.com/ucscGenomeBrowser/kent/releases).
:::
- Mapear as sequências utilizando diamond:
```bash=
for f in ./input_file_*; do
emapper.py -m diamond --no_annot --no_file_comments --cpu 16 -i $f -o $f;
done
```
2. Ortologia e anotação funcional
- Concatenando todos os arquivos resultantes da etapa anterior
```bash=
cat *.emapper.seed_orthologs > input_file.emapper.seed_orthologs
```
- Busca por ortólogos e anotação
```bash=
emapper.py --annotate_hits_table input_file.emapper.seed_orthologs --no_file_comments -o output_file --cpu 16
```
## Enriquecimento Funcional de termos do GO
### Criação do diretório de saída
```bash=
mkdir -p goatools
```
### Preparo do arquivo das associações
O primeiro arquivo será o goatools/**association**, que irá conter todas as relações entre os IDs possíveis de estarem no conjunto que será avaliado o enriquecimento dos termos. Neste caso serão todas as proteínas (proteoma) de *E. coli*.
O preparo do arquivo goatools/**association** começa a partir do resultado de anotação com eggNOG-mapper vamos selecionar a primeira (ID) e a sétima coluna (GOs):
```bash=
cut -f 1,7 output_file.emapper.annotations > emapper_gos.txt
```
Devemos formatar o arquivo resultante (emapper_gos.txt) para conter um único ID associado a um único GO por linha:
```bash=
splitteR.R x="./emapper_gos.txt" --by.x="," col.x="GO" colnames.x="ID,GO" --noh.x --out=./emapper_go.txt
```
Então, agregar novamente por ID, porém, separando com ";" os GOs:
```bash=
aggregatoR.R --x="emapper_go.txt" --sep=";" --by.x="ID" --out=./emapper_gos_ok.txt
```
Remover os cabeçalhos:
```bash=
grep -P -v '^ID\tGO' ./emapper_gos_ok.txt > ./goatools/association
```
pode ser feito também com *tail*:
```bash=
tail -n +2 emapper_gos_ok.txt > ./goatools/association
```
Abaixo o resultado das primeiras linhas do arquivo ./goatools/**association**:
### Montando arquivo de entrada para o WEGO
O arquivo a seguir contém todos os IDs de genes de *E. coli* com anotação GO no formato de entrada nativo do [WEGO](http://wego.genomics.org.cn/):
```bash=
tail -n +2 emapper_gos_ok.txt | perl -F"\t" -lane 'my @GO=split(/;/, $F[1]); print join("\t", $F[0], @GO);' > WEGO.txt
```
Para utilizá-lo, basta acessar o site http://wego.genomics.org.cn/.
..., selecionar o tipo de formato ("Native Format") do arquivo (WEGO.txt), carregá-lo e submeter.
Para obter o clássico gráfico de barras com a proporção de termos do GO, clique no botão **Graph**:
::: info
Para adicionar mais de uma lista a fim de comparar a proporção de termos do GO, é necessário selecionar mais de um arquivo contendo a lista desejada, antes de submeter.
:::
### Preparo do arquivo da população
O segundo arquivo será o goatools/**population**, que irá conter todos os IDs possíveis de estarem no conjunto que será avaliado o enriquecimento dos termos. Neste caso serão todas as proteínas (proteoma) de *E. coli*, ou seja, de toda a população de proteínas que possivelmente podem estar no arquivo de estudo (goatools/**study**).
```bash=
grep '^>' ./UP000000625_83333.fasta | sed 's/ .*$//' | sed 's/^>//' > ./goatools/population
```
Download **go-basic.obo**:
```bash=
wget --timestamping http://geneontology.org/ontology/go-basic.obo 2> /dev/null
```
### Preparo do arquivo do estudo
O arquivo goatools/**study** deve conter os IDs das proteínas em que se deseja investigar se há enriquecimento de termos do GO.
:::warning
É importante que todo ID presente no arquivo de estudo (goatools/**study**) também esteja presente no arquivo com a população (goatools/**population**).
:::
Abaixo, vamos apenas fazer um teste, considerando alguns IDs de proteínas que estão associadas com alguns GOs da lista de GOs mais frequentes. Essa prática serve apenas para esse efeito de teste.
- Verificar os GOs mais frequentes:
```bash=
grep -v '^ID' emapper_go.txt | cut -f 2 | nsort | uniq -c | sed 's/^ *//' | sed 's/ /\t/' | nsort -t$'\t' -k1nr,1nr | more
```
- Utilizar um dos IDs (Ex: **GO:0003824** - [catalytic activity](http://amigo.geneontology.org/amigo/term/GO:0003824)) e capturar os IDs relacionados para compor o estudo:
```bash=
grep 'GO:0003824' emapper_go.txt | cut -f 1 > ./goatools/study
```
:::info
Os primeiros IDs são de termos MUITO abrangentes:
GO:0008150 - [biological_process](http://amigo.geneontology.org/amigo/term/GO:0008150)
GO:0003674 - [molecular_function](http://amigo.geneontology.org/amigo/term/GO:0003674)
GO:0005575 - [cellular_component](http://amigo.geneontology.org/amigo/term/GO:0005575)
GO:0005623 - [cell](http://amigo.geneontology.org/amigo/term/GO:0005623)
GO:0044464 - [cell part](http://amigo.geneontology.org/amigo/term/GO:0044464)
GO:0009987 - [cellular process](http://amigo.geneontology.org/amigo/term/GO:0009987)
GO:0008152 - [metabolic process](http://amigo.geneontology.org/amigo/term/GO:0008152)
:::
### Enriquecimento funcional
```bash=
find_enrichment.py ./goatools/study \
./goatools/population \
./goatools/association \
--outfile=./goatools/goea_all.xlsx,./goatools/goea_all.tsv \
--pvalcalc=fisher_scipy_stats \
--obo=./go-basic.obo \
--pval=1 \
--alpha=0.05
```
O resultado estará nos arquivos ./goatools/**goea_all.xlsx** e ./goatools/**./goatools/goea_all.tsv**:
Observando o resultado final, certamente terá como um dos itens mais enriquecidos o termo **catalytic activity** (GO:0003824):
```bash=
cut -f 1,2,4,5,6,7,8,9,10,11,12,13 ./goatools/goea_all.tsv | nsort -t$'\t' -k6n,6n | head -5
```
Sign in with Wallet
Connect another wallet