Automatización del análisis de qPCR: calcula ΔΔCt y el cambio de expresión sin esfuerzo

Al finalizar este tema

Podrás crear tu propia secuencia de comandos que, a partir de un archivo CSV con los resultados de un experimento de qPCR, calcule automáticamente el ΔΔCt, el cambio de expresión y los intervalos de confianza, combinando los conceptos de pandas groupby, propagación de errores y entrada/salida de archivos CSV que aprendiste en clase. La tarea de calcular manualmente cada celda en Excel se reemplazará por una simple llamada a una función.

Este artículo es un ejemplo didáctico. El análisis de qPCR es una tarea común en los laboratorios de biología molecular, por lo que lo hemos utilizado como tema.

"Un solo error en una celda puede arruinarlo todo" — Los riesgos de Excel

Después de realizar un experimento de qPCR, obtendrás un archivo CSV como este.

sample,gene,ct
control_rep1,BRCA1,25.3
control_rep1,GAPDH,18.2
control_rep2,BRCA1,25.5
control_rep2,GAPDH,18.4
treated_rep1,BRCA1,22.1
treated_rep1,GAPDH,18.3
treated_rep2,BRCA1,22.3
treated_rep2,GAPDH,18.5
...

Lo que quieren hacer es lo siguiente:

1. Para cada muestra, ΔCt = Ct del gen de interés − Ct del gen de referencia (GAPDH)
2. Calcular la ΔCt media para cada condición entre las réplicas
3. ΔΔCt = ΔCt media del grupo de tratamiento − ΔCt media del grupo de control
4. Cambio relativo = 2^(-ΔΔCt)
5. Calcular la desviación estándar del cambio relativo mediante la propagación de errores

Están familiarizados con cómo hacerlo en Excel: VLOOKUP, AVERAGEIF, seleccionar celdas manualmente, arrastrar fórmulas, mover los resultados a una nueva hoja y volver a calcular...

Pero a menudo ocurren cosas como estas:

• Se selecciona una celda incorrecta, por lo que los datos del grupo de control se incluyen en el cálculo del grupo de tratamiento.
• El primer experimento funciona bien, pero en el segundo experimento, el orden de las columnas cambia, por lo que las fórmulas no funcionan.
• Se elimina una de las tres réplicas porque se considera un valor atípico, pero no se registra cuál se eliminó, por lo que no se puede reproducir.
• El cambio relativo en la presentación y el cambio relativo en la figura del artículo son ligeramente diferentes, pero no se puede encontrar dónde se produjo la discrepancia.

Estos errores han sido la causa de la retractación de varios artículos. Una sola celda puede arruinar un artículo.

En este artículo, crearemos una canalización que elimine la posibilidad de cometer estos errores. Simplemente introduzca el archivo CSV y obtendrá los resultados. No se requiere intervención manual en el proceso.

Primero, veamos el producto final (ejecutemos la caja negra primero)

Así es como se utilizará la herramienta que vamos a crear.

result = analyze_qpcr(
    csv_path="qpcr_data.csv",
    reference_gene="GAPDH",
    control_condition="control",
)
print(result)

=== Resultados del análisis qPCR ===
condition   gene     mean_ddct   fold_change   fc_std
control     BRCA1    0.00        1.00          0.15
control     TP53     0.00        1.00          0.11
treated     BRCA1    -3.05       8.28          0.72
treated     TP53     -1.72       3.29          0.31

Archivo guardado: qpcr_result.csv (reproducible)

Cada vez que se utilizan los mismos datos experimentales, se obtienen los mismos resultados. Solo hay que introducir la ruta del archivo, y el resto lo hace el código automáticamente. No hay margen para errores.

¿De qué componentes está compuesto este programa (desglose de componentes)?

Pipeline automatizado de análisis qPCR
   ┌──────────────────────────────────────────────────┐
   │ [Entrada] Cargar CSV ───── pieza: entrada/salida CSV│ ← Se proporciona completa
   │              │                                      │
   │              ▼                                      │
   │ [Paso 1] Agrupar por condición y gen + media         │
   │        pieza: pandas groupby                         │ ← La construyes tú ★
   │              │                                      │
   │              ▼                                      │
   │ [Paso 2] Cálculo vectorial de ΔCt, ΔΔCt y fold change│
   │        pieza: aritmética vectorizada                 │ ← La construyes tú ★
   │              │                                      │
   │              ▼                                      │
   │ [Paso 3] Propagación de errores                      │
   │        pieza: regla de combinación de desviaciones   │ ← La construyes tú ★
   │              │                                      │
   │              ▼                                      │
   │ [Salida] CSV de resultados + gráfico                 │ ← Se proporciona completa
   └──────────────────────────────────────────────────┘

📌 Si no estás familiarizado con estos conceptos (enlace en la parte superior)

• pandas groupby — dividir, aplicar, combinar
• Propagación de errores — reglas para combinar desviaciones estándar
• CSV — formato estándar para datos tabulares

Los nuevos conceptos que crearemos son groupby, cálculo vectorial y propagación de errores. La entrada/salida CSV y los gráficos se proporcionan como herramientas completas. Son solo tres, lo que está dentro de los límites de la capacidad cognitiva.

Paso 1: Preparación de datos (proporcionado completo)

Creemos datos que simulen resultados reales de qPCR. En la práctica, estos datos se obtienen en formato CSV desde el equipo.

import pandas as pd
import numpy as np

# Números aleatorios reproducibles para simular ruido experimental
rng = np.random.default_rng(42)

def make_qpcr_data():
    rows = []
    # Escenario: frente al control, BRCA1 aumenta ~8 veces y TP53 ~3 veces en treated
    scenarios = {
        ("control", "BRCA1"): 25.0,
        ("control", "TP53"):  23.5,
        ("control", "GAPDH"): 18.3,
        ("treated", "BRCA1"): 22.0,   # Ct disminuye 3 = aumento aproximado de 8 veces
        ("treated", "TP53"):  21.8,   # Ct disminuye 1,7 = aumento aproximado de 3 veces
        ("treated", "GAPDH"): 18.4,   # El gen de referencia apenas cambia
    }
    for (condition, gene), base_ct in scenarios.items():
        for rep in range(1, 4):  # 3 replicates
            ct = base_ct + rng.normal(0, 0.15)  # Variación técnica
            rows.append({
                "sample": f"{condition}_rep{rep}",
                "condition": condition,
                "gene": gene,
                "ct": round(ct, 3),
            })
    return pd.DataFrame(rows)

df = make_qpcr_data()

# Verificación de la forma de los datos
assert len(df) == 18                                # 2 conditions × 3 genes × 3 replicates
assert set(df.columns) == {"sample", "condition", "gene", "ct"}
assert set(df["gene"].unique()) == {"BRCA1", "TP53", "GAPDH"}
assert (df["ct"] > 15).all() and (df["ct"] < 30).all()
print(df.head(6))

Esta tabla representa el formato de los datos originales del experimento de qPCR. Ahora, lo utilizaremos como datos de entrada.

Paso 2: Calcular el promedio por condición y gen ★ (pandas groupby)

✍️ Sección para completar manualmente. Componente = pandas groupby. Objetivo: calcular el promedio y la desviación estándar de Ct para cada combinación de (condición, gen) entre las réplicas.

Dado que tenemos 3 réplicas, debemos extraer el promedio y la desviación estándar de estas 3 réplicas. Esto debe hacerse para cada combinación de (condición, gen). Si se implementa de forma sencilla, se obtendrá un bucle for anidado como este.

def summarize_naive(df):
    result = []
    for condition in df["condition"].unique():
        for gene in df["gene"].unique():
            subset = df[(df["condition"] == condition) & (df["gene"] == gene)]
            result.append({
                "condition": condition,
                "gene": gene,
                "mean_ct": subset["ct"].mean(),
                "std_ct": subset["ct"].std(),
            })
    return pd.DataFrame(result)

summary_naive = summarize_naive(df)
assert len(summary_naive) == 6   # 2 × 3

Esto sí que da una respuesta. Pero si hay 10 condiciones y 100 genes, el doble bucle for repetirá el filtrado 1000 veces. Si los datos aumentan, se vuelve más lento.

pandas groupby reemplaza este patrón con una sola línea.

🔎 ¿Qué es groupby? (Diapositiva — pandas-groupby) "Agrupa las filas del mismo grupo, aplica una función a cada grupo y combina los resultados" (split-apply-combine). Es exactamente el mismo concepto que GROUP BY en SQL. Se procesa en vector sin usar bucles for, por lo que es mucho más rápido.

def summarize_ct(df):
    return (
        df.groupby(["condition", "gene"])["ct"]
          .agg(["mean", "std", "count"])
          .rename(columns={"mean": "mean_ct", "std": "std_ct", "count": "n"})
          .reset_index()
    )

summary = summarize_ct(df)
print(summary)

# Verificación: coincide con la versión ingenua
naive = summarize_naive(df).sort_values(["condition", "gene"]).reset_index(drop=True)
smart = summary.sort_values(["condition", "gene"]).reset_index(drop=True)
assert np.allclose(smart["mean_ct"].values, naive["mean_ct"].values)
assert (smart["n"] == 3).all()   # Confirma tres réplicas

El mismo resultado, en una sola línea. Esta es la potencia de groupby.

🤔 Indicación para la autoevaluación ¿Por qué es necesario calcular count a partir de .agg(["mean", "std", "count"])? ¿En la práctica, se puede garantizar que el número de réplicas sea siempre 3 en la tabla? (Pista: si se eliminan los valores atípicos o hay pocillos con errores de carga, el número de réplicas puede variar entre los grupos).

Paso 3: Calcular el vector ΔCt · ΔΔCt · cambio de pliegue ★

✍️ Sección para completar manualmente. Componente = aritmética vectorial. Objetivo: rellenar las columnas ΔCt/ΔΔCt/cambio de pliegue de la tabla resumen sin usar un bucle for.

Actualmente, nuestra tabla resumen tiene este aspecto.

condition  gene    mean_ct  std_ct
control    BRCA1   25.00    0.12
control    GAPDH   18.30    0.10
control    TP53    23.50    0.11
treated    BRCA1   22.00    0.14
treated    GAPDH   18.40    0.13
treated    TP53    21.80    0.12

Aquí debemos calcular ΔCt = mean_ct - mean_ct(GAPDH, mismas condiciones). Parece complicado porque debemos encontrar el GAPDH de cada fila y restarlo. Sin embargo, en pandas, podemos hacerlo todo de una vez con merge.

def add_delta_ct(summary, reference_gene="GAPDH"):
    # Extrae solo el gen de referencia
    ref = (
        summary[summary["gene"] == reference_gene]
        [["condition", "mean_ct", "std_ct"]]
        .rename(columns={"mean_ct": "mean_ct_ref", "std_ct": "std_ct_ref"})
    )
    # Une por condición
    merged = summary.merge(ref, on="condition", how="left")
    merged["dct"] = merged["mean_ct"] - merged["mean_ct_ref"]
    # Propaga la desviación estándar: sqrt(std_target^2 + std_ref^2)
    merged["dct_std"] = np.sqrt(merged["std_ct"]**2 + merged["std_ct_ref"]**2)
    return merged.drop(columns=["mean_ct_ref", "std_ct_ref"])

with_dct = add_delta_ct(summary)
print(with_dct[["condition", "gene", "mean_ct", "dct", "dct_std"]])

# Verificación: el ΔCt del gen de referencia es 0 porque se resta a sí mismo
assert np.allclose(
    with_dct[with_dct["gene"] == "GAPDH"]["dct"].values, 0.0
)
# ΔCt es el valor del gen de interés menos GAPDH
control_brca1_dct = with_dct[
    (with_dct["condition"] == "control") & (with_dct["gene"] == "BRCA1")
]["dct"].iloc[0]
# 25,0 - 18,3 ≈ 6,7, teniendo en cuenta el ruido experimental
assert 6.0 < control_brca1_dct < 7.5

Ahora, ΔΔCt es, para cada gen, ΔCt del grupo tratado − ΔCt del grupo de control. Se utiliza el mismo patrón de combinación.

def add_ddct_and_fold(with_dct, control_condition="control"):
    # Extrae solo el grupo de control
    ctrl = (
        with_dct[with_dct["condition"] == control_condition]
        [["gene", "dct", "dct_std"]]
        .rename(columns={"dct": "dct_ctrl", "dct_std": "dct_std_ctrl"})
    )
    merged = with_dct.merge(ctrl, on="gene", how="left")
    merged["ddct"] = merged["dct"] - merged["dct_ctrl"]
    # Propagación del error: supone independientes los dos ΔCt
    merged["ddct_std"] = np.sqrt(merged["dct_std"]**2 + merged["dct_std_ctrl"]**2)
    # fold change = 2^(-ΔΔCt)
    merged["fold_change"] = 2 ** (-merged["ddct"])
    # Propagación: sd de fc = fc * ln(2) * ddct_std, derivada de la transformación log
    merged["fc_std"] = merged["fold_change"] * np.log(2) * merged["ddct_std"]
    return merged.drop(columns=["dct_ctrl", "dct_std_ctrl"])

final = add_ddct_and_fold(with_dct)
print(final[["condition", "gene", "ddct", "fold_change", "fc_std"]])

# Verificación: en control, ΔΔCt es 0 y fold change es 1
ctrl_rows = final[final["condition"] == "control"]
assert np.allclose(ctrl_rows["ddct"].values, 0.0)
assert np.allclose(ctrl_rows["fold_change"].values, 1.0)
# Por diseño, BRCA1 en treated debe estar cerca de 8 veces
treated_brca1_fc = final[
    (final["condition"] == "treated") & (final["gene"] == "BRCA1")
]["fold_change"].iloc[0]
assert 5 < treated_brca1_fc < 12

Este código no contiene ningún bucle for. Los cálculos se aplican de forma vectorial a todo el dataframe.

🤔 Indicación de autoexplicación Explique por qué la fórmula de la desviación estándar del cambio de pliegue es fc * ln(2) * ddct_std. Se deriva de la diferenciación de la función f(x) = 2^(-x) y de la fórmula general de propagación de errores (σ_f ≈ |f'| · σ_x). (Pista: d/dx[2^(-x)] = -ln(2) · 2^(-x))

Paso 4: Profundizar en el principio de propagación de errores ★ (error-propagation)

Ya hemos utilizado la propagación de errores dos veces. Vamos a analizarlo explícitamente.

🔎 ¿Qué es la propagación de errores? (Diagrama — error-propagation) Supongamos que tenemos dos mediciones, A y B, con desviaciones estándar σA y σB, respectivamente. La suma o resta de estas dos mediciones se propaga como σ = √(σA² + σB²) (asumiendo independencia). En el caso de la multiplicación o división, las desviaciones estándar relativas se combinan. Para una función f(A), la desviación estándar se aproxima como σ ≈ |f'(A)| · σA. Estas tres reglas son fundamentales para trabajar con datos experimentales.

Apliquemos este principio a nuestro flujo de trabajo:

• ΔCt = Ct(target) - Ct(GAPDH) → Resta → σ_ΔCt = √(σ_target² + σ_GAPDH²)
• ΔΔCt = ΔCt(treated) - ΔCt(control) → Resta → σ_ΔΔCt = √(σ_treated² + σ_control²)
• Cambio de pliegue = 2^(-ΔΔCt) → Función no lineal → σ_fc ≈ |d/dx[2^(-x)]| · σ_ΔΔCt = fc · ln(2) · σ_ΔΔCt

# Comprobación manual
# Supuesto: σ_target=0,15, σ_GAPDH=0,10
# ΔCt σ = sqrt(0.15^2 + 0.10^2) ≈ 0.180
manual_dct_std = np.sqrt(0.15**2 + 0.10**2)
assert abs(manual_dct_std - 0.180) < 0.01
# Después, σ_ΔΔCt = sqrt(0,180^2 + 0,180^2) ≈ 0,255
manual_ddct_std = np.sqrt(manual_dct_std**2 + manual_dct_std**2)
assert abs(manual_ddct_std - 0.255) < 0.01
# Si fc=8, σ_fc = 8 * ln(2) * 0,255 ≈ 1,41
manual_fc_std = 8 * np.log(2) * 0.255
assert 1.3 < manual_fc_std < 1.5

Este cálculo manual debe tener una magnitud similar a los resultados de nuestro flujo de trabajo. El hecho de que la propagación de errores pueda verificarse manualmente es la fortaleza de este principio.

Uniendo las piezas: la función del flujo de trabajo completo

Combinamos los tres componentes en una sola función.

def analyze_qpcr(
    df: pd.DataFrame,
    reference_gene: str = "GAPDH",
    control_condition: str = "control",
) -> pd.DataFrame:
    """Datos qPCR sin procesar → tabla de ΔΔCt, fold change y propagación de errores."""
    # 1. Resumen por condición y gen
    summary = summarize_ct(df)
    # 2. ΔCt
    with_dct = add_delta_ct(summary, reference_gene=reference_gene)
    # 3. ΔΔCt + fold change
    result = add_ddct_and_fold(with_dct, control_condition=control_condition)
    # 4. Selecciona solo las columnas que se mostrarán
    return result[[
        "condition", "gene", "n", "mean_ct", "std_ct",
        "dct", "dct_std", "ddct", "ddct_std", "fold_change", "fc_std",
    ]].sort_values(["condition", "gene"]).reset_index(drop=True)

result = analyze_qpcr(df)
print(result)

# Control: fold change = 1 y desviación cercana a 0
ctrl = result[result["condition"] == "control"]
assert np.allclose(ctrl["fold_change"].values, 1.0)
# BRCA1 tratado: aproximadamente 8 veces por diseño
treated_brca1_fc = result[
    (result["condition"] == "treated") & (result["gene"] == "BRCA1")
]["fold_change"].iloc[0]
assert 5 < treated_brca1_fc < 12

Esta herramienta es una implementación del método de Livak (ΔΔCt, 2001) que estamos analizando hoy. Las herramientas prácticas (qbase+, PrimePCR Analysis, etc.) simplemente añaden características adicionales, como la corrección de la eficiencia, múltiples genes de referencia y la detección automática de valores atípicos. La estructura básica es la misma que la que acabas de crear.

Análisis detallado del rendimiento: ¿por qué es tan rápido?

import time

# Simulación grande: 100 genes, 20 condiciones y 4 réplicas
def big_data(n_genes=100, n_conditions=20, n_reps=4):
    genes = [f"G{i}" for i in range(n_genes)]
    conditions = ["control"] + [f"cond_{i}" for i in range(1, n_conditions)]
    rows = []
    for c in conditions:
        for g in genes + ["GAPDH"]:
            for r in range(n_reps):
                rows.append({
                    "sample": f"{c}_rep{r}",
                    "condition": c,
                    "gene": g,
                    "ct": 20 + rng.normal(0, 0.2),
                })
    return pd.DataFrame(rows)

df_big = big_data()
t0 = time.time()
result_big = analyze_qpcr(df_big)
elapsed = time.time() - t0
print(f"100 genes · 20 condiciones · 4 réplicas → {elapsed*1000:.1f} ms")
assert elapsed < 2.0  # Debe terminar en pocos segundos

Aunque el número de filas sea de miles, el procesamiento se realiza en milisegundos. Esto se debe a que se utiliza groupby + cálculo vectorial en lugar de bucles for.

¿Qué pasa si hacemos lo mismo en Excel? Tendríamos que crear una hoja de cálculo para cada gen, y volver a establecer las celdas de referencia para cada condición... Esto llevaría horas o incluso días. Además, si una sola celda es incorrecta, los resultados se verán silenciosamente afectados.

Hay otros caminos (reflexión multipaso)

• Corrección de eficiencia (método de Pfaffl): El método de Livak asume que la eficiencia de la PCR es del 100% (exactamente el doble en cada ciclo). En realidad, la eficiencia varía ligeramente según el cebador y las condiciones. El método de Pfaffl tiene en cuenta la eficiencia real de cada gen (E, entre 1,8 y 2,0). Cuándo y qué: cribado/estudio piloto = Livak / resultados críticos/gráficos para publicaciones = Pfaffl.
• Múltiples genes de referencia: Depender de un solo gen, como GAPDH, es arriesgado. El propio GAPDH puede variar según las condiciones. El promedio geométrico de múltiples genes de referencia (GeNorm, NormFinder) es más seguro.
• Detección automática de valores atípicos por pocillo: Si uno de los tres valores repetidos difiere significativamente de los otros dos (por ejemplo, >0,5 de diferencia en Ct), se marca automáticamente como un valor atípico. Si lo añadimos a nuestro flujo de trabajo, la reproducibilidad mejorará considerablemente.
• DataFrame vs. SQL: Hemos utilizado pandas, pero si los datos experimentales son muy grandes (por ejemplo, el archivo de experimentos de todo el laboratorio), es mejor migrar a DuckDB o PostgreSQL y procesarlos con SQL. El algoritmo sigue siendo el mismo.

Clave: "groupby en lugar de bucles for, aritmética vectorial en lugar de cálculos manuales, archivos en lugar de flujos". Estos tres principios determinan la reproducibilidad del análisis experimental. Lo que acabas de crear es la materialización de estos principios.

Próximos pasos (enlaces de salida en la parte inferior)

• Por qué la aritmética vectorial es mucho más rápida que los bucles for → Vectorización
• Reglas detalladas de propagación de errores (multiplicación, división, no lineal) → Principio de propagación de errores
• Para seguir con la visualización de los resultados → Aplicación Panel de mapas de calor de RNA-seq

Ponlo en práctica (problema independiente)

1. Filtro de valores atípicos: Añade un preprocesamiento que elimine automáticamente los valores repetidos que difieren en 0,5 Ct o más de los otros dos valores repetidos. También incluye en la tabla de resultados una columna que indique cuántos valores repetidos quedan después del filtrado.
2. Genes de referencia múltiples: Cree una versión que reciba varios genes reference_gene (por ejemplo, ["GAPDH", "ACTB"]) y utilice su media geométrica de Ct como referencia.
3. Gráfico de cambio de expresión: Utilice matplotlib para crear un gráfico de barras del cambio de expresión + barras de error (fc_std). Diferencie el grupo de control (cambio de expresión = 1) en gris y los grupos de tratamiento en color.
4. Desafío: Extensión de Pfaffl: Cree una versión que reciba la eficiencia E de cada gen como una columna CSV y calcule (E_target^-ΔCt_target) / (E_ref^-ΔCt_ref).

Resumen

Hemos resuelto el problema de "extraer el cambio de expresión y los intervalos de confianza de los resultados de qPCR" en tres partes:

• pandas groupby nos permitió obtener el promedio por condición y gen sin usar bucles.
• La aritmética vectorial calculó ΔCt, ΔΔCt y el cambio de expresión completo de una sola vez en el DataFrame.
• La propagación de errores transfirió con precisión la desviación de los datos originales al intervalo de confianza de los resultados.

Un cálculo que tardaba un día en Excel ahora se realiza con una sola llamada de función. Los errores que hacían que los datos de cambio de expresión de las presentaciones y los artículos fueran diferentes desaparecen. El riesgo de que un error en una sola celda destruya un artículo se elimina al integrarlo en la canalización.

Este artículo es un ejemplo educativo general. Las herramientas de análisis de qPCR prácticas (qbase+, PrimePCR Analysis, etc.) incluyen corrección de eficiencia, referencias múltiples, detección automática de valores atípicos y pruebas estadísticas GLM. Puede agregar estas funciones a esta estructura básica o confiar en herramientas validadas.