Lesson 17 of 20

커버리지 — 평균이 감추는 것

<code>samtools depth</code> · GC 편향 · 코드랩 3단계 끝

10강에서 깊이가 무엇인지 봤습니다. 이 강은 그것을 실제로 재는 이야기입니다.

samtools depth -a tumor.sorted.bam     # 자리마다 한 줄: 염색체 · 위치 · 깊이

samtools coverage tumor.sorted.bam      # 염색체마다 요약 한 줄

-a 가 중요합니다. 없으면 리드가 하나라도 있는 자리만 나옵니다. 깊이 0 인 자리 — 즉 가장 알고 싶은 자리 — 가 표에서 사라집니다.

그래서 보고서에는 평균 대신 "X× 이상으로 덮인 표적의 비율" 을 적습니다.

Python Runs in your browser

# 평균만 보면 놓치는 것 — 같은 평균, 다른 자료
import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(11)
n = 20000

flat  = rng.poisson(60, n)                                   # 고르다
gcbad = np.concatenate([rng.poisson(85, n//2), rng.poisson(35, n//2)])  # GC 편향
hole  = np.concatenate([rng.poisson(63, int(n*.9)), np.zeros(int(n*.1), int)])  # 구멍

fig, ax = plt.subplots(1, 3, figsize=(9.6, 2.8), sharey=True)
for a, d, t in zip(ax, (flat, gcbad, hole), ("even", "GC bias", "10% gaps")):
    a.hist(d, bins=50); a.set_title(f"{t}  mean={d.mean():.0f}x")
    a.set_xlabel("depth")
plt.tight_layout(); plt.show()

print(f"{'자료':<12}{'평균':>7}{'중앙값':>8}{'10x 미만':>10}{'0x':>8}")
for d, t in zip((flat, gcbad, hole), ("고르다", "GC 편향", "구멍")):
    print(f"{t:<12}{d.mean():>7.0f}{np.median(d):>8.0f}"
          f"{(d<10).mean():>10.1%}{(d==0).mean():>8.1%}")
print()
print("평균은 셋 다 60x 근처입니다. 쓸 수 있는 자료는 전혀 다릅니다.")

실무에서 실제로 보고하는 값

지표	무엇	왜
평균 깊이	전체 평균	비용을 정한다. 품질은 못 말한다
중앙값 깊이	가운데 값	극단값에 안 흔들린다. 평균보다 낫다
C10 · C20 · C30	10× / 20× / 30× 이상으로 덮인 표적의 비율	이것이 진짜 지표. "C20 = 95%" 처럼 쓴다
균일성	평균의 0.2배 이상인 자리의 비율	편향을 잡는다
0× 구간	아예 안 덮인 자리	여기는 결과를 말할 수 없다

임상 보고서에는 0× 구간 목록이 붙습니다. "이 exon 은 못 봤습니다" 를 적어 두지 않으면, 읽는 사람이 "변이 없음" 으로 오해합니다.

왜 고르지 않은가 — GC 편향

PCR 은 GC 가 아주 높거나 아주 낮은 조각을 잘 못 늘립니다. GC 가 높으면 가닥이 잘 안 풀리고, 낮으면 프라이머가 잘 안 붙습니다.

결과는 종 모양입니다 — GC 40~50% 근처는 잘 덮이고, 양 끝은 얕아집니다.

여기서 6강의 이야기가 돌아옵니다. Dog10K 조립이 CanFam3.1 대비 특히 많이 채운 것이 GC 가 높은 구간이었습니다. 그 구간에 프로모터와 첫 exon 이 몰려 있습니다.

읽기도 어렵고, 참조에도 없었다

두 어려움이 같은 자리에 겹칩니다.

가장 중요한 자리가 가장 얕습니다. 그래서 GC 가 높은 표적을 노리는 패널은 별도의 프로토콜을 씁니다.

Python Runs in your browser

# GC 편향의 모양 — 그리고 그것이 어느 유전체 요소를 때리나
import numpy as np
import matplotlib.pyplot as plt

gc = np.linspace(0.15, 0.85, 200)
rel = np.exp(-((gc - 0.44) ** 2) / (2 * 0.11 ** 2))    # 종 모양

plt.figure(figsize=(7.2, 3.2))
plt.plot(gc * 100, rel * 60, lw=2.2)
plt.axhline(20, ls="--", lw=1.2)
plt.xlabel("GC content (%)"); plt.ylabel("expected depth (x)")
plt.tight_layout(); plt.show()

print(f"{'구간':<22}{'대표 GC':>9}{'기대 깊이':>11}")
for name, g in [("유전자 사이 (AT rich)", 0.36), ("전체 평균", 0.41),
                ("코딩 exon", 0.52), ("CpG 섬 · 프로모터", 0.66),
                ("첫 exon (GC rich)", 0.72)]:
    d = 60 * float(np.exp(-((g - .44) ** 2) / (2 * .11 ** 2)))
    flag = "  ← 위험" if d < 20 else ""
    print(f"{name:<22}{g:>8.0%}{d:>11.0f}x{flag}")

노트북에서 커버리지를 보는 자리

코드랩 3단계 끝에서 samtools depth 로 위치별 깊이를 그립니다. 종양이 정상보다 두 배쯤 두껍게 나오는 것이 보입니다.

여기서 눈으로 확인할 것 두 가지입니다.

📏

양 끝이 얕다

20 kb 창의 가장자리는 리드가 절반만 걸친다. 가장자리 효과

📊

비율이 맞다

종양 : 정상 ≈ 2 : 1. 심으라고 준 커버리지가 실제로 들어왔는가

가장자리 효과 때문에 변이가 창의 끝에 있으면 놓칩니다. 노트북이 드라이버를 창 가운데에 심는 이유입니다.

그래서 표적 구간을 잡을 때 양쪽에 여유(padding) 를 100~200 bp 붙입니다.

Python Runs in your browser

# 가장자리 효과 — 창 끝은 왜 얕은가
import numpy as np
import matplotlib.pyplot as plt

W, read_len = 2000, 150
starts = np.arange(0, W - read_len)        # 리드가 시작할 수 있는 자리
depth = np.zeros(W)
for s in starts:
    depth[s:s + read_len] += 1
depth = depth / depth.max() * 60

plt.figure(figsize=(7.2, 2.8))
plt.plot(depth, lw=1.8)
plt.xlabel("position in window (bp)"); plt.ylabel("relative depth (x)")
plt.tight_layout(); plt.show()

for p in (0, 50, 149, 150, 1000, W - 150, W - 1):
    print(f"pos {p:>5}   depth {depth[p]:>5.1f}x")
print()
print("리드 길이만큼(150 bp)의 양 끝이 삼각형으로 깎입니다.")
print("실제 유전체에서는 염색체 끝과 조립 구멍 가장자리가 전부 이렇게 생깁니다.")

🤔 Check yourself

samtools depth 를 -a 없이 돌리면 무엇이 빠집니까?

Answer -a 없이는 리드가 하나라도 있는 자리만 출력됩니다. 안 덮인 자리가 표에서 사라지므로, 평균만 계산하면 실제보다 좋아 보입니다. 임상 보고서에 0× 구간 목록을 따로 붙이는 이유도 같습니다.

🤔 Check yourself

평균 깊이 60× 인 두 시료 가운데 어느 쪽이 더 쓸 만한지 가리려면 무엇을 봐야 합니까?

🤔 Check yourself

GC 편향이 특히 아픈 이유를 개 유전체 맥락에서 고르면?