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

광고

“plink를 이용한 GWAS 분석 및 Manhattan plot 만들기”의 9개의 생각

  1. 포스팅 감사합니다! 덕분에 많은 도움 받고 있습니다.
    질문이 있는데요. 맨하탄 플랏에서 유의한 SNP은 어떻게 추출하나요?
    위의 예에서는 data의 P값이 가장 작은 값이 y축에 높게 위치하는 것 같은데, min 값 등으로 추출하나요?

    Liked by 1명

    응답

    1. 아래 threshold 보다 작은 값 만족시키는 snp을 추출한 후에 개별 snp중에 생물학적으로 의미가 있는 녀석을 찾아내야합니다. false positive가 많기 때문이지요.

      좋아요

      응답

  2. 또 유의한 SNP의 기준은 어떻게 정하는지 알려주시면 감사하겠습니다! (그림에서 1개는 독보적이지만 다른 것들도 유의하다고 판단할 수 있는 기준 같은 것들이 있는지 궁금합니다)

    좋아요

    응답

    1. 유의한 P값의 threshold는 정의하기 나름입니다. 다만, 일반적으로 통계에서 0.05를 쓰는 것처럼 Genome-wide significant level 5*10^-8을 많이 씁니다. 0.05에 대한 multiple comparison correction을 bon ferroni 방법으로 적용시키면 대강 이러한 값이 나옵니다.

      Liked by 1명

      응답

      1. 아! 정말 감사합니다. 혼자 공부하고 있어서 긴가민가하는 부분이 많았는데, 많은 도움 되었습니다 ㅜㅜ

        한가지만 더 여쭤봐도 될까요?
        cohort study에서 심한 데이터 불균형 문제로, 기준에 따른 case-control로 나눠서 분석해 보고 있는데요.
        plink 문서에 따르면,
        PED file에서 affection column을 수정하거나 –pheno 옵션을 사용해서 association(chi-square) 분석을 할 수 있는 것으로 알고 있습니다. 저는 bfile로 바꾸고 –pheno 옵션을 사용해서 분석하고 있는데요.
        당연하게도 qq-plot을 그렸을 때, case 집단과 control집단, 합친 집단의 모양이 다르게 나오고 유의한 SNP도 다릅니다.

        –pheno 옵션을 사용해서 (1, 0으로 코딩된 phenotype.txt 파일 이용) 분석하는 경우,
        1) case+control이 합쳐진 fam 파일을 이용하는 것이 맞는지,
        2) case, control, case+control의 유의한 SNP이 다를 때, 영향을 준다고 해석할 수 있는 SNP은 case를 보는 것이 맞는지?
        의견을 부탁드려도 괜찮을까요?

        좋아요

      2. 디테일한 내용들은 데이터를 직접 핸들링하는 사람 외에는 다른 사람이 직접 데이터를 보지 않는한 설명만으로 알아듣기가 힘듭니다. 제가 느끼기에 질문 내용이 아직 정확하게 분석의 의도와 방향을 잘 모른다는 인상을 받습니다. 충분히 더 공부하시고 다양하게 분석을 하여 의미 있는 결과를 도출하시길 빌겠습니다.

        좋아요

  4. 안녕하세요? 간결하게 핵심작으로 잘 설명해주셔서 감사드립니다. 그런데
    “분석을 위한 Input Phenotype data를 준비하는 과정도 중요한데, 다음 페이지에서 자세하게 소개가 되어 있습니다” 라고 되어 있는데 다음페이지를 찾을수가 없네요. 혹시 가르쳐 주실 수 있나요?

    좋아요

    응답

답글 남기기

아래 항목을 채우거나 오른쪽 아이콘 중 하나를 클릭하여 로그 인 하세요:

Gravatar
WordPress.com 로고

WordPress.com의 계정을 사용하여 댓글을 남깁니다. ( 로그아웃 /  변경 )

Facebook 사진

Facebook의 계정을 사용하여 댓글을 남깁니다. ( 로그아웃 /  변경 )

취소

%s에 연결하는 중