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.

 

약물 유전체 연구가 어려운 이유

저는 작년 2월부터 1년 반정도의 기간을 약물유전체 연구를 하고 있습니다. 지도 교수님이신 이민구 교수님과 다양한 약물 반응에 대한 유전적 바이오 마커를 발굴하는 연구를 하고 있는데, 생각보다 쉽지가 않고, 좋은 결과가 나오지 않고 있습니다. 그에 비해 최근에 선천성 기형의 일종인 두개골 조기유합증이라는 희귀질환에 대해 성형외과 및 신경외과와 공동연구에도 참여하고 있는데, 많은 환자들의 유전적 원인들을 확인할 수 있었습니다. 유전적 소인과 형질 간에는 어떠한 연관이 있는 것일까요? 이번 글은 흔히 말하는 Common diseaseRare disease 의 차이와 더불어, 지난 1년반정도의 기간을 약물 유전체 연구를 하며 느낀 점들과 왜 약물 유전체 연구가 어려운지에 대해서 정리해보고자 합니다.

기본적으로 약물 유전체 연구는 크게 여러 사람들이 동일한 약물을 먹었을 때 혈중 유효 농도가 다양하게 나타나는 것에서, 어떤 유전적 차이가 이러한 약물 대사에 기인하는지부작용 발생 유무의 위험도를 예측할 수 있는 유전적 바이오마커가 있는지에 관심을 갖춰 연구되고 있습니다.

관련 포스팅 > 약물 유전학은 왜 정밀의료에서 중요한가?

slide_42

I. 약물 반응은 복합 형질 (Complex trait)이다 : 기본적으로 약물의 대사 과정에는 다양한 약물 효소가 관련합니다. 또한 약물이 흡수되어 배출되기까지의 대사 과정 (ADME) 또는 약동학 (Pharmacokinetics) 과정에는 다양한 요소들이 관여하기 때문에, 한 두가지 유전적 소인이 형질에 결정적 차이를 나타내기 어렵습니다. 복합 형질로 가장 많이 연구되는 질병 중 하나가 2형 당뇨병 (Type 2 Diabetes mellitus; T2DM)인데, 당뇨병 발생의 원인과 그 유전적 요인에 대해서 많은 연구가 진행되었지만 여전히 속 시원한 유전적 원인에 대해서는 알지 못하고 있습니다. 특히 이러한 복합 형질에서 발굴된 유전적 마커들은 형질의 차이에 기여하는 정도가 매우 작아서, 대부분의 효과 크기 (Effect size)가 매우 작습니다. 그래서 그나마 연구가 잘되고 결과가 잘 나오는 것은 효과 크기가 매우 큰 한 두가지의 유전적 인자가 약물의 부작용 발생 유무에 영향을 미치는 경우입니다.

II. 약물 반응의 측정 자체가 어렵다 : 체내 약물 대사능에 영향을 주는 유전적 인자를 확인하고자 하는 연구의 경우, 일단 환자에서 해당 약물 농도 측정 자체가 매우 어렵습니다. 현실적으로 환자들에게는 의사들이 체중이나 대사능 등을 고려하여 약을 처방하기 때문에 복용한 약물의 양도 간격도 전부 달라지게 되며, 약물 농도라는 것도 매우 변동성이 심하기 때문에 언제 채혈하였는지, 다른 약과 함께 복용하였는지 (drug-drug interaction), 음주 & 흡연 여부, 성별 등 다양한 요소에 영향을 받게 됩니다. 기본적으로 이러한 요소들에 대한 명확한 통제가 어렵고, 보정을 한다고 하더라도 그 측정 약물 농도가 명확하게 그 사람의 약물 대사능을 대변하지도 못합니다. 즉, 처음부터 얻어지는 정보 자체에 매우 큰 변동성이 있기 때문에 해당 데이터와 유전적 정보 간의 연관성을 찾으려고 해도, 그 영향이 명확하게 큰 경우가 아니면 연관성을 찾기가 매우 어렵습니다.

III. 약물 대사 경로에는 다양한 대체자가 존재한다. : 이 세상에는 정말로 다양한 약물이 존재합니다. 기본적으로 약물은 간에서 대사되어 신장을 통해 배설된다고 알려져 있습니다만, 약물 개별로 보면 어떤 약물이 정확하게 어떠한 효소에 의해 대사되어 어떠한 형태로 배설되는지, 명확하게 알려져 있는 약물은 그리 많지 않습니다. 희귀 질환의 경우에는 생명에 필수적인 역할을 하는 어떠한 유전자에 문제가 생겨서 바로 질환으로 나타나는 경우가 많습니다. 이는 해당 유전자가 만들어내는 단백질이 중요한 역할을 하고, 다른 유전자가 대신 기능을 해주지 못하기 때문입니다. 반면에 약물 유전자가 만들어내는 약물 효소의 종류는 워낙 다양해서 한 두가지 효소에 문제가 생긴다고 하더라도, 비슷한 다른 효소가 이러한 역할을 대신해주게 됩니다. 그리고 대사 경로 자체가 한가지 방향으로만 정해져 있는 것이 아니라, 어떠한 길이 막히면 다른 길로 돌아갈 수 있는 대체 경로가 존재하게 됩니다. 즉, 약물 대사능은 한가지 유전자와의 1:1 대응이 아니라, 다수의 효소들이 관여하여 복합적으로 나타나기 때문에 동시에 고려해야할 요소들이 많아지게 됩니다. 이를 유전학적으로 나타내보면 다음과 같습니다.

  • A number of isoforms (e.g. Cytochrome P450 family, GST family)
  • Many different transcription mode in a single gene: alternative splicing

 

IV. 연구 방법의 한계 : 유전적 바이오 마커 발굴의 연구 방법으로 많이 사용하고 있는 것이 SNP array chip 또는 NGS를 통한 시퀀싱입니다. SNP array는 주로 GWAS 연구에 사용하기 때문에 인구집단에 흔하게 존재하는 common variant 연구에 사용하고, NGS 시퀀싱은 유전자의 개별 변이까지 모두 확인하기 때문에 rare variant 발굴에 사용하게 됩니다. 그러나 두 연구 방법 모두 한계가 있습니다. 앞에서 언급한 것처럼 복합형질에서 common variant는 그 효과 크기에 대부분 매우 작기 때문에 GWAS 연구로는 새로운 마커의 발굴이 쉽지 않은 편입니다. 반면 Rare variant 발굴에 유리한 NGS 방법으로는 rare variant를 발굴하여도 그 변이의 해석이 쉽지 않고, 더불어 통계적으로 의미 있는 결과를 얻기 위해 필요한 n수가 매우 커서 현실적으로 연구가 어렵게 됩니다.

관련 포스팅 >

[유전자칩 비교] SNP array와 array CGH, 그리고 한국인칩

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

유전자 변이의 해석: 대용량 기능 검사의 필요성

위에서 언급한 여러가지 이유들로 인해, 약물 유전체 연구는 정말 어려운 분야인 것 같습니다. 하지만 다른 한편으로는 정밀의료 분야의 발전으로 가장 많은 사람들이 혜택을 볼 수 있는 분야도 약물과 관련된 분야이기 때문에, 그만큼 의미가 크다고 할 수 있겠습니다. 이러한 여러가지 어려운 점에도 불구하고, 열심히 노력하고 있는 연구자들이 함께 좋은 연구가 나올 수 있었으면 하는 바람입니다.