13강 / 20강

SAM · BAM — 플래그와 CIGAR

정렬 결과를 실제로 열어 보는 일 · <code>samtools view</code>

정렬이 끝나면 SAM 이 나옵니다. 탭으로 나뉜 표이고, 리드 하나가 한 줄입니다. 그것을 압축하고 색인한 것이 BAM 입니다. 내용은 같습니다.

#	이름	무엇
1	QNAME	리드 이름. 짝인 두 리드가 같은 이름을 갖는다
2	FLAG	비트 열두 개를 더한 정수. 이 강의의 절반
3·4	RNAME · POS	어느 염색체, 몇 번 자리 (1부터)
5	MAPQ	이 자리를 얼마나 확신하는가 — 14강
6	CIGAR	어떻게 붙었는가. 이 강의의 나머지 절반
7·8·9	RNEXT · PNEXT · TLEN	짝은 어디 있나, 조각 전체 길이 — 15강
10·11	SEQ · QUAL	FASTQ 의 둘째·넷째 줄이 그대로 온 것

samtools view -f 4 는 '4 비트가 켜진 것만', -F 1024 는 '1024 비트가 꺼진 것만' 입니다.

Python 브라우저에서 실행

# FLAG 는 비트를 더한 것 — 되돌려 읽어 본다
BITS = [
    (1,    "짝이 있는 리드 (paired)"),
    (2,    "짝이 제대로 붙음 (proper pair)"),
    (4,    "이 리드가 안 붙음 (unmapped)"),
    (8,    "짝이 안 붙음 (mate unmapped)"),
    (16,   "이 리드가 역방향 (reverse)"),
    (32,   "짝이 역방향 (mate reverse)"),
    (64,   "첫 번째 리드 = R1"),
    (128,  "두 번째 리드 = R2"),
    (256,  "이차 정렬 (secondary)"),
    (512,  "품질 검사 실패"),
    (1024, "PCR/광학 중복  ← 16강"),
    (2048, "보조 정렬 (supplementary)"),
]

def explain(flag):
    return [name for bit, name in BITS if flag & bit]

for f in (99, 147, 1024 + 99, 4, 2048 + 16):
    print(f"FLAG {f}")
    for line in explain(f):
        print("   ·", line)
    print()

print("자주 보는 값:")
print("  99  = 1+2+32+64    정방향 R1, 짝이 제대로 붙음")
print(" 147  = 1+2+16+128   역방향 R2, 짝이 제대로 붙음")
print("  4                  안 붙은 리드")

CIGAR — 리드가 참조에 어떻게 붙었나

150M 이면 150글자가 그대로 붙었다는 뜻입니다. 글자와 숫자를 번갈아 적습니다.

글자	뜻	참조 좌표가	리드 좌표가
M	정렬됨 (일치 또는 불일치)	나아간다	나아간다
I	삽입 — 리드에만 있다	멈춘다	나아간다
D	결실 — 참조에만 있다	나아간다	멈춘다
S	soft clip — 끝을 안 쓰고 남겨 둠	멈춘다	나아간다
H	hard clip — 끝을 아예 잘라 버림	멈춘다	멈춘다
N	건너뜀 — RNA 에서 인트론	나아간다	멈춘다

M 이 일치를 뜻하지 않는다는 것이 첫 함정입니다. 150M 인 리드에 변이가 다섯 개 있을 수도 있습니다. 진짜 일치·불일치를 구별하는 = 와 X 가 있지만 거의 안 쓰입니다.

75M75S 처럼 절반만 붙은 리드가 한자리에 몰리면 구조변이의 신호입니다. ITD 도 이렇게 보일 수 있습니다.

Python 브라우저에서 실행

# CIGAR 를 풀어 참조 위 길이와 리드 위 길이를 따로 센다
import re

def parse(cig):
    return [(int(n), op) for n, op in re.findall(r"(\d+)([MIDNSHP=X])", cig)]

REF = set("MDN=X")       # 참조 좌표를 소비하는 연산
QRY = set("MIS=X")       # 리드 좌표를 소비하는 연산

def spans(cig):
    ops = parse(cig)
    return (sum(n for n, o in ops if o in REF),
            sum(n for n, o in ops if o in QRY))

cases = [
    ("150M",          "그대로 붙음"),
    ("100M30I20M",    "30 bp 삽입 — ITD 가 이렇게 보인다"),
    ("100M30D20M",    "30 bp 결실"),
    ("20S130M",       "앞 20 bp 를 안 씀 — 어댑터나 poly-G 일 수 있다"),
    ("75M75S",        "절반만 붙음 — 구조변이 의심"),
]
print(f"{'CIGAR':<14}{'참조 위':>8}{'리드 위':>8}   메모")
for cig, memo in cases:
    r, q = spans(cig)
    print(f"{cig:<14}{r:>8}{q:>8}   {memo}")
print()
print("리드 위 길이는 언제나 150 이어야 합니다 (H 가 없다면).")
print("참조 위 길이는 삽입·결실에 따라 달라집니다 — 그래서 END 를 따로 계산해야 합니다.")

samtools 로 실제로 하는 일 넷

             samtools view -b  aln.sam > aln.bam    # SAM → BAM (압축)

             samtools sort -o aln.sorted.bam aln.bam  # 좌표 순서로 줄 세우기

             samtools index aln.sorted.bam          # .bai 색인 — 임의 접근

             samtools depth -a aln.sorted.bam        # 자리별 깊이 — 17강

정렬(sort)이 필수인 이유는 7강의 .fai 와 같습니다 — 줄이 좌표 순서로 서 있어야 색인을 만들 수 있고, 색인이 있어야 "13번 염색체 47 Mb 근처만" 을 바로 꺼낼 수 있습니다.

정렬하지 않은 BAM 에 index 를 걸면 오류가 납니다. 다행히 이건 시끄러운 사고입니다.

실제 .bai 는 구간을 계층적으로 묶은 R-tree 비슷한 구조지만, 원리는 이 이진 탐색과 같습니다.

Python 브라우저에서 실행

# 왜 정렬(sort)이 필요한가 — 색인이 가능한 형태
reads = [("13", 47_258_400), ("13", 47_258_010), ("13", 47_259_900),
         ("13", 47_258_120), ("13", 47_259_100)]

print("정렬 전 (정렬기가 뱉은 순서 = 읽어 온 순서):")
for r in reads:
    print("  ", r[0], f"{r[1]:,}")

print()
print("47,258,000 ~ 47,258,200 구간을 꺼내려면 → 파일을 처음부터 끝까지 훑어야 함")
print()

reads_sorted = sorted(reads, key=lambda r: (r[0], r[1]))
print("정렬 후:")
for r in reads_sorted:
    print("  ", r[0], f"{r[1]:,}")
print()
import bisect
pos = [r[1] for r in reads_sorted]
lo = bisect.bisect_left(pos, 47_258_000)
hi = bisect.bisect_right(pos, 47_258_200)
print(f"이제는 이진 탐색 두 번: 인덱스 {lo} ~ {hi} 만 읽으면 됩니다")
print("→ 이것이 .bai 색인이 하는 일입니다")

노트북 3단계를 다시 읽으면

amap = f"{TOOLS['bwa']} mem -t 4 -R '{rg}' kit_locus.fa {r1} {r2}"

sh(f"{amap} 2>/dev/null | {samtools} sort -o {sample}.sorted.bam -")

sh(f"{samtools} index {sample}.sorted.bam")

파이프(|)로 이어져 있는 것이 요점입니다. 정렬기의 출력을 디스크에 쓰지 않고 바로 sort 에 넘깁니다.

실제 WGS 에서 중간 SAM 파일은 수백 GB 입니다. 그것을 쓰고 다시 읽는 것만으로 몇 시간이 갑니다. 파이프가 그 일을 없앱니다.

마지막의 - 가 "표준입력에서 읽어라" 입니다.

🤔 확인 문제

SAM 의 CIGAR 에서 M 의 뜻으로 정확한 것은?

해설 M 은 alignment match 이지 sequence match 가 아닙니다. 150M 인 리드에 SNV 가 다섯 개 있을 수 있습니다. 진짜 일치·불일치를 구별하는 =·X 가 규격에 있지만 실무 도구는 거의 쓰지 않습니다.

🤔 확인 문제

samtools index 를 걸기 전에 반드시 해야 하는 것은?

🤔 확인 문제

정렬기 출력을 파일로 쓰지 않고 파이프로 바로 sort 에 넘기는 이유는?

해설 실제 WGS 에서 비압축 SAM 은 원본 FASTQ 보다도 큽니다. 파이프로 이으면 디스크 왕복이 사라집니다. 노트북 3단계 끝의 - 가 '표준입력에서 읽어라' 입니다.