Blast és l'eina principal de l'NCBI per comparar una seqüència d’ADN o proteïna amb altres seqüències de diverses bases de dades

Introducció

Basic Local Alignment Search Tool (BLAST) és l'eina principal de l'NCBI per comparar una seqüència d’ADN o proteïna amb altres seqüències de diverses bases de dades.

La cerca BLAST és una de les maneres fonamentals d'aprendre sobre una proteïna o un gen: la cerca revela quines seqüències relacionades estan presents en el mateix organisme i en altres organismes.

Blast+

Els programes BLAST hostetjats al lloc web de NCBI són extremadament populars.

Com a alternativa, també pots utilitzar BLAST+, un conjunt d'executables de línia d'ordres per a ús local. Aquesta alternativa és molt interessant, ja que així pots configurar bases de dades personalitzades, fer cerques massives (és a dir, utilitzar un nombre molt gran de consultes), implementar estratègies de cerca complexes mitjançant scripts personalitzats o utilitzar el teu clúster d'ordinadors per millorar el rendiment.

En aquest enllaç tens el manual: BLAST Command Line Applications User Manual.

Entorno de trabajo

Per treballar amb blast utilitzarem la imatge ncbi\blast.

Descarrega l'script blast.sh:

$ wget https://gitlab.com/xtec/bio-blast/-/raw/main/blast.sh
$ chmod u+x blast.sh

Executa l'script - el primer cop tardará més temps en executar-se perquè ha de descarregar la imatge:

$ ./blast.sh

...
Status: Downloaded newer image for ncbi/blast:latest
root@8e54f06bb10e:/blast#

L’script també crea unes carpetes compartides entre el contenidor i la teva màquina:

$ ls blast
...

Per sortir del contenidor executa exit:

$ exit

Crea una base de dades personalitzada

BLAST et permet crear bases de dades personalitzades

Per realitzar la consulta utilitzarem una base de dades amb 7 seqüències de proteïnes del tauró nodrissa.

Primer hem de descarregar els fitxers fasta de la base de dades protein mitjançant Entrez:

$ ./blast.sh 
$ efetch -db protein -format fasta -id Q90523,P80049,P83981,P83982,P83983,P83977,P83984,P83985,P27950 > fasta/nurse-shark-proteins.fasta

Pots veure que les seqüències s'han guardat en el fitxer nurse-shack-proteins.fasta:

$ cat fasta/nurse-shark-proteins.fasta | head -n 1

>sp|P27950.1|NDK_GINCI RecName: Full=Nucleoside diphosphate kinase; Short=NDK; Short=NDP kinase

A continuació crea una base de dades BLAST amb el nom nurse-shark-proteins:

$ makeblastdb -in fasta/nurse-shark-proteins.fasta -dbtype prot -parse_seqids -out blastdb_custom/nurse-shark-proteins -title "Nurse shark proteins" -taxid 7801 -blastdb_version 5

Building a new DB, current time: 03/18/2024 17:08:17
New DB name:   /blast/blastdb_custom/nurse-shark-proteins
New DB title:  Nurse shark proteins
Sequence type: Protein
Keep MBits: T
Maximum file size: 3000000000B
Adding sequences from FASTA; added 7 sequences in 0.00194597 seconds.

Com que els fitxers fasta són pocs i petits, la base de dades es construeix en un moment.

Pots veure els fitxers que formen la base de dades:

$ ls blastdb_custom/

nurse-shark-proteins.pdb  nurse-shark-proteins.pin  nurse-shark-proteins.pog  nurse-shark-proteins.pot  nurse-shark-proteins.ptf
nurse-shark-proteins.phr  nurse-shark-proteins.pjs  nurse-shark-proteins.pos  nurse-shark-proteins.psq  nurse-shark-proteins.pto

Pots verificar la base de dades BLAST que acabes de crear executar l'ordre següent per mostrar les accessions, la longitud de la seqüència i el nom comú de les seqüències a la base de dades.

$ blastdbcmd -entry all -db nurse-shark-proteins -outfmt "%a %l %T"

P27950.1 151 7801
P80049.1 132 7801
P83977.1 95 7801
P83981.1 53 7801
P83984.1 190 7801
P83985.1 195 7801
Q90523.1 106 7801

Consultar la base de dades

Per realitzar una consulta, primer hem d’obtenir la seqüència de consulta.

En el nostre exemple serà la proteína P01349 del tauró tigre de sorra que es va seqüenciar el 1986:

$ efetch -db protein -format fasta -id P01349 > queries/P01349.fasta

Mira el contingut del fitxer fasta que acabes de baixar:

$ cat queries/P01349.fasta

>sp|P01349.2|RELX_CARTA RecName: Full=Relaxin; Contains: RecName: Full=Relaxin B chain; Contains: RecName: Full=Relaxin A chain
QLCGRGFIRAIIFACGGSRWATSPAMSIKCCIYGCTKKDISVLC

A continuació pots executar la consulta blastp:

$ blastp -query queries/P01349.fasta -db nurse-shark-proteins

...
>P80049.1 RecName: Full=Fatty acid-binding protein, liver; AltName: Full=Liver-type 
fatty acid-binding protein; Short=L-FABP
Length=132

 Score = 14.2 bits (25),  Expect = 0.96, Method: Compositional matrix adjust.
 Identities = 3/9 (33%), Positives = 6/9 (67%), Gaps = 0/9 (0%)

Query  2    LCGRGFIRA  10
            +C R ++R 
Sbjct  123  VCTREYVRE  131
...

El resultat d’aquest consulta és la identificació de la seqüència de proteïnes P80049.1 com una coincidència amb una puntuació de 14,2 i un valor E de 0,96.

Per un anàlisi posterior pots utilitzar l'opció -out per guardar la sortida en un fitxer.

$ blastp -query queries/P01349.fasta -db nurse-shark-proteins  -out results/blastp.out

Protein Data Bank

Les bases de dades personals són útils quan estàs treballant amb seqüències que no són públiques, però en la majoria dels casos es treballa amb bases de dades públiques.

Primer has de descarregar la base de dades d'aminoàcids Protein Data Bank (pdbaa)

El fitxer update_blastdb.pl` és un script de Perl que proporciona l’NCBI que ja està incorporat i configurat a les variables d’entorn.

L’has d’executar dins de la carpeta blastdb!

Fem servir l'opció -source gcp perquè descarregui les dades desde Google Cloud enlloc del servidor FTP del NCBI:

$ cd blastdb
$ update_blastdb.pl -source gcp pdbaa
Connected to GCP
$ cd ..

Pots veure els fitxers que composen la base de dades s’han descarregat a la carpeta blastdb:

$ ls blastdb

pdbaa.pdb  pdbaa.phr  pdbaa.pin  pdbaa.pog  pdbaa.pos  pdbaa.pot  pdbaa.ppd  pdbaa.ppi  pdbaa.psq  pdbaa.ptf  pdbaa.pto

A continuació podem consultar la seqüència P01349.fasta contra la base de dades PDB tal com es mostra a continuació:

$ blastp -query queries/P01349.fasta -db pdbaa
...
Sequences producing significant alignments:                          (Bits)  Value

2FHW_B Chain B, Relaxin 3 (Prorelaxin H3) (Insulin-like peptide I...  38.1    1e-05
2H8B_B Chain B, Insulin-like 3 [synthetic construct]                  34.7    3e-04
2K6U_B Chain B, Insulin-like 3 B chain [synthetic construct]          30.8    0.010
...

En aquest cas hi ha moltes més coincidències i amb un score més alt.

Bases de datos

Bases de dades

Tal com has vist abans update_blastdb.pl és un script en perl que et permet gestionar les bases de dades.

Executa update_blastdb.pl sense arguments perquè et mostri el menú d’ajuda, i poder veure totes les bases de dades disponibles

$ update_blastdb.pl --showall
18S_fungal_sequences
Betacoronavirus
...

Algoritme

Els algorismes de semblança global optimitzen l'alineament global de dues seqüències, i són els més adequats per trobar coincidències que consisteixen en trams llargs de poca similitud.

En canvi, els algorismes de similitud local com BLAST identifiquen alineaments relativament curts, i són útils per a la cerca en bases de dades perquè moltes seqüències de consulta (“query”) tenen dominis, llocs actius o altres motius que tenen regions locals de semblança amb altres proteïnes ( però no globals),

Normalment, les bases de dades també tenen fragments d'ADN i seqüències de proteïnes que es poden alinear localment a una consulta.

Llista

En la primera fase, l'algorisme BLASTP compila una llista de "paraules" de longitud fixa w que es deriven de la seqüència de consulta.

Per a les cerques de proteïnes, la mida de la paraula normalment té un valor predeterminat de 3 aminoàcids, encara que es pot modificar. Com que hi ha 20 aminoàcids, hi ha 203 = 8000 paraules possibles. Això és així perquè es tracta d’un cas de variacions amb repetició (és important l’ordre de cada aminoàcid i es poden repetir).

Per exemple si la seqüència “query” és la “human beta globin” (NP_000509.1) la seqüència incloura aquests aminoàcdis …VTALWGKVNVD… i algunes de les paraules derivades de la seqüència de consulta seran VTA TAL ALW ``LWG WGK GKV KVN VNV NVD.

Per tant, si una seqüència té 100 aminoàcids tindrem 98 paraules de consulta de mida 3.

A continuació, per cada paraula de consulta es genera una llista de paraules objectiu (“target”) de la mateixa mida amb la seva puntuació de similitud derivada d’una matriu com pot ser BLOSUM62.

Aquesta llista està ordenada pel valor de similitud entre la paraula de consulta i la paraula target.

Per exemple, per a la paraula de consulta LWG es pot generar aquesta llista:

LWG 4+11+6=21	LWA  4+11+0=15
IWG 2+11+6=19	LYG  4+ 2+6=12
MWG 2+11+6=19	--------------
VWG 1+11+6=18	LFG  4+ 1+6=11
FWG 0+11+6=17	FWS  0+11+0=11
AWG 0+11+6=17	AWS -1+11+0=10
LWS 4+11+0=15	CWS -1+11+0=10
LWN 4+11+0=15	IWC  2+11-3=10

No totes les paraules es fan servir, sinó aquelles que estiguin per sobre d’un llindar establert T (“threshold”). En el nostre exemple hem establert el llindar a 12. Si reduïm el valor de T augmentem la sensibilitat de la cerca (tindrem més “hits” perquè tindrem més paraules target) a canvi de que la cerca trigui més temps en executar-se.

Scan y extensiones

A continuació, per cada paraula de consulta es busca a l'índex de la base de dades totes les paraules target associades a la paraula de consulta amb una puntuació superior al llindar T. Si una paraula target és de mida 3, i suposem que tots els aminoàcids tenen la mateixa freqüència d’aparició (que sabem que no és cert), una paraula target tindrà 1 “hit” cada 8000 seqüències de la base de dades. Si una seqüència te un “hit”, es comença a fer una alieament en ambdues direccions tal com es mostra en aquest exemple amb l’alpha globin:

LTPEEKSAVTA[LWG]KV--NVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPDAVMGNPKV HBB
L+P +K+ V A[ WG]KV + E G EAL R+ + +P T+ +F F D G+ +V
LSPADKTNVKA[AWG]KVGAHAGEYGAEALERMFLSFPTTKTYFPHF------DLSHGSAQV HBA

Tal com pots veure a l’exemple, primer l’extensió es fa sense buits i després amb buits si la puntuació( “score”) ha arribat a un llindar, per obtenir parells de segments d'alta puntuació (HSP). Si una puntuació HSP supera una puntuació de tall S particular, aquest alineament local passa a formar part del resultat. Si durant l’extensió la puntuació cau per sota d'un límit l’alineament es descarta la seqüència.

Traceback

Un cop es tenen totes les seqüències amb un puntuació HSP superior al valor de tall S, s’assignen les ubicacions d'insercions, supressions i desajustos de les seqüències seleccionades. També s’apliquen unes estadístiques que descriurem més endavant.

BLASTN

Per a BLASTN, la primera fase és lleugerament diferent perquè l'algoritme exigeix ​​coincidències exactes de paraules (no hi ha matrius de similitut en l’alineament de seqüències d’ADN).

La mida predeterminada de la paraula és 11 – es pot ajustar-la a valors de 7 o 15. Reduir la longitud de la paraula té el mateix efecte que reduir la puntuació llindar en consultes blastp.

Threshold

A continuació es demostrar l'efecte de diferents nivells de llindar en una cerca BLASTP canviant el paràmetre T del seu valor predeterminat (11) a un rang d'altres valors en una consulta de betaglobina humana:

Python

Amb Python pots interactuar directament amb les eines BLAST+ tal com s’explica a Docker.

A continuació tens un exemple:

import os
from python_on_whales import docker

CWD = os.getcwd()

def blast(command):
    output = docker.run(
        image="ncbi/blast",
        remove=True,
        volumes=[(f"{CWD}/blast/blastdb", "/blast/blastdb"), (f"{CWD}/blast/blastdb_custom", "/blast/blastdb_custom"),
                 (f"{CWD}/blast/fasta", "/blast/fasta"), (f"{CWD}/blast/queries", "/blast/queries")],
        command=command)
    return output


output = blast(["efetch", "-db", "protein", "-format", "fasta", "-id", "NP_000509"])

print(output)

Pots veure que la funció blast et retorna el fitxer fasta de l’hemoglobina beta:

>NP_000509.1 hemoglobin subunit beta [Homo sapiens]
MVHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPDAVMGNPKVKAHGKKVLG
AFSDGLAHLDNLKGTFATLSELHCDKLHVDPENFRLLGNVLVCVLAHHFGKEFTPPVQAAYQKVVAGVAN
ALAHKYH

Al realitzar una consulta amb blastp pots utilitzar l’argument -outfmt per especificar un format que puguis utilitzar directament amb Python.

A continuació es mostren tots els formats disponibles en JSON (amb la comanda blastp -help pots veure els formats disponibles):

$ blastp --help
...
	12 = Seqalign (JSON),
	13 = Multiple-file BLAST JSON,
	15 = Single-file BLAST JSON,
...

Per analitzar un resultat, pots indicar que la sortida sigui en BLAST JSON tal com es mostra a continuació:

record = blast(["blastp", "-query", "/blast/queries/P01349.fasta", "-db", "pdbaa", "-outfmt", "15"])
record = json.loads(record)
hits = record["BlastOutput2"][0]["report"]["results"]["search"]["hits"]
print(json.dumps(hits[0], indent=4))

Pots veure com tens accés a tota la informació al convertir l’string JSON en un diccionari:`

{
	"num": 1,
	"description": [
    	{
        	"id": "pdb|2FHW|B",
        	"accession": "2FHW_B",
        	"title": "Chain B, Relaxin 3 (Prorelaxin H3) (Insulin-like peptide INSL7) (Insulin-like peptide 7) [synthetic construct]",
...

Cerca BLAST

A continuació es descriu els diferents passos per fer una consulta.

Especificar la seqüència d'interès

Una cerca BLAST comença amb la selecció d'una seqüència d'ADN o proteïna a analitzar.

La seqüència es pot passar directament – per exemple exemple en el format FASTA – o mitjançant el número d'accés – per exemple, un número d'identificació RefSeq o GenBank (GI).

La cerca BLAST també et permet seleccionar un subconjunt d'una seqüència de consulta sencera, com ara una regió o un domini d'interès.

Seleccionar el programa

La família de programes NCBI BLAST inclou cinc programes principals:

  1. El programa BLASTP (protein) compara una seqüència de consulta d'aminoàcids amb una base de dades de seqüències de proteïnes. Tingues en compte que per a aquest tipus de cerca hi ha paràmetres opcionals que són especialment rellevants per a les cerques de proteïnes, com ara l'elecció de diverses matrius de puntuació PAM i BLOSUM.

  2. El programa BLASTN (nucleotid) s'utilitza per comparar una seqüència de consulta de nucleòtids amb una base de dades de seqüències de nucleòtids.

    Ja saps que qualsevol seqüència d'ADN es pot transcriure i traduir en sis marcs de lectura potencials (tres a la cadena superior i tres a la cadena inferior). Pert tant, donada una seqüència d’ADN, tenim 6 seqüències proteïnes potencials.

Els altres tres programes utilitzen aquest fet per realitzar alineaments proteïna-proteïna.

  1. El programa BLASTX compara una seqüència de consulta de nucleòtids, traduïda a 6 proteïnes potencials, amb una base de dades de seqüències de proteïnes.

    Per tant, si tens una seqüència d'ADN i vols saber si codifica alguna proteïna pots fer una cerca blastx.

  2. El programa TBLASTN compara una seqüència de consulta de proteïnes, traduïda a 6 seqüencies d’ADN, amb una base de dades de seqüències de nucleòtids.

    Per tant, si tens una base de dades d'ADN i vols saber si quines seqüències de la base de dades codifiquen una proteïna, pots fer una cerca tblastn.

  3. El programa TBLASTX compara una seqüència de consulta de nucleòtids, traduïda a 6 proteïnes, que a la seva vegada es tradueix cada proteïna a 6 seqüencies d’ADN, amb una base de seqüències de nucleòtids.

    Per exemple, si tens una seqüència d'ADN sense coincidències de base de dades evidents i vols saber si codifica una proteïna amb coincidències de base de dades distants i estadísticament significatives en una base de dades d'etiquetes de seqüències expressades, una cerca tblastx és més sensible que blastn i, per tant, útil per revelar gens que codifiquen proteïnes homòlogues a la teva consulta.

Selecció d'una base de dades

Les bases de dades disponibles per a la cerca BLAST es mostren a cada pàgina BLAST.

  1. Per a les cerques de bases de dades de proteïnes (BLASTP i BLASTX), l'opció predeterminada és la base de dades no redundant (nr). Consisteix en els registres de proteïnes combinats de GenBank, el Protein Data Bank (PDB), SwissProt, PIR i PRF. Una altra opció és cercar només proteïnes Refeq.

  2. Per a les cerques de bases de dades d'ADN (BLASTN, TBLASTN, TBLASTX) l'opció predeterminada és cercar a la base de dades de nucleòtids nr/nt. Això inclou seqüències de nucleòtids de GenBank, EMBL, DDBJ, PDB i RefSeq. No obstant això, la base de dades nr no té registres de l'EST, el lloc etiquetat per seqüències (STS), la seqüència de genoma sencer (WGS), la seqüència d'enquesta del genoma (GSS), l'assemblatge del transcriptoma (TSA), les patents o la seqüència genòmica d'alt rendiment. bases de dades (HTGS).

Altres opcions que s'utilitzen habitualment inclouen la base de dades genòmica més transcripció humana (o ratolí) o la base de dades EST.

Les bases de dades nr es deriven fusionant diverses bases de dades principals de proteïnes o ADN. Aquestes bases de dades sovint contenen seqüències idèntiques. Generalment, la base de dades nr conserva només una d'aquestes seqüències, juntament amb diversos números d'accés. (Encara que dues seqüències de la base de dades nr semblen ser idèntiques, almenys haurien de tenir alguna diferència subtil.) Les bases de dades nr són sovint els llocs preferits per cercar la majoria de seqüències disponibles.

Selecció de paràmetres de cerca opcionals

TODO

Revisar taules:

https://www.ncbi.nlm.nih.gov/books/NBK279684/

NCBI-BLAST API

La API de URL común NCBI-BLAST te permite realizar búsquedas de forma remota. Soporta lo mismo comandos en el servidor web NCBI y en la instalación de un proveedor de nube. Esto te permite ejecutar búsquedas en el servidor web del NCBI o en un proveedor de nube sin tener que modificar el código.

Pràctica

El NCBI té una API REST per fer consultes BLAST+, però amb un límit de 1 consulta cada 10 segons.

Necessites fer cerques BLAST de gran volum? Teniu dades de seqüència propietat per cercar i no podeu utilitzar el lloc web de NCBI BLAST? Voleu evitar la feina de muntar i mantenir BLAST autònom?

NCBI ofereix ElasticBLAST que pot executar les vostres cerques al núvol.

https://blast.ncbi.nlm.nih.gov/doc/blast-help/cloudblast.html

Una altra opció és que tu mateix construeixis un aplicació web que permeti fer les consultes en un servidor flask.

Es tracta que no només provis els exemples proporcionats en aquestes activitats sinó que usis el teu nou servidor local per altres BLAST que consideris interessants.

Ja sigui cercant exemples de BLAST ja realitzats de fonts fiables o bé experimentant amb gens i organismes que et resultin interessants.

Algunes fonts que pots provar si vols:

Millores

Pots crear una interfície web.