2. 세로로 쌓기 — bcftools mpileup · Dr. Son's Canine Mast Cell Tumour Therapeutics, First Steps II

BAM 안의 자료는 리드마다 한 줄입니다. 가로로 누워 있습니다.

그런데 변이를 찾으려면 물어야 할 것은 이겁니다 — "이 자리를 덮은 리드들이 무엇이라고 말하는가."

그러려면 자료를 세로로 세워야 합니다. 그 일을 하는 것이 mpileup 입니다.

참조    … A C G A T G C …

read1 … A C G A T G C

read2    C G T T G C A

read3 … A C G T T G C

read4 … A C G T T G

                  ↑ 이 자리: A 1개, T 3개

리드 다섯 개 가운데 세 개가 이 자리를 덮었습니다. DP = 3 입니다.

Python Runs in your browser

# pileup 을 손으로 만들어 본다
ref   = "ACGATGCTTAGGCTA"
reads = [
    ("read1",  0, "ACGATGCTTA"),
    ("read2",  2, "GTTGCATAGG"),
    ("read3",  1, "CGTTGCTTAG"),
    ("read4",  0, "ACGTTGCTTA"),
    ("read5",  3, "TTGCTTAGGC"),
]
POS = 3          # 0부터 세는 좌표. 참조는 여기서 'A'

print("참조   ", ref)
for name, start, seq in reads:
    print(f"{name:<7}{' ' * start}{seq}")
print(f"{'':<7}{' ' * POS}^  ← 우리가 묻는 자리 (참조 {ref[POS]})")
print()

col = []
for name, start, seq in reads:
    off = POS - start
    if 0 <= off < len(seq):
        col.append((name, seq[off]))

from collections import Counter
cnt = Counter(b for _, b in col)
print("이 자리를 덮은 리드:", len(col), "개")
for b, n in cnt.most_common():
    tag = "참조" if b == ref[POS] else "변이"
    print(f"   {b}  {n} 개   {tag}")
print()
print("이것이 pileup 한 칸입니다. mpileup 은 이 일을 유전체의 모든 자리에서 합니다.")

리드를 쌓고 문턱을 옮겨 보세요. 어느 리드가 세어지고 어느 리드가 빠지는지 나옵니다.

Loading…

노트북이 실제로 부르는 줄

bcftools mpileup -f kit_locus.fa -a FORMAT/AD,FORMAT/DP -Q15 -q10 tumor.sorted.bam

옵션	무엇	왜
`-f kit_locus.fa`	참조 FASTA	참조 염기를 알아야 '다르다'를 말할 수 있다. `.fai` 가 있어야 한다 (I 권 7강)
`-Q 15`	염기 품질 문턱	Q15 미만인 염기는 세지 않는다. 오류율 3% 이상인 글자
`-q 10`	MAPQ 문턱	자리를 못 믿는 리드를 뺀다 (I 권 14강). MAPQ 0 이 여기서 걸린다
`-a FORMAT/AD,FORMAT/DP`	주석 추가	대립유전자별 리드 수(AD)와 깊이(DP)를 출력에 남긴다. 없으면 VAF 를 계산할 수 없다

두 문턱이 서로 다른 것을 잰다는 것을 놓치기 쉽습니다. -Q 는 글자를 못 믿는 것, -q 는 자리를 못 믿는 것입니다.

문턱 두 개로 VAF 가 움직입니다. 어느 쪽이 옳은 값인지는 정해져 있지 않습니다 — 11강에서 다시 봅니다.

Python Runs in your browser

# 두 문턱이 무엇을 거르나 — 같은 자리, 다른 리드들
reads = [
    # (이름, 이 자리의 염기, 염기품질 Q, 리드의 MAPQ)
    ("r01", "A", 37, 60), ("r02", "A", 35, 60), ("r03", "A", 33, 60),
    ("r04", "T", 36, 60), ("r05", "T", 38, 60), ("r06", "T", 34, 60),
    ("r07", "T", 12, 60),   # 글자를 못 믿는다
    ("r08", "T",  8, 60),   # 글자를 못 믿는다
    ("r09", "T", 36,  0),   # 자리를 못 믿는다 — 반복 구간
    ("r10", "T", 35,  3),   # 자리를 못 믿는다
    ("r11", "G", 30, 60),   # 진짜 오류 하나
]
Q_MIN, MAPQ_MIN = 15, 10

kept, dropped = [], []
for name, b, q, mq in reads:
    if q < Q_MIN:      dropped.append((name, b, "염기품질 Q%d < %d" % (q, Q_MIN)))
    elif mq < MAPQ_MIN: dropped.append((name, b, "MAPQ %d < %d" % (mq, MAPQ_MIN)))
    else:               kept.append((name, b))

from collections import Counter
print("걸러진 리드")
for n, b, why in dropped:
    print(f"   {n}  {b}   {why}")
print()
c = Counter(b for _, b in kept)
dp = len(kept)
alt = c.get("T", 0)
print(f"남은 깊이 DP = {dp}")
print(f"   A {c.get('A',0)} 개 (참조)   T {alt} 개   G {c.get('G',0)} 개")
print(f"   VAF(T) = {alt}/{dp} = {alt/dp:.3f}")
print()
print("문턱을 안 걸었다면 VAF(T) = %d/%d = %.3f 였습니다." %
      (sum(1 for r in reads if r[1] == "T"), len(reads),
       sum(1 for r in reads if r[1] == "T") / len(reads)))

인델(삽입·결실)은 왜 어려운가

점돌연변이는 세로 한 칸만 보면 됩니다. 인델은 그렇지 않습니다.

점돌연변이

칸이 어긋나지 않는다. 그냥 세면 된다

인델

리드마다 정렬이 다르게 될 수 있다. 같은 인델이 서로 다른 자리에 적힌다 — 특히 반복 구간에서

AAAAA 안에서 A 하나가 빠지면, 어느 A 가 빠졌는지 알 방법이 없습니다. 정렬기는 아무 A 나 고릅니다. 그러면 같은 인델이 리드마다 다른 좌표로 적히고, 세로로 세었을 때 흩어져 문턱을 못 넘습니다.

해결책이 둘입니다 — 좌측 정렬(left-align) 로 표기를 통일하고, 국소 재정렬(indel realignment) 로 리드들을 다시 맞춥니다. bcftools 는 mpileup 안에서 어느 정도 자동으로 합니다(BAQ).

🤔 Check yourself

-Q 15 와 -q 10 이 각각 거르는 것은?

🤔 Check yourself

-a FORMAT/AD,FORMAT/DP 를 빼면 무슨 일이 생깁니까?

🤔 Check yourself

반복 구간(AAAAA) 안의 인델이 어려운 이유는?