Guía de ejercicios básicos de R

HOLA MUNDO

print("Hola MUNDO")

para obtener ayuda de cualquier función usamos

?print

lo podemos usar como calculadora

1+1+2+3+5+8
8-13
13*21
34*34
55/34 * sqrt(89)

Para guardar alguno de estos resultados necesitamos asignarlos a un “objeto”. Para esto podemos usar el símbolo “=” o “<-”

x = 2+1 #x <- 2+1
x

Se pueden hacer operaciones directamente con los objetos

otro_numero = 55/34 * sqrt(89)
otro_numero
otro_numero*x

para eliminar un objeto usamos rm

rm(otro_numero)
otro_numero

Para crear un objeto con más de un elemento usamos c()

Ventas = c(1,1,2,3,5,8,13)
Ventas

R posee un método recursivo que nos permite realizar una operación sobre el objeto múltiples veces

Ventas + 5

También podemos usar R con funciones. Una función es una palabra seguida por paréntesis word()

sum(Ventas)
mean(Ventas)
max(Ventas)
min(Ventas)
range(Ventas)
unique(Ventas)

Para ordenar usamos sort()

sort_increase=sort(Ventas)

OBJETOS Y VECTORES

frutas = c("manzana" , "banana" , "sandía" , "ananá" , "limón" , "frutilla" , "naranja")

para asignar nombres a elementos en un objeto existe la función names()

names(Ventas) = frutas
Ventas

¿Qué pasó?

El objeto Ventas es un vector. El vector posee múltiples elementos en una determinada posición que se indica mediante el empleo de los corchetes: Ventas[posición]

Ventas[3]

Múltiples series de elementos pueden ser obtenidos al mismo tiempo. En este caso necesitamos usar c()

Ventas[c(1,3:7)]

Y si por el contrario, queremos eliminar esos elementos:

Ventas[-c(1,5)]

EJERCICIO: ¿CUál es la media de frutas vendidas sin contar “frutilla”? Utilizar lo aprendido. El resultado debe ser 4.166667

OPERADORES RELACIONALES

A veces queremos saber características de los números. Por ejemplo: El número de tiburones con longitud mayor a 2 metros. Esto se puede responder con operadores relacionales (más grande que, más largo que, igual a, etc.)

Con el objeto ventas. ¿Qué frutas se vendieron más de 5 veces?

Ventas>5
Ventas == 5

En el caso anterior usamos “==” para evitar asignar 5 al objeto Ventas. Como estas comparaciones nos dan vectores lógicos TRUE/FALSE o 1/0. Esto puede resultar muy útil para hacer algunos cálculos:

sum(Ventas > 5)
Ventas[Ventas > 5]

EJERCICIO: ¿Cómo sabemos que fruta se vendió más caro que la media?

AHORA VAMOS A HACER GRÁFICOS!!!

R es muy popular por la calidad y variedad de gráficas que puede generar. Una vez que los datos están en un objeto, graficarlo es muy simple. Por ejemplo: grafiquemos las frutas vendidas

plot(Ventas)

Si queremos otro tipo de gráficas

barplot(Ventas)
pie(Ventas)
hist(Ventas)

Podemos ver todas las opciones para cada tipo de gráfico usando “?”, pero sigamos con ejemplos

hist(Ventas, col="skyblue", main="Histograma en azúl", ylab="Frutas vendidas")

Existen varias maneras de guardar gráficas. Por ejemplo: como pdf, pero primero definimos la carpeta de trabajo

getwd()
setwd("./curso/")
pdf("histogram.pdf")
hist(Ventas, col="skyblue", main="Histograma en azúl", ylab="Frutas vendidas")
dev.off()

LECTURA DE ARCHIVOS Y MANIPULACIÓN DE TABLAS

Cuando ya tenemos valores en un archivo de texto, podemos importarlo a R, para editarlo. Buscar archivo “file1.txt”" en el directorio de trabajo. Deben primero descargar desde el drive:

setwd("./<Especificar_directorio/")
getwd()
list.files()

sin preocuparnos por el formato del archivo podemos revisarlo con scan. PUEDEN USAR TAB para AUTOCOMPLETAR

num = scan("../file1.txt")
num

usemos un ejemplo más realista con el archivo “tab.txt”. Descargar del drive

tab = read.table("tab.txt")
head(tab)

existe una forma rápida de visualizar este tipo de datos

?boxplot
boxplot(x = tab, col =c("red","green","blue"))

hemos visto previamente que usando “[ ]” podemos restringir los elementos del vector. Podemos hacer lo mismo con tablas pero en este caso necesitamos agregar una coma para referirnos a [lineas,columnas]

tab[c(1,6,10), c(1,3)]

EJEMPLOS DE TRABAJO EN DATOS DE EXPRESIÓN

La tabla del archivo “expr.txt” del drive contiene el valor de expresión génica de 10000 genes entre dos condiciones, con tres réplicas cada uno

expr = read.table("expr.txt")
head(expr)
dim(expr)
boxplot(expr, col=rainbow(6), main="Boxplot, datos crudos")

¿Cómo se ve esa gráfica? ¿Qué pasa si graficamos la distribución mediante un histograma?

hist(expr[,1], breaks=50, main=1, col="skyblue", xlab="gene expression")

Pero acá sólo obtenemos el histograma para los datos de la primera columna, ¿Cómo hacemos para graficar todas las columnas?. Necesitamos usar un control de estructura, al emplear “n1:n2”, crea una secuencia de n1 a n5 aumentando de a 1

for (i in 1:5) {
    print(i)
  }

Volvamos a los datos, ¿cuantos histogramas necesitamos? Si no nos acordamos las dimensiones podemos usar dim(), nrow() o ncol()

ncol(expr)

Ahora si, para producir 6 histogramas usamos:

for (i in 1:ncol(expr)) {
    hist(expr[,i], breaks=50, main=i, col="skyblue", xlab="gene expression")
  }

Esto nos grafica los 6 pero cada uno por separado, para poder graficarlos todos juntos, necesitamos establecer cómo queremos la estructura de los gráficos. La estructura por defecto es mfrow=c(1,1), pero nosotros queremos 2 líneas y 3 columnas. [ mfrow: number of Multiple Figures (use ROW-wise)]

pdf("multiplot.pdf")
par(mfrow = c(3,2))
for (i in 1:ncol(expr)) {
  hist(expr[,i], breaks=50, main=i, col="skyblue", xlab="gene expression")
}
dev.off()

¿Cómo se ven esos gráficos? ¿Es similar la distribución de los datos? Si asumimos una distribución idéntica, con algunos valores que se desvían, podemos normalizar usando el valor de la mediana

exprNorm = expr
for (i in 1:ncol(exprNorm)) {
    exprNorm[,i] = exprNorm[,i] - median(exprNorm[,i])
  }

¿Funcionó? Hagamos una imagen para responder

par(mfrow=c(1,2))
boxplot(expr, col=rainbow(6), main="Raw data")
boxplot(exprNorm, col=rainbow(6), main="Normalized data")

COMPARACIÓN DE VALORES ENTRE CONDICIONES

Este experimento es muy simple porque sólo tenemos dos condiciones, por lo que queremos determinar cuales son los genes expresados diferencialente al comparar la condicion “ko”y “wt”.¿Cómo lo podemos hacer? Primero, observamos los datos

head(exprNorm)

1)Necesitamos determinar el valor promedio en ko y wt

2)Comparar el valor promedio de cada gen en ko y wt

Para lo primero hacemos una evaluación con sólo algunos datos usando head() y cambiamos de la estructura de data.frame a matrix

mini=head(exprNorm,3)
mini = as.matrix(head(exprNorm, 3))
mini
rownames(mini)

rownames nos regresa los nombres de los renglones

for (i in rownames(mini)) {
  koMean = mean(mini[i,1:3])
  wtMean = mean(mini[i,4:6])
  print("#############")
  print(i)
  print(koMean)
  print(wtMean)
}

for() es fácil de usar pero, el resultado puede ser difícil de asignar a valores,por lo que no es muy eficiente. Una alternativa eficiente es apply(). Esto nos permite aplicar una función sobre muchos elementos en un sólo paso

apply(mini[,1:3], 1, mean)
apply(mini[,4:6], 1, mean)

apply() necesita de 3 elementos

1. Una tabla (como matriz)

2. Especificar las dimensiones (1=rows, 2=columns)

3. La función a aplicar sobre los datos

EJERCICIO: Usar apply() para obtener la expresión promedio por condición (6 condiciones), no por gen con todos los datos en exprNorm para todos los genes, obtener un vector que contenga la diferencia entre el valor promedio de ko y wt.

Crear el histograma correspondiente

¿Es lo que esperábamos?

exprNorm=as.matrix(exprNorm)
meanko=apply(exprNorm[,1:3],1,mean)
meanwt=apply(exprNorm[,4:6],1,mean)
diff=meanwt-meanko
head(diff)
hist(diff)

ahora obtuvimos las diferencias, pero ¿son significativas?

El objetivo principal cuando analizamos expresión es detectar los cambios asociados a la expresión génica. Sin embargo grandes cambios no es indicativo de “más confiable”. Necesitamos considerar la variación de los valores dentro de las muestras de la misma condición. Para esto se han desarrollado distintos tests y métodos estadísticos. Lo que queremos verificar mediante el test estadístico es si los cambios obtenidos ocurrieron al azar (hipótesis nula). Si encontramos que la probabilidad de que el cambio sea al azar es muy bajo, rechazaremos la hipótesis nula y diremos que el cambio es estadísticamente significativo

Usaremos como ejemplo T-test o T-Student. En este caso asumimos una distribución normal, T-test se ha usado para el análisis de microarreglos. Trabajaremos primero con la tabla mini, sobre el primer gen

mini[1,1:3]
mini[1,4:6]
t.test(mini[1,1:3], mini[1,4:6])

podemos ver que nos da un p-value=0.06. ¿Esto es realmente significativo? Analicemos más genes con la herramienta for() y sólo imprimimos pvalue

for (i in rownames(mini)) {
    t = t.test(mini[i,1:3], mini[i,4:6])
    print(i)
    print(t$p.value) #Imprimir sólo el valor de 't'
  }

ahora analicemos todos los genes, filtrando aquellos que no superen el umbral de p-value de 1

exprNorm$p.value = 1

para crear una nueva columna en exprNorm (se llenará de 1)

for (i in rownames(exprNorm)) {
  exprNorm[i,"p.value"] = t.test(exprNorm[i,1:3], exprNorm[i,4:6])$p.value
}
head(exprNorm)
pvals = exprNorm$p.value
names(pvals) = rownames(exprNorm)
head(pvals)

EJERCICIO:¿Cómo obtienen los pvals más chicos? ¿y los más grandes?¿Cuántos pvals son menores a un umbral, por ejemplo 0.05 o 0.001?

pvals_0.05=pvals[pvals<0.05]
length(pvals)
length(pvals_0.05)
pvals_0.001=pvals[pvals<0.001]
length(pvals_0.001)
head(pvals_0.001)

¿Cómo interpretan este resultado?

¿Creen que hay algún problema con el análisis hasta ahora?

INSTALAR PAQUETES REQUERIDOS para la práctica de mañana