BLAST ejecuta alineamientos locales entre una secuencia “problema” o query, generalmente no muy larga, y un conjunto grande de secuencias (la base de datos) entre las que esperamos encontrar alguna secuencia homóloga a la query. Cada secuencia de la base de datos que pueda alinearse con la query se la llama subject en el alineamiento. Los objetivos de esta práctica son:
Identificar proteínas ortólogas entre dos especies con BLASTP.
Realizar un PSI-BLAST contra la base de datos Swissprot.
Determinar el origen de las lecturas que no pudieron ser mapeadas a su genoma de referencia, usando BLASTN remoto.
Este guion se dividirá en tres ejercicios con diferentes objetivos.
Ejercicio 1
Nota: Para realizar este ejercicio primero debemos descargar desde uniprot los proteomas de Buchnera aphidicola (107806) y Blochmannia ocreatus (251538). Estas dos son especies bacterianas endosimbiontes. Los guardamos como dos archivos FASTA (están descomprimidos ambos): ../../data/buchnera.fa y ../../data/blochmannia.fa.
Una estrategia habitual para identificar genes ortólogos es comprobar si cada uno es, respecto del otro, su mejor resultado de BLAST. Lo que en inglés llaman reciprocal best hits. Aunque se trata de un método heurístico, los resultados son bastante fiables. Vamos a practicarlo entre dos proteomas que debes haber descargado, y que primero debemos enmascarar e indexar para BLAST.
Enmascarar e indexar los proteomas
En primer lugar enmascaramos e indexeamos los datos. Primero lo hacemos con Buchnera:
pwdif[!-e ../../data/buchnera.paa ];then# Creamos un archivo para enmascarar zonas de baja complejidad:if[!-e buchnera_seg.asnb ];thensegmasker-in ../../data/buchnera.fasta \-infmt fasta \-parse_seqids\-outfmt maskinfo_asn1_bin \-out ../../data/buchnera_seg.asnbfi# Y indexamos el archivo FASTA para BLAST:makeblastdb-in ../../data/buchnera.fasta \-input_type fasta \-dbtype prot \-parse_seqids\-mask_data ../../data/buchnera_seg.asnb \-out ../../data/buchnera \-title"Proteoma de Buchnera"fi
/home/alicia/Website/results/2026-03-24
Ahora lo repetimos con Blochmannia:
pwdif[!-e ../../data/blochmania.paa ];then# Creamos un archivo para enmascarar zonas de baja complejidad:if[!-e blochmania_seg.asnb ];thensegmasker-in ../../data/blochmania.fasta \-infmt fasta \-parse_seqids\-outfmt maskinfo_asn1_bin \-out ../../data/blochmania_seg.asnbfi# Y indexamos el archivo FASTA para BLAST:makeblastdb-in ../../data/blochmania.fasta \-input_type fasta \-dbtype prot \-parse_seqids\-mask_data ../../data/blochmania_seg.asnb \-out ../../data/blochmania \-title"Proteoma de Blochmannia"fi
/home/alicia/Website/results/2026-03-24
Realizar las búsquedas de BLASTP
A continuación realizamos la búsqueda BLASTP. El objetivo es encontrar la proteína de Blochmannia ocreatus que más se parece a cada una de las proteínas de Buchnera aphidicola, y viceversa. Entonces, seleccionaremos las parejas de proteínas que sean recíprocamente las más similares una de otra, por ser más probablemente ortólogas. Otras relaciones de semejanza no recíprocas pueden ser debidas a la paralogía.
En primer lugar lo hacemos de Buchnera aphidicola hacia Blochmannia ocreatus, para ello:
Warning: [blastp] Examining 5 or more matches is recommended
NOTA:Revisa los argumentos utilizados, prueba otros si quieres y asegúrate de que son adecuados. Puedes consultar la ayuda con blastp -help en el terminal.
Comparación de resultados
Buscamos comparar los resultados. Para ello lo hacemos de esta forma:
NOTA:Echa un vistazo a estas tablas de resultados. ¿Sabrías identificar las parejas de proteínas que son recíprocamente la que más se parece una a la otra?
Vemos que hay 433 coincidencias (he escrito “coincidencias” en la consola y obtengo 433 resultados).
¡OJO! Cuando utilizamos secuencias de proteínas para estudiar la taxonomía de las especies buscamos copias ortologas. Una forma de identificar los ortologos es descargando ambos proteomas, hacemos un BLAST y podemos encontrar más de un hit. Tenemos dos listas de resultados al realizar el BLAST. Cuando tenemos esto en R, utilizamos las herramientas para comparar y ver si entre ambas listas qué parejas son recíprocas (son el mejor hit de ambas).
Usamos la función match para encontrar en la primera tabla de proteína de Blochmania la que coincide en Buchnera. Queremos ver si cuando buscamos en Buchnera también aparece la mimsa en Blochmania. Por eso usamos la función reverse.
NOTA:¿Qué datos te gustaría representar gráficamente para ganar confianza en que los resultados son correctos y tienen sentido?
Quizá nos interesaría representar la longitud de la proteína de ambas especies.
Ejercicio 2
Nota: Para realizar este ejercicio descargamos la base de datos de Swissprot para BLAST. Los archivos descargados y descomprimidos deben guardarse en la carpeta ../../data. Por tanto, los comandos siguientes deben ejecutarse desde esa carpeta, en el terminal de Bash.
El PSI-BLAST nos permite identificar semejanzas sútiles, como las que existen entre miembros alejados de una superfamilia proteica. Para comprobarlo vamos a realizar una búsqueda a partir de cualquiera de las proteínas siguientes: Q9CC91, P39615, B1AJI8, P47343, Q5UPT2, Q6UDG6, A1AFY9, O59825 y P30314. Todas ellas están en Swissprot y pertenecen a una misma superfamilia. Sin embargo, un simple BLASTP con cualquiera de ellas no encuentra ninguna de las otras1. Comprobemos si con un PSI-BLAST lo conseguimos. A continuación lo intento a partir de P39615. Repite los pasos usando otra proteína de la lista como query y compara los resultados.
# Como Swissprot debería estar ya descargada en `../../data`, podemos sacar# de ahí la secuencia de la proteína escogida:blastdbcmd-db ../../data/swissprot -entry P39615 -out P39615.fa# Y ahora la usamos como query en psi-blast:psiblast-db ../../data/swissprot \-query P39615.fa \-out P39615.blastout.txt \-evalue 0.1 \-outfmt 7 \-inclusion_ethresh 0.001 \-num_iterations 7 \-matrix BLOSUM62
Warning: [psiblast] Query_1 P39615.1 RecName:.. : Composition-based score adjustment conditioned on sequence properties and unconditional composition-based score adjustment is not supported with PSSMs, resetting to default value of standard composition-based statistics
Lo probamos con otra proteína de la superfamilia, en este caso B1AJI8 (he copiado el mismo código que el anterior pero lo he ajustado a B1AJI8):
# Como Swissprot debería estar ya descargada en `../../data`, podemos sacar# de ahí la secuencia de la proteína escogida:blastdbcmd-db ../../data/swissprot -entry B1AJI8 -out B1AJI8.fa# Y ahora la usamos como query en psi-blast:psiblast-db ../../data/swissprot \-query B1AJI8.fa \-out B1AJI8.blastout.txt \-evalue 0.1 \-outfmt 7 \-inclusion_ethresh 0.001 \-num_iterations 7 \-matrix BLOSUM62
Warning: [psiblast] Query_1 B1AJI8.1 RecName:.. : Composition-based score adjustment conditioned on sequence properties and unconditional composition-based score adjustment is not supported with PSSMs, resetting to default value of standard composition-based statistics
NOTA:Consulta la ayuda de psiblast y asegúrate de entender para qué sirven los argumentos que hemos usado. Por ejemplo, ¿cuál es la diferencia entre -evalue y -inclusion_ethresh?¿Cómo harías que la búsqueda fuera más sensible todavía?
Inclusion_ethresh sirve para definir el valor umbral de e-valor, mientras que el e-valor nos indica el número de hits aleatorios esperados por casualidad al buscar en una base datos. Para aumentar la sensibilidad de la búsqueda se repite el psi-blast sucesivamente hasta llegar a una convergencia.
Ejercicio 3
Descargamos mediante el enlace del guion del Aula Virtual el archivo comprimido contaminacion.fa.gz. No debe abrirse ni descomprimirse. Sólo guardadlo en la carpeta de datos tal como está. En adelante supondré que podéis acceder al archivo mediante ../../data/contaminacion.fa.gz desde la carpeta de la práctica.
El archivo ../../data/contaminacion.fa.gz contiene 33174 lecturas procedentes de un experimento de secuenciación del insecto Aphis aurantii, con Ion Torrent. (SRR13258790). Concretamente son las lecturas que no han podido mapearse al genoma de referencia (GCA_055615905) de este insecto. Hay dos motivos principales por los que muchos fragmentos no pueden mapearse al genoma de referencia: o bien el genoma de referencia es incompleto, o bien los fragmentos no pertenecen a esa especie. La secuenciación no distingue al individuo de interés de sus simbiontes, parásitos o presas. Los fragmentos de ADN secuenciados que no proceden del individuo diana se consideran contaminación. El objetivo de este ejercicio es determinar la procedencia de estos fragmentos mediante BLASTN remoto a la base de datos nr (no redundante). Lo que nos interesa es sólo el origen taxonómico de la secuencia más parecida a cada una de esas lecturas.
Selección aleatoria de 20 queries
Si enviamos 33174 queries a los servidores del NCBI, o bien tardarán muchas horas en devolver los resultados o simplemente interrumpirán el servicio. Cada persona escogerá 20 secuencias al azar mediante las órdenes siguientes, y después compartiremos los resultados para tener una muestra más representativa.
library('Biostrings')
Loading required package: BiocGenerics
Loading required package: generics
Attaching package: 'generics'
The following objects are masked from 'package:base':
as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
setequal, union
Attaching package: 'BiocGenerics'
The following objects are masked from 'package:stats':
IQR, mad, sd, var, xtabs
La orden siguiente pide una sola secuencia (la más parecida) y alguna información sobre el alineamiento, para poder garantizar la calidad del resultado.
La única información taxonómica que podemos obtener de las secuencias encontradas cuando realizamos un BLAST remoto es el identificador taxonómico. Esta restricción tiene por objetivo agilizar el servicio remoto. Ahora es asunto nuestro saber a qué especies corresponden los identificadores, de manera automática. Hay varias opciones. Para no tener que descargar la base de datos taxonómica, intentaremos recuperar la información completa con una nueva consulta remota. Esta vez, con los programas datasets y dataformat de NCBI, que ya conocéis.
New version of client (18.26.0) available at https://ftp.ncbi.nlm.nih.gov/pub/datasets/command-line/LATEST/linux-amd64/datasets.
Nota
Puedes importar tanto los resultados del BLAST como el informe taxonómico a una sesión de R y determinar que fracción de las lecturas seleccionadas parecen proceder de un hemímptero y qué fracción proceden de otros organismos. Tienen sentido los resultados? Encuentras alguna secuencia de origen vegetal?
# Extraemos la columna 9 (TaxID), ordenamos y quitamos duplicadoscut-f 9 resultado.txt |sort|uniq> taxids_lista.txt# Comprobamos qué IDs tenemos (deberías ver números como 1470492, 80765, etc.)cat taxids_lista.txt
1470492
44670
80765
# Consultamos la taxonomía y la guardamos en formato TSV/home/alicia/bin/datasets summary taxonomy taxon --inputfile taxids_lista.txt --as-json-lines|\/home/alicia/bin/dataformat tsv taxonomy --template tax-summary > taxreport.txt# Verificamos que se haya escrito correctamentehead taxreport.txt
New version of client (18.26.0) available at https://ftp.ncbi.nlm.nih.gov/pub/datasets/command-line/LATEST/linux-amd64/datasets.
Query Taxid Tax name Authority Rank Basionym Basionym authority Curator common name Has type material Group name Domain/Realm name Domain/Realm taxid Kingdom name Kingdom taxid Phylum name Phylum taxid Class name Class taxid Order name Order taxid Family name Family taxid Genus name Genus taxid Species name Species taxid Scientific name is formal
1470492 1470492 Buchnera aphidicola (Aphis aurantii) FORMA_SPECIALIS no enterobacteria Bacteria 2 Pseudomonadati 3379134 Pseudomonadota 1224 Gammaproteobacteria 1236 Enterobacterales 91347 Erwiniaceae 1903409 Buchnera 32199 Buchnera aphidicola 9 FALSE
44670 44670 Metopolophium dirhodum (Walker, 1849) SPECIES rose-grain aphid no aphids Eukaryota 2759 Metazoa 33208 Arthropoda 6656 Insecta 50557 Hemiptera 7524 Aphididae 27482 Metopolophium 44669 Metopolophium dirhodum 44670 TRUE
80765 80765 Aphis gossypii Glover, 1877 SPECIES cotton aphid no aphids Eukaryota 2759 Metazoa 33208 Arthropoda 6656 Insecta 50557 Hemiptera 7524 Aphididae 27482 Aphis 80764 Aphis gossypii 80765 TRUE
# Leemos los resultados del BLASTblast <-read.table('resultado.txt',col.names =c('qaccver', 'saccver', 'length', 'qlen','qstart', 'qend', 'evalue', 'bitscore', 'staxid'))# Leemos el informe taxonómico (usamos quote="" para evitar que se cuelgue)taxreport <-read.table('taxreport.txt', sep ='\t', header =TRUE, fill =TRUE, quote ="")# Unimos las tablas para asignar el Orden y el Reino a cada hitblast$orden <- taxreport[match(blast$staxid, taxreport$Taxid), 'Order.name']blast$reino <- taxreport[match(blast$staxid, taxreport$Taxid), 'Kingdom.name']# CÁLCULOS # 1. Fracción de Hemípteroshemipteros <-sum(blast$orden =='Hemiptera', na.rm =TRUE)frac_hemipteros <- hemipteros /nrow(blast)# 2. Fracción de Plantas (Viridiplantae)plantas <-sum(blast$reino =='Viridiplantae', na.rm =TRUE)frac_plantas <- plantas /nrow(blast)# --- RESULTADOS ---cat('Fracción hemípteros:', frac_hemipteros, '\n')