유용한 VCF manipulation commands

시퀀싱 결과 최종 파일은 변이들에 대한 정보를 담고 있는 VCF 파일입니다. 여러 가지 조건에 맞게 변이들을 Filtering하고, 입맛에 맞게 관련 정보들을 Annotation하는 것은 후속 단계에서 의미 있는 변이들만 남기는데 있어서 거의 필수적이라고 할 수 있습니다. 이번 포스팅은 VCF 파일을 조작하는데 있어 유용한 bcftools의 명령어들을 정리하려고 합니다.

[관련 포스팅 보기]

• NGS DNA-SEQ PIPELINE: GATK BEST PRACTICE CODE – PART3. VCF MANIPULATION
• [계속 UPDATE 예정] 자주 쓰는 LINUX 명령어 및 프로그램 관련 자료
• ANNOVAR: POPULATION FREQUENCY, IN-SILICO PREDICTION TOOL 및 기타 DATABASE 활용

[bcftools Manual 페이지 바로가기] https://samtools.github.io/bcftools/bcftools.html

VCF 파일 sorting 후 압축 및 index 파일 생성

vcf-sort data.vcf |bgzip -c >data.sorted.vcf.gz
tabix -p vcf data.sorted.vcf.gz

VCF 파일 내 chromosome 표시 방법 변경하기 (chr1 <-> 1)

awk '{gsub(/^chr/,""); print}' your.vcf > no_chr.vcf
awk '{if($0 !~ /^#/) print "chr"$0; else print $0}' no_chr.vcf > with_chr.vcf

awk, cut, sed 명령어는 매우 유용하게 쓰일 수 있는 linux command이다. 자세한 내용은 아래를 참고하자.

• AWK command in Unix/Linux with examples
• cut command in Linux/Unix with examples
• sed command in Linux/Unix with examples

Merging multiple VCF files into a single VCF file

bcftools merge *.vcf.gz -Oz -o merged.vcf.gz # merge multiple samples into one
bcftools concat *.vcf.gz -Oz -o merged.vcf.gz # merge all chromosomes into one

• How to merge 7000 VCF files with bcftools merge?
• 샘플별로 쪼개져있는 vcf 파일을 하나의 파일로 합칠 때: bcftools merge
• Chromosome 단위로 쪼개져있는 vcf 파일을 하나의 파일로 합칠 때: bcftools concat

Splitting VCF into individual samples or chromosomes

# split to samples
for sample in `bcftools query -l multisample.vcf.gz`; do
   bcftools view -c1 -Oz -s $sample -o $sample.vcf.gz multisample.vcf.gz
done

# split to chromosomes
for i in {1..22}; do
   bcftools view -r ${i} file.vcf.gz -o file.chr${i}.vcf.gz
done

VCF 파일 변이 갯수 확인하기

bcftools view -H data.vcf.gz | wc -l

Allele frequency로 filtering 하기

bcftools view -q 0.01 data.vcf.gz -Oz -o filtered.common.vcf.gz # keep only common variants
bcftools view -Q 0.01 data.vcf.gz -Oz -o filtered.rare.vcf.gz #keep only rare variants
bcftools view -c 100 data.vcf.gz -Oz -o filtered.common.vcf.gz # keep AC >= 100
bcftools view -C 100 data.vcf.gz -Oz -o filtered.rare.vcf.gz #keep AC <= 100

-q, –min-af minimum allele frequency (INFO/AC / INFO/AN) of sites to be printed.

-Q, –max-af maximum allele frequency (INFO/AC / INFO/AN) of sites to be printed.

-c, –min-ac minimum allele count (INFO/AC) of sites to be printed.

-C, –max-ac maximum allele count (INFO/AC) of sites to be printed.

BED 파일을 이용하여 특정 Target 영역 추출하기

vcftools --vcf input.vcf --bed target.bed --out output_prefix --recode --keep-INFO-all

[참고 사이트]

https://plink.readthedocs.io/en/latest/