Ejercicio 1

Supongamos que vamos a estudiar un gen vegetal llamado “factor de transcripción de estrés por calor A-2” (HSFA2) [1]. Por lo tanto, se mide este nivel de transcripción génica (variable de interes) en dos(2) condicioneso tratamientos (un factor): normal y estresado por calor.Para este ejercicio 1, se reviso en literatura la distribución de la expresión en la población de las hojas de Arabidopsis thaliana para simularar los datos: muestra de 1000 plántulas de Arabidopsis en condiciones normales y en condiciones de estrés por calor.
Supongase que esudios han demostrados que la expresión tiene una distribución normal en esta población con un valor medio de expresion en condicionews normales de \(\mu_1 =2\) y para condiciones de stress los estudios indican que es el doble es decir; \(\mu_2=4\). Adicionalmente, la variabilidad esta alrededor de \(0.01\) y \(0.25\) para el tratamiento 1 y 2 respectivamente(\(\sigma_1=0.1;\sigma_2=0.5\)).
** Los datos se crearan en el formato tibble y para esto hay que instalar el paquete “tidyverse”**

library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.2.3
## Warning: package 'ggplot2' was built under R version 4.2.3
## Warning: package 'tibble' was built under R version 4.2.3
## Warning: package 'tidyr' was built under R version 4.2.3
## Warning: package 'lubridate' was built under R version 4.2.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.0     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
set.seed(1234) # semilla para obtnere el mismo resultado en cada run
#datos del trtamiento 1
xp_normal <- tibble( expresion = rnorm(n = 1000, mean = 2, sd = 0.1),
  condicion = "normal"  )
head(xp_normal)
## # A tibble: 6 × 2
##   expresion condicion
##       <dbl> <chr>    
## 1      1.88 normal   
## 2      2.03 normal   
## 3      2.11 normal   
## 4      1.77 normal   
## 5      2.04 normal   
## 6      2.05 normal
#datos del tratamiento 2
xp_stress <- tibble(expresion = rnorm(n = 1000, mean = 4, sd = 0.5),
                         condicion = "estres por calor")
head(xp_stress)
## # A tibble: 6 × 2
##   expresion condicion       
##       <dbl> <chr>           
## 1      3.40 estres por calor
## 2      4.15 estres por calor
## 3      3.23 estres por calor
## 4      4.32 estres por calor
## 5      4.35 estres por calor
## 6      3.05 estres por calor
#Todo en un tibble
xp = bind_rows(xp_normal, xp_stress) 
head(xp)
## # A tibble: 6 × 2
##   expresion condicion
##       <dbl> <chr>    
## 1      1.88 normal   
## 2      2.03 normal   
## 3      2.11 normal   
## 4      1.77 normal   
## 5      2.04 normal   
## 6      2.05 normal

1.1 Genere una muestra para cada tratamiento de tamaño \(n=500\) y \(\mu_1 =3\)
Se puede visualizar la distribución de la expresión para cada tratamiento y observar, que implicar que \(\sigma_1=0.1 <\sigma_2=0.5\) utilizando el paquete “ggplot2”.

grafica1 <- ggplot(xp, aes(x = expresion, fill = condicion)) +
  ggtitle("Distribucion de la expresion HSF2A ") +
  geom_density(color = "black", alpha = 0.5) +
  scale_x_continuous(labels = waiver())
grafica1

1.2 ¿Como se interpreta el efecto sobre la expresión de que la varainza en condiciones normales es menor que en condiciones de stress?
1.3 Realice una grafica si \(\sigma_1=1.5 <\sigma_2=1.5\).
1.4 Realice e interprete los resultados descriptivos para cada grupo de tratamiento (utilice la función summary)

summary(xp_normal)
##    expresion      condicion        
##  Min.   :1.660   Length:1000       
##  1st Qu.:1.933   Class :character  
##  Median :1.996   Mode  :character  
##  Mean   :1.997                     
##  3rd Qu.:2.062                     
##  Max.   :2.320

No obstante, que los datos fueron generados en forma aleatoria con la distribución normal vamos a aplicar uan prueba para este caso.

Pruebas de normalidad

Queremos responder porque se selecciona una u otra pruebas de normalidad: Kolmogorov-Smirnov, Shapiro-Wilk, Anderson-Darling, Ryan-Joiner, Crammer-von Misses, Ji cuadrado. Las pruebas sirven para determinar la normalidad de los datos [2][3].
Las hipótesis que se contrastan en este caso, para cada una las pruebas señaladas anteriormente, son:

  • H0: Los datos siguen una distribución normal.

  • H1: Los datos no siguen una distribución normal.
    Cabe señalar que las pruebas de Anderson-Darling y Kolmogórov-Smirnov se basan en la función de distribución empírica, entre tanto, la prueba de Ryan-Joiner similar a la prueba de Shapiro-Wilk se basa en regresión y correlación y la prueba Cramer-von Mises utiliza los momentos. Las tres pruebas tienden a ser adecuadas para identificar una distribución no normal cuando la distribución es asimétrica y las tres pruebas distinguen menos la distribución normal cuando la distribución subyacente es una distribución t y la no normalidad se debe a la curtosis. Por lo general, entre las pruebas que se basan en la función de distribución empírica, la prueba de Anderson-Darling tiende a ser más efectiva para detectar desviaciones en las colas de la distribución. Generalmente, si la desviación de la normalidad en las colas es el problema principal, se recomienda el uso de la prueba de Anderson-Darling como primera opción. Es decir, las pruebas de Anderson-Darling y Ryan-Joiner tienen una potencia similar para detectar la ausencia de normalidad, mientras la prueba de Kolmogórov-Smirnov tiene una potencia menor para obtener explicaciones con respecto a la normalidad de los datos.
    ### Packages

library(nortest) ###REALIZA 10 PRUEBAS DE NORMALIDAD###
library(moments) ###REALIZA 1 PRUEBA DE NORMALIDAD###

Nivel de Significancia

El nivel de significancia que se trabajará es de 0.05. \(\alpha=0.05\)

Criterio de Decisión

Si $p < $ Se rechaza Ho; Si \(p >= \alpha\) No se rechaza Ho.

###Prueba de Anderson-Darling 
ad.test(xp_normal$expresion)
## 
##  Anderson-Darling normality test
## 
## data:  xp_normal$expresion
## A = 0.57857, p-value = 0.1323
###Prueba de Cramer-von Mises: útil para pequeñas muestras 
cvm.test(xp_normal$expresion)
## 
##  Cramer-von Mises normality test
## 
## data:  xp_normal$expresion
## W = 0.067511, p-value = 0.3012
###Pruena de Lilliefors (Kolmogorov-Smirnov) para una distribución especifica
lillie.test(xp_normal$expresion)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  xp_normal$expresion
## D = 0.021443, p-value = 0.3213
###Prueba de Pearson chi-square: corresponde a una prueba de bondad de ajuste
pearson.test(xp_normal$expresion)
## 
##  Pearson chi-square normality test
## 
## data:  xp_normal$expresion
## P = 23.744, p-value = 0.7415
###Prueba de Shapiro-Wilk: Es más poderosa cuando se compara con otras pruebas de normalidad cuando la muestra es pequeña.
shapiro.test(xp_normal$expresion)
## 
##  Shapiro-Wilk normality test
## 
## data:  xp_normal$expresion
## W = 0.99737, p-value = 0.1053

Como podemos ver todas las pruebas no rechazan \(H_0\) No obstante, la prueba de Shapiro- Wilk tiene el valor p más pequeño.
2.1 Genere una muestra aleatoria de la normal de tamaño n=5000 con un \(\alpha=0.1\) para comprobar la normalidad de los datos de expresión.

Pruebas de Hipotesis

Supongamos que queremos compara o contranstar los valores de expresión entre la condición normal vs bajo stress, es decir:
H0: La expresión promedio del gen HSFA2 es la misma en condiciones normales que bajo estrés por calor en mis plántulas de Arabidopsis”. H1: “La expresión promedio del gen HSFA2 es diferente bajo estrés por calor en comparación con las condiciones normales en mis plántulas de Arabidopsis”.
Teniendo en cuenta que la variable es cuantitativa, datos gaussianos y que queremos comparar dos grupos; se puede realizar una prueba t de Student para probar el hipótesis nula H0 (no hay diferencia entre los medios normales y de estrés por calor).

t.test(x = xp_normal$expresion, 
       y = xp_stress$expresion, 
       alternative = "two.sided" ,      
       var.equal = FALSE,#recordar que la varianzas no son iguales    
       conf.level = 0.95) 
## 
##  Welch Two Sample t-test
## 
## data:  xp_normal$expresion and xp_stress$expresion
## t = -126.96, df = 1081.4, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.040977 -1.978850
## sample estimates:
## mean of x mean of y 
##  1.997340  4.007254

El resultado de la prueba t de Student arroja: \(H_0\) se rechaza lo que significa que hay diferencias significativas en la expresión media del gen HSFA2 entre las condiciones normales y con stress.
2.2 Realice un prueba t con una muestra de \(n=100\) para cada grupo con \(\sigma_1=\sigma_2=0.5\) y valor de media iguales.
2.3 Interprete el intervalo de confianza

Errores de tipo I y tipo II

Los errores de tipo I se producen cuando el H0 Hipótesis (nula) se rechaza erróneamente a favor de una hipótesis alternativa.
Errores de una prueba ###Análisis de potencia: ¿cuántas muestras necesito?

Para calcular la potencia estadística \(1-\beta\): la probabilidad de que la prueba (en este caso \(t\) encuentre una diferencia estadísticamente significativa entre las condiciones normal y bajo stress, en función del tamaño de la verdadera diferencia entre esas dos poblaciones (efecto) es decir, cuando realmente hay diferencias.Para esto, necesitamos:

Error de tipo I: controlado por el valor \(\alpha\). A menudo se establece en 0.01 (1%) o 0.001 (0.1%) en experimentos de RNA-seq.
Error de tipo II: controlado por el valor \(\beta\)).Debe establecerse en 70 u 80% para detectar el 70 u 80% de los genes expresados diferencialmente. El número de réplicas biológicas puede ser difícil de alcanzar en la práctica para los experimentos de RNA-seq.
Tamaño del efecto: Por ejemplo, si desea investigar genes que difieren entre tratamientos con una diferencia de su media de 2, entonces el tamaño del efecto es igual a 2 (ver presentación de la clase).
Si tenemos el valor de \(n\) entonces, se evalua la potencia (aposteriori) pero sino, podemos calcular el tamaño de la muestra con base en una potencia deseada (apriori). Vamos a hacer uso del paquete de R (ejemplo de uso aquí).pwr.

library("pwr")
## Warning: package 'pwr' was built under R version 4.2.3
pwr.t.test(d = 2,
           power = .8,
           sig.level = .05,
           type = "two.sample",
           alternative = "two.sided")
## 
##      Two-sample t test power calculation 
## 
##               n = 5.089995
##               d = 2
##       sig.level = 0.05
##           power = 0.8
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

Lo que significa que el tamaño de la muestra indicado es 6 aproximadamente para condiciones normales y 5 muestras en condiciones de stress.

Análisis de potencia para RNA-seq

Los experimentos de RNA-seq a menudo sufren de un bajo poder estadístico. Un bajo poder estadístico refleja un error de tipo II [4]. Hay que recordar que baja potencia implica que faltan “efectos verdaderos” y hay muchos falsos negativos: genes que deberías haber sido llamados diferenciales pero no lo son. También afecta el valor predictivo positivo de sus hallazgos, que es la probabilidad de que un gen diferencial (p < 0.01) sea un hallazgo verdadero real. Por lo tanto, el poder está vinculado a la capacidad de reclamar “hallazgos verdaderos”.

pwr.t.test(n=5:50,#probar con tamaños de muestra de 5 hasta 50
           d = 1,
                      sig.level = .05,
           type = "two.sample",
           alternative = "two.sided")
## 
##      Two-sample t test power calculation 
## 
##               n = 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50
##               d = 1
##       sig.level = 0.05
##           power = 0.2862955, 0.3473537, 0.4058070, 0.4612388, 0.5133625, 0.5620066, 0.6070978, 0.6486426, 0.6867106, 0.7214198, 0.7529230, 0.7813978, 0.8070367, 0.8300406, 0.8506124, 0.8689530, 0.8852578, 0.8997137, 0.9124984, 0.9237780, 0.9337077, 0.9424303, 0.9500773, 0.9567684, 0.9626126, 0.9677083, 0.9721439, 0.9759988, 0.9793441, 0.9822429, 0.9847512, 0.9869189, 0.9887897, 0.9904023, 0.9917908, 0.9929848, 0.9940105, 0.9948906, 0.9956451, 0.9962912, 0.9968439, 0.9973164, 0.9977198, 0.9980641, 0.9983575, 0.9986074
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

3.1 ¿cómo interpreta estos resultados?
3.2 ¿cuál es el tamaño de muestra a seleccionar?
3.3 De un ejemplo desde su area de interes para evaluar el tamaño de la muestra. 3.4 ¿Cuál es el tamaño de la muestra si queremos verificar si la media de la expresión en condiciones normales es menor que bajo stress?

Referencias

[1] Nishizawa A, Yabuta Y, Yoshida E, Maruta T, Yoshimura K, Shigeoka S. Arabidopsis heat shock transcription factor A2 as a key regulator in response to several types of environmental stress. Plant J. 2006 Nov;48(4):535-47. doi: 10.1111/j.1365-313X.2006.02889.x. Epub 2006 Oct 19. PMID: 17059409

[2] IBM. (2020). SPSS Statistics. https://www.ibm.com/products/spss-statistics

[3] Jensen, W., & Alexander, M. (2016). Statistics for Engineering and the Sciences. Journal of Quality Technology, 48(3), 297–299. https://www.tandfonline.com/doi/abs/10.1080/00224065.2016.11918168.
[4] Button, K., Ioannidis, J., Mokrysz, C. et al. Fallo de alimentación: por qué el pequeño tamaño de la muestra socava la fiabilidad de la neurociencia. Nat Rev Neurosci 14, 365–376 (2013). https://doi.org/10.1038/nrn3475.