Perturbación génica in silico de modelos fundacionales de células individuales: Vista previa en 3 minutos del cribado de 200 genes esenciales en K562

Los experimentos de knockout génico con CRISPR requieren de 2 a 4 semanas, desde el cultivo celular y el diseño de gARN (capítulo 04), hasta la transducción lentiviral, la selección, la secuenciación y el análisis. Además, si no se observa un fenotipo interesante en la práctica, todo ese esfuerzo se pierde. Los modelos fundacionales de células individuales, como Geneformer y scGPT, abordan este problema desde una perspectiva fundamentalmente nueva. Simulan el knockout reduciendo artificialmente el rango de expresión de un gen específico en el espacio latente aprendido (disminución de rango) y realizan un cribado preliminar de candidatos experimentales midiendo cuánto se modifican los embeddings celulares (desplazamiento coseno). En este capítulo, desarrollamos a fondo una canalización de perturbación in silico utilizando como ejemplo 200 genes esenciales en K562 (una línea celular de leucemia mieloide crónica humana) y validamos su concordancia con los datos de Perturb-seq obtenidos en experimentos de laboratorio de Replogle, 2022.

📚 Recomendación de capítulos previos (muy recomendable)

Este capítulo es un análisis avanzado de la IA y la biología. Antes de comenzar, se recomienda encarecidamente que escuche y vea primero los siguientes capítulos de DryBench:

• DryBench ai-native #3: Transformadores y embeddings
• DryBench ai-native #5: Mecanismos de atención
• DryBench ai-native #7: Ingeniería de prompts

Si accede a este capítulo sin haber visto los anteriores, comenzará directamente con el código práctico sin una nueva explicación sobre la manipulación del espacio de embeddings del transformador, el enmascaramiento de la atención y el ajuste del rango de los tokens génicos, por lo que le resultará difícil seguirlo.

Ya aprendimos esto en DryBench

En DryBench ai-native #3, aprendimos que los embeddings de los transformadores constituyen un espacio latente aprendido; en #5, aprendimos que la atención aprende las relaciones entre elementos arbitrarios dentro de una secuencia; y en #7, aprendimos que es posible dirigir la distribución de salida de un LLM mediante prompts.

Los modelos fundacionales de células individuales aplican estos principios a los datos de expresión génica. Cada célula se representa como una "secuencia de rangos de expresión génica", y el transformador aprende los patrones de esta secuencia. Una vez entrenado, si se elimina un gen específico de la secuencia (simulación de knockout) o se fuerza su movimiento hacia arriba (simulación de sobreexpresión), la forma en que reaccionan los embeddings del modelo se correlaciona sorprendentemente con los resultados experimentales obtenidos en el laboratorio. Este capítulo convierte esa observación en una herramienta práctica de cribado.

Definición del problema de nivel avanzado

Escenario de I+D en la práctica: cribado de genes esenciales en K562

En el estudio Perturb-seq de Replogle (2022) [1], se midió experimentalmente una biblioteca CRISPR a escala genómica (11.258 KO génicos) en células K562. De este conjunto, se seleccionaron 200 genes esenciales:

• En el laboratorio (estándar): Biblioteca de 200 gRNA → 4 a 8 semanas.
• In silico (este artículo): KO in silico de 200 genes en Geneformer 95M → 3 a 10 minutos de cálculo + validación en el laboratorio solo de los 20 mejores. Reducción del tiempo y el costo por un factor de 20.
• Validación: Cálculo de la tasa de reproducción de los top-K utilizando los datos experimentales de Replogle (disponibles en SRA y GEO).

Indicadores objetivo:

• Perturbación in silico de 200 genes en menos de 30 minutos (con GPU de centro de datos de 24 GB de VRAM).
• Tasa de reproducción de los top-20 genes superior al 60 % en comparación con los datos experimentales de Perturb-seq.
• Visualización de las trayectorias de transición del estado celular por gen (desplazamiento en UMAP).
• Reproducibilidad de la ejecución: fijación de la semilla, versión del modelo y conjunto de datos.

Espectro de enfoques existentes

• Basados en correlación (ingenuo): Redes de correlación de expresión génica (WGCNA · GENIE3). Sin causalidad, poder predictivo limitado.
• CausalPath: Inferencia de redes de regulación génica. Válido solo en condiciones específicas.
• Trayectorias de aprendizaje profundo (PAGA · SCENIC): Aprendizaje de trayectorias celulares. Limitado para la simulación de perturbaciones.
• scGPT (Nature Methods 2024): Primer modelo fundacional con SOTA. Soporta perturbación in silico [2].
• Geneformer (Nature 2023): Codificación rank-value, predicción de redes génicas [3].
• scFoundation (2024): Escala de 100 millones de células [4].
• UCE (Universal Cell Embedding, 2024): Transversal a múltiples especies y tejidos [5].
• Nicheformer (2024): Modelo fundacional de espacio celular único.

Este artículo utiliza la línea base de pesos 95M de Geneformer + un conjunto de modelos scGPT.

Stack de herramientas e infraestructura requerida

Requisitos de infraestructura:

• Ejecución de perturbaciones in silico: GPU de centro de datos con 24 GB+ de VRAM (Geneformer 95M fp16). El cribado paralelo a gran escala requiere una estación de trabajo de centro de datos con 80 GB+ de VRAM (el acceso general no está disponible; se recomienda el acceso bajo demanda en la nube).
• Experimentos a pequeña escala: Una GPU de gama alta para consumidores (RTX 4090 de 24 GB) permite procesar lotes pequeños.
• RAM de al menos 32 GB (datos K562, incrustaciones, matrices de resultados).
• Disco: los pesos de Geneformer ocupan aproximadamente 4 GB, el conjunto de datos K562 ocupa aproximadamente 5 GB y Replogle Perturb-seq ocupa aproximadamente 10 GB (subconjunto).

Costo estimado para la reproducción por parte del estudiante: Si no se dispone de una GPU de alto rendimiento local, consulte las tarifas por hora de las instancias bajo demanda en la nube. El cribado de 200 genes tarda aproximadamente entre 30 y 60 minutos.

Implementación práctica del flujo de trabajo

Flujo completo:

Paso 1. Carga y control de calidad de los datos de células individuales de K562

Utilice el conjunto de datos Replogle 2022 o el subconjunto de K562 de CELLxGENE.

from pathlib import Path

import scanpy as sc
import anndata as ad
import numpy as np


def load_and_qc(
    h5ad_path: Path,
    min_genes_per_cell: int = 500,
    max_genes_per_cell: int = 8000,
    max_pct_mito: float = 15.0,
    max_pct_ribo: float = 50.0,
) -> ad.AnnData:
    """QC estándar del conjunto de datos K562."""
    adata = sc.read_h5ad(h5ad_path)
    print(f"Carga: células {adata.n_obs}, genes {adata.n_vars}")

    # Genes de mitocondria y ribosoma como etiquetas
    adata.var["mt"] = adata.var_names.str.startswith("MT-")
    adata.var["ribo"] = adata.var_names.str.startswith(("RPS", "RPL"))
    sc.pp.calculate_qc_metrics(adata, qc_vars=["mt", "ribo"], inplace=True)

    # Filtrado
    sc.pp.filter_cells(adata, min_genes=min_genes_per_cell)
    sc.pp.filter_cells(adata, max_genes=max_genes_per_cell)
    adata = adata[adata.obs["pct_counts_mt"] < max_pct_mito, :].copy()
    adata = adata[adata.obs["pct_counts_ribo"] < max_pct_ribo, :].copy()

    # Normalización · Transformación logarítmica
    sc.pp.normalize_total(adata, target_sum=1e4)
    sc.pp.log1p(adata)

    # HVG · Escalado (Geneformer usa recuentos crudos, por lo que se mantiene una capa separada)
    if "raw_counts" not in adata.layers:
        adata.layers["raw_counts"] = adata.X.copy()  # Para la entrada de Geneformer

    print(f"Después del QC: células {adata.n_obs}, genes {adata.n_vars}")
    return adata

Paso 2. Tokenización de Geneformer (codificación por rango)

Geneformer ordena los genes según su nivel de expresión en cada célula y los representa como una secuencia. Se necesita el ID de gen de Ensembl.

import torch


class GeneformerEmbedder:
    """Envoltura del modelo preentrenado Geneformer de 95M.

    Utiliza EmbExtractor y TranscriptomeTokenizer del repositorio oficial de Geneformer.
    """

    MODEL_VARIANTS = {
        "gf-6L-30M-i2048": "Geneformer 30M small, context 2048",
        "gf-12L-95M-i2048": "Geneformer 95M standard, context 2048",
        "gf-12L-95M-i4096": "Geneformer 95M extended, context 4096",
    }

    def __init__(
        self,
        device: str = "cuda",
        model_variant: str = "gf-12L-95M-i2048",
    ):
        # Geneformer utiliza su propio tokenizador + la estructura BertForMaskedLM
        # Para el código en producción, consulte los ejemplos de helical o del repositorio oficial
        from geneformer import TranscriptomeTokenizer, EmbExtractor
        self.tokenizer = TranscriptomeTokenizer(
            custom_attr_name_dict={"cell_type": "cell_type", "target_gene": "target_gene"},
            nproc=4,
        )
        self.extractor = EmbExtractor(
            model_type="Pretrained",
            num_classes=0,
            emb_mode="cell",       # Incrustación celular (similar al token CLS)
            filter_data=None,
            max_ncells=None,
            emb_layer=-1,
            forward_batch_size=8,
            nproc=4,
        )
        self.device = device
        self.model_variant = model_variant

    def tokenize_adata(self, adata: ad.AnnData, output_dir: Path) -> Path:
        """AnnData → Archivo de tokens de Geneformer (.dataset).

        ID de gen Ensembl (adata.var["ensembl_id"]) obligatorio.
        """
        if "ensembl_id" not in adata.var.columns:
            # Mapeo de símbolo → Ensembl (mygene · pyensembl)
            raise ValueError("Se requiere adata.var['ensembl_id']. Mapear previamente con mygene o pyensembl.")

        output_dir.mkdir(parents=True, exist_ok=True)
        # Conversión anndata → loom y tokenización (flujo de trabajo oficial de Geneformer)
        loom_path = output_dir / "cells.loom"
        adata.write_loom(loom_path)
        self.tokenizer.tokenize_data(
            data_directory=str(output_dir),
            output_directory=str(output_dir),
            output_prefix="tokenized",
            file_format="loom",
        )
        return output_dir / "tokenized.dataset"

    def extract_baseline_embeddings(self, tokenized_path: Path, output_dir: Path) -> np.ndarray:
        """Incrustaciones de referencia celular (antes de la perturbación)."""
        embs = self.extractor.extract_embs(
            model_directory=self.model_variant,
            input_data_file=str(tokenized_path),
            output_directory=str(output_dir),
            output_prefix="baseline_embs",
        )
        return np.asarray(embs)

Paso 3. Eliminación in silico (protocolo de descenso de rango)

Idea principal: desplazar el gen objetivo al último lugar en la secuencia de clasificación de la expresión génica celular, simulando así un efecto de eliminación real. Geneformer proporciona la API oficial InSilicoPerturber.

def in_silico_knockout(
    embedder: GeneformerEmbedder,
    tokenized_path: Path,
    target_gene_ensembl: str,
    output_dir: Path,
    max_ncells: int = 100,
) -> np.ndarray:
    """Knockout in silico de un gen específico → incrustaciones de células perturbadas.

    API oficial InSilicoPerturber de Geneformer [3].
    """
    from geneformer import InSilicoPerturber

    perturber = InSilicoPerturber(
        perturb_type="delete",             # delete: eliminación del gen, overexpress: desplazamiento hacia arriba
        perturb_rank_shift=None,            # None para delete
        genes_to_perturb=[target_gene_ensembl],
        combos=0,                            # Gen único
        anchor_gene=None,
        model_type="Pretrained",
        num_classes=0,
        emb_mode="cell",                     # Observar cambios en las incrustaciones celulares
        cell_emb_style="mean_pool",
        filter_data=None,
        cell_states_to_model=None,
        max_ncells=max_ncells,               # Subconjunto de células (ahorro computacional)
        emb_layer=-1,
        forward_batch_size=8,
        nproc=4,
    )

    perturbed_embs = perturber.perturb(
        model_directory=embedder.model_variant,
        input_data_file=str(tokenized_path),
        output_directory=str(output_dir),
        output_prefix=f"perturbed_{target_gene_ensembl}",
    )
    return np.asarray(perturbed_embs)

Paso 4. Cuantificación del desplazamiento del coseno

Cuantificar el cambio en el estado celular mediante la distancia coseno entre las representaciones vectoriales de las células antes y después de la inactivación génica.

def cosine_shift(baseline_embs: np.ndarray, perturbed_embs: np.ndarray) -> np.ndarray:
    """Desplazamiento coseno por célula. forma: (N_células,) devuelto."""
    baseline_norm = baseline_embs / (np.linalg.norm(baseline_embs, axis=1, keepdims=True) + 1e-8)
    perturbed_norm = perturbed_embs / (np.linalg.norm(perturbed_embs, axis=1, keepdims=True) + 1e-8)
    cos_sim = np.sum(baseline_norm * perturbed_norm, axis=1)
    return 1.0 - cos_sim  # distancia


def euclidean_shift(baseline_embs: np.ndarray, perturbed_embs: np.ndarray) -> np.ndarray:
    """Distancia euclidiana (indicador auxiliar)."""
    return np.linalg.norm(perturbed_embs - baseline_embs, axis=1)


def screening_by_shift(
    embedder: GeneformerEmbedder,
    tokenized_path: Path,
    baseline_embs: np.ndarray,
    target_genes: list[str],       # Lista de Ensembl ID
    output_dir: Path,
    max_ncells_per_gene: int = 100,
) -> dict[str, dict]:
    """Múltiples genes KO in silico secuencial → estadísticas de desplazamiento."""
    results = {}
    for i, gene in enumerate(target_genes):
        if i % 20 == 0:
            print(f"[{i}/{len(target_genes)}] {gene}")
        try:
            perturbed = in_silico_knockout(
                embedder, tokenized_path, gene, output_dir,
                max_ncells=max_ncells_per_gene,
            )
            cos = cosine_shift(baseline_embs[:len(perturbed)], perturbed)
            euc = euclidean_shift(baseline_embs[:len(perturbed)], perturbed)
            results[gene] = {
                "mean_cosine_shift": float(np.mean(cos)),
                "median_cosine_shift": float(np.median(cos)),
                "std_cosine_shift": float(np.std(cos)),
                "mean_euclidean_shift": float(np.mean(euc)),
                "n_cells": int(len(perturbed)),
            }
        except Exception as e:
            results[gene] = {"error": str(e)}

    return dict(sorted(
        results.items(),
        key=lambda x: -(x[1].get("mean_cosine_shift", 0.0) if "error" not in x[1] else 0.0),
    ))

Paso 5. Visualización de la trayectoria UMAP

Visualizar la dirección del cambio en el estado celular mediante la proyección conjunta de los embeddings celulares antes y después del knockout en un gráfico UMAP.

import umap
import matplotlib.pyplot as plt


def visualize_perturbation_trajectory(
    baseline_embs: np.ndarray,
    perturbed_embs: np.ndarray,
    gene_symbol: str,
    output_path: str = "perturbation_umap.png",
    n_arrows: int = 20,
) -> None:
    """Visualización del movimiento del estado celular antes y después del knock-out en UMAP."""
    combined = np.vstack([baseline_embs, perturbed_embs])
    reducer = umap.UMAP(n_components=2, n_neighbors=15, min_dist=0.1, random_state=42)
    embedding_2d = reducer.fit_transform(combined)

    n_cells = len(baseline_embs)
    baseline_2d = embedding_2d[:n_cells]
    perturbed_2d = embedding_2d[n_cells:]

    fig, ax = plt.subplots(figsize=(12, 8))
    ax.scatter(baseline_2d[:, 0], baseline_2d[:, 1],
               color="lightgray", alpha=0.5, s=20, label="baseline (WT)")
    ax.scatter(perturbed_2d[:, 0], perturbed_2d[:, 1],
               color="crimson", alpha=0.7, s=30, label=f"KO({gene_symbol})")

    # Flechas de movimiento de células de muestra
    n_arrows_actual = min(n_arrows, n_cells)
    idx = np.random.choice(n_cells, n_arrows_actual, replace=False)
    for i in idx:
        ax.annotate("", xy=perturbed_2d[i], xytext=baseline_2d[i],
                    arrowprops=dict(arrowstyle="->", color="black", alpha=0.4, lw=0.8))

    ax.set_title(f"Movimiento del estado celular: {gene_symbol} knockout in silico\n(n_células={n_células}, flechas={n_flechas_actual})")
    ax.set_xlabel("UMAP1")
    ax.set_ylabel("UMAP2")
    ax.legend()
    plt.tight_layout()
    plt.savefig(output_path, dpi=150)
    plt.close()

Paso 6. Validación de los datos obtenidos en el laboratorio de Replogle Perturb-seq

Descargar los datos de Perturb-seq a escala genómica de Replogle, 2022 [1], desde GEO y compararlos con las predicciones in silico.

import pandas as pd
from scipy.stats import spearmanr


def load_replogle_ground_truth(
    replogle_h5ad_path: Path,
    baseline_ctrl_label: str = "non-targeting",
) -> pd.DataFrame:
    \"\"\"Cargar efectos KO por gen de Perturb-seq en Replogle 2022.

    Descarga y procesamiento de la biblioteca esencial de K562 desde GEO GSE168191.
    Cuantificación del desplazamiento en el estado celular de cada target_gene en comparación con el control no dirigible.
    """
    adata = sc.read_h5ad(replogle_h5ad_path)
    # se requiere la columna target_gene
    if "target_gene" not in adata.obs.columns:
        raise ValueError("adata.obs['target_gene'] obligatorio")

    # Incrustación de células de control (ej: PCA)
    if "X_pca" not in adata.obsm:
        sc.pp.pca(adata, n_comps=50)

    ctrl_mask = adata.obs["target_gene"] == baseline_ctrl_label
    ctrl_centroid = adata.obsm["X_pca"][ctrl_mask].mean(axis=0)

    # Distancia entre el centroide de cada célula objetivo y el centroide de control
    effects = []
    for target in adata.obs["target_gene"].unique():
        if target == baseline_ctrl_label:
            continue
        target_mask = adata.obs["target_gene"] == target
        target_centroid = adata.obsm["X_pca"][target_mask].mean(axis=0)
        # Distancia euclidiana y distancia coseno entre centroides
        eucl_dist = float(np.linalg.norm(target_centroid - ctrl_centroid))
        cos_sim = float(np.dot(target_centroid, ctrl_centroid) /
                        (np.linalg.norm(target_centroid) * np.linalg.norm(ctrl_centroid) + 1e-8))
        effects.append({
            "target_gene": target,
            "measured_shift": eucl_dist,
            "measured_cosine": 1.0 - cos_sim,
            "n_cells": int(target_mask.sum()),
        })
    return pd.DataFrame(effects)


def compare_with_wetlab(
    in_silico_results: dict[str, dict],
    wet_lab_df: pd.DataFrame,
    top_k: int = 20,
    gene_symbol_mapping: dict[str, str] | None = None,  # Ensembl → symbol
) -> dict:
    """Comparación de predicción in silico con medición experimental."""
    is_records = []
    for ensembl, stats in in_silico_results.items():
        if "error" in stats:
            continue
        symbol = gene_symbol_mapping.get(ensembl, ensembl) if gene_symbol_mapping else ensembl
        is_records.append({
            "target_gene": symbol,
            "silico_shift": stats["mean_cosine_shift"],
        })
    is_df = pd.DataFrame(is_records)
    merged = is_df.merge(wet_lab_df, on="target_gene", how="inner")
    print(f"Genes comunes: {len(merged)}")

    if len(merged) < 3:
        return {"error": "Insuficientes genes comunes"}

    # Correlación de Spearman (tamaño del desplazamiento)
    rho, pval = spearmanr(merged["silico_shift"], merged["measured_cosine"])

    # Tasa de reproducción Top-K
    is_topk = set(merged.nlargest(top_k, "silico_shift")["target_gene"])
    wl_topk = set(merged.nlargest(top_k, "measured_cosine")["target_gene"])
    recall_at_k = len(is_topk & wl_topk) / min(top_k, len(merged))

    # Precision at top-K
    precision_at_k = len(is_topk & wl_topk) / len(is_topk)

    return {
        "n_common_genes": len(merged),
        "spearman_rho": float(rho),
        "spearman_p": float(pval),
        f"recall_at_top{top_k}": float(recall_at_k),
        f"precision_at_top{top_k}": float(precision_at_k),
        "is_topk_genes": list(is_topk),
        "wl_topk_genes": list(wl_topk),
        "overlap_genes": list(is_topk & wl_topk),
    }

Paso 7. Flujo de trabajo integrado · cribado de genes esenciales en K562

K562_ESSENTIAL_GENES_EXAMPLE = [
    # Ejemplo. En la práctica, seleccionar de DepMap · MAGeCK · base de datos de genes esenciales.
    "ENSG00000141510",  # TP53
    "ENSG00000012048",  # BRCA1
    "ENSG00000139618",  # BRCA2
    "ENSG00000186092",  # MYC (ejemplo)
    # ... 200 elementos
]


def full_perturbation_pipeline(
    k562_h5ad_path: Path,
    essential_gene_ensembls: list[str],
    replogle_ref_path: Path | None,
    output_dir: Path,
    device: str = "cuda",
) -> dict:
    """K562 gene esencial in silico screening · Validación con Perturb-seq."""
    output_dir.mkdir(parents=True, exist_ok=True)

    print("[1/6] Control de calidad de datos K562")
    adata = load_and_qc(k562_h5ad_path)

    print("[2/6] Cargando Geneformer")
    embedder = GeneformerEmbedder(device=device)

    print("[3/6] Tokenización · incrustación de línea base")
    tokenized = embedder.tokenize_adata(adata, output_dir / "tokens")
    baseline = embedder.extract_baseline_embeddings(tokenized, output_dir / "baseline")
    print(f"      baseline shape={baseline.shape}")

    print(f"[4/6] Screening in silico de {len(essential_gene_ensembls)} genes KO")
    silico_results = screening_by_shift(
        embedder, tokenized, baseline, essential_gene_ensembls,
        output_dir / "perturbations", max_ncells_per_gene=50,
    )
    import json
    with open(output_dir / "silico_shifts.json", "w") as f:
        json.dump(silico_results, f, indent=2, ensure_ascii=False)

    print("[5/6] Visualización de trayectorias UMAP de los genes principales")
    top_genes = [g for g in list(silico_results.keys())[:5] if "error" not in silico_results[g]]
    for gene in top_genes:
        perturbed_dir = output_dir / "perturbations" / f"perturbed_{gene}"
        # ... recargar o almacenar en caché los resultados perturbados desde el embedder
        pass  # (En producción, utiliza caching para visualizar sin recalcular)

    print("[6/6] Comparación con datos reales de Perturb-seq")
    validation = {}
    if replogle_ref_path:
        wet_lab_df = load_replogle_ground_truth(replogle_ref_path)
        validation = compare_with_wetlab(silico_results, wet_lab_df, top_k=20)
        with open(output_dir / "validation.json", "w") as f:
            json.dump(validation, f, indent=2, ensure_ascii=False)
        print(f"      Spearman ρ={validation.get('spearman_rho', 0):.3f}, "
              f"Recall@20={validation.get('recall_at_top20', 0):.2f}")

    return {"silico_ranking": silico_results, "validation": validation}

Rendimiento, costo y casos de fallo conocidos

Referencias de rendimiento (citación de evaluaciones comparativas públicas)

Costo estimado de reproducción para el usuario

• Costo de API: 0 (al usar GPU local).
• Para usuarios sin GPU de alto rendimiento, consulte la tarifa por hora bajo demanda en la nube (Geneformer 95M es posible con 24 GB de VRAM).
• Análisis de ~200 genes: aproximadamente 30～60 minutos.
• Descarga de datos: pesos de Geneformer 4 GB + subconjunto Perturb-seq de Replogle 5～10 GB.

5 casos de fallo conocidos (recopilación comunitaria y bibliográfica)

1. Discrepancia entre KO in silico y resultados de laboratorio (en tipos celulares específicos) Síntoma: Predicciones de KO in silico con valores sin sentido en tipos celulares no presentes en los datos de preentrenamiento (por ejemplo, subtipos específicos de neuronas cerebrales). Causa: Limitaciones de cobertura del preentrenamiento del modelo fundacional. Prevención: (a) Ajuste fino con datos que incluyen el tipo celular objetivo, (b) Ensemble de múltiples modelos fundacionales (Geneformer + scGPT), (c) Filtrado por confianza de predicción (si la magnitud del desplazamiento del embedding es muy pequeña), (d) K562 está incluido en el preentrenamiento y es relativamente estable. Fuente: Kedzierska et al. "Assessing the limits of zero-shot foundation models in single-cell biology." bioRxiv 2023 [6].

2. Errores de mapeo de IDs de Ensembl de genes Síntomas: Fallo en el mapeo del símbolo génico (TP53) y el ID de Ensembl (ENSG00000141510) → la KO in silico se aplica al gen incorrecto. Causa: Los símbolos génicos tienen alias; Geneformer está entrenado en una versión específica de Ensembl (por ejemplo, GRCh38.p13). Solución: (a) Mejorar la precisión del mapeo con pyensembl o mygene.info y especificar la versión; (b) fijar la versión de Ensembl; (c) registrar en el log y omitir los genes cuyo mapeo falle; (d) verificar la versión más reciente del reentrenamiento de Geneformer. Fuente: Issues de GitHub de Geneformer [7].

3. El efecto de lote interfiere con el desplazamiento in silico Síntomas: Al analizar células de diferentes lotes conjuntamente, el propio efecto de lote se interpreta como un desplazamiento, distorsionando el efecto génico. Causa: El control de calidad previo y la normalización no eliminan completamente el efecto de lote. Solución: (a) Corrección del efecto de lote (Harmony, scVI, Scanorama); (b) realizar la KO in silico únicamente dentro de un solo lote; (c) utilizar el desplazamiento de los genes de control (genes neutrales) como ruido de referencia para evaluar la relación señal-ruido; (d) verificar la reproducibilidad de los valores de desplazamiento en múltiples lotes. Fuente: Peidli S et al. "scPerturb: datos de perturbación de células individuales armonizados". Nat Methods 2024 [8].

4. Discrepancia entre el protocolo de caída de rango y la KO real Síntomas: La simulación de caída de rango difiere de los resultados reales de KO con gRNA-Cas9. Causa: La caída de rango implica simplemente un cambio en el orden de expresión, mientras que la KO real refleja proteínas, complejos, cascadas posteriores y dinámicas temporales. Solución: (a) Interpretar como señal fuerte únicamente los genes con gran desplazamiento; (b) probar varios modos de perturbación de Geneformer (eliminación, sobreexpresión, reducción de la expresión); (c) utilizar los resultados solo para la generación de hipótesis y realizar una validación experimental obligatoria; (d) cuantificar la consistencia mediante el conjunto de referencia Perturb-seq. Fuente: Discusión de Theodoris et al. Nature 2023 [3]; Kedzierska 2023 [6].

5. Carga por descarga y almacenamiento de conjuntos de datos a gran escala Síntomas: El conjunto de datos completo de Perturb-seq de Replogle 2022 ocupa decenas de GB, y el subconjunto esencial de K562 también entre 5 y 10 GB. Causa: Los datos de recuento sin procesar de células individuales, aunque sean matrices dispersas, tienen un gran tamaño. Solución: (a) utilizar subconjuntos (solo el subconjunto de 200 genes esenciales), (b) descargar solo las muestras necesarias de GEO, (c) almacenamiento comprimido (h5ad + zstd), (d) recursos en la nube académica (por ejemplo, CZI Chan Zuckerberg Initiative BioHub). Fuente: GEO GSE168191, Replogle, 2022 [1].

Ideas para ampliar

• Simulación de sobreexpresión: simular la sobreexpresión aumentando el rango en lugar de disminuirlo (factores de transcripción, genes supresores de tumores, etc.).
• Combinaciones de múltiples genes: inactivación simultánea de dos genes (cribado de candidatos para letalidad sintética).
• Predicción de la respuesta a fármacos: comparar el cambio celular tras la inactivación del objetivo del fármaco con datos reales de tratamiento farmacológico → reposicionamiento de fármacos (relacionado con los capítulos 07 y 08).
• Trasplante entre especies: aplicar patrones aprendidos en datos de ratón a datos humanos (utilizando UCE).
• Integración de scRNA-seq en vivo: transmisión en tiempo real de datos experimentales → comparación en tiempo real con predicciones in silico.
• Predicción detallada de la red HIF y las vías de señalización: combinaciones de inactivación de vías de señalización específicas.

Próximo capítulo

• Capítulo 04 crispr-guide-scoring: Diseño de gRNA para los genes seleccionados mediante el cribado in silico (introducción en el laboratorio de experimentos con K562).
• Capítulo 07 drug-target-gnn: Utilización de las predicciones de inactivación génica para la priorización de dianas farmacológicas.
• Capítulo 13 protein-design-multimodal: Diseño de proteínas activadoras artificiales en lugar de genes inactivados.
• Capítulo 14 bio-mcp-agent: Exposición de las perturbaciones in silico mediante la herramienta MCP → ejecución autónoma de "extraer los genes con mayor cambio en K562".

Referencias

1. Replogle JM, Saunders RA, Pogson AN, et al. "Mapping information-rich genotype-phenotype landscapes with genome-scale Perturb-seq." Cell 2022. https://www.cell.com/cell/fulltext/S0092-8674(22)00597-9 · GEO GSE168191
2. Cui H, Wang C, Maan H, et al. "scGPT: toward building a foundation model for single-cell multi-omics using generative AI." Nature Methods 2024. https://www.nature.com/articles/s41592-024-02201-0
3. Theodoris CV, Xiao L, Chopra A, et al. "Transfer learning enables predictions in network biology (Geneformer)." Nature 2023. https://www.nature.com/articles/s41586-023-06139-9
4. Hao M, Gong J, Zeng X, et al. "Large-scale foundation model on single-cell transcriptomics (scFoundation)." Nature Methods 2024. https://www.nature.com/articles/s41592-024-02305-7
5. Rosen Y, Roohani Y, Agarwal A, et al. "Universal Cell Embeddings: A Foundation Model for Cell Biology (UCE)." bioRxiv 2024. https://www.biorxiv.org/content/10.1101/2023.11.28.568918v2
6. Kedzierska KZ, Crawford L, Amini AP, Lu AX. "Assessing the limits of zero-shot foundation models in single-cell biology." bioRxiv 2023. https://www.biorxiv.org/content/10.1101/2023.10.16.561085
7. Geneformer GitHub Issues: https://huggingface.co/ctheodoris/Geneformer/discussions
8. Peidli S, Green TD, Shen C, et al. "scPerturb: harmonized single-cell perturbation data." Nature Methods 2024. https://www.nature.com/articles/s41592-023-02144-y
9. Geneformer HuggingFace: https://huggingface.co/ctheodoris/Geneformer
10. scGPT GitHub: https://github.com/bowang-lab/scGPT
11. CELLxGENE (CZI): https://cellxgene.cziscience.com/
12. scanpy: https://scanpy.readthedocs.io/
13. anndata: https://anndata.readthedocs.io/
14. Human Cell Atlas: https://www.humancellatlas.org/
15. UMAP-learn: https://umap-learn.readthedocs.io/
16. Helical SDK (integración de Geneformer, scGPT y UCE): https://github.com/helicalAI/helical
17. DepMap (base de datos de genes esenciales): https://depmap.org/
18. MAGeCK (análisis de cribado CRISPR): https://sourceforge.net/p/mageck/wiki/Home/
19. Corrección de lotes con Harmony: https://github.com/immunogenomics/harmony
20. scVI (corrección de lotes y modelo generativo): https://scvi-tools.org/