El problema de la profundidad variable por célula

En S27 se han conservado únicamente las células limpias. Al intentar comparar estas células entre sí, nos encontramos con el primer obstáculo: el número total de UMI varía según la célula. Algunas células tienen 5.000 recuentos, mientras que otras solo 500. La mayor parte de esta diferencia no se debe a factores biológicos, sino a ruido técnico —qué tan bien se capturó cada célula dentro de la gota y con qué eficiencia ocurrió la transcripción inversa—.

Si un gen aparece 50 veces en una célula con profundidad 5.000 y 5 veces en una célula con profundidad 500, la proporción relativa es la misma. Es necesario eliminar esta diferencia técnica de profundidad para poder comparar las células de manera justa. Posteriormente, de entre los 30.000 genes, debemos seleccionar únicamente aquellos que realmente ayudan a distinguir las células. Estos dos pasos de preprocesamiento —la normalización y la selección de HVG— son el tema de esta sección.

Normalización estándar: log1p

El método más utilizado consta de tres pasos:

1. Corrección de profundidad: dividir los recuentos de cada célula por su suma total y escalar a un tamaño objetivo (por ejemplo, 1e4 = CP10K). "Unificar a 10.000 recuentos por célula".
2. Transformación logarítmica: aplicar $\log(x+1)$. +1(pseudocount) evita la divergencia al calcular el logaritmo de cero.
3. (Opcional) Escalado por gen (z-score).

La razón para aplicar el logaritmo es fundamental. Los recuentos de expresión están dominados por unos pocos genes altamente expresados. El logaritmo comprime esta cola larga, reduciendo la diferencia de escala entre valores grandes y pequeños (estabilización de la varianza) y transformando diferencias multiplicativas en aditivas, lo que permite que los análisis lineales posteriores (como PCA) funcionen correctamente.

Este método es rápido y robusto, por lo que se utiliza como valor predeterminado. Sin embargo, no es perfecto. En genes de baja expresión, log1p no logra eliminar completamente el efecto de la profundidad.

Un enfoque más sofisticado: residual de Pearson (sctransform)

sctransform adopta un planteamiento diferente. Modela cada gen mediante una regresión binomial negativa (usando la profundidad celular como variable explicativa) y expresa cuánto se desvía el valor observado del valor esperado del modelo mediante el residual de Pearson.

$\mu_{ci}$ es el conteo esperado que refleja la profundidad, y el denominador es la desviación estándar de la distribución binomial negativa (¡con una estructura similar a $\mu + \alpha\mu^2$ de S19!). De esta manera, el efecto de la profundidad se elimina explícitamente en el modelo, y la varianza de los genes altamente y poco expresados se estabiliza uniformemente. El cálculo es más costoso, pero es robusto cuando se deben distinguir tipos celulares sutiles. Scanpy lo ofrece mediante sc.experimental.pp.normalize_pearson_residuals.

HVG — De 30.000 a 2.000

Aunque tras la normalización quedan 30.000 genes, la mayoría son genes de mantenimiento (housekeeping) que se expresan de manera similar en todas las células y no permiten distinguir entre ellas. Lo que separa los tipos celulares es un pequeño número de genes con alta variabilidad.

La trampa principal: no basta con seleccionar "los genes con mayor varianza". Los genes con alta expresión tienen naturalmente una varianza mayor (relación media-varianza). Por lo tanto, se deben seleccionar aquellos cuya varianza sea mayor a la esperada en relación con su media. highly_variable_genes de Scanpy divide los genes en intervalos (bins) según su expresión media y luego selecciona aquellos con una varianza estandarizada (dispersion) alta dentro de cada intervalo. Por lo general, se conservan los 2.000 superiores.

Esto constituye la propia selección de características (feature selection). Al conservar solo un pequeño subconjunto discriminativo entre decenas de miles de características, se mitiga la maldición de la dimensionalidad y se mejora la relación señal-ruido en el PCA y el clustering posteriores.

Práctica con Scanpy (PBMC 3k)

Continuamos desde el objeto adata que ha pasado el QC de S27.

import scanpy as sc
# adata = Estado completado del QC de S27

# Copia de seguridad de los conteos crudos (para poder revertir)
adata.layers["counts"] = adata.X.copy()

# Normalización estándar
sc.pp.normalize_total(adata, target_sum=1e4)   # Corrección de profundidad
sc.pp.log1p(adata)                              # Transformación logarítmica
adata.raw = adata                               # Conservación de los valores normalizados por log

# Selección de HVG (top 2000)
sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor="seurat")
sc.pl.highly_variable_genes(adata)              # Visualización del resultado de la selección

# Conservar solo los HVG
adata = adata[:, adata.var.highly_variable].copy()
print(adata)                                    # Verificar si n_vars se redujo a 2000

En el gráfico de dispersión generado por sc.pl.highly_variable_genes, los puntos marcados en negro corresponden a los HVG seleccionados. Observemos visualmente cómo se seleccionan los genes que sobresalen por encima de la línea de tendencia media-varianza. Estos genes destacados (por ejemplo, las familias de marcadores en el caso de PBMC) serán los protagonistas en la posterior identificación de tipos celulares.

Es importante tener el hábito de conservar adata.raw y layers["counts"], ya que la visualización de genes marcadores (S30) o el análisis de expresión diferencial requieren alternar entre valores previos y posteriores a la normalización.

Mapeo CS

• Selección de características (feature selection): La selección de HVG consiste en conservar un pequeño número discriminante entre decenas de miles de características, siendo esta una etapa fundamental del preprocesamiento en machine learning.
• Estabilización de la varianza: log1p y los residuos de Pearson son transformaciones que igualan la varianza dependiente de la media. Cumplen el mismo propósito que las transformaciones de estabilización de la varianza en estadística (como la de Anscombe, etc.).
• Normalización = invariancia de escala: La corrección de la profundidad de secuenciación consiste en ajustar las observaciones a una escala común, lo cual está directamente relacionado con la estandarización de datos.
• Mitigación de la maldición de la dimensionalidad: Al reducir las características, se aumenta la fiabilidad de los cálculos de distancia posteriores (PCA·kNN).

Errores frecuentes

• Realizar la selección de HVG antes de la normalización: El efecto de la profundidad contamina la varianza, lo que provoca la selección de genes incorrectos. Debemos respetar el orden: normalización → HVG.
• No realizar un respaldo de los conteos crudos (raw counts): Los datos crudos son necesarios posteriormente para la expresión diferencial y las pruebas de marcadores. Es esencial seguir el hábito de layers["counts"].
• Realizar la expresión diferencial utilizando valores escalados (z-score): Los valores escalados son para visualización y PCA. La DGE debe realizarse con los valores log-normalizados (adata.raw).
• Aumentar el número de HVG de forma arbitraria: Lo habitual es alrededor de 2,000. Si hay demasiados, se introduce ruido; si hay muy pocos, se pierde la señal de los tipos celulares raros.

Para profundizar más

El texto principal es una descripción reconstruida directamente por BPD. Para profundizar, utilicemos los siguientes materiales de referencia:

• Libro web gratuito de referencia: Single-cell best practices (sc-best-practices.org) — Los capítulos sobre normalización y selección de características son la referencia definitiva para la comparación de métodos.
• Artículo original (sctransform): Hafemeister & Satija (2019), Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression, Genome Biology 20:296.
• Artículo original (scran pooling): Lun, Bach, Marioni (2016), Pooling across cells to normalize single-cell RNA sequencing data with many zero counts, Genome Biology 17:75.
• Formación de EMBL-EBI — Single-Cell Transcriptomics: Normalisation & Feature Selection (con subtítulos). Proporciona una intuición práctica sobre los parámetros.

Ahora se dispone de una matriz de 2.000 genes × células limpias. Sin embargo, las 2.000 dimensiones siguen siendo inobservables para el ser humano. En el siguiente episodio S29, reduciremos esta alta dimensionalidad mediante PCA, t-SNE y UMAP para desplegar la estructura que forman las células en un mapa bidimensional.