[카테고리:] bioinformatics

구조 변이 annotation tool: AnnotSV

유전체 정보로 부터 임상적으로 중요한 변이를 검출하기 위해서는 NGS 시퀀싱 기기의 read 정보로부터 변이 검출까지의 파이프 라인 못지 않게, 얻어낸 수많은 변이로 부터 병인에 중요한 역할을 할 것으로 생각되는 후보 변이를 필터링하는 전략이 매우 중요하게 됩니다. 따라서 적절한 데이터 베이스로 부터 Annotation을 하는 과정은 매우 중요하게 되는데, 이번 포스팅은 다양한 변이 중에서 구조 변이 (Structural variation; SV)을 대상으로 Annotation을 할 수 있는 도구 중에 하나인 AnnotSV에 대해서 소개하고자 합니다. AnnotSV는 이전에 소개했던 Annovar의 CNV (copy number variant) 버젼에 해당한다고 볼 수 있습니다.

[관련 포스팅 보기]

AnnotSV는 다양한 구조 변이의 Annotation 기능 뿐만 아니라, ACMG (American College of Medical Genetics)에서 권장하는 구조 변이의 판독 기준에 따라서, 해당 변이의 중요도를 5가지 카테고리로 구분해줍니다. Input으로는 bed 파일 또는 vcf 파일을 받으며, 다양한 유전자, 조절 인자, 기존에 알려진 병적 변이, 질병과의 연관성 등을 기준으로 ACMG class를 보고해줍니다. 위 그림은 AnnotSV의 이러한 분석 과정을 보여주고 있습니다.

[bed 파일의 기본 구조] bed 파일은 1) 염색체 번호 (Chromosome), 2) 시작 지점 (Start), 3) 끝 지점 (End)의 3가지 기본적인 정보를 토대로 유전체 내의 특정 범위에 대한 정보를 제공해줍니다.

구조 변이 (CNV)의 경우, 적은 수의 염기 변이 (SNV)보다 short-read sequencing을 이용하는 경우, 기술적으로 검출하는 해상도의 한계가 있으며 (deletion보다 duplication 검출이 어려움. 충분한 Depth와 Supporting read가 확보되어야 하며, 이 때문에 translocation도 검출이 어려움.) 변이의 해석도 더 어려운 경우가 많습니다. 따라서, 적절한 한계점을 인지하고 적절한 분석 방법론을 적용하는 것이 중요하며, 현재도 많은 부분들이 현재 진행형으로 연구가 되고 있는 분야입니다.

다만, 최근 ACMG에서 구조 변이의 임상적 해석을 위한 Criteria를 제시해주어, 많은 부분 임상적으로 활용이 가능해진 부분이 있습니다. (아래 참고 논문: Riggs, Erin Rooney, et al. Genetics in Medicine 22.2 (2020): 245-257) 그동안 구조 변이의 해석에 여러가지 어려운 점들이 많았는데, 최근 이 쪽 분야도 많은 툴들과 방법론 들이 개발되고 있는 것 같습니다. 그런 점에서 AnnotSV는 구조 변이를 연구하고 해석하는 입장에서 매우 유용한 툴임이 분명합니다.

[References]

AnnotSV Github: https://github.com/lgmgeo/AnnotSV

AnnotSV Homepage: https://www.lbgi.fr/AnnotSV/

Geoffroy, Véronique, et al. “AnnotSV: an integrated tool for structural variations annotation.” Bioinformatics 34.20 (2018): 3572-3574.

Geoffroy, Véronique, et al. “AnnotSV and knotAnnotSV: a web server for human structural variations annotations, ranking and analysis.” Nucleic Acids Research (2021).

Riggs, Erin Rooney, et al. “Technical standards for the interpretation and reporting of constitutional copy-number variants: a joint consensus recommendation of the American College of Medical Genetics and Genomics (ACMG) and the Clinical Genome Resource (ClinGen).” Genetics in Medicine 22.2 (2020): 245-257.

AlphaFold를 이용한 단백질 구조 예측

최근 Alphafold의 소스 코드와 전체 단백질 예측 구조 데이터 베이스가 공개되면서, 정말 다양한 방면에 활용이 가능할 것으로 생각되고 있습니다. 유전체 연구자의 입장에서는 Alphafold를 이용하면 수 많은 변이들의 in-silico structural prediction이 가능하기에, 과거에 존재하던 variant prediction tool들 보다 훨씬 강력하고 정확한 예측 도구가 나타날 것으로 생각됩니다. 이번 포스팅은 공개된 Alphafold2의 소스 코드와 데이터 베이스 자료를 정리해보고자 합니다.

[Alphafold2 모델 모식도]

Alphafold2의 경우, 기존에 실험적으로 17% 정도만 알려져있던 단백질의 3차원적인 구조를 딥러닝 기법을 활용하여, 존재하는 거의 모든 단백질에 대해서 높은 정확도로 예측하여 해당 구조를 데이터베이스로 제공하고 있습니다. 이는 과거에 알지 못했던 많은 단백질의 3차원 구조들에 대한 높은 예측을 제공하기 때문에, 완전히 새로운 필드를 열었다고 해도 과언이 아닙니다. 특히나 단백질의 3차원 모델화가 가능하다는 것은 이 모델을 활용하여 추후에 파생되는 약물 상호 작용이나 변이의 효과에 대한 예측도 가능하다는 이야기가 되기 때문에, 그 파급력은 더욱 클 것으로 생각됩니다.

현재 AlphaFold Protein Structure Database는 약 2만개의 단백질에 대한 정상 구조를 예측되어 제공하고 있지만, 개인적으로는 인간의 유전체에 존재하는 수 많은 변이 (Missense & Indel variants) 에 대해, mutant protein 들의 예측도 포함하는 Database가 update가 되길 기대해 봅니다. 현재까지의 아미노산 서열의 변화에 대한 많은 예측 Tool이 개발되었지만, 임상적으로 해석하고 활용하기에는 많은 한계가 있었는데, 최종 변이 단백질의 구조를 높은 정확도로 예측할 수 있다면, 그 해석과 활용도 용이해지기 때문입니다. 단백질 하나에 존재할 수 있는 아미노산 변이의 종류는 수만가지가 넘을 것이기 때문에 (아미노산 20종 x 단백질 서열 크기 500~2,000개 아미노산), 이러한 Database를 구축하기 위해서 컴퓨터가 연산해야할 데이터량이 무지막지하겠지만, 결국에는 이러한 in-silico saturated mutagenesis proteome Database도 추후에 구축되지 않을까 예상합니다. (비트 코인 채굴하는 자원이면, 금방 될 듯 합니다…)

[관련 포스팅 보기]

[References]

Deepmind Blog 소개 글: https://deepmind.com/blog/article/putting-the-power-of-alphafold-into-the-worlds-hands

Open source github 자료: https://github.com/deepmind/alphafold/

AlphaFold Protein Structure Database: https://alphafold.ebi.ac.uk/

Jumper, John, et al. “Highly accurate protein structure prediction with AlphaFold.” Nature (2021): 1-11.

Jumper, J., Evans, R., Pritzel, A. et al. “Highly accurate protein structure prediction with AlphaFold”. Nature (2021).

유전체 데이터 분석을 위한 서버 운영 (feat. harvard e2)

저는 현재 하버드 부속 연구기관 중 하나인 BCH (Boston Children’s Hospital)의 연구 서버 (Researching computing) “E2“를 이용하여, WGS 분석을 진행하고 있습니다. E2는 대략 1000개 이상의 연산용 cpu를 가지고, BCH 소속 많은 연구실의 연구자들의 Bioinformatics 관련 분석을 처리해주는 서버입니다. Genome sequencing 관련 분석은 Target 또는 Exome sequencing과 다르게 그 엄청난 데이터 사이즈로 인하여, 연구용 서버를 이용하여 많은 수의 샘플을 처리하는데 많은 Computational Resource (=cpu&memory)를 요구하게 됩니다. 이번 포스팅에서는 추후 효율적인 연구용 서버의 활용에 참고하기 위해, 하버드 부속 연구 기관에서 다양한 Omics 관련 분석을 위한 연구용 서버를 어떻게 운영하는지에 대해서, 기록으로 남겨두고자 합니다.

[e2 서버의 대략적인 구성과 스펙]

연산 서버와 데이터 저장 서버

기본적으로 연산 서버와 데이터 저장 서버가 나뉘어져 있고, 연산용 CPU 서버는 중앙 관리자가 관리를 맡아서 연구자들이 공유하여 사용할 수 있도록 돕고 있다. 연구를 지원하는 연구 지원 파트에 서버 관리만 전담으로 하는 인력들이 존재한다. (이들의 역할이 서버 유지에 매우 중요함.) 반면, 데이터 저장 (Storage) 공간은 연구실별로 사용료를 지불하고 저장 공간을 구매해서 사용한다. 연산용 서버는 여러 연구실에서 기기를 기증을 받은 후에, 이들을 모두 묶어서 하나의 커다란 Cluster를 구축하였다. 기증자에게는 High Priority (우선권)를 주어 해당 서버에 작업을 던지는 경우, 다른 일반 사용자보다 먼저 작업이 돌아가도록 구축되어 있다. 이렇게 모든 Computational resource를 하나의 Cluster로 구축하면, 더 빠르고 집중적으로 자원을 배분하면서 사용할 수 있기 때문에, 효율적으로 서버를 사용하는데 유리하다. 서버의 연산 작업은 24시간 내내 돌아가는 것이 아니기 때문에, 컴퓨팅 파워를 집약적으로 사용하면서, 동시에 기증자에게는 Advantage를 주어, 운영에 적합한 정책으로 생각된다.

Software Package

기본적으로 서버에 설치되어 있는 프로그램들은 연구를 진행하면서, 흔하게 사용되는 tool들을 위주로 Biogrid라는 Package로 구성되어 있다. 이는 Harvard에서 BioGrids consortium을 구축하여 관리하고 있는데, 대략 300 여개의 Genomics and Bioinformatics toolsets을 포함하고 있다. 해당 패키지는 매일 업데이트 되고, 서버에서는 해당 패키지를 로딩하면 바로 쓸 수 있게 구성되어 있어서, 프로그램을 따로 설치할 필요가 없다. 역시나 BioGrids consortium에는 연구용 프로그램 패키지만을 지원하고 관리하는 인력이 따로 존재한다!

[Dependency 관련 문제] Biogrids 패키지에 속해있는 프로그램들의 경우는 따로 dependency를 고려하지 않아도 되나, 추가적인 프로그램들을 구동하는데 필요한 환경 (Dependency)을 구축할 경우, Docker 또는 Singularity와 같은 Container를 이용하여 Handling을 하도록 되어있다.

[관련 페이지] BioGrids Consortium : https://biogrids.org/

[관련 페이지] Singularity Documentation : https://sylabs.io/guides/2.6/user-guide/index.html

Job Scheduler: Slurm

연산 작업의 매니저로는 Slurm을 이용하여, 관리하고 있다. 모든 사용자들은 Slurm을 통해서만 작업을 제출할 수 있으며, 사용자 정책 (User Policy)을 따라서 사용하도록 교육을 받는다. 사용자에게 할당된 기본적인 디스크는 아래와 같다.

[관련 페이지] Slurm 관련 문서 보기 : https://slurm.schedmd.com/documentation.html

On Demand

추가적으로 E2는 Web-page를 통해 GUI 기반으로 서버에 접근할 수 있는 (Interactive Session) On Demand 서비스를 구축하고 있다. 특히, Jupyter-Notebook, Rstudio, MATLAB, IGV Browser와 같이 범용적이고 흔하게 쓰이는 프로그램의 경우에는, 사용자들이 바로 개인 노트북 등을 통해 서버에 접근하고 작업을 제출하는 용도로 널리 쓰이고 있다.

[NGS DNA-SEQ] Functional Equivalence pipeline: CROMWELL, WDL

gnomAD, TOPMed 등 대규모 유전체 코호트들이 만들어지면서, 여기서 생산된 데이터를 이용하는데 중요한 문제가 부각되었는데, 바로 분석 결과 간의 재현성과 호환성에 있었습니다. 즉, 연구자가 GATK Best practice를 이용하여 Exome 또는 Genome 시퀀싱 분석을 진행하더라도 어떤 설정과 파라미터 값을 넣느냐에 따라, 최종 검출 변이의 결과가 달라지고, 이것은 연구 결과 간의 재현성의 측면에서 매우 중요한 문제가 되었던 것이지요. 참고논문에서 진행한 실험 결과를 보면, 동일한 샘플로 생산된 FASTQ 파일을 서로 다른 5개의 기관에 보내 각각의 파이프 라인으로 분석한 결과, Call된 변이들 간에 많은 차이가 있었다고 보고하고 있습니다.

[관련 포스팅 보기]

따라서, 점점 늘어나는 유전체 데이터만큼 유전체 분석 파이프 라인을 하나의 표준화된 파이프 라인으로 통합하는 것이 매우 중요해졌고, 그 결과 개발된 것이 “Functional Equivalence (FE)” Pipeline (기능적으로 동등한 파이프라인) 입니다. 사실 NGS 분석을 하는 사용자의 입장에서는 Input만 넣고, Output만 나오면 편한데, 그동안 개발된 툴들은 이를 모두 아우르는 것이 아니라, 그때 그때마다 필요한 부분들을 개발했기 때문에, 분석 파이프 라인도 이제야 어느 정도 성숙 단계에 이르렀다고 할 수 있습니다. 따라서 최근의 대규모 유전체 컨소시엄들은 모두 “Functional Equivalence ” Pipeline 을 통해 생산된 유전체 데이터를 생산하는 것으로 채택하고 있습니다. (그래서 저도 functional equivalent한 결과를 얻기 위해서 최근에 새롭게 공부를 하게 되었습니다.)

[ Functional Equivalence Pipeline Overview]

FE 파이프라인을 제공하기 위해서, Broad Institute의 개발진들은 WDL (Workflow Description Langauge)과 Cromwell이라고 하는 프로그래밍 언어를 개발하는데, 하나의 파이프라인을 패키지로 묶은 WDL 파일을 만들고, 이를 Cromwell이라는 프로그램으로 구동시켜주는 원리라고 합니다. 사실 사용자의 입장에서는 과거에 개별 프로그램을 설치하고, 개별 스텝을 따로 돌려야했다면, 지금은 이러한 것들이 모두 하나의 패키지 형태로 제공되기에 더욱 편해졌다고 할 수 있습니다 (?).

[Cromwell 페이지 바로 가기] https://cromwell.readthedocs.io/en/stable/

[WARP 페이지 바로가기] https://broadinstitute.github.io/warp/docs/get-started/

Cromwell의 로고: 꼬마돼지 베이브와 스타트렉에 출연한 배우 James Cromwell을 오마주한 로고라고 합니다 🙂

병원 검사실을 운영하는데, 검사 장비와 보고 방법을 표준화하는 것은 매우 중요한데, 유전체 분석 파이프 라인에도 이제야 이러한 개념이 들어왔다는 점이 반갑습니다. (분석 파이프 라인 하나도 이렇게 표준화하기가 어렵습니다.) Genome의 경우에는 처리해야할 데이터의 크기가 워낙 방대하기때문에 더욱 어려운 점이 있는 것 같습니다. 현재 이쪽 분야도 많은 Computational Scienctist들이 뛰어들어서 개발을 진행하고 있는 중이기 때문에, 추후에 더 User-friendly하고 간편한 파이프 라인이 개발되어 제공되기를 기대해 봅니다. (점차 대세는 클라우드로 옮겨가지 않을까 합니다?)

[References]

Regier, Allison A., et al. “Functional equivalence of genome sequencing analysis pipelines enables harmonized variant calling across human genetics projects.” Nature communications 9.1 (2018): 1-8.

유용한 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/

Protein Structure Visualization: RCSB Protein Data Bank (PDB) & PyMOL

단백질을 구성하는 아미노산의 서열을 안다고 하더라도 단백질의 3차원적인 구조 (Protein 3D structure)를 예측하는 것은 매우 어려운 일입니다. 최근 구글 딥마인드에서 딥러닝을 이용한 AlphaFold를 개발해서, 3차원적인 구조를 성공적으로 예측한 논문 (아래 reference) 이 발표된 바 있지만, 전통적으로 가장 확실하게 구조를 확인하는 방법은 X-ray 회절을 이용한 Crystallography 방법 입니다.

많은 연구를 통해서, 이러한 3차원 단백 구조를 보관하고 있는 data bank가 있는데, 오늘 이용할 RCSB Protein Data Bank (PDB)가 그 곳입니다.  최근 논문을 쓰다보니, 단백질의 치환 변이를 시각화한 Figure를 만드려다 보니, PyMOL을 이용하게 되어, 해당 과정을 블로그 포스팅으로 남깁니다.

 

I. PDB를 통한 Protein Struture 정보 확인

PDB

RCSB PDB 바로가기 > https://www.rcsb.org/

 가장 먼저 할 일은 위의 database에 해당하는 단백의 PDB 파일을 받거나, PDB ID를 확인하는 것입니다. 사람의 경우, Homo Sapiens인지를 확인하고, 관심있는 유전자 또는 단백의 이름을 입력한 후에 검색을 하면, 여러개의 결과가 display됩니다. 단백질 결정의 구조를 나타낸 것이기 때문에, 어떠한 Ligand를 가지고 실험을 하였는가에 따라서, 구조가 조금씩 달라지게 됩니다. 가장 적절한 것을 선택하여, 구조를 확인하거나, PDB 파일을 다운로드 받아서 활용하게 됩니다.

 

II. PyMOL 설치 및 PDB 파일 불러오기

위의 PDB에서 받은 .pdb 파일은 다양한 툴로 분석 및 시각화가 가능한데, 가장 널리 쓰이는 프로그램은 파이썬 기반으로 개발된 PyMOL입니다. 최근에 무료에서 유료로 전환되면서, 접근성이 떨어지긴 하였지만, 설치 후 바로 30일간 무료 체험판으로 이용 가능합니다. 윈도우에서 바로 이용 가능하며, anaconda 등을 이용하여 Python 기반으로도 이용 가능합니다.

PyMOL 다운로드 바로가기 > https://pymol.org/

 

III. PyMOL을 이용한 시각화

위의 링크를 통해, PyMOL을 설치하면, 프로그램에서 PDB 파일을 불러올 수 있습니다. 여러가지 명령어와 패키지 등을 활용하면, 단백질 구조에 대한 고급 분석이 가능하다고 하는데, 저는 가장 기본적인 기능만을 활용해서, 유전체 분석을 통해 확인된 변이들의 위치를 시각화하는데 이용 해보도록 하겠습니다. (논문용 Figure…) 더 자세한 세부 내용은 링크에서 찾으시길 바랍니다.

초보자를 위한 PyMOL 시작하기 > https://pymolwiki.org/index.php/Practical_Pymol_for_Beginners

기타 PyMOL의 활용법을 위한 PyMOL Wiki 바로가기 > https://pymolwiki.org/

자주 사용하는 유용한 PyMOL 명령어 > http://pymol.sourceforge.net/newman/user/S0210start_cmds.html

CYP1A1

저는 약물 유전자 중 하나인 CYP1A1 유전자에서 발견된 변이들의 위치를 아래의 Command들을 이용하여, Cartoon View로 표시해 보았습니다.  위의 그림은 최종적으로 완성한 그림입니다.

 

[References]

Senior, Andrew W., et al. “Improved protein structure prediction using potentials from deep learning.” Nature (2020): 1-5.

Rare variant association analysis: SKAT, SKAT-O, Burden test

관련 포스팅 보기>

전장 유전체 연관 분석, GWAS란 무엇인가?

Annovar: Population frequency, in-silico prediction tool 및 기타 database 활용

plink를 이용한 GWAS 분석 및 Manhattan plot 만들기

plink를 이용한 GWAS 분석에 대한 지난 포스팅에 이어서, 이번에는 SKAT을 이용한 rare variant 분석에 대한 포스팅을 정리해보고자 합니다. 유전자에 존재하는 변이(Variant)는 인구 집단 내의 분포 빈도 (Allele frequency)에 따라서, 흔한 변이 (Common variant)와 희귀 변이 (Rare variant)로 구분됩니다. 전장 유전체 연관 분석 (Genome-wide association study, GWAS)은 이러한 변이 중에서 일반적으로 인구 집단 내에 변이 빈도 분포가 5% 이상인 흔한 변이를 대상으로 하게 되는데, 대부분의 흔한 변이들은 유전학적 선택압 (Selective pressure)이 작은 변이들이기 때문에, 효과 크기 (Effect size)가 작은 것들이 대부분입니다.

 

41586_2009_Article_BFnature08494_Fig1_HTML

따라서 Common variant를 이용한 GWAS 분석으로, 설명이 안되는 유전력을 갖는 질환에 대해서 “Missing heritability“라는 용어가 등장하게 되었습니다. 이후의 후속 연구들에서는, Missing heritability를 설명하기 위해서, 변이의 빈도가 매우 드물지만, 효과 크기가 훨씬 큰 Rare variant 들이 조명을 받게 됩니다.

Missing heritability에 관한 Nature 사설 보기>

Maher, B. Personal genomes: The case of the missing heritability. Nature 456, 18–21 (2008)

Missing heritability

다만 Rare variant의 경우는, 변이 자체의 빈도가 매우 드물기 때문에 통계적으로 변이의 영향을 비교하기 위해서는 매우 큰 수의 표본이 필요한 한계가 있습니다. 이러한 한계를 극복하기 위해서, 같은 유전자 내에 존재하는 희귀 변이를 유전자 또는 특정 구역의 단위로 집합 시켜서 분석하는 방법이 제안되었는데, 이것이 Gene-level aggregation test 입니다. 즉, GWAS에서는 개별 SNP의 효과를 변이 단위로 분석이 진행되었다면 (Single variant association test), rare variant의 경우는 여러 개의 variant를 Gene 단위로 그룹화하여, Gene의 효과를 비교 분석 (Multiple variant association test)하는 분석을 진행하게 됩니다.

gene aggregation

이와 같은 분석에서 중요한 점은 ‘어떠한 기준으로 변이를 유전자 단위로 그룹화 할 것인가?’ 하는 문제가 발생하게 되는데, 일반적으로 변이 빈도의 threshold 설정 (MAF 5% vs 1%), 변이의 Functional classification 등을 사용자가 임의로 설정하여 분석을 진행하게 됩니다. 더불어, 각 희귀 변이의 효과들이 동일하지 않기 때문에 개별 변이의 효과를 보정해주는 방법으로 Rare variant allele frequency, In-silico prediction score 등을 이용하게 됩니다. 아래는 이러한 희귀 변이 분석 방법과 각 분석 방법의 장, 단점, 그리고 분석 software를 정리한 표입니다.

rare variant association

위 표의 희귀 변이를 위한 다양한 분석 도구 중에서 오늘은 R package로 개발된 대표적인 분석 도구로 사용되는 Sequence Kernal Association Test (SKAT)에 대해서, 간단히 정리하고 마치도록 하겠습니다. Burden test의 경우는 변이들의 효과 방향이 일정하고, 대부분이 원인 변이이 경우에 더 강력한 통계 검정 방법이고, SKAT의 경우는 각 변이들의 효과 방향이 제각각이거나, 원인 변이의 비율이 적은 경우에 더 강력한 통계 검정 방법입니다.  SKAT-O는 이 둘의 효과를 통계적으로 최적화하는 optimized rho value를 계산하여, 양쪽의 장단점을 모두 이용하는 방법입니다.

마지막으로 ‘SKAT’ R package를 이용을 위한 code를 공유합니다. 분석에 앞서, 이전 plink 분석 포스팅에서 언급했던 Bed, Bim, Fam 파일 및 phenotype data를 포함하는 PED 파일이 input data로 준비되어야 합니다. 마지막으로 Rare variant 변이들은 어떤 기준으로 그룹화할지에 대한 정보를 지정한 SetID 파일을 생성하여야 합니다. 자세한 내용은 SKAT 패키지의 manual을 참조하시기 바랍니다.


library(SKAT)

setwd("/plink_result")
file_name <- "my_plink"

File.Bed <- paste(file_name,".bed",sep = "")
File.Bim <- paste(file_name,".bim",sep = "")
File.Fam <- paste(file_name,".fam",sep = "")
File.SSD <- paste(file_name,".SSD",sep = "")
File.Info <- paste(file_name,".info",sep = "")
File.SetID <- paste(file_name,".SetID",sep = "")

Generate_SSD_SetID(File.Bed, File.Bim, File.Fam, File.SetID, File.SSD, File.Info)

SSD.INFO <- Open_SSD(File.SSD, File.Info)
FAM <- Read_Plink_FAM(File.Fam, Is.binary=FALSE)
COV <- Read_Plink_FAM_Cov(File.Fam,"pheno_data.ped", Is.binary=FALSE, flag1=0, cov_header=TRUE)

# continuous phenotype
obj <- SKAT_Null_Model(y ~ covariates, out_type="C")
# dichotomous phenotype
obj <- SKAT_Null_Model(y ~ covariates, out_type="D")

# SKAT
out.skat <- SKAT.SSD.All(SSD.INFO, obj)
#SKAT-O
out.skato <- SKAT.SSD.All(SSD.INFO, obj, method="optimal")
#Burden test
out.burden <- SKAT.SSD.All(SSD.INFO, obj, r.corr=1)

data <- cbind(out.skat$results,out.skato$results,out.burden$results)
write.table(data,"SKAT_results.txt",col.names=T,row.names=F,quote=F,sep="\t")

Close_SSD()

 

 

[References]

‘SKAT’ R package

SKAT GitHub: https://github.com/leeshawn/SKAT

Lee, Seunggeung, et al. “Rare-variant association analysis: study designs and statistical tests.” The American Journal of Human Genetics 95.1 (2014): 5-23.

Lee, Seunggeun, et al. “Optimal unified approach for rare-variant association testing with application to small-sample case-control whole-exome sequencing studies.” The American Journal of Human Genetics 91.2 (2012): 224-237.

Bansal, Vikas, et al. “Statistical analysis strategies for association studies involving rare variants.” Nature Reviews Genetics 11.11 (2010): 773-785.

 

plink를 이용한 GWAS 분석 및 Manhattan plot 만들기

유전체 연구에 있어서 연구 디자인 (Study Design)과 형질 (Phenotype) 은 매우 중요합니다. 특히 약물 유전체 연구에 있어서의 관심 형질은 체내 약물 농도나 대사능, 부작용의 발생 여부 등이기 때문에 더욱더 정확한 표현형을 수집하기가 어려운 점이 있습니다. 최근에 논문을 쓰면서, GWAS (Genome-wide association study)를 돌리고, Manhattan plot을 그릴 일이 있어서, 관련 분석 과정을 정리해볼까 합니다.

관련 포스팅 보기>

 

I. PLINK

1-s2.0-S0002929707613524-gr4_lrg

대부분의 GWAS는 SNP array를 이용하여, 대표 유전자 마커를 이용한 표현형 연관성 연구로 진행이 되는데, 이때 주로 사용하는 Tool이 plink입니다. (해당 tool이 논문으로 나온게 2007년이니까 벌써 10년도 넘은 소프트웨어입니다..) 그러나 아직도 쓰이고 있다는 건, 그만큼 많은 연구자들이 쓴다는 것이고, 대표적인 소프트웨어라고 할 수 있습니다. (1.9 버젼이 나온 이후, 2.0 버젼을 베타 테스트하고 있다고 한지도 꽤 오래 되었는데, 그 이후 업데이트가 매우 느리게 진행되고 있는 것이 단점입니다.) 물론, BI tool 답게 많은 경쟁 소프트웨어들이 나왔는데 (ex. EPACTS), 아직도 대부분의 논문에서 plink를 쓰는 것을 보면, 대부분의 분석을 하는데 plink만 있어도 크게 무리가 없기 때문이 아닐까 합니다. plink의 사용법은 plink 홈페이지 (PLINK: Whole genome data analysis toolset)의 tutorial page에 매우 자세하게 소개가 되어 있어서, 그때 그때 필요한 내용들을 찾아서 쓰면 됩니다.

plink 다운로드 및 설치>

<VCF 파일 압축 및 인덱싱>

bgzip -c [myvcf.vcf] > [myvcf.vcf.gz]
tabix -p vcf -f [myvcf.vcf.gz]

<PED, MAP 파일 or BED, BIM, FAM 파일 만들기>

plink 실행을 위해서는 PED & MAP file 또는 binary 형식으로 변환된 BED, BIM, FAM file이 필요합니다. 일반적으로 SNP array 데이터를 생산하면 만들어주기 때문에 따로 준비할 필요는 없습니다. 간혹 NGS로 생산된 시퀀싱 데이터로 plink로 실행하고 싶은 경우, vcf 파일을 위의 형식으로 변환하면 좋은데, 아래와 같은 command가 유용합니다.

plink --noweb --vcf [myvcf.vcf.gz] --recode --out myplink
plink --noweb --vcf [myvcf.vcf.gz] --recode --make-bed --out myplink

<PLINK 파일 기본 QC>

plink --file myplink --missing-genotype N --make-bed --mind 0.05 --maf 0.05 --geno 0.1 --hwe 1e-6 --recode --out myplink.QC

missing genotype 여부, genotyping calling rate, minor allele feqeuncy, HWE (Hardy-Weinberg equilibrium) cut-off 기준으로 이를 위반하는 SNP들은 모두 날려버리는 quality control 과정입니다.

<plink를 이용한 연관 분석>

plink를 이용한 연관 분석은 통계 모형에 기반하기 때문에 우선적으로 어떤 모델을 이용하여, 어떻게 분석을 할지를 고려해야 합니다. Genetic inheritance mode (Additive, Dominant, Recessive)와 분석 형질이 Dichotomous trait인지 Continuous trait 인지에 따라서 Case-control, linear regression, logistic regression model 등을 적용할 수 있습니다. 더불어, 보정을 위한 공변량(covariate)으로 무엇을 선택할 것인지도 중요합니다.

분석을 위한 Input Phenotype data를 준비하는 과정도 중요한데, 다음 페이지에서 자세하게 소개가 되어 있습니다.

[분석을 위한 command]

plink --noweb --bfile [mydata] --[additive/dominant/recessive] --[assoc/linear/logistic] --pheno [phenotype_file] --pheno-name [phenotype_name] --covar [covariates_file] --covar-name [covariates_name] --out [result_file]

위의 command에 적절한 inheritance mode [additive/dominant/recessive]와 분석 모델 [assoc/linear/logistic]을 골라서, 분석을 실행하면 됩니다. 특정 SNP에 대한 Conditioning을 원하는 경우, –cond [SNP ID]를 추가합니다.

위의 분석 과정을 거치면, 모든 SNP 위치에 대한 Beta 및 P value가 계산됩니다. Beta는 해당 SNP의 Effect size를 나타내는 통계량이고, P value는 해당 SNP의 통계적 유의도를 의미합니다. Manhattan plot은 일반적으로 여기서 계산된 P value에 -log를 취한 형태로 그리게 됩니다.

II. Manhattan Plot 그리기

Manhattan Plot을 그리는 방법도 다양하지만, 여기서는 제일 간편한 qqman R package를 이용하도록 하겠습니다. 자세한 option은 아래 Reference의 자료들을 참고 바랍니다.

library(qqman)

## plink 결과 파일 불러오기
data &amp;amp;lt;- read.table("plink_result", header = T, stringsAsFactors=F)

## Manhattan plot 그리기
manhattan(data, main = "Manhattan Plot", ylim = c(0, 40), cex = 0.8, cex.axis = 0.9, col = c("grey", "skyblue"))

## QQ plot 그리기
qq(data$P)
GWAS

위의 패키지를 이용하면, 위와 같은 Manhattan plot을 손쉽게 만들 수 있습니다.

III. HaploView

마지막으로, SNP 정보의 linkage 여부에 따른 LD block의 시각화를 위한 Haploview에 대해서 간단히 정리하고, 포스팅을 마치도록 하겠습니다.

Haploview 4.2 Download

plink --noweb --bfile [mydata] --extract [Gene_SNP_list] --recodeHV --out [Gene_haploview]

plink의 위의 command를 이용하여, Haploview를 원하는 SNP의 list에 대해 ped 및 info 파일을 생성합니다. 이를 HaploView 프로그램을 통해 loading해주면, 생성된 LD block과  계산된Haplotype 조회가 가능합니다. 아래 그림은 HaploView를 이용하여, 생성된 LD block 입니다.

Figure S5

[References]

PLINK: Whole genome data analysis toolset

Purcell, Shaun, et al. “PLINK: a tool set for whole-genome association and population-based linkage analyses.” The American journal of human genetics 81.3 (2007): 559-575.

Chang, Christopher C., et al. “Second-generation PLINK: rising to the challenge of larger and richer datasets.” Gigascience 4.1 (2015): s13742-015.

qqman R package GitHub

Manhattan plot in R: a review

NGS DNA-seq pipeline: GATK Best Practice Code – Part3. Vcf manipulation

앞선 포스팅의 두가지 과정을 거치고 나서 생성된 VCF 파일을 이용하면 드디어 분석 가능한 변이들을 확인할 수 있습니다. 그러나 실제로 이 데이터를 열어보면, 지저분하고 활용하기 위해서는 어느 정도 가공이 필요합니다. 그래서 이번 포스팅은 VCF Filter를 적용하여 분석을 위한 변이들을 정제하고, 분석에 참조하기 위한 Annotation 작업을 위한 Code까지 정리해보도록 하겠습니다.

관련 포스팅 보기>

NGS 분석 파이프 라인의 이해: GATK Best Practice

NGS DNA-seq pipeline: GATK Best Practice Code – Part1. Fastq to Bam

NGS DNA-seq pipeline: GATK Best Practice Code – Part2. Bam to Vcf

vcf_3

[일반적인 VCF 파일의 구조]

IV. Variant Filtration and Annotation : Genotype Refinement Workflow

VCF call을 하고 나면, 많은 변이 데이터가 생산되는데, 실제로 데이터의 퀄리티와 오류를 보정하는 작업이 필요합니다. 이 부분은 아직까지 확립된 최적의 Work Flow는 없지만, GATK에서는 Genotype Refinement Workflow를 제공하고 있습니다. 여기서는 Call된 genotype이 인구집단이나 가계도 정보에 근거하여 적절한 Call인지를 판단하게 됩니다. 가계도의 경우에는 양쪽 부모 정보가 모두 (Trio) 필요합니다. 아래에서는 이 Workflow를 따라 Code를 작성해보도록 하겠습니다.

Genotype refinement

1. CalculateGenotypePosteriors

VCF call을 통해 생성된 genotype을 기반으로 이번에는 거꾸로 Genotype의 posterior를 계산해줍니다. 이 결과를 토대로 변이 call을 판정해서, Genotype의 quality를 계산합니다. 이 과정에는 인구 집단이나 가계도의 Prior 정보를 이용하는데, 가계도 정보를 이용하기 위해서는 PED 파일 생성이 필요합니다.

PED 파일 생성 정보는 아래 링크를 참고하시기 바랍니다.

http://zzz.bwh.harvard.edu/plink/data.shtml

gatk --java-options 'Xmx16g' 'Xms8g' CalculateGenotypePosteriors -V [cohort.recal.vcf] \
-ped family.ped --skip-population-priors -O [cohort.GP.vcf]


2. Filter Low Quality Genotypes


위의 과정에서 계산된 GQ 값을 기준으로 Filter 조건을 걸어서 지저분한 변이 call들을 제외 시킬 수 있습니다. 아직까지 표준화된 Filter 조건은 없기 때문에, 사용자가 적절하게 필터 조건을 설정해주어야 합니다. 여기서 선언해주는 조건은 흔히 JEXL expression이라고 부르는 방법으로 작성을 해야한다고 하는데, 여기서는 널리 쓰이는 Filter option인 GQ < 20 을 걸도록 하겠습니다.  필요에 따라 수치를 변경시키거나, 조건을 바꾸면 됩니다. 자세한 내용은 아래에서 찾아보시기 바랍니다.


article about using JEXL expressions


https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_filters_VariantFiltration.php


마지막으로 filter조건에 부합하는 변이만 남기기 위해, GATK의 SelectVariants 이용하도록 하겠습니다.


gatk --java-options 'Xmx16g' 'Xms8g' VariantFiltration -R [hg19_reference.fa] -V [cohort.GP.vcf] \
--genotype-filter-expression "GQ<20" --genotype-filter-name "lowGQ" -O [cohort.filter.vcf]
gatk --java-options 'Xmx16g' 'Xms8g' SelectVariants -V [cohort.filter.vcf] --exclude-filtered true --max-nocoll-fraction 0.1 -O [cohort.filtered.vcf]


3. VCF annotation using Annovar


이제 어느 정도 정제가 된 변이 데이터셋이 생성되었습니다. 이제 불러온 변이들에 맞는 annotation 작업을 진행하면 분석을 진행할 수가 있습니다. 최근에 GATK에서 Funcotator를 개발하여, annotation이 가능해졌는데 아직까지는 Annovar를 이용하는게 대세이므로, Annotation 작업은 Annovar를 이용하도록 하겠습니다.


관련 포스팅 보기>


Annovar: Population frequency, in-silico prediction tool 및 기타 database 활용


Annovar 홈페이지를 통해 다운로드 신청을 하면,  등록된 메일로 Kai Wang에게서 자동으로 다운로드 url을 받을 수 있습니다. Annotation을 원하는 Database의 크기가 크기 때문에, 처음에는 DB 다운로드에 시간이 많이 소요가 될 수 있습니다.


Annovar를 설치하고, 아래 부위에 원하는 DB의 이름을 찾아서 넣으면 다운로드가 진행됩니다.


annotate_variation.pl -buildver hg19 -downdb -webfrom annovar [refGene] humandb


Annotation은 아래 code를 통해서 진행합니다. 여기서는 refGene, cytoBand, avsnp150, dbnsfp35a, exac03, gnomad_exome, clinvar_20190305 DB를 이용하여 annotation 작업을 진행하도록 하겠습니다. 각각의 DB는 Gene-based, Region-based, Filter-based annotation에 해당하는 카테고리가 있고, 이를 operation 옵션에 g, r, f를 이용하여 나타내줍니다.


table_annovar.pl [cohort.filtered.vcf] humandb/ -buildver hg19 -out [myanno] -protocol refGene,cytoBand,avsnp150,dbnsfp35a,exac03,gnomad_exome,clinvar_20190305 -operation g,r,f,f,f,f,f -nastring . --vcfinput --remove --thread 8


[References]


VCFtools를 이용한 VCF filtering: 위에서는 GATK를 이용하였지만, 과거에 많이 사용하였던 VCFtools를 이용하는 방법도 있습니다. 자세한 내용은 아래를 참조하시기 바랍니다.


http://www.ddocent.com/filtering/


GriffithLab에서 소개하고 있는 VCF filter 과정: 해당 내용은 GATK4로 넘어오면서 코드가 약간 달라졌습니다. 다만 전반적인 과정은 유사하므로, 참고하시기 바랍니다.


https://pmbio.org/module-04-germline/0004/02/02/Germline_SnvIndel_FilteringAnnotationReview/

NGS DNA-seq pipeline: GATK Best Practice Code – Part2. Bam to Vcf

지난 글에 이어서, 생성된 Bam 파일로부터 변이들을 읽어 들이고, haplotype call을 하는 두번째 파트의 code를 정리해보겠습니다. 아래 코드는 GATK 4.1.3 버젼을 기반으로 작성되었습니다. GATK 버젼에 따라서 조금씩 Tool과 명령어에 차이가 있습니다. 전반적인 흐름은 아래와 같습니다.

관련 포스팅 보기>

NGS DNA-seq pipeline: GATK Best Practice Code – Part1. Fastq to Bam

NGS 분석 파이프 라인의 이해: GATK Best Practice

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

1.jpg

III. Germline short variant discovery : Bam to Vcf

1. HaplotypeCaller

https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_haplotypecaller_HaplotypeCaller.php

해당 위치의 변이와 read 갯수등을 바탕으로 하여, haplotype을 생성해주는 과정입니다. 일반적으로 우리가 heterozygote인지 homozygote인지에 해당하는 변이의 zygosity를 이 과정에서 생성한다고 볼 수 있습니다. 앞선 과정에서 BQSR을 통해 최종적으로 생성된 Bam 파일을 이용하여, haplotype call을 진행합니다.

gatk -- java-options 'Xmx16g' 'Xms8g' HaplotypeCaller&nbsp;-R [hg19_reference.fa] -I [sample01.final.bam] -ERC GVCF [-L targets.interval_list] \ -O [sample01.vcf] --genotyping-mode DISCOVERY --standard-min-confidence-threshold-for-calling 30 -ploidy 2 --output-mode EMIT_VARIANTS_ONLY


-L option에는 타겟 영역에 해당하는 BED 파일을 넣어주면 해당 위치에서만 Call이 일어나게 됩니다. Whole genome이 아닌 Target panel이나 Exome panel의 경우, 타겟 유전자들로 구성된 BED 파일을 지정해줍니다.


2. GenomicsDBImport


https://github.com/GenomicsDB/GenomicsDB/wiki


이전의 GATK와 달라진 부분입니다. 예전에 여러개의 샘플로 구성된 VCF를 Genotype GVCF로 합쳤는데, 추후의 편의성 및 연산 속도를 고려하여 GenomicDB를 구축하는 과정이 추가되었습니다. 이 단계를 건너 뛰어도 좋지만, 일반적으로 이렇게 DB를 구축하고 나면 여러모로 좋다고 소개하고 있네요. 여러개의 샘플로 구성된 환자 코호트의 경우는 이렇게 DB를 구축하고 하나로 합쳐주는 과정이 있는게 좋을 것 같고, Single sample이라면 이 과정은 건너뛰어도 좋을 것 같습니다.


gatk -- java-options 'Xmx16g' 'Xms8g' GenomicsDBImport&nbsp;-V [vcf file_list] \ [-L targets.interval_list] --genomicsdb-workspace-path \ [/genomicDB] --merge-input-intervals true --tmp-dir=/tmp&nbsp;


[vcf file_list] 부분에는 -V sample01.vcf -V sample02.vcf -V sample03.vcf 와 같이 하나의 cohort로 구성하고자 하는 여러개의 샘플을 하나로 쭉 나열해주면 됩니다.


3. GenotypeGVCFs


위에서 GenomicDB를 구축하고, 여러 개의 샘플로 구성된 VCF가 있었다면, 이들을 Cohort 단위로 묶어서 하나의 VCF로 만들어 주는 과정입니다.


gatk -- java-options 'Xmx16g' 'Xms8g' GenotypeGVCFs -R [hg19_reference.fa] -V gendb://genomicDB --tmp-dir=/tmp -O [cohort.vcf]


4. VariantRecalibrator, ApplyRecalibration


https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_vqsr_VariantRecalibrator.php


https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_vqsr_ApplyVQSR.php


Bam을 구성할 때, 각 염기별로 recalibration을 진행했던 것과 비슷하게 call된 변이들에서도 recalibration을 진행해서 QC를 하는 부분입니다. 머신러닝을 이용하여, 기존 DB로 부터 변이들을 학습 시키고 걸러내는 방법이라고 하는데, 이를 위해서 hapmap, omni, 100G, dbsnp database를 다운 받아야합니다. 최근에는 CNN (convolutional neural network)에 기반한 모델을 수립하여, 테스트 중이라고 하는데 아직까지는 Beta 버젼에 머물고 있어 기존 코드를 이용합니다. 추후 이 부분은 CNN 기반 Recalibration code로 바뀔 가능성이 있습니다.


gatk -- java-options 'Xmx16g' 'Xms8g' VariantRecalibrator -R [hg19_reference.fa] -V [cohort.vcf] \
--resource:hapmap,known=false,training=true,truth=trueprior=15.0 hapmap_3.3.hg19.sites.vcf.gz \
--resource:omni,known=false,training=true,truth=false,prior=12.0 1000G_omni2.5.hg19.sites.vcf.gz \
--resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G_phase1.snps.high_confidence.hg19.vcf.gz \
--resource:dbsnp,known=true,training=false,truth=false,prior=2.0 Homo_sapiens_assembly19.dbsnp138.vcf.gz \
-an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR \
-mode SNP -O [cohort.recal.vcf] --tranches-file [output.tranches] --rscript-file [output.plots.R]

gatk ApplyVQSR -R [hg19_reference.fa] -V [cohort.vcf] -O [cohort.final.vcf] \
--truth-sensitivity-filter-level 99.0 --tranches-file [output.tranches] --recal-file [cohort.recal.vcf] -mode SNP 


위의 과정을 거치면, 최종적으로 변이들이 call 되어 분석이 가능해집니다. 마지막 파트는 Annotation 과정인데, 해당 코드는 마지막 Part 3 포스팅에서 다루도록 하겠습니다.


[Reference]


GATK Best Practice for Germline short variant discovery (SNPs + Indels)


: https://software.broadinstitute.org/gatk/best-practices/workflow?id=11145


How to run VQSR from GATK


https://software.broadinstitute.org/gatk/documentation/article?id=2805