Preparação para docking com o DOCK 6
Published:
Preparação do sistema para Docking usando DOCK6
Organização
Crie alguns diretórios:
001.initial_files
002.lig_prep
003.rec_prep
004.complex_prep
005.gmx_prep
006.dockprep
007.spheres
008.grid
Salve o arquivo .pdb
e todas as informações iniciais importantes na pasta 001.initial_files
. A pasta 002.lig_prep
será usada para armazenar todas as estruturas com átomos de hidrogênio e cargas parciais para uma equilibração em dinâmica molecular e o experimento de docking.
Preparação de ligante e cofator, se existir
Use o UCSF Chimera para isolar a estrutura do receptor, do ligante e do cofator como
${PDB}.pro.noH.pdb
, para a proteína${PDB}.lig.pdb
, para o ligante${PDB}.cof.pdb
, para o cofator, caso exista.
Salve esses arquivos em 001.initial_files
Ainda usando o Chimera, adicione átomos de hidrogênio em todas as estruturas e as renomeie de acordo:
${PDB}.lig.withH.pdb
, para o ligante${PDB}.cof.withH.pdb
, para o cofator, caso exista.
Observação: Caso não haja cofator, renomeie o arquivo do receptor proteico para ${PDB}.rec.noH.pdb
Crie arquivos .mol2
para o ligante usando o programa Antechamber do pacote AmberTools:
antechamber -fi pdb -fo mol2 -c bcc -at sybyl -i ${PDB}.lig.withH.pdb -o ${PDB}.lig.charged.mol2 -rn LIG -nc ${carga}
Observe que ${carga}
deve ser substituída pela carga do ligante e o mesmo procedimento deve ser feito para o cofator, caso exista. O script fix_charges.py
pode ser usado para consertar pequenas eventuais falhas de ajuste do programa sqm
associado ao Antechamber.
python ${path_to_script}/fix_charges.py -m ${PDB}.lig.charged.mol2
Transfira o arquivo ${PDB}.lig.charged.mol2
para a pasta 002.lig_prep
.
Importante! A ferramenta dockprep
do UCSF Chimera pode determinar a carga total do ligante e preparar os arquivos da simulação, mas não é recomendado para trabalhos com vistas a publicações. Se for usá-lo, use para dados preliminares.
Preparação do receptor
O preparo do receptor exige uma etapa de minimização, equilibração NVT e minimização com a presença do ligante (e do cofator) para que as cadeias laterais se acomodem apropriadamente. Quando o PDB em questão tem resolução < 2 angström, somente a minimização basta. Recomenda-se que estruturas obtidas de cryo-EM passem pelas três etapas de preparo.
Entre no diretório 003.rec_prep
e crie um arquivo para ser lido pelo tleap
:
cat <<EOF > ${PDB}.rec.leap.in
set default PBradii mbondi2
source oldff/leaprc.ff14SB
source leaprc.DNA.bsc1
loadamberparams frcmod.ions234lm_126_tip3p
loadamberparams frcmod.ions1lm_126_tip3p
loadamberparams frcmod.tip3p
REC = loadpdb ../001.initial_files/${PDB}.rec.noH.pdb
saveamberparm REC ${PDB}.rec.prmtop ${PDB}.rec.inpcrd
charge REC
quit
EOF
Ao executar o tleap
, serão gerados dois arquivos, ${PDB}.rec.prmtop
and ${PDB}.rec.inpcrd
:
tleap -f ${PDB}.rec.leap.in
Leia atentamente as mensagens de erro, pois elas indicam modificações que você precisará fazer no arquivo .pdb
.
Atenção! Se o tleap
reclamar algo como:
/Users/guilherme/miniconda3/envs/py37/bin/teLeap: Error!
Could not find bond parameter for: SH - SH
ele está tendo problemas para diferenciar cisteínas que fazem ligações dissulfídicas de cisteínas que não fazem tais ligações. Você deve renomear as cisteínas que fazem ligações S-S de CYS
para CYX
. Para descobrir quais cisteínas formam tal ligação, use o UCSF Chimera, selecione os resíduos CYS e mostre as cadeias laterais.
Outra mensagem que pode aparecer e deve ser lidada é a presença de avisos (“Warnings”) do tipo:
/Users/guilherme/miniconda3/envs/py37/bin/teLeap: Warning!
There is a bond of 11.508708 angstroms between:
Isso é facilmente resolvido com a adição de uma linha TER
entre os resíduos em questão no arquivo .pdb
. Em alguns casos é necessário adicionar grupos neutros (ACE
e NME
) nos terminais, mas deve-se evitá-los a menos que seja necessário (quando estamos interessados em estudar um fragmento específico, por exemplo).
tleap
informa a carga total da proteína em um aviso e no final da sua execução:
/Users/guilherme/miniconda3/envs/py37/bin/teLeap: Warning!
The unperturbed charge of the unit (1.000000) is not zero.
e
Total unperturbed charge: 1.000000
Total perturbed charge: 1.000000
Quit
Atenção! É importante guardar o .pdb
original porque o tleap
renumera todos os resíduos! Esta etapa é necessária para conhecermos as características principais da proteína e consertar quaisquer problemas diretamente no .pdb
.
Preparação do complexo
Entre em 004.complex_prep
na pasta do projeto. De forma semelhante à preparação do receptor, para preparar o complexo, precisamos criar um input para o tleap
:
cat <<EOF > com.leap.in
set default PBradii mbondi2
source oldff/leaprc.ff14SB
source leaprc.gaff2
loadamberparams frcmod.ions234lm_126_tip3p
loadamberparams frcmod.ions1lm_126_tip3p
loadamberparams frcmod.tip3p
REC = loadpdb ../001.initial_files/${PDB}.rec.noH.pdb
loadamberparams ${PDB}.lig.antechamber.frcmod
LIG = loadmol2 ${PDB}.lig.antechamber.mol2
COM = combine { REC LIG }
saveamberparm LIG ${PDB}.lig.prmtop ${PDB}.lig.ori.inpcrd
saveamberparm REC ${PDB}.rec.prmtop ${PDB}.rec.ori.inpcrd
saveamberparm COM ${PDB}.com.prmtop ${PDB}.com.ori.inpcrd
solvcomplex= combine { REC LIG }
solvateoct solvcomplex TIP3PBOX 12.0
addions solvcomplex Cl- 0
addions solvcomplex Na+ 0
savepdb solvcomplex ${PDB}.com.solv.pdb
charge solvcomplex
check solvcomplex
saveamberparm solvcomplex ${PDB}.com.solv.prmtop ${PDB}.com.solv.inpcrd
quit
EOF
Observe, entretanto, que devemos fazer algumas modificações com relação ao ligante. Precisamos de um arquivo .frcmod
e uma versão do .mol2
que não use tipos atômicos SYBYL, mas tipos compatíveis com o Amber. Para isso:
antechamber -i ../002.lig_prep/${PDB}.lig.charged.mol2 -fi mol2 -o ${PDB}.lig.antechamber.mol2 -fo mol2 -dr n
parmchk2 -i ${PDB}.lig.antechamber.mol2 -f mol2 -o ${PDB}.lig.antechamber.frcmod
Se houver cofator:
cat <<EOF > com.leap.in
set default PBradii mbondi2
source oldff/leaprc.ff14SB
source leaprc.gaff2
loadamberparams frcmod.ions234lm_126_tip3p
loadamberparams frcmod.ions1lm_126_tip3p
loadamberparams frcmod.tip3p
PRO = loadpdb ../001.initial_files/${PDB}.pro.noH.pdb
loadamberparams ${PDB}.lig.antechamber.frcmod
LIG = loadmol2 ${PDB}.lig.antechamber.mol2
loadamberparams ${PDB}.cof.antechamber.frcmod
loadamberprep ${PDB}.cof.antechamber.prep
COF = loadpdb ${PDB}.cof.antechamber.pdb
REC = combine { PRO COF }
saveamberparm COF ${PDB}.cof.prmtop ${PDB}.cof.ori.inpcrd
COM = combine { REC LIG }
saveamberparm LIG ${PDB}.lig.prmtop ${PDB}.lig.ori.inpcrd
saveamberparm PRO ${PDB}.pro.prmtop ${PDB}.pro.ori.inpcrd
saveamberparm REC ${PDB}.rec.prmtop ${PDB}.rec.ori.inpcrd
saveamberparm COM ${PDB}.com.prmtop ${PDB}.com.ori.inpcrd
solvcomplex= combine { REC LIG }
solvateoct solvcomplex TIP3PBOX 12.0
addions solvcomplex Cl- 0
addions solvcomplex Na+ 0
savepdb solvcomplex ${PDB}.com.solv.pdb
charge solvcomplex
check solvcomplex
saveamberparm solvcomplex ${PDB}.com.solv.prmtop ${PDB}.com.solv.inpcrd
quit
EOF
e também:
antechamber -i ../002.lig_prep/${PDB}.cof.charged.mol2 -fi mol2 -o ${PDB}.cof.antechamber.prep -fo prepi
antechamber -i ../002.lig_prep/${PDB}.cof.charged.mol2 -fi mol2 -o ${PDB}.cof.antechamber.pdb -fo pdb
parmchk2 -i ${PDB}.cof.antechamber.prep -f prepi -o ${PDB}.cof.antechamber.frcmod
Tendo feito todos esses passos, execute o tleap
e encontrará o seu sistema solvatado e neutralizado em dois arquivos, um .prmtop
e um .inpcrd
que contém respectivamente os parâmetros dos campos de força e as coordenadas de cada átomo do sistema. Use o script amber_to_gmx.py
para criar os arquivos de entrada para o GROMACS, ${PDB}.gro
e ${PDB}.top
:
python ${path_to_script}/amber_to_gmx.py -p ${PDB}.com.solv.prmtop -c ${PDB}.com.solv.inpcrd
Copie os arquivos .gro
e .top
para a pasta 005.gmx_prep
.
Equilibração pré-docking
Para criar os arquivos .mdp
da minimização e da dinâmica, execute o script em 005.gmx_prep
:
bash ${path_to_script}/create_gromacs_mdps.dockprep.sh
No caso de simulações contendo proteínas e ligantes, alguns cuidados a mais devem ser tomados. O primeiro estágio da equilibração, por exemplo, exige que a proteína e o ligante tenham restrições de movimento para que as moléculas de água que envelopam o sistema se assentem. O GROMACS exige que essas restrições sejam incluídas nos arquivos .top
após a definição dos últimos parâmetros de campo de força de cada molécula. A diretiva:
#ifdef POSRES
#include "${nome_do_arquivo_de_restrição}.itp"
#endif
normalmente é incluída no final da definição de cada espécie, antes da nova seção moleculetype
das topologias. Para criar os arquivos .itp
, executamos o seguinte programa do pacote GROMACS:
gmx genrestr -f ${PDB}.gro -o posres.itp
É costume selecionar a opção Protein-H
, isto é, todos os átomos da proteína que não sejam hidrogênio. Outra opção para quem quiser garantir a acomodação das cadeias laterais é restringir somente o esqueleto (opção Backbone
). Para restringir o movimento de átomos pesados do ligante é necessário primeiramente criar um arquivo de índice que indica somente os átomos que não são hidrogênio:
gmx make_ndx -f ${PDB}.gro -o index.ndx
No programa gmx make_ndx
você escolhe a opção do ligante:
> ${número do ligante} & ! a H*
> q
onde ${número do ligante}
é uma das opções impressas na tela. Esse comando selecionará somente os átomos que não sejam hidrogênio do ligante e os distinguirá no arquivo index_lig.ndx
. Você pode criar as restrições para esses átomos usando o gmx genrestr
:
gmx genrestr -f ${PDB}.gro -n index.ndx -o posres_lig.itp
selecionando a opção que você gerou (algo como LIG_&_!H*
. Isso gerará o arquivo desejado. Caso haja cofatores, é possível considerá-los parte da proteína, basta fazer as modificações textuais adequadas nos arquivos .ndx
gerados.
Importante! Os arquivos .itp
precisam ser numerados a partir de 1 e não do número do .gro
do sistema. Se for uma subestrutura da proteína, às vezes pode ser necessário renumerar usando um script de python específico:
python ${path_to_script}/fix_itp.py -i posres.itp
Simulações de dinâmica molecular
Após colocar as cláusulas de restrição no arquivo .top
, execute o GROMACS:
gmx grompp -f 01.minimization.mdp -p ${PDB}.top -c ${PDB}.gro -o 01.minimization.tpr
gmx mdrun -deffnm 01.minimization
gmx grompp -f 02.equil_nvt.mdp -p ${PDB}.top -c 01.minimization.gro -o 02.equil_nvt.tpr -r 01.minimization.gro
gmx mdrun -deffnm 02.equil_nvt
gmx grompp -f 03.minimization.mdp -p ${PDB}.top -c 02.equil_nvt.gro -o 03.minimization.tpr
gmx mdrun -deffnm 03.minimization
Extrair configuração final do resultado do GROMACS
O arquivo 03.minimization.gro
conterá a estrutura minimizada pronta para o docking. Precisamos, primeiramente, centralizar a proteína usando o comando:
gmx trjconv -f 03.minimization.gro -s 03.minimization.tpr -o centered.gro -pbc mol
Abra o arquivo centered.gro
no UCSF Chimera ou UCSF ChimeraX, elimine íons desnecessários e todas as moléculas de água. Salve o complexo proteína e ligante como ${PDB}.ready.pdb
e transfira o arquivo para a pasta 06.dockprep
.
Preparando os arquivos para o uso pelo DOCK6
Entre na pasta 06.dockprep
e abra sando o UCSF Chimera, abra ${PDB}.ready.pdb
. Primeiramente vá em:
Select > Residue > LIG
para selecionar a molécula de ligante. Apague-a com:
Actions > Atoms/Bonds > delete
Remanescerá a estrutura do receptor. Salve-a como ${PDB}.rec.withH.pdb
. Apague todos os átomos de hidrogênio por meio da sequência:
Select > Chemistry > element > H
Actions > Atoms/Bonds > delete
e salve a estrutura restante como ${PDB}.rec.noH.pdb
. Reinicie o Chimera (File > Close Session
) e reabra ${PDB}.ready.pdb
. Desta vez, selecione o ligante conforme as instruções acimae inverta a seleção:
Select > Invert (all models)
Ficará claro que a proteína está selecionada e o ligante não. Apague a proteína da forma exposta acima. O ligante original do complexo servirá de referência para as moléculas ancoradas ou criadas pelo DOCK6 e deve ser preparado como tal. Para isso:
Tools > Structure Editing > Dock Prep
Uma janela se abrirá com diversas opções que se aplicam a proteínas e ligantes. No caso do ligante, se aplicam somente as opções Add hydrogens
, Add charges
e Write Mol2 file
. Clique em OK
e outras janelas se abrirão. Se a sua molécula já tiver hidrogênios, nada acontecerá ao apertar OK
na janela Add Hydrogens for Dock Prep
. Caso não tenha, hidrogênios serão adicionados. A janela seguinte se refere à adição de cargas parciais atômicas. Escolha AM1-BCC e pressione OK
. Uma nova janela aparecerá para confirmar a carga do ligante. Você deve saber qual é a carga da molécula e confira se é o que aparece sob Net Charge
. Aperte OK
e um cálculo quântico terá início para o cálculo das cargas.
Ao terminar, salve o arquivo como ${PDB}.lig.charged.mol2
e use o script em python para corrigir ajustes defeituosos:
python ${path_to_script}/fix_charges.py -m ${PDB}.lig.charged.mol2
Para preparar o receptor, abra ${PDB}.rec.withH.pdb
e faça:
Tools > Structure Editing > Dock Prep
Aperte OK
em todas as janelas e salve o arquivo como ${PDB}.rec.charged.mol2
.
Determinando o sítio de ligação
Use o UCSF Chimera para abrir ${PDB}.rec.noH.pdb
. Mostre a superfície da proteína por meio de:
Actions > Surface > show
Salve a superfície por meio da seguinte sequência:
Tools > Structure Editing > Write DMS
Salve o arquivo como ${PDB}.rec.dms
na pasta 007.spheres
. Entre nessa pasta.
O sítio de ligação é representado por esferas. Para criá-las, temos que usar o programa sphgen
, normalmente instalado junto ao dock6
. Criar as esferas demanda criar um arquivo de input para o sphgen
:
cat <<EOF > INSPH
${PDB}.rec.dms
R
X
0.0
4.0
1.4
${PDB}.rec.sph
EOF
A primeira linha do arquivo INSPH gerado pelo comando acima é a superfície usada como input para a criação das esferas. R
indica que as esferas serão criadas fora da superfície do receptor, X
indica que todos os pontos da superfície serão explorados. 0.0
é a distância em angstrom entre as esferas e a superfície; 4.0
é o raio máximo das esferas e 1.4
é o raio mínimo em angstroms. A última linha de INSPH é o nome do arquivo contendo todas as esferas geradas pelo sphgen
. Para gerar as esferas:
sphgen -i INSPH -o OUTSPH
O arquivo ${PDB}.rec.sph
pode ser aberto no UCSF Chimera. Ele contém esferas ao redor de todo o receptor, o que não é interessante. Precisamos selecionar as esferas que cobrem o ligante original, ${PDB}.lig.charged.mol2
. Para isso usamos outro programa acessório instalado junto ao dock6
, o sphere_selector
:
sphere_selector ${PDB}.rec.sph ../006.dockprep/${PDB}.lig.charged.mol2 8.0
Observe que o caminho para a estrutura do ligante deve ser informada ao sphere_selector
assim como uma distância em angstrom. No caso acima, todas as esferas em até 8.0 angstrom do ligante são selecionadas como parte do sítio de ligação.
Gerando a caixa e o grid de energia potencial
Entre na pasta 008.grid
e crie o arquivo showbox.in
com o comando abaixo: Entre na pasta 008.grid
e crie o arquivo showbox.in
com o comando abaixo:
cat <<EOF > showbox.in
Y
8.0
../007.spheres/selected_spheres.sph
1
${PDB}.box.pdb
EOF
A primeira linha do arquivo showbox.in
significa que desejamos gerar uma caixa. A segunda linha determina que cada lado da caixa deve ter comprimento de 8 angstrom. A terceira linha diz que as esferas em selected_spheres.sph
devem estar no centro da caixa. A última linha determina o nome do arquivo de saída. A caixa é gerada por meio do comando:
showbox < showbox.in
Em posse da caixa, ${PDB}.box.pdb
, podemos gerar o grid que será usado para calcular as interações entre as diversas poses e o receptor. Crie o arquivo:
cat <<EOF > grid.in
compute_grids yes
grid_spacing 0.4
output_molecule no
contact_score no
energy_score yes
energy_cutoff_distance 9999
atom_model a
attractive_exponent 6
repulsive_exponent 12
distance_dielectric yes
dielectric_factor 4
allow_non_integral_charges yes
bump_filter yes
bump_overlap 0.75
receptor_file ../006.dockprep/${PDB}.rec.charged.mol2
box_file ${PDB}.box.pdb
vdw_definition_file ${caminho_para_a_pasta_do_dock6}/parameters/vdw_de_novo.defn
score_grid_prefix grid
EOF
O grid é criado pelo comando:
grid -i grid.in -o grid.out
o programa gerará três arquivos, grid.out
que contém as informações da execução do programa, grid.ngr
e grid.bmp
, binários contendo o grid de energia.
Próximos passos
Pronto! Agora você já pode fazer suas simulações com o dock6
. Para virtual screening ou de novo design, veja tutorial específico.