관련 포스팅 보기>
Annovar: Population frequency, in-silico prediction tool 및 기타 database 활용
plink를 이용한 GWAS 분석에 대한 지난 포스팅에 이어서, 이번에는 SKAT을 이용한 rare variant 분석에 대한 포스팅을 정리해보고자 합니다. 유전자에 존재하는 변이(Variant)는 인구 집단 내의 분포 빈도 (Allele frequency)에 따라서, 흔한 변이 (Common variant)와 희귀 변이 (Rare variant)로 구분됩니다. 전장 유전체 연관 분석 (Genome-wide association study, GWAS)은 이러한 변이 중에서 일반적으로 인구 집단 내에 변이 빈도 분포가 5% 이상인 흔한 변이를 대상으로 하게 되는데, 대부분의 흔한 변이들은 유전학적 선택압 (Selective pressure)이 작은 변이들이기 때문에, 효과 크기 (Effect size)가 작은 것들이 대부분입니다.
따라서 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)
다만 Rare variant의 경우는, 변이 자체의 빈도가 매우 드물기 때문에 통계적으로 변이의 영향을 비교하기 위해서는 매우 큰 수의 표본이 필요한 한계가 있습니다. 이러한 한계를 극복하기 위해서, 같은 유전자 내에 존재하는 희귀 변이를 유전자 또는 특정 구역의 단위로 집합 시켜서 분석하는 방법이 제안되었는데, 이것이 Gene-level aggregation test 입니다. 즉, GWAS에서는 개별 SNP의 효과를 변이 단위로 분석이 진행되었다면 (Single variant association test), rare variant의 경우는 여러 개의 variant를 Gene 단위로 그룹화하여, Gene의 효과를 비교 분석 (Multiple variant association test)하는 분석을 진행하게 됩니다.
이와 같은 분석에서 중요한 점은 ‘어떠한 기준으로 변이를 유전자 단위로 그룹화 할 것인가?’ 하는 문제가 발생하게 되는데, 일반적으로 변이 빈도의 threshold 설정 (MAF 5% vs 1%), 변이의 Functional classification 등을 사용자가 임의로 설정하여 분석을 진행하게 됩니다. 더불어, 각 희귀 변이의 효과들이 동일하지 않기 때문에 개별 변이의 효과를 보정해주는 방법으로 Rare variant allele frequency, In-silico prediction score 등을 이용하게 됩니다. 아래는 이러한 희귀 변이 분석 방법과 각 분석 방법의 장, 단점, 그리고 분석 software를 정리한 표입니다.
위 표의 희귀 변이를 위한 다양한 분석 도구 중에서 오늘은 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 GitHub: https://github.com/leeshawn/SKAT
두마디 정밀의료 