【Linux】遺伝統計ソフトPLINKでGWAS

#はじめに

最近、遺伝統計ソフトのPLINKを使い始めた。
この本を参考にPLINKの使い方を勉強。

ゼロから実践する　遺伝統計学セミナー

この本はWindows用に書かれているので、Macでのやり方を備忘録的に書いていく。
今回は、PLINKでゲノムワイド関連解析（GWAS）を行う。

#####参考
PLINKの基本的な使い方は以前に投稿済み
【Linux】遺伝統計ソフトPLINKを使ってみた

#使うデータ

####bed|bim|fam形式ファイル

• SNP.bed
• SNP.bim
• SNP.fam

####表現型ファイル

• phenotype1.txt

このファイルは、ケースコントロールのデータファイルで、ケースは2、コントロールは1となっていて、それぞれにファミリーID、サンプルIDが割り当てられていて、SNP.famファイルのファミリーID、サンプルIDと紐付いている。

これらのファイルを作業ディレクトリに格納する。

#PLINKを起動

作業ディレクトリを指定して、

作業ディレクトリの指定

$ cd /作業ディレクトリのパス/

PLINK（PLINK実行ファイル）を作業ディレクトリに移動して起動する。

PLINK起動

$ ./plink

#SNPのフィルタリング

--bfileで、SNPファイル（bed|bim|fam形式）を読み込んで、--make-bedでフィルタリング後のデータを新しいbed|bim|fam形式ファイルとして作成。
--outで出力するファイルの名前をSNP_QCに指定する。

今回は以下の条件のSNPを除外する。
--maf 0.05：マイナーアレル頻度が5%以下（通常のGWASでは1%あるいは0.5%以下が一般的）
--hwe 0.000001：ハーディー・ワインベルグ平衡の検定でp値が10^-6以下
--indep-pairwise 100 5 0.8：連鎖平衡係数r2値が0.8以上

SNPの除去

$ ./plink --bfile SNP --out SNP_QC --maf 0.05 --hwe 0.000001 --indep-pairwise 100 5 0.8 --make-bed

実行結果

PLINK v1.90b6.16 64-bit (19 Feb 2020)          www.cog-genomics.org/plink/1.9/
(C) 2005-2020 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to SNP_QC.log.
Options in effect:
  --bfile SNP
  --hwe 0.000001
  --indep-pairwise 100 5 0.8
  --maf 0.05
  --make-bed
  --out SNP_QC

16384 MB RAM detected; reserving 8192 MB for main workspace.
8830185 variants loaded from .bim file.
381 people (178 males, 203 females) loaded from .fam.
381 phenotype values loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 381 founders and 0 nonfounders present.
Calculating allele frequencies... done.
--hwe: 14399 variants removed due to Hardy-Weinberg exact test.
2692226 variants removed due to minor allele threshold(s)
(--maf/--max-maf/--mac/--max-mac).
6123560 variants and 381 people pass filters and QC.
Among remaining phenotypes, 0 are cases and 381 are controls.
--make-bed to SNP_QC.bed + SNP_QC.bim + SNP_QC.fam ... done.
Pruned 366094 variants from chromosome 1, leaving 103214.
Pruned 403320 variants from chromosome 2, leaving 106161.
Pruned 339719 variants from chromosome 3, leaving 85861.
Pruned 353504 variants from chromosome 4, leaving 86916.
Pruned 310443 variants from chromosome 5, leaving 79003.
Pruned 322217 variants from chromosome 6, leaving 81854.
Pruned 285601 variants from chromosome 7, leaving 83237.
Pruned 264127 variants from chromosome 8, leaving 71155.
Pruned 208275 variants from chromosome 9, leaving 76275.
Pruned 245016 variants from chromosome 10, leaving 67270.
Pruned 241219 variants from chromosome 11, leaving 63223.
Pruned 226971 variants from chromosome 12, leaving 62436.
Pruned 177571 variants from chromosome 13, leaving 45339.
Pruned 154730 variants from chromosome 14, leaving 44160.
Pruned 134104 variants from chromosome 15, leaving 45299.
Pruned 142197 variants from chromosome 16, leaving 50213.
Pruned 124358 variants from chromosome 17, leaving 41954.
Pruned 134574 variants from chromosome 18, leaving 38435.
Pruned 105678 variants from chromosome 19, leaving 37131.
Pruned 102672 variants from chromosome 20, leaving 31013.
Pruned 68848 variants from chromosome 21, leaving 23361.
Pruned 63204 variants from chromosome 22, leaving 25608.
Pruning complete.  4774442 of 6123560 variants removed.
Marker lists written to SNP_QC.prune.in and SNP_QC.prune.out .

作業ディレクトリに以下のファイルが出力されているのを確認する。

• SNP_QC.bed
• SNP_QC.bim
• SNP_QC.fam

#ゲノムワイド関連解析（GWAS）

以下のコマンドを使ってGWAS解析を行う。
--pheno：GWASに使う表現型ファイル（今回はphenotype1.txt）を入力
--logistic：ロジスティック回帰を実施
--ci 0.95：95%信頼区間を出力

GWAS（ロジスティック回帰分析）

$ ./plink --bfile SNP_QC --out SNP_QC_Pheno1 --pheno phenotype1.txt --logistic --ci 0.95

GWASの実行結果

PLINK v1.90b6.16 64-bit (19 Feb 2020)          www.cog-genomics.org/plink/1.9/
(C) 2005-2020 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to SNP_QC_Pheno1.log.
Options in effect:
  --bfile SNP_QC
  --ci 0.95
  --logistic
  --out SNP_QC_Pheno1
  --pheno phenotype1.txt

16384 MB RAM detected; reserving 8192 MB for main workspace.
6123560 variants loaded from .bim file.
381 people (178 males, 203 females) loaded from .fam.
381 phenotype values present after --pheno.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 381 founders and 0 nonfounders present.
Calculating allele frequencies... done.
6123560 variants and 381 people pass filters and QC.
Among remaining phenotypes, 188 are cases and 193 are controls.
Writing logistic model association results to SNP_QC_Pheno1.assoc.logistic ...
done.

作業ディレクトリにSNP_QC_Pheno1.assoc.logisticという名前のファイルが出力されているのを確認して、テキストエディタで開く。
1列目が染色体番号、2列目がSNP ID、3列目が染色体位置、12列目がp値となっている。

#AWKによる要素の抽出

GWAS結果から、AWKコマンドで「染色体番号」、「SNP ID」、「染色体位置」、「p値」の列を取り出す。
AWKコマンドで、入力ファイルをSNP_QC_Pheno1.assoc.logistic、出力ファイルをテキストファイルとしたSNP_QC_Pheno1.assoc.logistic.P.txtとする。
AWKでは' 'で区切って、この中にコマンドを書いて実行する。
{print $2"\t"$1"\t"$3"\t"$12}によって、
「2列目[SNP ID]　1列目[染色体番号]　3列目[染色体位置]　4列目[p値]」というデータフレームになる。
"\t"はタブで区切るというコマンド。
>によって指定したテキストファイルとして出力する。

AWKコマンドでGWAS結果から要素を抽出し、テキストファイルを出力

$ awk '{print $2"\t"$1"\t"$3"\t"$12}' SNP_QC_Pheno1.assoc.logistic > SNP_QC_Pheno1.assoc.logistic.P.txt

出力ファイルをCSVファイルにすることもできる。

AWKコマンドでGWAS結果から要素を抽出し、CSVファイルを出力

$ awk '{print $2"\t"$1"\t"$3"\t"$12}' SNP_QC_Pheno1.assoc.logistic > SNP_QC_Pheno1.assoc.logistic.P.csv

このGWAS結果を使ってマンハッタンプロットを描く。
マンハッタンプロットの書き方は以前に投稿済みのため、今回は省略。

#####参考
Rでのマンハッタンプロットの書き方も以前に投稿済み
【R】GWAS結果でマンハッタンプロットを描いてみた
【R】GWAS結果でマンハッタンプロットを描いてみた2