search
LoginSignup
2

More than 1 year has passed since last update.

posted at

updated at

【R】GWAS結果でマンハッタンプロットを描いてみた

はじめに

仕事で統計遺伝学が必要になりそうなので、GWAS結果からマンハッタンプロットを描いてみた。

備忘録的なやつ。

使用するデータ

NIHのGRASPデータベースに様々なGWAS結果が落ちている。

GRASP database

そのうちの一つを使用してみる。

Yang et al.

必要なパッケージ

qqmanパッケージとdata.tableパッケージをインストールしておく。

qqmanパッケージはマンハッタンプロットの作図に必要。

GWAS結果は一般にファイルサイズが大きいため、高速にデータを読み込むためfread()関数を使用する。

fread()関数はdata.tableパッケージに入っている。

library(qqman)
library(data.table)

データの読み込み

fread()関数でデータを読み込む。

DF <- fread("Yang_YoungOnsetHypertension_Illumina550_SBAS_rawpv.csv",
               stringsAsFactors = F, header = T)

str()関数でデータの内容を確認する。

実行結果。

> DF <- fread("Yang_YoungOnsetHypertension_Illumina550_SBAS_rawpv.csv",
+                stringsAsFactors = F, header = T)
> str(DF)
Classes data.table and 'data.frame':  475157 obs. of  8 variables:
 $ dbSNP_RS_ID      : chr  "rs3094315" "rs12562034" "rs3934834" "rs9442372" ...
 $ Chromosome       : int  1 1 1 1 1 1 1 1 1 1 ...
 $ Physical_position: int  742429 758311 995669 1008567 1011278 1020428 1021403 1038818 1039813 1051029 ...
 $ Symbol           : chr  "FLJ22639" "FLJ22639" "rs3934834" "C1orf159" ...
 $ CLR_N_pv         : num  0.977 0.706 0.651 0.673 0.731 ...
 $ CLR_N_BMI_pv     : num  0.357 0.298 0.348 0.844 0.429 ...
 $ CLR_C_pv         : num  0.692 0.438 0.513 1 0.527 ...
 $ CLR_C_BMI_pv     : num  0.114 0.0697 0.255 0.8723 0.2363 ...
 - attr(*, ".internal.selfref")=<externalptr>

dbSNP_RS_ID:SNPs ID
Chromosome:染色体番号
Physical_position:染色体位置
CLR_N_pvCLR_N_BMI_pvCLR_C_pvCLR_C_BMI_pv:p値(BMI調整あり/なし、条件付きロジスティック回帰)らしい...

マンハッタンプロットの描画

いよいよ描画。

qqmanパッケージのmanhattan()関数を使う。

#マンハッタンプロットの描画
manhattan(DF, chr = "Chromosome", 
          bp = "Physical_position", 
          p = "CLR_C_BMI_pv", 
          snp = "dbSNP_RS_ID") 

引数の説明。

chrには染色体番号のカラムを指定する。 X、Y、またはMT染色体がある場合は、これらの23、24、25などの番号に置換する。デフォルトはPLINKの「CHR」。

bpには染色体位置のカラムを指定する。デフォルトはPLINKの「BP」。この列は数値でなきゃダメ。

pはp値のカラム。デフォルトはPLINKの「P」。この列ももちろん数値でなきゃダメ。
※今回は「dbSNP_RS_ID」とした。

snpはSNP名のカラム(rs番号)を指定する。デフォルトはPLINKの「SNP」。このカラムは文字列でなきゃダメ。

描画結果

Rplot.png

青色のラインはp値が1.0×10-5のレベルであり、suggestive level(有意ではないが、なんらかの関連がありそうなレベル)を示している。

ここにはないけど、赤色のラインはp値が5.0×10-8のレベルであり、GWASでな有意レベルを示している。

最終スクリプト

#Rをキレイにしておく
rm(list = ls())

#ライブラリの読み込み
library(qqman)
library(data.table)

#データの読み込みと内容の確認
DF <- fread("Yang_YoungOnsetHypertension_Illumina550_SBAS_rawpv.csv",
               stringsAsFactors = F, header = T)
str(DF)

#マンハッタンプロットの描画
manhattan(DF, chr = "Chromosome", #染色体番号のカラムを指定する。
          bp = "Physical_position", #染色体位置のカラムを指定する。
          p = "CLR_C_BMI_pv", #p値のカラムを指定する。
          snp = "dbSNP_RS_ID") #SNP名のカラム(rs番号)を指定する。

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
What you can do with signing up
2