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

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

관련 포스팅 보기>

 

I. PLINK

1-s2.0-S0002929707613524-gr4_lrg

대부분의 GWASSNP 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

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

어제 정신과 전문의 친구와 점심을 먹었습니다. 제가 병원 연구실에서 유전체 연구를 하는 것을 듣고, 함께 연구할 아이디어에 대해서 이야기를 나누자고 만났는데, 안타깝게도 GWAS에 대한 개념이 전혀 없더군요. 지금은 바야흐로 GWAS의 시대입니다. 그래서 이번 글을 유전학 전공자가 아닌 일반인도 쉽게 이해할 수 있도록 전장 유전체 연관 분석 (Genome Wide Association Study; GWAS)의 개념과 연구 방법론에 대해서 글을 써 보고자 합니다.

저는 새로운 개념을 배울 때 항상 그 이름이 의미하는 바를 이해하려고 노력합니다. 그런 의미에서 GWAS라는 이름부터 파헤쳐보겠습니다.

Genome Wide = 전장 유전체 : 모든 유전체 위치에 대해서,

Association Study = 연관 분석: 관심을 가진 형질(Target phenotype)연관성을 갖는 유전적 위치를 찾는다.

GWAS
[GWAS 분석 방법의 개념] 일반적으로 Case (관심 형질을 가진 집단; 환자군)Control (형질을 갖지 않는 집단; 정상군)의 유전 정보를 서로 비교하여, case에서 더 많은 빈도를 갖는 = 연관성을 가진 유전자를 찾게 됩니다.
앞선 글에서 최근의 유전학 연구는 각 유전자 위치와 관련된 형질을 밝혀 그 발현 기전을 이해하는데 집중되고 있다는 말을 했습니다. GWAS는 그러한 유전자와 연관된 형질을 찾는 하나의 탐색 (Exploratory) 방법을 말합니다. 사실 무수히 많은 형질이 어떤 유전자와 관련되어 있는지 실험적으로 찾아내는 것은 정말로 어렵습니다. GWAS는 모든 유전자 위치에 대해 연관성의 정도를 분석하기 때문에, 관심있는 형질 또는 질환에 1차적으로 관련되어 있는 후보 유전자를 찾아내는 데 매우 유용한 탐색 도구 (screening method)가 됩니다.

GWAS는 일반적으로 Case (관심 형질을 가진 집단; 환자군)Control (형질을 갖지 않는 집단; 정상군)의 두 집단의 유전 정보를 얻은 후에 서로 비교하여, case에서 더 많은 빈도를 갖는, 즉 연관성을 가진 유전자를 찾게 됩니다. 한 가지 중요한 내용은 GWAS에서 찾아낸 유전자라 하더라도, 그것이 항상 원인 유전자는 아니라는 점입니다. 즉 GWAS는 인과 관계를 찾는 것이 아니라 우연히 연관되어 나타나는 유전자들의 후보를 찾는 과정입니다.

따라서 일반적으로 연구는 GWAS를 통한 후보 유전자 탐색 > 그리고 이 후에 더 많은 환자군에서 확인 (replication cohort) > 동물 & 세포 실험에서 생물학적 입증의 결과를 거쳐 최종적으로 유전자-형질의 관계를 밝히는 과정으로 진행됩니다. 이러한 GWAS 연구의 역사도 10년이 넘었습니다. GWAS는 강력한 툴 임에 틀림이 없지만, 그 원리가 통계적 연관성 분석에 기인하는 것이기 때문에 그 한계점도 분명히 인식하는 것이 중요합니다. 따라서 명확한 Case와 Control군을 확보하고, 통계적으로 분석이 가능한 충분한 수의 환자수를 확보하는 점도 중요하겠습니다. 하지만 현실에서는 이것이 쉽지만은 않죠.

linkage_disequilibrium
우리는 부모로 부터 두 쌍의 염색체 쌍 (상동 염색체)을 물려 받아 무작위적으로 재조합이 일어나게 됩니다. 그러나 유전자 재조합은 덩어리로 일어나기 때문에, 서로 거리가 가까운 유전자 위치 끼리는 유전형이 섞이지 않고 모자이크 패턴으로 함께 이동하게 되며, 이러한 하나의 덩어리를 일반적으로 LD block이라고 부릅니다.

더불어 GWAS 분석 방법을 이해하는 데 중요한 개념이 있는데, 흔히 LD (Linkage Disequilibrium)라고 부르는 ‘연관 비평형’ 입니다. 우리는 부모로부터 한 쌍씩 유전자를 물려받게 되는데, 생식 세포는 분열되면서 같은 세포 내에서도 끊임없이 유전형의 재배열이 일어납니다. 그러나 유전자 재조합은 덩어리로 일어나기 때문에, 서로 거리가 가까운 유전자 위치 끼리는 유전형이 섞이지 않고 모자이크 패턴으로 함께 이동하게 되며, 이러한 하나의 덩어리를 일반적으로 ‘LD block’이라고 부릅니다. 같은 LD block에 포함된 위치에 대해서는 연관성 분석을 하게 되면, 동일한 연관성을 보인 p 값을 보이기 됩니다. LD block의 존재는 다음과 같이 4가지를 시사합니다.

  1. GWAS 분석은 30억쌍의 모든 염기 서열에 대해서 할 필요가 없다. 같은 LD block에서 대표적인 하나의 마커만 이용해도 된다. > 분석 위치의 수가 축소화 됩니다.
  2. GWAS 연관 분석으로 후보 위치를 찾았다 하더라도, 정확한 원인 유전자의 위치는 LD block 내에 존재한 다른 위치일 수 있다. > GWAS로 찾아낸 후보 위치 근처의 유전형을 상세하게 다시 살펴봐야 하는 이유입니다.
  3. GWAS에 흔히 이용되는 Manhattan plot (맨하탄 플롯)에서 시그널이 하나의 탑처럼 주위에서 모두 높게 나오는 이유가 됩니다.
  4. 흔히 Imputation이라고 부르는 과정을 통해, 같은 LD block 내의 검사하지 않은 부위의 유전형도 추정이 가능해집니다.

Manhattan_Plot
[Manhattan plot] GWAS 분석 결과 의미 있는 시그널이 마치 맨하탄 가에 위치한 고층 빌딩들처럼 나온다고 하여 맨하탄 플롯 이라는 이름이 붙게 되었습니다.
마지막으로 GWAS에 관한 글은 GWAS catalog를 소개하면서 마치도록 하겠습니다. 지금 까지 무수히 많은 형질에 대한 GWAS 연구가 진행되었고, 최근에는 일반적인 형질에 대해서 UK biobank에 유전 정보와 형질이 공개되면서, 많은 부분 형질과 유전형 간의 GWAS 연구 및 관계가 드러나는 중입니다. 하지만 다시 한번 강조하면, 연관성과 인과 관계는 다릅니다. 따라서 확실한 생물학적 메카니즘으로 이를 설명하기 위해서는 후속 연구가 중요하게 됩니다. 이러한 GWAS 연구를 통해 형질과 유전자 위치의 관계가 명확하게 드러난 데이터를 모아 놓은 것이 GWAS catalog입니다. GWAS catalog는 지금도 계속 업데이트 되는 중이며, 나중에는 많은 질병과 유전병에 대해서 정보가 추가되기를 기대합니다.

아래 유튜브 자료에 GWAS catalog에 관한 내용이 잘 소개되어 있어 참고하면 좋을 것 같습니다.

[References]

Bush, William S., and Jason H. Moore. “Genome-wide association studies.” PLoS computational biology 8.12 (2012).

MacArthur, Jacqueline, et al. “The new NHGRI-EBI Catalog of published genome-wide association studies (GWAS Catalog).” Nucleic acids research 45.D1 (2017): D896-D901.

http://www.ebi.ac.uk/gwas/