Try   HackMD

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 a partir do portal do UniProt.

Link direto para a página do UniProt: Proteoma de E. coli no UniProt

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/training/online/course/uniprot-exploring-protein-sequence-and-functional/when-use-uniprot-guided-examples/download-p

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.

Para efeito de treinamento, vamos obter o transcriptoma desta mesma E. coli a partir do genoma, utilizando a ferramenta gffread:

Obtendo o genoma a partir do NCBI:

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:

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:

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
    ​​​​TransDecoder.LongOrfs -t transcriptome.fa
    ​​​​​​  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
    ​​​​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

    Atenção:

    1. O banco de dados para blastp antes deve ser obtido e indexado:
    ​​​​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
    1. Mude o número de threads de acordo com o servidor disponível;
    ​​​​hmmscan --cpu 10 --domtblout ./pfam.domtblout ./Pfam-A.hmm ./transcriptome.fa.transdecoder_dir/longest_orfs.pep

    Obtenha e prepare o banco de dados do Pfam para execução de buscas HMM:

    ​​​​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):
    ​​​​TransDecoder.Predict -t transcriptome.fa \ ​​​​--retain_pfam_hits ./pfam.domtblout \ ​​​​--retain_blastp_hits ./blastp.outfmt6
    ​​​​​​  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 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:
faSplit sequence ./UP000000625_83333.fasta 4 input_file_

Para utilizar o programa faSplit, deve compilar as ferramentas do UCSC-GenomeBrowser.

  • Mapear as sequências utilizando diamond:
for f in ./input_file_*; do emapper.py -m diamond --no_annot --no_file_comments --cpu 16 -i $f -o $f; done
  1. Ortologia e anotação funcional
  • Concatenando todos os arquivos resultantes da etapa anterior
cat *.emapper.seed_orthologs > input_file.emapper.seed_orthologs
  • Busca por ortólogos e anotação
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

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

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:

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:

aggregatoR.R --x="emapper_go.txt" --sep=";" --by.x="ID" --out=./emapper_gos_ok.txt

Remover os cabeçalhos:

grep -P -v '^ID\tGO' ./emapper_gos_ok.txt > ./goatools/association

pode ser feito também com tail:

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:

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:

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

grep '^>' ./UP000000625_83333.fasta | sed 's/ .*$//' | sed 's/^>//' > ./goatools/population

Download go-basic.obo:

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.

É 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:
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) e capturar os IDs relacionados para compor o estudo:
grep 'GO:0003824' emapper_go.txt | cut -f 1 > ./goatools/study

Os primeiros IDs são de termos MUITO abrangentes:
GO:0008150 - biological_process
GO:0003674 - molecular_function
GO:0005575 - cellular_component
GO:0005623 - cell
GO:0044464 - cell part
GO:0009987 - cellular process
GO:0008152 - metabolic process

Enriquecimento funcional

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

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