12강 / 20강

44,000 read 를 290 유전자에 붙인다

붙이는 함수를 짜고 · 어느 유전자인지 정할 수 없는 read 를 만난다

재료가 다 모였습니다. 색인(10강)과 역상보(11강)를 합치면 붙이는 함수가 됩니다.

고등학생붙이는 함수를 직접 짠다

이제 조각들을 모은다. 44,000개의 read 를 290개 유전자에 배정하는 함수를 여덟 줄쯤 으로 짠다.

그 여덟 줄이 하는 일

줄	하는 일
1~2	read 를 k-mer 로 자름
3~4	각 k-mer 를 색인에서 조회
5~6	유전자별로 맞은 횟수를 셈
7~8	가장 많이 맞은 유전자를 고름

salmon 이나 featureCounts 가 하는 일의 뼈대 가 이것이다. 그 도구들은 여기에 정교함을 얹었을 뿐, 골격은 이 여덟 줄이다.

복잡한 도구도 뼈대는 단순하다 — 그 뼈대를 아는 것이 도구를 아는 것이다

여덟 줄입니다. 이것이 우리의 정렬기 전부예요.

Python 브라우저에서 실행

import bu

def assign(seq, idx, k=31):
    "read 하나를 유전자 하나에 붙인다. 못 붙이면 None."
    for p in range(len(seq) - k + 1):          # ① 정방향으로 훑는다
        g = idx.get(seq[p:p+k])
        if g:
            return g
    r = bu.rc(seq)
    for p in range(len(r) - k + 1):            # ② 없으면 역상보로 훑는다
        g = idx.get(r[p:p+k])
        if g:
            return g
    return None                                 # ③ 그래도 없으면 포기

idx = bu.index()
for h, s, q in bu.reads("WT_1", 6):
    print(f"{assign(s, idx) or '—':<12}{s[:56]}")

카운트 표의 첫 열이 만들어졌습니다.

Python 브라우저에서 실행

import bu, time
from collections import Counter

idx = bu.index()
t = time.time()
c = Counter()
n = 0
for _, s, _q in bu.reads("WT_1"):
    n += 1
    g = bu.assign(s, idx)
    if g:
        c[g] += 1

print(f"{n:,} read 를 붙이는 데 {time.time()-t:.1f} 초")
print(f"붙은 read : {sum(c.values()):,}  ({sum(c.values())/n*100:.2f} %)")
print(f"나온 유전자 : {len(c)} / {len(bu.panel())}")
print()
print("── 가장 많이 붙은 유전자 12개 ──")
for g, k in c.most_common(12):
    print(f"  {g:<12}{k:>5}  {'▇' * (k // 4)}")

고등학생정할 수 없는 read 를 만난다

돌려 보면 반드시 나오는 문제가 있다. 두 유전자에 똑같이 잘 붙는 read.

왜 생기나

원인	예
유전자족(paralog)	진화적으로 복제된 형제 유전자
공통 도메인	여러 단백질이 공유하는 기능 부위
반복 서열	Alu, LINE 같은 것
같은 유전자의 다른 전사체	공통 엑손

이건 우리 코드의 결함이 아니다. 유전체가 원래 그렇게 생겼다. 76글자만 봐서는 구별이 안 되는 자리가 실제로 존재한다.

모호함은 알고리즘의 실패가 아니라 자료의 성질이다

이름을 보면 알 수 있는 것

상위 목록이 Rpl·Rps(리보솜 단백질), Eef1a1(번역 인자), Actb(액틴), Hsp90ab1(샤페론)로 채워집니다.

이것은 세포가 가장 부지런히 만드는 것들입니다. 어떤 조직, 어떤 실험이든 상위 목록은 대개 이렇습니다. 목록이 이렇지 않다면 그것이 오히려 신호예요 — 뭔가 잘못됐거나, 아주 특별한 세포이거나.

이것을 multi-mapping 이라고 합니다. 진짜 도구들은 여기서 훨씬 정교합니다.

Python 브라우저에서 실행

import bu

# 어느 유전자인지 정할 수 없는 read
idx = bu.index(); K = bu.K
amb = 0
examples = []
for _, s, _q in bu.reads("WT_1", 5000):
    gs = set()
    for p in range(len(s) - K + 1):
        for probe in (s[p:p+K], bu.rc(s)[p:p+K]):
            g = idx.get(probe)
            if g:
                gs.add(g)
    if len(gs) > 1:
        amb += 1
        if len(examples) < 6:
            examples.append((sorted(gs), s[:40]))

print(f"두 개 이상의 유전자에 걸린 read : {amb} / 5,000  ({amb/50:.2f} %)")
print()
for gs, s in examples:
    print(f"  {', '.join(gs):<34}{s}")
print()
print("우리 함수는 '먼저 걸린 것' 을 그냥 씁니다 — 정직하지 않은 처리입니다.")

고등학생다중 매핑을 처리하는 네 가지 방법

방법	내용	문제
버린다	모호한 read 는 세지 않음	유전자족 발현을 과소평가
모두에 센다	후보 전부에 1 씩	총합이 read 수보다 커짐
나눠 센다	후보가 n 개면 1/n 씩	실제 비율을 무시
EM 으로 추정	고유 read 로 비율을 추정한 뒤 그 비율대로 배분	계산이 무거움

마지막 방법이 왜 나은가

유전자 A 와 B 가 헷갈린다고 하자. 그런데 A 에만 고유하게 붙는 read 가 900개, B 에만 900개 대 100개라면? 헷갈리는 read 도 9:1 로 나누는 것이 합리적 이다.

이것이 salmon·kallisto·RSEM 이 쓰는 EM 알고리즘의 발상이다. 확실한 것으로 불확실한 것을 배분한다. /ml 에서 배우는 EM 과 같은 논리다.

모르는 것을 아는 것의 비율로 나눈다 — 통계적 배분의 기본 발상

이제 표가 생겼습니다. 논문이 GEO 에 올린 것과 같은 모양이죠.

Python 브라우저에서 실행

import bu, time

# 일곱 표본을 모두 붙여 카운트 표를 만든다
t = time.time()
T = bu.table()
print(f"{T.shape[0]} 유전자 × {T.shape[1]} 표본   ({time.time()-t:.1f} 초)")
print()
print(T.loc[["Actb", "Ddx54", "Rpl13", "Hsp90ab1", "Vim"]])
print()
print("표본별 총 카운트")
print(T.sum().to_string())

🎮 붙이기 실험실

read 를 하나씩 넘기며 어디에 어떻게 붙는지 따라가 보세요.

고등학생우리가 대충 넘어간 곳

이 강의의 여덟 줄이 진짜 도구와 다른 점을 정직하게 적는다.

우리가 안 한 것	영향
다중 매핑을 EM 으로 배분	유전자족 발현량 부정확
이어붙이기(엑손 경계) 처리	엑손을 넘는 read 를 놓침 (18강)
불일치 허용	변이·오류가 있는 read 를 놓침
가닥 방향 엄격 적용	일부 read 를 잘못 세거나 놓침
유전자 길이 보정	16강에서 다룸

그래도 이 표를 만드는 이유

완벽하지 않은 표를 만들고 어디가 왜 부족한지 아는 것 이, 완벽한 도구를 블랙박스로 쓰는 것보다 낫다. 14강에서 공개 표와 비교해 우리 표가 얼마나 맞는지 재고, 17강에서 애초에 얼마나 맞을 수 있었는지 를 계산한다.

부족함의 목록을 쓸 수 있으면 그 결과는 쓸 만하다

🤔 확인 문제

붙는 비율이 26% 밖에 안 되는 가장 큰 이유는?

🤔 확인 문제

multi-mapping read 를 '먼저 걸린 유전자에 준다' 는 처리의 문제는?

🤔 확인 문제

카운트 상위가 Rpl·Rps·Eef1a1·Actb 로 채워진 것은?