Agrupamiento de expresión — Encuentra grupos de coexpresión de genes mediante correlación y ordenación

Al finalizar este tema

Podrás implementar tu propio agrupamiento jerárquico para encontrar grupos de genes que se mueven juntos (coexpresión) en datos de expresión génica en múltiples muestras, combinando numpy, correlación y ordenación aprendidos en los libros de texto. Comprenderás los principios del funcionamiento interno de scikit-learn o scipy.cluster a través del código.

Este artículo es un ejemplo didáctico. Para el agrupamiento en la práctica, utiliza pandas/scanpy/scikit-learn. Aquí, abordaremos con precisión los conceptos internos de estas herramientas.

"¿Cuáles de estos 5000 genes se mueven juntos?" — No se puede ver con los ojos

Tienes datos de RNA-seq. Una matriz de expresión de 5000 genes × 30 muestras. Planteas las siguientes preguntas:

• ¿Qué genes se activan y desactivan juntos?
• ¿A qué función biológica corresponden estos grupos de coexpresión?

Es imposible comparar visualmente los vectores de 30 dimensiones de 5000 genes. Debemos encontrar y ordenar las relaciones automáticamente.

Enfoque ingenuo: Dibuja el perfil de expresión para cada gen. 5000 gráficos. De todos modos, no se verán las relaciones.

Enfoque real: Hay dos ejes:

1. Correlación: Cuantifica cuánto se mueven juntos los perfiles de expresión de dos genes. Correlación de Pearson.
2. Agrupamiento jerárquico: Agrupa y fusiona los pares más similares. El resultado es similar a un dendrograma (árbol genealógico).

Al combinar ambos, los grupos de genes que se mueven juntos se mostrarán en una estructura jerárquica natural.

De la caja negra a los componentes

Componente 1: Matriz de expresión y correlación

import numpy as np


def load_expression_matrix(path: str) -> tuple[np.ndarray, list[str], list[str]]:
    """
    Devuelve: (matrix, gene_names, sample_names)
    matrix shape: (n_genes, n_samples)
    """
    with open(path) as f:
        header = f.readline().strip().split("\t")
        sample_names = header[1:]

        gene_names = []
        rows = []
        for line in f:
            parts = line.strip().split("\t")
            gene_names.append(parts[0])
            rows.append([float(x) for x in parts[1:]])

    return np.array(rows), gene_names, sample_names


expr, genes, samples = load_expression_matrix("expression.tsv")
print(f"Genes: {len(genes)}, muestras: {len(samples)}")
print(f"Forma de la matriz: {expr.shape}")

Ahora, la matriz de correlación. Cada elemento representa la correlación de Pearson entre dos genes.

def compute_correlation_matrix(expr: np.ndarray) -> np.ndarray:
    """
    expr shape: (n_genes, n_samples)
    Devuelve: matriz de correlación (n_genes, n_genes)
    """
    return np.corrcoef(expr)


corr = compute_correlation_matrix(expr)
print(f"corr shape: {corr.shape}")  # (5000, 5000)
print(f"Diagonal: {corr[0, 0]}")  # 1.0 (consigo mismo)

np.corrcoef trata cada fila como una variable. Es decir, si expr tiene la forma (gen, muestra), el resultado es una correlación gen-gen.

Parte 2: Matriz de distancia

En el agrupamiento, se utiliza la distancia en lugar de la similitud. La correlación es una medida de similitud que varía de -1 a +1, por lo que se transforma en distancia mediante la siguiente conversión.

def correlation_to_distance(corr: np.ndarray) -> np.ndarray:
    """
    distancia = 1 - correlación.
    Correlación positiva fuerte → distancia pequeña.
    Correlación negativa (movimiento opuesto) → distancia grande.
    """
    return 1.0 - corr


dist = correlation_to_distance(corr)

Nota: Existen otras definiciones. sqrt(2 * (1 - corr)) es una definición consistente con la distancia euclidiana. Seleccione según el propósito.

Parte 3: Agrupamiento jerárquico (Enlace simple)

Enfoque aglomerativo: Inicialmente, cada gen es un clúster separado. Se repite el proceso de fusionar los dos clústeres más cercanos hasta que quede un solo clúster.

Es necesario definir la distancia entre dos clústeres. El enlace simple es la definición más sencilla: la distancia mínima entre pares de elementos de los dos clústeres.

def hierarchical_clustering_single(dist: np.ndarray) -> list[tuple[int, int, float]]:
    """
    Devuelve eventos de fusión [(cluster_a, cluster_b, merge_distance), ...].
    Tras cada fusión, el id nuevo es el número original de clústeres más el índice del evento.
    """
    n = dist.shape[0]
    active_clusters = {i: [i] for i in range(n)}
    cluster_distances = {(i, j): dist[i, j] for i in range(n) for j in range(i + 1, n)}

    events: list[tuple[int, int, float]] = []
    next_id = n

    while len(active_clusters) > 1:
        best_pair = min(cluster_distances, key=cluster_distances.get)
        i, j = best_pair
        merge_dist = cluster_distances[best_pair]
        events.append((i, j, merge_dist))

        new_cluster = active_clusters[i] + active_clusters[j]
        del active_clusters[i]
        del active_clusters[j]
        active_clusters[next_id] = new_cluster

        # Distancias entre el clúster nuevo y los restantes (enlace simple = mínimo)
        new_distances = {}
        for existing_id in active_clusters:
            if existing_id == next_id:
                continue
            existing = active_clusters[existing_id]
            min_d = min(
                dist[a, b] for a in new_cluster for b in existing
            )
            new_distances[(min(existing_id, next_id), max(existing_id, next_id))] = min_d

        cluster_distances = {
            k: v for k, v in cluster_distances.items()
            if i not in k and j not in k
        }
        cluster_distances.update(new_distances)

        next_id += 1

    return events

Complejidad temporal: Con una implementación ingenua, es O(n³). Con 5.000 genes, esto sería difícil de manejar, y en la práctica se utiliza un algoritmo de O(n² log n). Este tutorial tiene como objetivo comprender el concepto, por lo que se utiliza una implementación ingenua.

Parte 4: Extracción de clústeres de los resultados

Extraiga un número específico de clústeres de la lista de eventos de fusión.

def cut_dendrogram(events: list[tuple[int, int, float]], n_leaves: int, num_clusters: int):
    parent = list(range(n_leaves + len(events)))

    def find(x):
        while parent[x] != x:
            parent[x] = parent[parent[x]]
            x = parent[x]
        return x

    def union(x, y, new_id):
        rx, ry = find(x), find(y)
        parent[rx] = new_id
        parent[ry] = new_id

    # Ignora las últimas (num_clusters - 1) fusiones y aplica union a las anteriores
    events_to_apply = events[:-num_clusters + 1] if num_clusters > 1 else events

    next_id = n_leaves
    for a, b, _ in events_to_apply:
        union(a, b, next_id)
        next_id += 1

    clusters = {}
    for leaf in range(n_leaves):
        root = find(leaf)
        clusters.setdefault(root, []).append(leaf)

    return list(clusters.values())


cluster_lists = cut_dendrogram(events, n_leaves=len(genes), num_clusters=10)
for i, members in enumerate(cluster_lists):
    print(f"Clúster {i}: {len(members)} genes")

Union-Find: estructura de datos que permite realizar operaciones de unión y búsqueda en tiempo O(α(n)) ≈ O(1). Es útil incluso en conjuntos de datos grandes.

Ordenación y mapa de calor ordenado

El resultado visual del agrupamiento se observa en el mapa de calor reordenado. Los genes de cada clúster se ordenan de manera que formen filas adyacentes.

def cluster_order(cluster_lists: list[list[int]]) -> list[int]:
    """Mantiene el orden original dentro de cada clúster y ordena por clúster."""
    result = []
    for members in cluster_lists:
        result.extend(sorted(members))
    return result


new_order = cluster_order(cluster_lists)
reordered_expr = expr[new_order]

Al graficar esta matriz reordenada como un mapa de calor, cada clúster aparece como un bloque distinto.

import matplotlib.pyplot as plt


def plot_heatmap(matrix: np.ndarray, ax=None) -> None:
    if ax is None:
        _, ax = plt.subplots(figsize=(6, 8))
    im = ax.imshow(matrix, cmap="RdBu_r", aspect="auto")
    plt.colorbar(im, ax=ax)


fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 8))
plot_heatmap(expr, ax1)
ax1.set_title("Original")
plot_heatmap(reordered_expr, ax2)
ax2.set_title("Clustered")
plt.tight_layout()
plt.show()

Esta comparación antes y después de la reorganización es una visualización clave que demuestra el poder del agrupamiento.

Desvanecimiento: tres espacios en blanco para que los completes

Espacio en blanco 1: Enlace completo / Enlace promedio

Implementa enlace completo (la distancia máxima entre dos pares de elementos de clúster) o enlace promedio (la distancia promedio) en lugar del enlace único.

def hierarchical_clustering_complete(dist: np.ndarray) -> list[tuple[int, int, float]]:
    # Casi idéntico al enlace simple: sustituye min por max
    # ...
    # Cálculo de la distancia entre clústeres:
    # TODO: sustituir min(dist[a, b] for ...) por max(dist[a, b] for ...)
    pass

Pista: Parametrizar solo una función. linkage: Callable[[list[float]], float] = min or max or (lambda xs: sum(xs) / len(xs)).

Espacio en blanco 2: Gen representativo del clúster

En cada clúster, el gen que tiene la correlación promedio más alta con los demás miembros = gen representativo.

def find_hub_genes(
    corr: np.ndarray,
    clusters: list[list[int]],
    gene_names: list[str]
) -> list[tuple[str, float]]:
    """Gen representativo de cada clúster y su correlación media."""
    hubs = []
    for members in clusters:
        # TODO: calcular la correlación media de cada miembro con los demás
        # Seleccionar como hub el gen con el valor más alto
        pass
    return hubs

Pista:

best_score = -1
best_gene = None
for m in members:
    score = np.mean([corr[m, other] for other in members if other != m])
    if score > best_score:
        best_score = score
        best_gene = gene_names[m]

Espacio en blanco 3: anotación funcional (integración con herramientas externas)

Asigna la lista de genes del clúster a vías KEGG o términos GO. Aquí, la asignación se carga desde un archivo CSV.

def enrich_clusters(
    clusters: list[list[str]],
    gene_to_pathway_csv: str
) -> dict:
    """
    Las tres vías más frecuentes en cada clúster.
    """
    # TODO: cargar CSV (columnas gene, pathway)
    # Calcular la frecuencia de cada vía en cada clúster
    # Devolver las tres vías principales de cada clúster
    pass

Pista: from collections import Counter; counter = Counter(); for g in members: counter.update(gene_pathways.get(g, [])).

Reflexiones: diferencias con las herramientas de coexpresión del mundo real

scipy.cluster.hierarchy: Las herramientas del mundo real utilizan la función linkage de este módulo. Es un algoritmo en C con una complejidad de O(n² log n). Es más de 100 veces más rápido que tu implementación ingenua de O(n³).

WGCNA: Es una herramienta estándar para redes de coexpresión génica. Utiliza conceptos sofisticados como la matriz de adyacencia firmada, el umbral suave y la preservación de módulos. Tu herramienta es una versión simplificada de esto.

scanpy: Es una pila estándar para el análisis de ARN de células individuales. El agrupamiento suele realizarse con el algoritmo de Leiden, un enfoque basado en grafos, que es diferente del agrupamiento jerárquico que has utilizado.

Efecto de lote: Los datos de expresión del mundo real suelen provenir de diferentes lotes experimentales. El efecto de lote puede contaminar la señal de coexpresión real. Herramientas como Combat y Harmony resuelven este problema.

Corte dinámico de árboles: Tu cut_dendrogram especifica el número de clústeres. Las herramientas de WGCNA del mundo real utilizan un algoritmo que corta el árbol dinámicamente en función de la forma del dendrograma (corte dinámico de árboles).

Proyectos de ampliación

1. Reimplementación con scipy: Reemplaza tu agrupamiento ingenuo con scipy.cluster.hierarchy.linkage y compara el rendimiento.

2. Visualización del dendrograma: Visualiza el árbol con scipy.cluster.hierarchy.dendrogram. Convierte tu lista de eventos al formato de scipy.

3. Enriquecimiento de términos GO: Realiza el enriquecimiento de términos GO para cada clúster con gseapy.

4. Panel de Streamlit: Crea un panel en el que los usuarios puedan ajustar el número de clústeres mediante un control deslizante, y el mapa de calor y los genes representativos se actualicen en tiempo real.

Mapa de los componentes de este ejemplo

• [F] numpy: Manipulación de la matriz de expresión y la matriz de correlación. Utilizado en el mundo real np.corrcoef.
• [F] Correlación: Se utiliza la correlación de Pearson como medida cuantitativa de la similitud entre dos genes.
• [F] Ordenación: Reordenación de los clústeres para crear bloques en el mapa de calor. Extracción de los resultados con Union-Find.
• [W] matplotlib: visualización de mapas de calor (se proporciona el script completo).

[F] = usted lo implementa / [W] = se proporciona el código completo.