numpy·scipy로 생물 데이터 수치·통계 — 배열 연산·가설검정
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 scipyimport 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 = 10GeneD는 처리군 평균(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(또는 statsmodelsmultipletests)로 보정한 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 편의 분류·교차검증으로 넘기면, "어떤 유전자 조합이 군을 가르는가"까지 이어집니다.
• (데이터 정제) pandas로 생물 데이터 정제하기 — 검정에 넣을 발현 행렬을 분석 직전까지 다듬는 전 단계
• (R로 같은 검정) R로 기초 통계 검정하기 — t검정·상관·비모수 검정을 R로, 이 글과 1:1 비교
• (개념 토대) p-value와 FDR, 다중검정보정이란? — 왜 adjusted p값을 쓰는지, FDR의 직관
• (머신러닝) scikit-learn으로 생물 데이터 머신러닝 — 여기서 만든 배열을 특징으로 삼아 분류·교차검증
References
- Harris, C. R. et al. (2020). Array programming with NumPy. Nature, 585, 357–362. (numpy 원논문)
- Virtanen, P. et al. (2020). SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17, 261–272. (scipy 원논문)
- NumPy developers. (2026). NumPy v2.5 documentation — Broadcasting · Indexing. numpy.org/broadcasting
- SciPy developers. (2026). scipy.stats.false_discovery_control — Benjamini-Hochberg·Benjamini-Yekutieli FDR 보정(scipy 1.11+). false_discovery_control docs
- SciPy developers. (2026). scipy.stats.ttest_ind · mannwhitneyu — 독립표본 검정 API. ttest_ind docs
- 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
'Code Bench > Python for Bio' 카테고리의 다른 글
| gseapy로 GO·GSEA 농축분석 — 파이썬으로 유전자셋 분석 (0) | 2026.07.09 |
|---|---|
| pydeseq2로 RNA-seq 차등발현 — R 없이 파이썬만으로 DESeq2 (카운트→볼케이노) (1) | 2026.07.07 |
| scikit-learn으로 생물 데이터 머신러닝 — 분류·교차검증·ROC (0) | 2026.06.29 |
| matplotlib·seaborn으로 생물 데이터 시각화 — 발현 데이터 플롯 실전 (0) | 2026.06.25 |
| scanpy로 단일세포 RNA-seq 분석하기 — QC부터 UMAP·클러스터까지 (Python for Bio) (0) | 2026.06.23 |
