Preparación

En primer lugar, cargamos las librerías necesarias.

if (!require(readr)) install.packages('readr', dependencies = T)
## Loading required package: readr
library(readr)

Ejercicio 1: Enfoque FDR

Cargamos la base de datos que necesitaremos para el primer ejercicio.

LEUCE <- read_rds("LEUCE.RDS")
LEUCECL <- read_rds("LEUCECL.RDS")

En primer lugar, fusionamos los datos de LEUCE con los de LEUCECL

x <- LEUCE
y <- as.numeric(LEUCECL)

Al observar las dimensiones de x, vemos que la tabla está girada ya que nos presenta los datos de 38 genes de 3051 pacientes, por lo que la convertimos.

dim(x)
## [1] 3051   38
x <- data.frame(t(x))

Ahora ya tenemos las observaciones de 3051 genes de 38 pacientes, de los cuales 27 padecen leucemia linfoblástica aguda (ALL) y 11 tienen leucemia mieloide aguda (AML) (o y 1 respectivamente en la table(y))

dim(x)
## [1]   38 3051
table(y)
## y
##  0  1 
## 27 11

La realización de una prueba estándar de dos muestras en el gen número 18 gen produce un estadístico de -2.352069 y un valor p asociado de 0.02425779, lo que sugiere una modesta evidencia de una diferencia en los niveles de expresión media entre los dos tipos de leucemia.

x <- as.matrix(x)
x1 <- x[which(y == 0), ]
x2 <- x[which(y == 1), ]
n1 <- nrow(x1)
n2 <- nrow(x2)
t.out <- t.test(x1[, 18], x2[, 18], var.equal = TRUE)
TT <- t.out$statistic
TT
##         t 
## -2.352069
t.out$p.value
## [1] 0.02425779

Sin embargo, este p-valor p se basa en la suposición de que bajo la hipótesis nula de que no hay diferencia entre los dos grupos, la estadística de la prueba sigue una distribución con 27 + 11 - 2 = 36 grados de libertad. En lugar de utilizar esta distribución nula teórica, podemos dividir aleatoriamente a los 38 pacientes en dos grupos de 20 y 28, y calcular una nueva estadística de prueba. Bajo la hipótesis nula de no diferencia entre los grupos, esta nueva estadística de prueba debería tener la misma distribución que la original. Repetir este proceso 500 veces nos permite aproximarnos a la distribución nula del estadístico de la prueba. Calculamos la fracción de veces que nuestro estadístico de prueba observado supera al estadístico de prueba obtenido mediante el remuestreo.

set.seed(1)
B <- 500
Tbs <- rep(NA, B)
for (b in 1:B) {
  dat <- sample(c(x1[, 18], x2[, 18]))
  Tbs[b] <- t.test(dat[1:n1], dat[(n1 + 1):(n1 + n2)],
                   var.equal = TRUE
  )$statistic
}
mean((abs(Tbs) >= abs(TT)))
## [1] 0.016

Esta fracción, 0,016, es nuestro p-valor basado en el remuestreo. Es bastante inferior al p-valor p de 0.02425779 obtenido utilizando la distribución nula teórica.

Al igual que en el manual, podemos trazar un histograma de los estadísticos de prueba basados en el remuestreo.

hist(Tbs, breaks = 30, xlim = c(-3.5, 3.5), main = "",
    xlab = "Distribución nula de la estadística de prueba", col = 3)
lines(seq(-3.5, 3.5, len = 500),
    dt(seq(-3.5, 3.5, len = 500),
      df = (n1 + n2 - 2)
    ) * 1000, col = 2, lwd = 3)
abline(v = TT, col = 4, lwd = 2)
text(TT + 0.5, 350, paste("T = ", round(TT, 4), sep = ""),
    col = 4)

Por último, implementamos el enfoque FDR de remuestreo de plug-in descrito en el manual. Para su realización, consideraremos un subconjunto aleatorio de 100 genes. Para cada gen, primero calculamos la estadística de prueba observada y, a continuación, producimos 500 estadísticas de prueba remuestreadas. Como podemos obsevar, en este caso el p-valor es muy superior.

m <- 100
set.seed(1)
B <- 500
index <- sample(ncol(x1), m)
Ts <- rep(NA, m)
Ts.star <- matrix(NA, ncol = m, nrow = B)
for (j in 1:m) {
  k <- index[j]
  Ts[j] <- t.test(x1[, k], x2[, k],
                  var.equal = TRUE
  )$statistic
  for (b in 1:B) {
    dat <- sample(c(x1[, k], x2[, k]))
    Ts.star[b, j] <- t.test(dat[1:n1],
                            dat[(n1 + 1):(n1 + n2)], var.equal = TRUE
    )$statistic
  }
}

mean ((abs(Ts) >= abs(TT)))
## [1] 0.28

Ejercicio 2: Hipótesis nulas rechazadas, falsos positivos y FDR estimado

A continuación, calcularemos el número de hipótesis nulas rechazadas R, el número estimado de falsos positivos v, y el FDR estimado, para un rango de valores de umbral c. Los valores del umbral los elegimos utilizando los valores absolutos de los estadísticos de prueba de los 100 genes usados anteriormente.

cs <- sort(abs(Ts))
FDRs <- Rs <- Vs <- rep(NA, m)
for (j in 1:m) {R <- sum(abs(Ts) >= cs[j])
  V <- sum(abs(Ts.star) >= cs[j]) / B
  Rs[j] <- R
  Vs[j] <- V
  FDRs[j] <- V / R}

Ejercicio 3: Genes que serán rechazados con el FDR controlado

Ahora, para cualquier FDR dado, podemos encontrar los genes que serán rechazados. Por ejemplo, con el FDR controlado a 0,1, rechazamos 32 de las 100 hipótesis nulas. Por término medio, cabría esperar que tres de estos genes (es decir, el 10% de 32) fueran falsos descubrimientos. Con un FDR de 0,2, podemos rechazar la hipótesis nula de 38 genes, de los que esperamos que unos seis sean falsos descubrimientos. El índice variable es necesario en este caso, ya que restringimos nuestro análisis a sólo 100 genes seleccionados al azar.

max(Rs[FDRs <= .1])
## [1] 32
sort(index[abs(Ts) >= min(cs[FDRs < .1])])
##  [1]   84  248  270  316  330  526  557  733  858  990 1069 1075 1134 1426 1516
## [16] 1640 1696 1742 1766 1948 1974 2087 2300 2347 2374 2430 2483 2667 2668 2903
## [31] 2922 2937
max(Rs[FDRs <= .2])
## [1] 38
sort(index[abs(Ts) >= min(cs[FDRs < .2])])
##  [1]   84  248  270  316  330  471  526  557  733  858  930  990 1069 1075 1134
## [16] 1395 1426 1516 1639 1640 1696 1742 1766 1826 1948 1974 2087 2300 2347 2374
## [31] 2430 2483 2580 2667 2668 2903 2922 2937

Podemos graficar la tasa estimada de falsos descubrimientos frente al número de hipótesis nulas rechazadas, para 100 genes seleccionados aleatoriamente.

plot(Rs, FDRs, xlab = "Número de rechazos", type = "l",
    ylab = "Tasa de falsos descubrimientos", col = 4, lwd = 3)