Categorias
Artigos Introdução à Programação para Bioinformática com Biopython

BLAST

Capítulo 11

Este conteúdo faz parte do livro “Introdução à Programação para Bioinformática com Biopython“. Você pode adquirir a versão impressa desse livro aqui ou a versão para Kindle aqui. Para nos citar, consulte este link.

BLAST (Basic Local Alignment Search Tool ou Ferramenta de Busca de Alinhamento Local Básico) é uma ferramenta de comparação entre sequências primárias (nucleotídeos ou aminoácidos) que utiliza um algoritmo de programação dinâmica para alinhamento local. BLAST apresenta uma metodologia de grande eficiência para comparação par-a-par entre sequências, quanto comparação de uma sequência com uma grande base de dados.

A versão mais recente do BLAST foi introduzida em 2009 pelo NCBI. A suíte de aplicativos stand-alone BLAST+ foi construída com a linguagem C++, o que permitiu uma maior performance e uma melhoria nos recursos das aplicações. Assim, a suíte apresenta cinco aplicativos: blastn, blastp, blastx, tblastn e tblastx.

blastn

Compara uma sequência de nucleotídeos contra uma base de dados de sequência de nucleotídeos.

blastp

Compara uma sequência de aminoácidos contra uma base de dados de sequência de proteínas.

blastx

Compara uma sequência de nucleotídeos traduzida, em todos os quadros de leitura, contra uma base de dados de sequência de proteínas.

tblastn

Compara uma sequência de proteína contra uma base de dados de sequência de nucleotídeos traduzida dinamicamente em todos os quadros de leitura.

tblastx

Compara os seis quadros de leitura traduzidos de uma sequência de nucleotídeos contra os seis quadros de leitura traduzidos de uma base de dados de sequência de nucleotídeos.

Em Bioinformática, BLAST é uma ferramenta tão difundida que poderíamos afirmar que BLAST é a ferramenta favorita de 11 em cada 10 bioinformatas (BLAST S2). Brincadeiras a parte, Biopython apresenta o pacote Bio.Blast que permite o uso do BLAST tanto com bases de dados locais quanto online.

Executando BLAST pela internet

BLAST pode ser executado online pelo Website oficial no NCBI <http://blast.ncbi.nlm.nih.gov>, entretanto o módulo Bio.Blast possui o método qblast que permite a execução de buscas através de scripts em Python.

A sintaxe do método é: “NCBIWWW.qblast(parâmetros)”. Sendo os parâmetros:

program

blastn, blastp, blastx, tblastn ou tblastx.

database

Base de dados a qual será feita a busca (ex.: nr).

sequence

Sequência que será buscada.

descriptions

Número de descrições que serão exibidas (valor padrão: 500).

alignments

Número de alinhamentos que serão exibidos (valor padrão: 500).

expect

Valor de cutoff esperado (valor padrão: 10.0).

matrix_name

Especificar uma matriz de pontuação (PAM30, PAM70, BLOSUM80, BLOSUM45).

filter

Valor “nome” desativa filtragem.

format_type

“HTML”, “Text”, “ASN.1” ou “XML” (valor padrão: “XML”).

hitlist_size

Número de resultados exibidos (valor padrão: 50).

megablast

TRUE/FALSE: Ativa o algoritmo MEGABLAST (exclusivo para blastn).

service

plain, psi, phi, rpsblast, megablast.

Para testarmos as buscas na base de dados nt do NCBI com BLAST através do método qblast( ), primeiro crie um arquivo chamado exemplo.fasta.

>beta-glucosidase_parcial_mRNA
CGTGTGCTTGCGCTACTGCTTGGTCTCACTGTGTGAA
CAACTGCTGCCGCCGCCTGCCGGCGGTGTGCTC

Agora vamos criar um script chamado blast.py (script c11_s1.py no GitHub).

from Bio.Blast import NCBIWWW
from Bio import SeqIO

arquivo = SeqIO.read("exemplo.fasta", format="fasta")

print("Iniciando busca no NCBIWWW..." )

resultado = NCBIWWW.qblast("blastn", "nt", \
arquivo.seq, format_type="Text")

print(resultado.read())

print("Busca concluida." )

Se você obteve ao final da execução a mensagem “Busca concluida”, sua consulta foi realizada com sucesso. Se o script travou com a mensagem “Iniciando busca no NCBIWWW…”, seu script provavelmente obteve algum problema ao conectar-se a banco de dados do NCBI. Recomendamos que aguarde alguns minutos. Se ainda não obtiver nenhum resultado tente interromper a execução do script (comando CONTROL + C) e executá-lo novamente (se estiver utilizando o Sublime Text pode ser necessário finalizar o aplicativo e abri-lo novamente). Lembre-se que é necessária conexão com a internet. Agora vamos entender melhor o que foi realizado nesse script.

No início do código, importamos o módulo NCBIWWW do pacote Bio.Blast e também o pacote SeqIO. O arquivo FASTA foi lido e aplicado à variável arquivo. O método qblast( ) recebeu como parâmetros o programa blastn, a base de dados nt, a sequência (seq) armazenada na variável arquivo e o tipo de formato a qual seria exibido o resultado: “Text”. O resultado da consulta foi armazenado na variável resultado e foi impresso através do comando “print resultado.read( )”. Observe que inserimos uma mensagem antes da execução do comando qblast( ), indicando o início da busca, e outra após a execução, avisando sobre o término. A impressão de mensagens na tela é importante para detectarmos em qual ponto do código há maior tempo de execução. Por se tratar de um comando que faz uma consulta remota à base de dados do NCBI, não podemos definir com precisão qual o tempo total de execução. Em nossos testes, obtivemos consultas realizadas entre 9 s e 267 s, entretanto esse número irá variar de acordo com a conexão com a internet do usuário.

O formato Text nos dá uma visualização simples dos resultados, entretanto podemos obter uma visualização mais eficiente utilizando o formato XML. Entretanto, a análise de arquivos no formato XML é um pouco mais complexa.

Analisando resultados de BLAST no formato XML

Vamos alterar o script anterior para obter como resultado um arquivo no formato XML. Por fim, vamos ler esse arquivo XML e imprimi-lo formatado na tela.

from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML
from Bio import SeqIO

arquivo = SeqIO.read("exemplo.fasta", format="fasta")

print("Iniciando busca no NCBIWWW..." )

resultado = NCBIWWW.qblast("blastn", "nt", \
arquivo.seq, format_type="XML")

print("Busca concluida. Salvando resultados..." )

saida = open("blast_resultado.xml", "w")
saida.write(resultado.read())
saida.close() 

# Lendo arquivos XML
arquivo_xml = open("blast_resultado.xml","r")
dados = NCBIXML.parse(arquivo_xml)
item = next(dados)

i = 1

for a in item.alignments:
	for hsp in a.hsps:
		print('Alinhamento',i)
		print('Sequencia: '+a.title)
		print('Tamanho: ',a.length)
		print('Score: ',hsp.score)
		print('Gaps: ',hsp.gaps)
		print(hsp.query )
		print(hsp.match)
		print(hsp.sbjct)
		print("\n")
		i+=1

print("Executado com sucesso." )

Para a análise de arquivos em XML foi necessário importar o módulo NCBIXML. Esse módulo possui dois métodos: read para leitura de resultados de BLAST em arquivos XML e parse para análise iterativa de resultados. Em nosso exemplo anterior, após o fechamento do arquivo XML utilizamos a função open( ) para reabrí-lo novamente. Você deve estar pensando que essa não é a melhor maneira de realizar esta tarefa, e tem razão, afinal por que fecharíamos um arquivo e na próxima linha o abriríamos novamente? O programa perde em desempenho, entretanto essa é a forma mais didática de explicar como analisar arquivos XML.

Após a abertura do arquivo XML os comandos “dados = NCBIXML.parse(arquivo_xml)” e “item = next(dados)” permitem a ler o arquivo XML iterativamente e receber o próximo item em uma iteração, respectivamente. É declarada uma variável i que recebe o valor 1. Essa variável será utilizada como contador na impressão dos alinhamentos. Em seguida, dois laços de repetição percorrem os itens da lista de alinhamentos e para cada alinhamento é analisado o conteúdo de hsp (high-scoring segment pairs ou pares de segmentos de alta pontuação), onde estão contidos os resultados de pontuações de alinhamento (score), falhas no alinhamento (gaps), sequências buscadas (query), sequências correspondentes nos bancos de dados (sbjct) e as partes semelhantes entre sequências (match).

Se tudo ocorrer sem erros, o usuário recebe na tela mensagem: “Executado com sucesso.”.

Realizando BLAST local

Executar o BLAST localmente permite que suas buscas sejam realizadas com menor tempo. Ainda permite que você construa bases de dados personalizadas. Para executar o BLAST local, Biopython requer que a suíte de aplicativos BLAST+ esteja instalada. A última versão da suíte BLAST+ pode ser obtida em: <ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/> e tem suporte para sistemas operacionais Windows, Linux e MacOS.

Biopython trata o BLAST como um Wrapper (uma ferramenta que pode ser acessada através de outra ferramenta), e fornece métodos de acesso através do módulo Applications do pacote Blast. Para isso, Biopython permite a construção de linhas de comando que executarão o BLAST, e por fim, essas serão executadas através do Python.

Como exemplo, vamos comparar duas sequências de nucleotídeos. Faça uma cópia do arquivo exemplo.fasta, criado anteriormente (ou obtê-lo no diretório exemplos do GitHub). Pode chamá-lo de exemplo2.fasta. Faça algumas modificações no arquivo, como por exemplo, insira um número dois no cabeçalho e altere o quinto nucleotídeo de “T” para “A”. Agora vamos compará-los com BLAST. Veja:

from Bio.Blast.Applications import *
 
comando_blastn = NcbiblastnCommandline( \
query="exemplo.fasta", subject="exemplo2.fasta", \
outfmt=0, out="out.txt")
print(comando_blastn)
 
# Executando
stdout, stderr = comando_blastn()
 
# Abrindo resultado
blast_result = open("out.txt","r")
 
linhas = blast_result.read()
print(linhas)

É possível determinar se o BLAST foi executado com sucesso se você obtiver como resultado algo similar a:

/usr/local/ncbi/blast/bin/blastn -out out.txt -outfmt 0 -query exemplo.fasta -subject exemplo2.fasta

BLASTN 2.2.27+

Query= beta-glucosidase_parcial_mRNA

...

Entretanto, se você obtiver como resultado: “blastn: command not found”, provavelmente o comando blastn não está declarado na variável $PATH. Você pode corrigir esse problema inserindo no seu script uma variável com localização do executável blastn e declarando essa variável no método NcbiblastnCommandline( ).

Em geral, BLAST+ é instalado nos diretórios:

  • /usr/bin/blastn” no Linux;
  • /usr/local/ncbi/blast/bin/blastn” no MacOS;
  • C:\Program Files\NCBI\blast-2.2.31+\bin\blastn” no Windows.

Entretanto, o endereço dos diretórios pode variar de acordo com a versão do BLAST ou do Sistema Operacional. Você pode tentar encontrá-lo via terminal utilizando os seguintes comandos:

  • whereis blastn” ou “locate blastn” no Linux;
  • mdfind “kMDItemFSName=blastn”’ no MacOS;
  • no Windows abra um diretório qualquer, pressione F3 e busque por “blastn”.

Por exemplo, no MacOS o script correto ficaria desta forma:

from Bio.Blast.Applications import *
 
blastn = "/usr/local/ncbi/blast/bin/blastn" 
 
comando_blastn = NcbiblastnCommandline(cmd=blastn, \
query="exemplo.fasta", subject="exemplo2.fasta", \
outfmt=0, out="out.txt")
print(comando_blastn)
 
# Executando
stdout, stderr = comando_blastn()
 
# Abrindo resultado
blast_result = open("out.txt","r")
 
linhas = blast_result.read()
print(linhas)

Caso deseje executar esse script em um sistema operacional diferente, altere o valor da variável blastn.

Nesse script gostaríamos de destacar a maneira com a qual a linha de comando foi executada. A classe NcbiblastnCommandline( ) inicialmente gera a linha de comando. Assim, a execução só é realizada com a chamada do objeto criado, e quando isso ocorre duas variáveis recebem os resultados dessa chamada: (i) stdout, que recebe a saída da execução, e (ii) stderr, que recebe o log de erros.

Parâmetros do BLAST

Parâmetros do BLAST podem ser aplicados a classe NcbiblastnCommandline( ). Veja alguns parâmetros abaixo:

db

Nome do banco de dados BLAST. Input: string; valor padrão: none (nenhum).

query

Nome de um arquivo para busca. Input: string; valor padrão: requer o nome de um arquivo.

out

Nome do arquivo onde os resultados de BLAST serão salvos. Input: string; valor padrão: stdout.

evalue

E-value. Input: real; valor padrão: 10.0.

subject

Arquivo com sequências onde será realizada a busca. Input: string; valor padrão: none.

num_threads

Número de threads (CPUs) usadas em buscas pelo BLAST. Input: inteiro; valor padrão: 1.

word_size

Tamanho mínimo do fragmento de sequência usado na busca. Input: valor inteiro; valores padrões: (i) blastn: 11; (ii) blastp: 3; (iii) blastx: 3; (iv) tblastn: 3; (v) tblastx: 3; (vi) rpsblast: 40; (vii) megablast: 28;

matrix

Define a matriz de pontuação usada no alinhamento (usado em alinhamento de proteínas). Input: string; valor padrão: “BLOSSUM 62”.

outfmt

Opções de visualização de alinhamentos. Pode variar de 0 a 14. Input: string; valor padrão: “0”.

0 = par a par (textual),

1 = query-anchored mostrando identidades,

2 = query-anchored sem identidades,

3 = flat query-anchored (mostra identidade),

4 = flat query-anchored (sem identidades),

5 = XML,

6 = tabular,

7 = tabular com linhas comentadas,

8 = Textual ASN.1,

9 = Binário ASN.1,

10 = Valores separados por vírgula,

11 = Formato de arquivo BLAST (ASN.1),

12 = Saída JSON Seqalign,

13 = Saída JSON Blast,

14 = Saída XML2 Blast

Resultados de BLAST em formato tabular

Até agora, analisamos resultados de BLAST em formato textual par a par e em formato XML, entretanto o formato tabular (6) pode auxiliar na visualização em paralelo de diversos resultados, podendo ser aberto através de editores de planilhas, como MS Office Excel ou o Calc, ou facilmente analisado através de scripts. Por padrão o formato tabular exibe os seguintes resultados: ‘qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore‘, entretanto outras opções podem ser adicionadas. Vamos entender o que cada um desses elementos representa:

qseqid

Nome da sequência usada na busca (query). Pode ser utilizada mais de uma sequência (arquivo Multi-FASTA).

qlen

Tamanho da sequência usada na busca (query).

sseqid

Nome da sequência a qual é realizada a busca (subject). Pode ser utilizada mais de uma sequência (arquivo Multi-FASTA).

slen

Tamanho da sequência a qual é realizada a busca (subject).

qstart

Posição de início do alinhamento na query.

qend

Posição de término do alinhamento na query.

sstart

Posição de início do alinhamento no subject.

send

Posição de término do alinhamento no subject.

evalue

Valor mínimo esperado para o evalue.

score

Valor de pontuação.

length

Indica tamanho do alinhamento.

pident

Indica o percentual de correspondências idênticas (nucleotídeos/aminoácidos).

nident

Indica o número de correspondências idênticas (nucleotídeos/aminoácidos).

mismatch

Indica o número de mismatches (nucleotídeos/aminoácidos diferentes em um alinhamento).

positive

Número de nucleotídeos/aminoácidos idênticos.

gapopen

Número de gaps abertos.

gaps

Número de gaps existentes.

ppos

Percentual de nucleotídeos/aminoácidos classificados como positivos.

bitscore

Pontuação bit.

No exemplo abaixo, vamos gravar em um arquivo os atributos: qlen, qstart, qend, slen, sstart, send, positive, gaps, mismatch e pident.

from Bio.Blast.Applications import *
 
comando_blastn = NcbiblastnCommandline(query= \
"example.fasta", subject="example2.fasta", \
outfmt="'6 qlen qstart qend slen sstart send \
positive gaps mismatch pident'", out="out.txt")
print(comando_blastn)
 
# Executando
stdout, stderr = comando_blastn()
 
# Abrindo resultado
blast_result = open("out.txt","r")
 
linhas = blast_result.read()
print(linhas)

Observe a forma a qual o parâmetro outfmt é inserido na classe NcbiblastnCommandline( ). O valor contido em outfmt é registrado com aspas duplas seguido de aspas simples. Isso ocorre pois Biopython requer aspas duplas no parâmetro, e registra o conteúdo da variável como do tipo inteiro. Ao inserir as aspas simples após as aspas duplas, forçamos Biopython a enviar o argumento como string. Isso só é possível devido ao fato de estamos utilizando um wrapper e do BLAST aceitar como argumentos valores inteiros e strings no campo outfmt.

Se você executou o script com os arquivos de exemplo criados anteriormente, você deve ter obtido um resultado similar a este:

blastn -out out.txt -outfmt '6 qlen qstart qend slen sstart send positive gaps mismatch pident' -query exemplo.fasta -subject exemplo2.fasta
70	1	70	70	1	70	69	0	1	98.57

Na primeira linha temos a impressão do comando criado. A seguir os resultados na respectiva ordem: tamanho da sequência query (70), posição primeiro nucleotídeo do alinhamento na query (1), posição do último nucleotídeo do alinhamento na query (70), tamanho da sequência subject (70), posição primeiro nucleotídeo do alinhamento no subject (1), posição do último nucleotídeo do alinhamento no subject (70), número de nucleotídeos idênticos (69), número de gaps (0), número de mismatches (1) e percentual de identidade (98,57%).

Outras classes do Bio.Blast.Applications

Nos exemplos demonstrados anteriormente, apresentamos alinhamentos de nucleotídeos com a ferramenta blastn, entretanto é possível realizar alinhamentos utilizando as outras ferramentas da suíte BLAST+. Para isso é necessário apenas algumas alterações nos métodos utilizados. Segue abaixo, algumas classes do módulo Applications do pacote Bio.Blast:

Ncbiblastp

Commandline( )

Permite alinhamentos de sequências de proteínas com outras sequências de proteínas (blastp).

Ncbiblastx

Commandline( )

Permite a comparação de uma sequência de nucleotídeos traduzidas, em todos os quadros de leitura, contra uma base de dados de sequências de proteínas (blastx).

Ncbitblastn

Commandline( )

Permite a comparação de uma sequência de proteína contra uma base de dados de sequência de nucleotídeos traduzida dinamicamente em todos os quadros de leitura (tblastn).

Ncbitblastx

Commandline( )

Permite a comparação de seis quadros de leitura traduzidos de uma sequência de nucleotídeos contra os seis quadros de leitura traduzidos de um banco de dados sequência de nucleotídeos (tblastx).

Ncbipsiblast

Commandline( )

Permite o uso do programa psiblast. Exemplo:

NcbipsiblastCommandline(cmd=’psiblast’, help=True)

Ncbirpsblast

Commandline( )

Permite o uso do programa rpsblast. Exemplo:

NcbirpsblastCommandline(cmd=rpsblast, help=True)

Ncbirpstblastn

Commandline( )

Permite o uso do programa rpstblastn. Exemplo:

NcbirpstblastnCommandline(cmd=rpstblastn, help=True)

Ncbideltablast

Commandline( )

Permite o uso do programa deltablast (para proteínas). Exemplo:

NcbideltablastCommandline(cmd=’deltablast’, query=’rosemary.pro’, db=’nr’, evalue=0.001, remote=True)

Ncbiblastforma

tterCommand

line( )

Permite a conversão de arquivos ASN.1 para outros formatos de saída BLAST. Exemplo:

NcbiblastformatterCommandline(cmd=’blast_formatter’, out=’exemplo.xml’, outfmt=5, archive=’exemplo.asn’)

Para mais informações sobre o pacote Bio.Blast, acesse: <http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc84>.

Construindo seu próprio banco de dados (opcional)

Em alguns casos, comparações entre uma sequência (query) e um grande quantidade de sequências (subject), como por exemplo um arquivo Multi-FASTA com centenas de milhares de sequências, pode demandar uma grande quantidade de tempo. Assim, torna-se viável a criação de bancos de dados que facilitem o trabalho de busca.

Infelizmente, Biopython não fornece métodos para criação de bancos de dados, entretanto é possível utilizar o comando makeblastdb da suíte BLAST+ para realizar essa tarefa.

Como exemplo, vamos converter nosso arquivo exemplo.fasta em um banco de dados (Linux e MacOS X; para windows acrescente o endereço completo do programa makeblastdb).

makeblastdb -in exemplo.fasta -dbtype 'nucl' -out exemplo_db -title exemplo_db

Você perceberá que três novos arquivos foram criados: exemplo_db.nhr, exemplo_db.nin e exemplo_db.nsq.

Agora vamos modificar o script usado anteriormente, substituindo o parâmetro subject por db e acrescentando o nome do banco de dados criado anteriormente: exemplo_db. Não se esqueça de alterar o nome da query, afinal criamos nossa base de dados com o arquivo exemplo.fasta, logo faz mais sentido compararmos o arquivo exemplo2.fasta com o banco de dados criado.

from Bio.Blast.Applications import *

comando_blastn = NcbiblastnCommandline( \
query="example2.fasta", db="exemplo_db", \
outfmt="'6 qlen qstart qend slen sstart send \
positive gaps mismatch pident'", out="out.txt")
print(comando_blastn)
 
# Executando
stdout, stderr = comando_blastn()
 
# Abrindo resultado
blast_result = open("out.txt","r")
 
linhas = blast_result.read()
print(linhas)

Você notou que o script pode ser executado mais rápido? É claro que não! Diferenças só serão perceptíveis com grandes bancos de dados.

Quer aprender mais? Conheça nossos cursos profissionalizantes à partir de R$19,99:

Livro Introdução à Programação para Bioinformática com Biopython

Capítulo 1
Introdução ao Python

Capítulo 2
Comandos condicionais

Capítulo 3
Laços de repetição

Capítulo 4
Trabalhando com strings

Capítulo 5
Listas

Capítulo 6
Manipulando arquivos

Capítulo 7
Funções

Capítulo 8
Princípios da orientação a objetos

Capítulo 9
Introdução ao Biopython

Capítulo 10
Sequências

Capítulo 11
BLAST

Capítulo 12
PDB

Capítulo 13
Visualização de dados em Python

Capítulo 14
Outras coisas interessantes que se pode fazer com Biopython

Capítulo 15
Hierarquia do Biopython

Por favor, nos cite:

MARIANO, D. C. B.; BARROSO, J. R. P. M. ; CORREIA, T. S. ; de MELO-MINARDI, R. C. . Introdução à Programação para Bioinformática com Biopython. 3. ed. North Charleston, SC (EUA): CreateSpace Independent Publishing Platform, 2015. v. 1. 230p .

Por Diego Mariano

Doutor em Bioinformática pela Universidade Federal de Minas Gerais com atuação na área de ciência de dados e aprendizado de máquina aplicados ao aperfeiçoamento de enzimas usadas na produção de biocombustíveis. Mestre em Bioinformática, também pela UFMG, atuando na área de desenvolvimento de sistemas Web para montagem de genomas. Atualmente realiza estágio pós-doutoral no Departamento de Ciência da Computação da UFMG com foco em desenvolvimento de sistemas Web para Bioinformática, análise exploratória e visualização de dados. Tem conhecimentos nas linguagens: PHP, JavaScript, Python, R, Perl, HTML, CSS e SQL.