BioPlayground

🧬
목록으로

BWA-MEM 실무: SAM/BAM 태깅과 다중 샘플 정렬 오케스트레이션

BWA-MEM은 왜 10년 넘게 표준일까요. 시드 앤 익스텐드 원리부터 SAM/BAM 태그의 실무적 의미, read group 설정, 다중 샘플 병렬 오케스트레이션까지 하나로 정리합니다.

입문
|
20분
|
검증 완료 (2026-07-22)
short-read alignmentBWASAM/BAM
진행률0/52 (0%)

왜 아직도 BWA-MEM인가요

2013년 Heng Li가 BWA-MEM을 발표한 지 12년이 넘었습니다. 그동안 Minimap2·Bowtie2·HISAT2·STAR가 나왔고, DeepVariant는 아예 딥러닝을 얹었습니다. 그런데도 임상 유전체 파이프라인, GATK4 Best Practices, gnomAD, UK Biobank, All of Us — 이 모두가 여전히 BWA-MEM으로 정렬합니다. 왜일까요.

이유는 단순합니다. 잘 붙는 리드는 다 잘 붙이고, 안 붙는 리드는 안 붙는다고 정직하게 말한다. 이 정직성이 뒤 파이프라인의 통계 가정을 만족시킵니다. 오늘은 이 안정성이 어떤 설계 결정에서 오는지, 실무에서 무엇을 신경 써야 하는지 짚어봅시다.

이전 편들과의 연결 — Micro 티어의 재활용

BWA-MEM 안에서 우리가 M 티어에서 배운 것들이 이렇게 결합됩니다.

  • M13～M15: 참조 게놈을 접미사 배열(SA) → BWT → FM-index로 미리 인덱싱.
  • M10～M11: 리드를 k-mer 시드로 자르고 FM-index로 빠른 후보 위치 찾기 (BLAST 정신).
  • M06～M08: 후보 위치 주변에서 아핀 갭 페널티 국소 정렬 (Smith-Waterman 변형).

즉 BWA-MEM은 우리가 이미 유도한 세 알고리즘의 실무 조립체입니다. 이 관점을 갖고 옵션을 보면 훨씬 편해집니다.

인덱싱 — 한 번 하고 재사용

참조 게놈(hg38이면 3.1GB fa)을 다음 명령으로 인덱싱합니다.

bash
bwa index -a bwtsw hg38.fa

이 한 번의 작업이 5～6GB 인덱스 파일(.bwt, .pac, .ann, .amb, .sa)을 생성합니다. 인간 게놈이면 대략 30～60분. 한 번 만들면 모든 프로젝트에서 재활용하니 팀 서버의 공용 폴더에 두는 게 관습입니다.

BWA-MEM2를 쓴다면 별도 인덱스가 필요합니다.

bash
bwa-mem2 index hg38.fa

BWA-MEM2는 원조 BWA-MEM과 정확히 같은 결과를 내지만 벡터 명령어(AVX-512)를 사용해 1.7～3배 빠릅니다. 최근 클라우드 라이브러리는 대부분 BWA-MEM2 우선입니다.

기본 실행 — 한 샘플

가장 단순한 실무 커맨드입니다.

bash
bwa-mem2 mem -t 16 \
  -R "@RG\tID:S1_lane1\tSM:SAMPLE01\tLB:lib1\tPL:ILLUMINA\tPU:unit1" \
  hg38.fa clean_R1.fq.gz clean_R2.fq.gz \
  | samtools sort -@ 4 -o SAMPLE01.bam -
samtools index SAMPLE01.bam

세 부분으로 나눠 봅시다.

  • bwa-mem2 mem: 정렬 자체. -t는 스레드 수, -R은 read group 헤더.
  • samtools sort: SAM 스트림을 좌표 정렬된 BAM으로 변환.
  • samtools index: .bai 인덱스 생성 (IGV·GATK가 요구).

파이프(|)로 이어서 SAM을 디스크에 안 쓰는 게 실무 관습입니다. 100GB SAM을 파일로 두면 저장 낭비고 I/O 낭비입니다.

Read Group — 파이프라인의 이름표

-R 옵션의 그 문자열이 뭘까요. @RG 태그는 각 리드가 어느 실험에서 왔는지 알려주는 메타데이터입니다.

필드의미실무 예시
ID유일 식별자 (레인·플로우셀 단위)S1_lane1
SM샘플명 (환자·조직)SAMPLE01
LB라이브러리 (같은 SM에서 여러 라이브러리 가능)lib1
PL플랫폼ILLUMINA / PACBIO
PU플랫폼 단위 (실무에서 종종 생략 가능하나 GATK가 선호)run1_lane1

중요: GATK MarkDuplicates가 LB로 라이브러리 단위 중복을 판정합니다. HaplotypeCaller는 SM으로 샘플별 변이를 부릅니다. Read group을 대충 넣으면 뒤 파이프라인이 조용히 잘못된 결과를 냅니다. 임상 파이프라인 사고 절반이 이 read group 실수에서 옵니다.

SAM/BAM 태그 — 무엇을 어디서 확인할까요

정렬 결과 한 줄을 열어봅시다.

text
SRR..._1  99  chr1  100501  60  150M  =  100601  250  ACGT...  IIII...  NM:i:2  MD:Z:75A20T53  AS:i:145  XS:i:87

필드 순서(왼쪽에서): QNAME · FLAG · RNAME · POS · MAPQ · CIGAR · RNEXT · PNEXT · TLEN · SEQ · QUAL · 태그들.

실무에서 매번 보는 필드는 세 가지입니다.

1. MAPQ (mapping quality)

  • 정수 0～60. -10 log10(잘못된 위치 확률)이라는 정신은 Phred와 같습니다.
  • BWA-MEM 규약: 60 = 유일 최선 정렬, 0 = 여러 위치에 똑같이 잘 붙음 (multi-mapper).
  • GATK BQSR/HaplotypeCaller가 기본 MAPQ 20 이상만 사용합니다. MAPQ가 낮은 리드는 뒷단이 조용히 버립니다.

2. CIGAR

  • 리드가 참조에 어떻게 붙었는지 문자열로 표현.
  • M(매치·불일치 공통), I(리드 삽입), D(리드 삭제), S(soft-clip, 리드에 남되 정렬에서 제외), H(hard-clip, 리드 제거).
  • 예: 100M50S는 앞 100bp만 정렬되고 뒤 50bp는 어댑터 잔재로 잘렸다는 뜻.

3. FLAG

  • 12비트 비트마스크. 리드가 paired인지, first인지, unmapped인지, duplicate인지 등을 한 정수에 담습니다.
  • 자주 쓰는 값: 4=unmapped, 1024=PCR duplicate, 256=secondary, 2048=supplementary.
  • 팁: broadinstitute.github.io/picard/explain-flags.html 에서 정수를 넣으면 뜻을 풀어줍니다.

이 세 필드만 편안하게 읽을 수 있으면 파이프라인 디버깅의 80%가 해결됩니다.

Supplementary / chimeric alignment — 롱리드 시대의 함정

BWA-MEM의 특수 결과 중 실무 사고 자주 내는 게 2048 플래그(supplementary alignment)입니다.

  • 하나의 리드가 참조에서 두 곳에 나눠 붙습니다. 앞부분은 chr1, 뒷부분은 chr7에.
  • BWA는 이걸 "chimeric alignment"로 판정하고 리드를 두 줄로 뱉습니다. 하나는 primary, 다른 하나는 supplementary(2048 플래그).
  • 이 두 줄을 각각 별개 리드로 세어버리면 커버리지가 두 배로 뜨거나, 구조변이(융합 유전자) 신호를 놓칩니다.

실무 규칙: 커버리지 계산 전에 primary만 필터 (samtools view -F 2304, 즉 secondary+supplementary 제외). 반대로 SV 검출은 이 supplementary가 핵심 증거입니다. 용도가 다르면 필터도 다르다를 기억합시다.

옵션 서너 개만 진짜로 알아둡시다

BWA-MEM 매뉴얼은 옵션 30개 이상이지만 실무에서 만지는 건 셋뿐입니다.

  • -t N: 스레드. CPU 코어 수만큼.
  • -M: split 리드를 secondary로 마킹 (예전 GATK 호환). 최신 GATK4는 필요 없지만, 팀의 옛 파이프라인이면 유지.
  • -K 100000000: 배치 크기 고정. 스레드 수를 바꿔도 결과가 동일하도록 보장 (재현성).

정렬 정확도 자체는 기본값이 이미 극도로 잘 조율되어 있어 손대지 않습니다.

다중 샘플 오케스트레이션 — Bash 한 장

30개 샘플을 어떻게 돌릴까요. Snakemake/Nextflow를 배우기 전에는 이런 Bash 루프로 시작합니다.

bash
#!/bin/bash
set -euo pipefail
REF=/data/ref/hg38.fa
OUT=/data/aligned


for R1 in fastq/*_R1.fq.gz; do
    SAMPLE=$(basename "$R1" _R1.fq.gz)
    R2=fastq/${SAMPLE}_R2.fq.gz
    RG="@RG\tID:${SAMPLE}\tSM:${SAMPLE}\tLB:${SAMPLE}_lib\tPL:ILLUMINA"


    echo "[$(date +%T)] Aligning $SAMPLE"
    bwa-mem2 mem -t 16 -K 100000000 -R "$RG" "$REF" "$R1" "$R2" \
        | samtools sort -@ 4 -o "$OUT/${SAMPLE}.bam" -
    samtools index "$OUT/${SAMPLE}.bam"
done

이 스크립트에서 배울 것 세 가지.

  1. set -euo pipefail: 하나 실패하면 즉시 중단. 파이프라인은 침묵의 실패를 가장 경계합니다.
  2. 샘플명을 파일명에서 자동 추출.
  3. Read group을 샘플별로 자동 세팅.

이 뼈대에서 GNU Parallel이나 Slurm으로 옮겨가는 게 다음 단계입니다.

손 계산 — MAPQ의 감

BWA-MEM이 MAPQ를 어떻게 매기는지 대충 감을 잡아봅시다. 두 후보 위치의 정렬 점수가 AS=145, XS=87이라고 하면(SAM 태그 실제 예시), 최선과 차선의 차이가 58점. BWA는 이 차이를 확률로 변환해서 대략 MAPQ 60을 부여합니다. 반대로 AS=145, XS=143이면 두 위치가 거의 동점이라 MAPQ가 3～5 정도로 뚝 떨어집니다. 이게 반복서열이나 유사 유전자 부위에서 발생합니다.

  • 함의: MAPQ 60인 리드는 통계적으로 자신 있게 붙은 리드입니다.
  • 함의: MAPQ 0은 "여러 곳에 똑같이 잘 붙어서 어디라고 확정할 수 없음"입니다. 정보로는 유용하지만 변이 콜링에서는 대개 제외됩니다.

CS 매핑 — 시드 조회 + 국소 DP

  • FM-index 조회: 리드의 k-mer 시드를 참조에서 빠르게 찾는 부분. DryBench "FM-index 기본" 편(fm-index-fundamentals) 그대로입니다.
  • 국소 DP: 시드 확장 단계. Smith-Waterman(M07)과 아핀 갭(M08)의 산업용 튜닝.
  • 배치 처리: -K 100000000은 재현성을 위한 배치 크기 고정. 병렬 처리 시 결과 유일성을 지키는 CS 관례.

마무리

BWA-MEM은 12년째 표준입니다. 그 이유가 알고리즘의 세련됨보다 안정된 태그·MAPQ·chimeric 처리에 있음을 오늘 확인했습니다. 다음 편(S04)에서는 RNA-seq처럼 스플라이스가 있는 리드에 대해 STAR/HISAT2가 어떤 다른 결정을 내리는지 볼 겁니다. 그 뒤 S05에서 samtools로 BAM/CRAM을 정리하고 S06부터 GATK4 파이프라인에 진입합니다.

더 깊게 파고 싶다면

  • Li H. (2013), Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv:1303.3997. BWA-MEM 원 논문 — 시드 앤 익스텐드의 실무 튜닝이 자세합니다.
  • Vasimuddin M. et al. (2019), Efficient architecture-aware acceleration of BWA-MEM for multicore systems. IPDPS. BWA-MEM2의 SIMD 벡터화 논문.
  • Broad Institute BroadE 유튜브 — Read Alignment with BWA 30분 강의. 자막 완비. GATK4 파이프라인 관점에서 read group 실수 사례가 나옵니다.
  • Heng Li 개인 블로그: https://lh3.github.io/. 저자의 BWA-MEM 옵션 결정 후일담이 자주 올라옵니다.
  • SAMtools 공식 문서 (https://www.htslib.org/doc/): SAM/BAM/CRAM 규격의 원전.

한 샘플을 손으로 정렬해봐야 태그의 감이 옵니다. Galaxy EU에서 소형 FASTQ 하나 골라 BWA-MEM → samtools sort → IGV 시각화까지 눌러보고 오면 다음 편이 훨씬 수월합니다.