BLAST

En este capitulo vamos a explorar el uso de blast en local, es decir, instalado en nuestras máquinas. Lo primero que teneis que hacer, si no lo habeis hecho ya, es seguir las instrucciones para instalar Blast en vuestra máquina, podeis encontrarlo aqui.

En este tema solo se van a dar unas pinceladas de como funciona Blast y sus parámetro. Para el resto tendreis que visitar la documentación y su manual https://www.ncbi.nlm.nih.gov/books/NBK279690/ (esta disponible para descargar en PDF).

Homologia de secuencias

En biología, dos genes se definen como homólogos si derivan de un ancestro común. Pero la homología no termina ahí, ya que dos secuencias pueden ser homólogas ortólogas o parálogas, tal y como podemos observar en la Figura1.


Figura 1. Evolución de los genes de las globinas

Ortólogos Secuencias similares en dos os mas organismos diferentes que provienen de un proceso de especiación. Se espera que tengan la misma función biológica.

parálogos Secuecuencias similares dentro de un mismo organismo que proceden de un proceso de duplicación génica. No tienen porque tener la misma función


Alineamiento de secuencias

Un alineamiento de secuencias no es mas que el reordenamiento de dos cadenas de caracteres para conseguir el mayor numero de matches posibles entre ambas. El mejor alineamiento posible es aquel que muestra un mayor nivel de correspondencia entre ambas secuencias, y el menor numero de diferencias o discrepancias. Dos secuencias similares tienen una alta probabilidad de tener una misma función, pero siempre se cumple (recordar parálogos). Términos importantes:


identidad: El mismo aminoácido o nucleotido en la misma posicion en dos secuencias.

similaridad: Numero de positivos entre ambas secuencias, es decir, sustituciones de aminoácidos o nucleótidos con las mismas propiedades bioquímica en dos secuencias.


>Importante similaridad no es igual que homología

Si las secuencias son mayores de 100 nucleótidos o 100 aminoácidos, puedes decir que dos proteinas son homólogas si al menos el 25 % de los aminoácidos son idénticos. En el caso del DNA se requiere un porcentage de identidad del 70 %.

Hasta la fecha miles de genomas de diferentes especies se han secuenciado y están disponibles para que trabajemos con ellas, revelando información valiosa sobre sus genes y su función correspondiente. Para asignar funciones a los genes de nuevos genomas recien secuenciados se utilizan herramientas que buscan y miden los niveles de similitud entre millones de secuencias, ahí es donde entran en juegos herramientas como BLAST y otros alineadores de secuencias.


BLAST

BLAST es el acrónimo de Basic Local Alignment Search Tool. es el algoritmo mas empleado por el NCBI. es un programa informáJco de alineamiento de secuencias de tipo local, ya sea de ADN, ARN o de proteínas. El programa es capaz de comparar una secuencia problema (query) contra una gran canJdad de secuencias que se encuentren en una base de datos (subject).


Figura 3. Portada libro sobre BLAST.

Tipos de Blast

Podemos usar todos los tipos de blast con blast local, es decir, blastx, blastp, blastn y tblastn al instalarnos blast+ de NCBI en local. Recordad que:



  • blastn: nucleotide > nucleotide (nucleotide blast)
  • blastp : protein > protein (protein blast)
  • tblastn : protein > translated nucleotide
  • blastx: translated nucleotide > protein

Syntaxis general de blast

Cada vez que utilicemos blast usaremos un comando diferente segun el tipo de blast de los listados arriba vayamos a usar, es decir (blastp, blastn, blastx o tblastn). Además, es obligatorio indicar el archivo fasta que contiene la secuencia query (-query), el archivo con la/s secuencia/s subjects (-subject) e indicar el output si queremos guardar el resultado en un archivo con la opción –out

~$ [t]blast[pnx] -query fasta -subject fasta  -out result.out

query: secuencia problema fasta

subject: archivo fasta con la/s secuencia/s que enfrentaremos a nuestra secuencia problema

out: nombre del archivo de salida


veamos un ejemplo con blastp. Para ello usaremos los archivos fasta PBP1_staphylococcus.fasta como query PBP1.fasta como subject. Estos archivos contienen proteinas PBPs que son proteinas de union a la penicilina de diferentes bacterias.

~$ blastp -query PBP1_staphylococcus.fasta -subject PBP1.fasta -out result.out

...
Query= ASD49649.1 penicillin binding protein 1 [Staphylococcus aureus]

Length=744
                                                                    Score     E
Sequences producing significant alignments:                         (Bits) Value

ANW82078.1 penicillin binding protein 1 [Staphylococcus aureus]     1522    0.0  
ASD49626.1 penicillin binding protein 1 [Staphylococcus aureus]     1522    0.0  
ASD49630.1 penicillin binding protein 1 [Staphylococcus aureus]     1522    0.0  
ASD49632.1 penicillin binding protein 1 [Staphylococcus aureus]     1522    0.0  
ASD49633.1 penicillin binding protein 1 [Staphylococcus aureus]     1522    0.0  
ASD49643.1 penicillin binding protein 1 [Staphylococcus aureus]     1522    0.0  
ASD49646.1 penicillin binding protein 1 [Staphylococcus aureus]     1522    0.0  
ASD49647.1 penicillin binding protein 1 [Staphylococcus aureus]     1522    0.0  

Se obtiene un archivo con una estructura de datos muy similar a la que vemos cuando usamos la version web, sin embargo, el problema es que estas estructura es muy dificil de parsear. Las opciones de salida de datos que permite Blast son extensas y aqui solo vamos a mostrar algunas, os recomiendo que mireis el manual, si qureis profundizar en el tema.

En la Tabla 1 teneis algunas de las opciones disponibles comunes a todos los tipos diferentes de Blast.


Tabla 1. Opciones generales de Blast

optiontypedefault value description and notes
dbstringnoneBLAST database name.
querystringstdinQuery file name.
query_locstringnoneLocation on the query sequence (Format: start-stop)
outstringstdoutOutput file name
evaluereal10Expect value (E) for saving hits
subjectstringnoneFile with subject sequence(s) to search.
show_gisflagN/AShow NCBI GIs in report.
num_descriptionsinteger500Show one-line descriptions for this number of database sequences.
num_alignmentsinteger250Show alignments for this number of database sequences.
htmlflagN/AProduce HTML output
giliststringnoneRestrict search of database to GI’s listed in this file. Local searches only.
negative_giliststringnoneRestrict search of database to everything except the GI’s listed in this file. Local searches only.
entrez_querystringnoneRestrict search with the given Entrez query. Remote searches only.
best_hit_overhangrealnoneBest Hit algorithm overhang value (recommended value: 0.1)
best_hit_score_edgerealnoneBest Hit algorithm score edge value (recommended value: 0.1)
dbsizeintegernoneEffective size of the database
searchspintegernoneEffective length of the search space
import_search_strategystringnoneSearch strategy file to read.
export_search_strategystringnoneRecord search strategy to this file.
parse_deflinesflagN/AParse query and subject bar delimited sequence identifiers (e.g., gi|129295).
num_threadsinteger1Number of threads (CPUs) to use in blast search.
remoteflagN/AExecute search on NCBI servers?


output format

Si utilizamos la opción -outfmt podremos usar diferentes estilos de salida de datos (18 en total), nosotros usaremos la numero 6 (tabular), veamos que aspecto tiene:


~$ blastp -query PBP1_staphylococcus.fasta -subject PBP1.fasta -outfmt "6"

ASD49649.1 ANW82078.1 100.000 744 0 0 1 744 1 744 0.0 1522
ASD49649.1 ASD49626.1 100.000 744 0 0 1 744 1 744 0.0 1522
ASD49649.1 ASD49630.1 100.000 744 0 0 1 744 1 744 0.0 1522
ASD49649.1 ASD49632.1 100.000 744 0 0 1 744 1 744 0.0 1522
ASD49649.1 ASD49633.1 100.000 744 0 0 1 744 1 744 0.0 1522
ASD49649.1 ASD49643.1 100.000 744 0 0 1 744 1 744 0.0 1522
ASD49649.1 ASD49646.1 100.000 744 0 0 1 744 1 744 0.0 1522
ASD49649.1 ASD49647.1 100.000 744 0 0 1 744 1 744 0.0 1522
ASD49649.1 ASD49649.1 100.000 744 0 0 1 744 1 744 0.0 1522
ASD49649.1 ASD49651.1 100.000 744 0 0 1 744 1 744 0.0 1522
...

Pero quizas a veces necesitemos menos información. Utilizando la opción -outfmt tambien podemos especificar que campos son los que queremos mostrar, en este ejemplo, mostraremos las siguientes opciones que son de las mas importantes en cualquier analisis de blast:

  • qseqid: query seq ID
  • sseqid: subject seq ID
  • pident: Percentage of identical matches
  • length: Alignment length
  • mismatch: Number of mismatches
  • gapopen: Number of gap openings
  • qstart: Start of alignment in query
  • qend: End of alignment in query
  • sstart: Start of alignment in subject
  • send: End of alignment in subject
  • evalue: Expect value (E) for saving hits
  • bitscore: Bit score

como podeis ver ahora obtenemos un formato en columna mucho mas apetecible de ser parseado. las opciones mostradas por default en la opción 6 son por orden:

opciones que son de las mas usadas en cualquier análisis de blast:

qseqid

sseqid

qcovs

pident

evalue

~$ blastp -query PBP1_staphylococcus.fasta -subject PBP1.fasta -outfmt "6 qseqid sseqid qcovs pident evalue"
 ASD49649.1    ANW82078.1  100 100.000 0.0
 ASD49649.1    ASD49626.1  100 100.000 0.0
 ASD49649.1    ASD49630.1  100 100.000 0.0
 ASD49649.1    ASD49632.1  100 100.000 0.0
 ASD49649.1    ASD49633.1  100 100.000 0.0
 ASD49649.1    ASD49643.1  100 100.000 0.0
 ASD49649.1    ASD49646.1  100 100.000 0.0
 ASD49649.1    ASD49647.1  100 100.000 0.0
 ASD49649.1    ASD49649.1  100 100.000 0.0
 ASD49649.1    ASD49651.1  100 100.000 0.0
 ASD49649.1    ANW82077.1  100 99.866  0.0
 ASD49649.1    ASD49623.1  100 99.866  0.0
 …

tambien podemos controlar el valor de cut off del e-value por defecto es 10, pero no deberíamos de considerar resultados que no estén por debajo de 1xe-5, a menos que sepamos muy bien lo que estamos haciendo.

Por ejemplo veamos que pasa si filtramos por e-value de forma restrictiva utilizando la opción -evalue:

~$ blastp -query PBP1_staphylococcus.fasta -subject PBP1.fasta -evalue "0.0000001" -outfmt "6 qseqid sseqid qcovs pident evalue"

vereis que al ser mas restrictivo a nivel de e-value hemos obtenido una cantidad menor de hits.


Bases de datos de Blast

En vez de analizar mediante Blast un archivo Fasta, podemos generar o descargar de NCBI bases de datos de BLAST con secuencias de ADN o proteínas. La ventaja de usar bases de datos, es que la velocidad de analisis y comparación de secuencias de BLAST aumentará notablemente, especialmente cuanto utilizemos un conjunto de secuencias subject del orden de cientos de miles o millones.


Crear una base de datos en blast es muy sencillo:

~$ makeblastdb –in mydb.fsa –dbtype "nucl/prot" –out mydb

in: el multifasta con las secuencias que formarán la base de datos.

dbtype: Indica el tipo de datos que van a formar la base de datos, proteinas (prot) o nucleótidos (nucl).

out: el nombre que le pondremos a la base de datos


Hagamos un ejemplo con nuestras secuencias PBPs, vamos a crear una base de datos que se llamará:

~$ makeblastdb -in PBP1.fasta -dbtype 'prot' -out PBP1

Building a new DB, current time: 04/24/2020 09:50:15
New DB name:   /Users/JGL/Dropbox/programacion_bioinformatica/Tema_5_blast/PBP1
New DB title: PBP1.fasta
Sequence type: Protein
Keep MBits: T
Maximum file size: 1000000000B
Adding sequences from FASTA; added 54 sequences in 0.00509191 seconds.

Ahora hechemos un vistazo a los archivos que se han creado:

-rw-r--r--@ 1 JGL  staff   6.8K Apr 24 09:50 PBP1.phr
-rw-r--r--@ 1 JGL staff   504B Apr 24 09:50 PBP1.pin
-rw-r--r--@ 1 JGL staff   34K Apr 24 09:50 PBP1.psq

Lo mas interesante es guardar las bases de datos en un directorio, con el mismo nombre. En este caso por ejemplo PBP1/.

Como podemos observar las bases de datos de Blast contienen tres clases de archivos, .phr .pin y .psq.

Una vez creada nuestra base de datos, hacer un blast frente a la misma es muy fácil, solo tenemos que

~$ blastp -query PBP1_staphylococcus.fasta -db PBP1/PBP1 -outfmt "6 qseqid sseqid qcovs pident evalue"

Tambien podeis bajaros bases de datos de NCBI y utlizarlas para hacer análisis en local, de aqui podeis bajaros la que querais y explorarlas:

ftp://ftp.ncbi.nlm.nih.gov/blast/db/



Podeis copias el link address de cualquiera de ellas y descargadlas en un directorio con wget o curl. Utilizando el comando Blastdbcmd podréis hacer análisis de las bases datos que os descarguéis, ¡os animo a explorarlo!.