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.

 

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

오늘 포스팅은 Annovar에 관한 내용으로 작성해볼까 합니다. 언젠가 한번은 정리할 필요가 생각하고 있던 내용인데, 이제서야 차례가 왔네요. 시퀀싱 기술이 보편화되면서, 대부분 표준화된 파이프라인을 이용하여 시퀀싱 raw data로부터 변이를 검출하는 Up-stream analysis는 대부분의 업체에서 대동소이한 결과를  주고 있습니다. 연구자의 입장에서 사실 더 중요한 것은 이 정보를 어떻게 효과적으로 이용할까 하는 부분에 있고, 그래서 Bam file 또는 VCF file에서 시작하는 Down-stream analysis가 더 중요하다고 볼 수 있습니다. 그 첫 단계로 필요한 적절한 정보를 활용하여 주석을 다는 과정이 있고, 이를 우리는 Annotation (주석 달기) 과정이라고 합니다. 그리고 이 과정에서 우리는 대부분 Annovar를 사용하게 됩니다. 물론 업체에 의뢰하면 대부분 기본적인 Annotation이 끝난 파일도 전달을 받게 되는데, 사실 필요 없는 내용이 잔뜩 달려서 파일의 용량만 무지막지하게 커진다거나, 정작 필요한 내용이 빠진 경우도 종종 발생합니다. 그래서 이번에 다룰 내용은 주석 달기의 각 항목에 대한 의미와 주요 활용 항목에 대해서 정리해보겠습니다.

관련 포스팅 보기>

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

Annovar는 Perl script 기반으로 짜여 있어서, linux에서 명령어를 주면 자동으로 annotation을 달아주게 되어있습니다. 이때 몇가지 option에 따라서 원하는 내용대로 다양한 Database에서 관련 내용을 annotation 할 수가 있습니다. 더 많은 내용을 추가함에 따라서 연산 시간과 최종 파일 용량은 커지고, 가끔은 이것이 오히려 독이 되기도 합니다. 따라서 원하는 내용을 입맛에 따라 최적화하여 활용하는 것이 필요하게 됩니다. 아래 그림은 Annovar를 이용하여, 주석을 다는 과정의 전반적인 흐름을 보여주고 있습니다. 간단히, “Database 파일을 통째로 다운 받아서, Perl script 기반의 annotate_variation.pl을 실행하면, 그에 맞는 내용을 찾아서 주석으로 달아준다.” 정도로 말할 수 있겠습니다.

ANNOVAR_main_package_workflows.svg.png

Annovar에서 제공하는 Database 목록 보기

위 페이지에서는 다운로드 가능한 Database의 목록들을 보여주고 있습니다. 필요한 항목의 이름과 업데이트 날짜 등을 참고하여, DB를 다운로드하고 적절하게 활용하기 바랍니다.

annotate_variation.pl -buildver hg19 -downdb -webfrom annovar [위 목록의 Database 이름] humandb/

 


I. Gene-based Annotation

VCF 파일에서 가장 핵심 정보는 사실 몇개 없습니다. 정확하게는 5개 정보만 있어도 충분한데, “몇번째 염색체 (Chr)의 몇번째 염기 위치 (Position)가 원래 무엇인데 (Ref) 무엇 (Alt)으로 바뀌었다.” (흔히 이러한 형태의 파일을 MAF format이라고 합니다.) 이 말을 하기 위해서, 부수적인 내용들이 잔뜩 달려 있는 셈이죠. 그러나 이 정보만을 이용해서는 사람들이 알아들을 수가 없습니다. 우리는 “어떤 유전자의 몇번째 엑손 영역의 몇번째 아미노산이 무엇으로 바뀌었다“와 같은 정보가 필요하죠. 그런데 사실 아직도 유전자의 정확한 정의와 위치에 대해서는 100% 밝혀지지 않았습니다. 따라서 이 유전자라고 하는 부분도 다양한 데이터 베이스가 존재하게 됩니다. 대표적으로 RefGene, UCSC/Ensemble Gene, Known Gene, CCDS 등등의 데이터 베이스가 있습니다. 따라서 여기서 어떤 데이터 베이스를 이용하여, annotation을 하냐에 따라서 출력이 달라집니다. (하지만 사실 큰 차이는 없습니다. 대부분의 업체에서는 2~3개의 데이터 베이스를 이용하여 annotation을 해주는데, 대부분 내용이 중복되어 용량만 커짐… 그래서 개인적으로는 그냥 RefGene만으로도 충분합니다.)

<가장 핵심적인 VCF 파일의 정보> = MAF format

Chromosome : Position (Start_End) : Reference sequence > Alternative sequence

위의 database를 이용하면, 위의 정보가 어떤 유전자에 속하고, 해당 유전자에서 어떤 기능을 하는 어떤 부위의 변이인지, 기능적으로 변화가 있는지 없는지 등에 대한 기본적인 정보를 제공해주게 됩니다. 어떻게 보면 가장 핵심적인 정보를 추가하는 부분이라고 할 수 있습니다.

 

II. Filter-based Annotation

사실 변이를 Genome Browser에서 찾는 가장 빠른 방법은 rsID를 이용하는 것입니다. 그런 점에서 dbSNP 또는 avSNP의 rsID를 주석으로 달아놓는 것은 활용도가 높습니다. rsID는 변이 보고가 점점 늘어남에 따라서 계속 갱신되고 있는데, 가장 최근 database는 avSNP 151 버젼이지만, 보편적으로 아직까지는 avSNP 147 버젼을 사용하고 있는 것 같습니다.

관련 포스팅 보기>

바이오 연구자를 위한 Genome Browser 비교 및 활용

아미노산 치환의 효과 예측: In silico tool의 원리와 종류

변이빈도와 효과 크기

사실 변이의 생물학적 의미를 해석하는 과정에서 정말로 중요한 내용은 일반 인구 집단에서 얼마나 흔하게 존재하는가?에 있다고 할 수 있습니다. 그러한 의미에서 Population Frequency 정보는 변이의 의미를 파악하는데 매우 중요한 단서를 제공해 줍니다. 이와 같은 Population Frequency는 다양한 집단에서 시퀀싱을 시행하여, 그 빈도를 계산한 다양한 데이터 베이스에 기반하게 되는데, 대표적으로 EXAC, gnomAD, 1000Genome, ESP database 등이 있습니다. 이때 중요한 점은 인구 집단이 얼마나 대표성을 띄는가?에 있다고 볼 수 있는데, 빈도를 계산한 인구 집단의 크기 (n수)와 인종적 배경 (Ethnicity) 등이 특히 중요합니다. 더불어, Rare variant의 경우는 더욱더 인종에 따른 차이가 크기 때문에 일반적으로는 한국인의 경우는 EAS (East Asian population) 정보를 이용하지만, 해당 데이터 베이스가 충분히 한국인을 대표하지 못하는 경우가 많습니다. (위의 데이터 베이스에서는 gnomAD가 가장 n 수가 가장 크기 때문에 저는 주로 EXAC과 gnomAD의 EAS 인구 집단 정보를 활용하고 있습니다.)

사실 데이터 용량이 가장 뻥튀기 되는 부분이 in-silico prediction tool 부분에 있습니다. 아미노산 치환 효과를 예측해주는 tool은 100가지가 넘게 있는데, 각 tool 들이 예측해주는 정보를 얼마나 포함할 것인가에 따라서 추가되는 정보도 달라집니다. 크게 얼마나 진화적으로 보존된 지역인지 또는 아미노산 구조적으로 변화를 유발하는지 등에 기반한 알고리즘으로 개발되어, SIFT, PolyPhen과 같은 고전적 tool부터, GERP, SiPhy, MutationTaster, FATHMM, MetaSVM, CADD, DANN과 같은 다양한 tool이 존재합니다. (역시나 이쪽도 절대 지존은 없기 때문에 경우에 따라 적절하게 활용합니다. 보통 Voting Method, 즉 10개 또는 선택된 갯수의 서로 다른 알고리즘 중에서 몇개가 Deleterious로 예측하는가? 와 같은 방식으로 많은 경우 활용하게 됩니다. 저는 개인적으로 Ensemble 기반의 tool을 주로 활용하고 있습니다.)

마지막으로 임상 정보들을 annotation 하는 database가 존재하는데, 대표적으로 OMIM, HGMD, ClinVar, COSMIC 등등이 있습니다. 희귀 유전 질환에 대해서 연구를 한다면, OMIM이나 HGMD, Cancer 관련 연구를 한다면, COSMIC database 정보를 annotation 하는 것이 도움이 될 수 있습니다. 그러나 사실 이렇게 annotation을 덕지덕지 붙이다보면 파일 크기가 엄청나게 불어나게 됩니다. 현재까지 대부분의 연구자들은 Coding region의 Functional variant에만 관심이 있기 때문에, 1차적으로 Gene-based annotation 후 exon 영역의 functional variant만 filtering하고나서, 해당 변이들에 대해서 annotation 하는 방법이 시간과 데이터를 절약하는 방법이 될 수 있습니다.

 

III. 기타 annotation 방법

관련 포스팅 보기>

암유전체 분석: Driver mutation prediction tools

위의 annovar를 이용하는 방법은 linux 기반의 서버를 통해서 대용량으로 실행하는 방법입니다. 그러나 서버를 구축하지 못하거나, linux를 친숙하게 이용하지 못하는 경우에는 그러면 어떻게 annotation을 하는가? 에 대한 문제가 발생합니다. 이를 위해서 다양한 Web 기반의 annotation tool 들이 존재하게 됩니다. 가장 대표적인 것이 wANNOVAR입니다. 기타 cancer를 다룬다면, Oncotator 또는 Cancer Genome Interpreter도 대안이 될 수 있습니다. 그러나 역시 이러한 tool들은 Annovar에 비해서 자유도는 떨어지기 때문에 기능에 제약이 있다는 단점이 있습니다. 마지막으로 R을 활용하여, annotation이 가능한 몇가지 package들이 개발되어 있습니다. 대표적인 package로는 MAFtools, VariantAnnotation 등이 있으나, 역시 기능이 AnnoVar에 비하면 제한적입니다. 그러나, 소수의 변이에 대해서 빠르게 annotation이 필요한 경우라면 이러한 도구들도 적절하게 활용하는게 도움이 될 수 있습니다. 더 자세한 정보는 아래 github를 활용하시기 바랍니다.

 


 

References>

Yang, Hui, and Kai Wang. “Genomic variant annotation and prioritization with ANNOVAR and wANNOVAR.” Nature protocols 10.10 (2015): 1556.

ANNOVAR: http://annovar.openbioinformatics.org/en/latest/

wANNOVAR: https://github.com/WGLab/doc-ANNOVAR

Oncotator: https://github.com/broadinstitute/oncotator

VariantAnnotation: https://github.com/Bioconductor/VariantAnnotation

MAFtools: https://github.com/PoisonAlien/maftools

Web resources for Bioinformatics database: https://netbiolab.org/w/Web_Resources

NGS 타깃 시퀀싱 패널 검사의 분석 및 해석시 고려할 사항

지난 포스팅에서는 임상의의 입장에서 NGS 검사를 통한 변이의 해석을 이해하기 위한 기초적인 내용들을 언급했습니다.

<관련 포스팅 보기>

임상의를 위한 NGS 레포트 해석의 이해

NGS 검사: Whole Genome & Exome, Targeted Sequencing 비교

그러나, 언급된 내용들이 기초적이고 많이 부족하다고 느껴, 이번 포스팅에서는 타깃 시퀀싱 패널의 분석시 고려할 내용 및 팁 등을 언급해 보고자 합니다.

1) 검사 데이터의 Quality check: 사실 환자를 보는 의사의 입장에서는 최종 결과만 확인하기 때문에 가장 간과하기 쉬운 단계입니다. 그러나 본인이 직접 데이터를 분석하고 해석하고자 한다면, 가장 중요한 단계입니다. 확인해야할 여러 파라미터들이 있습니다만, 가장 기본적으로 target region의 coverage 및 depth를 확인해야합니다. 이는 우리가 검사하고자 하는 영역을 타깃 시퀀싱 패널이 얼마나 잘 디자인되어 검출하는지를 나타내주는 지표입니다. 경우에 따라 다르지만, 일반적으로 100X 이상의 depth로 원하는 영역의 99% 이상 커버(100X over target ratio > 99%)한다면 디자인이 매우 잘된 것으로 생각할 수 있습니다.

doc-1024x216

1
매우 잘 디자인된 타깃 시퀀싱 패널의 성능 결과 예: 평균 depth도 매우 많고, 원하는 영역을 골고루 잘 커버하는 것을 알 수 있다.

2) 검출된 변이가 true signal인가? false positive 인가?: 위와 같이 잘 디자인된 시퀀싱 패널이라 하더라도, micro-insertion 또는 deletion이 일어나는 경우에는 reference 패널과 read의 시퀀스 차이가 크기 때문에 엉뚱한 곳에 가서 read가 붙는 일이 일어날 수 있습니다. 그 결과 엉뚱한 위양성 변이가 검출되기도 합니다.

따라서 임상적으로 중요한 변이가 검출된 경우에는 (특히 frameshift mutation), 실제로 해당 변이를 IGV 와 같은 genome viewer를 통해서 직접 눈으로 확인해보는 것이 중요합니다.

2
국내의 모 NGS 검사 업체에서 전달받은 결과에서 보고한 한 환자의 검체에서 무더기로 검출된 변이. 검출된 변이의 빈도가 1% 정도로 매우 적기 때문에 noise signal로 판단할 수 있었지만, 이렇게 작은 빈도라 하더라도 somatic mutation을 타깃으로 하는 cancer panel이었다고 하면 이야기는 완전히 달라집니다.
Mutation-visualization-in-the-Integrative-Genomics-Viewer-IGV-browser-showing
경우에 따라서는 위와 같이 IGV를 통해 실제로 deletion 된 영역이 있고, read들이 올바르게 달라붙었는지를 눈으로 직접 확인해보는 것도 필요합니다.

3) 검사 목적에 따른 변이의 filter 전략: 타깃 시퀀싱 패널은 크게 2가지 유전 질환과 관련하여 생식 세포 돌연변이 (germ-line mutation) 또는 de novo mutation을 검출하거나 암 환자에서 체세포 돌연변이 (somatic mutation)을 검출할 목적으로 디자인됩니다. 이 두 가지는 구분하여 NGS 검사에서 검출된 변이를 적절하게 필터링하는 전략이 필요합니다.

<관련 포스팅 보기> 유전학 중요개념 정리: Germline vs. Somatic mutation

Germ-line의 경우에는 부모로 부터 한쌍씩 유전형을 물려받기 때문에 검출되는 변이의 상대 빈도는 ~50% 또는 ~100%일 수 밖에 없습니다. 반면, 암 세포의 경우에는 다양한 변이들이 섞여 있기 때문에 (tumor heterogeneity) 다양한 상대 빈도로 검출이 됩니다. 따라서, 변이를 필터링할 경우, 이러한 점을 염두해 두고 환자들에서 의미있는 병적 변이들을 검출하게 됩니다.

4) 집단 내 변이 빈도에 따른 filter 전략: 매우 드문 희귀 유전 질환의 변이를 검출하고자 하는 경우, 해당 변이의 집단 내 변이 빈도에 따라 필터링하는 전략도 유효합니다. 1000 Genome project 또는 ExAC과 같은 유전체 database는 인구 집단에서 해당 변이의 빈도에 대한 정보를 제공해주기 때문에, 이미 알려진 변이 빈도를 기반으로 인구 집단에서 흔하게 존재하는 변이(1% 이상)는 필터링하고 남은 변이들을 대상으로 임상적 평가를 확인하는 것이 좋습니다.

변이빈도와 효과 크기
일반적으로 희귀 유전질환의 경우에는 집단 내 변이 빈도가 매우 작고, 효과 크기가 큰 변이에 의해 유발되는 것으로 생각되고 있기 때문에, NGS 검사를 통해 검출된 흔한 변이들은 크게 임상적인 의미가 없는 경우가 많습니다.

5) 유전체 데이터 베이스를 활용한 Clinical annotation: 최근 다양한 생명정보학 및 유전체 툴들이 개발되어 검출된 변이의 특성 및 정보들을 자동으로 처리해주고 있습니다. 이러한 툴들을 보조적으로 잘 활용하면 변이 판독에 소요되는 시간을 최소화하면서, 동시에 효과적으로 판정할 수 있습니다. 그러나, 경우에 따라서 이러한 판정들이 항상 옳은 것은 아니기 때문에 변이의 판정 및 판독에는 다양한 경험을 바탕으로 한 전문가의 수기 판독이 필요하게 됩니다.

가장 대표적으로 많이 사용되는 툴은 Annovar이며, annovar 내에서도 다양한 툴들이 존재하기 때문에, 사용자가 필요에 따라 이를 선택할 수 있습니다.

Annovar 홈페이지 방문하기