유용한 VCF manipulation commands

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

[관련 포스팅 보기]

[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이다. 자세한 내용은 아래를 참고하자.

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

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/

[계속 Update 예정] 자주 쓰는 linux 명령어 및 프로그램 관련 자료

[Linux 자주 쓰는 명령어]

  • ls: 현재 위치 display
  • cd /folder : 이동
  • rm -rf /folder : 하위 폴더 포함 전체 삭제
  • mkdir /abc: abc 디렉토리 생성
  • chmod -R 777 ‘file1’: file 1 권한 변경 / 읽기: 4, 쓰기: 2, 실행: 1
  • cat, grep: cat /BiO/project/example.txt | grep “test”  내용 중 특정 문자열을 포함하는 줄만을 불러옴.
  • screen -ls : screen list 표시 / -S ‘screen1’ : screen1 생성 / -r ‘screen1’ : screen1 재입장 / Ctrl + A + D: screen 나가기
  • pkill screen : screen 모두 강제 종료
  • gzip ‘file1’ : file1 압축하기 / gzip -d ‘file1’ : file1 압축 풀기
  • tar 압축하기 : tar -cvf [파일명.tar] [폴더명]
  • tar 압축풀기: tar -xvf [파일명.tar]
  • tar.gz 압축하기 및 풀기: tar -zcvf / tar -zxvf : option에 z 추가
  • 여러개 gz 파일 합치기: cat file1.gz file2.gz file3.gz > allfiles.gz
  • cat /proc/cpuinfo | grep name: 서버 CPU 모델명 확인
  • nproc: 서버 코어 갯수 출력
  • df -h: 시스템 디스크 공간 확인
  • free -h: 현재 메모리 상황 확인
  • sar -r 10: 10초 간격으로 메모리 사용 현황 표시
  • rsync [file1] [path] -r –progress : 하위 폴더 포함 복사 & 복사 진행 상황 함께 보기
  • sudo: super user do – root 권한으로 실행
  • apt-get upgrade / apt-get update / apt-get install : 우분투 관련 패키지 업데이트 및 설치
  • 파일명 변경: mv or rename 이용
  • mv [파일명1] [파일명2]
  • rename [변경전문자] [변경후문자] [대상파일]: rename은 파일패턴문자(?, * 등)를 인식하기 때문에 여러 개의 파일에 동시 사용이 가능. Regex 정규 표현식을 이용해 사용 (링크 참조)

ex> file1~file9.png > file1~file9.jpg 확장자 변경
>rename png jpg file?.png

https://qgp9.github.io/blog/2016/03/24/rename-batch-rename-files

 

[Linux PATH 관련 설정]

  • echo $PATH: 현재 걸려있는 PATH 모두 확인
  • export PATH=”$PATH:/path/to/dir”: 경로 디렉토리를 일시적으로 추가

경로 디렉토리를 영구적으로 추가하기 위해서는, vi 에디터를 이용하여 ~/.profile에 위 경로를 추가하여 수정해주어야한다. (우분투 16의 경우, /etc/environment 에 등록하여도 됨.) 저장하고 나가기”:wq”
>sudo vi /etc/environment
export PATH="$PATH:/path/to/dir

이후에 해당 내용을 적용시키려면, source code를 실행시킨다.
>source /etc/environment

[java 관련 설정]

GATK 4를 돌리기 위해서는 JAVA 1.8 이상의 버젼이 필요한데, 최근에 오라클에서 유료로 전환되면서 업데이트가 힘들어졌다. 다행히 GATK는 무료로 배포되고 있는 Open JDK를 지원하기 때문에, Open JDK로 업데이트를 할 경우에 무리 없이 실행이 가능하다. 다음은 Ubuntu에서 Open JDK의 하나인 Zulu-8으로 업데이트를 위한 실행 명령어이다.

>sudo apt-key adv --keyserver hkp://keyserver.ubuntu.com:80 --recv-keys 0x219BD9C9
>sudo apt-add-repository 'deb http://repos.azulsystems.com/ubuntu stable main'
>apt-get update
>apt-get install zulu-8

위와 같이 zulu를 설치 후에, 기본 자바 홈을 vi 에디터를 이용하여, 설치 폴더로 변경해주어야 한다.

> vi /etc/profile 실행, 기존 자바 홈 설정 문구를 지우고, 아래 문구를 추가해준다. “:wq”
export JAVA_HOME=/usr/lib/jvm/zulu-8-amd64

>apt-get remove oracle*
>apt-get autoremove --purge
>apt-get autoclean

마지막으로 Java version을 확인하여, 설치가 제대로 되었는지 본다.
>> java -version
openjdk version "1.8.0_181"
OpenJDK Runtime Environment (Zulu 8.31.0.1-linux64) (build 1.8.0_181-b02)
OpenJDK 64-Bit Server VM (Zulu 8.31.0.1-linux64) (build 25.181-b02, mixed mode)

[Java option에 따른 의미]

>C:\java –X
–Xms : set initial Java heap size
–Xmx : set maximum Java heap size
–Xss : set java thread stack size

[conda, pip SSL error 관련]

> conda config --set ssl_verify False
> pip --trusted-host pypi.org --trusted-host files.pythonhosted.org install pip

회사 네트워크 보안 등의 설정으로 방화벽이 걸려있는 경우, SSL 인증서 문제가 발생하면서 업데이트가 안된다. 통신 에러 또는 방화벽 차단에 의해 conda update가 안되는 경우에는,  위의 명령어를 이용, ssl_verify 옵션을 꺼준다. pip의 경우는 –trusted-host 부분에 다운 받을 저장소의 주소를 등록하여 준다. 매번 이런 명령어를 쓰는 것이 귀찮으면, config 파일에 아래와 같이 등록하면 된다.

>sudo vi /etc/pip.conf
[global]
trusted-host = pypi.org
files.pythonhosted.org

[Samtools 이용 bam file header 변경]

for file in *.bam; do filename=`echo $file | cut -d "." -f 1`; samtools view -H $file | sed -e 's/SN:\([0-9XY]\)/SN:chr\1/' -e 's/SN:MT/SN:chrM/' | samtools reheader - $file > ${filename}_chr.bam; done

Last updated 2020-08-13