VCF와 BED, 명령어로 어떻게 다룰까요? 변이·구간 파일 실전

📌 이 글에서 알 수 있는 것

• VCF 한 줄(CHROM~FORMAT, GT·DP·GQ·AD) 읽는 법
• bcftools view·query·norm·filter·isec·stats의 옵션별 뜻
• BED와 VCF 좌표가 하나씩 어긋나는 이유와 변환법
• bedtools sort·merge·intersect·coverage, IGV로 변이 눈으로 확인하기

🧬 VCF 한 줄 읽기

VCF(Variant Call Format)는 변이 호출 결과를 담는 탭 구분 텍스트입니다. ##로 시작하는 메타 정보 줄(첫 줄은 ##fileformat), #CHROM으로 시작하는 헤더 줄, 그리고 변이마다 한 줄이 이어집니다¹. 예시는 변이 호출기가 만든 sample.vcf, ref.fa, targets.bed·genes.bed, 7편의 sample.markdup.bam을 가정합니다(호출 원리는 3편).

#CHROM POS  ID REF ALT  QUAL FILTER INFO  FORMAT      sample1
chr1   1000 .  A   G    50   PASS   DP=32 GT:DP:GQ:AD 0/1:32:99:15,17
chr1   2000 .  C   T,G  40   PASS   DP=28 GT:DP:GQ:AD 1/2:28:60:2,12,14

• POS: 1부터 셉니다. REF·ALT: 참조·대체 염기. indel은 앞 염기 하나를 붙여 씁니다.
• QUAL: ALT 판정이 틀렸을 확률의 Phred 값. FILTER: 모든 필터 통과면 PASS, 실패하면 필터 이름, 필터를 안 돌렸으면 .입니다.
• GT: 0은 REF, 1은 첫째 ALT, 2는 둘째 ALT. /는 phase 없음, |는 phase 있음, ./.는 결측. DP는 리드 깊이, GQ는 유전자형 품질, AD는 대립유전자별 리드 수(REF부터)입니다¹.
• 둘째 줄처럼 ALT가 여러 개면 multi-allelic 레코드입니다.

🗜️ 압축과 색인

bgzip -c sample.vcf > sample.vcf.gz
tabix -p vcf sample.vcf.gz
bcftools index -s sample.vcf.gz

• bgzip -c: 작은 블록(BGZF) 단위로 압축해 표준출력으로 보내고 원본은 남깁니다. gzip과 호환됩니다².
• tabix -p vcf: .tbi 색인을 만듭니다. 입력은 위치순 정렬과 bgzip 압축이 필요합니다. TBI는 512 Mbp까지 되고, 더 긴 염색체는 CSI가 필요합니다³.
• bcftools index: 기본은 CSI(-t면 TBI). -s는 색인에서 염색체별 레코드 수를 바로 읽습니다⁴.
• 자주 하는 실수: gzip으로 압축해 놓고 색인이 안 된다고 헤맵니다.

이하 bcftools 옵션은 1.25 매뉴얼 기준입니다⁴. 여러 VCF를 동시에 읽는 명령은 모두 색인이 필요합니다.

🔎 view·query: 거르고 표로 만들기

bcftools view -r chr1:1000000-2000000 -s sample1,sample2 \
  -f PASS -i 'QUAL>=30 && INFO/DP>=10' -Oz -o sub.vcf.gz -W sample.vcf.gz

• -r: 영역(1-based). POS가 영역 밖이라도 겹치는 indel은 포함됩니다. 여러 영역은 -R regions.bed로 주는데, 확장자가 .bed면 0-based BED로 읽습니다.
• -s: 여러 시료 VCF에서 시료 선택, ^sample3처럼 ^를 붙이면 제외. AC=0 줄도 남으니 -c 1을 더하세요.
• -f PASS: FILTER가 PASS인 줄만. 필터를 안 돌린 .까지 남기려면 -f .,PASS.
• -i·-e: 포함·제외 조건식. INFO/DP, FMT/DP, GT=“het”, TYPE=“snp” 등을 씁니다. 태그 이름은 대소문자를 구분합니다.
• -Oz -o: 압축 VCF로 저장(파이프 사이는 -Ou가 빠름). -W는 출력 색인을 함께 만듭니다.
• 자주 하는 실수: 조건식을 작은따옴표로 감싸지 않아 셸이 &&·>를 가로챕니다.

bcftools query -H -f '%CHROM\t%POS\t%REF\t%ALT[\t%SAMPLE=%GT:%DP]\n' sub.vcf.gz > table.tsv

• -f: 출력 형식. %열이름, %INFO/태그, [ ] 안은 시료마다 반복됩니다. -H는 열 이름 줄을 붙입니다.
• 자주 하는 실수: FORMAT 태그(%GT 등)를 [ ] 밖에 써서 오류가 납니다.

📐 norm: 비교 전에 정규화하기

같은 변이도 VCF에 여러 방식으로 적힐 수 있습니다. TCACACAG에서 CA 하나가 빠진 결실은 POS 100에 TCA→T로도, POS 104에 ACA→A로도 적을 수 있고, 결과 서열은 같습니다. 표준 정규화는 multi-allelic 분리, 불필요한 염기 다듬기, 왼쪽 이동입니다. 다만 인접한 복합 변이는 정규화 뒤에도 표현이 하나로 정해지지 않습니다⁵.

bcftools norm -m -both -f ref.fa -Oz -o sample.norm.vcf.gz -W sample.vcf.gz

• -m -both: multi-allelic을 biallelic 여러 줄로 나눕니다. 1/2는 1/0·0/1이 되어 het처럼 보이니 접합성 판정에 주의하세요.
• -f ref.fa: 참조 FASTA를 줘야 indel 왼쪽 정렬이 됩니다. -c는 REF가 참조와 다를 때 동작(기본 e 중단, w 경고, x 제외)을 정합니다.
• 자주 하는 실수: 호출에 쓴 것과 다른 판본의 FASTA를 주는 것입니다.

64개 phased 콜셋을 다시 채점한 2023년 연구는, 기존 평가가 특정 표현 방식에 유리하게 치우쳐 파이프라인 순위를 잘못 보여 줄 수 있다고 보고했습니다⁵. 임상 보고에 쓰는 HGVS(변이 이름 표기 규칙)는 VCF와 달리 전사체 3′ 쪽(+ 가닥 유전자면 오른쪽)으로 옮겨 적으니, 두 표기를 맞출 때 주의하세요⁵.

🚦 filter·isec·stats

bcftools filter -s LowQual -e 'QUAL<30 || INFO/DP<10' -Oz -o sample.flt.vcf.gz sample.norm.vcf.gz

• -e: 조건을 만족하는 줄을 걸러냅니다. -s LowQual을 주면 지우지 않고 FILTER 열에 LowQual을 표시합니다(soft filter). 기본은 기존 FILTER 값을 덮어쓰니, 이전 표시도 남기려면 -m +를 줍니다.
• 자주 하는 실수: 처음부터 지워 버려, 나중에 기준을 바꾸거나 걸러진 이유를 볼 수 없게 됩니다.

bcftools isec -p isec_out callerA.norm.vcf.gz callerB.norm.vcf.gz

• -p: 결과 폴더. 0000은 A에만, 0001은 B에만, 0002·0003은 공통(각각 A·B 기준) 레코드이고, README.txt에 설명이 남습니다⁶.
• 기본(-c none)은 REF·ALT가 같아야 같은 변이로 봅니다. -c all은 위치만 봅니다. -n=2 -w1은 공통 레코드를 A 기준으로 하나만 씁니다.
• 입력은 두 호출기(A·B) 결과를 각각 norm 한 파일입니다. 자주 하는 실수: 정규화하지 않은 VCF를 넣는 것입니다.

bcftools stats -f PASS -s - sample.flt.vcf.gz > sample.vchk
plot-vcfstats -p stats_plots sample.vchk

• stats: 변이 종류·Ts/Tv(전이/전환 비)·QUAL·깊이 분포 등. -f PASS는 LowQual 등 필터 표시가 붙은 줄을 빼고, -s -는 시료별 통계를 더합니다.
• plot-vcfstats -p: 그림 폴더. matplotlib이 필요하고, pdflatex가 없으면 -P로 PDF를 생략합니다.

📏 BED: 0부터 세는 반열림 구간

BED의 start는 0부터 세고 end는 포함하지 않습니다(반열림). start=9, end=20은 10~20번 염기입니다⁷.

bcftools query -i 'FILTER="PASS"' -f '%CHROM\t%POS0\t%END\t%REF>%ALT\n' sample.flt.vcf.gz > variants.bed

• %POS0: 0-based 시작, %END: REF가 끝나는 위치. 그대로 BED가 됩니다⁴. indel은 앞 염기(padding)까지 구간에 들어가니 엑손 경계 판정에 주의하세요.
• 자주 하는 실수: awk로 POS를 그대로 start·end에 넣어 구간이 한 칸 밀립니다. 또 -i 'FILTER="PASS"'를 빼면 LowQual 줄까지 BED에 들어갑니다.

🧰 bedtools: 구간 연산

이하 옵션은 bedtools 문서(2.31.0) 기준입니다⁸.

bedtools sort -i targets.bed > targets.sorted.bed
bedtools merge -i targets.sorted.bed -d 0 > targets.merged.bed

• sort: 염색체, 시작 위치 순으로 정렬. -faidx ref.fa.fai면 참조 순서를 따릅니다.
• merge: 겹치거나 맞닿은 구간을 합칩니다. -d는 합칠 최대 간격(기본 0). 정렬된 입력이 필요합니다.
• 자주 하는 실수: 정렬하지 않은 BED를 merge에 넣는 것입니다.

bedtools intersect -a variants.bed -b genes.bed -wa -wb > var_gene.tsv
bedtools intersect -a variants.bed -b targets.merged.bed -u > on_target.bed
bedtools intersect -a variants.bed -b targets.merged.bed -v > off_target.bed

• -a·-b: A의 구간마다 B와의 겹침을 찾습니다. BED·VCF·BAM·GFF를 받습니다.
• -wa -wb: A와 겹친 B 줄을 나란히 출력(어느 유전자에 걸렸나). -u: 겹침이 있으면 A를 한 번만. -v: 겹치지 않는 A만.
• -f 0.5처럼 최소 겹침 비율을 줄 수 있고(기본 1 bp), 큰 파일은 -sorted가 메모리를 아낍니다.
• 자주 하는 실수: -u 없이 실행해 여러 B와 겹친 변이가 중복으로 세어집니다.

bedtools coverage -a targets.merged.bed -b sample.markdup.bam > target_cov.tsv
bedtools getfasta -fi ref.fa -bed targets.merged.bed -fo targets.fa

• coverage: A 구간마다 겹친 리드 수, 덮인 염기 수, 구간 길이, 덮인 비율을 붙입니다. -hist는 깊이 분포. 중복·MAPQ 0 리드도 세니, 깊이 판정은 7편의 samtools coverage(중복 제외)와 맞춰 보세요.
• getfasta: 구간 서열을 뽑습니다. FASTA 헤더와 BED 염색체 이름이 정확히 같아야 합니다.
• 자주 하는 실수: 2.24.0 이전 문서대로 -a·-b를 뒤바꿔 넣는 것입니다.

👀 IGV로 눈으로 확인하기

IGV(Integrative Genomics Viewer)는 정렬과 변이를 보는 무료 뷰어입니다. 데스크톱은 File > Load from File로 BAM과 VCF를 열고, BAM은 같은 폴더에 .bai 색인이 있어야 합니다. 참조와 다른 염기는 색으로, 품질이 낮으면 흐리게 보입니다. 기본 30 kb보다 넓게 보면 리드가 숨겨지니 변이 위치로 확대하세요⁹. VCF 트랙에서 회색은 참조와 같음, 진한 파랑은 heterozygous, 청록은 homozygous variant입니다¹⁰.

볼 것은 대체 염기가 양쪽 가닥 리드에 고르게 있는지, 리드 끝에만 몰렸는지, MAPQ가 낮은 리드가 많은 반복 영역인지입니다. IGV의 genome 판본이 VCF와 같아야 합니다.

⚠️ 한계와 주의

• chr 표기: bcftools는 “chr20"과 “20"을 다른 이름으로 봅니다⁴. 결과가 조용히 비면 bcftools annotate --rename-chrs map.txt로 맞추세요.
• 참조 판본: GRCh37과 GRCh38은 좌표가 다릅니다. 다시 매핑할 수 없으면 BCFtools/liftover처럼 indel·multi-allelic을 지원하는 도구를 검토합니다¹¹.
• 정렬·색인: merge·-sorted·tabix는 정렬된 입력이, -r·isec은 색인이 필요합니다.
• 정규화 없는 비교: 같은 변이가 두 개로 세어집니다. 정규화한 뒤에도 복합 변이는 hap.py 같은 서열 기반 비교(3편)가 낫습니다⁵.

🏛️ 지금 어디까지 왔나

• 도구 판본(2026-11-06): bcftools·htslib 1.25(2026-10-27 배포), bedtools 2.31.1(2023-11, 문서는 2.31.0), IGV 데스크톱 2.19.x입니다. samtools·BCFtools는 Bioconda로 100만 회 넘게 설치됐습니다(2021)¹².
• 형식 표준: VCF 명세는 hts-specs에서 관리되며, 이 글은 2026-02-25 판 VCFv4.5를 확인했습니다¹.

🧪 연구자 노트

• 좌표 기준(0부터인지 1부터인지)이나 참조 유전체 판본(GRCh37·GRCh38)이 맞지 않으면, 오류 없이 엉뚱한 결과가 나옵니다. 이 두 가지를 먼저 확인하세요.
• 변이 필터는 PASS 표시만 믿지 말고, 깊이와 품질 분포를 직접 본 뒤 기준을 정하세요.

🙋 그래서 나에게는? — 실습 전 체크리스트

• bcftools --version, bedtools --version을 기록했다
• VCF는 bgzip 압축 후 색인했다
• VCF·BED·BAM·FASTA의 염색체 이름(chr1/1)과 판본이 같다
• 비교 전 norm -m -both -f ref.fa로 정규화했다
• 필터는 soft filter로 남기고 기준을 기록했다
• VCF→BED 변환에 %POS0을 썼다
• BED는 정렬 후 merge·intersect했다
• 핵심 변이는 IGV에서 리드로 확인했다

⚠️ 이 글은 의학적 조언이 아니며, 특정 제품을 권하지 않습니다. 실제 옵션은 사용하는 도구 버전의 매뉴얼과 기관 SOP, 자체 검증 결과를 우선하세요.

❓ 자주 묻는 질문

Q. bcftools index와 tabix 중 무엇을 써야 하나요? A. 둘 다 됩니다. bcftools index는 기본 CSI, tabix는 TBI를 만들고, 512 Mbp보다 긴 염색체가 있으면 CSI를 쓰세요.

Q. isec에서 공통 변이가 너무 적게 나옵니다. A. 먼저 두 파일을 같은 참조로 norm 했는지, chr 표기가 같은지 보세요. 그래도 차이가 크면 복합 변이 표현 차이일 수 있어 서열 기반 비교 도구로 다시 확인합니다.

Q. BED를 엑셀로 열어 고쳐도 되나요? A. 피하세요. 좌표가 숫자 서식으로 바뀌거나 탭이 깨질 수 있습니다. 고친 뒤에는 bedtools sort로 다시 정렬하세요.

📚 오믹스 기술 시리즈

• 1편 WGS 플랫폼 비교
• 2편 WGS 커버리지·깊이 설계
• 3편 WGS 데이터 처리와 QC
• 4편 WGS vs WES vs 패널
• 5편 리눅스 명령줄 기본기
• 6편 분석 환경 만들기
• 7편 FASTQ·BAM 실전
• 8편 VCF·BED 실전 (이 글)
• 9편 워크플로 관리
• 10편 R과 Python
• 11편 RNA-seq 실험 설계
• 12편 RNA-seq 데이터 분석
• 13편 small RNA-seq 실험 설계
• 14편 small RNA-seq 데이터 분석

참고