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.
Protein Data Bank (PDB) é um banco de dados que armazena estruturas tridimensionais de proteínas. O PDB é um banco de dados fundamental nas áreas de biologia estrutural, sendo o alicerce de uma sub-área da Bioinformática conhecida como Bioinformática Estrutural.
O PDB é mantido por uma organização internacional chamada Worldwide Protein Data Bank (wwPDB) e pode ser acessado em: <http://www.wwpdb.org/>. A base de dados também pode ser acessada através das organizações PDBe (http://www.pdbe.org), PDBj (http://www.pdbj.org) ou RCSB (http://www.rcsb.org).
O formato de armazenamento PDB apresenta um arquivo textual, o qual cada linha representa um determinado recurso, como por exemplo, linhas iniciadas com ATOM representam átomos específicos de uma proteína. Abaixo vemos como exemplo um fragmento de um arquivo PDB (código de identificação 1BGA).
HEADER GLYCOSIDASE 04-APR-97 1BGA
ATOM 1 N THR A 2 61.990 84.851 70.908 1.00 32.14 N
ATOM 2 CA THR A 2 62.828 85.531 69.934 1.00 20.14 C
ATOM 3 C THR A 2 62.799 85.087 68.476 1.00 17.99 C
Para mais informações acesse a documentação de formatos de arquivos PDB em: <http://www.wwpdb.org/documentation/file-format>.
Exemplo de arquivo PDB (código 1BGA) – extraído do PDBe (http://www.ebi.ac.uk/pdbe/node/1).
Pacote PDB
Em 2003, Thomas Hamelryck e Bernard Manderick desenvolveram um conjunto de classes para analisar arquivos e estruturas do PDB em Python (pacote Bio.PDB). Eles definiram que uma proteína poderia ser representada através de cinco classes hierárquicas: (i) estrutura; (ii) modelo; (iii) cadeia; (iv) resíduo; e (v) átomo. Sendo possível obter os objetos filhos a partir de seus pais utilizando um identificador:
- modelo (model) = estrutura[0] (structure[0])
- cadeia (chain) = modelo[‘A’] (model[‘A’])
- resíduo (residue) = cadeia[1] (chain[1])
- átomo (atom) = resíduo[‘CA’] (residue[‘CA’])
Eles afirmam que essa estrutura de dados é altamente eficaz para acessar dados em arquivos PDB. Nessa hierarquia, cada modelo pode ter um ou mais estruturas (apesar de modelos cristalográficos obtidos por cristalografia de raio-x apresentarem apenas uma estrutura, os modelos obtidos por ressonância magnética podem apresentar várias); cada estrutura pode ter uma ou mais cadeias, que são representadas através de caracteres (A, B, C, etc.); cada cadeia apresenta diversos resíduos, que podem ser acessados através de sua posição; e cada resíduo apresenta átomos, que são acessados a partir de um código de dois caracteres (ex.: “CA” representa “carbono alfa”).
O pacote PDB requer o pacote NumPy. Para mais informações acesse: <http://www.numpy.org>.
Baixando arquivos PDB
Antes de analisarmos arquivos PDB é necessário baixá-los de bancos de dados públicos. Como exemplo utilizaremos a proteína Beta-glicosidase A do organismo Bacillus polymyxa (código PDB 1BGA). A maneira mais comum de fazer o download é através de um dos Websites que fornecem acesso aos dados do PDB. Por exemplo, você pode acessar <http://www.rcsb.org/pdb/files/1bga.pdb> e fazer o download do arquivo PDB no formato texto PDB File (Text).
Entretanto, Biopython também fornce métodos de download através da classe PDBList e do método retrieve_pdb_file( ).
from Bio.PDB import *
# Fazendo download da estrutura 1BGA
pdb = PDBList()
pdb.retrieve_pdb_file('1BGA')
Ao executar esse código a seguinte mensagem será exibida: “Downloading PDB structure ‘1BGA’…”. Em seguida, uma pasta chamada bg será criada e dentro dela haverá um arquivo chamado “pdb1bga.ent”. Agora vamos analisar esse arquivo.
Primeiro, o renomeie para 1BGA.pdb. Agora crie um novo script chamado “pdb_parser.py” no mesmo diretório do arquivo PDB.
# Importando Bio.PDB
from Bio.PDB import *
# Criando um objeto PDBParser
parser = PDBParser()
# Declarando a estrutura
estrutura = parser.get_structure('BGA', '1BGA.pdb')
No código acima, inicialmente importamos o pacote Bio.PDB, em seguida criamos um objeto PDBParser, necessário para análise do arquivo PDB e, por fim, declaramos a estrutura do arquivo 1BGA.pdb com o nome BGA.
Obtendo informações do cabeçalho
A classe PDBParser fornece métodos para obtenção de amplos aspectos de arquivos PDB, como por exemplo cabeçalhos. Podemos obter as informações de cabeçalhos de um arquivo PDB e armazená-las em uma variável chamada “cabecalho” utilizando a função get_header( ). Entretanto, se você imprimir essa variável obterá um texto no formato json similar a este:
{'structure_method': 'x-ray diffraction', 'head': 'glycosidase', 'journal': 'AUTH J. SANZ-APARICIO, J. A. HERMOSO, M. MARTINEZ-RIPOLL
...
Nessa variável estão armazenados: o método de determinação da estrutura, informações de publicações sobre a proteína, quantidade de cadeias, palavras-chave, organismo, gene, dentre outras informações. Como pode ser percebido, a quantidade de dados armazenada no cabeçalho de um arquivo PDB é imensa. Entretanto é possível acessá-las individualmente.
from Bio.PDB import *
parser = PDBParser()
estrutura = parser.get_structure('BGA', '1BGA.pdb')
cabecalho = parser.get_header()
# Metodo de determinacao da estrutura
metodo = estrutura.header['structure_method']
# Resolucao
resolucao = estrutura.header['resolution']
print("Metodo: ", metodo)
print("Resolucao: ", resolucao)
# Metodo: x-ray diffraction
# Resolucao: 2.4
No exemplo anterior extraímos o método de determinação da estrutura tridimensional da proteína, além da resolução. Apesar de criarmos a variável cabecalho com todas as informações, obtemos as informações diretamente da variável estrutura através das chamadas de estrutura.header[‘structure_method’] para determinação do método, e estrutura.header[‘resolution’] para determinação da resolução.
Navegando na estrutura do arquivo PDB
Agora aprenderemos a navegar pela estrutura de um arquivo PDB. Crie um novo arquivo chamado “navegando.py” e insira o código abaixo:
from Bio.PDB import *
parser = PDBParser()
estrutura = parser.get_structure('BGA', '1BGA.pdb')
for modelo in estrutura:
print("Modelo: ", modelo.id)
for cadeia in modelo:
print("\t - Cadeia: ", cadeia.id)
for residuo in cadeia:
print("\t\t - Residuo: ",residuo.resname, \
"(",residuo.id[1],")" )
for atomo in residuo:
print("\t\t\t - Atomo:",atomo.name, \
"-> Coordenadas: ( X:",atomo.coord[0], \
"- Y:",atomo.coord[1],"- Z:", \
atomo.coord[2],")" )
Ao executar o script acima descrito, você obterá uma lista com todos os modelos, cadeias, resíduos e átomos da estrutura. Agora vamos analisar individualmente algumas partes desse código. No trecho:
for modelo in estrutura:
print("Modelo: ", modelo.id)
Um laço for percorre todos os modelos presentes na estrutura declarada. Como nossa estrutura possui apenas um modelo, o comando print exibirá na tela “Modelo: 0”. No trecho:
for cadeia in modelo:
print("\t - Cadeia: ", cadeia.id)
Para cada modelo existente, um outro laço for busca todas as cadeias presentes. A estrutura “1BGA.pdb” possui quatro cadeias nomeadas de A a D. Assim, será impresso em pontos distintos: “- Cadeia: A”, “- Cadeia: B”, “- Cadeia: C” e “- Cadeia: D”. Observe que inserimos no comando print os caracteres “\t”. Esses caracteres indicam que antes do texto será impresso um espaçamento tabular, que será utilizado para facilitar a visualização dos resultados. Em:
for residuo in cadeia:
print("\t\t - Residuo: ",residuo.resname, \
"(",residuo.id[1],")" )
Para cada cadeia, todos os resíduos serão impressos. Observe que duas variáveis são impressas: (i) nome do resíduo (residuo.resname); e (ii) posição do resíduo (residuo.id[1]). Observe que para aperfeiçoar a visualização dos resultados, desta vez foram inseridos dois espaçamentos (\t\t) para aumentar o recuo e diferenciar a impressão de resíduos da impressão de cadeias, além de parênteses para separar o nome do átomo de sua posição. Por fim, em:
for atomo in residuo:
print("\t\t\t - Atomo:",atomo.name, \
"-> Coordenadas: ( X:",atomo.coord[0], \
"- Y:",atomo.coord[1], \
"- Z:",atomo.coord[2],")" )
Para cada átomo presente em cada resíduo é impresso o nome do átomo (atomo.name), além das coordenadas espaciais X (atomo.coord[0]), Y (atomo.coord[1]) e Z (atomo.coord[2]).
Para cada resíduo, você obterá uma impressão com formatação similar a esta:
- Residuo: GLY ( 318 )
- Atomo: N -> Coordenadas: (X: 21.219 - Y: 33.363 - Z: 83.027)
- Atomo: CA -> Coordenadas: (X: 22.569 - Y: 33.738 - Z: 83.411)
- Atomo: C -> Coordenadas: (X: 23.646 - Y: 33.292 - Z: 82.434)
- Atomo: O -> Coordenadas: (X: 24.817 - Y: 33.589 - Z: 82.64)
Removendo heteroátomos
No exemplo anterior você deve ter notado que diversos resíduos são registrados como “HOH”.
- Residuo: HOH ( 602 )
- Atomo: O -> Coordenadas: (X: 24.992 - Y: 10.822 - Z: 67.761)
Heteroátomos são átomos que, na maioria dos casos, não pertencem necessariamente à estrutura da proteína, entretanto como suas coordenadas foram determinadas durante a obtenção da estrutura, eles são registrados nos arquivos PDB.
Em geral, a maior parte dos heteroátomos são moléculas de água (“HOH”). Assim, em alguns casos convém removê-los. Para isso, vamos alterar o script criado anteriormente e inserir um comando condicional que valida se um resíduo é ou não um heteroátomo. Veja como ficaria:
from Bio.PDB import *
parser = PDBParser()
estrutura = parser.get_structure('BGA', '1BGA.pdb')
for modelo in estrutura:
print("Modelo: ", modelo.id)
for cadeia in modelo:
print("\t - Cadeia: ", cadeia.id)
for residuo in cadeia:
if residuo.id[0] != ' ':
cadeia.detach_child( \
residuo.id)
else:
print("\t\t - Residuo: ",\
residuo.resname,"(", \
residuo.id[1],")")
for atomo in residuo:
print(\
"\t\t\t - Atomo:", atomo.name, \
"-> Coordenadas: ( X:",atomo.coord[0], \
"- Y:",atomo.coord[1], "- Z:",atomo.coord[2],")" )
Nesse exemplo, o comando “cadeia.detach_child(residuo.id)” remove os resíduos que são heteroátomos (determinados por “if residuo.id[0] != ‘ ‘:”).
Renumerando resíduos
Se o posicionamento de resíduos apresentar variações não contínuas, é possível reescrever as numerações. Para isso, insira um contador chamado i, que receberá o valor 1 antes do laço que lê cada resíduo. Em seguida, reescreva a numeração do resíduo com “residue.id = (‘ ‘, i, ‘ ‘)” e incremente o contador com “ i += 1”. Reescreva o código anterior, inserindo os trechos abaixo entre as linhas: “print “\t – Cadeia: “, cadeia.id” e “for atomo in residuo:”. As modificações corrigirão erros de numeração em resíduos.
# [...]
i = 1
for residuo in cadeia:
if residuo.id[0] != ' ':
cadeia.detach_child(residuo.id)
else:
residue.id = (' ', i, ' ')
i += 1
# [...]
Salvando as alterações realizadas
Apesar de termos feito diversas modificações em arquivos PDB, até o momento não salvamos nenhuma delas. Você pode salvar as alterações usando o módulo PDBIO (como exemplo instanciaremos um objeto chamado w), além dos métodos set_structure, que recebe um objeto com uma estrutura de um arquivo PDB, e save, que recebe como argumento o nome do novo arquivo PDB.
w = PDBIO()
w.set_structure(estrutura)
w.save('nova_1BGA.pdb')
Um novo arquivo chamado “nova_1BGA.pdb” será gravado no mesmo diretório em que está o script.
Extraindo outras informações de átomos e resíduos
Além do que foi apresentado até agora, Biopython ainda apresenta uma série de métodos para extração de informações de átomos e resíduos.
# ATOMOS
a.get_name() # nome do atomo
a.get_id() # id
a.get_coord() # coordenadas atomicas
a.get_vector() # coordenadas como vetor
a.get_bfactor() # fator B
a.get_occupancy() # ocupancia
a.get_altloc() # localizacao alternativa
a.get_sigatm() # parametros atomicos
a.get_anisou() # fator B anisotropico
a.get_fullname() # nome completo do atomo
# RESIDUOS
r.get_resname() # retorna o nome do residuo
r.has_id(name) # testa se ha certo atomo
Extraindo a estrutura primária
Podemos extrair a estrutura primária de cadeias de uma proteína usando a classe CaPPBuilder().
Como exemplo, salvaremos a sequência da proteína 1BGA no formato FASTA em um arquivo chamado 1BGA.fasta. Veja como ficaria o script:
from Bio.PDB import *
parser = PDBParser()
peptideos = CaPPBuilder()
estrutura = parser.get_structure('BGA', '1BGA.pdb')
# Criando um arquivo FASTA
w = open("1BGA.fasta","w")
for cadeia in estrutura[0]:
cadeia_atual = cadeia.id
for seq in peptideos.build_peptides(cadeia):
seq_atual = seq.get_sequence()
tamanho_seq_atual = len(seq_atual)
seq_fasta = \
">Cadeia_%s_tamanho_%d\n%s\n" \
%(cadeia_atual,tamanho_seq_atual,seq_atual)
print(seq_fasta)
w.write(seq_fasta)
w.close()
Abrimos o arquivo FASTA na linha “w = open(“1BGA.fasta”,”w”)”. A seguir, iniciamos um laço para percorrer todas as cadeias da primeira estrutura em “for cadeia in estrutura[0]:”; obtemos o nome da cadeia atual; iniciamos um novo laço para obter a estrutura primária em “for seq in peptideos.build_peptides(cadeia):”; obtemos a estrutura primária e o tamanho da sequência; e por fim, formatamos e gravamos os dados da variável “seq_fasta” no arquivo 1BGA.fasta.
Medindo a distância entre dois átomos
Medir a distância entre átomos é importante para diversas tarefas em Bioinformática Estrutural, como por exemplo, o cálculo de contatos. Para isso, Biopython utiliza as coordenadas geográficas do átomo no espaço (X, Y e Z) e calcula a distância euclidiana entre eles.
A seguir, aprenderemos a calcular a distância entre carbonos alfa de dois diferentes resíduos. Entretanto, primeiro aprenderemos a acessar um resíduo específico dentro de uma cadeia.
Pode-se obter um resíduo específico da seguinte forma:
residuo = estrutura[modelo][cadeia][posicao_residuo]
O mesmo é válido para encontrar um átomo específico de um resíduo:
atomo = estrutura[modelo][cadeia][posicao_residuo][atomo]
Com essa informação já é possível calcular de maneira simples a distância entre dois átomos, dada em ångström (Å). Um ångström equivale a dez elevado a menos dez metros, ou seja, 1 Å = 0,0000000001 m.
Agora, veja como calcular a distância entre dois átomos:
from Bio.PDB import *
parser = PDBParser()
estrutura = parser.get_structure('BGA', '1BGA.pdb')
residuo_100_A = estrutura[0]['A'][100]['CA']
residuo_101_A = estrutura[0]['A'][101]['CA']
distancia = residuo_100_A - residuo_101_A
print(distancia)
# 3.79421
Nesse exemplo, é calculada a distância entre os carbonos alfa dos resíduos 100 e 101 da cadeia A (a qual a distância é 3,79421 Å). Biopython fornece a distância euclidiana através de uma simples operação de subtração entre os objetos que armazenam os átomos que se quer saber a distância.
Construindo uma matriz de distância todos contra todos
Agora vamos calcular a distância de todos contra todos os átomos de uma das cadeias e construir uma matriz de distância. Armazenaremos nossa matriz de distância em um arquivo tabular que pode ser acessado através de qualquer editor de planilhas.
from Bio.PDB import *
parser = PDBParser()
estrutura = parser.get_structure('BGA', '1BGA.pdb')
# Resultado sera gravado num arquivo tabular
w = open("matriz_distancia.txt","w")
# Preenchendo a primeira linha
for residuo in estrutura[0]['A']:
if is_aa(residuo):
rid = residuo.id[1]
rname = \
Polypeptide.three_to_one(residuo.resname)
rnameid = "%s%s" %(rid,rname)
w.write("\t")
w.write(rnameid)
w.write("\n")
# Preenchendo a matriz
# comparacao todos contra todos - necessario 2 lacos
for residuo_C1 in estrutura[0]['A']:
if is_aa(residuo_C1):
# Gravamos a primeira coluna
rid = residuo_C1.id[1]
rname = \
Polypeptide.three_to_one(residuo_C1.resname)
rnameid = "%s%s" %(rid,rname)
w.write(rnameid)
w.write("\t")
# Segundo laco
for residuo_C2 in estrutura[0]['A']:
if is_aa(residuo_C2):
# id e carbono alfa
rid1 = residuo_C1.id[1]
rid2 = residuo_C2.id[1]
ca1 = \
estrutura[0]['A'][rid1]['CA']
ca2 = \
estrutura[0]['A'][rid2]['CA']
# calculando distancias \
entre CA
distancia = ca1 - ca2
# imprimindo resultados
w.write(str(distancia))
w.write("\t")
w.write("\n")
w.close()
Ao executar o script acima nada será exibido na tela, entretanto um arquivo chamado matriz_distancia.txt será criado no diretório. Você pode importar esse arquivo para um editor de planilhas (configure tabulações como separadores e marque o tipo de célula como texto para não ter problema com casas decimais). Assim, sua tabela deve apresentar centenas de resultados. Abaixo uma pequena fração da matriz de distância:
T2 | I3 | F4 | Q5 | |
T2 | 0.0 | 3.80041 | 6.55643 | 9.83362 |
I3 | 3.80041 | 0.0 | 3.81633 | 6.36557 |
F4 | 6.55643 | 3.81633 | 0.0 | 3.79293 |
Q5 | 9.83362 | 6.36557 | 3.79293 | 0.0 |
Nesse exemplo vemos que o carbono alfa da treonina na posição dois está a ~3,8 Å do carbono alfa da isoleucina na posição três e a ~9,83 Å do carbono alfa da glutamina na posição cinco. Matrizes de distância entre resíduos podem ter diversas aplicações em Bioinformática Estrutural, entretanto não entraremos em detalhes quanto a seu uso e sim em compreender melhor o código presente no script. Para isso, vamos decompor o script em pequenas partes e analisá-las individualmente.
from Bio.PDB import *
parser = PDBParser()
estrutura = parser.get_structure('BGA', '1BGA.pdb')
# Resultado sera gravado num arquivo tabular
w = open("matriz_distancia.txt","w")
Nesse trecho não temos muitas novidades: todo pacote Bio.PDB é importado, a variável parser recebe um objeto instanciado, que será utilizada para analisar o arquivo PDB e armazená-lo na variável estrutura. Por fim, um novo arquivo chamado matriz_distancia.txt é criado.
# Preenchendo a primeira linha
for residuo in estrutura[0]['A']:
if is_aa(residuo):
rid = residuo.id[1]
rname = \
Polypeptide.three_to_one(residuo.resname)
rnameid = "%s%s" %(rid,rname)
w.write("\t")
w.write(rnameid)
w.write("\n")
O trecho do código acima preenche apenas a primeira linha do arquivo matriz_distancia.txt, inserindo o cabeçalho com o código de uma letra dos aminoácidos e sua posição. Na segunda linha, declaramos um laço que varre a cadeia A do primeiro modelo da proteína e busca todos os resíduos. Nesse exemplo escolhemos trabalhar apenas com a cadeia A. Em seguida temos o uso de uma função ainda não descrita: “is_aa( )”. Essa função determina se o resíduo é um aminoácido ou um heteroátomo. Essa é uma maneira diferente de diferenciar aminoácidos de heteroátomos.
Se for determinado que o resíduo é um aminoácido, a variável rid recebe a posição dele e a variável rname chama a função Polypeptide.three_to_one( ) para converter o código do aminoácido de três caracteres presente em residuo.resname para o código de aminoácidos de um caractere. A variável rnameid une os dados presentes em rid e rname. A seguir, é impresso no arquivo uma tabulação (\t), seguido do conteúdo da variável rnameid. Ao final do laço é impresso uma quebra de linha (\n).
# Preenchendo a matriz
# comparacao todos contra todos - necessario 2 lacos
for residuo_C1 in estrutura[0]['A']:
if is_aa(residuo_C1):
# Gravamos a primeira coluna
rid = residuo_C1.id[1]
rname = \
Polypeptide.three_to_one(residuo_C1.resname)
rnameid = "%s%s" %(rid,rname)
w.write(rnameid)
w.write("\t")
Nessa parte do código iniciamos o preenchimento da matriz. É criado um novo laço que percorre todos os resíduos (chamamos de residuo_C1 para diferenciar do laço anterior), e é feita a validação se o resíduo é um aminoácido. O trecho a seguir é similar ao código anterior, entretanto nesse caso ele será utilizado para gravar os aminoácidos e posições na primeira coluna.
# Segundo laco
for residuo_C2 in estrutura[0]['A']:
if is_aa(residuo_C2):
# id e carbono alfa
rid1 = residuo_C1.id[1]
rid2 = residuo_C2.id[1]
ca1 = estrutura[0]['A'][rid1]['CA']
ca2 = estrutura[0]['A'][rid2]['CA']
# calculando distancias entre CA
distancia = ca1 - ca2
# imprimindo resultados
w.write(str(distancia))
w.write("\t")
w.write("\n")
No segundo laço serão feitas as comparações de distância de todos contra todos. No começo do trecho, o laço é iniciado (a variável agora é chamada de residuo_C2) e, mais uma vez, o resíduo é validado. A posição é recebida nas variáveis rid1 (para residuo_C1) e rid2 (para residuo_C2), e a posição será utilizada para determinar as posições dos carbonos alfa nas variáveis ca1 e ca2. O cálculo da distância é feito no trecho “distancia = ca1 – ca2”, e esse valor é impresso devidamente formatado no arquivo nas linhas seguintes. Observe que foi utilizado a função str( ) na gravação da variável distancia. Essa função força a gravação da variável como um string (indicado para evitar problemas com gravação de variáveis flutuantes).
Após o laço, o arquivo é fechado com segurança (w.close( )) e poderá ser aberto por outra aplicação.
Aprendendo mais sobre arquivos PDB
Para mais funções acesse a documentação oficial do módulo Bio.PDB em <http://biopython.org/DIST/docs/api/Bio.PDB-module.html> ou a página com as perguntas mais frequentes do sobre o módulo em <http://biopython.org/wiki/The_Biopython_Structural_Bioinformatics_FAQ>.
Quer aprender mais? Conheça nossos cursos profissionalizantes à partir de R$19,99:
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