프로테오믹스 분석 플랫폼을 직접 만들어봤다

BioAI Market이란

BioAI Market은 웹 브라우저에서 프로테오믹스 데이터를 업로드하고, QC → DE → Pathway 분석을 원클릭으로 실행할 수 있는 플랫폼이다. R과 Python 코딩 없이도 분석이 가능하도록 만드는 게 목표였다.

기존에 프로테오믹스 분석을 하려면 R에서 limma 스크립트를 짜거나, Perseus 같은 데스크톱 소프트웨어를 쓰거나, Python으로 직접 파이프라인을 구축해야 했다. 연구자들이 이 과정에서 겪는 고통을 줄이고 싶었다.

분석 파이프라인 설계

전체 파이프라인은 4단계로 구성된다:

[전처리] → [QC] → [DE] → [Pathway]

전처리: Filtering → Log2 Transform → Median Normalization → (선택) Imputation
QC:     PCA, CV(Coefficient of Variation), Sample Correlation
DE:     limma / Welch's t-test / ANOVA + Tukey HSD
Pathway: GO (Biological Process, Molecular Function, Cellular Component) + KEGG

전처리 — 순서가 중요하다

프로테오믹스 raw 데이터는 결측치(missing values)가 많고, intensity 값의 분포가 skewed되어 있다. 전처리 순서를 잘못 잡으면 결과가 완전히 달라진다.

1. Filtering    — 그룹별 valid value 기준으로 결측이 많은 단백질 제거
2. Log2 Transform — 0을 NA로 바꾼 뒤 log2
3. Median Normalization — 샘플 간 중앙값을 맞춤
4. (선택) Imputation — log2 스케일에서, 결측 원인(MNAR/MAR)에 맞는 방법으로

처음에는 Log2 → Imputation 순서로 했다가, log2(0) = -Infinity 문제에 부딪혔다. 그때는 «imputation을 먼저 하면 된다»고 정리했는데, 다시 보니 그건 잘못된 해법이었다. 맞는 해법은 0을 NA로 바꾼 뒤 log2를 하고, 필요하면 그다음에 log 스케일에서 대체하는 것이다. raw intensity에서 «최솟값의 1/2»로 채우면 분산도 인위적으로 줄어든다. (2026-10 수정)

DE 방법 자동 선택 — 결국 limma

처음에는 사용자가 통계를 몰라도 되도록, 정규성(Shapiro-Wilk)과 등분산성(Levene) 검정 결과로 DE 방법을 자동 추천하는 함수를 만들었다. 지금 보면 이 설계는 잘못됐다. 시료 하나(첫 열)로 정규성을 보고, 모든 그룹을 합쳐서 등분산을 봤으며, 그룹당 3~4개 값으로는 이런 검정이 거의 아무것도 말해주지 못한다. 실제로 함수는 거의 항상 limma를 돌려줬다. (2026-10 수정: 코드를 내리고 설명으로 바꿨다)

결론은 단순하다. 2그룹이든 3그룹 이상이든 limma(3그룹 이상은 moderated F-test + contrast). Empirical Bayes moderated t-test는 소규모 샘플에서 분산 추정이 안정적이라, 프로테오믹스처럼 반복 수가 적은 실험에 맞다.

ANOVA 벡터화 — 다시 재 보니

3그룹 이상 비교에서 ANOVA를 돌려야 하는데, 처음에는 순진하게 per-protein 루프를 돌렸다:

# ❌ Before: per-protein 루프
from scipy.stats import f_oneway

results = []
for protein_idx in range(n_proteins):  # 5000+ proteins
    groups_data = []
    for group in unique_groups:         # 9 groups
        mask = group_labels == group
        groups_data.append(data[mask, protein_idx])

    f_stat, p_value = f_oneway(*groups_data)
    results.append({"protein": proteins[protein_idx], "f": f_stat, "p": p_value})

5000개 단백질 × 9그룹에서 이 루프가 병목이었다.

Numpy 행렬 연산으로 벡터화했다:

# ✅ After: numpy 벡터화
import numpy as np

def vectorized_anova(data, group_labels):
    """
    data: (n_samples, n_proteins) array
    group_labels: (n_samples,) array
    """
    unique_groups = np.unique(group_labels)
    n_groups = len(unique_groups)
    n_total = data.shape[0]
    grand_mean = data.mean(axis=0)  # (n_proteins,)

    # SSB (Sum of Squares Between)
    ssb = np.zeros(data.shape[1])
    # SSW (Sum of Squares Within)
    ssw = np.zeros(data.shape[1])

    for group in unique_groups:
        mask = group_labels == group
        n_g = mask.sum()
        group_mean = data[mask].mean(axis=0)
        ssb += n_g * (group_mean - grand_mean) ** 2
        ssw += ((data[mask] - group_mean) ** 2).sum(axis=0)

    df_between = n_groups - 1
    df_within = n_total - n_groups

    msb = ssb / df_between
    msw = ssw / df_within
    f_stats = msb / msw

    from scipy.stats import f as f_dist
    p_values = f_dist.sf(f_stats, df_between, df_within)   # 1 - cdf 는 아주 작은 p값에서 0으로 뭉개진다

    return f_stats, p_values

(2026-10 수정) 이 글의 첫 판에는 «215초 → 5초»라고 적었는데, 위 두 코드만 그대로 다시 재 보니 그 숫자가 나오지 않았다. 5,000 단백질 × 9그룹 × 3반복 합성 데이터, i9-14900K, 5회 중앙값:

Before (per-protein f_oneway 루프): 1.61초
After  (numpy 벡터화):              0.5밀리초  (F값은 두 방법이 소수점 14자리까지 같음)

벡터화가 극적으로 빠른 건 맞다. 다만 처음 적은 215초에는 이 코드 말고 다른 단계(데이터 로딩, 모든 단백질에 대한 사후검정 등)가 섞여 있었던 것으로 보인다. 그 측정이 정확히 무엇을 포함했는지는 확인하지 못해서, 숫자는 위 재측정으로 바꿨다.

Tukey HSD — significant만 처리

ANOVA에서 유의미한 단백질이 나오면 어느 그룹 간 차이인지 Tukey HSD로 다중 비교를 한다. 하지만 5000개 전부에 대해 Tukey를 돌리면 또 느려지니, BH(Benjamini-Hochberg) correction 후 significant한 단백질만 처리했다:

from statsmodels.stats.multicomp import pairwise_tukeyhsd
from statsmodels.stats.multitest import multipletests

# BH correction
_, adj_pvalues, _, _ = multipletests(p_values, method='fdr_bh')

# Significant 단백질만 Tukey HSD
significant_mask = adj_pvalues < 0.05
significant_indices = np.where(significant_mask)[0]

tukey_results = {}
for idx in significant_indices:
    result = pairwise_tukeyhsd(data[:, idx], group_labels, alpha=0.05)
    tukey_results[proteins[idx]] = result

멀티그룹 UI — 드래그앤드롭의 고통

기술적으로 가장 까다로웠던 건 UI였다. 사용자가 N개 그룹을 동적으로 추가/삭제하고, 샘플을 드래그앤드롭으로 그룹에 배정하는 인터페이스를 만들어야 했다.

처음에는 그룹 이름으로 매칭을 시도했다:

// ❌ 이름 매칭 — 실패
// "Control_1", "Control_2" → "Control" 그룹
// "Treatment_A_1", "Treatment_A_2" → "Treatment_A" 그룹
// 근데 "WT_young_rep1" 같은 이름은? 규칙이 없다

샘플 이름의 네이밍 규칙이 연구실마다 다 달라서 자동 매칭이 불가능했다. 결국 위치 기반 매핑으로 전환했다. 사용자가 CSV 컬럼 순서대로 그룹을 지정하는 방식:

// ✅ 위치 기반 매핑
// 컬럼 순서: [S1, S2, S3, S4, S5, S6]
// 그룹 매핑: [Control, Control, Control, Treatment, Treatment, Treatment]
// → 앞 3개가 Control, 뒤 3개가 Treatment

그리고 비교 설정도 위치 기반으로:

interface Comparison {
  control: number  // 그룹 인덱스
  treatment: number  // 그룹 인덱스
}

// 그룹 0 vs 그룹 1 → {control: 0, treatment: 1}

3종 생물 지원 — Human, Mouse, Rat

Pathway 분석에서 organism을 지원해야 한다. GO/KEGG 분석에 사용하는 annotation DB가 종마다 다르기 때문이다:

ORGANISM_CONFIG = {
    "human": {
        "orgDb": "org.Hs.eg.db",
        "kegg": "hsa",
        "name": "Homo sapiens"
    },
    "mouse": {
        "orgDb": "org.Mm.eg.db",
        "kegg": "mmu",
        "name": "Mus musculus"
    },
    "rat": {
        "orgDb": "org.Rn.eg.db",
        "kegg": "rno",
        "name": "Rattus norvegicus"
    }
}

처음에는 Human만 지원했는데, "Mouse 데이터도 분석할 수 있나요?"라는 피드백을 받고 추가했다. annotation DB가 각각 수백 MB라서 Docker 이미지가 4GB를 넘긴 이유이기도 하다.

현재 상태와 남은 과제

BioAI Market은 현재 아래 기능이 동작한다:

• ✅ CSV 업로드 및 전처리
• ✅ QC 분석 (PCA, CV, Correlation)
• ✅ DE 분석 (limma, Welch, ANOVA + Tukey)
• ✅ Pathway 분석 (GO + KEGG)
• ✅ AI 결과 해석 (Ollama)
• ✅ 바이오마커 시맨틱 검색
• ✅ 리포트 자동 생성

남은 과제:

• TMT/iTRAQ 등 다양한 정량 방식 지원
• GSEA(Gene Set Enrichment Analysis) 추가
• 배치 분석 (여러 데이터셋 동시 처리)
• 논문 figure 수준의 시각화

직접 만들어보면서 느낀 건, 프로테오믹스 분석이 "코드 몇 줄"이 아니라 수많은 통계적 판단과 예외 처리의 집합체라는 것이다. 그래도 만드는 과정 자체가 공부가 많이 됐다.

---

참고 링크:

• limma Bioconductor 패키지
• clusterProfiler — Yu Lab
• Perseus 소프트웨어