miércoles, 23 de febrero de 2011

Descarga de secuencias Nucleotidicas de GenBank desde R

El post anterior manejaba distancias geneticas desde secuencias nucleotidicas; pero como y donde se obtienen estas secuencias?:


Una forma de obtener las secuencias es bajandolas una por una! desde el GenBank. Pero esto definitivamente es super tedioso!!!, ya que de hecho navegar por la pagina del GenBank es complicado como tal sí uno no sabe como buscar y en que enlaces entrar, y cuando uno necesita un gran numero de secuencias lleva mucho pero mucho tiempo decargarlas una por una!!!
Gracias a (Dios, Alá, Yahveh o Jehová, Zeus, Achamán, Ahura Mazda, Kamisama, Cao Dái, Waheguru, Shang Di, Elohim, cualquiera que sea!! o definitivamente a Emmanuel Paradis) existe ape en R y nos permite bajar secuencias del GenBank mucho mas sencillo y de una forma muy rapida por muchas que sean las secuencias a bajar... 


Para ello lo primero que debemos hacer es construir una matriz de formato .csv que contenga el codigo de cada secuencia que vayamos a bajar, los nombres de todos los taxa y el nombre del gen de las secuencias que vayamos a utilizar!!: (Para este ejemplo utilizaremos una matriz en formato .csv con 29 Especies de Canidos representativos de todo el mundo y  las secuencias que se bajaran pertenecen al mismo gen (TRSP) para todo los taxa).
Entonces la apariencia de la matriz de datos construida en formato .csv debe der asi:



"cod","taxon","gen"
AY609086,Alopex_lagopus,TRSP
AY609087,Atelocynus_microtis,TRSP
AY609088,Canis_adustus,TRSP
AY609089,Cuon_alpinus,TRSP
AY609091,Canis_aureus,TRSP
AY609092,Chrysocyon_brachyurus,TRSP
AY609094,Canis_latrans,TRSP
AY609095,Canis_lupus,TRSP
AY609097,Canis_mesomelas,TRSP
AY609098,Cerdocyon_thous,TRSP
AY609084,Canis_familiaris,TRSP
AY609101,Pseudalopex_griseus,TRSP
AY609104,Pseudalopex_gymnocercus,TRSP
AY609100,Pseudalopex_sechurae,TRSP
AY609106,Fennecus_zerda,TRSP
AY609108,Lycaon_pictus,TRSP
AY609109,Nycteruetes_procyonoides,TRSP
AY609112,Otocyon_megalotis,TRSP
AY609115,Speothos_venaticus,TRSP
AY609117,Urocyon_cinereoargenteus,TRSP
AY609120,Vulpes_corsac,TRSP
AY609121,Vulpes_macrotis,TRSP
AY609122,Vulpes_vulpes,TRSP
AY609123,Ursus_americanus,TRSP
DQ093097,Ailuropoda_melanoleuca,TRSP
DQ093098,Mirounga_augustirostris,TRSP
DQ093099,Odobenus_rosmarus,TRSP
DQ093100,Enhydra_lutris,TRSP
DQ093101,Lontra_longicaudis,TRSP
DQ093102,Procyon_lotor,TRSP

Donde la primera fila muestra que la tabla contiene tres columnas con sus respectivos nombres: "cod" (codigo), "taxon" y "gen".
Despues de tener el archivo formato.csv llamado "canidos.csv"que acabamos de crear y ponerlo en el directorio sobre el cual estamos trabajando en R, Iniciamos con la descarga de secuencias:


>library(ape)
#(Suponiendo que ya hemos descargado el paquete "ape" se debe llamar el paquete a R con el comando anterior.)


>berraco <- read.csv("canidos.csv",header=T)
# (Con este comando leemos la tabla creada anteriormente y usamos el vector de codigos para hacer del objeto secuencias, con el argumento "header=T" o "header=TRUE" le indicamos que la primera fila de nuestra tabla son los nombre de cada columna, de lo contrario utilizamos "header=F" o "header=FALSE")


>seqcanidae <- read.GenBank(berraco$cod)
#(leemos los codigos del GenBank del objeto antes creado "berraco" y leemos la columna "cod")


>nombres <-mat.or.vec(length(berraco$taxon),1)
#(Tomamos los nombres del objeto "berraco" y consturimos un vector de nombres para las secuencias de los codigos que acabamos de leer)

>nombres <-paste(berraco$taxon)
#(Unimos los nombres que esten espaciados)

>names(seqcanidae) <- nombres
#(Se le asiganan los nombres a las secuencias del objeto "seqcanidae" con el vector de nombres que creamos)

>write.dna(seqcanidae, "canidae.fas", format = "fasta")
#(Bajamos las secuencias finalmente y las escribimos en formato fasta con el nombre de "canidae.fas", este archivo quedara guardado en el directorio sobre el cual estemos trabajando)


>write.nexus.data(seqcanidae, "canidae.nex", format = "nexus")
#(Tambien podemos escribir las secuencias en otor tipo de formatos como el formato "nexus" pero esto ya se hace depende de los requerimientos del usuario. Otro tipo de formatos mirar el manual de ape.)

No siendo mas....
Hasta la proxima....

miércoles, 9 de febrero de 2011

Cálculo de distancias nucleotídicas en R.

Como preámbulo nos entusiasma comunicarles que recientemente R para chichombianos ha sido aceptado como parte de la comunidad de R blogglers!!! Para quienes no lo conocen, es una central de blogs que incorpora información publicada exclusivamente sobre programación en R acerca de infinidad de tópicos. Por tanto los invitamos a visitar la pagina oficial y unirse a esta comunidad de entusiastas y locos por R! http://www.r-bloggers.com/
Para variar de tema, y dado que el objetivo de este blog es abordar las bondades de utilizar R y sus aplicaciones en la solución de múltiples necesidades que se presentan en diferentes áreas del conocimiento, en esta ocasión el turno es para los interesados en filogenia/evolución/bioinformática y la temática a tratar es sobre el "cálculo de de distancias desde secuencias nucleotídicas" es un script corto y sencillo, asi que la transicion sera "suave".
La definición general de distancia nucleotídica hace referencia al grado de diferenciación en la composición nucleotídica entre dos unidades (secuencias moleculares) y puede reflejar el grado de similaridad ó divergencia genética actual, en donde pequeñas distancias genéticas podrían indicar relaciones cercanas mientras grandes distancias señalan relaciones distantes. Ver wiki: http://es.wikipedia.org/wiki/Distancia_genética
Aunque se ve simple y bonito, la interpretación de las distancias genéticas debe hacerse con cautela dado que puede verse ensombrecida por diversos factores (reversiones, eventos de recombinacion, definición de un modelo evolutivo, cuantificación del error en procesos de amplificación y secuenciación, etc.). Bajo este contexto, tal vez sea mejor ver el calculo de estas distancias como una herramienta exploratoria de los datos la cual debe complementarse con análisis mas profundos.
Iniciemos...

>install.packages ("ape")
# (Generalmente, ape es uno de los paquetes por defecto incorporados en R pero de no ser el caso es posible descargarlo visitando http://cran.r-project.org/web/packages/ape/index.html). Allí podrán encontrar los archivos comprimidos para los sistemas operativos más comunes, además de toda la información referente al paquete, incluyendo el manual).
>library(ape)
# Lo siguiente que se debe hacer es cargar el paquete ape (Analyses of Phylogenetics and Evolution), sí no lo conocen es un paquete realmente útil que provee funciones para leer, escribir, graficar y manipular filogenias, realizar análisis comparativos, analisis de diversificación y macroevolución, computar distancias entre muchas otras aplicaciones.)
>DNA <- read.dna("Genc.fas", format="fasta", skip=0, nlines = 0, comment.char="#", seq.names = NULL, as.character = FALSE, as.matrix = NULL)

# (Obviamente es necesario leer los datos, en este caso con la función "read.dna" desarrollada por Emmanuel Paradis. Con el argumento "Genc.fas" me refiero al nombre del archivo que contiene las secuencias nucleotidicas alineadas)

# (El formato puede ser "interleaved" o "sequential", es decir, sí las secuencias se organizan en bloques/múltiples lineas de x tamaño o sí son continuas siguiendo una línea.También es posible especificar formatos "clustal", o "fasta")
# (Asignandole un valor de cero al argumento "skip" le estoy señalando que lea el archivo desde el principio, sin omitir nunguna línea. Es necesario, especificar que los comentarios estan identificados con el simbolo #)
# (Con el argumento seq.names los nombres son especificados a partir del archivo leido y con "as.character" las secuencias son recuperadas como un objeto de clase "DNAbin" el cual permite que estas sean manipulables.)
# (Sí el formato elegido es el fasta y sí estas tienen la misma longitud, tengo la opción de recobrar las secuencias en una matriz con el argumento "as.matrix". Para este caso, no se utilizará por ende sera "NULL")
### Para el cálculo de las distancias nucleotídicas, es posible especificar o no un modelo de evolutivo:  
>DIST <- dist.dna(DNA, model = "raw", variance = TRUE, gamma = FALSE, pairwise.deletion = FALSE, base.freq = NULL, as.matrix = FALSE)
# (La función dist.dna computa una matriz de distancias por pares a partir del vector "DNA" previamente creado, el cual contiene las secuencias, en este ejemplo no se especifica un modelo evolutivo, pero es posible definir un modelo bajo el cual se asume que evolucionan las secuencias (ver manual APE) 
# (El argumento "variance = TRUE" señala que se computaran las varianzas de las distancias. Es posible además, asignar un valor al parámetro gamma para hacer una corrección sobre las distancias, si el argumento es "gamma = FALSE" no se aplicará la corrección) 
# (Sí existen sitios con datos perdidos en la comparación de las secuencias es posible eliminarlos o no, con el argumento "pairwise.deletion = TRUE/FALSE". Asi mismo es factible definir con el argumento "base.freq" si hay una proporción especifica en la frecuencia de bases (generalmente de acuerdo al modelo especificado), o si estas deben ser estimadas a partir de los datos.
# (Con "as.matrix" recupero los resultados como una matriz (TRUE), o como un objeto de clase dist (FALSE), ver manual.)

### Computación de distancias genéticas:

>GEN <- dist.gene(denv1.fas, method = "pairwise", pairwise.deletion = FALSE, variance = TRUE)

# (Esta función es mas general que la anterior y acepta diferentes tipos de datos como alelos, haplotipos, SNP y secuencias de DNA. La ventaja de este método es que el usuario puede especificar el método que quiere usar para computar las distancias, esto es comparación por pares "pairwise", a partir del numero de loci en que se diferencian dos secuencias o de acuerdo a un porcentaje "percentage" donde la distancia se divide en el número de loci. En todos los casos es posible estimar la varianza asociada.)

### Con esta función se pueden indicar los sitios polimórficos en el set de datos y de esta forma puedo inferir grosso modo cuantos sitios son los que posiblemente generan las diferencias entre las secuencias, y por tanto las distancias.

>y <- seg.sites("denv1.fas")
>y
>length(y)

### Adicionalmente, puedo obtener un resúmen de las distancias genéticas y graficarlas, en este caso sin generar una salida

>summary(DIST)
>plot(DIST)

### Sí deseo puedo generar una tabla que contenga las distancias calculadas en el vector DIST
>write.table (DIST, file = "distancias.csv", append = TRUE, sep = " ", na = "NA", dec = ".", row.names = TRUE, col.names = TRUE)

# (En este caso, el formato de salida es .csv en el que las columnas se separan por comas y las filas por saltos de línea.)

# (Si hay stios perdidos en los datos estos son denominados como NA y los valores decimales se indican por medio de punto.)

# (Con "append = true" los datos son anexados en el archivo y sep = "" las columnas son separadas por comillas)

#(Finalmente, si se desea asignar como encabezado nombre a las filas y a las columnas, con row.names y col.names = TRUE es posible, de esta forma se obtiene una tabla tabulada y ordenada en un archivo de texto perfecto para ser importado en una hoja de cálculo.)
# ## Para mayor informacion acerca de los argumentos que pueden defirnirse, visite: http://stat.ethz.ch/R-manual/R-patched/library/utils/html/write.table.html
Eso es todo! buena vibra!!! :)
 
Por Susana Ortiz-Baez.