¿Por qué la transformación logarítmica de los datos de RNA-seq reduce la cantidad de varianza explicada en el PCA?

Question

Estoy ejecutando un PCA en un conjunto de datos con 2k filas y 36k columnas. Me he dado cuenta de que cuando transformo los datos en logaritmos tengo que pedir más componentes principales durante el ACP para conseguir la misma cantidad de varianza explicada (ver imágenes adjuntas).

¿Cuál es la razón subyacente? Tengo la intuición ingenua de que la transformación logarítmica "comprime" los datos, por lo que la cantidad de varianza explicada para los diferentes componentes principales se está "homogeneizando". ¿Existe una explicación rigurosa para esto o he cometido un error?

Answer 1

Basándome en el tamaño de su conjunto de datos, sospecho que está trabajando con los datos de RNA-seq de una sola célula.

Si es así, puedo confirmar su observación: con los datos de scRNA-seq, las varianzas explicadas por el PCA después de la log-transformación son típicamente mucho más bajas que antes. Aquí hay una réplica de su hallazgo con el Tasic et al. 2016 conjunto de datos que tengo a mano:

Aquí he utilizado $\log(x+1)$ debido a los ceros exactos. Obsérvese que los datos transformados logarítmicamente producen varianzas explicadas aproximadamente similares a las de los datos estandarizados (cuando cada variable está centrada y escalada para tener una varianza unitaria).

La razón es que las diferentes variables (genes) tienen varianzas MUY diferentes. Los datos de RNA-seq son, en última instancia, recuentos de moléculas de RNA, y la varianza crece monotónicamente con la media (piense en la distribución de Poisson). Así que los genes que se expresan mucho tendrán una varianza alta, mientras que los genes que apenas se expresan o se detectan en absoluto, tendrán una varianza casi nula:

Sin ninguna transformación, hay un gen que solo explica más del 40% de la varianza (es decir, su varianza es superior al 40% de la varianza total). En este conjunto de datos, resulta ser este gen: https://en.wikipedia.org/wiki/Neuropeptide_Y que está muy expresado (valores RPKM superiores a 100000) en algunas células y tiene una expresión nula en otras. Cuando se realiza el PCA sobre los datos brutos, el PC1 coincidirá básicamente con este único gen.

Esto es similar a lo que ocurre en la respuesta aceptada a ¿PCA sobre correlación o covarianza? :

Obsérvese que el ACP sobre la covarianza está dominado por run800m y javelin : PC1 es casi igual a run800m (y explica el 82% de la varianza) y PC2 es casi igual a javelin (juntos explican el 97%). El ACP sobre la correlación es mucho más informativo y revela cierta estructura en los datos y las relaciones entre las variables (pero hay que tener en cuenta que las varianzas explicadas bajan al 64% y al 71%).

Actualización

En los comentarios, @An-old-man-in-the-sea sacó a relucir el tema de la transformaciones estabilizadoras de la varianza . Los recuentos de RNA-seq suelen modelarse con una distribución binomial negativa que tiene la siguiente relación media-varianza: $$V(\mu) = \mu + \frac{1}{r}\mu^2.$$

Si despreciamos el primer término (lo que tiene cierto sentido bajo la suposición de que los genes altamente expresados llevan la mayor información para el PCA), entonces la relación media-varianza restante se convierte en cuadrática, lo que coincide con la distribución log-normal y tiene el logaritmo como transformación estabilizadora de la varianza: $$\int\frac{1}{\sqrt{\mu^2}}d\mu=\log(\mu).$$

Alternativamente, las pequeñas potencias $\mu^\alpha$ con, por ejemplo $\alpha\approx 0.1$ más o menos también puede tener sentido y se estabilizaría la varianza para $V(\mu)=\mu^{2-2\alpha}$ Así que algo entre la relación lineal y la cuadrática de la media-varianza.

Otra opción es utilizar $$\int\frac{1}{\sqrt{\mu+\frac{1}{r}\mu^2}}d\mu=\operatorname{arsinh}\Big(\frac{x}{r}\Big)=\log\Big(\sqrt{\frac{\mu}{r}}+\sqrt{\frac{\mu^2}{r^2}+1}\Big),$$ posiblemente con la corrección de Anscombe como $$\operatorname{arsinh}\Big(\frac{x+3/8}{r-3/4}\Big).$$ Está claro que para los grandes $\mu$ todas estas fórmulas se reducen a $\log(x)$ .

Ver Harrison, 2015, la transformación estabilizadora de la varianza de 1948 de Anscombe para la distribución binomial negativa se adapta bien a los datos de expresión de RNA-Seq y Anscombe, 1948, The Transformation of Poisson, Binomial and Negative-Binomial Data .

Answer 2

El PCA sólo utiliza el álgebra lineal para encontrar los "mejores" componentes para explicar la varianza. Véase esta respuesta para una explicación más larga sobre el PCA. Por lo tanto, si su conjunto de datos ya tiene relaciones lineales entre sus variables, obtendrá los mejores modelos lineales sin ninguna transformación no lineal en los datos brutos.

Puede volver a formular su pregunta utilizando exp o tanh o cualquier función no lineal en lugar de log y se podría esperar un efecto similar. Básicamente estás empeorando las relaciones lineales entre las variables cuando aplicas el logaritmo (o cualquier función) a una variable. Normalmente es bueno tener alguna base teórica para logaritmizar una variable.