Cómo modelar proteínas con colabfold en tu GPU local

Hola,

hoy explicaré cómo he configurado ColabFold para ejecutarlo en hardware local, en concreto en una máquina con Ubuntu 20.04 que tiene una CPU Xeon CascadeLake Silver 4210R y una tarjeta gráfica NVIDIA RTX 3090. Puedes leer más sobre AlphaFold y ColabFold aquí o en este vídeo.

1) Necesité actualizar cuda, en concreto con la versión 11.8, algo que hice como se explica aquí:

wget https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2004/x86_64/cuda-ubuntu2004.pin
sudo mv cuda-ubuntu2004.pin /etc/apt/preferences.d/cuda-repository-pin-600
sudo apt-key adv --fetch-keys https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2004/x86_64/3bf863cc.pub
sudo add-apt-repository "deb http://developer.download.nvidia.com/compute/cuda/repos/ubuntu2004/x86_64/ /"
sudo apt update
sudo apt install cuda-toolkit-11-8

2) Tras reinicar, actualicé la variable de ambiente $PATH añadiendo estas líneas a mi fichero .bashrc:

export PATH=/usr/local/cuda/bin:$PATH

export LD_LIBRARY_PATH=/usr/local/cuda/lib64:$LD_LIBRARY_PATH

3) Seguí las instrucciones para Linux en https://github.com/YoshitakaMo/localcolabfold?tab=readme-ov-file#for-linux . En mi caso tardó unos pocos minutos y sumó 15G al disco duro. 

4) Probé que todo funciona con un fichero FASTA qee contiene varias secuencias, guardando los resultados en la carpeta multi/ :

colabfold_batch test.multi.faa multi/

Ahora resumo los resultados que obtuve:

• Por defecto colabfold_batch se conecta a https://api.colabfold.com para hacer búsquedas de secuencias similares y construir alineamientos múltiples (MSA) en un formato similar a FASTA que se llama a3m. Por tanto esa parte del trabajo no se hace localmente y tendrás que usarla con medida. Si quieres saber qué versión de las bases de datos de secuencias de ColabFold estás usando puedes consultar https://github.com/sokrypton/ColabFold/wiki/MSA-Server-Database-History
• Las primeras secuencias que usé para construir modelos en formato PDB tenían entre 114 y 162 resíduos y tardaban un par de minutos, pego aquí el log: 
• 2024-02-26 13:05:56,639 Running colabfold 1.5.5 (d36504fad856a0e1df511c5b0434957707030319)
  2024-02-26 13:05:56,862 Running on GPU
  2024-02-26 13:05:57,354 Found 5 citations for tools or databases
  2024-02-26 13:05:57,355 Query 1/29: test1 (length 114)
  2024-02-26 13:05:58,348 Sleeping for 6s. Reason: PENDING
  2024-02-26 13:06:05,308 Sleeping for 10s. Reason: RUNNING
  2024-02-26 13:06:30,822 Padding length to 124
  2024-02-26 13:06:58,791 alphafold2_ptm_model_1_seed_000 recycle=0 pLDDT=67.9 pTM=0.31
  2024-02-26 13:07:00,321 alphafold2_ptm_model_1_seed_000 recycle=1 pLDDT=68.8 pTM=0.329 tol=9.09
  2024-02-26 13:07:01,845 alphafold2_ptm_model_1_seed_000 recycle=2 pLDDT=69.7 pTM=0.358 tol=2.28
  2024-02-26 13:07:03,373 alphafold2_ptm_model_1_seed_000 recycle=3 pLDDT=69.8 pTM=0.367 tol=3.04
  2024-02-26 13:07:03,374 alphafold2_ptm_model_1_seed_000 took 32.6s (3 recycles)
  2024-02-26 13:07:04,871 alphafold2_ptm_model_2_seed_000 recycle=0 pLDDT=71.2 pTM=0.308
  2024-02-26 13:07:06,323 alphafold2_ptm_model_2_seed_000 recycle=1 pLDDT=71.6 pTM=0.346 tol=2.14
  2024-02-26 13:07:07,848 alphafold2_ptm_model_2_seed_000 recycle=2 pLDDT=71.7 pTM=0.358 tol=2.38
  2024-02-26 13:07:09,345 alphafold2_ptm_model_2_seed_000 recycle=3 pLDDT=71.8 pTM=0.365 tol=1.31
  2024-02-26 13:07:09,346 alphafold2_ptm_model_2_seed_000 took 5.9s (3 recycles)
  2024-02-26 13:07:10,984 alphafold2_ptm_model_3_seed_000 recycle=0 pLDDT=68.1 pTM=0.298
  2024-02-26 13:07:12,529 alphafold2_ptm_model_3_seed_000 recycle=1 pLDDT=68.6 pTM=0.34 tol=4.11
  2024-02-26 13:07:13,992 alphafold2_ptm_model_3_seed_000 recycle=2 pLDDT=69.2 pTM=0.36 tol=2.49
  2024-02-26 13:07:15,484 alphafold2_ptm_model_3_seed_000 recycle=3 pLDDT=68.8 pTM=0.367 tol=1.67
  2024-02-26 13:07:15,485 alphafold2_ptm_model_3_seed_000 took 6.1s (3 recycles)
  2024-02-26 13:07:16,987 alphafold2_ptm_model_4_seed_000 recycle=0 pLDDT=66.1 pTM=0.289
  2024-02-26 13:07:18,435 alphafold2_ptm_model_4_seed_000 recycle=1 pLDDT=66.8 pTM=0.283 tol=5.61
  2024-02-26 13:07:19,933 alphafold2_ptm_model_4_seed_000 recycle=2 pLDDT=67.7 pTM=0.298 tol=1.03
  2024-02-26 13:07:21,444 alphafold2_ptm_model_4_seed_000 recycle=3 pLDDT=67.9 pTM=0.318 tol=2.04
  2024-02-26 13:07:21,445 alphafold2_ptm_model_4_seed_000 took 5.9s (3 recycles)
  2024-02-26 13:07:22,931 alphafold2_ptm_model_5_seed_000 recycle=0 pLDDT=66.8 pTM=0.322
  2024-02-26 13:07:24,403 alphafold2_ptm_model_5_seed_000 recycle=1 pLDDT=68.2 pTM=0.345 tol=9.46
  2024-02-26 13:07:25,860 alphafold2_ptm_model_5_seed_000 recycle=2 pLDDT=68.8 pTM=0.354 tol=2.3
  2024-02-26 13:07:27,342 alphafold2_ptm_model_5_seed_000 recycle=3 pLDDT=69.4 pTM=0.358 tol=1.58
  2024-02-26 13:07:27,342 alphafold2_ptm_model_5_seed_000 took 5.9s (3 recycles)
  2024-02-26 13:07:27,369 reranking models by 'plddt' metric
  2024-02-26 13:07:27,369 rank_001_alphafold2_ptm_model_2_seed_000 pLDDT=71.8 pTM=0.365
  2024-02-26 13:07:27,369 rank_002_alphafold2_ptm_model_1_seed_000 pLDDT=69.8 pTM=0.367
  2024-02-26 13:07:27,370 rank_003_alphafold2_ptm_model_5_seed_000 pLDDT=69.4 pTM=0.358
  2024-02-26 13:07:27,370 rank_004_alphafold2_ptm_model_3_seed_000 pLDDT=68.8 pTM=0.367
  2024-02-26 13:07:27,370 rank_005_alphafold2_ptm_model_4_seed_000 pLDDT=67.9 pTM=0.318
  2024-02-26 13:07:28,679 Query 2/29: test2 (length 120)
  2024-02-26 13:07:29,695 Sleeping for 9s. Reason: PENDING
  2024-02-26 13:07:39,667 Sleeping for 9s. Reason: PENDING
  2024-02-26 13:07:49,628 Sleeping for 6s. Reason: PENDING
  2024-02-26 13:07:56,610 Sleeping for 6s. Reason: PENDING
  2024-02-26 13:08:03,608 Sleeping for 5s. Reason: PENDING
  2024-02-26 13:08:09,564 Sleeping for 6s. Reason: PENDING
  2024-02-26 13:08:16,534 Sleeping for 7s. Reason: PENDING
  2024-02-26 13:08:24,518 Sleeping for 5s. Reason: PENDING
  2024-02-26 13:08:30,471 Sleeping for 7s. Reason: PENDING
  2024-02-26 13:08:38,498 Sleeping for 5s. Reason: PENDING
  2024-02-26 13:08:44,459 Sleeping for 6s. Reason: PENDING
  2024-02-26 13:08:51,412 Sleeping for 9s. Reason: PENDING
  2024-02-26 13:09:01,412 Sleeping for 9s. Reason: PENDING
  2024-02-26 13:09:11,370 Sleeping for 8s. Reason: PENDING
  2024-02-26 13:09:20,337 Sleeping for 8s. Reason: PENDING
  2024-02-26 13:09:29,316 Sleeping for 6s. Reason: RUNNING
  2024-02-26 13:09:39,703 Padding length to 124
  2024-02-26 13:09:41,194 alphafold2_ptm_model_1_seed_000 recycle=0 pLDDT=73.9 pTM=0.55
  2024-02-26 13:09:42,664 alphafold2_ptm_model_1_seed_000 recycle=1 pLDDT=73.8 pTM=0.549 tol=3.08
  2024-02-26 13:09:44,110 alphafold2_ptm_model_1_seed_000 recycle=2 pLDDT=73.6 pTM=0.549 tol=1.59
  2024-02-26 13:09:45,593 alphafold2_ptm_model_1_seed_000 recycle=3 pLDDT=74.4 pTM=0.555 tol=1.67
  2024-02-26 13:09:45,593 alphafold2_ptm_model_1_seed_000 took 5.9s (3 recycles)
  2024-02-26 13:09:47,073 alphafold2_ptm_model_2_seed_000 recycle=0 pLDDT=76.7 pTM=0.565
  2024-02-26 13:09:48,523 alphafold2_ptm_model_2_seed_000 recycle=1 pLDDT=77.1 pTM=0.57 tol=0.571
  2024-02-26 13:09:49,977 alphafold2_ptm_model_2_seed_000 recycle=2 pLDDT=76.7 pTM=0.569 tol=0.958
  2024-02-26 13:09:51,421 alphafold2_ptm_model_2_seed_000 recycle=3 pLDDT=76.9 pTM=0.572 tol=0.881
  2024-02-26 13:09:51,421 alphafold2_ptm_model_2_seed_000 took 5.8s (3 recycles)
  2024-02-26 13:09:52,877 alphafold2_ptm_model_3_seed_000 recycle=0 pLDDT=75.6 pTM=0.542
  2024-02-26 13:09:54,315 alphafold2_ptm_model_3_seed_000 recycle=1 pLDDT=75.9 pTM=0.548 tol=1.52
  2024-02-26 13:09:55,763 alphafold2_ptm_model_3_seed_000 recycle=2 pLDDT=75.9 pTM=0.552 tol=1.69
  2024-02-26 13:09:57,218 alphafold2_ptm_model_3_seed_000 recycle=3 pLDDT=75.8 pTM=0.555 tol=0.883
  2024-02-26 13:09:57,219 alphafold2_ptm_model_3_seed_000 took 5.8s (3 recycles)
  2024-02-26 13:09:58,705 alphafold2_ptm_model_4_seed_000 recycle=0 pLDDT=73.9 pTM=0.56
  2024-02-26 13:10:00,177 alphafold2_ptm_model_4_seed_000 recycle=1 pLDDT=75.1 pTM=0.57 tol=2.2
  2024-02-26 13:10:01,620 alphafold2_ptm_model_4_seed_000 recycle=2 pLDDT=75.4 pTM=0.571 tol=1.78
  2024-02-26 13:10:03,076 alphafold2_ptm_model_4_seed_000 recycle=3 pLDDT=75.7 pTM=0.575 tol=2.04
  2024-02-26 13:10:03,077 alphafold2_ptm_model_4_seed_000 took 5.8s (3 recycles)
  2024-02-26 13:10:04,572 alphafold2_ptm_model_5_seed_000 recycle=0 pLDDT=75.2 pTM=0.573
  2024-02-26 13:10:06,026 alphafold2_ptm_model_5_seed_000 recycle=1 pLDDT=76.2 pTM=0.585 tol=2.12
  2024-02-26 13:10:07,498 alphafold2_ptm_model_5_seed_000 recycle=2 pLDDT=76.2 pTM=0.587 tol=1.44
  2024-02-26 13:10:08,958 alphafold2_ptm_model_5_seed_000 recycle=3 pLDDT=76.6 pTM=0.589 tol=1.21
  2024-02-26 13:10:08,959 alphafold2_ptm_model_5_seed_000 took 5.9s (3 recycles)
  2024-02-26 13:10:08,986 reranking models by 'plddt' metric
  2024-02-26 13:10:08,987 rank_001_alphafold2_ptm_model_2_seed_000 pLDDT=76.9 pTM=0.572
  2024-02-26 13:10:08,987 rank_002_alphafold2_ptm_model_5_seed_000 pLDDT=76.6 pTM=0.589
  2024-02-26 13:10:08,987 rank_003_alphafold2_ptm_model_3_seed_000 pLDDT=75.8 pTM=0.555
  2024-02-26 13:10:08,987 rank_004_alphafold2_ptm_model_4_seed_000 pLDDT=75.7 pTM=0.575
  2024-02-26 13:10:08,987 rank_005_alphafold2_ptm_model_1_seed_000 pLDDT=74.4 pTM=0.555
  2024-02-26 13:10:10,274 Query 3/29: test3 (length 162)
  2024-02-26 13:10:11,241 Sleeping for 8s. Reason: PENDING
  2024-02-26 13:10:20,230 Sleeping for 10s. Reason: PENDING
  2024-02-26 13:10:31,195 Sleeping for 5s. Reason: RUNNING
  2024-02-26 13:10:37,194 Sleeping for 6s. Reason: RUNNING
  2024-02-26 13:10:44,153 Sleeping for 9s. Reason: RUNNING
  2024-02-26 13:10:54,142 Sleeping for 10s. Reason: RUNNING
  2024-02-26 13:11:05,109 Sleeping for 8s. Reason: RUNNING
  2024-02-26 13:11:14,082 Sleeping for 6s. Reason: RUNNING
  2024-02-26 13:11:21,030 Sleeping for 8s. Reason: RUNNING
  2024-02-26 13:11:30,005 Sleeping for 9s. Reason: RUNNING
  2024-02-26 13:11:39,984 Sleeping for 7s. Reason: RUNNING
  2024-02-26 13:11:47,941 Sleeping for 10s. Reason: RUNNING
  2024-02-26 13:11:58,903 Sleeping for 9s. Reason: RUNNING
  2024-02-26 13:12:08,881 Sleeping for 5s. Reason: RUNNING
  2024-02-26 13:12:14,891 Sleeping for 9s. Reason: RUNNING
  2024-02-26 13:12:32,470 Padding length to 172
  2024-02-26 13:13:00,100 alphafold2_ptm_model_1_seed_000 recycle=0 pLDDT=62.9 pTM=0.433
  2024-02-26 13:13:02,186 alphafold2_ptm_model_1_seed_000 recycle=1 pLDDT=63.4 pTM=0.433 tol=8.27
  2024-02-26 13:13:04,282 alphafold2_ptm_model_1_seed_000 recycle=2 pLDDT=64.1 pTM=0.431 tol=8.02
  2024-02-26 13:13:06,403 alphafold2_ptm_model_1_seed_000 recycle=3 pLDDT=63.8 pTM=0.427 tol=8.51
  2024-02-26 13:13:06,404 alphafold2_ptm_model_1_seed_000 took 33.9s (3 recycles)
  2024-02-26 13:13:08,535 alphafold2_ptm_model_2_seed_000 recycle=0 pLDDT=60.2 pTM=0.417
  2024-02-26 13:13:10,637 alphafold2_ptm_model_2_seed_000 recycle=1 pLDDT=61 pTM=0.423 tol=6.09
  2024-02-26 13:13:12,742 alphafold2_ptm_model_2_seed_000 recycle=2 pLDDT=61.4 pTM=0.428 tol=3.33
  2024-02-26 13:13:14,846 alphafold2_ptm_model_2_seed_000 recycle=3 pLDDT=61.2 pTM=0.425 tol=1.8
  2024-02-26 13:13:14,846 alphafold2_ptm_model_2_seed_000 took 8.4s (3 recycles)
  2024-02-26 13:13:16,979 alphafold2_ptm_model_3_seed_000 recycle=0 pLDDT=62 pTM=0.425
  2024-02-26 13:13:19,099 alphafold2_ptm_model_3_seed_000 recycle=1 pLDDT=62.3 pTM=0.43 tol=7.21
  2024-02-26 13:13:21,197 alphafold2_ptm_model_3_seed_000 recycle=2 pLDDT=61.9 pTM=0.426 tol=4.32
  2024-02-26 13:13:23,303 alphafold2_ptm_model_3_seed_000 recycle=3 pLDDT=62.1 pTM=0.427 tol=5.17
  2024-02-26 13:13:23,304 alphafold2_ptm_model_3_seed_000 took 8.4s (3 recycles)
  2024-02-26 13:13:25,461 alphafold2_ptm_model_4_seed_000 recycle=0 pLDDT=60.5 pTM=0.418
  2024-02-26 13:13:27,552 alphafold2_ptm_model_4_seed_000 recycle=1 pLDDT=60.8 pTM=0.417 tol=9.52
  2024-02-26 13:13:29,658 alphafold2_ptm_model_4_seed_000 recycle=2 pLDDT=60.3 pTM=0.41 tol=9.23
  2024-02-26 13:13:31,749 alphafold2_ptm_model_4_seed_000 recycle=3 pLDDT=60.5 pTM=0.411 tol=6.08
  2024-02-26 13:13:31,750 alphafold2_ptm_model_4_seed_000 took 8.4s (3 recycles)
  2024-02-26 13:13:33,905 alphafold2_ptm_model_5_seed_000 recycle=0 pLDDT=59.9 pTM=0.416
  2024-02-26 13:13:36,038 alphafold2_ptm_model_5_seed_000 recycle=1 pLDDT=60.1 pTM=0.415 tol=9.96
  2024-02-26 13:13:38,154 alphafold2_ptm_model_5_seed_000 recycle=2 pLDDT=59.7 pTM=0.409 tol=3.89
  2024-02-26 13:13:40,252 alphafold2_ptm_model_5_seed_000 recycle=3 pLDDT=59.4 pTM=0.415 tol=11.4
  2024-02-26 13:13:40,253 alphafold2_ptm_model_5_seed_000 took 8.5s (3 recycles)
  2024-02-26 13:13:40,294 reranking models by 'plddt' metric
  2024-02-26 13:13:40,294 rank_001_alphafold2_ptm_model_1_seed_000 pLDDT=63.8 pTM=0.427
  2024-02-26 13:13:40,294 rank_002_alphafold2_ptm_model_3_seed_000 pLDDT=62.1 pTM=0.427
  2024-02-26 13:13:40,294 rank_003_alphafold2_ptm_model_2_seed_000 pLDDT=61.2 pTM=0.425
  2024-02-26 13:13:40,295 rank_004_alphafold2_ptm_model_4_seed_000 pLDDT=60.5 pTM=0.411
  2024-02-26 13:13:40,295 rank_005_alphafold2_ptm_model_5_seed_000 pLDDT=59.4 pTM=0.415
  
• Como ves el propio script espera cuando el servidor remoto está ocupado.
• Para cada secuencia problema obtienes figuras como éstas:
• 
  
  
  
  
  
  
  
  
  
  

Hasta pronto,

Bruno

PD Cuando acabes de instalar deberías tener algo similar en tu fichero $HOME/.bashrc:

export PATH=/usr/local/cuda/bin:$PATH

export LD_LIBRARY_PATH=/usr/local/cuda/lib64:$LD_LIBRARY_PATH

export PATH="$HOME/colabfold/colabfold-conda/bin:$PATH"

Cuando no vayas a usar colabfold comenta estas líneas para usar perl y python del sistema

PD2 Me comentan colegas de ULiverpool que haciendo 800-900 MSAs al día en https://api.colabfold.com no han tenido problemas 

PD3 Si no reconoce tu GPU mira posibles soluciones en https://github.com/YoshitakaMo/localcolabfold/issues/210

PD4: Puedes bloquear la versión de CUDA que hayas instalado con algo como:

 sudo apt-mark hold cuda-toolkit-11-8

Hito 200K del Protein Data Bank

Hola,

a punto de empezar el congreso PAG30, empezamos el año con una buena noticia: 

 

Pego aquí la reseña oficial del Protein Data Bank: 

Date: Wed, 11 Jan 2023 09:55:41 -0500
Subject: pdb-l: PDB Reaches a New Milestone: 200,000+ Entries

With this week's update, the PDB archive contains a record 200,069 
entries. The archive passed 150,000 structures in 2019 and 
100,000 structures in 2014. 

Established in 1971, this central, public archive has reached this 
critical milestone thanks to the efforts of structural biologists 
throughout the world who contribute their experimentally-determined 
protein and nucleic acid structure data.

wwPDB data centers support online access to three-dimensional structures 
of biological macromolecules that help researchers understand many 
facets of biomedicine, agriculture, and ecology, from protein synthesis 
to health and disease to biological energy. Many milestones have been 
reached since the archive released the 100,000th structure in 2014. PDB 
data have been seminal in understanding SARS-CoV-2, and provided the 
foundation for the development of AI/ML techniques for predicting 
protein structure. The 50th anniversary of the PDB was celebrated 
throughout 2021 <https://www.wwpdb.org/pdb50>.

Today, the archive is quite large, containing more than 3,000,000 files 
related to these PDB entries that require more than 1086 Gbytes of 
storage. PDB structures contain more than 1.8 billion non-hydrogen atoms.


Function follows form
In the 1950s, scientists had their first direct look at the structures 
of proteins and DNA at the atomic level. Determination of these early 
three-dimensional structures by X-ray crystallography ushered in a new 
era in biology-one driven by the intimate link between form and 
biological function. As the value of archiving and sharing these data 
were quickly recognized by the scientific community, the Protein Data 
Bank (PDB) was established as the first open access digital resource in 
all of biology by an international collaboration in 1971 with data 
centers located in the US and the UK.

Among the first structures deposited in the PDB were those of myoglobin 
and hemoglobin, two oxygen-binding molecules whose structures were 
elucidated by Chemistry Nobel Laureates John Kendrew and Max Perutz. 
With this week's regular update, the PDB welcomes 266 new structures 
into the archive. These structures join others vital to drug discovery, 
bioinformatics and education.

The PDB is growing rapidly, increasing in size ~13% since 2011. In 2022, 
an average of 275 new structures were released to the scientific 
community each week. The resource is accessed hundreds of millions of 
times annually by researchers, students, and educators intent on 
exploring how different proteins are related to one another, to clarify 
fundamental biological mechanisms and discover new medicines.

Twenty Years of Collaboration
Since its inception, the PDB has been a community-driven enterprise, 
evolving into a mission critical international resource for biological 
research. The wwPDB partnership was established in July 2003 with PDBe, 
PDBj, and RCSB PDB. Today, the collaboration includes partners BMRB 
(joined in 2006) and EMDB (2021).

The wwPDB ensures that these valuable PDB data are securely stored, 
expertly managed, and made freely available for the benefit of 
scientists and educators around the globe. wwPDB data centers work 
closely with community experts to define deposition and annotation 
policies, resolve data representation issues, and implement community 
validation standards. In addition, the wwPDB works to raise the profile 
of structural biology with increasingly broad audiences.

Each structure submitted to the archive is carefully curated by wwPDB 
staff before release. New depositions are checked and enhanced with 
value-added annotations and linked with other important biological data 
to ensure that PDB structures are discoverable and interpretable by 
users with a wide range of backgrounds and interests.

wwPDB eagerly awaits the next 100,000 structures and the invaluable 
knowledge these new data will bring.

Hasta pronto,

Bruno

Alphafold2: estructuras de proteínas para todos los públicos

Hola,

el jueves pasado escribía aquí sobre RoseTTAFold, que ponía al alcance de los programadores una potente herramienta para predecir la estructura de cualquier proteína que tenga suficientes secuencias parecidas en los archivos. Unas horas más tarde me enteré de que su principal competidor, AlphaFold2, del que hemos escrito aquí, aquí y aquí , había hecho también sus deberes y publicaba sus modelos neuronales para predecir el plegamiento de secuencias de aminoácidos. Hay muchos hilos en Twitter, como éstos: 1 y 2 .

En un artículo en Nature repasan sus resultados en CASP14, como su notable precisión de 1.5 Å de RMSD para todos los átomos modelados, y muestran que esos resultados excelentes se sostienen en gran medida en un conjunto de 3144 estructuras recientemente publicadas en el Protein Data Bank tras descartar las estructuras conocidas con % identidad de secuencia > 40. El siguiente vídeo resume el proceso de predicción de una secuencia con múltiples dominios y 863 aminoácidos. Recuerdo que la única información de partida es la secuencia de aminoácidos y un alineamiento múltiple de esa secuencia con otras secuencias no redundantes de proteínas homólogas:

 

Fuente: https://www.nature.com/articles/s41586-021-03819-2

Alphafold2 no es perfecto, ya que en la validación contra el PDB hay un en torno a un 25% de predicciones con valores de RMSD elevados, pero en la mayor parte de los casos es excelente:

Fuente: https://www.nature.com/articles/s41586-021-03819-2

Una de las fortalezas del método es que es capaz de estimar bien la calidad los modelos que genera (ver Figura 2b de https://www.nature.com/articles/s41586-021-03828-1). Entre las limitaciones reseñadas por los autores del trabajo, destacan:

1) La calidad de las predicciones cae de manera significativa si el alineamiento múltiple de secuencias homólogas a la de interés tiene una profundidad < 30 (leer más aquí). 

2) Las estructuas de Alphafold son de menor calidad para dominios poco compactos, con pocos contactos, como los dominios puente. Esto contrasta con las buenas prestaciones observadas (3) para predecir estructuras cuaternarias de homómeros.

Para los que queráis probarlo hay varias opciones:

i) El contenedor Docker descrito en https://github.com/deepmind/alphafold que ojo, requiere 2.2TB de espacio si instalas todas las bases de datos. 

ii) Un cuaderno Colab con un predictor simplificado en https://colab.research.google.com/github/deepmind/alphafold/blob/main/notebooks/AlphaFold.ipynb

iii) Las predicciones ya disponibles en UniProt para un total de 21 especies (https://alphafold.ebi.ac.uk/download), incluyendo 4 plantas (Arabidopsis thaliana, soja, arroz y maíz). Puedes leer más en 4, 5 y 6 y 7.

Hasta pronto,

Bruno

progreso en la predicción estructural de proteínas

Hola,
hace unos meses contaba aquí el algoritmo AlphaFold para plegar proteínas por predicción de distancias entre residuos, que había escuchado de boca de uno de sus creadores. Hoy me he encontrado con la evaluación oficial de estructura terciaria del experimento CASP13, donde AlphaFold se destacó como mejor grupo predictor. La conclusión se resume en esta figura:

Fuente: https://onlinelibrary.wiley.com/doi/10.1002/prot.25787

Parece seguro decir que en CASP13 ha habido un salto en la calidad de las predicciones respecto a ediciones previas, a pesar de que la dificultad en esta edición es comparable a la anterior (Tabla 1 del artículo de Abriata et al de la figura). Los evaluadores achacanla mejoría precisamente a que más allá de predecir contactos, algunos de los mejores predictores, como A7D (AlphaFold), MULTICOM o RaptorX han empezado a predecir directamente distancias entre residuos, algo para la cual hacen falta alineamientos múltiples de secuencia muy profundos. Hasta luego,
Bruno

modelado comparativo de proteínas multidominio

Hola,
en muchas ocasiones el modelado por homología o comparativo es  la única manera que tenemos trabajar con la estructura de una proteína que todavía no está en el Protein Data Bank. De hecho muchos artículos han sido publicados con figuras construidas sobre este tipo de modelos porque ayudan a comprender y poner en contexto tridimensional los resultados.

Interfaz entre dos monómeros modelada por homología, tomada de https://science.sciencemag.org/content/364/6445/1095.

Sin embargo, casi todas las herramientas que existen para modelar proteínas se han centrado históricamente en modelar dominios de proteína uno a uno, cuando la realidad es que muchas proteínas contienen varios dominios. Precisamente para modelar las conformaciones de este tipo de proteínas ha sido publicado recientemente https://zhanglab.ccmb.med.umich.edu/DEMO.

Diagrama de flujo de DEMO, tomado de https://www.pnas.org/content/116/32/15930.

Con la ayuda de DEMO podrás ensamblar dominios previamente modelados de dos en dos. El algoritmo consulta una colección no redundante de estructuras multidominio y optimiza las orientaciones entre dominios, además de que puede usar datos experimentales (cross-linking y crioEM) para guiar el proceso.

Un saludo, Bruno

Introducción a AlphaFold

Hola,
vengo de escuchar a Andrew W Senior en el auditorio Kendrew del EMBL-EBI hablar de cómo funciona AlphaFold para predecir la estructura de proteínas. Como titulaba en Diciembre, todavía no sabemos plegar proteínas, pero cada vez tenemos más información en el Protein Data Bank (del orden de 150K estructuras) para poder predecir distancias entre residuos y ángulos diedros.

Fuente: https://deepmind.com/blog/alphafold

Trataré de resumir lo que nos ha contado antes de que salga el artículo, pero ojo, ha advertido que el código no piensan publicarlo por ahora. Su sistema comprende varios módulos que se ejecutan secuencialmente:

1) Extración de correlaciones evolutivas entre residuos de una secuencia de aminoácidos en forma de modelo de Potts (CCMpred) a partir de perfiles de secuencias homólogas obtenidas con PSIBLAST o HHblits, como hace HHpred.

2) Red neuronal profunda, con al menos dos variantes:
2.1) predicción de distancias reales entre C-betas, no contactos, a partir de histogramas precalculados en el rango de 2 a 22 Ansgtrom. En esto se parecen a RaptorX.
2.2) predicción de ángulos diedros phi y psi

NOTA1: Las predicciones de distancias les permiten asignar estructura secundaria con una precisión Q3  del 84% usando las estructuras de CASP11.

NOTA2: Parten la matriz de contactos en submatrices (crops) que comprenden siempre un trozo de la diagonal principal y esa manera, y con la estructura secundaria, propagan localmente las restricciones de distancias entre residuos cercanos.

3) Diferenciación del potencial de distancias/ángulos por métodos de minimización de gradientes. Lo hacen partiendo la secuencia en dominios o con la secuencia entre y observan que la segunda manera funciona mejor.

4) Relajan el esqueleto obtenido y le añaden cadenas laterales completas con Rosetta. En sus manos esto no mejora el modelo de manera significativa.

Hasta luego,
Bruno

estructuras del PDB parecidas en secuencia (REST)

Hola,
hace unos días un usuario preguntaba en la lista de usuarios del Protein Data Bank (PDB) cómo usar la interfaz de servicios REST, documentada en https://www.rcsb.org/pdb/software/rest.do, para hacer búsquedas BLAST.

Mientras otro usuario compartía un script escrito en Python3, llamado sequenceSimilarity.py, que requiere mmtf-pyspark para hacer consultas PSI-BLAST en tiempo real contra el PDB, a mi me llamó la atención la simplicidad del servicio sequenceCluster, que para cualquier cadena de una estructura depositada en el PDB permite obtener el clúster de secuencias con un porcentaje de identidad controlado por el usuario:

https://www.rcsb.org/pdb/rest/sequenceCluster?cluster=70&structureId=9ant.A

Esto devuelve una lista de cadenas de estructuras PDB en formato XML que podemos procesar por ejemplo como se explica en

https://developer.atlassian.com/server/fisheye-crucible/writing-a-rest-client-in-perl

Hasta luego, buen verano,
Bruno

growth of protein-DNA complexes in the Protein Data Bank

Hi,
while checking the update logs of our good old 3D-footprint, a database of DNA-binding protein structures updated weekly from the Protein Data Bank, I found a folder with logs starting Februrary, 2009. The plot below shows how the number of non-redundant complexes, filtered in terms of protein sequence identity, has doubled in just a decade:

The nr95 bundle can be downloaded in PDB format at
http://maya.ccg.unam.mx/tfmodeller/get_library.cgi

Other related files are available at:
http://floresta.eead.csic.es/3dfootprint/download.html

cheers,
Bruno

formato PDBx/mmCIF del Protein Data Bank

Hola,
aunque el formato PDB sea probablemente todavía el estándar de facto para intercambiar información de macromoléculas biológicas, se ha quedado pequeño para representar grandes complejos moleculares y por esa razón desde 2014 el formato oficial del Protein Data Bank es el PDBx/mmCIF. Hay además otro formato en liza, el PDBML, pero tiene la desventaja de ocupar mucho más espacio.

Complejo proteína-ADN 1LFU, tomado de http://www.rcsb.org/pdb/explore/explore.do?structureId=1lfu

Repasemos cómo pasan a ser los campos ATOM de un fichero PDB, los más relevantes. Veamos el siguiente ejemplo en el vetusto formato PDB, tomado de la estructura 1LFU, con columnas de ancho fijo:

ATOM      1  O5'  DA C 100      31.258  -2.296  76.212  1.00 81.62           O  
ATOM      2  C5'  DA C 100      29.867  -2.121  76.367  1.00 69.89           C  
ATOM      3  C4'  DA C 100      28.980  -3.049  77.172  1.00 67.21           C  
ATOM      4  O4'  DA C 100      29.376  -3.145  78.557  1.00 64.58           O  
ATOM      5  C3'  DA C 100      27.626  -2.376  77.196  1.00 64.41           C  
ATOM      6  O3'  DA C 100      26.569  -3.309  77.165  1.00 66.18           O  
ATOM      7  C2'  DA C 100      27.647  -1.527  78.451  1.00 63.85           C  
ATOM      8  C1'  DA C 100      28.739  -2.123  79.322  1.00 56.01           C  
ATOM      9  N9   DA C 100      29.771  -1.142  79.635  1.00 49.13           N  
ATOM     10  C8   DA C 100      30.533  -0.428  78.740  1.00 48.58           C  
ATOM     11  N7   DA C 100      31.429   0.348  79.306  1.00 43.14           N  
ATOM     12  C5   DA C 100      31.218   0.141  80.664  1.00 40.35           C  
ATOM     13  C6   DA C 100      31.837   0.679  81.794  1.00 42.42           C  
ATOM     14  N6   DA C 100      32.826   1.571  81.750  1.00 48.24           N  
ATOM     15  N1   DA C 100      31.393   0.262  82.998  1.00 42.81           N  
ATOM     16  C2   DA C 100      30.397  -0.626  83.046  1.00 45.83           C  
ATOM     17  N3   DA C 100      29.734  -1.195  82.054  1.00 41.14           N  
ATOM     18  C4   DA C 100      30.197  -0.765  80.875  1.00 41.62           C  
ATOM     19  P    DG C 101      25.116  -2.785  76.764  1.00 78.33           P  
ATOM     20  OP1  DG C 101      24.278  -3.977  76.501  1.00 74.57           O

En formato PDBx/mmCIF esto mismo se expresa así, con columnas separadas con espacios en blanco y cabeceras auto-explicativas, incluyendo una para la carga:

loop_
_atom_site.group_PDB 
_atom_site.id 
_atom_site.type_symbol 
_atom_site.label_atom_id 
_atom_site.label_alt_id 
_atom_site.label_comp_id 
_atom_site.label_asym_id 
_atom_site.label_entity_id 
_atom_site.label_seq_id 
_atom_site.pdbx_PDB_ins_code 
_atom_site.Cartn_x 
_atom_site.Cartn_y 
_atom_site.Cartn_z 
_atom_site.occupancy 
_atom_site.B_iso_or_equiv 
_atom_site.pdbx_formal_charge 
_atom_site.auth_seq_id 
_atom_site.auth_comp_id 
_atom_site.auth_asym_id 
_atom_site.auth_atom_id 
_atom_site.pdbx_PDB_model_num 
ATOM   1    O  "O5'" . DA  A 1 1  ? 31.258 -2.296 76.212  1.00 81.62  ? 100 DA  C "O5'" 1 
ATOM   2    C  "C5'" . DA  A 1 1  ? 29.867 -2.121 76.367  1.00 69.89  ? 100 DA  C "C5'" 1 
ATOM   3    C  "C4'" . DA  A 1 1  ? 28.980 -3.049 77.172  1.00 67.21  ? 100 DA  C "C4'" 1 
ATOM   4    O  "O4'" . DA  A 1 1  ? 29.376 -3.145 78.557  1.00 64.58  ? 100 DA  C "O4'" 1 
ATOM   5    C  "C3'" . DA  A 1 1  ? 27.626 -2.376 77.196  1.00 64.41  ? 100 DA  C "C3'" 1 
ATOM   6    O  "O3'" . DA  A 1 1  ? 26.569 -3.309 77.165  1.00 66.18  ? 100 DA  C "O3'" 1 
ATOM   7    C  "C2'" . DA  A 1 1  ? 27.647 -1.527 78.451  1.00 63.85  ? 100 DA  C "C2'" 1 
ATOM   8    C  "C1'" . DA  A 1 1  ? 28.739 -2.123 79.322  1.00 56.01  ? 100 DA  C "C1'" 1 
ATOM   9    N  N9    . DA  A 1 1  ? 29.771 -1.142 79.635  1.00 49.13  ? 100 DA  C N9    1 
ATOM   10   C  C8    . DA  A 1 1  ? 30.533 -0.428 78.740  1.00 48.58  ? 100 DA  C C8    1 
ATOM   11   N  N7    . DA  A 1 1  ? 31.429 0.348  79.306  1.00 43.14  ? 100 DA  C N7    1 
ATOM   12   C  C5    . DA  A 1 1  ? 31.218 0.141  80.664  1.00 40.35  ? 100 DA  C C5    1 
ATOM   13   C  C6    . DA  A 1 1  ? 31.837 0.679  81.794  1.00 42.42  ? 100 DA  C C6    1 
ATOM   14   N  N6    . DA  A 1 1  ? 32.826 1.571  81.750  1.00 48.24  ? 100 DA  C N6    1 
ATOM   15   N  N1    . DA  A 1 1  ? 31.393 0.262  82.998  1.00 42.81  ? 100 DA  C N1    1 
ATOM   16   C  C2    . DA  A 1 1  ? 30.397 -0.626 83.046  1.00 45.83  ? 100 DA  C C2    1 
ATOM   17   N  N3    . DA  A 1 1  ? 29.734 -1.195 82.054  1.00 41.14  ? 100 DA  C N3    1 
ATOM   18   C  C4    . DA  A 1 1  ? 30.197 -0.765 80.875  1.00 41.62  ? 100 DA  C C4    1 
ATOM   19   P  P     . DG  A 1 2  ? 25.116 -2.785 76.764  1.00 78.33  ? 101 DG  C P     1 
ATOM   20   O  OP1   . DG  A 1 2  ? 24.278 -3.977 76.501  1.00 74.57  ? 101 DG  C OP1   1

La diferencia más notable, además de que las columnas no son de ancho fijo, es que cada átomo tiene ahora nuevas etiquetas (label) identificativas, además de las asignadas originalmente por los autores. En el ejemplo, la desoxiadenina 100 para a llevar el número 1. Otra diferencia notable es que en estructuras de NMR, como este caso, el número de modelo se indica también para cada átomo (entity_id).

Hasta luego,
Bruno

trabajo en PDB Europa

The Protein Data Bank in Europe (PDBe) is seeking to recruit a 
structural analyst/programmer to join the team at the European 
Bioinformatics Institute located on the Wellcome Trust Genome Campus 
near Cambridge in the UK.

For further information and to apply, please see the full EMBL-EBI job 
listing at http://bit.ly/1jvNjIO
To discuss this post informally, please contact PDBe's Dr Sameer 
Velankar (sameer@ebi.ac.uk).

About the position:

Established in 1994, CATH and SCOP are the world?s most comprehensive 
resources classifying protein-domain structures into evolutionary 
superfamilies. They are currently being combined in a collaborative 
project - Genome3D - that aims to provide predicted 3D structures for 
sequences assigned to SCOP and CATH superfamilies. Combining SCOP and 
CATH based predictions allows us to identify more accurately regions 
that agree between the two methods. We aim to provide these 
Genome3D-predicted structures via the PDBe resource. To improve the 
assessment of the reliability of the predictions it is necessary to 
develop a mapping between SCOP and CATH and to remove any conflicts.

We are now seeking to recruit a structural analyst/programmer to assist 
with this task. The structural analyst/programmer will also build an 
automatic pipeline to generate putative domain assignments (from CATH) 
for new PDB structures, prior to classification in SCOP or CATH. The 
structural analyst/programmer will be based at PDBe, and will be jointly 
supervised by Prof. Christine Orengo (UCL, London) and Dr Alexey Muzin 
(MRC-LMB, Cambridge), together with Prof. Gerard Kleywegt and Dr Sameer 
Velankar at PDBe.

About PDBe:

The Protein Data Bank in Europe (PDBe; pdbe.org) is part of the 
Worldwide Protein Data Bank organisation (wwPDB; wwpdb.org), which 
maintains the global archive of 3D structural data on biomacromolecules. 
The PDBe team also maintains a number of databases that support 
deposition and advanced search services for structural biologists and 
the wider scientific community. The team consists of an international 
and inter-disciplinary mix of professionals (scientists and IT specialists).

Best regards,
Gary.


-- 
Gary Battle
Protein Data Bank in Europe (PDBe)
EMBL-EBI
Wellcome Trust Genome Campus
Hinxton, Cambridge
CB10 1SD

http://www.facebook.com/proteindatabank
http://twitter.com/PDBeurope

TFcompare - a tool for structural alignment of DNA binding protein complexes

I want to introduce you the new bioinformatic contribution of our lab to the science world: TFcompare (http://floresta.eead.csic.es/tfcompare/)

TFcompare is a tool for structural alignment of DNA motifs and protein domains from DNA binding protein complexes in Protein Data Bank. 

The TFcompare algorithm calculates structural alignments between three dimensional structures of two DNA-protein complexes. The most interesting feature of TFcompare when compared with other methods is that it extracts individual protein domains and their recognized DNA sequences, aligning them separately and returning not only the structure superposition but the DNA sequence superposition too. In this way we can compare single domain affinity for different DNA sequences in DNA-protein complexes, especially transcription factors and their recognized cis elements.

The working schema of TFcompare is the following:

TFcompare takes as input two PDB identifiers. Structures from PDB are retrieved automatically and Pfam domains contacting DNA are calculated and trimmed from the original structure. Then all the domains from the first structure are aligned to all the domains from the second in several steps:

1. The program MAMMOTH performs the structural alignment.
2. The produced transformation matrices are applied to the coordinates of the DNA binding sites in order to derive the equivalent cis element superpositions.
3. Root-mean-squared deviations of superposed coordinates are calculated with beta-carbon atoms (proteins) and with N9 (purines) and N1 (pyrimidines) atoms (DNA).
4. Structural alignments are scored in terms of i) the number of identical superposed nucleotides (DNA Score 1-0) and ii) the sum of N9 and N1 atom pairs within 3.5 Å (DNA Score).

We can take as example the alignment of 1D5Y and 1BL0 structures, both are bacterial proteins with Helix-turn-helix (HTH) protein domains binding DNA.

We obtain the following results:

Results are ordered by structural similarity (RMSD), from both protein domain and DNA. In green colour are showed the alignment of similar structures (protein RMSD <=5.0 Å and DNA RMSD <= 3.5 Å) and in red colour dissimilar ones. 1D5Y contains two protein chains with HTH domains contacting DNA (trimmed domain structures are 1d5y_A1 and 1d5y_C1). 1BL0 have two HTH domains in its unique chain A (1bl0_A1 and 1bl0_A2). Structural alignment results show how 1bl0_A1 superposes very well with 1d5y_A1 and 1d5y_C1 (green colour). When DNA sequences recognized by these domains are aligned, they show a DNA motif conserved with three common nucleotides ‘CAC’. However, 1bl0_A2 superpositions (red colour) are not as good as previous ones, and DNA motif ‘CAC’ is not preserved when checking the resulting DNA alignment.

Each row contains an alignment of a pair of DNA binding domains, showing a picture of their structures before and after superposition. DNA alignment is also shown.

PDB files with aligned structures can be downloaded by left-clicking on the domain names and DNA sequences. Opening them with PDB viewer software (Pymol for ex.) is possible to visualize the resulting superposition after structural alignment.

Results column headers and their meaning:
Pair: Pair number
Domain_Query: PDB name, chain and domain number of the Query
Domain_Sbjct: PDB name, chain and domain number of the Sbjct
DNA_Query: DNA site recognized by the Query domain
DNA_Sbjct: DNA site recognized by the Query domain
Similar: 1 if both protein domains and DNA sites are below RMSD thresholds, 5.0 A and 3.5 A
respectively
DNA_Alignment: DNA sites structurally aligned
DNA_Aligned: Number of aligned nucleotides
DNA_Score_1-0: Number of identical nucleotides
DNA_Score: Structural alignment score
DNA_RMSD: RMSD of the structurally aligned DNA sites
PROT_RMSD: RMSD of the structurally aligned protein domains
3D_Alignment: 3D Visualization of aligned structures

Felices Navidades bioinformáticas!!!

Ya que estamos en Navidad os dejo nuestra felicitación navideña especialmente pensada para los bioinformáticos... con algo de inspiración en estructura tridimensional de proteínas y dianas terapéuticas.

Adjunto la imagen a alta resolución para que podáis intentar adivinar con qué estructura del PDB se han generado las bolas (proteínas) y adornos (ligandos) que decoran el árbol, dejar en los comentarios vuestras respuestas...

Y por supuesto... FELIZ NAVIDAD Y PRÓSPERO AÑO 2012!!!

English version:

In this Christmas time we want to share our Christmas card especially designed for bioinformaticians ... with some inspiration from three-dimensional structure of proteins and therapeutic targets.

It is a high resolution image so you can try to guess from which PDB structure have been extracted the balls (proteins) and ornaments (ligands) that decorate the tree, give your answers in the comment section ...

And of course... MERRY CHRISTMAS AND HAPPY NEW YEAR 2012!!!

Sincronizando el Protein Data Bank

Hola,
para los que trabajamos de manera habitual con archivos del Protein Data Bank es muy conveniente tener acceso a la colección de estructuras en formato PDB. Ojo, necesitarás al menos unos 14Gb de disco, pero tendrás al  alcance de tu mano más de 77 mil estructuras de proteínas (y sus ligandos). La sección de descargas del PDB ofrece varias opciones, pero en mi opinión, si necesitas acceder desde tus programas a los archivos PDB sin supervisión, lo más eficiente es tener una copia local, lo que en la jerga se llama mirror o sitio espejo. Como el PDB se va actualizando semanalmente y no sólo se añaden entradas nuevas, si no que se van corrigiendo o completando estructuras antiguas, el proceso de sincronización no es tan sencillo como pudiera parecer. Es aquí donde nos salva la vida el software rsync, parte integrante de cualquier Linux moderno, que permite hacer estas actualizaciones de manera eficiente, y que fácilmente podemos añadir a un cron semanal:

 #!/usr/bin/perl -w   
 use strict;  
   
 my $PDBSERVER = 'rsync.wwpdb.org';  
 my $PDBPORT  = 33444;  
 my $LOCALPDBDIR = '/path/to/my/PDB/mirror/';  
 my $LOGFILE   = '/path/to/PDBmirror.log';  
   
 my $rsync_command = "/usr/bin/rsync -rlpt -v -z --delete --exclude=biounit ".  
      "--exclude=monomers --exclude=status --exclude=obsolete --exclude=models ".  
      "--exclude=mmCIF --exclude=nmr_restraints* --exclude=structure_factors " .  
   "--exclude=XML --exclude=XML-extatom --exclude=XML-noatom ".  
      "--port=$PDBPORT $PDBSERVER\::ftp_data/ $LOCALPDBDIR";  
   
 open(LOG,">$LOGFILE") || die "# $0 : cannot write to $LOGFILE\n";  
 open(RSYNC,"$rsync_command 2>&1 |") || die "# $0 : cannot execute $rsync_command\n";  
 while(<RSYNC>)  
 {  
   print LOG;   
   if(/rsync: failed to connect to/)  
   {  
       die "# $0 : error:\n$_\n";  
   }  
 }  
 close(RSYNC);  
 close(LOG);  

Un saludo, Bruno

Extraer coordenadas de átomos en un fichero PDB

Los ficheros PDB (Protein Data Bank) contienen las coordenadas espaciales de los átomos de proteínas cuya estructura está resuelta por técnicas de rayos X o resonancia magnética nuclear (RMN). En una entrada de blog de Bruno podéis encontrar una muy buena explicación de estas técnicas realizada por el periódico El País.

Las estructuras de proteínas en formato PDB contienen mucha información que habitualmente no nos interesa para nuestros experimentos e interfiere con muchos programas que únicamente necesitan las coordenadas espaciales de los átomos. Además un fichero PDB normalmente contiene datos de múltiples moléculas y en diferentes posiciones en el cristal. Por ello es muy conveniente extraer de los ficheros PDB únicamente la información estructural que vamos a utilizar y descartar el resto.

En el post anterior vimos como extraer líneas de un texto o archivo. En la entrada de hoy veremos como reusar ese código para extraer los átomos deseados de un fichero PDB.

En el ejemplo se indica el fichero PDB original en la variable $pdb_file, en nuestro caso será la estructura de la famosa insulina que se puede descargar aquí (2INS.pdb). Y en el array @desired_coords le indicaremos el tipo de átomos que queremos extraer en un formato muy similar al usado por PyMOL.

   
 # File with 3D coords of insulin  
 my $pdb_file = "2INS.pdb";  
   
 # Extract coords of the alpha carbons of the chain A, and alpha and beta carbons of histidines of chain B  
 my @desired_coords = ('a/*/ca/', 'b/*/ca|cb/his');  
   
 open(PDBFILE,$pdb_file);  
 my $pdb_content = join('',<PDBFILE>);  
 close PDBFILE;  
   
 my $pdb_coords = join('',extract_pdb_coords($pdb_content,\@desired_coords));  
   
 open(PDBFILE,">$pdb_file.out");  
 print PDBFILE $pdb_coords;  
 close PDBFILE;  
   
 # Extract desired PDB coordinates from PDB entries  
 sub extract_pdb_coords {  
   
      my ($pdb_content, $types) = @_;  
   
      my $pdb_data;  
      my $patterns;  
      my @pattern_parts = ('chain','res_number','res_type','res_name');  
      foreach my $type (@{$types}) {  
           my $count = 0;  
           my %pattern;  
           while ($type =~ /([^\/]+)/g){  
                $type = $'; # Text to the right of the match  
                $pattern{$pattern_parts[$count]} = $1;  
                $count++;  
           }  
           push(@{$patterns},\%pattern);  
      }  
   
      foreach my $pattern (@{$patterns}){  
           my ($chain,$res_number,$res_type,$res_name);  
           if (!defined($pattern->{'chain'}) || !$pattern->{'chain'} || $pattern->{'chain'} eq '*'){  
                $chain = '\w{1}';  
           } else {  
                $chain = uc($pattern->{'chain'});  
           }  
           $pattern->{'res_number'} =~ s/\s+//;  
           if (!defined($pattern->{'res_number'}) || !$pattern->{'res_number'} || $pattern->{'res_number'} eq '*'){  
                $res_number = '\d+';  
           } elsif ($pattern->{'res_number'} =~ /^(\d+)$/) {  
                $res_number = $1;  
           } elsif ($pattern->{'res_number'} =~ /^(\d+)-(\d+)$/) {  
                $res_number = '('.join('|', $1 .. $2).')';  
           } elsif ($pattern->{'res_number'} =~ /\d,/) {  
                $res_number = '('.join('|', split(",",$pattern->{'res_number'})).')';  
           }  
           if (!defined($pattern->{'res_type'}) || !$pattern->{'res_type'} || $pattern->{'res_type'} eq '*'){  
                $res_type = '[\w\d]+';  
           } else {  
                $res_type = uc($pattern->{'res_type'});  
           }  
           if (!defined($pattern->{'res_name'}) || !$pattern->{'res_name'} || $pattern->{'res_name'} eq '*'){  
                $res_name = '\w{3}';  
           } else {  
                $res_name = uc($pattern->{'res_name'});  
           }  
           $pattern = '/(ATOM|HETATM)\s+\d+\s+('.$res_type.')\s+('.$res_name.')\s('.$chain.')\s+('.$res_number.')\s+.+/';  
      }  
   
      my @pdb_data_lines = extract_lines_from_text($pdb_content, $patterns);  
      if (@pdb_data_lines){  
           $pdb_data = join("\n",@pdb_data_lines);  
      }  
   
      return $pdb_data;  
   
 }  
   
 # Extract lines from text with the desired patterns  
 sub extract_lines_from_text {  
   
      my ($text, $patterns) = @_;  
   
      my @data;  
      my @lines = split("\n",$text);  
   
      foreach my $line (@lines){  
           foreach my $pattern (@{$patterns}){  
                if ($pattern =~ /^\/(.+)\/$/){  
                     if ($line =~ /$1/){  
                          push(@data,$line);  
                          last;  
                     }  
                } else {  
                     if ($line =~ /\Q$pattern\E/){  
                          push(@data,$line);  
                          last;  
                     }  
                }  
           }  
      }  
   
      return @data;  
   
 }  

El Protein Data Bank en El País

Buenos días,
antes de seguir con las estructuras de datos, que dejamos para otro día, hoy me gustaría destacar que recientemente ha salido en el diario El País, uno de los más leídos en español, un artículo donde se resumían los últimos avances en las técnicas de resolución de estructuras moleculares, que poco a poco parece que se van encaminando al estudio de moléculas individuales. El artículo original, muy interesante,  es Destellos brillantes y ultracortos iluminarán la nueva biología estructural.
La verdad me sorprendió gratamente encontrarme en la prensa generalista con información reciente y relevante de acerca del Protein Data Bank (PDB), el recurso sobre el que se construye la Bioinformática Estructural, pero más me sorprendí al descubrir que en realidad el PDB había sido objeto de al menos otros 3 artículos en el mismo diario en años recientes, todos ellos firmados por Cele Abad Zapatero:
2004) La revolución de los rayos X
2007) ¿Morirá de éxito la biología estructural?
2010) Medio siglo de las primeras estructuras de proteínas

Si no habéis usado nunca el PDB, el enlace principal es http://www.rcsb.org,  y una entrada típica, como la 1le8, contiene imágenes como ésta:

Por supuesto podéis acceder a él desde vuestros programas Perl, por ejemplo con el módulo WWW::PDB. En el curso de Algoritmos en Bioinformática Estructural hay varios ejemplos de usos de archivos de coordenadas en formato PDB.

Un saludo y que vaya bien la semana,
Bruno

Cómo insertar una molécula interactiva en una página web con Jmol

Jmol es un visor de estructuras químicas y moléculas en 3D, está programado en Java y su código es libre.

Jmol se puede descargar aquí, es aconsejable descargar la última versión estable y existen 2 opciones, descargar todo el código (sólo para programadores) o descargar sólo los archivos binarios para ejecutarlo.

Una vez descargado, hay que descomprimirlo (si se ha descargado el archivo .tar.gz, ejecutar: tar -xzf archivo.tar.gz). Jmol se inicia ejecutando jmol.bat si usamos Windows y jmol.sh si usamos  Linux. Ni que decir que Jmol necesita que tengamos una versión de Java Runtime instalada en nuestro sistema operativo.

Una vez iniciado Jmol, el programa nos ofrece un interfaz gráfico donde podemos abrir archivos de moléculas fácilmente y visualizarlas. En el siguiente enlace hay diversos tutoriales para manejar Jmol: http://wiki.jmol.org/index.php/Jmol_Tutorials.

Jmol ofrece una interesante característica, nos permite insertar fácilmente estructuras de moléculas en páginas web y visualizarlas de forma interactiva.
Para ello hay que copiar la carpeta de Jmol en nuestro servidor e insertar el siguiente código dentro de nuestra web:

 <html>  
      <head>  
           <script type="text/javascript" src="./jmol/Jmol.js"></script>  
      </head>  
      <body>  
           <script type="text/javascript">  
                jmolInitialize("./jmol", "JmolAppletSigned.jar");  
                jmolApplet(800, "load 1je8_AB.pdb; spacefill off; wireframe off; select all; cartoon; color structure;");  
           </script>  
      </body>  
 </html>  

En primer lugar se indica en la sección <head> de la página web la ruta a la librería JavaScript de Jmol: "<script type="text/javascript" src="./jmol/Jmol.js"></script>", en nuestro caso la librería está en el directorio 'jmol' en la misma ruta que nuestra página web.

En segundo lugar se inicializala librería anterior (dentro de la sección </body> de la página web): "jmolInitialize("./jmol", "JmolAppletSigned.jar");". El segundo parámetro (JmolAppletSigned.jar) sólo es necesario si el servidor web es nuestro propio ordenador.

Finalmente se ejecuta la miniaplicación (applet) indicando el tamaño de la molécula y los comandos de Jmol que queremos ejecutar antes de mostrar la molécula: "jmolApplet(800, "load 1je8_AB.pdb; spacefill off; wireframe off; select all; cartoon; color structure;");". En este caso se ha borrado la visualización previa de la molécula y se ha representado en modo 'cartoon' (utilizado para visualizar fácilmente estructuras secundarias de proteínas).

Así es como se vería en nuestro navegador: