Con la base de datos BLAST que preparamos en el Taller 2 ya estamos listos para hacer búsquedas locales. Vamos a utilizar las secuencias de las proteinas encontradas en el suero de leche de cabras (Taller 1). ,Al pie de esta página hay una copia del archivo fasta que generamos (simil.cabras.fasta).
Las proteínas halladas en la leche se determinaron a través de similitudes con proteínas que en la mayoría de los casos no se corresponden a proteínas de cabra, sino a proteínas de otros animales. En el Taller 1 ya habíamos extraído estas secuencias (simil.cabras.fasta), y también hicimos un análisis del origen taxonómico de esas proteínas.
Mediante un análisis BLAST vamos a recuperar las proteínas homólogas de cabra.
En primer lugar hacemos un BLAST que genere una salida html, para explorar que características tienen los hits que podemos encontrar, los valores E y scores característicos y las coberturas.
Desde la linea de comando de linux corrermos:
blastp -query simil.cabras.fasta -db cabra.protein -evalue 1e-50 -html -num_threads 2 > test1.html
Los programas BLAST tiene muchos argumentos. Aqui pueden encontrar explicaciones y recomendaciones detalladas de configuración y uso. Y aqui una ayuda específica para los argumentos de las variantes de BLAST (blastp, blastn, blastx, etc.)
Veamos los argumentos que usamos en nuestra corrida:
-query: es el archivo donde tenemos las secuencias que queremos analizar
-db: el nombre de la base de datos BLAST a consultar (ver Taller 2)
-evalue: el umbral de los valores E (no nos interesan hits con alineamientos que tengan un E-value mayor que el especificado aqui)
-html: queremos una salida en formato html
-num_threads: número de núcleos que queremos usar para el análisis
Abramos el archivo test1.html y analicemos las siguientes cuestiones:
¿BLAST encontró hits para todas las secuencias?
¿Qué características tienen las consultas para las cuales no hubo hits?
¿Cómo son las coberturas de los hits?
¿Hay un número máximo de hits de muy buena calidad (alto score, buena cobertura, pocas no coincidencias)?
¿Qué pasa con las secuencias consultas que en Uniprot están etiquetadas como fragmento?
A partir de esto podemos plantear una segunda búsqueda exploratoria:
blastp -query simil.cabras.fasta -db cabra.protein -evalue 0.01 -html -num_threads 2 -num_descriptions 4 -num_alignments 4 > test2.html
Y ya estamos listos para correr el análisis BLAST definitivo
blastp -query simil.cabras.fasta -db cabra.protein -evalue 0.01 -max_target_seqs 4 -outfmt '6 qseqid sseqid pident ppos length mismatch gapopen qstart qend sstart send evalue bitscore qlen slen' -num_threads 2 > blastp.cabras.csv
Ya estamos listos para obtener estos resultados en forma de tabla para analizarlos en R. Por lo tanto, creamos un proyecto que podemos llamar "taller_3". Además copiamos a esa carpeta el archivo blastp.cabras.csv. Este archivo contiene en forma tabular todos los datos del análisis.
Este archivo csv no tiene encabezado. Los nombres de las variables son:
qseqid: ID de la consulta
sseqid: ID del hit
pident: % de coincidencias
ppos: % de positivos
length: largo del alineamiento
mismatch: cantidad de no coincidencias
gapopen: cantidad de inicios de gaps
qstart: inicio del alineamiento (coordenadas de la consulta)
qend: fin del alineamiento (coordenadas de la consulta)
sstart: inicio del alineamiento (coordenadas del hit)
send: fin del alineamiento (coordenadas del hit)
evalue: valor E
bitscore: score
qlen: largo de la consulta
slen: largo del hit
qcovs: cobertura de la consulta
Comenzamos el análisis en R leyendo el archivo con la salida del análisis BLAST, y le agregamos los nombres a las variables.
blast.cabras <- read.table("blastp.cabras.csv", header=F, stringsAsFactors = F)names(blast.cabras) <- c("qseqid", "sseqid", "pident", "ppos", "length", "mismatch", "gapopen", "qstart", "qend", "sstart", "send", "evalue", "bitscore", "qlen", "slen", "qcovs")
Luego verificamos el paso anterior
head(blast.cabras, 4)
## qseqid sseqid pident ## 1 tr|A0A0D3RIW8|A0A0D3RIW8_BUBBU gi|926694174|ref|XP_005681778.2| 89.47 ## 2 tr|A0A0D3RIW8|A0A0D3RIW8_BUBBU gi|926694176|ref|XP_013820153.1| 89.47 ## 3 tr|A0A0D5CBR5|A0A0D5CBR5_BUBBU gi|926695622|ref|XP_013820690.1| 94.92 ## 4 tr|A0A0D5CBR5|A0A0D5CBR5_BUBBU gi|926730042|ref|XP_005701550.2| 94.77 ## ppos length mismatch gapopen qstart qend sstart send evalue bitscore ## 1 92.11 76 6 1 1 76 182 255 4e-37 127 ## 2 92.11 76 6 1 1 76 182 255 4e-37 127 ## 3 97.43 1516 76 1 1 1515 1 1516 0e+00 2958 ## 4 97.87 1223 64 0 439 1661 1 1223 0e+00 2425 ## qlen slen qcovs ## 1 76 257 100 ## 2 76 257 100 ## 3 1661 1516 91 ## 4 1661 1223 74
dim(blast.cabras)
## [1] 1231 16
Para cada consulta BLAST devolvió varios hits. Al revisar la salida en formato HTML vimos que hay uno o más hits con alta similitud (hay proteínas codificadas por genes duplicados) Y después la calidad de los alineamientos baja rápidamente.
Podríamos establecer como criterio inicial de homología que consideraremos como homólogas a aquellas proteínas de cabra que tienen al menos un 79% de aminoácidos idénticos, 85% de aminoácidos similares (idénticos o positivos) y que la cobertura es de al menos un 85%. Es decir, tratándose de especies relacionadas es un criterio relativamente laxo.
Para extraer estos datos podemos utilizar las herramientas de consulta de R base que ya conocemos, Pero existe un paquete de R, dplyr que facilita el trabajo de filtrado y selección y que gradualmente se está convirtiendo en la herramienta estándar para filtrar, agrupar y en general manipular dataframes, especialmente si son grandes. Además funciona con tuberías (introducidas en el paquete magrittr), que en R se escriben como “%>%”
Como hicimos otras veces, descargamos el paquete y lo instalamos la primera vez que lo vamos a usar, y en las siguientes sesiones smplemente lo cargamos en la biblioteca de paquetes de la sesión.
# install.packages("dplyr")library(dplyr)
## ## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats': ## ## filter, lag
## The following objects are masked from 'package:base': ## ## intersect, setdiff, setequal, union
¿Cómo funciona dplyr?
Por ejemplo, para ver los IDs únicos de nuestras consultas en R estándar podríamos hacer:
unique( blast.cabras$qseqid )
Con dplyr sería
distinct(select(blast.cabras, qseqid))
Y con dplyr y pipes
blast.cabras %>% select(qseqid) %>% distinct() %>% dim()
En este último caso, el comando es más largo, pero es más claro que la opción con dplyr sin pipes. Definitivamente es más largo que la versión de R estándar, pero en seguida veremos que esto cambia cuando queremos hacer consultas más complejas
Por ejemplo supongamos que queremos calcular cuantos hits o hsp mayores de 500 aminoácidos obtuvimos por cada secuencia en nuestra consulta:
# R base: table( blast.cabras$qseqid[blast.cabras$qlen >= 500] ) # Con dplyr: blast.cabras %>% filter(qlen >= 500) %>% select(qseqid) %>% table()
En este caso el comando también es más largo, pero mucho más claro, ya que estamos separando nuestro operación de partes elementales. Primero filtramos, luego seleccionamos que campo queremos analizar y luego realizamos la tabla.
Para leer una introducción a las operaciones que se pueden hacer con dplyr, pueden consultar este tutorial
Ahora vamos a extraer los hits que cumplen con nuestro criterio de homología:
seqs.ok <- blast.cabras %>% filter(pident >=79, ppos >= 85, qcovs >= 95) dim(seqs.ok)
## [1] 319 16
length( unique(seqs.ok$qseqid) )
## [1] 212
Pudimos encontrar entonces al menos una potencial proteína de cabra homóloga para 212 proteínas reportadas en el análisis proteómico. A continuación habría que revisar si estas secuencias de cabra realmente se pueden considerar homólogas. Es posible que en los casos donde hay más de un hit, solo aquella proteína con los mejores indicadores sea la verdadero homóloga.
También deberíamos analizar que sucede con las consultas que no presentaron una similitud suficiente como para declarar una homóloga. Estas consultas son aquellas que no están presentes en el dataframe seqs.ok:
seqs.revisar <- !(blast.cabras$qseqid %in% seqs.ok$qseqid)seqs.no.ok <- blast.cabras[seqs.revisar,]
Para completar esta tarea, o como complemento, nos podría interesar guardar estos dos dataframes en archivos csv separados para después leerlos en Excel.
write.csv(seqs.ok, "candidatos.csv", row.names = F)write.csv(seqs.no.ok, "para.revisar.csv", row.names = F)
También seguramente vamos a necesitar los números de acceso de las secuencias de cabra para poder recuperar las secuencias fasta. En particular necesitamos los números de acceso Genbank, y sin el número de versión
accesos.crudos <- seqs.ok$sseqidaccesos.lista <- strsplit(accesos.crudos, "\\|")accesos.ok <- sapply(accesos.lista, function(x) x[4])accesos.ok <- sub("\\.\\d+$", "", accesos.ok)write(accesos.ok, "accesos.txt", sep="\\n")
El archivo accesos.txt lo podemos para conseguir las secuencias en formato fasta o genbank desde la herramienta BatchEntrez del NCBI.
Y también desde la linea de comandos de linux, interrogando a la base de datos BLAST que creamos anteriormente. Una aclaración, si el archivo accesos.txt se creó en Windows para que sirva para recuperar secuencias desde linux hace falta convertir el formato con el programa dos2linux.
sudo apt-get install dos2unix
dos2unix accesos.txt
blastdbcmd -db cabra.protein -entry_batch accesos.txt > proteinas_cabra_ok.fasta