Introdução ao BioPython (parte II): utilizando o BLAST com o Python

Olá pessoal! Hoje daremos continuidade à nossa série de posts sobre BioPython. No nosso primeiro post falamos sobre os a manipulação se sequências, dados de anotação, e como podemos realizar o parsing e criação de novos records com uso dos módulos Bio.Seq, Bio.SeqRecord, Bio.SeqFeature e Bio.SeqIO. Agora abordaremos o módulo Bio.Blast, que fornece uma interface para uso da ferramenta BLAST pelo Python.

Como vimos na nossa série de posts sobre o BLAST, esta é uma ferramenta que permite a realização de busca por similaridade em bancos de dados de sequências biológicas, retornando aquelas (hits) que possuem maior score e significação (e-value). É uma ferramenta amplamente utilizada em bioinformática e que pode ser utilizada tanto através de uma interface web, como a disponibilizada pelo NCBI, quando através da linha de comando.

Caso tenha dúvidas quanto a instalação do pacote BioPython, só conferir os exemplos no primeiro post.

Obtendo os dados

Para este tutorial vamos utilizar dois arquivos FASTA: um contendo um conjunto de sequências de proteínas (dataset) e outro contendo apenas uma sequência (query), que podem ser baixados nos links abaixo.

• dataset.fasta
• query.fasta

Instalando o BLAST

Para este tutorial é importante que tenhamos o BLAST instalado. Para isso, podemos utilizar o seguinte comando no Linux:

$ sudo apt install ncbi-blast+

Criando um banco de dados para o BLAST

Para que possamos utilizar o BLAST é necessário antes de mais nada criar um banco de dados a partir do nosso arquivos de sequência. Este processo é chamado “formatação”, e é realizado por uma ferramenta chamada makeblastdb. Para que possamos realizar essa formatação com o BioPython é necessário utilizar a classe NcbimakeblastdbCommandline disponível no módulo Bio.Blast.Applications.

Neste exemplo, estamos passando para o construtor da classe 3 parâmetros:

• input_file: o arquivos de sequências que queremos formatar;
• dbtype: o tipo de banco de dados (prot é usado pois é um banco de proteínas, caso fosse de nucleotídeos seria nucl);
• out: prefixo do banco de dados a ser criado.

Ao executarmos, este código criará 3 arquivos: database.phr, database.pin e database.psq. Juntos, estes arquivos formam o banco de dados, que podem ser referenciado apenas pelo seu prefixo.

Executando uma busca com o BLAST

Agora para realizarmos a comparação da sequência presente no arquivo query.fasta com o banco de dados, tento em vista que tanto as sequência a ser consultado como o banco de dados são de proteínas, é necessário utilizar a classe NcbiblastpCommandline, do mesmo módulo, visto que esta serve de wrapper para o programa BLASTp. Para o construtor dela vamos passar quatro parâmetros:

• query: nome do arquivo que queremos comparar com o banco de dados;
• db: prefixo do banco de dados;
• out: nome do arquivo de saída;
• outfmt: formato de saída desejado. Neste caso vamos utilizar o formato 5, que gera um arquivo XML. Para mais informações sobre os formatos de saída, consulte a documentação do NCBI-BLAST+.

Ao executarmos este código produziremos o arquivo results.xml.

Processando os resultados

Para processarmos os resultados XML produzidos pelos BLAST podemos utilizar o função parse do módulo Bio.Blast.NCBIXML. Esta função recebe um handler de arquivo em modo de leitura e produz um generator, que a cada iteração retorna um o resultado (record) para uma das sequências no arquivo de consulta (query). Cada record possui um ou mais alinhamentos, sendo cada alinhamento referente a uma das sequências que bateu como similar na consulta. Os alinhamento por sua vez possuem listados todos os seus HSPs, que são as regiões onde o BLAST realizou o alinhamento local.

• Os records guardam dados sobre as queries os resultados da consulta, incluindo uma lista de alinhamentos (alignments), que são os hits;
• Cada alignment possui uma ou mais HSPs;
• Os dados de identidade / similaridade, e-value, bitscore, dentre outras coisas, ficam guardados nos HSPs.

Caso queiramos ler o conteúdo do resultado, listando cada HSPs e o seu respectivo e-value, podemos usar o seguinte código:

Próximos passos

O BioPython fornece wrappers e parsers para várias ferramentas além do BLAST, o que facilita em muito a automatização de tarefas em bioinformática, sobretudo em projetos de análise de genomas, transcriptomas e metagenomas. Abordaremos mais sobre estas ferramentas em posts futuros! :D