Lesson 12 of 20

BWA-MEM 과 minimap2 — 언제 무엇을

코드랩 3단계 · <code>ALIGNER</code> 분기가 하는 일

워크숍 노트북 0단계에 이런 줄이 있습니다.

for n in ("bwa", "minimap2", "samtools", "bcftools"):

    p = shutil.which(n)

그리고 3단계에서 ALIGNER 가 "bwa" 인지 아닌지로 갈라집니다. 왜 둘을 다 챙겨 두었을까요.

답은 단순합니다 — 설치가 되는 쪽을 씁니다. Colab 에서는 apt 로 둘 다 되지만, 관리자 권한이 없는 서버에서는 하나만 있을 수 있습니다. 그리고 둘의 결과가 거의 같기 때문에 이런 분기가 가능합니다.

둘의 차이

	BWA-MEM	minimap2
나온 해	2013	2018
만든 사람	둘 다 Heng Li — 같은 사람입니다
노리는 리드	짧은 리드 (70 bp ~ 1 kb)	긴 리드 (ONT · PacBio) 와 짧은 리드 둘 다
씨앗 방식	FM-index 로 정확한 최대 일치	minimizer — 해시로 성긴 씨앗
속도 (짧은 리드)	기준	비슷하거나 조금 빠르다
스플라이스 인식	없다	`-x splice` 로 RNA 도 붙인다
DNA 짧은 리드 표준	사실상 표준. 대부분의 논문이 이것	쓸 수 있지만 관행이 아니다

같은 사람이 만들었고, 나중 것이 앞의 것을 대체하려고 만든 것이 아닙니다. 다른 문제를 풀려고 만든 것입니다.

minimizer 가 무엇인가 — 성기게 골라 잡는 씨앗

BWA 는 리드의 모든 자리에서 씨앗을 찾습니다. 정확하지만 무겁습니다.

minimap2 는 다르게 합니다. 창을 하나 잡고, 그 안의 k-mer 가운데 해시값이 가장 작은 것 하나만 씨앗으로 씁니다. 이것이 minimizer 입니다.

요점은 이겁니다 — 같은 서열이면 같은 minimizer 가 뽑힙니다. 그래서 씨앗을 10분의 1 로 줄여도 두 서열이 겹치면 대개 같은 씨앗을 공유합니다.

긴 리드는 오류가 많아 정확한 일치를 찾기 어렵습니다. 대신 길어서 성긴 씨앗으로도 충분합니다. minimap2 가 긴 리드에 강한 이유입니다.

이 성김이 속도의 전부입니다. 그리고 같은 것이 뽑힌다는 성질이 정확도를 지킵니다.

Python Runs in your browser

# minimizer 를 손으로 뽑아 본다
def h(s):
    "간단한 해시 — 실제 구현은 다르지만 원리는 같다"
    v = 0
    for c in s:
        v = (v * 131 + ord(c)) % 100003
    return v

seq = "ACGTTTCAGGATCCAAGGTTACGGATCC"
k, w = 5, 6          # k-mer 길이 5, 창 크기 6

mins = []
for i in range(len(seq) - k - w + 2):
    win = [(h(seq[j:j+k]), j) for j in range(i, i + w)]
    mins.append(min(win))

uniq = sorted(set(mins), key=lambda x: x[1])
print("서열       :", seq)
print("k-mer 총수 :", len(seq) - k + 1)
print("minimizer  :", len(uniq), "개")
print()
for hv, pos in uniq:
    print(f"  pos {pos:>3}  {seq[pos:pos+k]}  hash {hv:>6}")
print()
print(f"씨앗이 {len(seq)-k+1} 개에서 {len(uniq)} 개로 줄었습니다.")
print("같은 서열이면 같은 것이 뽑히므로, 겹치는 두 서열은 여전히 씨앗을 공유합니다.")

실제로 부르는 명령 — 노트북 3단계

BWA 쪽

bwa index kit_locus.fa

bwa mem -t 4 -R '@RG\tID:tumor\tSM:tumor\tPL:ILLUMINA' kit_locus.fa tumor_R1.fastq.gz tumor_R2.fastq.gz

minimap2 쪽

minimap2 -ax sr -t 4 -R '@RG\tID:tumor\tSM:tumor\tPL:ILLUMINA' kit_locus.fa tumor_R1.fastq.gz tumor_R2.fastq.gz

토막	뜻
`-ax sr`	minimap2 에게 짧은 리드(short read) 라고 알려 주는 프리셋. `-a` 는 SAM 으로 출력
`-t 4`	실 4개. 정렬은 리드마다 독립이라 코어 수만큼 거의 그대로 빨라진다
`-R '@RG…'`	읽기 그룹. 이 리드들이 어느 시료에서 왔는지 BAM 에 새긴다

출력의 칸 사이는 진짜 탭입니다. 셸에서는 -R '@RG\tID:…' 처럼 역슬래시-t 두 글자로 적고, bwa 가 그것을 탭으로 바꿔 넣습니다.

Python Runs in your browser

# 읽기 그룹 문자열을 만들어 본다 — 어느 칸이 무엇인가
def rg(sample, lane=1, lib="lib1", platform="ILLUMINA"):
    return "\t".join([
        "@RG",
        f"ID:{sample}.{lane}",     # 이 정렬 실행을 가리키는 고유 이름
        f"SM:{sample}",            # 시료 — 변이 부르기가 보는 것
        f"LB:{lib}",               # 라이브러리 — 중복 표시가 보는 것
        f"PU:{sample}.{lane}.unit",# 플랫폼 유닛 — 플로우셀+레인
        f"PL:{platform}",          # 플랫폼
    ])

for s in ("tumor", "normal"):
    print(rg(s))
print()
print("SM 이 시료를 가릅니다  → II 권의 뺄셈이 이것으로 두 시료를 구별합니다")
print("LB 가 라이브러리를 가릅니다 → 16강의 중복 표시가 이것을 봅니다")
print()
print("같은 시료를 여러 레인에서 읽었다면 ID 는 달라도 SM 은 같아야 합니다.")

🤔 Check yourself

BWA-MEM 과 minimap2 의 관계로 옳은 것은?

Answer 둘 다 Heng Li 가 만들었습니다. BWA-MEM 은 짧은 리드 DNA 정렬의 사실상 표준이고, minimap2 는 오류가 많은 긴 리드를 위해 minimizer 방식으로 만들어졌습니다. minimap2 도 -x sr 로 짧은 리드를 잘 붙이므로 노트북이 대체품으로 쓸 수 있는 것입니다.

🤔 Check yourself

minimizer 방식이 씨앗을 크게 줄이고도 민감도를 지키는 이유는?

🤔 Check yourself

정렬할 때 -R '@RG\tID:…\tSM:tumor…' 를 빠뜨리면?

Answer SM: 이 시료 이름이고 변이 부르기가 이것으로 종양과 정상을 구별합니다. GATK 는 아예 거부하고, bcftools 는 파일 이름으로 넘어가지만 병합 시 섞입니다. 정렬을 다시 돌리는 것은 비싸므로 처음에 붙입니다.