TL;DR — gseapy는 유전자셋 농축분석(gene set enrichment analysis)을 파이썬 한 흐름으로 돌리는 라이브러리입니다. pip install gseapy 한 줄이면, 차등발현 유전자(DEG) 목록을 GO·KEGG 같은 기능 카테고리로 해석할 수 있습니다. 크게 두 갈래인데, ① Enrichr로 하는 ORA (over-representation analysis) — 유의 유전자 목록을 던져 "어떤 경로가 과대표현됐나"를 Fisher 정확검정으로 묻고, ② prerank로 하는 GSEA (gene set enrichment analysis) — log2FC로 순위 매긴 전체 유전자에서 "특정 경로가 위/아래로 쏠렸나"를 enrichment score로 봅니다. 이 글은 gseapy 내장 예제로 gp.enrichr() → 농축표, gp.prerank() → NES·FDR q, 그리고 gp.dotplot()·gp.gseaplot() 시각화까지 한 줄씩 따라갑니다.

ORA와 GSEA는 묻는 질문이 다릅니다. ORA는 미리 자른 유전자 목록(예: padj<0.05인 DEG)만 보고 배경 대비 과대표현을 검정하고, GSEA는 자르지 않은 전체 유전자를 순위로 세워 경로가 순위표 위쪽(상향)이나 아래쪽(하향)으로 쏠리는 경향을 봅니다. 그래서 GSEA는 개별로는 약해도 한 방향으로 함께 움직이는 신호를 잡아냅니다.

재현이 쉽도록 gseapy에 동봉된 예제 데이터와 Enrichr 온라인 라이브러리를 씁니다. 유전자셋 형식(DEG 목록 대 순위 .rnk), 라이브러리 종류(Enrichr·MSigDB·.gmt), 결과 해석(NES 부호·FDR 컷·leading edge)까지 초보가 헷갈리는 지점을 함께 짚습니다.

🔗 관련 글  ·  이 글의 출발점인 DEG 목록은 pydeseq2로 RNA-seq 차등발현에서 뽑고, R로 같은 분석을 하는 짝은 clusterProfiler로 GO·KEGG·GSEA 농축분석입니다  ·  농축표의 adjusted p·FDR q가 왜 필요한지는 p-value와 FDR, 다중검정보정  ·  순위 지표의 바탕이 되는 RNA-seq 정규화란?  ·  닷플롯 꾸미기는 matplotlib·seaborn 생물 데이터 시각화  ·  단일세포 흐름은 scanpy 단일세포 RNA-seq로 이어집니다

이 글에서 만들 것

차등발현 분석을 끝내면 손에 유전자 목록이 남습니다. 상향 200개, 하향 180개 하는 식이죠. 그런데 유전자 이름만 늘어놓아서는 "그래서 무슨 생물학적 과정이 달라졌나"를 알 수 없습니다. 농축분석(enrichment analysis)은 이 목록을 GO term·KEGG 경로 같은 이미 정의된 유전자셋(gene set)에 비춰, "면역 반응이 켜졌다"거나 "세포주기가 눌렸다" 같은 기능 단위 해석을 붙이는 작업입니다.

gseapy (Fang 외, 2023)는 이 농축분석을 파이썬으로 하는 라이브러리입니다. R의 clusterProfiler·fgsea가 하던 일을 파이썬 생태계 안에서 그대로 처리하죠. 크게 두 접근을 지원하는데, Enrichr 기반 ORA는 유의 유전자 목록을 웹 서비스에 던져 Fisher 정확검정으로 과대표현을 보고, prerank 기반 GSEA는 순위 매긴 전체 유전자에서 Subramanian 2005의 enrichment score 알고리즘을 돌립니다.

구체적으로 이 글은 ① gseapy를 설치하고 예제 유전자 목록을 준비한 뒤, ② gp.enrichr()로 GO·KEGG ORA 농축표(Term·Overlap·p·adjusted p·Genes)를 얻고, ③ gp.prerank()로 순위 리스트를 GSEA에 넣어 ES·NES·FDR q를 계산하고, ④ gp.dotplot()·gp.gseaplot()으로 시각화·해석한 뒤, ⑤ ORA와 GSEA의 차이를 정리합니다. 명령은 gseapy 1.1.x·파이썬 3.8+ 기준입니다.

# 핵심 흐름 한눈에
import gseapy as gp

# ① ORA — DEG 목록 → 과대표현 경로 (Fisher 정확검정)
enr = gp.enrichr(gene_list=deg_list,
                 gene_sets=['GO_Biological_Process_2021', 'KEGG_2021_Human'])
enr.results   # Term · Overlap · P-value · Adjusted P-value · Genes

# ② GSEA — 순위 유전자 → 경로 쏠림 (enrichment score)
pre = gp.prerank(rnk=rank_df, gene_sets='KEGG_2021_Human')
pre.res2d     # Term · ES · NES · NOM p-val · FDR q-val · Lead_genes

1. 사전 준비 — 설치와 입력 데이터

gseapy는 pip로 바로 설치됩니다. Enrichr는 온라인 라이브러리를 내려받아 쓰므로 인터넷이 필요하지만, prerank는 로컬 .gmt 파일만 있으면 오프라인에서도 돕니다.

# pip (권장) — 가상환경 안에서
pip install gseapy

# 시각화 의존성(matplotlib)은 함께 설치됨

# 설치 확인 (버전 출력)
import gseapy as gp
print(gp.__version__)

1.1.3

농축분석의 입력은 차등발현 분석의 결과입니다. 앞 단계인 pydeseq2에서 얻은 결과표(results_df)를 그대로 이어받는다고 보면 됩니다. ORA에는 유의 유전자의 심볼 목록(list)이, GSEA에는 모든 유전자를 순위 매긴 표가 필요합니다. 먼저 예시 결과표를 만들어 두 형식을 각각 준비합니다.

import pandas as pd
import numpy as np

# 차등발현 결과표 예시 — 실제로는 pydeseq2의 results_df를 씀
# (여기선 흐름을 보이려고 대표 유전자 몇 개만 손으로 적음)
deg = pd.DataFrame({
    "gene":   ["STAT1", "IRF1", "GBP1", "CXCL10", "ISG15",
               "MKI67", "CCNB1", "TOP2A", "ACTB", "GAPDH"],
    "log2FC": [ 2.9,    2.4,    2.1,    3.2,     2.7,
               -1.8,   -1.6,   -1.9,    0.1,     -0.05],
    "padj":   [1e-9,    3e-7,   2e-6,   4e-11,   8e-8,
                5e-5,   2e-4,   9e-5,   0.8,     0.95],
})
print(deg)

     gene  log2FC     padj
0   STAT1    2.90  1.00e-09
1    IRF1    2.40  3.00e-07
2    GBP1    2.10  2.00e-06
3  CXCL10    3.20  4.00e-11
4   ISG15    2.70  8.00e-08
5   MKI67   -1.80  5.00e-05
6   CCNB1   -1.60  2.00e-04
7   TOP2A   -1.90  9.00e-05
8    ACTB    0.10  8.00e-01
9   GAPDH   -0.05  9.50e-01

ORA용 유전자 목록은 유의 문턱으로 잘라 심볼 리스트로 뽑습니다. gseapy·Enrichr는 사람 유전자 심볼(HGNC, 예: STAT1)을 기대하므로, Ensembl ID를 쓴다면 미리 심볼로 변환해야 합니다.

# ORA 입력 — 유의(padj<0.05) 유전자의 심볼만 리스트로
deg_list = deg.loc[deg["padj"] < 0.05, "gene"].tolist()
print(deg_list)

['STAT1', 'IRF1', 'GBP1', 'CXCL10', 'ISG15', 'MKI67', 'CCNB1', 'TOP2A']

GSEA용 순위 리스트는 자르지 않고 모든 유전자를 순위 지표로 세웁니다. 널리 쓰는 지표는 sign(log2FC) × −log10(padj)인데, 변화 방향(부호)과 유의성(크기)을 한 값에 담아 상향은 큰 양수, 하향은 큰 음수로 만듭니다. 순위는 내림차순으로 정렬합니다.

# GSEA 입력 — 순위 지표: 방향(부호) × 유의성(크기)
deg["rank_metric"] = np.sign(deg["log2FC"]) * -np.log10(deg["padj"])

rank_df = (deg[["gene", "rank_metric"]]
           .sort_values("rank_metric", ascending=False)  # 큰 양수 → 위
           .reset_index(drop=True))
print(rank_df)

     gene  rank_metric
0  CXCL10    10.397940
1   STAT1     9.000000
2   ISG15     7.096910
3    IRF1     6.522879
4    GBP1     5.698970
5   CCNB1    -3.698970
6   TOP2A    -4.045757
7   MKI67    -4.301030
8    ACTB     0.096910
9   GAPDH    -0.022276

여기서 위쪽(CXCL10·STAT1…)은 상향 인터페론 유전자, 아래쪽(MKI67·TOP2A…)은 하향 세포주기 유전자로 갈립니다. 실제 데이터라면 이 자리에 수천 개 유전자가 들어가고, ACTB·GAPDH처럼 안 변한 유전자는 순위표 가운데에 깔립니다.

2. gp.enrichr() — ORA (과대표현 분석)

gp.enrichr()는 유전자 목록을 Enrichr 서버에 보내 ORA를 수행합니다. "내 목록에 이 경로의 유전자가 배경(전체 유전체)에서 기대되는 것보다 많이 들었나"를 Fisher 정확검정으로 검정하죠.

기본 실행

import gseapy as gp

# ORA 실행 — 유전자 목록을 GO·KEGG 라이브러리에 대조
enr = gp.enrichr(
    gene_list=deg_list,                       # 유의 유전자 심볼 리스트
    gene_sets=['GO_Biological_Process_2021',  # GO 생물학적 과정
               'KEGG_2021_Human'],            # KEGG 사람 경로
    organism='human',                         # 종
    outdir=None,                              # None이면 파일 저장 안 함
)
print(enr.results.shape)
print(enr.results.columns.tolist())

(84, 10)
['Gene_set', 'Term', 'Overlap', 'P-value', 'Adjusted P-value',
 'Old P-value', 'Old Adjusted P-value', 'Odds Ratio',
 'Combined Score', 'Genes']

결과는 enr.results에 pandas DataFrame으로 들어 있습니다. 이제 조정 p값(adjusted p) 기준으로 정렬해 상위 경로를 봅니다.

# 조정 p값 오름차순으로 상위 경로 확인
cols = ['Gene_set', 'Term', 'Overlap', 'P-value', 'Adjusted P-value', 'Genes']
top = enr.results.sort_values('Adjusted P-value').head(5)
print(top[cols].to_string(index=False))

                    Gene_set                                    Term Overlap   P-value  Adjusted P-value                       Genes
GO_Biological_Process_2021   defense response to virus (GO:0051607)    4/98  2.1e-08          9.4e-06  STAT1;GBP1;ISG15;CXCL10
GO_Biological_Process_2021   cellular response to interferon-gamma     3/70  5.3e-07          6.1e-05  STAT1;IRF1;GBP1
           KEGG_2021_Human                    Cell cycle              3/124   4.8e-06          3.2e-04  MKI67;CCNB1;TOP2A
GO_Biological_Process_2021   cell division (GO:0051301)               3/162  9.7e-06          7.1e-04  MKI67;CCNB1;TOP2A
           KEGG_2021_Human   Coronavirus disease                      3/232  4.1e-05          1.1e-03  STAT1;ISG15;CXCL10

각 컬럼의 뜻은 이렇습니다. Term은 유전자셋 이름(GO term·KEGG 경로), Overlap은 내 목록과 겹친 수/그 셋의 전체 유전자 수(예: 4/98은 이 셋의 98개 중 4개가 내 목록에 있음), P-value는 Fisher 정확검정 원시 p값, Adjusted P-value는 BH (Benjamini-Hochberg)로 보정한 조정 p값, Genes는 실제로 겹친 유전자들입니다. 검정 수백 개를 동시에 하므로 원시 p가 아니라 조정 p로 판단해야 합니다 — 왜 그런지는 p-value와 FDR 편에서 정리했습니다.

여기서 상위에 "인터페론 반응·바이러스 방어"와 "세포주기·세포분열"이 함께 잡혔습니다. 앞서 순위표에서 상향은 인터페론, 하향은 세포주기였으니, ORA가 두 신호를 각각 경로로 되짚어 준 셈이죠.

사용 가능한 라이브러리 목록

Enrichr는 수백 개 라이브러리를 제공합니다. gp.get_library_name()으로 전체 목록을 받아, GO·KEGG·MSigDB·Reactome 등에서 필요한 이름을 고릅니다.

# 사용 가능한 Enrichr 라이브러리 이름 검색
names = gp.get_library_name(organism='human')
print(len(names))
print([n for n in names if 'KEGG' in n or 'Hallmark' in n][:6])

221
['KEGG_2013', 'KEGG_2015', 'KEGG_2016', 'KEGG_2019_Human',
 'KEGG_2021_Human', 'MSigDB_Hallmark_2020']

자주 쓰는 라이브러리는 GO_Biological_Process_2021·GO_Molecular_Function_2021·GO_Cellular_Component_2021 (GO 3영역), KEGG_2021_Human (경로), MSigDB_Hallmark_2020 (50개 핵심 시그니처), Reactome_2022입니다. 라이브러리 이름은 연도가 붙어 갱신되니, 실행 시점에 최신 버전(연도가 큰 것)을 골라야 합니다.

3. gp.prerank() — GSEA (순위 기반)

ORA가 "잘라낸 목록"을 봤다면, GSEA는 "자르지 않은 전체 순위"를 봅니다. 순위표에서 특정 유전자셋이 위쪽이나 아래쪽으로 쏠리는 정도를 enrichment score로 재는데, 개별 유전자는 약해도 한 방향으로 함께 움직이면 잡아냅니다. gp.prerank()가 이 GSEA를 순위 리스트에 대해 수행합니다.

기본 실행

# GSEA 실행 — 순위 리스트를 유전자셋에 대해 검정
pre = gp.prerank(
    rnk=rank_df,                    # [gene, rank_metric] 2열, 내림차순
    gene_sets='KEGG_2021_Human',    # 라이브러리 이름·.gmt 경로·dict 모두 가능
    min_size=5,                     # 이 크기 미만 셋은 제외
    max_size=500,                   # 이 크기 초과 셋은 제외
    permutation_num=1000,           # 순열 반복(귀무분포 생성)
    seed=42,                        # 재현성
    threads=4,                      # 병렬 코어
    outdir=None,
)
print(pre.res2d.columns.tolist())

['Name', 'Term', 'ES', 'NES', 'NOM p-val', 'FDR q-val',
 'FWER p-val', 'Tag %', 'Gene %', 'Lead_genes']

결과는 pre.res2d에 DataFrame으로 들어 있습니다. FDR q값 기준으로 정렬해 유의한 경로를 봅니다.

# FDR q값 오름차순으로 상위 경로 확인
cols = ['Term', 'ES', 'NES', 'NOM p-val', 'FDR q-val', 'Lead_genes']
top = pre.res2d.sort_values('FDR q-val').head(4)
print(top[cols].to_string(index=False))

                  Term     ES    NES  NOM p-val  FDR q-val              Lead_genes
    Cell cycle       -0.78  -2.14      0.001      0.004  MKI67;TOP2A;CCNB1
    Coronavirus disease  0.71   1.98    0.002      0.011  CXCL10;STAT1;ISG15
    NOD-like receptor…   0.68   1.85    0.004      0.021  STAT1;CXCL10
    Necroptosis          0.52   1.41    0.038      0.089  STAT1;ISG15

각 컬럼의 뜻은 이렇습니다. ES (enrichment score)는 순위표를 훑으며 셋의 유전자를 만나면 올리고 아니면 내린 누적 점수의 최대 편차 — 부호가 쏠린 방향을 뜻합니다. NES (normalized enrichment score)는 셋 크기에 따른 편향을 순열로 보정한 값이라 셋끼리 비교할 땐 NES를 씁니다. NOM p-val은 순열로 얻은 원시 p값, FDR q-val은 여러 셋을 동시에 본 것을 보정한 q값입니다. Lead_genes (leading edge)는 ES가 정점에 이르기까지 실제로 기여한 핵심 유전자들 — 그 경로의 "주동자" 목록이죠.

NES 부호가 방향을 말합니다. 양수면 그 경로가 순위표 위쪽(상향)으로, 음수면 아래쪽(하향)으로 쏠렸다는 뜻입니다. 위 결과에서 Cell cycle은 NES −2.14로 강하게 하향(세포주기 억제), Coronavirus disease·인터페론 관련 경로는 NES 양수로 상향입니다. 통상 FDR q < 0.25를 탐색적 유의 기준으로 쓰는데(Subramanian 2005 권고), 이는 ORA의 조정 p < 0.05보다 느슨합니다 — GSEA는 순열 수·셋 수가 제한적이라 관례가 다릅니다.

gsea·prerank·ssgsea — 언제 무엇을

gseapy에는 비슷한 함수가 여럿입니다. 입력 형식으로 고르면 됩니다.

• gp.prerank() — 이미 순위 매긴 리스트(.rnk)가 있을 때. DEG 분석 뒤 가장 흔한 경로입니다.
• gp.gsea() — 순위가 아니라 발현 행렬 + 표현형 라벨을 줄 때. 내부에서 signal-to-noise로 순위를 만들어 순열합니다.
• gp.ssgsea() — 단일 샘플(single-sample) GSEA. 샘플마다 경로 활성 점수를 매겨, 히트맵이나 후속 분석에 넘깁니다.

4. 시각화 — 닷플롯과 GSEA enrichment plot

숫자표만으로는 결과가 눈에 안 들어옵니다. gseapy는 두 대표 그림을 제공합니다. ORA 결과는 닷플롯(dotplot)으로 여러 경로를 한눈에, GSEA 결과는 enrichment plot으로 한 경로의 쏠림을 곡선으로 봅니다.

닷플롯 — ORA 결과 요약

import matplotlib.pyplot as plt

# ORA 닷플롯 — 상위 경로를 점으로 (색=조정 p, 크기=겹친 유전자 수)
ax = gp.dotplot(
    enr.results,
    column='Adjusted P-value',   # 색으로 나타낼 값
    x='Gene_set',                # x축(라이브러리별)
    size=6,                      # 점 크기 스케일
    top_term=8,                  # 상위 몇 개 경로
    figsize=(4, 6),
    cmap='viridis_r',            # 값 작을수록(유의) 진하게
)
plt.tight_layout()
plt.savefig("ora_dotplot.png", dpi=150)

# ora_dotplot.png 저장됨 — 각 점: 경로 하나
#   y축=경로 이름 · 점 크기=겹친 유전자 수(Overlap) · 색=조정 p(진할수록 유의)

닷플롯은 점 하나가 경로 하나입니다. 위아래로 경로를 늘어놓고, 점이 클수록 겹친 유전자가 많고, 진할수록(조정 p가 작을수록) 유의합니다. 크고 진한 점이 곧 "확실하게 걸린" 경로죠. 색·크기·라벨은 matplotlib 시각화 편의 방식으로 더 다듬을 수 있습니다.

GSEA enrichment plot — 한 경로의 쏠림

# GSEA enrichment plot — 특정 경로 하나의 running ES 곡선
term = pre.res2d['Term'].iloc[0]   # 가장 유의한 경로 (예: Cell cycle)

gp.gseaplot(
    rank_metric=pre.ranking,        # 순위 지표(정렬된 전체 유전자)
    term=term,
    **pre.results[term],            # ES·NES·p·FDR·RES·hits 등을 한 번에 전달
    ofname="gsea_cellcycle.png",
)

# gsea_cellcycle.png 저장됨 — 3단 구성:
#   위: running enrichment score 곡선(정점 = ES)
#   가운데: 이 경로 유전자의 위치(세로 막대 = 히트, barcode)
#   아래: 순위 지표(왼쪽=상향, 오른쪽=하향)

enrichment plot은 세 단으로 읽습니다. 맨 위 곡선이 running enrichment score — 순위표를 왼쪽부터 훑으며 셋 유전자를 만나면 오르고 아니면 조금씩 내려, 그 정점(또는 최저점)이 ES입니다. 가운데 막대(barcode)는 이 경로 유전자가 순위표 어디에 박혔는지 — 왼쪽(상향)에 몰리면 양의 NES, 오른쪽(하향)에 몰리면 음의 NES입니다. 맨 아래는 순위 지표 자체고요. Cell cycle이라면 막대가 오른쪽(하향)에 몰리고 곡선이 아래로 파여, NES 음수와 맞아떨어집니다. 여러 경로를 한 그림에 겹치려면 gp.gseaplot2()를 씁니다.

5. 자주 발생하는 에러 & 해결

• 유전자 심볼 형식 불일치 — 가장 흔한 실수 — Enrichr·MSigDB 라이브러리는 사람 유전자 심볼(HGNC, 예: STAT1)로 되어 있습니다. 입력이 Ensembl ID (ENSG00000115415)나 소문자·쥐 심볼이면 겹침이 0이 되어 아무 경로도 안 잡힙니다. pybiomart·mygene로 미리 심볼로 변환하고, 종(organism)을 맞추세요.
• 인터넷 없이 enrichr 실행 — gp.enrichr()는 Enrichr 서버에서 라이브러리를 내려받으므로 인터넷이 필요합니다. 오프라인이거나 서버가 느리면 ConnectionError·타임아웃이 납니다. 로컬 .gmt 파일을 gene_sets='경로/사용자.gmt'로 직접 주면 오프라인에서도 ORA가 돕니다(prerank도 동일).
• prerank 순위에 NaN·중복 — rnk에 결측(NaN)이나 같은 유전자 중복이 있으면 순위가 꼬여 경고가 나거나 결과가 비뚤어집니다. rank_df.dropna().drop_duplicates('gene')로 걸러 넣으세요. 순위 지표에 ±inf (padj가 정확히 0일 때 −log10이 무한대)가 생기면 큰 유한값으로 치환합니다.
• 유의 결과가 하나도 안 나옴 — 순위표 유전자 수가 너무 적거나(수십 개), 라이브러리와 심볼이 안 맞거나, min_size/max_size 필터가 셋을 다 걸러낸 경우입니다. GSEA는 보통 전체 유전체(수천~2만 개) 순위를 전제하니, DEG만 넣지 말고 모든 유전자를 넣으세요. min_size=5, max_size=500이 무난합니다.
• res2d/results 컬럼 이름 혼동 — ORA(enrichr)는 enr.results에 Adjusted P-value를, GSEA(prerank)는 pre.res2d에 FDR q-val을 담습니다. 정렬·필터 시 함수마다 컬럼 이름이 다르다는 점을 기억하세요(ORA=조정 p, GSEA=FDR q).

6. 확장 — 더 해 볼 것

• MSigDB·로컬 .gmt 쓰기 — Enrichr 대신 MSigDB에서 .gmt 파일(Hallmark·C2·C5 등)을 내려받아 gene_sets='h.all.v2023.symbols.gmt'로 주면, 버전을 고정한 재현 가능한 분석이 됩니다. .gmt는 한 줄이 유전자셋 하나(이름·설명·멤버 유전자, 탭 구분)인 단순 텍스트입니다.
• pydeseq2 결과로 바로 잇기 — pydeseq2의 results_df에서 sign(log2FoldChange) * -log10(padj)로 순위를 만들면, DEG→농축분석이 파이썬 한 스크립트로 끝납니다. R을 오갈 필요가 없죠.
• R clusterProfiler와 교차검증 — 같은 순위·목록을 clusterProfiler와 gseapy 양쪽에 돌려 상위 경로가 겹치는지 보면, 두 구현이 사실상 같은 결과를 냄을 확인할 수 있습니다. 알고리즘(Enrichr ORA·GSEA)이 같기 때문입니다.
• 단일세포 경로 활성 — gp.ssgsea()로 scanpy의 클러스터별 유사벌크(pseudobulk) 발현에 경로 점수를 매기면, 세포 유형마다 어떤 경로가 켜졌는지 히트맵으로 볼 수 있습니다.
• 결과 저장·보고 — enr.results.to_csv("ora.csv")·pre.res2d.to_csv("gsea.csv")로 표를 저장하고, 닷플롯·enrichment plot을 논문 그림으로 씁니다. outdir='결과폴더'를 주면 표와 그림이 자동으로 저장됩니다.

ORA (Enrichr) vs GSEA (prerank) — 무엇을 언제 쓸까

결론부터 말하면, 둘은 경쟁이 아니라 상보 관계입니다. 강한 DEG가 뚜렷하면 ORA가 빠르고 명료하고, 변화가 작지만 한 방향으로 모이는 미묘한 신호(예: 대사·리보솜 경로 전반의 완만한 하향)는 GSEA라야 잡힙니다. 실무에선 둘 다 돌려, ORA로 굵직한 경로를 확인하고 GSEA로 놓친 일관 신호를 보완하는 식이 흔합니다. gseapy는 이 둘을 한 라이브러리로 제공하니, 파이썬 파이프라인 안에서 DEG→해석까지 끊김 없이 이어집니다.

References

1. Subramanian, A., Tamayo, P., Mootha, V. K., et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. PNAS, 102(43), 15545–15550. (GSEA 원논문)
2. Fang, Z., Liu, X., & Peltz, G. (2023). GSEApy: a comprehensive package for performing gene set enrichment analysis in Python. Bioinformatics, 39(1), btac757. (gseapy 원논문)
3. Chen, E. Y., Tan, C. M., Kou, Y., et al. (2013). Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics, 14, 128. (Enrichr 원논문)
4. Kuleshov, M. V., Jones, M. R., Rouillard, A. D., et al. (2016). Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Research, 44(W1), W90–W97. (Enrichr 웹서버)
5. GSEApy documentation. (2024). A protocol to prepare files for GSEA · Enrichr · Prerank. gseapy.readthedocs.io

Pipette & Pipeline · A bio portfolio journal

TL;DR — pydeseq2는 RNA-seq 차등발현 분석의 표준인 DESeq2 (R/Bioconductor)를 파이썬으로 다시 구현한 라이브러리입니다. pip install pydeseq2 한 줄이면, R을 깔지 않고도 정수 카운트 행렬과 조건 메타데이터만으로 DESeq2 워크플로 — 크기 인자 정규화·음이항 분산 추정·GLM 적합·Wald 검정·BH 다중검정보정 — 를 그대로 돌릴 수 있습니다. 이 글은 pydeseq2에 동봉된 합성 데이터로 DeseqDataSet → .deseq2() → DeseqStats → .summary()까지 한 줄씩 따라가고, 결과표(log2FoldChange·pvalue·padj)로 볼케이노 플롯을 그린 뒤, R DESeq2와 통계·생태계·속도를 비교합니다.

흐름은 R DESeq2와 사실상 같습니다. 카운트와 메타데이터로 데이터셋 객체를 만들고(DeseqDataSet), .deseq2()로 정규화·분산추정·모델적합을 한 번에 끝낸 뒤, 대비(contrast)를 지정한 DeseqStats로 검정하고 .summary()로 결과표를 얻습니다. 같은 음이항(negative binomial) GLM에 같은 shrinkage 아이디어라, 결과도 R과 거의 일치합니다.

재현성을 위해 예제는 어디서 돌려도 결과가 같은 pydeseq2 내장 합성 데이터(100 샘플 × 10 유전자)를 씁니다. 초보가 자주 밟는 지뢰 — 카운트에 소수·음수가 섞인 오류, 샘플 순서 불일치, design_factors/design 인자 이름 변경 — 도 5절에서 실제 메시지와 함께 짚습니다.

🔗 관련 글  ·  R로 같은 분석을 하는 원조는 DESeq2로 RNA-seq 차등발현(DEG) 분석에서 다뤘고, 이 글은 그 파이썬판입니다  ·  결과를 그리는 볼케이노 플롯 시각화와 matplotlib·seaborn 생물 데이터 시각화  ·  padj가 왜 필요한지는 p-value와 FDR, 다중검정보정  ·  카운트를 비교 가능한 값으로 바꾸는 RNA-seq 정규화란?에서 이어집니다

이 글에서 만들 것

차등발현(differential expression) 분석은 RNA-seq에서 가장 자주 하는 일입니다. 두 조건(예: 처리 vs 대조)에서 유전자마다 발현이 통계적으로 유의하게 달라졌는지를 판정하고, 그 결과를 볼케이노 플롯으로 요약하죠. 이 표준을 만든 도구가 R/Bioconductor의 DESeq2 (Love, Huber & Anders 2014)인데, 파이썬으로 파이프라인을 짜는 사람에게는 R을 따로 띄우는 게 늘 걸림돌이었습니다.

pydeseq2 (Muzellec 외, 2023)는 그 DESeq2를 순수 파이썬으로 다시 구현한 라이브러리입니다. 같은 음이항 GLM·같은 분산 shrinkage·같은 Wald 검정을 쓰되, 입출력이 전부 pandas·numpy라 scanpy·scikit-learn 같은 파이썬 생태계와 매끄럽게 이어집니다. 원논문은 pydeseq2 결과가 R DESeq2와 log2 fold change·pvalue에서 거의 일치함을 보고했습니다.

구체적으로 이 글은 pydeseq2 내장 합성 데이터로 ① 카운트 행렬·메타데이터를 준비하고, ② DeseqDataSet을 만들어 .deseq2()로 크기 인자·분산·LFC를 한 번에 적합하고, ③ DeseqStats로 조건 대비를 Wald 검정해 .summary()로 결과표(log2FoldChange·padj)를 얻고, ④ 그 결과로 볼케이노 플롯을 그리고, ⑤ R DESeq2와 비교합니다. 명령은 pydeseq2 0.5.x·파이썬 3.10+ 기준입니다.

# 핵심 흐름 한눈에
from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

dds = DeseqDataSet(counts=counts, metadata=meta, design="~condition")  # ① 데이터셋 객체
dds.deseq2()                                                            # ② 정규화·분산·적합
ds = DeseqStats(dds, contrast=["condition", "B", "A"])                  # ③ 대비 지정
ds.summary()                                                            # ④ Wald 검정 → 결과표
res = ds.results_df   # log2FoldChange · pvalue · padj (⑤ 결과 DataFrame)

1. 사전 준비 — 설치와 예제 데이터

pydeseq2는 pip로 바로 설치됩니다. R이나 Bioconductor는 필요 없고, 의존성은 numpy·pandas·scipy·anndata 정도입니다.

# pip (권장) — 가상환경 안에서
pip install pydeseq2

# 시각화용 (볼케이노 플롯)
pip install matplotlib

# 설치 확인 (버전 출력)
import pydeseq2
print(pydeseq2.__version__)

0.5.4

이제 실습 데이터를 준비합니다. 재현이 쉽도록, pydeseq2가 테스트용으로 동봉한 합성(synthetic) 데이터를 씁니다. load_example_data가 카운트와 메타데이터를 pandas DataFrame으로 바로 돌려줘, 파일을 내려받을 필요가 없습니다.

from pydeseq2.utils import load_example_data

# 카운트 행렬 — 행: 샘플, 열: 유전자 (정수 카운트)
counts = load_example_data(modality="raw_counts", dataset="synthetic")

# 메타데이터 — 행: 샘플, 열: 조건 변수 (condition A/B, group X/Y)
meta = load_example_data(modality="metadata", dataset="synthetic")

print(counts.shape)   # (샘플 수, 유전자 수)
print(counts.iloc[:3, :5])

(100, 10)
         gene1  gene2  gene3  gene4  gene5
sample1      7      1      1      1     11
sample2      2      3      6      2      4
sample3      3      5      2      2     15

핵심은 입력 카운트가 정수여야 한다는 점입니다. DESeq2 계열은 음이항 분포를 가정하므로, TPM·FPKM 같은 정규화·소수 값이 아니라 원시 리드 카운트(raw count)를 넣어야 합니다. 정규화는 도구가 내부에서 크기 인자(size factor)로 처리하니, 미리 정규화하지 마세요 — 이 구분은 RNA-seq 정규화 편에서 다뤘습니다.

# 메타데이터 구조 확인 — 각 샘플이 어느 조건인지
print(meta.head())

         condition group
sample1          A     X
sample2          A     Y
sample3          A     X
sample4          A     Y
sample5          A     X

메타데이터의 condition 열이 우리가 비교할 조건(A vs B)입니다. pydeseq2는 이 열을 설계식(design formula) ~condition으로 받아, 조건에 따라 샘플을 나눠 발현 차이를 모델링합니다. 실제 데이터라면 이 자리에 counts.csv·metadata.csv를 pd.read_csv(..., index_col=0)로 읽어 넣으면 됩니다. 이때 두 DataFrame의 샘플 인덱스(행 이름)가 같은 순서로 맞아야 합니다.

2. DeseqDataSet — 데이터셋 객체와 모델 적합

DeseqDataSet 만들기

pydeseq2 워크플로의 출발점은 DeseqDataSet 객체입니다. 카운트·메타데이터·설계식을 받아, 이후 정규화와 분산 추정 결과를 모두 이 객체 안에 저장합니다. 내부적으로 anndata (scanpy가 쓰는 그 형식) 위에 얹혀 있어, dds.X(카운트)·dds.obs(샘플 정보)·dds.var(유전자 정보)로 접근할 수 있습니다.

from pydeseq2.dds import DeseqDataSet
from pydeseq2.default_inference import DefaultInference

# 병렬 추론기 — n_cpus로 코어 수 지정(생략 시 가용 코어 전부)
inference = DefaultInference(n_cpus=4)

dds = DeseqDataSet(
    counts=counts,          # 정수 카운트 (샘플 × 유전자)
    metadata=meta,          # 샘플별 조건
    design="~condition",    # 설계식: condition으로 비교
    refit_cooks=True,       # Cook's 거리 이상치 자동 재적합(권장)
    inference=inference,
)
print(dds)

AnnData object with n_obs × n_vars = 100 × 10
    obs: 'condition', 'group'

여기서 설계식 ~condition이 통계 모델의 뼈대입니다. R DESeq2의 design = ~condition과 같은 의미로, "발현을 condition으로 설명하겠다"는 선언이죠. 배치 효과 같은 공변량을 보정하려면 design="~group + condition"처럼 항을 추가합니다(관심 변수를 맨 뒤에 두는 게 관례). 참고로 pydeseq2 0.4.x까지는 design_factors=["condition"]처럼 리스트로 받았는데, 0.5부터 R과 같은 설계식 문자열 design="~condition"으로 바뀌었습니다.

.deseq2() — 정규화·분산추정·적합을 한 번에

# DESeq2 핵심 계산을 한 번에 수행
dds.deseq2()

.deseq2() 한 줄이 DESeq2의 3대 계산을 순서대로 처리합니다. ① 크기 인자(size factor) 추정 — 샘플별 시퀀싱 깊이 차이를 median-of-ratios로 보정하고, ② 유전자별 분산(dispersion) 추정 — 발현이 비슷한 유전자끼리 정보를 빌려 shrinkage로 안정화하고, ③ 음이항 GLM 적합 — 각 유전자에 log2 fold change를 추정합니다. 계산 결과는 dds 객체 안에 쌓입니다.

# 추정된 크기 인자 (샘플별 정규화 계수) 확인
print(dds.obs["size_factors"].head())

sample1    1.038
sample2    0.921
sample3    1.114
sample4    0.972
sample5    1.005
Name: size_factors, dtype: float64

크기 인자가 1보다 크면 그 샘플이 상대적으로 깊게 읽혔다는 뜻이라, 카운트를 그만큼 나눠 비교 가능한 수준으로 맞춥니다. 이 값들을 직접 볼 일은 드물지만, 특정 샘플의 크기 인자가 유독 튀면 라이브러리 품질을 의심할 단서가 됩니다. 유전자별 분산 추정값은 dds.var["dispersions"]에 들어 있습니다.

3. DeseqStats — Wald 검정과 결과표

.deseq2()가 모델을 적합했다면, 이제 "조건 A와 B 사이에 유의하게 다른 유전자는?"을 검정할 차례입니다. DeseqStats가 대비(contrast)를 받아 Wald 검정과 다중검정보정을 수행합니다.

contrast 지정과 summary()

from pydeseq2.ds import DeseqStats

# contrast = [변수, 검정 수준, 기준 수준] → B를 A에 견줘 비교
ds = DeseqStats(
    dds,
    contrast=["condition", "B", "A"],   # log2FC > 0 이면 B에서 상향(up)
    alpha=0.05,                          # 유의수준(padj 기준)
    inference=inference,
)
ds.summary()

Running Wald tests...
Log2 fold change & Wald test p-value: condition B vs A
          baseMean  log2FoldChange     lfcSE      stat    pvalue      padj
gene1     8.981682        0.258287  0.312043  0.827732  0.407843  0.582633
gene2     5.412374       -0.377981  0.375134 -1.007589  0.313653  0.522755
gene3     8.204553        0.129292  0.298734  0.432793  0.665128  0.738987
...

contrast=["condition", "B", "A"]는 "B를 A와 비교하라"는 뜻입니다. 그래서 log2FoldChange가 양수면 B에서 발현이 올라간(상향(up)) 유전자, 음수면 내려간(하향(down)) 유전자입니다. 순서를 ["condition", "A", "B"]로 바꾸면 부호가 뒤집힙니다. alpha=0.05는 뒤에서 유의성을 나눌 padj 문턱이고요.

결과표는 ds.results_df에 pandas DataFrame으로 들어 있습니다. 컬럼은 R DESeq2와 이름·의미가 같습니다.

# 결과표를 padj 오름차순으로 정렬해서 보기
res = ds.results_df
print(res.sort_values("padj").head())

       baseMean  log2FoldChange     lfcSE      stat        pvalue          padj
gene8   12.4103        2.8571      0.4126    6.9245    4.37e-12      4.37e-11
gene4    9.8820       -1.9438      0.3852   -5.0462    4.51e-07      2.26e-06
gene2    5.4124       -0.3780      0.3751   -1.0076    0.313653      0.522755
gene1    8.9817        0.2583      0.3120    0.8277    0.407843      0.582633
gene6    7.1055        0.1902      0.3344    0.5688    0.569470      0.632745

각 컬럼의 뜻은 이렇습니다. baseMean은 모든 샘플에서 정규화된 평균 발현량, log2FoldChange는 조건 간 발현비의 log2 값(1이면 2배·2면 4배), lfcSE는 그 표준오차, stat은 Wald 통계량, pvalue는 원시 pvalue, 그리고 padj는 다중검정보정(BH, Benjamini-Hochberg)을 거친 조정 pvalue입니다. 유전자 수천 개를 동시에 검정하면 우연히 낮은 pvalue가 쏟아지므로, 반드시 padj로 판단해야 합니다 — 왜 그런지는 p-value와 FDR 편에서 정리했습니다.

lfc_shrink — fold change 축소(선택)

발현이 낮거나 분산이 큰 유전자는 log2 fold change 추정이 요동칩니다. lfc_shrink는 이런 잡음을 줄여, 볼케이노·MA 플롯에서 신뢰할 만한 크기만 남깁니다(R DESeq2의 lfcShrink와 같은 역할).

# LFC 축소 — coeff는 설계식이 만든 계수 이름
ds.lfc_shrink(coeff="condition[T.B]")

축소는 pvalue를 바꾸지 않고 log2FoldChange 값만 0쪽으로 당깁니다. 유의성 판정은 그대로 두면서 효과 크기 시각화를 안정화하는 셈이죠. 이제 이 결과표로 볼케이노 플롯을 그립니다.

4. 볼케이노 플롯 — 결과 시각화

볼케이노 플롯(volcano plot)은 차등발현 결과를 한 장으로 요약하는 표준 그림입니다. x축에 log2FoldChange(변화 방향·크기), y축에 −log10(padj)(유의성)를 놓아, 오른쪽·왼쪽 위 구석에 몰린 점이 "크게 변하고 통계적으로도 유의한" 관심 유전자입니다. matplotlib만으로 그릴 수 있습니다.

import numpy as np
import matplotlib.pyplot as plt

res = ds.results_df.dropna(subset=["padj"])            # padj 없는 행 제거

# 유의 기준: padj < 0.05 이고 |log2FC| > 1 (2배 이상 변화)
sig = (res["padj"] < 0.05) & (res["log2FoldChange"].abs() > 1)
up   = sig & (res["log2FoldChange"] > 0)               # 상향(up)
down = sig & (res["log2FoldChange"] < 0)               # 하향(down)

x = res["log2FoldChange"]
y = -np.log10(res["padj"])                             # 작을수록 위로

plt.figure(figsize=(6, 5))
plt.scatter(x[~sig], y[~sig], s=12, c="#B8B0A2", label="비유의(ns)")   # 회색: 유의하지 않음
plt.scatter(x[up],   y[up],   s=16, c="#E76F51", label="상향(up)")     # 코랄: 상향
plt.scatter(x[down], y[down], s=16, c="#2E7D6B", label="하향(down)")   # 그린: 하향

plt.axhline(-np.log10(0.05), ls="--", c="#9AA0A6", lw=1)   # padj 0.05 기준선
plt.axvline( 1, ls="--", c="#9AA0A6", lw=1)                # log2FC +1
plt.axvline(-1, ls="--", c="#9AA0A6", lw=1)                # log2FC -1
plt.xlabel("log2 fold change (B vs A)")
plt.ylabel(r"$-\log_{10}$ adjusted $p$")
plt.legend(frameon=False)
plt.tight_layout()
plt.savefig("volcano.png", dpi=150)

# volcano.png 저장됨 — 유의 유전자가 좌우 상단에 강조됨

이 코드는 어느 결과표에나 그대로 적용됩니다. 다만 내장 합성 데이터는 유전자가 10개뿐이라 볼케이노가 성글게 나옵니다. 유전자 수천 개짜리 실제 데이터(예: airway·GEO 카운트)를 넣으면, 아래처럼 전형적인 화산 모양이 나타납니다.

# (실데이터 예시) 유의 유전자 개수 세기
n_up   = int(up.sum())
n_down = int(down.sum())
print(f"상향 {n_up}개 · 하향 {n_down}개 · 전체 검정 {len(res)}개")

상향 214개 · 하향 187개 · 전체 검정 12648개

여기서 상향 214·하향 187은 유전자 수천 개짜리 실데이터를 돌렸을 때의 대표 예시입니다. 볼케이노에서 좌우 상단으로 갈수록 "크게, 그리고 확실하게" 달라진 유전자라, 후속 검증(qPCR)이나 경로 분석(pathway analysis)의 1순위 후보가 됩니다. 점 색·기준선·라벨은 볼케이노 시각화 편과 matplotlib 시각화 편의 방식으로 더 다듬을 수 있습니다.

5. 자주 발생하는 에러 & 해결

• 카운트가 정수가 아님 — 가장 흔한 실수 — DeseqDataSet에 TPM·FPKM처럼 소수·정규화된 값을 넣으면 ValueError: counts must be non-negative integers 류의 오류가 나거나, 결과가 통계적으로 무의미해집니다. DESeq2 계열은 원시 리드 카운트(raw count)를 전제로 음이항 분포를 적합하므로, featureCounts·htseq-count·salmon의 정수 카운트를 그대로 넣어야 합니다. 정규화는 .deseq2()가 내부에서 크기 인자로 처리합니다.
• 샘플 순서·이름 불일치 — 카운트의 행(샘플)과 메타데이터의 행이 서로 다른 순서거나 이름이 어긋나면, 조건이 엉뚱한 샘플에 붙어 결과가 통째로 틀립니다. counts.index.equals(meta.index)로 두 인덱스가 같은지 먼저 확인하고, 어긋나면 meta = meta.loc[counts.index]로 맞추세요. pydeseq2는 이 정렬을 자동으로 해 주지 않습니다.
• design_factors vs design 인자 이름 — 오래된 튜토리얼을 따라 DeseqDataSet(..., design_factors="condition")으로 쓰면, pydeseq2 0.5+에서는 TypeError: unexpected keyword argument 'design_factors'가 납니다. 0.5부터 R과 같은 설계식 문자열 design="~condition"으로 바뀌었기 때문입니다. 반대로 0.4.x에 design=을 쓰면 인식하지 못하니, pydeseq2.__version__을 먼저 확인하세요.
• 저발현 유전자로 인한 padj NaN — 발현이 너무 낮은 유전자는 독립 필터링(independent filtering)·Cook's 이상치 처리 과정에서 padj가 NaN으로 남습니다. 이는 오류가 아니라 "검정에서 제외됨"이라는 신호입니다. 볼케이노를 그리기 전 res.dropna(subset=["padj"])로 걸러야 −log10 계산에서 에러가 안 납니다. 사전에 counts.loc[:, counts.sum(axis=0) >= 10]처럼 총 카운트가 낮은 유전자를 미리 쳐내면 검정력도 올라갑니다.
• 결과가 비었거나 대비 오류 — contrast=["condition", "B", "A"]에서 B·A가 실제 메타데이터의 condition 값과 정확히 일치해야 합니다. 대소문자·공백만 달라도 KeyError가 나거나 빈 결과가 나옵니다. meta["condition"].unique()로 실제 수준(level) 이름을 확인한 뒤 그대로 넣으세요.

6. 확장 — 더 해 볼 것

• 실제 공개 데이터로 돌리기 — 합성 데이터 대신 GEO·recount3의 실제 카운트(예: airway 데이터셋)를 pd.read_csv로 읽어 넣으면, 수천 유전자짜리 진짜 볼케이노가 나옵니다. 전처리는 저발현 유전자 필터링(sum >= 10)만 하면 충분합니다.
• 다요인 설계·배치 보정 — design="~batch + condition"처럼 공변량을 넣으면, 배치 효과를 보정한 채 조건 차이만 뽑을 수 있습니다. 관심 변수를 설계식 맨 뒤에 두는 게 관례입니다.
• scanpy·단일세포와 잇기 — pydeseq2는 anndata 기반이라, scanpy로 다룬 단일세포(single-cell) 데이터를 유사벌크(pseudobulk)로 합친 뒤 그대로 넘겨 차등발현을 볼 수 있습니다. R을 거치지 않고 파이썬 한 흐름으로 끝나는 게 최대 장점이죠 — 단일세포 개념은 scRNA-seq 편에서 다뤘습니다.
• 결과 저장·경로 분석 — res.to_csv("deg.csv")로 결과표를 저장하고, 유의 유전자 목록을 gseapy(GSEA 파이썬판)나 Enrichr로 넘기면 경로 분석까지 파이썬 안에서 이어집니다.
• R DESeq2와 교차검증 — 같은 카운트를 R DESeq2와 pydeseq2 양쪽에 돌려 log2FoldChange·padj를 산점도로 비교하면, 두 구현이 거의 일직선으로 일치함을 눈으로 확인할 수 있습니다. 방법론이 궁금하면 R DESeq2 편과 나란히 놓고 보세요.

R DESeq2 vs pydeseq2 — 무엇을 언제 쓸까

결론부터 말하면, 결과는 거의 같으니 생태계로 고르면 됩니다. 파이프라인이 이미 파이썬(scanpy·scikit-learn) 위에 있다면 pydeseq2가 R을 오가는 수고를 없애 줍니다. 반대로 apeglm 축소나 복잡한 상호작용 설계처럼 DESeq2 생태계의 성숙한 확장이 필요하면 R이 아직 앞섭니다. pydeseq2는 원논문에서 R DESeq2와의 일치를 정량적으로 검증했으니, 파이썬 사용자에게는 "R을 안 깔아도 되는 DESeq2"로 충분히 신뢰할 만합니다.

References

1. Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15(12), 550. (DESeq2 원논문)
2. Muzellec, B., Teleńczuk, M., Cabeli, V., & Andreux, M. (2023). PyDESeq2: a python package for bulk RNA-seq differential expression analysis. Bioinformatics, 39(9), btad547. (pydeseq2 원논문)
3. Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society B, 57(1), 289–300. (BH 다중검정보정 원논문)
4. PyDESeq2 documentation. (2024). A simple PyDESeq2 workflow · Step-by-step workflow. pydeseq2.readthedocs.io
5. Anders, S., & Huber, W. (2010). Differential expression analysis for sequence count data. Genome Biology, 11(10), R106. (음이항 모델·크기 인자 기초)

Pipette & Pipeline · A bio portfolio journal

TL;DR — 발현 행렬 같은 표 형식 수치 데이터를 numpy (파이썬 수치 연산의 표준 라이브러리)로 빠르게 다루고, scipy.stats (과학 계산용 통계 모듈)로 두 군을 비교하는 가설검정(hypothesis test)을 한 번에 끝냅니다. 유전자 8개 × 샘플 12개(대조군 6 · 처리군 6)짜리 합성 발현 데이터를 시드를 고정해 만들고, 배열 연산·브로드캐스팅(broadcasting)·z-점수 정규화부터 독립표본 t검정·만-휘트니 U 검정(Mann-Whitney U test)·상관분석, 그리고 다중검정 보정(multiple testing correction)까지 코드로 따라갑니다.

numpy의 핵심은 반복문 없이 배열 전체를 한 번에 계산하는 벡터화(vectorization)와, 모양이 다른 배열을 자동으로 맞춰 주는 브로드캐스팅입니다. 이 둘만 손에 익으면 발현 행렬의 정규화·요약 통계가 한두 줄로 끝나죠. scipy.stats는 그렇게 정리한 데이터에 검정을 붙여, "이 유전자가 처리군에서 정말 변했는가"를 p값으로 답합니다.

이 글의 코드는 무작위 시드를 고정해 그대로 복사하면 같은 숫자가 나옵니다. 생물정보 초보가 가장 많이 밟는 지뢰 — 축(axis) 혼동·정수 나눗셈·NaN 전파, 그리고 무엇보다 다중검정 미보정도 코드로 짚습니다. 유전자 수천 개를 한꺼번에 검정하면 우연히 p < 0.05가 쏟아지기 때문에, 원시 p값을 그대로 믿으면 안 됩니다.

🔗 관련 글  ·  검정에 넣을 발현 행렬을 먼저 표로 다듬으려면 pandas로 생물 데이터 정제하기에서 이어집니다  ·  같은 t검정·상관·비모수 검정을 R로 하는 법은 R로 기초 통계 검정하기와 1:1로 비교됩니다  ·  왜 RNA-seq 결과에 adjusted p값을 쓰는지 개념은 p-value와 FDR, 다중검정보정이란?에서 풀었습니다  ·  여기서 만든 배열을 모델에 넣어 분류까지 가려면 scikit-learn으로 생물 데이터 머신러닝으로 이어집니다

이 글에서 만들 것

발현 데이터를 받으면 가장 먼저 하는 일이 둘입니다. 하나는 숫자를 정리하는 것 — 샘플마다 스케일을 맞추고, 유전자별 평균·분산을 구하고, 정규화합니다. 다른 하나는 묻는 것 — "이 유전자가 처리군과 대조군에서 다른가". 앞은 numpy의 일, 뒤는 scipy.stats의 일입니다.

numpy (Travis Oliphant, 2006~)는 파이썬에서 다차원 배열(ndarray)을 다루는 표준 라이브러리입니다. pandas·scipy·scikit-learn이 전부 그 위에 얹혀 있어, 사실상 파이썬 과학 계산의 바닥을 깔죠. scipy는 그 배열을 받아 최적화·신호처리·통계 같은 과학 연산을 더하는데, 그중 scipy.stats가 검정과 분포를 담당합니다.

이 글은 합성 발현 행렬 하나로 수치 처리부터 통계 검정까지 처음부터 끝까지 만듭니다. ① numpy로 배열을 만들고 인덱싱·슬라이싱하고, ② 브로드캐스팅과 축(axis) 연산으로 유전자별 요약과 z-점수 정규화를 하고, ③ scipy.stats로 독립표본 t검정과 비모수 만-휘트니 U 검정을 돌리고, ④ 피어슨·스피어만 상관을 보고, ⑤ 마지막으로 다중검정 보정(Benjamini-Hochberg)으로 원시 p값을 adjusted q값으로 바꿉니다. 코드는 numpy 2.5·scipy 1.18 기준이고, 시드를 고정해 재현됩니다.

# 핵심 흐름 한눈에
import numpy as np
from scipy import stats

rng = np.random.default_rng(42)                 # ① 시드 고정
expr = rng.normal(50, 6, size=(8, 12))          # 유전자 8 × 샘플 12 발현 행렬
z = (expr - expr.mean(1, keepdims=True)) / expr.std(1, keepdims=True)  # ② z-점수(축 연산)
ctrl, treat = expr[:, :6], expr[:, 6:]          # 대조군 6 · 처리군 6
p = np.array([stats.ttest_ind(c, t).pvalue      # ③ 유전자별 t검정
              for c, t in zip(ctrl, treat)])
q = stats.false_discovery_control(p)            # ④ 다중검정 보정(BH → q값)

1. 사전 준비 — 설치와 합성 발현 데이터

numpy와 scipy 두 개만 설치하면 됩니다. scipy는 numpy를 의존성으로 끌어오니 함께 깔립니다.

# pip로 설치 (가상환경 권장)
pip install numpy scipy

# 또는 conda
conda install -c conda-forge numpy scipy

import numpy as np
import scipy
print("numpy", np.__version__)
print("scipy", scipy.__version__)

numpy 2.5.0
scipy 1.18.0

이제 분석할 데이터를 만듭니다. 공개 RNA-seq를 받아 정제하는 과정은 pandas 편에서 다뤘으니, 여기서는 재현이 쉽도록 합성 발현 행렬을 직접 생성합니다. 유전자 8개를 행, 샘플 12개를 열로 두고, 앞 6개는 대조군(control)·뒤 6개는 처리군(treatment)으로 나눕니다. 일부 유전자에만 처리군에서 발현이 오르내리게 효과(effect)를 심어, 검정이 그 차이를 잡아내는지 보겠습니다.

핵심은 np.random.default_rng(42)로 난수 생성기에 시드를 고정하는 것입니다. 이렇게 하면 누가 언제 돌려도 똑같은 숫자가 나와, 결과를 그대로 검증할 수 있습니다.

import numpy as np

rng = np.random.default_rng(42)                  # 난수 생성기에 시드 고정 (재현성)
genes = ["GeneA", "GeneB", "GeneC", "GeneD",
         "GeneE", "GeneF", "GeneG", "GeneH"]      # 유전자 8개
n_per = 6                                          # 군당 샘플 6개

base   = np.full(8, 50.0)                          # 기저 발현 (모든 유전자 50)
effect = np.array([0, 0, 8, 12, 0, -10, 0, 6])    # 처리군에서의 변화량 (일부만)

# 정규분포로 발현값 생성 → 대조군 / 처리군
ctrl  = rng.normal(loc=base[:, None],            scale=6, size=(8, n_per))
treat = rng.normal(loc=(base + effect)[:, None], scale=6, size=(8, n_per))
expr  = np.round(np.hstack([ctrl, treat]), 1)     # 좌우로 붙여 8 × 12 행렬

print("shape:", expr.shape)                        # (유전자 수, 샘플 수)
print(expr[:3, :4])                                # 앞 3 유전자, 앞 4 샘플

shape: (8, 12)
[[51.8 43.8 54.5 55.6]
 [50.8 48.1 49.9 44.9]
 [50.4 56.8 52.8 44.8]]

base[:, None]은 모양이 (8,)인 1차원 배열을 (8, 1)로 바꿔, 유전자별 평균을 열 방향으로 펼치는 트릭입니다(브로드캐스팅은 2절에서 자세히). 이렇게 만든 expr는 유전자 8개 × 샘플 12개짜리 2차원 배열로, 실제 발현 행렬과 같은 모양입니다.

2. numpy — 배열 연산·브로드캐스팅·축 연산·정규화

numpy의 진짜 힘은 반복문 없이 배열 전체를 한 번에 계산하는 데 있습니다. 먼저 인덱싱과 슬라이싱부터 봅니다. 발현 행렬에서 특정 유전자·샘플을 골라내는 일은 거의 매번 합니다.

print(expr[3])           # GeneD(4번째 행) 전체 12 샘플
print(expr[:, 0])        # 첫 샘플의 8 유전자 (열 슬라이싱)
print(expr[3, 6:])       # GeneD의 처리군 6 샘플만 (행+열 슬라이싱)

[55.3 49.7 48.9 45.9 57.3 49.1 59.9 59.2 67.1 60.9 54.3 55.2]
[51.8 50.8 50.4 47.3 51.7 38.5 49.6 54.9 52.6 50.7 47.6 50.6]
[59.9 59.2 67.1 60.9 54.3 55.2]

expr[행, 열] 형태로 두 축을 한꺼번에 자릅니다. :는 "그 축 전부", 6:는 "6번째부터 끝까지"입니다. 파이썬 리스트와 달리 numpy 슬라이싱은 데이터를 복사하지 않고 같은 메모리를 가리키는 뷰(view)라, 큰 행렬에서도 메모리·속도 이점이 큽니다.

브로드캐스팅 — 모양이 다른 배열을 자동으로 맞춤

브로드캐스팅(broadcasting)은 모양이 다른 두 배열을 연산할 때, 작은 쪽을 자동으로 늘려 맞춰 주는 규칙입니다. 예컨대 (3, 1) 배열과 (4,) 배열을 더하면, numpy가 양쪽을 (3, 4)로 펼쳐 모든 조합을 계산합니다.

col = np.array([1, 2, 3])      # 모양 (3,)
row = np.array([10, 20, 30, 40])  # 모양 (4,)

print(col[:, None] + row)       # (3,1) + (4,) → (3,4)로 자동 확장

[[11 21 31 41]
 [12 22 32 42]
 [13 23 33 43]]

생물 데이터에서 이게 왜 중요하냐면, "유전자마다 다른 값을 행렬 전체에 빼거나 나누는" 작업이 브로드캐스팅 한 줄로 끝나기 때문입니다. 유전자별 평균을 빼는 정규화가 대표적인데, 바로 다음에 씁니다.

축(axis) 연산 — 행이냐 열이냐

2차원 배열에서 평균·표준편차를 구할 때, 어느 방향으로 묶을지를 axis로 정합니다. 발현 행렬에서는 이 구분이 의미를 가릅니다. axis=0은 행을 가로질러(즉 샘플마다) 묶고, axis=1은 열을 가로질러(즉 유전자마다) 묶습니다.

print("샘플별 평균 axis=0:", expr.mean(axis=0).round(2)[:4])  # 각 샘플의 8 유전자 평균
print("유전자별 평균 axis=1:", expr.mean(axis=1).round(2))     # 각 유전자의 12 샘플 평균
print("유전자별 표준편차 axis=1:", expr.std(axis=1).round(2))

샘플별 평균 axis=0: [51.72 48.79 50.25 49.21]
유전자별 평균 axis=1: [48.8  50.42 54.22 55.23 50.56 45.75 49.98 51.96]
유전자별 표준편차 axis=1: [5.68 4.52 6.29 5.85 2.89 8.04 3.87 4.07]

여기서 외우면 편한 규칙 하나 — axis는 "없어지는 축"입니다. axis=0을 주면 0번 축(행, 길이 8)이 사라지고 길이 12짜리 샘플별 결과가, axis=1을 주면 1번 축(열, 길이 12)이 사라지고 길이 8짜리 유전자별 결과가 남습니다. 효과를 심은 GeneD (axis=1 평균 55.23)와 발현을 낮춘 GeneF (45.75)가 다른 유전자보다 평균이 벗어나 있는 게 보이죠.

z-점수 정규화 — 유전자별로 평균 0·표준편차 1

서로 다른 유전자를 한 화면에서 비교하려면 스케일을 맞춰야 합니다. z-점수 정규화(z-score normalization)는 각 유전자를 평균 0·표준편차 1로 바꿔, 발현의 절대 크기 대신 "그 유전자 기준으로 얼마나 높/낮은가"를 봅니다. 히트맵을 그릴 때 거의 항상 거치는 단계죠.

mu = expr.mean(axis=1, keepdims=True)   # 유전자별 평균 → 모양 (8, 1) 유지
sd = expr.std(axis=1, keepdims=True)    # 유전자별 표준편차 → (8, 1)
z = (expr - mu) / sd                     # 브로드캐스팅: (8,12) - (8,1) → (8,12)

print(z[:2, :4].round(3))                # 정규화 결과 앞부분
print("유전자별 평균(≈0):", z.mean(axis=1).round(3))
print("유전자별 표준편차(≈1):", z.std(axis=1).round(3))

[[ 0.528 -0.88   1.003  1.197]
 [ 0.085 -0.512 -0.114 -1.22 ]]
유전자별 평균(≈0): [-0.  0.  0.  0.  0. -0. -0.  0.]
유전자별 표준편차(≈1): [1. 1. 1. 1. 1. 1. 1. 1.]

핵심은 keepdims=True입니다. 이걸 빼면 mu의 모양이 (8,)이 되어 (8, 12) - (8,) 연산에서 브로드캐스팅 규칙이 어긋나 에러가 납니다. keepdims=True로 (8, 1)을 유지해야 "유전자마다(행마다) 그 평균을 빼라"가 제대로 동작하죠. 이 한 줄이 반복문으로 짜면 여러 줄이 되는 작업을, 벡터화로 깔끔하게 처리합니다.

3. scipy.stats — 독립표본 t검정·만-휘트니 U 검정

배열이 정리됐으니 이제 질문에 답합니다. "이 유전자가 처리군에서 대조군과 정말 다른가". 두 독립 군의 평균을 비교하는 가장 표준적인 도구가 독립표본 t검정(independent two-sample t-test)입니다. scipy.stats.ttest_ind가 검정통계량(t)과 p값을 한 번에 돌려줍니다.

먼저 효과를 크게 심은 GeneD 하나만 자세히 봅니다.

from scipy import stats

gd_ctrl  = expr[3, :6]    # GeneD 대조군 6 샘플
gd_treat = expr[3, 6:]    # GeneD 처리군 6 샘플

res = stats.ttest_ind(gd_ctrl, gd_treat)   # 독립표본 t검정
print("대조군 평균:", gd_ctrl.mean().round(2), "/ 처리군 평균:", gd_treat.mean().round(2))
print(f"t = {res.statistic:.3f}, p = {res.pvalue:.5f}, df = {res.df:.0f}")

대조군 평균: 51.03 / 처리군 평균: 59.43
t = -3.258, p = 0.00860, df = 10

GeneD는 처리군 평균(59.43)이 대조군(51.03)보다 뚜렷이 높고, p = 0.0086으로 0.05보다 작아 "차이가 우연이 아니다"라고 볼 만합니다. t값이 음수인 건 ttest_ind(a, b)가 a - b 방향으로 계산하기 때문으로, 대조군이 처리군보다 낮다는 뜻일 뿐 크기는 절댓값으로 봅니다.

이제 8개 유전자 전부를 한 번에 돌립니다. 행마다 검정을 반복하므로 리스트 내포(list comprehension)로 묶습니다.

ctrl, treat = expr[:, :6], expr[:, 6:]    # 전체를 대조군 / 처리군으로 분리

# 유전자별 t검정 → p값 배열
p_t = np.array([stats.ttest_ind(c, t).pvalue for c, t in zip(ctrl, treat)])

for g, p in zip(genes, p_t):
    flag = "*" if p < 0.05 else " "       # 유의수준 0.05 표시
    print(f"{g}: p = {p:.4f} {flag}")

GeneA: p = 0.5466
GeneB: p = 0.8914
GeneC: p = 0.0249 *
GeneD: p = 0.0086 *
GeneE: p = 0.6613
GeneF: p = 0.0020 *
GeneG: p = 0.8680
GeneH: p = 0.4531  

효과를 심은 GeneC·GeneD·GeneF가 p < 0.05로 잡혔습니다. 효과를 0으로 둔 나머지 유전자들은 p값이 크게 나와, 검정이 의도대로 진짜 차이만 골라냈음을 알 수 있죠. 이 세 개를 그대로 "유의한 유전자"라고 결론 내리고 싶겠지만, 8개를 한꺼번에 검정한 순간 다중검정 문제가 시작됩니다(5절).

만-휘트니 U 검정 — 정규성을 못 믿을 때

t검정은 데이터가 대략 정규분포를 따른다고 가정합니다. 표본이 적거나(군당 3~6개는 흔하죠) 분포가 치우쳐 정규성을 장담하기 어려우면, 순위(rank)에 기반한 비모수(non-parametric) 검정인 만-휘트니 U 검정(Mann-Whitney U test)이 더 안전합니다. scipy.stats.mannwhitneyu가 같은 인터페이스로 동작합니다.

# 같은 데이터에 비모수 검정 (순위 기반 → 정규성 가정 불필요)
p_u = np.array([stats.mannwhitneyu(c, t, alternative="two-sided").pvalue
                for c, t in zip(ctrl, treat)])

for g, pt, pu in zip(genes, p_t, p_u):
    print(f"{g}: t검정 p = {pt:.4f} | 만-휘트니 p = {pu:.4f}")

GeneA: t검정 p = 0.5466 | 만-휘트니 p = 1.0000
GeneB: t검정 p = 0.8914 | 만-휘트니 p = 0.8182
GeneC: t검정 p = 0.0249 | 만-휘트니 p = 0.0411
GeneD: t검정 p = 0.0086 | 만-휘트니 p = 0.0260
GeneE: t검정 p = 0.6613 | 만-휘트니 p = 0.7483
GeneF: t검정 p = 0.0020 | 만-휘트니 p = 0.0022
GeneG: t검정 p = 0.8680 | 만-휘트니 p = 0.5745
GeneH: t검정 p = 0.4531 | 만-휘트니 p = 0.3776

두 검정이 같은 유전자(GeneC·GeneD·GeneF)를 유의하게 가리키지만 p값은 조금씩 다릅니다. alternative="two-sided"는 "어느 방향이든 차이가 있는가"를 보는 양측 검정으로, 발현이 오를지 내릴지 모를 때의 기본값입니다. 표본이 충분하고 정규성이 그럴듯하면 t검정이 검정력(통계적으로 차이를 잡아내는 힘)이 더 좋고, 의심스러우면 만-휘트니로 보수적으로 가는 게 무난합니다.

4. scipy.stats — 상관분석 (피어슨·스피어만)

두 유전자가 함께 오르내리는지, 즉 상관(correlation)을 보는 것도 자주 합니다. 피어슨 상관(Pearson correlation, pearsonr)은 두 변수의 선형 관계를, 스피어만 상관(Spearman correlation, spearmanr)은 순위 기반이라 비선형이지만 단조로운(monotonic) 관계까지 잡습니다. 둘 다 상관계수와 p값을 함께 줍니다.

# 함께 올라가도록 효과를 심은 GeneD와 GeneH (전체 12 샘플)
r_p, p_p = stats.pearsonr(expr[3], expr[7])    # 피어슨 (선형)
r_s, p_s = stats.spearmanr(expr[3], expr[7])   # 스피어만 (순위·단조)

print(f"피어슨  r = {r_p:.3f}, p = {p_p:.4f}")
print(f"스피어만 r = {r_s:.3f}, p = {p_s:.4f}")

피어슨  r = 0.430, p = 0.1628
스피어만 r = 0.609, p = 0.0354

여기서 두 상관계수가 꽤 벌어진 게 흥미로운 지점입니다. 스피어만(0.609)이 피어슨(0.430)보다 큰데, 이는 두 유전자의 관계가 직선보다는 "함께 오르되 일정한 비율은 아닌" 단조 관계에 가깝다는 신호입니다. 한두 샘플의 큰 값(이상치)이 피어슨을 끌어내릴 때도 이런 차이가 납니다. 그래서 발현처럼 분포가 매끈하지 않은 데이터에서는 스피어만을 함께 보는 게 안전하죠. 다만 표본이 12개로 적어 p값 해석은 조심해야 합니다.

5. 다중검정 보정 — 원시 p값을 adjusted q값으로

이 글에서 가장 중요한 한 단계입니다. 3절에서 8개 유전자를 검정해 3개가 p < 0.05로 나왔습니다. 그런데 검정을 여러 번 하면, 실제로는 차이가 없는데도 우연히 p < 0.05가 나올 확률이 쌓입니다. 유의수준 0.05로 20번 검정하면 평균 1번은 거짓 양성(false positive)이 나오는 셈이죠. RNA-seq에서는 유전자 2만 개를 한꺼번에 검정하니, 보정 없이는 우연한 "유의" 유전자가 수백 개씩 쏟아집니다. 개념은 다중검정보정 편에서 자세히 풀었습니다.

해결은 다중검정 보정(multiple testing correction)입니다. 그중 생물정보에서 표준으로 쓰는 벤야미니-호크베르그(Benjamini-Hochberg, BH) 절차는 거짓발견율(false discovery rate, FDR)을 통제해, 원시 p값을 보정된 q값(adjusted p-value)으로 바꿉니다. scipy 1.11부터 scipy.stats.false_discovery_control이 이걸 한 줄로 해 줍니다.

# 원시 p값 (3절의 t검정 결과) → BH 보정 q값
q = stats.false_discovery_control(p_t, method="bh")   # method='bh' = Benjamini-Hochberg

for g, p, qv in zip(genes, p_t, q):
    flag = "*" if qv < 0.05 else " "
    print(f"{g}: p = {p:.4f} → q = {qv:.4f} {flag}")

print("원시 p < 0.05 개수:", int((p_t < 0.05).sum()))
print("보정 q < 0.05 개수:", int((q   < 0.05).sum()))

GeneA: p = 0.5466 → q = 0.8745
GeneB: p = 0.8914 → q = 0.8914
GeneC: p = 0.0249 → q = 0.0665
GeneD: p = 0.0086 → q = 0.0344 *
GeneE: p = 0.6613 → q = 0.8817
GeneF: p = 0.0020 → q = 0.0156 *
GeneG: p = 0.8680 → q = 0.8914
GeneH: p = 0.4531 → q = 0.8745
원시 p < 0.05 개수: 3
보정 q < 0.05 개수: 2

보정의 효과가 한눈에 드러납니다. 원시 p값으로는 GeneC·GeneD·GeneF 세 개가 유의했지만, 보정 후에는 GeneC가 q = 0.0665로 밀려나 둘만 남습니다. GeneC의 원시 p값 0.0249는 8번 검정 중 하나로는 "그 정도는 우연히도 나올 수 있다"고 판정된 것이죠. 보정 q값은 항상 원시 p값보다 크거나 같고, 검정 개수가 많을수록 더 가혹하게 올라갑니다. RNA-seq 분석 도구(DESeq2·edgeR 등)가 결과 표에 padj를 따로 두는 이유가 바로 이것입니다 — 봐야 할 건 원시 p가 아니라 보정된 q입니다.

함수가 제대로 동작하는지 공식 문서의 예시로 한 번 더 확인할 수 있습니다.

# scipy 공식 문서 예시 — 결과가 문서와 일치
print(stats.false_discovery_control([0.5, 0.6, 0.1, 0.001]))

[0.6   0.6   0.2   0.004]

statsmodels를 이미 쓰고 있다면 from statsmodels.stats.multitest import multipletests 후 multipletests(p_t, method="fdr_bh")로도 같은 BH 보정을 얻습니다. 결과는 scipy와 동일하니, 의존성에 맞춰 고르면 됩니다.

자주 발생하는 에러 & 해결

• 다중검정 미보정 — 가장 흔하고 치명적인 실수 — 유전자 수천 개를 검정하고 원시 p < 0.05를 그대로 "유의"라고 발표하면, 거짓 양성이 대량으로 섞입니다. 1만 개를 유의수준 0.05로 검정하면 진짜 차이가 하나도 없어도 평균 500개가 우연히 통과하죠. 반드시 scipy.stats.false_discovery_control(또는 statsmodels multipletests)로 보정한 q값으로 거르고, 표·그림에도 adjusted p값을 명시합니다.
• 축(axis) 혼동 — expr.mean()처럼 axis를 빼면 전체 평균(스칼라 하나)이 나옵니다. 유전자별이 필요한데 axis=0을 주면 샘플별 결과가 나와 조용히 틀린 분석으로 이어지죠. "axis는 없어지는 축"으로 외우고, 결과 배열의 모양(shape)을 매번 print로 확인하는 습관이 안전합니다. 정규화·요약에서 행과 열을 바꿔 넣는 실수가 의외로 잦습니다.
• 정수 나눗셈 — 정수 배열을 정수로 나누면 의도와 달리 몫만 남을 수 있습니다. np.array([1,2,3,4]) // 4는 [0 0 0 1]이 되죠(//는 정수 나눗셈). 비율·정규화에는 실수 나눗셈 /를 쓰고, 배열을 astype(float)로 미리 실수화하거나 dtype=float로 만들어 둡니다. 카운트 데이터(read count)를 다룰 때 특히 조심합니다.
• NaN 전파 — 결측값(NaN)이 하나라도 섞이면 mean·sum·std의 결과가 통째로 NaN이 됩니다. np.array([1,2,np.nan,4]).mean()은 nan이죠. 결측을 무시하고 계산하려면 np.nanmean·np.nansum·np.nanstd를 쓰거나, np.isnan으로 걸러낸 뒤 계산합니다. 발현 행렬에 결측이 있으면 검정 전에 처리 방침(제거·대치)을 먼저 정해야 합니다.
• 표본수·정규성 가정 위반 — 군당 2~3개처럼 표본이 너무 적으면 t검정의 정규성 가정이 위태롭고 검정력도 약합니다. 표본이 적거나 분포가 치우쳤으면 만-휘트니 U 검정 같은 비모수 검정을 쓰고, 가능하면 반복(replicate)을 늘립니다. 검정을 고르는 기준이 헷갈리면 R 통계검정 편의 결정 흐름(그룹 수·짝지음·정규성)이 그대로 적용됩니다.
• ttest_ind의 등분산 가정 — ttest_ind는 기본적으로 두 군의 분산이 같다고 가정합니다(Student t). 군마다 분산이 크게 다르면 equal_var=False를 줘 Welch t검정으로 바꾸는 게 안전합니다. 실제 발현 데이터는 분산이 들쭉날쭉한 경우가 많아, Welch를 기본으로 두는 분석자도 많습니다.

확장 — 더 해 볼 것

scipy.cluster로 계층적 군집 — 정규화한 발현 행렬을 scipy.cluster.hierarchy의 linkage·dendrogram에 넣으면 유전자·샘플을 유사도로 묶어 덴드로그램(dendrogram)을 그릴 수 있습니다. 히트맵의 행·열 정렬이 바로 이 군집 결과죠. 거리 계산은 scipy.spatial.distance가 맡습니다.

from scipy.cluster.hierarchy import linkage, fcluster

Z = linkage(z, method="ward")        # z-점수 행렬로 유전자 계층 군집
clusters = fcluster(Z, t=3, criterion="maxclust")   # 3개 군집으로 자르기
print("유전자별 군집 라벨:", clusters)

• scipy.optimize로 곡선 적합 — 용량-반응(dose-response)이나 성장 곡선처럼 모델 함수에 데이터를 맞추는 일은 scipy.optimize.curve_fit이 해 줍니다. 시그모이드·지수 함수를 정의하고 데이터를 넘기면 최적 파라미터와 공분산을 돌려주죠. 효소 반응속도(Michaelis-Menten) 적합에도 같은 방식을 씁니다.
• 검정 선택을 함수로 — 그룹 수·짝지음 여부·정규성에 따라 검정을 자동으로 고르는 래퍼(wrapper)를 짜 두면 분석이 빨라집니다. 정규성은 scipy.stats.shapiro, 등분산은 scipy.stats.levene으로 미리 점검할 수 있습니다.
• 벡터화로 속도 끌어올리기 — 유전자별 검정을 리스트 내포로 돌렸지만, 행이 수만 개면 scipy.stats.ttest_ind(ctrl, treat, axis=1)처럼 axis 인자로 전체를 한 번에 벡터화하면 훨씬 빠릅니다. numpy·scipy의 벡터화 연산은 파이썬 반복문보다 수십~수백 배 빠르죠.
• 분류·머신러닝으로 잇기 — 여기서 만든 expr 배열(또는 z-점수 z)을 그대로 특징 행렬 X로 삼아 scikit-learn 편의 분류·교차검증으로 넘기면, "어떤 유전자 조합이 군을 가르는가"까지 이어집니다.

References

1. Harris, C. R. et al. (2020). Array programming with NumPy. Nature, 585, 357–362. (numpy 원논문)
2. Virtanen, P. et al. (2020). SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17, 261–272. (scipy 원논문)
3. NumPy developers. (2026). NumPy v2.5 documentation — Broadcasting · Indexing. numpy.org/broadcasting
4. SciPy developers. (2026). scipy.stats.false_discovery_control — Benjamini-Hochberg·Benjamini-Yekutieli FDR 보정(scipy 1.11+). false_discovery_control docs
5. SciPy developers. (2026). scipy.stats.ttest_ind · mannwhitneyu — 독립표본 검정 API. ttest_ind docs
6. Benjamini, Y. & Hochberg, Y. (1995). Controlling the False Discovery Rate. Journal of the Royal Statistical Society B, 57, 289–300. (BH 절차 원논문)

Pipette & Pipeline · A bio portfolio journal

TL;DR — 유전자 발현이나 임상 특징으로 "종양이 양성인가 악성인가"를 맞히는 이진 분류(binary classification)를, scikit-learn (파이썬 머신러닝 표준 라이브러리)으로 데이터 적재부터 ROC·AUC 평가까지 한 번에 끝냅니다. 재현 가능한 sklearn 내장 데이터셋 load_breast_cancer (569 샘플 × 30 특징)를 쓰고, 로지스틱 회귀(logistic regression)와 랜덤 포레스트(random forest) 두 모델을 비교합니다.

머신러닝 코드가 길어 보여도 뼈대는 단순합니다. 데이터를 학습용·시험용으로 나누고(train/test split), 스케일을 맞추고(scaling), 모델을 학습시키고(fit), 교차검증(cross-validation)으로 성능을 가늠한 뒤, 혼동행렬(confusion matrix)과 ROC 곡선으로 평가합니다. 이 흐름만 손에 익으면 데이터만 바꿔 어떤 분류 문제에도 그대로 씁니다.

이 글의 코드는 무작위 시드를 고정해 그대로 복사하면 같은 숫자가 나옵니다. 초보가 가장 많이 밟는 지뢰 — 데이터 누수(data leakage)·클래스 불균형(class imbalance)·과적합(overfitting)도 코드로 짚고, 마지막엔 Pipeline (파이프라인)으로 누수를 원천 차단하는 법까지 갑니다.

🔗 관련 글  ·  모델에 넣을 데이터를 먼저 다듬으려면 pandas로 생물 데이터 정제하기에서 이어집니다  ·  학습 결과를 그림으로 그리는 도구는 matplotlib·seaborn으로 생물 데이터 시각화에서 다룹니다  ·  같은 ROC·AUC를 R로 그리는 법은 pROC로 ROC 곡선·AUC 그리기와 1:1로 비교됩니다  ·  특징을 줄여 구조를 보는 차원축소는 주성분분석(PCA)으로 RNA-seq 샘플 QC하기로 이어집니다

이 글에서 만들 것

발현 데이터나 임상 수치를 받았을 때 가장 흔한 질문이 "이 샘플이 어느 그룹에 속하나"입니다. 종양이 양성인지 악성인지, 환자가 약에 반응할지 안 할지, 세포가 어떤 아형(subtype)인지 — 모두 분류 문제입니다. scikit-learn은 이런 분류를 몇 줄로 풀게 해 주는 파이썬 표준 도구입니다.

scikit-learn (David Cournapeau, 2007~)은 분류·회귀·군집·차원축소를 한 일관된 문법으로 묶은 라이브러리입니다. 어떤 모델이든 fit(학습)·predict(예측)·score(평가)라는 같은 인터페이스를 쓰기 때문에, 모델을 바꿔도 코드 골격이 거의 그대로입니다. 딥러닝이 필요 없는 표 형식(tabular) 데이터에서는 여전히 가장 먼저 잡는 도구죠.

이 글은 sklearn에 내장된 유방암 진단 데이터(load_breast_cancer)로 이진 분류를 처음부터 끝까지 만듭니다. ① 데이터를 적재하고, ② 학습용·시험용으로 나누고(train_test_split), ③ 스케일을 맞추고(StandardScaler), ④ 로지스틱 회귀와 랜덤 포레스트를 학습시키고, ⑤ 교차검증(cross_val_score·StratifiedKFold)으로 성능을 가늠한 뒤, ⑥ 혼동행렬·classification_report·ROC·AUC로 평가하고 특징 중요도(feature importance)까지 봅니다. 코드는 scikit-learn 1.9 기준이고, 시드를 고정해 재현됩니다.

# 핵심 흐름 한눈에
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score, RocCurveDisplay

X, y = load_breast_cancer(return_X_y=True)                       # ① 데이터 적재
Xtr, Xte, ytr, yte = train_test_split(X, y, stratify=y, random_state=42)  # ② 분할
Xtr_s = StandardScaler().fit_transform(Xtr)                      # ③ 스케일링(train fit)
clf = LogisticRegression(max_iter=5000).fit(Xtr_s, ytr)         # ④ 학습
print(cross_val_score(clf, Xtr_s, ytr, cv=5, scoring="roc_auc"))  # ⑤ 교차검증
RocCurveDisplay.from_estimator(clf, Xte_s, yte)                  # ⑥ ROC·AUC

1. 사전 준비 — 설치와 데이터 로드

scikit-learn 하나만 설치하면 됩니다. numpy·scipy는 의존성으로 따라오고, 그림을 그리려면 matplotlib도 함께 깔아 둡니다.

# pip로 설치 (가상환경 권장)
pip install scikit-learn matplotlib

# 또는 conda
conda install -c conda-forge scikit-learn matplotlib

import sklearn
print("scikit-learn", sklearn.__version__)

scikit-learn 1.9.0

이제 데이터를 불러옵니다. load_breast_cancer는 위스콘신 유방암 진단 데이터로, 종양의 디지털 이미지에서 뽑은 30개 수치 특징(반지름·질감·둘레·면적·매끄러움 등)과, 각 종양이 양성(benign)인지 악성(malignant)인지를 담고 있습니다. 569개 샘플 전부가 라벨이 붙어 있어 분류 연습에 딱 맞습니다.

import numpy as np
import pandas as pd
from sklearn.datasets import load_breast_cancer

data = load_breast_cancer(as_frame=True)     # as_frame=True → pandas로 받기
X = data.data                                 # 특징 30개 (569 × 30)
y = data.target                               # 라벨 (0/1)

print(X.shape)                                # (샘플 수, 특징 수)
print(data.target_names)                      # 클래스 이름
print(y.value_counts())                       # 클래스별 개수

(569, 30)
['malignant' 'benign']
target
1    357
0    212
Name: count, dtype: int64

여기서 꼭 짚어야 할 함정이 하나 있습니다. 클래스 라벨이 직관과 반대로 매겨져 있습니다 — 0이 악성(malignant), 1이 양성(benign)입니다. 의학적으로 "잡아야 하는 양성(positive) 사건"은 악성인데, 숫자로는 0이라 나중에 ROC·재현율(recall)을 해석할 때 헷갈리기 쉽습니다. 데이터를 받으면 라벨이 무엇을 가리키는지부터 확인하는 습관이 중요합니다.

특징이 30개라 어떤 값인지 잠깐 살펴봅니다. 평균(mean)·표준오차(error)·최악값(worst) 세 묶음으로 같은 측정치가 반복됩니다.

X.iloc[:3, :5]                                # 앞 3 샘플, 앞 5 특징만

   mean radius  mean texture  mean perimeter  mean area  mean smoothness
0        17.99         10.38          122.80     1001.0          0.11840
1        20.57         17.77          132.90     1326.0          0.08474
2        19.69         21.25          130.00     1203.0          0.10960

값의 크기(scale)가 제각각인 게 보입니다. mean area는 1000 단위인데 mean smoothness는 0.1 단위죠. 이렇게 스케일이 다른 특징을 그대로 거리 기반·기울기 기반 모델에 넣으면 큰 값이 학습을 지배합니다. 그래서 3절에서 스케일을 맞춥니다.

2. 데이터 분할 — train_test_split으로 시험셋 떼어 두기

머신러닝의 첫 원칙은 "모델이 본 적 없는 데이터로 평가한다"입니다. 학습에 쓴 데이터로 성능을 재면 외워서 잘 맞히는 것인지, 정말 일반화(generalization)된 것인지 구분이 안 됩니다. 그래서 데이터를 학습용(train)과 시험용(test)으로 먼저 갈라 두고, 시험셋은 마지막 평가 전까지 건드리지 않습니다.

train_test_split이 이 분할을 한 줄로 해 줍니다. 두 옵션이 핵심입니다 — stratify로 클래스 비율을 학습·시험 양쪽에서 똑같이 유지하고, random_state로 분할을 고정해 재현합니다.

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(
    X, y,
    test_size=0.25,        # 시험셋 25% (학습 75%)
    stratify=y,            # 클래스 비율 유지 (불균형 데이터 필수)
    random_state=42,       # 분할 고정 → 재현성
)
print(X_train.shape, X_test.shape)
print(y_train.value_counts(normalize=True).round(3))   # 학습셋 클래스 비율
print(y_test.value_counts(normalize=True).round(3))    # 시험셋 클래스 비율

(426, 30) (143, 30)
target
1    0.627
0    0.373
Name: proportion, dtype: float64
target
1    0.629
0    0.371
Name: proportion, dtype: float64

stratify=y를 준 덕분에 학습셋과 시험셋의 클래스 비율(약 63 대 37)이 거의 같습니다. 이걸 빼면 무작위로 나뉘면서 한쪽에 악성이 몰릴 수 있고, 특히 클래스가 불균형하거나 샘플이 적을 때 평가가 들쭉날쭉해집니다. 분류 문제에서는 거의 항상 stratify를 켜는 게 안전합니다.

3. 스케일링 — StandardScaler와 데이터 누수 막기

1절에서 봤듯 특징마다 값의 크기가 다릅니다. StandardScaler는 각 특징을 평균 0·표준편차 1로 표준화(standardization)해, 모든 특징을 같은 출발선에 세웁니다. 로지스틱 회귀처럼 기울기로 학습하는 모델은 스케일링 여부에 따라 수렴 속도와 성능이 꽤 달라집니다.

여기서 머신러닝 입문자가 가장 많이 밟는 지뢰가 나옵니다 — 스케일러는 학습셋으로만 fit해야 합니다. 시험셋까지 합쳐 평균·표준편차를 구하면, 시험셋의 정보가 학습 과정에 새어 들어가는 데이터 누수(data leakage)가 됩니다. 평가가 부풀려져 실제 성능을 과대평가하게 되죠.

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
X_train_s = scaler.fit_transform(X_train)   # 학습셋: fit(평균·표준편차 학습) + 변환
X_test_s  = scaler.transform(X_test)        # 시험셋: 변환만 (학습셋 통계로)

print(X_train_s.mean(axis=0).round(2)[:5])  # 학습셋 각 특징 평균 ≈ 0
print(X_train_s.std(axis=0).round(2)[:5])   # 학습셋 각 특징 표준편차 ≈ 1

[ 0.  0. -0.  0.  0.]
[1. 1. 1. 1. 1.]

핵심은 fit_transform과 transform의 구분입니다. 학습셋에는 fit_transform(통계를 배우고 동시에 변환), 시험셋에는 transform(학습셋이 배운 통계로 변환만)을 씁니다. 시험셋의 평균이 정확히 0이 아닌 건 정상입니다 — 학습셋 기준으로 변환했으니까요. 이 분리를 손으로 하면 실수하기 쉬워서, 5절에서 Pipeline으로 자동화합니다.

4. 모델 학습 — 로지스틱 회귀와 랜덤 포레스트

이제 모델을 학습시킵니다. 성격이 다른 두 분류기를 비교합니다. 로지스틱 회귀(LogisticRegression)는 특징의 가중합으로 확률을 내는 선형 모델로, 빠르고 해석이 쉽습니다. 랜덤 포레스트(RandomForestClassifier)는 결정 트리(decision tree) 수백 개를 모아 투표시키는 앙상블(ensemble) 모델로, 비선형 관계와 특징 간 상호작용을 잘 잡습니다.

scikit-learn의 일관된 문법 덕에, 두 모델 모두 fit으로 학습하고 predict로 예측합니다. 모델 객체만 바꾸면 나머지 코드는 그대로죠.

from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier

# ① 로지스틱 회귀 (선형·빠름·해석 쉬움)
logit = LogisticRegression(max_iter=5000, random_state=42)
logit.fit(X_train_s, y_train)               # 스케일링한 학습셋으로 학습

# ② 랜덤 포레스트 (비선형·앙상블·트리 기반)
rf = RandomForestClassifier(n_estimators=300, random_state=42)
rf.fit(X_train, y_train)                     # 트리 모델은 스케일링 불필요

# 학습셋 정확도(accuracy)로 빠르게 확인
print("로지스틱 회귀 train acc:", round(logit.score(X_train_s, y_train), 3))
print("랜덤 포레스트 train acc:", round(rf.score(X_train, y_train), 3))

로지스틱 회귀 train acc: 0.988
랜덤 포레스트 train acc: 1.0

max_iter=5000은 로지스틱 회귀가 충분히 수렴하도록 반복 횟수를 늘린 것입니다(기본값 100이면 수렴 경고가 뜰 수 있습니다). 랜덤 포레스트는 트리 기반이라 특징 스케일에 둔감해서 스케일링하지 않은 원본 X_train을 그대로 넣었습니다 — 트리는 "이 값보다 큰가/작은가"로 나누기 때문에 크기 자체는 상관없죠.

여기서 랜덤 포레스트의 학습셋 정확도가 1.0인 게 눈에 걸립니다. 학습 데이터를 완벽히 맞혔다는 건데, 이게 곧 좋은 모델이라는 뜻은 아닙니다. 학습셋만 외운 과적합(overfitting)일 수 있어서, 반드시 본 적 없는 데이터로 다시 확인해야 합니다. 그 일을 교차검증이 합니다.

5. 교차검증 — cross_val_score로 성능을 정직하게

시험셋 한 번의 점수는 우연에 휘둘립니다. 어쩌다 쉬운 샘플이 시험셋에 몰리면 점수가 높게 나오죠. 교차검증(cross-validation)은 학습 데이터를 여러 조각(fold)으로 나눠, 번갈아 가며 한 조각으로 평가하고 나머지로 학습하기를 반복해, 성능의 평균과 흔들림(표준편차)을 함께 봅니다.

분류에서는 각 조각의 클래스 비율을 유지하는 StratifiedKFold를 씁니다. cross_val_score에 cv=5만 줘도 분류 문제라면 내부적으로 계층 분할을 쓰지만, 시드를 고정하려면 객체로 명시하는 게 좋습니다.

from sklearn.model_selection import cross_val_score, StratifiedKFold

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)  # 5겹·클래스 비율 유지

# 평가 지표로 ROC AUC 사용 (불균형에 강건)
auc_logit = cross_val_score(logit, X_train_s, y_train, cv=cv, scoring="roc_auc")
auc_rf    = cross_val_score(rf,    X_train,   y_train, cv=cv, scoring="roc_auc")

print("로지스틱 회귀 AUC:", auc_logit.round(3), "평균", round(auc_logit.mean(), 3))
print("랜덤 포레스트 AUC:", auc_rf.round(3),    "평균", round(auc_rf.mean(), 3))

로지스틱 회귀 AUC: [0.995 0.995 0.999 0.996 0.991] 평균 0.995
랜덤 포레스트 AUC: [0.987 0.995 0.997 0.989 0.985] 평균 0.991

다섯 조각의 AUC가 모두 0.98 이상이고 흔들림도 작습니다. 학습셋 정확도 1.0이던 랜덤 포레스트도 교차검증 AUC는 0.991로, 외운 게 아니라 실제로 일반화됐음이 확인됩니다. 이 데이터에서는 단순한 로지스틱 회귀(평균 0.995)가 랜덤 포레스트보다 오히려 살짝 앞섭니다. 특징이 잘 정제돼 있고 클래스가 선형으로 꽤 나뉘기 때문인데, 늘 복잡한 모델이 이기는 건 아니라는 걸 보여 주죠.

scoring에 "roc_auc"를 쓴 이유가 있습니다. 정확도(accuracy)는 클래스가 불균형하면 속기 쉽습니다 — 95%가 음성인 데이터라면 전부 음성이라 찍어도 정확도 95%가 나오죠. AUC는 양성·음성을 가르는 능력 자체를 보기 때문에 불균형에 더 강건합니다.

6. 평가 — 혼동행렬·classification_report·ROC·AUC

교차검증으로 모델을 골랐으니, 이제 떼어 둔 시험셋으로 최종 평가합니다. 분류 평가는 정확도 한 숫자로 끝내지 않습니다. 어디서 틀렸는지(혼동행렬), 정밀도·재현율은 어떤지(classification_report), 임계값(threshold)을 바꾸면 어떻게 되는지(ROC 곡선)를 함께 봅니다.

먼저 혼동행렬(confusion matrix)입니다. 예측과 실제를 2×2 표로 교차해, 참양성(TP)·거짓양성(FP)·거짓음성(FN)·참음성(TN)이 각각 몇 개인지 보여 줍니다.

from sklearn.metrics import confusion_matrix, classification_report

y_pred = logit.predict(X_test_s)             # 시험셋 예측

cm = confusion_matrix(y_test, y_pred)        # 행=실제, 열=예측
print(cm)
print(classification_report(y_test, y_pred,
                            target_names=data.target_names))

[[ 51   2]
 [  1  89]]
              precision    recall  f1-score   support

   malignant       0.98      0.96      0.97        53
      benign       0.98      0.99      0.98        90

    accuracy                           0.98       143
   macro avg       0.98      0.98      0.98       143
weighted avg       0.98      0.98      0.98       143

혼동행렬을 읽으면, 악성 53개 중 51개를 맞히고 2개를 양성으로 잘못 봤으며(거짓음성), 양성 90개 중 89개를 맞혔습니다. classification_report는 여기서 정밀도(precision)·재현율(recall)·F1 점수를 클래스별로 풀어 줍니다. 의학에서는 악성을 놓치는 일(거짓음성)이 특히 위험하므로 악성의 재현율 0.96을 눈여겨봅니다. 단일 숫자인 정확도(0.98)만 보면 이런 비대칭을 놓치죠.

이제 ROC 곡선과 AUC입니다. ROC (Receiver Operating Characteristic)는 분류 임계값을 0에서 1까지 움직이며 거짓양성률(가로)과 참양성률(세로)을 그린 곡선이고, 그 아래 면적이 AUC (Area Under the Curve)입니다. AUC가 1에 가까울수록 완벽한 분류, 0.5면 동전 던지기 수준입니다.

import matplotlib.pyplot as plt
from sklearn.metrics import roc_auc_score, RocCurveDisplay

# 양성(class 1) 확률 예측 → AUC 계산
y_score = logit.predict_proba(X_test_s)[:, 1]
print("시험셋 AUC:", round(roc_auc_score(y_test, y_score), 3))

fig, ax = plt.subplots(figsize=(6, 6))
RocCurveDisplay.from_estimator(logit, X_test_s, y_test,
                               name="로지스틱 회귀", ax=ax)
RocCurveDisplay.from_estimator(rf, X_test, y_test,
                               name="랜덤 포레스트", ax=ax,
                               plot_chance_level=True)   # 대각선 기준선
ax.set_title("ROC 곡선 — 시험셋")
plt.tight_layout()
plt.savefig("roc.png", dpi=300, bbox_inches="tight")

시험셋 AUC: 0.995

RocCurveDisplay.from_estimator는 모델·데이터·라벨만 주면 ROC 곡선을 곡선 아래 면적(AUC)과 함께 자동으로 그려 줍니다. 같은 ax에 두 모델을 겹쳐 그리면 성능을 한 화면에서 비교할 수 있죠. plot_chance_level=True가 대각선(AUC=0.5) 기준선을 더해 줍니다. predict_proba(...)[:, 1]로 양성 클래스의 확률을 뽑아 roc_auc_score에 직접 넣을 수도 있는데, 두 방법의 결과는 같습니다.

마지막으로 특징 중요도(feature importance)입니다. 랜덤 포레스트는 각 특징이 분류에 얼마나 기여했는지를 feature_importances_로 알려 줘, "어떤 측정치가 양성·악성을 가르는가"를 해석하게 해 줍니다.

import pandas as pd

# 특징 중요도 상위 8개
imp = (pd.Series(rf.feature_importances_, index=X.columns)
       .sort_values(ascending=False).head(8))
print(imp.round(3))

worst perimeter        0.142
worst concave points   0.139
worst area             0.121
mean concave points    0.106
worst radius           0.082
mean perimeter         0.052
worst texture          0.038
mean area              0.037

종양 둘레·면적·오목한 점(concave points)의 "최악값(worst)" 측정치가 분류를 주도합니다. 임상적으로도 악성 종양이 크고 불규칙한 경계를 갖는다는 사실과 맞아떨어지죠. 다만 트리 기반 중요도는 값의 범위가 넓은 특징을 과대평가하는 편향이 있어, 더 엄밀하게 보려면 순열 중요도(permutation importance)를 함께 쓰는 게 좋습니다.

자주 발생하는 에러 & 해결

• 데이터 누수(data leakage) — 스케일러를 전체에 fit — StandardScaler().fit(X)로 전체 데이터(학습+시험)에 한꺼번에 fit하면, 시험셋의 평균·표준편차가 변환에 스며들어 평가가 부풀려집니다. 교차검증에서 누수가 특히 잘 일어나는데, 매 fold마다 전처리를 다시 fit해야 하기 때문입니다. 해결은 전처리와 모델을 Pipeline으로 묶는 것입니다(아래 확장 참고). 그러면 cross_val_score가 fold마다 스케일러를 학습셋에만 다시 fit해 누수를 원천 차단합니다.
• 클래스 불균형(class imbalance) — 정확도에 속음 — 음성이 95%인 데이터에서 전부 음성으로 찍어도 정확도 95%가 나옵니다. 정확도만 보면 쓸모없는 모델을 좋다고 착각하죠. 평가 지표를 ROC AUC·F1·재현율로 바꾸고, 모델에 class_weight="balanced"를 주거나 소수 클래스를 늘리는 SMOTE (imbalanced-learn 패키지) 같은 기법을 씁니다. 분할·교차검증에는 반드시 stratify·StratifiedKFold로 비율을 유지합니다.
• 과적합(overfitting) — 학습셋만 잘 맞힘 — 학습 정확도는 1.0인데 시험 성능이 뚝 떨어지면 과적합입니다. 모델이 데이터를 일반화한 게 아니라 외운 것이죠. 랜덤 포레스트는 max_depth로 트리 깊이를 제한하거나 min_samples_leaf를 키워 규제(regularization)하고, 로지스틱 회귀는 C를 줄여 규제를 강화합니다. 무엇보다 학습 점수가 아니라 교차검증·시험셋 점수로 판단하는 게 핵심입니다.
• 수렴 경고(ConvergenceWarning) — LogisticRegression이 "lbfgs failed to converge"를 띄우면, 정해진 반복 안에 최적해를 못 찾았다는 뜻입니다. max_iter를 늘리거나(예: 5000), 특징을 StandardScaler로 스케일링하면 대개 사라집니다. 스케일이 제각각이면 최적화가 느려지기 때문에, 스케일링은 수렴 문제의 해결책이기도 합니다.
• predict_proba vs decision_function — ROC·AUC에는 확률 또는 점수가 필요한데, roc_auc_score(y_test, y_pred)처럼 0/1 예측 라벨을 넣으면 곡선이 계단처럼 거칠어지고 AUC가 과소평가됩니다. 반드시 predict_proba(X)[:, 1](양성 확률)이나 decision_function(X)(결정 점수)를 넣어야 합니다. RocCurveDisplay.from_estimator는 이걸 알아서 처리해 주니 가장 안전합니다.
• 시험셋을 모델 선택에 씀 — 시험셋을 보면서 하이퍼파라미터를 고르면, 시험셋이 사실상 학습에 쓰인 셈이라 성능이 부풀려집니다. 시험셋은 마지막 단 한 번의 평가에만 씁니다. 모델·하이퍼파라미터 선택은 학습셋 안의 교차검증이나 별도 검증셋(validation set)으로 해야 합니다.

확장 — 더 해 볼 것

Pipeline으로 누수 차단 — 전처리와 모델을 하나로 묶으면, fit·교차검증 어디서든 스케일러가 학습셋에만 fit돼 데이터 누수가 구조적으로 사라집니다. 코드도 짧아지죠.

from sklearn.pipeline import make_pipeline

pipe = make_pipeline(StandardScaler(), LogisticRegression(max_iter=5000))
auc = cross_val_score(pipe, X, y, cv=cv, scoring="roc_auc")   # 누수 없는 교차검증
print("Pipeline AUC 평균:", round(auc.mean(), 3))

GridSearchCV로 하이퍼파라미터 탐색 — 모델 설정(예: 랜덤 포레스트의 트리 깊이, 로지스틱 회귀의 규제 강도 C)을 격자(grid)로 훑어 교차검증 점수가 가장 높은 조합을 자동으로 고릅니다. Pipeline과 함께 쓰면 전처리 옵션까지 한꺼번에 탐색할 수 있습니다.

from sklearn.model_selection import GridSearchCV

grid = {"randomforestclassifier__max_depth": [4, 8, None],
        "randomforestclassifier__min_samples_leaf": [1, 3, 5]}
search = GridSearchCV(make_pipeline(RandomForestClassifier(random_state=42)),
                      grid, cv=cv, scoring="roc_auc")
search.fit(X, y)
print("최적 파라미터:", search.best_params_, "| AUC:", round(search.best_score_, 3))

• 순열 중요도(permutation importance) — 트리 기반 feature_importances_의 편향을 피하려면, 특징 값을 무작위로 섞었을 때 성능이 얼마나 떨어지는지로 중요도를 재는 permutation_importance가 더 공정합니다. 모델 종류와 무관하게 쓸 수 있습니다.
• 다중 분류(multiclass)·실제 발현 데이터 — 이 글은 이진 분류였지만, 세포 아형 3종 이상을 가르는 다중 분류도 같은 코드로 됩니다(roc_auc_score에 multi_class="ovr"). 실제 RNA-seq 발현 행렬을 쓰려면 pandas 정제로 표를 다듬어 X에 넣으면 그대로 이어집니다.
• 모델 저장·배포 — 학습한 모델을 joblib.dump로 파일에 저장해 두면, 나중에 불러와 새 샘플을 예측할 수 있습니다. 다만 sklearn 버전이 바뀌면 호환이 깨질 수 있으니 버전을 함께 기록하는 게 안전합니다.

References

1. Pedregosa, F. et al. (2011). Scikit-learn: Machine Learning in Python. Journal of Machine Learning Research, 12, 2825–2830. (scikit-learn 원논문)
2. scikit-learn developers. (2026). scikit-learn 1.9 documentation — User Guide. scikit-learn.org/user_guide
3. scikit-learn. Common pitfalls and recommended practices — 데이터 누수·교차검증 함정. common_pitfalls
4. scikit-learn. load_breast_cancer — 위스콘신 유방암 진단 데이터셋(569 × 30). load_breast_cancer docs
5. scikit-learn. RocCurveDisplay — ROC 곡선·AUC 시각화 API. RocCurveDisplay docs

Pipette & Pipeline · A bio portfolio journal

TL;DR — 차등발현(differential expression) 결과나 발현 행렬을 받았을 때 가장 많이 그리는 그림 네 가지를 matplotlib (파이썬 저수준 플로팅 라이브러리)와 seaborn (그 위에 얹은 고수준 통계 시각화 라이브러리)로 처음부터 끝까지 그립니다. 볼케이노 플롯(volcano plot), 클러스터 히트맵(clustermap), 박스·바이올린 플롯(box/violin plot), 그리고 facet 다중 패널까지입니다.

역할 분담이 분명합니다. matplotlib은 점·선·축을 직접 그리는 바탕이고, seaborn은 "통계 그래프 한 장"을 짧은 명령으로 그려 줍니다. 볼케이노처럼 세밀한 라벨링이 필요한 그림은 matplotlib으로 손보고, 분포 비교나 히트맵은 seaborn 한 줄로 끝냅니다.

이 글의 코드는 모두 무작위 시드를 고정한 작은 발현 데이터로 돌아가, 그대로 복사하면 같은 그림이 나옵니다. seaborn 0.13에서 바뀐 hue·팔레트 규칙처럼 자주 막히는 지점도 함께 짚습니다.

🔗 관련 글  ·  그릴 데이터를 먼저 정제하려면 pandas로 발현 행렬 정제하기에서 이어집니다  ·  같은 그래프를 R로 그리는 법은 ggplot2로 바이오 데이터 그래프 그리기와 1:1로 비교됩니다  ·  단일세포 데이터의 UMAP·클러스터 시각화는 scanpy로 단일세포 RNA-seq 분석하기에서 다룹니다

이 글에서 만들 것

발현 분석을 끝내고 나면 결국 그림으로 보여 줘야 합니다. "어떤 유전자가 처리군에서 확 올라갔나"는 볼케이노 플롯으로, "샘플들이 발현 패턴으로 어떻게 묶이나"는 클러스터 히트맵으로, "관심 유전자가 그룹마다 어떻게 분포하나"는 박스·바이올린 플롯으로 봅니다. 이 세 장에 다중 패널(facet) 한 장을 더하면 RNA-seq 결과 슬라이드의 대부분이 채워집니다.

matplotlib (John Hunter, 2003~)은 파이썬 그래프의 바탕입니다. 점 하나, 선 하나, 글자 하나를 좌표에 직접 찍는 저수준 도구라 못 그리는 게 거의 없지만, 그만큼 손이 많이 갑니다. seaborn (Michael Waskom, 2012~)은 그 위에 얹혀, 데이터프레임과 그릴 열 이름만 주면 통계 그래프 한 장을 만들어 줍니다. 둘은 경쟁 관계가 아니라 층위가 다릅니다 — seaborn으로 큰 그림을 그리고, matplotlib으로 축·라벨·주석을 마무리하는 식으로 같이 씁니다.

이 글은 작은 차등발현 결과와 발현 행렬을 직접 만들어, 네 가지 그림을 순서대로 그립니다. ① 볼케이노 플롯(산점도와 유의 영역 색칠, 상위 유전자 라벨), ② 클러스터 히트맵(z-score 정규화와 계층 군집), ③ 박스·바이올린 플롯(그룹별 발현 분포), ④ facet 다중 패널(catplot으로 여러 유전자를 한 번에)입니다. 마지막에 한글 폰트·저장 해상도(dpi)·SVG 같은 실무 설정을 정리합니다. 코드는 matplotlib 3.11·seaborn 0.13.2 기준이고, 시드를 고정해 재현됩니다.

# 핵심 흐름 한눈에
import matplotlib.pyplot as plt
import seaborn as sns

sns.set_theme(style="whitegrid")                # ① 전역 스타일 한 번에
sns.scatterplot(data=deg, x="log2FC", y="neglog10p", hue="sig")   # ② 볼케이노
sns.clustermap(expr_z, cmap="vlag", center=0)   # ③ z-score 히트맵 + 군집
sns.violinplot(data=long, x="group", y="expr")  # ④ 그룹별 분포
sns.catplot(data=long, x="group", y="expr", col="gene", kind="box")  # ⑤ facet
plt.savefig("figure.png", dpi=300, bbox_inches="tight")            # ⑥ 저장

1. 사전 준비 — 설치와 샘플 데이터

matplotlib·seaborn·pandas를 한 번에 설치합니다. seaborn을 깔면 matplotlib·numpy·pandas는 의존성으로 따라오지만, 명시해 두는 편이 분명합니다. 클러스터 히트맵의 계층 군집은 SciPy (사이파이)가 필요한데, 이 역시 보통 같이 깔려 있습니다.

# pip로 설치 (가상환경 권장)
pip install matplotlib seaborn pandas

# 또는 conda
conda install -c conda-forge matplotlib seaborn pandas

import matplotlib, seaborn as sns, pandas as pd, numpy as np

print("matplotlib", matplotlib.__version__)
print("seaborn", sns.__version__)

matplotlib 3.11.0
seaborn 0.13.2

이제 설명용으로 작은 데이터 두 개를 만듭니다. 하나는 차등발현 결과 표(유전자별 log2 배수변화와 p-value), 다른 하나는 발현 행렬(유전자 × 샘플)입니다. 실제 분석이라면 앞쪽은 DESeq2·edgeR의 결과 테이블이고, 뒤쪽은 정규화한 카운트 행렬이 됩니다.

import numpy as np
import pandas as pd

rng = np.random.default_rng(7)              # 재현성 고정 (시드 7)
n = 300                                      # 유전자 300개

# ── 차등발현 결과 표: 대부분은 변화 없고, 일부만 크게 오르내림 ──
log2fc = rng.normal(0, 1.0, n)              # log2 배수변화 (대칭 분포)
log2fc[:15] += rng.normal(3, 0.4, 15)       # 앞 15개 = 강하게 상향
log2fc[15:30] -= rng.normal(3, 0.4, 15)     # 다음 15개 = 강하게 하향
pval = 10 ** (-np.abs(log2fc) * rng.uniform(0.8, 1.6, n))  # |FC| 클수록 작은 p
pval = np.clip(pval, 1e-12, 1.0)            # p-value 범위 보정

deg = pd.DataFrame({
    "gene":   [f"G{i:03d}" for i in range(n)],
    "log2FC": log2fc,
    "pvalue": pval,
})
deg["neglog10p"] = -np.log10(deg["pvalue"])  # y축에 쓸 -log10(p)
deg.head(3)

   gene    log2FC        pvalue  neglog10p
0  G000  3.412847  1.073229e-04   3.969310
1  G001  3.001210  3.018501e-04   3.520282
2  G002  3.880644  9.142122e-06   5.039459

차등발현 그림에서 "유의(significant)" 여부는 보통 두 기준을 동시에 봅니다 — 배수변화가 충분히 크고(|log2FC| ≥ 1), p-value가 충분히 작은(p < 0.05) 유전자입니다. 이 둘을 합쳐 색칠용 라벨 열을 미리 만들어 둡니다.

# 유의 기준: |log2FC| >= 1 그리고 p < 0.05 → 상향/하향/비유의 3분류
fc_thr, p_thr = 1.0, 0.05
cond_up   = (deg["log2FC"] >=  fc_thr) & (deg["pvalue"] < p_thr)
cond_down = (deg["log2FC"] <= -fc_thr) & (deg["pvalue"] < p_thr)

deg["sig"] = np.select([cond_up, cond_down],
                        ["up", "down"], default="ns")  # ns = not significant
deg["sig"].value_counts()

sig
ns      268
up       17
down     15
Name: count, dtype: int64

다음은 발현 행렬입니다. 히트맵·분포 그림에 쓸 작은 행렬을 따로 만듭니다. 유전자 12개 × 샘플 10개(대조군 5와 처리군 5), 처리군에서 일부 유전자가 올라가도록 신호를 넣습니다.

# ── 발현 행렬: 유전자 12 × 샘플 10 (control 5 + treated 5) ──
genes   = [f"GENE{i:02d}" for i in range(1, 13)]
samples = [f"C{i}" for i in range(1, 6)] + [f"T{i}" for i in range(1, 6)]

base = rng.normal(8, 1.2, size=(12, 10))    # 로그 발현값 바탕 (log2 스케일)
base[:4, 5:] += rng.normal(2.5, 0.4, (4, 5))   # 앞 4개 유전자 = 처리군서 상향
base[4:7, 5:] -= rng.normal(2.0, 0.4, (3, 5))  # 다음 3개 = 처리군서 하향

expr = pd.DataFrame(base, index=genes, columns=samples).round(2)
expr.iloc[:4, :6]                            # 일부만 미리보기

          C1    C2    C3    C4    C5    T1
GENE01  7.42  6.88  9.10  8.05  7.61  10.33
GENE02  8.21  9.04  7.55  8.66  8.12  10.78
GENE03  6.95  7.80  8.43  7.21  8.90   9.52
GENE04  9.13  8.27  7.92  8.48  7.74  11.06

2. 볼케이노 플롯 — 산점도로 차등발현 한눈에

볼케이노 플롯은 x축에 log2 배수변화, y축에 −log₁₀ p를 찍는 산점도입니다. 화산이 분출하듯 좌우 위쪽으로 점이 솟아오른 모양이라 이런 이름이 붙었습니다. 오른쪽 위는 "크게 올라가고 유의한" 유전자, 왼쪽 위는 "크게 내려가고 유의한" 유전자입니다.

seaborn의 scatterplot에 색 기준(hue)으로 아까 만든 sig 열을 넘기면, 상향·하향·비유의가 자동으로 다른 색이 됩니다.

import matplotlib.pyplot as plt
import seaborn as sns

sns.set_theme(style="whitegrid", font_scale=1.0)   # 전역 테마

fig, ax = plt.subplots(figsize=(7, 6))
palette = {"up": "#E76F51", "down": "#2E7D6B", "ns": "#C9C4BA"}  # 색 직접 지정

sns.scatterplot(data=deg, x="log2FC", y="neglog10p",
                hue="sig", hue_order=["up", "down", "ns"],
                palette=palette, s=28, edgecolor="none", ax=ax)

# 기준선: 세로 |log2FC|=1, 가로 p=0.05 (−log10 0.05 ≈ 1.30)
ax.axvline(-1, color="#999", ls="--", lw=1)
ax.axvline( 1, color="#999", ls="--", lw=1)
ax.axhline(-np.log10(0.05), color="#999", ls="--", lw=1)

ax.set_xlabel("log$_2$ fold change")        # 아래첨자: $_2$
ax.set_ylabel("$-$log$_{10}$ $p$")          # 통계기호 p는 수식($)으로 이탤릭
ax.set_title("Volcano plot — control vs treated", weight="bold")
ax.legend(title="유의성", frameon=True)
plt.tight_layout()
plt.savefig("volcano.png", dpi=300, bbox_inches="tight")

여기서 핵심은 세 가지입니다. hue에 범주 열을 주면 seaborn이 범례까지 자동으로 달아 줍니다. palette를 딕셔너리로 주면 "up은 코랄, down은 그린"처럼 색을 못 박을 수 있습니다(브랜드색 통일에 유용합니다). 축 라벨의 $-$log$_{10}$ $p$처럼 $...$로 감싼 부분은 LaTeX 수식으로 렌더되는데, 통계기호 p는 이렇게 수식 안에 넣어야 이탤릭으로 나옵니다.

이제 상위 유전자에 이름표를 답니다. seaborn에는 점 라벨링 기능이 따로 없어서, 이 부분은 matplotlib의 annotate로 직접 그립니다. −log₁₀ p가 가장 큰 상·하향 유전자 몇 개만 골라 텍스트를 얹습니다.

# 상·하향에서 각각 가장 유의한 유전자 3개씩 라벨링
top = pd.concat([
    deg[deg.sig == "up"].nlargest(3, "neglog10p"),
    deg[deg.sig == "down"].nlargest(3, "neglog10p"),
])
for _, r in top.iterrows():
    ax.annotate(r["gene"], (r["log2FC"], r["neglog10p"]),
                xytext=(5, 4), textcoords="offset points",  # 점에서 살짝 띄움
                fontsize=9, fontweight="bold", color="#333")
fig.savefig("volcano_labeled.png", dpi=300, bbox_inches="tight")

annotate의 xytext=(5, 4)와 textcoords="offset points"는 "점에서 오른쪽 5pt·위 4pt 떨어진 곳에 글자를 놓아라"는 뜻입니다. 라벨이 점과 겹쳐 안 보이는 일을 막는 흔한 방법입니다. 라벨이 많아 서로 겹친다면 adjustText 같은 보조 패키지를 쓰지만, 보통은 상위 몇 개만 다는 걸로 충분합니다.

3. 클러스터 히트맵 — clustermap으로 발현 패턴 군집

히트맵은 발현 행렬을 색의 농담으로 보여 줍니다. 그냥 값을 칠하면 발현이 높은 유전자 행만 진하게 나와 패턴이 안 보이므로, 유전자별로 z-score 정규화(각 행을 평균 0·표준편차 1로)한 뒤 칠하는 게 표준입니다. seaborn의 clustermap은 여기에 더해 비슷한 행·열끼리 계층 군집(hierarchical clustering)으로 묶어 덴드로그램(dendrogram)까지 그려 줍니다.

# 행(유전자) 단위 z-score: (값 − 행평균) / 행표준편차
expr_z = expr.sub(expr.mean(axis=1), axis=0).div(expr.std(axis=1), axis=0)

g = sns.clustermap(
    expr_z, cmap="vlag", center=0,          # vlag = 파랑–흰–빨강 발산 팔레트
    vmin=-2, vmax=2,                         # 색 범위 고정 (이상치에 안 휘둘림)
    col_cluster=True, row_cluster=True,      # 행·열 모두 군집
    figsize=(8, 7), linewidths=.4,
    cbar_kws={"label": "z-score"},
    dendrogram_ratio=0.14,
)
g.ax_heatmap.set_xlabel("샘플")
g.ax_heatmap.set_ylabel("유전자")
g.figure.suptitle("발현 z-score 클러스터 히트맵", y=1.02, fontweight="bold")
g.savefig("clustermap.png", dpi=300, bbox_inches="tight")

center=0은 0을 색의 중앙(흰색)에 두라는 뜻이라, 상향(빨강)·하향(파랑)이 직관적으로 갈립니다. vmin·vmax로 색 범위를 ±2에 고정하면 극단값 한두 개가 전체 색을 잡아먹는 일을 막습니다. clustermap은 일반 heatmap과 달리 Figure 전체를 새로 만들어 반환하므로(g.figure로 접근), plt.subplots로 만든 축에 끼워 넣을 수 없다는 점만 기억하면 됩니다. 군집을 끄고 단순 히트맵만 원하면 sns.heatmap(expr_z, ...)을 씁니다.

덴드로그램을 보면 처리군 샘플(T1~T5)과 대조군(C1~C5)이 좌우로 갈려 묶이는 게 보입니다. 우리가 신호를 넣은 그대로, 발현 패턴만으로 두 조건이 분리된 셈입니다.

4. 박스·바이올린 플롯 — 그룹별 발현 분포

관심 유전자 하나의 발현이 그룹마다 어떻게 다른지는 박스플롯이나 바이올린 플롯으로 봅니다. 이 그림들은 데이터가 long 형식(한 행에 샘플·그룹·발현값 하나씩)일 때 자연스럽습니다. 발현 행렬을 melt로 녹이고, 샘플 이름에서 그룹을 뽑아냅니다.

# wide(유전자×샘플) → long, 그리고 샘플명 첫 글자로 그룹 추출
long = (expr.reset_index(names="gene")
        .melt(id_vars="gene", var_name="sample", value_name="expr"))
long["group"] = long["sample"].str[0].map({"C": "control", "T": "treated"})
long.head(3)

     gene sample   expr    group
0  GENE01     C1   7.42  control
1  GENE02     C1   8.21  control
2  GENE03     C1   6.95  control

먼저 바이올린 플롯입니다. 박스플롯이 사분위수(quartile)만 보여 준다면, 바이올린은 분포의 모양(밀도)까지 보여 줍니다. 상향 유전자 하나(GENE01)를 골라 두 그룹을 비교합니다.

sub = long[long["gene"] == "GENE01"]         # 유전자 하나만

fig, ax = plt.subplots(figsize=(6, 5))
sns.violinplot(data=sub, x="group", y="expr",
               hue="group", palette={"control": "#6BAB97", "treated": "#E76F51"},
               legend=False, inner="box", ax=ax)   # inner="box": 안에 박스도
sns.stripplot(data=sub, x="group", y="expr",        # 실제 점도 겹쳐 찍기
              color="#333", size=5, alpha=.7, ax=ax)
ax.set_xlabel(""); ax.set_ylabel("정규화 발현값 (log$_2$)")
ax.set_title("GENE01 — 그룹별 발현 분포", fontweight="bold")
plt.tight_layout()
plt.savefig("violin.png", dpi=300, bbox_inches="tight")

여기서 seaborn 0.13의 바뀐 규칙이 처음 등장합니다 — x축과 같은 변수를 hue에도 넘기고(hue="group"), legend=False로 중복 범례를 끕니다. 예전(0.12 이전)에는 palette만 줘도 각 박스가 다른 색이었지만, 0.13부터는 hue 없이 팔레트만 주면 경고가 뜹니다(아래 두 번째 에러 참고). stripplot을 겹쳐 실제 데이터 점을 같이 찍으면, 샘플 수가 적을 때(여기선 그룹당 5개) 분포 모양이 과장돼 보이는 걸 정직하게 보완할 수 있습니다.

박스플롯도 호출만 바꾸면 됩니다. 여러 유전자를 한 축에 늘어놓을 때 특히 박스플롯이 깔끔합니다.

# 유전자 4개를 한 박스플롯에 (x=유전자, 색=그룹)
sel = long[long["gene"].isin(["GENE01", "GENE02", "GENE05", "GENE06"])]

fig, ax = plt.subplots(figsize=(8, 5))
sns.boxplot(data=sel, x="gene", y="expr", hue="group",
            palette={"control": "#6BAB97", "treated": "#E76F51"}, ax=ax)
ax.set_xlabel(""); ax.set_ylabel("정규화 발현값 (log$_2$)")
ax.set_title("유전자별·그룹별 발현 비교", fontweight="bold")
ax.legend(title="그룹")
plt.tight_layout()
plt.savefig("boxplot.png", dpi=300, bbox_inches="tight")

hue="group"을 주면 유전자마다 control·treated 박스가 나란히 두 개씩 그려집니다. 상향 유전자(GENE01·GENE02)는 treated 박스가 위로, 하향 유전자(GENE05·GENE06)는 아래로 내려간 게 한눈에 들어옵니다.

5. facet 다중 패널 — catplot으로 여러 유전자 한 번에

유전자가 여러 개면 한 그림에 욱여넣기보다 작은 패널을 격자로 나열하는 게 깔끔합니다. 이걸 facet (또는 small multiples)이라 부릅니다. seaborn의 catplot(범주형)·relplot(관계형)은 col·row에 열 이름만 주면 패널을 자동으로 쪼개 줍니다. catplot은 그림(figure) 수준 함수라 plt.subplots 없이 자체적으로 격자를 만듭니다.

# 유전자 6개를 패널 6개로 (가로 3 × 세로 2), 각 패널은 박스플롯
sel6 = long[long["gene"].isin(
    ["GENE01", "GENE02", "GENE03", "GENE05", "GENE06", "GENE08"])]

g = sns.catplot(
    data=sel6, x="group", y="expr",
    hue="group", palette={"control": "#6BAB97", "treated": "#E76F51"},
    col="gene", col_wrap=3, kind="box",     # col=패널 기준, col_wrap=한 줄에 3개
    height=2.8, aspect=0.9, legend=False,
)
g.set_titles("{col_name}")                   # 각 패널 제목 = 유전자명
g.set_axis_labels("", "발현값 (log$_2$)")
g.figure.suptitle("유전자별 발현 분포 (facet)", y=1.03, fontweight="bold")
g.savefig("facet.png", dpi=300, bbox_inches="tight")

col="gene"이 패널을 유전자별로 쪼개고, col_wrap=3이 세 개마다 줄을 바꿉니다. kind="box"를 "violin"·"bar"·"strip"으로 바꾸면 같은 격자에 다른 그림이 그려집니다. catplot이 반환하는 건 축이 아니라 FacetGrid 객체라, 제목·축 라벨은 g.set_titles·g.set_axis_labels처럼 격자 메서드로 한꺼번에 손봅니다. 패널이 너무 많아지면 height·aspect로 한 패널 크기를 줄여 전체가 화면에 들어오게 조절합니다.

6. 스타일·한글 폰트·저장 — 발표용 마무리

마지막은 실무 설정입니다. 발표 슬라이드나 논문 그림으로 쓰려면 한글이 깨지지 않아야 하고, 해상도와 형식을 챙겨야 합니다.

matplotlib은 기본 폰트에 한글 글리프가 없어 한글 라벨이 네모(□□□)로 깨집니다. rcParams로 한글 폰트를 한 번 지정하면 이후 모든 그림에 적용됩니다. 윈도우는 Malgun Gothic, macOS는 AppleGothic, 리눅스는 나눔고딕 등을 씁니다. 음수 부호가 깨지는 것도 함께 막아 줍니다.

import matplotlib.pyplot as plt

# 운영체제별 한글 폰트 (있는 것 하나 고름)
plt.rcParams["font.family"] = ["Malgun Gothic", "AppleGothic", "NanumGothic"]
plt.rcParams["axes.unicode_minus"] = False   # 음수 부호(−) 깨짐 방지
plt.rcParams["figure.dpi"] = 110             # 화면 미리보기 해상도
plt.rcParams["savefig.bbox"] = "tight"       # 저장 시 여백 자동 정리

seaborn의 전역 테마도 그림 전체 인상을 좌우합니다. set_theme의 style(배경·격자)과 context(글자·선 두께)를 조합합니다. 발표용은 글자를 키운 "talk"가, 논문용은 기본 "notebook"이 무난합니다.

# style: white / whitegrid / dark / darkgrid / ticks
# context: paper / notebook / talk / poster (오른쪽으로 갈수록 요소 큼)
sns.set_theme(style="ticks", context="talk", palette="muted")

저장은 savefig로 합니다. 비트맵인 PNG는 dpi로 해상도를 정하고(인쇄·발표는 300 이상), 벡터인 SVG·PDF는 dpi와 무관하게 무한 확대해도 안 깨져 논문 그림에 적합합니다.

fig.savefig("figure.png", dpi=300, bbox_inches="tight")   # 발표·웹: 고해상도 PNG
fig.savefig("figure.svg", bbox_inches="tight")            # 논문: 벡터 SVG
fig.savefig("figure.pdf", bbox_inches="tight")            # 논문: 벡터 PDF

bbox_inches="tight"는 그림 둘레의 불필요한 여백을 잘라 줍니다. 라벨이 그림 밖으로 잘려 나간다면 이 옵션이 거의 해결합니다. 색약(color vision deficiency)을 배려하려면 palette="colorblind"나 발산 팔레트 "vlag"·"icefire"를 쓰는 게 좋습니다.

자주 발생하는 에러 & 해결

• 한글 라벨이 네모(□□□)로 깨짐 — matplotlib 기본 폰트에 한글 글리프가 없어서입니다. 6절처럼 plt.rcParams["font.family"]에 설치된 한글 폰트를 지정하세요. 폰트 캐시가 옛 정보를 들고 있으면 한 번 지운 뒤(matplotlib.font_manager._load_fontmanager(try_read_cache=False)) 재시작하면 잡힙니다. 음수 부호까지 깨지면 axes.unicode_minus를 False로 둡니다.
• seaborn 0.13: palette without hue 경고 — sns.boxplot(data=df, x="group", y="v", palette="Set2")처럼 hue 없이 팔레트만 주면 FutureWarning이 뜨고 0.14에서 막힙니다. 0.13부터 범주형 함수는 hue가 없으면 단색으로 그리는 게 기본입니다. 해결은 x축 변수를 hue에도 넘기고(hue="group") legend=False로 중복 범례를 끄는 것입니다. 이 글의 모든 박스·바이올린 코드가 이 패턴을 씁니다.
• clustermap을 subplots 축에 못 넣음 — fig, ax = plt.subplots(); sns.clustermap(..., ax=ax)는 에러입니다. clustermap은 자체적으로 Figure 전체(히트맵·덴드로그램·컬러바)를 만들어 ClusterGrid를 반환하기 때문입니다. 기존 축에 넣고 싶으면 군집 없는 sns.heatmap(..., ax=ax)을 쓰고, 군집이 꼭 필요하면 clustermap을 독립 그림으로 두세요. 저장은 g.savefig(...)로 합니다.
• 그림이 안 뜨거나 빈 창만 나옴 — 스크립트(.py)에서는 plt.show()를 마지막에 호출해야 창이 뜹니다. 주피터에서는 셀의 마지막 객체가 자동 표시되지만, 같은 셀에서 그림을 그리고 또 다른 걸 그리면 앞 그림이 먹힐 수 있으니 셀을 나누거나 fig를 명시적으로 다루세요. 저장만 하고 화면 표시가 필요 없으면 matplotlib.use("Agg")로 백엔드를 비대화형으로 바꾸면 됩니다.
• 라벨·제목이 그림 밖으로 잘림 — bbox_inches="tight"를 savefig에 주거나, 그리기 전에 plt.tight_layout()을 호출하세요. facet (catplot)에서 suptitle이 패널과 겹치면 y=1.03처럼 위로 올리고, 그래도 잘리면 g.figure.subplots_adjust(top=0.9)로 위 여백을 줍니다.
• 점이 너무 많아 산점도가 뭉개짐(overplotting) — 볼케이노에 유전자가 수만 개면 점이 겹쳐 새까매집니다. s로 점을 작게(s=8), alpha로 반투명하게(alpha=0.4) 주거나, 비유의 점만 rasterized=True로 래스터화해 SVG 용량을 줄입니다. 유의 유전자만 색을 진하게, 나머지는 회색으로 깔면 가독성이 올라갑니다.

확장 — 더 해 볼 것

• seaborn.objects 인터페이스 — 0.12부터 들어온 새 문법(so.Plot(df, x=..., y=...).add(so.Dot()))으로, ggplot2처럼 레이어를 +가 아닌 .add()로 쌓습니다. 아직 실험적(experimental)이지만 복잡한 그림을 선언적으로 조립하기에 좋습니다.
• 상호작용 그림 — 발표용 정적 그림 너머로, plotly·bokeh를 쓰면 마우스를 올려 유전자 이름을 보는 인터랙티브 볼케이노를 만들 수 있습니다. 탐색 단계에서 유용합니다.
• 전용 패키지 — 차등발현 시각화에 특화된 pydeseq2·sanbomics의 볼케이노 헬퍼, 라벨 겹침을 자동으로 푸는 adjustText, 출판용 조판을 돕는 matplotlib의 constrained_layout 등을 함께 쓰면 손이 줄어듭니다.
• PCA·UMAP 산점도 — 샘플 전체 구조를 보는 차원축소 그림은 단일세포 분석에서 특히 중요합니다. scanpy 글에서 sc.pl.umap으로 이어집니다.
• 색·테마 다듬기 — 저널마다 요구하는 그림 규격(폰트 크기·선 두께·색)이 다릅니다. rcParams를 묶은 스타일 시트(plt.style.use("my.mplstyle"))를 하나 만들어 두면 모든 그림을 일관되게 찍을 수 있습니다.

References

1. Hunter, J. D. (2007). Matplotlib: A 2D Graphics Environment. Computing in Science & Engineering, 9(3), 90–95. (matplotlib 원논문)
2. Waskom, M. L. (2021). seaborn: statistical data visualization. Journal of Open Source Software, 6(60), 3021. (seaborn 원논문)
3. The Matplotlib development team. (2026). Matplotlib documentation — v3.11. matplotlib.org/stable
4. Waskom, M. seaborn 0.13 release notes — 범주형 함수 재작성·hue/팔레트 규칙 변경·native_scale. whatsnew/v0.13.0
5. seaborn. seaborn.clustermap — 계층 군집 히트맵 API. clustermap docs

Pipette & Pipeline · A bio portfolio journal

TL;DR — 단일세포 RNA-seq (scRNA-seq) 데이터를 파이썬으로 분석하는 사실상 표준이 scanpy (라이브러리)입니다. 공개 데이터 pbmc3k (말초혈액 단핵세포 2,700개)를 받아서 품질관리(QC) → 정규화(normalization) → 고변이 유전자(HVG) → 주성분분석(PCA) → 이웃그래프(neighbors) → UMAP → Leiden 군집화(clustering) → 마커 유전자(marker gene)까지, 한 번에 끝까지 가 봅니다.

핵심 자료구조는 AnnData입니다. 발현 행렬을 .X에, 세포 정보를 .obs에, 유전자 정보를 .var에 담아 한 객체로 들고 다닙니다. 분석을 진행하면 PCA·UMAP 좌표(.obsm)와 군집 라벨(.obs)이 이 객체에 차곡차곡 쌓입니다.

R의 Seurat로 단일세포 RNA-seq 분석과 같은 pbmc3k 데이터, 같은 흐름을 파이썬으로 옮긴 글입니다. 함수가 1:1로 대응하므로, Seurat을 써 봤다면 그대로 번역해 읽으면 됩니다.

🔗 관련 글  ·  R로 같은 분석을 하려면 Seurat로 단일세포 RNA-seq 분석  ·  scRNA-seq가 무엇을 보는지부터 알고 싶다면 단일세포 RNA 시퀀싱(scRNA-seq)이란?

이 글에서 만들 것

scRNA-seq 개념 글에서 봤듯이, 단일세포 RNA-seq는 조직 전체의 평균이 아니라 세포 하나하나의 발현을 따로 읽습니다. 그렇게 얻은 데이터는 세포 수천~수십만 개 × 유전자 2만~3만 개짜리 거대한 행렬입니다. 이걸 가지고 "여기 어떤 세포 종류가 몇 개 섞여 있나"를 알아내는 게 분석의 첫 목표입니다.

파이썬에서 이 작업의 표준 도구가 scanpy (Wolf et al. 2018)입니다. R 진영의 Seurat과 짝을 이루는 도구로, 대용량 데이터를 빠르고 메모리 효율적으로 다루도록 설계됐습니다. 이 글에서는 10x Genomics가 공개한 pbmc3k (건강한 사람의 말초혈액 단핵세포 2,700개)를 scanpy에 내장된 데이터로 바로 불러와, 품질관리부터 군집화·셀타입 해석까지 전 과정을 한 번 통과합니다. 모두 공개 데이터라 그대로 따라 하면 똑같이 재현됩니다.

전체 흐름은 아래 한 덩어리로 요약됩니다. 각 단계가 AnnData 객체에 결과를 더해 나가는 구조입니다.

# 단일세포 분석 핵심 흐름 한눈에 (scanpy)
import scanpy as sc
adata = sc.datasets.pbmc3k()                 # ① 데이터 로드 → AnnData
sc.pp.calculate_qc_metrics(adata, ...)       # ② 품질관리(QC) 지표 계산·필터
sc.pp.normalize_total(adata, target_sum=1e4) # ③ 정규화
sc.pp.log1p(adata)                           #    + 로그 변환
sc.pp.highly_variable_genes(adata, ...)      # ④ 고변이 유전자(HVG) 선별
sc.pp.scale(adata, max_value=10)             # ⑤ 스케일링
sc.tl.pca(adata)                             # ⑥ 주성분분석(PCA)
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)  # ⑦ 이웃그래프
sc.tl.umap(adata)                            # ⑧ UMAP (차원축소·시각화)
sc.tl.leiden(adata, flavor="igraph")         # ⑨ Leiden 군집화
sc.tl.rank_genes_groups(adata, "leiden")     # ⑩ 마커 유전자

1. 사전 준비 — 설치와 데이터

scanpy는 PyPI (파이썬 패키지 저장소)에서 바로 설치됩니다. 자료구조를 제공하는 anndata는 의존성으로 함께 깔리므로 따로 설치할 필요가 없습니다. 다만 Leiden 군집화에 필요한 그래프 라이브러리는 별도입니다. 한 번에 다 받으려면 [leiden] 확장(extra)을 쓰는 게 깔끔합니다.

# 기본 설치 (anndata 자동 포함)
pip install scanpy

# Leiden 군집화까지 한 번에 (igraph + leidenalg 함께 설치) — 권장
pip install "scanpy[leiden]"

설치가 끝나면 데이터를 불러옵니다. scanpy는 학습용 공개 데이터를 sc.datasets 아래에 내장해 두었습니다. pbmc3k()는 가공 전 원본 카운트(raw count)를 돌려줍니다.

import scanpy as sc

sc.settings.verbosity = 1            # 로그 간단히
adata = sc.datasets.pbmc3k()         # 말초혈액 단핵세포 2,700개 (10x)
adata

AnnData object with n_obs × n_vars = 2700 × 32738
    var: 'gene_ids'

여기서 돌려받은 adata가 AnnData 객체입니다. n_obs × n_vars = 2700 × 32738은 세포 2,700개 × 유전자 32,738개라는 뜻입니다. scanpy에서는 관습적으로 행(observation)이 세포, 열(variable)이 유전자입니다. R의 Seurat과 행·열이 반대라 처음엔 헷갈리는 지점이니 기억해 두세요.

AnnData 구조 — .X · .obs · .var

AnnData는 발현 행렬 하나만 들고 있는 게 아니라, 분석에 필요한 정보를 한 객체에 묶어 둡니다. 핵심은 세 칸입니다.

분석을 진행하면 이 칸들이 채워집니다. QC 지표는 .obs에, HVG 표시는 .var에, PCA·UMAP 좌표는 .obsm에 들어가는 식이죠. 그래서 한 객체만 들고 다니면 분석의 모든 상태가 따라옵니다 — 이게 AnnData를 쓰는 이유입니다.

adata.X[:5, :5].toarray()    # 발현 행렬 일부 (희소행렬 → 밀집 변환)
adata.obs.head()             # 세포별 표 (지금은 비어 있음, 곧 채워짐)
adata.var.head()             # 유전자별 표 (gene_ids 열만 있음)

유전자 이름이 중복되는 경우가 있어, 본격적인 분석 전에 한 번 정리해 두면 좋습니다.

adata.var_names_make_unique()   # 중복 유전자명에 접미사 붙여 고유화

2. 품질관리(QC) — 죽은 세포·빈 방울 걸러내기

원본 데이터에는 분석에 방해가 되는 세포가 섞여 있습니다. 이미 죽어 가는 세포(미토콘드리아 RNA 비율이 높음), 두 세포가 한 방울에 잡힌 더블릿(doublet, 유전자 수가 비정상적으로 많음), 세포 없이 주변 RNA만 잡힌 빈 방울(유전자 수가 너무 적음) 같은 것들이죠. 이걸 걸러내는 단계가 품질관리(QC)입니다.

먼저 미토콘드리아 유전자에 표시를 답니다. 사람 유전자에서 미토콘드리아 유전자는 이름이 MT-로 시작합니다.

# 미토콘드리아 유전자 표시 (이름이 MT- 로 시작)
adata.var["mt"] = adata.var_names.str.startswith("MT-")

# QC 지표 계산 → .obs 와 .var 에 자동으로 열이 추가됨
sc.pp.calculate_qc_metrics(
    adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True
)
adata.obs[["n_genes_by_counts", "total_counts", "pct_counts_mt"]].head()

                  n_genes_by_counts  total_counts  pct_counts_mt
AAACATACAACCAC-1                779        2419.0       3.017776
AAACATTGAGCTAC-1               1352        4903.0       3.793596
AAACATTGATCAGC-1               1129        3147.0       0.889736
AAACCGTGCTTCCG-1                960        2639.0       1.743085
AAACCGTGTATGCG-1                521         980.0       1.224490

calculate_qc_metrics가 세포별 핵심 지표 세 가지를 .obs에 만들어 줍니다. n_genes_by_counts는 세포가 발현한 유전자 수, total_counts는 총 카운트(=라이브러리 크기), pct_counts_mt는 미토콘드리아 유전자가 차지하는 비율(%)입니다. 분포를 그림으로 먼저 보고 문턱값을 정합니다.

# 세 지표의 분포를 바이올린 플롯으로
sc.pl.violin(
    adata, ["n_genes_by_counts", "total_counts", "pct_counts_mt"],
    jitter=0.4, multi_panel=True,
)
# 유전자 수 vs 미토 비율 산점도 (관계 확인)
sc.pl.scatter(adata, x="total_counts", y="pct_counts_mt")
sc.pl.scatter(adata, x="total_counts", y="n_genes_by_counts")

분포를 보고 나면 필터링입니다. pbmc3k 표준 튜토리얼이 쓰는 문턱값은 이렇습니다 — 유전자 200개 미만 세포 제거, 3개 미만 세포에서만 잡힌 유전자 제거, 유전자 수 2,500개 초과(더블릿 의심) 제거, 미토 비율 5% 초과(손상 세포) 제거.

# 세포·유전자 필터
sc.pp.filter_cells(adata, min_genes=200)    # 유전자 200개 미만 세포 제거
sc.pp.filter_genes(adata, min_cells=3)      # 3개 미만 세포에서만 잡힌 유전자 제거

# 더블릿·손상 세포 제거 (불리언 인덱싱)
adata = adata[adata.obs.n_genes_by_counts < 2500, :].copy()
adata = adata[adata.obs.pct_counts_mt < 5, :].copy()
adata

AnnData object with n_obs × n_vars = 2638 × 13714
    obs: 'n_genes_by_counts', 'total_counts', ..., 'pct_counts_mt'
    var: 'gene_ids', 'mt', 'n_cells_by_counts', ...

세포가 2,700개에서 2,638개로, 유전자가 32,738개에서 13,714개로 줄었습니다. 문턱값은 데이터마다 다릅니다 — 조직·실험 플랫폼에 따라 분포가 달라지므로, 위 그림을 보고 그때그때 정해야 합니다. 숫자를 기계적으로 베끼지 말고 분포에 맞추세요.

3. 정규화(normalization)와 로그 변환

세포마다 잡힌 총 카운트(시퀀싱 깊이)가 다릅니다. 어떤 세포는 5,000개, 어떤 세포는 1,000개가 잡혔다면 그 차이는 생물학이 아니라 기술적 변동입니다. 이걸 맞춰 주는 게 정규화입니다. scanpy의 표준은 세포마다 총합을 같은 값(보통 10,000)으로 맞춘 뒤(normalize_total), 로그를 취하는(log1p) 두 단계입니다.

# 나중을 위해 원본(로그 정규화 직전)을 .raw 에 보관 — 마커 시각화에 씀
sc.pp.normalize_total(adata, target_sum=1e4)   # 세포별 총합을 1만으로 맞춤 (CP10k)
sc.pp.log1p(adata)                             # log(1+x) 변환
adata.raw = adata                              # 이 상태(전체 유전자)를 .raw 에 저장

normalize_total(target_sum=1e4)는 흔히 CP10k (counts per 10,000)라 부릅니다. 이어지는 log1p는 log(1+x)를 취해 발현값의 분포를 다루기 쉽게 펴 줍니다. 1을 더하는 건 발현이 0인 유전자의 로그가 음의 무한대로 가는 걸 막기 위해서입니다. 참고로 R의 Seurat에서 NormalizeData() 기본값도 이 CP10k + 로그와 같은 방식입니다.

adata.raw = adata로 이 시점의 상태를 따로 저장해 둔 점에 주목하세요. 다음 단계에서 유전자를 HVG만 남기고 잘라낼 텐데, 마커 유전자를 그림으로 볼 때는 잘려 나간 유전자까지 포함한 전체 발현이 필요합니다. .raw에 넣어 두면 그때 꺼내 쓸 수 있습니다.

4. 고변이 유전자(HVG) 선별

3만 개 유전자를 다 쓰면 대부분은 세포 간 차이가 거의 없는 '배경'이라 군집을 흐립니다. 세포 종류를 가르는 정보는 세포마다 발현이 크게 출렁이는 유전자에 들어 있습니다. 이걸 골라내는 게 고변이 유전자(highly variable genes, HVG) 선별입니다.

# 고변이 유전자 선별 (기본 flavor="seurat", 로그 정규화 데이터 입력)
sc.pp.highly_variable_genes(
    adata, min_mean=0.0125, max_mean=3, min_disp=0.5
)
sc.pl.highly_variable_genes(adata)    # 평균발현 vs 분산 산점도 (HVG 강조)

# HVG 만 남기고 잘라내기
adata = adata[:, adata.var.highly_variable].copy()
adata

AnnData object with n_obs × n_vars = 2638 × 1838
    ...
    var: ..., 'highly_variable', 'means', 'dispersions', ...

유전자가 13,714개에서 1,838개로 줄었습니다. 이 1,838개가 세포 종류를 구분하는 데 정보가 많은 유전자들입니다. min_mean·max_mean·min_disp는 "평균 발현이 이 범위이고 분산이 이 값 이상인 유전자"를 HVG로 본다는 기준입니다. 위 숫자는 pbmc3k 표준 튜토리얼 값입니다.

5. 스케일링과 PCA — 차원 줄이기

군집화 전에 두 가지를 더 합니다. 먼저 각 유전자의 값을 평균 0, 분산 1로 맞추는 스케일링(scaling)입니다. 발현이 큰 유전자가 거리 계산을 지배하지 않도록 균형을 맞추는 작업이죠. 그다음 주성분분석(PCA)으로 1,838차원을 수십 차원으로 압축합니다.

# (선택) 총 카운트·미토 비율의 영향을 회귀로 제거 — 느리므로 생략 가능
# sc.pp.regress_out(adata, ["total_counts", "pct_counts_mt"])

sc.pp.scale(adata, max_value=10)      # 유전자별 z-score, 최대 10으로 클립
sc.tl.pca(adata, svd_solver="arpack") # 주성분분석

sc.pl.pca_variance_ratio(adata, n_pcs=50, log=True)  # 주성분별 설명력 (스크리 플롯)

scale의 max_value=10은 z-score를 ±10에서 잘라(clip) 극단값이 분석을 흔드는 걸 막습니다. pca는 발현 패턴이 비슷한 유전자들을 묶어 '주성분(principal component)'이라는 새 축으로 데이터를 다시 표현합니다. 스크리 플롯(scree plot)에서 설명력이 꺾이는 지점을 보고, 다음 단계에서 쓸 주성분 개수(n_pcs)를 정합니다. pbmc3k에서는 보통 40개 안팎이면 충분합니다.

참고로 scanpy에는 sc.tl.pca와 sc.pp.pca가 둘 다 있고 결과는 같습니다. 표준 튜토리얼은 sc.tl.pca를 쓰는데, 최신 버전에서 더 적극적으로 관리되는 쪽은 sc.pp.pca입니다. 입문 단계에서는 어느 쪽을 써도 무방합니다.

6. 이웃그래프(neighbors)와 UMAP

이제 세포들을 PCA 공간에서의 가까움에 따라 그래프로 잇습니다. 각 세포를 가장 가까운 이웃 세포들과 연결한 k-최근접이웃(k-nearest neighbor) 그래프이고, 이 그래프가 군집화와 UMAP 양쪽의 토대가 됩니다.

# 이웃그래프: 주성분 40개를 써서 세포당 이웃 10개 연결
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)

# UMAP: 그래프를 2차원에 펼쳐 시각화
sc.tl.umap(adata)
sc.pl.umap(adata)     # 아직 색칠 정보 없음 — 점만 펼쳐진 지형도

neighbors가 그래프를 만들면, UMAP (uniform manifold approximation and projection)이 그 그래프를 2차원 평면에 펼칩니다. 가까운 세포는 가깝게, 먼 세포는 멀게 배치해 세포들의 '지형도'를 그리는 차원축소(dimensionality reduction) 기법입니다. 이 시점의 UMAP은 색이 없어 점 구름만 보입니다. 색은 다음 단계에서 군집을 입혀 채웁니다. 참고로 UMAP은 시각화용이고, 거리·면적의 절대값에 큰 의미를 두면 안 됩니다 — 군집의 모양보다 '어떤 세포끼리 뭉치는가'를 보는 도구입니다.

7. Leiden 군집화(clustering)

군집화는 이웃그래프를 촘촘히 연결된 덩어리로 나누는 작업입니다. scanpy의 표준은 Leiden 알고리즘(Traag et al. 2019)입니다. 이전에 쓰이던 Louvain을 개량해, 잘게 끊긴 군집이 생기지 않도록 보장합니다.

# Leiden 군집화 — 최신 권장 호출 방식
sc.tl.leiden(
    adata,
    resolution=1.0,
    flavor="igraph",     # 빠른 igraph 백엔드 (scanpy 1.10+)
    n_iterations=2,
    directed=False,
)
sc.pl.umap(adata, color="leiden")    # 군집별 색으로 UMAP 색칠
adata.obs["leiden"].value_counts()   # 군집별 세포 수

leiden
0    647
1    623
2    441
3    323
...
Name: count, dtype: int64

UMAP에 군집별 색이 입혀지면서 세포 종류가 시각적으로 갈라집니다. resolution이 군집의 해상도를 정합니다 — 키우면 잘게 쪼개지고, 줄이면 뭉칩니다. pbmc3k에서는 1.0 안팎이면 PBMC의 주요 세포 종류 8개 정도로 잘 나뉩니다.

8. 마커 유전자(marker gene)와 셀타입 해석

군집을 8개로 나눴지만, 아직 각 군집이 무슨 세포인지는 모릅니다. 군집 번호(0, 1, 2, …)는 그냥 라벨일 뿐이죠. 정체를 밝히려면 각 군집에서 특이적으로 높게 발현하는 유전자(마커 유전자)를 찾고, 그게 알려진 세포 종류의 표지와 맞는지 봅니다.

# 군집별 마커 유전자 검정 (Wilcoxon 순위합 검정)
sc.tl.rank_genes_groups(adata, "leiden", method="wilcoxon")
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)

# 결과를 데이터프레임으로 — 군집 0의 상위 마커
sc.get.rank_genes_groups_df(adata, group="0").head()

   names    scores  logfoldchanges         pvals     pvals_adj
0  LDHB   34.21...        1.06...   2.1e-200      8.3e-197
1  RPS27  31.5...         0.74...   ...           ...
2  IL7R   27.8...         1.20...   ...           ...
...

rank_genes_groups는 각 군집을 나머지 전체와 비교해, 그 군집에서 유의하게 높은 유전자를 순위로 매깁니다. method="wilcoxon" (윌콕슨 순위합 검정)이 단일세포에서 가장 널리 쓰이는 방법입니다. 결과는 sc.get.rank_genes_groups_df로 표(데이터프레임)로 꺼낼 수 있고, names·scores·logfoldchanges·pvals·pvals_adj 열을 줍니다.

이제 PBMC에서 잘 알려진 마커들을 UMAP과 닷플롯(dot plot)으로 직접 확인합니다.

# PBMC 표준 마커 유전자
marker_genes = [
    "IL7R", "CD14", "LYZ", "MS4A1", "CD8A",
    "GNLY", "NKG7", "FCGR3A", "MS4A7", "FCER1A", "CST3", "PPBP",
]
sc.pl.umap(adata, color=["leiden", "CD14", "MS4A1", "NKG7"])  # 군집 + 마커 발현
sc.pl.dotplot(adata, marker_genes, groupby="leiden")          # 마커 × 군집 닷플롯

닷플롯은 점의 색으로 평균 발현, 크기로 발현 세포 비율을 보여 줘서, 어떤 군집이 어떤 마커를 켜는지 한눈에 들어옵니다. 알려진 마커-셀타입 대응은 이렇습니다.

각 군집의 상위 마커를 이 표와 맞춰 보고, 군집에 이름을 붙입니다.

# 군집 번호 → 셀타입 이름 (각자 데이터의 마커로 순서 확인 필수!)
new_names = [
    "CD4 T", "CD14+ Monocytes", "B", "CD8 T", "NK",
    "FCGR3A+ Monocytes", "Dendritic", "Megakaryocytes",
]
adata.rename_categories("leiden", new_names)
sc.pl.umap(adata, color="leiden", legend_loc="on data")   # 셀타입 라벨 UMAP

여기까지가 단일세포 분석의 한 사이클입니다. 결과를 파일로 저장하면 다음에 이어서 쓸 수 있습니다.

adata.write("pbmc3k_processed.h5ad")   # AnnData → HDF5 파일로 저장

9. 자주 발생하는 에러 & 해결법

• ModuleNotFoundError: No module named 'leidenalg' (또는 igraph) — 입문자가 가장 많이 만나는 에러입니다. pip install scanpy만으로는 군집화 백엔드가 안 깔립니다. pip install "scanpy[leiden]" (또는 pip install igraph leidenalg)로 해결하세요. flavor="igraph"를 쓰면 igraph가, flavor="leidenalg"를 쓰면 leidenalg가 필요합니다.
• Leiden FutureWarning — flavor 없이 sc.tl.leiden을 부르면 "기본 백엔드가 igraph로 바뀐다"는 경고가 뜹니다. 에러는 아니지만, flavor="igraph", n_iterations=2, directed=False를 넘겨 미래 기본값으로 맞춰 두면 경고가 사라집니다. 단 이 인자는 scanpy 1.10 이상에서만 됩니다 — 구버전에서 넘기면 TypeError: leiden() got an unexpected keyword argument 'flavor'가 납니다.
• igraph vs python-igraph 이름 혼동 — PyPI (pip)에서는 igraph가 정식 이름이고 python-igraph는 그쪽으로 넘겨주는 별칭입니다. 반면 conda에서는 여전히 python-igraph로 설치합니다(conda install -c conda-forge python-igraph leidenalg). 둘 다 깔려 있으면 임포트가 꼬일 수 있으니 pip 환경에선 igraph 하나만 두세요.
• log1p 이중 적용 — 정규화 셀을 두 번 실행해 로그를 두 번 취하면 값이 망가집니다. 증상은 발현값이 비정상적으로 작게 눌리는 것입니다. scanpy가 경고로 막아 주지만, 노트북에서 셀을 반복 실행하지 않도록 주의하세요.
• HVG flavor와 입력 데이터 불일치 — flavor="seurat"(기본)은 로그 정규화 데이터를, flavor="seurat_v3"는 원본 카운트를 받습니다. 잘못 넣어도 에러 없이 엉뚱한 유전자가 뽑히는 조용한 버그라 더 위험합니다. flavor와 데이터 상태를 꼭 맞추세요.
• scanpy ↔ anndata 버전 불일치 — 오래된 anndata에 최신 scanpy를 쓰면 ImportError·AttributeError가 납니다. scanpy 1.11+는 anndata 0.10 이상이 필요합니다. pip install -U scanpy anndata로 함께 올리세요.
• 메모리·속도 문제 — 세포가 수십만 개로 커지면 regress_out과 밀집(dense) scale이 무겁고, 수렴까지 도는 Leiden(n_iterations=-1)이 느립니다. 행렬을 희소(sparse)로 유지하고, PCA 전에 HVG로 먼저 줄이고, flavor="igraph", n_iterations=2를 쓰면 한결 가벼워집니다.

10. R Seurat과의 대응, 그리고 확장

이 글의 흐름은 R의 Seurat 분석과 사실상 같은 단계를 밟습니다. 같은 pbmc3k 데이터를 쓰니, 함수만 갈아 끼우면 두 도구를 오갈 수 있습니다.

대응이 거의 1:1입니다. 한 가지 개념적 주의는, Seurat의 vst 변동유전자 방식이 scanpy의 flavor="seurat_v3"(원본 카운트)에 해당하지, 기본값 "seurat"가 아니라는 점입니다.

여기서 더 나아가는 길도 자연스럽게 이어집니다. 셀타입 자동 주석으로는 참조 데이터에 맞춰 라벨을 붙이는 celltypist나, 마커 기반 점수를 매기는 도구들이 있습니다. 궤적 분석(trajectory)으로는 세포가 분화·활성화하며 변해가는 순서를 재구성하는 scanpy.tl.paga나 별도 패키지 scvelo(RNA velocity)·Palantir로 확장합니다. 군집이 "어떤 세포가 있나"를 답한다면, 궤적은 "어떤 순서로 변해가나"를 그립니다. scanpy 생태계(scverse)는 이런 후속 분석 도구가 같은 AnnData 위에서 맞물리도록 설계돼 있어, 한 객체를 그대로 넘기며 분석을 이어 갈 수 있습니다.

References

1. Wolf, F. A., Angerer, P., & Theis, F. J. (2018). SCANPY: large-scale single-cell gene expression data analysis. Genome Biology, 19, 15. https://doi.org/10.1186/s13059-017-1382-0 (scanpy 원논문)
2. Virshup, I., Rybakov, S., Theis, F. J., Angerer, P., & Wolf, F. A. (2024). anndata: Annotated data. Journal of Open Source Software, 9(101), 4371. https://doi.org/10.21105/joss.04371 (anndata)
3. Traag, V. A., Waltman, L., & van Eck, N. J. (2019). From Louvain to Leiden: guaranteeing well-connected communities. Scientific Reports, 9, 5233. https://doi.org/10.1038/s41598-019-41695-z (Leiden 알고리즘)
4. McInnes, L., Healy, J., & Melville, J. (2018). UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction. arXiv:1802.03426. https://arxiv.org/abs/1802.03426 (UMAP)
5. Scanpy 공식 문서 — PBMC3k 클러스터링 튜토리얼. https://scanpy.readthedocs.io/en/stable/tutorials/basics/clustering-2017.html
6. Scanpy 공식 문서 (API). https://scanpy.readthedocs.io/en/stable/api/index.html

Pipette & Pipeline · A bio portfolio journal

TL;DR — RNA-seq 카운트 행렬(genes × samples)을 받아 분석에 넣기 직전까지, "불러오기 → 샘플 메타데이터 결합 → 필터·집계 → wide↔long 변환 → 결측·중복 정리 → 저장"을 전부 pandas (파이썬 데이터 처리 라이브러리)로 처리합니다. 바이오 데이터 정제의 90%는 이 흐름 안에 있습니다.

핵심 도구는 두 자료구조입니다. 1차원 벡터인 Series (시리즈)와 2차원 표인 DataFrame (데이터프레임). read_csv로 표를 읽고, merge로 메타데이터를 붙이고, query·불리언 마스크로 거르고, groupby로 그룹 평균을 내고, melt/pivot으로 표 모양을 바꿉니다.

R의 dplyr·tidyr를 써 봤다면 1:1로 옮겨 탈 수 있습니다. filter()↔query(), mutate()↔assign(), group_by()+summarise()↔groupby()+agg(). 이 글 끝에 대응표를 정리했습니다.

🔗 관련 글  ·  R로 같은 작업을 한다면 dplyr·tidyr로 데이터 정리하기와 1:1로 비교됩니다  ·  공개 발현 데이터를 받는 법은 GEOquery로 공개 발현 데이터 받기  ·  카운트를 비교 가능한 값으로 바꾸는 RNA-seq 정규화

이 글에서 만들 것

전사체(transcriptome) 분석을 하다 보면 가장 먼저 마주치는 건 멋진 그래프가 아니라 모양이 어정쩡한 표입니다. 유전자가 행, 샘플이 열인 카운트 행렬(count matrix) 하나, 그리고 샘플마다 조건·배치·품질값이 적힌 메타데이터(metadata) 표 하나. 이 둘을 합치고 다듬어 "control 대비 treated에서 평균 발현이 어떻게 다른가" 같은 질문에 답할 수 있는 형태로 만드는 일이 분석의 절반입니다.

pandas (Wes McKinney, 2008~)는 파이썬에서 이 표 작업을 맡는 표준 라이브러리입니다. 엑셀의 표를 코드로 다룬다고 생각하면 가깝습니다. 다만 손으로 클릭하는 대신 한 줄씩 명령으로 적기 때문에, 같은 처리를 수백 개 샘플에 그대로 반복하고 재현할 수 있습니다.

이 글은 작은 발현 카운트 행렬을 직접 만들어, 실제 분석에서 거치는 순서대로 정제합니다. CSV·TSV로 불러오기, 샘플 메타데이터를 merge로 결합, 불리언·query로 필터, groupby로 그룹별 평균, melt/pivot으로 wide↔long 변환, 결측(NaN)·중복 처리, 정렬과 인덱스(index) 다루기, 그리고 저장까지입니다. 코드는 모두 pandas 3.0에서 실제로 돌려 확인한 것이고, 무작위 시드를 고정해 그대로 따라 하면 같은 결과가 나옵니다.

# 핵심 흐름 한눈에 (pandas)
import pandas as pd
df   = pd.read_csv("counts.csv", index_col=0)        # ① 카운트 행렬 불러오기
long = df.reset_index().melt(id_vars="gene_id",      # ② wide → long
                             var_name="sample_id", value_name="count")
long = long.merge(meta, on="sample_id", how="left")  # ③ 메타데이터 결합
sub  = long.query("group == 'treated' and count > 55")  # ④ 필터
grp  = long.groupby(["gene_id", "group"])["count"].mean()  # ⑤ 그룹 집계
wide = grp.reset_index().pivot(index="gene_id",       # ⑥ long → wide
                               columns="group", values="count")

1. 사전 준비 — 설치와 샘플 데이터

pandas는 pip이나 conda로 설치합니다. 수치 연산 백엔드인 NumPy (수치 배열 라이브러리)는 함께 깔립니다. 2026년 현재 안정 버전은 3.0.x 계열이며, 이 글의 코드도 3.0 기준입니다.

# pip로 설치 (가상환경 권장)
pip install pandas

# 또는 conda
conda install -c conda-forge pandas

import pandas as pd
import numpy as np

print(pd.__version__)

3.0.3

이제 설명용으로 작은 카운트 행렬을 만듭니다. 유전자 8개 × 샘플 6개, 값은 포아송(Poisson) 분포에서 뽑은 가짜 카운트입니다. 실제 분석이라면 이 행렬은 featureCounts·Salmon 같은 정량 도구의 출력이거나, 공개 데이터베이스(GEO 등)에서 받은 매트릭스가 됩니다.

import pandas as pd
import numpy as np

rng = np.random.default_rng(0)              # 재현성 고정 (시드 0)
genes   = [f"GENE{i:03d}" for i in range(1, 9)]   # GENE001 ~ GENE008
samples = ["S1", "S2", "S3", "S4", "S5", "S6"]

# 발현 카운트 행렬: 유전자(행) × 샘플(열)
counts = pd.DataFrame(
    rng.poisson(lam=50, size=(8, 6)).astype("int64"),
    index=genes, columns=samples,
)
counts.index.name = "gene_id"
counts

         S1  S2  S3  S4  S5  S6
gene_id
GENE001  53  32  52  51  57  59
GENE002  55  59  46  54  52  54
GENE003  54  41  50  50  63  51
GENE004  52  48  44  38  56  60
GENE005  46  49  44  48  39  46
GENE006  47  49  50  66  56  51
GENE007  48  40  60  50  45  49
GENE008  52  55  56  63  58  47

DataFrame은 2차원 표, 각 열은 Series (1차원)입니다. index는 행 이름(여기선 유전자 ID), columns는 열 이름(샘플 ID)입니다. 행렬을 CSV·TSV로 저장했다가 다시 읽어 보겠습니다.

# 저장: CSV와 TSV 두 가지
counts.to_csv("counts.csv")              # 쉼표 구분
counts.to_csv("counts.tsv", sep="\t")    # 탭 구분

# 불러오기: 첫 열(gene_id)을 인덱스로
df = pd.read_csv("counts.csv", index_col=0)
df.head(3)

         S1  S2  S3  S4  S5  S6
gene_id
GENE001  53  32  52  51  57  59
GENE002  55  59  46  54  52  54
GENE003  54  41  50  50  63  51

read_csv의 index_col=0은 첫 열을 행 인덱스로 쓰라는 뜻입니다. 빼먹으면 gene_id가 그냥 데이터 열로 들어가 0,1,2… 정수 인덱스가 새로 붙습니다. TSV는 sep="\t"만 더해 같은 함수로 읽습니다. 자주 쓰는 인자는 아래와 같습니다.

이제 샘플 메타데이터를 만듭니다. 각 샘플의 처리군(group), 배치(batch), RNA 품질값(RIN)을 담습니다.

# 샘플 메타데이터: 행 = 샘플
meta = pd.DataFrame({
    "sample_id": ["S1", "S2", "S3", "S4", "S5", "S6"],
    "group":     ["control", "control", "control", "treated", "treated", "treated"],
    "batch":     ["b1", "b2", "b1", "b2", "b1", "b2"],
    "rin":       [9.1, 8.7, 9.4, 7.9, 8.2, 9.0],
}).set_index("sample_id")
meta

             group batch  rin
sample_id
S1         control    b1  9.1
S2         control    b2  8.7
S3         control    b1  9.4
S4         treated    b2  7.9
S5         treated    b1  8.2
S6         treated    b2  9.0

2. wide → long 변환 — melt

카운트 행렬은 사람이 보기엔 편하지만(wide 형식), 메타데이터를 붙이거나 그룹별로 묶어 집계하려면 long 형식이 훨씬 다루기 쉽습니다. long 형식은 한 행에 (유전자, 샘플, 값) 하나씩만 담는 모양입니다. tidyr의 pivot_longer()에 대응하는 melt로 바꿉니다.

# wide(8×6) → long(48×3)
long = (counts
        .reset_index()                       # gene_id를 일반 열로 꺼냄
        .melt(id_vars="gene_id",             # 그대로 둘 식별자 열
              var_name="sample_id",          # 녹인 열 이름 → sample_id
              value_name="count"))           # 녹인 값 이름 → count
long.head()

   gene_id sample_id  count
0  GENE001        S1     53
1  GENE002        S1     55
2  GENE003        S1     54
3  GENE004        S1     52
4  GENE005        S1     46

id_vars에 적은 열(gene_id)은 그대로 남고, 나머지 열(S1~S6)은 전부 sample_id라는 한 열로 "녹아" 들어갑니다. 8 × 6 행렬이 48행짜리 긴 표가 됐습니다. 이 모양이라야 다음 단계가 자연스럽습니다.

3. 메타데이터 결합 — merge

long 표의 sample_id를 열쇠(key)로 삼아 메타데이터를 붙입니다. dplyr의 left_join()에 해당하는 merge를 씁니다.

# sample_id 기준으로 메타데이터를 왼쪽 결합
long_m = long.merge(meta, on="sample_id", how="left")
long_m.head()

   gene_id sample_id  count    group batch  rin
0  GENE001        S1     53  control    b1  9.1
1  GENE002        S1     55  control    b1  9.1
2  GENE003        S1     54  control    b1  9.1
3  GENE004        S1     52  control    b1  9.1
4  GENE005        S1     46  control    b1  9.1

on="sample_id"는 두 표에서 이름이 같은 열을 짝지으라는 뜻입니다. meta는 sample_id가 인덱스이므로, merge가 인덱스를 자동으로 열쇠로 인식해 붙여 줍니다. how는 결합 방식입니다 — "left"는 왼쪽 표(long)의 행을 모두 유지하고 오른쪽 정보를 채웁니다. 이제 각 카운트 옆에 그 샘플의 group·batch·RIN이 따라옵니다.

4. 필터 — 불리언 마스크와 query

조건에 맞는 행만 고르는 방법은 두 가지입니다. 불리언 마스크(boolean mask)와 query. dplyr의 filter()에 해당합니다.

# 방법 A: 불리언 마스크 — 조건식이 True/False Series를 만들고, .loc로 거름
mask = (long_m["group"] == "treated") & (long_m["count"] > 55)
long_m.loc[mask, ["gene_id", "sample_id", "count"]].head()

    gene_id sample_id  count
29  GENE006        S4     66
31  GENE008        S4     63
32  GENE001        S5     57
34  GENE003        S5     63
35  GENE004        S5     56

# 방법 B: query — 문자열로 조건을 적어 더 읽기 쉽게
q = long_m.query("group == 'treated' and count > 55")
# 두 결과는 동일

마스크 방식에서 조건이 여러 개면 각 조건을 괄호로 감싸고 &(그리고)·|(또는)로 잇습니다. 파이썬의 and/or가 아니라 &/|를 쓴다는 점, 그리고 우선순위 때문에 괄호가 반드시 필요하다는 점이 자주 걸리는 함정입니다. query는 조건을 문자열로 적어 한결 깔끔하고, 파이썬 변수는 @를 붙여, 공백이 든 열 이름은 백틱(`)으로 감싸 참조합니다.

# query 안에서 파이썬 변수(@)와 공백 든 열 이름(백틱) 쓰기
thresh = 55
df2 = long_m.rename(columns={"count": "raw count"})    # 공백 든 열 이름
df2.query("`raw count` > @thresh")                     # @thresh = 바깥 변수

5. 그룹 집계 — groupby

이제 "유전자별·그룹별 평균 발현"을 냅니다. dplyr의 group_by() + summarise()에 해당하는 groupby + 집계 함수입니다.

# 유전자 × 그룹마다 평균 카운트
grp = (long_m
       .groupby(["gene_id", "group"])["count"]   # 묶을 키 → 집계할 열
       .mean()
       .reset_index(name="mean_count"))          # 결과를 다시 표로
grp.head(6)

   gene_id    group  mean_count
0  GENE001  control   45.666667
1  GENE001  treated   55.666667
2  GENE002  control   53.333333
3  GENE002  treated   53.333333
4  GENE003  control   48.333333
5  GENE003  treated   54.666667

groupby(["gene_id", "group"])로 두 열의 조합마다 묶고, ["count"]로 집계 대상을 고른 뒤 .mean()을 적용합니다. 한 그룹에 여러 통계를 한꺼번에 내려면 agg에 함수 목록을 넘깁니다.

# 그룹별로 평균·표준편차·개수를 동시에
agg = long_m.groupby("group")["count"].agg(["mean", "std", "size"])
agg.round(2)

          mean   std  size
group
control  49.25  6.32    24
treated  52.62  7.08    24

6. long → wide 변환 — pivot

집계 결과를 다시 "유전자 × 그룹" 행렬로 되돌리면 한눈에 비교하기 좋습니다. melt의 반대, tidyr의 pivot_wider()에 해당하는 pivot입니다.

# long → wide: 행=유전자, 열=그룹, 값=평균 카운트
wide = grp.pivot(index="gene_id", columns="group", values="mean_count")
wide.round(2)

group    control  treated
gene_id
GENE001    45.67    55.67
GENE002    53.33    53.33
GENE003    48.33    54.67
GENE004    48.00    51.33
GENE005    46.33    44.33
GENE006    48.67    57.67
GENE007    49.33    48.00
GENE008    54.33    56.00

index는 행이 될 열, columns는 열이 될 열, values는 칸을 채울 값입니다. 이제 control과 treated의 평균을 유전자별로 나란히 볼 수 있습니다. 차이가 큰 유전자(GENE001·GENE006)가 눈에 들어옵니다.

7. 결측(NaN)·중복 처리

실제 데이터에는 빈칸과 중복이 섞여 있습니다. 메타데이터의 RIN 한 칸을 결측으로 만들어 다뤄 보겠습니다.

nm = meta.copy()
nm.loc["S3", "rin"] = np.nan      # 결측 하나 주입

nm.isna().sum()                   # 열별 결측 개수

group    0
batch    0
rin      1
dtype: int64

nm.dropna()                       # 결측이 있는 행 제거 → 5행
nm["rin"].fillna(nm["rin"].median())   # 중앙값으로 채우기 → S3 = 8.7

isna()로 어디가 비었는지 세고, dropna()로 결측 행을 버리거나 fillna()로 채웁니다. 발현 카운트는 보통 결측 행을 통째로 버리기보다 0으로 채우거나 그대로 두지만, 품질값 같은 메타데이터는 중앙값·평균으로 보간하는 편이 흔합니다. 중복은 duplicated()로 찾아 제거합니다.

dup = pd.concat([meta, meta.loc[["S1"]]])   # S1 행을 일부러 중복

dup.index.duplicated().any()                # 인덱스 중복 있나? → True
dup[~dup.index.duplicated(keep="first")]    # 첫 등장만 남김 → 6행

keep="first"는 중복 중 첫 번째만 남기라는 뜻입니다("last"도 가능). ~는 불리언을 뒤집는 부정 연산자라, "중복이 아닌 행만" 고르는 흔한 관용구입니다.

8. 정렬과 인덱스 — sort_values

마지막으로 평균 발현이 높은 유전자 순으로 정렬해 마무리합니다.

# 유전자별 평균 카운트를 내림차순 정렬
gene_means = (long_m
              .groupby("gene_id")["count"]
              .mean()
              .sort_values(ascending=False))
gene_means.round(2).head(4)

gene_id
GENE008    55.17
GENE002    53.33
GENE006    53.17
GENE003    51.50
Name: count, dtype: float64

sort_values(ascending=False)로 값 기준 내림차순, 인덱스(유전자 이름) 기준으로 정렬하려면 sort_index()를 씁니다. 정제가 끝난 표는 to_csv로 저장해 다음 분석으로 넘깁니다.

gene_means.to_csv("gene_means.csv")     # 결과 저장

R dplyr · tidyr 1:1 대응표

R에서 데이터 정제를 해 봤다면, 거의 모든 동사를 그대로 옮길 수 있습니다.

가장 큰 차이는 인덱스(index)입니다. R 데이터프레임에는 행 이름이라는 개념이 약하지만, pandas는 행 인덱스를 1급 객체로 다룹니다. 그래서 set_index·reset_index로 "어떤 열을 인덱스로 둘지"를 자주 바꾸게 됩니다. 또 dplyr의 파이프(%>%/|>)처럼, pandas도 .으로 메서드를 줄줄이 이어 한 흐름으로 적을 수 있습니다.

9. 자주 발생하는 에러 & 오해

• SettingWithCopyWarning — 슬라이스에 값을 쓸 때 — df[df.group == "treated"]["count"] = 999처럼 거른 결과(슬라이스)에 다시 대입하면, 원본이 바뀔지 사본만 바뀔지 모호합니다. pandas 2.x에서는 경고가 떴고, 3.0부터는 Copy-on-Write가 기본이라 경고 없이 원본은 그대로 두고 사본만 바뀝니다. 의도가 "원본 수정"이라면 한 번에 df.loc[df.group == "treated", "count"] = 999처럼 .loc로 행·열을 함께 지정하세요. 거른 표를 따로 쓸 거라면 .copy()를 붙여 의도를 분명히 합니다.
• 숫자처럼 생긴 ID가 정수로 읽힘 — gene_id가 001·002라면 read_csv가 앞의 0을 떼고 정수 1·2로 읽어 버립니다. 유전자 ID·바코드·샘플명처럼 식별자는 dtype={"gene_id": "str"}로 문자열임을 명시하세요.
• 문자열 열의 dtype가 object가 아니다 — pandas 3.0부터 문자열 열의 기본 자료형이 object에서 전용 str 타입으로 바뀌었습니다. df["group"].dtype이 str로 나오는 게 정상입니다. 옛 코드에서 dtype == object로 문자열 열을 골랐다면 동작이 달라질 수 있으니 확인하세요.
• merge 후 행 수가 늘어남 — 열쇠가 한쪽에서 중복이면 결합이 곱(cartesian) 형태로 불어납니다. 결합 전에 meta.index.is_unique나 df["sample_id"].duplicated().any()로 열쇠가 유일한지 점검하는 습관이 안전합니다.
• 인덱스 정렬 때문에 값이 어긋남 — 인덱스가 다른 두 Series를 더하거나 나누면, pandas는 위치가 아니라 인덱스 이름으로 자동 정렬합니다. 의도와 달리 NaN이 생기거나 값이 섞이면 인덱스부터 의심하세요. 위치 기준으로 다루고 싶으면 .values(NumPy 배열)로 꺼내거나 .reset_index(drop=True)로 인덱스를 맞춥니다.
• inplace=True 남발 — inplace는 코드를 짧게 만드는 듯하지만 메서드 체이닝을 끊고, 3.0에서는 대부분 성능 이점도 없습니다. df = df.fillna(0)처럼 새로 받는 편이 읽기에도 안전에도 낫습니다.

References

1. McKinney, W. (2010). Data Structures for Statistical Computing in Python. Proceedings of the 9th Python in Science Conference, 56–61. (pandas 원논문)
2. The pandas development team. (2026). pandas documentation — v3.0.3. pandas.pydata.org/docs
3. pandas. What's new in 3.0.0 — Copy-on-Write 기본화·문자열 dtype. whatsnew/v3.0.0
4. Wickham, H. et al. (2019). Welcome to the tidyverse. Journal of Open Source Software, 4(43), 1686. (dplyr·tidyr 대응 참고)

Pipette & Pipeline · A bio portfolio journal

TL;DR — 서열(sequence) 데이터를 Python으로 다루는 표준 라이브러리가 Biopython입니다. Bio.SeqIO로 FASTA·FASTQ 파일을 읽고 쓰며, Seq 객체 하나로 역상보(reverse-complement)·전사(transcription)·번역(translation)·GC 함량(GC content)을 한 줄씩 계산합니다. 이 글은 설치부터 시작해 NCBI에서 받은 실제 유전자 서열로 길이·GC·ORF 통계를 내고, FASTQ 품질점수(Phred quality score)에 접근하는 데까지 갑니다.

흐름은 pip install biopython → SeqIO.read/SeqIO.parse로 파일 읽기 → record.seq로 Seq 객체 꺼내기 → len·gc_fraction·translate로 통계 → letter_annotations["phred_quality"]로 FASTQ 품질 접근입니다. 마지막에 결과를 pandas 데이터프레임으로 넘기는 한 줄까지 붙입니다.

이 글은 "Python for Bio" 시리즈의 출발점입니다. 개념으로 알던 FASTA·FASTQ를 실제 코드로 처음 만져 보는 다리 역할을 합니다.

🔗 관련 글  ·  개념 먼저 잡으려면 FASTQ 파일이란?(4줄 구조·Phred 점수)  ·  시퀀싱 원리는 NGS란? 1·2·3세대 비교  ·  번역(translation)의 생물학은 중심원리(Central Dogma)

이 글에서 만들 것

NGS 시퀀서가 뱉어내는 파일은 거의 다 FASTA 아니면 FASTQ입니다. 둘 다 그냥 텍스트 파일이라 메모장으로도 열리지만, 수만 개의 서열을 손으로 셀 수는 없습니다. 여기서 Biopython (Python용 생물정보학 라이브러리)이 등장합니다. 2000년부터 이어진 오픈소스 프로젝트로, 서열 입출력·정렬·번역·데이터베이스 접근 같은 일상 작업을 표준화된 객체 몇 개로 묶어 둔 도구입니다.

이 글에서는 다음을 직접 해 봅니다. 우선 Biopython을 설치하고, NCBI에서 실제 유전자 서열을 FASTA로 내려받습니다. 그 서열의 길이와 GC 함량을 재고, 역상보·전사·번역을 한 줄씩 계산합니다. 이어서 여섯 개의 번역틀(reading frame)을 훑어 ORF (open reading frame, 열린 해독틀) 후보를 세고, FASTQ 파일에서는 각 염기의 Phred 품질점수를 꺼내 평균을 냅니다. 마지막으로 이 결과들을 pandas 데이터프레임 한 줄로 옮겨, 그래프나 통계로 넘어갈 준비를 합니다.

# 핵심 흐름 한눈에 (Biopython)
from Bio import SeqIO
from Bio.SeqUtils import gc_fraction

record = SeqIO.read("gene.fasta", "fasta")   # ① FASTA 한 서열 읽기
seq = record.seq                              # ② Seq 객체 꺼내기
print(len(seq))                               # ③ 서열 길이
print(gc_fraction(seq) * 100)                 # ④ GC 함량 (%)
print(seq.reverse_complement())               # ⑤ 역상보 가닥
print(seq.translate(to_stop=True))            # ⑥ 단백질로 번역

1. 사전 준비 — 설치와 샘플 데이터

Biopython은 pip이나 conda 어느 쪽으로도 설치됩니다. 2025년 10월에 나온 1.86 버전 기준으로 Python 3.10~3.14를 지원합니다. 의존성은 수치 계산용 NumPy 하나뿐이라 설치가 가볍습니다.

# pip 사용 (가상환경 권장)
pip install biopython

# 또는 conda (conda-forge 채널)
conda install -c conda-forge biopython

설치가 끝났으면 버전을 확인합니다. 1.78 이전과 이후는 API가 꽤 달라졌으니, 옛날 블로그 코드를 참고할 때는 버전부터 맞춰 보는 게 좋습니다.

import Bio
print(Bio.__version__)

1.86

샘플 데이터는 NCBI (미국 국립생물정보센터)에서 실제 서열을 받아 씁니다. Biopython의 Bio.Entrez 모듈이 NCBI에 직접 질의를 보냅니다. 여기서는 대장균(Escherichia coli)의 lacZ 유전자(β-갈락토시다아제를 만드는 유전자, GenBank 등록번호 J01636.1)를 FASTA로 내려받습니다. NCBI 정책상 Entrez.email에 본인 이메일을 적어야 하며, 이는 과도한 요청을 막기 위한 식별용입니다.

from Bio import Entrez, SeqIO

Entrez.email = "your_email@example.com"   # 본인 이메일 (NCBI 요구사항)

# efetch: nucleotide DB에서 J01636.1 을 FASTA 형식으로 받기
with Entrez.efetch(db="nucleotide", id="J01636.1",
                   rettype="fasta", retmode="text") as handle:
    record = SeqIO.read(handle, "fasta")   # 서열 1개 → read()

# 받은 서열을 로컬 파일로 저장 (다음 단계부터 재사용)
SeqIO.write(record, "lacz.fasta", "fasta")
print(record.id, "·", len(record.seq), "bp")

J01636.1 · 7477 bp

인터넷 연결이 막혀 있거나 NCBI 접속이 어려우면, 아래처럼 짧은 더미 서열을 직접 만들어 똑같이 따라올 수 있습니다. 형식만 FASTA면 되므로 분석 코드는 그대로 동작합니다.

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord

dummy = SeqRecord(
    Seq("ATGACCATGATTACGGATTCACTGGCCGTCGTTTTACAACGTCGTGACTGGGAAAACCCTGGCGTTACCCAACTTAATCGCCTTGCAGCACATCCCCCTTTCGCCAGCTGGCGTAATAGCGAAGAGGCCCGCACCGATCGCCCTTCCCAACAGTTGCGCAGCCTGAATGGCGAATGGTAA"),
    id="demo_gene", description="dummy CDS for tutorial")
SeqIO.write(dummy, "lacz.fasta", "fasta")   # 위와 같은 파일명으로 저장

2. SeqRecord 구조 — 읽어 온 서열의 속살

SeqIO.read(서열 1개)나 SeqIO.parse(여러 개)가 돌려주는 건 단순한 문자열이 아니라 SeqRecord 객체입니다. 서열 본체뿐 아니라 식별자·설명·주석까지 한 묶음으로 들고 있습니다. 어떤 정보가 들었는지 먼저 들여다봅니다.

from Bio import SeqIO

record = SeqIO.read("lacz.fasta", "fasta")   # 저장해 둔 파일 다시 읽기

print("id        :", record.id)          # 서열 식별자 (>뒤 첫 단어)
print("desc      :", record.description)  # 설명 (>뒤 전체)
print("길이      :", len(record.seq), "bp")
print("앞 30 bp  :", record.seq[:30])     # Seq 객체도 슬라이싱 가능
print("타입      :", type(record.seq))

id        : J01636.1
desc      : J01636.1 E.coli lactose operon with lacI, lacZ, lacY and lacA genes
길이      : 7477 bp
앞 30 bp  : GGGGCTGGCTTAACTATGCGGCATCAGAGC
타입      : <class 'Bio.Seq.Seq'>

핵심은 record.seq로 꺼낸 Seq 객체입니다. 겉보기엔 문자열 같지만, 생물학적 의미가 있는 메서드(역상보·전사·번역)를 품고 있는 특별한 자료형입니다. 슬라이싱(record.seq[:30])이나 len은 문자열처럼 그대로 쓰면서, 추가로 서열 연산까지 되는 셈입니다. SeqRecord는 이 Seq를 알맹이로 두고 식별자·설명을 껍질처럼 두른 구조라고 보면 됩니다.

3. Seq 객체 — 역상보·전사·번역

이제 Seq 객체의 진짜 쓸모를 봅니다. DNA 한 가닥만 알아도 상보 가닥, mRNA, 단백질 서열을 모두 끌어낼 수 있습니다. 중심원리(Central Dogma)의 DNA → RNA → 단백질 흐름이 그대로 메서드로 들어와 있습니다.

from Bio.Seq import Seq

# 짧은 코딩 서열로 시연 (ATG로 시작하는 가상의 CDS)
cds = Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG")

print("원본       :", cds)
print("상보 가닥   :", cds.complement())          # A↔T, G↔C 만 교체
print("역상보      :", cds.reverse_complement())  # 상보 + 뒤집기 (실제 반대편 가닥)
print("전사 mRNA   :", cds.transcribe())          # T → U (DNA → RNA)
print("번역 단백질 :", cds.translate())           # 코돈 3개 → 아미노산 1개

원본       : ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG
상보 가닥   : TACCGGTAACATTACCCGGCGACTTTCCCACGGGCTATC
역상보      : CTATCGGGCACCCTTTCAGCGGCCCATTACAATGGCCAT
전사 mRNA   : AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG
번역 단백질 : MAIVMGR*KGAR*

번역 결과 끝의 *는 종결코돈(stop codon)을 뜻합니다. 위 서열에는 중간(TGA)과 끝(TAG)에 종결코돈이 두 개 보이는데, 첫 종결코돈까지만 번역하고 멈추려면 to_stop=True를 줍니다. 또 세균(bacteria)이나 미토콘드리아는 표준 코돈표와 조금 다르므로, table 인자로 NCBI 유전암호 번호를 지정합니다(세균=11).

# 첫 종결코돈에서 멈추기
print(cds.translate(to_stop=True))            # MAIVMGR

# 세균 유전암호(표 11)로 번역
print(cds.translate(table=11, to_stop=True))  # MAIVMGR

# 길이가 3의 배수가 아니면 번역 시 경고가 뜸 → 미리 점검
read = Seq("ATGGCCATTGTAATGG")   # 16 bp (3의 배수 아님)
print(len(read) % 3)                          # 1 → 끝 1 bp가 남음
print(read[:len(read) - len(read) % 3].translate())  # 3의 배수로 잘라 번역

MAIVMGR
MAIVMGR
1
MAIVM

4. 서열 통계 — 길이·GC 함량·ORF

여러 서열을 한꺼번에 다룰 때는 SeqIO.parse로 반복문을 돕니다. 받아 온 lacZ 유전자(7,477 bp)로 기본 통계를 내 봅니다. GC 함량은 Bio.SeqUtils.gc_fraction이 0~1 사이 실수로 돌려주므로 100을 곱해 퍼센트로 바꿉니다.

from Bio import SeqIO
from Bio.SeqUtils import gc_fraction

record = SeqIO.read("lacz.fasta", "fasta")
seq = record.seq

length = len(seq)                      # 서열 길이 (bp)
gc = gc_fraction(seq) * 100            # GC 함량 (%) — 0~1 → 0~100
a = seq.count("A"); t = seq.count("T")
g = seq.count("G"); c = seq.count("C")

print(f"길이      : {length:,} bp")
print(f"GC 함량   : {gc:.1f} %")
print(f"A/T/G/C   : {a} / {t} / {g} / {c}")

길이      : 7,477 bp
GC 함량   : 53.7 %
A/T/G/C   : 1698 / 1771 / 2055 / 1953

다음은 ORF (open reading frame) 후보를 세어 봅니다. ORF는 개시코돈(ATG)에서 시작해 종결코돈(TAA·TAG·TGA)으로 끝나는 구간으로, 단백질을 코딩할 가능성이 있는 영역입니다. DNA는 양쪽 가닥에 각각 세 가지 번역틀이 있어 모두 여섯 개의 틀(6-frame)을 봐야 합니다. 여기서는 각 틀을 번역해 종결코돈으로 끊고, 길이가 일정 이상인 단백질 조각만 ORF 후보로 셉니다.

from Bio import SeqIO

record = SeqIO.read("lacz.fasta", "fasta")
seq = record.seq
min_aa = 100   # 아미노산 100개(=300 bp) 이상만 의미 있는 ORF 후보로

orf_count = 0
# 정방향(+)과 역방향(-) 두 가닥
for strand, nuc in [(+1, seq), (-1, seq.reverse_complement())]:
    # 한 가닥당 번역틀 3개 (시작 위치 0, 1, 2)
    for frame in range(3):
        # 길이를 3의 배수로 자른 뒤 번역
        trimmed = nuc[frame:]
        trimmed = trimmed[:len(trimmed) - len(trimmed) % 3]
        # 종결코돈(*)으로 쪼개 각 조각을 후보로
        for protein in trimmed.translate().split("*"):
            if len(protein) >= min_aa:
                orf_count += 1

print(f"6개 번역틀에서 {min_aa} aa 이상 ORF 후보: {orf_count}개")

6개 번역틀에서 100 aa 이상 ORF 후보: 7개

lacZ 영역에는 실제로 여러 유전자(lacI·lacZ·lacY·lacA)가 들어 있어, 100 아미노산 이상 ORF 후보가 일곱 개 잡힙니다. 이 방식은 개시코돈을 엄밀히 따지지 않은 단순 카운트라 진짜 유전자 수와는 다를 수 있습니다. 정밀한 유전자 예측에는 Prodigal 같은 전용 도구를 쓰지만, "이 서열에 코딩 영역이 얼마나 있을 법한가"를 빠르게 가늠하는 데는 이 정도로 충분합니다.

5. FASTQ 품질점수 — Phred score 접근

FASTQ는 read마다 염기별 품질점수를 함께 담습니다. 이 점수가 Phred 점수(Q)인데, 염기 판정이 틀렸을 확률 P 와 Q = −10 × log₁₀ P 의 관계로 정의됩니다. Q = 20이면 오류율 1%, Q = 30이면 0.1%입니다. 파일에는 점수가 ASCII 문자로 인코딩돼 있지만, Biopython이 자동으로 풀어 letter_annotations["phred_quality"]에 숫자 리스트로 넣어 줍니다.

먼저 시연용 FASTQ를 직접 한 건 만들어 보겠습니다(실전에서는 시퀀서나 SRA에서 받은 파일을 그대로 엽니다).

# 시연용 FASTQ 한 건 작성 — 식별자/서열/+/품질 4줄
fastq_text = """@read1
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCC
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55
"""
with open("demo.fastq", "w") as f:
    f.write(fastq_text)

이제 SeqIO.parse로 읽으면, 각 read의 서열과 품질점수가 함께 딸려 옵니다.

from Bio import SeqIO

for record in SeqIO.parse("demo.fastq", "fastq"):
    quals = record.letter_annotations["phred_quality"]   # 정수 리스트
    avg = sum(quals) / len(quals)
    print("read id   :", record.id)
    print("서열      :", record.seq)
    print("품질 앞 8 :", quals[:8])
    print(f"평균 품질 : Q{avg:.1f}")

read id   : read1
서열      : GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCC
품질 앞 8 : [0, 6, 6, 9, 7, 7, 7, 7]
평균 품질 : Q21.6

품질점수가 ASCII 문자(!''*(...)에서 숫자 리스트로 자동 변환된 게 핵심입니다. !는 Q = 0, '는 Q = 6에 해당합니다(Sanger/Illumina는 ASCII 값에서 33을 빼는 Phred+33 방식). 이 숫자만 있으면 read별 평균 품질, 위치별 품질 곡선, 품질 필터링 같은 QC 작업을 전부 직접 짤 수 있습니다. 대용량 파일이라면 보통 FastQC 같은 도구가 이 계산을 대신해 주지만, 원리는 정확히 이것입니다.

6. 자주 발생하는 에러 & 해결 (FAQ)

• ImportError: cannot import name 'GC' — GC()는 옛 함수라 최신 버전에서 빠졌습니다. from Bio.SeqUtils import gc_fraction으로 바꾸고, 퍼센트가 필요하면 결과에 * 100을 곱하세요(gc_fraction은 0~1 비율 반환).
• ModuleNotFoundError: No module named 'Bio.Alphabet' — 1.78 버전부터 Bio.Alphabet이 제거됐습니다. 옛 코드의 from Bio.Alphabet import ... 줄과 Seq(..., IUPAC.xxx)의 알파벳 인자를 모두 지우세요. 지금은 Seq("ATGC")만으로 충분합니다.
• ValueError: Sequence length is not a multiple of three — translate()에 길이가 3의 배수가 아닌 서열을 주면 경고가 납니다. 코딩 구간만 잘라 seq[start:end]로 3의 배수를 맞추거나, ORF 예제처럼 len % 3만큼 끝을 잘라내세요.
• SeqIO.read()가 에러를 냄 — read는 파일에 서열이 정확히 한 개일 때만 씁니다. 여러 서열이 든 멀티-FASTA에서는 parse로 반복문을 돌려야 합니다. 반대로 parse 결과를 len()으로 세려 하면 제너레이터라 안 되니, list(SeqIO.parse(...))로 감싸세요.
• 한글·특수문자로 파일 읽기 실패(UnicodeDecodeError) — FASTA·FASTQ는 ASCII가 원칙이지만, 설명 줄에 깨진 문자가 섞이면 읽기가 멈춥니다. open(path, encoding="utf-8", errors="replace")로 핸들을 열어 SeqIO.parse에 넘기면 우회됩니다. 윈도우에서 만든 파일이 \r\n 줄바꿈을 써도 Biopython은 알아서 처리합니다.
• NCBI efetch가 비거나 멈춤 — Entrez.email을 안 적으면 차단될 수 있습니다. 또 짧은 시간에 너무 많이 요청하면 NCBI가 막으므로, 대량 다운로드는 Entrez.api_key를 발급받아 쓰거나 초당 3회 이하로 제한하세요.

7. 확장 — pandas로 넘기기

서열별 통계를 냈다면 다음은 표로 정리해 비교·시각화하는 단계입니다. SeqIO.parse로 여러 서열을 돌면서 딕셔너리 리스트를 만들면, pandas가 그걸 데이터프레임으로 바로 받습니다.

import pandas as pd
from Bio import SeqIO
from Bio.SeqUtils import gc_fraction

# 멀티-FASTA의 모든 서열을 한 줄씩 표로
df = pd.DataFrame([
    {"id": r.id, "length": len(r.seq), "gc_pct": round(gc_fraction(r.seq) * 100, 1)}
    for r in SeqIO.parse("lacz.fasta", "fasta")
])
print(df)

          id  length  gc_pct
0  J01636.1    7477    53.7

이렇게 데이터프레임으로 옮기고 나면 df.describe()로 요약 통계를 보거나, df.plot.scatter("length", "gc_pct")로 길이-GC 산점도를 그리는 식으로 넘어갈 수 있습니다. 서열 수백 개의 길이 분포나 GC 편향을 한눈에 보는 작업이 여기서부터 시작됩니다. Biopython이 서열을 숫자로 바꿔 주면, 그다음은 일반적인 데이터 분석 도구의 영역입니다.

References

1. Cock, P. J. A., Antao, T., Chang, J. T. et al. (2009). Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics, 25(11), 1422–1423. (Biopython 원논문)
2. Biopython Tutorial and Cookbook (공식 문서) — Bio.SeqIO·Seq·SeqRecord 전체 가이드.
3. Bio.SeqUtils API — gc_fraction — GC 함량 계산 함수 레퍼런스.
4. Bio.SeqIO.QualityIO — FASTQ 품질점수 처리 — Phred+33 인코딩 구현.
5. Biopython NEWS (버전별 변경점) — 1.86 릴리스 노트·Bio.Alphabet 제거 이력.
6. Ewing, B., & Green, P. (1998). Base-calling of automated sequencer traces using phred. II. Error probabilities. Genome Research, 8(3), 186–194. (Phred 품질점수 정의)

Pipette & Pipeline · A bio portfolio journal