Automatizando pipelines de bioinformática (parte II): anotação funcional com shell script

Nesta semana darei continuidade à série “automatizando pipelines de bioinformática”. No post anterior falei sobre shell scripts, mostrando como utilizar variáveis e estruturas lógicas (if , for e while) para a criação de programas simples que podem incluir também comandos do terminal. Neste segundo post, mostrarei um exemplo de uso do shell script em uma pipeline simples para realizar a anotação funcional de um genoma microbiano.

Anotação funcional

A anotação funcional consiste na identificação das funções desempenhadas por uma determinada biomolécula, geralmente nucleotídica ou proteica. Ao anotarmos um genoma, atribuímos um nome ao produto de cada gene presente (ex: “DNA Gyrase subunit B”), mas como esta nomenclatura pode variar de acordo o banco de dados que estamos utilizando, e em alguns casos inclusive ser pouco informativa (ex: “protease” ou “hypothetical protein”), é necessário haver formar padronizadas para atribuir informações sobre função molecular para estas moléculas.

Por conta disso, bancos de dados específicos para a organização de ontologias de funções e processos biológicos foram criados, sendo o Gene Ontology (GO) um dos principais. Este banco organiza os diferentes tipos de processos biológicos, funções molecular e localizações subcelulares aos quais uma determinada molécula pode estar associada em uma estrutura de grafo, o que permite diferente níveis de especificidade para estas anotações. A cada nodo no grafo é atribuído um ID único, chamado “termo GO”. No caso de bancos de dados curados manualmente, como o Swissprot, os termos GO relacionados à uma determinada molécula são selecionados a partir da literatura científica, mas é possível usar ferramentas automáticas que realizam a identificação dos termos mais prováveis para uma molécula a partir dos seus hits no BLAST. Esta estratégia é usada por ferramentas como o BLAST2GO, e neste tutorial mostrarei como criar um script que faz isso 😄

Estrutura do projeto

Para este projeto, usaremos a seguinte estrutura de diretórios:

project/
    genome/
        genome.fasta
    prodigal_out/
    emboss_out/
    blast_db/
        uniprot_sp.fasta
    blast_out/
    blast_ids/
    quickgo/
    pipeline.sh

Baixando os arquivos necessários

Os arquivos necessários para iniciarmos este projeto são:

• Genoma da Escherichia coli K-12. Este arquivo deve ficar na pasta genome/ com o nome genome.fasta.
• Arquivo FASTA contendo as sequências do Uniprot-Swissprot. Este arquivo deve ficar na pasta blast_db com o nome uniprot_sp.fasta.

Instalando dependências

Para este exemplo utilizaremos os seguintes programas:

• ncbi-blast+: pacote que implementa o algoritmo BLAST mantido pelo NCBI.
• prodigal: preditor de genes para genomas microbianos.
• emboss: pacote de utilitários para bioinformática.
• jq: utilitário para processamentos de arquivos JSON via linha de comando.

Estes programas podem ser instalados com o seguinte comando:

$ sudo apt install ncbi-blast+ prodigal emboss jq

Criando a pipeline

A nossa pipeline será implementada como um único script de shell com o nome pipeline.sh. Os passos serão executados serialmente, indo do processo de formatação do banco de dados do BLAST e predição das regiões de CDS até a compilação dos resultados em um arquivo final de report.

Formatando o banco de dados do BLAST

O primeiro passo é formatar o banco de dados da ferramenta usando o programa makeblastdb.

• -in: indica o arquivo de entrada.
• -out: indica o arquivo de saída.
• -dbtype: indica o tipo do banco de dados.

# formatar o banco de dados
makeblastdb -in blast_db/uniprot_sp.fasta \
            -out blast_db/uniprot_sp \
            -dbtype prot

Predizendo as CDS do genoma

A predição das regiões codificantes de DNA (coding DNA sequences, CDSs) é realizada com o programa prodigal.

• in: indica o arquivo de entrada.
• -a: indica o arquivo onde devem ser salvas as sequências das proteínas codificadas por cada CDS.

# executar a predição das proteínas codificadas pelo genoma
prodigal -i genome/genome.fasta \
         -a prodigal_out/proteome.fasta

Separando os produtos das CDS

Para os passos posteriormente é necessário que tenhamos uma proteína sendo analisada por vez. Para isso, usamos o utilitário seqretsplit do pacote EMBOSS.

• -sequence: indica o arquivo de entrada.
• -outseq: indica um formato que deve ser usado para o nome dos arquivos de saída.
• -osdirectory2: indica o diretório de saída.

# separar as proteínas em arquivos FASTA
seqretsplit -sequence prodigal_out/proteome.fasta \
            -outseq "*.format" \
            -osdirectory2 emboss_out/

Realizando o BLAST e extraindo os IDs dos melhores hits

Agora, vamos utilizar a estrutura for para repetir o comando blastp para cada uma das sequências de proteínas geradas pelo programa seqretsplit. As opções que utilizamos no BLAST são:

• -query : arquivo de entrada.
• -db : banco de dados que deve ser utilizado.
• -out : arquivo de saída.
• -outfmt : formato de saída (neste caso, 6 indica o formato tabular).
• -evalue : e-value máximo que deve ser aceito.
• -num_threads : número de threads que devem ser utilizada.

Após executarmos o BLAST para cada proteína, realizamos um segundo for para passar por cada resultado e selecionar apenas o Uniprot ID dos 5 primeiros hits de cada proteína. Para fazermos isso, utilizamos uma pipe composta pelos programas cat , head e cut.

Com o cat pegamos o conteúdo do arquivo de saída do BLAST e transferimos para o head, que seleciona apenas as 5 primeiras linhas. Depois, o cut é usado para filtrar apenas os dados da segunda coluna (com a opção -f2), que então é processada por uma segunda instância do cut que usa como delimitador a barra vertical (“|”), e seleciona apenas os valores que estão entre a primeira e segunda barra (-d"|" -f2). Por fim, estes resultados são redirecionados para a pasta blast_ids/.

# rodar o BLAST para cada proteína do proteoma
for protein_fasta in emboss_out/*; do
     blastp -query $protein_fasta \
            -db blast_db/swissprot \
            -out blast_out/$(basename $protein_fasta).tbl
            -outfmt 6 \
            -evalue 1e-100 \
            -num_threads 4
done
# selecionar o Uniprot ID dos 5 melhores hits 
# identificados pelo BLAST 
for blast_out in blast_out/*.tbl; do
    cat $blast_out | head -5 \
        | cut -f2 \
        | cut -d'|' -f2 \
        > blast_ids/$(basename $blast_out).ids
done

Identificando os termos do Gene Ontology

Para buscar os termos GO associados aos hits de cada proteína vamos utilizar a API do QuickGO, uma ferramenta que mapeia dados do Uniprot às anotações com Gene Ontology. Para isso, criamos um for que passa por cada arquivo na pasta blast_ids e lê os ids do Uniprot presentes no arquivo com um while.

Para cada id é usada ferramenta wget para obter os dados da API, que são retornados em formato JSON. Usamos a opção -O - --silent para que o conteúdo seja mostrado na tela, e não passado para um arquivo, e que nenhuma informação adicional (progresso) seja apresentada. Estas informações são passadas para o jq , que usa a query ".results[] | objects | .goId" para selecionar apenas os termos GO do JSON. Usamos o tr apenas para remover as aspas duplas que são retornadas junto com os termos, e então salvamos tudo em um arquivo com extensão .go.

# buscar no QuickGO os termos GO associados com cada ID do Uniprot
# e usar o utilitário "jq" para filtar o JSON retornado
QUICKGO_API=https://www.ebi.ac.uk/QuickGO/services/annotation/
for cds_blastid in blast_ids/*.ids; do
    while read uniprot_id; do
        wget -O - \
             --silent \
             $QUICKGO_API/search?geneProductId=$uniprot_id \
             | jq ".results[] | objects | .goId" \
             | tr -d '"' \
             > quickgo/$(basename $cds_blastid).go
    done < $cds_blastid
done

Integrando os resultados

Agora, que temos os termos GO referentes a cada proteína identificada pelo prodigal, precisamos apenas compilar todos os resultados em uma única tabela. Para isso, usamos um for que passa por cada arquivo na pasta quickgo/ e carrega o nome do gene identificado para a variável gene_id e a lista de termos GO para a variável go_terms, imprimindo na tela estes valores na forma de uma tabela. Note que para selecionar apenas o nome do gene identificado, usamos uma pipe composta pelos comandos echo e cut que é executada dentro de um $(...).

# integrar os dados das anotação funcional em um único arquivo
for cds_quickgoids in quickgo/*.go; do
    gene_id=$(echo $cds_quickgoids | cut -d'.' -f1,2)
    go_terms=$(cat quickgo/$gene_id.*)
    echo -e "gene_id\tgo_terms"
done > report.tbl

Ao final, teremos o seguinte arquivo sendo produzido:

obs: esta é apenas uma parte do arquivo original 😆

O script completo

O script contendo todas as etapas de análise está mostrado abaixo :)

Próximos passos

Neste caso a nossa pipeline de análise está codificada em um único arquivo, e caso precisemos re-executar alguma etapa, toda a pipeline deve ser executada novamente. Uma alternativa seria fragmentar cada e etapa em um script e chamar-los individualmente, e caso queiramos executar tudo com apenas um comando? Para resolver este problema, no próximo tutorial mostrarei como utilizar arquivos Makefile para estruturar as etapas de um projeto e organizar a execução de múltiplos scripts.