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.shExecuta 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.fastaPots 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 kinaseA 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.ptoPots 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 7801Consultar 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.fastaMira 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
QLCGRGFIRAIIFACGGSRWATSPAMSIKCCIYGCTKKDISVLCA 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.ptoA 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 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
...1.- Descarrega i descomprimeix la base de dades swissprot i realitza una consulta BLAST a la proteïna P01929 que guardi els resultats en un fitxer a la carpeta results.
update_blastdb.pl --source gcp --decompress swissprotefetch -db protein -format fasta -id P01929 > queries/P01929.fastablastp -query queries/P01929.fasta -db swissprot -out results/P01929_blast_results.txtAlguns resultats obtinguts:
Query= sp|P01929.1|HBA_LEOFU RecName: Full=Hemoglobin subunit alpha;
AltName: Full=Alpha-globin; AltName: Full=Hemoglobin alpha chain;
Contains: RecName: Full=Hemopressin
Length=141
Score E
Sequences producing significant alignments: (Bits) Value
P01929.1 RecName: Full=Hemoglobin subunit alpha; AltName: Full=Al... 282 1e-98
P63107.2 RecName: Full=Hemoglobin subunit alpha; AltName: Full=Al... 277 1e-96
P67817.2 RecName: Full=Hemoglobin subunit alpha; AltName: Full=Al... 276 2e-96
P67818.1 RecName: Full=Hemoglobin subunit alpha; AltName: Full=Al... 276 3e-96
P01923.1 RecName: Full=Hemoglobin subunit alpha; AltName: Full=Al... 275 7e-96Algoritme
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 20^3 = 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; que funcionarà fora del contenidor si has seguit totes les indicacions prèvies:
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:
$ python3 blastdemo.py
>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:
-
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.
-
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.
-
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. -
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.
-
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.
-
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.
-
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.
La API REST de l'NCBI et permet realitzar consultes BLAST+, però amb un límit de 1 consulta cada 10 segons.
- Necessites fer cerques BLAST de gran volum?
- Tens dades de seqüència propietaries per cercar i no pots 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.
Una altra opció és que tu mateix construeixis un aplicació web que permeti fer les consultes en un servidor.
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.