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.)
Eso es todo! buena vibra!!! :)
Por Susana Ortiz-Baez.
No hay comentarios:
Publicar un comentario