qPCR 분석 자동화 — ΔΔCt를 손 안 대고 fold change까지

이 토픽을 마치면

교과서에서 배운 pandas groupby, 오차 전파, CSV 입출력 을 엮어서, qPCR 실험 결과 CSV를 받아 자동으로 ΔΔCt·fold change·신뢰 구간까지 계산하는 파이프라인을 직접 만들 수 있습니다. Excel로 매번 셀을 손으로 잡아 계산하던 그 작업이 함수 호출 한 번으로 대체됩니다.

이 글은 교육용 일반 예제 입니다. qPCR 분석은 분자생물학 실험실에서 매일 하는 작업이라 소재로 삼았습니다.

"그 셀 하나만 잘못 잡아도" — Excel의 위험

qPCR 실험을 마치면 이런 CSV 파일이 나옵니다.

text

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
...

여러분이 하고 싶은 것은 이렇습니다.

각 샘플의 ΔCt = 관심 유전자 Ct − 기준 유전자(GAPDH) Ct
각 조건의 평균 ΔCt 를 replicate 간에 계산
ΔΔCt = 처리군 평균 ΔCt − 대조군 평균 ΔCt
Fold change = 2^(-ΔΔCt)
오차 전파 로 fold change의 표준편차

Excel로 이걸 하는 방법은 익숙합니다. VLOOKUP, AVERAGEIF, 셀을 손으로 잡고 수식 드래그, 결과를 새 시트에 옮기고, 또 계산하고...

그런데 자주 이런 일이 벌어집니다.

셀 하나만 잘못 잡아서 대조군 데이터가 처리군 계산에 들어감
첫 실험은 잘 됐는데 두 번째 실험에서 컬럼 순서가 달라져서 수식이 안 맞음
3반복 중 하나를 이상치라고 뺐는데, 어느 걸 뺐는지 기록이 없어 재현 안 됨
발표 자료의 fold change와 논문 그림의 fold change가 미묘하게 다른데 어디서 갈렸는지 찾을 수 없음

이런 오류가 실제로 논문 철회 사유가 된 사례가 여럿 있습니다. 셀 하나가 논문 하나를 부순다.

이 글에서는 그 실수의 여지를 없애는 파이프라인을 만듭니다. CSV를 넣으면 결과가 나온다. 중간에 사람 손이 안 들어간다.

먼저 완성품을 봅시다 (블랙박스 먼저 돌리기)

우리가 만들 도구는 이렇게 씁니다.

python

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

text

=== 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

파일 저장: qpcr_result.csv (재현 가능)

같은 실험 데이터에 대해 매번 동일한 결과 가 나옵니다. 파일 경로만 넣으면 나머지는 코드가 알아서 합니다. 실수가 개입할 자리가 없습니다.

이 도구는 어떤 부품으로 조립되는가 (부품 분해도)

text

qPCR 자동 분석 파이프라인
   ┌──────────────────────────────────────────────────┐
   │  [입력] CSV 로딩 ────────────── 부품: CSV 입출력    │  ← 완성 제공 (도구)
   │              │                                      │
   │              ▼                                      │
   │  [1단계] 조건·유전자별 그룹화 + 평균                   │
   │        부품: pandas groupby                          │  ← 직접 만든다 ★
   │              │                                      │
   │              ▼                                      │
   │  [2단계] ΔCt · ΔΔCt · fold change 벡터 계산          │
   │        부품: 벡터화 산술                              │  ← 직접 만든다 ★
   │              │                                      │
   │              ▼                                      │
   │  [3단계] 오차 전파                                    │
   │        부품: 표준편차 조합 규칙                        │  ← 직접 만든다 ★
   │              │                                      │
   │              ▼                                      │
   │  [출력] 결과 CSV + 그래프                             │  ← 완성 제공 (도구)
   └──────────────────────────────────────────────────┘

부품	어디서 배웠나	이 도구에서 하는 일
CSV 입출력	`csv-io`	실험 데이터 읽고 결과 저장
pandas groupby	`pandas-groupby`	조건·유전자별로 평균과 표준편차
오차 전파	`error-propagation`	ΔCt의 표준편차를 fold change의 표준편차로
matplotlib	`matplotlib-basics`	fold change와 오차 막대 그래프

📌 이 개념들 처음 본다면 (상단 진입 링크)

pandas groupby — split-apply-combine

오차 전파 — 표준편차의 조합 규칙

CSV — 표 데이터의 표준 형식

직접 만드는 새 개념은 groupby · 벡터 계산 · 오차 전파 3개. CSV I/O와 그래프는 도구라 완성본으로 드립니다. 딱 3개죠 — 인지 한계선 안입니다.

만들기 1단계 — 데이터 준비 (완성 제공)

실제 qPCR 결과를 흉내낸 데이터를 만들어봅시다. 실무에서는 이 데이터가 장비에서 CSV로 나옵니다.

python

import pandas as pd
import numpy as np

# 재현 가능한 난수 (실험 노이즈 흉내)
rng = np.random.default_rng(42)

def make_qpcr_data():
    rows = []
    # 진짜 시나리오: control 대비 treated에서 BRCA1이 약 8배, TP53이 약 3배 상승
    scenarios = {
        ("control", "BRCA1"): 25.0,
        ("control", "TP53"):  23.5,
        ("control", "GAPDH"): 18.3,
        ("treated", "BRCA1"): 22.0,   # Ct 3 감소 = 약 8배 상승
        ("treated", "TP53"):  21.8,   # Ct 1.7 감소 = 약 3배 상승
        ("treated", "GAPDH"): 18.4,   # 기준 유전자는 거의 안 변함
    }
    for (condition, gene), base_ct in scenarios.items():
        for rep in range(1, 4):  # 3 replicates
            ct = base_ct + rng.normal(0, 0.15)  # 기술 편차
            rows.append({
                "sample": f"{condition}_rep{rep}",
                "condition": condition,
                "gene": gene,
                "ct": round(ct, 3),
            })
    return pd.DataFrame(rows)

df = make_qpcr_data()

# 검증: 데이터 형태
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))

이 표가 qPCR 실험의 원 데이터 형태입니다. 이제 이것을 처리 대상으로 삼습니다.

만들기 2단계 — 조건·유전자별로 평균 내기 ★ (pandas groupby)

✍️ 직접 채우는 구간. 부품 = pandas groupby. 목표: 각 (조건, 유전자) 조합의 평균 Ct와 표준편차를 replicate 간에 계산.

3 replicates가 있으니 이 3개의 평균과 표준편차를 뽑아야 합니다. 이걸 각 (조건, 유전자) 조합마다 해야 합니다. 순진하게 짜면 이런 이중 for 문이 나옵니다.

python

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

이게 답을 주긴 합니다. 그런데 조건이 10개, 유전자가 100개면 이중 for는 1000번 필터링을 반복합니다. 데이터가 커지면 느려집니다.

pandas groupby 는 이 패턴을 한 줄로 대체합니다.

🔎 groupby란 (드로어 — pandas-groupby) "같은 그룹의 행을 모아, 각 그룹에 함수를 적용하고, 결과를 합친다"(split-apply-combine). SQL의 GROUP BY와 정확히 같은 개념입니다. for 문 없이 벡터로 처리되어 훨씬 빠릅니다.

python

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)

# 검증: 순진 버전과 결과 일치
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()   # 3 replicates 확인

같은 결과, 한 줄. 이것이 groupby의 힘입니다.

🤔 자기설명 프롬프트 .agg(["mean", "std", "count"])에서 count를 왜 굳이 계산할까요? 실무에서 replicate 개수가 표에서 늘 3이라 보장할 수 있나요? (힌트: 이상치를 뺐거나, 로딩 실패 웰이 있으면 그룹마다 n이 달라집니다.)

만들기 3단계 — ΔCt · ΔΔCt · fold change 벡터로 계산 ★

✍️ 직접 채우는 구간. 부품 = 벡터화 산술. 목표: 요약 표에 ΔCt/ΔΔCt/fold change 컬럼을 for 문 없이 채워 넣기.

지금 우리의 요약 표는 이렇게 생겼습니다.

text

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

여기서 ΔCt = mean_ct - mean_ct(GAPDH, 같은 조건) 을 계산해야 합니다. 각 행마다 그 조건의 GAPDH를 찾아 빼야 하니, 얼핏 복잡해 보입니다. 그런데 pandas에서는 merge 로 한 방에 됩니다.

python

def add_delta_ct(summary, reference_gene="GAPDH"):
    # 기준 유전자만 뽑기
    ref = (
        summary[summary["gene"] == reference_gene]
        [["condition", "mean_ct", "std_ct"]]
        .rename(columns={"mean_ct": "mean_ct_ref", "std_ct": "std_ct_ref"})
    )
    # 조건별로 붙이기
    merged = summary.merge(ref, on="condition", how="left")
    merged["dct"] = merged["mean_ct"] - merged["mean_ct_ref"]
    # 표준편차 전파: 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"]])

# 검증: 기준 유전자의 ΔCt는 자기 자신에서 뺀 것 = 0
assert np.allclose(
    with_dct[with_dct["gene"] == "GAPDH"]["dct"].values, 0.0
)
# ΔCt는 관심 유전자에서 GAPDH를 뺀 값
control_brca1_dct = with_dct[
    (with_dct["condition"] == "control") & (with_dct["gene"] == "BRCA1")
]["dct"].iloc[0]
# 25.0 - 18.3 ≈ 6.7 (실험 노이즈 감안)
assert 6.0 < control_brca1_dct < 7.5

이제 ΔΔCt 는 각 유전자에 대해 처리군 ΔCt − 대조군 ΔCt 입니다. 같은 merge 패턴을 씁니다.

python

def add_ddct_and_fold(with_dct, control_condition="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"]
    # 오차 전파: 두 Δ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"])
    # 오차 전파: fc의 sd = fc * ln(2) * ddct_std (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"]])

# 검증: 대조군의 ΔΔCt는 0, fold change는 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)
# treated의 BRCA1이 8배 근처여야 함 (설계상)
treated_brca1_fc = final[
    (final["condition"] == "treated") & (final["gene"] == "BRCA1")
]["fold_change"].iloc[0]
assert 5 < treated_brca1_fc < 12

이 코드에는 for 문이 하나도 없습니다. 계산이 데이터프레임 전체에 벡터로 적용됩니다.

🤔 자기설명 프롬프트 Fold change의 표준편차 공식이 fc * ln(2) * ddct_std인 이유를 짚어보세요. 함수 f(x) = 2^(-x)의 미분과 오차 전파의 일반 공식(σ_f ≈ |f'| · σ_x)에서 유도됩니다. (힌트: d/dx[2^(-x)] = -ln(2) · 2^(-x))

만들기 4단계 — 오차 전파 원리 짚기 ★ (error-propagation)

우리는 이미 오차 전파를 두 번 썼습니다. 명시적으로 짚어봅시다.

🔎 오차 전파란 (드로어 — error-propagation) 두 측정값 A와 B가 각각 표준편차 σA, σB를 가진다고 합시다. 이 둘의 합·차 는 σ = √(σA² + σB²)로 전파됩니다(독립 가정). 곱·나눗셈 은 상대 표준편차가 결합됩니다. 함수 f(A) 는 σ ≈ |f'(A)| · σA. 이 세 규칙이 우리가 실험 데이터를 다루는 핵심입니다.

이 원칙을 우리 파이프라인에 적용:

ΔCt = Ct(target) - Ct(GAPDH) → 차 → σ_ΔCt = √(σ_target² + σ_GAPDH²)
ΔΔCt = ΔCt(treated) - ΔCt(control) → 차 → σ_ΔΔCt = √(σ_treated² + σ_control²)
Fold change = 2^(-ΔΔCt) → 비선형 함수 → σ_fc ≈ |d/dx[2^(-x)]| · σ_ΔΔCt = fc · ln(2) · σ_ΔΔCt

python

# 손 계산으로 확인
# 가상: σ_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
# 다음, σ_ΔΔ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
# 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

이 손 계산이 우리 파이프라인 결과와 같은 규모여야 합니다. 오차 전파는 손으로도 검증 가능하다 는 사실이 이 원칙의 힘입니다.

부품을 하나로 — 완성된 파이프라인 함수

세 부품을 하나의 함수로 묶습니다.

python

def analyze_qpcr(
    df: pd.DataFrame,
    reference_gene: str = "GAPDH",
    control_condition: str = "control",
) -> pd.DataFrame:
    """qPCR 원 데이터 → ΔΔCt · fold change · 오차 전파 표."""
    # 1. 조건·유전자별 요약
    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. 표시할 컬럼만 선택
    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)

# 대조군: fold change = 1, 편차 = 0에 가까움
ctrl = result[result["condition"] == "control"]
assert np.allclose(ctrl["fold_change"].values, 1.0)
# 처리군 BRCA1: 설계상 약 8배
treated_brca1_fc = result[
    (result["condition"] == "treated") & (result["gene"] == "BRCA1")
]["fold_change"].iloc[0]
assert 5 < treated_brca1_fc < 12

이 도구가 오늘 우리가 이야기하는 Livak 방법(ΔΔCt method, 2001) 의 구현체입니다. 실무 도구(qbase+, PrimePCR Analysis 등)는 여기에 효율 보정, 다중 참조 유전자, 이상치 자동 탐지 등을 추가한 것뿐입니다. 뼈대는 여러분이 방금 만든 것과 같습니다.

성능 딥다이브 — 왜 이렇게 빠르나

python

import time

# 큰 데이터 시뮬레이션: 100개 유전자, 20개 조건, 4 replicates
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 · 조건 20 · rep 4 → {elapsed*1000:.1f}ms")
assert elapsed < 2.0  # 몇 초 안에 끝나야

행 수가 수천 개여도 밀리초 단위입니다. for 문 없이 groupby + 벡터 계산 이기 때문입니다.

Excel로 같은 크기를 해보면? 유전자별 시트 하나씩 만들고, 조건마다 참조 셀 다시 잡고... 몇 시간에서 하루가 걸립니다. 그리고 도중에 한 셀만 틀려도 결과가 조용히 오염됩니다.

다른 길도 있다 (멀티패스 성찰)

효율 보정(Pfaffl method): Livak 방법은 PCR 효율이 100%(매 사이클마다 정확히 2배)라고 가정합니다. 실제로는 프라이머·조건마다 효율이 조금씩 다릅니다. Pfaffl 방법은 각 유전자의 실제 효율(E, 1.8~2.0 사이)을 반영합니다. 언제 무엇: 스크리닝·pilot = Livak / 결정적 결과·논문 그림 = Pfaffl.
다중 참조 유전자: GAPDH 하나에 의존하는 것은 위험합니다. GAPDH 자체가 조건에 따라 변할 수 있어서요. geometric mean of multiple reference genes (GeNorm, NormFinder)가 더 안전합니다.
웰별 이상치 자동 탐지: 3 replicates 중 하나가 다른 둘과 크게 다르면(예: >0.5 Ct 차이) 자동으로 이상치 표시. 우리 파이프라인에 추가하면 재현성이 크게 올라갑니다.
DataFrame vs SQL: 우리는 pandas로 짰지만, 실험 데이터가 매우 크면(연구소 전체 실험 아카이브 등) DuckDB나 PostgreSQL로 옮겨서 SQL로 처리 하는 게 확장성이 낫습니다. 알고리즘은 그대로.

핵심: "for 문 대신 groupby, 손 계산 대신 벡터 산술, 흐름 대신 파일." 이 세 원칙이 실험 분석의 재현성을 결정합니다. 여러분이 방금 만든 것이 그 원칙의 실체입니다.

다음 단계로 (하단 출구 링크)

왜 벡터 산술이 for보다 압도적으로 빠른지 → 벡터화
오차 전파의 상세 규칙(곱·나눗셈·비선형) → 오차 전파 원리
결과 시각화로 이어가기 → 응용편 RNA-seq 히트맵 대시보드

직접 해보기 (독립 문제)

이상치 필터: 3 replicate 중 다른 두 개와 0.5 Ct 이상 차이나는 replicate를 자동으로 제거하는 전처리를 추가하세요. 필터 후 몇 개 replicate가 남았는지도 결과 표에 컬럼으로 넣으세요.
다중 참조 유전자: reference_gene 을 여러 개(예: ["GAPDH", "ACTB"])로 받아 이들의 기하 평균 Ct 를 참조로 삼는 버전을 만드세요.
fold change 그래프: matplotlib으로 fold change 막대 + 오차 막대(fc_std)를 그리세요. 대조군(fold=1)은 회색, 처리군은 색으로 구분.
도전 — Pfaffl 확장: 각 유전자의 효율 E를 CSV 컬럼으로 받고 (E_target^-ΔCt_target) / (E_ref^-ΔCt_ref) 로 계산하는 버전을 짜세요.

정리

우리는 "qPCR 결과에서 fold change와 신뢰 구간을 뽑는" 문제를, 세 개의 부품으로 정복했습니다.

pandas groupby 가 조건·유전자별 평균을 for 없이 얻게 해줬습니다.
벡터화 산술 이 ΔCt · ΔΔCt · fold change 전체를 데이터프레임 위에서 한 번에 계산했습니다.
오차 전파 가 원 데이터의 편차를 결과의 신뢰 구간으로 정확히 옮겼습니다.

Excel로 하루 걸리던 계산이 함수 호출 한 번이 됐습니다. 발표 자료와 논문의 fold change가 갈리는 사고도 사라집니다. 셀 한 개의 실수가 논문 하나를 부수던 위험 이 파이프라인 안으로 밀봉됩니다.

이 글은 일반 교육 예제 입니다. 실무 qPCR 분석 도구(qbase+, PrimePCR Analysis 등)는 여기에 효율 보정, 다중 참조, 자동 이상치 탐지, GLM 통계 검정 등이 추가된 것입니다. 그 상세한 버전은 이 뼈대 위에 여러분이 붙이거나, 검증된 도구에 맡기면 됩니다.