SAM

Guillermo Ayala Gallego

2025-05-21

Método SAM

Es la abreviatura de Significance Analysis of Microarrays.
El conjunto de muestras es \(\{1,\ldots,n\}\).
Empezamos con el caso en que pretendemos comparar dos grupos.
Para un gen dado las expresiones del grupo 1 (o del grupo 2) corresponden con los valores donde \(y_j = 1\) (respectivamente \(y_j = 2\)).
Tendremos \(J_1\) y \(J_2\) donde \(J_g\) son las muestras del grupo g con g= 1,2.
Tenemos \[ \bar{x}_{i g} = \frac{\sum_{j \in J_g} x_{ij}}{n_g} \] siendo \(n_g\) el número de muestras en \(J_g\).

Para el i-ésimo gen consideramos la diferencia relativa dada por \[ d_i = \frac{\bar{x}_{i 1} - \bar{x}_{i 2}}{s_i + s_0} \] donde \[ s_i = \sqrt{a \sum_{g=1}^2 \sum_{j \in J_g}(x_{ij} - \bar{x}_{i g})^2} \] y \[ a = \frac{1/n_1 + 1/n_2}{n_1 + n_2} \]
El estadístico \(d_i\) es el t-estadístico al que le modificamos el error estándar.
Si el valor \(s_i\) es muy pequeño entonces el valor \(d_i\) es muy grande.
Los distintos valores \(d_i\) no son comparables para los distintos genes, no son intercambiables.
Por ello se le suma un valor \(s_0\) que estabiliza la varianza (y nos permite comparar los t-valores entre si).
Este valor hay que estimarlo o calcularlo a partir de los propios datos.

El procedimiento

Calculamos los estadísticos \(d_i\) y los ordenamos: \[d_{(1)} \leq \ldots \leq d_{(N)}.\]
Realizamos \(B\) permutaciones aleatorias del vector \(y\) que nos indica el grupo de la columna.
Para cada permutación, \(b\), calculamos los nuevos valores \(d_i\) y los denotamos por \(d_i(b)\).
Ordenamos \[d_{(1)}(b) \leq \ldots \leq d_{(N)}(b).\]
Calculamos, para las \(B\) permutaciones, el valor medio de estos estadísticos ordenados. \[ \bar{d}_{(i)} = \sum_{b=1}^B \frac{d_{(i)}(b)}{B}. \]
Para un valor fijo positivo \(\Delta\) determinamos:
- Empezando en el origen y moviéndonos a la derecha determinamos el primer \(i = i_1\) tal que \[d_{(i)} -\bar{d}_{(i)} > \Delta.\] Todos los genes que verifican esta desigualdad son llamados significativamente positivos.
- Empezando en el origen y moviéndonos a la izquierda determinamos el
  primer \(i = i_2\) tal que \[\bar{d}_{(i)} - d_{(i)} > \Delta.\] Todos los genes que verifican esta desigualdad los llamamos significativamente negativos.
- Para cada \(\Delta\) definimos un punto de corte superior como \[ u(\Delta) = \min \{d_i: i \text{ es significativamente positivo} \}, \] y \[ l(\Delta) = \max \{d_i: i \text{ es significativamente negativo} \}. \]
Para distintos valores de \(\Delta\), calculamos:

El número total de genes significativos.
En la \(b\)-ésima permutación aleatoria hemos obtenido \(d_i(b)\) con \(i =1,\ldots, N\). Calculamos el número de genes falsamente llamados en la aleatorización \(b\)-ésima como \[ c_b = | \{i: d_i(b) > u(\Delta) \text{ o } d_i(b) < l(\Delta) \}\, \] Tendremos los valores \(c_b\) con \(b=1,\ldots, B\). Notemos que los valores \(d_i(b)\) han sido calculados bajo la hipótesis de no asociación, es cierta la hipótesis nula para cada gen. En consecuencia debieran de estar en el intervalo \([l(\Delta),u(\Delta)]\). Calculamos la mediana y el percentil de orden 0.90 de los valores \(c_b\) con \(b=1,\ldots, B\) y los denotamos como \(q_{0.5}(c)\) y \(q_{0.90}(c)\).

Estimamos \(\pi_0\), la proporción de genes sin expresión diferencial.

Calculamos los percentiles de orden 0.25 y 0.75, \(q_{0.25}\) y \(q_{0.75}\), de \(\{d_i(b): i = 1, \ldots,N; b = 1, \ldots,B\}\).
Calculamos \[ \hat{\pi}_0 = \frac{|\{d_i: d_i \in (q_{0.25},q_{0.75})\}|}{N/2} \] donde \(\{d_1, \ldots, d_N\}\) son los \(d\) valores originales.
Truncamos el valor de \(\hat{\pi}_0\) en 1, es decir, consideramos \[ \hat{\pi}_0 = \min\{\hat{\pi}_0,1\}. \]

Multiplicamos la mediana, \(q_{0.5}(c)\), y el percentil 0.9, \(q_{0.90}(c)\) por \(\hat{\pi}_0\).
Elegimos un valor de \(\Delta\) y determinamos el conjunto de genes significativos.
La tasa de falsos positivos FDR es estimada como el cociente entre la mediana, \(q_{0.5}(c)\), (o el percentil 0.90, \(q_{0.90}(c)\)) del número de genes falsamente llamados y el número de genes declarados significativos.
El q-valor: + Es pFDR cuando la lista de genes significativos está compuesta por este gen y todos los genes que son más significativos que este. + Se calcula buscando cuál es el valor más pequeño de \(\Delta\) para el cual este gen es declarado como significativo. + Si este valor es \(\Delta\) entonces la pFDR asociada a este valor es el q-valor asociado al gen.

Cálculo de \(s_0\)

Consideramos \(s = (s_1,\ldots,s_N)\) las desviaciones estándar muestrales.
Sea \(q_{\alpha}(s)\) el percentil de orden \(\alpha\) de los valores \(s_i\). Definimos \[ d_i^{(\alpha)} = \frac{\bar{x}_{i 1} - \bar{x}_{i 2}}{s_i + q_{\alpha}(s)}. \]
Consideramos los percentiles \(q_{j/100}(s)\) con \(j=1,\ldots,100\).
Para cada \(\alpha \in \{\) 0, 0.1, \(\ldots, 1\}\):
1. Se calcula \[ v_j = \frac{1}{0.64} MAD(\{d_i^{(\alpha)}: s_i \in [q_{j/100}(s), q_{(j+1)/100}(s)]\}) \] para \(j=1,\ldots,100\) donde MAD es la mediana de las desviaciones absolutas respecto de la mediana
2. Calculamos el valor \(CV(\alpha)\), el coeficiente de variación de los valores \(v_j\).
3. Se elige como valor de \(\alpha\): \(\hat{\alpha}\) que minimiza \(CV(\alpha)\).
Tomamos para \(s_0\): \(\hat{s}_0 = q_{\hat{\alpha}}(s)\).

samr

Cargamos el paquete.

pacman::p_load(samr)

Cargamos los datos.

data(golub,package="multtest")
fac0 = golub.cl +1 ## samr espero valores 1 y 2 para los grupos

Aplicamos el método: tenemos dos grupos independientes.

samfit = SAM(golub,fac0,resp.type="Two class unpaired",nperms=1000)

¿Qué valor de \(\Delta\) se ha utilizado?

samfit$del

    delta 
0.4607407

sigtabla = samfit$siggenes.table

Veamos gráficamente el resultado.

plot(samfit)

Podemos ver también una exploración de posibles valores para \(\Delta\) con

head(samfit$delta.table)

         delta # med false pos 90th perc false pos # called median FDR
[1,] 0.4607407       356.33661            503.3491     1689 0.21097490
[2,] 0.5289115       255.99049            390.6485     1506 0.16998041
[3,] 0.6017838       176.41495            293.0538     1339 0.13175127
[4,] 0.6793575       118.68895            211.5361     1188 0.09990653
[5,] 0.7616326        77.68732            149.0625     1071 0.07253718
[6,] 0.8486091        46.93609            102.5041      932 0.05036061
     90th perc FDR      cutlo     cuthi
[1,]     0.2980160 -0.8650472 0.9722594
[2,]     0.2593948 -0.9883173 1.1236393
[3,]     0.2188602 -1.1215032 1.2757362
[4,]     0.1780607 -1.2637244 1.4182065
[5,]     0.1391807 -1.4105379 1.5526948
[6,]     0.1099829 -1.5606679 1.7121587

Modificando FDR

Tenemos muchos genes significativos. Parece que podemos ser más exigentes con la tasa FDR.

samfit = SAM(golub,fac0,resp.type="Two class unpaired",nperms=1000,fdr.output=.001)
sigtabla = samfit$siggenes.table

La tabla resumen sería

sigtabla = samfit$siggenes.table

Otro tipo de covariables

Comparación de más de dos grupos

Tenemos \(K\) grupos a comparar (\(K >2\)) por lo que \(y_j \in \{1,\ldots,K\}\).
\(J_k\) son los índices de las observaciones en grupo \(k\).
\(N_k\) es el cardinal de \(J_k\) (\(\sum_{k=1}^K n_k = n\)).
Definimos \(d_i\) como \[ d_i = \frac{r_i}{s_i + s_0} \] con \[ r_i = \sqrt{\frac{\sum_{k=1}^K n_k}{\prod_{k=1}^K n_k} \sum_{k=1}^K n_k (\bar{x}_{i k} - \bar{x}_{i \cdot})^2} \] \[ s_i = \sqrt{\frac{1}{\sum_{k=1}^K (n_k-1)} \bigg ( \sum_{k=1}^K \frac{1}{n_k} \bigg ) \sum_{k=1}^K \sum_{j \in J_k} (x_{ij} - \bar{x}_{i k})^2}. \]

Covariable continua

Supongamos que los valores de \(y\) son valores numéricos.
Definimos \(d_i\) como \[ d_i = \frac{r_i}{s_i + s_0} \]
Tenemos: \[ r_i = \frac{\sum_{j=1}^n y_j (x_{ij} - \bar{x}_{i \cdot})}{\sum_{j=1}^n (y_j - \bar{y})^2} \] y \[ s_i = \frac{\hat{\sigma}_i}{ \sqrt{\sum_{j=1}^n (y_j - \bar{y})^2}} \] donde \[ \hat{\sigma}_i = \sqrt{\frac{\sum_{j=1}^n (x_{ij} - \hat{x}_{ij})^2}{n-2}}, \] y \[ \hat{x}_{ij} = \hat{\beta}_{i0} + r_i y_j, \] \[ \hat{\beta}_{i0} = \bar{x}_{i \cdot} - r_i \bar{y}. \]