16강 / 20강

긴 유전자가 read 를 더 받는다 — 그런데 나누면 더 나빠졌다

TPM · FPKM 이 무엇을 고치는가 · 그리고 언제 쓰면 안 되는가

mRNA 를 잘게 부수어 76bp 씩 읽습니다. 그러면 긴 mRNA 에서 조각이 더 많이 나옵니다.

유전자 A (1,000 염기) 분자 10개 → 조각 약 130개 유전자 B (4,000 염기) 분자 10개 → 조각 약 530개

같은 개수의 분자인데 카운트는 4배 차이가 납니다. 그래서 한 표본 안에서 유전자끼리 견주려면 길이로 나눠야 합니다. 그것이 TPM(transcripts per million)이 하는 일입니다.

RPK = count ÷ 길이(kb) → TPM = RPK ÷ ΣRPK × 10⁶

이 자료에서도 그럴까요? 재어 봅시다.

고등학생긴 유전자가 read 를 더 받는다

같은 양으로 발현되는 유전자 두 개가 있다. 하나는 1,000 글자, 하나는 5,000 글자. 어느 쪽이 read 를 더 많이 받을까?

긴 쪽이 다섯 배 받는다

RNA-seq 은 전사체를 조각내서 읽는다. 긴 전사체는 조각이 더 많이 나온다. 분자 수가 같아도 그렇다.

같은 길이의 실 열 가닥과, 다섯 배 긴 실 열 가닥을 각각 잘게 자르면 후자가 다섯 배 많은 토막을 낸다. 실의 개수는 같은데.

	분자 수	길이	read 수
유전자 A	100	1,000	기준
유전자 B	100	5,000	5배

read 수는 발현량 × 길이 에 비례한다 — 길이가 섞여 들어간다

기대는 양의 상관이었습니다. 자료는 반대로 말합니다.

Python 브라우저에서 실행

import bu
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

T = bu.table()
panel = bu.panel()
L = np.array([len(panel[g]) for g in T.index])
m = T["WT_1"].values.astype(float)

sel = m > 0
r = stats.pearsonr(np.log10(L[sel]), np.log2(m[sel]))[0]

fig, ax = plt.subplots(figsize=(6.2, 4.4))
ax.scatter(np.log10(L[sel]), np.log2(m[sel]), s=14, alpha=.5, c="#4c6ef5")
ax.set_xlabel("log10 전사체 길이"); ax.set_ylabel("log2 내 카운트")
ax.set_title(f"길이 대 카운트   r = {r:+.3f}")
plt.tight_layout()

print(f"전사체 길이 중앙값 : {int(np.median(L)):,} 염기")
print(f"길이와 카운트의 상관 : {r:+.3f}")
print()
print("양수가 나와야 할 것 같은데 음수입니다. 왜일까요?")

기대와 반대다

길이와 카운트의 상관이 −0.38 입니다. 긴 유전자가 덜 받았다는 뜻이죠. 길이 편향이 없어진 걸까요? 아닙니다. 더 센 것이 덮은 것입니다.

우리 패널의 최상위는 리보솜 단백질 유전자(Rpl·Rps)들인데, 이들은 짧고(500~1,000 염기) 엄청나게 많이 발현됩니다.

짐작으로 끝내지 말고 확인합시다 — 그 52개를 빼고 다시 재 보면 됩니다.

고등학생TPM 과 FPKM 이 고치려는 것

길이 편향을 없애려면 read 수를 길이로 나누면 된다. 이 발상에서 나온 것이 RPKM · FPKM · TPM 이다.

지표	계산 순서
FPKM	깊이로 먼저 나누고 → 길이로 나눔
TPM	길이로 먼저 나누고 → 합이 100만이 되게 맞춤

TPM 이 나은 이유

TPM 은 모든 시료에서 합이 정확히 100만 이다. 그래서 시료 사이 비교가 말이 된다. FPKM 은 합이 시료마다 달라서, 같은 FPKM 값이 시료마다 다른 뜻을 갖는다.

비율을 비교하려면 분모가 같아야 한다 는 당연한 원칙인데, FPKM 이 널리 퍼진 뒤에야 지적되었다. 지금도 옛 논문에는 FPKM 이 많다.

TPM 을 쓰라 — 분모가 고정되어 있기 때문이다

짐작을 확인하는 가장 싼 방법 — 의심되는 것을 빼고 다시 재기.

Python 브라우저에서 실행

import bu, re
import numpy as np
from scipy import stats

T = bu.table(); panel = bu.panel()
cpm = bu.cpm(bu.published())
L = np.array([len(panel[g]) for g in T.index], dtype=float)
m = T["WT_1"].values.astype(float)
p = cpm["WT_DDX54_1"].reindex(T.index).fillna(0).values
ribo = np.array([bool(re.match(r"Rp[ls]", g)) for g in T.index])

def lenr(v, sel):
    return stats.pearsonr(np.log10(L[sel]), np.log2(v[sel] + 1))[0]

print(f"리보솜 단백질 유전자 : {ribo.sum()}개")
print()
print(f"{'':<16}{'전체':>10}{'리보솜 뺀 것':>14}")
print("-" * 42)
print(f"{'내 카운트':<16}{lenr(m, m>0):>10.3f}{lenr(m, (m>0)&~ribo):>14.3f}")
print(f"{'공개 CPM':<16}{lenr(p, p>0):>10.3f}{lenr(p, (p>0)&~ribo):>14.3f}")
print()
print("52개를 빼니 −0.38 이 −0.15 로 줄었습니다 — 짐작이 맞았습니다.")
print("그리고 공개 표도 같은 방향입니다. 우리 방법의 문제가 아니라")
print("이 세포의 성질이라는 뜻이죠.")

정답지와의 비율이 어떤 변수를 따라 움직이면, 그것이 내 방법의 편향입니다.

Python 브라우저에서 실행

import bu
import numpy as np
from scipy import stats

# 그러면 우리 방법에는 길이 편향이 따로 있나?
# 공개 표를 기준으로 삼고, '내 값 ÷ 공개 값' 이 길이에 따라 변하는지 본다.
T = bu.table(); panel = bu.panel()
cpm = bu.cpm(bu.published())
L = np.array([len(panel[g]) for g in T.index], dtype=float)
m = T["WT_1"].values.astype(float)
p = cpm["WT_DDX54_1"].reindex(T.index).fillna(0).values

sel = (m > 0) & (p > 0)
ratio = np.log2((m[sel] / m.sum() * 1e6 + 1) / (p[sel] + 1))
r = stats.pearsonr(np.log10(L[sel]), ratio)[0]

print(f"(내 값 / 공개 값) 과 길이의 상관 : {r:+.3f}   (유전자 {sel.sum()}개)")
print()
print("음수입니다 — 우리는 긴 유전자를 공개 파이프라인보다 '적게' 셉니다.")
print("유전자마다 전사체를 하나만(가장 긴 protein_coding) 골랐기 때문입니다.")
print("긴 유전자일수록 실제 발현되는 이형질체와 어긋날 여지가 크죠.")

'맞는 보정' 을 했는데 결과가 나빠지는 경우 — 이유를 찾아야 합니다.

Python 브라우저에서 실행

import bu
import numpy as np
from scipy import stats

# 그러면 TPM 으로 바꾸면 공개 표와 더 잘 맞을까?
T = bu.table(); panel = bu.panel()
cpm = bu.cpm(bu.published())
L = np.array([len(panel[g]) for g in T.index], dtype=float)
m = T["WT_1"].values.astype(float)
truth = cpm["WT_DDX54_1"].reindex(T.index).fillna(0).values

# 14강과 같은 자로 잽니다 — log2(값+1) 끼리의 pearson.
rpk = m / (L / 1000)
tpm = rpk / rpk.sum() * m.sum()          # 총합을 원래 카운트와 맞춘다

print(f"원래 카운트 그대로 : r = {stats.pearsonr(np.log2(m+1), np.log2(truth+1))[0]:.3f}")
print(f"길이로 나눈 뒤     : r = {stats.pearsonr(np.log2(tpm+1), np.log2(truth+1))[0]:.3f}")
print(f"(순위만 보는 spearman : {stats.spearmanr(m, truth)[0]:.3f} → "
      f"{stats.spearmanr(tpm, truth)[0]:.3f})")
print()
print("길이로 나눴더니 오히려 나빠졌습니다. 순위로 봐도 마찬가지입니다.")

고등학생그런데 나누니까 더 나빠졌다

이 강의에서 실제로 길이로 나눠 보면, 공개 표와의 상관이 떨어진다. 고치려던 것이 더 나빠졌다.

왜 이런 일이

가능한 이유	설명
3' 편향	7강 — 꼬리 쪽만 읽혔으면 길이 비례가 성립 안 함
패널의 길이 정보 부정확	어느 전사체 길이를 썼나
비교 대상도 보정됨	이미 보정된 값에 또 보정
read 가 적어 잡음이 큼	나누기가 잡음을 증폭

첫 줄이 가장 유력하다

3' 편향이 심하면 read 는 길이에 비례하지 않는다. 꼬리 근처에서만 나오니 짧든 길든 비슷한 수가 나온다. 이런 자료에 길이 보정을 하면 있지도 않은 편향을 고치려다 새 편향을 만든다.

보정은 그 편향이 실제로 있을 때만 도움이 된다

🎮 길이 실험실

길이와 카운트의 관계를 발현 구간별로 갈라 보세요.

고등학생언제 길이 보정을 하고 언제 하지 않나

비교	길이 보정	이유
같은 유전자, 시료 A vs B	필요 없음	길이가 같으니 상쇄됨
같은 시료 안, 유전자 A vs B	필요	길이가 다름
DESeq2 · edgeR 입력	절대 하지 말 것	원시 카운트를 요구함
유전자 세트 분석 · 히트맵	대개 함	유전자끼리 견주므로

셋째 줄은 흔한 실수다

차등발현 분석 도구들은 카운트가 정수여야 통계 모형이 성립한다. 음이항 분포는 세어진 개수에 대한 모형이기 때문이다. TPM 을 넣으면 모형의 전제가 깨져서 p 값이 엉망이 된다.

도구가 요구하는 입력 형식에는 대개 통계적 이유 가 있다. "왜 이 형식이어야 하는가" 를 물으면 그 도구를 이해하게 된다.

정규화는 목적별로 다르다 — 한 번 보정한 값을 모든 곳에 쓰지 않는다

🤔 확인 문제

길이와 카운트의 상관이 음수로 나온 이유는?

🤔 확인 문제

공개 카운트 표와 견줄 때 TPM 으로 바꾸면 안 되는 이유는?

해설 비교하는 두 값은 같은 처리를 거쳐야 합니다. 파일 이름이 cnt 인 것이 힌트였죠 — 무엇을 받았는지 먼저 확인해야 합니다.

🤔 확인 문제

DESeq2·edgeR 가 TPM 이 아니라 raw count 를 요구하는 이유는?