14강 / 20강

카운트 표를 만들고 공개 표와 맞춰 본다

우리 손으로 만든 290 × 7 · 그리고 r = 0.76

12강에서 표를 만들었습니다. 이제 논문이 GEO 에 올린 표와 맞춰 봅니다. 같은 세포, 같은 read 에서 나온 표이니 맞아야 합니다.

단, 그대로 견주면 안 됩니다.

우리는 20,000 read, 저쪽은 2,700만 read — 규모가 다릅니다
우리는 290 유전자, 저쪽은 24,419 유전자 — 분모가 다릅니다

그래서 같은 290 유전자 안에서, 비율로 견줍니다.

고등학생드디어 표가 나온다

1강에서 "원자료는 표가 아니다" 로 시작했다. 열두 강을 지나 이제 표 를 만든다. 290 유전자 × 7 시료.

이 표 한 장에 들어 있는 결정들

결정	어느 강
트리밍 안 함	8강
k=31 색인	10강
가닥 처리 방식	11강
다중 매핑 처리	12강
이름 못 찾은 16개 제외	13강

논문에 실리는 표는 이 결정들을 보여주지 않는다. 표는 사실이 아니라 여러 판단의 결과물 이라는 것 — 이 강의의 첫 문장이 여기서 완결된다.

남의 표를 볼 때도 이 목록을 떠올릴 수 있어야 한다

규모를 맞추려면 내 카운트도 비율(ppm)로 바꿔야 합니다.

Python 브라우저에서 실행

import bu
import numpy as np, pandas as pd

T = bu.table()
pub = bu.published()
cpm = bu.cpm(pub)

g = T.index
mine = T["WT_1"].astype(float)
theirs = cpm["WT_DDX54_1"].reindex(g).fillna(0)

d = pd.DataFrame({"내 카운트": mine.astype(int),
                  "공개 CPM": theirs.round(0).astype(int)})
d["내 비율(ppm)"] = (mine / mine.sum() * 1e6).round(0).astype(int)
print(d.sort_values("공개 CPM", ascending=False).head(12).to_string())
print()
print(f"내 총 카운트 {int(mine.sum()):,} · 공개 표의 이 290 유전자 합 {theirs.sum():,.0f} CPM")

산점도가 대각선을 따라가면 맞는 것입니다. 흩어진 정도가 곧 오차예요.

Python 브라우저에서 실행

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

T = bu.table(); cpm = bu.cpm(bu.published())
g = T.index
mine = T["WT_1"].astype(float)
theirs = cpm["WT_DDX54_1"].reindex(g).fillna(0)

r = stats.pearsonr(np.log2(mine + 1), np.log2(theirs + 1))[0]
rho = stats.spearmanr(mine, theirs)[0]

fig, ax = plt.subplots(figsize=(5.8, 5.4))
ax.scatter(np.log2(theirs + 1), np.log2(mine + 1), s=14, alpha=.55, c="#4c6ef5")
ax.set_xlabel("공개 표  log2(CPM+1)"); ax.set_ylabel("내 표  log2(count+1)")
ax.set_title(f"290 유전자 · pearson r = {r:.3f}")
plt.tight_layout()

print(f"pearson  (log2) : {r:.3f}")
print(f"spearman (순위) : {rho:.3f}")
print(f"내 카운트가 0 인 유전자 : {int((mine==0).sum())} / {len(g)}")

고등학생공개 표와 맞춰 본다

논문 저자들이 GEO 에 올린 카운트 표가 있다. 우리가 만든 표와 견주면 우리가 얼마나 잘했는지 를 잴 수 있다.

무엇으로 재나 — 상관계수

같은 유전자에 대해 우리 값과 저자 값을 짝지어 놓고, 그 짝들이 얼마나 함께 움직이는지를 하나의 수로 요약한 것이 상관계수 r 이다.

r	뜻
1.0	완전히 같은 방향으로 움직임
0.7~0.9	강한 일치
0	아무 관계 없음

결과는 r = 0.76 이다. 좋은가? 나쁜가? 비교 대상이 없으면 답할 수 없다. 17강이 그 비교 대상을 만든다.

숫자 하나만으로는 아무 판단도 못 한다 — 기준이 있어야 한다

r = 0.76 — 좋은 건가 나쁜 건가

혼자서는 판단할 수 없는 숫자입니다. 견줄 기준이 필요합니다.

무엇과 무엇	r
공개 표의 WT_1 대 WT_2 (같은 290 유전자)	0.987
우리 표 대 공개 표 (WT_1)	0.762
우리 표 대 공개 표 (평균)	0.638

공개 표끼리는 0.987 로 거의 완벽합니다. 우리는 0.762 죠. 차이의 정체를 15강(깊이)과 17강(천장)에서 나눠 잽니다.

자를 하나로 못 박아 둡니다. 이 강의에서 "r" 은 언제나 pearson(log2(내 카운트+1), log2(공개 CPM+1)) 입니다. +1 을 어디에 어떻게 더하느냐로 값이 0.05 쯤 움직이므로, 같은 자로 재지 않은 숫자끼리 견주면 안 됩니다. 순위만 보는 spearman(0.732)은 이 선택에서 자유롭고, 늘 함께 적어 둡니다.

일곱 줄을 나란히 놓으면 깊이가 성적을 정한다는 것이 바로 보입니다.

Python 브라우저에서 실행

import bu
import numpy as np
from scipy import stats

# 일곱 표본 전부
T = bu.table(); cpm = bu.cpm(bu.published())
print(f"{'표본':<7}{'실은 read':>10}{'붙은 것':>9}{'pearson':>10}{'spearman':>10}")
print("-" * 48)
for s in bu.SAMPLES:
    m = T[s].astype(float)
    p = cpm[bu.PUB[s]].reindex(T.index).fillna(0)
    r = stats.pearsonr(np.log2(m + 1), np.log2(p + 1))[0]
    rho = stats.spearmanr(m, p)[0]
    print(f"{s:<7}{len(bu.reads(s)):>10,}{int(m.sum()):>9,}{r:>10.3f}{rho:>10.3f}")
print()
print("WT_1 만 유난히 높습니다 — read 를 5배 실었기 때문입니다.")
print("표본의 성질이 아니라 깊이의 차이예요. 다음 강에서 확인합니다.")

고등학생왜 완벽하게 안 맞나

같은 원자료인데 왜 r 이 1.0 이 아닌가? 이유가 여럿이다.

차이	우리	저자
참조	290 유전자 패널	전체 전사체
read 수	44,000	2,500만
다중 매핑	단순 처리	EM 배분
이어붙이기	처리 안 함	처리함
불일치 허용	없음	있음

이 중 가장 큰 것은

read 수 다. 2,500만 개에서 44,000 개를 뽑았으니 570분의 1이다. 발현량이 낮은 유전자는 우리 표본에 아예 안 들어왔을 수 있다. 이것은 알고리즘의 문제가 아니라 뽑기 잡음 이고, 아무리 잘 짜도 피할 수 없다.

17강에서 계산할 것이 정확히 이 부분이다 — 뽑기 잡음만으로 정해지는 상한.

차이의 목록을 만들고, 그중 피할 수 없는 것을 골라내는 것이 다음 단계다

🎮 표 대조기

유전자를 검색해 내 카운트와 공개 값을 나란히 보세요.

고등학생상관계수를 볼 때 조심할 것

r 하나로 판단하면 안 되는 이유가 몇 가지 있다.

함정	내용
큰 값이 지배	발현량이 큰 유전자 몇 개가 r 을 좌우
로그를 안 취함	발현량은 몇 자릿수 차이 — 로그가 보통 옳음
0 이 많음	안 잡힌 유전자가 많으면 r 이 부풀거나 꺼짐
그림을 안 봄	같은 r 이 전혀 다른 산점도에서 나옴

마지막 줄이 가장 중요하다

앤스콤의 4중주가 유명한 예다. 평균·분산·상관계수가 모두 같은 네 개의 자료가 그림으로는 완전히 다르게 생겼다. 요약 통계만 보고 그림을 안 보면 속는다.

이 강의에서는 반드시 산점도를 함께 그린다. /viz 와 /stat 이 같은 말을 한다.

r 을 보고했으면 산점도도 보여야 한다 — 둘은 한 벌이다

🤔 확인 문제

우리 카운트와 공개 CPM 을 견줄 때 같은 290 유전자 안에서 비교해야 하는 이유는?

🤔 확인 문제

WT_1 의 상관(0.762)이 다른 표본(0.57~0.66)보다 높은 이유는?

🤔 확인 문제

공개 표의 두 표본끼리 r = 0.987 이라는 사실이 알려 주는 것은?