Mapa de calor de RNA-seq: normalización, agrupamiento y visualización integrados

Al finalizar este tema

Podrás crear una herramienta que combine la normalización, el agrupamiento jerárquico y los mapas de calor de Seaborn aprendidos en los libros de texto, para que reciba automáticamente una matriz de recuentos de RNA-seq, la normalice, agrupe las muestras y los genes con patrones de expresión similares y genere un mapa de calor adecuado para incluir en publicaciones. Comprenderás que esa imagen atractiva es, en realidad, la superposición de tres capas.

Este texto es un ejemplo educativo general. Se ha utilizado RNA-seq como material porque es el método estándar en los estudios de expresión actuales.

"¿Por qué cambia la imagen si solo cambio el orden de las muestras?" — Las trampas de la normalización y el ordenamiento

Una vez finalizado el experimento de RNA-seq, se obtiene una matriz de recuentos como esta.

gene       tumor1  tumor2  tumor3  normal1  normal2  normal3
BRCA1      1520    1830    1650    340      420      380
TP53       850     920     780     650      710      680
MYC        4200    5100    4700    1200     1350     1280
GAPDH      18000   19500   17800   17200    17900    18300
ACTB       15000   16800   15500   14800    15600    15200
...

¿Qué pasaría si simplemente coloreáramos esta tabla como un mapa de calor?

• BRCA1 tiene valores que oscilan entre cientos y miles. Las celdas azul oscuro apenas son visibles.
• GAPDH · ACTB tienen valores del orden de decenas de miles. El mapa de calor completo se convierte en una franja roja brillante donde solo destacan estos dos genes.
• Al colocar el orden de las columnas de muestra según el tiempo experimental, los tumores y los normales están mezclados al azar, por lo que no se aprecia ningún patrón.

Al encontrarse con estos problemas, comprenderán que necesitan tres tipos de procesamiento diferentes:

1. Normalización: para hacer comparables los rangos de expresión entre genes y muestras.
2. Ordenamiento (agrupamiento): para colocar elementos similares juntos y que los patrones sean visibles.
3. Mapeo de colores: para convertir los valores en colores visualmente distinguibles.

Estas tres capas deben estar bien integradas para obtener una figura apta para publicación científica. A continuación, construiremos y ensamblaremos estas tres capas una por una.

Primero veamos el producto terminado (ejecutar la caja negra primero)

La herramienta que vamos a crear se usa así:

plot_rnaseq_heatmap(
    counts_df,
    top_n_variable=50,        # Solo los 50 genes con mayor variación
    cluster_samples=True,
    cluster_genes=True,
    save_path="heatmap.png",
)

Los resultados muestran lo siguiente.

=== Lectura del mapa de calor ===
1. Las tres muestras tumorales se agrupan automáticamente a la izquierda
2. Las tres muestras normales se agrupan automáticamente a la derecha
3. Los genes especialmente intensos en tumor se agrupan arriba
4. Los genes intensos en normal se agrupan abajo
5. Color = z-score (valor estandarizado según la variación entre muestras)
6. Dendrograma superior = jerarquía de muestras
7. Dendrograma izquierdo = jerarquía de genes

Archivo guardado: heatmap.png (300 dpi, apto para incluir en un artículo)

Con los mismos datos, siempre se obtiene la misma imagen. La reproducibilidad se logra mediante un único nombre de archivo.

¿De qué componentes está compuesta esta herramienta (diagrama de componentes)?

Pipeline de mapa de calor de RNA-seq
   ┌──────────────────────────────────────────────────┐
   │ [Entrada] Carga de matriz counts ─ pieza: pandas    │ ← Se proporciona completa
   │              │                                      │
   │              ▼                                      │
   │ [Paso 1] Normalización (CPM y z-score)                │
   │        pieza: normalization                           │ ← La construyes tú ★
   │              │                                      │
   │              ▼                                      │
   │ [Paso 2] Selección de genes con alta variación        │
   │        pieza: cálculo de varianza con pandas          │ ← Se proporciona completa
   │              │                                      │
   │              ▼                                      │
   │ [Paso 3] Clustering jerárquico de muestras y genes    │
   │        pieza: hierarchical clustering                 │ ← La construyes tú ★
   │              │                                      │
   │              ▼                                      │
   │ [Salida] seaborn clustermap (mapa + dendrograma)      │ ← Lo ensamblas tú ★
   └──────────────────────────────────────────────────┘

📌 Si es la primera vez que ves estos conceptos (enlaces de entrada arriba)

• Normalización: principio de unificación de escalas
• Agrupamiento jerárquico: reglas de crecimiento del dendrograma
• Mapa de calor de seaborn: asignación de valores de color

Los nuevos conceptos que crearás directamente son normalización, agrupamiento jerárquico y ensamblaje de mapas de calor (3 elementos). pandas es una herramienta, por lo que se proporciona como una solución completa. Solo 3 elementos, dentro del límite cognitivo.

Paso 1: Preparación de datos (se proporciona completo)

Creemos datos que imiten una matriz de conteos de RNA-seq real. En la práctica, estos datos provienen de herramientas como salmon, featureCounts y HTSeq, normalmente como archivos CSV.

import numpy as np
import pandas as pd

rng = np.random.default_rng(42)

def make_rnaseq_counts():
    genes = ["BRCA1", "TP53", "MYC", "KRAS", "EGFR", "PIK3CA", "PTEN", "APC",
             "GAPDH", "ACTB", "B2M", "HPRT1", "TBP"]  # Los últimos cinco son housekeeping
    tumor_samples = [f"tumor{i+1}" for i in range(3)]
    normal_samples = [f"normal{i+1}" for i in range(3)]
    samples = tumor_samples + normal_samples

    counts = np.zeros((len(genes), len(samples)), dtype=int)
    for gi, gene in enumerate(genes):
        # Housekeeping: similar en ambas condiciones
        base_tumor  = 1000 if gene in {"BRCA1", "MYC", "KRAS"} else 300
        base_normal = 300 if gene in {"BRCA1", "MYC", "KRAS"} else 300
        if gene in {"GAPDH", "ACTB", "B2M", "HPRT1", "TBP"}:
            base_tumor = base_normal = 15000  # Housekeeping
        for si, sample in enumerate(samples):
            base = base_tumor if sample.startswith("tumor") else base_normal
            counts[gi, si] = max(0, int(base * rng.lognormal(0, 0.15)))
    return pd.DataFrame(counts, index=genes, columns=samples)

counts = make_rnaseq_counts()
print(counts)

# Verificación
assert counts.shape == (13, 6)
assert (counts >= 0).all().all()
# Por diseño, BRCA1 debe ser mayor en tumor que en normal
assert counts.loc["BRCA1", "tumor1"] > counts.loc["BRCA1", "normal1"]
# Los genes housekeeping son similares en ambas condiciones
housekeeping_tumor = counts.loc["GAPDH", ["tumor1","tumor2","tumor3"]].mean()
housekeeping_normal = counts.loc["GAPDH", ["normal1","normal2","normal3"]].mean()
assert abs(housekeeping_tumor / housekeeping_normal - 1) < 0.3

Paso 2 de creación — Normalización CPM ★ (normalización: comparación entre muestras)

✍️ Sección para completar directamente. El componente es la normalización. Dado que el número total de lecturas varía entre cada muestra, no es posible comparar los recuentos absolutos. Se convierten a proporciones relativas.

La primera razón para la normalización es la diferencia en el número total de lecturas entre las muestras. Si una muestra ha sido secuenciada bien y produce 30 millones de lecturas, mientras que otra solo produce 20 millones, no se pueden comparar directamente los recuentos del mismo gen.

CPM (Recuentos por millón): Normalización del número total de lecturas a un millón.

CPM(gene, sample) = counts(gene, sample) / total_counts(sample) × 1,000,000

def to_cpm(counts):
    """Divide por el total de counts de cada muestra y escala a un millón."""
    library_sizes = counts.sum(axis=0)     # Counts totales de cada muestra
    cpm = counts.div(library_sizes, axis=1) * 1_000_000
    return cpm

cpm = to_cpm(counts)
print(cpm.round(1).head())

# Verificación: la suma de CPM de cada muestra es exactamente 1.000.000
assert np.allclose(cpm.sum(axis=0).values, 1_000_000)
# Comparación: CPM permite comparar lo que los counts sin procesar no permitían
tumor_avg = cpm[["tumor1","tumor2","tumor3"]].mean(axis=1)
normal_avg = cpm[["normal1","normal2","normal3"]].mean(axis=1)
# BRCA1 es efectivamente mayor en tumor
assert tumor_avg["BRCA1"] > normal_avg["BRCA1"]

Ahora podemos comparar entre muestras. Sin embargo, la comparación entre genes aún no es posible. El CPM de GAPDH sigue siendo 100 veces mayor que el CPM de BRCA1. La puntuación Z (z-score) resuelve este problema.

🔎 ¿Por qué el CPM es tan simple? (Draw — normalización) Es la forma más sencilla de normalización, que consiste en escalar a un tamaño estándar. En la práctica, se utilizan métodos más sofisticados, como TPM, TMM y la mediana de proporciones de DESeq2. La base conceptual es la misma: hacer que los datos sean comparables.

🤔 Pregunta para la reflexión ¿Por qué es más seguro usar CPM que usar directamente los recuentos brutos? Si el número total de lecturas de dos muestras es de 50 millones y 10 millones, ¿cómo se puede comparar la expresión del gen A utilizando los recuentos brutos? (Pista: no se puede).

Paso 3: Normalización con la puntuación Z (z-score) ★ (normalización: comparación entre genes)

✍️ Sección para completar directamente. Componente = normalización. Se estandariza cada gen con su propia media y desviación estándar para eliminar las diferencias de escala entre genes.

Para cada fila de genes:

z = (value - mean) / std

De esta manera, se expresa cuántas desviaciones estándar se aparta la expresión de cada gen de su valor medio. Esto permite comparar GAPDH y BRCA1 en la misma escala.

def to_zscore(cpm):
    """Z-score por gen tras transformación logarítmica, una práctica habitual."""
    log_cpm = np.log2(cpm + 1)                        # +1 evita log(0)
    means = log_cpm.mean(axis=1)
    stds = log_cpm.std(axis=1, ddof=0)
    # Los genes con desviación 0 (todos los valores iguales) se fijan de forma segura en 0
    stds = stds.replace(0, 1)
    return log_cpm.sub(means, axis=0).div(stds, axis=0)

zscore = to_zscore(cpm)
print(zscore.round(2))

# Verificación: para cada gen, media ≈ 0 y desviación ≈ 1
row_means = zscore.mean(axis=1)
row_stds = zscore.std(axis=1, ddof=0)
assert np.allclose(row_means.values, 0, atol=1e-9)
# Comprueba solo genes con variación, excluyendo la división por desviación 0
active_genes = cpm.std(axis=1) > 0
assert np.allclose(row_stds[active_genes].values, 1, atol=0.05)
# Por diseño, BRCA1 debe ser positivo en tumor y negativo en normal
assert (zscore.loc["BRCA1", ["tumor1","tumor2","tumor3"]] > 0).all()
assert (zscore.loc["BRCA1", ["normal1","normal2","normal3"]] < 0).all()

Ahora, todos los genes están en la misma escala (media 0, desviación estándar 1). Los colores del mapa de calor se distribuyen libremente para cada gen. Esta estandarización es la clave para lograr una mayor densidad de información en el mapa de calor.

🤔 Pregunta para la autoexplicación ¿Por qué aplicamos log2(cpm + 1) antes de calcular el puntaje Z? ¿Qué pasaría si simplemente aplicáramos el puntaje Z a los recuentos de millones de lecturas (CPM) sin procesar? (Pista: la expresión se distribuye de forma natural en una escala logarítmica, y los valores altos pueden dominar la desviación estándar como valores atípicos).

Paso 4 de la creación: seleccionar solo los genes con alta variabilidad (proporcionado)

Para que el mapa de calor sea informativo, debemos eliminar los genes que no muestran cambios. Los genes que siempre tienen el mismo valor, independientemente de las condiciones, solo dibujan líneas blancas y añaden ruido.

def select_top_variable(zscore, top_n=50):
    """Devuelve los top_n genes con mayor intervalo de z-score."""
    # Usa el rango (máx.-mín.) en vez de la varianza para priorizar cambios intensos en ambos sentidos
    variability = zscore.max(axis=1) - zscore.min(axis=1)
    top_genes = variability.sort_values(ascending=False).head(top_n).index
    return zscore.loc[top_genes]

top = select_top_variable(zscore, top_n=8)
print(top.index.tolist())

# Verificación: solo quedan ocho genes
assert len(top) == 8
# Los genes housekeeping, con pocos cambios, deben quedar fuera de los ocho primeros
housekeeping = {"GAPDH", "ACTB", "B2M", "HPRT1", "TBP"}
selected_housekeeping = set(top.index) & housekeeping
# Se excluye la mayoría; alguno puede aparecer por ruido
assert len(selected_housekeeping) <= 2

Paso 5 de creación: Clustering jerárquico ★ (clustering-hierarchical)

✍️ Sección para completar directamente. Componente = Clustering jerárquico. Calcula la distancia entre dos vectores y los agrupa secuencialmente según su proximidad para construir un árbol.

Idea: Dado que cada muestra (o gen) es un vector, se puede medir la distancia entre los vectores. Se agrupan los dos más cercanos en un grupo, y luego se trata a ese grupo como un único vector para continuar agrupando. El resultado es un dendrograma, un árbol jerárquico.

Este método está disponible de forma completa en scipy.

🔎 ¿Qué es el clustering jerárquico? (Draw — clustering-hierarchical) "Comienza con cada punto como su propio grupo → fusiona repetidamente los dos grupos más cercanos → hasta que quede uno solo." Al podar el árbol resultante, se obtiene el número deseado de clústeres. A diferencia de k-means, no es necesario definir previamente el número de clústeres.

from scipy.cluster.hierarchy import linkage, dendrogram

def cluster_order(data, axis="rows", method="average", metric="euclidean"):
    """Devuelve índices reordenados mediante clustering jerárquico."""
    matrix = data.values if axis == "rows" else data.values.T
    Z = linkage(matrix, method=method, metric=metric)
    # leaves_list contiene el orden óptimo de las hojas originales del árbol
    from scipy.cluster.hierarchy import leaves_list
    order = leaves_list(Z)
    if axis == "rows":
        return data.index[order]
    return data.columns[order]

# Reordena muestras y genes según sus respectivos clusters
sample_order = cluster_order(top, axis="cols")
gene_order   = cluster_order(top, axis="rows")

reordered = top.loc[gene_order, sample_order]
print("Orden de las muestras:", list(sample_order))
print("Orden de los genes:", list(gene_order))

# Verificación: la suma de la tabla reordenada coincide con la original
assert np.allclose(reordered.values.sum(), top.values.sum())
# Por diseño, las tres muestras tumorales y las tres normales deben agruparse por separado
sample_names = list(sample_order)
tumor_positions = [i for i, s in enumerate(sample_names) if s.startswith("tumor")]
normal_positions = [i for i, s in enumerate(sample_names) if s.startswith("normal")]
# Las muestras tumorales deben ocupar posiciones consecutivas
assert max(tumor_positions) - min(tumor_positions) == 2

Los tres tumores se agruparon automáticamente. Aunque no indicamos al algoritmo de agrupamiento que las muestras estuvieran etiquetadas como tumor/normal, los propios datos revelan ese grupo. Ese es el valor del agrupamiento: descubre la estructura sin etiquetas.

🤔 Autodescripción mediante indicaciones Seleccioné method="average" y metric="euclidean". ¿Qué pasaría si los cambiara por method="ward"? ¿Y si usara metric="correlation"? Cambie los valores y compare los resultados.

Paso 6: ensamblaje de seaborn.clustermap ★

La función clustermap de Seaborn genera simultáneamente el agrupamiento y el mapa de calor. Dado que ya conocemos los principios que aplicamos manualmente en los pasos anteriores, podemos comprender con precisión qué hace cada parámetro de esta función.

import matplotlib
matplotlib.use("Agg")   # Compatible con Pyodide y servidores
import matplotlib.pyplot as plt
import seaborn as sns

def plot_rnaseq_heatmap(counts_df, top_n=50, save_path=None):
    """Pipeline completo: normalización → máxima variación → clustermap."""
    cpm = to_cpm(counts_df)
    z = to_zscore(cpm)
    top = select_top_variable(z, top_n=top_n)

    g = sns.clustermap(
        top,
        cmap="RdBu_r",           # Rojo-azul, distingue claramente positivo y negativo
        center=0,                # Centra en 0 (la media)
        vmin=-2, vmax=2,         # Intervalo de color fijo
        method="average",
        metric="euclidean",
        figsize=(10, max(5, top_n * 0.15)),
        cbar_kws={"label": "z-score (log2 CPM)"},
    )
    if save_path:
        g.savefig(save_path, dpi=300, bbox_inches="tight")
    return g

# Ejecución en Colab/Jupyter
g = plot_rnaseq_heatmap(counts, top_n=8, save_path=None)
plt.close("all")

# Verificación: el grid devuelto es un objeto clustermap de seaborn
from seaborn.matrix import ClusterGrid
assert isinstance(g, ClusterGrid)
# Los datos reordenados tienen la forma esperada
assert g.data2d.shape == (8, 6)

Esta única función ensambla todas las etapas de normalización, agrupación y mapeo de colores que hemos creado. Si conoce los principios de cada capa, podrá ajustar los parámetros con confianza.

Unir las piezas: la clase de panel de control completa

Agrupamos todo el pipeline en una sola clase.

class RNAseqDashboard:
    def __init__(self, counts_df: pd.DataFrame):
        self.counts = counts_df
        self.cpm = to_cpm(counts_df)
        self.zscore = to_zscore(self.cpm)

    def summary(self) -> pd.DataFrame:
        return pd.DataFrame({
            "n_genes": [len(self.counts)],
            "n_samples": [self.counts.shape[1]],
            "total_reads_mean": [self.counts.sum(axis=0).mean()],
            "top_variable_gene": [
                (self.zscore.max(axis=1) - self.zscore.min(axis=1)).idxmax()
            ],
        })

    def plot(self, top_n=50, save_path=None):
        return plot_rnaseq_heatmap(self.counts, top_n=top_n, save_path=save_path)

dash = RNAseqDashboard(counts)
print(dash.summary())

# Verificación
assert dash.summary()["n_genes"].iloc[0] == 13
assert dash.summary()["n_samples"].iloc[0] == 6
# El gen de mayor variación debe mostrar una diferencia real entre condiciones: BRCA1, MYC o KRAS
assert dash.summary()["top_variable_gene"].iloc[0] in {"BRCA1", "MYC", "KRAS"}

Esta herramienta es una versión simplificada del mapa de calor de expresión de RNA-seq que se ve con frecuencia en los artículos científicos actuales. Las herramientas prácticas (como pheatmap en R o ComplexHeatmap) simplemente añaden anotaciones de color, funciones de distancia personalizadas y superposiciones de pruebas estadísticas a esta base.

Existen otros enfoques (reflexión multicanal)

• TPM vs. CPM: CPM solo normaliza según el tamaño de la muestra. TPM (Transcripts Per Million) también normaliza según la longitud del gen, lo que lo hace más preciso para las comparaciones entre genes. Cuándo usar cada uno: Solo para comparaciones entre muestras = CPM / Para comparaciones absolutas entre genes = TPM.
• Normalización de DESeq2: Si el objetivo es el análisis de expresión diferencial, se deben utilizar las normalizaciones median-of-ratios de DESeq2 o TMM (edgeR) en lugar de CPM/TPM, ya que se ajustan mejor a los supuestos de las pruebas estadísticas.
• Alternativas a k-means: El agrupamiento jerárquico requiere una memoria de O(n²), lo que lo hace inviable cuando hay decenas de miles de genes. En esos casos, se puede sustituir por k-means o UMAP + HDBSCAN. El objetivo de descubrir estructuras sin etiquetas sigue siendo el mismo.
• Selección de la función de distancia: euclidean es sensible al valor absoluto de los datos. La correlación se centra únicamente en la forma del patrón (ascensos y descensos). Tras la normalización con la puntuación z, ambos métodos son similares, pero conocer las alternativas permite elegir según la situación.
• Riesgo con tamaños de muestra pequeños: Si hay muy pocas muestras, como 3+3, el agrupamiento es inestable debido al ruido. Al incluir este tipo de mapas de calor en un artículo, se debe indicar en la interpretación de los resultados que el tamaño de la muestra no es grande.

Punto clave: "La normalización permite la comparabilidad, el agrupamiento descubre la estructura y el mapa de calor capta la atención". Si se conoce el papel de cada una de estas tres capas, el proceso desde los datos originales hasta el descubrimiento será reproducible. Lo que han creado es la materialización de ese proceso.

Siguiente paso (enlace de salida en la parte inferior)

• Más detalles sobre los diversos métodos de normalización → Principios de normalización
• Por qué el agrupamiento jerárquico encuentra estructuras sin etiquetas → Agrupamiento jerárquico
• Convertir los resultados de qPCR creados anteriormente en un mapa de calor → Automatización del análisis de qPCR

Pruébelo usted mismo (problema independiente)

1. Extensión TPM: Escriba el to_tpm(counts, gene_lengths) que normaliza a TPM en lugar de CPM, recibiendo además la información de la longitud del gen (archivo CSV con las longitudes).
2. Columna de anotación: En clustermap, agregue una barra de colores para las condiciones tumor/normal utilizando el parámetro col_colors. (Sugerencia: cree y pase una serie de colores para cada condición).
3. Solo genes seleccionados: Agregue una opción para que el usuario especifique directamente su lista de genes de interés (en lugar de usar la opción top_n).
4. Desafío: Mapa de calor de correlación: Calcule la matriz de correlación gen-gen (correlación de Pearson) y genere un mapa de calor agrupado (clustermap). Utilícelo para descubrir patrones de coexpresión.

Resumen

Hemos abordado el problema de "crear mapas de calor informativos a partir de los resultados de RNA-seq" mediante tres componentes clave.

• La normalización permitió la comparación entre genes y muestras.
• El agrupamiento jerárquico descubrió automáticamente la estructura de los grupos tumor/normal sin necesidad de etiquetas previas.
• seaborn clustermap ensambló los resultados en un formato visual adecuado para su inclusión en publicaciones científicas.

El "misterio" de esos hermosos mapas de calor que se ven en los artículos científicos se ha resuelto: son simplemente el resultado de tres capas que cumplen su función. Con el flujo de trabajo (pipeline) que ha creado, puede descubrir su propia historia en sus propios datos.

Este texto es un ejemplo educativo general. El análisis práctico de RNA-seq (DESeq2, edgeR, Bioconductor, etc.) incluye pruebas estadísticas, corrección de efectos de lote y análisis de enriquecimiento de GO, entre otros elementos. La versión detallada se puede añadir a esta estructura básica o delegar en bibliotecas validadas.