Extraer números de acceso de un archivo Excel. Analizar el origen taxonómico de los accesos. Descargar secuencias desde Uniprot.
En este paper se investiga el perfil de proteínas presentes en el suero de la leche de razas griegas de oveja y de cabra.
Los autores utilizaron técnicas de proteómica para aislar, analizar y caracterizar unas 500 proteínas. Para este taller nos va a interesar un grupo más chico, que son aquellas proteínas que estaban compartidas por las dos razas de cabra analizadas (algunas proteínas eran exclusivas de cada raza).
Antes de comenzar a trabajar con R, tenemos que crear un proyecto RStudio nuevo.
Los datos con los que vamos a trabajar están en un archivo Excel que los autores depositaron como material suplementario (mmc8.xlsx). Este archivo lo podemos descargar directamente desde el link con cualquier navegador, aunque también lo podemos hacer desde R con la funcióndownload.file(). Esta función tiene dos argumentos obligatorios, la dirección (url) del archivo, y que nombre le queremos dar al descargarlo.
download.file( "http://www.sciencedirect.com/science/MiamiMultiMediaURL/1-s2.0-S1874391916301294/1-s2.0-S1874391916301294-mmc8.xlsx/276834/html/S1874391916301294/a8e2c42ea7eccd0415326d59e0e8f957/mmc8.xlsx", destfile = "mmc8.xlsx")
En cualquiera de los dos casos (descarga desde el navegador o desde R), el archivo tiene que quedar en la carpeta del proyecto. Podría ir en otra carpeta, pero en ese caso debemos indicar el directorio cada vez que lo necesitemos; y además, no es muy prolijo.
Es posible leer archivos Excel desde R, pero para eso necesitamos un paquete. La instalación de paquetes se puede hacer desde la barra de menú de RStudio (Tools > Install Packages…), o desde la línea de comando, con el comando install.packages(). La instalación se hace una sola vez, o al menos hasta que se distribuya una nueva versión del paquete. Luego, en cada sesión tenemos que cargar el o los paquetes que necesitemos a la biblioteca de paquetes disponibles para la sesión. Esto se hace con el comando library(). El motivo de esto es que podemos llegar a tener decenas o cientos de paquetes diferentes instalados en nuestra computadora, y si todos estuvieran disponibles en memoria aun cuando no los necesitamos, sería muy ineficiente, y un gasto de memoria.
# Correr una vez# install.packages("readxl")library(readxl)
Antes de proceder a leer el archivo vamos a abrirlo en Excel para ver que información contiene y cómo está organizada.
La planilla tiene dos hojas “Sheet1” y “Sheet3”. Sólo nos interesa la primera. Alli vemos que hay una primera fila con el título de la tabla. Luego, la tabla de datos tiene dos columnas, una para los números de acceso y otra para la descripción. La primera fila de la tabla contiene los nombres de los campos y a continuación están los datos. Esa es toda la información que necesitamos para leer el archivo y cargar los datos a un dataframe de R.
ids_crudos <- read_excel("mmc8.xlsx", sheet = 1, col_names = T, skip = 1)head(ids_crudos, 2)
## Accession ## 1 A0A0D3RIW8 ## 2 A0A0D5CBR5 ## Description ## 1 Beta-casein (Fragment) OS=Bubalus bubalis GN=CSN2 PE=4 SV=1 - [A0A0D3RIW8_BUBBU] ## 2 Complement component 3 OS=Bubalus bubalis GN=C3 PE=2 SV=1 - [A0A0D5CBR5_BUBBU]
¿Para qué sirven los argumentos que le pasamos a la función read_excel()?
El primer argumento es el nombre del archivo.
sheet: número de hoja queremos leer (el número según el orden en la planilla)
col_names: ¿la primera columna de la tabla contiene los nombres de las variables (TRUE/FALSE)
skip Número de filas al principio de la planilla que hay que saltear antes de encontrar la tabla (default = 0)
En el vector ids_crudos$Accession tenemos los números de acceso de las proteínas en la base de datos UNIPROT (http://www.uniprot.org).
En la otra columna del dataframe, ids_crudos$Description, tenemos la descripción de esos números de acceso.
head( ids_crudos$Description, 10 )
## [1] "Beta-casein (Fragment) OS=Bubalus bubalis GN=CSN2 PE=4 SV=1 - [A0A0D3RIW8_BUBBU]" ## [2] "Complement component 3 OS=Bubalus bubalis GN=C3 PE=2 SV=1 - [A0A0D5CBR5_BUBBU]" ## [3] "Myostatin (Fragment) OS=Ovis aries PE=3 SV=1 - [A0FEQ6_SHEEP]" ## [4] "L-lactate dehydrogenase OS=Bos mutus grunniens PE=2 SV=1 - [A0FH35_BOSMU]" ## [5] "Lipoprotein lipase OS=Capra hircus GN=LPL PE=2 SV=1 - [A0FI82_CAPHI]" ## [6] "Smooth muscle and non-muscle myosin alkali light chain peptide 6 (Fragment) OS=Bos taurus GN=MYL6 PE=2 SV=1 - [A1XEA7_BOVIN]" ## [7] "Xanthine oxidoreductase OS=Capra hircus PE=2 SV=1 - [A1YZ34_CAPHI]" ## [8] "Serpin A3-6 OS=Bos taurus GN=SERPINA3-6 PE=3 SV=1 - [SPA36_BOVIN]" ## [9] "Butyrophilin subfamily 1 member A1 OS=Capra hircus GN=BTN1A1 PE=2 SV=1 - [A3EY52_CAPHI]" ## [10] "Lactoperoxidase OS=Capra hircus GN=LPO PE=1 SV=1 - [A3F9D6_CAPHI]"
¿Por qué estas proteínas no pertenecen a cabra (Capra hircus)? Esto se debe a que la base de datos de fragmentos de CG-masa que usaron los autores para identificar las proteínas no incluía muchas proteínas de cabra, y el buscador devolvió los números de acceso de las proteínas más similares.
¿Será posible armar un listado de las especies a las que pertenecen estas proteínas? Si, por supuesto!
Primero observemos que al menos las primeras lineas de descripción tienen una clave seguida por un signo igual, “OS=”, seguida del género y especie.
¿Este texto “OS=” estará presente en todos los registros?
Importante: en lo que sigue vamos a trabajar con expresiones regulares. En el curso le vamos a dedicar un tiempo exclusivo a este tema
Hagamos primero una prueba:
texto1 <- "texto sin clave de especie"texto2 <- "texto que incluye OS= en la mitad"grepl("OS=", texto1)
## [1] FALSE
grepl("OS=", texto2)
## [1] TRUE
Entonces podemos aprovechar esta función para hacer esto:
desc_os <- grepl("OS=", ids_crudos$Description)table(desc_os)
## desc_os ## TRUE ## 257
Bien! Todas las líneas de descripción tienen una clave “OS”. Ahora hay que extraerlas todas, junto con las dos palabras que siguen, que son el género y especie.
Ni grep() ni grepl() nos van a servir aqui. Pero hay otra función que si, regexpr().
Empecemos probando con la primera línea:
regexpr("OS=", ids_crudos$Description[1])
## [1] 24 ## attr(,"match.length") ## [1] 3 ## attr(,"useBytes") ## [1] TRUE
match_1 <- regexpr("OS=", ids_crudos$Description[1])str(match_1)
## atomic [1:1] 24 ## - attr(*, "match.length")= int 3 ## - attr(*, "useBytes")= logi TRUE
match_1 es un vector de tamaño 1 que indica donde empieza la coincidencia y contiene dos atributos extra. Nos interesa el primero,match.length. Veamos:
largo_match_1 <- attr(match_1, "match.length")substr(ids_crudos$Description[1], start = match_1, stop = match_1 + largo_match_1 - 1)
## [1] "OS="
Funciona, pero nos están faltando las cadenas de caracteres de género y especie.
match_2 <- regexpr("OS=\\w*\\s\\w*", ids_crudos$Description[1])largo_match_2 <- attr(match_2, "match.length")substr(ids_crudos$Description[1], start = match_2, stop = match_2 + largo_match_2 - 1)
## [1] "OS=Bubalus bubalis"
Ya tenemos casi todo lo que precisamos. En el código que sigue aprovechamos la capacidad e R para trabajar con vectores y procesamos todas las líneas de descripción juntas.
match_os <- regexpr("OS=\\w*\\s\\w*", ids_crudos$Description)largo_match_os <- attr(match_os, "match.length")os <- substr(ids_crudos$Description, match_os, match_os + largo_match_os - 1)head(os)
## [1] "OS=Bubalus bubalis" "OS=Bubalus bubalis" "OS=Ovis aries" ## [4] "OS=Bos mutus" "OS=Capra hircus" "OS=Bos taurus"
Solo nos resta sacar el texto “OS=” de los strings y armar una tabla.
os <- sub("OS=", "", os)length(os)
## [1] 257
length(ids_crudos$Description)
## [1] 257
table(os)
## os ## Bos indicus Bos mutus Bos taurus ## 2 33 65 ## Bubalus bubalis Capra hircus Cephalophus dorsalis ## 6 31 1 ## Damaliscus pygargus Nanger dama Oreotragus oreotragus ## 1 1 1 ## Ovis aries Ovis vignei ## 115 1
# O más prolijoas.data.frame( table(os) )
## os Freq ## 1 Bos indicus 2 ## 2 Bos mutus 33 ## 3 Bos taurus 65 ## 4 Bubalus bubalis 6 ## 5 Capra hircus 31 ## 6 Cephalophus dorsalis 1 ## 7 Damaliscus pygargus 1 ## 8 Nanger dama 1 ## 9 Oreotragus oreotragus 1 ## 10 Ovis aries 115 ## 11 Ovis vignei 1
Ahora ya sabemos a qué especies y con que abundancia pertenecen las proteínas que resultaron más similares a las que se obtuvieron en el suero de cabras.
Ademas de saber el orígen y la descripción de las proteínas que resultaron más similares a las encontradas en cabra, nos podría interesar disponer de las secuencias fasta de todas ellas.
Uniprot utiliza unas reglas relativamente sencillas para generar URLs con las cuales recuperar accesos específicos y en formatos predeterminados. Estas reglas se pueden consultar aqui.
Por ejemplo, para descargar una secuencias en formato FASTA el URL genérico es:
http://www.uniprot.org/uniprot/[acc_id].fasta
Donde acc_id es el número de acceso de la secuencia que queremos descargar. Nosotros tenemos que descargar 257, por lo que necesitamos implementar alguna estrategia para obtener estas secuencias sin tener que hacerlo de una por vez.
Pero vamos a ir por partes. Primero vamos a construir una función que descargue una única secuencia. Como argumento tenemos que pasarle el número de acceso.
leer.fasta <- function(num.acceso){# Esta es una función para descargar desde Uniprot una secuencia# en formato fasta a partir del número de acceso base <- "http://www.uniprot.org/uniprot/" tipo.fasta <- ".fasta" secuencia.fasta <- paste0(num.acceso, tipo.fasta) secuencia.url <- paste0(base, secuencia.fasta) download.file( secuencia.url, secuencia.fasta, quiet = T)}
Probemos esta función:
leer.fasta(ids_crudos$Accession[29])
Atención: Esta función va a generar un mensaje de error y crear un archivo fasta de cero bytes si el número de acceso no existe. Más adelante vamos a ver cómo manejar estas situaciones. Para este análisis, esto no es un problema.
Ahora necesitamos construir un loop para descargar todas las secuencias y guardarlas en un solo archivo, simil.cabras.fasta.
file.create("simil.cabras.fasta")
## [1] TRUE
for(accs in ids_crudos$Accession) { leer.fasta( accs ) # print(paste("leyendo:",accs)) file.append("simil.cabras.fasta", paste0(accs,".fasta") ) file.remove(paste0(accs,".fasta"))}
Ya tenemos las 257 secuencias fasta de las proteínas de cabra o de otras especies emparentadas con la cabra y que son similares a las detectadas en el suero de leche de cabra.