29.4 유전체·단일세포 데이터 · AI for Science 로드맵

단일세포 RNA 시퀀싱(single-cell RNA-seq)은 세포 하나하나에서 수천 개 유전자의 RNA 분자 수를 센다. 결과는 세포 × 유전자의 거대한 개수 행렬이다. 이 행렬은 ML 연구자에게 익숙한 "고차원 표 데이터"처럼 보이지만, 세 가지 점에서 다르다. 값이 유한한 표본에서 센 정수라 0이 매우 많고, 세포마다 측정 깊이가 달라 값의 크기를 그대로 비교할 수 없으며, 실험 배치마다 체계적인 차이가 있다.

이 장은 PCA·군집·시각화(16.5장)를 다시 가르치지 않는다. 대신 그 도구를 개수 데이터에 적용하기 전에 해야 하는 일(정규화, 변환, 유전자 선택)과, 적용한 결과를 해석할 때 빠지기 쉬운 함정(깊이 차이가 만든 축, 배치가 만든 군집, 보정이 지운 생물학)을 다룬다. 측정의 기초(시퀀싱 깊이, 기술적·생물학적 반복, 배치 효과)는 29.1장에서 보았다.

학습 목표

고차원 개수(count) 데이터를 전처리한다
배치 효과와 생물학적 신호를 구분한다

데이터의 특성

왜 필요한가. 전처리의 모든 단계는 데이터가 어떻게 생성되었는가에 대한 가정이다. 생성 과정을 모르면 어떤 단계를 왜 하는지, 언제 그 가정이 깨지는지 판단할 수 없다.

개수 데이터. 세포 $i$ 에서 유전자 $g$ 의 RNA 분자를 센 값을 $x_{ig} \in \{0, 1, 2, \dots\}$ 라 한다. 세포 하나에서 포획해 읽는 분자는 그 세포에 있던 RNA의 작은 일부다. 그래서 $x_{ig}$ 는 세포의 실제 발현 수준 $\mu_{ig}$ 에 비례하는 평균을 가진 확률변수로 보는 것이 자연스럽다.

x_{ig} \sim \mathrm{Poisson}(s_i \, \mu_{ig}) \quad \text{또는} \quad x_{ig} \sim \mathrm{NB}\big(s_i \, \mu_{ig},\ \theta\big).

여기서 $s_i$ 는 세포마다 다른 포획 효율·시퀀싱 깊이를 나타내는 크기 인자(size factor)다. 같은 상태의 세포라도 $s_i$ 가 다르면 모든 유전자의 개수가 함께 커지거나 작아진다. 음이항 분포(negative binomial, NB)는 Poisson보다 분산이 큰 개수 분포로, 평균이 $m$ 일 때 분산이 $m + m^2/\theta$ 다(12.6장). 세포 사이의 생물학적 변동이 Poisson 표본 추출 위에 더해지면 이런 과분산이 생긴다.

희소성. 평균이 $m$ 인 Poisson 변수가 0일 확률은 $e^{-m}$ 이다. 평균 0.1개인 유전자는 세포의 약 90%에서 0으로 관측된다. 단일세포 데이터의 0 대부분은 "그 유전자가 꺼져 있다"가 아니라 "표본에 잡히지 않았다"로 설명된다. 0을 별도의 "탈락"(dropout) 과정으로 모델링해야 하는지는 측정 기술과 데이터에 따라 다르므로, 먼저 개수 분포가 예측하는 0의 비율과 관측을 비교해 본다.

세 가지 세포 유형으로 된 합성 데이터를 만든다. 세포 600개, 유전자 500개이고, 유형마다 marker 유전자 20개의 발현이 8배 높다. 세포마다 크기 인자는 로그정규분포를 따르고, 개수는 음이항 분포( $\theta = 2$ )에서 뽑는다. 유형 2에만 반응하는 유전자 20개(100–119번)도 정의해 두고, 배치 효과 절에서 쓴다.

import numpy as np

rng = np.random.default_rng(0)
G, THETA = 500, 2.0                                    # 유전자 수, 음이항 분포의 모양 매개변수

def simulate(n_per_type, response=1.0):
    types = np.repeat(np.arange(3), n_per_type)        # (N,) 세포 유형 0, 1, 2
    mu = np.tile(base, (3, 1))                         # (3, 500) 유형별 평균 발현
    for t in range(3):
        mu[t, 20 * t:20 * t + 20] *= 8                 # 유형마다 marker 유전자 20개
    mu[2, 100:120] *= response                         # 유형 2만 반응하는 유전자(예: 처리 조건)
    size = rng.lognormal(0.0, 0.6, len(types))         # (N,) 세포마다 포획 효율
    lam = size[:, None] * mu[types]                    # (N, 500)
    return rng.poisson(rng.gamma(THETA, lam / THETA)), types   # 감마–Poisson = 음이항

base = rng.lognormal(0.0, 1.5, G)                     # (500,) 유전자별 기본 발현 수준
counts, types = simulate([200, 200, 200])              # (600, 500)
lib = counts.sum(1)                                    # (600,) 세포별 총 개수
lam_hat = lib[:, None] * counts.sum(0)[None, :] / counts.sum()   # (600, 500) 단순 기대값
print(f"0의 비율 {np.mean(counts == 0):.2f}")
print(f"세포별 총 개수: 최소 {lib.min()}, 중앙값 {int(np.median(lib))}, 최대 {lib.max()}")
print(f"Poisson이 예측하는 0의 비율 {np.mean(np.exp(-lam_hat)):.2f}")

0의 비율 0.45
세포별 총 개수: 최소 267, 중앙값 1691, 최대 12306
Poisson이 예측하는 0의 비율 0.39

0이 행렬의 45%다. 세포별 총 개수는 최소와 최대가 40배 넘게 차이 난다. 생물학이 아니라 크기 인자 때문이다. 세포·유전자별 기대값으로 계산한 Poisson의 0 비율(0.39)은 관측(0.45)보다 작다. 차이는 음이항 분포의 과분산에서 온다. 0의 대부분은 개수 분포만으로 설명되고, 나머지도 과분산으로 설명된다. 이 데이터에는 별도의 탈락 과정이 없다.

정규화. 크기 인자의 차이를 없애는 가장 단순한 방법은 세포마다 총 개수로 나누고 공통의 값(흔히 $10^4$ )을 곱하는 것이다. 이어서 $\log(1 + x)$ 로 변환한다.

\tilde{x}_{ig} = \log\!\Big(1 + 10^4 \cdot \frac{x_{ig}}{\sum_{h} x_{ih}}\Big).

로그는 곱셈적인 차이(2배, 8배)를 덧셈적인 차이로 바꾸고, 큰 값에 끌려가는 분산을 줄인다. 1을 더하는 이유는 0 때문이다. $\log 0$ 은 정의되지 않는다. 1을 더하는 것은 편의를 위한 선택이라, 작은 값 근처에서는 변환의 모양이 이 상수에 따라 달라진다는 점도 기억한다.

고변이 유전자. 유전자 수천 개 중 대부분은 세포 유형과 무관하게 비슷하거나, 너무 적게 발현되어 잡음뿐이다. 그래서 평균에 비해 분산이 큰 유전자(고변이 유전자, highly variable genes, HVG)를 수백–수천 개 골라 이후 분석에 쓴다. 평균 대비 분산을 쓰는 이유는 개수 데이터에서 분산이 평균과 함께 커지기 때문이다. 분산만으로 고르면 많이 발현된 유전자만 뽑힌다.

차원 축소와 군집

왜 필요한가. 정규화한 행렬도 수백 차원이다. 세포 유형이나 상태를 찾으려면 차원을 줄이고 비슷한 세포를 묶는다. 도구는 16.5장의 PCA, 군집, 시각화 그대로다. 달라지는 것은 입력이 무엇이냐다.

PCA는 정규화된 데이터에. 원 개수에 PCA를 적용하면 분산이 가장 큰 방향, 즉 크기 인자의 방향이 첫 주성분이 된다. 정규화·log1p·HVG를 거친 뒤의 PCA와 비교한다. 두 경우 모두 PCA 10차원 위에서 k-means로 세 군집을 찾고, 군집과 실제 세포 유형의 일치를 adjusted Rand index(ARI, 1이면 완전 일치, 0이면 우연 수준)로 잰다.

import warnings
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA
from sklearn.metrics import adjusted_rand_score as ari
warnings.filterwarnings("ignore", module="threadpoolctl")      # 브라우저 실행 환경의 무해한 경고만 끈다

def normalize(C):                                      # 세포마다 총합 1만으로 맞춘 뒤 log1p
    return np.log1p(C / C.sum(1, keepdims=True) * 1e4)          # (N, G)

def hvg(X, n=100):                                     # 평균 대비 분산이 큰 유전자 n개의 번호
    m, v = X.mean(0), X.var(0)
    return np.argsort(v / (m + 1e-3))[-n:]                       # (n,)

def embed_and_cluster(M, k=3):
    Z = PCA(10, random_state=0).fit_transform(M - M.mean(0))     # (N, 10)
    return Z, KMeans(k, n_init=10, random_state=0).fit_predict(Z)

norm = normalize(counts)                               # (600, 500)
for name, M in [("원 개수", counts.astype(float)), ("정규화+log1p+HVG", norm[:, hvg(norm)])]:
    Z, lab = embed_and_cluster(M)
    r = np.corrcoef(np.abs(Z[:, 0]), lib)[0, 1]
    print(f"{name}: |PC1|과 총 개수의 상관 {r:.2f}, 군집과 세포 유형의 ARI {ari(types, lab):.2f}")

원 개수: |PC1|과 총 개수의 상관 0.58, 군집과 세포 유형의 ARI 0.02
정규화+log1p+HVG: |PC1|과 총 개수의 상관 -0.03, 군집과 세포 유형의 ARI 0.83

원 개수의 첫 주성분은 세포별 총 개수와 0.58의 상관을 보이고, 군집은 세포 유형과 거의 무관하다(ARI 0.02). 군집이 찾은 것은 "많이 읽힌 세포와 적게 읽힌 세포"다. 정규화 뒤에는 첫 주성분과 총 개수의 상관이 사라지고 군집이 세포 유형을 대부분 되찾는다(ARI 0.83). 실제 분석에서도 군집이나 주성분이 총 개수, 검출된 유전자 수, 미토콘드리아 유전자 비율 같은 품질 지표와 상관을 보이면, 생물학이 아니라 측정을 본 것은 아닌지 먼저 의심한다.

그래프 기반 군집. 실제 단일세포 분석에서는 k-means보다 그래프 기반 군집을 흔히 쓴다. PCA 공간에서 세포마다 k-최근접 이웃 그래프를 만들고, 그래프에서 촘촘하게 연결된 공동체를 찾는다(Louvain, Leiden 같은 공동체 탐지 알고리즘). 군집 수를 미리 정하지 않아도 되지만, 해상도 매개변수가 군집의 수를 정한다. 군집 수는 데이터가 아니라 분석자가 고른 값이라는 점은 같다(16.5장).

시각화의 해석 한계. t-SNE와 UMAP 그림은 이웃 관계를 보여 주도록 설계되었다. 그림에서 군집 사이의 거리, 군집의 크기와 밀도는 원래 공간의 거리·크기·밀도를 반영한다는 보장이 없다(16.5장). 그림에서 떨어져 보이는 섬 하나를 새로운 세포 유형으로 보고하기 전에, PCA 공간의 군집 결과, marker 유전자의 발현, 그리고 다음 절의 배치 정보를 함께 확인한다.

배치 효과

왜 필요한가. 단일세포 연구는 여러 시료, 여러 날, 여러 기관의 데이터를 합친다. 배치마다 시약, 장비, 처리 시간이 다르고, 이 차이가 모든 세포에 체계적으로 들어간다(29.1장). 데이터를 합치면 세포가 유형이 아니라 배치로 묶인다.

보정과 과보정. 배치 보정(batch correction) 또는 데이터 통합(integration)은 배치의 차이를 제거해 같은 유형의 세포가 배치와 무관하게 함께 묶이게 하는 작업이다. 문제는 배치 차이 안에 생물학이 섞여 있을 수 있다는 점이다. 두 배치의 세포 구성이 다르거나, 한 배치에만 있는 조건(처리, 질병)이 있으면, 배치의 평균 차이는 배치 효과와 생물학적 차이의 합이다. 이를 통째로 제거하면 생물학적 차이까지 지운다. 이것을 과보정(over-correction)이라 한다.

배치 2를 만든다. 세포 구성이 다르고(유형 0이 적다), 유형 2 세포는 처리에 반응해 반응 유전자 20개의 발현이 4배 높다. 여기에 로그 공간에서 유전자마다 더해지는 배치 효과를 넣는다. 실제 배치 효과는 이보다 복잡하지만, 보정의 원리를 보기에는 충분한 단순화다. 세 가지를 비교한다. 보정하지 않은 데이터, 배치마다 전체 평균을 빼는 보정, 그리고 두 배치에 모두 있고 처리에 반응하지 않는 유형 1 세포의 평균으로만 배치 차이를 추정해 빼는 보정이다. 마지막 방법은 세포 유형을 marker로 확인할 수 있다고 가정한다. 비교 기준은 배치 효과를 넣기 전의 참값이다.

counts2, types2 = simulate([40, 200, 160], response=4.0)   # 배치 2: 구성이 다르고 유형 2가 반응
T = np.concatenate([types, types2])                    # (1000,) 세포 유형
B = np.repeat([0, 1], [600, 400])                      # (1000,) 배치
truth = normalize(np.vstack([counts, counts2]))        # (1000, 500) 배치 효과가 없을 때
M = truth.copy()
M[B == 1] += rng.normal(0, 0.8, G)                     # 로그 공간에서 유전자마다 더해지는 배치 효과
genes = hvg(M)                                         # (100,) 군집에 쓸 유전자는 한 번만 고른다

def center(M, ref):                                    # 배치마다 ref 세포의 평균을 뺀다
    out = M.copy()
    for b in (0, 1):
        out[B == b] -= M[(B == b) & ref].mean(0)
    return out

def response(X, t):                                    # 유형 t에서 반응 유전자 점수의 배치 간 차이
    r = X[:, 100:120].mean(1)
    return r[(T == t) & (B == 1)].mean() - r[(T == t) & (B == 0)].mean()

for name, X in [("참값(배치 효과 없음)", truth), ("보정 전", M),
                ("배치 평균 빼기", center(M, np.ones(len(B), bool))),
                ("공유 유형 1로 보정", center(M, T == 1))]:
    _, lab = embed_and_cluster(X[:, genes])
    print(f"{name}: ARI(유형) {ari(T, lab):.2f} ARI(배치) {ari(B, lab):.2f} | "
          f"반응 차이 유형0 {response(X, 0):+.2f} 유형2 {response(X, 2):+.2f}")

참값(배치 효과 없음): ARI(유형) 0.86 ARI(배치) 0.03 | 반응 차이 유형0 -0.07 유형2 +1.23
보정 전: ARI(유형) 0.25 ARI(배치) 0.65 | 반응 차이 유형0 +0.16 유형2 +1.46
배치 평균 빼기: ARI(유형) 0.84 ARI(배치) 0.01 | 반응 차이 유형0 -0.54 유형2 +0.76
공유 유형 1로 보정: ARI(유형) 0.87 ARI(배치) 0.03 | 반응 차이 유형0 -0.07 유형2 +1.23

보정하지 않으면 군집이 배치를 따라간다(ARI(배치) 0.65). 배치 평균을 빼면 군집은 세포 유형과 잘 맞고(ARI 0.84) 배치와는 무관해진다. 군집 지표만 보면 성공이다. 그러나 반응 유전자를 보면 문제가 드러난다. 참값에서 유형 2의 반응은 +1.23이고 유형 0은 반응이 없다(−0.07). 배치 평균을 빼자 유형 2의 반응은 +0.76으로 줄었고, 반응하지 않은 유형 0에 −0.54의 가짜 반응이 생겼다. 배치 2의 평균에는 유형 2의 반응이 섞여 있었고, 그것을 모든 세포에서 빼 버렸기 때문이다. 공유 유형으로 추정한 보정은 군집도 맞추고 반응도 참값 그대로 남긴다.

교훈은 세 가지다.

보정의 평가는 두 방향이다. 배치가 섞였는가(배치 혼합)와, 생물학적 차이가 남았는가(생물학 보존)를 함께 본다. 이 장의 실험처럼 군집 지표는 과보정을 드러내지 못할 수 있다.
보정은 무엇을 같다고 가정하느냐에 달려 있다. 배치 평균 빼기는 "두 배치의 구성과 조건이 같다"고 가정한다. 공유 세포 집단을 기준으로 하는 방법은 "그 집단은 배치 효과 외에는 같다"고 가정한다. 실제로 쓰이는 통합 방법들은 대개 배치 사이에서 서로 대응하는 세포를 찾아 그 차이로 보정하는데, 대응을 잘못 찾으면 서로 다른 유형을 하나로 합친다.
조건과 배치가 완전히 겹치면 보정할 수 없다. 처리군을 모두 한 배치에서 측정했다면, 어떤 보정도 처리 효과와 배치 효과를 나누지 못한다(29.1장). 이 경우 처리 효과에 대한 결론은 보정 전후 어느 쪽에서도 내릴 수 없다.

모델링

왜 필요한가. 정규화와 log1p는 편리하지만 임시방편이다. 크기 인자, 과분산, 배치를 확률 모델 안에 넣으면 각 요소의 역할이 분명해지고, 불확실성도 함께 다룰 수 있다.

확률 모델(개수 분포). 개수 자체를 음이항 분포로 모델링하고, 크기 인자를 평균에 곱하는 오프셋(offset)으로, 배치와 조건을 평균의 설명변수로 넣는다(12.6장의 일반화 선형 모델).

x_{ig} \sim \mathrm{NB}(m_{ig}, \theta_g), \qquad \log m_{ig} = \log s_i + \beta_{g,0} + \beta_{g,\text{batch}(i)} + \beta_{g,\text{cond}(i)}.

이 형태에서는 배치 효과와 조건 효과를 같은 식에서 동시에 추정하므로, 설계가 허락하는 한 둘을 분리할 수 있고, 겹쳐 있으면 추정 불가능성이 그대로 드러난다. 과분산 $\theta$ 는 무엇을 담는가. 한 세포 유형의 세포만 모아 평균–분산 관계 $v \approx m + m^2/\theta$ 를 맞춰 본다.

X0 = counts[types == 0]                                # (200, 500) 한 유형의 세포
m, v = X0.mean(0), X0.var(0)
keep = m > 0.5
inv_theta = np.sum((v - m)[keep] * m[keep] ** 2) / np.sum(m[keep] ** 4)   # v ≈ m + m²/θ 최소제곱
print(f"분산/평균 중앙값 {np.median(v[keep] / m[keep]):.1f}  (Poisson이면 1)")
print(f"추정 θ {1 / inv_theta:.2f}  (시뮬레이션의 θ = {THETA})")

분산/평균 중앙값 3.1  (Poisson이면 1)
추정 θ 1.20  (시뮬레이션의 θ = 2.0)

분산은 평균의 3.1배로 Poisson보다 훨씬 크다. 추정한 $\theta$ 는 1.20으로 시뮬레이션의 2.0보다 작다(과분산이 더 크다). 크기 인자 $s_i$ 의 세포 간 변동도 분산을 키우기 때문이다. 크기 인자를 모델에 넣지 않으면 그 변동이 과분산으로 흡수된다. 과분산은 생물학적 변동, 기술적 변동, 모델에서 빠진 요인이 합쳐진 값이다.

표현학습. 개수 분포를 likelihood로 쓰는 VAE(19.3장)로 세포마다 저차원 잠재변수를 학습하는 접근도 널리 쓰인다. 디코더가 음이항 분포의 매개변수를 내고, 크기 인자와 배치를 조건으로 넣어 잠재변수가 배치와 무관한 세포 상태를 담도록 유도한다. 이 경우에도 배치 정보를 넣는 방식이 과보정 여부를 정하므로, 앞 절의 두 방향 평가가 그대로 필요하다.

세포 상태 예측. 세포의 상태(유형, 처리 반응, 분화 단계)를 예측하는 지도학습도 같은 데이터 위에서 한다. 이때 독립성의 단위는 세포가 아니라 시료(개체, 기증자)다. 같은 기증자의 세포는 유전적 배경과 처리 과정을 공유하므로, 세포를 무작위로 나누면 기증자 정보가 누수된다(15.4장). "새 기증자에게 일반화한다"는 주장은 기증자 단위로 나눈 평가로만 뒷받침된다.

실습

두 배치의 단일세포 데이터를 통합하고, 보정 전후의 군집과 알려진 세포 유형의 대응을 비교한다. 실제 데이터 분석에는 Scanpy(Python)나 Seurat(R) 같은 단일세포 분석 도구를 쓸 수 있다(집필 시점 예시. 이 책의 실습 환경에는 들어 있지 않으므로 각 도구의 설치 안내를 확인한다). 단계별 권장 사항은 Luecken·Theis(2019)의 튜토리얼이 정리했다.

세포 유형 주석이 있는 공개 데이터에서 같은 조직을 다른 배치(기관, 기술, 날짜)로 측정한 두 데이터셋을 고른다. 데이터의 판과 배치 정보를 기록한다.
품질 관리(총 개수·검출 유전자 수가 극단적인 세포 제거)를 하고, 기준과 제거된 세포 수를 기록한다.
정규화 → log1p → 고변이 유전자 → PCA → 군집을 보정 없이 수행하고, 군집을 세포 유형과 배치에 각각 대응시킨 ARI를 계산한다.
배치 보정을 하나 이상 적용하고 같은 지표를 계산한다.
생물학 보존을 확인한다. 두 배치에 모두 있는 세포 유형마다 알려진 marker 유전자의 발현 차이가 보정 전후에 유지되는지, 한 배치에만 있는 세포 유형이 다른 유형에 합쳐지지 않았는지 본다.
결과를 "배치 혼합"과 "생물학 보존"의 두 열로 된 표로 보고한다.

정규화 함수를 직접 구현한다. 총 개수가 0인 세포처럼 실제 데이터에서 생기는 예외도 처리한다.

실패 사례

배치 효과를 생물학적 차이로 해석하는 것. 군집 하나가 거의 한 배치의 세포로만 이루어져 있으면, 새로운 세포 유형이라고 보고하기 전에 배치를 의심한다.
원 개수에 PCA를 적용하는 것. 첫 주성분이 크기 인자가 되고, 군집이 시퀀싱 깊이를 따라간다.
0을 "발현되지 않음"으로 해석하는 것. 평균이 작은 유전자는 발현되어 있어도 대부분의 세포에서 0으로 관측된다.
보정 성공을 배치 혼합만으로 판단하는 것. 배치 평균을 빼는 보정은 군집 지표를 좋게 만들면서 조건의 효과를 줄이고 다른 세포 유형에 가짜 효과를 만들 수 있다.
세포를 독립 표본으로 취급하는 것. 같은 기증자의 세포 수천 개를 독립 표본으로 세면 통계 검정의 p-값과 예측 모델의 성능이 모두 부풀려진다.

통과 기준

두 배치 데이터에 대해 보정 전후의 군집–세포 유형 ARI와 군집–배치 ARI를 보고한다.
두 배치에 공통인 세포 유형에서 marker 유전자 발현 차이가 보정 후에도 유지되는지 수치로 확인하고, 과보정 여부를 근거와 함께 쓴다.
연습문제 normalize-counts의 테스트를 통과한다.

참고 문헌

Luecken, Theis (2019). Current best practices in single-cell RNA-seq analysis: a tutorial. Molecular Systems Biology, 15(6), e8746. — 품질 관리·정규화·배치 보정·특징 선택·차원 축소의 단계별 권장 사항
Leek 외 (2010). Tackling the widespread and critical impact of batch effects in high-throughput data. Nature Reviews Genetics, 11(10), 733–739. — 배치 효과가 결과 변수와 겹칠 때의 위험