abc.r

abc.R

# Script de R: Introduccion a R
# Goku

# Solamente si estoy en la FACEN
Sys.setenv(http_proxy="http://laboratorio:laboratorio2017@172.16.1.1:3128/")
Sys.setenv(https_proxy="https://laboratorio:laboratorio2017@172.16.1.1:3128/")

## Trabajando con Secuencias de DNA

# Una de las cosas mas basicas qeu necesitas saber hacer con R para bioinformatica es
#leer y escribir secuencias de DNA o proteinas:
  #Necesitaras leer las secuencias antes de hacer cualquier analisis

# Puedes estar trabajando con secuencias ya establecidas en las bases de datos,
#o con nuevas secuencias tuyas.

# Hay algunos formatos básicos y estandares para secuencias, ej.:
  #Single sequence formats: FASTA, ASN.1, GCG, etc.
  #Multiple sequence formats: FASTA, Clustal, MSF, etc.

#Por lejos el formato mas simple y comun es FASTA, que, como puedes ver, es compatible 
#secuencias simples y secuencias múltiples -
#así que comencemos aprendiendo cómo importar archivos FASTA en R.

# El formato FASTA es muy simple: es un archivo de texto plano donde la primera línea comienza con ">" (la
#linea de descripcion) y contiene cualquier nombre o identificador que pertenezca a la secuencia.

# La (s) siguiente (s) línea (s) son la secuencia misma, ADN o proteína, en el código IUPAC estándar.

#Si el archivo FASTA contiene más de una secuencia, la siguiente secuencia se indica con una nueva
# línea de descripción: esto se puede repetir para tantas secuencias como se desee.

# Instalaciones (instalar solo si se usa fuera de biolinux.ova)
#install.packages("seqinr")
#install.packages("ape")
#If R version error (ape or others...)
#http://scottsfarley.com/research/cloudcomputing/2016/07/19/Updating-R-on-Debian.html
#https://stackoverflow.com/questions/10255082/installing-r-from-cran-ubuntu-repository-no-public-key-error
#install.packages("rentrez")
  #Si da error instalar en terminal linux:
    #sudo apt-get -y build-dep libcurl4-gnutls-dev
    #sudo apt-get -y install libcurl4-gnutls-dev
    #sudo apt-get install libxml2-dev
#install.packages("XML")
#install.packages("rentrez")
  #Si da error instalar en consola linux:
    #sudo apt-get install lib32z1-dev
    #sudo apt-get install libxml2-dev
    #sudo apt-get install libssl-dev
    #sudo apt-get install libcurl4-nss-dev


# Librerias
library(seqinr)
library(ape)
library(rentrez)
library(ggplot2)


# Lectura de un archivo local FASTA
# Empecemos leyendo un archivo FASTA en R usando seqinr library.

# Ahora que le hemos dicho a R que vamos a usar la biblioteca seqinr, 
#podemos escribir un pequeño fragmento de código para cargar un archivo FASTA.

# Descargar archivo
#https://prod-edxapp.edx-cdn.org/assets/courseware/v1/6ebb427f5f53343cbb7af97945366d02/asset-v1:USMx+BIF003x+1T2018+type@asset+block/cox1.fasta

cox1 <- read.fasta(file = "cox1.fasta", seqtype = "AA")
# Si todo salio bien, ahora deberías tener una variable llamada "cox1" en la ventana de su entorno.

# Abra esa variable haciendo clic en la flecha a la izquierda del nombre y le dirá que tiene 4
#objetos (secuencias) en esta variable, y al desplazarse por la lista, le dirá que son de tipo
  #SeqFastaAA (Amino Acid).

# Si hubiéramos omitido la parte "seqtype =" AA "" del código, read.fasta
#hubiera dejado de forma predeterminada a DNA.
#note que no es lo suficientemente "inteligente" como para detectar 
#automáticamente el tipo de secuencia; debes especificarla.

# Podemos dividir las cuatro secuencias en secuencias individuales también.

# Escribe en la consola
length(cox1) 
# y ejecutar - le dirá que hay cuatro secuencias en la variable.

# Para ver solo el primero podes escribir 
cox1[1]
# y ejecutar, en tu script podrias decir 
seq1 <- cox1[1] 
#– hagamos eso ahora

## Recuperar una secuencia de genbank como un objeto binario

# Antes de ir mucho más lejos, aprendamos cómo importar 
# una secuencia directamente desde una base de datos en lugar de desde un archivo local
# Esta vez vamos a recuperar una secuencia de DNA de Genbank. 

# Tenga en cuenta que hay muchos paquetes diferentes que nos permiten conectarnos 
# a bases de datos: APE es solo uno de los más fáciles.

#Vamos a agregar esta linea a tu script:
AB003468 <- read.GenBank("AB003468", as.character = "TRUE")

# Si ejecuta esto, verá aparecer una nueva variable (AB003468) en su entorno: esta es la
  # secuencia que acaba de descargar de Genbank.

# Abrir la variable haciendo clic en la flecha a la izquierda de su nombre le mostrará los atributos, es
  # un vector de clonación.

# Este archivo está actualmente en formato binario, necesitamos guardarlo como FASTA y
# leerlo antes de poder hacer cualquier estadística.

##Guardar una secuencia de genbank en formato FASTA
# Ahora podemos hacer eso usando este comando de APE y agregándolo a su script, luego ejecútelo:
write.dna(AB003468, file ="AB003468.fasta", format = "fasta", append = FALSE, nbcol = 6, colsep = " ", colw = 10)

# Nota IMPORTANTE:
# Asegúrese de revisar su directorio de trabajo abriendo el archivo haciendo clic en el panel ARCHIVOS para ver cómo se ve.
# Note cómo los diversos comandos de formato hacen que el archivo se vea, 
#ej 6 columnas, con un ancho de columna de 10, separadas por un espacio.

# También podemos recuperar secuencias de Genbank usando un paquete creado por NCBI llamado rentrez.

# El paquete rentrez nos permite recuperar secuencias de otras formas además del número de acceso
  # y no nos limita al formato binario que hace APE.

# Por ejemplo, podes buscar secuencias usando entrez_search:
entrez_search(db="nucleotide", term="human superoxide dismutase")

# Probemos eso en su consola y también puede recuperar toda la lista de ID que devuelve esta búsqueda,
# si lo desea, luego pasar esos ID a entrez_fetch; por ahora, usemos el comando APE ya que es más fácil.

# Nota IMPORTANTE:
# Para obtener más información sobre rentrez y acceder a un tutorial completo, 
#puede hacer clic en el enlace de rentrez.
#https://cran.r-project.org/web/packages/rentrez/vignettes/rentrez_tutorial.html


# Strip first sequence of AB003468 into just characters
# Primero necesitaremos convertir nuestra secuencia en una cadena simple para el análisis.
# Usemos la secuencia AB[[1]]003468.
CloningVector <- AB003468

#Esto le dice a seqinr que quite la primera (y única) secuencia
# de la variable AB003468 a los caracteres básicos (sin encabezado).

# #Haga un recuento de nucleótidos basico en CloningVector
# Ahora hagamos un simple conteo de nucleótidos y ejecutemos:
count <- count(CloningVector,1) 

# Esto recupera el recuento total de cada nucleótido (A, C, G y T) y lo almacena en el recuento de variables.
# Vamos a escribir
count
# en la consola y ejecutemos para ver que pasa.

# El comando de conteo requiere el objeto en este caso, CloningVector, así como un valor numérico -
# en este caso le dimos 1, lo que resultó en sequin que nos muestra todas las bases posibles.

# Podemos ver fácilmente todas las combinaciones de dinucleótidos usando combinaciones 
#de 2 o trinucleótidos usando un 3:count(CloningVector,2) #or 
count(CloningVector,3)

# Nota IMPORTANTE:
#Esto se llama análisis de "conteo de palabras", y se puede hacer para cualquier tamaño de "palabra".
# Otra métrica importante que se puede obtener fácilmente es el contenido de GC de una secuencia.
# El contenido de GC nos dice cuál es la proporción de residuos de G (guaninas) y C (citosinas) en comparación con A
# (adenina) y residuos de T (timidina) en una secuencia

# Esto es importante porque las regiones de codificación tienden a ser más altas en GC.

#El contenido de CG también afecta la temperatura de "melting" del ADN.
# Calcular el contenido de GC de CloningVector
GC <- GC(CloningVector)

#GC aparece ahora en el panel de Enviroment para CloningVector, el contenido de GC es 0.51,
# lo que significa que la relación GC a AT es bastante parejo en esta secuencia.

# Sin embargo, el valor de GC es un valor global, y si el contenido de GC varía en el transcurso de una secuencia 
# Si estuviera analizando un genoma completo, sería
# mucho más interesante ver cómo el contenido de GC cambia a lo largo de una secuencia.

# Para hacer eso tendremos que hacer algo un poco más sofisticado, calculando el contenido del GC por
  # "ventanas", o regiones cortas, de la secuencia.

#Luego podemos trazar el cambio en el contenido de GC sobre estas ventanas y ver si hay algún cambio interesante.

# Comencemos por calcular el contenido del GC en un tamaño de ventana de 200.
# Calcular ventanas de tamaño 200 para nuestra secuencia
# Usaremos una variable llamada GCwindow para dividir la secuencia en trozos de 200.
GCwindow <- seq(1, length(CloningVector)-200, by = 200)

# Esto usa el comando seq para encontrar los "puntos de quiebre" de CloningVector en trozos de 200

# Si escribis
GCwindow 
# en la consola y lo ejecutas veras los valores 1, 201, 401, etc..

# Estos son los puntos de inicio de ventana para nuestro análisis: 1 a 200, 201 a 400, etc.

# Nota IMPORTANTE:
# Asegúrese de comprender por qué tomamos la longitud de CloningVector y restamos 200 de ella.
# Si no tiene sentido, intente ejecutar este código sin el -200 y vea qué valores obtiene en su lugar
  # y comprobar si tienen sentido.


#Encuentra la longitud del vector GCwindow
# Ahora que tenemos nuestros puntos de inicio de ventana, podemos usarlos en un ciclo FOR para "dividir" 
#la secuencia. Primero, necesitamos saber con cuántos trozos estamos tratando agregando esta línea:
n <- length(GCwindow)

# Ahora creamos un nuevo vector en blanco con el mismo número de espacios en blanco que podemos 
#llenar con nuestros valores de GC:
Chunks <- numeric(n)

#El numérico asegura que obtenemos el número bruto para n.

# Finalmente podemos agregar el ciclo FOR.

# Primero, hagamos que imprima el contenido del GC en la consola para asegurarnos de tener el ciclo correcto.

# FOR loop para calcular GC por fragmento
# El loop FOR consta de estas líneas.
for (i in 1:n) {
  chunk <- CloningVector[GCwindow[i]:(GCwindow[i]+199)]
  chunkGC <- GC(chunk)
  print(chunkGC)
  Chunks[i] <- chunkGC
}

# Debería ver que la consola enumera una serie de valores de GC, uno para cada fragmento 
#(ventana) en nuestro análisis.

# Pasemos por esas líneas una a la vez antes de dar el último paso.

# La primera línea es fácil - para (i en 1: n) {- esto inicia el ciclo FOR y se asegura de que 
# se repita tantas veces como n, que es la cantidad de fragmentos o ventanas que tenemos.

# La segunda línea - fragmento <- CloningVector [GCwindow [i] :( GCwindow [i] +199)] - es crítica.

# Esta línea toma CloningVector por una longitud de i a i + 199 y subconjuntos de la secuencia, 
# estableciendo el valor del fragmento en esas bases:
  
# Para la primera iteración del ciclo FOR, i es 1, por lo que significa que el fragmento se establece 
# en las primeras 200 bases de CloningVector.

# Para la próxima iteración, i es 2, que va al segundo valor en GCwindow, que es 201, por lo que 
# las bases 201 a (201 + 199) o 400 se configuran en bloque.

# El fragmento de línea GC <- GC (fragmento) toma el valor actual del fragmento y calcula el valor del GC.

# La siguiente línea simplemente imprime el valor actual de GC (chunkGC).

# Luego, el "}" le dice al ciclo que vuelva a la parte superior y haga la siguiente iteración.

# El resto de este análisis es fácil. Agreguemos una línea que llene nuestra variable "en blanco", Chunks, 
# con el valor apropiado de chunkGC en cada paso.

# Como i se usa para representar cada contador de pasos, también podemos usar eso para decirle a R qué valor en Chunks completar.
# Agregar esta línea después del comando print (chunkGC):
#Chunks[i] <- chunkGC

# Ahora ejecute el código y los valores seguirán imprimiéndose en la consola, 
# pero también se agregarán al vector Chunks.
# Veámoslo en el Enviroment o escribe
Chunks 
#en la consola y ejecutalo.

# Plotting our GC Windows
# El último paso de este análisis particular será trazar el resultado.

# Al trazar el resultado, podemos ver gráficamente cómo el contenido del GC cambia en el espacio de la secuencia.

# Ya tenemos todos los datos que necesitamos, solo tenemos que construir el comando plot - 
# vamos a mantenerlo simple y solo usar las funciones integradas en R.
plot(GCwindow,Chunks,type="b",xlab="Nucleotide start position",ylab="GC content")

# Esto le dice a R que genere un diagrama usando GCwindow como nuestro eje X (recuerde que las ventanas son los tramos de 200 pb) y 
# los trozos como nuestro eje Y (donde los trozos son el valor GC real).

# El "tipo" le dice a R que use puntos conectados en lugar de un simple diagrama de dispersión.

# Cuando ejecutas esta línea, deberías ver que la trama aparece en tu panel PLOTS.

# Observe cómo los ejes X e Y tienen etiquetas agradables, ésos fueron generados por las partes xlab y ylab del comando.

# Este resultado te dice que, si bien el contenido del GC es 0.5 en general, cambia drásticamente 
# a lo largo de la secuencia si usamos una ventana de 200 pb.

# Si hubiéramos usado un tamaño de ventana diferente, ¿el gráfico habría sido aún más diferente?

# La respuesta es sí"!

# Sería útil mirar los mismos datos en diferentes tamaños de ventana, p. 100 o 50 (o quizás ventanas más grandes).

# Parece que podría ser "más fácil" simplemente copiar y pegar el código y volver a escribirlo para diferentes tamaños de ventana.

# Esto sería ineficiente, además el código sería menos reutilizable, digamos que quieres analizar 
# una secuencia diferente en el futuro, tal vez quieras tamaños de ventana completamente diferentes para esa, 
# si simplemente editaras el código, terminarías editando una y otra vez para cada análisis.

# Lo más inteligente es crear una pieza de código reutilizable.
# Hay varias formas de hacerlo, pero la más útil es escribir una función R personalizada.

# Una función personalizada en R es una forma de definir un fragmento de código que luego puede invocar,
# al igual que puede invocar cualquier otra función R.

# La mayor diferencia es que sus funciones personalizadas no estarán "integradas", lo que significa que cada vez 
# que quiera usar una función personalizada deberá incluirla en el script R que lo llame.

# Esto sigue siendo Good Thing™ - es fácil copiar y pegar, o almacenar una "biblioteca" de funciones en un archivo de texto.

# Si su función es realmente útil, incluso podría enviarla como un paquete que 
# otros usuarios podrían descargar e instalar.

# Vamos a escribir una función personalizada en R que calculará y trazará contenido de GC dado solo una secuencia de entrada 
# y el tamaño de ventana deseado.

# Función personalizada para el trazado de ventanas de GC
slidingwindowGCplot <- function(windowsize, inputseq) #Esta línea le dice a R que está definiendo una función personalizada 
  #llamada slidingwindowGCplot que aceptará dos entradas: lo primero que se le pasará será el tamaño 
  #de la ventana (windowsize), y la segunda será la secuencia real para analizar (inputseq) .
{# En la siguiente línea agrega un "corchete abierto" - {
  GCwindow <- seq(1, length(inputseq)-windowsize, by = windowsize) # Observe cómo hemos sustituido 
  # la variable windowsize (una de nuestras entradas a la función) por las 200 que teníamos 
  # en la secuencia de comandos "codificada" e inputseq para la secuencia que estaba
  # codificada previamente (CloningVector)
  
  # Encontrar la longitud de la ventana de GC
  n <- length(GCwindow) # Hacer un vector de la misma longitud, pero "en blanco" para que podamos llenar
  Chunks <- numeric(n) #Esas dos líneas a continuación no necesitaban modificaciones para 
  #incluirlas en nuestra función.
  for (i in 1:n) {
    chunk <- inputseq[GCwindow[i]:(GCwindow[i]+windowsize-1)] # Reemplazamos los valores codificados para la secuencia y el tamaño de la ventana con nuestras variables
    chunkGC <- GC(chunk)
    print(chunkGC)
    Chunks[i] <- chunkGC
  }
  plot(GCwindow, Chunks, type="b", xlab = "Nucleotide start position", ylab = "GC content")
} 

# Llamadas de diferentes tamanos de ventana
slidingwindowGCplot(200,CloningVector) #Esto le dice a R que use su función personalizada slidingwindowGCplot con un windowsize de 200 y CloningVector como la secuencia.
slidingwindowGCplot(175,CloningVector)

#Set slliding windowsize
slidingwindowGCplot(200,CloningVector)
#Custom Function for GC Window Plotting w/ main=paste
slidingwindowGCplot <- function(windowsize, inputseq)
{
  GCwindow <- seq(1, length(inputseq)-windowsize, by = windowsize)
  #Encontrar la longitud de GCwindow
  n <- length(GCwindow)
  #Haga un vector de la misma longitud, pero "en blanco" para que lo rellenemos
  Chunks <- numeric(n)
  for (i in 1:n) {
    chunk <- inputseq[GCwindow[i]:(GCwindow[i]+windowsize-1)]
    chunkGC <- GC(chunk)
    print(chunkGC)
    Chunks[i] <- chunkGC
  }
  plot(GCwindow, Chunks, type="b", xlab = "Nucleotide start position", ylab = "GC content", main = paste("GC Plot with windowsize ", windowsize))
}

# El comando pegar nos permite combinar varios elementos en la misma línea de título,
# pero tenemos que combinar los elementos entre paréntesis, por ejemplo: main = pegar ("título", variable)

# Asegúrese de que esto tenga sentido para usted antes de continuar.

# En el futuro, si desea utilizar esta función personalizada de GC, puede copiar y pegar 
# la función en un nuevo guión R y luego invocarlo en el guión mismo o en la consola escribiendo 
# el nombre de la función y pasando el parámetros. Por ejemplo: slidingwindowGCplot (200, CloningVector).