samtools·bcftools 실전 — BAM 정렬·인덱싱과 VCF 변이 필터링
TL;DR — NGS 파이프라인에서 정렬 결과(SAM/BAM)와 변이 결과(VCF)를 다루는 두 표준 도구가
samtools(SAM/BAM/CRAM 처리)와bcftools(VCF/BCF 처리)입니다. 이 글은 samtools에 동봉된 예제(toy.sam·toy.fa)로, SAM→정렬된 BAM→인덱스(.bai)→커버리지 확인→bcftools mpileup+call로 변이 호출→QUAL·DP로 필터링까지 한 줄씩 따라갑니다. 모든 명령·플래그에 한글 주석을 달고,flagstat·coverage출력 예시까지 그대로 보여 재현 가능하게 구성했습니다.samtools의 핵심 흐름은 단순합니다.
view로 형식을 바꾸거나 FLAG·MAPQ로 거르고,sort로 좌표순 정렬하고,index로 임의 좌표 접근을 위한 색인을 만든 뒤,flagstat·coverage·depth로 정렬 품질과 커버리지를 확인합니다. bcftools는 그 BAM에서mpileup+call로 변이를 뽑아(GATK의 가벼운 대안),filter표현식으로 저품질 변이를 걷어냅니다.재현성을 위해 예제는 어디서 돌려도 결과가 같은
toy.sam을 씁니다. 초보가 가장 많이 밟는 지뢰 — 정렬 안 된 BAM을 인덱싱하다 나는 오류,bcftools mpileup의-f참조 누락,@SQ헤더 불일치 — 도 4절에서 실제 에러 메시지와 함께 짚습니다.
🔗 관련 글 · SAM/BAM의 11개 필드·FLAG·CIGAR를 아직 안 봤다면 SAM/BAM 파일이란?을 먼저 읽는 걸 권합니다 · 변이 결과 포맷은 VCF 파일이란?에서 8개 컬럼·INFO·GT를 정리했습니다 · BAM을 만들기 직전 단계인 정렬은 BWA-MEM·STAR로 리드 정렬하기에서, 더 무거운 변이 호출기는 GATK으로 변이 호출하기에서 이어집니다
이 글에서 만들 것
시퀀싱 파이프라인은 대체로 이렇게 흐릅니다. FASTQ (원시 리드)를 품질검사·트리밍하고, 참조 유전체에 정렬해 SAM/BAM을 얻고, 그 정렬을 정리한 뒤 변이를 호출해 VCF를 만듭니다. 이 가운데 정렬 결과를 다루는 도구가 samtools, 변이 결과를 다루는 도구가 bcftools입니다. 정렬 자체(FASTQ→SAM)는 BWA·STAR 편에서 다뤘으니, 이 글은 그다음 — SAM을 정렬된 BAM으로 만들고, 커버리지를 확인하고, 변이를 뽑아 필터링하는 구간 — 을 처음부터 끝까지 만듭니다.
samtools (Heng Li 외, 2009~)는 SAM (Sequence Alignment/Map)·BAM (SAM의 이진 압축본)·CRAM을 읽고 쓰고 거르는 표준 도구입니다. bcftools는 그 짝으로, 변이 정보 포맷인 VCF (Variant Call Format)·BCF (VCF의 이진본)를 다룹니다. 둘 다 htslib이라는 같은 C 라이브러리 위에 얹혀 있어, 2021년 GigaScience 논문 "Twelve years of SAMtools and BCFtools"로 함께 정리됐습니다.
구체적으로 이 글은 samtools 예제 파일 하나로 ① view로 SAM을 BAM으로 바꾸고 FLAG·MAPQ로 거르고, ② sort로 좌표순 정렬한 뒤 index로 .bai 색인을 만들고, ③ flagstat·coverage·depth로 정렬 통계와 커버리지를 보고, ④ faidx로 참조를 색인한 뒤 bcftools mpileup+call로 변이를 호출하고, ⑤ bcftools filter·norm·query로 변이를 걸러 원하는 컬럼만 뽑습니다. 명령은 samtools·bcftools 1.21+ 기준입니다.
# 핵심 흐름 한눈에 (toy 예제 기준)
samtools view -b toy.sam | samtools sort -o toy.sorted.bam # ① SAM → 정렬된 BAM
samtools index toy.sorted.bam # ② .bai 색인 생성
samtools flagstat toy.sorted.bam # ③ 정렬 통계
samtools faidx toy.fa # ④ 참조 .fai 색인
bcftools mpileup -f toy.fa toy.sorted.bam \
| bcftools call -mv -Oz -o calls.vcf.gz # ⑤ 변이 호출 (VCF)
bcftools filter -e 'QUAL<20 || INFO/DP<10' calls.vcf.gz # ⑥ 저품질 변이 필터
1. 사전 준비 — 설치와 예제 데이터
samtools와 bcftools를 설치합니다. 대부분 환경에서 conda (bioconda 채널)가 가장 편하고, 데비안·우분투 계열이면 apt로도 됩니다.
# conda (bioconda) — 권장, 버전 고정이 쉬움
conda install -c bioconda -c conda-forge samtools bcftools
# 또는 데비안/우분투
sudo apt-get update && sudo apt-get install samtools bcftools# 설치 확인 (버전 출력)
samtools --version | head -1
bcftools --version | head -1samtools 1.21
bcftools 1.21이제 실습 데이터를 받습니다. 재현이 쉽도록, 어디서 돌려도 결과가 같은 samtools 공식 예제 toy.sam·toy.fa를 씁니다. 참조 서열 2개(ref 45 bp·ref2 40 bp)와 정렬 레코드 12개(ref에 r001~r004 6개·ref2에 x1~x6 6개)로 이루어진 아주 작은 데이터입니다.
# samtools 저장소에서 예제 파일 내려받기
wget https://raw.githubusercontent.com/samtools/samtools/develop/examples/toy.sam
wget https://raw.githubusercontent.com/samtools/samtools/develop/examples/toy.fa# 내용 확인 — 헤더(@SQ 2줄) + 정렬 레코드 12개
cat toy.sam@SQ SN:ref LN:45
@SQ SN:ref2 LN:40
r001 163 ref 7 30 8M4I4M1D3M = 37 39 TTAGATAAAGAGGATACTG * XX:B:S,12561,2,20,112
r002 0 ref 9 30 1S2I6M1P1I1P1I4M2I * 0 0 AAAAGATAAGGGATAAA *
r003 0 ref 9 30 5H6M * 0 0 AGCTAA *
r004 0 ref 16 30 6M14N1I5M * 0 0 ATAGCTCTCAGC *
r003 16 ref 29 30 6H5M * 0 0 TAGGC *
r001 83 ref 37 30 9M = 7 -39 CAGCGCCAT *
x1 0 ref2 1 30 20M * 0 0 aggttttataaaacaaataa ????????????????????
x2 0 ref2 2 30 21M * 0 0 ggttttataaaacaaataatt ?????????????????????
x3 0 ref2 6 30 9M4I13M * 0 0 ttataaaacAAATaattaagtctaca ??????????????????????????
x4 0 ref2 10 30 25M * 0 0 CaaaTaattaagtctacagagcaac ?????????????????????????
x5 0 ref2 12 30 24M * 0 0 aaTaattaagtctacagagcaact ????????????????????????
x6 0 ref2 14 30 23M * 0 0 Taattaagtctacagagcaacta ???????????????????????각 줄은 리드 하나의 정렬입니다. 왼쪽부터 리드 이름(QNAME)·FLAG·참조 이름(RNAME)·위치(POS)·MAPQ·CIGAR 순이고, 이 11개 필드의 의미는 SAM/BAM 파일 편에서 자세히 풀었습니다. 여기서 눈여겨볼 건 두 번째 열 FLAG입니다. FLAG는 여러 상태를 비트로 겹쳐 담은 정수인데, 이걸 이해하면 samtools view의 필터가 한결 명확해집니다.
2. samtools view — 형식 변환·FLAG·MAPQ 필터
SAM FLAG — 상태를 비트로 겹쳐 담은 정수
FLAG는 "이 리드가 짝이 있나·정방향인가·중복인가" 같은 여러 예/아니오를 2의 거듭제곱 비트로 합친 값입니다. 예컨대 FLAG가 163이면 1 (짝 있음) + 2 (제대로 짝지음) + 32 (짝이 역방향) + 128 (두 번째 리드)의 합이죠. 자주 쓰는 비트는 이렇습니다.
| 값(10진) | 16진 | 의미 |
|---|---|---|
| 1 | 0x1 | 리드에 짝이 있음(paired) |
| 4 | 0x4 | 정렬 안 됨(unmapped) |
| 16 | 0x10 | 역가닥에 정렬(reverse strand) |
| 256 | 0x100 | 이차 정렬(secondary) |
| 1024 | 0x400 | PCR·광학 중복(duplicate) |
| 2048 | 0x800 | 보충 정렬(supplementary) |
숫자를 외울 필요는 없습니다. samtools flags 서브명령이 비트↔이름을 바꿔 주거든요.
# FLAG 값 → 어떤 비트가 켜졌는지 이름으로 풀기
samtools flags 0x9040x904 2308 UNMAP,SECONDARY,SUPPLEMENTARY0x904(=2308)는 unmapped·secondary·supplementary 세 비트를 합친 값입니다. 이걸 -F(제외)에 넣으면 "정렬 실패·이차·보충 정렬을 빼고 진짜 일차 정렬만" 남길 수 있죠. 실전에서 자주 쓰는 관용구입니다.
view — 형식 변환과 필터
samtools view는 SAM↔BAM 형식 변환과 필터를 동시에 하는 만능 명령입니다. 자주 쓰는 옵션은 -b(BAM 출력), -h(헤더 포함), -q(MAPQ 문턱), -f(비트 요구), -F(비트 제외), -c(개수만)입니다.
# SAM → BAM 변환 (-b), 헤더 포함(-h는 BAM에선 기본 포함)
samtools view -b toy.sam -o toy.bam
# MAPQ ≥ 20 이고, unmapped·secondary·supplementary(0x904) 제외한 리드 개수만 세기
samtools view -c -q 20 -F 0x904 toy.bam12toy 데이터의 정렬 12개는 모두 MAPQ 30이고 이차·보충 정렬이 없어, 이 조건에 12개가 다 걸립니다. 문턱을 -q 60으로 올리면 0개가 되죠. -c는 레코드를 출력하지 않고 개수만 돌려줘, "필터 조건에 몇 개가 걸리나"를 빠르게 볼 때 씁니다. -q 20은 MAPQ (정렬 신뢰도, 값이 클수록 그 위치가 유일할 가능성이 높음)가 20 미만인 리드를 버립니다 — MAPQ 개념은 read alignment 편에서 다뤘습니다.
# 특정 좌표 구간만 잘라내기 (색인이 있어야 동작 — 3절에서 index 생성 후)
samtools view -b toy.sorted.bam ref:1-20 -o region.bamview에 참조이름:시작-끝 형식을 주면 그 구간에 걸친 리드만 뽑습니다. 다만 이 임의 좌표 접근은 BAM이 좌표순으로 정렬되고 색인(.bai)이 있을 때만 동작합니다. 그래서 다음 순서가 sort → index입니다.
3. samtools sort·index·flagstat·coverage — 정렬·색인·통계
sort — 좌표순 정렬
정렬기(BWA 등)가 내놓은 BAM은 리드가 읽힌 순서 그대로라, 참조 좌표순이 아닙니다. 대부분의 하류 도구(변이 호출·색인·구간 추출)는 좌표순 정렬을 전제하므로, samtools sort로 먼저 정렬합니다.
# 좌표순 정렬 (기본값) → toy.sorted.bam
samtools sort toy.bam -o toy.sorted.bamsort는 기본적으로 참조 위치(leftmost coordinate)순으로 정렬합니다. 리드 이름순이 필요하면(예: fixmate·markdup 전처리) -n을 줍니다. 큰 파일이면 -@로 스레드를, -m으로 스레드당 메모리를 늘려 속도를 올릴 수 있습니다.
index — 임의 좌표 접근을 위한 색인
samtools index는 정렬된 BAM에 .bai 색인을 붙여, 파일 전체를 읽지 않고도 특정 좌표로 바로 점프할 수 있게 합니다. IGV 같은 뷰어나 구간 추출이 이 색인에 기댑니다.
# .bai 색인 생성 → toy.sorted.bam.bai
samtools index toy.sorted.bam여기서 흔한 함정 하나. 정렬 안 된 BAM을 인덱싱하려 하면 오류가 납니다.
# (오류 예시) 정렬 안 된 toy.bam을 그대로 인덱싱하면
samtools index toy.bamsamtools index: "toy.bam" is not sorted or is malformed. Coordinate-sorted files are required for indexing.색인은 좌표순 정렬을 반드시 전제하므로, index 전에 sort가 와야 합니다. 이 순서는 4절에서 다시 정리합니다.
flagstat — 정렬 통계 한눈에
samtools flagstat은 FLAG를 집계해 "전체 몇 개 중 몇 개가 정렬됐고, 짝지음 비율은 얼마인가"를 한 화면에 보여 줍니다. 실제 대규모 BAM에서는 이런 모양입니다(paired-end, 약 100만 리드 예시).
samtools flagstat sample.sorted.bam1005678 + 0 in total (QC-passed reads + QC-failed reads)
1000000 + 0 primary
3210 + 0 secondary
2468 + 0 supplementary
18234 + 0 duplicates
18234 + 0 primary duplicates
998765 + 0 mapped (99.31% : N/A)
993087 + 0 primary mapped (99.31% : N/A)
1000000 + 0 paired in sequencing
500000 + 0 read1
500000 + 0 read2
985432 + 0 properly paired (98.54% : N/A)
991200 + 0 with itself and mate mapped
1887 + 0 singletons (0.19% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)각 줄은 QC통과 + QC실패 형식입니다. 가장 먼저 보는 값은 mapped 비율(정렬률, 여기선 99.31%)과 properly paired 비율(제대로 짝지어 정렬된 비율, 98.54%)입니다. 정렬률이 유난히 낮거나 properly paired가 낮으면 참조 불일치·오염·라이브러리 문제를 의심합니다. duplicates가 많으면 PCR 중복 표시(markdup)를 점검하죠.
coverage·depth — 얼마나 깊게 읽혔나
변이를 믿으려면 그 자리가 충분히 깊게 읽혔는지(커버리지)를 봐야 합니다. samtools coverage는 참조 서열마다 요약을, samtools depth는 위치마다 깊이를 줍니다. toy 데이터는 두 참조 모두에 리드가 걸려 있어(ref r001~r004·ref2 x1~x6) 둘 다 요약이 나옵니다.
# 참조 서열별 커버리지 요약 (한 줄에 서열 하나)
samtools coverage toy.sorted.bam#rname startpos endpos numreads covbases coverage meandepth meanbaseq meanmapq
ref 1 45 6 39 86.67 2.0 40 30
ref2 1 40 6 36 90 3.32 40 30coverage 컬럼은 그 서열에서 1× 이상 덮인 염기 비율(%), meandepth는 평균 깊이입니다. 실제 데이터라면 이 표에서 커버리지가 낮은 구간을 먼저 살펴, 변이를 믿어도 되는 자리인지 판단합니다. -m을 주면 같은 정보를 ASCII 히스토그램으로 그려 줍니다.
# 위치별 깊이 (CHROM POS DEPTH), -a는 깊이 0 위치도 포함
samtools depth toy.sorted.bam | head -5ref 7 1
ref 8 1
ref 9 2
ref 10 2
ref 11 2depth는 위치마다 세 컬럼(참조 이름·1-based 위치·깊이)을 내놓습니다. 특정 유전자·엑손의 깊이 프로파일을 그리거나, 최소 깊이 미달 구간을 찾을 때 씁니다.
4. bcftools — mpileup·call로 변이 호출하고 filter로 거르기
이제 정렬된 BAM에서 변이를 뽑습니다. bcftools는 두 단계로 나뉩니다. mpileup이 각 위치의 염기 증거를 모아 유전형 가능도(genotype likelihood)를 계산하고, call이 그걸 받아 실제 변이를 판정합니다.
faidx — 참조 색인 (변이 호출의 전제)
bcftools mpileup은 참조 FASTA를 -f로 반드시 받아야 합니다. 참조가 있어야 어느 자리가 참조와 다른지 알 수 있으니까요. 그리고 그 참조는 samtools faidx로 미리 색인(.fai)해 둬야 임의 좌표를 빠르게 읽습니다.
# 참조 FASTA 색인 → toy.fa.fai 생성
samtools faidx toy.fa
cat toy.fa.fairef 45 5 45 46
ref2 40 57 40 41.fai는 다섯 컬럼(이름·길이·바이트 오프셋·줄당 염기수·줄당 바이트수)입니다. toy.fa는 서열이 한 줄이라 줄당 염기수가 길이와 같죠. 이 색인 덕에 samtools faidx toy.fa ref:1-10처럼 참조의 특정 구간만 즉시 뽑을 수도 있습니다.
mpileup + call — 변이 호출
# mpileup: 위치별 염기 증거 → 유전형 가능도(VCF/BCF)
# -f: 참조 FASTA (필수)
# call: 가능도 → 변이 판정
# -m: 다대립 호출기(권장 기본) · -v: 변이 자리만 출력 · -Oz: gzip VCF 출력
bcftools mpileup -f toy.fa toy.sorted.bam \
| bcftools call -mv -Oz -o calls.vcf.gztoy 데이터는 리드가 적어 변이가 거의 안 나오지만, 명령의 흐름과 출력 구조는 실제와 같습니다. 아래는 실제 소형 BAM (예: 1000 Genomes 저커버리지 BAM의 한 구간)에서 얻은 것과 같은 형태의 VCF입니다(대표 예시).
# 결과 VCF 훑어보기 (헤더 제외 데이터 줄 앞부분)
bcftools view calls.vcf.gz | grep -v '^##' | head -4#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT sample
20 1001420 . G A 221.4 . DP=34;MQ=59 GT:PL 1/1:255,102,0
20 1004201 . C T 18.7 . DP=6;MQ=42 GT:PL 0/1:52,0,71
20 1008833 . T TA 164.2 . INDEL;DP=27 GT:PL 0/1:198,0,201VCF의 8개 고정 컬럼(CHROM·POS·ID·REF·ALT·QUAL·FILTER·INFO)에 이어, FORMAT과 샘플별 값이 붙습니다. QUAL은 변이가 진짜일 확률의 Phred 점수(클수록 신뢰), INFO의 DP는 그 자리 총 깊이, FORMAT의 GT는 유전형(0/1 이형접합·1/1 동형접합 변이)입니다. 이 컬럼들의 의미는 VCF 파일 편에서 정리했습니다. 위에서 두 번째 변이(POS 1004201)는 QUAL 18.7·DP 6으로, 다음 단계에서 걸러질 후보입니다.
filter — 표현식으로 저품질 변이 제거
bcftools filter는 -e(제외)·-i(포함) 표현식으로 변이를 거릅니다. 표현식에는 QUAL·INFO 필드를 이름으로 참조할 수 있습니다.
# QUAL이 20 미만 '또는' 깊이(DP)가 10 미만인 변이를 제외
bcftools filter -e 'QUAL<20 || INFO/DP<10' calls.vcf.gz -Oz -o filtered.vcf.gz
# 걸러진 뒤 남은 변이 확인
bcftools view filtered.vcf.gz | grep -v '^##'#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT sample
20 1001420 . G A 221.4 . DP=34;MQ=59 GT:PL 1/1:255,102,0
20 1008833 . T TA 164.2 . INDEL;DP=27 GT:PL 0/1:198,0,201QUAL<20 || INFO/DP<10 표현식이 QUAL 18.7·DP 6이던 POS 1004201을 정확히 걸러, 신뢰할 만한 변이 둘만 남았습니다. 실전에서는 변이를 삭제하는 대신, -s LowQual로 FILTER 컬럼에 이름표만 붙여(soft filter) 나중에 되살릴 수 있게 하는 방식도 많이 씁니다.
# (대안) 삭제 대신 FILTER 컬럼에 'LowQual' 표시만 (soft filter)
bcftools filter -s LowQual -e 'QUAL<20 || INFO/DP<10' calls.vcf.gz -Oz -o soft.vcf.gznorm·stats·query — 정규화·요약·추출
변이를 비교·병합하려면 표현을 통일해야 합니다. bcftools norm은 인델을 왼쪽 정렬(left-align)하고, 한 줄에 뭉친 다대립 변이를 -m -로 이대립씩 쪼갭니다.
# 인델 정규화(left-align) + 다대립 → 이대립 분해 (-m -)
bcftools norm -f toy.fa -m - filtered.vcf.gz -Oz -o norm.vcf.gz# 변이 요약 통계 (SNP·인델 개수 등)
bcftools stats norm.vcf.gz | grep '^SN'SN 0 number of samples: 1
SN 0 number of records: 2
SN 0 number of SNPs: 1
SN 0 number of indels: 1
SN 0 number of multiallelic sites: 0bcftools query는 원하는 컬럼만 뽑아 표로 만들 때 씁니다. INFO·고정 컬럼은 %필드로, 샘플별 FORMAT 값은 대괄호 [...]로 감쌉니다.
# 원하는 컬럼만 탭 구분으로 추출 ([%GT]는 샘플별 유전형)
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t[%GT]\n' norm.vcf.gz20 1001420 G A 221.4 1/1
20 1008833 T TA 164.2 0/1\t는 탭, \n은 줄바꿈입니다. 이렇게 뽑은 표는 그대로 R·pandas로 넘겨 하류 분석에 쓰기 좋습니다. query는 VCF를 사람이·스크립트가 다루기 쉬운 평범한 표로 바꿔 주는, 파이프라인 마지막의 단골 도구입니다.
QUAL<20 || DP<10 같은 표현식으로 저품질 변이를 걸러 신뢰할 만한 콜셋만 남긴다.자주 발생하는 에러 & 해결
- 정렬 안 된 BAM 인덱싱 오류 — 가장 흔한 함정 —
samtools index가"... is not sorted or is malformed. Coordinate-sorted files are required for indexing."를 뱉으면, BAM이 좌표순이 아니라는 뜻입니다. 정렬기 출력이나 이름순 정렬(-n) 결과를 그대로 인덱싱하면 이 오류가 납니다. 반드시samtools sort(좌표순)를 먼저 돌린 뒤index합니다. 순서는 view → sort → index로 기억하세요. bcftools mpileup의-f참조 누락 —mpileup에-f ref.fa를 빼면 변이를 계산할 기준이 없어 실패합니다. 참조 FASTA는 필수이며, 그 참조는samtools faidx ref.fa로 미리 색인(.fai)돼 있어야 합니다(없으면 mpileup이 자동 생성하려 하지만, 쓰기 권한이 없는 경로면 실패). 정렬에 쓴 참조와 같은 FASTA를 써야 좌표가 맞습니다.@SQ헤더 불일치 — 참조가 다를 때 — BAM 헤더의@SQ(참조 이름·길이)와 변이 호출에 준 참조 FASTA가 다르면, 좌표가 어긋나거나the sequence ... not found류의 오류가 납니다. 정렬·변이 호출·주석 달기까지 같은 참조 유전체 버전(예: GRCh38)을 일관되게 써야 합니다.chr1vs1처럼 이름 규칙만 달라도 문제가 되니, 헤더를samtools view -H로 먼저 확인하세요.- 인덱스가 낡음(BAM보다 오래됨) — BAM을 다시 만들었는데
.bai가 옛것이면 뷰어·구간 추출이 어긋납니다. BAM을 갱신할 때마다samtools index를 다시 돌려 색인을 최신으로 유지합니다. 마찬가지로.fa를 바꾸면.fai도 다시 만들어야 합니다. - 큰 염색체 색인 —
.csi가 필요할 때 — 참조 서열 하나가 512 Mbp를 넘으면 기본.bai색인이 못 담습니다. 이때는samtools index -c로.csi색인을 만듭니다. 사람 유전체는 대개.bai로 충분하지만, 일부 식물·양서류 유전체에서 마주칩니다. - 원격 BAM 구간 추출이 느리거나 실패 — samtools는 http·https·ftp·s3 URL을 직접 읽어
samtools view -b URL chr:start-end로 원격 BAM의 한 구간만 잘라올 수 있습니다. 다만 그 URL 옆에.bai색인이 함께 있어야 하고, 공개 FTP는 느리거나 경로가 바뀌기도 합니다. 재현이 중요하면 이 글처럼 작은 로컬 예제로 먼저 검증한 뒤 실데이터로 옮기는 편이 안전합니다.
확장 — 더 해 볼 것
markdup으로 PCR 중복 표시 — 라이브러리 증폭 과정에서 생긴 중복 리드는 변이 근거를 부풀립니다. samtools markdup으로 표시(또는 제거)하는데, 전처리 순서가 정해져 있습니다.
samtools sort -n -o namesort.bam input.bam # 1) 이름순 정렬
samtools fixmate -m namesort.bam fixmate.bam # 2) 짝 정보 태그 부착(-m 필수)
samtools sort -o positionsort.bam fixmate.bam # 3) 다시 좌표순 정렬
samtools markdup positionsort.bam markdup.bam # 4) 중복 표시 (-r이면 제거)- CRAM으로 저장 공간 줄이기 — BAM보다 더 압축되는 CRAM은 참조를 알아야 복원됩니다.
samtools view -T ref.fa -C in.bam -o out.cram으로 변환하면 대용량 저장에서 용량을 크게 아낄 수 있습니다. 참조만 잘 관리하면 무손실입니다. - 여러 샘플 합동 호출(joint calling) —
bcftools mpileup -f ref.fa a.bam b.bam c.bam | bcftools call -mv처럼 BAM 여러 개를 한 번에 넘기면, 샘플들을 함께 보고 변이를 호출합니다. 코호트 분석의 기본이죠. 더 정교한 GVCF 기반 합동 호출은 GATK 편의 방식과 비교됩니다. - VCF 병합·교집합 —
bcftools merge로 여러 샘플 VCF를 하나로 합치고,bcftools isec으로 두 콜셋의 교집합·차집합을 구할 수 있습니다. 서로 다른 호출기(bcftools vs GATK vs DeepVariant) 결과를 대조할 때 유용합니다. 변이 호출기별 차이는 변이 호출 개념 편에서 다뤘습니다. - 파이프라인으로 묶기 — 위 명령들을 셸 스크립트나 Snakemake·Nextflow 규칙으로 엮으면, FASTQ부터 필터된 VCF까지 재현 가능한 한 흐름으로 자동화할 수 있습니다. 각 단계 출력에 로그·체크섬을 남기면 디버깅이 쉬워집니다.
• (포맷 이해) SAM/BAM 파일이란? — 이 글에서 다룬 BAM의 11개 필드·FLAG·CIGAR를 개념부터
• (변이 포맷) VCF 파일이란? — 필터한 VCF의 8컬럼·INFO·GT를 완전 정리
• (직전 단계) BWA-MEM·STAR로 리드 정렬하기 — FASTQ를 정렬해 이 글의 입력인 SAM을 만드는 단계
• (더 무거운 호출기) GATK으로 변이 호출하기 — bcftools의 대안, BAM에서 VCF까지의 표준 파이프라인
References
- Li, H., Handsaker, B., Wysoker, A., et al. (2009). The Sequence Alignment/Map format and SAMtools. Bioinformatics, 25(16), 2078–2079. (samtools 원논문)
- Danecek, P., Bonfield, J. K., Liddle, J., et al. (2021). Twelve years of SAMtools and BCFtools. GigaScience, 10(2), giab008. (samtools·bcftools 통합 논문)
- Danecek, P., Auton, A., Abecasis, G., et al. (2011). The variant call format and VCFtools. Bioinformatics, 27(15), 2156–2158. (VCF 포맷 원논문)
- HTSlib project. (2024). samtools manual — view · sort · index · flagstat · coverage · faidx · markdup. htslib.org/doc/samtools
- HTSlib project. (2024). bcftools manual — mpileup · call · filter · norm · stats · query. htslib.org/doc/bcftools
- SAMtools/hts-specs. (2024). Sequence Alignment/Map (SAM) · Variant Call Format (VCF) 명세. hts-specs
Pipette & Pipeline · A bio portfolio journal
'Code Bench > NGS Pipelines' 카테고리의 다른 글
| Snakemake 재현 파이프라인 — rule·wildcard·DAG로 NGS 워크플로 자동화 (0) | 2026.07.10 |
|---|---|
| salmon으로 전사체 정량 — 유사정렬·RNA-seq 업스트림 (0) | 2026.07.08 |
| GATK으로 변이 호출하기 — BAM에서 VCF까지 (NGS Pipelines 3) (0) | 2026.06.24 |
| BWA-MEM·STAR로 리드 정렬하기 — FASTQ에서 BAM까지 (NGS Pipelines 2) (0) | 2026.06.24 |
| FastQC·fastp로 시퀀싱 리드 품질검사·트리밍 (NGS Pipelines 1) (0) | 2026.06.24 |
