BWA-MEM·STAR로 리드 정렬하기 — FASTQ에서 BAM까지 (NGS Pipelines 2)

TL;DR — 시퀀서가 뱉어낸 수백만 개의 짧은 리드(read)는 그 자체로는 위치 정보가 없습니다. 각 리드가 참조유전체(reference genome)의 어디에서 온 것인지 찾아 붙이는 과정이 정렬(alignment)이고, 그 결과가 BAM (정렬 결과의 이진 포맷) 파일입니다. 이 글은 conda로 bwa·STAR·samtools를 설치하고, 작은 대장균 유전체로 인덱싱→정렬→정렬(sort)→인덱스→통계까지 직접 돌려 봅니다.

DNA (WGS·WES) 리드는 bwa mem으로, RNA-seq 리드는 인트론을 건너뛰어 붙일 줄 아는 STAR로 정렬합니다. 둘 다 결과는 SAM/BAM으로 모이고, samtools로 정렬·압축·인덱싱한 뒤 samtools flagstat으로 매핑률을 확인합니다.

정렬은 변이 호출(variant calling)·발현 정량 같은 거의 모든 후속 분석의 입구입니다. 여기서 만든 정렬 BAM이 다음 단계의 입력이 됩니다.

🔗 관련 글  ·  정렬의 개념·MAPQ는 서열 정렬(read alignment)이란?  ·  결과 파일의 11개 필드·FLAG·CIGAR는 SAM/BAM 파일이란?  ·  시퀀싱 세대 차이는 NGS란? 1·2·3세대 정리

이 글에서 만들 것

NGS 데이터 분석의 첫 갈림길은 정렬입니다. 시퀀서는 유전체를 잘게 쪼개 짧은 조각(리드)으로 읽어 내는데, 이 조각들은 순서도 위치도 뒤섞인 채로 FASTQ 파일에 담겨 나옵니다. 각 리드가 원래 유전체의 어느 좌표에서 왔는지를 되찾아 줄을 세우는 일, 그게 정렬입니다. 비유하자면 산산조각 난 책을 페이지 번호 없이 받아 들고, 원본(참조유전체)과 대조해 가며 각 조각이 몇 쪽 몇 번째 줄에 들어가는지 찾아 꽂는 작업입니다.

이 글에서는 공개·무료 데이터로 처음부터 끝까지 재현 가능한 정렬 파이프라인을 돌립니다. 대장균(E. coli) 유전체를 참조로 인덱싱하고, 시뮬레이션 리드를 bwa mem으로 정렬해 SAM (정렬 결과의 텍스트 포맷) 파일을 얻은 뒤, samtools로 BAM 변환·좌표 정렬·인덱싱을 거쳐 samtools flagstat으로 매핑 통계를 읽습니다. RNA-seq용 splice-aware 정렬 도구인 STAR는 왜·언제 쓰는지, 인덱스 생성과 정렬 명령을 따로 다룹니다. 마지막에 흔한 에러를 정리합니다.

# 핵심 흐름 한눈에 (DNA 정렬)
conda install -c bioconda bwa samtools          # ① 도구 설치
bwa index ref.fa                                # ② 참조유전체 인덱싱
bwa mem -t 8 ref.fa R1.fq.gz R2.fq.gz \
  | samtools sort -@ 8 -o sorted.bam -          # ③ 정렬 → 좌표정렬 BAM
samtools index sorted.bam                       # ④ BAM 인덱스(.bai)
samtools flagstat sorted.bam                    # ⑤ 매핑 통계

1. 사전 준비 — 설치·참조유전체·샘플 리드

도구는 conda(또는 더 빠른 mamba)로 한 번에 설치합니다. bioconda 채널에 정렬 3대 도구가 모두 올라와 있습니다.

# bioconda 채널 설정 (한 번만)
conda config --add channels bioconda
conda config --add channels conda-forge
conda config --set channel_priority strict

# 정렬 환경 한 번에 생성
conda create -n align -c conda-forge -c bioconda \
  bwa bwa-mem2 star samtools wgsim
conda activate align

2026년 6월 기준 현행 버전은 bwa 0.7.19, bwa-mem2 2.3, STAR 2.7.11b, samtools 1.23입니다. wgsim은 테스트 리드를 만들 시뮬레이터입니다. 예전엔 samtools에 딸려 있었지만 지금은 별도 도구(lh3/wgsim)라 따로 설치해야 합니다.

참조유전체는 대장균 K-12 MG1655 (~4.6 Mb)를 씁니다. 사람 유전체(~3.1 Gb)는 인덱싱에만 한참 걸리고 RAM도 많이 먹으니, 학습·재현용으로는 작은 세균 유전체가 제격입니다.

# 대장균 K-12 MG1655 참조유전체 (ASM584v2) 내려받기
curl -L -o ecoli.fa.gz \
  https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz
gunzip ecoli.fa.gz                              # → ecoli.fa (4,641,652 bp)

# 페어드엔드(paired-end) 리드 10만 쌍 시뮬레이션 (150 bp × 2)
wgsim -N 100000 -1 150 -2 150 ecoli.fa R1.fq R2.fq
gzip R1.fq R2.fq                                # → R1.fq.gz, R2.fq.gz

💡 실제 데이터라면 — wgsim 대신 시퀀싱 시설에서 받은 FASTQ를, ecoli.fa 대신 분석 대상 종의 참조유전체(사람이면 GRCh38)를 씁니다. 흐름과 명령은 똑같습니다. 정렬 전에 품질관리(QC)·어댑터 제거를 거치는 게 보통입니다.

2. 참조유전체 인덱싱 — bwa index

정렬은 리드 하나하나를 4.6 Mb(사람이면 31억 bp) 전체와 일일이 비교하지 않습니다. 그러면 너무 느리니까요. 대신 참조유전체를 미리 색인(index) 구조로 만들어 두고, 그 색인을 빠르게 뒤져 후보 위치를 좁힙니다. bwa는 BWT (Burrows-Wheeler Transform) 기반 FM-index를 만듭니다.

# 참조유전체 인덱스 생성 (한 번만 — 결과는 재사용)
bwa index ecoli.fa
# → ecoli.fa.amb, ecoli.fa.ann, ecoli.fa.bwt, ecoli.fa.pac, ecoli.fa.sa 생성

대장균이면 1초 안에 끝나지만, 사람 유전체는 한 시간 가까이 걸립니다. 인덱스는 한 번 만들면 같은 참조에 대해 계속 재사용하므로, 보통 분석 시작 때 한 번만 돌립니다.

📌 bwa-mem2를 쓴다면 — bwa-mem2는 classic bwa와 결과가 동일하면서 매핑 속도가 1.3~3.1배 빠른 후속 버전입니다. 다만 인덱스를 만들 때 참조 크기의 28배에 달하는 RAM이 필요합니다(사람 유전체 인덱싱은 큰 메모리 노드에서). 인덱스 파일 형식도 달라서 bwa 인덱스와 호환되지 않습니다. 명령은 bwa-mem2 index ecoli.fa처럼 거의 같습니다.

3. BWA-MEM으로 DNA 리드 정렬 — FASTQ → SAM → BAM

이제 리드를 정렬합니다. bwa mem은 짧은 리드(70 bp~수백 bp)를 참조에 붙이는 사실상의 표준 알고리즘으로, WGS (전장유전체)·WES (전장엑솜) 같은 DNA 시퀀싱에 씁니다. 핵심 옵션은 스레드 수(-t)와 리드 그룹(-R)입니다.

# BWA-MEM 정렬: 페어드엔드 FASTQ → SAM
bwa mem -t 8 \
  -R '@RG\tID:sample1\tSM:sample1\tLB:lib1\tPL:ILLUMINA\tPU:unit1' \
  ecoli.fa R1.fq.gz R2.fq.gz > aln.sam

-R로 넣는 리드 그룹(read group, @RG)은 "이 리드들이 어느 샘플·라이브러리·레인에서 왔는지"를 꼬리표로 붙이는 정보입니다. GATK 같은 변이 호출 도구가 이 태그를 요구하므로, 정렬 단계에서 미리 박아 두면 나중에 다시 손볼 일이 없습니다. 작은따옴표로 감싸야 셸이 \t를 그대로 넘기고, bwa가 이를 실제 탭으로 바꿔 줍니다. 주요 필드는 ID(식별자)·SM(샘플)·LB(라이브러리)·PL(플랫폼)·PU(플랫폼 유닛)입니다.

정렬 결과 SAM은 사람이 읽을 수 있는 텍스트지만 용량이 큽니다. 이를 이진 포맷인 BAM으로 바꾸고, 좌표 순으로 정렬하고, 인덱스(.bai)를 만들어야 뒤따르는 도구들이 빠르게 특정 영역만 꺼내 볼 수 있습니다. samtools가 이 셋을 담당합니다.

# SAM → BAM 변환
samtools view -b -o aln.bam aln.sam            # -b = BAM으로 출력

# 좌표 순 정렬 (sort)
samtools sort -@ 8 -o sorted.bam aln.bam       # -@ = 스레드 수

# BAM 인덱스 생성 → sorted.bam.bai
samtools index sorted.bam

실무에서는 SAM 중간 파일을 디스크에 쓰지 않고 파이프로 바로 흘려보내는 방식을 더 많이 씁니다. 한 줄로 정렬→좌표정렬→압축까지 끝납니다.

# 권장: 중간 SAM 없이 파이프로 한 번에
bwa mem -t 8 \
  -R '@RG\tID:sample1\tSM:sample1\tLB:lib1\tPL:ILLUMINA\tPU:unit1' \
  ecoli.fa R1.fq.gz R2.fq.gz \
  | samtools sort -@ 8 -o sorted.bam -          # 끝의 '-' = 표준입력
samtools index sorted.bam

⚠️ 예전 튜토리얼이 쓰던 samtools view -bS의 S는 지금은 무시됩니다(입력 포맷을 자동 감지). -b만 주면 됩니다.

4. 정렬 통계 — samtools flagstat과 MAPQ

정렬이 잘 됐는지는 samtools flagstat으로 확인합니다. FLAG 비트를 집계해 전체 리드 수, 매핑된 비율, 제대로 짝지어진(properly paired) 비율, 중복(duplicates) 등을 한눈에 보여 줍니다.

samtools flagstat sorted.bam

200000 + 0 in total (QC-passed reads + QC-failed reads)
200000 + 0 primary
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
199876 + 0 mapped (99.94% : N/A)
199876 + 0 primary mapped (99.94% : N/A)
200000 + 0 paired in sequencing
100000 + 0 read1
100000 + 0 read2
199654 + 0 properly paired (99.83% : N/A)
199810 + 0 with itself and mate mapped
66 + 0 singletons (0.03% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

읽는 법은 이렇습니다. mapped 비율은 전체 리드 중 어디든 붙은 비율로, 좋은 데이터·올바른 참조면 보통 95% 이상입니다. 이 값이 낮으면 참조유전체를 잘못 골랐거나 오염(contamination)을 의심합니다. properly paired는 두 짝(R1·R2)이 기대한 방향·거리로 함께 붙은 비율로, 라이브러리 품질의 지표입니다. secondary/supplementary는 한 리드가 여러 곳에 붙은 보조 정렬, duplicates는 PCR 증폭 등으로 생긴 중복인데 여기서는 0으로 나옵니다(중복 표시를 아직 안 했으니까요).

각 리드에는 MAPQ (mapping quality) 값이 붙습니다. 정렬 위치가 틀렸을 확률을 Phred 척도로 옮긴 값으로, MAPQ = −10 × log₁₀(P(위치가 틀림))입니다. MAPQ 30이면 1/1,000, 40이면 1/10,000 확률로 틀렸다는 뜻이죠. bwa mem은 최대 60까지 매기며, 어디에 붙여야 할지 애매한(여러 곳에 똑같이 잘 붙는) 리드는 MAPQ 0을 받습니다. 다중 매핑 리드를 걸러 내려면 samtools view의 -q 옵션을 씁니다.

# MAPQ 30 미만 리드 제거 (다중 매핑·저신뢰 정렬 걸러내기)
samtools view -b -q 30 sorted.bam > filtered.bam

5. STAR로 RNA-seq 리드 정렬 — splice-aware 정렬

RNA-seq 리드를 bwa mem으로 정렬하면 안 됩니다. 성숙한 mRNA에는 인트론(intron)이 잘려 나가고 엑손(exon)만 이어져 있어서, 한 리드가 유전체 상에서 멀리 떨어진 두 엑손에 걸쳐 있을 수 있기 때문입니다. bwa는 이런 인트론 건너뛰기(splice)를 모릅니다. 그래서 RNA-seq에는 splice-aware 정렬 도구인 STAR를 씁니다. STAR는 인트론을 가로지른 리드를 끊어서 양쪽 엑손에 정확히 붙입니다.

STAR도 먼저 참조유전체 인덱스를 만듭니다. 이때 유전자 주석(annotation) GTF 파일을 함께 줘서 알려진 스플라이스 접합부 위치를 미리 알려 줍니다.

# ① STAR 인덱스 생성 (genomeGenerate)
STAR --runMode genomeGenerate \
     --genomeDir ./star_index \
     --genomeFastaFiles ecoli.fa \
     --sjdbGTFfile annotation.gtf \
     --sjdbOverhang 149 \
     --genomeSAindexNbases 10 \
     --runThreadN 8

--sjdbOverhang은 접합부 양쪽에 둘 서열 길이로, 리드 길이에서 1을 뺀 값이 이상적입니다(150 bp 리드면 149). --genomeSAindexNbases는 작은 유전체에서 반드시 줄여 줘야 하는 값입니다. 사람 유전체면 기본값 14로 두지만, 작은 유전체에 14를 쓰면 STAR가 경고를 내고 정렬 단계에서 세그폴트(seg-fault)가 날 수 있습니다. 권장값은 min(14, log₂(L)/2 − 1)로, 대장균(4.64 Mb)이면 10입니다(내림). 세균처럼 스플라이싱이 없는 유전체라면 GTF 관련 옵션은 빼도 됩니다.

# ② STAR 정렬 (alignReads) — 좌표 정렬 BAM으로 바로 출력
STAR --runMode alignReads \
     --genomeDir ./star_index \
     --readFilesIn R1.fq.gz R2.fq.gz \
     --readFilesCommand zcat \
     --outSAMtype BAM SortedByCoordinate \
     --outSAMattrRGline ID:sample1 SM:sample1 LB:lib1 PL:ILLUMINA PU:unit1 \
     --runThreadN 8 \
     --outFileNamePrefix sample1_
# → sample1_Aligned.sortedByCoord.out.bam 생성
samtools index sample1_Aligned.sortedByCoord.out.bam

--outSAMtype BAM SortedByCoordinate를 주면 STAR가 좌표 정렬된 BAM을 바로 뱉어 주므로 samtools sort를 따로 돌릴 필요가 없습니다. gzip된 FASTQ는 --readFilesCommand zcat으로 풀어 읽습니다(단, genomeGenerate에 넣는 참조 FASTA는 압축을 풀어 둬야 합니다). 한 가지 함정은 리드 그룹 문법이 bwa와 다르다는 점입니다. STAR의 --outSAMattrRGline은 @RG 접두어도 \t도 없이 공백으로 태그를 구분합니다. 복사·붙여넣기로 가장 자주 틀리는 부분이니 주의하세요.

그림 1. 정렬의 개념 — 위치가 뒤섞인 FASTQ 리드를 참조유전체에 대조해 원래 좌표에 줄 세우고, 결과를 BAM으로 저장한다.

그림 2. BWA-MEM vs STAR — DNA(WGS·WES)는 BWA-MEM으로, 인트론을 건너뛰어야 하는 RNA-seq는 splice-aware인 STAR로 정렬한다.

6. 자주 발생하는 에러 & 해결법

• fail to locate the index files (bwa) — 인덱스를 안 만들었거나 참조 경로(prefix)가 틀렸습니다. bwa index ref.fa를 먼저 돌리고, bwa mem에 같은 prefix를 주세요. bwa는 ref.fa.{amb,ann,bwt,pac,sa}를 찾습니다. bwa 인덱스와 bwa-mem2 인덱스는 파일 형식이 달라 서로 못 씁니다.
• STAR genomeGenerate가 Killed (메모리 부족) — 사람 유전체 인덱스 생성에는 32 GB 이상의 RAM이 필요합니다. 부족하면 OS가 프로세스를 죽이는데, STAR 에러 없이 그냥 Killed만 떠서 헷갈립니다(dmesg로 확인). RAM이 넉넉한 노드에서 돌리거나 --limitGenomeGenerateRAM으로 상한을 거세요.
• STAR 작은 유전체 경고·세그폴트 — 작은 유전체에 --genomeSAindexNbases 기본값(14)을 쓰면 경고가 뜨고 정렬에서 세그폴트가 날 수 있습니다. min(14, log₂(L)/2 − 1)로 줄이세요(대장균이면 10).
• samtools sort 메모리·임시공간 부족 — sort는 스레드당 기본 768 MiB를 쓰므로 -@를 키우면 총 메모리가 함께 늡니다(-@ 16 -m 4G는 OOM 위험). 임시 파일은 -T로 여유 있는 볼륨을 가리키게 하세요: samtools sort -@ 8 -m 1G -T /scratch/srt -o sorted.bam in.bam.
• 리드 그룹 누락 → GATK 에러 — @RG가 없으면 GATK가 변이 호출 단계에서 멈춥니다. 정렬 때 -R(또는 STAR --outSAMattrRGline)로 넣거나, 사후에 picard AddOrReplaceReadGroups로 보강하세요.
• bwa-mem2 index가 RAM을 너무 먹음 — 인덱스 생성에는 참조 크기의 28배 RAM이 듭니다(매핑 자체는 ~10 GB로 가볍습니다). 큰 RAM 노드에서 한 번 만들어 두고 매핑은 어디서든 하면 됩니다.

7. 확장 / 다음 단계

정렬 BAM은 거의 모든 후속 분석의 출발점입니다. DNA 정렬 BAM에서는 변이 호출(variant calling)로 넘어갑니다. 중복 표시(MarkDuplicates)와 염기 품질 재보정(BQSR)을 거친 뒤 GATK·DeepVariant로 SNP·INDEL을 찾아 VCF 파일로 내보내는 흐름입니다. 정렬 단계에서 리드 그룹을 미리 넣어 둔 이유가 여기서 드러납니다 — 변이 호출 도구가 그 태그를 요구하니까요. 개념과 도구 비교는 변이 호출(variant calling)이란?에 정리해 두었습니다.

RNA-seq 정렬 BAM에서는 발현 정량으로 이어집니다. 유전자별로 붙은 리드 수를 세거나(featureCounts·HTSeq), 전사체 수준에서 정량해(salmon) 차등발현 분석으로 넘어갑니다.

그림 3. samtools flagstat 통계 읽기 — 매핑률·properly paired·중복 등 핵심 지표를 막대로 본다.

🎓 다음 학습
• (개념 복습) 서열 정렬(read alignment)이란? — BWA·minimap2·MAPQ의 원리
• (결과 포맷) SAM/BAM 파일이란? — 11개 필드·FLAG·CIGAR 완전 정리
• (다음 단계) 변이 호출(variant calling)이란? — BAM에서 VCF로, GATK·DeepVariant
• (시퀀싱 기초) NGS란? 1·2·3세대 시퀀싱 — Sanger·Illumina·Nanopore

References

1. Li, H. & Durbin, R. (2009). Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics, 25(14), 1754–1760. doi:10.1093/bioinformatics/btp324 (BWA 원논문)
2. Li, H. (2013). Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv:1303.3997. arXiv:1303.3997 (BWA-MEM)
3. Vasimuddin, M. et al. (2019). Efficient architecture-aware acceleration of BWA-MEM for multicore systems. IEEE IPDPS. bwa-mem2 GitHub
4. Dobin, A. et al. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics, 29(1), 15–21. doi:10.1093/bioinformatics/bts635 (STAR 원논문)
5. Danecek, P. et al. (2021). Twelve years of SAMtools and BCFtools. GigaScience, 10(2), giab008. doi:10.1093/gigascience/giab008 (samtools)
6. SAMtools / Htslib 공식 문서. htslib.org/doc · STAR 매뉴얼. github.com/alexdobin/STAR

Pipette & Pipeline · A bio portfolio journal