はじめに
最近、遺伝統計ソフトのPLINKを使い始めた。
この本を参考にPLINKの使い方を勉強。
この本はWindows用に書かれているので、Macでのやり方を備忘録的に書いていく。
今回は、PLINKでeQTL解析をやってみた。
参考
PLINKの基本的な使い方、GWAS解析方法は以前に投稿済み
【Linux】遺伝統計ソフトPLINKを使ってみた
【Linux】遺伝統計ソフトPLINKでGWAS
使うデータ
bed|bim|fam形式ファイル
- SNP_QC.bed
- SNP_QC.bim
- SNP_QC.fam
SNPをフィルタリング済みのファイル
表現型ファイル
- Exp_BLK.txt
BLK遺伝子の発現量データ
eQTL解析
eQTL解析では、遺伝子発現量という量的形質に対して解析するため、線形解析を行う。
使うコマンド一覧。
--pheno：GWASに使う表現型ファイル（今回はExp_BLK.txt）を入力
--linear：線形回帰を実施
--ci 0.95：95%信頼区間を出力
$ ./plink --bfile SNP_QC --out SNP_QC_Exp_BLK --pheno Exp_BLK.txt --linear --ci 0.95
作業ディレクトリにSNP_QC_Exp_BLK.assoc.linearという名前のファイルが出力されているのを確認して、テキストエディタで開く。
1列目が染色体番号、2列目がSNP ID、3列目が染色体位置、12列目がp値となっている。
AWKによる要素の抽出
GWAS結果から、AWKコマンドで「染色体番号」、「SNP ID」、「染色体位置」、「p値」の列を取り出す。
AWKコマンドで、入力ファイルをSNP_QC_Exp_BLK.assoc.linear、出力ファイルをテキストファイルとしたSNP_QC_Exp_BLK.assoc.linear.P.txtとする。
AWKでは' 'で区切って、この中にコマンドを書いて実行する。
{print $2"\t"$1"\t"$3"\t"$12}によって、
「2列目[SNP ID] 1列目[染色体番号] 3列目[染色体位置] 4列目[p値]」というデータフレームになる。
"\t"はタブで区切るというコマンド。
>によって指定したテキストファイルとして出力する。
$ awk '{print $2"\t"$1"\t"$3"\t"$12}' SNP_QC_Exp_BLK.assoc.linear > SNP_QC_Exp_BLK.assoc.linear.P.txt
このGWAS結果を使ってマンハッタンプロットを描く。
マンハッタンプロットの描画
マンハッタンプロットを描いてみたらこうなった。
SNPの特定
eQTL効果を示したSNPを抜き出してみる。
awk '$4<=10^-12 {print $0}' 1KG_EUR_QC_Exp_BLK.assoc.linear.P.txt
出力
rs13255193 8 11309192 4.539e-13
rs13257831 8 11332964 6.545e-13
rs2736345 8 11352485 1.707e-14
rs1478898 8 11395079 6.882e-16
rs2244894 8 11448659 1.497e-13
rs2244648 8 11450422 2.068e-14
rs13273172 8 11461111 1.188e-14
p値が最小のSNPはrs1478898であった。