アメリエフの技術ブログ

Amelieff Staff Blog

PLINKの使い方1

これまでQTL解析のお話をしてきました。
今日はPLINKを使って基本的なSNPのデータ解析(GWAS:ゲノムワイド関連解析)を行います。

1.サンプルデータを用意します
NCBI Gene Expression Omnibus (GEO)に登録されているデータが整理されている「遺伝子発現バンク(GEO)目次」という大変便利なサイトがあります。
では、イルミナ社が公開しているHuman610-Quad v1.0ビーズチップでHapMapサンプルを解析したデータをダウンロードして使用します。
CHBとJPT(75サンプル)をコントロールCEU(73サンプル)をケースにします。

2.MAPファイルを作成します
MAPファイルの名前はhapmap.mapにしました。
マーカー数が約62万なので、62万行4列の縦長の形になります。

3.PEDファイルを作成します
まずasia.pedおよびceu.pedのふたつのPEDファイルを作成します。
affection statusカラムには、コントロールは「1」、ケースは「2」を入力します。つまり、asia.pedの6列目はすべて「1」、ceu.pedの6列目はすべて「2」となります。
ふたつのファイルを合体させると、縦がサンプル数・横がSNP数の、
148行62万列の横長の形になります。
(イルミナ社のゲノムブラウザーをご使用の方は、pedファイル作成のプラグインがございますのでご活用ください)

cat asia.ped ceu.ped > hapmap.ped


拡張子より前の部分はMAPファイルと同じ名前にする必要があるので、PEDファイルの名前はhapmap.pedにしました。

ふたつのファイルを作成するときに大切なのは、MAPファイルとPEDファイルのSNPの順番を一致させることです。

同じ要領で、HapMapのデータをコントロール、自分の実験データをケースとしてGWASをすることも可能です。

次のステップです。

4.MAPとPEDの書式を確認します
plink --file hapmap

画面に、マーカー数・サンプル数・ケースとコントロールそれぞれの数・性別などの情報が表示されますので、意図どおりに認識されているかチェックしてください。

5.PEDファイルとMAPファイルから、binary PEDファイルを作成します
plink --file hapmap --make-bed --out hapmap

6.いよいよ解析です
・Basic association analysis
plink --bfile hapmap --assoc --out as1

as1.assocというファイルが出来上がります。
カラム数は9つで
Chromosome
SNP identifier
Code for allele 1 (the minor, rare allele based on the entire sample frequencies)
The frequency of this variant in cases
The frequency of this variant in controls
Code for the other allele
The chi-squared statistic for this test (1 df)
The asymptotic significance value for this test
The odds ratio for this test
です。

・補正をするための「adjust」(ボンフェローニ補正など)
plink --bfile hapmap --assoc --adjust --out as2

カラム数は10で
Chromosome
SNP identifier
Unadjusted, asymptotic significance value
Genomic control adjusted significance value
Bonferroni adjusted significance value
Holm step-down adjusted significance value
Sidak single-step adjusted significance value
Sidak step-down adjusted significance value
Benjamini & Hochberg (1995) step-up FDR control
Benjamini & Yekutieli (2001) step-up FDR control
です。
補正前(as2.assoc)と補正後(as2.assoc.adjust)のふたつのファイルが作成されます。

【参考】
・遺伝子発現バンク(GEO)目次
http://lifesciencedb.jp/geo/

PLINK
http://pngu.mgh.harvard.edu/~purcell/plink/index.shtml